SND@LHC Software
Loading...
Searching...
No Matches
Functions | Variables
study_muEloss Namespace Reference

Functions

 run ()
 
 makePlot (f, book=True)
 
 readChain ()
 
 NA62 ()
 
 makeSummaryPlot ()
 

Variables

str mcEngine = "TGeant4"
 
int runnr = 1
 
int nev = 1000000
 
dict setup = {}
 
 s = sys.argv[1]
 
dict thickness = setup[s]['thickness']
 
dict material = setup[s]['material']
 
dict momentum = setup[s]['momentum']
 
dict maxTheta = setup[s]['maxTheta']
 
bool checkOverlap = True
 
bool storeOnlyMuons = True
 
str outFile = "msc"+s+".root"
 
int theSeed = 0
 
float ecut = 0.0
 
dict h = {}
 

Function Documentation

◆ makePlot()

study_muEloss.makePlot (   f,
  book = True 
)

Definition at line 116 of file study_muEloss.py.

116def makePlot(f,book=True):
117# print interaction and radiation length of target
118 sGeo=ROOT.gGeoManager
119 if sGeo:
120 v = sGeo.FindVolumeFast('target')
121 m = v.GetMaterial()
122 length = v.GetShape().GetDZ()*2
123 print("Material:",m.GetName(),'total interaction length=',length/m.GetIntLen(),'total rad length=',length/m.GetRadLen())
124 else:
125 density= 2.413
126 length= 125.0
127 print("Use predefined values:",density,length)
128 if book:
129 ut.bookHist(h,'theta','scattering angle '+str(momentum)+'GeV/c;{Theta}(rad)',500,0,maxTheta)
130 ut.bookHist(h,'eloss','rel energy loss as function of momentum GeV/c',100,0,maxTheta,10000,0.,1.)
131 ut.bookHist(h,'elossRaw','energy loss as function of momentum GeV/c',100,0,maxTheta, 10000,0.,100.)
132 sTree = f.cbmsim
133 for n in range(sTree.GetEntries()):
134 rc = sTree.GetEvent(n)
135 Ein = sTree.MCTrack[0].GetEnergy()
136 M = sTree.MCTrack[0].GetMass()
137 Eloss = 0
138 for aHit in sTree.vetoPoint:
139 Eloss+=aHit.GetEnergyLoss()
140 print(Ein,Eloss/Ein)
141 rc = h['eloss'].Fill(Ein,Eloss/Ein)
142 rc = h['elossRaw'].Fill(Ein,Eloss)
143 ut.bookCanvas(h,key=s,title=s,nx=900,ny=600,cx=1,cy=1)
144 tc = h[s].cd(1)
145 if s=="NA62":
146 h['eloss'].Draw()
147 h['95'] = h['eloss'].ProjectionX('95',96,100)
148 h['95'].Sumw2()
149 h['0'] = h['eloss'].ProjectionX('0',1,100)
150 h['0'].Sumw2()
151 rc = h['95'].Divide(h['0'] )
152 h['95'].Draw()
153 h['meanEloss'] = h['elossRaw'].ProjectionX()
154 for n in range(1,h['elossRaw'].GetNbinsX()+1):
155 tmp = h['elossRaw'].ProjectionY('tmp',n,n)
156 eloss = tmp.GetMean()
157 h['meanEloss'].SetBinContent(n,eloss/density/length*1000)
158 h['meanEloss'].SetTitle('mean energy loss MeV cm2 / g')
159 h['meanEloss'].Draw()
160 elif s=="ATLAS":
161 h['eloss'].Draw()
162 h['>eloss']=h['eloss'].ProjectionY().Clone('>eloss')
163 cum = 0
164 N = float(h['>eloss'].GetEntries())
165 for n in range(h['>eloss'].GetNbinsX(),0,-1):
166 cum+=h['>eloss'].GetBinContent(n)
167 h['>eloss'].SetBinContent(n,cum/N)
168 print("Ethreshold event fraction in %")
169 for E in [15.,20.,30.,50.,80.]:
170 n = h['>eloss'].FindBin(E/350.)
171 print(" %5.0F %5.2F "%(E,h['>eloss'].GetBinContent(n)*100))
172 else:
173 tc.SetLogy(1)
174 h['theta_100']=h['theta'].Clone('theta_100')
175 h['theta_100']=h['theta'].Rebin(5)
176 h['theta_100'].Scale(1./h['theta_100'].GetMaximum())
177 h['theta_100'].Draw()
178 h[s].Print(s+'.png')
179 h[s].Print(s+'.root')
180 f.Write(h['theta'].GetName())
181 f.Write(h['theta_100'].GetName())
182

◆ makeSummaryPlot()

study_muEloss.makeSummaryPlot ( )

Definition at line 209 of file study_muEloss.py.

209def makeSummaryPlot():
210# using data in /mnt/hgfs/microDisk/Data/eloss/eloss_sum.root
211# krypton total interaction length= 1.97246306079 total rad length= 26.5231000393
212 pdg={10.0:1.914,14.0:1.978,20.0:2.055,30.0:2.164,40.0:2.263,80.0:2.630,100.:2.810,140.:3.170,200.:3.720,277.:4.420,300.:4.631,400.:5.561}
213 h['Gpdg'] = ROOT.TGraph(len(pdg))
214 Gpdg = h['Gpdg']
215 Gpdg.SetMarkerColor(ROOT.kRed)
216 Gpdg.SetMarkerStyle(20)
217 keys = sorted(pdg.keys())
218 for n in range(len(keys)):
219 Gpdg.SetPoint(n,keys[n],pdg[keys[n]])
220 density= 2.413
221 length= 125.0
222 ut.readHists(h,"/mnt/hgfs/microDisk/Data/eloss/eloss_sum.root")
223 ut.readHists(h,"/mnt/hgfs/microDisk/Data/eloss/eloss_withRaw.root")
224 ut.bookCanvas(h,key='summary',title=" ",nx=1200,ny=600,cx=2,cy=1)
225 tc = h['summary'].cd(1)
226 h['0'] = h['eloss'].ProjectionX('0',1,h['eloss'].GetNbinsY())
227 h['0'].Sumw2()
228 NA62()
229 for t in [93,95]:
230 h[t] = h['eloss'].ProjectionX(str(t),int(h['eloss'].GetNbinsY()*t/100.),h['eloss'].GetNbinsY())
231 h[t].Sumw2()
232 h[t].SetStats(0)
233 h[t].SetMarkerStyle(24)
234 rc = h[t].Divide(h['0'] )
235 h[t].Rebin(2)
236 h[t].Scale(1./2.)
237 if t!=93:
238 h[t].SetMarkerColor(ROOT.kBlue)
239 h[t].Draw('same')
240 else:
241 h[t].SetMaximum(1E-5)
242 h[t].SetMarkerColor(ROOT.kMagenta)
243 h[t].SetXTitle('incoming muon momentum [GeV/c]')
244 h[t].SetYTitle('prob #DeltaE>X%')
245 h[t].SetTitle('')
246 h[t].Draw()
247 h['NA62'].Draw('sameP')
248 h['lg'] = ROOT.TLegend(0.53,0.79,0.98,0.94)
249 h['lg'].AddEntry(h['NA62'],'NA62 measurement >95%','PL')
250 h['lg'].AddEntry(h[95],'FairShip >95%','PL')
251 h['lg'].AddEntry(h[93],'FairShip >93%','PL')
252 h['lg'].Draw()
253 tc = h['summary'].cd(2)
254 h['meanEloss'] = h['elossRaw'].ProjectionX()
255 for n in range(1,h['elossRaw'].GetNbinsX()+1):
256 tmp = h['elossRaw'].ProjectionY('tmp',n,n)
257 eloss = tmp.GetMean()
258 h['meanEloss'].SetBinContent(n,eloss/density/length*1000)
259 h['meanEloss'].SetBinError(n,0)
260 h['meanEloss'].SetTitle('mean energy loss MeV cm^{2}/g')
261 h['meanEloss'].SetStats(0)
262 h['meanEloss'].SetMaximum(7.)
263 h['meanEloss'].SetXTitle('incoming muon momentum [GeV/c]')
264 h['meanEloss'].SetYTitle('mean energy loss [MeV cm^[2]]/g')
265 h['meanEloss'].SetTitle('')
266 h['meanEloss'].Draw()
267 Gpdg.Draw('sameP')
268 h['lg2'] = ROOT.TLegend(0.53,0.79,0.98,0.94)
269 h['lg2'].AddEntry(h['Gpdg'],'muon dE/dx, PDG ','PL')
270 h['lg2'].AddEntry(h['meanEloss'],'energy deposited in krypton, FairShip','PL')
271 h['lg2'].Draw()
272 h['summary'].Print('catastrophicEnergyLoss.png')

◆ NA62()

study_muEloss.NA62 ( )

Definition at line 190 of file study_muEloss.py.

190def NA62():
191 na62Points = open('NA62.points')
192 allPoints = na62Points.readlines()
193 N = int((len(allPoints)-1)/3.)
194 h['NA62']=ROOT.TGraphErrors(N)
195 for l in range(N):
196 tmp = allPoints[3*l].split(',')
197 x=float(tmp[0])
198 y=float(tmp[1].replace('\n',''))
199 tmp = allPoints[3*l+1].split(',')
200 y1=float(tmp[1].replace('\n',''))
201 tmp = allPoints[3*l+2].split(',')
202 y2=float(tmp[1].replace('\n',''))
203 h['NA62'].SetPoint(l,x,y*1E-6)
204 h['NA62'].SetPointError(l,0,abs(y1-y2)/2.*1E-6)
205 h['NA62'].SetLineColor(ROOT.kRed)
206 h['NA62'].SetMarkerColor(ROOT.kRed)
207 h['NA62'].SetMarkerStyle(20)
208

◆ readChain()

study_muEloss.readChain ( )

Definition at line 183 of file study_muEloss.py.

183def readChain():
184 tmp = "/mnt/hgfs/microDisk/Data/mscNA62_X.root"
185 for i in [0,1]:
186 f = ROOT.TFile(tmp.replace('X',str(i)))
187 if i==1: makePlot(f,False)
188 else: makePlot(f)
189

◆ run()

study_muEloss.run ( )

Definition at line 48 of file study_muEloss.py.

48def run():
49# -------------------------------------------------------------------
50 ROOT.gRandom.SetSeed(theSeed) # this should be propagated via ROOT to Pythia8 and Geant4VMC
51 shipRoot_conf.configure() # load basic libraries, prepare atexit for python
52 # ship_geo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = 10, tankDesign = 5, muShieldDesign = 7, nuTauTargetDesign=1)
53# -----Timer--------------------------------------------------------
54 timer = ROOT.TStopwatch()
55 timer.Start()
56# -----Create simulation run----------------------------------------
57 gFairBaseContFact = ROOT.FairBaseContFact() # required by change to FairBaseContFact to avoid TList::Clear errors
58 run = ROOT.FairRunSim()
59 run.SetName(mcEngine) # Transport engine
60 if nev==0: run.SetOutputFile("dummy.root")
61 else: run.SetOutputFile(outFile) # Output file
62 run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C
63 rtdb = run.GetRuntimeDb()
64# -----Materials----------------------------------------------
65 run.SetMaterials("media.geo")
66# -----Create geometry----------------------------------------------
67 cave= ROOT.ShipCave("CAVE")
68 cave.SetGeometryFileName("cave.geo")
69 run.AddModule(cave)
70#
71 target = ROOT.simpleTarget()
72 material, thickness, 0
73#
74 target.SetEnergyCut(ecut*u.GeV)
75 if storeOnlyMuons: target.SetOnlyMuons()
76 target.SetParameters(material,thickness,0.)
77 run.AddModule(target)
78#
79 primGen = ROOT.FairPrimaryGenerator()
80 myPgun = ROOT.FairBoxGenerator(13,1) # pdg id and multiplicity
81 if s=="NA62": myPgun.SetPRange(momentum,maxTheta)
82 else: myPgun.SetPRange(momentum-0.01,momentum+0.01)
83 myPgun.SetPhiRange(0,0) # // Azimuth angle range [degree]
84 myPgun.SetThetaRange(0,0) # // Polar angle in lab system range [degree]
85 myPgun.SetXYZ(0.*u.cm, 0.*u.cm, -1.*u.mm - (thickness) )
86 primGen.AddGenerator(myPgun)
87#
88 run.SetGenerator(primGen)
89# -----Initialize simulation run------------------------------------
90 run.Init()
91 if nev==0: return
92 gMC = ROOT.TVirtualMC.GetMC()
93
94 fStack = gMC.GetStack()
95 fStack.SetMinPoints(1)
96 fStack.SetEnergyCut(-1.)
97
98# -----Start run----------------------------------------------------
99 print("run for ",nev,"events")
100 run.Run(nev)
101
102# -----Start Analysis---------------
103 ROOT.gROOT.ProcessLine('#include "Geant4/G4EmParameters.hh"')
104 emP = ROOT.G4EmParameters.Instance()
105 emP.Dump()
106 h['f']= ROOT.gROOT.GetListOfFiles()[0].GetName()
107# -----Finish-------------------------------------------------------
108 timer.Stop()
109 rtime = timer.RealTime()
110 ctime = timer.CpuTime()
111 print(' ')
112 print("Macro finished succesfully.")
113 print("Output file is ", outFile)
114 print("Real time ",rtime, " s, CPU time ",ctime,"s")
115
shipRoot_conf.configure
configure(darkphoton=None)
Definition shipRoot_conf.py:67

Variable Documentation

◆ checkOverlap

bool study_muEloss.checkOverlap = True

Definition at line 38 of file study_muEloss.py.

◆ ecut

float study_muEloss.ecut = 0.0

Definition at line 43 of file study_muEloss.py.

◆ h

dict study_muEloss.h = {}

Definition at line 46 of file study_muEloss.py.

◆ material

dict study_muEloss.material = setup[s]['material']

Definition at line 34 of file study_muEloss.py.

◆ maxTheta

dict study_muEloss.maxTheta = setup[s]['maxTheta']

Definition at line 36 of file study_muEloss.py.

◆ mcEngine

str study_muEloss.mcEngine = "TGeant4"

Definition at line 10 of file study_muEloss.py.

◆ momentum

dict study_muEloss.momentum = setup[s]['momentum']

Definition at line 35 of file study_muEloss.py.

◆ nev

int study_muEloss.nev = 1000000

Definition at line 12 of file study_muEloss.py.

◆ outFile

str study_muEloss.outFile = "msc"+s+".root"

Definition at line 41 of file study_muEloss.py.

◆ runnr

int study_muEloss.runnr = 1

Definition at line 11 of file study_muEloss.py.

◆ s

study_muEloss.s = sys.argv[1]

Definition at line 32 of file study_muEloss.py.

◆ setup

dict study_muEloss.setup = {}

Definition at line 14 of file study_muEloss.py.

◆ storeOnlyMuons

bool study_muEloss.storeOnlyMuons = True

Definition at line 39 of file study_muEloss.py.

◆ theSeed

int study_muEloss.theSeed = 0

Definition at line 42 of file study_muEloss.py.

◆ thickness

dict study_muEloss.thickness = setup[s]['thickness']

Definition at line 33 of file study_muEloss.py.