File indexing completed on 2024-04-06 12:30:35
0001 #include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h"
0002 #include <CLHEP/Random/RandGaussQ.h>
0003
0004 void GflashKaonPlusShowerProfile::loadParameters() {
0005 double einc = theShowino->getEnergy();
0006 Gflash3Vector position = theShowino->getPositionAtShower();
0007 int showerType = theShowino->getShowerType();
0008
0009
0010 double energyMeanHcal = 0.0;
0011 double energySigmaHcal = 0.0;
0012
0013 double r1 = 0.0;
0014 double r2 = 0.0;
0015
0016 if (showerType == 0 || showerType == 1 || showerType == 4 || showerType == 5) {
0017
0018 double energyRho = fTanh(einc, Gflash::kplus_correl_hadem);
0019
0020 if (showerType == 0 || showerType == 1) {
0021 do {
0022 r1 = CLHEP::RandGaussQ::shoot();
0023
0024 energyScale[Gflash::kESPM] =
0025 einc * (fTanh(einc, Gflash::kplus_emscale[0]) + fTanh(einc, Gflash::kplus_emscale[1]) * r1);
0026
0027
0028 energyMeanHcal = (fTanh(einc, Gflash::kplus_hadscale[0]) +
0029 fTanh(einc, Gflash::kplus_hadscale[1]) *
0030 depthScale(position.getRho(), Gflash::RFrontCrystalEB, Gflash::LengthCrystalEB));
0031 energySigmaHcal = (fTanh(einc, Gflash::kplus_hadscale[2]) +
0032 fTanh(einc, Gflash::kplus_hadscale[3]) *
0033 depthScale(position.getRho(), Gflash::RFrontCrystalEB, Gflash::LengthCrystalEB));
0034
0035 r2 = CLHEP::RandGaussQ::shoot();
0036 energyScale[Gflash::kHB] =
0037 exp(energyMeanHcal + energySigmaHcal * (energyRho * r1 + sqrt(1.0 - energyRho * energyRho) * r2));
0038 } while (energyScale[Gflash::kESPM] < 0 || energyScale[Gflash::kHB] > einc * 1.5);
0039 } else {
0040 do {
0041 r1 = CLHEP::RandGaussQ::shoot();
0042 energyScale[Gflash::kENCA] =
0043 einc * (fTanh(einc, Gflash::kplus_emscale[0]) + fTanh(einc, Gflash::kplus_emscale[1]) * r1);
0044
0045
0046 energyMeanHcal = (fTanh(einc, Gflash::kplus_hadscale[0]) +
0047 fTanh(einc, Gflash::kplus_hadscale[1]) *
0048 depthScale(std::fabs(position.getZ()), Gflash::ZFrontCrystalEE, Gflash::LengthCrystalEE));
0049 energySigmaHcal =
0050 (fTanh(einc, Gflash::kplus_hadscale[2]) +
0051 fTanh(einc, Gflash::kplus_hadscale[3]) *
0052 depthScale(std::fabs(position.getZ()), Gflash::ZFrontCrystalEE, Gflash::LengthCrystalEE));
0053 r2 = CLHEP::RandGaussQ::shoot();
0054 energyScale[Gflash::kHE] =
0055 exp(energyMeanHcal + energySigmaHcal * (energyRho * r1 + sqrt(1.0 - energyRho * energyRho) * r2));
0056 } while (energyScale[Gflash::kENCA] < 0 || energyScale[Gflash::kHE] > einc * 1.5);
0057 }
0058 } else if (showerType == 2 || showerType == 3 || showerType == 6 || showerType == 7) {
0059
0060
0061 energyMeanHcal = fTanh(einc, Gflash::kplus_hadscale[4]);
0062 energySigmaHcal = fTanh(einc, Gflash::kplus_hadscale[5]);
0063
0064 double gap_corr = einc * fTanh(einc, Gflash::kplus_hadscale[6]);
0065
0066 if (showerType == 2 || showerType == 3) {
0067 energyScale[Gflash::kESPM] = 0.0;
0068
0069 do {
0070 r1 = CLHEP::RandGaussQ::shoot();
0071 energyScale[Gflash::kHB] = exp(energyMeanHcal + energySigmaHcal * r1);
0072 } while (energyScale[Gflash::kHB] > einc * 1.5);
0073
0074 if (showerType == 2) {
0075 energyScale[Gflash::kHE] = std::max(
0076 0.0, energyScale[Gflash::kHB] - gap_corr * depthScale(position.getRho(), Gflash::Rmin[Gflash::kHB], 28.));
0077 }
0078 } else if (showerType == 6 || showerType == 7) {
0079 energyScale[Gflash::kENCA] = 0.0;
0080
0081 do {
0082 r1 = CLHEP::RandGaussQ::shoot();
0083 energyMeanHcal += std::log(1.0 - fTanh(einc, Gflash::kplus_hadscale[7]));
0084 energyScale[Gflash::kHE] = exp(energyMeanHcal + energySigmaHcal * r1);
0085 } while (energyScale[Gflash::kHE] > einc * 1.5);
0086
0087 if (showerType == 6) {
0088 energyScale[Gflash::kHE] =
0089 std::max(0.0,
0090 energyScale[Gflash::kHE] -
0091 gap_corr * depthScale(std::fabs(position.getZ()), Gflash::Zmin[Gflash::kHE], 66.));
0092 }
0093 }
0094 }
0095
0096
0097
0098
0099 double *rhoHcal = new double[2 * Gflash::NPar];
0100 double *correlationVectorHcal = new double[Gflash::NPar * (Gflash::NPar + 1) / 2];
0101
0102
0103
0104 if (showerType > 3) {
0105 showerType -= 4;
0106 }
0107
0108 if (showerType == 0)
0109 showerType = 1;
0110
0111
0112 for (int i = 0; i < 2 * Gflash::NPar; i++) {
0113 rhoHcal[i] = fTanh(einc, Gflash::kplus_rho[i + showerType * 2 * Gflash::NPar]);
0114 }
0115
0116 getFluctuationVector(rhoHcal, correlationVectorHcal);
0117
0118 double normalZ[Gflash::NPar];
0119 for (int i = 0; i < Gflash::NPar; i++)
0120 normalZ[i] = CLHEP::RandGaussQ::shoot();
0121
0122 for (int i = 0; i < Gflash::NPar; i++) {
0123 double correlationSum = 0.0;
0124
0125 for (int j = 0; j < i + 1; j++) {
0126 correlationSum += correlationVectorHcal[i * (i + 1) / 2 + j] * normalZ[j];
0127 }
0128 longHcal[i] = fTanh(einc, Gflash::kplus_par[i + showerType * Gflash::NPar]) +
0129 fTanh(einc, Gflash::kplus_par[i + (4 + showerType) * Gflash::NPar]) * correlationSum;
0130 }
0131 delete[] rhoHcal;
0132 delete[] correlationVectorHcal;
0133
0134
0135
0136 for (int i = 0; i < Gflash::Nrpar; i++) {
0137 lateralPar[Gflash::kHB][i] = fTanh(einc, Gflash::kplus_rpar[i + showerType * Gflash::Nrpar]);
0138 lateralPar[Gflash::kHE][i] = lateralPar[Gflash::kHB][i];
0139 }
0140
0141
0142
0143
0144 if (showerType == 1) {
0145 double *rhoEcal = new double[2 * Gflash::NPar];
0146 double *correlationVectorEcal = new double[Gflash::NPar * (Gflash::NPar + 1) / 2];
0147 for (int i = 0; i < 2 * Gflash::NPar; i++)
0148 rhoEcal[i] = fTanh(einc, Gflash::kplus_rho[i]);
0149
0150 getFluctuationVector(rhoEcal, correlationVectorEcal);
0151
0152 for (int i = 0; i < Gflash::NPar; i++)
0153 normalZ[i] = CLHEP::RandGaussQ::shoot();
0154 for (int i = 0; i < Gflash::NPar; i++) {
0155 double correlationSum = 0.0;
0156
0157 for (int j = 0; j < i + 1; j++) {
0158 correlationSum += correlationVectorEcal[i * (i + 1) / 2 + j] * normalZ[j];
0159 }
0160 longEcal[i] = fTanh(einc, Gflash::kplus_par[i]) +
0161 0.5 * fTanh(einc, Gflash::kplus_par[i + 4 * Gflash::NPar]) * correlationSum;
0162 }
0163
0164 delete[] rhoEcal;
0165 delete[] correlationVectorEcal;
0166
0167
0168
0169 for (int i = 0; i < Gflash::Nrpar; i++) {
0170 lateralPar[Gflash::kESPM][i] = fTanh(einc, Gflash::kplus_rpar[i]);
0171 lateralPar[Gflash::kENCA][i] = lateralPar[Gflash::kESPM][i];
0172 }
0173 }
0174 }