GFlash/src/GflashPiKShowerProfile.cc

0001 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
0002 #include <CLHEP/Random/RandGaussQ.h>
0003
0004 void GflashPiKShowerProfile::loadParameters() {
0005   double einc = theShowino->getEnergy();
0006   Gflash3Vector position = theShowino->getPositionAtShower();
0007   int showerType = theShowino->getShowerType();
0008
0009   // energy scale
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     //@@@ energy dependent energyRho based on tuning with testbeam data
0018     double energyRho = fTanh(einc, Gflash::pion_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::pion_emscale[0]) + fTanh(einc, Gflash::pion_emscale[1]) * r1);
0026
0027         // LogNormal mean and sigma of Hcal energy
0028         energyMeanHcal = (fTanh(einc, Gflash::pion_hadscale[0]) +
0029                           fTanh(einc, Gflash::pion_hadscale[1]) *
0030                               depthScale(position.getRho(), Gflash::RFrontCrystalEB, Gflash::LengthCrystalEB));
0031         energySigmaHcal = (fTanh(einc, Gflash::pion_hadscale[2]) +
0032                            fTanh(einc, Gflash::pion_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::pion_emscale[0]) + fTanh(einc, Gflash::pion_emscale[1]) * r1);
0044
0045         //@@@extend depthScale for HE
0046         energyMeanHcal = (fTanh(einc, Gflash::pion_hadscale[0]) +
0047                           fTanh(einc, Gflash::pion_hadscale[1]) *
0048                               depthScale(std::fabs(position.getZ()), Gflash::ZFrontCrystalEE, Gflash::LengthCrystalEE));
0049         energySigmaHcal =
0050             (fTanh(einc, Gflash::pion_hadscale[2]) +
0051              fTanh(einc, Gflash::pion_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     // Hcal response for mip-like pions (mip)
0060
0061     energyMeanHcal = fTanh(einc, Gflash::pion_hadscale[4]);
0062     energySigmaHcal = fTanh(einc, Gflash::pion_hadscale[5]);
0063
0064     double gap_corr = einc * fTanh(einc, Gflash::pion_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::pion_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   // parameters for the longitudinal profiles
0097   //@@@check longitudinal profiles of endcaps for possible variations
0098
0099   double *rhoHcal = new double[2 * Gflash::NPar];
0100   double *correlationVectorHcal = new double[Gflash::NPar * (Gflash::NPar + 1) / 2];
0101
0102   //@@@until we have a separate parameterization for Endcap
0103
0104   if (showerType > 3) {
0105     showerType -= 4;
0106   }
0107   // no separate parameterization before crystal
0108   if (showerType == 0)
0109     showerType = 1;
0110
0111   // Hcal parameters are always needed regardless of showerType
0112   for (int i = 0; i < 2 * Gflash::NPar; i++) {
0113     rhoHcal[i] = fTanh(einc, Gflash::pion_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::pion_par[i + showerType * Gflash::NPar]) +
0129                   fTanh(einc, Gflash::pion_par[i + (4 + showerType) * Gflash::NPar]) * correlationSum;
0130   }
0131   delete[] rhoHcal;
0132   delete[] correlationVectorHcal;
0133
0134   // lateral parameters for Hcal
0135
0136   for (int i = 0; i < Gflash::Nrpar; i++) {
0137     lateralPar[Gflash::kHB][i] = fTanh(einc, Gflash::pion_rpar[i + showerType * Gflash::Nrpar]);
0138     lateralPar[Gflash::kHE][i] = lateralPar[Gflash::kHB][i];
0139   }
0140
0141   // Ecal parameters are needed if and only if the shower starts inside the
0142   // crystal
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::pion_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] =
0161           fTanh(einc, Gflash::pion_par[i]) + 0.5 * fTanh(einc, Gflash::pion_par[i + 4 * Gflash::NPar]) * correlationSum;
0162     }
0163
0164     delete[] rhoEcal;
0165     delete[] correlationVectorEcal;
0166
0167     // lateral parameters for Ecal
0168
0169     for (int i = 0; i < Gflash::Nrpar; i++) {
0170       lateralPar[Gflash::kESPM][i] = fTanh(einc, Gflash::pion_rpar[i]);
0171       lateralPar[Gflash::kENCA][i] = lateralPar[Gflash::kESPM][i];
0172     }
0173   }
0174 }