GflashHadronShowerProfile Class Reference

#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

(

const edm::ParameterSet &

parSet

)

Definition at line 17 of file GflashHadronShowerProfile.cc.

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

                                                                                  : theParSet(parSet) {
  theBField = parSet.getParameter<double>("bField");
  theGflashHcalOuter = parSet.getParameter<bool>("GflashHcalOuter");

  theShowino = new GflashShowino();
  theHisto = GflashHistogram::instance();
}

GflashHadronShowerProfile::~GflashHadronShowerProfile ( )

virtual

Definition at line 25 of file GflashHadronShowerProfile.cc.

References theShowino.

                                                      {
  if (theShowino)
    delete theShowino;
}

Member Function Documentation

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

protected

Definition at line 487 of file GflashHadronShowerProfile.cc.

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

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

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

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

protected

Definition at line 396 of file GflashHadronShowerProfile.cc.

References mps_fire::i, gen::k, and mathSSE::sqrt().

Referenced by getFluctuationVector(), and getGflashHitList().

                                                                                            {
  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];
    }
  }
}

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

protected

Definition at line 480 of file GflashHadronShowerProfile.cc.

References jets_cff::func, and cmsBatch::log.

Referenced by getGflashHitList().

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

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

protected

Definition at line 473 of file GflashHadronShowerProfile.cc.

References jets_cff::func, and cmsBatch::log.

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

                                                                      {
  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;
}

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

protected

Definition at line 494 of file GflashHadronShowerProfile.cc.

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

Referenced by getGflashHitList(), and twoGammaProfile().

                                                                                                               {
  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;
}

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

protected

Definition at line 358 of file GflashHadronShowerProfile.cc.

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

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

                                                                                                   {
  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;
}

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

inline

Definition at line 33 of file GflashHadronShowerProfile.h.

References alpha, pfBoostedDoubleSVAK8TagInfos_cfi::beta, egammaForCoreTracking_cff::depth, depthScale(), doCholeskyReduction(), fLnE1(), fTanh(), gammaProfile(), getFluctuationVector(), getNumberOfSpots(), hoProfile(), locateHitPosition(), longitudinalProfile(), medianLateralArm(), point, setEnergyScale(), theGflashHitList, and twoGammaProfile().

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

{ return theGflashHitList; };

GflashShowino* GflashHadronShowerProfile::getGflashShowino ( )

inline

Definition at line 32 of file GflashHadronShowerProfile.h.

References theShowino.

Referenced by CalorimetryManager::HDShowerSimulation().

{ return theShowino; }

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

protected

Definition at line 426 of file GflashHadronShowerProfile.cc.

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

Referenced by getGflashHitList(), and hadronicParameterization().

                                                                              {
  // 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;
}

void GflashHadronShowerProfile::hadronicParameterization

(

)

Definition at line 36 of file GflashHadronShowerProfile.cc.

References GflashShowino::addEnergyDeposited(), Gflash::divisionStep, MillePedeFileConverter_cfg::e, energyScale, Gflash::EtaMax, fTanh(), Gflash::getCalorimeterNumber(), GflashShowino::getEnergy(), GflashShowino::getEnergyDeposited(), GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getGlobalTime(), GflashShowino::getHelix(), getNumberOfSpots(), GflashShowino::getPathLength(), GflashTrajectory::getPathLengthAtRhoEquals(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPosition(), GflashShowino::getPositionAtShower(), GflashShowino::getStepLengthToOut(), GeV, Gflash::ho_nonzero, hoProfile(), Gflash::kHB, Gflash::kHE, Gflash::kHO, Gflash::kNULL, locateHitPosition(), cmsBatch::log, longitudinalProfile(), SiStripPI::max, medianLateralArm(), Gflash::MinEnergyCutOffForHO, Gflash::Rmax, Gflash::Rmin, GflashHit::setEnergy(), GflashShowino::setPathLength(), GflashHit::setPosition(), GflashHit::setTime(), theGflashHcalOuter, theGflashHitList, theShowino, GflashShowino::updateShowino(), and Gflash::Zmax.

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

                                                         {
  // 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;
}

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

protected

Definition at line 346 of file GflashHadronShowerProfile.cc.

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

Referenced by getGflashHitList(), and hadronicParameterization().

                                                                              {
  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;
}

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

Definition at line 30 of file GflashHadronShowerProfile.cc.

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

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

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

void GflashHadronShowerProfile::loadParameters ( )

virtual

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

Definition at line 247 of file GflashHadronShowerProfile.cc.

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

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

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

protected

Definition at line 272 of file GflashHadronShowerProfile.cc.

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

Referenced by getGflashHitList(), and hadronicParameterization().

                                                                                                          {
  // 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;
}

double GflashHadronShowerProfile::longitudinalProfile ( )

protected

Definition at line 299 of file GflashHadronShowerProfile.cc.

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 getGflashHitList(), and hadronicParameterization().

                                                      {
  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;
}

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

protected

Definition at line 252 of file GflashHadronShowerProfile.cc.

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

Referenced by getGflashHitList(), and hadronicParameterization().

                                                                                                        {
  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;
}

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

protected

Referenced by getGflashHitList().

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

protected

Definition at line 505 of file GflashHadronShowerProfile.cc.

References egammaForCoreTracking_cff::depth, JetChargeProducer_cfi::exp, gammaProfile(), Gflash::intLength, SiStripPI::max, min(), and Gflash::radLength.

Referenced by getGflashHitList(), and longitudinalProfile().

                                                                                                               {
  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;
}

Member Data Documentation

double GflashHadronShowerProfile::averageSpotEnergy[Gflash::kNumberCalorimeter]

protected

Definition at line 63 of file GflashHadronShowerProfile.h.

double GflashHadronShowerProfile::energyScale[Gflash::kNumberCalorimeter]

protected

Definition at line 62 of file GflashHadronShowerProfile.h.

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

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

protected

Definition at line 66 of file GflashHadronShowerProfile.h.

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

double GflashHadronShowerProfile::longEcal[Gflash::NPar]

protected

Definition at line 64 of file GflashHadronShowerProfile.h.

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

double GflashHadronShowerProfile::longHcal[Gflash::NPar]

protected

Definition at line 65 of file GflashHadronShowerProfile.h.

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

double GflashHadronShowerProfile::theBField

protected

Definition at line 56 of file GflashHadronShowerProfile.h.

Referenced by GflashHadronShowerProfile(), and initialize().

bool GflashHadronShowerProfile::theGflashHcalOuter

protected

Definition at line 57 of file GflashHadronShowerProfile.h.

Referenced by GflashHadronShowerProfile(), and hadronicParameterization().

std::vector<GflashHit> GflashHadronShowerProfile::theGflashHitList

protected

Definition at line 68 of file GflashHadronShowerProfile.h.

Referenced by getGflashHitList(), and hadronicParameterization().

GflashHistogram* GflashHadronShowerProfile::theHisto

protected

Definition at line 60 of file GflashHadronShowerProfile.h.

Referenced by GflashHadronShowerProfile(), and locateHitPosition().

edm::ParameterSet GflashHadronShowerProfile::theParSet

protected

Definition at line 55 of file GflashHadronShowerProfile.h.

GflashShowino* GflashHadronShowerProfile::theShowino

protected

Definition at line 59 of file GflashHadronShowerProfile.h.

Referenced by getGflashShowino(), getNumberOfSpots(), GflashHadronShowerProfile(), hadronicParameterization(), hoProfile(), initialize(), GflashProtonShowerProfile::loadParameters(), GflashKaonMinusShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), GflashAntiProtonShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), longitudinalProfile(), and ~GflashHadronShowerProfile().