1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | #include "TLorentzVector.h" |
14 | #include "TLorentzRotation.h" |
15 | #include "TRandom3.h" |
16 | |
17 | #include "particleType.h" |
18 | |
19 | #include "AMPTOOLS_MCGEN/GammaZToXZ.h" |
20 | #include "UTILITIES/BeamProperties.h" |
21 | |
22 | GammaZToXZ::GammaZToXZ( float massX, TString beamConfigFile, Double_t Bslope) : |
23 | m_target( 0, 0, 0, ParticleMass(Pb208)), |
24 | m_slope( Bslope ), |
25 | m_recoil (ParticleMass(Pb208)), |
26 | m_childMass( 0 ) |
27 | { |
28 | kMproton=ParticleMass(Proton); |
29 | kMneutron=ParticleMass(Neutron); |
30 | |
31 | |
32 | kMZ = ParticleMass(Pb208); |
33 | kMPion = ParticleMass(PiPlus); |
34 | kMPi0 = ParticleMass(Pi0); |
35 | kMKaon = ParticleMass(KPlus); |
36 | |
37 | m_childMass.push_back( massX ); |
38 | |
39 | |
40 | BeamProperties beamProp(beamConfigFile); |
41 | cobrem_vs_E = (TH1D*)beamProp.GetFlux(); |
42 | cobrem_vs_E->GetName(); |
43 | } |
44 | |
45 | Kinematics* |
46 | GammaZToXZ::generate(){ |
47 | |
48 | |
49 | const double kPi = 3.14159; |
50 | |
51 | |
52 | double EtaMass=m_childMass[0]; |
53 | double recoilMass = m_recoil; |
54 | |
55 | |
56 | double beamE = cobrem_vs_E->GetRandom(); |
57 | TLorentzVector beam; |
58 | beam.SetPxPyPzE(0,0,beamE,beamE); |
59 | TLorentzVector cm = beam + m_target; |
60 | |
61 | TLorentzRotation lab2cmBoost( -cm.BoostVector() ); |
62 | TLorentzRotation cm2labBoost( cm.BoostVector() ); |
63 | |
64 | TLorentzVector cmCM = lab2cmBoost * cm; |
65 | |
66 | double cmEnergy = ( lab2cmBoost * cm ).E(); |
67 | double beamMomCM = cmMomentum( cmEnergy, beam.M(), m_target.M() ); |
68 | double EtaMomCM = cmMomentum( cmEnergy, EtaMass, recoilMass ); |
69 | |
70 | |
71 | |
72 | |
73 | |
74 | |
75 | |
76 | double t, tMaxkin, tMin, tMax; |
77 | |
78 | |
79 | tMaxkin = 4. * beamMomCM * EtaMomCM; |
80 | tMax = 0.2; |
81 | |
82 | tMin = abs(pow(EtaMass,4)/(2*cmEnergy*cmEnergy) - (beamMomCM-EtaMomCM)*(beamMomCM-EtaMomCM)); |
83 | |
84 | double Irandom = gRandom->Uniform(); |
85 | |
86 | t = -log( Irandom*(exp(-m_slope*tMax) - exp(-m_slope*tMin)) + exp(-m_slope*tMin))/m_slope; |
87 | |
88 | TVector3 EtaMom3CM; |
89 | double thetaCM = 2.*sqrt((t- tMin)/tMaxkin); |
90 | |
91 | double phiCM = random( -kPi, kPi ); |
92 | |
93 | EtaMom3CM.SetMagThetaPhi( EtaMomCM, thetaCM, phiCM); |
94 | |
95 | TLorentzVector EtaMom4CM( EtaMom3CM, sqrt( EtaMom3CM.Mag2() + EtaMass * EtaMass ) ); |
96 | TLorentzVector recoilMom4CM = cmCM - EtaMom4CM; |
97 | |
98 | |
99 | |
100 | TLorentzVector EtaMom4Lab = cm2labBoost * EtaMom4CM; |
101 | TLorentzVector recoil = cm2labBoost * recoilMom4CM; |
102 | |
103 | double trecoil = 2*recoil.M()*(recoil.E()-recoil.M()); |
| Value stored to 'trecoil' during its initialization is never read |
104 | |
105 | |
106 | |
107 | vector< TLorentzVector > allPart; |
108 | allPart.push_back( beam ); |
109 | allPart.push_back( EtaMom4Lab ); |
110 | allPart.push_back( recoil ); |
111 | |
112 | |
113 | |
114 | |
115 | |
116 | |
117 | double genWeight = 1.; |
118 | |
119 | return new Kinematics( allPart, genWeight ); |
120 | } |
121 | double |
122 | GammaZToXZ::cmMomentum( double M, double m1, double m2 ) const { |
123 | |
124 | |
125 | |
126 | double num1 = ( M * M - ( m1 + m2 ) * ( m1 + m2 ) ); |
127 | double num2 = ( M * M - ( m1 - m2 ) * ( m1 - m2 ) ); |
128 | |
129 | return( sqrt( num1 * num2 ) / ( 2 * M ) ); |
130 | } |
131 | |
132 | double |
133 | GammaZToXZ::random( double low, double hi ) const { |
134 | |
135 | return( ( hi - low ) * gRandom->Uniform() + low ); |
136 | } |
137 | |