Changes in modules/PhotonConversions.cc [341014c:bdfe7f2] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/PhotonConversions.cc
r341014c rbdfe7f2 17 17 */ 18 18 19 19 20 /** \class PhotonConversions 20 21 * … … 29 30 30 31 #include "classes/DelphesClasses.h" 32 #include "classes/DelphesFactory.h" 31 33 #include "classes/DelphesCylindricalFormula.h" 32 #include "classes/DelphesFactory.h" 33 34 35 #include "ExRootAnalysis/ExRootResult.h" 36 #include "ExRootAnalysis/ExRootFilter.h" 34 37 #include "ExRootAnalysis/ExRootClassifier.h" 35 #include "ExRootAnalysis/ExRootFilter.h" 36 #include "ExRootAnalysis/ExRootResult.h" 37 38 39 #include "TMath.h" 40 #include "TF1.h" 41 #include "TString.h" 42 #include "TFormula.h" 43 #include "TRandom3.h" 44 #include "TObjArray.h" 38 45 #include "TDatabasePDG.h" 39 #include "TF1.h"40 #include "TFormula.h"41 46 #include "TLorentzVector.h" 42 #include "TMath.h"43 #include "TObjArray.h"44 #include "TRandom3.h"45 #include "TString.h"46 47 #include "TVector3.h" 47 48 48 49 #include <algorithm> 50 #include <stdexcept> 49 51 #include <iostream> 50 52 #include <sstream> 51 #include <stdexcept>52 53 53 54 using namespace std; … … 58 59 fItInputArray(0), fConversionMap(0), fDecayXsec(0) 59 60 { 60 fDecayXsec = new TF1("decayXsec", 61 fDecayXsec = new TF1("decayXsec","1.0 - 4.0/3.0 * x * (1.0 - x)", 0.0, 1.0); 61 62 fConversionMap = new DelphesCylindricalFormula; 62 63 } … … 73 74 { 74 75 fRadius = GetDouble("Radius", 1.0); 75 fRadius2 = fRadius *fRadius;76 fRadius2 = fRadius*fRadius; 76 77 77 78 fHalfLength = GetDouble("HalfLength", 3.0); … … 121 122 122 123 fItInputArray->Reset(); 123 while((candidate = static_cast<Candidate 124 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 124 125 { 125 126 … … 132 133 candidatePosition = candidate->Position; 133 134 candidateMomentum = candidate->Momentum; 134 x = candidatePosition.X() *1.0E-3;135 y = candidatePosition.Y() *1.0E-3;136 z = candidatePosition.Z() *1.0E-3;135 x = candidatePosition.X()*1.0E-3; 136 y = candidatePosition.Y()*1.0E-3; 137 z = candidatePosition.Z()*1.0E-3; 137 138 138 139 // check that particle position is inside the cylinder … … 151 152 152 153 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0 153 tmp = px * y - py *x;154 discr2 = pt2 * fRadius2 - tmp *tmp;154 tmp = px*y - py*x; 155 discr2 = pt2*fRadius2 - tmp*tmp; 155 156 156 157 if(discr2 < 0.0) … … 160 161 } 161 162 162 tmp = px * x + py *y;163 tmp = px*x + py*y; 163 164 discr = TMath::Sqrt(discr2); 164 t1 = (-tmp + discr) /pt2;165 t2 = (-tmp - discr) /pt2;165 t1 = (-tmp + discr)/pt2; 166 t2 = (-tmp - discr)/pt2; 166 167 t = (t1 < 0.0) ? t2 : t1; 167 168 168 z_t = z + pz *t;169 z_t = z + pz*t; 169 170 if(TMath::Abs(z_t) > fHalfLength) 170 171 { … … 175 176 176 177 // final position 177 x_t = x + px * t; 178 y_t = y + py * t; 179 z_t = z + pz * t; 180 181 r_t = TMath::Sqrt(x_t * x_t + y_t * y_t + z_t * z_t); 178 x_t = x + px*t; 179 y_t = y + py*t; 180 z_t = z + pz*t; 181 182 r_t = TMath::Sqrt(x_t*x_t + y_t*y_t + z_t*z_t); 183 182 184 183 185 // here starts conversion code 184 nsteps = Int_t(r_t /fStep);186 nsteps = Int_t(r_t/fStep); 185 187 186 188 x_i = x; … … 188 190 z_i = z; 189 191 190 dt = t /nsteps;192 dt = t/nsteps; 191 193 192 194 converted = false; … … 194 196 for(i = 0; i < nsteps; ++i) 195 197 { 196 x_i += px *dt;197 y_i += py *dt;198 z_i += pz *dt;199 pos_i.SetXYZ(x_i, y_i,z_i);198 x_i += px*dt; 199 y_i += py*dt; 200 z_i += pz*dt; 201 pos_i.SetXYZ(x_i,y_i,z_i); 200 202 201 203 // convert photon position into cylindrical coordinates, cylindrical r,phi,z !! 202 204 203 r_i = TMath::Sqrt(x_i * x_i + y_i *y_i);205 r_i = TMath::Sqrt(x_i*x_i + y_i*y_i); 204 206 phi_i = pos_i.Phi(); 205 207 … … 208 210 209 211 // convert into conversion probability 210 p_conv = 1 - TMath::Exp(-7.0 / 9.0 * fStep *rate);212 p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate); 211 213 212 214 // case conversion occurs … … 219 221 x2 = 1 - x1; 220 222 221 ep = static_cast<Candidate 222 em = static_cast<Candidate 223 224 ep->Position.SetXYZT(x_i * 1.0E3, y_i * 1.0E3, z_i * 1.0E3, candidatePosition.T() + nsteps * dt * e *1.0E3);225 em->Position.SetXYZT(x_i * 1.0E3, y_i * 1.0E3, z_i * 1.0E3, candidatePosition.T() + nsteps * dt * e *1.0E3);226 227 ep->Momentum.SetPtEtaPhiE(x1 * pt, eta, phi, x1 *e);228 em->Momentum.SetPtEtaPhiE(x2 * pt, eta, phi, x2 *e);223 ep = static_cast<Candidate*>(candidate->Clone()); 224 em = static_cast<Candidate*>(candidate->Clone()); 225 226 ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 227 em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 228 229 ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e); 230 em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e); 229 231 230 232 ep->PID = -11; … … 249 251 250 252 //------------------------------------------------------------------------------ 253
Note:
See TracChangeset
for help on using the changeset viewer.