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::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  theBField = parSet.getParameter<double>("bField");
22  theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
23  theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");
24 
26 
27  theShowino = new GflashShowino();
28 
30 }
static GflashHistogram * instance()
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
GflashHistogram * theHisto
Gflash::CalorimeterNumber jCalorimeter
edm::ParameterSet theParSet

◆ ~GflashEMShowerProfile()

GflashEMShowerProfile::~GflashEMShowerProfile ( )

Definition at line 32 of file GflashEMShowerProfile.cc.

References theShowino.

32  {
33  if (theShowino)
34  delete theShowino;
35 }

Member Function Documentation

◆ getDistanceToOut()

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

Definition at line 279 of file GflashEMShowerProfile.cc.

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

Referenced by parameterization().

279  {
280  double stepLengthLeft = 0.0;
281  if (kCalor == Gflash::kESPM) {
284  } else if (kCalor == Gflash::kENCA) {
285  double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
286  stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign * Gflash::Zmax[Gflash::kENCA]) -
288  }
289  return stepLengthLeft;
290 }
const double Zmax[kNumberCalorimeter]
Gflash3Vector & getPosition()
Definition: GflashShowino.h:37
GflashTrajectory * getHelix()
Definition: GflashShowino.h:33
float getEta(float r, float z)
Definition: Hit.h:38
double getPathLengthAtRhoEquals(double rho) const
double getPathLengthAtShower()
Definition: GflashShowino.h:29
double getPathLengthAtZ(double z) const
const double Rmax[kNumberCalorimeter]

◆ getGflashHitList()

std::vector<GflashHit>& GflashEMShowerProfile::getGflashHitList ( )
inline

Definition at line 32 of file GflashEMShowerProfile.h.

References theGflashHitList.

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

32 { return theGflashHitList; };
std::vector< GflashHit > theGflashHitList

◆ getGflashShowino()

GflashShowino* GflashEMShowerProfile::getGflashShowino ( )
inline

Definition at line 31 of file GflashEMShowerProfile.h.

References theShowino.

31 { return theShowino; }

◆ initialize()

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

Definition at line 37 of file GflashEMShowerProfile.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, hcalRecHitTable_cff::energy, GflashShowino::initialize(), position, theBField, and theShowino.

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

38  {
39  theShowino->initialize(showerType, energy, globalTime, charge, position, momentum, theBField);
40 }
static int position[264][3]
Definition: ReadPGInfo.cc:289
void initialize(int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum, double magneticField)

◆ locateHitPosition()

Gflash3Vector GflashEMShowerProfile::locateHitPosition ( GflashTrajectoryPoint point,
double  rCore,
double  rTail,
double  probability,
double &  rShower 
)
private

Definition at line 292 of file GflashEMShowerProfile.cc.

References funct::cos(), jCalorimeter, point, position, Gflash::rMoliere, funct::sin(), mathSSE::sqrt(), testProducerWithPsetDescEmpty_cfi::u1, and MetAnalyzer::u2.

Referenced by parameterization().

293  {
294  double u1 = CLHEP::HepUniformRand();
295  double u2 = CLHEP::HepUniformRand();
296  double rInRM = 0.0;
297 
298  if (u1 < probability) {
299  rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
300  } else {
301  rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
302  }
303 
304  rShower = rInRM * Gflash::rMoliere[jCalorimeter];
305 
306  // Uniform & random rotation of spot along the azimuthal angle
307  double azimuthalAngle = CLHEP::twopi * CLHEP::HepUniformRand();
308 
309  // actual spot position by adding a radial vector to a trajectoryPoint
310  Gflash3Vector position = point.getPosition() + rShower * std::cos(azimuthalAngle) * point.getOrthogonalUnitVector() +
311  rShower * std::sin(azimuthalAngle) * point.getCrossUnitVector();
312 
313  return position;
314 }
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:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
CLHEP::Hep3Vector Gflash3Vector
Definition: Gflash3Vector.h:6
static int position[264][3]
Definition: ReadPGInfo.cc:289
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5

◆ parameterization()

void GflashEMShowerProfile::parameterization ( )

Definition at line 42 of file GflashEMShowerProfile.cc.

References funct::abs(), simBeamSpotPI::alpha, HLT_2023v12_cff::beta, Gflash::criticalEnergy, l1tTrackerHTMiss_cfi::deltaZ, 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, hcalRecHitTable_cff::energy, JetChargeProducer_cfi::exp, Gflash::getCalorimeterNumber(), getDistanceToOut(), GflashShowino::getEnergy(), GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getGlobalTime(), GflashShowino::getHelix(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPositionAtShower(), GflashHistogram::getStoreFlag(), jCalorimeter, Gflash::kENCA, Gflash::kESPM, Gflash::kHE, locateHitPosition(), dqm-mbProfile::log, listHistos::logY, SiStripPI::max, SiStripPI::min, LaserDQM_cfg::p1, SiStripOfflineCRack_cfg::p2, chargedHadronTrackResolutionFilter_cfi::p3, conifer::pow(), Gflash::radLength, rho, Gflash::rMoliere, mathSSE::sqrt(), metsig::tau, theEnergyScaleEB, theEnergyScaleEE, theGflashHitList, theHisto, theShowino, tmax, y, Gflash::Z, testProducerWithPsetDescEmpty_cfi::z2, and Gflash::Zmax.

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

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

Member Data Documentation

◆ jCalorimeter

Gflash::CalorimeterNumber GflashEMShowerProfile::jCalorimeter
private

◆ theBField

double GflashEMShowerProfile::theBField
private

Definition at line 43 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and initialize().

◆ theEnergyScaleEB

double GflashEMShowerProfile::theEnergyScaleEB
private

Definition at line 44 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

◆ theEnergyScaleEE

double GflashEMShowerProfile::theEnergyScaleEE
private

Definition at line 45 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

◆ theGflashHitList

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

Definition at line 49 of file GflashEMShowerProfile.h.

Referenced by getGflashHitList(), and parameterization().

◆ theHisto

GflashHistogram* GflashEMShowerProfile::theHisto
private

Definition at line 48 of file GflashEMShowerProfile.h.

Referenced by GflashEMShowerProfile(), and parameterization().

◆ theParSet

edm::ParameterSet GflashEMShowerProfile::theParSet
private

Definition at line 42 of file GflashEMShowerProfile.h.

◆ theShowino

GflashShowino* GflashEMShowerProfile::theShowino
private