CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimGeneral/GFlash/src/GflashShowino.cc

Go to the documentation of this file.
00001 #include "SimGeneral/GFlash/interface/GflashShowino.h"
00002 #include <CLHEP/Random/Randomize.h>
00003 
00004 GflashShowino::GflashShowino() 
00005   : theShowerType(-1), theEnergy(0), thePositionAtShower(0,0,0), thePathLengthAtShower(0), 
00006     thePathLengthOnEcal(0), theStepLengthToHcal(0), theStepLengthToOut(0),  
00007     thePathLength(0), theGlobalTime(0), thePosition(0,0,0), theEnergyDeposited(0)
00008 {
00009   theHelix =  new GflashTrajectory;;
00010 }
00011 
00012 GflashShowino::~GflashShowino() 
00013 {
00014   delete theHelix;
00015 }
00016 
00017 void GflashShowino::initialize(int showerType, double energy, double globalTime, double charge, 
00018                                Gflash3Vector &position, Gflash3Vector &momentum,
00019                                double magneticField) {
00020 
00021   theEnergy = energy;
00022   theGlobalTime = globalTime;
00023   theEnergyDeposited = 0.0;
00024 
00025   // inside the magnetic field (tesla unit);
00026   theHelix->initializeTrajectory(momentum,position,charge,magneticField);
00027 
00028   if(showerType<100) {
00029     thePositionAtShower = position;  
00030     thePosition = thePositionAtShower;
00031     theShowerType = showerType;
00032 
00033   }
00034   else {
00035     //this input is from FastSimulation
00036     //1. simulate the shower starting position
00037     thePositionAtShower = simulateFirstInteractionPoint(showerType,position);
00038     thePosition = thePositionAtShower;
00039 
00040     //2. find shower type depending on where is the shower starting point
00041     theShowerType = Gflash::findShowerType(thePositionAtShower);
00042 
00043   }  
00044 
00045   evaluateLengths();
00046 
00047 }
00048 
00049 void GflashShowino::updateShowino(double deltaStep)
00050 {
00051   thePathLength += deltaStep;
00052   //trajectory point of showino along the shower depth at the pathLength
00053   GflashTrajectoryPoint trajectoryShowino;
00054   theHelix->getGflashTrajectoryPoint(trajectoryShowino,thePathLength);
00055 
00056   thePosition = trajectoryShowino.getPosition();
00057 
00058   theGlobalTime +=  (theEnergy/100.0)*deltaStep/30.0; //@@@calculate exact time change in nsec
00059 }
00060 
00061 void GflashShowino::evaluateLengths() {
00062   //thePathLengthAtShower: path Length from the origin to the shower starting point in cm
00063   //theStepLengthToOut: the total path length from the starting point of
00064   //                    shower to the maximum distance inside paramerized envelopes
00065   double eta = thePosition.getEta();
00066   
00067   if (std::fabs(eta) < Gflash::EtaMax[Gflash::kESPM]) {
00068     thePathLengthOnEcal   = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
00069     thePathLengthAtShower = theHelix->getPathLengthAtRhoEquals(thePosition.getRho());
00070     double pathLengthAtHcalBack = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
00071     if (pathLengthAtHcalBack>0) {
00072       theStepLengthToOut  = std::min(300.,pathLengthAtHcalBack - thePathLengthAtShower);
00073     }
00074     else {
00075       theStepLengthToOut  = 200.;
00076     }
00077     theStepLengthToHcal   = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]) - thePathLengthAtShower;
00078 
00079   }
00080   else if (std::fabs(eta) < Gflash::EtaMax[Gflash::kENCA] ) {
00081     double zsign = (eta > 0) ? 1.0 : -1.0 ; 
00082     thePathLengthOnEcal   = theHelix->getPathLengthAtZ(zsign*Gflash::ZFrontCrystalEE);
00083     thePathLengthAtShower = theHelix->getPathLengthAtZ(thePosition.getZ());
00084     theStepLengthToOut    = std::min(300.,theHelix->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kHE]) - thePathLengthAtShower);
00085     theStepLengthToHcal   = theHelix->getPathLengthAtZ(zsign*Gflash::Zmin[Gflash::kHE]) - thePathLengthAtShower;
00086   }
00087   else { 
00088     //@@@extend for HF later
00089     theStepLengthToOut = 200.0;
00090   }
00091   
00092   thePathLength = thePathLengthAtShower;
00093 
00094 }
00095 
00096 Gflash3Vector& GflashShowino::simulateFirstInteractionPoint(int fastSimShowerType, Gflash3Vector &position) {
00097 
00098   //determine the shower starting point (ssp):
00099   //the position at the entrance + the mean free path till the inelastic interaction inside calo
00100 
00101   double depthAtShower = 0.0;
00102 
00103   //set thePathLengthOnEcal, the pathLength at the reference (r=123.8 for barrel and z=304.5 for endcap)
00104     
00105   //effective interaction length fitter to ssp from Geant4
00106   double effectiveLambda = 0.0;
00107   if(theEnergy > 0.0 &&  theEnergy < 15) {
00108     effectiveLambda = 24.6+2.6*std::tanh(3.0*(std::log(theEnergy)-1.43));
00109   }
00110   else {
00111     effectiveLambda = 28.4+1.20*std::tanh(1.5*(std::log(theEnergy)-4.3));
00112   }
00113   //fraction before the crystal, but inside Ecal
00114   //  double frac_ssp1 = 1.5196e-01+1.3300e-01*tanh(-4.6971e-01*(std::log(theEnergy)+2.4162e+00));
00115   //fraction after the crystal, but before Hcal
00116   double frac_ssp2 = 2.8310e+00+2.6766e+00*tanh(-4.8068e-01*(std::log(theEnergy)+3.4857e+00));
00117 
00118   if(fastSimShowerType == 100 ) { //fastTrack.onEcal() == 1
00119 
00120     //    double rhoTemp = Gflash::ROffCrystalEB + Gflash::LengthCrystalEB*std::sin(position.getTheta());
00121     double rhoTemp = Gflash::LengthCrystalEB*std::sin(position.getTheta());
00122     thePathLengthOnEcal  = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
00123     double pathLengthAt2 = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB + rhoTemp );
00124     double pathLengthAt3 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]);
00125 
00127     /*
00128     if(CLHEP::HepUniformRand() < frac_ssp1 ) {
00129       depthAtShower = (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
00130     }
00131     else {
00132     */
00133     //inside the crystal
00134     //      depthAtShower = (pathLengthAt1-thePathLengthOnEcal) - effectiveLambda*log(CLHEP::HepUniformRand());
00135       depthAtShower = - effectiveLambda*log(CLHEP::HepUniformRand());
00136       //after the crystal
00137       if(depthAtShower > (pathLengthAt2 - thePathLengthOnEcal) ) {
00138         //before Hcal
00139         if(CLHEP::HepUniformRand() < frac_ssp2 ) {
00140           depthAtShower = (pathLengthAt2 - thePathLengthOnEcal) + (pathLengthAt3-pathLengthAt2)*CLHEP::HepUniformRand();
00141         }
00142         //inside Hcal
00143         else {
00144           depthAtShower = (pathLengthAt3 - thePathLengthOnEcal) - effectiveLambda*log(CLHEP::HepUniformRand());
00145           //check whether the shower starts beyond HB
00146           double pathLengthAt4 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
00147           if(depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
00148             depthAtShower = (pathLengthAt4 - pathLengthAt3)*CLHEP::HepUniformRand();
00149           }
00150         }
00151       }
00152       //    }
00153   }
00154   else if(fastSimShowerType == 101 ) { //fastTrack.onEcal() == 2 
00155 
00156     double zTemp = Gflash::LengthCrystalEE;
00157     double zsign = (position.getEta() > 0) ? 1.0 : -1.0 ; 
00158 
00159     thePathLengthOnEcal  = theHelix->getPathLengthAtZ(zsign*Gflash::ZFrontCrystalEE);
00160     double pathLengthAt2 = theHelix->getPathLengthAtZ(zsign*(Gflash::ZFrontCrystalEE + zTemp));
00161     double pathLengthAt3 = theHelix->getPathLengthAtZ(zsign*Gflash::Zmin[Gflash::kHE]);
00162       
00163     /*
00164     if(CLHEP::HepUniformRand() <  frac_ssp1 ) {
00165       depthAtShower = (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
00166     }
00167     else {
00168     */
00169     //      depthAtShower = (pathLengthAt1-thePathLengthOnEcal)-effectiveLambda*std::log(CLHEP::HepUniformRand());
00170       depthAtShower = -effectiveLambda*std::log(CLHEP::HepUniformRand());
00171 
00172       if(depthAtShower > (pathLengthAt2 - thePathLengthOnEcal) ) {
00173         if(CLHEP::HepUniformRand() < frac_ssp2 ) {
00174           depthAtShower = (pathLengthAt2 - thePathLengthOnEcal) +(pathLengthAt3 - pathLengthAt2)*CLHEP::HepUniformRand();
00175         }
00176         else {
00177           depthAtShower = (pathLengthAt3 - thePathLengthOnEcal)-effectiveLambda*std::log(CLHEP::HepUniformRand());
00178           //shower starts beyond HE
00179           double pathLengthAt4 = theHelix->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kHE]);
00180           if(depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
00181             depthAtShower = (pathLengthAt4 - pathLengthAt3)*CLHEP::HepUniformRand();
00182           }
00183         }
00184       }
00185       //    }
00186   }
00187   else {
00188     depthAtShower = 0.0; 
00189   }
00190 
00191   double pathLength  = thePathLengthOnEcal + depthAtShower;
00192   GflashTrajectoryPoint trajectoryPoint;
00193   theHelix->getGflashTrajectoryPoint(trajectoryPoint,pathLength);
00194 
00195   //return the initial showino position at the shower starting position
00196   return trajectoryPoint.getPosition();
00197 
00198 }