Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:28 

0001 #include "GeneratorInterface/CosmicMuonGenerator/interface/SingleParticleEvent.h"
0002 
0003 void SingleParticleEvent::create(
0004     int id, double px, double py, double pz, double e, double m, double vx, double vy, double vz, double t0) {
0005   ID = ID_in = id;
0006   Px = Px_in = px;
0007   Py = Py_in = py;
0008   Pz = Pz_in = pz;
0009   E = E_in = e;
0010   M = M_in = m;
0011   Vx = Vx_in = vx;
0012   Vy = Vy_in = vy;
0013   Vz = Vz_in = vz;
0014   T0 = T0_in = t0;
0015   HitTarget = false;
0016 }
0017 
0018 void SingleParticleEvent::propagate(double ElossScaleFac,
0019                                     double RadiusTarget,
0020                                     double Z_DistTarget,
0021                                     double Z_CentrTarget,
0022                                     bool TrackerOnly,
0023                                     bool MTCCHalf) {
0024   MTCC = MTCCHalf;  //need to know this boolean in absVzTmp()
0025   // calculated propagation direction
0026   dX = Px / absmom();
0027   dY = Py / absmom();
0028   dZ = Pz / absmom();
0029   // propagate with decreasing step size
0030   tmpVx = Vx;
0031   tmpVy = Vy;
0032   tmpVz = Vz;
0033   double RadiusTargetEff = RadiusTarget;
0034   double Z_DistTargetEff = Z_DistTarget;
0035   double Z_CentrTargetEff = Z_CentrTarget;
0036   if (TrackerOnly == true) {
0037     RadiusTargetEff = RadiusTracker;
0038     Z_DistTargetEff = Z_DistTracker;
0039   }
0040   HitTarget = true;
0041   if (HitTarget == true) {
0042     HitTarget = false;
0043     double stepSize = MinStepSize * 100000.;
0044     double acceptR = RadiusTargetEff + stepSize;
0045     double acceptZ = Z_DistTargetEff + stepSize;
0046     bool continuePropagation = true;
0047     while (continuePropagation) {
0048       //if (tmpVy < -acceptR) continuePropagation = false;
0049       if (dY < 0. && tmpVy < -acceptR)
0050         continuePropagation = false;
0051       if (dY >= 0. && tmpVy > acceptR)
0052         continuePropagation = false;
0053       //if (absVzTmp() < acceptZ && rVxyTmp() < acceptR){
0054       if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0055         HitTarget = true;
0056         continuePropagation = false;
0057       }
0058       if (continuePropagation)
0059         updateTmp(stepSize);
0060     }
0061   }
0062   if (HitTarget == true) {
0063     HitTarget = false;
0064     double stepSize = MinStepSize * 10000.;
0065     double acceptR = RadiusTargetEff + stepSize;
0066     double acceptZ = Z_DistTargetEff + stepSize;
0067     bool continuePropagation = true;
0068     while (continuePropagation) {
0069       //if (tmpVy < -acceptR) continuePropagation = false;
0070       if (dY < 0. && tmpVy < -acceptR)
0071         continuePropagation = false;
0072       if (dY >= 0. && tmpVy > acceptR)
0073         continuePropagation = false;
0074       //if (absVzTmp() < acceptZ && rVxyTmp() < acceptR){
0075       if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0076         HitTarget = true;
0077         continuePropagation = false;
0078       }
0079       if (continuePropagation)
0080         updateTmp(stepSize);
0081     }
0082   }
0083   if (HitTarget == true) {
0084     HitTarget = false;
0085     double stepSize = MinStepSize * 1000.;
0086     double acceptR = RadiusTargetEff + stepSize;
0087     double acceptZ = Z_DistTargetEff + stepSize;
0088     bool continuePropagation = true;
0089     while (continuePropagation) {
0090       //if (tmpVy < -acceptR) continuePropagation = false;
0091       if (dY < 0. && tmpVy < -acceptR)
0092         continuePropagation = false;
0093       if (dY >= 0. && tmpVy > acceptR)
0094         continuePropagation = false;
0095       //if (absVzTmp() < acceptZ && rVxyTmp() < acceptR){
0096       if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0097         HitTarget = true;
0098         continuePropagation = false;
0099       }
0100       if (continuePropagation)
0101         updateTmp(stepSize);
0102     }
0103   }
0104   if (HitTarget == true) {
0105     HitTarget = false;
0106     double stepSize = MinStepSize * 100.;
0107     double acceptR = RadiusTargetEff + stepSize;
0108     double acceptZ = Z_DistTargetEff + stepSize;
0109     bool continuePropagation = true;
0110     while (continuePropagation) {
0111       //if (tmpVy < -acceptR) continuePropagation = false;
0112       if (dY < 0. && tmpVy < -acceptR)
0113         continuePropagation = false;
0114       if (dY >= 0. && tmpVy > acceptR)
0115         continuePropagation = false;
0116       //if (absVzTmp() < acceptZ && rVxyTmp() < acceptR){
0117       if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0118         HitTarget = true;
0119         continuePropagation = false;
0120       }
0121       if (continuePropagation)
0122         updateTmp(stepSize);
0123     }
0124   }
0125   if (HitTarget == true) {
0126     HitTarget = false;
0127     double stepSize = MinStepSize * 10.;
0128     double acceptR = RadiusTargetEff + stepSize;
0129     double acceptZ = Z_DistTargetEff + stepSize;
0130     bool continuePropagation = true;
0131     while (continuePropagation) {
0132       //if (tmpVy < -acceptR) continuePropagation = false;
0133       if (dY < 0. && tmpVy < -acceptR)
0134         continuePropagation = false;
0135       if (dY >= 0. && tmpVy > acceptR)
0136         continuePropagation = false;
0137       //if (absVzTmp() < acceptZ && rVxyTmp() < acceptR){
0138       if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0139         HitTarget = true;
0140         continuePropagation = false;
0141       }
0142       if (continuePropagation)
0143         updateTmp(stepSize);
0144     }
0145   }
0146   if (HitTarget == true) {
0147     HitTarget = false;
0148     double stepSize = MinStepSize * 1.;
0149     double acceptR = RadiusTargetEff + stepSize;
0150     double acceptZ = Z_DistTargetEff + stepSize;
0151     bool continuePropagation = true;
0152     while (continuePropagation) {
0153       //if (tmpVy < -acceptR) continuePropagation = false;
0154       if (dY < 0. && tmpVy < -acceptR)
0155         continuePropagation = false;
0156       if (dY >= 0. && tmpVy > acceptR)
0157         continuePropagation = false;
0158       //if (0 < absVzTmp()){ //only check for MTCC setup in last step of propagation, need fine stepSize
0159       if (absVzTmp() < acceptZ && rVxyTmp() < acceptR) {
0160         if (std::fabs(tmpVz - Z_CentrTargetEff) < acceptZ && rVxyTmp() < acceptR) {
0161           HitTarget = true;
0162           continuePropagation = false;
0163         }
0164       }
0165       if (continuePropagation)
0166         updateTmp(stepSize);
0167     }
0168   }
0169   // actual propagation + energy loss
0170   if (HitTarget == true) {
0171     HitTarget = false;
0172     //int nAir = 0; int nWall = 0; int nRock = 0; int nClay = 0; int nPlug = 0;
0173     int nMat[6] = {0, 0, 0, 0, 0, 0};
0174     double stepSize = MinStepSize * 1.;  // actual step size
0175     double acceptR = RadiusCMS + stepSize;
0176     double acceptZ = Z_DistCMS + stepSize;
0177     if (TrackerOnly == true) {
0178       acceptR = RadiusTracker + stepSize;
0179       acceptZ = Z_DistTracker + stepSize;
0180     }
0181     bool continuePropagation = true;
0182     while (continuePropagation) {
0183       //if (Vy < -acceptR) continuePropagation = false;
0184       if (dY < 0. && tmpVy < -acceptR)
0185         continuePropagation = false;
0186       if (dY >= 0. && tmpVy > acceptR)
0187         continuePropagation = false;
0188       //if (absVz() < acceptZ && rVxy() < acceptR){
0189       if (std::fabs(Vz - Z_CentrTargetEff) < acceptZ && rVxy() < acceptR) {
0190         HitTarget = true;
0191         continuePropagation = false;
0192       }
0193       if (continuePropagation)
0194         update(stepSize);
0195 
0196       int Mat = inMat(Vx, Vy, Vz, PlugVx, PlugVz, ClayWidth);
0197 
0198       nMat[Mat]++;
0199     }
0200 
0201     if (HitTarget) {
0202       double lPlug = double(nMat[Plug]) * stepSize;
0203       double lWall = double(nMat[Wall]) * stepSize;
0204       double lAir = double(nMat[Air]) * stepSize;
0205       double lClay = double(nMat[Clay]) * stepSize;
0206       double lRock = double(nMat[Rock]) * stepSize;
0207       //double lUnknown = double(nMat[Unknown])*stepSize;
0208 
0209       double waterEquivalents =
0210           (lAir * RhoAir + lWall * RhoWall + lRock * RhoRock + lClay * RhoClay + lPlug * RhoPlug) * ElossScaleFac /
0211           10.;  // [g cm^-2]
0212       subtractEloss(waterEquivalents);
0213       if (E < MuonMass)
0214         HitTarget = false;  // muon stopped in the material around the target
0215     }
0216   }
0217   // end of propagation part
0218 }
0219 
0220 void SingleParticleEvent::update(double stepSize) {
0221   Vx += stepSize * dX;
0222   Vy += stepSize * dY;
0223   Vz += stepSize * dZ;
0224 }
0225 
0226 void SingleParticleEvent::updateTmp(double stepSize) {
0227   tmpVx += stepSize * dX;
0228   tmpVy += stepSize * dY;
0229   tmpVz += stepSize * dZ;
0230 }
0231 
0232 void SingleParticleEvent::subtractEloss(double waterEquivalents) {
0233   double L10E = log10(E);
0234   // parameters for standard rock (PDG 2004, page 230)
0235   double A = (1.91514 + 0.254957 * L10E) / 1000.;                              // a [GeV g^-1 cm^2]
0236   double B = (0.379763 + 1.69516 * L10E - 0.175026 * L10E * L10E) / 1000000.;  // b [g^-1 cm^2]
0237   double EPS = A / B;                                                          // epsilon [GeV]
0238   E = (E + EPS) * exp(-B * waterEquivalents) - EPS;                            // updated energy
0239   double oldAbsMom = absmom();
0240   double newAbsMom = sqrt(E * E - MuonMass * MuonMass);
0241   Px = Px * newAbsMom / oldAbsMom;  // updated px
0242   Py = Py * newAbsMom / oldAbsMom;  // updated py
0243   Pz = Pz * newAbsMom / oldAbsMom;  // updated pz
0244 }
0245 
0246 double SingleParticleEvent::Eloss(double waterEquivalents, double Energy) {
0247   double L10E = log10(Energy);
0248   // parameters for standard rock (PDG 2004, page 230)
0249   double A = (1.91514 + 0.254957 * L10E) / 1000.;                              // a [GeV g^-1 cm^2]
0250   double B = (0.379763 + 1.69516 * L10E - 0.175026 * L10E * L10E) / 1000000.;  // b [g^-1 cm^2]
0251   double EPS = A / B;                                                          // epsilon [GeV]
0252   double newEnergy = (Energy + EPS) * exp(-B * waterEquivalents) - EPS;        // updated energy
0253   double EnergyLoss = Energy - newEnergy;
0254   return EnergyLoss;
0255 }
0256 
0257 void SingleParticleEvent::setEug(double Eug) { E_ug = Eug; }
0258 
0259 double SingleParticleEvent::Eug() { return E_ug; }
0260 
0261 double SingleParticleEvent::deltaEmin(double E_sf) {
0262   double dE = Eloss(waterEquivalents, E_sf);
0263   return E_ug - (E_sf - dE);
0264 }
0265 
0266 void SingleParticleEvent::SurfProj(double Vx_in,
0267                                    double Vy_in,
0268                                    double Vz_in,
0269                                    double Px_in,
0270                                    double Py_in,
0271                                    double Pz_in,
0272                                    double& Vx_up,
0273                                    double& Vy_up,
0274                                    double& Vz_up) {
0275   //determine vertex of muon at Surface (+PlugWidth)
0276   double dy = Vy_in - (SurfaceOfEarth + PlugWidth);
0277   Vy_up = Vy_in - dy;
0278   Vx_up = Vx_in - dy * Px_in / Py_in;
0279   Vz_up = Vz_in - dy * Pz_in / Py_in;
0280   if (Debug)
0281     std::cout << "Vx_up=" << Vx_up << " Vy_up=" << Vy_up << " Vz_up=" << Vz_up << std::endl;
0282 }
0283 
0284 double SingleParticleEvent::absVzTmp() {
0285   if (MTCC == true) {
0286     return tmpVz;  //need sign to be sure muon hits half of CMS with MTCC setup
0287   } else {
0288     return std::fabs(tmpVz);
0289   }
0290 }
0291 
0292 double SingleParticleEvent::rVxyTmp() { return sqrt(tmpVx * tmpVx + tmpVy * tmpVy); }
0293 
0294 bool SingleParticleEvent::hitTarget() { return HitTarget; }
0295 
0296 int SingleParticleEvent::id_in() { return ID_in; }
0297 
0298 double SingleParticleEvent::px_in() { return Px_in; }
0299 
0300 double SingleParticleEvent::py_in() { return Py_in; }
0301 
0302 double SingleParticleEvent::pz_in() { return Pz_in; }
0303 
0304 double SingleParticleEvent::e_in() { return E_in; }
0305 
0306 double SingleParticleEvent::m_in() { return M_in; }
0307 
0308 double SingleParticleEvent::vx_in() { return Vx_in; }
0309 
0310 double SingleParticleEvent::vy_in() { return Vy_in; }
0311 
0312 double SingleParticleEvent::vz_in() { return Vz_in; }
0313 
0314 double SingleParticleEvent::t0_in() { return T0_in; }
0315 
0316 int SingleParticleEvent::id() { return ID; }
0317 
0318 double SingleParticleEvent::px() { return Px; }
0319 
0320 double SingleParticleEvent::py() { return Py; }
0321 
0322 double SingleParticleEvent::pz() { return Pz; }
0323 
0324 double SingleParticleEvent::e() { return E; }
0325 
0326 double SingleParticleEvent::m() { return M; }
0327 
0328 double SingleParticleEvent::vx() { return Vx; }
0329 
0330 double SingleParticleEvent::vy() { return Vy; }
0331 
0332 double SingleParticleEvent::vz() { return Vz; }
0333 
0334 double SingleParticleEvent::t0() { return T0; }
0335 
0336 double SingleParticleEvent::WaterEquivalents() { return waterEquivalents; }
0337 
0338 double SingleParticleEvent::phi() {
0339   double phiXZ = atan2(Px, Pz);
0340   if (phiXZ < 0.)
0341     phiXZ = phiXZ + TwoPi;
0342   return phiXZ;
0343 }
0344 
0345 double SingleParticleEvent::theta() { return atan2(sqrt(Px * Px + Pz * Pz), -Py); }
0346 
0347 double SingleParticleEvent::absmom() { return sqrt(Px * Px + Py * Py + Pz * Pz); }
0348 
0349 double SingleParticleEvent::absVz() { return std::fabs(Vz); }
0350 
0351 double SingleParticleEvent::rVxy() { return sqrt(Vx * Vx + Vy * Vy); }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!