GFlash/src/GflashHadronShowerProfile.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
0004 #include "SimGeneral/GFlash/interface/GflashHit.h"
0005 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
0006
0007 #include <CLHEP/GenericFunctions/IncompleteGamma.hh>
0008 #include <CLHEP/GenericFunctions/LogGamma.hh>
0009 #include <CLHEP/Random/RandGaussQ.h>
0010 #include <CLHEP/Random/RandPoissonQ.h>
0011 #include <CLHEP/Random/Randomize.h>
0012 #include <CLHEP/Units/PhysicalConstants.h>
0013 #include <CLHEP/Units/SystemOfUnits.h>
0014
0015 #include <cmath>
0016
0017 GflashHadronShowerProfile::GflashHadronShowerProfile(const edm::ParameterSet &parSet) : theParSet(parSet) {
0018   theBField = parSet.getParameter<double>("bField");
0019   theGflashHcalOuter = parSet.getParameter<bool>("GflashHcalOuter");
0020
0021   theShowino = new GflashShowino();
0022   theHisto = GflashHistogram::instance();
0023 }
0024
0025 GflashHadronShowerProfile::~GflashHadronShowerProfile() {
0026   if (theShowino)
0027     delete theShowino;
0028 }
0029
0030 void GflashHadronShowerProfile::initialize(
0031     int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum) {
0032   // initialize GflashShowino for this track
0033   theShowino->initialize(showerType, energy, globalTime, charge, position, momentum, theBField);
0034 }
0035
0036 void GflashHadronShowerProfile::hadronicParameterization() {
0037   // The skeleton of this method is based on the fortran code gfshow.F
0038   // originally written by S. Peters and G. Grindhammer (also see NIM A290
0039   // (1990) 469-488), but longitudinal parameterizations of hadron showers are
0040   // significantly modified for the CMS calorimeter
0041
0042   // unit convention: energy in [GeV] and length in [cm]
0043   // intrinsic properties of hadronic showers (lateral shower profile)
0044
0045   // The step size of showino along the helix trajectory in cm unit
0046   double showerDepthR50 = 0.0;
0047   bool firstHcalHit = true;
0048
0049   // initial valuses that will be changed as the shower developes
0050   double stepLengthLeft = theShowino->getStepLengthToOut();
0051
0052   double deltaStep = 0.0;
0053   double showerDepth = 0.0;
0054
0055   Gflash::CalorimeterNumber whichCalor = Gflash::kNULL;
0056
0057   theGflashHitList.clear();
0058
0059   GflashHit aHit;
0060
0061   while (stepLengthLeft > 0.0) {
0062     // update shower depth and stepLengthLeft
0063     if (stepLengthLeft < Gflash::divisionStep) {
0064       deltaStep = stepLengthLeft;
0065       stepLengthLeft = 0.0;
0066     } else {
0067       deltaStep = Gflash::divisionStep;
0068       stepLengthLeft -= deltaStep;
0069     }
0070
0071     showerDepth += deltaStep;
0072     showerDepthR50 += deltaStep;
0073
0074     // update GflashShowino
0075     theShowino->updateShowino(deltaStep);
0076
0077     // evaluate energy in this deltaStep along the longitudinal shower profile
0078     double heightProfile = 0.;
0079     double deltaEnergy = 0.;
0080
0081     whichCalor = Gflash::getCalorimeterNumber(theShowino->getPosition());
0082
0083     // skip if Showino is outside envelopes
0084     if (whichCalor == Gflash::kNULL)
0085       continue;
0086
0087     heightProfile = longitudinalProfile();
0088
0089     // skip if the delta energy for this step will be very small
0090     if (heightProfile < 1.00e-08)
0091       continue;
0092
0093     // get energy deposition for this step
0094     deltaEnergy = heightProfile * Gflash::divisionStep * energyScale[whichCalor];
0095     theShowino->addEnergyDeposited(deltaEnergy);
0096
0097     // apply sampling fluctuation if showino is inside the sampling calorimeter
0098     double hadronicFraction = 1.0;
0099     double fluctuatedEnergy = deltaEnergy;
0100
0101     int nSpotsInStep =
0102         std::max(1, static_cast<int>(getNumberOfSpots(whichCalor) * (deltaEnergy / energyScale[whichCalor])));
0103     double sampleSpotEnergy = hadronicFraction * fluctuatedEnergy / nSpotsInStep;
0104
0105     // Lateral shape and fluctuations
0106
0107     if ((whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) && firstHcalHit) {
0108       firstHcalHit = false;
0109       // reset the showerDepth used in the lateral parameterization inside Hcal
0110       showerDepthR50 = Gflash::divisionStep;
0111     }
0112
0113     // evaluate the fluctuated median of the lateral distribution, R50
0114     double R50 = medianLateralArm(showerDepthR50, whichCalor);
0115
0116     double hitEnergy = sampleSpotEnergy * CLHEP::GeV;
0117     double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond;
0118
0119     Gflash::CalorimeterNumber hitCalor = Gflash::kNULL;
0120
0121     for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
0122       // Compute global position of generated spots with taking into account
0123       // magnetic field Divide deltaStep into nSpotsInStep and give a spot a
0124       // global position
0125       double incrementPath =
0126           theShowino->getPathLength() + (deltaStep / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0127
0128       // trajectoryPoint give a spot an imaginary point along the shower
0129       // development
0130       GflashTrajectoryPoint trajectoryPoint;
0131       theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, incrementPath);
0132
0133       Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, R50);
0134
0135       hitCalor = Gflash::getCalorimeterNumber(hitPosition);
0136
0137       if (hitCalor == Gflash::kNULL)
0138         continue;
0139
0140       hitPosition *= CLHEP::cm;
0141
0142       aHit.setTime(hitTime);
0143       aHit.setEnergy(hitEnergy);
0144       aHit.setPosition(hitPosition);
0145       theGflashHitList.push_back(aHit);
0146
0147     }  // end of for spot iteration
0148   }    // end of while for longitudinal integration
0149
0150   // HO parameterization
0151
0152   if (theShowino->getEnergy() > Gflash::MinEnergyCutOffForHO &&
0153       fabs(theShowino->getPositionAtShower().pseudoRapidity()) < Gflash::EtaMax[Gflash::kHO] && theGflashHcalOuter) {
0154     // non zero ho fraction to simulate based on geant4
0155     double nonzeroProb = 0.7 * fTanh(theShowino->getEnergy(), Gflash::ho_nonzero);
0156     double r0 = CLHEP::HepUniformRand();
0157     double leftoverE = theShowino->getEnergy() - theShowino->getEnergyDeposited();
0158
0159     //@@@ nonzeroProb is not random - need further correlation for non-zero HO
0160     // energy
0161
0162     if (r0 < nonzeroProb && leftoverE > 0.0) {
0163       // starting path Length and stepLengthLeft
0164       double pathLength = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHO] - 10);
0165       stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO] + 10) - pathLength;
0166       showerDepth = pathLength - theShowino->getPathLengthAtShower();
0167
0168       theShowino->setPathLength(pathLength);
0169
0170       double pathLengthx = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO]);
0171       double pathLengthy = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
0172
0173       while (stepLengthLeft > 0.0) {
0174         // update shower depth and stepLengthLeft
0175         if (stepLengthLeft < Gflash::divisionStep) {
0176           deltaStep = stepLengthLeft;
0177           stepLengthLeft = 0.0;
0178         } else {
0179           deltaStep = Gflash::divisionStep;
0180           stepLengthLeft -= deltaStep;
0181         }
0182
0183         showerDepth += deltaStep;
0184
0185         // update GflashShowino
0186         theShowino->updateShowino(deltaStep);
0187
0188         // evaluate energy in this deltaStep along the longitudinal shower
0189         // profile
0190         double heightProfile = 0.;
0191         double deltaEnergy = 0.;
0192
0193         double hoScale = leftoverE * (pathLengthx - pathLengthy) / (pathLengthx - theShowino->getPathLengthAtShower());
0194         double refDepth =
0195             theShowino->getPathLength() - theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
0196
0197         if (refDepth > 0) {
0198           heightProfile = hoProfile(theShowino->getPathLength(), refDepth);
0199           deltaEnergy = heightProfile * Gflash::divisionStep * hoScale;
0200         }
0201
0202         int nSpotsInStep = std::max(
0203             50, static_cast<int>((160. + 40 * CLHEP::RandGaussQ::shoot()) * std::log(theShowino->getEnergy()) + 50.));
0204
0205         double hoFraction = 1.00;
0206         double poissonProb = CLHEP::RandPoissonQ::shoot(1.0);
0207
0208         double fluctuatedEnergy = deltaEnergy * poissonProb;
0209         double sampleSpotEnergy = hoFraction * fluctuatedEnergy / nSpotsInStep;
0210
0211         // Lateral shape and fluctuations
0212
0213         // evaluate the fluctuated median of the lateral distribution, R50
0214         double R50 = medianLateralArm(showerDepth, Gflash::kHB);
0215
0216         double hitEnergy = sampleSpotEnergy * CLHEP::GeV;
0217         double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond;
0218
0219         for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
0220           double incrementPath =
0221               theShowino->getPathLength() + (deltaStep / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0222
0223           // trajectoryPoint give a spot an imaginary point along the shower
0224           // development
0225           GflashTrajectoryPoint trajectoryPoint;
0226           theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, incrementPath);
0227
0228           Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, R50);
0229           hitPosition *= CLHEP::cm;
0230
0231           if (std::fabs(hitPosition.getZ() / CLHEP::cm) > Gflash::Zmax[Gflash::kHO])
0232             continue;
0233
0234           aHit.setTime(hitTime);
0235           aHit.setEnergy(hitEnergy);
0236           aHit.setPosition(hitPosition);
0237           theGflashHitList.push_back(aHit);
0238
0239         }  // end of for HO spot iteration
0240       }    // end of while for HO longitudinal integration
0241     }
0242   }
0243
0244   //  delete theGflashNavigator;
0245 }
0246
0247 void GflashHadronShowerProfile::loadParameters() {
0248   edm::LogInfo("SimGeneralGFlash") << "GflashHadronShowerProfile::loadParameters() "
0249                                    << "should be implimented for each particle type";
0250 }
0251
0252 double GflashHadronShowerProfile::medianLateralArm(double showerDepthR50, Gflash::CalorimeterNumber kCalor) {
0253   double lateralArm = 0.0;
0254   if (kCalor != Gflash::kNULL) {
0255     double showerDepthR50X = std::min(showerDepthR50 / 22.4, Gflash::maxShowerDepthforR50);
0256     double R50 = lateralPar[kCalor][0] + std::max(0.0, lateralPar[kCalor][1]) * showerDepthR50X;
0257     double varinanceR50 = std::pow((lateralPar[kCalor][2] + lateralPar[kCalor][3] * showerDepthR50X) * R50, 2);
0258
0259     // Simulation of lognormal distribution
0260
0261     if (R50 > 0) {
0262       double sigmaSq = std::log(varinanceR50 / (R50 * R50) + 1.0);
0263       double sigmaR50 = std::sqrt(sigmaSq);
0264       double meanR50 = std::log(R50) - (sigmaSq / 2.);
0265
0266       lateralArm = std::exp(meanR50 + sigmaR50 * CLHEP::RandGaussQ::shoot());
0267     }
0268   }
0269   return lateralArm;
0270 }
0271
0272 Gflash3Vector GflashHadronShowerProfile::locateHitPosition(GflashTrajectoryPoint &point, double lateralArm) {
0273   // Smearing in r according to f(r)= 2.*r*R50**2/(r**2+R50**2)**2
0274   double rnunif = CLHEP::HepUniformRand();
0275   double rxPDF = std::sqrt(rnunif / (1. - rnunif));
0276   double rShower = lateralArm * rxPDF;
0277
0278   // rShower within maxLateralArmforR50
0279   rShower = std::min(Gflash::maxLateralArmforR50, rShower);
0280
0281   // Uniform smearing in phi
0282   double azimuthalAngle = CLHEP::twopi * CLHEP::HepUniformRand();
0283
0284   // actual spot position by adding a radial vector to a trajectoryPoint
0285   Gflash3Vector position = point.getPosition() + rShower * std::cos(azimuthalAngle) * point.getOrthogonalUnitVector() +
0286                            rShower * std::sin(azimuthalAngle) * point.getCrossUnitVector();
0287
0288   //@@@debugging histograms
0289   if (theHisto->getStoreFlag()) {
0290     theHisto->rshower->Fill(rShower);
0291     theHisto->lateralx->Fill(rShower * std::cos(azimuthalAngle));
0292     theHisto->lateraly->Fill(rShower * std::sin(azimuthalAngle));
0293   }
0294   return position;
0295 }
0296
0297 // double GflashHadronShowerProfile::longitudinalProfile(double showerDepth,
0298 // double pathLength) {
0299 double GflashHadronShowerProfile::longitudinalProfile() {
0300   double heightProfile = 0.0;
0301
0302   Gflash3Vector pos = theShowino->getPosition();
0303   int showerType = theShowino->getShowerType();
0304   double showerDepth = theShowino->getDepth();
0305   double transDepth = theShowino->getStepLengthToHcal();
0306
0307   // Energy in a delta step (dz) = (energy to
0308   // deposite)*[Gamma(z+dz)-Gamma(z)]*dz where the incomplete Gamma function
0309   // gives an intergrate probability of the longitudinal shower up to the shower
0310   // depth (z). Instead, we use approximated energy; energy in dz = (energy to
0311   // deposite)*gamma(z)*dz where gamma is the Gamma-distributed probability
0312   // function
0313
0314   Gflash::CalorimeterNumber whichCalor = Gflash::getCalorimeterNumber(pos);
0315
0316   if (showerType == 0 || showerType == 4) {
0317     double shiftDepth = theShowino->getPathLengthOnEcal() - theShowino->getPathLengthAtShower();
0318     if (shiftDepth > 0) {
0319       heightProfile = twoGammaProfile(longEcal, showerDepth - shiftDepth, whichCalor);
0320     } else {
0321       heightProfile = 0.;
0322       //      std::cout << "negative shiftDepth for showerType 0 " << shiftDepth
0323       //      << std::endl;
0324     }
0325   } else if (showerType == 1 || showerType == 5) {
0326     if (whichCalor == Gflash::kESPM || whichCalor == Gflash::kENCA) {
0327       heightProfile = twoGammaProfile(longEcal, showerDepth, whichCalor);
0328     } else if (whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) {
0329       heightProfile = twoGammaProfile(longHcal, showerDepth - transDepth, whichCalor);
0330     } else
0331       heightProfile = 0.;
0332   } else if (showerType == 2 || showerType == 6) {
0333     // two gammas between crystal and Hcal
0334     if ((showerDepth - transDepth) > 0.0) {
0335       heightProfile = twoGammaProfile(longHcal, showerDepth - transDepth, Gflash::kHB);
0336     } else
0337       heightProfile = 0.;
0338   } else if (showerType == 3 || showerType == 7) {
0339     // two gammas inside Hcal
0340     heightProfile = twoGammaProfile(longHcal, showerDepth, Gflash::kHB);
0341   }
0342
0343   return heightProfile;
0344 }
0345
0346 double GflashHadronShowerProfile::hoProfile(double pathLength, double refDepth) {
0347   double heightProfile = 0;
0348
0349   GflashTrajectoryPoint tempPoint;
0350   theShowino->getHelix()->getGflashTrajectoryPoint(tempPoint, pathLength);
0351
0352   double dint = 1.4 * Gflash::intLength[Gflash::kHO] * std::sin(tempPoint.getPosition().getTheta());
0353   heightProfile = std::exp(-1.0 * refDepth / dint);
0354
0355   return heightProfile;
0356 }
0357
0358 void GflashHadronShowerProfile::getFluctuationVector(double *lowTriangle, double *correlationVector) {
0359   const int dim = Gflash::NPar;
0360
0361   double **xr = new double *[dim];
0362   double **xrho = new double *[dim];
0363
0364   for (int j = 0; j < dim; j++) {
0365     xr[j] = new double[dim];
0366     xrho[j] = new double[dim];
0367   }
0368
0369   for (int i = 0; i < dim; i++) {
0370     for (int j = 0; j < i + 1; j++) {
0371       if (j == i)
0372         xrho[i][j] = 1.0;
0373       else {
0374         xrho[i][j] = lowTriangle[i * (i - 1) / 2 + j];
0375         xrho[j][i] = xrho[i][j];
0376       }
0377     }
0378   }
0379
0380   doCholeskyReduction(xrho, xr, dim);
0381
0382   for (int i = 0; i < dim; i++) {
0383     for (int j = 0; j < i + 1; j++) {
0384       correlationVector[i * (i + 1) / 2 + j] = xr[i][j];
0385     }
0386   }
0387
0388   for (int j = 0; j < dim; j++)
0389     delete[] xr[j];
0390   delete[] xr;
0391   for (int j = 0; j < dim; j++)
0392     delete[] xrho[j];
0393   delete[] xrho;
0394 }
0395
0396 void GflashHadronShowerProfile::doCholeskyReduction(double **vv, double **cc, const int ndim) {
0397   double sumCjkSquare;
0398   double vjjLess;
0399   double sumCikjk;
0400
0401   cc[0][0] = std::sqrt(vv[0][0]);
0402
0403   for (int j = 1; j < ndim; j++) {
0404     cc[j][0] = vv[j][0] / cc[0][0];
0405   }
0406
0407   for (int j = 1; j < ndim; j++) {
0408     sumCjkSquare = 0.0;
0409     for (int k = 0; k < j; k++)
0410       sumCjkSquare += cc[j][k] * cc[j][k];
0411
0412     vjjLess = vv[j][j] - sumCjkSquare;
0413
0414     // check for the case that vjjLess is negative
0415     cc[j][j] = std::sqrt(std::fabs(vjjLess));
0416
0417     for (int i = j + 1; i < ndim; i++) {
0418       sumCikjk = 0.;
0419       for (int k = 0; k < j; k++)
0420         sumCikjk += cc[i][k] * cc[j][k];
0421       cc[i][j] = (vv[i][j] - sumCikjk) / cc[j][j];
0422     }
0423   }
0424 }
0425
0426 int GflashHadronShowerProfile::getNumberOfSpots(Gflash::CalorimeterNumber kCalor) {
0427   // generator number of spots: energy dependent Gamma distribution of Nspots
0428   // based on Geant4 replacing old parameterization of H1, int numberOfSpots =
0429   // std::max( 50, static_cast<int>(80.*std::log(einc)+50.));
0430
0431   double einc = theShowino->getEnergy();
0432   int showerType = theShowino->getShowerType();
0433
0434   int numberOfSpots = 0;
0435   double nmean = 0.0;
0436   double nsigma = 0.0;
0437
0438   if (showerType == 0 || showerType == 1 || showerType == 4 || showerType == 5) {
0439     if (kCalor == Gflash::kESPM || kCalor == Gflash::kENCA) {
0440       nmean = 10000 + 5000 * log(einc);
0441       nsigma = 1000;
0442     }
0443     if (kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
0444       nmean = 5000 + 2500 * log(einc);
0445       nsigma = 500;
0446     }
0447   } else if (showerType == 2 || showerType == 3 || showerType == 6 || showerType == 7) {
0448     if (kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
0449       nmean = 5000 + 2500 * log(einc);
0450       nsigma = 500;
0451     } else {
0452       nmean = 10000;
0453       nsigma = 1000;
0454     }
0455   }
0456   //@@@need correlation and individual fluctuation on alphaNspots and betaNspots
0457   // here: evaluating covariance should be straight forward since the
0458   // distribution is 'one' Gamma
0459
0460   numberOfSpots = std::max(500, static_cast<int>(nmean + nsigma * CLHEP::RandGaussQ::shoot()));
0461
0462   // until we optimize the reduction scale in the number of Nspots
0463
0464   if (kCalor == Gflash::kESPM || kCalor == Gflash::kENCA) {
0465     numberOfSpots = static_cast<int>(numberOfSpots / 100);
0466   } else {
0467     numberOfSpots = static_cast<int>(numberOfSpots / 3.0);
0468   }
0469
0470   return numberOfSpots;
0471 }
0472
0473 double GflashHadronShowerProfile::fTanh(double einc, const double *par) {
0474   double func = 0.0;
0475   if (einc > 0.0)
0476     func = par[0] + par[1] * std::tanh(par[2] * (std::log(einc) - par[3])) + par[4] * std::log(einc);
0477   return func;
0478 }
0479
0480 double GflashHadronShowerProfile::fLnE1(double einc, const double *par) {
0481   double func = 0.0;
0482   if (einc > 0.0)
0483     func = par[0] + par[1] * std::log(einc);
0484   return func;
0485 }
0486
0487 double GflashHadronShowerProfile::depthScale(double ssp, double ssp0, double length) {
0488   double func = 0.0;
0489   if (length > 0.0)
0490     func = std::pow((ssp - ssp0) / length, 2.0);
0491   return func;
0492 }
0493
0494 double GflashHadronShowerProfile::gammaProfile(double alpha, double beta, double showerDepth, double lengthUnit) {
0495   double gamma = 0.0;
0496   //  if(alpha > 0 && beta > 0 && lengthUnit > 0) {
0497   if (showerDepth > 0.0) {
0498     Genfun::LogGamma lgam;
0499     double x = showerDepth * (beta / lengthUnit);
0500     gamma = (beta / lengthUnit) * std::pow(x, alpha - 1.0) * std::exp(-x) / std::exp(lgam(alpha));
0501   }
0502   return gamma;
0503 }
0504
0505 double GflashHadronShowerProfile::twoGammaProfile(double *longPar, double depth, Gflash::CalorimeterNumber kIndex) {
0506   double twoGamma = 0.0;
0507
0508   longPar[0] = std::min(1.0, longPar[0]);
0509   longPar[0] = std::max(0.0, longPar[0]);
0510
0511   if (longPar[3] > 4.0 || longPar[4] > 4.0) {
0512     double rfactor = 2.0 / std::max(longPar[3], longPar[4]);
0513     longPar[3] = rfactor * (longPar[3] + 1.0);
0514     longPar[4] *= rfactor;
0515   }
0516
0517   twoGamma = longPar[0] * gammaProfile(exp(longPar[1]), exp(longPar[2]), depth, Gflash::radLength[kIndex]) +
0518              (1 - longPar[0]) * gammaProfile(exp(longPar[3]), exp(longPar[4]), depth, Gflash::intLength[kIndex]);
0519   return twoGamma;
0520 }