#include <GflashEMShowerProfile.h>

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

GflashShowino *	getGflashShowino ()

	GflashEMShowerProfile (const 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 ( const edm::ParameterSet & parSet )

Definition at line 20 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 32 of file GflashEMShowerProfile.cc.

References theShowino.

                                               {
   if (theShowino)
     delete theShowino;
 }

Member Function Documentation

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

private

Definition at line 279 of file GflashEMShowerProfile.cc.

References mkfit::getEta(), 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 32 of file GflashEMShowerProfile.h.

References theGflashHitList.

Referenced by GflashEMShowerModel::makeHits(), and GFlashEMShowerModel::makeHits().

32 { return theGflashHitList; };

GflashEMShowerProfile::theGflashHitList

std::vector< GflashHit > theGflashHitList

Definition: GflashEMShowerProfile.h:49

GflashShowino* GflashEMShowerProfile::getGflashShowino ( )

inline

Definition at line 31 of file GflashEMShowerProfile.h.

References theShowino.

31 { return theShowino; }

GflashEMShowerProfile::theShowino

GflashShowino * theShowino

Definition: GflashEMShowerProfile.h:47

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

Definition at line 37 of file GflashEMShowerProfile.cc.

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

Referenced by GflashEMShowerModel::DoIt(), and 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 292 of file GflashEMShowerProfile.cc.

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

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 42 of file GflashEMShowerProfile.cc.

Referenced by GflashEMShowerModel::DoIt(), and 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;
 }