ParticleGun_FastCalo_ChargeFlip_Config.py 2.85 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#! -*- python -*-
evgenConfig.description = "Single particle gun for FastCaloSim event generation"
evgenConfig.keywords = ["singleParticle",]
evgenConfig.generators = ["ParticleGun"]
evgenConfig.contact = ["david.sosa@cern.ch"]

import ParticleGun as PG
import ROOT

class MyParticleSampler(PG.ParticleSampler):
    def __init__(self,energy,eta,pid,shift_z=0):
        self.pid = pid
        self.shift_z = shift_z
        pdg_table = ROOT.TDatabasePDG.Instance()
        mass = pdg_table.GetParticle(self.pid()).Mass()
        self.mom1 = PG.EEtaMPhiSampler(energy=energy,eta=eta,mass=mass)

    def shoot(self):
        pid = self.pid()
        
        shift_z = self.shift_z

        mom = self.mom1.shoot()
        pos_temp = mom.Vect().Unit()
        
        # Would it hit the barrel, or the endcap?
        if abs(pos_temp.Z())/3550.<pos_temp.Perp()/1148.: # Hit the barrel!
            pos_temp *= 1148./pos_temp.Perp()
        else: # Hit the endcap!
            pos_temp *= 3550./abs(pos_temp.Z())

        # Shift position of vector in the Z direction
        pos_temp_2 = ROOT.TVector3()
        pos_temp_2.SetXYZ(pos_temp.X(), pos_temp.Y(), pos_temp.Z()+shift_z)
        pos_temp_2 *= 1. / pos_temp_2.Mag(); # reduce magnitude of vector

        # recalculate; Would it hit the barrel, or the endcap?
        if abs(pos_temp_2.Z())/3550.<pos_temp_2.Perp()/1148.:
            pos_temp_2 *= 1148./pos_temp_2.Perp()
        else:
            pos_temp_2 *= 3550./abs(pos_temp_2.Z())

        pos = ROOT.TLorentzVector(pos_temp_2.X(),pos_temp_2.Y(),pos_temp_2.Z(), pos_temp_2.Mag())
        
        #print "pid ",pid
        
        return [ PG.SampledParticle( pid , mom , pos ) ]

49
50
myE = float(jofile.split('_E')[1].split('_')[0])
myZV = float(jofile.split('_')[-1].split('.py')[0].replace("m","-"))
51

52
myPDGID = jofile.split('_pid')[1].split('_')[0].replace('n','-')
53
54
55
56
myPDGID = int(float(myPDGID.replace('p','')))

eta_li = []

57
58
59
60
61
if "disj" in jofile:
    myLowEta1  = 0.01*float(jofile.split('eta_')[1].split('_')[0].replace('m','-'))
    myLowEta2  = 0.01*float(jofile.split('eta_')[1].split('_')[1].replace('m','-'))
    myHighEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[2].replace('m','-'))
    myHighEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[3].replace('m','-'))
62
63
64
    eta_li.extend([myLowEta1,myLowEta2,myHighEta1,myHighEta2])

else:
65
66
    myLowEta  = 0.01*float(jofile.split('eta')[1].split('_')[0].replace('m','-'))
    myHighEta = 0.01*float(jofile.split('eta')[1].split('_')[1].replace('m','-'))
67
68
69
70
71
72
73
74
75
76
77
78
    eta_li.extend([myLowEta,myHighEta])


print "================ SETTTINGS ================="
print ("energy = ", myE)
print ("eta = ", eta_li)
print ("pid = ", myPDGID)
print ("shift_z = ", myZV)
print "============================================"

genSeq += PG.ParticleGun()
genSeq.ParticleGun.sampler = MyParticleSampler(energy=myE,eta=eta_li,pid=(myPDGID,myPDGID),shift_z=myZV) #unmixed