#include <GflashEMShowerProfile.h>

Public Member Functions
std::vector< GflashHit > &	getGflashHitList ()

GflashShowino *	getGflashShowino ()

	GflashEMShowerProfile (edm::ParameterSet parSet)

void	initialize (int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum)

void	parameterization ()

	~GflashEMShowerProfile ()

Private Member Functions
double	getDistanceToOut (Gflash::CalorimeterNumber kCalor)

Gflash3Vector	locateHitPosition (GflashTrajectoryPoint &point, double rCore, double rTail, double probability, double &rShower)

Private Attributes
Gflash::CalorimeterNumber	jCalorimeter

double	theBField

double	theEnergyScaleEB

double	theEnergyScaleEE

std::vector< GflashHit >	theGflashHitList

GflashHistogram *	theHisto

edm::ParameterSet	theParSet

GflashShowino *	theShowino

Detailed Description

Definition at line 15 of file GflashEMShowerProfile.h.

Constructor & Destructor Documentation

GflashEMShowerProfile::GflashEMShowerProfile ( edm::ParameterSet parSet )

Definition at line 21 of file GflashEMShowerProfile.cc.

References edm::ParameterSet::getParameter(), GflashHistogram::instance(), jCalorimeter, Gflash::kNULL, theBField, theEnergyScaleEB, theEnergyScaleEE, theHisto, and theShowino.

                                                                    : theParSet(parSet)
 {
   theBField = parSet.getParameter<double>("bField");
   theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
   theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");
 
   jCalorimeter = Gflash::kNULL;
 
   theShowino = new GflashShowino();
 
   theHisto = GflashHistogram::instance();
 
 }

GflashEMShowerProfile::~GflashEMShowerProfile ( )

Definition at line 36 of file GflashEMShowerProfile.cc.

References theShowino.

 {
   if(theShowino) delete theShowino;
 }

Member Function Documentation

double GflashEMShowerProfile::getDistanceToOut ( Gflash::CalorimeterNumber kCalor )

private

Definition at line 281 of file GflashEMShowerProfile.cc.

References GflashShowino::getHelix(), GflashTrajectory::getPathLengthAtRhoEquals(), GflashShowino::getPathLengthAtShower(), GflashTrajectory::getPathLengthAtZ(), GflashShowino::getPosition(), Gflash::kENCA, Gflash::kESPM, Gflash::Rmax, theShowino, and Gflash::Zmax.

Referenced by parameterization().

                                                                              {
 
   double stepLengthLeft = 0.0;
   if(kCalor == Gflash::kESPM ) {
     stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kESPM])
                    - theShowino->getPathLengthAtShower();
   }
   else if (kCalor == Gflash::kENCA) {
     double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
     stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kENCA])
                    - theShowino->getPathLengthAtShower();
   }
   return stepLengthLeft;
 
 }

std::vector<GflashHit>& GflashEMShowerProfile::getGflashHitList ( )

inline

Definition at line 30 of file GflashEMShowerProfile.h.

References theGflashHitList.

Referenced by GflashEMShowerModel::makeHits().

30 {return theGflashHitList;};

GflashEMShowerProfile::theGflashHitList

std::vector< GflashHit > theGflashHitList

Definition: GflashEMShowerProfile.h:48

GflashShowino* GflashEMShowerProfile::getGflashShowino ( )

inline

Definition at line 29 of file GflashEMShowerProfile.h.

References theShowino.

29 { return theShowino; }

GflashEMShowerProfile::theShowino

GflashShowino * theShowino

Definition: GflashEMShowerProfile.h:46

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

Definition at line 41 of file GflashEMShowerProfile.cc.

References GflashShowino::initialize(), theBField, and theShowino.

Referenced by GflashEMShowerModel::DoIt().

 {
   theShowino->initialize(showerType, energy, globalTime, charge,
                          position, momentum, theBField);
 }

Gflash3Vector GflashEMShowerProfile::locateHitPosition	(	GflashTrajectoryPoint &	point,
		double	rCore,
		double	rTail,
		double	probability,
		double &	rShower
	)

private

Definition at line 297 of file GflashEMShowerProfile.cc.

References funct::cos(), GflashTrajectoryPoint::getCrossUnitVector(), GflashTrajectoryPoint::getOrthogonalUnitVector(), GflashTrajectoryPoint::getPosition(), jCalorimeter, position, Gflash::rMoliere, funct::sin(), and mathSSE::sqrt().

Referenced by parameterization().

 {
   double u1 = CLHEP::HepUniformRand();
   double u2 = CLHEP::HepUniformRand();
   double rInRM = 0.0;
   
   if (u1 < probability ) {
     rInRM = rCore* std::sqrt( u2/(1.0-u2) );
   }
   else {
     rInRM = rTail * std::sqrt( u2/(1.0-u2) );
   }
   
   rShower =  rInRM * Gflash::rMoliere[jCalorimeter];
 
   // Uniform & random rotation of spot along the azimuthal angle
   double azimuthalAngle = CLHEP::twopi*CLHEP::HepUniformRand();
 
   // actual spot position by adding a radial vector to a trajectoryPoint
   Gflash3Vector position = point.getPosition() +
     rShower*std::cos(azimuthalAngle)*point.getOrthogonalUnitVector() +
     rShower*std::sin(azimuthalAngle)*point.getCrossUnitVector();
   
   return position;
 }

void GflashEMShowerProfile::parameterization ( )

Definition at line 48 of file GflashEMShowerProfile.cc.

Referenced by GflashEMShowerModel::DoIt().

 {
   // This part of code is copied from the original GFlash Fortran code.
   // reference : hep-ex/0001020v1
   // The units used in Geant4 internally are in mm, MeV.
   // For simplicity, the units here are in cm, GeV.
 
   const double energyCutoff     = 0.01; 
   const int    maxNumberOfSpots = 100000;
 
   double incomingEnergy   = theShowino->getEnergy();
   Gflash3Vector showerStartingPosition = theShowino->getPositionAtShower();
 
   //find the calorimeter at the shower starting point
   jCalorimeter = Gflash::getCalorimeterNumber(showerStartingPosition);
 
   double logEinc = std::log(incomingEnergy);
   double y = incomingEnergy / Gflash::criticalEnergy; // y = E/Ec, criticalEnergy is in GeV
   double logY = std::log(y);
 
   // Total number of spots are not yet optimized.
   double nSpots = 93.0 * std::log(Gflash::Z[jCalorimeter]) * std::pow(incomingEnergy,0.876);
 
   //path Length from the origin to the shower starting point in cm
   double pathLength0 = theShowino->getPathLengthAtShower();
   double pathLength = pathLength0; // this will grow along the shower development
 
   //--- 2.2  Fix intrinsic properties of em. showers.
 
   double fluctuatedTmax = std::log(logY - 0.7157);
   double fluctuatedAlpha = std::log(0.7996 +(0.4581 + 1.8628/Gflash::Z[jCalorimeter])*logY);
 
   double sigmaTmax = 1.0/( -1.4  + 1.26 * logY);
   double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
   double rho = 0.705  - 0.023 * logY;
   double sqrtPL = std::sqrt((1.0+rho)/2.0);
   double sqrtLE = std::sqrt((1.0-rho)/2.0);
 
   double norm1 =  CLHEP::RandGaussQ::shoot();
   double norm2 =  CLHEP::RandGaussQ::shoot();
   double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
   double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
 
   // tmax, alpha, beta : parameters of gamma distribution
   double tmax = std::exp(tempTmax);
   double alpha = std::exp(tempAlpha);
   double beta = std::max(0.0,(alpha - 1.0)/tmax);
  
   // spot fluctuations are added to tmax, alpha, beta
   double averageTmax = logY-0.858;
   double averageAlpha = 0.21+(0.492+2.38/Gflash::Z[jCalorimeter])*logY;
   double spotTmax  = averageTmax * (0.698 + .00212*Gflash::Z[jCalorimeter]);
   double spotAlpha = averageAlpha * (0.639 + .00334*Gflash::Z[jCalorimeter]);
   double spotBeta = std::max(0.0,(spotAlpha-1.0)/spotTmax);
 
    if(theHisto->getStoreFlag()) {
     theHisto->em_incE->Fill(incomingEnergy);
     theHisto->em_ssp_rho->Fill(showerStartingPosition.rho());
     theHisto->em_ssp_z->Fill(showerStartingPosition.z());
   }
 
   //  parameters for lateral distribution and fluctuation
   double z1=0.0251+0.00319*logEinc;
   double z2=0.1162-0.000381*Gflash::Z[jCalorimeter];
 
   double k1=0.659 - 0.00309 * Gflash::Z[jCalorimeter];
   double k2=0.645;
   double k3=-2.59;
   double k4=0.3585+ 0.0421*logEinc;
 
   double p1=2.623 -0.00094*Gflash::Z[jCalorimeter];
   double p2=0.401 +0.00187*Gflash::Z[jCalorimeter];
   double p3=1.313 -0.0686*logEinc;
 
   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : e25Scale = 1.006 for EB with ecalNotContainment = 1.0.
   // Now e25Scale is a configurable parameter with default ecalNotContainment which is 0.97 for EB and 0.975 for EE.
   // So if ecalNotContainment constants are to be changed in the future, e25Scale should be changed accordingly.
   double e25Scale = 1.0;
   if(jCalorimeter == Gflash::kESPM) e25Scale = theEnergyScaleEB;
   else if(jCalorimeter == Gflash::kENCA) e25Scale = theEnergyScaleEE;
 
   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : p1 *= 0.965
   p1 *= 0.965;
 
   // preparation of longitudinal integration
   double stepLengthLeft = getDistanceToOut(jCalorimeter);
 
   int    nSpots_sd = 0; // count total number of spots in SD
   double zInX0 = 0.0; // shower depth in X0 unit
   double deltaZInX0 = 0.0; // segment of depth in X0 unit
   double deltaZ = 0.0; // segment of depth in cm
   double stepLengthLeftInX0 = 0.0; // step length left in X0 unit
 
   const double divisionStepInX0 = 0.1; //step size in X0 unit
   double energy = incomingEnergy; // energy left in GeV
 
   Genfun::IncompleteGamma gammaDist;
 
   double energyInGamma = 0.0; // energy in a specific depth(z) according to Gamma distribution
   double preEnergyInGamma = 0.0; // energy calculated in a previous depth
   double sigmaInGamma  = 0.; // sigma of energy in a specific depth(z) according to Gamma distribution
   double preSigmaInGamma = 0.0; // sigma of energy in a previous depth
 
   //energy segment in Gamma distribution of shower in each step  
   double deltaEnergy =0.0 ; // energy in deltaZ
   int spotCounter = 0; // keep track of number of spots generated
 
   //step increment along the shower direction
   double deltaStep = 0.0;
 
   theGflashHitList.clear();
 
   // loop for longitudinal integration
   while(energy > 0.0 && stepLengthLeft > 0.0) { 
 
     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
 
     if ( stepLengthLeftInX0 < divisionStepInX0 ) {
       deltaZInX0 = stepLengthLeftInX0;
       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
       stepLengthLeft = 0.0;
     }
     else {
       deltaZInX0 = divisionStepInX0;
       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
       stepLengthLeft -= deltaZ;
     }
 
     zInX0 += deltaZInX0;
 
     int nSpotsInStep = 0;
 
     if ( energy > energyCutoff  ) {
       preEnergyInGamma  = energyInGamma;
       gammaDist.a().setValue(alpha);  //alpha
       energyInGamma = gammaDist(beta*zInX0); //beta
       double energyInDeltaZ  = energyInGamma - preEnergyInGamma;
       deltaEnergy   = std::min(energy,incomingEnergy*energyInDeltaZ);
  
       preSigmaInGamma  = sigmaInGamma;
       gammaDist.a().setValue(spotAlpha);  //alpha spot
       sigmaInGamma = gammaDist(spotBeta*zInX0); //beta spot
       nSpotsInStep = std::max(1,int(nSpots * (sigmaInGamma - preSigmaInGamma)));
     }
     else {
       deltaEnergy = energy;
       preSigmaInGamma  = sigmaInGamma;
       nSpotsInStep = std::max(1,int(nSpots * (1.0 - preSigmaInGamma)));
     }
 
     if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy = energy;
 
     energy  -= deltaEnergy;
 
     if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
       nSpotsInStep = maxNumberOfSpots - spotCounter;
       if ( nSpotsInStep < 1 ) { // @@ check
     edm::LogInfo("SimGeneralGFlash") << "GflashEMShowerProfile: Too Many Spots ";
         edm::LogInfo("SimGeneralGFlash") << " - break to regenerate nSpotsInStep ";
         break;
       }
     }
 
     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
     deltaStep  += 0.5*deltaZ;
     pathLength += deltaStep;
     deltaStep   =  0.5*deltaZ;
 
     // lateral shape and fluctuations for  homogenous calo.
     double tScale = tmax *alpha/(alpha-1.0) * (std::exp(fluctuatedAlpha)-1.0)/std::exp(fluctuatedAlpha);
     double tau = std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
     double rCore = z1 + z2 * tau; 
     double rTail = k1 *( std::exp(k3*(tau-k2)) + std::exp(k4*(tau-k2))); // @@ check RT3 sign
     double p23 = (p2 - tau)/p3;
     double probabilityWeight = p1 *  std::exp( p23 - std::exp(p23) );
 
     // Deposition of spots according to lateral distr.
     // Apply absolute energy scale
     // Convert into MeV unit
     double hitEnergy = deltaEnergy / nSpotsInStep * e25Scale * CLHEP::GeV;
     double hitTime = theShowino->getGlobalTime()*CLHEP::nanosecond + (pathLength - pathLength0)/30.0;
 
     GflashHit aHit;
 
     for (int ispot = 0 ;  ispot < nSpotsInStep ; ispot++) {
       spotCounter++;
 
       // Compute global position of generated spots with taking into account magnetic field
       // Divide deltaZ into nSpotsInStep and give a spot a global position
       double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
 
       // trajectoryPoint give a spot an imaginary point along the shower development
       GflashTrajectoryPoint trajectoryPoint;
       theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint,pathLength+incrementPath);
 
       double rShower = 0.0;
       Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint,rCore,rTail,probabilityWeight,rShower);
 
       // Convert into mm unit
       hitPosition *= CLHEP::cm;
 
       if( std::fabs(hitPosition.getZ()/CLHEP::cm) > Gflash::Zmax[Gflash::kHE]) continue;
 
       // put energy and position to a Hit
       aHit.setTime(hitTime);
       aHit.setEnergy(hitEnergy);
       aHit.setPosition(hitPosition);
       theGflashHitList.push_back(aHit);
 
       double zInX0_spot = std::abs(pathLength+incrementPath - pathLength0)/Gflash::radLength[jCalorimeter];
 
       nSpots_sd++;
 
       // for histogramming      
       if(theHisto->getStoreFlag()) {
     theHisto->em_long->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
     theHisto->em_lateral->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
     theHisto->em_long_sd->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
     theHisto->em_lateral_sd->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
       }
       
     } // end of for spot iteration
 
   } // end of while for longitudinal integration
 
   if(theHisto->getStoreFlag()) {
     theHisto->em_nSpots_sd->Fill(nSpots_sd);
   }
 
   //  delete theGflashNavigator;
 
 }