CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
GflashEMShowerProfile Class Reference

#include <GflashEMShowerProfile.h>

Public Member Functions

std::vector< GflashHit > & getGflashHitList ()
 
GflashShowinogetGflashShowino ()
 
 GflashEMShowerProfile (const edm::ParameterSet &parSet)
 
void initialize (int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum)
 
void parameterization ()
 
 ~GflashEMShowerProfile ()
 

Private Member Functions

double getDistanceToOut (Gflash::CalorimeterNumber kCalor)
 
Gflash3Vector locateHitPosition (GflashTrajectoryPoint &point, double rCore, double rTail, double probability, double &rShower)
 

Private Attributes

Gflash::CalorimeterNumber jCalorimeter
 
double theBField
 
double theEnergyScaleEB
 
double theEnergyScaleEE
 
std::vector< GflashHittheGflashHitList
 
GflashHistogramtheHisto
 
edm::ParameterSet theParSet
 
GflashShowinotheShowino
 

Detailed Description

Definition at line 15 of file GflashEMShowerProfile.h.

Constructor & Destructor Documentation

GflashEMShowerProfile::GflashEMShowerProfile ( const edm::ParameterSet parSet)

Definition at line 20 of file GflashEMShowerProfile.cc.

References edm::ParameterSet::getParameter(), GflashHistogram::instance(), jCalorimeter, Gflash::kNULL, theBField, theEnergyScaleEB, theEnergyScaleEE, theHisto, and theShowino.

20  : theParSet(parSet)
21 {
22  theBField = parSet.getParameter<double>("bField");
23  theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
24  theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");
25 
27 
28  theShowino = new GflashShowino();
29 
31 
32 }
T getParameter(std::string const &) const
static GflashHistogram * instance()
GflashHistogram * theHisto
Gflash::CalorimeterNumber jCalorimeter
edm::ParameterSet theParSet
GflashEMShowerProfile::~GflashEMShowerProfile ( )

Definition at line 35 of file GflashEMShowerProfile.cc.

References theShowino.

36 {
37  if(theShowino) delete theShowino;
38 }

Member Function Documentation

double GflashEMShowerProfile::getDistanceToOut ( Gflash::CalorimeterNumber  kCalor)
private

Definition at line 280 of file GflashEMShowerProfile.cc.

References getEta(), GflashShowino::getHelix(), GflashTrajectory::getPathLengthAtRhoEquals(), GflashShowino::getPathLengthAtShower(), GflashTrajectory::getPathLengthAtZ(), GflashShowino::getPosition(), Gflash::kENCA, Gflash::kESPM, Gflash::Rmax, theShowino, and Gflash::Zmax.

Referenced by getGflashHitList(), and parameterization().

280  {
281 
282  double stepLengthLeft = 0.0;
283  if(kCalor == Gflash::kESPM ) {
286  }
287  else if (kCalor == Gflash::kENCA) {
288  double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
289  stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign*Gflash::Zmax[Gflash::kENCA])
291  }
292  return stepLengthLeft;
293 
294 }
double getPathLengthAtZ(double z) const
const double Zmax[kNumberCalorimeter]
Gflash3Vector & getPosition()
Definition: GflashShowino.h:34
GflashTrajectory * getHelix()
Definition: GflashShowino.h:30
static float getEta(float eta, int bits=5)
double getPathLengthAtShower()
Definition: GflashShowino.h:26
const double Rmax[kNumberCalorimeter]
double getPathLengthAtRhoEquals(double rho) const
std::vector<GflashHit>& GflashEMShowerProfile::getGflashHitList ( )
inline
GflashShowino* GflashEMShowerProfile::getGflashShowino ( )
inline

Definition at line 29 of file GflashEMShowerProfile.h.

References theShowino.

29 { return theShowino; }
void GflashEMShowerProfile::initialize ( int  showerType,
double  energy,
double  globalTime,
double  charge,
Gflash3Vector position,
Gflash3Vector momentum 
)

Definition at line 40 of file GflashEMShowerProfile.cc.

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

Referenced by GflashEMShowerModel::DoIt(), and GFlashEMShowerModel::DoIt().

42 {
43  theShowino->initialize(showerType, energy, globalTime, charge,
44  position, momentum, theBField);
45 }
static int position[264][3]
Definition: ReadPGInfo.cc:509
void initialize(int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum, double magneticField)
Gflash3Vector GflashEMShowerProfile::locateHitPosition ( GflashTrajectoryPoint point,
double  rCore,
double  rTail,
double  probability,
double &  rShower 
)
private

Definition at line 296 of file GflashEMShowerProfile.cc.

References funct::cos(), GflashTrajectoryPoint::getCrossUnitVector(), GflashTrajectoryPoint::getOrthogonalUnitVector(), GflashTrajectoryPoint::getPosition(), jCalorimeter, position, Gflash::rMoliere, funct::sin(), mathSSE::sqrt(), MetAnalyzer::u1, and MetAnalyzer::u2.

Referenced by getGflashHitList(), and parameterization().

298 {
299  double u1 = CLHEP::HepUniformRand();
300  double u2 = CLHEP::HepUniformRand();
301  double rInRM = 0.0;
302 
303  if (u1 < probability ) {
304  rInRM = rCore* std::sqrt( u2/(1.0-u2) );
305  }
306  else {
307  rInRM = rTail * std::sqrt( u2/(1.0-u2) );
308  }
309 
310  rShower = rInRM * Gflash::rMoliere[jCalorimeter];
311 
312  // Uniform & random rotation of spot along the azimuthal angle
313  double azimuthalAngle = CLHEP::twopi*CLHEP::HepUniformRand();
314 
315  // actual spot position by adding a radial vector to a trajectoryPoint
317  rShower*std::cos(azimuthalAngle)*point.getOrthogonalUnitVector() +
318  rShower*std::sin(azimuthalAngle)*point.getCrossUnitVector();
319 
320  return position;
321 }
Gflash3Vector getCrossUnitVector()
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Gflash::CalorimeterNumber jCalorimeter
const double rMoliere[kNumberCalorimeter]
T sqrt(T t)
Definition: SSEVec.h:18
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Gflash3Vector getOrthogonalUnitVector()
CLHEP::Hep3Vector Gflash3Vector
Definition: Gflash3Vector.h:6
Gflash3Vector & getPosition()
static int position[264][3]
Definition: ReadPGInfo.cc:509
void GflashEMShowerProfile::parameterization ( )

Definition at line 47 of file GflashEMShowerProfile.cc.

References funct::abs(), alpha, beta, Gflash::criticalEnergy, GflashHistogram::em_incE, GflashHistogram::em_lateral, GflashHistogram::em_lateral_sd, GflashHistogram::em_long, GflashHistogram::em_long_sd, GflashHistogram::em_nSpots_sd, GflashHistogram::em_ssp_rho, GflashHistogram::em_ssp_z, JetChargeProducer_cfi::exp, Gflash::getCalorimeterNumber(), getDistanceToOut(), GflashShowino::getEnergy(), GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getGlobalTime(), GflashShowino::getHelix(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPositionAtShower(), GflashHistogram::getStoreFlag(), GeV, jCalorimeter, Gflash::kENCA, Gflash::kESPM, Gflash::kHE, locateHitPosition(), cmsBatch::log, listHistos::logY, hpstanc_transforms::max, min(), p1, p2, p3, funct::pow(), Gflash::radLength, rho, Gflash::rMoliere, mathSSE::sqrt(), metsig::tau, theEnergyScaleEB, theEnergyScaleEE, theGflashHitList, theHisto, theShowino, tmax, y, Gflash::Z, and Gflash::Zmax.

Referenced by GflashEMShowerModel::DoIt(), and GFlashEMShowerModel::DoIt().

48 {
49  // This part of code is copied from the original GFlash Fortran code.
50  // reference : hep-ex/0001020v1
51  // The units used in Geant4 internally are in mm, MeV.
52  // For simplicity, the units here are in cm, GeV.
53 
54  const double energyCutoff = 0.01;
55  const int maxNumberOfSpots = 100000;
56 
57  double incomingEnergy = theShowino->getEnergy();
58  Gflash3Vector showerStartingPosition = theShowino->getPositionAtShower();
59 
60  //find the calorimeter at the shower starting point
61  jCalorimeter = Gflash::getCalorimeterNumber(showerStartingPosition);
62 
63  double logEinc = std::log(incomingEnergy);
64  double y = incomingEnergy / Gflash::criticalEnergy; // y = E/Ec, criticalEnergy is in GeV
65  double logY = std::log(y);
66 
67  // Total number of spots are not yet optimized.
68  double nSpots = 93.0 * std::log(Gflash::Z[jCalorimeter]) * std::pow(incomingEnergy,0.876);
69 
70  //path Length from the origin to the shower starting point in cm
71  double pathLength0 = theShowino->getPathLengthAtShower();
72  double pathLength = pathLength0; // this will grow along the shower development
73 
74  //--- 2.2 Fix intrinsic properties of em. showers.
75 
76  double fluctuatedTmax = std::log(logY - 0.7157);
77  double fluctuatedAlpha = std::log(0.7996 +(0.4581 + 1.8628/Gflash::Z[jCalorimeter])*logY);
78 
79  double sigmaTmax = 1.0/( -1.4 + 1.26 * logY);
80  double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
81  double rho = 0.705 - 0.023 * logY;
82  double sqrtPL = std::sqrt((1.0+rho)/2.0);
83  double sqrtLE = std::sqrt((1.0-rho)/2.0);
84 
85  double norm1 = CLHEP::RandGaussQ::shoot();
86  double norm2 = CLHEP::RandGaussQ::shoot();
87  double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
88  double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
89 
90  // tmax, alpha, beta : parameters of gamma distribution
91  double tmax = std::exp(tempTmax);
92  double alpha = std::exp(tempAlpha);
93  double beta = std::max(0.0,(alpha - 1.0)/tmax);
94 
95  // spot fluctuations are added to tmax, alpha, beta
96  double averageTmax = logY-0.858;
97  double averageAlpha = 0.21+(0.492+2.38/Gflash::Z[jCalorimeter])*logY;
98  double spotTmax = averageTmax * (0.698 + .00212*Gflash::Z[jCalorimeter]);
99  double spotAlpha = averageAlpha * (0.639 + .00334*Gflash::Z[jCalorimeter]);
100  double spotBeta = std::max(0.0,(spotAlpha-1.0)/spotTmax);
101 
102  if(theHisto->getStoreFlag()) {
103  theHisto->em_incE->Fill(incomingEnergy);
104  theHisto->em_ssp_rho->Fill(showerStartingPosition.rho());
105  theHisto->em_ssp_z->Fill(showerStartingPosition.z());
106  }
107 
108  // parameters for lateral distribution and fluctuation
109  double z1=0.0251+0.00319*logEinc;
110  double z2=0.1162-0.000381*Gflash::Z[jCalorimeter];
111 
112  double k1=0.659 - 0.00309 * Gflash::Z[jCalorimeter];
113  double k2=0.645;
114  double k3=-2.59;
115  double k4=0.3585+ 0.0421*logEinc;
116 
117  double p1=2.623 -0.00094*Gflash::Z[jCalorimeter];
118  double p2=0.401 +0.00187*Gflash::Z[jCalorimeter];
119  double p3=1.313 -0.0686*logEinc;
120 
121  // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : e25Scale = 1.006 for EB with ecalNotContainment = 1.0.
122  // Now e25Scale is a configurable parameter with default ecalNotContainment which is 0.97 for EB and 0.975 for EE.
123  // So if ecalNotContainment constants are to be changed in the future, e25Scale should be changed accordingly.
124  double e25Scale = 1.0;
125  if(jCalorimeter == Gflash::kESPM) e25Scale = theEnergyScaleEB;
126  else if(jCalorimeter == Gflash::kENCA) e25Scale = theEnergyScaleEE;
127 
128  // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : p1 *= 0.965
129  p1 *= 0.965;
130 
131  // preparation of longitudinal integration
132  double stepLengthLeft = getDistanceToOut(jCalorimeter);
133 
134  int nSpots_sd = 0; // count total number of spots in SD
135  double zInX0 = 0.0; // shower depth in X0 unit
136  double deltaZInX0 = 0.0; // segment of depth in X0 unit
137  double deltaZ = 0.0; // segment of depth in cm
138  double stepLengthLeftInX0 = 0.0; // step length left in X0 unit
139 
140  const double divisionStepInX0 = 0.1; //step size in X0 unit
141  double energy = incomingEnergy; // energy left in GeV
142 
143  Genfun::IncompleteGamma gammaDist;
144 
145  double energyInGamma = 0.0; // energy in a specific depth(z) according to Gamma distribution
146  double preEnergyInGamma = 0.0; // energy calculated in a previous depth
147  double sigmaInGamma = 0.; // sigma of energy in a specific depth(z) according to Gamma distribution
148  double preSigmaInGamma = 0.0; // sigma of energy in a previous depth
149 
150  //energy segment in Gamma distribution of shower in each step
151  double deltaEnergy =0.0 ; // energy in deltaZ
152  int spotCounter = 0; // keep track of number of spots generated
153 
154  //step increment along the shower direction
155  double deltaStep = 0.0;
156 
157  theGflashHitList.clear();
158 
159  // loop for longitudinal integration
160  while(energy > 0.0 && stepLengthLeft > 0.0) {
161 
162  stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
163 
164  if ( stepLengthLeftInX0 < divisionStepInX0 ) {
165  deltaZInX0 = stepLengthLeftInX0;
166  deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
167  stepLengthLeft = 0.0;
168  }
169  else {
170  deltaZInX0 = divisionStepInX0;
171  deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
172  stepLengthLeft -= deltaZ;
173  }
174 
175  zInX0 += deltaZInX0;
176 
177  int nSpotsInStep = 0;
178 
179  if ( energy > energyCutoff ) {
180  preEnergyInGamma = energyInGamma;
181  gammaDist.a().setValue(alpha); //alpha
182  energyInGamma = gammaDist(beta*zInX0); //beta
183  double energyInDeltaZ = energyInGamma - preEnergyInGamma;
184  deltaEnergy = std::min(energy,incomingEnergy*energyInDeltaZ);
185 
186  preSigmaInGamma = sigmaInGamma;
187  gammaDist.a().setValue(spotAlpha); //alpha spot
188  sigmaInGamma = gammaDist(spotBeta*zInX0); //beta spot
189  nSpotsInStep = std::max(1,int(nSpots * (sigmaInGamma - preSigmaInGamma)));
190  }
191  else {
192  deltaEnergy = energy;
193  preSigmaInGamma = sigmaInGamma;
194  nSpotsInStep = std::max(1,int(nSpots * (1.0 - preSigmaInGamma)));
195  }
196 
197  if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy = energy;
198 
199  energy -= deltaEnergy;
200 
201  if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
202  nSpotsInStep = maxNumberOfSpots - spotCounter;
203  if ( nSpotsInStep < 1 ) { // @@ check
204  edm::LogInfo("SimGeneralGFlash") << "GflashEMShowerProfile: Too Many Spots ";
205  edm::LogInfo("SimGeneralGFlash") << " - break to regenerate nSpotsInStep ";
206  break;
207  }
208  }
209 
210  // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
211  deltaStep += 0.5*deltaZ;
212  pathLength += deltaStep;
213  deltaStep = 0.5*deltaZ;
214 
215  // lateral shape and fluctuations for homogenous calo.
216  double tScale = tmax *alpha/(alpha-1.0) * (std::exp(fluctuatedAlpha)-1.0)/std::exp(fluctuatedAlpha);
217  double tau = std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
218  double rCore = z1 + z2 * tau;
219  double rTail = k1 *( std::exp(k3*(tau-k2)) + std::exp(k4*(tau-k2))); // @@ check RT3 sign
220  double p23 = (p2 - tau)/p3;
221  double probabilityWeight = p1 * std::exp( p23 - std::exp(p23) );
222 
223  // Deposition of spots according to lateral distr.
224  // Apply absolute energy scale
225  // Convert into MeV unit
226  double hitEnergy = deltaEnergy / nSpotsInStep * e25Scale * CLHEP::GeV;
227  double hitTime = theShowino->getGlobalTime()*CLHEP::nanosecond + (pathLength - pathLength0)/30.0;
228 
229  GflashHit aHit;
230 
231  for (int ispot = 0 ; ispot < nSpotsInStep ; ispot++) {
232  spotCounter++;
233 
234  // Compute global position of generated spots with taking into account magnetic field
235  // Divide deltaZ into nSpotsInStep and give a spot a global position
236  double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
237 
238  // trajectoryPoint give a spot an imaginary point along the shower development
239  GflashTrajectoryPoint trajectoryPoint;
240  theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint,pathLength+incrementPath);
241 
242  double rShower = 0.0;
243  Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint,rCore,rTail,probabilityWeight,rShower);
244 
245  // Convert into mm unit
246  hitPosition *= CLHEP::cm;
247 
248  if( std::fabs(hitPosition.getZ()/CLHEP::cm) > Gflash::Zmax[Gflash::kHE]) continue;
249 
250  // put energy and position to a Hit
251  aHit.setTime(hitTime);
252  aHit.setEnergy(hitEnergy);
253  aHit.setPosition(hitPosition);
254  theGflashHitList.push_back(aHit);
255 
256  double zInX0_spot = std::abs(pathLength+incrementPath - pathLength0)/Gflash::radLength[jCalorimeter];
257 
258  nSpots_sd++;
259 
260  // for histogramming
261  if(theHisto->getStoreFlag()) {
262  theHisto->em_long->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
263  theHisto->em_lateral->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
264  theHisto->em_long_sd->Fill(zInX0_spot,hitEnergy/CLHEP::GeV);
265  theHisto->em_lateral_sd->Fill(zInX0_spot,rShower/Gflash::rMoliere[jCalorimeter],hitEnergy/CLHEP::GeV);
266  }
267 
268  } // end of for spot iteration
269 
270  } // end of while for longitudinal integration
271 
272  if(theHisto->getStoreFlag()) {
273  theHisto->em_nSpots_sd->Fill(nSpots_sd);
274  }
275 
276  // delete theGflashNavigator;
277 
278 }
const double Z[kNumberCalorimeter]
const double beta
float alpha
Definition: AMPTWrapper.h:95
const double GeV
Definition: MathUtil.h:16
const double Zmax[kNumberCalorimeter]
GflashTrajectory * getHelix()
Definition: GflashShowino.h:30
double getEnergy()
Definition: GflashShowino.h:24
double getDistanceToOut(Gflash::CalorimeterNumber kCalor)
GflashHistogram * theHisto
std::vector< GflashHit > theGflashHitList
Gflash::CalorimeterNumber jCalorimeter
const double rMoliere[kNumberCalorimeter]
CalorimeterNumber getCalorimeterNumber(const Gflash3Vector &position)
double getPathLengthAtShower()
Definition: GflashShowino.h:26
T sqrt(T t)
Definition: SSEVec.h:18
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const double criticalEnergy
T min(T a, T b)
Definition: MathUtil.h:58
Gflash3Vector & getPositionAtShower()
Definition: GflashShowino.h:27
double p2[4]
Definition: TauolaWrapper.h:90
void getGflashTrajectoryPoint(GflashTrajectoryPoint &point, double s) const
Gflash3Vector locateHitPosition(GflashTrajectoryPoint &point, double rCore, double rTail, double probability, double &rShower)
static const double tmax[3]
const double radLength[kNumberCalorimeter]
CLHEP::Hep3Vector Gflash3Vector
Definition: Gflash3Vector.h:6
double getGlobalTime()
Definition: GflashShowino.h:32
double p1[4]
Definition: TauolaWrapper.h:89
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double p3[4]
Definition: TauolaWrapper.h:91

Member Data Documentation

Gflash::CalorimeterNumber GflashEMShowerProfile::jCalorimeter
private
double GflashEMShowerProfile::theBField
private

Definition at line 42 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and initialize().

double GflashEMShowerProfile::theEnergyScaleEB
private

Definition at line 43 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

double GflashEMShowerProfile::theEnergyScaleEE
private

Definition at line 44 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

std::vector<GflashHit> GflashEMShowerProfile::theGflashHitList
private

Definition at line 48 of file GflashEMShowerProfile.h.

Referenced by getGflashHitList(), and parameterization().

GflashHistogram* GflashEMShowerProfile::theHisto
private

Definition at line 47 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

edm::ParameterSet GflashEMShowerProfile::theParSet
private

Definition at line 41 of file GflashEMShowerProfile.h.

GflashShowino* GflashEMShowerProfile::theShowino
private