SND@LHC Software
Loading...
Searching...
No Matches
Functions | Variables
Survey-MufiScifi Namespace Reference

Functions

 pyExit ()
 
 map2Dict (aHit, T='GetAllSignals', mask=True)
 
 systemAndOrientation (s, plane)
 
 getAverageZpositions ()
 
 smallSiPMchannel (i)
 
 fit_langau (hist, o, bmin, bmax, opt='')
 
 twoLangaufun (x, par)
 
 langaufun (x, par)
 
 myPrint (tc, name, withRootFile=True)
 
 makeAnimation (histname, j0=1, j1=2, animated=True, findMinMax=True, lim=50)
 
 Scifi_hitMaps (Nev=options.nEvents)
 
 Mufi_hitMaps (Nev=options.nEvents)
 
 smallVsLargeSiPMs (Nev=-1)
 
 makeIndividualPlots (run=options.runNumber)
 
 makeQDCcorHTML (run=options.runNumber)
 
 makeLogVersion (run)
 
 eventTime (Nev=options.nEvents)
 
 TimeStudy (Nev=options.nEvents, withDisplay=False)
 
 beamSpot ()
 
 USshower (Nev=options.nEvents)
 
 USDSEnergy (Nev=options.nEvents, nproc=15)
 
 Mufi_Efficiency (Nev=options.nEvents, optionTrack=options.trackType, withReco='True', NbinsRes=100, X=10.)
 
 mips (readHists=True, option=0)
 
 mipsAfterBurner (readhisto=True)
 
 plotMips (readhisto=True)
 
 plotMipsGainRatio ()
 
 plotRMS (readHists=True, t0_channel=4)
 
 plotMipsTimeResT0 (readHists=True, animation=False)
 
 plotMipsTimeRes (readHists=True, animation=False)
 
 weightedY (aHist)
 
 makeDetPrints ()
 
 drawArea (s=3, p=0, opt='', color=ROOT.kGreen)
 
 plotSignalsUS ()
 
 analyze_EfficiencyAndResiduals (readHists=False, mode='S', local=True, zoom=False)
 
 TimeCalibrationNtuple (Nev=options.nEvents, nStart=0)
 
 analyzeTimings (c='', vsignal=15., args={}, opt=0)
 
 scifiTimings ()
 
 ScifiChannel_residuals ()
 
 SciFiPropSpeed ()
 
 Scifi_xPos (detID)
 
 Scifi_slopes (Nev=options.nEvents)
 
 mergeScifiSignals (hstore)
 
 mergeMuFilterSignals (hstore, tag='signalT', system='', fname=None)
 
 signalZoom (smax)
 
 scifi_beamSpot ()
 
 TwoTrackFinder (nstart=0, Nev=-1, sMin=10, dClMin=7, minDistance=1.5, sepDistance=0.5, debug=False)
 
 Scifi_residuals (Nev=options.nEvents, NbinsRes=100, xmin=-2000., alignPar=False, unbiased=True, minEnergy=False, remakeClusters=options.remakeScifiClusters, nproc=1)
 
 printScifi_residuals (tag='v0')
 
 printScifiAlignment ()
 
 minimizeAlignScifi (first=True, level=1, migrad=False)
 
 FCN (npar, gin, f, par, iflag)
 
 plotsForCollabMeeting ()
 
 testReversChannelMapping (p=None)
 

Variables

 Tkey = ROOT.std.vector('TString')()
 
 Ikey = ROOT.std.vector('int')()
 
 Value = ROOT.std.vector('float')()
 
dict h = {}
 
 parser = ArgumentParser()
 
 dest
 
 help
 
 type
 
 int
 
 required
 
 False
 
 default
 
 str
 
 None
 
 options = parser.parse_args()
 
 runNr = str(options.runNumber).zfill(6)
 
 path = options.path
 
int partitions = 0
 
 dirlist = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls "+options.path,shell=True) )
 
str data = "sndsw_raw_"+runNr+".root"
 
 geo
 
 MuFilter = geo.modules['MuFilter']
 
 Scifi = geo.modules['Scifi']
 
 nav = ROOT.gGeoManager.GetCurrentNavigator()
 
 nScifi = Scifi.GetConfParI("Scifi/nscifi")
 
 nMats = Scifi.GetConfParI("Scifi/nmats")
 
 nVeto = MuFilter.GetConfParI("MuFilter/NVetoPlanes")
 
 A
 
 B
 
 locA
 
 locB
 
 globA
 
 globB
 
 latex = ROOT.TLatex()
 
dict systemAndPlanes
 
dict systemAndBars
 
dict systemAndChannels
 
dict sdict = {1:'Veto',2:'US',3:'DS'}
 
 freq = u.snd_freq
 
 TDC2ns = u.snd_TDC2ns
 
 zPos = getAverageZpositions()
 
 eventChain = ROOT.TChain('rawConv')
 
 f = ROOT.TFile.Open(options.fname)
 
 rc = eventChain.GetEvent(0)
 
 run = ROOT.FairRunAna()
 
 ioman = ROOT.FairRootManager.Instance()
 
 outFile = ROOT.TMemFile('dummy','CREATE')
 
 source = ROOT.FairFileSource(eventChain.GetCurrentFile())
 
 sink = ROOT.FairRootFileSink(outFile)
 
bool houghTransform = False
 
 muon_reco_task = SndlhcMuonReco.MuonReco()
 
 trackTask = SndlhcTracking.Tracking()
 
 xrdb = ROOT.FairRuntimeDb.instance()
 
 eventTree
 
 OT = ioman.GetSink().GetOutTree()
 
 Reco_MuonTracks = trackTask.fittedTracks
 
 Cluster_Scifi = trackTask.clusScifi
 
list barColor
 
 tmp = options.command.split(';')
 
str command = tmp[0]+"("
 

Function Documentation

◆ analyze_EfficiencyAndResiduals()

Survey-MufiScifi.analyze_EfficiencyAndResiduals (   readHists = False,
  mode = 'S',
  local = True,
  zoom = False 
)

Definition at line 2670 of file Survey-MufiScifi.py.

2670def analyze_EfficiencyAndResiduals(readHists=False,mode='S',local=True,zoom=False):
2671 if readHists: ut.readHists(h,'MuFilterEff_run'+str(options.runNumber)+'.root',withProj=False)
2672# start with tracks
2673 ut.bookCanvas(h,'Tracks','',900,1600,1,2)
2674 h['Tracks'].cd(1)
2675 h['tracksChi2Ndof'].Draw('colz')
2676 h['Tracks'].cd(2)
2677 h['trackslxy'].Draw('colz')
2678
2679 # params for the double symmetrical Crystal Ball - symetric core and tails
2680 a1 = ROOT.RooRealVar("alpha", "alpha", 2, 1e-3, 10.)
2681 p1 = ROOT.RooRealVar("n", "n", 1, 1e-3, 10.)
2682
2683 if local: res = 'res'
2684 else: res='gres'
2685 for proj in ['X','Y']:
2686 if mode == "S":
2687 if proj=="Y": ut.bookCanvas(h,'EVetoUS'+proj,'',1200,2000,2,7)
2688 if proj=="Y": ut.bookCanvas(h,'EDS'+proj,'',1400,2100,2,3)
2689 if proj=="X":
2690 ut.bookCanvas(h,'EDS'+proj,'',1000,2000,2,4)
2691 ut.bookCanvas(h,'EVeto'+proj,'',800,400,2,1)
2692 else:
2693 ut.bookCanvas(h,'EVetoUS'+proj,'',1200,2000,4,8)
2694 ut.bookCanvas(h,'EDS'+proj,'',1200,2000,3,7)
2695 k = 1
2696 residualsAndSigma = {}
2697 for s in [1,2,3]:
2698 if mode == "S" and s==2 and proj=="X": continue
2699 if s==3: k=1
2700 sy = sdict[s]
2701 if s<3: tname = 'EVetoUS'+proj # changed for veto3 in the plane loop below
2702 else: tname = 'EDS'+proj
2703 for plane in range(systemAndPlanes[s]):
2704 if mode == "S" and s==1 and proj=="X" and plane<2: continue
2705 if s==1 and plane==2:
2706 if proj=="Y": continue# vertical Veto3
2707 tname = 'EVeto'+proj
2708 if s==3 and proj=="X" and (plane%2==0 and plane<5): continue
2709 if s==3 and proj=="Y" and (plane%2==1 or plane==6): continue
2710 tc = h[tname].cd(k)
2711 key = sy+str(s*10+plane)
2712 h['track_' + key].Draw('colz')
2713 k+=1
2714 for side in ['L','R','S']:
2715 if side=='S' and (s==3 or mode=='S') : continue
2716 if side=='R' and (mode=='S') : continue
2717 tc = h[tname].cd(k)
2718 k+=1
2719 if mode=="S":
2720 hname = res+proj+'_'+key
2721 else:
2722 hname = res+proj+'_'+sy+side+str(s*10+plane)
2723 resH = h[hname].ProjectionX(hname+'_proj')
2724 if zoom:
2725 if s==3: resH.GetXaxis().SetRangeUser(-4.,4.)
2726 else: resH.GetXaxis().SetRangeUser(-10.,10.)
2727 resH.Draw()
2728 binw = resH.GetBinWidth(1)
2729 if resH.GetEntries()>10:
2730 #fit using a double sided CB shape with symmetric tails
2731 x = ROOT.RooRealVar("x", "residual [cm]", -10., 10.)
2732 mu = ROOT.RooRealVar("mu", "mu", resH.GetMean(), resH.GetMean()-1, resH.GetMean()+1)
2733 width = ROOT.RooRealVar("width", "width", resH.GetRMS(), 1e-5, resH.GetRMS()+2)
2734 dcbPdf = ROOT.RooCrystalBall("dcbPdf", "DoubleSidedCB", x, mu, width, a1, p1, True)
2735 roofit_hist = ROOT.RooDataHist("roofit_hist", "roofit_hist", ROOT.RooArgList(x), ROOT.RooFit.Import(resH))
2736 fitResult = dcbPdf.fitTo(roofit_hist, ROOT.RooFit.PrintLevel(-1), ROOT.RooFit.Save(True))
2737 pl = x.frame()
2738 roofit_hist.plotOn(pl)
2739 dcbPdf.plotOn(pl)
2740 pl.Draw()
2741 latex.DrawLatexNDC(0.13,0.8, "#mu = {:.3f} +/- {:.3f}".format(mu.getVal(), mu.getError()))
2742 latex.DrawLatexNDC(0.13,0.7,"#sigma = {:.3f} +/- {:.3f}".format(width.getVal(), width.getError()))
2743 '''
2744 # changed to CB fit function so commenting the modified Gaussian fit
2745 # Keeping it however since Thomas was able to get resolutions ~ bar_size/sqrt(12) in 2022.
2746 # For the 2024 data that is not the case, so swapped with CB
2747 myGauss = ROOT.TF1('gauss','abs([0])*'+str(binw)+'/(abs([2])*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+abs([3])',4)
2748 myGauss.SetParameter(0,resH.GetEntries())
2749 myGauss.SetParameter(1,0)
2750 myGauss.SetParameter(2,2.)
2751 rc = resH.Fit(myGauss,'SL','',-15.,15.)
2752 fitResult = rc.Get()
2753 '''
2754 if not (s*10+plane) in residualsAndSigma:
2755 residualsAndSigma[s*10+plane] = [mu.getVal()/u.mm,mu.getError()/u.mm,abs(width.getVal()/u.cm),width.getError()/u.cm]
2756 #residualsAndSigma[s*10+plane] = [fitResult.Parameter(1)/u.mm,fitResult.ParError(1)/u.mm,abs(fitResult.Parameter(2)/u.cm),fitResult.ParError(2)/u.cm]
2757 tc.Update()
2758 '''
2759 stats = resH.FindObject('stats')
2760 stats.SetOptFit(111111)
2761 stats.SetX1NDC(0.63)
2762 stats.SetY1NDC(0.25)
2763 stats.SetX2NDC(0.98)
2764 stats.SetY2NDC(0.94)
2765 tracks = h['track_'+key].GetEntries()
2766 if tracks >0:
2767 eff = fitResult.Parameter(0)/tracks
2768 effErr = fitResult.ParError(0)/tracks
2769 latex.DrawLatexNDC(0.2,0.8,'eff=%5.2F+/-%5.2F%%'%(eff,effErr))
2770 '''
2771 if 'EVetoUS'+proj in h: myPrint(h['EVetoUS'+proj],'EVetoUS'+proj)
2772 if 'EVeto'+proj in h: myPrint(h['EVeto'+proj],'EVeto'+proj)
2773 if 'EDS'+proj in h: myPrint(h['EDS'+proj],'EDS'+proj)
2774# summary plot of residuals
2775 h['gResiduals'], h['gSigma'] = ROOT.TGraphErrors(), ROOT.TGraphErrors()
2776 k = 0
2777 for sp in residualsAndSigma:
2778 h['gResiduals'].SetPoint(k,sp,residualsAndSigma[sp][0])
2779 h['gResiduals'].SetPointError(k,0.5,residualsAndSigma[sp][1])
2780 h['gSigma'].SetPoint(k,sp,residualsAndSigma[sp][2])
2781 h['gSigma'].SetPointError(k,0.5,residualsAndSigma[sp][3])
2782 k+=1
2783 h['gResiduals'].SetLineWidth(2)
2784 h['gSigma'].SetLineWidth(2)
2785 ut.bookCanvas(h,'MufiRes'+proj,'Mufi residuals',750,750,1,1)
2786 tc = h['MufiRes'+proj].cd()
2787 h['gResiduals'].SetTitle(';station; offset [mm]')
2788 h['gResiduals'].Draw("AP")
2789 ut.bookCanvas(h,'MufiSigm','Mufi sigma',750,750,1,1)
2790 tc = h['MufiSigm'].cd()
2791 h['gSigma'].SetTitle(';station; #sigma [cm]')
2792 h['gSigma'].Draw("AP")
2793 myPrint(h['MufiRes'+proj],'MufiResiduals'+proj)
2794 myPrint(h['MufiSigm'],'MufiSigma'+proj)
2795
2796 ut.bookCanvas(h,'TVE','',800,900,1,nVeto)
2797 ut.bookCanvas(h,'TUS','',800,2000,1,5)
2798 ut.bookCanvas(h,'TDS','',800,1200,1,3)
2799 latex.SetTextColor(ROOT.kRed)
2800
2801 S = {1:'TVE',2:'TUS',3:'TDS'}
2802 for s in S:
2803 sy = sdict[s]
2804 tname = S[s]
2805 k=1
2806 for plane in range(systemAndPlanes[s]):
2807 if s==1 and proj=="X" and plane<2: continue
2808 if s==1 and plane==2:
2809 if proj=="Y": continue# vertical Veto3
2810 if s==3 and proj=="X" and (plane%2==0 and plane<5): continue
2811 if s==3 and proj=="Y" and (plane%2==1 or plane==6): continue
2812 tc = h[tname].cd(k)
2813 k+=1
2814 key = sy+str(s*10+plane)
2815 hist = h['dtLRvsX_'+key]
2816 if hist.GetSumOfWeights()==0: continue
2817# find beam spot
2818 tmp = h['track_'+key].ProjectionX()
2819 rc = tmp.Fit('gaus','S')
2820 res = rc.Get()
2821 fmin = res.Parameter(1) - 3*res.Parameter(2)
2822 fmax = res.Parameter(1) + 3*res.Parameter(2)
2823 hist.SetStats(0)
2824 xmin = max( hist.GetXaxis().GetBinCenter(1), fmin)
2825 xmax = min( hist.GetXaxis().GetBinCenter(hist.GetNbinsX()), fmax)
2826 hist.GetXaxis().SetRangeUser(xmin,xmax)
2827 hist.GetYaxis().SetRange(1,hist.GetNbinsY())
2828 rc = hist.ProjectionY().Fit('gaus','S')
2829 res = rc.Get()
2830 ymin = max( hist.GetYaxis().GetBinCenter(1), -5*res.Parameter(2))
2831 ymax = min( hist.GetYaxis().GetBinCenter(hist.GetNbinsY()), 5*res.Parameter(2))
2832 hist.GetYaxis().SetRangeUser(ymin,ymax)
2833 hist.Draw('colz')
2834# get time x correlation, X = m*dt + b
2835 h['gdtLRvsX_'+key] = ROOT.TGraph()
2836 g = h['gdtLRvsX_'+key]
2837 xproj = hist.ProjectionX('tmpx')
2838 if xproj.GetSumOfWeights()==0: continue
2839 np = 0
2840 Lmin = hist.GetSumOfWeights() * 0.005
2841 for nx in range(hist.GetXaxis().FindBin(xmin),hist.GetXaxis().FindBin(xmax)):
2842 tmp = hist.ProjectionY('tmp',nx,nx)
2843 if tmp.GetEntries()<Lmin:continue
2844 X = hist.GetXaxis().GetBinCenter(nx)
2845 rc = tmp.Fit('gaus','NQS')
2846 res = rc.Get()
2847 if not res: dt = tmp.GetMean()
2848 else: dt = res.Parameter(1)
2849 g.SetPoint(np,X,dt)
2850 np+=1
2851 g.SetLineColor(ROOT.kRed)
2852 g.SetLineWidth(2)
2853 g.Draw('same')
2854 rc = g.Fit('pol1','SQ','',xmin,xmax)
2855 g.GetFunction('pol1').SetLineColor(ROOT.kGreen)
2856 g.GetFunction('pol1').SetLineWidth(2)
2857 result = rc.Get()
2858 if not result: continue
2859 m = 1./result.Parameter(1)
2860 b = -result.Parameter(0) * m
2861 sign = '+'
2862 if b<0: sign = '-'
2863 txt = 'X = %5.1F #frac{cm}{ns} #times #frac{1}{2}#Deltat %s %5.1F '%(m*2,sign,abs(b))
2864 latex.SetTextSize(0.07)
2865 latex.DrawLatexNDC(0.2,0.8,txt)
2866 myPrint(h['TVE'],'dTvsX_Veto')
2867 myPrint(h['TUS'],'dTvsX_US')
2868 myPrint(h['TDS'],'dTvsX_DS')
2869
2870# timing relative to Scifi tracks
2871 colors = {1:ROOT.kGreen,2:ROOT.kRed,3:ROOT.kBlue}
2872 ut.bookCanvas(h,'relT','',800,1000,1,7)
2873 for s in systemAndPlanes:
2874 for plane in range(systemAndPlanes[s]):
2875 hname = 'atLRvsX_'+sdict[s]+str(s*10+plane)
2876 h[hname+'_proj'] = h[hname].ProjectionY(hname+'_proj')
2877 h[hname+'_proj'].SetLineColor(colors[s])
2878 h[hname+'_proj'].SetStats(0)
2879 h[hname+'_proj'].SetTitle('; #Deltat [ns]')
2880 if s==1 and plane == 0: h[hname+'_proj'] .Draw()
2881 else: h[hname+'_proj'] .Draw('same')
2882
2883
2884#for VETO
2885 ut.bookCanvas(h,'veto','',800,1600,1,7)
2886 s=1
2887 for plane in range(systemAndPlanes[s]):
2888 if plane<2: hname = 'resBarY_'+sdict[s]+str(s*10+plane)
2889 elif plane==2: hname = 'resBarX_'+sdict[s]+str(s*10+plane)
2890 for nbar in range(1,h[hname].GetNbinsY()+1):
2891 vname = 'veto'+str(s*10+plane)+'_'+str(nbar)
2892 h[vname] = h[hname].ProjectionX(vname,nbar+1,nbar+1)
2893 binw = h[vname].GetBinWidth(1)
2894 myGauss = ROOT.TF1('gauss','abs([0])*'+str(binw)+'/(abs([2])*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+abs([3])',4)
2895 myGauss.SetParameter(0,resH.GetEntries())
2896 myGauss.SetParameter(1,0)
2897 myGauss.SetParameter(2,2.)
2898 rc = h[vname].Fit(myGauss,'SL','',-15.,15.)
2899 fitResult = rc.Get()
2900
2901# Scifi specific code

◆ analyzeTimings()

Survey-MufiScifi.analyzeTimings (   c = '',
  vsignal = 15.,
  args = {},
  opt = 0 
)

Definition at line 2969 of file Survey-MufiScifi.py.

2969def analyzeTimings(c='',vsignal = 15., args={},opt=0):
2970 print('-->',c)
2971 if c=='':
2972 analyzeTimings(c='findMax',vsignal = vsignal)
2973 new_list = sorted(h['mult'].items(), key=lambda x: x[1], reverse=True)
2974 args = {"detID0":new_list[0][0]}
2975 print(args)
2976 analyzeTimings(c='firstIteration',vsignal = vsignal, args=args,opt=0)
2977 h['dTtwin'].Draw()
2978 return
2979 fNtuple = ROOT.TFile("scifi_timing_"+str(options.runNumber)+".root")
2980 tTimes = fNtuple.tTimes
2981 if c=="findMax":
2982 h['mult'] = {}
2983 for nt in tTimes:
2984 for nHit in range(nt.nChannels):
2985 detID = nt.detIDs[nHit]
2986 if not detID in h['mult']: h['mult'][detID] = 0
2987 h['mult'][detID] += 1
2988 new_list = sorted(h['mult'].items(), key=lambda x: x[1], reverse=True)
2989 detID0 = new_list[0][0] # args = {"detID0":new_list[0][0]}
2990 if c=="firstIteration":
2991 ut.bookHist(h,'dTtwin','dt between neighbours;[ns]',100,-5,5)
2992 ut.bookHist(h,'dTtwinMax','dt between neighbours; [ns]',100,-5,5)
2993 for x in h:
2994 if x[0]=='c':h[x].Reset()
2995
2996 for nt in tTimes:
2997 tag = False
2998 for nHit in range(nt.nChannels):
2999 if nt.detIDs[nHit] == args["detID0"]:
3000 tag = nHit
3001 break
3002 if not tag: continue
3003 for nHit in range(nt.nChannels):
3004 detID = nt.detIDs[nHit]
3005 dLL = (nt.dL[nHit] - nt.dL[tag]) / vsignal # light propagation in cm/ns
3006 dtdc = nt.tdc[nHit] - nt.tdc[tag]
3007 if opt>0:
3008 hname = 'c'+str(detID)
3009 if not hname in h: ut.bookHist(h,hname,";dt [ns]",200,-10.,10.)
3010 rc = h[hname].Fill(dtdc-dLL)
3011 for nHit2 in range(nHit+1,nt.nChannels):
3012 detID2 = nt.detIDs[nHit2]
3013 if abs(detID-detID2)==1:
3014 dt = nt.tdc[nHit] - nt.tdc[nHit2]
3015 if detID2>detID: dt=-dt
3016 if nt.qdc[nHit] < 15 or nt.qdc[nHit2] < 15: continue
3017 rc = h['dTtwin'].Fill(dt)
3018 if detID == args["detID0"]: rc = h['dTtwinMax'].Fill(dt)
3019 ut.bookCanvas(h,'scifidT','',1200,900,1,1)
3020 tc = h['scifidT'].cd()
3021 h['dTtwin'].Draw()
3022

◆ beamSpot()

Survey-MufiScifi.beamSpot ( )

Definition at line 1041 of file Survey-MufiScifi.py.

1041def beamSpot():
1042 trackTask.ExecuteTask()
1043 Xbar = -10
1044 Ybar = -10
1045 for aTrack in Reco_MuonTracks:
1046 state = aTrack.getFittedState()
1047 pos = state.getPos()
1048 rc = h['bs'].Fill(pos.x(),pos.y())
1049 points = aTrack.getPoints()
1050 keys = ROOT.std.vector('int')()
1051 detIDs = ROOT.std.vector('int')()
1052 ROOT.fixRoot(points, detIDs,keys,True)
1053 for k in range(keys.size()):
1054 # m = p.getRawMeasurement()
1055 detID =detIDs[k] # m.getDetId()
1056 key = keys[k] # m.getHitId()//1000 # for mufi
1057 aHit = eventTree.Digi_MuFilterHits[key]
1058 if aHit.GetDetectorID() != detID: continue # not a Mufi hit
1059 s = detID//10000
1060 l = (detID%10000)//1000 # plane number
1061 bar = (detID%1000)
1062 if s>2:
1063 l=2*l
1064 if bar>59:
1065 bar=bar-60
1066 if l<6: l+=1
1067 if s==3 and l%2==0: Ybar=bar
1068 if s==3 and l%2==1: Xbar=bar
1069 nSiPMs = aHit.GetnSiPMs()
1070 nSides = aHit.GetnSides()
1071 for p in range(nSides):
1072 c=bar*nSiPMs*nSides + p*nSiPMs
1073 for i in range(nSiPMs):
1074 signal = aHit.GetSignal(i+p*nSiPMs)
1075 if signal > 0:
1076 rc = h['Tsig_'+str(s)+str(l)].Fill(signal)
1077 mom = state.getMom()
1078 slopeY= mom.X()/mom.Z()
1079 slopeX= mom.Y()/mom.Z()
1080 h['slopes'].Fill(slopeX,slopeY)
1081 if not Ybar<0 and not Xbar<0 and abs(slopeY)<0.01: rc = h['bsDS'].Fill(Xbar,Ybar)
1082

◆ drawArea()

Survey-MufiScifi.drawArea (   s = 3,
  p = 0,
  opt = '',
  color = ROOT.kGreen 
)

Definition at line 2598 of file Survey-MufiScifi.py.

2598def drawArea(s=3,p=0,opt='',color=ROOT.kGreen):
2599 ut.bookCanvas(h,'area','',900,900,1,1)
2600 hname = 'track_'+sdict[s]+str(s)+str(p)
2601 barw = 1.0
2602 if s<3: barw=3.
2603 h[hname].SetStats(0)
2604 h[hname].SetTitle('')
2605 h[hname].Draw('colz'+opt)
2606 mufi = geo.modules['MuFilter']
2607 if opt=='': h['lines'] = []
2608 lines = h['lines']
2609 latex.SetTextColor(ROOT.kRed)
2610 latex.SetTextSize(0.03)
2611 for bar in range(systemAndBars[s]):
2612 mufi.GetPosition(s*10000+p*1000+bar,A,B)
2613 barName = str(bar)
2614 botL = ROOT.TVector3(A.X(),A.Y()-barw,0)
2615 botR = ROOT.TVector3(B.X(),B.Y()-barw,0)
2616 topL = ROOT.TVector3(A.X(),A.Y()+barw,0)
2617 topR = ROOT.TVector3(B.X(),B.Y()+barw,0)
2618 lines.append(ROOT.TLine(topL.X(),topL.Y(),topR.X(),topR.Y()))
2619 lines.append(ROOT.TLine(botL.X(),botL.Y(),botR.X(),botR.Y()))
2620 lines.append(ROOT.TLine(botL.X(),botL.Y(),topL.X(),topL.Y()))
2621 lines.append(ROOT.TLine(botR.X(),botR.Y(),topR.X(),topR.Y()))
2622 if s<3 or (s==3 and bar%5==0): latex.DrawLatex(botR.X()-5.,botR.Y()+0.3,barName)
2623 N = len(lines)
2624 for n in range(N-4,N):
2625 l=lines[n]
2626 l.SetLineColor(color)
2627 l.SetLineWidth(4)
2628 for l in lines: l.Draw('same')
2629 latex.DrawLatexNDC(0.2,0.2,"Right")
2630 latex.DrawLatexNDC(0.8,0.2,"Left")
2631 if s<3: return
2632 h['linesV'] = []
2633 linesV = h['linesV']
2634 for bar in range(systemAndBars[s]):
2635 mufi.GetPosition(s*10000+(p+1)*1000+bar+60,A,B)
2636 topL = ROOT.TVector3(A.X()+barw,A.Y(),0)
2637 botL = ROOT.TVector3(B.X()+barw,B.Y(),0)
2638 topR = ROOT.TVector3(A.X()-barw,A.Y(),0)
2639 botR = ROOT.TVector3(B.X()-barw,B.Y(),0)
2640 linesV.append(ROOT.TLine(topL.X(),topL.Y(),topR.X(),topR.Y()))
2641 linesV.append(ROOT.TLine(topR.X(),topR.Y(),botR.X(),botR.Y()))
2642 linesV.append(ROOT.TLine(botR.X(),botR.Y(),botL.X(),botL.Y()))
2643 linesV.append(ROOT.TLine(botL.X(),botL.Y(),topL.X(),topL.Y()))
2644 barName = str(bar)
2645 if bar%5==0: latex.DrawLatex(topR.X(),topR.Y()+2,barName)
2646 for l in linesV:
2647 l.SetLineColor(ROOT.kRed)
2648 l.SetLineWidth(4)
2649 l.Draw('same')
2650

◆ eventTime()

Survey-MufiScifi.eventTime (   Nev = options.nEvents)

Definition at line 899 of file Survey-MufiScifi.py.

899def eventTime(Nev=options.nEvents):
900 if Nev < 0 : Nev = eventTree.GetEntries()
901 ut.bookHist(h,'Etime','delta event time; dt [s]',100,0.0,1.)
902 ut.bookHist(h,'EtimeZ','delta event time; dt [ns]',1000,0.0,10000.)
903 ut.bookCanvas(h,'T',' ',1024,3*768,1,3)
904
905# need to make extra gymnastiques since absolute time is missing
906 Ntinter = []
907 N = 0
908 for f in eventTree.GetListOfFiles():
909 dN = f.GetEntries()
910 rc = eventTree.GetEvent(N)
911 t0 = eventTree.EventHeader.GetEventTime()/freq
912 rc = eventTree.GetEvent(N+dN-1)
913 tmax = eventTree.EventHeader.GetEventTime()/freq
914 Ntinter.append([t0,tmax])
915 N+=dN
916
917 Tduration = 0
918 for x in Ntinter:
919 Tduration += (x[1]-x[0])
920 tsep = 3600.
921 t0 = Ntinter[0][0]
922 tmax = Tduration+(tsep*(len(Ntinter)-1))
923 nbins = 1000
924 yunit = "events per %5.0F s"%( (tmax-t0)/nbins)
925 if 'time' in h: h.pop('time').Delete()
926 ut.bookHist(h,'time','elapsed time; t [s];'+yunit,nbins,0,tmax-t0)
927
928 N=-1
929 Tprev = 0
930 Toffset = 0
931 for event in eventTree:
932 N+=1
933 if N>Nev: break
934 T = event.EventHeader.GetEventTime()
935 dT = T-Tprev
936 if N>0 and dT >0:
937 rc = h['Etime'].Fill( dT/freq )
938 rc = h['EtimeZ'].Fill( dT*1E9/freq )
939 rc = h['time'].Fill( (T+Toffset)/freq-t0 )
940 elif dT<0:
941 Toffset+=tsep*freq+Tprev
942 rc = h['time'].Fill( (T+Toffset)/freq-t0 )
943 else: rc = h['time'].Fill( T/freq-t0 ) # very first event
944 Tprev = T
945
946 tc = h['T'].cd(1)
947 h['time'].SetStats(0)
948 h['time'].Draw()
949 tend = 0
950 for x in Ntinter:
951 tend += x[1]+tsep/2.
952 m = str(int(tend))
953 h['line'+m]=ROOT.TLine(tend,0,tend,h['time'].GetMaximum())
954 h['line'+m].SetLineColor(ROOT.kRed)
955 h['line'+m].Draw()
956 tend += tsep/2.
957 tc = h['T'].cd(2)
958 tc.SetLogy(1)
959 h['EtimeZ'].Draw()
960 rc = h['EtimeZ'].Fit('expo','S','',0.,250.)
961 h['T'].Update()
962 stats = h['EtimeZ'].FindObject('stats')
963 stats.SetOptFit(1111111)
964 tc = h['T'].cd(3)
965 tc.SetLogy(1)
966 h['Etime'].Draw()
967 rc = h['Etime'].Fit('expo','S')
968 h['T'].Update()
969 stats = h['Etime'].FindObject('stats')
970 stats.SetOptFit(1111111)
971 h['T'].Update()
972 myPrint(h['T'],'time')
973

◆ FCN()

Survey-MufiScifi.FCN (   npar,
  gin,
  f,
  par,
  iflag 
)

Definition at line 3959 of file Survey-MufiScifi.py.

3959def FCN(npar, gin, f, par, iflag):
3960#calculate chisquare
3961 h['iter']+=1
3962 chisq = 0
3963 h['alignPar'] = {}
3964 theArray = [h['iter'],0]
3965 for p in range(h['npar']):
3966 if h['npar'] == 10:
3967 s = p//2+1
3968 o = p%2
3969 name = "Scifi/LocD"+str(s*10+o)
3970 if h['npar'] > 29 and h['npar'] != 32:
3971 if p<30:
3972 s = p//6+1
3973 o = (p%6)//3
3974 m = p%3
3975 name = "Scifi/LocM"+str(s*100+o*10+m)
3976 else:
3977 s = (p-30)//3 + 1
3978 if p%3 == 0: name = "Scifi/RotPhiS"+str(s)
3979 if p%3 == 1: name = "Scifi/RotPsiS"+str(s)
3980 if p%3 == 2: name = "Scifi/RotThetaS"+str(s)
3981 if h['npar'] == 32:
3982 if p<8:
3983 s = p//2+1
3984 o = p%2
3985 m = 0
3986 name = "Scifi/LocM"+str(s*100+o*10+m)
3987 else:
3988 s = (p-8)//6 + 1
3989 o = ((p-8)//3)%2
3990 if p%3 == 2: name = "Scifi/RotPhiS"+str(s)+str(o)
3991 if p%3 == 0: name = "Scifi/RotPsiS"+str(s)+str(o)
3992 if p%3 == 1: name = "Scifi/RotThetaS"+str(s)+str(o)
3993 h['alignPar'][name] = par[p]
3994 print('minuit %s %6.4F'%(name,par[p]), p)
3995 if h['level']==1 and p<nScifi*nMat*2: # station alignment
3996 if m>0: h['alignPar'][name] = par[p-m]
3997 else:
3998 h['alignPar'][name] = par[p]
3999 theArray.append(par[p])
4000 X = Scifi_residuals(Nev=h['Nevents'],NbinsRes=100,xmin=h['xmin'],alignPar=h['alignPar'],nproc=10)
4001 projs = {1:'V',0:'H'}
4002 for s in range(1,5):
4003 sigma = 15
4004 if s==1 or s==nScifi: sigma = 35
4005 for o in range(2):
4006 for p in projs:
4007 proj = projs[p]
4008 for x in ['X','Y']:
4009 print('res'+x+proj+'_Scifi'+str(s*10+o))
4010 hist = h['res'+x+proj+'_Scifi'+str(s*10+o)]
4011 nbins = hist.GetNbinsY()
4012 for k in range(1,nbins+1):
4013 tmp = hist.ProjectionX('tmp',k,k)
4014 chisq += (tmp.GetMean()/sigma)**2
4015 print('chisq=',chisq,iflag,h['iter'])
4016 f.value = chisq
4017 theArray[1]=chisq
4018 theTuple = array('f',theArray)
4019 h['evol'].Fill(theTuple)
4020 return

◆ fit_langau()

Survey-MufiScifi.fit_langau (   hist,
  o,
  bmin,
  bmax,
  opt = '' 
)

Definition at line 195 of file Survey-MufiScifi.py.

195def fit_langau(hist,o,bmin,bmax,opt=''):
196 if opt == 2:
197 params = {0:'Width(scale)',1:'mostProbable',2:'norm',3:'sigma',4:'N2'}
198 F = ROOT.TF1('langau',langaufun,0,200,len(params))
199 else:
200 params = {0:'Width(scale)',1:'mostProbable',2:'norm',3:'sigma'}
201 F = ROOT.TF1('langau',twoLangaufun,0,200,len(params))
202 for p in params: F.SetParName(p,params[p])
203 rc = hist.Fit('landau','S'+o,'',bmin,bmax)
204 res = rc.Get()
205 if not res: return res
206 F.SetParameter(2,res.Parameter(0))
207 F.SetParameter(1,res.Parameter(1))
208 F.SetParameter(0,res.Parameter(2))
209 F.SetParameter(3,res.Parameter(2))
210 F.SetParameter(4,0)
211 F.SetParLimits(0,0,100)
212 F.SetParLimits(1,0,100)
213 F.SetParLimits(3,0,10)
214 rc = hist.Fit(F,'S'+o,'',bmin,bmax)
215 res = rc.Get()
216 return res
217

◆ getAverageZpositions()

Survey-MufiScifi.getAverageZpositions ( )

Definition at line 168 of file Survey-MufiScifi.py.

168def getAverageZpositions():
169 zPos={'MuFilter':{},'Scifi':{}}
170 for s in systemAndPlanes:
171 for plane in range(systemAndPlanes[s]):
172 bar = 4
173 p = plane
174 if s==3 and (plane%2==0 or plane==7):
175 bar = 90
176 p = plane//2
177 elif s==3 and plane%2==1:
178 bar = 30
179 p = plane//2
180 MuFilter.GetPosition(s*10000+p*1000+bar,A,B)
181 zPos['MuFilter'][s*10+plane] = (A.Z()+B.Z())/2.
182 for s in range(1,nScifi+1):
183 mat = 1
184 sipm = 1
185 channel = 64
186 for o in range(2):
187 Scifi.GetPosition(channel+1000*sipm+10000*mat+100000*o+1000000*s,A,B)
188 zPos['Scifi'][s*10+o] = (A.Z()+B.Z())/2.
189 return zPos
190
MuFilter::GetPosition
void GetPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition MuFilter.cxx:639
Scifi::GetPosition
void GetPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition Scifi.cxx:664

◆ langaufun()

Survey-MufiScifi.langaufun (   x,
  par 
)

Definition at line 224 of file Survey-MufiScifi.py.

224def langaufun(x,par):
225 #Fit parameters:
226 #par[0]=Width (scale) parameter of Landau density
227 #par[1]=Most Probable (MP, location) parameter of Landau density
228 #par[2]=Total area (integral -inf to inf, normalization constant)
229 #par[3]=Width (sigma) of convoluted Gaussian function
230 #
231 #In the Landau distribution (represented by the CERNLIB approximation),
232 #the maximum is located at x=-0.22278298 with the location parameter=0.
233 #This shift is corrected within this function, so that the actual
234 #maximum is identical to the MP parameter.
235#
236 # Numeric constants
237 invsq2pi = 0.3989422804014 # (2 pi)^(-1/2)
238 mpshift = -0.22278298 # Landau maximum location
239#
240 # Control constants
241 np = 100.0 # number of convolution steps
242 sc = 5.0 # convolution extends to +-sc Gaussian sigmas
243#
244 # Variables
245 summe = 0.0
246#
247 # MP shift correction
248 mpc = par[1] - mpshift * par[0]
249#
250 # Range of convolution integral
251 xlow = max(0,x[0] - sc * par[3])
252 xupp = x[0] + sc * par[3]
253#
254 step = (xupp-xlow) / np
255#
256 # Convolution integral of Landau and Gaussian by sum
257 i=1.0
258 if par[0]==0 or par[3]==0: return 9999
259 while i<=np/2:
260 i+=1
261 xx = xlow + (i-.5) * step
262 fland = ROOT.TMath.Landau(xx,mpc,par[0]) / par[0]
263 summe += fland * ROOT.TMath.Gaus(x[0],xx,par[3])
264#
265 xx = xupp - (i-.5) * step
266 fland = ROOT.TMath.Landau(xx,mpc,par[0]) / par[0]
267 summe += fland * ROOT.TMath.Gaus(x[0],xx,par[3])
268#
269 return (par[2] * step * summe * invsq2pi / par[3])
270

◆ makeAnimation()

Survey-MufiScifi.makeAnimation (   histname,
  j0 = 1,
  j1 = 2,
  animated = True,
  findMinMax = True,
  lim = 50 
)

Definition at line 277 of file Survey-MufiScifi.py.

277def makeAnimation(histname,j0=1,j1=2,animated=True, findMinMax=True, lim = 50):
278 ut.bookCanvas(h,'ani','',900,800,1,1)
279 tc = h['ani'].cd()
280 jmin,jmax = j0,j1
281 if findMinMax:
282 jmin,jmax = 0,0
283 for j in range(j0,j1):
284 tmp = histname.replace('$$$',str(j))
285 if tmp in h:
286 if h[tmp].GetEntries()>lim:
287 jmin = j
288 print(j,tmp,h[tmp].GetEntries())
289 break
290 for j in range(j1,j0,-1):
291 tmp = histname.replace('$$$',str(j))
292 if tmp in h:
293 if h[tmp].GetEntries()>lim:
294 jmax = j
295 break
296 for j in range(jmin,jmax):
297 tmp = histname.replace('$$$',str(j))
298 h[tmp].Draw()
299 tc.Update()
300 stats = h[tmp].FindObject('stats')
301 stats.SetOptFit(1111111)
302 h[tmp].Draw()
303 if animated:
304 h['ani'].Print('picAni'+str(j)+'.png')
305 else:
306 rc = input("hit return for next event or q for quit: ")
307 if rc=='q': break
308 if animated and jmax > jmin:
309 os.system("convert -delay 60 -loop 0 picAni*.png "+histname+".gif")
310 os.system('rm picAni*.png')
311
312
313
314# initialize
315

◆ makeDetPrints()

Survey-MufiScifi.makeDetPrints ( )

Definition at line 2594 of file Survey-MufiScifi.py.

2594def makeDetPrints():
2595 for s in range(1,4):
2596 drawArea(s,p=0,opt='',color=ROOT.kGreen)
2597 myPrint(h['area'],'areaDet_'+sdict[s])

◆ makeIndividualPlots()

Survey-MufiScifi.makeIndividualPlots (   run = options.runNumber)

Definition at line 844 of file Survey-MufiScifi.py.

844def makeIndividualPlots(run=options.runNumber):
845 ut.bookCanvas(h,'dummy','',900,800,1,1)
846 if not "run"+str(run) in os.listdir(): os.system("mkdir run"+str(run))
847 for l in range(5):
848 for side in ['L','R']:
849 f=ROOT.TFile('QDCcor'+side+str(l)+'-run'+str(run)+'.root')
850 tcanv = f.Get('cor'+side+str(l))
851 for pad in tcanv.GetListOfPrimitives():
852 if not hasattr(pad,"GetListOfPrimitives"): continue
853 for aHist in pad.GetListOfPrimitives():
854 if not aHist.ClassName() == 'TH2D': continue
855 hname = aHist.GetName()
856 tmp = hname.split('_')
857 bar = tmp[1][2]
858 pname = 'corUS'+str(l)+'-'+str(bar)+side+'_'+tmp[0][3:]
859 aHist.SetDirectory(ROOT.gROOT)
860 ROOT.gROOT.cd()
861 tc=h['dummy'].cd()
862 tc.SetLogz(1)
863 aHist.Draw('colz')
864 tc.Update()
865 stats = aHist.FindObject('stats')
866 stats.SetOptStat(11)
867 stats.SetX1NDC(0.15)
868 stats.SetY1NDC(0.75)
869 stats.SetX2NDC(0.35)
870 stats.SetY2NDC(0.88)
871 h['dummy'].Update()
872 tc.Print('run'+str(run)+'/'+pname+'.png')
873 tc.Print('run'+str(run)+'/'+pname+'.root')
874 #os.system("convert -delay 120 -loop 0 run"+str(run)+"/corUS*.png corUS-"+str(run)+".gif")
875

◆ makeLogVersion()

Survey-MufiScifi.makeLogVersion (   run)

Definition at line 886 of file Survey-MufiScifi.py.

886def makeLogVersion(run):
887 for l in range(5):
888 for side in ['L','R']:
889 fname = 'QDCcor'+side+str(l)+'-run'+str(run)
890 f=ROOT.TFile('QDCcor'+side+str(l)+'-run'+str(run)+'.root')
891 c = 'cor'+side+str(l)
892 h['X'] = f.Get(c).Clone(c)
893 for pad in h['X'].GetListOfPrimitives():
894 pad.SetLogz(1)
895 h['X'].Draw()
896 h['X'].Print(fname+'.pdf')
897
898

◆ makeQDCcorHTML()

Survey-MufiScifi.makeQDCcorHTML (   run = options.runNumber)

Definition at line 876 of file Survey-MufiScifi.py.

876def makeQDCcorHTML(run=options.runNumber):
877 F = ROOT.TFile.Open('QDCcorrelations-run'+str(run)+'.root','recreate')
878 for l in range(5):
879 for side in ['L','R']:
880 key = side+str(l)
881 f=ROOT.TFile('QDCcor'+key+'-run'+str(run)+'.root')
882 tcanv = f.Get('cor'+key).Clone()
883 F.mkdir(key)
884 F.cd(key)
885 tc.Write()

◆ map2Dict()

Survey-MufiScifi.map2Dict (   aHit,
  T = 'GetAllSignals',
  mask = True 
)

Definition at line 67 of file Survey-MufiScifi.py.

67def map2Dict(aHit,T='GetAllSignals',mask=True):
68 if T=="SumOfSignals":
69 key = Tkey
70 elif T=="GetAllSignals" or T=="GetAllTimes":
71 key = Ikey
72 else:
73 print('use case not known',T)
74 1/0
75 key.clear()
76 Value.clear()
77 if T=="GetAllTimes": ROOT.fixRootT(aHit,key,Value,mask)
78 else: ROOT.fixRoot(aHit,key,Value,mask)
79 theDict = {}
80 for k in range(key.size()):
81 if T=="SumOfSignals": theDict[key[k].Data()] = Value[k]
82 else: theDict[key[k]] = Value[k]
83 return theDict
84

◆ mergeMuFilterSignals()

Survey-MufiScifi.mergeMuFilterSignals (   hstore,
  tag = 'signalT',
  system = '',
  fname = None 
)

Definition at line 3271 of file Survey-MufiScifi.py.

3271def mergeMuFilterSignals(hstore,tag='signalT',system='',fname=None):
3272 if fname: F=ROOT.TFile(fname)
3273 ROOT.gROOT.cd()
3274 if system == '': s = [1,2,3]
3275 else: s = [system]
3276 for x in s:
3277 for side in ['L','R']:
3278 hname = tag+side+'_'+sdict[x]
3279 xname = tag+side+'_'+sdict[x]+str(x*10)+'_'+'0'
3280 if F:
3281 hstore[xname] = F.Get(xname).Clone(xname)
3282 hstore[hname] = hstore[tag+side+'_'+sdict[x]+str(x*10)+'_'+'0'].Clone(hname)
3283 hstore[hname].Reset()
3284 hstore[hname].SetTitle('QDC '+sdict[x]+'; QDC [a.u.]; d [cm]')
3285 for l in range(systemAndPlanes[x]):
3286 for bar in range(systemAndBars[x]):
3287 xname = tag+side+'_'+sdict[x]+str(x*10+l)+'_'+str(bar)
3288 if F:
3289 hstore[xname] = F.Get(xname).Clone(xname)
3290 hstore[hname].Add(hstore[xname])
3291 print('summary histogram:',hname)
3292

◆ mergeScifiSignals()

Survey-MufiScifi.mergeScifiSignals (   hstore)

Definition at line 3265 of file Survey-MufiScifi.py.

3265def mergeScifiSignals(hstore):
3266 ut.bookHist(hstore,'signalAll','signal all mat',150,0.0,150.)
3267 for mat in range(30):
3268 hstore['signalAll'].Add(hstore['sig_'+str(mat)])
3269 hstore['signalAll'].Scale(1./hstore['signalAll'].GetSumOfWeights())
3270

◆ minimizeAlignScifi()

Survey-MufiScifi.minimizeAlignScifi (   first = True,
  level = 1,
  migrad = False 
)

Definition at line 3744 of file Survey-MufiScifi.py.

3744def minimizeAlignScifi(first=True,level=1,migrad=False):
3745 eventTree.SetBranchStatus("Cluster_Scifi",0)
3746 global trackTask
3747 trackTask = SndlhcTracking.Tracking()
3748 trackTask.SetName('simpleTracking')
3749 trackTask.Init()
3750 trackTask.makeScifiClusters = True
3751 trackTask.clusScifi = ROOT.TObjArray(100);
3752
3753 npar = nScifi*nMats*2 + nScifi*3
3754 if nMats ==1: npar += nScifi*3
3755 vstart = array('d',[0]*npar)
3756 h['Nevents'] = 500000
3757 if first:
3758 h['xmin'] =-5000.
3759 X = Scifi_residuals(Nev=10000,NbinsRes=100,xmin=h['xmin'])
3760 for m in range(nMats):
3761 vstart[m] = -X["Scifi/LocD10"][0]
3762 vstart[3+m] = -X["Scifi/LocD11"][0]
3763 err = 1000.
3764 h['npar'] = 10
3765 else:
3766 if level==1:
3767 h['xmin'] =-2000.
3768 h['npar'] = 30 + 15
3769 h['level'] = 1 # station alignment
3770
3771 # take relative mat alignment from H6
3772
3773 vstart[0] = 0 # s1
3774 vstart[1] = 0
3775 vstart[2] = 0
3776 vstart[3] = 0
3777 vstart[4] = -44.7
3778 vstart[5] = 0
3779 vstart[6] = 0 # s2
3780 vstart[7] = 0
3781 vstart[8] = 0
3782 vstart[9] = 0
3783 vstart[10] = 0
3784 vstart[11] = 0.
3785 vstart[12] = 0. # s3
3786 vstart[13] = 0.
3787 vstart[14] = 0.
3788 vstart[15] = 0.
3789 vstart[16] = -92.4
3790 vstart[17] = 0
3791 vstart[18] = 0 # s3
3792 vstart[19] = 28.0
3793 vstart[20] = 168.9
3794 vstart[21] = 0.
3795 vstart[22] = 25.5
3796 vstart[23] = -21.0
3797 vstart[24] = 0 # s4
3798 vstart[25] = -46.7
3799 vstart[26] = 46.9
3800 vstart[27] = 0
3801 vstart[28] = 523
3802 vstart[29] = 169.6
3803# rotation, three angles / station
3804 vstart[30] = 0.0000 # s1
3805 vstart[31] = 0.0
3806 vstart[32] = 0.0
3807 vstart[33] = 0.0000 # s2
3808 vstart[34] = 0.0
3809 vstart[35] = 0.0000
3810 vstart[36] = 0.0000 # s3
3811 vstart[37] = 0.0
3812 vstart[38] = 0.0000
3813 vstart[39] = 0.0000 # s4
3814 vstart[40] = 0.0
3815 vstart[41] = 0.0000
3816 vstart[42] = 0.0000 # s5
3817 vstart[43] = 0.0
3818 vstart[44] = 0.0000
3819
3820 err = 500.
3821 if level==2:
3822 h['level'] = 2 # mat alignment
3823 for i in range (npar):
3824 vstart[i]=0.0
3825
3826 err = 20.
3827 h['xmin'] =-2000.
3828 h['npar'] = npar
3829 if level==3:
3830 h['Nevents'] = 100000
3831 p = 0
3832 for s in range(1, nScifi+1):
3833 for o in range(2):
3834 for m in range(nMats):
3835 x = "Scifi/LocM"+str(s*100+o*10+m)
3836 vstart[p] = Scifi.GetConfParF(x)
3837 p+=1
3838 err = 25.
3839 h['xmin'] =-2000.
3840 h['npar'] = 30
3841
3842 ierflg = ctypes.c_int(0)
3843 gMinuit = ROOT.TMinuit(npar) # https://root.cern.ch/download/minuit.pdf
3844 gMinuit.SetFCN(FCN)
3845 gMinuit.SetErrorDef(1.0)
3846
3847 p=0
3848 variables = "n:chi2:"
3849 if nMats ==1 : plane = ['0', '1']
3850 else: plane = ['']
3851 for s in range(1, nScifi+1):
3852 for o in range(2):
3853 name = "Scifi/LocM"
3854 for m in range(nMats):
3855 gMinuit.mnparm(p, name+str(s*100+o*10+m), vstart[p], err, 0.,0.,ierflg)
3856 variables +=name.replace('/','')+str(s*100+o*10+m)+":"
3857 p+=1
3858 for s in range(1, nScifi+1):
3859 for pln in plane:
3860 for r in ["RotPhiS","RotPsiS", "RotThetaS"]:
3861 name = "Scifi/"+r+str(s)+pln
3862 gMinuit.mnparm(p, name, vstart[p], err/100, 0.,0.,ierflg)
3863 variables +=name.replace('/','')+":"
3864 p+=1
3865 if level == 1:
3866 p=0
3867 for s in range(1, nScifi+1):
3868 for o in range(2):
3869 for m in range(nMats):
3870 if s==1 or s==5 : gMinuit.FixParameter(p)
3871 elif m==1 or m==2 : gMinuit.FixParameter(p)
3872 p+=1
3873 for s in range(1, nScifi+1):
3874 for a in range(nMats):
3875 gMinuit.FixParameter(p)
3876 p+=1
3877
3878 if level == 2:
3879 p=0
3880 for s in range(1,nScifi+1):
3881 for o in range(2):
3882 for m in range(nMats):
3883 if nMats == 1:
3884 # fix one plane in one station
3885 # since for some stations in the H8 2023 testbeam
3886 # there are rotations btw planes that are accounted for
3887 if s==2 and o==0 : gMinuit.FixParameter(p)
3888 else:
3889 if s==3 or s==4 : gMinuit.FixParameter(p)
3890 p+=1
3891 for s in range(1, nScifi+1):
3892 if nMats == 1:
3893 for o in range(2):
3894 for a in range(3):
3895 if a==0 or a==2 or (s==2 and o==0): gMinuit.FixParameter(p)
3896 p+=1
3897 else:
3898 for a in range(3):
3899 if a==0 or a==2 or s==3 : gMinuit.FixParameter(p)
3900 p+=1
3901
3902 h['evol'] = ROOT.TNtuple("nt","evolution",variables[0:len(variables)-2])
3903 if level == 3:
3904# station 1 H
3905 gMinuit.FixParameter(0)
3906 gMinuit.FixParameter(1)
3907 gMinuit.FixParameter(2)
3908# station 1 V
3909 gMinuit.FixParameter(3)
3910 #gMinuit.FixParameter(4)
3911 gMinuit.FixParameter(5)
3912# station 2 H
3913 gMinuit.FixParameter(6)
3914 gMinuit.FixParameter(7)
3915 gMinuit.FixParameter(8)
3916# station 2 V
3917 gMinuit.FixParameter(9)
3918 gMinuit.FixParameter(10)
3919 gMinuit.FixParameter(11)
3920# station 3 H
3921 gMinuit.FixParameter(12)
3922 gMinuit.FixParameter(13)
3923 gMinuit.FixParameter(14)
3924# station 3 V
3925 gMinuit.FixParameter(15)
3926 #gMinuit.FixParameter(16)
3927 gMinuit.FixParameter(17)
3928# station 4 H
3929 #gMinuit.FixParameter(18)
3930 gMinuit.FixParameter(19)
3931 #gMinuit.FixParameter(20)
3932# station 4 V
3933 #gMinuit.FixParameter(21)
3934 #gMinuit.FixParameter(22)
3935 #gMinuit.FixParameter(23)
3936# station 5 H
3937 gMinuit.FixParameter(24)
3938 #gMinuit.FixParameter(25)
3939 #gMinuit.FixParameter(26)
3940# station 5 V
3941 #gMinuit.FixParameter(27)
3942 #gMinuit.FixParameter(28)
3943 #gMinuit.FixParameter(29)
3944
3945 h['iter'] = 0
3946 strat = array('d',[0])
3947 gMinuit.mnexcm("SET STR",strat,1,ierflg) # 0 faster, 2 more reliable
3948 gMinuit.mnexcm("SIMPLEX",vstart,npar,ierflg)
3949 if migrad: gMinuit.mnexcm("MIGRAD",vstart,npar,ierflg)
3950
3951 p = 'C2'
3952 ut.bookCanvas(h,p,p,750,750,1,1)
3953 tc = h[p].cd()
3954 h['evol'].Draw("chi2:n",'','*')
3955
3956 h[p].SaveAs("chi2ndf.root","root")
3957 h[p].SaveAs("chi2ndf.pdf","pdf")
3958
Scifi::GetConfParF
Float_t GetConfParF(TString name)
Definition Scifi.h:50

◆ mips()

Survey-MufiScifi.mips (   readHists = True,
  option = 0 
)

Definition at line 1732 of file Survey-MufiScifi.py.

1732def mips(readHists=True,option=0):
1733# plot mean sipm channel fired vs X
1734 if readHists:
1735 for x in h:
1736 if hasattr(x,'Reset'): x.Reset()
1737 ut.readHists(h,'MuFilterEff_run'+str(options.runNumber)+'.root',withProj=False)
1738 s = 2
1739 for plane in range(5):
1740 for bar in range(10):
1741 for p in ['L','R']:
1742 for T in ['T','']:
1743 if T=='T': nloop = 1
1744 else: nloop = systemAndChannels[s][1]+systemAndChannels[s][0]
1745 for i in range(nloop):
1746 if s==2 and smallSiPMchannel(i): continue
1747 name = 'signal'+T+p+'_US'+str(s*10+plane)+'_'+str(bar)
1748 if nloop>1: name += '-c'+str(i)
1749 histo = h[name]
1750 for x in ['M','C','W','S','E']:
1751 h[x+name]=histo.ProjectionY(x+name)
1752 h[x+name].Reset()
1753 h[x+name].SetTitle(histo.GetName()+';distance [cm]')
1754 for n in range(1,h['M'+name].GetNbinsX()+1):
1755 h[name+'-X'+str(n)] = h[name].ProjectionX(name+'-X'+str(n),n,n)
1756 tmp = h[name+'-X'+str(n)]
1757 h['C'+name].SetBinContent(n,-1)
1758 h['E'+name].SetBinContent(n,tmp.GetEntries())
1759 if tmp.GetEntries()>50:
1760 print('Fit ',options.runNumber,name,n)
1761 if T=='T': bmin,bmax = 10,tmp.GetMaximumBin()
1762 else: bmin,bmax = 0,30
1763 for k in range(tmp.GetMaximumBin(),1,-1):
1764 if tmp.GetBinContent(k)<2:
1765 bmin = k
1766 break
1767 for k in range(tmp.GetMaximumBin(),tmp.GetNbinsX()):
1768 if tmp.GetBinContent(k) + tmp.GetBinContent(k+1)<2:
1769 bmax = k
1770 break
1771 res = fit_langau(tmp,'LQ',0.8*tmp.GetBinCenter(bmin),1.5*tmp.GetBinCenter(bmax))
1772 if not res: continue
1773 if not res.IsValid(): continue
1774 h['C'+name].SetBinContent(n,res.Chi2()/res.Ndf())
1775 h['M'+name].SetBinContent(n,res.Parameter(1))
1776 h['M'+name].SetBinError(n,res.ParError(1))
1777 h['W'+name].SetBinContent(n,res.Parameter(0))
1778 h['W'+name].SetBinError(n,res.ParError(0))
1779 h['S'+name].SetBinContent(n,res.Parameter(3))
1780 h['S'+name].SetBinError(n,res.ParError(3))
1781 name = 'nSiPMs'+p+'_US'+str(s*10+plane)+'_'+str(bar)
1782 histo = h[name]
1783 x='N'
1784 h[x+name]=histo.ProjectionY(x+name)
1785 h[x+name].Reset()
1786 h[x+name].SetTitle(histo.GetName()+';distance [cm]')
1787 for n in range(1,h[x+name].GetNbinsX()+1):
1788 tmp = h[name].ProjectionX('tmp',n,n)
1789 if tmp.GetEntries()<10: continue
1790 h[x+name].SetBinContent(n,tmp.GetMean())
1791 h[x+name].SetBinError(n,tmp.GetRMS())
1792
1793# make DS
1794# add all bars in one plane
1795 s=3
1796 for plane in range(4):
1797 for p in ['L','R']:
1798 name = 'signal'+p+'_DS'+str(s*10+plane)
1799 h[name]=h[name+'_'+str(0)].Clone(name)
1800 for bar in range(60):
1801 h[name].Add(h[name+'_'+str(bar)])
1802#
1803 histo = h[name]
1804 for x in ['M','W','S']:
1805 h[x+name]=histo.ProjectionY(x+name)
1806 h[x+name].Reset()
1807 h[x+name].SetTitle(histo.GetName()+';distance [cm]')
1808 for n in range(1,h['M'+name].GetNbinsX()+1):
1809 h[name+'-X'+str(n)] = h[name].ProjectionX(name+'-X'+str(n),n,n)
1810 tmp = h[name+'-X'+str(n)]
1811 if tmp.GetEntries()>50:
1812 print('Fit ',options.runNumber,name,n)
1813 bmin,bmax = 0,80
1814 for k in range(tmp.GetMaximumBin(),1,-1):
1815 if tmp.GetBinContent(k)<2:
1816 bmin = k
1817 break
1818 res = fit_langau(tmp,'LQ',0.8*tmp.GetBinCenter(bmin),1.5*tmp.GetBinCenter(bmax))
1819 if not res: continue
1820 if not res.IsValid(): continue
1821 h['M'+name].SetBinContent(n,res.Parameter(1))
1822 h['M'+name].SetBinError(n,res.ParError(1))
1823 h['W'+name].SetBinContent(n,res.Parameter(0))
1824 h['W'+name].SetBinError(n,res.ParError(0))
1825 h['S'+name].SetBinContent(n,res.Parameter(3))
1826 h['S'+name].SetBinError(n,res.ParError(3))
1827
1828 ut.writeHists(h,'LandauFits_run'+str(options.runNumber)+'.root')

◆ mipsAfterBurner()

Survey-MufiScifi.mipsAfterBurner (   readhisto = True)

Definition at line 1829 of file Survey-MufiScifi.py.

1829def mipsAfterBurner(readhisto=True):
1830 if readhisto: ut.readHists(h,'LandauFits_run'+str(options.runNumber)+'.root',withProj=False)
1831 s=2
1832 for plane in range(5):
1833 for p in ['L','R']:
1834 for bar in range(10):
1835 name = 'signalT'+p+'_US'+str(s*10+plane)+'_'+str(bar)
1836 for nb in range(1,h['M'+name].GetNbinsX()+1):
1837 tmp = h[name+'-X'+str(nb)]
1838 N = tmp.GetEntries()
1839 if N>50 and h['M'+name].GetBinContent(nb)==0:
1840 fg = tmp.GetFunction('langau')
1841 bmin,bmax = fg.GetXmin(),fg.GetXmax()
1842 res = fit_langau(tmp,'LQ',bmin,bmax)
1843 if res.IsValid():
1844 print(name,nb,N,'fixed')
1845 h['M'+name].SetBinContent(nb,res.Parameter(1))
1846 h['M'+name].SetBinError(nb,res.ParError(1))
1847 h['W'+name].SetBinContent(nb,res.Parameter(0))
1848 h['W'+name].SetBinError(nb,res.ParError(0))
1849 h['S'+name].SetBinContent(nb,res.Parameter(3))
1850 h['S'+name].SetBinError(nb,res.ParError(3))
1851 else:
1852 print(name,nb,N,'problem')

◆ Mufi_Efficiency()

Survey-MufiScifi.Mufi_Efficiency (   Nev = options.nEvents,
  optionTrack = options.trackType,
  withReco = 'True',
  NbinsRes = 100,
  X = 10. 
)

Definition at line 1298 of file Survey-MufiScifi.py.

1298def Mufi_Efficiency(Nev=options.nEvents,optionTrack=options.trackType,withReco='True',NbinsRes=100,X=10.):
1299
1300 projs = {1:'Y',0:'X'}
1301 v = Scifi.GetConfParF("Scifi/signalSpeed") #signal propagation in fibre
1302 for s in range(1,nScifi+1):
1303 for p in projs:
1304 ut.bookHist(h,'dtScifivsX_'+str(s)+projs[p],'dt vs x track '+projs[p]+";X [cm]; dt [ns]",100,-10.,40.,260,-8.,5.)
1305 ut.bookHist(h,'clN_'+str(s)+projs[p],'cluster size '+projs[p],10,-0.5,9.5)
1306 ut.bookHist(h,'dtScifivsdL_'+str(s),'dt vs dL '+str(s)+";X [cm]; dt [ns]",100,-40.,0.,200,-5.,5.)
1307
1308 for s in systemAndPlanes:
1309 for l in range(systemAndPlanes[s]):
1310 ut.bookHist(h,'dtLRvsX_'+sdict[s]+str(s*10+l),'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",100,-70.,30.,260,-8.,5.)
1311 ut.bookHist(h,'atLRvsX_'+sdict[s]+str(s*10+l),'mean time - T0track vs x '+str(s*10+l)+";X [cm]; dt [ns]",20,-70.,30.,250,-10.,15.0)
1312 ut.bookHist(h,'VetoatLRvsX_'+sdict[s]+str(s*10+l),'mean time - T0Veto vs x '+str(s*10+l)+";X [cm]; dt [ns]",20,-70.,30.,250,-10.,15.0)
1313
1314 scale = 1.
1315 if s==3: scale = 0.4
1316 for side in ['','L','R','S']:
1317 for proj in ['X','Y']:
1318 xmin = -X*NbinsRes/100. * scale
1319 xmax = -xmin
1320 ut.bookHist(h,'res'+proj+'_'+sdict[s]+side+str(s*10+l),'residual '+proj+str(s*10+l)+
1321 '; residual [cm]; local position [cm]',NbinsRes,xmin,xmax,100,-100.,100.)
1322 ut.bookHist(h,'gres'+proj+'_'+sdict[s]+side+str(s*10+l),'residual '+proj+str(s*10+l)+
1323 '; residual [cm]; global position [cm]',NbinsRes,xmin,xmax,100,-100.,100.)
1324 if side=='S': continue
1325 if side=='':
1326 if s==1:
1327 ut.bookHist(h,'resBar'+proj+'_'+sdict[s]+str(s*10+l),'residual '+proj+str(s*10+l)+
1328 '; residual [cm]; bar number',NbinsRes,xmin,xmax,7,-0.5,6.5)
1329 if side=='':
1330 ut.bookHist(h,'track_'+sdict[s]+str(s*10+l),'track x/y '+str(s*10+l)+';X [cm];Y [cm]',100,-90.,10.,100,-20.,80.)
1331 ut.bookHist(h,'locBarPos_'+sdict[s]+str(s*10+l),'bar sizes;X [cm];Y [cm]',100,-100,100,100,-100,100)
1332 ut.bookHist(h,'locEx_'+sdict[s]+str(s*10+l),'loc track pos;X [cm];Y [cm]',100,-100,100,100,-100,100)
1333 for bar in range(systemAndBars[s]):
1334 key = sdict[s]+str(s*10+l)+'_'+str(bar)
1335 if side=="":
1336 ut.bookHist(h,'dtLRvsX_'+key,'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",100,-70.,30.,260,-8.,5.)
1337 ut.bookHist(h,'dtF1LRvsX_'+key,'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",100,-70.,30.,260,-8.,5.)
1338 ut.bookHist(h,'dtfastLRvsX_'+key,'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",100,-70.,30.,260,-8.,5.)
1339 ut.bookHist(h,'atLRvsX_'+key,'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",100,-70.,30.,260,-8.,5.)
1340 else:
1341 ut.bookHist(h,'nSiPMs'+side+'_'+key,'#sipms',16,-0.5,15.5,20,0.,100.)
1342 ut.bookHist(h,'tvsX'+side+'_'+key,"t-t0 vs x track;X [cm]; dt [ns]",100,-70.,30.,200,-12.,12.)
1343 ut.bookHist(h,'tFastvsX'+side+'_'+key,"t-t0 vs x track;X [cm]; dt [ns]",100,-70.,30.,200,-12.,12.)
1344 for i in range(systemAndChannels[s][1]+systemAndChannels[s][0]):
1345 if s==2 and smallSiPMchannel(i):
1346 ut.bookHist(h,'signalS'+side+'_'+key+'-c'+str(i),'signal',200,0.,100.,20,0.,100.)
1347 else:
1348 ut.bookHist(h,'signal'+side+'_'+key+'-c'+str(i),'signal',200,0.,100.,20,0.,100.)
1349 if s==3: continue
1350 ut.bookHist(h,'sigmaTDC'+side+'_'+key+'-c'+str(i),'rms TDC ;dt [ns]',200,-10.0,10.0)
1351 ut.bookHist(h,'TDCcalib'+side+'_'+key+'-c'+str(i),'rms TDC ;dt [ns]',200,-10.0,10.0)
1352 ut.bookHist(h,'sigmaQDC'+side+'_'+key+'-c'+str(i),'rms QDC ; QDC ',200,-50.0,50.)
1353 ut.bookHist(h,'tvsX'+side+'_'+key+'-c'+str(i),"t-t0 vs x track;X [cm]; dt [ns]",100,-70.,30.,200,-12.,12.)
1354
1355 ut.bookHist(h,'signalT'+side+'_'+key,'signal',400,0.,400.,20,0.,100.)
1356 ut.bookHist(h,'signalTS'+side+'_'+key,'signal',400,0.,400.,20,0.,100.)
1357 ut.bookHist(h,'signal'+side+'_'+key,'signal',200,0.,100.,20,0.,100.)
1358
1359 ut.bookHist(h,'resVETOY_1','channel vs residual 1',NbinsRes,xmin,xmax,112,-0.5,111.5)
1360 ut.bookHist(h,'resVETOY_2','channel vs residual 2',NbinsRes,xmin,xmax,112,-0.5,111.5)
1361 ut.bookHist(h,'resVETOX_3','channel vs residual 3',NbinsRes,xmin,xmax,112,-0.5,111.5)
1362
1363 ut.bookHist(h,'trackslxy','track direction',200,-0.1,0.1,200,-0.1,0.1)
1364 ut.bookHist(h,'trackslxy_badChi2','track direction',200,-0.1,0.1,200,-0.1,0.1)
1365 ut.bookHist(h,'tracksChi2Ndof','chi2/ndof',100,0.0,100.,10,-0.5,9.5)
1366 ut.bookHist(h,'NdofvsNMeas','ndof Nmeas',20,-0.5,19.5,20,-0.5,19.5)
1367
1368 if Nev < 0 : Nev = eventTree.GetEntries()
1369 N=0
1370 for event in eventTree:
1371 rc = eventTree.GetEvent(N)
1372 N+=1
1373 if N>Nev: break
1374 if withReco:
1375 for aTrack in Reco_MuonTracks: aTrack.Delete()
1376 Reco_MuonTracks.Clear()
1377 if optionTrack=='DS': rc = trackTask.ExecuteTask("DS")
1378 else : rc = trackTask.ExecuteTask("Scifi")
1379 if not Reco_MuonTracks.GetEntries()==1: continue
1380 theTrack = Reco_MuonTracks[0]
1381 if not theTrack.getFitStatus().isFitConverged() and optionTrack!='DS': # for H8 where only two planes / proj were avaiable
1382 continue
1383# only take horizontal tracks
1384 state = theTrack.getFittedState(0)
1385 pos = state.getPos()
1386 mom = state.getMom()
1387 slopeX= mom.X()/mom.Z()
1388 slopeY= mom.Y()/mom.Z()
1389 if abs(slopeX)>0.25: continue # 4cm distance, 250mrad = 1cm
1390 if abs(slopeY)>0.1: continue
1391
1392# now extrapolate to US and check for hits.
1393 fitStatus = theTrack.getFitStatus()
1394 chi2Ndof = fitStatus.getChi2()/(fitStatus.getNdf()+1E-10)
1395 rc = h['tracksChi2Ndof'].Fill(chi2Ndof,fitStatus.getNdf())
1396 rc = h['NdofvsNMeas'].Fill(fitStatus.getNdf(),theTrack.getNumPointsWithMeasurement())
1397# map clusters to hit keys
1398 DetID2Key={}
1399 if not hasattr(event,"Cluster_Scifi"): # or remakeClusters
1400 if hasattr(trackTask,"clusScifi"):
1401 trackTask.clusScifi.Clear()
1402 trackTask.scifiCluster()
1403 clusters = trackTask.clusScifi
1404 else:
1405 clusters = event.Cluster_Scifi
1406
1407 for aCluster in clusters:
1408 for nHit in range(event.Digi_ScifiHits.GetEntries()):
1409 if event.Digi_ScifiHits[nHit].GetDetectorID()==aCluster.GetFirst():
1410 DetID2Key[aCluster.GetFirst()] = nHit
1411 for aCluster in clusters:
1412 detID = aCluster.GetFirst()
1413 s = int(detID/1000000)
1414 p= int(detID/100000)%10
1415 rc = h['clN_'+str(s)+projs[p]].Fill(aCluster.GetN())
1416
1417 if chi2Ndof> 9 and optionTrack!='DS':
1418 rc = h['trackslxy_badChi2'].Fill(mom.x()/mom.Mag(),mom.y()/mom.Mag())
1419 continue
1420 rc = h['trackslxy'].Fill(mom.x()/mom.Mag(),mom.y()/mom.Mag())
1421# get T0 from Track
1422 if optionTrack=='DS':
1423 # define T0 using mean TDC L/R of the horizontal planes
1424 T0track = 0
1425 Z0track = -1
1426 T0s = []
1427 for nM in range(theTrack.getNumPointsWithMeasurement()):
1428 M = theTrack.getPointWithMeasurement(nM)
1429 W = M.getRawMeasurement()
1430 detID = W.getDetId()
1431 hkey = W.getHitId()
1432 aHit = event.Digi_MuFilterHits[ hkey ]
1433 if aHit.isVertical(): continue
1434 allTDCs = map2Dict(aHit,'GetAllTimes')
1435 if (0 in allTDCs) and (1 in allTDCs) : T0s.append( (allTDCs[0]+allTDCs[1])*TDC2ns/2. )
1436 if len(T0s)==2:
1437 lam = (zPos['MuFilter'][32]-pos.z())/mom.z()
1438 xEx,yEx = lam*mom.x(),pos.y()+lam*mom.y()
1439 # need track length
1440 L = ROOT.TMath.Sqrt( (lam*mom.x())**2 + (lam*mom.y())**2 + (zPos['MuFilter'][32]-pos.z())**2)
1441 ToF = L / u.speedOfLight
1442 T0track = (T0s[0]+T0s[1]-ToF)/2.
1443 Z0track = pos[2]
1444 deltaZ02 = (zPos['MuFilter'][32]-zPos['MuFilter'][30])
1445 # print('debug 1',L,ToF,T0track,Z0track,deltaZ02)
1446 else:
1447 M = theTrack.getPointWithMeasurement(0)
1448 W = M.getRawMeasurement()
1449 detID = W.getDetId()
1450 aHit = event.Digi_ScifiHits[ DetID2Key[detID] ]
1451 Scifi.GetSiPMPosition(detID,A,B)
1452 X = B-pos
1453 L0 = X.Mag()/v
1454 # need to correct for signal propagation along fibre
1455 clkey = W.getHitId()
1456 aCl = clusters[clkey]
1457 T0track = aCl.GetTime() - L0
1458 TZero = aCl.GetTime()
1459 Z0track = pos[2]
1460 times = {}
1461 for nM in range(theTrack.getNumPointsWithMeasurement()):
1462 state = theTrack.getFittedState(nM)
1463 posM = state.getPos()
1464 M = theTrack.getPointWithMeasurement(nM)
1465 W = M.getRawMeasurement()
1466 detID = W.getDetId()
1467 clkey = W.getHitId()
1468 aCl = clusters[clkey]
1469 aHit = event.Digi_ScifiHits[ DetID2Key[detID] ]
1470 Scifi.GetSiPMPosition(detID,A,B)
1471 if aHit.isVertical(): X = B-posM
1472 else: X = A-posM
1473 L = X.Mag()/v
1474 # need to correct for signal propagation along fibre
1475 dT = aCl.GetTime() - L - T0track - (posM[2] -Z0track)/u.speedOfLight
1476 ss = str(aHit.GetStation())
1477 prj = 'X'
1478 l = posM[0]
1479 if aHit.isVertical():
1480 prj='Y'
1481 l = posM[1]
1482 rc = h['dtScifivsX_'+ss+prj].Fill(X.Mag(),dT)
1483 times[ss+prj]=[aCl.GetTime(),L*v,detID,l]
1484 for s in range(1,nScifi+1):
1485 if str(s)+'X' in times and str(s)+'Y' in times:
1486 deltaT = times[str(s)+'X'][0] - times[str(s)+'Y'][0]
1487 deltaL = times[str(s)+'X'][1] - times[str(s)+'Y'][1]
1488 rc = h['dtScifivsdL_'+str(s)].Fill(deltaL,deltaT)
1489 deltaZ02 = 40. # to be fixed
1490 ToF = 1.
1491 #print(detID,aHit.GetDetectorID(),aHit.GetTime()*TDC2ns-TZero,dT,L,aHit.GetTime()*TDC2ns - L,T0)
1492
1493 muHits = {}
1494 for s in systemAndPlanes:
1495 for p in range(systemAndPlanes[s]):
1496 if s == 3:
1497 muHits[s*10+2*p]=[]
1498 muHits[s*10+2*p+1]=[]
1499 else:
1500 muHits[s*10+p]=[]
1501 for p in range(systemAndPlanes[s]): muHits[s*10+p]=[]
1502 for aHit in event.Digi_MuFilterHits:
1503 if not aHit.isValid(): continue
1504 s = aHit.GetDetectorID()//10000
1505 p = (aHit.GetDetectorID()//1000)%10
1506 bar = (aHit.GetDetectorID()%1000)%60
1507 plane = s*10+p
1508 if s==3:
1509 if aHit.isVertical():
1510 plane = s*10+2*p+1
1511 if p==3: plane = s*10+2*p
1512 else: plane = s*10+2*p
1513 muHits[plane].append(aHit)
1514
1515# get T0 from VETO
1516 s = 1
1517 Z0Veto = zPos['MuFilter'][1*10+0]
1518 dZ = zPos['MuFilter'][1*10+1] - zPos['MuFilter'][1*10+0]
1519 avT = {}
1520 for p in range(systemAndPlanes[s]-1):# exclude the vertical Veto3
1521 plane = s*10+p
1522 if len(muHits[plane])!=1: continue
1523 aHit = muHits[plane][0]
1524# check if hit within track extrapolation
1525 zEx = zPos['MuFilter'][s*10+plane]
1526 lam = (zEx-pos.z())/mom.z()
1527 xEx,yEx = pos.x()+lam*mom.x(),pos.y()+lam*mom.y()
1528 detID = aHit.GetDetectorID()
1529 MuFilter.GetPosition(detID,A,B)
1530 D = (A[1]+B[1])/2. - yEx
1531 if abs(D)>5: continue
1532 avT[plane] = aHit.GetImpactT()
1533 T0Veto = -999
1534 if len(avT)==2:
1535 T0Veto = (avT[10]+(avT[11]-dZ/u.speedOfLight))/2.
1536
1537 vetoHits = {0:[],1:[], 2:[]}
1538 for s in sdict:
1539 name = str(s)
1540 for plane in range(systemAndPlanes[s]):
1541 zEx = zPos['MuFilter'][s*10+plane]
1542 lam = (zEx-pos.z())/mom.z()
1543 xEx,yEx = pos.x()+lam*mom.x(),pos.y()+lam*mom.y()
1544 # tag with station close by
1545 if plane ==0: tag = 1
1546 else: tag = plane -1
1547 tagged = False
1548 for aHit in muHits[s*10+tag]:
1549 detID = aHit.GetDetectorID()
1550 MuFilter.GetPosition(detID,A,B)
1551 if aHit.isVertical() : D = (A[0]+B[0])/2. - xEx
1552 else: D = (A[1]+B[1])/2. - yEx
1553 if abs(D)<5: tagged = True
1554 #if not tagged: continue
1555 rc = h['track_'+sdict[s]+str(s*10+plane)].Fill(xEx,yEx)
1556 for aHit in muHits[s*10+plane]:
1557 detID = aHit.GetDetectorID()
1558 bar = (detID%1000)%60
1559 nSiPMs = aHit.GetnSiPMs()
1560 nSides = aHit.GetnSides()
1561 MuFilter.GetPosition(detID,globA,globB)
1562 MuFilter.GetLocalPosition(detID,locA,locB)
1563 globEx = array('d',[xEx,yEx,zEx])
1564 locEx = array('d',[0,0,0])
1565 nav.MasterToLocal(globEx,locEx)
1566 locPos = 0.5* (locA+locB)
1567 globPos = 0.5 * (globA+globB)
1568 dy = locPos[1] - locEx[1]
1569 dx = locPos[0] - locEx[0]
1570 gdy = globPos[1] - globEx[1]
1571 gdx = globPos[0] - globEx[0]
1572 rc = h['locBarPos_'+sdict[s]+str(s*10+plane)].Fill( locPos[0],locPos[1])
1573 rc = h['locEx_'+sdict[s]+str(s*10+plane)].Fill( locEx[0],locEx[1])
1574 rc = h['resY_'+sdict[s]+str(s*10+plane)].Fill(dy,locEx[0])
1575 rc = h['resX_'+sdict[s]+str(s*10+plane)].Fill(dx,locEx[1])
1576 rc = h['gresY_'+sdict[s]+str(s*10+plane)].Fill(gdy,globEx[0])
1577 rc = h['gresX_'+sdict[s]+str(s*10+plane)].Fill(gdx,globEx[1])
1578 S = map2Dict(aHit,'GetAllSignals')
1579 # check for signal in left / right or small sipm
1580 left,right,smallL,smallR,Sleft,Sright,SleftS,SrightS = 0,0,0,0,0,0,0,0
1581 if s==1:
1582 if plane<2:
1583 vetoHits[plane].append( [gdy,bar] )
1584 rc = h['resBarY_'+sdict[s]+str(s*10+plane)].Fill(gdy,bar)
1585 elif plane==2:
1586 vetoHits[plane].append( [gdx,bar] )
1587 rc = h['resBarX_'+sdict[s]+str(s*10+plane)].Fill(gdx,bar)
1588 for x in S:
1589 if s==1:
1590 if plane<2:
1591 nc = x + 2*nSiPMs*bar
1592 h['resVETOY_'+str(plane+1)].Fill(dy,nc)
1593 elif plane==2:
1594 nc = x + nSiPMs*bar
1595 h['resVETOX_'+str(plane+1)].Fill(dx,nc)
1596 if x<nSiPMs:
1597 if s==2 and smallSiPMchannel(x): smallL+=1
1598 else: left+=1
1599 else:
1600 if s==2 and smallSiPMchannel(x): smallR+=1
1601 else: right+=1
1602 if left>0:
1603 rc = h['resY_'+sdict[s]+'L'+str(s*10+plane)].Fill(dy,locEx[1])
1604 rc = h['resX_'+sdict[s]+'L'+str(s*10+plane)].Fill(dx,locEx[0])
1605 rc = h['gresY_'+sdict[s]+'L'+str(s*10+plane)].Fill(gdy,globEx[1])
1606 rc = h['gresX_'+sdict[s]+'L'+str(s*10+plane)].Fill(gdx,globEx[0])
1607 if right>0:
1608 rc = h['resY_'+sdict[s]+'R'+str(s*10+plane)].Fill(dy,locEx[1])
1609 rc = h['resX_'+sdict[s]+'R'+str(s*10+plane)].Fill(dx,locEx[0])
1610 rc = h['gresY_'+sdict[s]+'R'+str(s*10+plane)].Fill(gdy,globEx[1])
1611 rc = h['gresX_'+sdict[s]+'R'+str(s*10+plane)].Fill(gdx,globEx[0])
1612 if s==2 and (smallL>0 or smallR>0):
1613 rc = h['resY_'+sdict[s]+'S'+str(s*10+plane)].Fill(dy,locEx[1])
1614 rc = h['resX_'+sdict[s]+'S'+str(s*10+plane)].Fill(dx,locEx[0])
1615 rc = h['gresY_'+sdict[s]+'S'+str(s*10+plane)].Fill(gdy,globEx[1])
1616 rc = h['gresX_'+sdict[s]+'S'+str(s*10+plane)].Fill(gdx,globEx[0])
1617 dist = abs(dy)
1618 if aHit.isVertical() : dist = abs(dx)
1619 if dist<3.0: # check channels
1620 if aHit.isVertical():
1621 dL = locA[1]- locEx[1]
1622 dR = locEx[1] - locB[1]
1623 else:
1624 dR = locA[0] - locEx[0]
1625 dL = locEx[0] - locB[0]
1626 barName = sdict[s]+str(s*10+plane)+'_'+str(bar)
1627 rc = h['nSiPMsL_'+barName].Fill(left,dL)
1628 rc = h['nSiPMsR_'+barName].Fill(right,dR)
1629 for x in S:
1630 qcd = S[x]
1631 if x<nSiPMs:
1632 if s==2 and smallSiPMchannel(x):
1633 rc = h['signalSL_'+barName+'-c'+str(x)].Fill(qcd,dL)
1634 SleftS+=qcd
1635 else:
1636 rc = h['signalL_'+barName+'-c'+str(x)].Fill(qcd,dL)
1637 Sleft+=qcd
1638 else:
1639 if s==2 and smallSiPMchannel(x):
1640 rc = h['signalSR_'+barName+'-c'+str(x-nSiPMs)].Fill(qcd,dR)
1641 SrightS+=qcd
1642 else:
1643 rc = h['signalR_'+barName+'-c'+str(x-nSiPMs)].Fill(qcd,dR)
1644 Sright+=qcd
1645 rc = h['signalTL_'+barName].Fill(Sleft,dL)
1646 rc = h['signalTR_'+barName].Fill(Sright,dR)
1647 rc = h['signalTSL_'+barName].Fill(SleftS,dL)
1648 rc = h['signalTSR_'+barName].Fill(SrightS,dR)
1649
1650# look at delta time vs track X, works only for horizontal planes.
1651 allTDCs = map2Dict(aHit,'GetAllTimes')
1652 if not aHit.isVertical():
1653 dt = aHit.GetDeltaT()
1654 dtF = aHit.GetFastDeltaT()
1655 mtVeto = aHit.GetImpactT() - T0Veto - (globPos[2] - Z0Veto)/u.speedOfLight
1656 h['dtLRvsX_'+sdict[s]+str(s*10+plane)].Fill(xEx,dt*TDC2ns)
1657 h['dtLRvsX_'+barName].Fill(xEx,dt*TDC2ns)
1658 if (1 in allTDCs) and (9 in allTDCs): h['dtF1LRvsX_'+barName].Fill(xEx,(allTDCs[1]-allTDCs[9])*TDC2ns)
1659 h['dtfastLRvsX_'+barName].Fill(xEx,dtF*TDC2ns)
1660 if Z0track>0:
1661 tcor = (globPos[2] - Z0track)/deltaZ02 * ToF
1662 mtTrack = aHit.GetImpactT() - T0track - tcor
1663 h['atLRvsX_'+sdict[s]+str(s*10+plane)].Fill(xEx,mtTrack)
1664 h['atLRvsX_'+barName].Fill(xEx,mtTrack)
1665 if s <3:
1666 for i in allTDCs:
1667 dt = allTDCs[i]*TDC2ns-T0track - tcor
1668 #print('debug 2',tcor,dt)
1669 if i<nSiPMs: h['tvsXL_'+barName+'-c'+str(i)].Fill(xEx,dt)
1670 else:
1671 h['tvsXR_'+barName+'-c'+str(i-nSiPMs)].Fill(xEx,dt)
1672
1673 h['VetoatLRvsX_'+sdict[s]+str(s*10+plane)].Fill(xEx,mtVeto)
1674# QDC/TDC channel variations
1675 if s==3: continue
1676 meanL,meanR,nL,nR=0,0,0,0
1677 t0Left = 999
1678 t0Right = 999
1679 for i in allTDCs:
1680 if i==4: t0Left = allTDCs[i]
1681 if i==12: t0Right = allTDCs[i]
1682
1683 for i in allTDCs:
1684 if s==2 and smallSiPMchannel(i): continue
1685 if i < nSiPMs: # left side
1686 nL+=1
1687 meanL+=allTDCs[i]
1688 else:
1689 nR+=1
1690 meanR+=allTDCs[i]
1691 for i in allTDCs:
1692 if s==2 and smallSiPMchannel(i): continue
1693 if i<nSiPMs and nL>0:
1694 key = sdict[s]+str(s*10+plane)+'_'+str(bar)+'-c'+str(i)
1695 rc = h['sigmaTDCL_'+key].Fill( (allTDCs[i]-meanL/nL)*TDC2ns )
1696 if t0Left<900: rc = h['TDCcalibL_'+key].Fill( (allTDCs[i]-t0Left)*TDC2ns )
1697 elif not(i<nSiPMs) and nR>0:
1698 key = sdict[s]+str(s*10+plane)+'_'+str(bar)+'-c'+str(i-nSiPMs)
1699 rc = h['sigmaTDCR_'+key].Fill((allTDCs[i]-meanR/nR)*TDC2ns )
1700 if t0Right<900: rc = h['TDCcalibR_'+key].Fill( (allTDCs[i]-t0Right)*TDC2ns )
1701
1702 meanL,meanR,nL,nR=0,0,0,0
1703 for i in S:
1704 if s==2 and smallSiPMchannel(i): continue
1705 if i < nSiPMs: # left side
1706 nL+=1
1707 meanL+=S[i]
1708 else:
1709 nR+=1
1710 meanR+=S[i]
1711 for i in S:
1712 if s==2 and smallSiPMchannel(i): continue
1713 if i<nSiPMs and nL>0:
1714 key = sdict[s]+str(s*10+plane)+'_'+str(bar)+'-c'+str(i)
1715 rc = h['sigmaQDCL_'+key].Fill( (S[i]-meanL/nL) )
1716 elif not(i<nSiPMs) and nR>0:
1717 key = sdict[s]+str(s*10+plane)+'_'+str(bar)+'-c'+str(i-nSiPMs)
1718 rc = h['sigmaQDCR_'+key].Fill((S[i]-meanR/nR) )
1719
1720 for s in [1,2]:
1721 for l in range(systemAndPlanes[s]):
1722 for side in ['L','R']:
1723 for bar in range(systemAndBars[s]):
1724 key = 'sigmaTDC'+side+'_'+sdict[s]+str(s*10+l)+'_'+str(bar)
1725 h[key]=h[key+'-c0'].Clone(key)
1726 h[key].Reset()
1727 for i in range(systemAndChannels[s][1]+systemAndChannels[s][0]):
1728 h[key].Add(h[key+'-c'+str(i)])
1729
1730 ut.writeHists(h,'MuFilterEff_run'+str(options.runNumber)+'.root')
1731
MuFilter::GetLocalPosition
void GetLocalPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition MuFilter.cxx:563
Scifi::GetSiPMPosition
void GetSiPMPosition(Int_t SiPMChan, TVector3 &A, TVector3 &B)
Definition Scifi.cxx:727

◆ Mufi_hitMaps()

Survey-MufiScifi.Mufi_hitMaps (   Nev = options.nEvents)

Definition at line 473 of file Survey-MufiScifi.py.

473def Mufi_hitMaps(Nev = options.nEvents):
474 # veto system 2 layers with 7 bars and 8 sipm channels on both ends
475 # US system 5 layers with 10 bars and 8 sipm channels on both ends
476 # DS system horizontal(3) planes, 60 bars, readout on both sides, single channel
477 # vertical(4) planes, 60 bar, readout on top, single channel
478 for s in systemAndPlanes:
479 for l in range(systemAndPlanes[s]):
480 ut.bookHist(h,'hitmult_'+str(s*10+l),'hit mult / plane '+str(s*10+l),61,-0.5,60.5)
481 ut.bookHist(h,'hit_'+str(s*10+l),'channel map / plane '+str(s*10+l),160,-0.5,159.5)
482 if s==3: ut.bookHist(h,'bar_'+str(s*10+l),'hit map / plane '+str(s*10+l),60,-0.5,59.5)
483 else: ut.bookHist(h,'bar_'+str(s*10+l),'hit map / plane '+str(s*10+l),10,-0.5,9.5)
484 ut.bookHist(h,'sig_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
485 if s==2: ut.bookHist(h,'sigS_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
486 ut.bookHist(h,'sigL_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
487 ut.bookHist(h,'sigR_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
488 ut.bookHist(h,'Tsig_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
489 ut.bookHist(h,'occ_'+str(s*10+l),'channel occupancy '+str(s*10+l),100,0.0,200.)
490 ut.bookHist(h,'occTag_'+str(s*10+l),'channel occupancy '+str(s*10+l),100,0.0,200.)
491
492 ut.bookHist(h,'leftvsright_1','Veto hits in left / right',10,-0.5,9.5,10,-0.5,9.5)
493 ut.bookHist(h,'leftvsright_2','US hits in left / right',10,-0.5,9.5,10,-0.5,9.5)
494 ut.bookHist(h,'leftvsright_3','DS hits in left / right',2,-0.5,1.5,2,-0.5,1.5)
495 ut.bookHist(h,'leftvsright_signal_1','Veto signal in left / right',100,-0.5,200.,100,-0.5,200.)
496 ut.bookHist(h,'leftvsright_signal_2','US signal in left / right',100,-0.5,200.,100,-0.5,200.)
497 ut.bookHist(h,'leftvsright_signal_3','DS signal in left / right',100,-0.5,200.,100,-0.5,200.)
498
499 ut.bookHist(h,'dtime','delta event time; dt [ns]',100,0.0,1000.)
500 ut.bookHist(h,'dtimeu','delta event time; dt [us]',100,0.0,1000.)
501 ut.bookHist(h,'dtimem','delta event time; dt [ms]',100,0.0,1000.)
502
503 ut.bookHist(h,'bs','beam spot',100,-100.,10.,100,0.,80.)
504 ut.bookHist(h,'bsDS','beam spot',60,-0.5,59.5,60,-0.5,59.5)
505 ut.bookHist(h,'slopes','track slopes',100,-0.1,0.1,100,-0.1,0.1)
506
507 for s in range(1,4):
508 ut.bookHist(h,sdict[s]+'Mult','QDCs vs nr hits',200,0.,800.,200,0.,300.)
509
510 N=-1
511 Tprev = 0
512 if Nev < 0 : Nev = eventTree.GetEntries()
513 eventTree.GetEvent(0)
514 t0 = eventTree.EventHeader.GetEventTime()/freq
515 listOfHits = {1:[],2:[],3:[]}
516 mult = {}
517 for event in eventTree:
518 N+=1
519 if N%options.heartBeat == 0: print('event ',N,' ',time.ctime())
520 if N>Nev: break
521 withX = False
522 planes = {}
523 listOfHits[1].clear()
524 listOfHits[2].clear()
525 listOfHits[3].clear()
526 for s in systemAndPlanes:
527 for l in range(systemAndPlanes[s]): mult[s*10+l]=0
528 for aHit in event.Digi_MuFilterHits:
529 if not aHit.isValid(): continue
530 detID = aHit.GetDetectorID()
531 if aHit.isVertical(): withX = True
532 s = detID//10000
533 l = (detID%10000)//1000 # plane number
534 bar = (detID%1000)
535 if s>2:
536 l=2*l
537 if bar>59:
538 bar=bar-60
539 if l<6: l+=1
540 mult[s*10+l]+=1
541 key = s*100+l
542 if not key in planes: planes[key] = {}
543 sumSignal = map2Dict(aHit,'SumOfSignals')
544 planes[key][bar] = [sumSignal['SumL'],sumSignal['SumR']]
545 nSiPMs = aHit.GetnSiPMs()
546 nSides = aHit.GetnSides()
547# check left/right
548 allChannels = map2Dict(aHit,'GetAllSignals')
549 for c in allChannels:
550 listOfHits[s].append(allChannels[c])
551 if nSides==2:
552 Nleft = 0
553 Nright = 0
554 Sleft = 0
555 Sright = 0
556 for c in allChannels:
557 if nSiPMs > c: # left side
558 Nleft+=1
559 Sleft+=allChannels[c]
560 else:
561 Nright+=1
562 Sright+=allChannels[c]
563 rc = h['leftvsright_'+str(s)].Fill(Nleft,Nright)
564 rc = h['leftvsright_signal_'+str(s)].Fill(Sleft,Sright)
565#
566 for c in allChannels:
567 channel = bar*nSiPMs*nSides + c
568 rc = h['hit_'+str(s)+str(l)].Fill( int(channel))
569 rc = h['bar_'+str(s)+str(l)].Fill(bar)
570 if s==2 and smallSiPMchannel(c) : rc = h['sigS_'+str(s)+str(l)].Fill(allChannels[c])
571 elif c<nSiPMs: rc = h['sigL_'+str(s)+str(l)].Fill(allChannels[c])
572 else : rc = h['sigR_'+str(s)+str(l)].Fill(allChannels[c])
573 rc = h['sig_'+str(s)+str(l)].Fill(allChannels[c])
574 allChannels.clear()
575#
576 for s in listOfHits:
577 nhits = len(listOfHits[s])
578 qcdsum = 0
579 for i in range(nhits):
580 rc = h[sdict[s]+'Mult'].Fill(nhits, listOfHits[s][i])
581 for s in systemAndPlanes:
582 for l in range(systemAndPlanes[s]):
583 rc = h['hitmult_'+str(s*10+l)].Fill(mult[s*10+l])
584
585 maxOneBar = True
586 for key in planes:
587 if len(planes[key]) > 2: maxOneBar = False
588 if withX and maxOneBar: beamSpot()
589
590 S = {1:[1800,800,2,1],2:[1800,1500,2,3],3:[1800,1800,2,4]}
591 for s in S:
592 ut.bookCanvas(h,'hitmaps' +str(s),' ',S[s][0],S[s][1],S[s][2],S[s][3])
593 ut.bookCanvas(h,'barmaps'+str(s),' ',S[s][0],S[s][1],S[s][2],S[s][3])
594 ut.bookCanvas(h,'signal' +str(s),' ',S[s][0],S[s][1],S[s][2],S[s][3])
595 ut.bookCanvas(h,'Tsignal' +str(s),' ',S[s][0],S[s][1],S[s][2],S[s][3])
596
597 for l in range(systemAndPlanes[s]):
598 n = l+1
599 if s==3 and n==7: n=8
600 tc = h['hitmaps'+str(s)].cd(n)
601 tag = str(s)+str(l)
602 h['hit_'+tag].Draw()
603 tc = h['barmaps'+str(s)].cd(n)
604 h['bar_'+tag].Draw()
605 tc = h['signal'+str(s)].cd(n)
606 h['sig_'+tag].Draw()
607 tc = h['Tsignal'+str(s)].cd(n)
608 h['Tsig_'+tag].Draw()
609
610 ut.bookCanvas(h,'hitmult','hit multiplicities per plane',2000,1600,4,3)
611 k=1
612 for s in systemAndPlanes:
613 for l in range(systemAndPlanes[s]):
614 tc = h['hitmult'].cd(k)
615 tc.SetLogy(1)
616 k+=1
617 rc = h['hitmult_'+str(s*10+l)].Draw()
618
619 ut.bookCanvas(h,'VETO',' ',1200,1800,1,2)
620 for l in range(2):
621 tc = h['VETO'].cd(l+1)
622 hname = 'hit_'+str(1)+str(l)
623 h[hname].SetStats(0)
624 h[hname].Draw()
625 for n in range(7):
626 x = (n+1)*16-0.5
627 y = h['hit_'+str(1)+str(l)].GetMaximum()
628 lname = 'L'+str(n)+hname
629 h[lname] = ROOT.TLine(x,0,x,y)
630 h[lname].SetLineColor(ROOT.kRed)
631 h[lname].SetLineStyle(9)
632 h[lname].Draw('same')
633
634 ut.bookCanvas(h,'USBars',' ',1200,900,1,1)
635 colours = {0:ROOT.kOrange,1:ROOT.kRed,2:ROOT.kGreen,3:ROOT.kBlue,4:ROOT.kMagenta,5:ROOT.kCyan,
636 6:ROOT.kAzure,7:ROOT.kPink,8:ROOT.kSpring}
637 for i in range(5):
638 h['bar_2'+str(i)].SetLineColor(colours[i])
639 h['bar_2'+str(i)].SetLineWidth(2)
640 h['bar_2'+str(i)].SetStats(0)
641 h['bar_20'].Draw()
642 h['bar_21'].Draw('same')
643 h['bar_22'].Draw('same')
644 h['bar_23'].Draw('same')
645 h['bar_24'].Draw('same')
646 h['lbar2']=ROOT.TLegend(0.6,0.6,0.99,0.99)
647 for i in range(5):
648 h['lbar2'].AddEntry(h['bar_2'+str(i)],'plane '+str(i+1),"f")
649 h['lbar2'].Draw()
650 for i in range(7):
651 h['hit_3'+str(i)].SetLineColor(colours[i])
652 h['hit_3'+str(i)].SetLineWidth(2)
653 h['hit_3'+str(i)].SetStats(0)
654 h['hit_30'].Draw()
655 for i in range(1,7):
656 h['hit_3'+str(i)].Draw('same')
657 h['lbar3']=ROOT.TLegend(0.6,0.6,0.99,0.99)
658 for i in range(7):
659 h['lbar3'].AddEntry(h['hit_3'+str(i)],'plane '+str(i+1),"f")
660 h['lbar3'].Draw()
661
662 ut.bookCanvas(h,'LR',' ',1800,900,3,2)
663 h['LR'].cd(1)
664 h['leftvsright_'+str(1)].Draw('textBox')
665 h['LR'].cd(2)
666 h['leftvsright_'+str(2)].Draw('textBox')
667 h['LR'].cd(3)
668 h['leftvsright_'+str(3)].Draw('textBox')
669 h['LR'].cd(4)
670 h['leftvsright_signal_1'].SetMaximum(h['leftvsright_signal_1'].GetBinContent(10,10))
671 h['leftvsright_signal_2'].SetMaximum(h['leftvsright_signal_2'].GetBinContent(10,10))
672 h['leftvsright_signal_3'].SetMaximum(h['leftvsright_signal_3'].GetBinContent(10,10))
673 h['leftvsright_signal_'+str(1)].Draw('colz')
674 h['LR'].cd(5)
675 h['leftvsright_signal_'+str(2)].Draw('colz')
676 h['LR'].cd(6)
677 h['leftvsright_signal_'+str(3)].Draw('colz')
678
679 ut.bookCanvas(h,'LRinEff',' ',1800,450,3,1)
680 for s in range(1,4):
681 h['lLRinEff'+str(s)]=ROOT.TLegend(0.6,0.54,0.99,0.93)
682 name = 'leftvsright_signal_'+str(s)
683 h[name+'0Y'] = h[name].ProjectionY(name+'0Y',1,1)
684 h[name+'0X'] = h[name].ProjectionX(name+'0X',1,1)
685 h[name+'1X'] = h[name].ProjectionY(name+'1Y')
686 h[name+'1Y'] = h[name].ProjectionX(name+'1X')
687 tc = h['LRinEff'].cd(s)
688 tc.SetLogy()
689 h[name+'0X'].SetStats(0)
690 h[name+'0Y'].SetStats(0)
691 h[name+'1X'].SetStats(0)
692 h[name+'1Y'].SetStats(0)
693 h[name+'0X'].SetLineColor(ROOT.kRed)
694 h[name+'0Y'].SetLineColor(ROOT.kGreen)
695 h[name+'1X'].SetLineColor(ROOT.kMagenta)
696 h[name+'1Y'].SetLineColor(ROOT.kCyan)
697 h[name+'0X'].SetMaximum(max(h[name+'1X'].GetMaximum(),h[name+'1Y'].GetMaximum()))
698 h[name+'0X'].Draw()
699 h[name+'0Y'].Draw('same')
700 h[name+'1X'].Draw('same')
701 h[name+'1Y'].Draw('same')
702 # Fill(Sleft,Sright)
703 h['lLRinEff'+str(s)].AddEntry(h[name+'0X'],'left with no signal right',"f")
704 h['lLRinEff'+str(s)].AddEntry(h[name+'0Y'],'right with no signal left',"f")
705 h['lLRinEff'+str(s)].AddEntry(h[name+'1X'],'left all',"f")
706 h['lLRinEff'+str(s)].AddEntry(h[name+'1Y'],'right all',"f")
707 h['lLRinEff'+str(s)].Draw()
708
709 ut.bookCanvas(h,'signalUSVeto',' ',1200,1600,3,7)
710 s = 1
711 l = 1
712 for plane in range(2):
713 for side in ['L','R','S']:
714 tc = h['signalUSVeto'].cd(l)
715 l+=1
716 if side=='S': continue
717 h['sig'+side+'_'+str( s*10+plane)].Draw()
718 s=2
719 for plane in range(5):
720 for side in ['L','R','S']:
721 tc = h['signalUSVeto'].cd(l)
722 l+=1
723 h['sig'+side+'_'+str( s*10+plane)].Draw()
724 ut.bookCanvas(h,'signalDS',' ',900,1600,2,7)
725 s = 3
726 l = 1
727 for plane in range(7):
728 for side in ['L','R']:
729 tc = h['signalDS'].cd(l)
730 l+=1
731 h['sig'+side+'_'+str( s*10+plane)].Draw()
732
733
734 for canvas in ['signalUSVeto','LR','USBars']:
735 h[canvas].Update()
736 myPrint(h[canvas],canvas)
737 for canvas in ['hitmaps','barmaps','signal','Tsignal']:
738 for s in range(1,4):
739 h[canvas+str(s)].Update()
740 myPrint(h[canvas+str(s)],canvas+sdict[s])

◆ myPrint()

Survey-MufiScifi.myPrint (   tc,
  name,
  withRootFile = True 
)

Definition at line 271 of file Survey-MufiScifi.py.

271def myPrint(tc,name,withRootFile=True):
272 tc.Update()
273 tc.Print(name+'-run'+str(options.runNumber)+'.png')
274 tc.Print(name+'-run'+str(options.runNumber)+'.pdf')
275 if withRootFile: tc.Print(name+'-run'+str(options.runNumber)+'.root')
276

◆ plotMips()

Survey-MufiScifi.plotMips (   readhisto = True)

Definition at line 1856 of file Survey-MufiScifi.py.

1856def plotMips(readhisto=True):
1857 if readhisto: ut.readHists(h,'LandauFits_run'+str(options.runNumber)+'.root',withProj=False)
1858 langau = ROOT.TF1('Langau',langaufun,0.,100.,4)
1859 ut.bookCanvas(h,'CsignalT','',1200,2000,2,5)
1860 ut.bookCanvas(h,'WsignalT','',1200,2000,2,5)
1861 ut.bookCanvas(h,'SsignalT','',1200,2000,2,5)
1862 ut.bookCanvas(h,'Msignal1','',900,600,1,1)
1863 ut.bookCanvas(h,'MsignalT','',1200,2000,2,5)
1864 ut.bookCanvas(h,'MsignalT1','',900,600,1,1)
1865 ut.bookCanvas(h,'NnSiPMs1','',900,600,1,1)
1866 ut.bookCanvas(h,'NnSiPMs','',1200,2000,2,5)
1867#example plots:
1868 for T in ['T']:
1869 name = 'signal'+T+'L_US22_5-X10'
1870 tc = h['Msignal'+T+'1'].cd()
1871 h[name].Draw()
1872 tc.Update()
1873 stats = h[name].FindObject('stats')
1874 stats.SetOptFit(1111111)
1875 stats.SetX1NDC(0.51)
1876 stats.SetY1NDC(0.54)
1877 stats.SetX2NDC(0.98)
1878 stats.SetY2NDC(0.94)
1879 tc.Update()
1880 myPrint(h['Msignal'+T+'1'],'signal'+T+'1')
1881 tc = h['NnSiPMs1'].cd(1)
1882 h['nSiPMs1'] = h['nSiPMsL_US21_5'].ProjectionX('nSiPMs1',1,9)
1883 h['nSiPMs1'].Draw('hist')
1884 tc.Update()
1885 myPrint(h['NnSiPMs1'],'nSiPMs1')
1886
1887 s=2
1888# make summary of mpv and att length
1889 ut.bookCanvas(h,'Mpvs','',2400,1200,1,2)
1890 ut.bookCanvas(h,'attLs','',2400,1200,1,2)
1891 tc = h['Mpvs'].cd(1)
1892 h['attLengthL'],h['attLengthR'],h['attLengthLC'],h['attLengthRC']= ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph()
1893 h['mpv0L'],h['mpv0R'],h['mpv0LC'],h['mpv0RC'] = ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph()
1894 for t in ['T','']:
1895 for p in ['L','R']:
1896 nk = 0
1897 for plane in range(5):
1898 for bar in range(10):
1899 if t=='T': nloop = 1
1900 else: nloop = systemAndChannels[s][1]+systemAndChannels[s][0]
1901 for i in range(nloop):
1902 if s==2 and smallSiPMchannel(i): continue
1903 name = 'Msignal'+t+p+'_US'+str(s*10+plane)+'_'+str(bar)
1904 chi2N = 'Csignal'+t+p+'_US'+str(s*10+plane)+'_'+str(bar)
1905 if nloop>1:
1906 name += '-c'+str(i)
1907 chi2N += '-c'+str(i)
1908 # find fit limits
1909 limits = []
1910 for i in range(1,h[chi2N].GetNbinsX()+1):
1911 funL = h[name[1:]+'-X'+str(i)].GetFunction('langau')
1912 if not funL: continue
1913 if not funL.GetParameter(1)>0: continue
1914 mpvRelErr = funL.GetParError(1)/funL.GetParameter(1)*ROOT.TMath.Sqrt(h[name[1:]+'-X'+str(i)].GetEntries())
1915 if mpvRelErr < 0.1:
1916 print('exclude bin ',name,i,mpvRelErr,h[name[1:]+'-X'+str(i)].GetEntries())
1917 continue
1918 if len(limits)==0 and h[chi2N].GetBinContent(i)>0 and h[chi2N].GetBinContent(i)<9:
1919 limits.append(h[chi2N].GetBinCenter(i))
1920 break
1921 for i in range(h[chi2N].GetNbinsX(),1,-1):
1922 funL = h[name[1:]+'-X'+str(i)].GetFunction('langau')
1923 if not funL: continue
1924 if not funL.GetParameter(1)>0: continue
1925 mpvRelErr = funL.GetParError(1)/funL.GetParameter(1)*ROOT.TMath.Sqrt(h[name[1:]+'-X'+str(i)].GetEntries())
1926 if mpvRelErr < 0.1:
1927 print('exclude bin ',name,i,mpvRelErr,h[name[1:]+'-X'+str(i)].GetEntries())
1928 continue
1929 if len(limits)==1 and h[chi2N].GetBinContent(i)>0 and h[chi2N].GetBinContent(i)<9:
1930 limits.append(h[chi2N].GetBinCenter(i))
1931 break
1932 if len(limits)<2:
1933 X = nk+0.5
1934 tag = p
1935 if nloop>1:
1936 tag = p+'C'
1937 X = nk/6.
1938 h['mpv0'+tag].SetPoint(nk,X,-1)
1939 h['attLength'+tag].SetPoint(nk,X,0)
1940 nk+=1
1941 continue
1942 rc = h[name].Fit('expo','SQ','',limits[0],limits[1])
1943 fun = h[name].GetFunction('expo')
1944 fun.SetLineColor(barColor[bar%10])
1945 res = rc.Get()
1946 tag = p
1947 X = nk+0.5
1948 val = fun.Eval(40)
1949 if nloop>1:
1950 tag = p+'C'
1951 X = nk/6.
1952 val = fun.Eval(0)
1953 h['mpv0'+tag].SetPoint(nk,X,val)
1954 h['attLength'+tag].SetPoint(nk,X,-1./res.Parameter(1))
1955 nk+=1
1956
1957 ut.bookHist(h,'hMpvs',';stations with bars; MPV [QCD]',101,-0.5,50.5)
1958 ut.bookHist(h,'hAttl',';plane number; att length [cm]',101,-0.5,50.5)
1959 params = {'hMpvs':['Mpvs',20,'mpv0',0.], 'hAttl':['attLs',0.,'attLength',-2000.]}
1960 side = {'L':[ROOT.kBlue,''], 'R':[ROOT.kRed,'same']}
1961 for P in params:
1962 h[P].SetStats(0)
1963 h[P].SetMaximum(params[P][1])
1964 h[P].SetMinimum(params[P][3])
1965 X = h[P].GetXaxis()
1966 X.SetTitleFont(42)
1967 X.SetNdivisions(51)
1968 for tag in ['','C']:
1969 tc = h[params[P][0]].cd(1)
1970 if tag == 'C': tc = h[params[P][0]].cd(2)
1971 tc.SetGridx(1)
1972 h[P].DrawClone()
1973 for S in side:
1974 h[params[P][2]+S+tag].SetMarkerStyle(92)
1975 if tag=='C':
1976 h[params[P][2]+S+tag].SetMarkerSize(1)
1977 h[params[P][2]+S+tag].SetMarkerStyle(71)
1978 else: h[params[P][2]+S+tag].SetMarkerSize(2)
1979 h[params[P][2]+S+tag].SetMarkerColor(side[S][0])
1980 h[params[P][2]+S+tag].Draw('P'+side[S][1])
1981
1982 myPrint(h['Mpvs'],'Mpvs')
1983 myPrint(h['attLs'],'AttLength')
1984
1985 k=1
1986 tMax = {'CsignalT':0,'MsignalT':0,'NnSiPMs':0,'WsignalT':0,'SsignalT':0}
1987 for plane in range(5):
1988 for p in ['L','R']:
1989 for t in tMax:
1990 for bar in range(10):
1991 name = t+p+'_US'+str(s*10+plane)+'_'+str(bar)
1992 nmax = h[name].GetBinContent(h[name].GetMaximumBin())
1993 err = h[name].GetBinError(h[name].GetMaximumBin())
1994 if err/nmax > 0.2: continue
1995 if nmax > tMax[t]: tMax[t]=nmax
1996 for plane in range(5):
1997 for p in ['L','R']:
1998 for t in tMax:
1999 tc = h[t].cd(k)
2000 for bar in range(10):
2001 color = barColor[bar%10]
2002 name = t+p+'_US'+str(s*10+plane)+'_'+str(bar)
2003 h[name].SetStats(0)
2004 h[name].SetMarkerStyle(34)
2005 h[name].SetMarkerColor(color)
2006 h[name].SetLineColor(color)
2007 dropt = ''
2008 if len(tc.GetListOfPrimitives())>0: dropt='same'
2009 h[name].SetMinimum(0.)
2010 h[name].SetMaximum(tMax[t]*1.1)
2011 if name.find('CsignalT')==0:
2012 h[name].SetTitle('station '+p+' '+str(plane+1)+';chi2/nDoF')
2013 if name.find('MsignalT')==0:
2014 h[name].SetTitle('QDC sum station '+p+' '+str(plane+1)+';d [cm];QDC')
2015 elif name.find('SsignalT')==0:
2016 h[name].SetTitle('sigma '+p+' '+str(plane+1)+';d [cm];QDC')
2017 elif name.find('WsignalT')==0:
2018 h[name].SetTitle('width '+p+' '+str(plane+1)+';d [cm];QDC')
2019 if t== 'MsignalT': rc = h[name].Draw('P'+dropt)
2020 else: rc = h[name].Draw('PL'+dropt)
2021 if t == 'NnSiPMs' :
2022 tc = h['NnSiPMs1'].cd(1)
2023 h[name].SetMinimum(3)
2024 if k>1: dropt = 'same'
2025 h[name].Draw(dropt)
2026 k+=1
2027 for t in ['MsignalT','NnSiPMs','WsignalT','SsignalT']: myPrint(h[t],t)
2028 myPrint(h['NnSiPMs1'],'nSiPMsAll')
2029
2030# DS
2031 ut.bookCanvas(h,'Msignal1','',900,600,1,1)
2032#example plots:
2033 tc = h['Msignal1'].cd()
2034 name = 'signalR_DS31-X11'
2035 h[name].Draw()
2036 tc.Update()
2037 stats = h[name].FindObject('stats')
2038 stats.SetOptFit(1111111)
2039 stats.SetX1NDC(0.51)
2040 stats.SetY1NDC(0.54)
2041 stats.SetX2NDC(0.98)
2042 stats.SetY2NDC(0.94)
2043 tc.Update()
2044 myPrint(h['Msignal1'],'signalDS1')
2045
2046 ut.bookCanvas(h,'MDsignal2d','',900,600,2,4)
2047 k=1
2048 for plane in range(4):
2049 for p in ['L','R']:
2050 tc = h['MDsignal2d'].cd(k)
2051 h['signal'+p+'_DS3'+str(plane)].SetTitle('; QDC ;d [cm]')
2052 h['signal'+p+'_DS3'+str(plane)].Draw('colz')
2053 k+=1
2054 myPrint(h['MDsignal2d'],'MDsignal2d')
2055 ut.bookCanvas(h,'MDSsignal','',1200,900,1,1)
2056 plColor={0:ROOT.kBlue,1:ROOT.kGreen,2:ROOT.kCyan,3:ROOT.kOrange}
2057 pMarker={'L':22,'R':23}
2058 slopes={}
2059 for p in ['L','R']:
2060 for plane in range(4):
2061 hist = h['Msignal'+p+'_DS3'+str(plane)]
2062 if hist.GetEntries()<1: continue
2063 hist.SetStats(0)
2064 hist.SetTitle(";distance [cm]")
2065 hist.SetMarkerStyle(pMarker[p])
2066 hist.SetMarkerColor(plColor[plane])
2067 hist.SetLineColor(plColor[plane])
2068 if p=='R' or plane%2==1: rc=hist.Fit('expo','S','',0,60)
2069 if p=='L' and plane%2==0: rc=hist.Fit('expo','S','',20,80)
2070 res = rc.Get()
2071 if res: slopes[str(plane)+p]=[res.Parameter(1),res.ParError(1)]
2072 hist.GetFunction('expo').SetLineColor(plColor[plane])
2073
2074 h['MsignalR_DS30'].SetMaximum(60)
2075 h['MsignalR_DS30'].Draw()
2076 y=0.12
2077 for p in ['L','R']:
2078 for plane in range(4):
2079 hist = h['Msignal'+p+'_DS3'+str(plane)]
2080 if hist.GetEntries()<1: continue
2081 hist.Draw('same')
2082 pp = p
2083 if plane%2==1: pp='T'
2084 txt = "%i%s (%5.1F+/-%5.1F)cm"%(plane,pp,-1./slopes[str(plane)+p][0],slopes[str(plane)+p][1]/slopes[str(plane)+p][0]**2)
2085 latex.DrawLatexNDC(0.3,0+y,txt)
2086 y += 0.05
2087 myPrint(h['MDSsignal'],'DSattenuation')

◆ plotMipsGainRatio()

Survey-MufiScifi.plotMipsGainRatio ( )

Definition at line 2088 of file Survey-MufiScifi.py.

2088def plotMipsGainRatio():
2089 t = 'MsignalT'
2090 for run in [89,90,91]:
2091 h[run] = {}
2092 f = ROOT.TFile('LandauFits_run'+str(run)+'.root')
2093 for plane in range(5):
2094 for p in ['L','R']:
2095 for bar in range(10):
2096 name = t+p+'_US'+str(s*10+plane)+'_'+str(bar)
2097 h[run][name] = f.FindObjectAny(name).Clone()
2098 h[run][name].SetDirectory(ROOT.gROOT)
2099 nRef = 90 # 90 is with gain 2.5 and stable
2100 for run in [89,91]:
2101 ut.bookCanvas(h,'MsignalN'+str(run),'',1200,2000,2,5)
2102 k = 1
2103 for plane in range(5):
2104 for p in ['L','R']:
2105 nMax = 0
2106 for bar in range(10):
2107 name = t+p+'_US'+str(s*10+plane)+'_'+str(bar)
2108 h[run]['N'+name] = h[run][name].Clone('N'+name)
2109 h[run]['N'+name].Divide(h[nRef][name])
2110 tc = h['MsignalN'+str(run)].cd(k)
2111 color = barColor[bar]
2112 h[run]['N'+name].SetStats(0)
2113 h[run]['N'+name].SetMarkerStyle(34)
2114 h[run]['N'+name].SetMarkerColor(color)
2115 h[run]['N'+name].SetLineColor(color)
2116 M = h[run]['N'+name].GetBinContent(h[run]['N'+name].GetMaximumBin())
2117 if M>nMax and M<20: nMax = M
2118 for bar in range(10):
2119 name = t+p+'_US'+str(s*10+plane)+'_'+str(bar)
2120 h[run]['N'+name].SetMaximum(1.2*nMax)
2121 h[run]['N'+name].SetMinimum(0)
2122 opt = 'P'
2123 if bar>0: opt+='same'
2124 h[run]['N'+name].Draw(opt)
2125 k+=1
2126 h['MsignalN'+str(run)].Print('MPVgainRatio_'+str(run)+'_'+str(nRef)+'.png')
2127

◆ plotMipsTimeRes()

Survey-MufiScifi.plotMipsTimeRes (   readHists = True,
  animation = False 
)

Definition at line 2430 of file Survey-MufiScifi.py.

2430def plotMipsTimeRes(readHists=True,animation=False):
2431 maxRes = {1:0.6,2:1.0,3:0.6}
2432 if readHists:
2433 for x in h:
2434 if hasattr(x,'Reset'): x.Reset()
2435 ut.readHists(h,'MuFilterEff_run'+str(options.runNumber)+'.root',withProj=False)
2436 ut.bookCanvas(h,'dummy','',900,800,1,1)
2437 S = [1,2,3]
2438 if options.runNumber<100: S = [2,3]
2439 h['results']={}
2440 results = h['results']
2441 for mode in ['','fast','F1']:
2442 results[mode] = {}
2443 for s in S:
2444 results[mode][s] = {}
2445 if s==1: ut.bookCanvas(h,mode+'TimeLR-'+str(s),'',1800,800,2,1)
2446 if s==2: ut.bookCanvas(h,mode+'TimeLR-'+str(s),'',1800,1500,2,3)
2447 if s==3: ut.bookCanvas(h,mode+'TimeLR-'+str(s),'',900,1500,1,3)
2448 nbars = 0
2449 for plane in range(systemAndPlanes[s]):
2450 if systemAndOrientation(s,plane) == "vertical": continue
2451 nbars += systemAndBars[s]
2452 ut.bookHist(h,mode+'tsigma_'+str(s),';'+str(systemAndBars[s])+'*station + bar number;#sigma_{t} [ns]',nbars,-0.5,nbars-0.5)
2453 k = 1
2454 for plane in range(systemAndPlanes[s]):
2455 if systemAndOrientation(s,plane) == "vertical": continue
2456 results[mode][s][plane]={}
2457 for bar in range(systemAndBars[s]):
2458 results[mode][s][plane][bar]={}
2459 tc = h['dummy'].cd()
2460 name0 = 'dt'+mode+'LRvsX_'+sdict[s]+str(s*10+plane)+'_'+str(bar)
2461 name = 'Rdt'+mode+'LRvsX_'+sdict[s]+str(s*10+plane)+'_'+str(bar)
2462 tname = mode+'TimeRes_'+sdict[s]+str(s*10+plane)+'_'+str(bar)
2463 h[name] = h[name0].Clone(name)
2464 h[tname] = h[name].ProjectionX(tname)
2465 h[tname].Reset()
2466 h['g_'+name] = ROOT.TGraphErrors()
2467 g = h['g_'+name]
2468 np = 0
2469 xax = h[tname].GetXaxis()
2470 yax = h[tname].GetYaxis()
2471 yax.SetTitle('#sigma_{t} [ns]')
2472 ymin = h[name].GetYaxis().GetBinCenter(1)
2473 ymax = h[name].GetYaxis().GetBinCenter(h[name].GetYaxis().GetNbins())
2474 minContent = 50
2475 for j in range(1,xax.GetNbins()+1):
2476 jname = name+'-X'+str(j)
2477 h[jname] = h[name].ProjectionY(jname,j,j)
2478 if h[jname].GetSumOfWeights()<minContent: continue
2479 xx = h[jname].GetBinCenter(h[jname].GetMaximumBin())
2480 xmin = max(ymin,xx-2.)
2481 xmax = min(ymax,xx+2.)
2482 rc = h[jname].Fit('gaus','SQL','',xmin,xmax)
2483 res = rc.Get()
2484 h[tname].SetBinContent(j,res.Parameter(2))
2485 h[tname].SetBinError(j,res.ParError(2))
2486 g.SetPoint(np,xax.GetBinCenter(j),res.Parameter(1))
2487 g.SetPointError(np,xax.GetBinCenter(j),res.ParError(1))
2488 np+=1
2489 A = weightedY(h[tname])
2490 if A[1]<100:
2491 h[mode+'tsigma_'+str(s)].SetBinContent(plane*systemAndBars[s]+bar,A[0])
2492 h[mode+'tsigma_'+str(s)].SetBinError(plane*systemAndBars[s]+bar,A[1])
2493 rc = g.Fit('pol1','SQ')
2494 res = rc.Get()
2495 if not res: continue
2496 results[mode][s][plane][bar]['velocity'] = 2./res.Parameter(1)
2497
2498 if bar==0: opt=''
2499 else: opt='same'
2500 color = barColor[bar%10]
2501 h[tname].SetMaximum(maxRes[s])
2502 h[tname].SetStats(0)
2503 h[tname].SetMarkerStyle(34)
2504 h[tname].SetMarkerColor(color)
2505 h[tname].SetLineColor(color)
2506 tc = h[mode+'TimeLR-'+str(s)].cd(k)
2507 h[tname].Draw('P'+opt)
2508 h[tname].Draw('HISTLSAME')
2509 k+=1
2510
2511 ut.bookCanvas(h,mode+'dummy'+sdict[s],'',900,800,1,1)
2512 tc = h[mode+'dummy'+sdict[s]].cd()
2513 tc.SetGridx(1)
2514 name = mode+'tsigma_'+str(s)
2515 h[name].SetStats(0)
2516 h[name].SetMarkerColor(ROOT.kBlue)
2517 h[name].SetMarkerStyle(34)
2518 h[name].SetMarkerSize(2)
2519 X = h[name].GetXaxis()
2520 X.SetTitleFont(42)
2521 X.SetNdivisions(505)
2522 h[name].Draw('P')
2523 myPrint(h[mode+'dummy'+sdict[s]],mode+'TimeResol-'+sdict[s])
2524 myPrint(h[mode+'TimeLR-'+str(s)],mode+'TimeResol-'+sdict[s]+'vsX')
2525
2526 if animation:
2527 tc = h['dummy'].cd()
2528 for s in S:
2529 for l in range(systemAndPlanes[s]):
2530 if systemAndOrientation(s,l) == "vertical": continue
2531 if l==4 and s==3: continue # only 2 station for H8
2532 for bar in range(systemAndBars[s]):
2533 name = 'dtLRvsX_'+sdict[s]+str(s*10+l)+'_'+str(bar)
2534 if h[name].GetSumOfWeights()>50: h[name].Draw('colz')
2535 myPrint(h['dummy'],str(l)+'-'+str(bar).zfill(2)+'Xanimat',withRootFile=False)
2536 os.system("convert -delay 120 -loop 0 *Xanimat-run"+str(options.runNumber)+".png TDCvsX_"+sdict[s]+"-animat.gif")
2537 os.system("rm *Xanimat-run*.png")
2538
2539 mode = ''
2540 for s in S:
2541 h['galldtV_'+str(s)] = ROOT.TGraph()
2542 if s==3: ut.bookHist(h,'alldtV_'+str(s),';bar;[cm/ns]',180,0.,180.)
2543 if s==2: ut.bookHist(h,'alldtV_'+str(s),';bar;[cm/ns]',800,0.,800.)
2544 if s==1: ut.bookHist(h,'alldtV_'+str(s),';bar;[cm/ns]',15,0.,15.)
2545 ut.bookHist(h,'alldtVdis_'+str(s),';[cm/ns]',100,10.,20.)
2546 np = 0
2547 nx = 0
2548 for plane in range(systemAndPlanes[s]):
2549 if systemAndOrientation(s,plane) == "vertical": continue
2550 for bar in range(systemAndBars[s]):
2551 X = results[mode][s][plane][bar]
2552 if 'velocity' in X:
2553 V = X['velocity']
2554 if abs(V)>1:
2555 h['galldtV_'+str(s)].SetPoint(nx,np,abs(V))
2556 nx+=1
2557 rc = h['alldtVdis_'+str(s)].Fill(abs(V))
2558 np+=1
2559 ut.bookCanvas(h,'sdtV'+str(s),'',1600,800,1,1)
2560 tc = h['sdtV'+str(s)].cd()
2561 h['alldtVdis_'+str(s)].SetStats(0)
2562 h['alldtV_'+str(s)].SetStats(0)
2563 rc = h['alldtVdis_'+str(s)].Fit('gaus','SL','',10,18)
2564 res = rc.Get()
2565 txt = ROOT.TLatex()
2566 txt.DrawLatexNDC(0.3,0.5,"v=%5.2F cm/ns"%(res.Parameter(1)))
2567 myPrint(h['sdtV'+str(s)],sdict[s]+'dTsignalSpeedSum')
2568 hist = h['alldtV_'+str(s)]
2569 hist.SetMinimum(10)
2570 hist.SetMaximum(18)
2571 hist.Reset()
2572 hist.Draw()
2573 h['galldtV_'+str(s)].SetMarkerStyle(29)
2574 h['galldtV_'+str(s)].SetMarkerSize(1.5)
2575 h['galldtV_'+str(s)].SetMarkerColor(ROOT.kBlue)
2576 h['galldtV_'+str(s)].Draw('P')
2577 myPrint(h['sdtV'+str(s)],sdict[s]+'dTsignalSpeed')
2578

◆ plotMipsTimeResT0()

Survey-MufiScifi.plotMipsTimeResT0 (   readHists = True,
  animation = False 
)

Definition at line 2234 of file Survey-MufiScifi.py.

2234def plotMipsTimeResT0(readHists=True,animation=False):
2235 if readHists:
2236 for x in h:
2237 if hasattr(x,'Reset'): x.Reset()
2238 ut.readHists(h,'MuFilterEff_run'+str(options.runNumber)+'.root',withProj=False)
2239 ut.bookCanvas(h,'dummy','',900,800,1,1)
2240 S = [1,2]
2241 if options.runNumber<100: S = [2]
2242 for mode in ['t0']:
2243 for s in S:
2244 for side in ['L','R']:
2245 if s==1: ut.bookCanvas(h,mode+'Time'+side+'-'+str(s),'',1800,800,2,1)
2246 if s==2: ut.bookCanvas(h,mode+'Time'+side+'-'+str(s),'',1800,1500,2,3)
2247 nbars = 0
2248 for plane in range(systemAndPlanes[s]):
2249 if systemAndOrientation(s,plane) == "vertical": continue
2250 nbars += systemAndBars[s]
2251 ut.bookHist(h,mode+'tsigma_'+str(s),';'+str(systemAndBars[s])+'*station + bar number;#sigma_{t} [ns]',nbars,-0.5,nbars-0.5)
2252
2253 h['resultsT0'] = {}
2254 resultsT0 = h['resultsT0']
2255
2256 tc = h['dummy'].cd()
2257 mode = 't0'
2258 for s in S:
2259 resultsT0[s]={}
2260 for l in range(systemAndPlanes[s]):
2261 resultsT0[s][l]={}
2262 for side in ['L','R']:
2263 resultsT0[s][l][side]={}
2264 for bar in range(systemAndBars[s]):
2265 resultsT0[s][l][side][bar]={}
2266 for sipm in range(systemAndChannels[s][0]+systemAndChannels[s][1]):
2267 resultsT0[s][l][side][bar][sipm] = {}
2268 key = sdict[s]+str(s*10+l)+'_'+str(bar)
2269 name = 'tvsX'+side+'_'+key+'-c'+str(sipm)
2270 tname = 'T'+name
2271 h[tname] = h[name].ProjectionX(tname)
2272 h[tname].Reset()
2273 h['g_'+name] = ROOT.TGraphErrors()
2274 g = h['g_'+name]
2275 np = 0
2276 xax = h[tname].GetXaxis()
2277 yax = h[tname].GetYaxis()
2278 yax.SetTitle('#sigma_{t} [ns]')
2279 minContent = 50
2280 for j in range(1,xax.GetNbins()+1):
2281 jname = name+'-X'+str(j)
2282 h[jname] = h[name].ProjectionY(jname,j,j)
2283 if h[jname].GetEntries()<minContent: continue
2284 rc = h[jname].Fit('gaus','SQ0')
2285 res = rc.Get()
2286 if res:
2287 h[tname].SetBinContent(j,res.Parameter(2))
2288 h[tname].SetBinError(j,res.ParError(2))
2289 g.SetPoint(np,xax.GetBinCenter(j),res.Parameter(1))
2290 g.SetPointError(np,xax.GetBinCenter(j),res.ParError(1))
2291 np+=1
2292 A = weightedY(h[tname])
2293 if A[1]<100:
2294 resultsT0[s][l][side][bar][sipm]['meanSigma'] = A[0]
2295 resultsT0[s][l][side][bar][sipm]['errorSigma'] = A[1]
2296 rc = g.Fit('pol1','SQ')
2297 res = rc.Get()
2298 if not res: continue
2299 resultsT0[s][l][side][bar][sipm]['velocity'] = 1./res.Parameter(1)
2300#
2301 h['sideLines'] = []
2302 for s in resultsT0:
2303 ut.bookHist(h,'allV_'+str(s),';sipm channel;[cm/ns]',800,0.,800.)
2304 ut.bookHist(h,'allVdis_'+str(s),';[cm/ns]',100,-20.,20.)
2305 h['gallV_'+str(s)] = ROOT.TGraph()
2306 np = 0
2307 nx = 0
2308 for l in resultsT0[s]:
2309 for side in resultsT0[s][l]:
2310 for bar in resultsT0[s][l][side]:
2311 for sipm in resultsT0[s][l][side][bar]:
2312 X = resultsT0[s][l][side][bar][sipm]
2313 if 'velocity' in X:
2314 V = X['velocity']
2315 rc = h['allV_'+str(s)].Fill(np,V)
2316 if abs(V)>1:
2317 h['gallV_'+str(s)].SetPoint(nx,np,abs(V))
2318 nx+=1
2319 rc = h['allVdis_'+str(s)].Fill(V)
2320 np+=1
2321 ut.bookCanvas(h,'sV'+str(s),'',1600,800,1,1)
2322 tc = h['sV'+str(s)].cd()
2323 h['allVdis_'+str(s)].SetStats(0)
2324 h['allV_'+str(s)].SetStats(0)
2325 rc = h['allVdis_'+str(s)].Fit('gaus','S','',-18,-8)
2326 resL = rc.Get()
2327 rc = h['allVdis_'+str(s)].Fit('gaus','S+','',8,18)
2328 resR = rc.Get()
2329 txt = ROOT.TLatex()
2330 txt.DrawLatexNDC(0.3,0.5,"v=%5.2F cm/ns"%(resL.Parameter(1)))
2331 txt.DrawLatexNDC(0.3,0.4,"v=+%5.2F cm/ns"%(resR.Parameter(1)))
2332 myPrint(h['sV'+str(s)],sdict[s]+'signalSpeedSum')
2333 hist = h['allV_'+str(s)]
2334 hist.SetMinimum(10)
2335 hist.SetMaximum(18)
2336 hist.Reset()
2337 hist.Draw()
2338 h['gallV_'+str(s)].SetMarkerStyle(29)
2339 h['gallV_'+str(s)].SetMarkerSize(1.5)
2340 h['gallV_'+str(s)].SetMarkerColor(ROOT.kBlue)
2341 h['gallV_'+str(s)].Draw('P')
2342 if s == 2:
2343 xchan = 0
2344 for l in range(systemAndPlanes[s]):
2345 for side in ['L','R']:
2346 xchan+=80
2347 h['sideLines'].append(ROOT.TLine(xchan,hist.GetMinimum(),xchan,hist.GetMaximum()) )
2348 L = h['sideLines'][len(h['sideLines'])-1]
2349 L.SetLineColor(ROOT.kGreen)
2350 L.Draw()
2351
2352 myPrint(h['sV'+str(s)],sdict[s]+'signalSpeedChannel')
2353
2354 for s in resultsT0:
2355 h['gallSigm_'+str(s)] = ROOT.TGraphErrors()
2356 ut.bookHist(h,'allSigm_'+str(s),';SiPM channel ; [ns]',800,0.,800.)
2357 hist = h['allSigm_'+str(s)]
2358 ut.bookCanvas(h,'sS'+str(s),'',1600,800,1,1)
2359 np = 0
2360 nx = 0
2361 for l in resultsT0[s]:
2362 for side in resultsT0[s][l]:
2363 for bar in resultsT0[s][l][side]:
2364 for sipm in resultsT0[s][l][side][bar]:
2365 X = resultsT0[s][l][side][bar][sipm]
2366 if 'meanSigma' in X:
2367 V = X['meanSigma']
2368 E = X['errorSigma']
2369 h['gallSigm_'+str(s)].SetPoint(nx,np,V)
2370 h['gallSigm_'+str(s)].SetPointError(nx,0.5,E)
2371 nx+=1
2372 np+=1
2373# printout
2374 for s in resultsT0:
2375 for l in resultsT0[s]:
2376 for bar in resultsT0[s][l][side]:
2377 for side in resultsT0[s][l]:
2378 txt = '%s%i %i %s:'%(sdict[s],10*s+l,bar,side)
2379 for sipm in resultsT0[s][l][side][bar]:
2380 X = resultsT0[s][l][side][bar][sipm]
2381 V = 0
2382 E = 0
2383 if 'meanSigma' in X:
2384 V = X['meanSigma']
2385 E = X['errorSigma']
2386 txt += ' %5.0F+/-%5.0F '%(V*1000,E*1000)
2387 print(txt)
2388
2389 tc = h['sS'+str(s)].cd()
2390 hist.SetMaximum(1.5)
2391 hist.SetMinimum(0.)
2392 hist.SetStats(0)
2393 hist.Draw()
2394 h['gallSigm_'+str(s)].SetMarkerStyle(29)
2395 h['gallSigm_'+str(s)].SetMarkerSize(1.5)
2396 h['gallSigm_'+str(s)].SetMarkerColor(ROOT.kBlue)
2397 h['gallSigm_'+str(s)].Draw('P')
2398 if s == 2:
2399 xchan = 0
2400 for l in range(systemAndPlanes[s]):
2401 for side in ['L','R']:
2402 xchan+=80
2403 h['sideLines'].append(ROOT.TLine(xchan,0.,xchan,hist.GetMaximum()) )
2404 L = h['sideLines'][len(h['sideLines'])-1]
2405 L.SetLineColor(ROOT.kGreen)
2406 L.Draw()
2407 myPrint(h['sS'+str(s)],sdict[s]+'SigmaChannel')
2408
2409 if animation:
2410 tc = h['dummy'].cd()
2411 for s in S:
2412 for l in range(systemAndPlanes[s]):
2413 for bar in range(systemAndBars[s]):
2414 for side in ['L','R']:
2415 txt=str(s*10+l)+side+"_"+str(bar)
2416 for sipm in range(systemAndChannels[s][0]+systemAndChannels[s][1]):
2417 if s==2 and smallSiPMchannel(sipm): continue
2418 X = resultsT0[s][side][bar][sipm]
2419 if len(X)==2: txt+=" "+"%5.2F"%(X['meanSigma'] )
2420 else: txt+=" "+"%5.2F"%(-1)
2421 key = sdict[s]+str(s*10+l)+'_'+str(bar)
2422 name = 'tvsX'+side+'_'+key+'-c'+str(sipm)
2423 tname = 'T'+name
2424 h[name].Draw('colz')
2425 myPrint(h['dummy'],name+'X',withRootFile=False)
2426 print(txt, tname)
2427 os.system("convert -delay 120 -loop 0 *X.png TDCallChannels_"+sdict[s]+".gif")
2428 os.system("rm *X.png")
2429

◆ plotRMS()

Survey-MufiScifi.plotRMS (   readHists = True,
  t0_channel = 4 
)

Definition at line 2128 of file Survey-MufiScifi.py.

2128def plotRMS(readHists=True,t0_channel=4):
2129 if readHists:
2130 ut.readHists(h,options.path+"MuFilterEff_run"+str(r)+".root",withProj=False)
2131 for s in [1,2]:
2132 h['tdcCalib'+sdict[s]] = {}
2133 for l in range(systemAndPlanes[s]):
2134 for bar in range(systemAndBars[s]):
2135 for side in ['L','R']:
2136 key = sdict[s]+str(s*10+l)+'_'+str(bar)
2137 for i in range(systemAndChannels[s][1]+systemAndChannels[s][0]):
2138 if i==t0_channel: continue
2139 j = side+'_'+key+'-c'+str(i)
2140 x = 'TDCcalib'+j
2141 # get TDC calibration constants:
2142 h['tdcCalib'+sdict[s]][j] =[0,0,0,0]
2143 if h[x].GetEntries()>10:
2144 rc = h[x].Fit('gaus','SNQ')
2145 rc = h[x].Fit('gaus','SNQ')
2146 if rc:
2147 res = rc.Get()
2148 h['tdcCalib'+sdict[s]][j] =[res.Parameter(1),res.ParError(1),res.Parameter(2),res.ParError(2)]
2149#
2150 for l in range(systemAndPlanes[s]):
2151 tag = sdict[s]+str(l)
2152 ut.bookCanvas(h,'sigmaTDC_'+tag,'TDC RMS',2000,600,systemAndBars[s],2)
2153 ut.bookCanvas(h, 'TDCcalib_'+tag,'TDC TDCi-T0',2000,600,systemAndBars[s],2)
2154 ut.bookCanvas(h,'sigmaQDC_'+tag,'QDC RMS',2000,600,systemAndBars[s],2)
2155 k=1
2156 for bar in range(systemAndBars[s]):
2157 for side in ['L','R']:
2158 key = sdict[s]+str(10*s+l)+'_'+str(bar)
2159 for i in range(systemAndChannels[s][1]+systemAndChannels[s][0]):
2160 j = side+'_'+key+'-c'+str(i)
2161 for x in ['sigmaTDC_'+tag,'TDCcalib_'+tag,'sigmaQDC_'+tag]:
2162 if x.find('calib')>0 and (i==t0_channel): continue
2163 tc=h[x].cd(k)
2164 opt='same'
2165 if i==0 and x.find('calib')<0: opt=""
2166 if i==1 and x.find('calib')>0: opt=""
2167 aHist = x.split('_')[0]+j
2168 h[aHist].SetLineColor(barColor[i])
2169 if x.find('QDC')>0: h[aHist].GetXaxis().SetRangeUser(-15,15)
2170 else: h[aHist].GetXaxis().SetRangeUser(-5,5)
2171 h[aHist].Draw(opt)
2172 h[aHist].Draw(opt+'HIST')
2173 k+=1
2174 myPrint(h['sigmaTDC'+tag],'TDCrms'+tag)
2175 myPrint(h['TDCcalib'+tag],'TDCcalib'+tag)
2176 myPrint(h['sigmaQDC'+tag],'QDCrms'+tag)
2177
2178 ut.bookHist(h,'TDCcalibMean_'+tag,';SiPM channel ; [ns]',160,0.,160.)
2179 ut.bookHist(h,'TDCcalibSigma_'+tag,';SiPM channel ; [ns]',160,0.,160.)
2180 h['gTDCcalibMean_'+tag]=ROOT.TGraphErrors()
2181 h['gTDCcalibSigma_'+tag]=ROOT.TGraphErrors()
2182
2183 x = 'TDCcalib'+sdict[s]
2184 tmp = h['tdcCalib'+sdict[s]][x]
2185 side = 0
2186 if x[0]=='R': side = 1
2187 l = x[5:6]
2188 bar = x[7:8]
2189 c = x[10:11]
2190 xbin = int(bar)*16+side*8+int(c)
2191 tag = sdict[s]+"_"+str(l)
2192 rc = h['TDCcalibMean_'+tag].SetBinContent(xbin+1,tmp[0])
2193 rc = h['TDCcalibMean_'+tag].SetBinError(xbin+1,tmp[1])
2194 rc = h['TDCcalibSigma_'+tag].SetBinContent(xbin+1,tmp[2])
2195 rc = h['TDCcalibSigma_'+tag].SetBinError(xbin+1,tmp[3])
2196 #print("%s %5.2F+/-%5.2F %5.2F+/-%5.2F"%(x,tmp[0],tmp[1],tmp[2],tmp[3]))
2197 ut.bookCanvas(h,'tTDCcalib'+sdict[s],'TDC calib',2400,1800,2,5)
2198 for l in range(systemAndPlanes[s]):
2199 tag = sdict[s]+"_"+str(l)
2200 aHistS = h['TDCcalibSigma_'+tag]
2201 aHistM = h['TDCcalibMean_'+tag]
2202 k=0
2203 for i in range(1,aHistS.GetNbinsX()):
2204 if aHistS.GetBinContent(i)>0:
2205 h['gTDCcalibSigma_'+tag].SetPoint(k,i-1,aHistS.GetBinContent(i))
2206 h['gTDCcalibSigma_'+tag].SetPointError(k,0.5,aHistS.GetBinError(i))
2207 h['gTDCcalibMean_'+tag].SetPoint(k,i-1,aHistM.GetBinContent(i))
2208 h['gTDCcalibMean_'+tag].SetPointError(k,0.5,aHistM.GetBinError(i))
2209 k+=1
2210
2211 planeColor = {0:ROOT.kBlue,1:ROOT.kRed,2:ROOT.kGreen,3:ROOT.kCyan,4:ROOT.kMagenta}
2212 for l in range(systemAndPlanes[s]):
2213 tag = sdict[s]+"_"+str(l)
2214 tc = h['tTDCcalib'].cd(2*l+1)
2215 aHistS = h['TDCcalibSigma_'+tag]
2216 aHistS.Reset()
2217 aHistS.SetMaximum(2.0)
2218 aHistS.Draw()
2219 h['gTDCcalibSigma_'+tag].SetLineColor(planeColor[l])
2220 h['gTDCcalibSigma_'+tag].Draw('same')
2221
2222 for l in range(systemAndPlanes[s]):
2223 tag = sdict[s]+"_"+str(l)
2224 tc = h['tTDCcalib'+sdict[s]]+cd(2*l+2)
2225 aHistM = h['TDCcalibMean_'+tag]
2226 aHistM.Reset()
2227 aHistM.SetMaximum(2.0)
2228 aHistM.SetMinimum(-2.0)
2229 aHistM.Draw()
2230 h['gTDCcalibMean_'+tag].SetLineColor(planeColor[l])
2231 h['gTDCcalibMean_'+tag].Draw('same')
2232 myPrint(h['tTDCcalib'+sdict[s]],'TDCcalibration'+sdict[s])
2233

◆ plotsForCollabMeeting()

Survey-MufiScifi.plotsForCollabMeeting ( )

Definition at line 4021 of file Survey-MufiScifi.py.

4021def plotsForCollabMeeting():
4022 drawArea(s=3,p=0,opt='',color=ROOT.kGreen)
4023 h['track_DS30'].Draw('colzsame')
4024 myPrint(h['area'],'Extrap2DS30')
4025 name = 'track_DS30_projx'
4026 h[name] = h['track_DS30'].ProjectionX(name)
4027 h[name] .Fit('gaus')
4028 h['area'].Update()
4029 stats = h[name].FindObject('stats')
4030 stats.SetOptStat(1000000000)
4031 stats.SetOptFit(1111111)
4032 stats.SetX1NDC(0.62)
4033 stats.SetY1NDC(0.66)
4034 stats.SetX2NDC(0.98)
4035 stats.SetY2NDC(0.94)
4036 h['area'].Update()
4037 myPrint(h['area'],'Extrap2DS30_projx')
4038#
4039 drawArea(s=1,p=0,opt='',color=ROOT.kGreen)
4040 myPrint(h['area'],'Extrap2Veto10')
4041#
4042 ut.bookCanvas(h,'DSXRes','',2000,600,4,1)
4043 L = {1:31,2:33,3:35,4:36}
4044 for n in L:
4045 tc = h['DSXRes'].cd(n)
4046 histo = ROOT.gROOT.FindObject('resX_DS'+str(L[n])+'_proj')
4047 histo.SetTitle('DS'+str(L[n])+';X [cm]')
4048 histo.Draw()
4049 tc.Update()
4050 stats = histo.FindObject('stats')
4051 stats.SetOptStat(1000000000)
4052 myPrint(h['DSXRes'],'resDSX')
4053 ut.bookCanvas(h,'DSYRes','',1500,600,3,1)
4054 L = {1:30,2:32,3:34}
4055 for n in L:
4056 tc = h['DSYRes'].cd(n)
4057 histo = ROOT.gROOT.FindObject('resY_DS'+str(L[n])+'_proj')
4058 histo.SetTitle('DS'+str(L[n])+';Y [cm]')
4059 histo.Draw()
4060 tc.Update()
4061 stats = histo.FindObject('stats')
4062 stats.SetOptStat(1000000000)
4063 if n==4:
4064 rc = histo.Fit('gaus','SQ','',-2.,2.)
4065 myPrint(h['DSYRes'],'resDSY')
4066#
4067 ut.bookCanvas(h,'USYRes','',2000,600,5,1)
4068 L = {1:20,2:21,3:22,4:23,5:24}
4069 for n in L:
4070 tc = h['USYRes'].cd(n)
4071 histo = ROOT.gROOT.FindObject('resY_US'+str(L[n])+'_proj')
4072 histo.SetTitle('US'+str(L[n])+';Y [cm]')
4073 histo.Draw()
4074 tc.Update()
4075 if n==3 or n==4:
4076 histo.SetStats(0)
4077 else:
4078 stats = histo.FindObject('stats')
4079 stats.SetOptStat(1000000000)
4080 stats.SetX1NDC(0.54)
4081 stats.SetY1NDC(0.77)
4082 stats.SetX2NDC(0.98)
4083 stats.SetY2NDC(0.94)
4084 myPrint(h['USYRes'],'resUSY')
4085#
4086 ut.bookCanvas(h,'VEYRes','',800,600,2,1)
4087 L = {1:10,2:11}
4088 for n in L:
4089 tc = h['VEYRes'].cd(n)
4090 histo = ROOT.gROOT.FindObject('resY_Veto'+str(L[n])+'_proj')
4091 histo.SetTitle('Veto'+str(L[n])+';Y [cm]')
4092 histo.Draw()
4093 tc.Update()
4094 stats = histo.FindObject('stats')
4095 stats.SetOptStat(1000000000)
4096 stats.SetX1NDC(0.54)
4097 stats.SetY1NDC(0.77)
4098 stats.SetX2NDC(0.98)
4099 stats.SetY2NDC(0.94)
4100 tc.Update()
4101 myPrint(h['VEYRes'],'resVetoY')
4102#
4103 S = {1:'TVE',2:'TUS',3:'TDS'}
4104 proj = 'X'
4105 for s in S:
4106 sy = sdict[s]
4107 tname = S[s]
4108 k=1
4109 for plane in range(systemAndPlanes[s]):
4110 if s==3 and (plane%2==1 or plane>5): continue
4111 tc = h[tname].cd(k)
4112 k+=1
4113 key = sy+str(s*10+plane)
4114 hist = h['dtLRvsX_'+key]
4115 hist.SetTitle(key+';X [cm];#Deltat [ns]')
4116 hist.Draw('colz')
4117 if hist.GetSumOfWeights()==0: continue
4118 g = h['gdtLRvsX_'+key]
4119 g.Draw('same')
4120 xmin = hist.GetXaxis().GetXmin()
4121 xmax = hist.GetXaxis().GetXmax()
4122 rc = g.Fit('pol1','SQ','',xmin,xmax)
4123 g.GetFunction('pol1').SetLineColor(ROOT.kGreen)
4124 g.GetFunction('pol1').SetLineWidth(2)
4125 result = rc.Get()
4126 if not result: continue
4127 m = 1./result.Parameter(1)
4128 b = -result.Parameter(0) * m
4129 sign = '+'
4130 if b<0: sign = '-'
4131 txt = 'X = %5.1F #frac{cm}{ns} #times #frac{1}{2}#Deltat %s %5.1F '%(m*2,sign,abs(b))
4132 latex.SetTextSize(0.07)
4133 latex.DrawLatexNDC(0.2,0.8,txt)
4134 myPrint(h['TVE'],'dTvsX_Veto')
4135 myPrint(h['TUS'],'dTvsX_US')
4136 myPrint(h['TDS'],'dTvsX_DS')
4137#

◆ plotSignalsUS()

Survey-MufiScifi.plotSignalsUS ( )

Definition at line 2651 of file Survey-MufiScifi.py.

2651def plotSignalsUS():
2652 ut.bookCanvas(h,'signalUS','',1200,1800,4,5)
2653 k = 1
2654 s = 2
2655 for plane in range(systemAndPlanes[s]):
2656 for side in ['L','R']:
2657 for sipm in ['','S']:
2658 tc=h['signalUS'].cd(k)
2659 k+=1
2660 name = 'signalT'+sipm+side+'_'+sdict[s]+str(s*10+plane)
2661 h[name] = h[name+'_0'].ProjectionX(name)
2662 for bar in range(1,systemAndBars[s]):
2663 h[name].Add(h[name+'_'+str(bar)].ProjectionX())
2664 X = h[name].GetXaxis()
2665 X.SetRangeUser(1.5,X.GetXmax())
2666 h[name].Draw()
2667 myPrint(h['signalUS'],'signalPlane'+sdict[s])
2668
2669

◆ printScifi_residuals()

Survey-MufiScifi.printScifi_residuals (   tag = 'v0')

Definition at line 3695 of file Survey-MufiScifi.py.

3695def printScifi_residuals(tag='v0'):
3696 P = {'':'','X':'colz','Y':'colz','C':'colz'}
3697 for proj in P:
3698 myPrint(h['scifiRes'+proj],tag+'-scifiRes'+proj)
3699 for p in h['globalPos']:
3700 myPrint(h[p],tag+'-scifiRes'+p)
3701 for p in h['globalPosM']:
3702 myPrint(h[p+'M'],tag+'-scifiResM'+p)
3703# make report about alignment constants
3704 for s in range(1, nScifi+1):
3705 for o in range(2):
3706 mean = []
3707 for m in range(nMats):
3708 x = "Scifi/LocM"+str(s*100+o*10+m)
3709 mean.append(Scifi.GetConfParF(x)/u.um)
3710 XM = sum(mean)/nMats
3711 print("mean value %8.2F" %XM, " delta s:", ["%8.2F" %(mean[i]-XM) for i in range(nMats)])
3712# for H6 beam
3713 h6 = False
3714 if h6:
3715 h['trackSlopes'].GetYaxis().SetRangeUser(-20,20)
3716 h['trackSlopes'].GetXaxis().SetRangeUser(-40,0)
3717 ut.bookCanvas(h,'beamSpot','track slopes',750,750,1,1)
3718 tc = h['beamSpot'].cd()
3719 h['trackSlopes'].Draw('colz')
3720 myPrint(h['beamSpot'],tag+'-beamSpot')
3721

◆ printScifiAlignment()

Survey-MufiScifi.printScifiAlignment ( )

Definition at line 3722 of file Survey-MufiScifi.py.

3722def printScifiAlignment():
3723 if nMats ==1 : plane = ['0', '1']
3724 if nMats ==3 : plane = ['']
3725 for s in range(1, nScifi+1):
3726 for o in range(2):
3727 p = 10*s+o
3728 if nMats ==1: txt = " c.Scifi.LocM"+str(p)+"0 = "
3729 if nMats ==3: txt = " c.Scifi.LocM"+str(p)+"0,c.Scifi.LocM"+str(p)+"1,c.Scifi.LocM"+str(p)+"2 = "
3730 for m in range(nMats):
3731 x = "Scifi/LocM"+str(s*100+o*10+m)
3732 txt += "%8.2F*u.um,"%(Scifi.GetConfParF(x)/u.um)
3733 l = len(txt)
3734 print(txt[:l-1])
3735 for s in range(1, nScifi+1):
3736 for p in plane:
3737 txt = " c.Scifi.RotPhiS"+str(s)+p+",c.Scifi.RotPsiS"+str(s)+p+",c.Scifi.RotThetaS"+str(s)+p+" = "
3738 for a in ["RotPhiS","RotPsiS","RotThetaS"]:
3739 x = "Scifi/"+a+str(s)+p
3740 txt += "%8.2F*u.mrad,"%(Scifi.GetConfParF(x)/u.mrad)
3741 l = len(txt)
3742 print(txt[:l-1])
3743

◆ pyExit()

Survey-MufiScifi.pyExit ( )

Definition at line 58 of file Survey-MufiScifi.py.

58def pyExit():
59 print("Make suicide until solution found for freezing")
60 os.system('kill '+str(os.getpid()))

◆ scifi_beamSpot()

Survey-MufiScifi.scifi_beamSpot ( )

Definition at line 3299 of file Survey-MufiScifi.py.

3299def scifi_beamSpot():
3300 A,B = ROOT.TVector3(),ROOT.TVector3()
3301 ut.bookHist(h,'bs','beam spot',100,-100.,10.,100,0.,80.)
3302 for event in eventTree:
3303 xMean = 0
3304 yMean = 0
3305 w=0
3306 for d in event.Digi_ScifiHits:
3307 detID = d.GetDetectorID()
3308 s = int(detID/1000000)
3309 modules['Scifi'].GetSiPMPosition(detID,A,B)
3310 vertical = int(detID/100000)%10==1
3311 if vertical: xMean+=A[0]
3312 else: yMean+=A[1]
3313 w+=1
3314 rc = h['bs'].Fill(xMean/w,yMean/w)
3315

◆ Scifi_hitMaps()

Survey-MufiScifi.Scifi_hitMaps (   Nev = options.nEvents)

Definition at line 405 of file Survey-MufiScifi.py.

405def Scifi_hitMaps(Nev=options.nEvents):
406 Scifi = geo.modules['Scifi']
407 nScifi = Scifi.GetConfParI("Scifi/nscifi")
408 nMats = Scifi.GetConfParI("Scifi/nmats")
409 nMats = Scifi.GetConfParI("Scifi/nmats")
410 A=ROOT.TVector3()
411 B=ROOT.TVector3()
412
413 for s in range(nScifi*2):
414 ut.bookHist(h,'posX_'+str(s),'x',2000,-100.,100.)
415 ut.bookHist(h,'posY_'+str(s),'y',2000,-100.,100.)
416 if s%2==1: ut.bookHist(h,'mult_'+str(s),'mult vertical station '+str(s//2+1),100,-0.5,99.5)
417 else: ut.bookHist(h,'mult_'+str(s),'mult horizontal station '+str(s//2+1),100,-0.5,99.5)
418 for mat in range(nScifi*nMats*2):
419 ut.bookHist(h,'mat_'+str(mat),'hit map / mat',512,-0.5,511.5)
420 ut.bookHist(h,'sig_'+str(mat),'signal / mat',200,-50.0,150.)
421 ut.bookHist(h,'tdc_'+str(mat),'tdc / mat',200,-1.,4.)
422 N=-1
423 if Nev < 0 : Nev = eventTree.GetEntries()
424 for event in eventTree:
425 N+=1
426 if N%options.heartBeat == 0: print('event ',N,' ',time.ctime())
427 if N>Nev: break
428 mult = [0]*10
429 for aHit in event.Digi_ScifiHits:
430 if not aHit.isValid(): continue
431 X = Scifi_xPos(aHit.GetDetectorID())
432 rc = h['mat_'+str(X[0]*3+X[1])].Fill(X[2])
433 rc = h['sig_'+str(X[0]*3+X[1])].Fill(aHit.GetSignal(0))
434 rc = h['tdc_'+str(X[0]*3+X[1])].Fill(aHit.GetTime(0))
435 Scifi.GetSiPMPosition(aHit.GetDetectorID(),A,B)
436 if aHit.isVertical(): rc = h['posX_'+str(X[0])].Fill(A[0])
437 else: rc = h['posY_'+str(X[0])].Fill(A[1])
438 mult[X[0]]+=1
439 for s in range(10):
440 rc = h['mult_'+str(s)].Fill(mult[s])
441
442 ut.bookCanvas(h,'hitmaps',' ',1024,768,6,5)
443 ut.bookCanvas(h,'signal',' ',1024,768,6,5)
444 ut.bookCanvas(h,'tdc',' ',1024,768,6,5)
445 for mat in range(30):
446 tc = h['hitmaps'].cd(mat+1)
447 A = h['mat_'+str(mat)].GetSumOfWeights()/512.
448 if h['mat_'+str(mat)].GetMaximum()>10*A: h['mat_'+str(mat)].SetMaximum(10*A)
449 h['mat_'+str(mat)].Draw()
450 tc = h['signal'].cd(mat+1)
451 h['sig_'+str(mat)].Draw()
452 tc = h['tdc'].cd(mat+1)
453 h['tdc_'+str(mat)].Draw()
454
455 ut.bookCanvas(h,'positions',' ',2048,768,5,2)
456 ut.bookCanvas(h,'mult',' ',2048,768,5,2)
457 for s in range(5):
458 tc = h['positions'].cd(s+1)
459 h['posY_'+str(2*s)].Draw()
460 tc = h['positions'].cd(s+6)
461 h['posX_'+str(2*s+1)].Draw()
462 tc = h['mult'].cd(s+1)
463 tc.SetLogy(1)
464 h['mult_'+str(2*s)].Draw()
465 tc = h['mult'].cd(s+6)
466 tc.SetLogy(1)
467 h['mult_'+str(2*s+1)].Draw()
468
469 for canvas in ['hitmaps','signal','mult']:
470 h[canvas].Update()
471 myPrint(h[canvas],"Scifi-"+canvas)
472
Scifi::GetConfParI
Int_t GetConfParI(TString name)
Definition Scifi.h:51

◆ Scifi_residuals()

Survey-MufiScifi.Scifi_residuals (   Nev = options.nEvents,
  NbinsRes = 100,
  xmin = -2000.,
  alignPar = False,
  unbiased = True,
  minEnergy = False,
  remakeClusters = options.remakeScifiClusters,
  nproc = 1 
)

Definition at line 3420 of file Survey-MufiScifi.py.

3420def Scifi_residuals(Nev=options.nEvents,NbinsRes=100,xmin=-2000.,alignPar=False,unbiased = True,minEnergy=False,remakeClusters=options.remakeScifiClusters,nproc=1):
3421 projs = {1:'V',0:'H'}
3422 scifi = geo.modules['Scifi']
3423 nScifi = scifi.GetConfParI("Scifi/nscifi")
3424 nMats = scifi.GetConfParI("Scifi/nmats")
3425
3426 for s in range(1, nScifi+1):
3427 for o in range(2):
3428 for p in projs:
3429 proj = projs[p]
3430 xmax = -xmin
3431 ut.bookHist(h,'res'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax)
3432 ut.bookHist(h,'resX'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,100,-50.,0.)
3433 ut.bookHist(h,'resY'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,100,10.,60.)
3434 ut.bookHist(h,'resC'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,128*4*nMats+10,-0.5,128*4*nMats-0.5+10)
3435 ut.bookHist(h,'track_Scifi'+str(s*10+o),'track x/y '+str(s*10+o),80,-70.,10.,80,0.,80.)
3436 ut.bookHist(h,'trackChi2/ndof','track chi2/ndof vs ndof',100,0,100,20,0,20)
3437 ut.bookHist(h,'trackSlopes','track slope; x [mrad]; y [mrad]',1000,-100,100,1000,-100,100)
3438 parallelToZ = ROOT.TVector3(0., 0., 1.)
3439
3440 if Nev < 0 : Nev = eventTree.GetEntries()
3441# set alignment parameters
3442 if alignPar:
3443 for x in alignPar:
3444 if not x.find('Rot')>0:
3445 Scifi.SetConfPar(x,alignPar[x]*u.um)
3446 else:
3447 Scifi.SetConfPar(x,alignPar[x]*u.mrad)
3448
3449 process = []
3450 for iproc in range(nproc):
3451 try:
3452 if nproc>1: pid = os.fork()
3453 else: pid = 0
3454 except OSError:
3455 print("Could not create a child process")
3456 if pid!=0:
3457 process.append(pid)
3458 else:
3459 dN = Nev//nproc
3460 nstart = iproc*dN
3461 nstop = nstart + dN
3462 if iproc==(nproc-1): nstop = Nev
3463 for k in range(nstart,nstop):
3464 event = eventTree
3465 eventTree.GetEvent(k)
3466 if minEnergy:
3467 trackTypes = [0,0,0,0,0]
3468 for aTrack in Reco_MuonTracks:
3469 if aTrack.GetUniqueID()==1:
3470 trackTypes[1]+=1
3471 fstate = aTrack.getFittedState()
3472 mom = fstate.getMom()
3473 if abs(mom.X()/mom.Z())<0.1:
3474 # check for muon hits
3475 mult = {}
3476 for i in range(30,40): mult[i]=0
3477 for aHit in eventTree.Digi_MuFilterHits:
3478 if not aHit.isValid(): continue
3479 detID = aHit.GetDetectorID()
3480 s = detID//10000
3481 if s<3: continue
3482 l = (detID%10000)//1000 # plane number
3483 bar = (detID%1000)
3484 if s>2:
3485 l=2*l
3486 if bar>59:
3487 bar=bar-60
3488 if l<6: l+=1
3489 mult[s*10+l]+=1
3490 if mult[35]==1 and mult[34]==1: trackTypes[4]+=1 # from IP1 and high momentum
3491
3492 if aTrack.GetUniqueID()==3: trackTypes[3]+=1
3493 if not (trackTypes[4]==1): continue
3494
3495# required to bypass memory leak for files with tracks
3496 if event.GetBranch("Reco_MuonTracks"):
3497 for aTrack in event.Reco_MuonTracks: aTrack.Delete()
3498
3499 if not hasattr(event,"Cluster_Scifi") or alignPar or remakeClusters:
3500 if hasattr(trackTask,"clusScifi"):
3501 trackTask.clusScifi.Clear()
3502 trackTask.scifiCluster()
3503 clusters = trackTask.clusScifi
3504 else:
3505 clusters = event.Cluster_Scifi
3506
3507 sortedClusters={}
3508 for aCl in clusters:
3509 so = aCl.GetFirst()//100000
3510 if not so in sortedClusters: sortedClusters[so]=[]
3511 sortedClusters[so].append(aCl)
3512
3513# select events with clusters in each plane
3514 if len(sortedClusters)<2*nScifi: continue
3515 goodEvent = True
3516 for s in sortedClusters:
3517 if len(sortedClusters[s])>1: goodEvent=False
3518 if not goodEvent: continue
3519 validTrack = True
3520 for s in range(1, nScifi+1):
3521# build trackCandidate
3522 hitlist = {}
3523 if unbiased or s==1:
3524 k=0
3525 for so in sortedClusters:
3526 if so//10 == s and unbiased: continue
3527 for x in sortedClusters[so]:
3528 hitlist[k] = x
3529 k+=1
3530 theTrack = trackTask.fitTrack(hitlist)
3531 if not hasattr(theTrack,"getFittedState"):
3532 validTrack = False
3533 continue
3534 fitStatus = theTrack.getFitStatus()
3535 if not fitStatus.isFitConverged() or theTrack.getNumPointsWithMeasurement()<2*(nScifi-1):
3536 theTrack.Delete()
3537 validTrack = False
3538 continue
3539 rc = h['trackChi2/ndof'].Fill(fitStatus.getChi2()/(fitStatus.getNdf()+1E-10),fitStatus.getNdf() )
3540 fstate = theTrack.getFittedState()
3541 mom = fstate.getMom()
3542 rc = h['trackSlopes'].Fill(mom.X()/mom.Z()*1000,mom.Y()/mom.Z()*1000)
3543# check residuals
3544 if not validTrack and not unbiased: break
3545# test plane
3546 for o in range(2):
3547 testPlane = s*10+o
3548 z = zPos['Scifi'][testPlane]
3549 rep = ROOT.genfit.RKTrackRep(13)
3550 state = ROOT.genfit.StateOnPlane(rep)
3551# find closest track state
3552 mClose = 0
3553 mZmin = 999999.
3554 for m in range(0,theTrack.getNumPointsWithMeasurement()):
3555 st = theTrack.getFittedState(m)
3556 Pos = st.getPos()
3557 if abs(z-Pos.z())<mZmin:
3558 mZmin = abs(z-Pos.z())
3559 mClose = m
3560 if mZmin>10000:
3561 print("something wrong here with measurements",mClose,mZmin,theTrack.getNumPointsWithMeasurement())
3562 continue
3563 fstate = theTrack.getFittedState(mClose)
3564 pos,mom = fstate.getPos(),fstate.getMom()
3565 rep.setPosMom(state,pos,mom)
3566 NewPosition = ROOT.TVector3(0., 0., z) # assumes that plane in global coordinates is perpendicular to z-axis, which is not true for TI18 geometry.
3567 rep.extrapolateToPlane(state, NewPosition, parallelToZ )
3568 pos = state.getPos()
3569 xEx,yEx = pos.x(),pos.y()
3570 rc = h['track_Scifi'+str(testPlane)].Fill(xEx,yEx)
3571 for aCl in sortedClusters[testPlane]:
3572 aCl.GetPosition(A,B)
3573 detID = aCl.GetFirst()
3574 channel = detID%1000 + ((detID%10000)//1000)*128 + (detID%100000//10000)*512
3575# calculate DOCA
3576 pq = A-pos
3577 uCrossv= (B-A).Cross(mom)
3578 doca = pq.Dot(uCrossv)/uCrossv.Mag()
3579 rc = h['resC'+projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,channel)
3580 rc = h['res'+projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um)
3581 rc = h['resX'+projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,xEx)
3582 rc = h['resY'+projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,yEx)
3583
3584 if unbiased: theTrack.Delete()
3585 if not unbiased and validTrack: theTrack.Delete()
3586
3587 if nproc>1:
3588 ut.writeHists(h,'tmp_'+str(iproc))
3589 exit(0)
3590
3591 if nproc>1:
3592 while process:
3593 pid,exit_code = os.wait()
3594 if pid == 0: time.sleep(100)
3595 else:
3596 process.remove(pid)
3597 for i in range(nproc):
3598 ut.readHists(h,'tmp_'+str(i))
3599
3600# analysis and plots
3601 P = {'':'','X':'colz','Y':'colz','C':'colz'}
3602 h['globalPos'] = {'meanH':ROOT.TGraphErrors(),'sigmaH':ROOT.TGraphErrors(),'meanV':ROOT.TGraphErrors(),'sigmaV':ROOT.TGraphErrors()}
3603 h['globalPosM'] = {'meanH':ROOT.TGraphErrors(),'sigmaH':ROOT.TGraphErrors(),'meanV':ROOT.TGraphErrors(),'sigmaV':ROOT.TGraphErrors()}
3604 globalPos = h['globalPos']
3605 # params for the double symmetrical Crystal Ball - symetric core and tails
3606 x = ROOT.RooRealVar("x", "residual [cm]", xmin, xmax)
3607 a1 = ROOT.RooRealVar("alpha", "alpha", 2, 1e-3, 10.)
3608 p1 = ROOT.RooRealVar("n", "n", 1, 1e-3, 10.)
3609 for proj in P:
3610 ut.bookCanvas(h,'scifiRes'+proj,'',1600,1900,2,5)
3611 k=1
3612 j = {0:0,1:0}
3613 for s in range(1, nScifi+1):
3614 for o in range(2):
3615 so = s*10+o
3616 tc = h['scifiRes'+proj].cd(k)
3617 k+=1
3618 hname = 'res'+proj+projs[o]+'_Scifi'+str(so)
3619 h[hname].Draw(P[proj])
3620 if proj == '':
3621 #fit using a double sided CB shape with symmetric tails
3622 mu = ROOT.RooRealVar("mu", "mu", h[hname].GetMean(), h[hname].GetMean()-200, h[hname].GetMean()+200)
3623 width = ROOT.RooRealVar("width", "width", h[hname].GetRMS(), 1e-5, h[hname].GetRMS()+1e+3)
3624 Par = {'mean':mu,'sigma':width}
3625 dcbPdf = ROOT.RooCrystalBall("dcbPdf", "DoubleSidedCB", x, mu, width, a1, p1, True)
3626 roofit_hist = ROOT.RooDataHist("roofit_hist", "roofit_hist", ROOT.RooArgList(x), ROOT.RooFit.Import(h[hname]))
3627 fitResult = dcbPdf.fitTo(roofit_hist, ROOT.RooFit.PrintLevel(-1), ROOT.RooFit.Save(True))
3628 pl = x.frame()
3629 roofit_hist.plotOn(pl)
3630 dcbPdf.plotOn(pl)
3631 pl.Draw()
3632 latex.DrawLatexNDC(0.13,0.8, "#mu = {:.3f} +/- {:.3f}".format(mu.getVal(), mu.getError()))
3633 latex.DrawLatexNDC(0.13,0.7,"#sigma = {:.3f} +/- {:.3f}".format(width.getVal(), width.getError()))
3634 tc.Update()
3635 for p in Par:
3636 globalPos[p+projs[o]].SetPoint(s-1,s,Par[p].getVal())
3637 globalPos[p+projs[o]].SetPointError(s-1,0.5,Par[p].getError())
3638 globalPos[p+projs[o]].SetMarkerStyle(21)
3639 globalPos[p+projs[o]].SetMarkerColor(ROOT.kBlue)
3640 if proj == 'C':
3641 for m in range(nMats):
3642 h[hname+str(m)] = h[hname].ProjectionX(hname+str(m),m*512,m*512+512)
3643 mu = ROOT.RooRealVar("mu", "residual [cm]", h[hname+str(m)].GetMean(), h[hname+str(m)].GetMean()-200, h[hname+str(m)].GetMean()+200)
3644 width = ROOT.RooRealVar("width", "width", h[hname+str(m)].GetRMS(), 1e-5, h[hname+str(m)].GetRMS()+1e+3)
3645 dcbPdf = ROOT.RooCrystalBall("dcbPdf", "DoubleSidedCB", x, mu, width, a1, p1, True)
3646 roofit_hist = ROOT.RooDataHist("roofit_hist", "roofit_hist", ROOT.RooArgList(x), ROOT.RooFit.Import(h[hname+str(m)]))
3647 fitResult = dcbPdf.fitTo(roofit_hist, ROOT.RooFit.PrintLevel(-1), ROOT.RooFit.Save(True))
3648 if not fitResult: continue
3649 pl = x.frame()
3650 roofit_hist.plotOn(pl)
3651 dcbPdf.plotOn(pl)
3652 pl.Draw()
3653 for p in Par:
3654 h['globalPosM'][p+projs[o]].SetPoint(j[o], s*10+m,Par[p].getVal())
3655 h['globalPosM'][p+projs[o]].SetPointError(j[o],0.5,Par[p].getError())
3656 j[o]+=1
3657 h['globalPosM'][p+projs[o]].SetMarkerStyle(21)
3658 h['globalPosM'][p+projs[o]].SetMarkerColor(ROOT.kBlue)
3659
3660 S = ctypes.c_double()
3661 M = ctypes.c_double()
3662 alignResult = {}
3663 for p in globalPos:
3664 ut.bookCanvas(h,p,p,750,750,1,1)
3665 tc = h[p].cd()
3666 if p.find('mean')<0: globalPos[p].SetTitle(p+';station; #sigma [#mum]')
3667 else: globalPos[p].SetTitle(p+';station; offset [#mum]')
3668 globalPos[p].Draw("ALP")
3669 if p.find('mean')==0:
3670 for n in range(globalPos[p].GetN()):
3671 rc = globalPos[p].GetPoint(n,S,M)
3672 E = globalPos[p].GetErrorY(n)
3673 print("station %i: offset %s = %5.2F um sigma = %5.2F um"%(S.value,p[4:5],M.value,E))
3674 s = int(S.value*10)
3675 if p[4:5] == "V": s+=1
3676 alignResult["Scifi/LocD"+str(s)] = [M.value,E]
3677
3678 for p in h['globalPosM']:
3679 ut.bookCanvas(h,p+'M',p,750,750,1,1)
3680 tc = h[p+'M'].cd()
3681 if p.find('mean')<0: h['globalPosM'][p].SetTitle(p+';mat ; #sigma [#mum]')
3682 else: h['globalPosM'][p].SetTitle(p+';mat ; offset [#mum]')
3683 h['globalPosM'][p].Draw("ALP")
3684 if p.find('mean')==0:
3685 for n in range(h['globalPosM'][p].GetN()):
3686 rc = h['globalPosM'][p].GetPoint(n,S,M)
3687 E = h['globalPosM'][p].GetErrorY(n)
3688 print("station %i: offset %s = %5.2F um sigma = %5.2F um"%(S.value,p[4:5],M.value,E))
3689 s = int(S.value*10)
3690 if p[4:5] == "V": s+=1
3691 alignResult["Scifi/LocM"+str(s)] = [M.value,E]
3692
3693 return alignResult
3694
Scifi::SetConfPar
void SetConfPar(TString name, Float_t value)
Definition Scifi.h:46

◆ Scifi_slopes()

Survey-MufiScifi.Scifi_slopes (   Nev = options.nEvents)

Definition at line 3117 of file Survey-MufiScifi.py.

3117def Scifi_slopes(Nev=options.nEvents):
3118 A,B = ROOT.TVector3(),ROOT.TVector3()
3119 ut.bookHist(h,'slopesX','slope diffs',1000,-1.0,1.0)
3120 ut.bookHist(h,'slopesY','slope diffs',1000,-1.0,1.0)
3121 ut.bookHist(h,'clX','cluster size',10,0.5,10.5)
3122 ut.bookHist(h,'clY','cluster size',10,0.5,10.5)
3123# assuming cosmics make straight line
3124 if Nev < 0 : Nev = eventTree.GetEntries()
3125 for event in eventTree:
3126 if Nev<0: break
3127 Nev=Nev-1
3128 clusters = []
3129 hitDict = {}
3130 for k in range(event.Digi_ScifiHits.GetEntries()):
3131 d = event.Digi_ScifiHits[k]
3132 if not d.isValid(): continue
3133 hitDict[d.GetDetectorID()] = k
3134 hitList = list(hitDict.keys())
3135 if len(hitList)>0:
3136 hitList.sort()
3137 tmp = [ hitList[0] ]
3138 cprev = hitList[0]
3139 ncl = 0
3140 last = len(hitList)-1
3141 hitlist = ROOT.std.vector("sndScifiHit*")()
3142 for i in range(len(hitList)):
3143 if i==0 and len(hitList)>1: continue
3144 c=hitList[i]
3145 if (c-cprev)==1:
3146 tmp.append(c)
3147 if (c-cprev)!=1 or c==hitList[last]:
3148 first = tmp[0]
3149 N = len(tmp)
3150 hitlist.clear()
3151 for aHit in tmp: hitlist.push_back( event.Digi_ScifiHits[hitDict[aHit]])
3152 aCluster = ROOT.sndCluster(first,N,hitlist,Scifi)
3153 clusters.append(aCluster)
3154 if c!=hitList[last]:
3155 ncl+=1
3156 tmp = [c]
3157 cprev = c
3158 xHits = {}
3159 yHits = {}
3160 for s in range(5):
3161 xHits[s]=[]
3162 yHits[s]=[]
3163 for aCl in clusters:
3164 aCl.GetPosition(A,B)
3165 vertical = int(aCl.GetFirst()/100000)%10==1
3166 s = int(aCl.GetFirst()/1000000)-1
3167 if vertical:
3168 xHits[s].append(ROOT.TVector3(A))
3169 rc = h['clX'].Fill(aCl.GetN())
3170 else:
3171 yHits[s].append(ROOT.TVector3(A))
3172 rc = h['clY'].Fill(aCl.GetN())
3173 proj = {'X':xHits,'Y':yHits}
3174 for p in proj:
3175 sls = []
3176 for s1 in range(0,5):
3177 if len(proj[p][s1]) !=1: continue
3178 cl1 = proj[p][s1][0]
3179 for s2 in range(s1+1,5):
3180 if len(proj[p][s2]) !=1: continue
3181 cl2 = proj[p][s2][0]
3182 dz = abs(cl1[2]-cl2[2])
3183 if dz < 5: continue
3184 dzRep = 1./dz
3185 m = dzRep*(cl2-cl1)
3186 sls.append( m )
3187 for ix1 in range(0,len(sls)-1):
3188 for ix2 in range(ix1+1,len(sls)):
3189 if p=="X": rc = h['slopes'+p].Fill( sls[ix2][0]-sls[ix1][0])
3190 if p=="Y": rc = h['slopes'+p].Fill( sls[ix2][1]-sls[ix1][1])
3191 ut.bookCanvas(h,'slopes',' ',1024,768,1,2)
3192 h['slopes'].cd(1)
3193 h['slopesX'].GetXaxis().SetRangeUser(-0.2,0.2)
3194 h['slopesX'].SetTitle('x projection; delta slope [rad]')
3195 h['slopesX'].Draw()
3196 h['slopesX'].Fit('gaus','S','',-0.02,0.02)
3197 h['slopes'].Update()
3198 stats = h['slopesX'].FindObject('stats')
3199 stats.SetOptFit(111)
3200 h['slopes'].cd(2)
3201 h['slopesY'].GetXaxis().SetRangeUser(-0.2,0.2)
3202 h['slopesY'].SetTitle('y projection; delta slope [rad]')
3203 h['slopesY'].Draw()
3204 h['slopesY'].Fit('gaus','S','',-0.02,0.02)
3205 h['slopes'].Update()
3206 stats = h['slopesY'].FindObject('stats')
3207 stats.SetOptFit(111)
3208 for event in eventTree:
3209 if Nev<0: break
3210 Nev=Nev-1
3211 clusters = []
3212 hitDict = {}
3213 for k in range(event.Digi_ScifiHits.GetEntries()):
3214 d = event.Digi_ScifiHits[k]
3215 if not d.isValid(): continue
3216 hitDict[d.GetDetectorID()] = k
3217 hitList = list(hitDict.keys())
3218 if len(hitList)>0:
3219 hitList.sort()
3220 tmp = [ hitList[0] ]
3221 cprev = hitList[0]
3222 ncl = 0
3223 last = len(hitList)-1
3224 hitlist = ROOT.std.vector("sndScifiHit*")()
3225 for i in range(len(hitList)):
3226 if i==0 and len(hitList)>1: continue
3227 c=hitList[i]
3228 if (c-cprev)==1:
3229 tmp.append(c)
3230 if (c-cprev)!=1 or c==hitList[last]:
3231 first = tmp[0]
3232 N = len(tmp)
3233 hitlist.clear()
3234 for aHit in tmp: hitlist.push_back( event.Digi_ScifiHits[hitDict[aHit]])
3235 aCluster = ROOT.sndCluster(first,N,hitlist,scifiDet)
3236 clusters.append(aCluster)
3237 if c!=hitList[last]:
3238 ncl+=1
3239 tmp = [c]
3240 cprev = c
3241 xHits = {}
3242 yHits = {}
3243 for s in range(5):
3244 xHits[s]=[]
3245 yHits[s]=[]
3246 for aCl in clusters:
3247 aCl.GetPosition(A,B)
3248 vertical = int(aCl.GetFirst()/100000)%10==1
3249 s = int(aCl.GetFirst()/1000000)-1
3250 if vertical:
3251 xHits[s].append(ROOT.TVector3(A))
3252 rc = h['clX'].Fill(aCl.GetN())
3253 else:
3254 yHits[s].append(ROOT.TVector3(A))
3255 rc = h['clY'].Fill(aCl.GetN())
3256 proj = {'X':xHits,'Y':yHits}
3257 for p in proj:
3258 sls = []
3259 for s1 in range(0,5):
3260 if len(proj[p][s1]) !=1: continue
3261 cl1 = proj[p][s1][0]
3262 for s2 in range(s1+1,5):
3263 if len(proj[p][s2]) !=1: continue
3264

◆ Scifi_xPos()

Survey-MufiScifi.Scifi_xPos (   detID)

Definition at line 3110 of file Survey-MufiScifi.py.

3110def Scifi_xPos(detID):
3111 orientation = (detID//100000)%10
3112 nStation = 2*(detID//1000000-1)+orientation
3113 mat = (detID%100000)//10000
3114 X = detID%1000+(detID%10000)//1000*128
3115 return [nStation,mat,X] # even numbers are Y (horizontal plane), odd numbers X (vertical plane)
3116

◆ ScifiChannel_residuals()

Survey-MufiScifi.ScifiChannel_residuals ( )

Definition at line 3036 of file Survey-MufiScifi.py.

3036def ScifiChannel_residuals():
3037 for s in range(1, nScifi+1):
3038 for o in range(2):
3039 ut.bookHist(h,'res_Scifi'+str(s*10+o),'residual '+str(s*10+o)+';channel; [#mum]',
3040 512*3,-0.5,512*3-0.5,100,-2500.,2500.)
3041 fNtuple = ROOT.TFile("scifi_timing.root")
3042 tTimes = fNtuple.tTimes
3043 for nt in tTimes:
3044 for nHit in range(nt.nChannels):
3045 detID = nt.detIDs[nHit]
3046 X = Scifi_xPos(detID)
3047 s = X[0]//2+1
3048 o = X[0]%2
3049 n = X[1]*512+X[2]
3050 rc = h['res_Scifi'+str(s*10+o)].Fill(n,nt.D[nHit]*10000)
3051

◆ SciFiPropSpeed()

Survey-MufiScifi.SciFiPropSpeed ( )

Definition at line 3052 of file Survey-MufiScifi.py.

3052def SciFiPropSpeed():
3053 ut.bookCanvas(h,'v','',1600,1200,5,1)
3054 for s in range(1,6):
3055 tc = h['v'].cd(s)
3056 hname = 'dtScifivsdL_'+str(s)
3057 hist = h[hname]
3058 if hist.GetSumOfWeights()==0: continue
3059# find beam spot
3060 tmp = hist.ProjectionX()
3061 rc = tmp.Fit('gaus','S')
3062 res = rc.Get()
3063 fmin = res.Parameter(1) - 3*res.Parameter(2)
3064 fmax = res.Parameter(1) + 3*res.Parameter(2)
3065 hist.SetStats(0)
3066 xmin = max( hist.GetXaxis().GetBinCenter(1), fmin)
3067 xmax = min( hist.GetXaxis().GetBinCenter(hist.GetNbinsX()), fmax)
3068 hist.GetXaxis().SetRangeUser(xmin,xmax)
3069 hist.GetYaxis().SetRange(1,hist.GetNbinsY())
3070 rc = hist.ProjectionY().Fit('gaus','S')
3071 res = rc.Get()
3072 ymin = max( hist.GetYaxis().GetBinCenter(1), -5*res.Parameter(2))
3073 ymax = min( hist.GetYaxis().GetBinCenter(hist.GetNbinsY()), 5*res.Parameter(2))
3074 hist.GetYaxis().SetRangeUser(ymin,ymax)
3075 hist.Draw('colz')
3076# get time x correlation, X = m*dt + b
3077 h['g'+hname] = ROOT.TGraph()
3078 g = h['g'+hname]
3079 xproj = hist.ProjectionX('tmpx')
3080 if xproj.GetSumOfWeights()==0: continue
3081 np = 0
3082 Lmin = hist.GetSumOfWeights() * 0.005
3083 for nx in range(hist.GetXaxis().FindBin(xmin),hist.GetXaxis().FindBin(xmax)):
3084 tmp = hist.ProjectionY('tmp',nx,nx)
3085 if tmp.GetEntries()<Lmin:continue
3086 X = hist.GetXaxis().GetBinCenter(nx)
3087 rc = tmp.Fit('gaus','NQS')
3088 res = rc.Get()
3089 if not res: dt = tmp.GetMean()
3090 else: dt = res.Parameter(1)
3091 g.SetPoint(np,X,dt)
3092 np+=1
3093 g.SetLineColor(ROOT.kRed)
3094 g.SetLineWidth(2)
3095 g.Draw('same')
3096 rc = g.Fit('pol1','SQ','',xmin,xmax)
3097 g.GetFunction('pol1').SetLineColor(ROOT.kGreen)
3098 g.GetFunction('pol1').SetLineWidth(2)
3099 result = rc.Get()
3100 if not result: continue
3101 m = 1./result.Parameter(1)
3102 b = -result.Parameter(0) * m
3103 sign = '+'
3104 if b<0: sign = '-'
3105 txt = 'X = %5.1F #frac{cm}{ns} #times #frac{1}{2}#Deltat %s %5.1F '%(m,sign,abs(b))
3106 latex.SetTextSize(0.07)
3107 latex.DrawLatexNDC(0.2,0.8,txt)
3108
3109# Scifi specific code

◆ scifiTimings()

Survey-MufiScifi.scifiTimings ( )

Definition at line 3023 of file Survey-MufiScifi.py.

3023def scifiTimings():
3024# plot TDC vs dL
3025 fNtuple = ROOT.TFile("scifi_timing_"+str(options.runNumber)+".root")
3026 tTimes = fNtuple.tTimes
3027 ut.bookHist(h,'dTtwin','dt between neighbours;[ns]',100,-5,5)
3028 for nt in tTimes:
3029 for nHit in range(nt.nChannels):
3030 detID = nt.detIDs[nHit]
3031 dLL = (nt.dL[nHit] - nt.dL[tag]) / vsignal # light propagation in cm/ns
3032 dtdc = nt.tdc[nHit] - nt.tdc[tag]
3033
3034
3035# h['c'+str(new_list[3][0])].Draw()

◆ signalZoom()

Survey-MufiScifi.signalZoom (   smax)

Definition at line 3293 of file Survey-MufiScifi.py.

3293def signalZoom(smax):
3294 for mat in range(30):
3295 h['sig_'+str(mat)].GetXaxis().SetRangeUser(0.,smax)
3296 tc = h['signal'].cd(mat+1)
3297 tc.Update()
3298

◆ smallSiPMchannel()

Survey-MufiScifi.smallSiPMchannel (   i)

Definition at line 191 of file Survey-MufiScifi.py.

191def smallSiPMchannel(i):
192 if i==2 or i==5 or i==10 or i==13: return True
193 else: return False
194

◆ smallVsLargeSiPMs()

Survey-MufiScifi.smallVsLargeSiPMs (   Nev = -1)

Definition at line 741 of file Survey-MufiScifi.py.

741def smallVsLargeSiPMs(Nev=-1):
742 S = 2
743 for l in range(systemAndPlanes[S]):
744 ut.bookHist(h,'SVSl_'+str(l),'QDC large vs small sum',200,0.,200.,200,0.,200.)
745 ut.bookHist(h,'sVSl_'+str(l),'QDC large vs small average',200,0.,200.,200,0.,200.)
746 for side in ['L','R']:
747 for i1 in range(7):
748 for i2 in range(i1+1,8):
749 tag=''
750 if S==2 and smallSiPMchannel(i1): tag = 's'+str(i1)
751 else: tag = 'l'+str(i1)
752 if S==2 and smallSiPMchannel(i2): tag += 's'+str(i2)
753 else: tag += 'l'+str(i2)
754 ut.bookHist(h,'cor'+tag+'_'+side+str(l),'QDC channel i vs channel j',200,0.,200.,200,0.,200.)
755 for bar in range(systemAndBars[S]):
756 ut.bookHist(h,'cor'+tag+'_'+side+str(l)+str(bar),'QDC channel i vs channel j',200,0.,200.,200,0.,200.)
757
758 N=-1
759 if Nev < 0 : Nev = eventTree.GetEntries()
760 for event in eventTree:
761 N+=1
762 if N%options.heartBeat == 0: print('event ',N,' ',time.ctime())
763 if N>Nev: break
764 for aHit in event.Digi_MuFilterHits:
765 if not aHit.isValid(): continue
766 detID = aHit.GetDetectorID()
767 s = detID//10000
768 bar = (detID%1000)
769 if s!=2: continue
770 l = (detID%10000)//1000 # plane number
771 sumL,sumS,SumL,SumS = 0,0,0,0
772 allChannels = map2Dict(aHit,'GetAllSignals',mask=False)
773 nS = 0
774 nL = 0
775 for c in allChannels:
776 if s==2 and smallSiPMchannel(c) :
777 sumS+= allChannels[c]
778 nS += 1
779 else:
780 sumL+= allChannels[c]
781 nL+=1
782 if nL>0: SumL=sumL/nL
783 if nS>0: SumS=sumS/nS
784 rc = h['sVSl_'+str(l)].Fill(SumS,SumL)
785 rc = h['SVSl_'+str(l)].Fill(sumS/4.,sumL/12.)
786#
787 for side in ['L','R']:
788 offset = 0
789 if side=='R': offset = 8
790 for i1 in range(offset,offset+7):
791 if not i1 in allChannels: continue
792 qdc1 = allChannels[i1]
793 for i2 in range(i1+1,offset+8):
794 if not (i2) in allChannels: continue
795 if s==2 and smallSiPMchannel(i1): tag = 's'+str(i1-offset)
796 else: tag = 'l'+str(i1-offset)
797 if s==2 and smallSiPMchannel(i2): tag += 's'+str(i2-offset)
798 else: tag += 'l'+str(i2-offset)
799 qdc2 = allChannels[i2]
800 rc = h['cor'+tag+'_'+side+str(l)].Fill(qdc1,qdc2)
801 rc = h['cor'+tag+'_'+side+str(l)+str(bar)].Fill(qdc1,qdc2)
802 allChannels.clear()
803
804 ut.bookCanvas(h,'TSL','',1800,1400,3,2)
805 ut.bookCanvas(h,'STSL','',1800,1400,3,2)
806 for l in range(systemAndPlanes[S]):
807 tc = h['TSL'].cd(l+1)
808 tc.SetLogz(1)
809 aHist = h['sVSl_'+str(l)]
810 aHist.SetTitle(';small SiPM QCD av:large SiPM QCD av')
811 nmax = aHist.GetBinContent(aHist.GetMaximumBin())
812 aHist.SetMaximum( 0.1*nmax )
813 tc = h['sVSl_'+str(l)].Draw('colz')
814 myPrint(h['TSL'],"largeSiPMvsSmallSiPM")
815 for l in range(systemAndPlanes[S]):
816 tc = h['STSL'].cd(l+1)
817 tc.SetLogz(1)
818 aHist = h['SVSl_'+str(l)]
819 aHist.SetTitle(';small SiPM QCD sum/2:large SiPM QCD sum/6')
820 nmax = aHist.GetBinContent(aHist.GetMaximumBin())
821 aHist.SetMaximum( 0.1*nmax )
822 tc = h['SVSl_'+str(l)].Draw('colz')
823 myPrint(h['STSL'],"SumlargeSiPMvsSmallSiPM")
824
825 for l in range(systemAndPlanes[S]):
826 for side in ['L','R']:
827 ut.bookCanvas(h,'cor'+side+str(l),'',1800,1400,7,4)
828 k=1
829 for i1 in range(7):
830 for i2 in range(i1+1,8):
831 tag=''
832 if S==2 and smallSiPMchannel(i1): tag = 's'+str(i1)
833 else: tag = 'l'+str(i1)
834 if S==2 and smallSiPMchannel(i2): tag += 's'+str(i2)
835 else: tag += 'l'+str(i2)
836 tc = h['cor'+side+str(l)].cd(k)
837 for bar in range(systemAndBars[S]):
838 if bar == 0: h['cor'+tag+'_'+side+str(l)+str(bar)].Draw('colz')
839 else: h['cor'+tag+'_'+side+str(l)+str(bar)].Draw('colzsame')
840 k+=1
841 myPrint(h['cor'+side+str(l)],'QDCcor'+side+str(l))
842
843

◆ systemAndOrientation()

Survey-MufiScifi.systemAndOrientation (   s,
  plane 
)

Definition at line 153 of file Survey-MufiScifi.py.

153def systemAndOrientation(s,plane):
154 if s==1 and plane==2: return "vertical"
155 if s==3 and (plane%2==1 or plane == 6): return "vertical"
156 return "horizontal"
157

◆ testReversChannelMapping()

Survey-MufiScifi.testReversChannelMapping (   p = None)

Definition at line 4138 of file Survey-MufiScifi.py.

4138def testReversChannelMapping(p=None):
4139 import reverseMapping
4140 R = reverseMapping.reversChannelMapping()
4141 if not p: p = options.path.replace("convertedData","raw_data")+"/data/run_"+runNr
4142 R.Init(p)
4143 for event in eventTree:
4144 for aHit in eventTree.Digi_MuFilterHits:
4145 allChannels = map2Dict(aHit,'GetAllSignals')
4146 for c in allChannels:
4147 print(R.daqChannel(aHit,c))
4148

◆ TimeCalibrationNtuple()

Survey-MufiScifi.TimeCalibrationNtuple (   Nev = options.nEvents,
  nStart = 0 
)

Definition at line 2902 of file Survey-MufiScifi.py.

2902def TimeCalibrationNtuple(Nev=options.nEvents,nStart=0):
2903 if Nev < 0 : Nev = eventTree.GetEntries()
2904 maxD = 100
2905 fNtuple = ROOT.TFile("scifi_timing_"+str(options.runNumber)+".root",'recreate')
2906 tTimes = ROOT.TTree('tTimes','Scifi time calib')
2907 nChannels = array('i',1*[0])
2908 fTime = array('f',1*[0])
2909 detIDs = array('i',maxD*[0])
2910 tdc = array('f',maxD*[0.])
2911 qdc = array('f',maxD*[0.])
2912 dL = array('f',maxD*[0.])
2913 D = array('f',maxD*[0.])
2914 tTimes.Branch('fTime',fTime,'fTime/F')
2915 tTimes.Branch('nChannels',nChannels,'nChannels/I')
2916 tTimes.Branch('detIDs',detIDs,'detIDs[nChannels]/I')
2917 tTimes.Branch('tdc',tdc,'tdc[nChannels]/F')
2918 tTimes.Branch('qdc',tdc,'qdc[nChannels]/F')
2919 tTimes.Branch('D',D,'D[nChannels]/F')
2920 tTimes.Branch('dL',dL,'dL[nChannels]/F')
2921
2922 for n_ in range(nStart,nStart+Nev):
2923 rc = eventTree.GetEvent(n_)
2924 event = eventTree
2925 n_+=1
2926 if n_%100000==0: print("now at event ",n_)
2927 theTrack = Scifi_track()
2928 if not hasattr(theTrack,"getFittedState"): continue
2929 status = theTrack.getFitStatus()
2930 if status.isFitConverged():
2931 DetID2Key={}
2932 for nHit in range(event.Digi_ScifiHits.GetEntries()):
2933 DetID2Key[event.Digi_ScifiHits[nHit].GetDetectorID()] = nHit
2934 nHit = 0
2935 for nM in range(theTrack.getNumPointsWithMeasurement()):
2936 state = theTrack.getFittedState(nM)
2937 posM = state.getPos()
2938 M = theTrack.getPointWithMeasurement(nM)
2939 W = M.getRawMeasurement()
2940 clkey = W.getHitId()
2941 aCl = eventTree.Cluster_Scifi[clkey]
2942 fTime[0] = aCl.GetTime()
2943 for nh in range(aCl.GetN()):
2944 detID = aCl.GetFirst() + nh
2945 aHit = event.Digi_ScifiHits[ DetID2Key[detID] ]
2946 geo.modules['Scifi'].GetSiPMPosition(detID,A,B)
2947 if aHit.isVertical(): X = B-posM
2948 else: X = A-posM
2949 dL[nHit] = X.Mag()
2950 detIDs[nHit] = detID
2951 tdc[nHit] = aHit.GetTime()*TDC2ns
2952 qdc[nHit] = aHit.GetSignal()
2953 N = B-A
2954 X = posM-A
2955 V = X.Cross(N)
2956 D[nHit] = 0
2957 if abs(V.Z())>0: D[nHit] = V.Mag()/N.Mag()*V.Z()/abs(V.Z())
2958 nHit += 1
2959 if nHit==maxD:
2960 print('too many hits:',n_)
2961 break
2962 if nHit==maxD: break
2963 nChannels[0] = nHit
2964 tTimes.Fill()
2965 theTrack.Delete()
2966 fNtuple.Write()
2967 fNtuple.Close()
2968

◆ TimeStudy()

Survey-MufiScifi.TimeStudy (   Nev = options.nEvents,
  withDisplay = False 
)

Definition at line 974 of file Survey-MufiScifi.py.

974def TimeStudy(Nev=options.nEvents,withDisplay=False):
975 if Nev < 0 : Nev = eventTree.GetEntries()
976 ut.bookHist(h,'Vetotime','time',1000,0.,50.)
977 ut.bookHist(h,'UStime','time',1000,0.,50.)
978 ut.bookHist(h,'DStime','time',1000,0.,50.)
979 ut.bookHist(h,'Stime','time',1000,0.,50.)
980 ut.bookHist(h,'SvsDStime','; mean Scifi time [ns];mean Mufi time [ns]',100,0.,50.,100,0.,50.)
981 ut.bookHist(h,'VEvsUStime','; mean US time [ns];mean VE time [ns]',100,0.,50.,100,0.,50.)
982 ut.bookCanvas(h,'T','',900,1200,1,2)
983 tc = h['T'].cd(1)
984 h['Vetotime'].SetLineColor(ROOT.kOrange)
985 h['UStime'].SetLineColor(ROOT.kGreen)
986 h['DStime'].SetLineColor(ROOT.kRed)
987 N=-1
988 for event in eventTree:
989 N+=1
990 if N>Nev: break
991 h['UStime'].Reset()
992 h['DStime'].Reset()
993 h['Vetotime'].Reset()
994 h['Stime'].Reset()
995 for aHit in eventTree.Digi_MuFilterHits:
996 T = aHit.GetAllTimes()
997 s = aHit.GetDetectorID()//10000
998 for x in T:
999 t = x.second*TDC2ns
1000 if t>0:
1001 if s==1: rc = h['Vetotime'].Fill(t)
1002 if s==2: rc = h['UStime'].Fill(t)
1003 if s==3: rc = h['DStime'].Fill(t)
1004 stations = {}
1005 for aHit in eventTree.Digi_ScifiHits:
1006 t = aHit.GetTime()*TDC2ns
1007 rc = h['Stime'].Fill(t)
1008 stations[aHit.GetDetectorID()//1000000] = 1
1009 if len(stations)>3:
1010 rc = h['SvsDStime'].Fill(h['Stime'].GetMean(),h['DStime'].GetMean())
1011 rc = h['VEvsUStime'].Fill(h['UStime'].GetMean(),h['Vetotime'].GetMean())
1012 if withDisplay:
1013 tc = h['T'].cd(1)
1014 h['UStime'].Draw()
1015 h['DStime'].Draw('same')
1016 h['Vetotime'].Draw('same')
1017 tc = h['T'].cd(2)
1018 h['Stime'].Draw()
1019 rc = input("hit return for next event or q for quit: ")
1020 if rc=='q': break
1021 tc = h['T'].cd(1)
1022 h['SvsDStime'].Draw('colz')
1023 tc = h['T'].cd(2)
1024 h['SvsDStime_mufi'] = h['SvsDStime'].ProjectionY('SvsDStime_mufi')
1025 h['SvsDStime_scifi'] = h['SvsDStime'].ProjectionX('SvsDStime_scifi')
1026 h['Vetime'] = h['VEvsUStime'].ProjectionY('Vetime')
1027 h['UStime'] = h['VEvsUStime'].ProjectionX('UStime')
1028 h['SvsDStime_mufi'].SetLineColor(ROOT.kRed)
1029 h['SvsDStime_scifi'].SetLineColor(ROOT.kGreen)
1030 h['UStime'].SetLineColor(ROOT.kBlue)
1031 h['Vetime'].SetLineColor(ROOT.kOrange)
1032 h['UStime'].SetStats(0)
1033 h['Vetime'].SetStats(0)
1034 h['SvsDStime_mufi'].SetStats(0)
1035 h['SvsDStime_scifi'].SetStats(0)
1036 h['SvsDStime_mufi'].Draw()
1037 h['SvsDStime_scifi'].Draw('same')
1038 h['UStime'].Draw('same')
1039 h['Vetime'].Draw('same')
1040

◆ twoLangaufun()

Survey-MufiScifi.twoLangaufun (   x,
  par 
)

Definition at line 218 of file Survey-MufiScifi.py.

218def twoLangaufun(x,par):
219 N1 = langaufun(x,par)
220 par2 = [par[0],par[1]*2,par[4],par[3]]
221 N2 = langaufun(x,par2)
222 return N1+abs(N2)
223

◆ TwoTrackFinder()

Survey-MufiScifi.TwoTrackFinder (   nstart = 0,
  Nev = -1,
  sMin = 10,
  dClMin = 7,
  minDistance = 1.5,
  sepDistance = 0.5,
  debug = False 
)

Definition at line 3316 of file Survey-MufiScifi.py.

3316def TwoTrackFinder(nstart=0,Nev=-1,sMin=10,dClMin=7,minDistance=1.5,sepDistance=0.5,debug=False):
3317 if Nev < 0 :
3318 nstop = eventTree.GetEntries() - nstart
3319 for ecounter in range(nstart,nstop):
3320 event = eventTree
3321 rc = eventTree.GetEvent(ecounter)
3322 E = eventTree.EventHeader
3323 if ecounter%1000000==0: print('still here',ecounter)
3324 trackTask.clusScifi.Clear()
3325 trackTask.scifiCluster()
3326 clusters = trackTask.clusScifi
3327 sortedClusters={}
3328 for aCl in clusters:
3329 so = aCl.GetFirst()//100000
3330 if not so in sortedClusters: sortedClusters[so]=[]
3331 sortedClusters[so].append(aCl)
3332 if len(sortedClusters)<sMin: continue
3333 M=0
3334 for x in sortedClusters:
3335 if len(sortedClusters[x]) == 2: M+=1
3336 if M < dClMin: continue
3337 seeds = {}
3338 S = [-1,-1]
3339 for o in range(0,2):
3340# same procedure for both projections
3341# take seeds from from first station with 2 clusters
3342 for s in range(1,6):
3343 x = 10*s+o
3344 if x in sortedClusters:
3345 if len(sortedClusters[x])==2:
3346 sortedClusters[x][0].GetPosition(A,B)
3347 if o%2==1: pos0 = (A[0]+B[0])/2
3348 else: pos0 = (A[1]+B[1])/2
3349 sortedClusters[x][1].GetPosition(A,B)
3350 if o%2==1: pos1 = (A[0]+B[0])/2
3351 else: pos1 = (A[1]+B[1])/2
3352 if abs(pos0-pos1) > minDistance:
3353 S[o] = s
3354 break
3355 if S[o]<0: break # no seed found
3356 seeds[o]={}
3357 k = -1
3358 for c in sortedClusters[S[o]*10+o]:
3359 k += 1
3360 c.GetPosition(A,B)
3361 if o%2==1: pos = (A[0]+B[0])/2
3362 else: pos = (A[1]+B[1])/2
3363 seeds[o][k] = [[c,pos]]
3364 if k!=1: continue
3365 if abs(seeds[o][0][0][1] - seeds[o][1][0][1]) < sepDistance: continue
3366 for s in range(1,6):
3367 if s==S[o]: continue
3368 for c in sortedClusters[s*10+o]:
3369 c.GetPosition(A,B)
3370 if o%2==1: pos = (A[0]+B[0])/2
3371 else: pos = (A[1]+B[1])/2
3372 for k in range(2):
3373 if abs(seeds[o][k][0][1] - pos) < sepDistance:
3374 seeds[o][k].append([c,pos])
3375 if S[0]<0 or S[1]<0:
3376 passed = False
3377 else:
3378 passed = True
3379 if debug:
3380 print('-----',E.GetEventNumber(),'------')
3381 clusPos = {}
3382 for x in sortedClusters:
3383 clusPos[x] = []
3384 for c in sortedClusters[x]:
3385 c.GetPosition(A,B)
3386 if x%2==1: pos = (A[0]+B[0])/2
3387 else: pos = (A[1]+B[1])/2
3388 clusPos[x].append(pos)
3389 for x in sortedClusters: print(x,clusPos[x])
3390
3391 for o in range(0,2):
3392 for k in range(2):
3393 if len(seeds[o][k])<3:
3394 passed = False
3395 break
3396 if passed:
3397 tracks = []
3398 for k in range(2):
3399 # arbitrarly combine X and Y of combination 0
3400 n = 0
3401 hitlist = {}
3402 for o in range(0,2):
3403 for X in seeds[o][k]:
3404 hitlist[n] = X[0]
3405 n+=1
3406 theTrack = trackTask.fitTrack(hitlist)
3407 if not hasattr(theTrack,"getFittedState"):
3408 validTrack = False
3409 continue
3410 fitStatus = theTrack.getFitStatus()
3411 if not fitStatus.isFitConverged():
3412 theTrack.Delete()
3413 else:
3414 tracks.append(theTrack)
3415
3416 if len(tracks)==2:
3417 print('another 2track event',ecounter,E.GetEventNumber(),E.GetRunId(),E.GetFillNumber())
3418 for t in tracks: t.Delete()
3419

◆ USDSEnergy()

Survey-MufiScifi.USDSEnergy (   Nev = options.nEvents,
  nproc = 15 
)

Definition at line 1166 of file Survey-MufiScifi.py.

1166def USDSEnergy(Nev=options.nEvents,nproc=15):
1167#
1168 for b in ['','lb1','lip1','lb2']:
1169 ut.bookHist(h,'energy23'+b,'tot energy in DS vs US; US tot.QDC;DS tot.QDC ', 300,0.,15000,250,0.,20000.)
1170 ut.bookHist(h,'Tenergy23'+b,'tot energy in DS vs US with scifi track; US tot.QDC;DS tot.QDC',300,0.,15000,250,0.,20000.)
1171 ut.bookHist(h,'Denergy23'+b,'tot energy in DS vs US with DS track; US tot.QDC;DS tot.QDC',300,0.,15000,250,0.,20000.)
1172#
1173 fg = ROOT.TFile.Open(options.path+'FSdict.root')
1174 pkl = Unpickler(fg)
1175 FSdict = pkl.load('FSdict')
1176 fg.Close()
1177
1178 if options.runNumber in FSdict: fsdict = FSdict[options.runNumber]
1179 else: fsdict = False
1180
1181 if Nev < 0 : Nev = eventTree.GetEntries()
1182 N=0
1183 process = []
1184 print('number of events',Nev)
1185 for i in range(nproc):
1186 try:
1187 pid = os.fork()
1188 except OSError:
1189 print("Could not create a child process")
1190 print('pid',pid,i)
1191 if pid!=0:
1192 process.append(pid)
1193 # print('append process to list',process)
1194 else:
1195 dN = Nev//nproc
1196 nstart = i*dN
1197 nstop = nstart + dN
1198 if i==(nproc-1): nstop = Nev
1199 print('start ',i,nstart,nstop)
1200 for k in range(nstart,nstop):
1201 event = eventTree
1202 eventTree.GetEvent(k)
1203# check for b1,b2,IP1,IP2
1204 if fsdict:
1205 bunchNumber = eventTree.EventHeader.GetEventTime()%(4*3564)//4
1206 nb1 = (3564 + bunchNumber - fsdict['phaseShift1'])%3564
1207 nb2 = (3564 + bunchNumber - fsdict['phaseShift1']- fsdict['phaseShift2'])%3564
1208 b1 = nb1 in fsdict['B1']
1209 b2 = nb2 in fsdict['B2']
1210 IP1 = False
1211 IP2 = False
1212 if b1:
1213 IP1 = fsdict['B1'][nb1]['IP1']
1214 if b2:
1215 IP2 = fsdict['B2'][nb2]['IP2']
1216 # look for isolated bunch types
1217 lb1 = False
1218 if b1 and not b2: lb1 = True
1219 lip1 = False
1220 if lb1 and IP1: lip1 = True
1221 lb2 = False
1222 if b2 and not b1: lb2 = True
1223
1224 N+=1
1225 if N>Nev: break
1226 for aTrack in Reco_MuonTracks: aTrack.Delete()
1227 Reco_MuonTracks.Clear()
1228 rc = trackTask.ExecuteTask("DSScifi")
1229 Etot = {2:0,3:0}
1230 for aHit in eventTree.Digi_MuFilterHits:
1231 if not aHit.isValid(): continue
1232 detID = aHit.GetDetectorID()
1233 s = aHit.GetDetectorID()//10000
1234 if s<2: continue
1235 S = map2Dict(aHit,'SumOfSignals')
1236 Etot[s]+=S['Sum']
1237 rc = h['energy23'].Fill(Etot[2],Etot[3])
1238 if lb1: rc = h['energy23lb1'].Fill(Etot[2],Etot[3])
1239 if lb2: rc = h['energy23lb2'].Fill(Etot[2],Etot[3])
1240 if lip1: rc = h['energy23lip1'].Fill(Etot[2],Etot[3])
1241 trackTypes = [0,0]
1242 for t in Reco_MuonTracks:
1243 for aTrack in Reco_MuonTracks:
1244 if aTrack.GetUniqueID()==1: trackTypes[0]+=1
1245 else: trackTypes[1]+=1
1246 if trackTypes[0]==1:
1247 rc = h['Tenergy23'].Fill(Etot[2],Etot[3])
1248 if lb1: rc = h['Tenergy23lb1'].Fill(Etot[2],Etot[3])
1249 if lb2: rc = h['Tenergy23lb2'].Fill(Etot[2],Etot[3])
1250 if lip1: rc = h['Tenergy23lip1'].Fill(Etot[2],Etot[3])
1251 if trackTypes[1]==1:
1252 rc = h['Denergy23'].Fill(Etot[2],Etot[3])
1253 if lb1: rc = h['Denergy23lb1'].Fill(Etot[2],Etot[3])
1254 if lb2: rc = h['Denergy23lb2'].Fill(Etot[2],Etot[3])
1255 if lip1: rc = h['Denergy23lip1'].Fill(Etot[2],Etot[3])
1256 if Etot[2]>5000 and Etot[3]>5000:
1257 print("HE",nstart,k,eventTree.GetChainOffset(),eventTree.GetCurrentFile().GetName(),Etot[2],Etot[3],lb1,lb2,lip1)
1258 if Etot[2]<800 and Etot[3]>2000 and trackTypes[0]==1:
1259 print("DIS",nstart,k,eventTree.GetChainOffset(),eventTree.GetCurrentFile().GetName(),Etot[2],Etot[3],lb1,lb2,lip1)
1260
1261 ut.writeHists(h,'tmp_'+str(i))
1262 exit(0)
1263#
1264 while process:
1265 pid,exit_code = os.wait()
1266 if pid == 0: time.sleep(100)
1267 else:
1268 print('child process has finished',pid,exit_code)
1269 process.remove(pid)
1270 for i in range(nproc):
1271 ut.readHists(h,'tmp_'+str(i))
1272 print('i am finished')
1273 ut.writeHists(h,"USDSEnergy-"+str(options.runNumber)+".root")
1274# plot
1275 for b in ['','lb1','lip1','lb2']:
1276 ut.bookCanvas(h,'tc1'+b,'USDS Energy',2400,1800,3,2)
1277 tc = h['tc1'+b].cd(1)
1278 tc.SetLogz()
1279 h['energy23'+b].Draw('colz')
1280 tc = h['tc1'+b].cd(2)
1281 tc.SetLogy()
1282 h['energy23'+b].ProjectionX().Draw()
1283 tc = h['tc1'+b].cd(3)
1284 tc.SetLogy()
1285 h['energy23'+b].ProjectionY().Draw()
1286 tc = h['tc1'+b].cd(4)
1287 tc.SetLogz()
1288 h['Tenergy23'+b].Draw('colz')
1289 tc = h['tc1'+b].cd(5)
1290 tc.SetLogy()
1291 h['Tenergy23'+b].ProjectionX().Draw()
1292 tc = h['tc1'+b].cd(6)
1293 tc.SetLogy()
1294 h['Tenergy23'+b].ProjectionY().Draw()
1295
1296 return
1297

◆ USshower()

Survey-MufiScifi.USshower (   Nev = options.nEvents)

Definition at line 1083 of file Survey-MufiScifi.py.

1083def USshower(Nev=options.nEvents):
1084 zUS0 = zPos['MuFilter'][20] - 10
1085 zUS4 = zPos['MuFilter'][24] + 10
1086 for x in ['','-small']:
1087 ut.bookHist(h,'shower'+x,'energy vs z',200,0.,10000.,20,zUS0,zUS4)
1088 ut.bookHist(h,'showerX'+x,'energy vs z',200,0.,10000.,20,zUS0,zUS4)
1089 ut.bookHist(h,'wshower'+x,'z weighted energy ',100,-300.,0.)
1090 ut.bookHist(h,'zyshower'+x,'y vs z weighted energy ',20,zUS0,zUS4,11,-0.5,10.5)
1091 for p in range(systemAndPlanes[2]):
1092 ut.bookHist(h,'SvsL'+str(p),'small vs large Sipms plane' +str(p)+';large [QCD]; small [QCD] ',100,0.1,250.,100,0.1,100.)
1093 if Nev < 0 : Nev = eventTree.GetEntries()
1094 N=0
1095 for event in eventTree:
1096 N+=1
1097 if N>Nev: break
1098 UShits = {}
1099 UShitsBar = {}
1100 for aHit in eventTree.Digi_MuFilterHits:
1101 if not aHit.isValid(): continue
1102 detID = aHit.GetDetectorID()
1103 s = aHit.GetDetectorID()//10000
1104 if s!=2: continue
1105 p = (aHit.GetDetectorID()//1000)%10
1106 S = map2Dict(aHit,'SumOfSignals')
1107 rc = h['SvsL'+str(p)].Fill(S['SumL'],S['SumS'])
1108 plane = (aHit.GetDetectorID()//1000)%10
1109 bar = aHit.GetDetectorID()%100
1110 if not plane in UShits:
1111 UShits[plane]=0
1112 UShitsBar[plane]={}
1113 UShits[100+plane]=0
1114 UShitsBar[100+plane]={}
1115 if not bar in UShitsBar[plane]:
1116 UShitsBar[plane][bar]=0
1117 UShitsBar[100+plane][bar]=0
1118 UShits[plane]+=S['Sum']
1119 UShitsBar[plane][bar]+=S['Sum']
1120 UShits[100+plane]+=S['SumS']
1121 UShitsBar[100+plane][bar]+=S['SumS']
1122 s = 2
1123 for plane in UShits:
1124 z = zPos['MuFilter'][s*10+plane%100]
1125 x = ''
1126 if plane > 99: x='-small'
1127 rc = h ['shower'+x].Fill(UShits[plane],z)
1128 if 0 in UShits:
1129 if UShits[0]>750: rc = h['showerX'+x].Fill(UShits[plane],z)
1130 rc = h ['wshower'+x].Fill(z,UShits[plane])
1131 for bar in UShitsBar[plane]:
1132 rc = h ['zyshower'+x].Fill(z,bar,UShitsBar[plane][bar])
1133 ut.bookCanvas(h,'lego','',900,1600,1,2)
1134 energy = {46:180,49:180,56:140,58:140,72:240,73:240,74:240,89:300,90:300,91:300}
1135 gain = {46:2.5,49:3.65,52:3.65,54:2.5,56:2.5,58:3.65,72:1.0,73:2.5,74:3.65,86:3.65,87:2.5,88:1.0,89:3.65,90:2.5,91:1.0}
1136 tc = h['lego'].cd(1)
1137 tc.SetPhi(-20.5)
1138 tc.SetTheta(21.1)
1139
1140 text = ""
1141 if options.runNumber in energy: text +="E=%5.1F"%(energy[options.runNumber])
1142 if options.runNumber in gain: text +=" with gain=%5.2F"%(gain[options.runNumber])
1143 h ['zyshower'].SetTitle(text+';z [cm]; y [bar number];QDC')
1144 h ['zyshower'].SetStats(0)
1145 h ['zyshower'].Draw('lego2')
1146 tc = h['lego'].cd(2)
1147 tc.SetPhi(-20.5)
1148 tc.SetTheta(21.1)
1149 h ['zyshower-small'].SetTitle('small sipms;z [cm]; y [bar number];QDC')
1150 h ['zyshower-small'].SetStats(0)
1151 h ['zyshower-small'].Draw('lego2')
1152 myPrint(h['lego'],'shower',withRootFile=True)
1153
1154 ut.bookCanvas(h,'CorLS','',900,1200,1,5)
1155 h['SvsL'] = h['SvsL0'].Clone('SvsL')
1156 h['SvsL'].SetTitle('small vs large Sipms all planes')
1157 for p in range(1,systemAndPlanes[2]):
1158 h['SvsL'].Add(h['SvsL'+str(p)])
1159 h['SvsL'].SetStats(0)
1160 for p in range(systemAndPlanes[2]):
1161 tc = h['CorLS'].cd(p+1)
1162 h['SvsL'+str(p)].SetStats(0)
1163 h['SvsL'+str(p)].Draw('colz')
1164 myPrint(h['CorLS'],'QCDsmallCSQCDlarge')
1165

◆ weightedY()

Survey-MufiScifi.weightedY (   aHist)

Definition at line 2579 of file Survey-MufiScifi.py.

2579def weightedY(aHist):
2580 xax = aHist.GetXaxis()
2581 nx = xax.GetNbins()
2582 W = 1E-20
2583 mean = 0
2584 for j in range(1,nx+1):
2585 N = aHist.GetBinContent(j)
2586 if not N>0: continue
2587 E = aHist.GetBinError(j)
2588 if E>0:
2589 mean+=N/E**2
2590 W+= 1/E**2
2591 return [mean/W,ROOT.TMath.Sqrt(1./W)]
2592
2593

Variable Documentation

◆ A

Survey-MufiScifi.A

Definition at line 141 of file Survey-MufiScifi.py.

◆ B

Survey-MufiScifi.B

Definition at line 141 of file Survey-MufiScifi.py.

◆ barColor

list Survey-MufiScifi.barColor
Initial value:
1= [ROOT.kRed,ROOT.kRed-7,ROOT.kMagenta,ROOT.kMagenta-6,ROOT.kBlue,ROOT.kBlue-9,
2 ROOT.kCyan,ROOT.kAzure-4,ROOT.kGreen,ROOT.kGreen+3]

Definition at line 1853 of file Survey-MufiScifi.py.

◆ Cluster_Scifi

Survey-MufiScifi.Cluster_Scifi = trackTask.clusScifi

Definition at line 368 of file Survey-MufiScifi.py.

◆ command

str Survey-MufiScifi.command = tmp[0]+"("

Definition at line 4156 of file Survey-MufiScifi.py.

◆ data

str Survey-MufiScifi.data = "sndsw_raw_"+runNr+".root"

Definition at line 109 of file Survey-MufiScifi.py.

◆ default

Survey-MufiScifi.default

Definition at line 91 of file Survey-MufiScifi.py.

◆ dest

Survey-MufiScifi.dest

Definition at line 90 of file Survey-MufiScifi.py.

◆ dirlist

Survey-MufiScifi.dirlist = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls "+options.path,shell=True) )

Definition at line 107 of file Survey-MufiScifi.py.

◆ eventChain

Survey-MufiScifi.eventChain = ROOT.TChain('rawConv')

Definition at line 319 of file Survey-MufiScifi.py.

◆ eventTree

Survey-MufiScifi.eventTree

Definition at line 363 of file Survey-MufiScifi.py.

◆ f

Survey-MufiScifi.f = ROOT.TFile.Open(options.fname)

Definition at line 328 of file Survey-MufiScifi.py.

◆ False

Survey-MufiScifi.False

Definition at line 91 of file Survey-MufiScifi.py.

◆ freq

Survey-MufiScifi.freq = u.snd_freq

Definition at line 163 of file Survey-MufiScifi.py.

◆ geo

Survey-MufiScifi.geo

Definition at line 130 of file Survey-MufiScifi.py.

◆ globA

Survey-MufiScifi.globA

Definition at line 141 of file Survey-MufiScifi.py.

◆ globB

Survey-MufiScifi.globB

Definition at line 141 of file Survey-MufiScifi.py.

◆ h

dict Survey-MufiScifi.h = {}

Definition at line 87 of file Survey-MufiScifi.py.

◆ help

Survey-MufiScifi.help

Definition at line 90 of file Survey-MufiScifi.py.

◆ houghTransform

bool Survey-MufiScifi.houghTransform = False

Definition at line 346 of file Survey-MufiScifi.py.

◆ Ikey

Survey-MufiScifi.Ikey = ROOT.std.vector('int')()

Definition at line 64 of file Survey-MufiScifi.py.

◆ int

Survey-MufiScifi.int

Definition at line 90 of file Survey-MufiScifi.py.

◆ ioman

Survey-MufiScifi.ioman = ROOT.FairRootManager.Instance()

Definition at line 334 of file Survey-MufiScifi.py.

◆ latex

Survey-MufiScifi.latex = ROOT.TLatex()

Definition at line 142 of file Survey-MufiScifi.py.

◆ locA

Survey-MufiScifi.locA

Definition at line 141 of file Survey-MufiScifi.py.

◆ locB

Survey-MufiScifi.locB

Definition at line 141 of file Survey-MufiScifi.py.

◆ MuFilter

Survey-MufiScifi.MuFilter = geo.modules['MuFilter']

Definition at line 132 of file Survey-MufiScifi.py.

◆ muon_reco_task

Survey-MufiScifi.muon_reco_task = SndlhcMuonReco.MuonReco()

Definition at line 348 of file Survey-MufiScifi.py.

◆ nav

Survey-MufiScifi.nav = ROOT.gGeoManager.GetCurrentNavigator()

Definition at line 134 of file Survey-MufiScifi.py.

◆ nMats

Survey-MufiScifi.nMats = Scifi.GetConfParI("Scifi/nmats")

Definition at line 137 of file Survey-MufiScifi.py.

◆ None

Survey-MufiScifi.None

Definition at line 92 of file Survey-MufiScifi.py.

◆ nScifi

Survey-MufiScifi.nScifi = Scifi.GetConfParI("Scifi/nscifi")

Definition at line 136 of file Survey-MufiScifi.py.

◆ nVeto

Survey-MufiScifi.nVeto = MuFilter.GetConfParI("MuFilter/NVetoPlanes")

Definition at line 139 of file Survey-MufiScifi.py.

◆ options

Survey-MufiScifi.options = parser.parse_args()

Definition at line 99 of file Survey-MufiScifi.py.

◆ OT

Survey-MufiScifi.OT = ioman.GetSink().GetOutTree()

Definition at line 366 of file Survey-MufiScifi.py.

◆ outFile

Survey-MufiScifi.outFile = ROOT.TMemFile('dummy','CREATE')

Definition at line 336 of file Survey-MufiScifi.py.

◆ parser

Survey-MufiScifi.parser = ArgumentParser()

Definition at line 89 of file Survey-MufiScifi.py.

◆ partitions

int Survey-MufiScifi.partitions = 0

Definition at line 103 of file Survey-MufiScifi.py.

◆ path

Survey-MufiScifi.path = options.path

Definition at line 102 of file Survey-MufiScifi.py.

◆ rc

Survey-MufiScifi.rc = eventChain.GetEvent(0)

Definition at line 332 of file Survey-MufiScifi.py.

◆ Reco_MuonTracks

Survey-MufiScifi.Reco_MuonTracks = trackTask.fittedTracks

Definition at line 367 of file Survey-MufiScifi.py.

◆ required

Survey-MufiScifi.required

Definition at line 90 of file Survey-MufiScifi.py.

◆ run

Survey-MufiScifi.run = ROOT.FairRunAna()

Definition at line 333 of file Survey-MufiScifi.py.

◆ runNr

Survey-MufiScifi.runNr = str(options.runNumber).zfill(6)

Definition at line 101 of file Survey-MufiScifi.py.

◆ Scifi

Survey-MufiScifi.Scifi = geo.modules['Scifi']

Definition at line 133 of file Survey-MufiScifi.py.

◆ sdict

dict Survey-MufiScifi.sdict = {1:'Veto',2:'US',3:'DS'}

Definition at line 161 of file Survey-MufiScifi.py.

◆ sink

Survey-MufiScifi.sink = ROOT.FairRootFileSink(outFile)

Definition at line 343 of file Survey-MufiScifi.py.

◆ source

Survey-MufiScifi.source = ROOT.FairFileSource(eventChain.GetCurrentFile())

Definition at line 337 of file Survey-MufiScifi.py.

◆ str

Survey-MufiScifi.str

Definition at line 92 of file Survey-MufiScifi.py.

◆ systemAndBars

dict Survey-MufiScifi.systemAndBars
Initial value:
1= {1:MuFilter.GetConfParI("MuFilter/NVetoBars"),
2 2:MuFilter.GetConfParI("MuFilter/NUpstreamBars"),
3 3:MuFilter.GetConfParI("MuFilter/NDownstreamBars")}
MuFilter::GetConfParI
Int_t GetConfParI(TString name)
Definition MuFilter.h:47

Definition at line 149 of file Survey-MufiScifi.py.

◆ systemAndChannels

dict Survey-MufiScifi.systemAndChannels
Initial value:
1= {1:[MuFilter.GetConfParI("MuFilter/VetonSiPMs"),0],
2 2:[MuFilter.GetConfParI("MuFilter/UpstreamnSiPMs")-2,2],
3 3:[MuFilter.GetConfParI("MuFilter/DownstreamnSiPMs"),0]}

Definition at line 158 of file Survey-MufiScifi.py.

◆ systemAndPlanes

dict Survey-MufiScifi.systemAndPlanes
Initial value:
1= {1:MuFilter.GetConfParI("MuFilter/NVetoPlanes"),
2 2:MuFilter.GetConfParI("MuFilter/NUpstreamPlanes"),
3 3:2*MuFilter.GetConfParI("MuFilter/NDownstreamPlanes")-1}

Definition at line 146 of file Survey-MufiScifi.py.

◆ TDC2ns

Survey-MufiScifi.TDC2ns = u.snd_TDC2ns

Definition at line 164 of file Survey-MufiScifi.py.

◆ Tkey

Survey-MufiScifi.Tkey = ROOT.std.vector('TString')()

Definition at line 63 of file Survey-MufiScifi.py.

◆ tmp

Survey-MufiScifi.tmp = options.command.split(';')

Definition at line 4150 of file Survey-MufiScifi.py.

◆ trackTask

Survey-MufiScifi.trackTask = SndlhcTracking.Tracking()

Definition at line 353 of file Survey-MufiScifi.py.

◆ type

Survey-MufiScifi.type

Definition at line 90 of file Survey-MufiScifi.py.

◆ Value

Survey-MufiScifi.Value = ROOT.std.vector('float')()

Definition at line 65 of file Survey-MufiScifi.py.

◆ xrdb

Survey-MufiScifi.xrdb = ROOT.FairRuntimeDb.instance()

Definition at line 358 of file Survey-MufiScifi.py.

◆ zPos

Survey-MufiScifi.zPos = getAverageZpositions()

Definition at line 316 of file Survey-MufiScifi.py.