#include <GflashHadronShowerProfile.h>

Inheritance diagram for GflashHadronShowerProfile:

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

GflashShowino *	getGflashShowino ()

	GflashHadronShowerProfile (const edm::ParameterSet &parSet)

void	hadronicParameterization ()

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

virtual void	loadParameters ()

virtual	~GflashHadronShowerProfile ()

Protected Member Functions
double	depthScale (double ssp, double ssp0, double length)

void	doCholeskyReduction (double cc, double vv, const int ndim)

double	fLnE1 (double einc, const double *par)

double	fTanh (double einc, const double *par)

double	gammaProfile (double alpha, double beta, double depth, double lengthUnit)

void	getFluctuationVector (double lowTriangle, double correlationVector)

int	getNumberOfSpots (Gflash::CalorimeterNumber kCalor)

double	hoProfile (double pathLength, double refDepth)

Gflash3Vector	locateHitPosition (GflashTrajectoryPoint &point, double lateralArm)

double	longitudinalProfile ()

double	medianLateralArm (double depth, Gflash::CalorimeterNumber kCalor)

void	setEnergyScale (double einc, const Gflash3Vector &ssp)

double	twoGammaProfile (double *par, double depth, Gflash::CalorimeterNumber kIndex)

Protected Attributes
double	averageSpotEnergy [Gflash::kNumberCalorimeter]

double	energyScale [Gflash::kNumberCalorimeter]

double	lateralPar [Gflash::kNumberCalorimeter][Gflash::Nrpar]

double	longEcal [Gflash::NPar]

double	longHcal [Gflash::NPar]

double	theBField

bool	theGflashHcalOuter

std::vector< GflashHit >	theGflashHitList

GflashHistogram *	theHisto

edm::ParameterSet	theParSet

GflashShowino *	theShowino

Detailed Description

Definition at line 15 of file GflashHadronShowerProfile.h.

Constructor & Destructor Documentation

◆ GflashHadronShowerProfile()

GflashHadronShowerProfile::GflashHadronShowerProfile ( const edm::ParameterSet & parSet )

Definition at line 17 of file GflashHadronShowerProfile.cc.

                                                                                   : theParSet(parSet) {
   theBField = parSet.getParameter<double>("bField");
   theGflashHcalOuter = parSet.getParameter<bool>("GflashHcalOuter");
  
   theShowino = new GflashShowino();
   theHisto = GflashHistogram::instance();
 }

References edm::ParameterSet::getParameter(), GflashHistogram::instance(), theBField, theGflashHcalOuter, theHisto, and theShowino.

◆ ~GflashHadronShowerProfile()

GflashHadronShowerProfile::~GflashHadronShowerProfile ( )

virtual

Definition at line 25 of file GflashHadronShowerProfile.cc.

                                                       {
   if (theShowino)
     delete theShowino;
 }

References theShowino.

Member Function Documentation

◆ depthScale()

double GflashHadronShowerProfile::depthScale	(	double	ssp,
		double	ssp0,
		double	length
	)

protected

Definition at line 487 of file GflashHadronShowerProfile.cc.

                                                                                    {
   double func = 0.0;
   if (length > 0.0)
     func = std::pow((ssp - ssp0) / length, 2.0);
   return func;
 }

References TrackCollections2monitor_cff::func, and funct::pow().

Referenced by GflashAntiProtonShowerProfile::loadParameters(), GflashProtonShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), and GflashKaonMinusShowerProfile::loadParameters().

◆ doCholeskyReduction()

void GflashHadronShowerProfile::doCholeskyReduction	(	double **	cc,
		double **	vv,
		const int	ndim
	)

protected

Definition at line 396 of file GflashHadronShowerProfile.cc.

                                                                                             {
   double sumCjkSquare;
   double vjjLess;
   double sumCikjk;
  
   cc[0][0] = std::sqrt(vv[0][0]);
  
   for (int j = 1; j < ndim; j++) {
     cc[j][0] = vv[j][0] / cc[0][0];
   }
  
   for (int j = 1; j < ndim; j++) {
     sumCjkSquare = 0.0;
     for (int k = 0; k < j; k++)
       sumCjkSquare += cc[j][k] * cc[j][k];
  
     vjjLess = vv[j][j] - sumCjkSquare;
  
     // check for the case that vjjLess is negative
     cc[j][j] = std::sqrt(std::fabs(vjjLess));
  
     for (int i = j + 1; i < ndim; i++) {
       sumCikjk = 0.;
       for (int k = 0; k < j; k++)
         sumCikjk += cc[i][k] * cc[j][k];
       cc[i][j] = (vv[i][j] - sumCikjk) / cc[j][j];
     }
   }
 }

References mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, and mathSSE::sqrt().

Referenced by getFluctuationVector().

◆ fLnE1()

double GflashHadronShowerProfile::fLnE1	(	double	einc,
		const double *	par
	)

protected

Definition at line 480 of file GflashHadronShowerProfile.cc.

                                                                       {
   double func = 0.0;
   if (einc > 0.0)
     func = par[0] + par[1] * std::log(einc);
   return func;
 }

References TrackCollections2monitor_cff::func, and dqm-mbProfile::log.

◆ fTanh()

double GflashHadronShowerProfile::fTanh	(	double	einc,
		const double *	par
	)

protected

Definition at line 473 of file GflashHadronShowerProfile.cc.

                                                                       {
   double func = 0.0;
   if (einc > 0.0)
     func = par[0] + par[1] * std::tanh(par[2] * (std::log(einc) - par[3])) + par[4] * std::log(einc);
   return func;
 }

References TrackCollections2monitor_cff::func, and dqm-mbProfile::log.

Referenced by hadronicParameterization(), GflashProtonShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), GflashKaonMinusShowerProfile::loadParameters(), and GflashAntiProtonShowerProfile::loadParameters().

◆ gammaProfile()

double GflashHadronShowerProfile::gammaProfile	(	double	alpha,
		double	beta,
		double	depth,
		double	lengthUnit
	)

protected

Definition at line 494 of file GflashHadronShowerProfile.cc.

                                                                                                                {
   double gamma = 0.0;
   //  if(alpha > 0 && beta > 0 && lengthUnit > 0) {
   if (showerDepth > 0.0) {
     Genfun::LogGamma lgam;
     double x = showerDepth * (beta / lengthUnit);
     gamma = (beta / lengthUnit) * std::pow(x, alpha - 1.0) * std::exp(-x) / std::exp(lgam(alpha));
   }
   return gamma;
 }

References alpha, HLT_FULL_cff::beta, JetChargeProducer_cfi::exp, CustomPhysics_cfi::gamma, funct::pow(), and x.

Referenced by twoGammaProfile().

◆ getFluctuationVector()

void GflashHadronShowerProfile::getFluctuationVector	(	double *	lowTriangle,
		double *	correlationVector
	)

protected

Definition at line 358 of file GflashHadronShowerProfile.cc.

                                                                                                    {
   const int dim = Gflash::NPar;
  
   double **xr = new double *[dim];
   double **xrho = new double *[dim];
  
   for (int j = 0; j < dim; j++) {
     xr[j] = new double[dim];
     xrho[j] = new double[dim];
   }
  
   for (int i = 0; i < dim; i++) {
     for (int j = 0; j < i + 1; j++) {
       if (j == i)
         xrho[i][j] = 1.0;
       else {
         xrho[i][j] = lowTriangle[i * (i - 1) / 2 + j];
         xrho[j][i] = xrho[i][j];
       }
     }
   }
  
   doCholeskyReduction(xrho, xr, dim);
  
   for (int i = 0; i < dim; i++) {
     for (int j = 0; j < i + 1; j++) {
       correlationVector[i * (i + 1) / 2 + j] = xr[i][j];
     }
   }
  
   for (int j = 0; j < dim; j++)
     delete[] xr[j];
   delete[] xr;
   for (int j = 0; j < dim; j++)
     delete[] xrho[j];
   delete[] xrho;
 }

References doCholeskyReduction(), mps_fire::i, dqmiolumiharvest::j, and Gflash::NPar.

Referenced by GflashAntiProtonShowerProfile::loadParameters(), GflashProtonShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), and GflashKaonMinusShowerProfile::loadParameters().

◆ getGflashHitList()

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

inline

Definition at line 33 of file GflashHadronShowerProfile.h.

33 { return theGflashHitList; };

References theGflashHitList.

Referenced by CalorimetryManager::HDShowerSimulation(), GflashHadronShowerModel::makeHits(), and GFlashHadronShowerModel::makeHits().

◆ getGflashShowino()

GflashShowino* GflashHadronShowerProfile::getGflashShowino ( )

inline

Definition at line 32 of file GflashHadronShowerProfile.h.

32 { return theShowino; }

References theShowino.

Referenced by CalorimetryManager::HDShowerSimulation().

◆ getNumberOfSpots()

int GflashHadronShowerProfile::getNumberOfSpots ( Gflash::CalorimeterNumber kCalor )

protected

Definition at line 426 of file GflashHadronShowerProfile.cc.

                                                                               {
   // generator number of spots: energy dependent Gamma distribution of Nspots
   // based on Geant4 replacing old parameterization of H1, int numberOfSpots =
   // std::max( 50, static_cast<int>(80.*std::log(einc)+50.));
  
   double einc = theShowino->getEnergy();
   int showerType = theShowino->getShowerType();
  
   int numberOfSpots = 0;
   double nmean = 0.0;
   double nsigma = 0.0;
  
   if (showerType == 0 || showerType == 1 || showerType == 4 || showerType == 5) {
     if (kCalor == Gflash::kESPM || kCalor == Gflash::kENCA) {
       nmean = 10000 + 5000 * log(einc);
       nsigma = 1000;
     }
     if (kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
       nmean = 5000 + 2500 * log(einc);
       nsigma = 500;
     }
   } else if (showerType == 2 || showerType == 3 || showerType == 6 || showerType == 7) {
     if (kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
       nmean = 5000 + 2500 * log(einc);
       nsigma = 500;
     } else {
       nmean = 10000;
       nsigma = 1000;
     }
   }
   //@@@need correlation and individual fluctuation on alphaNspots and betaNspots
   // here: evaluating covariance should be straight forward since the
   // distribution is 'one' Gamma
  
   numberOfSpots = std::max(500, static_cast<int>(nmean + nsigma * CLHEP::RandGaussQ::shoot()));
  
   // until we optimize the reduction scale in the number of Nspots
  
   if (kCalor == Gflash::kESPM || kCalor == Gflash::kENCA) {
     numberOfSpots = static_cast<int>(numberOfSpots / 100);
   } else {
     numberOfSpots = static_cast<int>(numberOfSpots / 3.0);
   }
  
   return numberOfSpots;
 }

References GflashShowino::getEnergy(), GflashShowino::getShowerType(), Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, dqm-mbProfile::log, SiStripPI::max, and theShowino.

Referenced by hadronicParameterization().

◆ hadronicParameterization()

void GflashHadronShowerProfile::hadronicParameterization ( )

Definition at line 36 of file GflashHadronShowerProfile.cc.

                                                          {
   // The skeleton of this method is based on the fortran code gfshow.F
   // originally written by S. Peters and G. Grindhammer (also see NIM A290
   // (1990) 469-488), but longitudinal parameterizations of hadron showers are
   // significantly modified for the CMS calorimeter
  
   // unit convention: energy in [GeV] and length in [cm]
   // intrinsic properties of hadronic showers (lateral shower profile)
  
   // The step size of showino along the helix trajectory in cm unit
   double showerDepthR50 = 0.0;
   bool firstHcalHit = true;
  
   // initial valuses that will be changed as the shower developes
   double stepLengthLeft = theShowino->getStepLengthToOut();
  
   double deltaStep = 0.0;
   double showerDepth = 0.0;
  
   Gflash::CalorimeterNumber whichCalor = Gflash::kNULL;
  
   theGflashHitList.clear();
  
   GflashHit aHit;
  
   while (stepLengthLeft > 0.0) {
     // update shower depth and stepLengthLeft
     if (stepLengthLeft < Gflash::divisionStep) {
       deltaStep = stepLengthLeft;
       stepLengthLeft = 0.0;
     } else {
       deltaStep = Gflash::divisionStep;
       stepLengthLeft -= deltaStep;
     }
  
     showerDepth += deltaStep;
     showerDepthR50 += deltaStep;
  
     // update GflashShowino
     theShowino->updateShowino(deltaStep);
  
     // evaluate energy in this deltaStep along the longitudinal shower profile
     double heightProfile = 0.;
     double deltaEnergy = 0.;
  
     whichCalor = Gflash::getCalorimeterNumber(theShowino->getPosition());
  
     // skip if Showino is outside envelopes
     if (whichCalor == Gflash::kNULL)
       continue;
  
     heightProfile = longitudinalProfile();
  
     // skip if the delta energy for this step will be very small
     if (heightProfile < 1.00e-08)
       continue;
  
     // get energy deposition for this step
     deltaEnergy = heightProfile * Gflash::divisionStep * energyScale[whichCalor];
     theShowino->addEnergyDeposited(deltaEnergy);
  
     // apply sampling fluctuation if showino is inside the sampling calorimeter
     double hadronicFraction = 1.0;
     double fluctuatedEnergy = deltaEnergy;
  
     int nSpotsInStep =
         std::max(1, static_cast<int>(getNumberOfSpots(whichCalor) * (deltaEnergy / energyScale[whichCalor])));
     double sampleSpotEnergy = hadronicFraction * fluctuatedEnergy / nSpotsInStep;
  
     // Lateral shape and fluctuations
  
     if ((whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) && firstHcalHit) {
       firstHcalHit = false;
       // reset the showerDepth used in the lateral parameterization inside Hcal
       showerDepthR50 = Gflash::divisionStep;
     }
  
     // evaluate the fluctuated median of the lateral distribution, R50
     double R50 = medianLateralArm(showerDepthR50, whichCalor);
  
     double hitEnergy = sampleSpotEnergy * CLHEP::GeV;
     double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond;
  
     Gflash::CalorimeterNumber hitCalor = Gflash::kNULL;
  
     for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
       // Compute global position of generated spots with taking into account
       // magnetic field Divide deltaStep into nSpotsInStep and give a spot a
       // global position
       double incrementPath =
           theShowino->getPathLength() + (deltaStep / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
  
       // trajectoryPoint give a spot an imaginary point along the shower
       // development
       GflashTrajectoryPoint trajectoryPoint;
       theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, incrementPath);
  
       Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, R50);
  
       hitCalor = Gflash::getCalorimeterNumber(hitPosition);
  
       if (hitCalor == Gflash::kNULL)
         continue;
  
       hitPosition *= CLHEP::cm;
  
       aHit.setTime(hitTime);
       aHit.setEnergy(hitEnergy);
       aHit.setPosition(hitPosition);
       theGflashHitList.push_back(aHit);
  
     }  // end of for spot iteration
   }    // end of while for longitudinal integration
  
   // HO parameterization
  
   if (theShowino->getEnergy() > Gflash::MinEnergyCutOffForHO &&
       fabs(theShowino->getPositionAtShower().pseudoRapidity()) < Gflash::EtaMax[Gflash::kHO] && theGflashHcalOuter) {
     // non zero ho fraction to simulate based on geant4
     double nonzeroProb = 0.7 * fTanh(theShowino->getEnergy(), Gflash::ho_nonzero);
     double r0 = CLHEP::HepUniformRand();
     double leftoverE = theShowino->getEnergy() - theShowino->getEnergyDeposited();
  
     //@@@ nonzeroProb is not random - need further correlation for non-zero HO
     // energy
  
     if (r0 < nonzeroProb && leftoverE > 0.0) {
       // starting path Length and stepLengthLeft
       double pathLength = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHO] - 10);
       stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO] + 10) - pathLength;
       showerDepth = pathLength - theShowino->getPathLengthAtShower();
  
       theShowino->setPathLength(pathLength);
  
       double pathLengthx = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO]);
       double pathLengthy = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
  
       while (stepLengthLeft > 0.0) {
         // update shower depth and stepLengthLeft
         if (stepLengthLeft < Gflash::divisionStep) {
           deltaStep = stepLengthLeft;
           stepLengthLeft = 0.0;
         } else {
           deltaStep = Gflash::divisionStep;
           stepLengthLeft -= deltaStep;
         }
  
         showerDepth += deltaStep;
  
         // update GflashShowino
         theShowino->updateShowino(deltaStep);
  
         // evaluate energy in this deltaStep along the longitudinal shower
         // profile
         double heightProfile = 0.;
         double deltaEnergy = 0.;
  
         double hoScale = leftoverE * (pathLengthx - pathLengthy) / (pathLengthx - theShowino->getPathLengthAtShower());
         double refDepth =
             theShowino->getPathLength() - theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
  
         if (refDepth > 0) {
           heightProfile = hoProfile(theShowino->getPathLength(), refDepth);
           deltaEnergy = heightProfile * Gflash::divisionStep * hoScale;
         }
  
         int nSpotsInStep = std::max(
             50, static_cast<int>((160. + 40 * CLHEP::RandGaussQ::shoot()) * std::log(theShowino->getEnergy()) + 50.));
  
         double hoFraction = 1.00;
         double poissonProb = CLHEP::RandPoissonQ::shoot(1.0);
  
         double fluctuatedEnergy = deltaEnergy * poissonProb;
         double sampleSpotEnergy = hoFraction * fluctuatedEnergy / nSpotsInStep;
  
         // Lateral shape and fluctuations
  
         // evaluate the fluctuated median of the lateral distribution, R50
         double R50 = medianLateralArm(showerDepth, Gflash::kHB);
  
         double hitEnergy = sampleSpotEnergy * CLHEP::GeV;
         double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond;
  
         for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
           double incrementPath =
               theShowino->getPathLength() + (deltaStep / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
  
           // trajectoryPoint give a spot an imaginary point along the shower
           // development
           GflashTrajectoryPoint trajectoryPoint;
           theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, incrementPath);
  
           Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, R50);
           hitPosition *= CLHEP::cm;
  
           if (std::fabs(hitPosition.getZ() / CLHEP::cm) > Gflash::Zmax[Gflash::kHO])
             continue;
  
           aHit.setTime(hitTime);
           aHit.setEnergy(hitEnergy);
           aHit.setPosition(hitPosition);
           theGflashHitList.push_back(aHit);
  
         }  // end of for HO spot iteration
       }    // end of while for HO longitudinal integration
     }
   }
  
   //  delete theGflashNavigator;
 }

Referenced by GflashHadronShowerModel::DoIt(), GFlashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

◆ hoProfile()

double GflashHadronShowerProfile::hoProfile	(	double	pathLength,
		double	refDepth
	)

protected

Definition at line 346 of file GflashHadronShowerProfile.cc.

                                                                               {
   double heightProfile = 0;
  
   GflashTrajectoryPoint tempPoint;
   theShowino->getHelix()->getGflashTrajectoryPoint(tempPoint, pathLength);
  
   double dint = 1.4 * Gflash::intLength[Gflash::kHO] * std::sin(tempPoint.getPosition().getTheta());
   heightProfile = std::exp(-1.0 * refDepth / dint);
  
   return heightProfile;
 }

References JetChargeProducer_cfi::exp, GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getHelix(), GflashTrajectoryPoint::getPosition(), Gflash::intLength, Gflash::kHO, funct::sin(), and theShowino.

Referenced by hadronicParameterization().

◆ initialize()

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

Definition at line 30 of file GflashHadronShowerProfile.cc.

                                                                                                                        {
   // initialize GflashShowino for this track
   theShowino->initialize(showerType, energy, globalTime, charge, position, momentum, theBField);
 }

References ALCARECOTkAlJpsiMuMu_cff::charge, HCALHighEnergyHPDFilter_cfi::energy, GflashShowino::initialize(), position, theBField, and theShowino.

Referenced by GflashHadronShowerModel::DoIt(), GFlashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

◆ loadParameters()

void GflashHadronShowerProfile::loadParameters ( )

virtual

Reimplemented in GflashAntiProtonShowerProfile, GflashKaonMinusShowerProfile, GflashKaonPlusShowerProfile, GflashPiKShowerProfile, and GflashProtonShowerProfile.

Definition at line 247 of file GflashHadronShowerProfile.cc.

                                                {
   edm::LogInfo("SimGeneralGFlash") << "GflashHadronShowerProfile::loadParameters() "
                                    << "should be implimented for each particle type";
 }

Referenced by GflashHadronShowerModel::DoIt(), GFlashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

◆ locateHitPosition()

Gflash3Vector GflashHadronShowerProfile::locateHitPosition	(	GflashTrajectoryPoint &	point,
		double	lateralArm
	)

protected

Definition at line 272 of file GflashHadronShowerProfile.cc.

                                                                                                           {
   // Smearing in r according to f(r)= 2.*r*R50**2/(r**2+R50**2)**2
   double rnunif = CLHEP::HepUniformRand();
   double rxPDF = std::sqrt(rnunif / (1. - rnunif));
   double rShower = lateralArm * rxPDF;
  
   // rShower within maxLateralArmforR50
   rShower = std::min(Gflash::maxLateralArmforR50, rShower);
  
   // Uniform smearing in phi
   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();
  
   //@@@debugging histograms
   if (theHisto->getStoreFlag()) {
     theHisto->rshower->Fill(rShower);
     theHisto->lateralx->Fill(rShower * std::cos(azimuthalAngle));
     theHisto->lateraly->Fill(rShower * std::sin(azimuthalAngle));
   }
   return position;
 }

References funct::cos(), GflashHistogram::getStoreFlag(), GflashHistogram::lateralx, GflashHistogram::lateraly, Gflash::maxLateralArmforR50, min(), point, position, GflashHistogram::rshower, funct::sin(), mathSSE::sqrt(), and theHisto.

Referenced by hadronicParameterization().

◆ longitudinalProfile()

double GflashHadronShowerProfile::longitudinalProfile ( )

protected

Definition at line 299 of file GflashHadronShowerProfile.cc.

                                                       {
   double heightProfile = 0.0;
  
   Gflash3Vector pos = theShowino->getPosition();
   int showerType = theShowino->getShowerType();
   double showerDepth = theShowino->getDepth();
   double transDepth = theShowino->getStepLengthToHcal();
  
   // Energy in a delta step (dz) = (energy to
   // deposite)*[Gamma(z+dz)-Gamma(z)]*dz where the incomplete Gamma function
   // gives an intergrate probability of the longitudinal shower up to the shower
   // depth (z). Instead, we use approximated energy; energy in dz = (energy to
   // deposite)*gamma(z)*dz where gamma is the Gamma-distributed probability
   // function
  
   Gflash::CalorimeterNumber whichCalor = Gflash::getCalorimeterNumber(pos);
  
   if (showerType == 0 || showerType == 4) {
     double shiftDepth = theShowino->getPathLengthOnEcal() - theShowino->getPathLengthAtShower();
     if (shiftDepth > 0) {
       heightProfile = twoGammaProfile(longEcal, showerDepth - shiftDepth, whichCalor);
     } else {
       heightProfile = 0.;
       //      std::cout << "negative shiftDepth for showerType 0 " << shiftDepth
       //      << std::endl;
     }
   } else if (showerType == 1 || showerType == 5) {
     if (whichCalor == Gflash::kESPM || whichCalor == Gflash::kENCA) {
       heightProfile = twoGammaProfile(longEcal, showerDepth, whichCalor);
     } else if (whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) {
       heightProfile = twoGammaProfile(longHcal, showerDepth - transDepth, whichCalor);
     } else
       heightProfile = 0.;
   } else if (showerType == 2 || showerType == 6) {
     // two gammas between crystal and Hcal
     if ((showerDepth - transDepth) > 0.0) {
       heightProfile = twoGammaProfile(longHcal, showerDepth - transDepth, Gflash::kHB);
     } else
       heightProfile = 0.;
   } else if (showerType == 3 || showerType == 7) {
     // two gammas inside Hcal
     heightProfile = twoGammaProfile(longHcal, showerDepth, Gflash::kHB);
   }
  
   return heightProfile;
 }

References Gflash::getCalorimeterNumber(), GflashShowino::getDepth(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPathLengthOnEcal(), GflashShowino::getPosition(), GflashShowino::getShowerType(), GflashShowino::getStepLengthToHcal(), Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, longEcal, longHcal, theShowino, and twoGammaProfile().

Referenced by hadronicParameterization().

◆ medianLateralArm()

double GflashHadronShowerProfile::medianLateralArm	(	double	depth,
		Gflash::CalorimeterNumber	kCalor
	)

protected

Definition at line 252 of file GflashHadronShowerProfile.cc.

                                                                                                         {
   double lateralArm = 0.0;
   if (kCalor != Gflash::kNULL) {
     double showerDepthR50X = std::min(showerDepthR50 / 22.4, Gflash::maxShowerDepthforR50);
     double R50 = lateralPar[kCalor][0] + std::max(0.0, lateralPar[kCalor][1]) * showerDepthR50X;
     double varinanceR50 = std::pow((lateralPar[kCalor][2] + lateralPar[kCalor][3] * showerDepthR50X) * R50, 2);
  
     // Simulation of lognormal distribution
  
     if (R50 > 0) {
       double sigmaSq = std::log(varinanceR50 / (R50 * R50) + 1.0);
       double sigmaR50 = std::sqrt(sigmaSq);
       double meanR50 = std::log(R50) - (sigmaSq / 2.);
  
       lateralArm = std::exp(meanR50 + sigmaR50 * CLHEP::RandGaussQ::shoot());
     }
   }
   return lateralArm;
 }

References JetChargeProducer_cfi::exp, Gflash::kNULL, lateralPar, dqm-mbProfile::log, SiStripPI::max, Gflash::maxShowerDepthforR50, min(), funct::pow(), and mathSSE::sqrt().

Referenced by hadronicParameterization().

◆ setEnergyScale()

void GflashHadronShowerProfile::setEnergyScale	(	double	einc,
		const Gflash3Vector &	ssp
	)

protected

◆ twoGammaProfile()

double GflashHadronShowerProfile::twoGammaProfile	(	double *	par,
		double	depth,
		Gflash::CalorimeterNumber	kIndex
	)

protected

Definition at line 505 of file GflashHadronShowerProfile.cc.

                                                                                                                {
   double twoGamma = 0.0;
  
   longPar[0] = std::min(1.0, longPar[0]);
   longPar[0] = std::max(0.0, longPar[0]);
  
   if (longPar[3] > 4.0 || longPar[4] > 4.0) {
     double rfactor = 2.0 / std::max(longPar[3], longPar[4]);
     longPar[3] = rfactor * (longPar[3] + 1.0);
     longPar[4] *= rfactor;
   }
  
   twoGamma = longPar[0] * gammaProfile(exp(longPar[1]), exp(longPar[2]), depth, Gflash::radLength[kIndex]) +
              (1 - longPar[0]) * gammaProfile(exp(longPar[3]), exp(longPar[4]), depth, Gflash::intLength[kIndex]);
   return twoGamma;
 }