HDShower.cc

Go to the documentation of this file.

//updated by Reza Goldouzian
//FastSimulation Headers
#include "FastSimulation/ShowerDevelopment/interface/HDShower.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"

#include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
#include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"

// CMSSW headers
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <cmath>

// number attempts for transverse distribution if exit on a spec. condition
#define infinity 10000
// debugging flag ( 0, 1, 2, 3)
#define debug 0

using namespace edm;

HDShower::HDShower(const RandomEngineAndDistribution* engine,
           HDShowerParametrization* myParam, 
           EcalHitMaker* myGrid,
           HcalHitMaker* myHcalHitMaker,
           int onECAL,
           double epart,
                   double pmip)
  : theParam(myParam), 
    theGrid(myGrid),
    theHcalHitMaker(myHcalHitMaker),
    onEcal(onECAL),
    e(epart),
//    pmip(pmip),
    random(engine)
{ 
  // To get an access to constants read in FASTCalorimeter
  //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();

  // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
  lossesOpt       = myParam->hsParameters()->getHDlossesOpt();
  nDepthSteps     = myParam->hsParameters()->getHDnDepthSteps();
  nTRsteps        = myParam->hsParameters()->getHDnTRsteps();
  transParam      = myParam->hsParameters()->getHDtransParam();
  eSpotSize       = myParam->hsParameters()->getHDeSpotSize();
  depthStep       = myParam->hsParameters()->getHDdepthStep();
  criticalEnergy  = myParam->hsParameters()->getHDcriticalEnergy();
  maxTRfactor     = myParam->hsParameters()->getHDmaxTRfactor();
  balanceEH       = myParam->hsParameters()->getHDbalanceEH();
  hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();

  // Special tr.size fluctuations 
  transParam *= (1. + random->flatShoot()); 

  // Special long. fluctuations
  hcalDepthFactor +=  0.05 * (2.* random->flatShoot() - 1.);

  transFactor = 1.;   // normally 1, in HF - might be smaller 
                      // to take into account
                      // a narrowness of the HF shower (Cherenkov light) 

  // simple protection ...
  if(e < 0) e = 0.;

  // Get the Famos Histos pointer
  //  myHistos = FamosHistos::instance();
  //  std::cout << " Hello FamosShower " << std::endl;
  
  theECALproperties = theParam->ecalProperties();
  theHCALproperties = theParam->hcalProperties();


  double emax = theParam->emax(); 
  double emid = theParam->emid(); 
  double emin = theParam->emin(); 
  double effective = e;

  if( e < emid ) {
    theParam->setCase(1);
    // avoid "underflow" below Emin (for parameters calculation only)
    if(e < emin) effective = emin; 
  } 
  else 
    theParam->setCase(2);
 
  // Avoid "overflow" beyond Emax (for parameters)
  if(effective > 0.5 * emax) {
    eSpotSize *= 2.5;
    if(effective > emax) {
      effective = emax; 
      eSpotSize *= 4.; 
      depthStep *= 2.;
      if(effective > 1000.)
      eSpotSize *= 2.; 
    }
  }

  if(debug == 2 )
    LogInfo("FastCalorimetry") <<  " HDShower : " << std::endl 
         << "       Energy   "  <<               e << std::endl     
         << "      lossesOpt "  <<       lossesOpt << std::endl  
         << "    nDepthSteps "  <<     nDepthSteps << std::endl
         << "       nTRsteps "  <<        nTRsteps << std::endl
         << "     transParam "  <<      transParam << std::endl
         << "      eSpotSize "  <<       eSpotSize << std::endl
         << " criticalEnergy "  <<  criticalEnergy << std::endl
         << "    maxTRfactor "  <<     maxTRfactor << std::endl
         << "      balanceEH "  <<       balanceEH << std::endl
         << "hcalDepthFactor "  << hcalDepthFactor << std::endl;


  double alpEM1 = theParam->alpe1();
  double alpEM2 = theParam->alpe2();

  double betEM1 = theParam->bete1();
  double betEM2 = theParam->bete2();

  double alpHD1 = theParam->alph1();
  double alpHD2 = theParam->alph2();

  double betHD1 = theParam->beth1();
  double betHD2 = theParam->beth2();

  double part1 = theParam->part1();
  double part2 = theParam->part2();

  aloge = std::log(effective);
 
  double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
  double aedep = std::log(edpar);

  if(debug == 2)
    LogInfo("FastCalorimetry") << " HDShower : " << std::endl
         << "     edpar " <<   edpar << "   aedep " << aedep << std::endl 
         << "    alpEM1 " <<  alpEM1 << std::endl  
         << "    alpEM2 " <<  alpEM2 << std::endl  
         << "    betEM1 " <<  betEM1 << std::endl  
         << "    betEM2 " <<  betEM2 << std::endl  
         << "    alpHD1 " <<  alpHD1 << std::endl  
         << "    alpHD2 " <<  alpHD2 << std::endl  
         << "    betHD1 " <<  betHD1 << std::endl  
         << "    betHD2 " <<  betHD2 << std::endl  
         << "     part1 " <<   part1 << std::endl  
         << "     part2 " <<   part2 << std::endl; 

  // private members to set
  theR1  = theParam->r1();
  theR2  = theParam->r2();
  theR3  = theParam->r3();

  alpEM  = alpEM1 + alpEM2 * aedep;
  tgamEM = tgamma(alpEM);
  betEM  = betEM1 - betEM2 * aedep;
  alpHD  = alpHD1 + alpHD2 * aedep;
  tgamHD = tgamma(alpHD);
  betHD  = betHD1 - betHD2 * aedep;
  part   = part1  -  part2 * aedep;
  if(part > 1.) part = 1.;          // protection - just in case of 

  if(debug  == 2 )
    LogInfo("FastCalorimetry") << " HDShower : " << std::endl 
         << "    alpEM " <<  alpEM << std::endl  
         << "   tgamEM " << tgamEM << std::endl
         << "    betEM " <<  betEM << std::endl
         << "    alpHD " <<  alpHD << std::endl
         << "   tgamHD " << tgamHD << std::endl
         << "    betHD " <<  betHD << std::endl
         << "     part " <<   part << std::endl;
       

  if(onECAL){
    lambdaEM = theParam->ecalProperties()->interactionLength();
    x0EM     = theParam->ecalProperties()->radLenIncm();
  } 
  else {
    lambdaEM = 0.;
    x0EM     = 0.;
  }
  lambdaHD = theParam->hcalProperties()->interactionLength();
  x0HD     = theParam->hcalProperties()->radLenIncm();

  if(debug == 2)
    LogInfo("FastCalorimetry") << " HDShower e " << e        << std::endl
         << "          x0EM = " << x0EM     << std::endl 
         << "          x0HD = " << x0HD     << std::endl 
         << "         lamEM = " << lambdaEM << std::endl
         << "         lamHD = " << lambdaHD << std::endl;


  // Starting point of the shower
  // try first with ECAL lambda  

  double sum1   = 0.;  // lambda path from the ECAL/HF entrance;  
  double sum2   = 0.;  // lambda path from the interaction point;
  double sum3   = 0.;  // x0     path from the interaction point;
  int  nsteps   = 0;   // full number of longitudinal steps (counter);

  int nmoresteps;      // how many longitudinal steps in addition to 
                       // one (if interaction happens there) in ECAL

  mip = 1;             // just to initiate particle as MIP in ECAL    

  if(e < criticalEnergy ) nmoresteps = 1;  
  else                    nmoresteps = nDepthSteps;
  
  depthECAL  = 0.;
  depthGAP   = 0.;
  depthGAPx0 = 0.; 
  if(onECAL ) {
    depthECAL  = theGrid->ecalTotalL0();         // ECAL depth segment
    //depthECAL  = theGrid->ecalTotalL0() + theGrid->ps1TotalL0() + theGrid->ps2TotalL0() + theGrid->ps2eeTotalL0(); //TEST: include preshower
    depthGAP   = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
    depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
  }
  
  depthHCAL   = theGrid->hcalTotalL0();    // HCAL depth segment
  depthToHCAL = depthECAL + depthGAP;    

  //---------------------------------------------------------------------------
  // Depth simulation & various protections, among them
  // if too deep - get flat random in the allowed region
  // if no HCAL material behind - force to deposit in ECAL
  double maxDepth    = depthToHCAL + depthHCAL - 1.1 * depthStep;
  double depthStart  = std::log(1./random->flatShoot()); // starting point lambda unts
  //std::cout<<"generated depth             "<<depthStart<<std::endl;
  if (pmip==1) {depthStart=depthToHCAL;}
  else {depthStart=depthStart*0.9*depthECAL/std::log(1./pmip);}
 // std::cout<<"modified  depth             "<<depthStart<<std::endl;


  if(e < emin) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl; 
    depthStart = 0.;
  }
 
  if(depthStart > maxDepth) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " << depthStart << std::endl;     
    depthStart = maxDepth *  random->flatShoot();
    if(depthStart < 0.) depthStart = 0.;
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart re-calculated = " << depthStart << std::endl; 
  }
  
  if(onECAL && e < emid) {
    if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.2) {
      depthStart = 0.5 * depthECAL * random->flatShoot();
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : small energy, " << " depthStart reduced to = " << depthStart << std::endl; 
    }
  }

  if(depthHCAL < depthStep) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " << depthHCAL << " depthStart -> forced to 0 !!!" << std::endl;
    depthStart = 0.;    
    nmoresteps = 0;
    
    if(depthECAL < depthStep) {
      nsteps = -1;
      LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - " << " particle is lost !!! " << std::endl; 
    }
  }

  if(debug)
    LogInfo("FastCalorimetry") << " FamosHDShower  depths(lam) - "  << std::endl 
         << "          ECAL = " << depthECAL  << std::endl
         << "           GAP = " << depthGAP   << std::endl
         << "          HCAL = " << depthHCAL  << std::endl
         << "  starting point = " << depthStart << std::endl; 

  if( onEcal ) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
    if(depthStart < depthECAL) {
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
      if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.1) {
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step" << std::endl;
        
        //  ECAL - one step
        nsteps++; 
        sum1   += depthECAL;             // at the end of step
        sum2   += depthECAL-depthStart; 
        sum3   += sum2 * lambdaEM / x0EM;
        lamtotal.push_back(sum1);
        lamdepth.push_back(sum2);
        lamcurr.push_back(lambdaEM);
        lamstep.push_back(depthECAL-depthStart);
        x0depth.push_back(sum3);
        x0curr.push_back(x0EM);
        detector.push_back(1);
        mip = 0;    
        
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : " << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;
        
        // Gap - no additional step after ECAL
        // just move further to HCAL over the gap
        sum1   += depthGAP;          
        sum2   += depthGAP; 
        sum3   += depthGAPx0;
      }
      else { // Just shift starting point to HCAL
        //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
        if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;

        depthStart = depthToHCAL;
        sum1 += depthStart;     
      }
    }
    else { // GAP or HCAL         
      if(depthStart >= depthECAL &&  depthStart < depthToHCAL ) {
        depthStart = depthToHCAL; // just a shift to HCAL for simplicity
      }
      sum1 += depthStart;
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
    }
  }
  else {   // Forward 
    if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
    sum1 += depthStart;
    transFactor = 0.5;   // makes narower tresverse size of shower     
  }
 
  for (int i = 0; i < nmoresteps ; i++) {
    sum1 += depthStep;
    if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; 
    sum2 += depthStep;
    sum3 += sum2 * lambdaHD / x0HD;
    lamtotal.push_back(sum1);
    lamdepth.push_back(sum2);
    lamcurr.push_back(lambdaHD);
    lamstep.push_back(depthStep);
    x0depth.push_back(sum3);
    x0curr.push_back(x0HD);
    detector.push_back(3);
    nsteps++;
  }

  // Make fractions of energy and transverse radii at each step 

  if(nsteps > 0) { makeSteps(nsteps); }

}

void HDShower::makeSteps(int nsteps) {

  double sumes = 0.;
  double sum   = 0.;
  std::vector<double> temp;

  if(debug)
    LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - " 
       << " nsteps required : " << nsteps << std::endl;

  int count = 0;
  for (int i = 0; i < nsteps; i++) {    
    
    double deplam = lamdepth[i] - 0.5 * lamstep[i];
    double depx0  = x0depth[i]  - 0.5 * lamstep[i] / x0curr[i]; 
    double     x = betEM * depx0;
    double     y = betHD * deplam;
    
    if(debug == 2)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps " 
                                  << " - step " << i
                  << "   depx0, x = " << depx0 << ", " << x 
                  << "   deplam, y = " << deplam << ", "
                  << y << std::endl;
    
    double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
          (x0curr[i] * tgamEM) + 
          (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
    
    // protection ...
    if(est < 0.) {
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " - negative step energy !!!" 
       << std::endl;
      est = 0.;
      break ; 
    }

    // for estimates only
    sum += est;
    int nPest = (int) (est * e / sum / eSpotSize) ;

    if(debug == 2)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " 
       <<  nPest << std::endl;

    if(nPest <= 1 && count !=0 ) break;

    // good step - to proceed

    temp.push_back(est);
    sumes += est;
    
    rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge))
               * deplam * transFactor); 
    count ++;
  }

  // fluctuations in ECAL and re-distribution of remaining energy in HCAL
  if(detector[0] == 1 && count > 1) {
    double oldECALenergy = temp[0];
    double oldHCALenergy = sumes - oldECALenergy ;
    double newECALenergy = 2. * sumes;
    for (int i = 0; newECALenergy > sumes && i < infinity; i++)
      newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy; 
     
    if(debug == 2)
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " ECAL fraction : old/new - "
       << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl;

    temp[0] = newECALenergy;
    double newHCALenergy = sumes - newECALenergy;
    double newHCALreweight =  newHCALenergy / oldHCALenergy;

    for (int i = 1; i < count; i++) {
      temp[i] *= newHCALreweight;
    }
    
  }
  
  // final re-normalization of the energy fractions  
  for (int i = 0; i < count ; i++) {
    eStep.push_back(temp[i] * e / sumes );
    nspots.push_back((int)(eStep[i]/eSpotSize)+1);

   if(debug)
     LogInfo("FastCalorimetry") 
       << " step " << i
       << "  det: " << detector[i]   
       << "  xO and lamdepth at the end of step = " 
       << x0depth[i] << " "  
       << lamdepth[i] << "   Estep func = " <<  eStep[i]
       << "   Rstep = " << rlamStep[i] << "  Nspots = " <<  nspots[i]
       << "  espot = " <<  eStep[i] / (double)nspots[i]
       << std::endl; 

  }

  // The only step is in ECAL - let's make the size bigger ...  
  if(count == 1 and detector[0] == 1) rlamStep[0] *= 2.;

  if(debug) {
    if(eStep[0] > 0.95 * e && detector[0] == 1) 
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - " << "ECAL energy = " << eStep[0]
       << " out of total = " << e << std::endl;  
  }

}

bool HDShower::compute() {
  
  //  TimeMe theT("FamosHDShower::compute");

  bool status = false;
  int numLongit = eStep.size();
  if(debug)
    LogInfo("FastCalorimetry") << " FamosHDShower::compute - " 
        << " N_long.steps required : " << numLongit << std::endl;

  if(numLongit > 0) {

    status = true;    
    // Prepare the trsanverse probability function
    std::vector<double> Fhist;
    std::vector<double> rhist; 
    for (int j = 0; j < nTRsteps + 1; j++) {
      rhist.push_back(maxTRfactor * j / nTRsteps );  
      Fhist.push_back(transProb(maxTRfactor,1.,rhist[j]));
      if(debug == 3) LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
    }
    
    // Longitudinal steps
    for (int i = 0; i < numLongit ; i++) {
      
      double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
      // vary the longitudinal profile if needed
      if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor;                     
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i] << "  currentDepthL0 = " << currentDepthL0 << std::endl;
      
      double maxTRsize   = maxTRfactor * rlamStep[i];     // in lambda units
      double rbinsize    = maxTRsize / nTRsteps; 
      double espot       = eStep[i] / (double)nspots[i];  // re-adjust espot

      if(espot > 2. || espot < 0. ) LogInfo("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " << espot << std::endl;

      int ecal = 0;
      if(detector[i] != 1) { 
        bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
    
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of " << " theHcalHitMaker->setDepth(currentDepthL0) is " 
                    << setHDdepth << std::endl;
    
        if(!setHDdepth) {
          currentDepthL0 -= lamstep[i];
          setHDdepth =  theHcalHitMaker->setDepth(currentDepthL0);
        }
        
        if(!setHDdepth) continue;

        theHcalHitMaker->setSpotEnergy(espot);
        
        //fill hcal longitudinal distribution histogram
      }
      else {
        ecal = 1;
        bool status = theGrid->getPads(currentDepthL0);   
    
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;
    
        if(!status) continue; 

        int ntry = nspots[i] * 10;  
        if( ntry >= infinity ) {  // use max allowed in case of too many spots
          nspots[i] = 0.5 * infinity;  
          espot *= 0.1 * (double)ntry / double(nspots[i]);
        }     
        else {
          espot *= 0.1;  // fine-grain energy spots in ECAL
                         // to avoid false ECAL clustering 
          nspots[i] = ntry;
        }

        theGrid->setSpotEnergy(espot);
        
        //fill ecal longitudinal distribution histogram
      }  
      
      // Transverse distribition
      int nok   = 0;                          // counter of OK  
      int count = 0;
      int inf   = infinity;
      if(lossesOpt) inf = nspots[i];          // if losses are enabled, otherwise
      // only OK points are counted ...
      if(nspots[i] > inf ) std::cout << " FamosHDShower::compute - at long.step " << i << "  too many spots required : "  <<  nspots[i] << " !!! " << std::endl;

      for (int j = 0; j < inf; j++) {
        if(nok == nspots[i]) break;
        count ++;
        
        double prob   = random->flatShoot();
        int index     = indexFinder(prob,Fhist);
        double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize; // in-bin  
        double phi = 2.*M_PI*random->flatShoot();
    
        if(debug == 2)  LogInfo("FastCalorimetry") << std::endl << " FamosHDShower::compute " << " r = " << radius 
                        << "    phi = " << phi << std::endl;
    
        bool result;
        if(ecal) {
          result = theGrid->addHit(radius,phi,0);
      
          if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theGrid->addHit result = " << result << std::endl;
          
          //fill ecal transverse distribution histogram
        }
        else {
          result = theHcalHitMaker->addHit(radius,phi,0); 
 
          if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theHcalHitMaker->addHit result = " << result << std::endl;
          
          //fill hcal transverse distribution histogram
        }    
        
        if(result) nok ++; 



      } // end of tranverse simulation
      
      
      if(count == infinity) { 
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute " << " maximum number of" 
                    << " transverse points " << count << " is used !!!" << std::endl; 
      }

      if(debug) LogInfo("FastCalorimetry")  << " FamosHDShower::compute " << " long.step No." 
                << i << "   Ntry, Nok = " << count << " " << nok << std::endl; 
    } // end of longitudinal steps
  } // end of no steps

  return status;

}

int HDShower::indexFinder(double x, const std::vector<double> & Fhist) {
  // binary search in the vector of doubles
  int size = Fhist.size();

  int curr = size / 2;
  int step = size / 4;
  int iter;
  int prevdir = 0; 
  int actudir = 0; 

  for (iter = 0; iter < size ; iter++) {

    if( curr >= size || curr < 1 )
      LogWarning("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " 
       << curr << " !!!" << std::endl;

    if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break;
    prevdir = actudir;
    if(x > Fhist[curr]) {actudir =  1;}
    else                {actudir = -1;}
    if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
    curr += actudir * step;
    if(curr > size) curr = size;
    else { if(curr < 1) {curr = 1;}}
 
    if(debug == 3)
      LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter 
       << " curr, F[curr-1], F[curr] = "
       << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;
    
  }

  if(debug == 3)
    LogInfo("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr-1
         << std::endl;


  return curr-1;
}