Reference Manual

Go to the documentation of this file.
00001 //
00002 // $Id: GflashEMShowerProfile.cc,v 1.5 2010/04/30 19:10:11 dwjang Exp $
00003 // initial setup : Soon Jun & Dongwook Jang
00004 // Translated from Fortran code.
00005 //
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 
00008 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
00009 #include "SimGeneral/GFlash/interface/GflashTrajectory.h"
00010 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
00011 #include "SimGeneral/GFlash/interface/GflashHit.h"
00012 #include "SimGeneral/GFlash/interface/Gflash3Vector.h"
00013 #include "SimGeneral/GFlash/interface/GflashHistogram.h"
00014 #include "SimGeneral/GFlash/interface/GflashShowino.h"
00015 
00016 #include <CLHEP/GenericFunctions/IncompleteGamma.hh>
00017 #include <CLHEP/Units/PhysicalConstants.h>
00018 #include <CLHEP/Random/Randomize.h>
00019 #include <CLHEP/Random/RandGaussQ.h>
00020 
00021 GflashEMShowerProfile::GflashEMShowerProfile(edm::ParameterSet parSet) : theParSet(parSet)
00022 {
00023   theBField = parSet.getParameter<double>("bField");
00024   theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
00025   theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");
00026 
00027   jCalorimeter = Gflash::kNULL;
00028 
00029   theShowino = new GflashShowino();
00030 
00031   theHisto = GflashHistogram::instance();
00032 
00033 }
00034 
00035 
00036 GflashEMShowerProfile::~GflashEMShowerProfile()
00037 {
00038   if(theShowino) delete theShowino;
00039 }
00040 
00041 void GflashEMShowerProfile::initialize(int showerType, double energy, double globalTime, double charge,
00042                                        Gflash3Vector &position,Gflash3Vector &momentum) 
00043 {
00044   theShowino->initialize(showerType, energy, globalTime, charge,
00045                          position, momentum, theBField);
00046 }
00047 
00048 void GflashEMShowerProfile::parameterization()
00049 {
00050   // This part of code is copied from the original GFlash Fortran code.
00051   // reference : hep-ex/0001020v1
00052   // The units used in Geant4 internally are in mm, MeV.
00053   // For simplicity, the units here are in cm, GeV.
00054 
00055   const double energyCutoff     = 0.01; 
00056   const int    maxNumberOfSpots = 100000;
00057 
00058   double incomingEnergy   = theShowino->getEnergy();
00059   Gflash3Vector showerStartingPosition = theShowino->getPositionAtShower();
00060 
00061   //find the calorimeter at the shower starting point
00062   jCalorimeter = Gflash::getCalorimeterNumber(showerStartingPosition);
00063 
00064   double logEinc = std::log(incomingEnergy);
00065   double y = incomingEnergy / Gflash::criticalEnergy; // y = E/Ec, criticalEnergy is in GeV
00066   double logY = std::log(y);
00067 
00068   // Total number of spots are not yet optimized.
00069   double nSpots = 93.0 * std::log(Gflash::Z[jCalorimeter]) * std::pow(incomingEnergy,0.876);
00070 
00071   //path Length from the origin to the shower starting point in cm
00072   double pathLength0 = theShowino->getPathLengthAtShower();
00073   double pathLength = pathLength0; // this will grow along the shower development
00074 
00075   //--- 2.2  Fix intrinsic properties of em. showers.
00076 
00077   double fluctuatedTmax = std::log(logY - 0.7157);
00078   double fluctuatedAlpha = std::log(0.7996 +(0.4581 + 1.8628/Gflash::Z[jCalorimeter])*logY);
00079 
00080   double sigmaTmax = 1.0/( -1.4  + 1.26 * logY);
00081   double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
00082   double rho = 0.705  - 0.023 * logY;
00083   double sqrtPL = std::sqrt((1.0+rho)/2.0);
00084   double sqrtLE = std::sqrt((1.0-rho)/2.0);
00085 
00086   double norm1 =  CLHEP::RandGaussQ::shoot();
00087   double norm2 =  CLHEP::RandGaussQ::shoot();
00088   double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
00089   double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
00090 
00091   // tmax, alpha, beta : parameters of gamma distribution
00092   double tmax = std::exp(tempTmax);
00093   double alpha = std::exp(tempAlpha);
00094   double beta = std::max(0.0,(alpha - 1.0)/tmax);
00095  
00096   // spot fluctuations are added to tmax, alpha, beta
00097   double averageTmax = logY-0.858;
00098   double averageAlpha = 0.21+(0.492+2.38/Gflash::Z[jCalorimeter])*logY;
00099   double spotTmax  = averageTmax * (0.698 + .00212*Gflash::Z[jCalorimeter]);
00100   double spotAlpha = averageAlpha * (0.639 + .00334*Gflash::Z[jCalorimeter]);
00101   double spotBeta = std::max(0.0,(spotAlpha-1.0)/spotTmax);
00102 
00103    if(theHisto->getStoreFlag()) {
00104     theHisto->em_incE->Fill(incomingEnergy);
00105     theHisto->em_ssp_rho->Fill(showerStartingPosition.rho());
00106     theHisto->em_ssp_z->Fill(showerStartingPosition.z());
00107   }
00108 
00109   //  parameters for lateral distribution and fluctuation
00110   double z1=0.0251+0.00319*logEinc;
00111   double z2=0.1162-0.000381*Gflash::Z[jCalorimeter];
00112 
00113   double k1=0.659 - 0.00309 * Gflash::Z[jCalorimeter];
00114   double k2=0.645;
00115   double k3=-2.59;
00116   double k4=0.3585+ 0.0421*logEinc;
00117 
00118   double p1=2.623 -0.00094*Gflash::Z[jCalorimeter];
00119   double p2=0.401 +0.00187*Gflash::Z[jCalorimeter];
00120   double p3=1.313 -0.0686*logEinc;
00121 
00122   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : e25Scale = 1.006 for EB with ecalNotContainment = 1.0.
00123   // Now e25Scale is a configurable parameter with default ecalNotContainment which is 0.97 for EB and 0.975 for EE.
00124   // So if ecalNotContainment constants are to be changed in the future, e25Scale should be changed accordingly.
00125   double e25Scale = 1.0;
00126   if(jCalorimeter == Gflash::kESPM) e25Scale = theEnergyScaleEB;
00127   else if(jCalorimeter == Gflash::kENCA) e25Scale = theEnergyScaleEE;
00128 
00129   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : p1 *= 0.965
00130   p1 *= 0.965;
00131 
00132   // preparation of longitudinal integration
00133   double stepLengthLeft = getDistanceToOut(jCalorimeter);
00134 
00135   int    nSpots_sd = 0; // count total number of spots in SD
00136   double zInX0 = 0.0; // shower depth in X0 unit
00137   double deltaZInX0 = 0.0; // segment of depth in X0 unit
00138   double deltaZ = 0.0; // segment of depth in cm
00139   double stepLengthLeftInX0 = 0.0; // step length left in X0 unit
00140 
00141   const double divisionStepInX0 = 0.1; //step size in X0 unit
00142   double energy = incomingEnergy; // energy left in GeV
00143 
00144   Genfun::IncompleteGamma gammaDist;
00145 
00146   double energyInGamma = 0.0; // energy in a specific depth(z) according to Gamma distribution
00147   double preEnergyInGamma = 0.0; // energy calculated in a previous depth
00148   double sigmaInGamma  = 0.; // sigma of energy in a specific depth(z) according to Gamma distribution
00149   double preSigmaInGamma = 0.0; // sigma of energy in a previous depth
00150 
00151   //energy segment in Gamma distribution of shower in each step  
00152   double deltaEnergy =0.0 ; // energy in deltaZ
00153   int spotCounter = 0; // keep track of number of spots generated
00154 
00155   //step increment along the shower direction
00156   double deltaStep = 0.0;
00157 
00158   theGflashHitList.clear();
00159 
00160   // loop for longitudinal integration
00161   while(energy > 0.0 && stepLengthLeft > 0.0) { 
00162 
00163     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
00164 
00165     if ( stepLengthLeftInX0 < divisionStepInX0 ) {
00166       deltaZInX0 = stepLengthLeftInX0;
00167       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
00168       stepLengthLeft = 0.0;
00169     }
00170     else {
00171       deltaZInX0 = divisionStepInX0;
00172       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
00173       stepLengthLeft -= deltaZ;
00174     }
00175 
00176     zInX0 += deltaZInX0;
00177 
00178     int nSpotsInStep = 0;
00179 
00180     if ( energy > energyCutoff  ) {
00181       preEnergyInGamma  = energyInGamma;
00182       gammaDist.a().setValue(alpha);  //alpha
00183       energyInGamma = gammaDist(beta*zInX0); //beta
00184       double energyInDeltaZ  = energyInGamma - preEnergyInGamma;
00185       deltaEnergy   = std::min(energy,incomingEnergy*energyInDeltaZ);
00186  
00187       preSigmaInGamma  = sigmaInGamma;
00188       gammaDist.a().setValue(spotAlpha);  //alpha spot
00189       sigmaInGamma = gammaDist(spotBeta*zInX0); //beta spot
00190       nSpotsInStep = std::max(1,int(nSpots * (sigmaInGamma - preSigmaInGamma)));
00191     }
00192     else {
00193       deltaEnergy = energy;
00194       preSigmaInGamma  = sigmaInGamma;
00195       nSpotsInStep = std::max(1,int(nSpots * (1.0 - preSigmaInGamma)));
00196     }
00197 
00198     if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy = energy;
00199 
00200     energy  -= deltaEnergy;
00201 
00202     if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
00203       nSpotsInStep = maxNumberOfSpots - spotCounter;
00204       if ( nSpotsInStep < 1 ) { // @@ check
00205         edm::LogInfo("SimGeneralGFlash") << "GflashEMShowerProfile: Too Many Spots ";
00206         edm::LogInfo("SimGeneralGFlash") << " - break to regenerate nSpotsInStep ";
00207         break;
00208       }
00209     }
00210 
00211     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
00212     deltaStep  += 0.5*deltaZ;
00213     pathLength += deltaStep;
00214     deltaStep   =  0.5*deltaZ;
00215 
00216     // lateral shape and fluctuations for  homogenous calo.
00217     double tScale = tmax *alpha/(alpha-1.0) * (std::exp(fluctuatedAlpha)-1.0)/std::exp(fluctuatedAlpha);
00218     double tau = std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
00219     double rCore = z1 + z2 * tau; 
00220     double rTail = k1 *( std::exp(k3*(tau-k2)) + std::exp(k4*(tau-k2))); // @@ check RT3 sign
00221     double p23 = (p2 - tau)/p3;
00222     double probabilityWeight = p1 *  std::exp( p23 - std::exp(p23) );
00223 
00224     // Deposition of spots according to lateral distr.
00225     // Apply absolute energy scale
00226     // Convert into MeV unit
00227     double hitEnergy = deltaEnergy / nSpotsInStep * e25Scale * CLHEP::GeV;
00228     double hitTime = theShowino->getGlobalTime()*CLHEP::nanosecond + (pathLength - pathLength0)/30.0;
00229 
00230     GflashHit aHit;
00231 
00232     for (int ispot = 0 ;  ispot < nSpotsInStep ; ispot++) {
00233       spotCounter++;
00234 
00235       // Compute global position of generated spots with taking into account magnetic field
00236       // Divide deltaZ into nSpotsInStep and give a spot a global position
00237       double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
00238 
00239       // trajectoryPoint give a spot an imaginary point along the shower development
00240       GflashTrajectoryPoint trajectoryPoint;
00241       theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint,pathLength+incrementPath);
00242 
00243       double rShower = 0.0;
00244       Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint,rCore,rTail,probabilityWeight,rShower);
00245 
00246       // Convert into mm unit
00247       hitPosition *= CLHEP::cm;
00248 
00249       if( std::fabs(hitPosition.getZ()/CLHEP::cm) > Gflash::Zmax[Gflash::kHE]) continue;
00250 
00251       // put energy and position to a Hit
00252       aHit.setTime(hitTime);
00253       aHit.setEnergy(hitEnergy);
00254       aHit.setPosition(hitPosition);
00255       theGflashHitList.push_back(aHit);
00256 
00257       double zInX0_spot = std::abs(pathLength+incrementPath - pathLength0)/Gflash::radLength[jCalorimeter];
00258 
00259       nSpots_sd++;
00260 
00261       // for histogramming      
00262       if(theHisto->getStoreFlag()) {
00263         theHisto->em_long->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
00264         theHisto->em_lateral->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
00265         theHisto->em_long_sd->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
00266         theHisto->em_lateral_sd->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
00267       }
00268       
00269     } // end of for spot iteration
00270 
00271   } // end of while for longitudinal integration
00272 
00273   if(theHisto->getStoreFlag()) {
00274     theHisto->em_nSpots_sd->Fill(nSpots_sd);
00275   }
00276 
00277   //  delete theGflashNavigator;
00278 
00279 }
00280 
00281 double GflashEMShowerProfile::getDistanceToOut(Gflash::CalorimeterNumber kCalor) {
00282 
00283   double stepLengthLeft = 0.0;
00284   if(kCalor == Gflash::kESPM ) {
00285     stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kESPM])
00286                    - theShowino->getPathLengthAtShower();
00287   }
00288   else if (kCalor == Gflash::kENCA) {
00289     double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
00290     stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kENCA])
00291                    - theShowino->getPathLengthAtShower();
00292   }
00293   return stepLengthLeft;
00294 
00295 }
00296 
00297 Gflash3Vector GflashEMShowerProfile::locateHitPosition(GflashTrajectoryPoint& point, 
00298                                      double rCore, double rTail, double probability,double &rShower)
00299 {
00300   double u1 = CLHEP::HepUniformRand();
00301   double u2 = CLHEP::HepUniformRand();
00302   double rInRM = 0.0;
00303   
00304   if (u1 < probability ) {
00305     rInRM = rCore* std::sqrt( u2/(1.0-u2) );
00306   }
00307   else {
00308     rInRM = rTail * std::sqrt( u2/(1.0-u2) );
00309   }
00310   
00311   rShower =  rInRM * Gflash::rMoliere[jCalorimeter];
00312 
00313   // Uniform & random rotation of spot along the azimuthal angle
00314   double azimuthalAngle = CLHEP::twopi*CLHEP::HepUniformRand();
00315 
00316   // actual spot position by adding a radial vector to a trajectoryPoint
00317   Gflash3Vector position = point.getPosition() +
00318     rShower*std::cos(azimuthalAngle)*point.getOrthogonalUnitVector() +
00319     rShower*std::sin(azimuthalAngle)*point.getCrossUnitVector();
00320   
00321   return position;
00322 }
/data/refman/pasoursint/CMSSW_5_3_1/src/SimGeneral/GFlash/src/GflashEMShowerProfile.cc
