SND@LHC Software
Loading...
Searching...
No Matches
ana_thermalNeutrons.py
Go to the documentation of this file.
1#!/usr/bin/env python
2import ROOT
3from decorators import *
4
5import os
6import shipunit as u
7from argparse import ArgumentParser
8
9parser = ArgumentParser()
10parser.add_argument("-r", "--run", dest="run", type=int, help="production sequence number", default=0)
11parser.add_argument('--material', dest='targetMaterial', type=str,help="target material BoronCarbide or Boratedpolyethylene", default="BoronCarbide")
12parser.add_argument('--thick', dest='targetLength', type=float,help="target thickness", default=0.1)
13parser.add_argument('-c', '--command', dest='command', help="command",default="none")
14parser.add_argument("-f", dest="inputFile", help="Input file", required=False, default=False)
15parser.add_argument("-g", dest="geoFile", help="geo file", required=False, default="geofile-thermNeutron_Boratedpolyethylene_coldbox.root")
16
17options = parser.parse_args()
18
19h={}
20import rootUtils as ut
21def count(hFile):
22 f=ROOT.TFile(hFile)
23 ut.bookHist(h,'E',';log10(Ekin [MeV])',1000,-12.,4.)
24 ut.bookHist(h,'Epassed',';log10(Ekin [MeV])',1000,-12.,4.)
25 ut.bookHist(h,'Esecondary',';log10(Ekin [MeV])',1000,-12.,4.)
26 ut.bookHist(h,'electrons',';log10(Ekin [MeV])',1000,-12.,4.,100,-0.5,99.5)
27 ut.bookHist(h,'photons',';log10(Ekin [MeV])',1000,-12.,4.,100,-0.5,99.5)
28 ut.bookHist(h,'Lab',' ; logE[MeV];absorption Length [cm]',1000,-12.,4.,100000,-0.1,100.)
29 ut.bookHist(h,'Labz','; logE[MeV];absorption Length [cm]',1000,-12.,4.,100000,-0.1,100.)
30 ut.bookHist(h,'xyz','', 100,-0.1,0.1,100,-0.1,0.1,200,-1.,1.)
31 ut.bookHist(h,'dxyz','',100,-0.1,0.1,100,-0.1,0.1,200,-1.,1.)
32 Npassed = 0
33 for sTree in f.cbmsim:
34 N = sTree.MCTrack[0]
35 Ekin = N.GetP()**2/(2*N.GetMass())*1000.
36 logEkin = ROOT.TMath.Log10(Ekin)
37 rc = h['E'].Fill(logEkin)
38 Nelectrons = {11:0,22:0}
39 for m in sTree.MCTrack:
40# only count electrons if produced in first cm, don't want neutrons to thermalize
41 if m.GetStartZ()>1: continue
42 for x in Nelectrons:
43 if abs(m.GetPdgCode())==x: Nelectrons[x]+=1
44 rc = h['electrons'].Fill(logEkin,Nelectrons[11])
45 rc = h['photons'].Fill(logEkin,Nelectrons[22])
46 for p in sTree.vetoPoint:
47 if not p.PdgCode()==2112: continue
48 if not p.GetTrackID()==0: rc = h['Esecondary'].Fill(ROOT.TMath.Log10(Ekin))
49 if p.GetDetectorID()==1 :
50 mean = ROOT.TVector3(p.GetX(),p.GetY(),p.GetZ())
51 end = p.LastPoint()
52 D = 2*(end-mean)
53 start = 2*mean-end
54 rc = h['Lab'].Fill(logEkin,D.Mag())
55 rc = h['Labz'].Fill(logEkin,D.Z())
56 rc = h['xyz'].Fill(start.X(),start.Y(),start.Z())
57 rc = h['dxyz'].Fill(end.X(),end.Y(),end.Z())
58 else:
59 Npassed+=1
60 rc = h['Epassed'].Fill(ROOT.TMath.Log10(Ekin))
61 h['Eff'] = ROOT.TEfficiency(h['Epassed'],h['E'])
62 h['Eff'].Draw()
63 print(Npassed)
64
65def beamGas():
66 f=open("/mnt/hgfs/microDisk/CERNBOX/SND@LHC/ThermalNeutrons/neutronsTI18.dat")
67 variables = "Eleft:Eright:N:sig"
68 fntuple = ROOT.TFile.Open('neutronsTI18.root', 'RECREATE')
69 nt = ROOT.TNtuple("nt","neutronGas",variables)
70 for l in f.readlines():
71 if not l.find('#')<0: continue
72 p = []
73 for x in l.split(' '):
74 if x == '' or x == '\n': continue
75 p.append(float(x))
76 nt.Fill(p[0]*1000.,p[1]*1000.,p[2]/1000.*10170.,p[3]) # / 1000 for MeV
77 # figure 12 of technical proposal
78 fntuple.cd()
79 nt.Write()
80 nt.Draw('log10(N*Eleft):log10(Eleft)','','box')
81
82def myPrint(tc,tname):
83 for z in ['.png','.pdf','.root']:
84 tc.Print(tname+z)
85
86materials = {'Boratedpolyethylene':ROOT.kGreen,'BoronCarbide':ROOT.kBlue,'Concrete':ROOT.kGray,'EmulsionAg109':ROOT.kOrange}
87
88def absorptionLength(plotOnly=True):
89 materials = {'Boratedpolyethylene':ROOT.kGreen,'BoronCarbide':ROOT.kBlue,'Concrete':ROOT.kGray,'EmulsionAg109':ROOT.kOrange,'H2O':ROOT.kCyan}
90 B10Parameter = {}
91 B10Parameter ['Boratedpolyethylene'] = 0.01*0.94 / (10.*1.672E-24)
92 B10Parameter ['BoronCarbide'] = 0.125*1.360 / (10.*1.672E-24)
93 Xsec = 0.61E-24
94 h['Fabsorp'] = {}
95 Fabsorp = h['Fabsorp']
96 for m in B10Parameter:
97 Fabsorp[m] = ROOT.TF1('fabs'+m,'1./([0]/sqrt(10**x)*[1])',-10.,4.)
98 Fabsorp[m].SetParameter(0,Xsec)
99 Fabsorp[m].SetParameter(1,B10Parameter[m])
100 if not plotOnly:
101 for material in materials:
102 for Erange in ['-14_-12','-12_-10','-10_-8','-8_-7','-7_-6','-6_-4','-4_-2']:
103 fname = "thermNeutron_"+material+"_100.0_"+Erange+"_0.root"
104 if not fname in os.listdir('.'): continue
105 count(fname)
106 h['Esecondary'+'_'+material+Erange] = h['Esecondary'].Clone('Esecondary'+'_'+material+Erange)
107 h['electrons'+'_'+material+Erange] = h['electrons'].Clone('electrons'+'_'+material+Erange)
108 h['photons'+'_'+material+Erange] = h['photons'].Clone('photons'+'_'+material+Erange)
109 for x in ['Lab','Labz']:
110 hname = x+'_'+material+Erange
111 h[hname] = h[x].ProfileX(hname,1,-1,'g')
112 hsum = x+'_'+material
113 if not hsum in h:
114 h[hsum] = h[hname].Clone(hsum)
115 h[hsum].SetLineColor(materials[material])
116 else: h[hsum].Add(h[hname])
117 ROOT.gROOT.cd()
118 ut.writeHists(h,'thermalNeutrons-histograms.root',plusCanvas=True)
119 else:
120 ut.readHists(h,'thermalNeutrons-histograms.root')
121# make stats about secondary neutrons:
122 for material in materials:
123 for Erange in ['-14_-12','-12_-10','-10_-8','-8_-7','-7_-6','-6_-4','-4_-2']:
124 histo = h['Esecondary'+'_'+material+Erange]
125 print("secondary neutrons for %s %s %5.2F%%"%(material,Erange,100.*histo.GetEntries()/h['Labz'+'_'+material+Erange].GetEntries()))
126 for X in ['electrons','photons']:
127 for material in materials:
128 h[X+'_'+material]=h[X].Clone(X+'_'+material)
129 for Erange in ['-14_-12','-12_-10','-10_-8','-8_-7','-7_-6','-6_-4','-4_-2']:
130 h[X+'_'+material].Add(h[X+'_'+material+Erange])
131 norm = h[X+'_'+material].ProjectionX('norm',1,101)
132 for p in [1,2,5,10]:
133 h[X+str(p)+'Percent_'+material] = h[X+'_'+material].ProjectionX(X+str(p)+'Percent_'+material,1,p)
134 h[X+str(p)+'Percent_'+material].Divide(norm)
135 ut.bookCanvas(h,'T'+X,'',1200,800,1,1)
136 h['T'+X].cd()
137 material = 'EmulsionAg109'
138 h[X+'1Percent_'+material].SetLineColor(ROOT.kRed)
139 h[X+'2Percent_'+material].SetLineColor(ROOT.kOrange)
140 h[X+'5Percent_'+material].SetLineColor(ROOT.kBlue)
141 h[X+'10Percent_'+material].SetLineColor(ROOT.kGreen)
142 h[X+'1Percent_'+material].SetTitle('')
143 h[X+'1Percent_'+material].GetYaxis().SetTitle('fraction of events with < N_{'+X+'}')
144 h[X+'1Percent_'+material].GetXaxis().SetRangeUser(-12.,1)
145 h[X+'1Percent_'+material].SetStats(0)
146 h[X+'1Percent_'+material].Draw('hist')
147 h['leg'+X]=ROOT.TLegend(0.63,0.25,0.88,0.40)
148 for p in [1,2,5,10]:
149 h[X+str(p)+'Percent_'+material].Draw('histsame')
150 rc = h['leg'+X].AddEntry(h[X+str(p)+'Percent_'+material],'N_{'+X+'}<'+str(p),'PL')
151 h['leg'+X].Draw('same')
152 myPrint(h['T'+X],'fracEveWith'+X)
153#
154 fntuple = ROOT.TFile.Open('neutronsTI18.root')
155 nt=fntuple.nt
156 ROOT.gROOT.cd()
157 tcanv = 'TFig12'
158 if tcanv in h: h.pop(tcanv)
159 ut.bookCanvas(h,tcanv,'',1200,800,1,1)
160 # figure 12 of technical proposal
161 nt.Draw('log10(N*Eleft):log10(Eleft)>>rates','','box')
162 h['rates']=ROOT.gROOT.FindObjectAny('rates').Clone('rates')
163
164 for x in ['Lab','Labz']:
165 tcanv = 'abs'+x
166 if tcanv in h: h.pop(tcanv)
167 ut.bookCanvas(h,tcanv,'',1200,800,1,1)
168 h['abs'+x].cd()
169 h['abs'+x].SetLogy()
170 hsum = x+'_Concrete'
171 h[hsum].GetXaxis().SetRangeUser(-12.,1)
172 h[hsum].SetMaximum(30.)
173 h[hsum].SetMinimum(0.004)
174 h[hsum].GetXaxis().SetTitle('logE[MeV]')
175 h[hsum].GetYaxis().SetTitle('absorption Length [cm]')
176 h[hsum].SetLineWidth(2)
177 h[hsum].Draw('hist')
178 h['leg'+x]=ROOT.TLegend(0.6,0.2,0.86,0.36)
179 for material in materials:
180 hsum = x+'_'+material
181 h[hsum].SetStats(0)
182 h[hsum].SetLineWidth(2)
183 h[hsum].Draw('histsame')
184 if material in Fabsorp: Fabsorp[material].Draw('same')
185 rc = h['leg'+x].AddEntry(h[hsum],material,'PL')
186 h['leg'+x].Draw('same')
187 myPrint(h['abs'+x],'AbsLength'+x)
188
189# folding with neutron rate
190 h['Fig12'] = ROOT.TGraph()
191 h['dE'] = ROOT.TGraph()
192 h['neutronRate'] = ROOT.TGraph()
193 h['dangerZone']=ROOT.TGraph()
194 h['dangerZone'].SetPoint(0,-5.6,0.)
195 h['dangerZone'].SetPoint(1,-5.6,1.E7)
196 h['dangerZone'].SetPoint(2,-5.1,1.E7)
197 h['dangerZone'].SetPoint(3,-5.1,0.)
198 h['dangerZone'].SetFillStyle(1001)
199 h['dangerZone'].SetFillColor(ROOT.kYellow)
200
201 n = 0
202 RateIntegrated = 0
203 RateIntegratedW = 0
204 for nt in fntuple.nt:
205 E = (nt.Eleft+nt.Eright)/2.
206 dE = nt.Eright - nt.Eleft
207 h['Fig12'].SetPoint(n,ROOT.TMath.Log10(E),ROOT.TMath.Log10(nt.N*E))
208 h['neutronRate'].SetPoint(n,E,nt.N)
209 h['dE'].SetPoint(n,E,dE)
210 RateIntegrated+=nt.N
211 RateIntegratedW+=nt.N*dE
212 n+=1
213
214 h['TFig12'].cd()
215 h['Fig12'].Draw('same')
216#
217 ut.bookHist(h,'Nr',';E [MeV];dn/dlogE [cm^{-2}y^{-1}] ',100,-12.,1.)
218 h['Nr'].SetMaximum(2.E8)
219 h['Nr'].SetMinimum(1.E-4)
220 h['Nr'].SetStats(0)
221 intRates = {}
222 thick = {0.0:ROOT.kBlue,0.5:ROOT.kOrange,1:ROOT.kRed,2:ROOT.kRed+2,5:ROOT.kRed-7,10:ROOT.kGreen}
223 for material in ['Boratedpolyethylene','BoronCarbide']:
224 intRates[material]={}
225 tcanv = 'ratesWith_'+material
226 if tcanv in h: h.pop(tcanv)
227 ut.bookCanvas(h,tcanv,material,1200,800,1,1)
228 h[tcanv].SetLogy(1)
229 h[tcanv].cd()
230 h['ratesWith_'+material].cd()
231 h['Nr'].Draw()
232 h['dangerZone'].Draw('sameF')
233 h['legthick'+material]=ROOT.TLegend(0.6,0.2,0.86,0.36)
234 for d in thick: # try different thicknesses
235 intRates[material][d]=0
236 absorbLength = h['Labz_'+material]
237 gname = 'neutronRate_'+material+'_'+str(d)
238 h[gname] = ROOT.TGraph()
239 h[gname].SetLineWidth(2)
240 h[gname].SetLineColor(thick[d])
241 for n in range(h['Fig12'].GetN()):
242 logE = h['Fig12'].GetX()[n]
243 R = ROOT.TMath.Power(10.,h['Fig12'].GetY()[n])
244 dE = h['dE'].GetPointY(n)
245 E = ROOT.TMath.Power(10.,logE)
246 absorpt = 0
247 if d==0.0:
248 h[gname].SetPoint(n,logE,R)
249 intRates[material][d] += dE*R/E
250 else:
251 L = absorbLength.GetBinContent(absorbLength.FindBin(logE))
252 Rnew = 0
253 if L>0:
254 absorpt = ROOT.TMath.Exp(-d/L)
255 Rnew = R * absorpt
256 h[gname].SetPoint(n,logE,Rnew)
257# E = (nt.Eleft+nt.Eright)/2. h['Fig12'].SetPoint(n,ROOT.TMath.Log10(E),ROOT.TMath.Log10(nt.N*E))
258 intRates[material][d] += dE * Rnew/E
259 h[gname].Draw('same')
260 rc = h['legthick'+material].AddEntry(h[gname],'thickness %3.1Fcm'%(d),'PL')
261 reduction = intRates[material][d]/intRates[material][0]
262 print('integrated rate with %s d=%3.1Fcm: %6.4G reduction factor=%6.2G'%(material,d,intRates[material][d],reduction) )
263# make integrated rates:
264 h['IUp-neutronRate_'+material+str(d)] = ROOT.TGraph()
265 h['IUp-neutronRateW_'+material+str(d)] = ROOT.TGraph()
266 h['IDown-neutronRate_'+material+str(d)] = ROOT.TGraph()
267 up = h['IUp-neutronRate_'+material+str(d)]
268 down = h['IDown-neutronRate_'+material+str(d)]
269 N = h[gname].GetN()
270 S = 0
271 for n in range(N):
272 logE = h[gname].GetX()[n]
273 dE = h['dE'].GetPointY(n)
274 Rnew = h[gname].GetY()[n]
275 E = ROOT.TMath.Power(10.,logE)
276 S += dE * Rnew/E
277 up.SetPoint(n,logE,S)
278# with damage function
279 W = 1 - h['electrons1Percent_EmulsionAg109'].GetBinContent(n)
280 h['IUp-neutronRateW_'+material+str(d)].SetPoint(n,logE,W*S)
281 S = up.GetY()[N-1]
282 for n in range(N):
283 logE = up.GetX()[n]
284 S -= up.GetY()[n]
285 down.SetPoint(n,logE,S)
286#
287 for Elimit in [100.,1.,0.1,0.01]: #MeV
288 g = h['IUp-neutronRate_'+material+str(d)]
289 gW = h['IUp-neutronRateW_'+material+str(d)]
290 N = g.GetN()
291 for n in range(N-1,1,-1):
292 if ROOT.TMath.Power(10,g.GetX()[n])<Elimit: break
293 reduction = g.GetY()[n]/h['IUp-neutronRate_'+material+str(0.0)].GetY()[N-1]
294 reductionW = gW.GetY()[n]/h['IUp-neutronRate_'+material+str(0.0)].GetY()[N-1]
295 print('integrated rate E < %5.3F with %s d=%3.1Fcm: %6.4G reduction factor=%6.2G | %6.4G reduction factor=%6.2G'%(
296 Elimit,material,d,g.GetY()[n],reduction,gW.GetY()[n],reductionW))
297
298 h['legthick'+material].Draw('same')
299 myPrint(h['ratesWith_'+material],'reducedRates_'+material)
300
301def printReactions():
302 for material in materials:
303 for Erange in ['-14_-12','-12_-10','-10_-8','-8_-7','-7_-6','-6_-4','-4_-2']:
304 print('reactions for material ',material, 'in energy range ',Erange)
305 fname = "thermNeutron_"+material+"_100.0_"+Erange+"_0.root"
306 reactions(fname,1)
307def reactions(hFile,distance=1E10):
308 f=ROOT.TFile(hFile)
309 elements = {1:'Hydrogen',2:'Helium',3:'Lithium',4:'Berylium',5:'Boron',7:'Nitrogen',6:'Carbon',
310 8:'Oxygen',11:'Sodium',12:'Magnesium ',14:'Silicon',15:'Phosphorus',16:'Sulfur',20:'Cadmium',
311 13:'Aluminum',26:'Iron',34:'Selenium',35:'Bromine',47:'Silver',53:"Iodine"}
312 ut.bookHist(h,'E',';log10(Ekin [MeV])',1000,-12.,4.)
313 stats = {}
314 n=-1
315 for sTree in f.cbmsim:
316 n+=1
317 daughters = []
318 for m in sTree.MCTrack:
319 if m.GetMotherId()==0: daughters.append(m)
320 nphot = 0
321 for d in daughters:
322 if d.GetPdgCode()==22: nphot+=1
323 isotopes =[]
324 for d in daughters:
325 if d.GetStartZ()>distance: continue
326 pid = d.GetPdgCode()
327 if pid> 1000000000: isotopes.append(pid)
328 if len(isotopes)==0: continue
329 isotopes.sort()
330 reaction=''
331 for x in isotopes: reaction+=str(x)+','
332 if not reaction in stats:
333 sTree.MCTrack.Dump()
334 stats[reaction]={}
335 if not nphot in stats[reaction]: stats[reaction][nphot]=0
336 stats[reaction][nphot]+=1
337 if n>50000: break
338 for reaction in stats:
339 ilist = ''
340 R=0
341 for c in reaction.split(','):
342 if len(c)<2: continue
343 x = int(c)
344 A = int((x/10)%1000)
345 Z = int((x/10000)%1000)
346 isoname = "unknown"
347 if Z in elements: isoname = elements[Z]
348 ilist+= " %s(%i,%i)"%(isoname,A,Z)
349 for nphot in stats[reaction]: R+=stats[reaction][nphot]
350 R = R/n
351 print('nuclear reaction n + %s: %5.3F%%'%(ilist,R*100))
352 return stats
353
354def absorptionLengthOLD():
355 Lfun = ROOT.TF1('Lfun','[0]*(10**x)**[1]',-9,-6)
356 Lfun.SetParameter(0,6.4)
357 Lfun.SetParameter(1,0.98)
358 hFiles = {"thermNeutron_BoronCarbide_X.XX_-E_-E_0.root":[0.08,0.3],"thermNeutron_Boratedpolyethylene_X.XX_0.root":[1.0,100.]}
359 # thermNeutron_BoronCarbide_4.0_-8_-7_0.root
360
361 Ls = {0.01:ROOT.kRed,0.1:ROOT.kOrange,0.04:ROOT.kCyan,0.4:ROOT.kBlue,4.0:ROOT.kMagenta}
362
363 for hFile in hFiles:
364 material = hFile.split('_')[1]
365 ut.bookCanvas(h,'absorb'+material,'',1200,800,1,1)
366 h['absorb'+material].cd()
367 for L in Ls:
368 l = str(L)
369 if L<3: tmp = hFile.replace("X.XX",l).replace(" _-E_-E","")
370 else: tmp = hFile.replace("X.XX",l).replace(" _-E_-E","_-8_-7")
371 count(tmp)
372 h['Eff_'+l]=h['Eff'].Clone('Eff_'+l)
373 h['L_'+l]=ROOT.TGraph()
374 h['L_'+l].SetLineColor(Ls[L])
375 h['Eff'].Draw()
376 g = h['Eff'].GetPaintedGraph()
377 for n in range(g.GetN()):
378 logE, p = g.GetPointX(n),g.GetPointY(n)
379 if p>0:
380 absL = -L/ROOT.TMath.Log(p)
381 else:
382 absL = 0
383 h['L_'+l].SetPoint(n,logE,absL)
384 ut.bookHist(h,'L',';logE; L [cm]',100,-9.,-6.)
385 h['L'].SetMaximum(hFiles[hFile][0])
386 h['L'].SetStats(0)
387 h['L'].Draw()
388 for L in Ls:
389 if L>hFiles[hFile][1]: continue
390 h['L_'+str(L)].Draw('same')
391 h['L_'+str(0.1)].Fit(Lfun)
392 myPrint(h['absorb'+material],'absorbLength'+material)
393
394def flukaRateIntegrated(save=False):
395 fntuple = ROOT.TFile.Open('neutronsTI18.root')
396 h['Fig12'] = ROOT.TGraph()
397 h['neutronRate'] = ROOT.TGraph()
398 h['N'] = ROOT.TGraph()
399 n = 0
400 for nt in fntuple.nt:
401 # nt.N = [cm/MeV]
402 E = (nt.Eleft+nt.Eright)/2.
403 dE = nt.Eright - nt.Eleft
404 h['Fig12'].SetPoint(n,ROOT.TMath.Log10(E),E*nt.N)
405 h['neutronRate'].SetPoint(n,ROOT.TMath.Log10(E),dE*nt.N)
406 h['N'].SetPoint(n,E,dE*nt.N)
407 n+=1
408 S= 0
409 h['neutronRate<E'] = h['neutronRate'].Clone('S')
410 for n in range(h['neutronRate'].GetN()-1,-1,-1):
411 S+=h['neutronRate'].GetY()[n]
412 h['neutronRate<E'].SetPoint(n,h['neutronRate'].GetX()[n],S)
413 h['neutronRate<E%']= h['neutronRate'].Clone('S%')
414 for n in range(h['neutronRate<E'].GetN()):
415 h['neutronRate<E%'].SetPoint(n,h['neutronRate<E'].GetX()[n],h['neutronRate<E'].GetY()[n]/h['neutronRate<E'].GetY()[0])
416 ut.bookHist(h,'Nr',';E [MeV];dn/dlogE [cm^{-2}y^{-1}] ',100,-12.,2.)
417 h['Nr'].SetMaximum(1.1)
418 h['Nr'].SetMinimum(0.)
419 h['Nr'].SetStats(0)
420 h['Nr'].SetTitle('; log10(E [MeV]();N(E>x)/total')
421 h['neutronRate<E%'].SetLineWidth(2)
422 h['Nr'].Draw()
423 h['neutronRate<E%'].Draw('same')
424 if save: ut.writeHists(h,'flukaNeutronRates')
425
426from array import array
427neutronMass = 939.565379/1000.
428hnorm = {}
429def coldBox(plotOnly=True,pas=''):
430
431 if 1>0:
432 tmp = options.inputFile.split('/')
433 gFile = "geofile-"+(tmp[len(tmp)-1].split('_coldbox')[0]+'_coldbox.root').replace('histos-','')
434 g = ROOT.TFile(gFile)
435 sGeo = g.FAIRGeom
436 ROOT.gROOT.cd()
437 vbox = sGeo.FindVolumeFast('vbox')
438 sbox = sGeo.FindVolumeFast('sensBox')
439 dX,dY,dZ = vbox.GetShape().GetDX(),vbox.GetShape().GetDY(),vbox.GetShape().GetDZ()
440 nav = sGeo.GetCurrentNavigator()
441 boundaries = {}
442 # find y boundaries
443 start = array('d',[0,200,0])
444 direction = array('d',[0,-1,0])
445 startnode = sGeo.InitTrack(start, direction)
446 length = c_double(200.)
447 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
448 boundaries['topIn'] = nav.GetCurrentPoint()[1]
449 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
450 boundaries['topSens'] = nav.GetCurrentPoint()[1]
451 start = array('d',[0,-200,0])
452 direction = array('d',[0,1,0])
453 startnode = sGeo.InitTrack(start, direction)
454 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
455 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
456 boundaries['botSens'] = nav.GetCurrentPoint()[1]
457 # find x boundaries
458 start = array('d',[-200,0,0])
459 direction = array('d',[1,0,0])
460 startnode = sGeo.InitTrack(start, direction)
461 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
462 boundaries['leftIn'] = nav.GetCurrentPoint()[0]
463 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
464 boundaries['leftSens'] = nav.GetCurrentPoint()[0]
465 start = array('d',[200,0,0])
466 direction = array('d',[-1,0,0])
467 startnode = sGeo.InitTrack(start, direction)
468 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
469 boundaries['rightIn'] = nav.GetCurrentPoint()[0]
470 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
471 boundaries['rightSens'] = nav.GetCurrentPoint()[0]
472 # find z boundaries
473 start = array('d',[0,0,-200])
474 direction = array('d',[0,0,1])
475 startnode = sGeo.InitTrack(start, direction)
476 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
477 boundaries['backIn'] = nav.GetCurrentPoint()[2]
478 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
479 boundaries['backSens'] = nav.GetCurrentPoint()[2]
480 start = array('d',[0,0,200])
481 direction = array('d',[0,0,-1])
482 startnode = sGeo.InitTrack(start, direction)
483 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
484 boundaries['frontIn'] = nav.GetCurrentPoint()[2]
485 node = sGeo.FindNextBoundaryAndStep(length, ROOT.kFALSE)
486 boundaries['frontSens'] = nav.GetCurrentPoint()[2]
487 Rin = {'':'','topIn':'Y','leftIn':'X','rightIn':'X','frontIn':'Z','backIn':'Z'}
488 Rsens = {'':'','topSens':'Y','botSens':'Y','leftSens':'X','rightSens':'X','frontSens':'Z','backSens':'Z'}
489 epsi = 0.1
490# side with holes is the front
491#
492 if not plotOnly:
493 # example file "root://eospublic.cern.ch:1094//eos/experiment/sndlhc/MonteCarlo/ThermalNeutrons//thermNeutron_Borated30polyethylene_10.0_coldbox_0-20M.root"
494 f = ROOT.TFile.Open(options.inputFile)
495 ROOT.gROOT.cd()
496 ut.bookHist(h,'start','start neutron;x [cm] ;y [cm] ;z [cm]',100,-200,200,100,-200,200,100,-200,200)
497 ut.bookHist(h,'startR','start neutron;R',100,0,200)
498 for o in ['_seco','_prim']:
499 for T in ['','cold','hot']:
500 ut.bookHist(h,'entry'+T+o,'entry neutron;x [cm] ;y [cm] ;z [cm]',100,-100,100,100,-100,100,100,-100,100)
501 ut.bookHist(h,'exit'+T+'-original'+o,'enter coldbox neutron;x [cm] ;y [cm] ;z [cm]',100,-100,100,100,-100,100,100,-100,100)
502 ut.bookHist(h,'exit'+T+o,'enter coldbox neutron;x [cm] ;y [cm] ;z [cm]',100,-100,100,100,-100,100,100,-100,100)
503 for r in Rin:
504 ut.bookHist(h,'EkinE'+r+o,'log10(Ekin)',100,-13.,0,100,-200.,200.)
505 ut.bookHist(h,'Ekin'+r+o,'log10(Ekin)',100,-13.,0,100,-200.,200.)
506 for r in Rsens:
507 ut.bookHist(h,'DF'+r+o,'travel distance',100,0.,200.)
508 ut.bookHist(h,'EkinF'+r+o,'log10(Ekin) vs distance',100,-13.,0,100,-200.,200.)
509 ut.bookHist(h,'EkinF-original'+r+o,'log10(Ekin) vs distance',100,-13.,0,100,-200.,200.)
510 ut.bookHist(h,'checkBox'+r+o,'enter coldbox neutron;x [cm] ;y [cm] ;z [cm]',100,-100,100,100,-100,100,100,-100,100)
511 ut.bookHist(h,'EkinG','log10(Ekin)',100,-13.,0.)
512 ut.bookHist(h,'EkinW','log10(Ekin)',100,-13.,0.)
513 ut.bookHist(h,'EkinWlin','Ekin',100,1E-9,100*1E-9)
514 ut.bookHist(h,'multS','mult veto points shielding',100,-0.5,99.5)
515 ut.bookHist(h,'multC','mult veto points inside',100,-0.5,99.5)
516 ut.bookHist(h,'multN','mult neutrons',100,-0.5,99.5)
517
518 flukaRateIntegrated()
519 Nsim = f.cbmsim.GetEntries()
520 for sTree in f.cbmsim:
521 rc=h['multN'].Fill(sTree.MCTrack.GetEntries())
522 neutron = sTree.MCTrack[0]
523 start = ROOT.TVector3(neutron.GetStartX(),neutron.GetStartY(),neutron.GetStartZ())
524 if start.y() < boundaries['botSens']: continue
525
526 P = ROOT.TVector3(neutron.GetPx(),neutron.GetPy(),neutron.GetPz())
527 Ekin = ROOT.TMath.Sqrt(P.Mag2()+neutronMass**2) - neutronMass
528 W = h['Fig12'].Eval(ROOT.TMath.Log10(Ekin*1000)) / Nsim
529 rc = h['start'].Fill(start.Z(),start.X(),start.Y(),W)
530 rc = h['EkinW'].Fill(ROOT.TMath.Log10(Ekin),W)
531 rc = h['EkinWlin'].Fill(Ekin,W)
532 rc = h['EkinG'].Fill(ROOT.TMath.Log10(Ekin))
533 rc = h['startR'].Fill(start.Mag(),W)
534 nC,nS=0,0
535 for p in sTree.vetoPoint:
536 if p.PdgCode()!=2112: continue
537 if p.GetDetectorID()==13: nC+=1
538 else: nS+=1
539 rc = h['multC'].Fill(nC)
540 rc = h['multS'].Fill(nS)
541 trajectory = {}
542 for p in sTree.vetoPoint:
543 if p.PdgCode()!=2112: continue
544 trackID = p.GetTrackID()
545 if not trackID in trajectory: trajectory[trackID]=[]
546 trajectory[trackID].append(p)
547 for trackID in trajectory:
548 firstSens = True
549 for p in trajectory[trackID]:
550 origin = '_seco'
551 if trackID==0: origin = '_prim'
552 lastPoint = p.LastPoint()
553 mPoint = ROOT.TVector3(p.GetX(),p.GetY(),p.GetZ())
554 firstPoint = 2*mPoint-lastPoint
555 D = lastPoint-firstPoint
556 TD = firstPoint - start
557 firstMom = ROOT.TVector3(p.GetPx(),p.GetPy(),p.GetPz())
558 Ekin_entry = ROOT.TMath.Sqrt(firstMom.Mag2()+neutronMass**2) - neutronMass
559 if p.GetDetectorID()==13 and firstSens: # first point inside coldbox
560 firstSens = False
561 rc = h['exit'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
562
563 if Ekin*1E9< 10: rc = h['exitcold-original'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W) # 10eV
564 else: rc = h['exithot-original'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
565
566 if Ekin_entry*1E9< 10: rc = h['exitcold'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W) # 10eV
567 else: rc = h['exithot'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
568 rc = h['DF'+origin].Fill(TD.Mag(),W)
569 rc = h['EkinF'+origin].Fill(ROOT.TMath.Log10(Ekin_entry),D.Mag(),W)
570 rc = h['EkinF-original'+origin].Fill(ROOT.TMath.Log10(Ekin),D.Mag(),W)
571 # find location Rsens = ['','topSens','botSens','leftSens','rightSens','frontSens','backSens']
572 found = False
573 for r in Rsens:
574 if r=='': continue
575 X = eval('firstPoint.'+Rsens[r]+'()')
576 if abs(X-boundaries[r]) < epsi:
577 # if r=='botSens' and p!=trajectory[trackID][0]: continue # to enter from bottom is without crossing shield. It is more complicated!
578 found = True
579 break
580 if not found:
581 txt = ''
582 for r in Rsens:
583 if r=='': continue
584 X = eval('firstPoint.'+Rsens[r]+'()')
585 txt+= " "+r+" "+str(X)
586 print("this should no happen", txt, boundaries)
587 for P in trajectory[trackID]: print(P.GetDetectorID(),P.LastPoint().X(),P.LastPoint().Y(),P.LastPoint().Z())
588 rc = h['DF'+r+origin].Fill(TD.Mag(),W)
589 rc = h['EkinF'+r+origin].Fill(ROOT.TMath.Log10(Ekin_entry),D.Mag(),W)
590 rc = h['EkinF-original'+r+origin].Fill(ROOT.TMath.Log10(Ekin),D.Mag(),W)
591 rc = h['checkBox'+r+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
592
593 if p.GetDetectorID()==1: # inside shielding
594# check that first point is further out.
595 if firstPoint.Mag()<lastPoint.Mag(): continue
596 rc = h['Ekin'+origin].Fill(ROOT.TMath.Log10(Ekin),D.Mag(),W)
597 rc = h['EkinE'+origin].Fill(ROOT.TMath.Log10(Ekin_entry),D.Mag(),W)
598 rc = h['entry'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
599 if Ekin*1E9< 10: rc = h['entrycold'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W) # 10eV
600 else: rc = h['entryhot'+origin].Fill(firstPoint.X(),firstPoint.Y(),firstPoint.Z(),W)
601 # find location
602 for r in Rin:
603 if r=='': continue
604 X = eval('firstPoint.'+Rin[r]+'()')
605 if abs(X-boundaries[r]) < epsi: break
606 rc = h['Ekin'+r+origin].Fill(ROOT.TMath.Log10(Ekin),X,W)
607 rc = h['EkinE'+r+origin].Fill(ROOT.TMath.Log10(Ekin_entry),X,W)
608
609 tmp = options.inputFile.split('/')
610 outFile = 'histos-'+ tmp[len(tmp)-1]
611# make sum of seco and prim
612 hkeys = list(h.keys())
613 for x in hkeys:
614 if x.find( '_seco')>0:
615 hname = x.replace('_seco','')
616 h[hname]=h[x].Clone(hname)
617 h[hname].Add(h[x.replace('_seco','_prim')])
618#
619 ut.writeHists(h,outFile)
620 print('finished making histograms ','histos-noConc-'+options.inputFile)
621
622 else:
623 ut.readHists(h,options.inputFile)
624 if pas!='': ut.readHists(hnorm,'histos-thermNeutron_vacuums_0.0001_coldbox_pass1.root')
625 else: ut.readHists(hnorm,options.inputFile)
626 tmp = options.inputFile.split('_')
627 material = tmp[1]+"_coldbox/"
628 thickness = "_"+tmp[2]
629
630 binning = {}
631 for c in ['entry','exit']:
632 binning[c]={}
633 for p in ['x','y','z']:
634 tmp = h[c].Project3D(p)
635 for imin in range(1,tmp.GetNbinsX()+1):
636 if tmp.GetBinContent(imin)>0: break
637 for imax in range(tmp.GetNbinsX(),0,-1):
638 if tmp.GetBinContent(imax)>0: break
639 binning[c][p]={'min':imin,'max':imax}
640
641#
642 ytop = {'axis':'Y','proj':'xz','entry':binning['entry']['y']['max'],'exit':h['exit'].GetYaxis().FindBin(boundaries['topSens']),'xdist':dZ,'ydist':dX}
643 ybot = {'axis':'Y','proj':'xz','entry':binning['entry']['y']['min'],'exit':h['exit'].GetYaxis().FindBin(boundaries['botSens']),'xdist':dZ,'ydist':dX}
644 xLeft = {'axis':'X','proj':'yz','entry':binning['entry']['x']['min'],'exit':h['exit'].GetXaxis().FindBin(boundaries['leftSens']),'xdist':dZ,'ydist':dY}
645 xRight = {'axis':'X','proj':'yz','entry':binning['entry']['x']['max'],'exit':h['exit'].GetXaxis().FindBin(boundaries['rightSens']),'xdist':dZ,'ydist':dY}
646 zMin = {'axis':'Z','proj':'yx','entry':binning['entry']['z']['min'],'exit':h['exit'].GetZaxis().FindBin(boundaries['backSens']),'xdist':dX,'ydist':dY}
647 zMax = {'axis':'Z','proj':'yx','entry':binning['entry']['z']['max'],'exit':h['exit'].GetXaxis().FindBin(boundaries['frontSens']),'xdist':dX,'ydist':dY}
648 print('boundaries',boundaries)
649
650 # make projections
651 projections = {'Top':ytop,'Bot':ybot,'Right':xRight,'Left':xLeft,'Front':zMax,'Back':zMin}
652 h['fluences'] = {}
653 for o in ['','_prim']:
654 for T in ['','cold','hot']:
655 for Z in ['entry','exit','exit-original']:
656 tmp = Z.split('-')
657 c = tmp[0]
658 x=''
659 if len(tmp)>1: x='-'+tmp[1]
660 case = c+T+x+o
661 ut.bookCanvas(h,'t'+case,case,1200,1800,2,3)
662 k=1
663 for p in projections:
664 h['t'+case].cd(k)
665 tmp = h[case]
666 axis = eval('tmp.Get'+projections[p]['axis']+'axis()')
667 if Z=='entry': axis.SetRange(projections[p][c]-1,projections[p][c]+1)
668 else:
669 axis.SetRange(projections[p][c],projections[p][c])
670 if p=='xBot':
671 xax = tmp.GetXaxis()
672 xax.SetRange(xax.FindBin(boundaries['leftSens'])+1,xax.FindBin(boundaries['rightSens'])-1)
673 zax = tmp.GetZaxis()
674 zax.SetRange(zax.FindBin(boundaries['backSens'])+1,zax.FindBin(boundaries['frontSens'])-1)
675 h[case+p] = tmp.Project3D(projections[p]['proj'])
676 h[case+p].SetName(case+p)
677 tmp.GetXaxis().SetRange(0,0)
678 tmp.GetYaxis().SetRange(0,0)
679 tmp.GetZaxis().SetRange(0,0)
680 h[case+p].SetStats(0)
681 h[case+p].SetMinimum(0)
682 h[case+p].SetTitle(p)
683 h[case+p].Draw('colz')
684# check uniformity:
685 h['X-'+case+p]=h[case+p].ProjectionX('X-'+case+p)
686 h['Y-'+case+p]=h[case+p].ProjectionY('Y-'+case+p)
687 k+=1
688 sqcm = projections[p]['xdist']* projections[p]['ydist']
689 entries = h[case+p].GetSumOfWeights()
690 X = entries/sqcm
691 if X > 100: txt = "%5.1F/cm^{2}"%(X)
692 else: txt = "%5.2F/cm^{2}"%(X)
693 h['fluences'][case+p] =X
694 L = ROOT.TLatex()
695 rc = L.DrawLatexNDC(0.2,0.85,txt)
696 h['X-'+case+p].Scale(1./ projections[p]['ydist'])
697 h['Y-'+case+p].Scale(1./ projections[p]['xdist'])
698 myPrint(h['t'+case],material+'t'+case+thickness)
699
700# make cross checks
701 ut.bookCanvas(h,'crosschecks','cross checks',900,600,1,1)
702 tc = h['crosschecks'].cd()
703 tc.SetLogy(1)
704 tc.SetGridx()
705 tc.SetGridy()
706 h['EkinW'].SetStats(0)
707 h['EkinW'].SetStats(0)
708 h['EkinW'].SetLineWidth(3)
709 h['EkinW'].SetTitle(';log(E) [GeV];dn/dlogE [cm^{-2}y^{-1}] ')
710 h['EkinW'].SetMaximum(1E8)
711 h['EkinW'].SetMinimum(1E2)
712 h['EkinW'].Draw()
713 myPrint(h['crosschecks'],material+'kinEnergy'+thickness )
714 tc.SetLogy(0)
715 k=-10
716 for p in projections:
717 for l in ['X-','Y-']:
718 case = l+'entry'+p
719 h[case].SetLineColor(ROOT.kRed+k)
720 h[case].SetStats(0)
721 if k<-9:
722 h[case].SetTitle('; x,y,z [cm]; N/L [cm^{-1}]')
723 tpl = ut.findMaximumAndMinimum(h[case])
724 h[case].SetMaximum(tpl[1]*1.5)
725 h[case].Draw()
726 else: h[case].Draw('same')
727 k+=1
728 myPrint(h['crosschecks'],material+'irradiationXYZ'+thickness )
729#
730 tc.SetLogy(1)
731# Ekin Rin = {'':'','topIn':'Y','leftIn':'X','rightIn':'X','frontIn':'Z','backIn':'Z'}
732# EkinF Rsens = {'':'','topSens':'Y','botSens':'Y','leftSens':'X','rightSens':'X','frontSens':'Z','backSens':'Z'}
733 ut.bookCanvas(h,'trej','rejections',1200,1800,2,3)
734 k=0
735 for r in ['','top','bot','right','left','front','back']:
736 rej = 'rej'+r
737 rejo = 'rejo'+r
738 if r=='':
739 norm = hnorm['Ekin'].ProjectionX('norm')
740 h[rej] = h['EkinF'].ProjectionX(rej)
741 h[rejo] = h['EkinF-original'].ProjectionX(rejo)
742 else:
743 k+=1
744 if r=='bot': norm = hnorm['EkintopIn'].ProjectionX('norm')
745 else: norm = hnorm['Ekin'+r+'In'].ProjectionX('norm')
746 h[rej] = h['EkinF'+r+'Sens'].ProjectionX(rej)
747 h[rejo] = h['EkinF-original'+r+'Sens'].ProjectionX(rejo)
748 h[rej].Divide(norm)
749 h[rejo].Divide(norm)
750 h[rej].SetStats(0)
751 h[rejo].SetStats(0)
752 h[rej].SetTitle(';log10(Ekin) GeV; rejection')
753 h[rej].SetLineColor(ROOT.kBlue)
754 h[rej].SetLineWidth(2)
755 h[rejo].SetLineWidth(2)
756 h[rejo].SetLineColor(ROOT.kGreen)
757 h[rej].GetXaxis().SetRangeUser(-13.,-1.)
758 h[rej].SetMaximum(1.2)
759 h[rej].SetMinimum(1.E-4)
760 if k==0: tc = h['crosschecks'].cd()
761 else:
762 tc = h['trej'].cd(k)
763 tc.SetGridx()
764 tc.SetGridy()
765 h[rej].SetMinimum(1.E-6)
766 h[rej].SetTitle(r[0].upper()+r[1:])
767 tc.SetLogy()
768 h[rej].Draw('hist')
769 h[rejo].Draw('histsame')
770 h['legR'+r]=ROOT.TLegend(0.11,0.74,0.72,0.82)
771 rc = h['legR'+r].AddEntry(h[rejo],'reduction as function of original E_{kin} when entering shield','PL')
772 rc = h['legR'+r].AddEntry(h[rej],'reduction as function of E_{kin} when leaving shield','PL')
773 h['legR'+r].Draw('same')
774 if k==0: myPrint(h['crosschecks'],material+'rejections'+thickness)
775 myPrint(h['trej'],material+'Listrejections'+thickness+pas)
776#
777 h['dangerZone']=ROOT.TGraph()
778 h['dangerZone'].SetPoint(0,-8.6,0.)
779 h['dangerZone'].SetPoint(1,-8.6,1.)
780 h['dangerZone'].SetPoint(2,-8.1,1.)
781 h['dangerZone'].SetPoint(3,-8.1,0.)
782 h['dangerZone'].SetFillStyle(1001)
783 h['dangerZone'].SetFillColor(ROOT.kYellow)
784 for r in ['rej','rejo']:
785 ut.bookCanvas(h,r+'ection2','rejectionRate',900,600,1,1)
786 h[r+'TopLeftRightBack'] = h[r+'top'].Clone(r+'TopLeftRightBack')
787 h[r+'TopLeftRightBack'].Add(h[r+'left'])
788 h[r+'TopLeftRightBack'].Add(h[r+'right'])
789 h[r+'TopLeftRightBack'].Add(h[r+'back'])
790 h[r+'TopLeftRightBack'].Scale(1./4.)
791 h[r+'TopLeftRightBack'].SetTitle('')
792 if r=='rej': h[r+'TopLeftRightBack'].GetXaxis().SetTitle('outgoing '+h[r+'TopLeftRightBack'].GetXaxis().GetTitle())
793 if r=='rejo': h[r+'TopLeftRightBack'].GetXaxis().SetTitle('incoming '+h[r+'TopLeftRightBack'].GetXaxis().GetTitle())
794 h[r+'TopLeftRightBack'].SetLineColor(ROOT.kGreen)
795 h[r+'bot'].SetLineColor(ROOT.kGray)
796 h[r+'front'].SetLineColor(ROOT.kRed)
797 h[r+'TopLeftRightBack'].SetMaximum(1.2)
798 h[r+'TopLeftRightBack'].SetMinimum(1.E-6)
799 tc = h[r+'ection2'].cd()
800 tc.SetGridx()
801 tc.SetGridy()
802 tc.SetLogy()
803 h[r+'TopLeftRightBack'].Draw()
804 h['dangerZone'].Draw('sameF')
805 h[r+'TopLeftRightBack'].Draw('same')
806 h[r+'bot'].Draw('same')
807 h[r+'front'].Draw('same')
808 h[r+'legR2']=ROOT.TLegend(0.55,0.30,0.86,0.45)
809 rc = h[r+'legR2'].AddEntry(h[r+'TopLeftRightBack'],'av. rejection top/left/right/back','PL')
810 rc = h[r+'legR2'].AddEntry(h[r+'bot'],'rejction bottom, only concrete','PL')
811 rc = h[r+'legR2'].AddEntry(h[r+'front'],'rejction front, with holes','PL')
812 h[r+'legR2'].Draw('same')
813 myPrint(h[r+'ection2'],material+'Sum '+r+'ections'+thickness+pas)
814
815 if pas=='': ut.writeHists(h,options.inputFile.replace('.root','_pass1.root'))
816
817# statistics:
818 for o in ['','-original']:
819 for T in ['cold','hot']:
820 norm = 0
821 exit = 0
822 concrete = 0
823 holes = 0
824 for p in projections:
825 if p!='Bot' and norm!='Front':
826 norm+=h['fluences']['entry'+T+p]
827 exit+=h['fluences']['exit'+T+o+p]
828 if p=='Bot': concrete = h['fluences']['exit'+T+o+p]
829 elif p=='Front': holes = h['fluences']['exit'+T+o+p]
830 norm = norm/float(4)
831 exit = exit/float(4)
832#! do not need average, need total.
833# incoming = 6 * average !
834 print('%s %s region: %5.2F concrete: %5.2F holes: %5.2F x permille total: %5.2F x permille'%(o,T,exit/norm*1000,concrete/norm*1000,holes/norm*1000,(4*exit+concrete+holes)/(6*norm)*1000))
835
836
837
838def drawNeutronGen(hname='start'):
839 ROOT.gStyle.SetPalette(ROOT.kRainBow)
840 pal = ROOT.TColor.GetPalette()
841 geoManager = ROOT.gGeoManager
842 geoManager.SetNsegments(12)
843 h['material'] = ROOT.TGeoMaterial("dummy")
844 h['medium'] = ROOT.TGeoMedium('dummy',1,h['material'])
845 h['top'] = geoManager.GetTopVolume()
846 top = h['top']
847 floor=geoManager.FindVolumeFast('vfloor')
848 floor.SetLineColor(ROOT.kGray)
849 box = geoManager.FindVolumeFast('vbox')
850 box.SetLineColor(ROOT.kBlue-9)
851 top.Draw('ogl')
852 nx,ny,nz = h[hname].GetXaxis().GetNbins(),h[hname].GetYaxis().GetNbins(),h[hname].GetZaxis().GetNbins()
853 hmax = pal.GetSize()/h[hname].GetMaximum()
854 h['x'] = geoManager.MakeBox('x',h['medium'] ,0.5,0.5,0.5)
855 h['x'].SetTransparency(50)
856 h['x'].SetLineColor(ROOT.kMagenta-6)
857 h['A'] = ROOT.TGeoVolumeAssembly("A")
858 n=0
859 h['T'] = []
860 for ix in range(1,nx+1,2):
861 for iy in range(1,ny+1,2):
862 for iz in range(1,nz+1,2):
863 C = h[hname].GetBinContent(ix,iy,iz)
864 if not C>0: continue
865 n+=1
866# h['start'].Fill(start.Z(),start.X(),start.Y(),W)
867 z,x,y = h[hname].GetXaxis().GetBinCenter(ix),h[hname].GetYaxis().GetBinCenter(iy),h[hname].GetZaxis().GetBinCenter(iz)
868 # kc = int(pal.GetAt(int(C*hmax)))
869 # print("color",int(C*hmax),kc)
870 h['T'].append(ROOT.TGeoTranslation(x,y, z))
871 h['A'].AddNode(h['x'], ix+1000*iy+1000000*iz ,h['T'][len(h['T'])-1])
872 s = str(ix+1000*iy+1000000*iz)
873 h['A'].AddNode(top,0)
874 h['A'].Draw('ogl')
875
876if options.command=='coldBox':
877 coldBox(plotOnly=False)
878if options.command=='coldBoxPlots':
879 coldBox()
880if options.command=='coldBoxPass1':
881 coldBox(pas='_pass1')
882
ana_thermalNeutrons.count
count(hFile)
Definition ana_thermalNeutrons.py:21
ana_thermalNeutrons.beamGas
beamGas()
Definition ana_thermalNeutrons.py:65
ana_thermalNeutrons.absorptionLength
absorptionLength(plotOnly=True)
Definition ana_thermalNeutrons.py:88
ana_thermalNeutrons.coldBox
coldBox(plotOnly=True, pas='')
Definition ana_thermalNeutrons.py:429
ana_thermalNeutrons.myPrint
myPrint(tc, tname)
Definition ana_thermalNeutrons.py:82
ana_thermalNeutrons.absorptionLengthOLD
absorptionLengthOLD()
Definition ana_thermalNeutrons.py:354
ana_thermalNeutrons.int
int
Definition ana_thermalNeutrons.py:10
ana_thermalNeutrons.float
float
Definition ana_thermalNeutrons.py:12
ana_thermalNeutrons.drawNeutronGen
drawNeutronGen(hname='start')
Definition ana_thermalNeutrons.py:838
ana_thermalNeutrons.flukaRateIntegrated
flukaRateIntegrated(save=False)
Definition ana_thermalNeutrons.py:394
ana_thermalNeutrons.reactions
reactions(hFile, distance=1E10)
Definition ana_thermalNeutrons.py:307
ana_thermalNeutrons.printReactions
printReactions()
Definition ana_thermalNeutrons.py:301
ana_thermalNeutrons.str
str
Definition ana_thermalNeutrons.py:11