ShowerDevelopment/src/HDRShower.cc

0001 //FastSimulation Headers
0002 #include "FastSimulation/ShowerDevelopment/interface/HDRShower.h"
0003 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0004 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0005 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0006
0007 //CMSSW headers
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 using namespace edm;
0011
0012 ////////////////////////////////////////////////////////////////////////////////
0013 // What's this? Doesn't seem to be needed. Maybe Geometry/CaloGeometry/interface/CaloCellGeometry.h?
0014 //#include "Calorimetry/CaloDetector/interface/CellGeometry.h"
0015
0016 // number attempts for transverse distribution if exit on a spec. condition
0017 #define infinity 5000
0018 // debugging flag ( 0, 1, 2, 3)
0019 #define debug 0
0020
0021 using namespace std;
0022
0023 HDRShower::HDRShower(const RandomEngineAndDistribution* engine,
0024                      HDShowerParametrization* myParam,
0025                      EcalHitMaker* myGrid,
0026                      HcalHitMaker* myHcalHitMaker,
0027                      int onECAL,
0028                      double epart)
0029     : theParam(myParam), theGrid(myGrid), theHcalHitMaker(myHcalHitMaker), onEcal(onECAL), e(epart), random(engine) {
0030   eHDspot = 0.2;
0031   EsCut = 0.050;
0032   EcalShift = 0.12;
0033   nthetaStep = 10;
0034   thetaStep = 0.5 * M_PI / nthetaStep;
0035
0036   if (e < 0)
0037     e = 0.;
0038   setFuncParam();
0039 }
0040
0041 bool HDRShower::computeShower() {
0042   if (onEcal) {
0043     depthECAL = theGrid->ecalTotalL0();        // ECAL depth segment
0044     depthGAP = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
0045   } else
0046     depthECAL = depthGAP = 0;
0047
0048   float depthHCAL = theGrid->hcalTotalL0();  // HCAL depth segment
0049
0050   //  maxDepth   = depthECAL + depthGAP + depthHCAL - 1.0;
0051   maxDepth = depthECAL + depthHCAL - 0.5;
0052   depthStart = log(1. / random->flatShoot());  // starting point lambda unts
0053
0054   if (depthStart > maxDepth) {
0055     depthStart = maxDepth * random->flatShoot();
0056     if (depthStart < 0.)
0057       depthStart = 0.;
0058   }
0059
0060   if (depthStart < EcalShift)
0061     depthStart = EcalShift;
0062
0063   decal = (depthECAL + depthStart) * 0.5;
0064   qstatus = false;
0065   if (decal < depthECAL) {
0066     qstatus = theGrid->getPads(decal);
0067     //    if(!qstatus)
0068     //      cout<<" depth rejected by getQuads(decal="<<decal<<") status="<<qstatus
0069     //    <<" depthECAL="<<depthECAL<<endl;
0070   }
0071
0072   thetaFunction(nthetaStep);
0073   int maxLoops = 10000;
0074   for (int itheta = 0; itheta < nthetaStep; itheta++) {
0075     float theta, es;
0076     for (int i = 0; i <= thetaSpots[itheta]; i++) {
0077       if (i == thetaSpots[itheta])
0078         es = elastspot[itheta];
0079       else
0080         es = eHDspot;
0081       for (int j = 0; j < maxLoops; j++) {
0082         theta = (itheta + random->flatShoot()) * thetaStep;
0083         if (setHit(es, theta))
0084           break;
0085       }
0086     }
0087   }
0088   return (true);
0089 }
0090
0091 bool HDRShower::setHit(float espot, float theta) {
0092   float phi = 2. * M_PI * random->flatShoot();  // temporary: 1st approximation
0093   float rshower = getR();                       // temporary: 1st approximation
0094
0095   float d = depthStart + rshower * cos(theta);
0096   if (d + depthGAP > maxDepth)
0097     return (false);
0098
0099   // Commented (F.B) to remove a warning. Not used anywhere ?
0100   //  bool inHcal = !onEcal || d>depthECAL || !qstatus;
0101   bool result = false;
0102   if (!onEcal || d > depthECAL || !qstatus) {  // in HCAL (HF or HB, HE)
0103     d += depthGAP;
0104     bool setHDdepth = theHcalHitMaker->setDepth(d);
0105     if (setHDdepth) {
0106       theHcalHitMaker->setSpotEnergy(espot);
0107       result = theHcalHitMaker->addHit(rshower * sin(theta), phi, 0);
0108     } else
0109       LogWarning("FastCalorimetry") << " setHit in HCAL failed d=" << d << " maxDepth=" << maxDepth << " onEcal'"
0110                                     << onEcal << endl;
0111   } else {
0112     //    bool status = theGrid->getQuads(d);
0113     theGrid->setSpotEnergy(espot);
0114     result = theGrid->addHit(rshower * sin(theta), phi, 0);
0115   }
0116   return (result);
0117 }
0118
0119 float HDRShower::getR() {
0120   float p = random->flatShoot();
0121   unsigned int i = 1;
0122   while (rpdf[i] < p && i < R_range - 1) {
0123     i++;
0124   }
0125   float r;
0126   float dr = rpdf[i] - rpdf[i - 1];
0127   if (dr != 0.0)
0128     r = (float(i) + (p - rpdf[i - 1]) / dr) / lambdaHD;
0129   else
0130     r = float(i) / lambdaHD;
0131   return (r);
0132 }
0133
0134 void HDRShower::thetaFunction(int nthetaStep) {
0135   unsigned int i = 0;
0136   while (EgridTable[i] < e && i < NEnergyScan - 1) {
0137     i++;
0138   }
0139
0140   float amean, asig, lambda1, lambda1sig, lam21, lam21sig;
0141   amean = Theta1amp[i];
0142   asig = Theta1ampSig[i];
0143   lambda1 = Theta1Lambda[i];
0144   lambda1sig = Theta1LambdaSig[i];
0145   lam21 = ThetaLam21[i];
0146   lam21sig = ThetaLam21Sig[i];
0147   if (i == 0)
0148     i = 1;  //extrapolation to the left
0149   float c = (e - EgridTable[i - 1]) / (EgridTable[i] - EgridTable[i - 1]);
0150
0151   amean += (Theta1amp[i] - Theta1amp[i - 1]) * c;
0152   asig += (Theta1ampSig[i] - Theta1ampSig[i - 1]) * c;
0153   lambda1 += (Theta1Lambda[i] - Theta1Lambda[i - 1]) * c;
0154   lambda1sig += (Theta1LambdaSig[i] - Theta1LambdaSig[i - 1]) * c;
0155   lam21 += (ThetaLam21[i] - ThetaLam21[i - 1]) * c;
0156   lam21sig += (ThetaLam21Sig[i] - ThetaLam21Sig[i - 1]) * c;
0157
0158   float a = exp(amean + asig * random->gaussShoot());
0159   float L1 = lambda1 + lambda1sig * random->gaussShoot();
0160   if (L1 < 0.02)
0161     L1 = 0.02;
0162   float L2 = L1 * (lam21 + lam21sig * random->gaussShoot());
0163
0164   vector<double> pdf;
0165   pdf.erase(pdf.begin(), pdf.end());
0166   thetaSpots.erase(thetaSpots.begin(), thetaSpots.end());
0167   elastspot.erase(elastspot.begin(), elastspot.end());
0168   double sum = 0;
0169   for (int it = 0; it < nthetaStep; it++) {
0170     float theta = it * thetaStep;
0171     float p = a * exp(L1 * theta) + exp(L2 * theta);
0172     sum += p;
0173     pdf.push_back(p);
0174   }
0175   float ntot = e / eHDspot;
0176   float esum = 0;
0177   for (int it = 0; it < nthetaStep; it++) {
0178     float fn = ntot * pdf[it] / sum;
0179     thetaSpots.push_back(int(fn));
0180     elastspot.push_back((fn - int(fn)) * eHDspot);
0181   }
0182
0183   for (int it = 0; it < nthetaStep; it++)
0184     if (elastspot[it] < EsCut) {
0185       esum += elastspot[it];
0186       elastspot[it] = 0;
0187     }
0188
0189   float en = esum / EsCut;
0190   int n = int(en);
0191   en = esum - n * EsCut;
0192
0193   for (int ie = 0; ie <= n; ie++) {
0194     int k = int(nthetaStep * random->flatShoot());
0195     if (k < 0 || k > nthetaStep - 1)
0196       k = k % nthetaStep;
0197     if (ie == n)
0198       elastspot[k] += en;
0199     else
0200       elastspot[k] += EsCut;
0201   }
0202 }
0203
0204 void HDRShower::setFuncParam() {
0205   lambdaHD = theParam->hcalProperties()->interactionLength();
0206   x0HD = theParam->hcalProperties()->radLenIncm();
0207   if (onEcal)
0208     lambdaEM = theParam->ecalProperties()->interactionLength();
0209   else
0210     lambdaEM = lambdaHD;
0211
0212   if (debug)
0213     LogDebug("FastCalorimetry") << "setFuncParam-> lambdaEM=" << lambdaEM << " lambdaHD=" << lambdaHD << endl;
0214
0215   float _EgridTable[NEnergyScan] = {10, 20, 30, 50, 100, 300, 500};
0216   float _Theta1amp[NEnergyScan] = {1.57, 2.05, 2.27, 2.52, 2.66, 2.76, 2.76};
0217   float _Theta1ampSig[NEnergyScan] = {2.40, 1.50, 1.25, 1.0, 0.8, 0.52, 0.52};
0218
0219   float _Theta1Lambda[NEnergyScan] = {0.086, 0.092, 0.88, 0.80, 0.0713, 0.0536, 0.0536};
0220   float _Theta1LambdaSig[NEnergyScan] = {0.038, 0.037, 0.027, 0.03, 0.023, 0.018, 0.018};
0221
0222   float _ThetaLam21[NEnergyScan] = {2.8, 2.44, 2.6, 2.77, 3.16, 3.56, 3.56};
0223   float _ThetaLam21Sig[NEnergyScan] = {1.8, 0.97, 0.87, 0.77, 0.7, 0.49, 0.49};
0224
0225   for (int i = 0; i < NEnergyScan; i++) {
0226     EgridTable[i] = _EgridTable[i];
0227     Theta1amp[i] = _Theta1amp[i];
0228     Theta1ampSig[i] = _Theta1ampSig[i];
0229     Theta1Lambda[i] = _Theta1Lambda[i];
0230     Theta1LambdaSig[i] = _Theta1LambdaSig[i];
0231     ThetaLam21[i] = _ThetaLam21[i];
0232     ThetaLam21Sig[i] = _ThetaLam21Sig[i];
0233   }
0234
0235 #define lambdafit 15.05
0236   float R_alfa = -0.0993 + 0.1114 * log(e);
0237   float R_p = 0.589191 + 0.0463392 * log(e);
0238   float R_beta_lam = (0.54134 - 0.00011148 * e) / 4.0 * lambdafit;  //was fitted in 4cmbin
0239   float LamOverX0 = lambdaHD / x0HD;                                // 10.52
0240   //  int R_range = 100; // 7 lambda
0241   //  rpdf.erase(rpdf.begin(),rpdf.end());
0242
0243   rpdf[0] = 0.;
0244   for (int i = 1; i < R_range; i++) {
0245     float x = (float(i)) / lambdaHD;
0246     float r = pow(x, R_alfa) * (R_p * exp(-R_beta_lam * x) + (1 - R_p) * exp(-LamOverX0 * R_beta_lam * x));
0247     rpdf[i] = r;
0248     //    rpdf.push_back(r);
0249   }
0250
0251   for (int i = 1; i < R_range; i++)
0252     rpdf[i] += rpdf[i - 1];
0253   for (int i = 0; i < R_range; i++)
0254     rpdf[i] /= rpdf[R_range - 1];
0255 }