Changeset c1ce3fe in git for modules/PhotonConversions.cc
- Timestamp:
- Jun 29, 2015, 10:28:05 PM (9 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- ff37d75
- Parents:
- 839deb7
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/PhotonConversions.cc
r839deb7 rc1ce3fe 57 57 58 58 PhotonConversions::PhotonConversions() : 59 f DecayXsec(0), fConversionMap(0), fItInputArray(0)59 fItInputArray(0), fConversionMap(0), fDecayXsec(0) 60 60 { 61 61 fDecayXsec = new TF1; … … 74 74 { 75 75 fRadius = GetDouble("Radius", 1.0); 76 fRadius2 = 77 76 fRadius2 = fRadius*fRadius; 77 78 78 fHalfLength = GetDouble("HalfLength", 3.0); 79 79 80 80 fEtaMax = GetDouble("EtaMax", 5.0); 81 81 fEtaMin = GetDouble("EtaMin", 2.0); 82 83 fStep = GetDouble("Step", 0.1); // in meters 84 85 fConversionMap->Compile(GetString("ConversionMap", "0.")); 86 87 fDecayXsec->Compile("1 - 4/3*x*(1-x)"); 88 fDecayXsec->SetRange(0.0,1.0); 89 82 83 fStep = GetDouble("Step", 0.1); // in meters 84 85 fConversionMap->Compile(GetString("ConversionMap", "0.0")); 86 87 #if ROOT_VERSION_CODE < ROOT_VERSION(6,04,00) 88 fDecayXsec->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)"); 89 #else 90 fDecayXsec->GetFormula()->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)"); 91 #endif 92 fDecayXsec->SetRange(0.0, 1.0); 93 90 94 // import array with output from filter/classifier module 91 95 … … 111 115 void PhotonConversions::Process() 112 116 { 113 Candidate *candidate, * mother, *ep, *em;117 Candidate *candidate, *ep, *em; 114 118 TLorentzVector candidatePosition, candidateMomentum; 115 119 TVector3 pos_i; 116 Double_t px, py, pz, pt, pt2, e, eta, phi , q;120 Double_t px, py, pz, pt, pt2, e, eta, phi; 117 121 Double_t x, y, z, t; 118 122 Double_t x_t, y_t, z_t, r_t; … … 123 127 Double_t rate, p_conv, x1, x2; 124 128 Bool_t converted; 125 126 const Double_t c_light = 2.99792458E8; 127 129 128 130 fItInputArray->Reset(); 129 131 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 130 132 { 131 132 if(candidate->PID != 22) fOutputArray->Add(candidate); 133 134 if(candidate->PID != 22) 135 { 136 fOutputArray->Add(candidate); 137 } 133 138 else 134 139 { … … 138 143 y = candidatePosition.Y()*1.0E-3; 139 144 z = candidatePosition.Z()*1.0E-3; 140 q = candidate->Charge;141 145 142 146 // check that particle position is inside the cylinder 143 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength) 144 { 145 continue; 146 } 147 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength) continue; 147 148 148 149 px = candidateMomentum.Px(); … … 154 155 phi = candidateMomentum.Phi(); 155 156 e = candidateMomentum.E(); 156 157 157 158 if(eta < fEtaMin || eta > fEtaMax) continue; 158 159 159 160 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0 160 161 tmp = px*y - py*x; 161 162 discr2 = pt2*fRadius2 - tmp*tmp; … … 182 183 183 184 // final position 184 185 185 x_t = x + px*t; 186 186 y_t = y + py*t; 187 187 z_t = z + pz*t; 188 188 189 189 r_t = TMath::Sqrt(x_t*x_t + y_t*y_t + z_t*z_t); 190 191 192 // here starts conversion code 190 191 192 // here starts conversion code 193 193 nsteps = Int_t(r_t/fStep); 194 194 195 195 x_i = x; 196 196 y_i = y; 197 197 z_i = z; 198 198 199 199 dt = t/nsteps; 200 200 201 201 converted = false; 202 203 for(i =0; i<nsteps; ++i)202 203 for(i = 0; i < nsteps; ++i) 204 204 { 205 206 x_i += px*dt; 207 y_i += py*dt; 208 z_i += pz*dt; 205 x_i += px*dt; 206 y_i += py*dt; 207 z_i += pz*dt; 209 208 pos_i.SetXYZ(x_i,y_i,z_i); 210 211 //convert photon position into cylindrical coordinates, cylindrical r,phi,z !!212 213 209 210 // convert photon position into cylindrical coordinates, cylindrical r,phi,z !! 211 212 r_i = TMath::Sqrt(x_i*x_i + y_i*y_i); 214 213 phi_i = pos_i.Phi(); 215 214 216 215 // read conversion rate/meter from card 217 rate = fConversionMap->Eval(r_i, phi_i, z_i); 218 219 // convert into conversion probability 220 p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate); 221 222 223 //case conversion occurs 224 if(gRandom->Uniform() < p_conv) 225 { 226 converted = true; 227 228 // generate x1 and x2, the fraction of the photon energy taken resp. by e+ and e- 229 x1 = fDecayXsec->GetRandom(); 230 x2 = 1 - x1; 231 232 mother = candidate; 233 ep = static_cast<Candidate*>(candidate->Clone()); 234 em = static_cast<Candidate*>(candidate->Clone()); 235 236 ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 237 em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 238 239 ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e); 240 em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e); 241 242 ep->PID = -11; 243 em->PID = 11; 244 245 ep->Charge = 1.0; 246 em->Charge = -1.0; 247 248 ep->IsFromConversion = 1; 249 em->IsFromConversion = 1; 250 251 fOutputArray->Add(em); 252 fOutputArray->Add(ep); 253 254 break; 255 256 } 257 216 rate = fConversionMap->Eval(r_i, phi_i, z_i); 217 218 // convert into conversion probability 219 p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate); 220 221 // case conversion occurs 222 if(gRandom->Uniform() < p_conv) 223 { 224 converted = true; 225 226 // generate x1 and x2, the fraction of the photon energy taken resp. by e+ and e- 227 x1 = fDecayXsec->GetRandom(); 228 x2 = 1 - x1; 229 230 ep = static_cast<Candidate*>(candidate->Clone()); 231 em = static_cast<Candidate*>(candidate->Clone()); 232 233 ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 234 em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3); 235 236 ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e); 237 em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e); 238 239 ep->PID = -11; 240 em->PID = 11; 241 242 ep->Charge = 1.0; 243 em->Charge = -1.0; 244 245 ep->IsFromConversion = 1; 246 em->IsFromConversion = 1; 247 248 fOutputArray->Add(em); 249 fOutputArray->Add(ep); 250 251 break; 252 } 258 253 } 259 260 if(!converted) fOutputArray->Add(candidate); 261 254 if(!converted) fOutputArray->Add(candidate); 262 255 } 263 256 }
Note:
See TracChangeset
for help on using the changeset viewer.