Bug Summary

File:libraries/AMPTOOLS_MCGEN/GammaZToXZ.cc
Location:line 103, column 8
Description:Value stored to 'trecoil' during its initialization is never read

Annotated Source Code

1/*
2 * GammaZToXZ.cc
3 * GlueXTools
4 *
5 * Copied from GammaPToXP. Elton 5/26/2020
6 *
7 * Created by Matthew Shepherd on 1/22/10.
8 * Copyright 2010 Home. All rights reserved.
9 *
10 * Copied from
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
22GammaZToXZ::GammaZToXZ( float massX, TString beamConfigFile, Double_t Bslope) :
23m_target( 0, 0, 0, ParticleMass(Pb208)),
24m_slope( Bslope ),
25m_recoil (ParticleMass(Pb208)),
26m_childMass( 0 )
27{
28 kMproton=ParticleMass(Proton);
29 kMneutron=ParticleMass(Neutron);
30 // kMZ = 108.; // mass of Sn116
31 // kMZ = 208.*0.931494; // use mass of Pb as it is in the particle table
32 kMZ = ParticleMass(Pb208); // use mass of Pb as it is in the particle table
33 kMPion = ParticleMass(PiPlus);
34 kMPi0 = ParticleMass(Pi0);
35 kMKaon = ParticleMass(KPlus);
36
37 m_childMass.push_back( massX );
38
39 // get beam properties from configuration file
40 BeamProperties beamProp(beamConfigFile);
41 cobrem_vs_E = (TH1D*)beamProp.GetFlux();
42 cobrem_vs_E->GetName();
43}
44
45Kinematics*
46GammaZToXZ::generate(){
47
48
49 const double kPi = 3.14159;
50
51// Double_t kMProton = ParticleMass(Proton);
52 double EtaMass=m_childMass[0];
53 double recoilMass = m_recoil;
54// Double_t masses[2] = {EtaMass,recoilMass};
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/*cout << "M=" << beam.M() << " beam="; beam.Print();
71 cout << "M=" << m_target.M() << " m_target="; m_target.Print();
72 cout << "M=" << cm.M() << " cm="; cm.Print();*/
73
74 // generate an exponetial t-slope between t1 and t2
75
76 double t, tMaxkin, tMin, tMax;
77 // generate the t-distribution. t is positive here (i.e. should be -t)
78
79 tMaxkin = 4. * beamMomCM * EtaMomCM;
80 tMax = 0.2; // restrict max to make more efficient for Primakoff generation (about 2. deg at 0.05 GeV-2)
81 // tMax = 1.; // restrict max to make more efficient for Primakoff generation
82 tMin = abs(pow(EtaMass,4)/(2*cmEnergy*cmEnergy) - (beamMomCM-EtaMomCM)*(beamMomCM-EtaMomCM)); // treating t as positive
83
84 double Irandom = gRandom->Uniform();
85 // generate random t with exponential between tMin and tMax
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); // acos( 1. - 2.*t/tMax ) -> use small angle approximation to avoid roundoff. For heavy target, CM=Lab. Opposite to target
90 // double thetaCM = acos( 1. - 2.*t/tMaxkin );
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
103double trecoil = 2*recoil.M()*(recoil.E()-recoil.M());
Value stored to 'trecoil' during its initialization is never read
104
105// cout << " tMin=" << tMin << " tMax=" << tMax << " t=" << t << " trecoil=" << trecoil << " m_slope=" << m_slope << " thetaCM=" << thetaCM*180/kPi << " phiCM=" << phiCM*180./kPi << endl;
106
107 vector< TLorentzVector > allPart;
108 allPart.push_back( beam );
109 allPart.push_back( EtaMom4Lab );
110 allPart.push_back( recoil );
111
112/*cout << "M=" << beam.M() << " beam="; beam.Print();
113 cout << "cmCM=" << cmCM.M() << " cmCM="; cmCM.Print();
114 cout << "M=" << recoil.M() << " recoil="; recoil.Print();
115 cout << " m_slope=" << m_slope << " M=" << EtaMom4Lab.M() << " EtaMom4Lab="; EtaMom4Lab.Print(); cout << endl;*/
116
117 double genWeight = 1.;
118
119 return new Kinematics( allPart, genWeight );
120}
121double
122GammaZToXZ::cmMomentum( double M, double m1, double m2 ) const {
123
124 // mini PDG Eq: 38.16
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
132double
133GammaZToXZ::random( double low, double hi ) const {
134
135 return( ( hi - low ) * gRandom->Uniform() + low );
136}
137