GflashShowino.cc

Go to the documentation of this file.

#include "SimGeneral/GFlash/interface/GflashShowino.h"
#include <CLHEP/Random/Randomize.h>

GflashShowino::GflashShowino() 
  : theShowerType(-1), theEnergy(0), thePositionAtShower(0,0,0), thePathLengthAtShower(0), 
    thePathLengthOnEcal(0), theStepLengthToHcal(0), theStepLengthToOut(0),  
    thePathLength(0), theGlobalTime(0), thePosition(0,0,0), theEnergyDeposited(0)
{
  theHelix =  new GflashTrajectory;;
}

GflashShowino::~GflashShowino() 
{
  delete theHelix;
}

void GflashShowino::initialize(int showerType, double energy, double globalTime, double charge, 
                   Gflash3Vector &position, Gflash3Vector &momentum,
                   double magneticField) {

  theEnergy = energy;
  theGlobalTime = globalTime;
  theEnergyDeposited = 0.0;

  // inside the magnetic field (tesla unit);
  theHelix->initializeTrajectory(momentum,position,charge,magneticField);

  if(showerType<100) {
    thePositionAtShower = position;  
    thePosition = thePositionAtShower;
    theShowerType = showerType;

  }
  else {
    //this input is from FastSimulation
    //1. simulate the shower starting position
    thePositionAtShower = simulateFirstInteractionPoint(showerType,position);
    thePosition = thePositionAtShower;

    //2. find shower type depending on where is the shower starting point
    theShowerType = Gflash::findShowerType(thePositionAtShower);

  }  

  evaluateLengths();

}

void GflashShowino::updateShowino(double deltaStep)
{
  thePathLength += deltaStep;
  //trajectory point of showino along the shower depth at the pathLength
  GflashTrajectoryPoint trajectoryShowino;
  theHelix->getGflashTrajectoryPoint(trajectoryShowino,thePathLength);

  thePosition = trajectoryShowino.getPosition();

  theGlobalTime +=  (theEnergy/100.0)*deltaStep/30.0; //@@@calculate exact time change in nsec
}

void GflashShowino::evaluateLengths() {
  //thePathLengthAtShower: path Length from the origin to the shower starting point in cm
  //theStepLengthToOut: the total path length from the starting point of
  //                    shower to the maximum distance inside paramerized envelopes
  double eta = thePosition.getEta();
  
  if (std::fabs(eta) < Gflash::EtaMax[Gflash::kESPM]) {
    thePathLengthOnEcal   = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
    thePathLengthAtShower = theHelix->getPathLengthAtRhoEquals(thePosition.getRho());
    double pathLengthAtHcalBack = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
    if (pathLengthAtHcalBack>0) {
      theStepLengthToOut  = std::min(300.,pathLengthAtHcalBack - thePathLengthAtShower);
    }
    else {
      theStepLengthToOut  = 200.;
    }
    theStepLengthToHcal   = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]) - thePathLengthAtShower;

  }
  else if (std::fabs(eta) < Gflash::EtaMax[Gflash::kENCA] ) {
    double zsign = (eta > 0) ? 1.0 : -1.0 ; 
    thePathLengthOnEcal   = theHelix->getPathLengthAtZ(zsign*Gflash::ZFrontCrystalEE);
    thePathLengthAtShower = theHelix->getPathLengthAtZ(thePosition.getZ());
    theStepLengthToOut    = std::min(300.,theHelix->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kHE]) - thePathLengthAtShower);
    theStepLengthToHcal   = theHelix->getPathLengthAtZ(zsign*Gflash::Zmin[Gflash::kHE]) - thePathLengthAtShower;
  }
  else { 
    //@@@extend for HF later
    theStepLengthToOut = 200.0;
  }
  
  thePathLength = thePathLengthAtShower;

}

Gflash3Vector& GflashShowino::simulateFirstInteractionPoint(int fastSimShowerType, Gflash3Vector &position) {

  //determine the shower starting point (ssp):
  //the position at the entrance + the mean free path till the inelastic interaction inside calo

  double depthAtShower = 0.0;

  //set thePathLengthOnEcal, the pathLength at the reference (r=123.8 for barrel and z=304.5 for endcap)
    
  //effective interaction length fitter to ssp from Geant4
  double effectiveLambda = 0.0;
  if(theEnergy > 0.0 &&  theEnergy < 15) {
    effectiveLambda = 24.6+2.6*std::tanh(3.0*(std::log(theEnergy)-1.43));
  }
  else {
    effectiveLambda = 28.4+1.20*std::tanh(1.5*(std::log(theEnergy)-4.3));
  }
  //fraction before the crystal, but inside Ecal
  //  double frac_ssp1 = 1.5196e-01+1.3300e-01*tanh(-4.6971e-01*(std::log(theEnergy)+2.4162e+00));
  //fraction after the crystal, but before Hcal
  double frac_ssp2 = 2.8310e+00+2.6766e+00*tanh(-4.8068e-01*(std::log(theEnergy)+3.4857e+00));

  if(fastSimShowerType == 100 ) { //fastTrack.onEcal() == 1

    //    double rhoTemp = Gflash::ROffCrystalEB + Gflash::LengthCrystalEB*std::sin(position.getTheta());
    double rhoTemp = Gflash::LengthCrystalEB*std::sin(position.getTheta());
    thePathLengthOnEcal  = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB);
    double pathLengthAt2 = theHelix->getPathLengthAtRhoEquals(Gflash::RFrontCrystalEB + rhoTemp );
    double pathLengthAt3 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]);

    /*
    if(CLHEP::HepUniformRand() < frac_ssp1 ) {
      depthAtShower = (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
    }
    else {
    */
    //inside the crystal
    //      depthAtShower = (pathLengthAt1-thePathLengthOnEcal) - effectiveLambda*log(CLHEP::HepUniformRand());
      depthAtShower = - effectiveLambda*log(CLHEP::HepUniformRand());
      //after the crystal
      if(depthAtShower > (pathLengthAt2 - thePathLengthOnEcal) ) {
    //before Hcal
    if(CLHEP::HepUniformRand() < frac_ssp2 ) {
      depthAtShower = (pathLengthAt2 - thePathLengthOnEcal) + (pathLengthAt3-pathLengthAt2)*CLHEP::HepUniformRand();
    }
    //inside Hcal
    else {
      depthAtShower = (pathLengthAt3 - thePathLengthOnEcal) - effectiveLambda*log(CLHEP::HepUniformRand());
      //check whether the shower starts beyond HB
      double pathLengthAt4 = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
      if(depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
        depthAtShower = (pathLengthAt4 - pathLengthAt3)*CLHEP::HepUniformRand();
      }
    }
      }
      //    }
  }
  else if(fastSimShowerType == 101 ) { //fastTrack.onEcal() == 2 

    double zTemp = Gflash::LengthCrystalEE;
    double zsign = (position.getEta() > 0) ? 1.0 : -1.0 ; 

    thePathLengthOnEcal  = theHelix->getPathLengthAtZ(zsign*Gflash::ZFrontCrystalEE);
    double pathLengthAt2 = theHelix->getPathLengthAtZ(zsign*(Gflash::ZFrontCrystalEE + zTemp));
    double pathLengthAt3 = theHelix->getPathLengthAtZ(zsign*Gflash::Zmin[Gflash::kHE]);
      
    /*
    if(CLHEP::HepUniformRand() <  frac_ssp1 ) {
      depthAtShower = (pathLengthAt1-thePathLengthOnEcal)*CLHEP::HepUniformRand();
    }
    else {
    */
    //      depthAtShower = (pathLengthAt1-thePathLengthOnEcal)-effectiveLambda*std::log(CLHEP::HepUniformRand());
      depthAtShower = -effectiveLambda*std::log(CLHEP::HepUniformRand());

      if(depthAtShower > (pathLengthAt2 - thePathLengthOnEcal) ) {
    if(CLHEP::HepUniformRand() < frac_ssp2 ) {
      depthAtShower = (pathLengthAt2 - thePathLengthOnEcal) +(pathLengthAt3 - pathLengthAt2)*CLHEP::HepUniformRand();
    }
    else {
      depthAtShower = (pathLengthAt3 - thePathLengthOnEcal)-effectiveLambda*std::log(CLHEP::HepUniformRand());
      //shower starts beyond HE
      double pathLengthAt4 = theHelix->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kHE]);
      if(depthAtShower > (pathLengthAt4 - thePathLengthOnEcal)) {
        depthAtShower = (pathLengthAt4 - pathLengthAt3)*CLHEP::HepUniformRand();
      }
    }
      }
      //    }
  }
  else {
    depthAtShower = 0.0; 
  }

  double pathLength  = thePathLengthOnEcal + depthAtShower;
  GflashTrajectoryPoint trajectoryPoint;
  theHelix->getGflashTrajectoryPoint(trajectoryPoint,pathLength);

  //return the initial showino position at the shower starting position
  return trajectoryPoint.getPosition();

}