SND@LHC Software
Loading...
Searching...
No Matches
makeALPACAEvents.py
Go to the documentation of this file.
1import os,sys,getopt
2import numpy,math
3import ROOT
4
5hGeV = 6.58211928*pow(10.,-16)* pow(10.,-9) # no units or it messes up!!
6c_light = 2.99792458e+10
7mres="1.0d0"
8gax="1d-7"
9nev="100"
10
11def ALPACAFormatting(s):
12 s=str(s)
13 if s.find('d')==-1:
14 s = numpy.format_float_scientific(float(s))
15 s= s.replace('e','d')
16 s= s.replace('+','')
17 return s
18
19
20def Ctau(mres,gax):
21 return c_light*hGeV*65.*math.pi/((float(gax)*float(gax))*(float(mres)*float(mres)*float(mres)))
22
23def Decaylength(e,p,ctau):
24 beta=p/e
25 gamma=e/math.sqrt(e*e-p*p)
26 return beta*gamma*ctau
27
28def Decayweight(Lmin,Lmax,Decaylength,LS):
29 return math.exp(-LS/Decaylength)*((Lmax-Lmin)/Decaylength)
30
31def inputWrite(mres,gax,nev,Lmin,Lmax):#need to apply Rmin & Rmax as well..
32 Mres = ALPACAFormatting(mres)
33 Gax = ALPACAFormatting(gax)
34 L = ALPACAFormatting(float(Lmax)-float(Lmin))
35 Lmin = ALPACAFormatting(Lmin)
36 Lmax = ALPACAFormatting(Lmax)
37 with open("input.DAT","w") as f:
38 f.write("****************************************************************************************\n")
39 f.write("*********** RE-RUN ./init IF FIRST FIVE PARAMETERS ARE CHANGED: **********************\n")
40 f.write("****************************************************************************************\n")
41 f.write("************************* Experimental parameters *************************************\n")
42 f.write("****************************************************************************************\n")
43 f.write("4d2 ! [ebeam] : Beam kinetic energy (GeV)\n")
44 f.write("'prot' ! [btype] : Beam type ('prot' = proton, 'elec' = electron)\n")
45 f.write("95.95d0 ! [aa] : Target mass number\n")
46 f.write("42d0 ! [az] : Target atomic number\n")
47 f.write("%s ! [lsh] : Shielding length (m)\n"%(Lmin))
48 f.write("%s ! [dvol] : Decay volume (m)\n"%(L))
49 f.write("1.5d-1 ! [rmin] : Inner radius (m)\n")
50 f.write("2.5d0 ! [rmax] : Outer radius (m)\n")
51 f.write("0.1d0 ! [dmin] : Minimum photon separation (m)\n")
52 f.write("1d0 ! [emin] : Minimum photon energy (GeV)\n")
53 f.write("1d-2 ! [athetamax] : Maximum ALP theta (when adecay is false)\n")
54 f.write("****************************************************************************************\n")
55 f.write("********************************** ALP parameters **************************************\n")
56 f.write("****************************************************************************************\n")
57 f.write("%s ! [mres] : ALP mass (GeV)\n"%(Mres))
58 f.write("%s ! [gax] : ALP photon coupling (GeV^-1)\n"%(Gax))
59 f.write("****************************************************************************************\n")
60 f.write("****************************** Output parameters **************************************\n")
61 f.write("****************************************************************************************\n")
62 f.write("'_%s_%s' ! [outtag] : For output file\n"%(mres,gax))
63 f.write("****************************************************************************************\n")
64 f.write("************************* Integration parameters ***************************************\n")
65 f.write("****************************************************************************************\n")
66 f.write("100000 ! [ncall] : Number of calls for preconditioning\n")
67 f.write("10 ! [itmx] : Number of iterations for preconditioning\n")
68 f.write("1d0 ! [prec] : Relative accuracy (in %) in main run\n")
69 f.write("100000 ! [ncall1] : Number of calls in first iteration\n")
70 f.write("100000 ! [inccall] : Number of increase calls per iteration\n")
71 f.write("50 ! [itend] : Maximum number of iterations\n")
72 f.write("6 ! [iseed] : Random number seed (integer > 0)\n")
73 f.write("****************************************************************************************\n")
74 f.write("******************************* Unweighted events **************************************\n")
75 f.write("****************************************************************************************\n")
76 f.write(".true. ! [genunw] : Generate unweighted events\n")
77 f.write("%s ! [nev] : Number of events ( < 1000000 recommended)\n"%(nev))
78 f.write("'hepevt' ! [erec] : Event record format ('lhe' = Les Houches, 'hepevt' = HEPEVT)\n")
79 f.write("****************************************************************************************\n")
80 f.write("******************************* Cuts ************************************************\n")
81 f.write(".true. ! [gencuts] : Generate cuts - [rmin], [rmax], [dmin], [emin]\n")
82 f.write(".true. ! [adecay] : Include ALP decay\n")
83 f.write("****************************************************************************************\n")
84 f.write("****************************************************************************************\n")
85 f.write("****************************************************************************************\n")
86
87def ntupleWrite(ctau,fn,Lmin,Lmax,startZ,endZ,SmearBeam,mres,gax,old):
88 Lmin = float(Lmin)
89 Lmax = float(Lmax)
90 startZ = float(startZ)
91 endZ = float(endZ)
92 SmearBeam = float(SmearBeam)
93 wdir = os.environ['ALPACABIN']
94 with open(wdir+"/outputs/output_"+str(mres)+"_"+str(gax)+".dat","r") as fp:
95 for line in fp:
96 if "Cross section" in line:
97 line=line.split()
98 xs=line[3]+"_"+line[5]
99 fn=fn.readlines()
100 inputFile=old+"/"+"alp_m"+str(mres)+"_g"+str(gax)+"_xs"+str(xs)+".root"
101 f=ROOT.TFile(inputFile,"recreate")
102 ntup=ROOT.TNtuple("MCTrack","Track Informations","event:track:pdg:px:py:pz:x:y:z:parent:decay:e:tof:w")
103 for i,j in enumerate(range(5,len(fn),9)):
104 LS = ROOT.gRandom.Uniform(Lmin*100.,Lmax*100.)
105 zinter = ROOT.gRandom.Uniform(startZ,endZ)
106 dx, dy = ROOT.gRandom.Uniform(-1,+1)*SmearBeam, ROOT.gRandom.Uniform(-1,+1)*SmearBeam
107 tr=fn[j].split()
108 dau1=fn[j+1].split()
109 dau2=fn[j+2].split()
110 px,py,pz=float(tr[7]),float(tr[8]),float(tr[9])
111 p = math.sqrt(px**2.+py**2.+pz**2.)
112 lam = LS/p
113 daux,dauy,dauz= dx+lam*px,dy+lam*py,zinter+lam*pz
114 dl=Decaylength(float(tr[10]),p,ctau)
115 w = Decayweight(Lmin*100.,Lmax*100.,dl,LS)
116 ntup.Fill(int(i),int(0),int(9900015),px,py,pz,dx,dy,zinter,int(-1),float(0),float(tr[10]),float(0),w)
117 #ntup.Fill(int(i),int(1),int(dau1[1]),float(dau1[7]),float(dau1[8]),float(dau1[9]),float(dau1[12])/10.+dx,float(dau1[13])/10.+dy,float(dau1[14])/10.+zinter,int(0),float(1),float(dau1[10]),float(dau1[15])/10./c_light,w)
118 #ntup.Fill(int(i),int(2),int(dau2[1]),float(dau2[7]),float(dau2[8]),float(dau2[9]),float(dau2[12])/10.+dx,float(dau2[13])/10.+dy,float(dau2[14])/10.+zinter,int(0),float(1),float(dau2[10]),float(dau2[15])/10./c_light,w)
119 ntup.Fill(int(i),int(1),int(dau1[1]),float(dau1[7]),float(dau1[8]),float(dau1[9]),daux,dauy,dauz,int(0),float(1),float(dau1[10]),float(dau1[15])/10./c_light,w)
120 ntup.Fill(int(i),int(2),int(dau2[1]),float(dau2[7]),float(dau2[8]),float(dau2[9]),daux,dauy,dauz,int(0),float(1),float(dau2[10]),float(dau2[15])/10./c_light,w)
121 ntup.Write()
122 f.Close()
123 return inputFile
124
125def runEvents(mres,gax,nev,Lmin,Lmax,startZ,endZ,SmearBeam):
126 print('ALPACA is starting for mass of {} GeV with photon coupling coeffiecient {} GeV^-1.'.format(mres,gax))
127 pdg = ROOT.TDatabasePDG.Instance()
128 pdg.AddParticle('a','ALP', mres, False, gax, 0., 'a', 9900015)
129 wdir = os.environ['ALPACABIN']
130 old=os.getcwd()
131 os.chdir(wdir)
132 inputWrite(mres,gax,nev,Lmin,Lmax)
133 rn="./alpaca < input.DAT"
134 os.system(rn)
135 print('ALPACA generated the events.')
136 pa="./evrecs/evrec_"+str(mres)+"_"+str(gax)+".dat"
137 fn=open(pa,"r")
138 ctau=Ctau(mres,gax)
139 print('Events are recording into a ntuple.')
140 inputFile = ntupleWrite(ctau,fn,Lmin,Lmax,startZ,endZ,SmearBeam,mres,gax,old)
141 os.chdir(old)
142 print('Ntuple is ready for the reading.')
143 return inputFile
144#runEvents()
makeALPACAEvents.Decayweight
Decayweight(Lmin, Lmax, Decaylength, LS)
Definition makeALPACAEvents.py:28
makeALPACAEvents.inputWrite
inputWrite(mres, gax, nev, Lmin, Lmax)
Definition makeALPACAEvents.py:31
makeALPACAEvents.Decaylength
Decaylength(e, p, ctau)
Definition makeALPACAEvents.py:23
makeALPACAEvents.ntupleWrite
ntupleWrite(ctau, fn, Lmin, Lmax, startZ, endZ, SmearBeam, mres, gax, old)
Definition makeALPACAEvents.py:87
makeALPACAEvents.runEvents
runEvents(mres, gax, nev, Lmin, Lmax, startZ, endZ, SmearBeam)
Definition makeALPACAEvents.py:125