GFlash/src/GflashShowino.cc

0001 #include "SimGeneral/GFlash/interface/GflashShowino.h"
0002 #include <CLHEP/Random/Randomize.h>
0003
0004 GflashShowino::GflashShowino()
0005     : theShowerType(-1),
0006       theEnergy(0),
0007       thePositionAtShower(0, 0, 0),
0008       thePathLengthAtShower(0),
0009       thePathLengthOnEcal(0),
0010       theStepLengthToHcal(0),
0011       theStepLengthToOut(0),
0012       thePathLength(0),
0013       theGlobalTime(0),
0014       thePosition(0, 0, 0),
0015       theEnergyDeposited(0) {
0016   theHelix = new GflashTrajectory;
0017 }
0018
0019 GflashShowino::~GflashShowino() { delete theHelix; }
0020
0021 void GflashShowino::initialize(int showerType,
0022                                double energy,
0023                                double globalTime,
0024                                double charge,
0025                                Gflash3Vector &position,
0026                                Gflash3Vector &momentum,
0027                                double magneticField) {
0028   theEnergy = energy;
0029   theGlobalTime = globalTime;
0030   theEnergyDeposited = 0.0;
0031
0032   // inside the magnetic field (tesla unit);
0033   theHelix->initializeTrajectory(momentum, position, charge, magneticField);
0034
0035   if (showerType < 100) {
0036     thePositionAtShower = position;
0037     thePosition = thePositionAtShower;
0038     theShowerType = showerType;
0039
0040   } else {
0041     // this input is from FastSimulation
0042     // 1. simulate the shower starting position
0043     thePositionAtShower = simulateFirstInteractionPoint(showerType, position);
0044     thePosition = thePositionAtShower;
0045
0046     // 2. find shower type depending on where is the shower starting point
0047     theShowerType = Gflash::findShowerType(thePositionAtShower);
0048   }
0049
0050   evaluateLengths();
0051 }
0052
0053 void GflashShowino::updateShowino(double deltaStep) {
0054   thePathLength += deltaStep;
0055   // trajectory point of showino along the shower depth at the pathLength
0056   GflashTrajectoryPoint trajectoryShowino;
0057   theHelix->getGflashTrajectoryPoint(trajectoryShowino, thePathLength);
0058
0059   thePosition = trajectoryShowino.getPosition();
0060
0061   theGlobalTime += (theEnergy / 100.0) * deltaStep / 30.0;  //@@@calculate exact time change in nsec
0062 }
0063
0064 void GflashShowino::evaluateLengths() {
0065   // thePathLengthAtShower: path Length from the origin to the shower starting
0066   // point in cm theStepLengthToOut: the total path length from the starting
0067   // point of
0068   //                    shower to the maximum distance inside paramerized
0069   //                    envelopes
0070   double eta = thePosition.getEta();
0071
0072   if (std::fabs(eta) < Gflash::EtaMax[Gflash::kESPM]) {
0073     thePathLengthOnEcal = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
0074     thePathLengthAtShower = theHelix->getPathLengthAtRhoEquals(thePosition.getRho());
0075     double pathLengthAtHcalBack = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
0076     if (pathLengthAtHcalBack > 0) {
0077       theStepLengthToOut = std::min(300., pathLengthAtHcalBack - thePathLengthAtShower);
0078     } else {
0079       theStepLengthToOut = 200.;
0080     }
0081     theStepLengthToHcal = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]) - thePathLengthAtShower;
0082
0083   } else if (std::fabs(eta) < Gflash::EtaMax[Gflash::kENCA]) {
0084     double zsign = (eta > 0) ? 1.0 : -1.0;
0085     thePathLengthOnEcal = theHelix->getPathLengthAtZ(zsign * Gflash::ZFrontCrystalEE);
0086     thePathLengthAtShower = theHelix->getPathLengthAtZ(thePosition.getZ());
0087     theStepLengthToOut =
0088         std::min(300., theHelix->getPathLengthAtZ(zsign * Gflash::Zmax[Gflash::kHE]) - thePathLengthAtShower);
0089     theStepLengthToHcal = theHelix->getPathLengthAtZ(zsign * Gflash::Zmin[Gflash::kHE]) - thePathLengthAtShower;
0090   } else {
0091     //@@@extend for HF later
0092     theStepLengthToOut = 200.0;
0093   }
0094
0095   thePathLength = thePathLengthAtShower;
0096 }
0097
0098 Gflash3Vector &GflashShowino::simulateFirstInteractionPoint(int fastSimShowerType, Gflash3Vector &position) {
0099   // determine the shower starting point (ssp):
0100   // the position at the entrance + the mean free path till the inelastic
0101   // interaction inside calo
0102
0103   double depthAtShower = 0.0;
0104
0105   // set thePathLengthOnEcal, the pathLength at the reference (r=123.8 for
0106   // barrel and z=304.5 for endcap)
0107
0108   // effective interaction length fitter to ssp from Geant4
0109   double effectiveLambda = 0.0;
0110   if (theEnergy > 0.0 && theEnergy < 15) {
0111     effectiveLambda = 24.6 + 2.6 * std::tanh(3.0 * (std::log(theEnergy) - 1.43));
0112   } else {
0113     effectiveLambda = 28.4 + 1.20 * std::tanh(1.5 * (std::log(theEnergy) - 4.3));
0114   }
0115   // fraction before the crystal, but inside Ecal
0116   //  double frac_ssp1
0117   //  = 1.5196e-01+1.3300e-01*tanh(-4.6971e-01*(std::log(theEnergy)+2.4162e+00));
0118   // fraction after the crystal, but before Hcal
0119   double frac_ssp2 = 2.8310e+00 + 2.6766e+00 * tanh(-4.8068e-01 * (std::log(theEnergy) + 3.4857e+00));
0120
0121   if (fastSimShowerType == 100) {  // fastTrack.onEcal() == 1
0122
0123     //    double rhoTemp = Gflash::ROffCrystalEB +
0124     //    Gflash::LengthCrystalEB*std::sin(position.getTheta());
0125     double rhoTemp = Gflash::LengthCrystalEB * std::sin(position.getTheta());
0126     thePathLengthOnEcal = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
0127     double pathLengthAt2 = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB + rhoTemp);
0128     double pathLengthAt3 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]);
0129
0130     /// fraction before the crystal, but inside Ecal
0131     /*
0132     if(CLHEP::HepUniformRand() < frac_ssp1 ) {
0133       depthAtShower =
0134     (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
0135     }
0136     else {
0137     */
0138     // inside the crystal
0139     //      depthAtShower = (pathLengthAt1-thePathLengthOnEcal) -
0140     //      effectiveLambda*log(CLHEP::HepUniformRand());
0141     depthAtShower = -effectiveLambda * log(CLHEP::HepUniformRand());
0142     // after the crystal
0143     if (depthAtShower > (pathLengthAt2 - thePathLengthOnEcal)) {
0144       // before Hcal
0145       if (CLHEP::HepUniformRand() < frac_ssp2) {
0146         depthAtShower =
0147             (pathLengthAt2 - thePathLengthOnEcal) + (pathLengthAt3 - pathLengthAt2) * CLHEP::HepUniformRand();
0148       }
0149       // inside Hcal
0150       else {
0151         depthAtShower = (pathLengthAt3 - thePathLengthOnEcal) - effectiveLambda * log(CLHEP::HepUniformRand());
0152         // check whether the shower starts beyond HB
0153         double pathLengthAt4 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
0154         if (depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
0155           depthAtShower = (pathLengthAt4 - pathLengthAt3) * CLHEP::HepUniformRand();
0156         }
0157       }
0158     }
0159     //    }
0160   } else if (fastSimShowerType == 101) {  // fastTrack.onEcal() == 2
0161
0162     double zTemp = Gflash::LengthCrystalEE;
0163     double zsign = (position.getEta() > 0) ? 1.0 : -1.0;
0164
0165     thePathLengthOnEcal = theHelix->getPathLengthAtZ(zsign * Gflash::ZFrontCrystalEE);
0166     double pathLengthAt2 = theHelix->getPathLengthAtZ(zsign * (Gflash::ZFrontCrystalEE + zTemp));
0167     double pathLengthAt3 = theHelix->getPathLengthAtZ(zsign * Gflash::Zmin[Gflash::kHE]);
0168
0169     /*
0170     if(CLHEP::HepUniformRand() <  frac_ssp1 ) {
0171       depthAtShower =
0172     (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
0173     }
0174     else {
0175     */
0176     //      depthAtShower =
0177     //      (pathLengthAt1-thePathLengthOnEcal)-effectiveLambda*std::log(CLHEP::HepUniformRand());
0178     depthAtShower = -effectiveLambda * std::log(CLHEP::HepUniformRand());
0179
0180     if (depthAtShower > (pathLengthAt2 - thePathLengthOnEcal)) {
0181       if (CLHEP::HepUniformRand() < frac_ssp2) {
0182         depthAtShower =
0183             (pathLengthAt2 - thePathLengthOnEcal) + (pathLengthAt3 - pathLengthAt2) * CLHEP::HepUniformRand();
0184       } else {
0185         depthAtShower = (pathLengthAt3 - thePathLengthOnEcal) - effectiveLambda * std::log(CLHEP::HepUniformRand());
0186         // shower starts beyond HE
0187         double pathLengthAt4 = theHelix->getPathLengthAtZ(zsign * Gflash::Zmax[Gflash::kHE]);
0188         if (depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
0189           depthAtShower = (pathLengthAt4 - pathLengthAt3) * CLHEP::HepUniformRand();
0190         }
0191       }
0192     }
0193     //    }
0194   } else {
0195     depthAtShower = 0.0;
0196   }
0197
0198   double pathLength = thePathLengthOnEcal + depthAtShower;
0199
0200   theHelix->getGflashTrajectoryPoint(theTrajectoryPoint, pathLength);
0201
0202   // return the initial showino position at the shower starting position
0203   return theTrajectoryPoint.getPosition();
0204 }