GflashEMShowerProfile.cc

Go to the documentation of this file.

//
// initial setup : Soon Jun & Dongwook Jang
// Translated from Fortran code.
//
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "SimGeneral/GFlash/interface/Gflash3Vector.h"
#include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
#include "SimGeneral/GFlash/interface/GflashHistogram.h"
#include "SimGeneral/GFlash/interface/GflashHit.h"
#include "SimGeneral/GFlash/interface/GflashShowino.h"
#include "SimGeneral/GFlash/interface/GflashTrajectory.h"
#include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"

#include <CLHEP/GenericFunctions/IncompleteGamma.hh>
#include <CLHEP/Random/RandGaussQ.h>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Units/PhysicalConstants.h>

GflashEMShowerProfile::GflashEMShowerProfile(const edm::ParameterSet &parSet) : theParSet(parSet) {
  theBField = parSet.getParameter<double>("bField");
  theEnergyScaleEB = parSet.getParameter<double>("energyScaleEB");
  theEnergyScaleEE = parSet.getParameter<double>("energyScaleEE");

  jCalorimeter = Gflash::kNULL;

  theShowino = new GflashShowino();

  theHisto = GflashHistogram::instance();
}

GflashEMShowerProfile::~GflashEMShowerProfile() {
  if (theShowino)
    delete theShowino;
}

void GflashEMShowerProfile::initialize(
    int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum) {
  theShowino->initialize(showerType, energy, globalTime, charge, position, momentum, theBField);
}

void GflashEMShowerProfile::parameterization() {
  // This part of code is copied from the original GFlash Fortran code.
  // reference : hep-ex/0001020v1
  // The units used in Geant4 internally are in mm, MeV.
  // For simplicity, the units here are in cm, GeV.

  const double energyCutoff = 0.01;
  const int maxNumberOfSpots = 100000;

  double incomingEnergy = theShowino->getEnergy();
  Gflash3Vector showerStartingPosition = theShowino->getPositionAtShower();

  // find the calorimeter at the shower starting point
  jCalorimeter = Gflash::getCalorimeterNumber(showerStartingPosition);

  double logEinc = std::log(incomingEnergy);
  double y = incomingEnergy / Gflash::criticalEnergy;  // y = E/Ec, criticalEnergy is in GeV
  double logY = std::log(y);

  // Total number of spots are not yet optimized.
  double nSpots = 93.0 * std::log(Gflash::Z[jCalorimeter]) * std::pow(incomingEnergy, 0.876);

  // path Length from the origin to the shower starting point in cm
  double pathLength0 = theShowino->getPathLengthAtShower();
  double pathLength = pathLength0;  // this will grow along the shower development

  //--- 2.2  Fix intrinsic properties of em. showers.

  double fluctuatedTmax = std::log(logY - 0.7157);
  double fluctuatedAlpha = std::log(0.7996 + (0.4581 + 1.8628 / Gflash::Z[jCalorimeter]) * logY);

  double sigmaTmax = 1.0 / (-1.4 + 1.26 * logY);
  double sigmaAlpha = 1.0 / (-0.58 + 0.86 * logY);
  double rho = 0.705 - 0.023 * logY;
  double sqrtPL = std::sqrt((1.0 + rho) / 2.0);
  double sqrtLE = std::sqrt((1.0 - rho) / 2.0);

  double norm1 = CLHEP::RandGaussQ::shoot();
  double norm2 = CLHEP::RandGaussQ::shoot();
  double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
  double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);

  // tmax, alpha, beta : parameters of gamma distribution
  double tmax = std::exp(tempTmax);
  double alpha = std::exp(tempAlpha);
  double beta = std::max(0.0, (alpha - 1.0) / tmax);

  // spot fluctuations are added to tmax, alpha, beta
  double averageTmax = logY - 0.858;
  double averageAlpha = 0.21 + (0.492 + 2.38 / Gflash::Z[jCalorimeter]) * logY;
  double spotTmax = averageTmax * (0.698 + .00212 * Gflash::Z[jCalorimeter]);
  double spotAlpha = averageAlpha * (0.639 + .00334 * Gflash::Z[jCalorimeter]);
  double spotBeta = std::max(0.0, (spotAlpha - 1.0) / spotTmax);

  if (theHisto->getStoreFlag()) {
    theHisto->em_incE->Fill(incomingEnergy);
    theHisto->em_ssp_rho->Fill(showerStartingPosition.rho());
    theHisto->em_ssp_z->Fill(showerStartingPosition.z());
  }

  //  parameters for lateral distribution and fluctuation
  double z1 = 0.0251 + 0.00319 * logEinc;
  double z2 = 0.1162 - 0.000381 * Gflash::Z[jCalorimeter];

  double k1 = 0.659 - 0.00309 * Gflash::Z[jCalorimeter];
  double k2 = 0.645;
  double k3 = -2.59;
  double k4 = 0.3585 + 0.0421 * logEinc;

  double p1 = 2.623 - 0.00094 * Gflash::Z[jCalorimeter];
  double p2 = 0.401 + 0.00187 * Gflash::Z[jCalorimeter];
  double p3 = 1.313 - 0.0686 * logEinc;

  // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : e25Scale = 1.006
  // for EB with ecalNotContainment = 1.0. Now e25Scale is a configurable
  // parameter with default ecalNotContainment which is 0.97 for EB and 0.975
  // for EE. So if ecalNotContainment constants are to be changed in the future,
  // e25Scale should be changed accordingly.
  double e25Scale = 1.0;
  if (jCalorimeter == Gflash::kESPM)
    e25Scale = theEnergyScaleEB;
  else if (jCalorimeter == Gflash::kENCA)
    e25Scale = theEnergyScaleEE;

  // @@@ dwjang, intial tuning by comparing 20-150GeV TB data : p1 *= 0.965
  p1 *= 0.965;

  // preparation of longitudinal integration
  double stepLengthLeft = getDistanceToOut(jCalorimeter);

  int nSpots_sd = 0;                // count total number of spots in SD
  double zInX0 = 0.0;               // shower depth in X0 unit
  double deltaZInX0 = 0.0;          // segment of depth in X0 unit
  double deltaZ = 0.0;              // segment of depth in cm
  double stepLengthLeftInX0 = 0.0;  // step length left in X0 unit

  const double divisionStepInX0 = 0.1;  // step size in X0 unit
  double energy = incomingEnergy;       // energy left in GeV

  Genfun::IncompleteGamma gammaDist;

  double energyInGamma = 0.0;     // energy in a specific depth(z) according to Gamma distribution
  double preEnergyInGamma = 0.0;  // energy calculated in a previous depth
  double sigmaInGamma = 0.;       // sigma of energy in a specific depth(z) according
                                  // to Gamma distribution
  double preSigmaInGamma = 0.0;   // sigma of energy in a previous depth

  // energy segment in Gamma distribution of shower in each step
  double deltaEnergy = 0.0;  // energy in deltaZ
  int spotCounter = 0;       // keep track of number of spots generated

  // step increment along the shower direction
  double deltaStep = 0.0;

  theGflashHitList.clear();

  // loop for longitudinal integration
  while (energy > 0.0 && stepLengthLeft > 0.0) {
    stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];

    if (stepLengthLeftInX0 < divisionStepInX0) {
      deltaZInX0 = stepLengthLeftInX0;
      deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
      stepLengthLeft = 0.0;
    } else {
      deltaZInX0 = divisionStepInX0;
      deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
      stepLengthLeft -= deltaZ;
    }

    zInX0 += deltaZInX0;

    int nSpotsInStep = 0;

    if (energy > energyCutoff) {
      preEnergyInGamma = energyInGamma;
      gammaDist.a().setValue(alpha);            // alpha
      energyInGamma = gammaDist(beta * zInX0);  // beta
      double energyInDeltaZ = energyInGamma - preEnergyInGamma;
      deltaEnergy = std::min(energy, incomingEnergy * energyInDeltaZ);

      preSigmaInGamma = sigmaInGamma;
      gammaDist.a().setValue(spotAlpha);           // alpha spot
      sigmaInGamma = gammaDist(spotBeta * zInX0);  // beta spot
      nSpotsInStep = std::max(1, int(nSpots * (sigmaInGamma - preSigmaInGamma)));
    } else {
      deltaEnergy = energy;
      preSigmaInGamma = sigmaInGamma;
      nSpotsInStep = std::max(1, int(nSpots * (1.0 - preSigmaInGamma)));
    }

    if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
      deltaEnergy = energy;

    energy -= deltaEnergy;

    if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
      nSpotsInStep = maxNumberOfSpots - spotCounter;
      if (nSpotsInStep < 1) {  // @@ check
        edm::LogInfo("SimGeneralGFlash") << "GflashEMShowerProfile: Too Many Spots ";
        edm::LogInfo("SimGeneralGFlash") << " - break to regenerate nSpotsInStep ";
        break;
      }
    }

    // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
    deltaStep += 0.5 * deltaZ;
    pathLength += deltaStep;
    deltaStep = 0.5 * deltaZ;

    // lateral shape and fluctuations for  homogenous calo.
    double tScale = tmax * alpha / (alpha - 1.0) * (std::exp(fluctuatedAlpha) - 1.0) / std::exp(fluctuatedAlpha);
    double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
    double rCore = z1 + z2 * tau;
    double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));  // @@ check RT3 sign
    double p23 = (p2 - tau) / p3;
    double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));

    // Deposition of spots according to lateral distr.
    // Apply absolute energy scale
    // Convert into MeV unit
    double hitEnergy = deltaEnergy / nSpotsInStep * e25Scale * CLHEP::GeV;
    double hitTime = theShowino->getGlobalTime() * CLHEP::nanosecond + (pathLength - pathLength0) / 30.0;

    GflashHit aHit;

    for (int ispot = 0; ispot < nSpotsInStep; ispot++) {
      spotCounter++;

      // Compute global position of generated spots with taking into account
      // magnetic field Divide deltaZ into nSpotsInStep and give a spot a global
      // position
      double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);

      // trajectoryPoint give a spot an imaginary point along the shower
      // development
      GflashTrajectoryPoint trajectoryPoint;
      theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);

      double rShower = 0.0;
      Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint, rCore, rTail, probabilityWeight, rShower);

      // Convert into mm unit
      hitPosition *= CLHEP::cm;

      if (std::fabs(hitPosition.getZ() / CLHEP::cm) > Gflash::Zmax[Gflash::kHE])
        continue;

      // put energy and position to a Hit
      aHit.setTime(hitTime);
      aHit.setEnergy(hitEnergy);
      aHit.setPosition(hitPosition);
      theGflashHitList.push_back(aHit);

      double zInX0_spot = std::abs(pathLength + incrementPath - pathLength0) / Gflash::radLength[jCalorimeter];

      nSpots_sd++;

      // for histogramming
      if (theHisto->getStoreFlag()) {
        theHisto->em_long->Fill(zInX0_spot, hitEnergy / CLHEP::GeV);
        theHisto->em_lateral->Fill(zInX0_spot, rShower / Gflash::rMoliere[jCalorimeter], hitEnergy / CLHEP::GeV);
        theHisto->em_long_sd->Fill(zInX0_spot, hitEnergy / CLHEP::GeV);
        theHisto->em_lateral_sd->Fill(zInX0_spot, rShower / Gflash::rMoliere[jCalorimeter], hitEnergy / CLHEP::GeV);
      }

    }  // end of for spot iteration

  }  // end of while for longitudinal integration

  if (theHisto->getStoreFlag()) {
    theHisto->em_nSpots_sd->Fill(nSpots_sd);
  }

  //  delete theGflashNavigator;
}

double GflashEMShowerProfile::getDistanceToOut(Gflash::CalorimeterNumber kCalor) {
  double stepLengthLeft = 0.0;
  if (kCalor == Gflash::kESPM) {
    stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kESPM]) -
                     theShowino->getPathLengthAtShower();
  } else if (kCalor == Gflash::kENCA) {
    double zsign = (theShowino->getPosition()).getEta() > 0 ? 1.0 : -1.0;
    stepLengthLeft = theShowino->getHelix()->getPathLengthAtZ(zsign * Gflash::Zmax[Gflash::kENCA]) -
                     theShowino->getPathLengthAtShower();
  }
  return stepLengthLeft;
}

Gflash3Vector GflashEMShowerProfile::locateHitPosition(
    GflashTrajectoryPoint &point, double rCore, double rTail, double probability, double &rShower) {
  double u1 = CLHEP::HepUniformRand();
  double u2 = CLHEP::HepUniformRand();
  double rInRM = 0.0;

  if (u1 < probability) {
    rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
  } else {
    rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
  }

  rShower = rInRM * Gflash::rMoliere[jCalorimeter];

  // Uniform & random rotation of spot along the azimuthal angle
  double azimuthalAngle = CLHEP::twopi * CLHEP::HepUniformRand();

  // actual spot position by adding a radial vector to a trajectoryPoint
  Gflash3Vector position = point.getPosition() + rShower * std::cos(azimuthalAngle) * point.getOrthogonalUnitVector() +
                           rShower * std::sin(azimuthalAngle) * point.getCrossUnitVector();

  return position;
}