GFlash/src/GflashEMShowerProfile.cc

0001 //
0002 // initial setup : Soon Jun & Dongwook Jang
0003 // Translated from Fortran code.
0004 //
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "SimGeneral/GFlash/interface/Gflash3Vector.h"
0008 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
0009 #include "SimGeneral/GFlash/interface/GflashHistogram.h"
0010 #include "SimGeneral/GFlash/interface/GflashHit.h"
0011 #include "SimGeneral/GFlash/interface/GflashShowino.h"
0012 #include "SimGeneral/GFlash/interface/GflashTrajectory.h"
0013 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
0014
0015 #include <CLHEP/GenericFunctions/IncompleteGamma.hh>
0016 #include <CLHEP/Random/RandGaussQ.h>
0017 #include <CLHEP/Random/Randomize.h>
0018 #include <CLHEP/Units/PhysicalConstants.h>
0019
0020 GflashEMShowerProfile::GflashEMShowerProfile(const edm::ParameterSet &parSet) : theParSet(parSet) {
0021   theBField = parSet.getParameter<double>("bField");
0022   theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
0023   theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");
0024
0025   jCalorimeter = Gflash::kNULL;
0026
0027   theShowino = new GflashShowino();
0028
0029   theHisto = GflashHistogram::instance();
0030 }
0031
0032 GflashEMShowerProfile::~GflashEMShowerProfile() {
0033   if (theShowino)
0034     delete theShowino;
0035 }
0036
0037 void GflashEMShowerProfile::initialize(
0038     int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum) {
0039   theShowino->initialize(showerType, energy, globalTime, charge, position, momentum, theBField);
0040 }
0041
0042 void GflashEMShowerProfile::parameterization() {
0043   // This part of code is copied from the original GFlash Fortran code.
0044   // reference : hep-ex/0001020v1
0045   // The units used in Geant4 internally are in mm, MeV.
0046   // For simplicity, the units here are in cm, GeV.
0047
0048   const double energyCutoff = 0.01;
0049   const int maxNumberOfSpots = 100000;
0050
0051   double incomingEnergy = theShowino->getEnergy();
0052   Gflash3Vector showerStartingPosition = theShowino->getPositionAtShower();
0053
0054   // find the calorimeter at the shower starting point
0055   jCalorimeter = Gflash::getCalorimeterNumber(showerStartingPosition);
0056
0057   double logEinc = std::log(incomingEnergy);
0058   double y = incomingEnergy / Gflash::criticalEnergy;  // y = E/Ec, criticalEnergy is in GeV
0059   double logY = std::log(y);
0060
0061   // Total number of spots are not yet optimized.
0062   double nSpots = 93.0 * std::log(Gflash::Z[jCalorimeter]) * std::pow(incomingEnergy, 0.876);
0063
0064   // path Length from the origin to the shower starting point in cm
0065   double pathLength0 = theShowino->getPathLengthAtShower();
0066   double pathLength = pathLength0;  // this will grow along the shower development
0067
0068   //--- 2.2  Fix intrinsic properties of em. showers.
0069
0070   double fluctuatedTmax = std::log(logY - 0.7157);
0071   double fluctuatedAlpha = std::log(0.7996 + (0.4581 + 1.8628 / Gflash::Z[jCalorimeter]) * logY);
0072
0073   double sigmaTmax = 1.0 / (-1.4 + 1.26 * logY);
0074   double sigmaAlpha = 1.0 / (-0.58 + 0.86 * logY);
0075   double rho = 0.705 - 0.023 * logY;
0076   double sqrtPL = std::sqrt((1.0 + rho) / 2.0);
0077   double sqrtLE = std::sqrt((1.0 - rho) / 2.0);
0078
0079   double norm1 = CLHEP::RandGaussQ::shoot();
0080   double norm2 = CLHEP::RandGaussQ::shoot();
0081   double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
0082   double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
0083
0084   // tmax, alpha, beta : parameters of gamma distribution
0085   double tmax = std::exp(tempTmax);
0086   double alpha = std::exp(tempAlpha);
0087   double beta = std::max(0.0, (alpha - 1.0) / tmax);
0088
0089   // spot fluctuations are added to tmax, alpha, beta
0090   double averageTmax = logY - 0.858;
0091   double averageAlpha = 0.21 + (0.492 + 2.38 / Gflash::Z[jCalorimeter]) * logY;
0092   double spotTmax = averageTmax * (0.698 + .00212 * Gflash::Z[jCalorimeter]);
0093   double spotAlpha = averageAlpha * (0.639 + .00334 * Gflash::Z[jCalorimeter]);
0094   double spotBeta = std::max(0.0, (spotAlpha - 1.0) / spotTmax);
0095
0096   if (theHisto->getStoreFlag()) {
0097     theHisto->em_incE->Fill(incomingEnergy);
0098     theHisto->em_ssp_rho->Fill(showerStartingPosition.rho());
0099     theHisto->em_ssp_z->Fill(showerStartingPosition.z());
0100   }
0101
0102   //  parameters for lateral distribution and fluctuation
0103   double z1 = 0.0251 + 0.00319 * logEinc;
0104   double z2 = 0.1162 - 0.000381 * Gflash::Z[jCalorimeter];
0105
0106   double k1 = 0.659 - 0.00309 * Gflash::Z[jCalorimeter];
0107   double k2 = 0.645;
0108   double k3 = -2.59;
0109   double k4 = 0.3585 + 0.0421 * logEinc;
0110
0111   double p1 = 2.623 - 0.00094 * Gflash::Z[jCalorimeter];
0112   double p2 = 0.401 + 0.00187 * Gflash::Z[jCalorimeter];
0113   double p3 = 1.313 - 0.0686 * logEinc;
0114
0115   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : e25Scale = 1.006
0116   // for EB with ecalNotContainment = 1.0. Now e25Scale is a configurable
0117   // parameter with default ecalNotContainment which is 0.97 for EB and 0.975
0118   // for EE. So if ecalNotContainment constants are to be changed in the future,
0119   // e25Scale should be changed accordingly.
0120   double e25Scale = 1.0;
0121   if (jCalorimeter == Gflash::kESPM)
0122     e25Scale = theEnergyScaleEB;
0123   else if (jCalorimeter == Gflash::kENCA)
0124     e25Scale = theEnergyScaleEE;
0125
0126   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : p1 *= 0.965
0127   p1 *= 0.965;
0128
0129   // preparation of longitudinal integration
0130   double stepLengthLeft = getDistanceToOut(jCalorimeter);
0131
0132   int nSpots_sd = 0;                // count total number of spots in SD
0133   double zInX0 = 0.0;               // shower depth in X0 unit
0134   double deltaZInX0 = 0.0;          // segment of depth in X0 unit
0135   double deltaZ = 0.0;              // segment of depth in cm
0136   double stepLengthLeftInX0 = 0.0;  // step length left in X0 unit
0137
0138   const double divisionStepInX0 = 0.1;  // step size in X0 unit
0139   double energy = incomingEnergy;       // energy left in GeV
0140
0141   Genfun::IncompleteGamma gammaDist;
0142
0143   double energyInGamma = 0.0;     // energy in a specific depth(z) according to Gamma distribution
0144   double preEnergyInGamma = 0.0;  // energy calculated in a previous depth
0145   double sigmaInGamma = 0.;       // sigma of energy in a specific depth(z) according
0146                                   // to Gamma distribution
0147   double preSigmaInGamma = 0.0;   // sigma of energy in a previous depth
0148
0149   // energy segment in Gamma distribution of shower in each step
0150   double deltaEnergy = 0.0;  // energy in deltaZ
0151   int spotCounter = 0;       // keep track of number of spots generated
0152
0153   // step increment along the shower direction
0154   double deltaStep = 0.0;
0155
0156   theGflashHitList.clear();
0157
0158   // loop for longitudinal integration
0159   while (energy > 0.0 && stepLengthLeft > 0.0) {
0160     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
0161
0162     if (stepLengthLeftInX0 < divisionStepInX0) {
0163       deltaZInX0 = stepLengthLeftInX0;
0164       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0165       stepLengthLeft = 0.0;
0166     } else {
0167       deltaZInX0 = divisionStepInX0;
0168       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0169       stepLengthLeft -= deltaZ;
0170     }
0171
0172     zInX0 += deltaZInX0;
0173
0174     int nSpotsInStep = 0;
0175
0176     if (energy > energyCutoff) {
0177       preEnergyInGamma = energyInGamma;
0178       gammaDist.a().setValue(alpha);            // alpha
0179       energyInGamma = gammaDist(beta * zInX0);  // beta
0180       double energyInDeltaZ = energyInGamma - preEnergyInGamma;
0181       deltaEnergy = std::min(energy, incomingEnergy * energyInDeltaZ);
0182
0183       preSigmaInGamma = sigmaInGamma;
0184       gammaDist.a().setValue(spotAlpha);           // alpha spot
0185       sigmaInGamma = gammaDist(spotBeta * zInX0);  // beta spot
0186       nSpotsInStep = std::max(1, int(nSpots * (sigmaInGamma - preSigmaInGamma)));
0187     } else {
0188       deltaEnergy = energy;
0189       preSigmaInGamma = sigmaInGamma;
0190       nSpotsInStep = std::max(1, int(nSpots * (1.0 - preSigmaInGamma)));
0191     }
0192
0193     if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
0194       deltaEnergy = energy;
0195
0196     energy -= deltaEnergy;
0197
0198     if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
0199       nSpotsInStep = maxNumberOfSpots - spotCounter;
0200       if (nSpotsInStep < 1) {  // @@ check
0201         edm::LogInfo("SimGeneralGFlash") << "GflashEMShowerProfile: Too Many Spots ";
0202         edm::LogInfo("SimGeneralGFlash") << " - break to regenerate nSpotsInStep ";
0203         break;
0204       }
0205     }
0206
0207     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
0208     deltaStep += 0.5 * deltaZ;
0209     pathLength += deltaStep;
0210     deltaStep = 0.5 * deltaZ;
0211
0212     // lateral shape and fluctuations for  homogenous calo.
0213     double tScale = tmax * alpha / (alpha - 1.0) * (std::exp(fluctuatedAlpha) - 1.0) / std::exp(fluctuatedAlpha);
0214     double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
0215     double rCore = z1 + z2 * tau;
0216     double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));  // @@ check RT3 sign
0217     double p23 = (p2 - tau) / p3;
0218     double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));
0219
0220     // Deposition of spots according to lateral distr.
0221     // Apply absolute energy scale
0222     // Convert into MeV unit
0223     double hitEnergy = deltaEnergy / nSpotsInStep * e25Scale * CLHEP::GeV;
0224     double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond + (pathLength - pathLength0) / 30.0;
0225
0226     GflashHit aHit;
0227
0228     for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
0229       spotCounter++;
0230
0231       // Compute global position of generated spots with taking into account
0232       // magnetic field Divide deltaZ into nSpotsInStep and give a spot a global
0233       // position
0234       double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0235
0236       // trajectoryPoint give a spot an imaginary point along the shower
0237       // development
0238       GflashTrajectoryPoint trajectoryPoint;
0239       theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
0240
0241       double rShower = 0.0;
0242       Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, rCore, rTail, probabilityWeight, rShower);
0243
0244       // Convert into mm unit
0245       hitPosition *= CLHEP::cm;
0246
0247       if (std::fabs(hitPosition.getZ() / CLHEP::cm) > Gflash::Zmax[Gflash::kHE])
0248         continue;
0249
0250       // put energy and position to a Hit
0251       aHit.setTime(hitTime);
0252       aHit.setEnergy(hitEnergy);
0253       aHit.setPosition(hitPosition);
0254       theGflashHitList.push_back(aHit);
0255
0256       double zInX0_spot = std::abs(pathLength + incrementPath - pathLength0) / Gflash::radLength[jCalorimeter];
0257
0258       nSpots_sd++;
0259
0260       // for histogramming
0261       if (theHisto->getStoreFlag()) {
0262         theHisto->em_long->Fill(zInX0_spot, hitEnergy / CLHEP::GeV);
0263         theHisto->em_lateral->Fill(zInX0_spot, rShower / Gflash::rMoliere[jCalorimeter], hitEnergy / CLHEP::GeV);
0264         theHisto->em_long_sd->Fill(zInX0_spot, hitEnergy / CLHEP::GeV);
0265         theHisto->em_lateral_sd->Fill(zInX0_spot, rShower / Gflash::rMoliere[jCalorimeter], hitEnergy / CLHEP::GeV);
0266       }
0267
0268     }  // end of for spot iteration
0269
0270   }  // end of while for longitudinal integration
0271
0272   if (theHisto->getStoreFlag()) {
0273     theHisto->em_nSpots_sd->Fill(nSpots_sd);
0274   }
0275
0276   //  delete theGflashNavigator;
0277 }
0278
0279 double GflashEMShowerProfile::getDistanceToOut(Gflash::CalorimeterNumber kCalor) {
0280   double stepLengthLeft = 0.0;
0281   if (kCalor == Gflash::kESPM) {
0282     stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kESPM]) -
0283                      theShowino->getPathLengthAtShower();
0284   } else if (kCalor == Gflash::kENCA) {
0285     double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
0286     stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign * Gflash::Zmax[Gflash::kENCA]) -
0287                      theShowino->getPathLengthAtShower();
0288   }
0289   return stepLengthLeft;
0290 }
0291
0292 Gflash3Vector GflashEMShowerProfile::locateHitPosition(
0293     GflashTrajectoryPoint &point, double rCore, double rTail, double probability, double &rShower) {
0294   double u1 = CLHEP::HepUniformRand();
0295   double u2 = CLHEP::HepUniformRand();
0296   double rInRM = 0.0;
0297
0298   if (u1 < probability) {
0299     rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
0300   } else {
0301     rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
0302   }
0303
0304   rShower = rInRM * Gflash::rMoliere[jCalorimeter];
0305
0306   // Uniform & random rotation of spot along the azimuthal angle
0307   double azimuthalAngle = CLHEP::twopi * CLHEP::HepUniformRand();
0308
0309   // actual spot position by adding a radial vector to a trajectoryPoint
0310   Gflash3Vector position = point.getPosition() + rShower * std::cos(azimuthalAngle) * point.getOrthogonalUnitVector() +
0311                            rShower * std::sin(azimuthalAngle) * point.getCrossUnitVector();
0312
0313   return position;
0314 }