SND@LHC Software
Loading...
Searching...
No Matches
compactingBackgroundProduction.py
Go to the documentation of this file.
1from __future__ import print_function
2import os,ROOT,sys,subprocess,pickle,time,datetime
3import rootUtils as ut
4
5pdg = ROOT.TDatabasePDG()
6charm = False
7
8# functions, PoT: statistics per file
9# TotStat: total statistics in given directory
10# GetGoodAndBadRuns: list of good and bad runs within given dates
11# addRuns: concatenate good runs
12# YandexProd(): specify start end date, select good runs, write to pickle file, concatenate runs requires start number !
13
14# prod = ''
15# globalPath="/eos/experiment/ship/skygrid/background-prod-2018/"
16# fnames = 'pythia8_Geant4_1_10.0.root'
17# ecut = '10.0'
18
19# final statistics: 65041000000 pot -> weight 768.75
20
21prod = '1GeV'
22globalPath="/eos/experiment/ship/skygrid/background-prod-2018-1gev/"
23fnames = 'pythia8_Geant4_1_1.0.root'
24ecut = '1.0'
25
26def GetGoodAndBadRuns(startDate,endDate):
27# find bad runs, try to recover
28 goodRuns=[]
29 badRuns=[]
30 fileName = "pythia8_Geant4_1_"+ecut+".root"
31 os.system("ls -l --full-time "+globalPath+" >inventory.lst")
32 f=open("inventory.lst")
33 runs = f.readlines()
34 f.close()
35 N=0
36 for r in runs:
37 tmp = r.split(' ')
38 if len(tmp) != 9:
39 print("wrong format",tmp)
40 continue
41 date = tmp[5].split('-')
42 time = tmp[6].split(':')
43 fileDate = datetime.datetime(int(date[0]),int(date[1]),int(date[2]),int(time[0]),int(time[1]),int(time[2].split('.')[0]))
44 if fileDate < startDate or fileDate>endDate: continue
45 theRun = int(tmp[8].replace('\n',''))
46 test = os.listdir(globalPath+str(theRun))
47 if len(test)<2:continue
48 N+=1
49 rc = os.system('grep -q "Number of events produced with activity after hadron absorber" '+globalPath+str(theRun)+"/stdout")
50 if rc!=0:
51 badRuns.append(theRun)
52 continue
53 for x in test:
54 if x.find('run_fixedTarget')>0:
55 test2 = os.listdir(globalPath+str(theRun)+'/'+x)
56 bad = True
57 tmpF = False
58 f = globalPath+str(theRun)+"/"+x+"/"+fileName
59 if fnames+"tmp" in test2:
60 tmpF = True
61 f = f+"tmp"
62 elif fnames in test2:
63 f = f
64 else:
65 badRuns.append(theRun)
66 continue
67 try:
68 t = ROOT.TFile.Open(os.environ["EOSSHIP"]+f)
69 if not t:
70 badRuns.append(theRun)
71 continue
72 if N%1000 == 0: print(N,t.GetName())
73 if t.ReadKeys()==0:
74 t.Close()
75 badRuns.append(theRun)
76 continue
77 if t.FindObjectAny('cbmsim'):
78 bad = False
79 goodRuns.append(f)
80 t.Close()
81 except:
82 badRuns.append(theRun)
83 continue
84 return goodRuns,badRuns
85
86def addRuns(goodRuns,Nstart=0):
87 N=0
88 while(1>0):
89 cmd = ""
90 for i in range(N,min(N+1000,len(goodRuns))):
91 cmd += " $EOSSHIP"+goodRuns[i]
92 tmpFile = "pythia8_Geant4_"+ecut+"_c"+str(N+Nstart)+".root"
93 rc = os.system("hadd -j 10 -O "+tmpFile + " " +cmd)
94 if rc != 0:
95 print("hadd failed, stop",N)
96 return
97 rc = os.system("xrdcp "+tmpFile+" $EOSSHIP/eos/experiment/ship/data/Mbias/background-prod-2018/"+tmpFile)
98 if rc != 0:
99 print("copy to EOS failed, stop",N,N+Nstart)
100 else:
101 rc = os.system("rm "+tmpFile)
102 N+=1000
103 if N > len(goodRuns): break
104def YandexProd():
105 #startDate = datetime.datetime(2018, 1, 1, 0, 0) # N=0
106 #endDate = datetime.datetime(2018, 2, 1, 0, 0)
107 #startDate = datetime.datetime(2018, 2, 1, 0, 0) # N=23000
108 #endDate = datetime.datetime(2018, 2, 6, 0, 0)
109 #startDate = datetime.datetime(2018, 2, 6, 0, 0) # N=33000
110 #endDate = datetime.datetime(2018, 2, 12, 0, 0)
111 #startDate = datetime.datetime(2018, 2, 12, 0, 0) # N=45000
112 #endDate = datetime.datetime(2018, 2, 15, 0, 0)
113 #startDate = datetime.datetime(2018, 2, 15, 0, 0) # N=51000
114 #endDate = datetime.datetime(2018, 2, 21, 0, 0)
115 #startDate = datetime.datetime(2018, 2, 21, 0, 0) # N=63000
116 #endDate = datetime.datetime(2018, 2, 2, 0, 0)
117 #startDate = datetime.datetime(2018, 1, 1, 0, 0) # start with 1 GeV production
118 #endDate = datetime.datetime(2018, 2, 27, 0, 0) # N=0
119 #startDate = datetime.datetime(2018, 2, 27, 0, 0) #
120 #endDate = datetime.datetime(2018, 3, 6, 0, 0) # N=3000
121 #startDate = datetime.datetime(2018, 3, 6, 0, 0) #
122 #endDate = datetime.datetime(2018, 3, 12, 0, 0) # N=9000
123 goodRuns,badRuns = GetGoodAndBadRuns(startDate,endDate)
124 pName = 'goodAndBadRuns'+prod+'_'+startDate.__str__().split(' ')[0]+'_'+endDate.__str__().split(' ')[0]+'.pkl'
125 fpi=open(pName,'w')
126 database = {}
127 database['goodruns']=goodRuns
128 database['badRuns']=badRuns
129 pickle.dump(database,fpi)
130 fpi.close()
131 fpi=open(pName)
132 database = pickle.load(fpi)
133 addRuns(database['goodruns'],20000) # next cycle
134
135def addAllHistograms():
136 h={}
137 ecut = '10.0'
138 Nmax=45000
139 path = os.environ['EOSSHIP']+"/eos/experiment/ship/data/Mbias/background-prod-2018/"
140 ut.bookCanvas(h,key='canvas',title='debug',nx=1600,ny=1200,cx=1,cy=1)
141 h['canvas'].SetLogy(1)
142 for i in range(0,Nmax,1000):
143 fname = "pythia8_Geant4_"+ecut+"_c"+str(i)+".root"
144 ut.readHists(h,path+fname)
145 if i==0:
146 h[1012].Draw()
147 for x in h.keys():
148 if h[x].GetName().find('proj')>0: rc = h.pop(x)
149 ut.writeHists(h,"pythia8_Geant4_"+ecut+"_c"+str(Nmax)+"-histos.root")
150
151def compactifyCascade(cycle):
152 ncpus = 20
153 path = "/afs/cern.ch/project/lbcern/vol1/truf/charm/"
154 cmd = ''
155 Ntot = 0
156 NperJob = 2000000
157 for i in range(cycle,cycle+ncpus):
158 fName = path+"run"+str(i)+"/Cascade-run"+str(i)+"-parp16-MSTP82-1-MSEL4.root"
159 f=open(path+"run"+str(i)+"/log"+str(i))
160 success = False
161 for l in f.readlines():
162 if not l.find('Macro finished succesfully')<0: success = True
163 if not success:
164 print("job not finished properly",fName)
165 continue
166 cmd += fName +" "
167 f.close()
168 Ntot+= NperJob
169 if cmd.find('root')<0:
170 print('no file found, exit')
171 else:
172 stat = str( int(Ntot/1E6))+'Mpot'
173 outFile = "Cascade-run"+str(cycle)+"-"+str(cycle+ncpus-1)+"-parp16-MSTP82-1-MSEL4-"+stat+".root"
174 rc = os.system("hadd -O "+outFile + " " +cmd)
175 rc = os.system("xrdcp "+outFile+" $EOSSHIP/eos/experiment/ship/data/Mbias/background-prod-2018/"+outFile)
176 if rc != 0:
177 print("copy to EOS failed, stop",outFile)
178 else:
179 rc = os.system("rm "+outFile)
180
181
182
183# some old stuff
184def compactify(charm):
185 globalPath="/afs/cern.ch/project/lbcern/vol2/truf/muonBackground"
186 ecut = '10.0'
187 if charm:
188 allDirs = os.listdir(globalPath+"/charm")
189 allFiles = []
190 for r in range(int(runMin),int(runMax)+1):
191# collect the 20 subdirectories connected to a run
192 nr = '9000'+str(r)
193 subruns = []
194 badFiles = []
195 for x in allDirs:
196 for k in range(20):
197 sr = "{0:0>2}".format(k)
198 if not x.find(nr+sr)<0:
199 if not x.find('log')<0: continue
200 fn = 'pythia8_Geant4_'+nr+sr+'_'+ecut+'.root'
201 aFile = globalPath+"/charm/"+x+'/'+fn
202 if os.path.exists(aFile):
203 try:
204 t = ROOT.TFile.Open(aFile)
205 if not t:
206 badFiles.append(d)
207 continue
208 if t.ReadKeys()==0:
209 badFiles.append(d)
210 continue
211 if t.FindObjectAny('cbmsim'):
212 subruns.append(aFile)
213 except:
214 badFiles.append(aFile)
215 continue
216 ldir = ' '
217 for x in subruns: ldir+= (x+" ")
218 outputr = 'pythia8_Geant4_charm_'+nr+'_'+ecut+'.root'
219 allFiles += (outputr+" ")
220 os.system('hadd ' + outputr +' '+ldir)
221 output = 'pythia8_Geant4_charm_'+str(runMin)+'_'+str(runMax)+'_'+ecut+'.root'
222 os.system('hadd ' + output +' '+allFiles)
223 makeHistos(output)
224 else:
225 output = 'pythia8_Geant4_'+str(runMin)+'-'+str(runMax)+'_'+ecut+'.root'
226 if not checkOnly:
227 ldir = ''
228 badFiles = []
229 for d in os.listdir(globalPath):
230 if d.find( 'run_fixedTarget') < 0 : continue
231 srun = d.split('run_fixedTarget_')[1]
232 run = int(srun)
233 if not run > runMax and not run < runMin:
234 f = globalPath+'/'+d+'/pythia8_Geant4_'+srun+'_'+ecut+'.root '
235 ftmp = f.replace('.root ','.roottmp')
236 if os.path.exists(ftmp):f=ftmp+' '
237 try:
238 t = ROOT.TFile.Open(f)
239 if not t:
240 badFiles.append(d)
241 continue
242 if t.ReadKeys()==0:
243 badFiles.append(d)
244 continue
245 if t.FindObjectAny('cbmsim'):
246 ldir+=f
247 except:
248 badFiles.append(d)
249 continue
250 os.system('hadd '+output+' '+ldir)
251 makeHistos(output)
252
253def makeHistos(rfile):
254 f=ROOT.TFile.Open(rfile)
255 sTree = f.cbmsim
256 nTot = 0
257 for k in f.GetListOfKeys():
258 if k.GetName() == 'FileHeader':
259 tmp = k.GetTitle().split('=')[1]
260 tmp2 = tmp.split('with')[0]
261 if tmp2.find('E')<0: nTot += int(tmp2)
262 else: nTot += float(tmp2)
263 print("POT = ",nTot," number of events:",sTree.GetEntries())
264# particle statistics
265 h={}
266 ut.bookHist(h,'pids','pid',19999,-9999.5,9999.5)
267 ut.bookHist(h,'test','muon p/pt',100,0.,400.,100,0.,5.)
268 diMuonDecays = [221, 223, 113, 331, 333]
269 pDict = {}
270 procDict = {}
271 for n in range(sTree.GetEntries()):
272 rc = sTree.GetEvent(n)
273 for p in sTree.vetoPoint:
274 t = sTree.MCTrack[p.GetTrackID()]
275 pid = t.GetPdgCode()
276 rc = h['pids'].Fill(pid)
277 if abs(pid)==13:
278 procID = t.GetProcName().Data()
279 mother = t.GetMotherId()
280 if not mother < 0:
281 moPid = sTree.MCTrack[mother].GetPdgCode()
282 name = pdg.GetParticle(moPid).GetName()
283 name = procID+' '+name
284 if name not in h: h[name]=h['test'].Clone(name)
285 rc=h[name].Fill(t.GetP(),t.GetPt())
286 if procID not in h: h[procID]=h['test'].Clone(procID)
287 rc=h[procID].Fill(t.GetP(),t.GetPt())
288 for x in h:
289 h[x].Scale(1./nTot)
290 tmp = rfile.split('/')
291 hname = tmp[len(tmp)-1].replace('pythia8_Geant4','Histos')
292 ut.writeHists(h,hname)
293
294hunbiased = {}
295hbiased = {}
296import operator
297def makePrintout():
298 # Histos_2000000-2001500_10.0.root
299 ut.readHists(hunbiased,'/media/microdisk/HNL/muonBackground/Histos_1000000-1000600_10.0.root')
300 ut.readHists(hbiased,'hadded_Histos_1_10.0.root')
301
302 unbiased = {}
303 biased = {}
304 for l in hunbiased:
305 unbiased[l]=hunbiased[l].GetSumOfWeights()
306 for l in hbiased:
307 if l.find('proj') < 0 and l.find('test') < 0 and l.find('pids') < 0 : biased[l]=hbiased[l].GetSumOfWeights()
308 p={}
309 for i in range( 1,hbiased['pids'].GetNbinsX() + 1 ):
310 c = hbiased['pids'].GetBinContent(i)
311 if c>0: p[int(hbiased['pids'].GetBinCenter(i))]=c
312
313 sorted_p = sorted(p.items(), key=operator.itemgetter(1))
314 for p in sorted_p:
315 print("%25s : %5.2G"%(pdg.GetParticle(p[0]).GetName(),float(p[1])))
316 sorted_pr = sorted(biased.items(), key=operator.itemgetter(1))
317 print("origin of muons")
318 for p in sorted_pr:
319 if not p[0].find('Hadronic inelastic')<0:
320 if len(p[0])>len( 'Hadronic inelastic' ): continue
321 denom = 0
322 if p[0] in unbiased: denom = unbiased[p[0]]
323 if denom >0:
324 fac = float(p[1])/denom
325 print("%40s : %5.2G %5.1F"%(p[0],float(p[1]),fac))
326 else:
327 print("%40s : %5.2G "%(p[0],float(p[1])))
328
329if len(sys.argv)>1:
330 runMin=sys.argv[1]
331 runMax=sys.argv[2]
332 if len(sys.argv)>3: charm=sys.argv[3]
333 compactify(charm)
334else:
335# production without boost factor
336#runMin = 1000000
337#runMax = 1000600
338# production with boost factor 100
339# runMin = 2001001
340# runMax = 2001500 # more does not work, maybe limit of hadd? 2001359
341# charm run_fixedTarget_9000119/pythia8_Geant4_9000219_10.0.root
342 print("methods to run: makeHistos(rfile),makePrintout()")
343
344# yandex compactification, Feb 8
345#... startDate = datetime.datetime(2018, 1, 1, 0, 0)
346#... endDate = datetime.datetime(2018, 2, 1, 0, 0)
347# ... goodAndBadRuns_2018-01-01_2018-02-01.pkl
348# hadd Source file 339: root://eospublic.cern.ch//eos/experiment/ship/skygrid/background-prod-2018/424121/80eefea19104_run_fixedTarget_1/pythia8_Geant4_1_10.0.root
349# hadd Target path: pythia8_Geant4_10.0_c22000.root:/
350# len(database['goodruns']) .22339
351# database['goodruns'][22338] '/eos/experiment/ship/skygrid/background-prod-2018/424121/80eefea19104_run_fixedTarget_1/pythia8_Geant4_1_10.0.root'
352#... startDate = datetime.datetime(2018, 2, 1, 0, 0)
353#... endDate = datetime.datetime(2018, 2, 6, 0, 0)
354
355# for charm
356# chicc=1.7e-3; //prob to produce primary ccbar pair/pot
357# 20*2000000 are equal to 23.5E9
358# for beauty
359# chibb=1.6e-7;
360
361# something went wrong Yandex2018Prod-23000.root onwards up to 43000, only have of yield, 44000 ok !!!
362def removeStupidFiles():
363 path = '/eos/experiment/ship/data/Mbias/background-prod-2018'
364 for f in os.listdir(path):
365 ff = path+'/'+f
366 x = os.path.getsize(ff)
367 if x < 500 : os.remove(ff)
368
369
370def check4DoubleRuns():
371 allRuns = ['goodAndBadRuns_2018-01-01_2018-02-01.pkl','goodAndBadRuns_2018-02-01_2018-02-06.pkl','goodAndBadRuns_2018-02-06_2018-02-12.pkl',\
372 'goodAndBadRuns_2018-02-12_2018-02-15.pkl','goodAndBadRuns_2018-02-15_2018-02-21.pkl','goodAndBadRuns_2018-02-21_2018-02-26.pkl']
373 Nruns=0
374 for x in allRuns:
375 fn = open(x)
376 dn = pickle.load(fn)
377 Nruns += len(dn['goodruns'])
378 print("Total number of runs:",Nruns)
379
380 for n in range( len(allRuns)-1 ):
381 fn = open(allRuns[n])
382 dn = pickle.load(fn)
383 for m in range( n+1, len(allRuns) ):
384 fm = open(allRuns[m])
385 dm = pickle.load(fm)
386 for rn in dn['goodruns']:
387 for rm in dm['goodruns']:
388 if rn == rm :
389 print("double entry found",rn,allRuns[n],allRuns[m])
390
391
392