Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:26:45

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   float esum = e;
0075   for (int itheta = 0; itheta < nthetaStep; itheta++) {
0076     float theta, es;
0077     for (int i = 0; i <= thetaSpots[itheta]; i++) {
0078       if (i == thetaSpots[itheta])
0079         es = elastspot[itheta];
0080       else
0081         es = eHDspot;
0082       float loops = 0;
0083       for (int j = 0; j < maxLoops; j++) {
0084         theta = (itheta + random->flatShoot()) * thetaStep;
0085         if (setHit(es, theta))
0086           break;
0087         loops++;
0088       }
0089       esum -= es;  // to check only
0090     }
0091   }
0092   return (true);
0093 }
0094 
0095 bool HDRShower::setHit(float espot, float theta) {
0096   float phi = 2. * M_PI * random->flatShoot();  // temporary: 1st approximation
0097   float rshower = getR();                       // temporary: 1st approximation
0098 
0099   float d = depthStart + rshower * cos(theta);
0100   if (d + depthGAP > maxDepth)
0101     return (false);
0102 
0103   // Commented (F.B) to remove a warning. Not used anywhere ?
0104   //  bool inHcal = !onEcal || d>depthECAL || !qstatus;
0105   bool result = false;
0106   if (!onEcal || d > depthECAL || !qstatus) {  // in HCAL (HF or HB, HE)
0107     d += depthGAP;
0108     bool setHDdepth = theHcalHitMaker->setDepth(d);
0109     if (setHDdepth) {
0110       theHcalHitMaker->setSpotEnergy(espot);
0111       result = theHcalHitMaker->addHit(rshower * sin(theta), phi, 0);
0112     } else
0113       LogWarning("FastCalorimetry") << " setHit in HCAL failed d=" << d << " maxDepth=" << maxDepth << " onEcal'"
0114                                     << onEcal << endl;
0115   } else {
0116     //    bool status = theGrid->getQuads(d);
0117     theGrid->setSpotEnergy(espot);
0118     result = theGrid->addHit(rshower * sin(theta), phi, 0);
0119   }
0120   return (result);
0121 }
0122 
0123 float HDRShower::getR() {
0124   float p = random->flatShoot();
0125   unsigned int i = 1;
0126   while (rpdf[i] < p && i < R_range - 1) {
0127     i++;
0128   }
0129   float r;
0130   float dr = rpdf[i] - rpdf[i - 1];
0131   if (dr != 0.0)
0132     r = (float(i) + (p - rpdf[i - 1]) / dr) / lambdaHD;
0133   else
0134     r = float(i) / lambdaHD;
0135   return (r);
0136 }
0137 
0138 void HDRShower::thetaFunction(int nthetaStep) {
0139   unsigned int i = 0;
0140   while (EgridTable[i] < e && i < NEnergyScan - 1) {
0141     i++;
0142   }
0143 
0144   float amean, asig, lambda1, lambda1sig, lam21, lam21sig;
0145   amean = Theta1amp[i];
0146   asig = Theta1ampSig[i];
0147   lambda1 = Theta1Lambda[i];
0148   lambda1sig = Theta1LambdaSig[i];
0149   lam21 = ThetaLam21[i];
0150   lam21sig = ThetaLam21Sig[i];
0151   if (i == 0)
0152     i = 1;  //extrapolation to the left
0153   float c = (e - EgridTable[i - 1]) / (EgridTable[i] - EgridTable[i - 1]);
0154 
0155   amean += (Theta1amp[i] - Theta1amp[i - 1]) * c;
0156   asig += (Theta1ampSig[i] - Theta1ampSig[i - 1]) * c;
0157   lambda1 += (Theta1Lambda[i] - Theta1Lambda[i - 1]) * c;
0158   lambda1sig += (Theta1LambdaSig[i] - Theta1LambdaSig[i - 1]) * c;
0159   lam21 += (ThetaLam21[i] - ThetaLam21[i - 1]) * c;
0160   lam21sig += (ThetaLam21Sig[i] - ThetaLam21Sig[i - 1]) * c;
0161 
0162   float a = exp(amean + asig * random->gaussShoot());
0163   float L1 = lambda1 + lambda1sig * random->gaussShoot();
0164   if (L1 < 0.02)
0165     L1 = 0.02;
0166   float L2 = L1 * (lam21 + lam21sig * random->gaussShoot());
0167 
0168   vector<double> pdf;
0169   pdf.erase(pdf.begin(), pdf.end());
0170   thetaSpots.erase(thetaSpots.begin(), thetaSpots.end());
0171   elastspot.erase(elastspot.begin(), elastspot.end());
0172   double sum = 0;
0173   for (int it = 0; it < nthetaStep; it++) {
0174     float theta = it * thetaStep;
0175     float p = a * exp(L1 * theta) + exp(L2 * theta);
0176     sum += p;
0177     pdf.push_back(p);
0178   }
0179   float ntot = e / eHDspot;
0180   float esum = 0;
0181   for (int it = 0; it < nthetaStep; it++) {
0182     float fn = ntot * pdf[it] / sum;
0183     thetaSpots.push_back(int(fn));
0184     elastspot.push_back((fn - int(fn)) * eHDspot);
0185   }
0186 
0187   for (int it = 0; it < nthetaStep; it++)
0188     if (elastspot[it] < EsCut) {
0189       esum += elastspot[it];
0190       elastspot[it] = 0;
0191     }
0192 
0193   float en = esum / EsCut;
0194   int n = int(en);
0195   en = esum - n * EsCut;
0196 
0197   for (int ie = 0; ie <= n; ie++) {
0198     int k = int(nthetaStep * random->flatShoot());
0199     if (k < 0 || k > nthetaStep - 1)
0200       k = k % nthetaStep;
0201     if (ie == n)
0202       elastspot[k] += en;
0203     else
0204       elastspot[k] += EsCut;
0205   }
0206 }
0207 
0208 void HDRShower::setFuncParam() {
0209   lambdaHD = theParam->hcalProperties()->interactionLength();
0210   x0HD = theParam->hcalProperties()->radLenIncm();
0211   if (onEcal)
0212     lambdaEM = theParam->ecalProperties()->interactionLength();
0213   else
0214     lambdaEM = lambdaHD;
0215 
0216   if (debug)
0217     LogDebug("FastCalorimetry") << "setFuncParam-> lambdaEM=" << lambdaEM << " lambdaHD=" << lambdaHD << endl;
0218 
0219   float _EgridTable[NEnergyScan] = {10, 20, 30, 50, 100, 300, 500};
0220   float _Theta1amp[NEnergyScan] = {1.57, 2.05, 2.27, 2.52, 2.66, 2.76, 2.76};
0221   float _Theta1ampSig[NEnergyScan] = {2.40, 1.50, 1.25, 1.0, 0.8, 0.52, 0.52};
0222 
0223   float _Theta1Lambda[NEnergyScan] = {0.086, 0.092, 0.88, 0.80, 0.0713, 0.0536, 0.0536};
0224   float _Theta1LambdaSig[NEnergyScan] = {0.038, 0.037, 0.027, 0.03, 0.023, 0.018, 0.018};
0225 
0226   float _ThetaLam21[NEnergyScan] = {2.8, 2.44, 2.6, 2.77, 3.16, 3.56, 3.56};
0227   float _ThetaLam21Sig[NEnergyScan] = {1.8, 0.97, 0.87, 0.77, 0.7, 0.49, 0.49};
0228 
0229   for (int i = 0; i < NEnergyScan; i++) {
0230     EgridTable[i] = _EgridTable[i];
0231     Theta1amp[i] = _Theta1amp[i];
0232     Theta1ampSig[i] = _Theta1ampSig[i];
0233     Theta1Lambda[i] = _Theta1Lambda[i];
0234     Theta1LambdaSig[i] = _Theta1LambdaSig[i];
0235     ThetaLam21[i] = _ThetaLam21[i];
0236     ThetaLam21Sig[i] = _ThetaLam21Sig[i];
0237   }
0238 
0239 #define lambdafit 15.05
0240   float R_alfa = -0.0993 + 0.1114 * log(e);
0241   float R_p = 0.589191 + 0.0463392 * log(e);
0242   float R_beta_lam = (0.54134 - 0.00011148 * e) / 4.0 * lambdafit;  //was fitted in 4cmbin
0243   float LamOverX0 = lambdaHD / x0HD;                                // 10.52
0244   //  int R_range = 100; // 7 lambda
0245   //  rpdf.erase(rpdf.begin(),rpdf.end());
0246 
0247   rpdf[0] = 0.;
0248   for (int i = 1; i < R_range; i++) {
0249     float x = (float(i)) / lambdaHD;
0250     float r = pow(x, R_alfa) * (R_p * exp(-R_beta_lam * x) + (1 - R_p) * exp(-LamOverX0 * R_beta_lam * x));
0251     rpdf[i] = r;
0252     //    rpdf.push_back(r);
0253   }
0254 
0255   for (int i = 1; i < R_range; i++)
0256     rpdf[i] += rpdf[i - 1];
0257   for (int i = 0; i < R_range; i++)
0258     rpdf[i] /= rpdf[R_range - 1];
0259 }