SND@LHC Software
Loading...
Searching...
No Matches
runPythia8.py
Go to the documentation of this file.
1import ROOT
2import rootUtils as ut
3from array import array
4theSeed = 0
5h = {}
6ROOT.gRandom.SetSeed(theSeed)
7ROOT.gSystem.Load("libpythia8")
8
9from argparse import ArgumentParser
10parser = ArgumentParser()
11parser.add_argument("-b", "--heartbeat", dest="heartbeat", type=int, help="progress report", default=10000)
12parser.add_argument("-n", "--pot", dest="NPoT", type=int, help="protons on target", default=1000000)
13parser.add_argument("-r", "--run", dest="run", type=int, help="production sequence number", default=0)
14parser.add_argument('-M', '--Emin', dest='Emin', type=float,help="cutOff", default=1.0)
15parser.add_argument('-E', '--Beam', dest='fMom', type=float,help="beam momentum", default=400.)
16parser.add_argument('-J', '--Jpsi-mainly', action='store_true', dest='JpsiMainly', default=False)
17parser.add_argument('-Y', '--DrellYan', action='store_true', dest='DrellYan', default=False)
18parser.add_argument('-P', '--PhotonCollision', action='store_true', dest='PhotonCollision', default=False)
19parser.add_argument('-C', '--charm', action='store_true', dest='charm', default=False)
20parser.add_argument('-X', '--PDFpSet',dest="PDFpSet", type=str, help="PDF pSet to use", default="13")
21parser.add_argument('-u', '--uOnly', action='store_true', dest='uOnly', default=False)
22
23# for lhapdf, -X LHAPDF6:cteq6
24
25options = parser.parse_args()
26print("start IGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNORE")
27X=ROOT.FixedTargetGenerator()
28print("end IGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNOREIGNORE")
29
30def yBeam(Mproton = 0.938272081, pbeam = 400.):
31 Ebeam = ROOT.TMath.Sqrt(pbeam**2+Mproton**2)
32 betaCM = pbeam / (Ebeam + Mproton)
33 y_beam = 0.5*ROOT.TMath.Log( (1+betaCM)/(1-betaCM)) # https://arxiv.org/pdf/1604.02651.pdf
34 return y_beam
35
36generators = {'p':ROOT.Pythia8.Pythia(),'n':ROOT.Pythia8.Pythia()}
37generators['p'].settings.mode("Beams:idB", 2212)
38generators['n'].settings.mode("Beams:idB", 2112)
39
40for g in generators:
41 ut.bookHist(h, 'xsec_'+g, ' total cross section',1,0.,1.)
42 ut.bookHist(h, 'M_'+g, ' N mu+mu-;M [GeV/c^{2}];y_{CM}',500,0.,10.,120,-3.,3.,100,0.,5.)
43 ut.bookHist(h, 'cosCS_'+g, ' N cosCS;cosCS;y_{CM}',100,-1.,1.,120,-3.,3.,100,0.,5.)
44 ut.bookHist(h, 'cosCSJpsi_'+g, ' N cosCS 2.9<M<4.5;cosCS;y_{CM}',100,-1.,1.,120,-3.,3.,100,0.,5.)
45 generators[g].settings.mode("Next:numberCount",options.heartbeat)
46 generators[g].settings.mode("Beams:idA", 2212)
47 generators[g].settings.mode("Beams:frameType", 2)
48 generators[g].settings.parm("Beams:eA",options.fMom)
49 generators[g].settings.parm("Beams:eB",0.)
50 generators[g].readString("PDF:pSet = "+options.PDFpSet)
51 if options.DrellYan:
52 generators[g].readString("WeakSingleBoson:ffbar2gmZ = on")
53 if options.uOnly:
54 generators[g].readString("23:onMode = off")
55 generators[g].readString("23:onIfAll = 13 13")
56 generators[g].readString("23:mMin = "+str(options.Emin))
57 generators[g].readString("PhaseSpace:mHatMin = "+str(options.Emin))
58 # generators[g].readString("BeamRemnants:primordialKThard = 0.8")
59 generators[g].readString("PartonLevel:FSR = on")
60 elif options.JpsiMainly:
61# use this for all onia productions
62 generators[g].readString("Onia:all(3S1) = on")
63 if options.uOnly:
64 generators[g].readString("443:onMode = off")
65 generators[g].readString("443:onIfAll = 13 13")
66 generators[g].readString("553:onMode = off")
67 generators[g].readString("553:onIfAll = 13 13")
68 elif options.PhotonCollision:
69 generators[g].readString("PhotonCollision:gmgm2mumu = on")
70 generators[g].readString("PhaseSpace:mHatMin = "+str(options.Emin))
71 elif options.charm:
72 generators[g].readString("HardQCD:hardccbar = on")
73 else:
74 generators[g].readString("SoftQCD:inelastic = on")
75 generators[g].init()
76
77rc = generators['p'].next()
78processes = generators['p'].info.codesHard()
79hname = 'pythia8_PDFpset'+options.PDFpSet+'_Emin'+str(options.Emin)+'_'+generators['p'].info.nameProc(processes[0])
80hname = hname.replace('*','star')
81hname = hname.replace('->','to')
82hname = hname.replace('/','')
83
84f = ROOT.TFile("ntuple-"+hname+".root","RECREATE")
85signal = ROOT.TNtuple("ntuple","ntuple","M:P:Pt:y:p1x:p1y:p1z:p2x:p2y:p2z:cosCS")
86
87timer = ROOT.TStopwatch()
88timer.Start()
89
90ntagged = {'p':0,'n':0}
91ybeam = yBeam(pbeam = options.fMom)
92for n in range(int(options.NPoT)):
93 for g in generators:
94 py = generators[g]
95 rc = py.next()
96 nmu = {}
97 for ii in range(1,py.event.size()):
98 if options.DrellYan and py.event[ii].id()!=23: continue
99 if options.PhotonCollision and py.event[ii].id()!=22: continue
100 for m in py.event.daughterList(ii):
101 if abs(py.event[m].id())==13: nmu[m]=ii
102 if len(nmu) == 2:
103 ntagged[g]+=1
104 ks = list(nmu)
105 if options.DrellYan:
106 Zstar = py.event[nmu[ks[0]]]
107 rc=h['M_'+g].Fill(Zstar.m(),py.event[nmu[ks[0]]].y()-ybeam,py.event[nmu[ks[0]]].pT())
108# what about polarization?
109 ii = nmu[ks[0]]
110 d0 = py.event.daughterList(ii)[0]
111 d1 = py.event.daughterList(ii)[1]
112 if py.event[d0].id() < 0:
113 nlep = py.event[d0]
114 nantilep = py.event[d1]
115 else:
116 nlep = py.event[d1]
117 nantilep = py.event[d0]
118 P1pl = nlep.e()+nlep.pz()
119 P2pl = nantilep.e()+nantilep.pz()
120 P1mi = nlep.e()-nlep.pz()
121 P2mi = nantilep.e()-nantilep.pz()
122 A = P1pl*P2mi-P2pl*P1mi
123 cosCS = Zstar.pz()/abs(Zstar.pz()) * 1./Zstar.m()/ROOT.TMath.Sqrt(Zstar.m2()+Zstar.pT()**2)*A
124 rc=h['cosCS_'+g].Fill(cosCS,Zstar.y()-ybeam,py.event[nmu[ks[0]]].pT())
125 if Zstar.m()>2.9 and Zstar.m()<4.5: rc=h['cosCSJpsi_'+g].Fill(cosCS,py.event[nmu[ks[0]]].y()-ybeam,py.event[nmu[ks[0]]].pT())
126 rc = signal.Fill(Zstar.m(),Zstar.pAbs(),Zstar.pT(),Zstar.y(),nlep.px(),nlep.py(),nlep.pz(),nantilep.px(),nantilep.py(),nantilep.pz(),cosCS)
127 if options.PhotonCollision:
128 M={}
129 k=0
130 for m in ks:
131 M[k]=ROOT.TLorentzVector(py.event[m].px(),py.event[m].py(),py.event[m].pz(),py.event[m].e())
132 k+=1
133 G = M[0]+M[1]
134 rc=h['M_'+g].Fill(G.M(),G.Rapidity()-ybeam,G.Pt())
135
136
137print(">>>> proton statistics, ntagged=",ntagged['p'])
138generators['p'].stat()
139print( ">>>> neutron statistics, ntagged=",ntagged['n'])
140generators['n'].stat()
141
142timer.Stop()
143rtime = timer.RealTime()
144ctime = timer.CpuTime()
145print("run ",options.run," PoT ",options.NPoT," Real time ",rtime, " s, CPU time ",ctime,"s")
146signal.Write()
147for g in generators:
148 sigma = generators[g].info.sigmaGen(processes[0])
149 h['xsec_'+g].SetBinContent(1,sigma)
150ut.writeHists(h,hname+'.root')
151
152def na50(online=True):
153 for g in generators:
154 if online:
155 processes = generators[g].info.codesHard()
156 name = generators[g].info.nameProc(processes[0])
157 sigma = generators[g].info.sigmaGen(processes[0])
158 else:
159 name = ''
160 sigma = h['xsec_'+g].GetBinContent(1)
161 yax = h['M_'+g].GetYaxis()
162 xax = h['M_'+g].GetXaxis()
163 Mmin = xax.FindBin(2.9)
164 Mmax = xax.FindBin(4.5)
165 Ymin = yax.FindBin(-0.425)
166 Ymax = yax.FindBin(0.575)
167 h['MA'] = h['M_'+g].ProjectionX('MA')
168 h['M'] = h['M_'+g].ProjectionX('M',Ymin,Ymax)
169 print("generator sigma mumu-ratio in-mass-range in-y-range")
170 print("%s %s %6.2F nbarn, %5.2F, %5.2G, %5.2F "%(g,name,sigma*1E6,\
171 float(h['MA'].GetEntries())/options.NPoT,\
172 h['MA'].Integral(Mmin,Mmax)/float(h['MA'].GetEntries()),\
173 h['M'].GetEntries()/float(h['MA'].GetEntries())))
174 fraction = h['M'].Integral(Mmin,Mmax)/options.NPoT
175 # multiply with 0.5 assuming no polarization -0.5 < cosCS < 0.5
176 print("cross section a la NA50 for : %s %5.2F pb"%(g,0.5*fraction*sigma*1E9))
177 # via cosCS
178 for g in generators:
179 if online:
180 processes = generators[g].info.codesHard()
181 name = generators[g].info.nameProc(processes[0])
182 sigma = generators[g].info.sigmaGen(processes[0])
183 else:
184 name = ''
185 sigma = h['xsec_'+g].GetBinContent(1)
186 yax = h['cosCSJpsi_'+g].GetYaxis()
187 xax = h['cosCSJpsi_'+g].GetXaxis()
188 Mmin = xax.FindBin(-0.5)
189 Mmax = xax.FindBin(0.5)
190 Ymin = yax.FindBin(-0.425)
191 Ymax = yax.FindBin(0.575)
192 h['MA'] = h['cosCSJpsi_'+g].ProjectionX('MA')
193 h['M'] = h['cosCSJpsi_'+g].ProjectionX('M',Ymin,Ymax)
194 print("generator sigma mumu-in-mass-range% cosCS in-y-range")
195 print("%s %s %6.2F nbarn, %5.2F, %5.2F, %5.2F "%(g,name,sigma*1E6,\
196 float(h['MA'].GetEntries())/options.NPoT*100.,\
197 h['MA'].Integral(Mmin,Mmax)/float(h['MA'].GetEntries()),\
198 h['M'].GetEntries()/float(h['MA'].GetEntries())))
199 fraction = h['M'].Integral(Mmin,Mmax)/options.NPoT
200 # taking polarization into account.
201 print("cross section a la NA50, -0.5<cosCS<0.5: %5.2F pb"%(fraction*sigma*1E9))
202
203def muflux():
204 Z_Mo = 96.
205 P_Mo = 42
206 fraction = {}
207 for g in generators:
208 processes = generators[g].info.codesHard()
209 name = generators[g].info.nameProc(processes[0])
210 sigma = generators[g].info.sigmaGen(processes[0])
211 yax = h['M_'+g].GetYaxis()
212 xax = h['M_'+g].GetXaxis()
213 Mmin = xax.FindBin(0.)
214 Mmax = xax.FindBin(5.)
215 Ymin = yax.FindBin(0.3)
216 Ymax = yax.FindBin(3.)
217 h['MA'] = h['M_'+g].ProjectionX('MA')
218 h['M'] = h['M_'+g].ProjectionX('M',Ymin,Ymax)
219 print(g,name,sigma,float(h['MA'].GetEntries())/options.NPoT,h['MA'].Integral(Mmin,Mmax)/float(h['MA'].GetEntries()),h['M'].GetEntries()/float(h['MA'].GetEntries()))
220 fraction[g] = h['M'].Integral(Mmin,Mmax)/options.NPoT
221 meanFraction = (fraction['p']*P_Mo+fraction['n']*(Z_Mo-P_Mo))/Z_Mo * sigma
222 print("cross section a la muflux: %5.2F pb"%(0.5*meanFraction*1E9))
223
224def debugging(g):
225 generators[g].settings.listAll()
runPythia8.yBeam
yBeam(Mproton=0.938272081, pbeam=400.)
Definition runPythia8.py:30
runPythia8.na50
na50(online=True)
Definition runPythia8.py:152