CMS 3D CMS Logo

GflashKaonPlusShowerProfile.cc
Go to the documentation of this file.
2 #include <CLHEP/Random/RandGaussQ.h>
3 
5  double einc = theShowino->getEnergy();
7  int showerType = theShowino->getShowerType();
8 
9  // energy scale
10  double energyMeanHcal = 0.0;
11  double energySigmaHcal = 0.0;
12 
13  double r1 = 0.0;
14  double r2 = 0.0;
15 
16  if (showerType == 0 || showerType == 1 || showerType == 4 || showerType == 5) {
17  //@@@ energy dependent energyRho based on tuning with testbeam data
18  double energyRho = fTanh(einc, Gflash::kplus_correl_hadem);
19 
20  if (showerType == 0 || showerType == 1) {
21  do {
22  r1 = CLHEP::RandGaussQ::shoot();
23 
25  einc * (fTanh(einc, Gflash::kplus_emscale[0]) + fTanh(einc, Gflash::kplus_emscale[1]) * r1);
26 
27  // LogNormal mean and sigma of Hcal energy
28  energyMeanHcal = (fTanh(einc, Gflash::kplus_hadscale[0]) +
29  fTanh(einc, Gflash::kplus_hadscale[1]) *
31  energySigmaHcal = (fTanh(einc, Gflash::kplus_hadscale[2]) +
32  fTanh(einc, Gflash::kplus_hadscale[3]) *
34 
35  r2 = CLHEP::RandGaussQ::shoot();
37  exp(energyMeanHcal + energySigmaHcal * (energyRho * r1 + sqrt(1.0 - energyRho * energyRho) * r2));
38  } while (energyScale[Gflash::kESPM] < 0 || energyScale[Gflash::kHB] > einc * 1.5);
39  } else {
40  do {
41  r1 = CLHEP::RandGaussQ::shoot();
43  einc * (fTanh(einc, Gflash::kplus_emscale[0]) + fTanh(einc, Gflash::kplus_emscale[1]) * r1);
44 
45  //@@@extend depthScale for HE
46  energyMeanHcal = (fTanh(einc, Gflash::kplus_hadscale[0]) +
47  fTanh(einc, Gflash::kplus_hadscale[1]) *
48  depthScale(std::fabs(position.getZ()), Gflash::ZFrontCrystalEE, Gflash::LengthCrystalEE));
49  energySigmaHcal =
50  (fTanh(einc, Gflash::kplus_hadscale[2]) +
51  fTanh(einc, Gflash::kplus_hadscale[3]) *
52  depthScale(std::fabs(position.getZ()), Gflash::ZFrontCrystalEE, Gflash::LengthCrystalEE));
53  r2 = CLHEP::RandGaussQ::shoot();
55  exp(energyMeanHcal + energySigmaHcal * (energyRho * r1 + sqrt(1.0 - energyRho * energyRho) * r2));
56  } while (energyScale[Gflash::kENCA] < 0 || energyScale[Gflash::kHE] > einc * 1.5);
57  }
58  } else if (showerType == 2 || showerType == 3 || showerType == 6 || showerType == 7) {
59  // Hcal response for mip-like pions (mip)
60 
61  energyMeanHcal = fTanh(einc, Gflash::kplus_hadscale[4]);
62  energySigmaHcal = fTanh(einc, Gflash::kplus_hadscale[5]);
63 
64  double gap_corr = einc * fTanh(einc, Gflash::kplus_hadscale[6]);
65 
66  if (showerType == 2 || showerType == 3) {
68 
69  do {
70  r1 = CLHEP::RandGaussQ::shoot();
71  energyScale[Gflash::kHB] = exp(energyMeanHcal + energySigmaHcal * r1);
72  } while (energyScale[Gflash::kHB] > einc * 1.5);
73 
74  if (showerType == 2) {
76  0.0, energyScale[Gflash::kHB] - gap_corr * depthScale(position.getRho(), Gflash::Rmin[Gflash::kHB], 28.));
77  }
78  } else if (showerType == 6 || showerType == 7) {
80 
81  do {
82  r1 = CLHEP::RandGaussQ::shoot();
83  energyMeanHcal += std::log(1.0 - fTanh(einc, Gflash::kplus_hadscale[7]));
84  energyScale[Gflash::kHE] = exp(energyMeanHcal + energySigmaHcal * r1);
85  } while (energyScale[Gflash::kHE] > einc * 1.5);
86 
87  if (showerType == 6) {
89  std::max(0.0,
91  gap_corr * depthScale(std::fabs(position.getZ()), Gflash::Zmin[Gflash::kHE], 66.));
92  }
93  }
94  }
95 
96  // parameters for the longitudinal profiles
97  //@@@check longitudinal profiles of endcaps for possible variations
98 
99  double *rhoHcal = new double[2 * Gflash::NPar];
100  double *correlationVectorHcal = new double[Gflash::NPar * (Gflash::NPar + 1) / 2];
101 
102  //@@@until we have a separate parameterization for Endcap
103 
104  if (showerType > 3) {
105  showerType -= 4;
106  }
107  // no separate parameterization before crystal
108  if (showerType == 0)
109  showerType = 1;
110 
111  // Hcal parameters are always needed regardless of showerType
112  for (int i = 0; i < 2 * Gflash::NPar; i++) {
113  rhoHcal[i] = fTanh(einc, Gflash::kplus_rho[i + showerType * 2 * Gflash::NPar]);
114  }
115 
116  getFluctuationVector(rhoHcal, correlationVectorHcal);
117 
118  double normalZ[Gflash::NPar];
119  for (int i = 0; i < Gflash::NPar; i++)
120  normalZ[i] = CLHEP::RandGaussQ::shoot();
121 
122  for (int i = 0; i < Gflash::NPar; i++) {
123  double correlationSum = 0.0;
124 
125  for (int j = 0; j < i + 1; j++) {
126  correlationSum += correlationVectorHcal[i * (i + 1) / 2 + j] * normalZ[j];
127  }
128  longHcal[i] = fTanh(einc, Gflash::kplus_par[i + showerType * Gflash::NPar]) +
129  fTanh(einc, Gflash::kplus_par[i + (4 + showerType) * Gflash::NPar]) * correlationSum;
130  }
131  delete[] rhoHcal;
132  delete[] correlationVectorHcal;
133 
134  // lateral parameters for Hcal
135 
136  for (int i = 0; i < Gflash::Nrpar; i++) {
137  lateralPar[Gflash::kHB][i] = fTanh(einc, Gflash::kplus_rpar[i + showerType * Gflash::Nrpar]);
139  }
140 
141  // Ecal parameters are needed if and only if the shower starts inside the
142  // crystal
143 
144  if (showerType == 1) {
145  double *rhoEcal = new double[2 * Gflash::NPar];
146  double *correlationVectorEcal = new double[Gflash::NPar * (Gflash::NPar + 1) / 2];
147  for (int i = 0; i < 2 * Gflash::NPar; i++)
148  rhoEcal[i] = fTanh(einc, Gflash::kplus_rho[i]);
149 
150  getFluctuationVector(rhoEcal, correlationVectorEcal);
151 
152  for (int i = 0; i < Gflash::NPar; i++)
153  normalZ[i] = CLHEP::RandGaussQ::shoot();
154  for (int i = 0; i < Gflash::NPar; i++) {
155  double correlationSum = 0.0;
156 
157  for (int j = 0; j < i + 1; j++) {
158  correlationSum += correlationVectorEcal[i * (i + 1) / 2 + j] * normalZ[j];
159  }
160  longEcal[i] = fTanh(einc, Gflash::kplus_par[i]) +
161  0.5 * fTanh(einc, Gflash::kplus_par[i + 4 * Gflash::NPar]) * correlationSum;
162  }
163 
164  delete[] rhoEcal;
165  delete[] correlationVectorEcal;
166 
167  // lateral parameters for Ecal
168 
169  for (int i = 0; i < Gflash::Nrpar; i++) {
172  }
173  }
174 }
const double ZFrontCrystalEE
double fTanh(double einc, const double *par)
const double kplus_hadscale[8][5]
double lateralPar[Gflash::kNumberCalorimeter][Gflash::Nrpar]
double getEnergy()
Definition: GflashShowino.h:27
int getShowerType()
Definition: GflashShowino.h:26
const double kplus_correl_hadem[5]
const int Nrpar
const double kplus_rho[8 *NPar][5]
T sqrt(T t)
Definition: SSEVec.h:18
const double LengthCrystalEE
const double LengthCrystalEB
Gflash3Vector & getPositionAtShower()
Definition: GflashShowino.h:30
const int NPar
const double kplus_rpar[4 *Nrpar][5]
const double RFrontCrystalEB
const double Zmin[kNumberCalorimeter]
CLHEP::Hep3Vector Gflash3Vector
Definition: Gflash3Vector.h:6
const double kplus_emscale[2][5]
const double kplus_par[8 *NPar][5]
static int position[264][3]
Definition: ReadPGInfo.cc:509
void getFluctuationVector(double *lowTriangle, double *correlationVector)
double depthScale(double ssp, double ssp0, double length)
const double Rmin[kNumberCalorimeter]
double energyScale[Gflash::kNumberCalorimeter]