HFGflash.cc

Go to the documentation of this file.

#include "SimG4CMS/Calo/interface/HFGflash.h"

#include "G4VPhysicalVolume.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4Navigator.hh"
#include "G4NavigationHistory.hh"
#include "CLHEP/Units/PhysicalConstants.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include "Randomize.hh"
#include "G4TransportationManager.hh"
#include "G4VPhysicalVolume.hh"
#include "G4LogicalVolume.hh"
#include "G4VSensitiveDetector.hh"
#include "G4EventManager.hh"
#include "G4SteppingManager.hh"
#include "G4FastTrack.hh"
#include "G4ParticleTable.hh"

#include "CLHEP/GenericFunctions/IncompleteGamma.hh"

#include "SimG4Core/Application/interface/SteppingAction.h"
#include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
#include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

#include <cmath>

//#define EDM_ML_DEBUG

HFGflash::HFGflash(edm::ParameterSet const& p) {
  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFGflash");
  theBField = m_HF.getUntrackedParameter<double>("BField", 3.8);
  theWatcherOn = m_HF.getUntrackedParameter<bool>("WatcherOn", true);
  theFillHisto = m_HF.getUntrackedParameter<bool>("FillHisto", true);
  edm::LogVerbatim("HFShower") << "HFGFlash:: Set B-Field to " << theBField << ", WatcherOn to " << theWatcherOn
                               << " and FillHisto to " << theFillHisto;

  theHelix = std::make_unique<GflashTrajectory>();
  theGflashStep = std::make_unique<G4Step>();
  theGflashNavigator = std::make_unique<G4Navigator>();
  theGflashTouchableHandle = std::make_unique<G4TouchableHistory>();

#ifdef EDM_ML_DEBUG
  if (theFillHisto) {
    edm::Service<TFileService> tfile;
    if (tfile.isAvailable()) {
      TFileDirectory showerDir = tfile->mkdir("GflashEMShowerProfile");

      em_incE = showerDir.make<TH1F>("em_incE", "Incoming energy (GeV)", 500, 0, 500.);
      em_ssp_rho =
          showerDir.make<TH1F>("em_ssp_rho", "Shower starting position;#rho (cm);Number of Events", 100, 100.0, 200.0);
      em_ssp_z =
          showerDir.make<TH1F>("em_ssp_z", "Shower starting position;z (cm);Number of Events", 2000, 0.0, 2000.0);
      em_long =
          showerDir.make<TH1F>("em_long", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
      em_lateral = showerDir.make<TH1F>("em_lateral", "Lateral Profile;Radiation Length;Moliere Radius", 100, 0.0, 5.0);
      em_2d = showerDir.make<TH2F>("em_2d",
                                   "Lateral Profile vs. Shower Depth;Radiation Length;Moliere Radius",
                                   800,
                                   800.0,
                                   1600.0,
                                   100,
                                   0.0,
                                   5.0);
      em_long_sd = showerDir.make<TH1F>("em_long_sd",
                                        "Longitudinal Profile in Sensitive Detector;Radiation Length;Number of Spots",
                                        800,
                                        800.0,
                                        1600.0);
      em_lateral_sd =
          showerDir.make<TH1F>("em_lateral_sd",
                               "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
                               100,
                               0.0,
                               5.0);
      em_2d_sd =
          showerDir.make<TH2F>("em_2d_sd",
                               "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
                               800,
                               800.0,
                               1600.0,
                               100,
                               0.0,
                               5.0);
      em_ze = showerDir.make<TH2F>(
          "em_ze", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 1.0);
      em_ratio = showerDir.make<TH2F>(
          "em_ratio", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 100.0);
      em_ratio_selected = showerDir.make<TH2F>("em_ratio_selected",
                                               "Profile vs. Energy;Radiation Length;Moliere Radius",
                                               800,
                                               800.0,
                                               1600.0,
                                               1000,
                                               0.0,
                                               100.0);
      em_nSpots_sd =
          showerDir.make<TH1F>("em_nSpots_sd",
                               "Number of Gflash Spots in Sensitive Detector;Number of Spots;Number of Events",
                               100,
                               0.0,
                               100);
      em_ze_ratio = showerDir.make<TH1F>("em_ze_ratio", "Ratio of Energy and Z Position", 1000, 0.0, 0.001);
    } else {
      theFillHisto = false;
      edm::LogVerbatim("HFShower") << "HFGFlash::No file is available for saving"
                                   << " histos so the flag is set to false";
    }
  }
#endif
  jCalorimeter = Gflash::kHF;
}

HFGflash::~HFGflash() {}

std::vector<HFGflash::Hit> HFGflash::gfParameterization(const G4Step* aStep, bool onlyLong) {
  double tempZCalo = 26;            // Gflash::Z[jCalorimeter]
  double hfcriticalEnergy = 0.021;  // Gflash::criticalEnergy

  std::vector<HFGflash::Hit> hit;
  HFGflash::Hit oneHit;

  auto const preStepPoint = aStep->GetPreStepPoint();
  auto const track = aStep->GetTrack();

  // This part of code is copied from the original GFlash Fortran code.
  // reference : hep-ex/0001020v1

  const G4double energyCutoff = 1;
  const G4int maxNumberOfSpots = 10000000;

  G4ThreeVector showerStartingPosition = track->GetPosition() / cm;
  G4ThreeVector showerMomentum = preStepPoint->GetMomentum() / GeV;
  jCalorimeter = Gflash::kHF;

  const G4double invgev = 1.0 / CLHEP::GeV;
  G4double energy = preStepPoint->GetTotalEnergy() * invgev;  // energy left in GeV
  G4double logEinc = std::log(energy);

  G4double y = energy / hfcriticalEnergy;  // y = E/Ec, criticalEnergy is in GeV
  G4double logY = std::log(y);

  G4double nSpots = 93.0 * std::log(tempZCalo) * energy;  // total number of spot due linearization
  if (energy < 1.6)
    nSpots = 140.4 * std::log(tempZCalo) * std::pow(energy, 0.876);

  //   // implementing magnetic field effects
  double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
  theHelix->initializeTrajectory(showerMomentum, showerStartingPosition, charge, theBField);

  G4double pathLength0 = theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
  G4double pathLength = pathLength0;  // this will grow along the shower development

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

  G4double fluctuatedTmax = std::log(logY - 0.7157);
  G4double fluctuatedAlpha = std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) * logY);

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

  G4double norm1 = G4RandGauss::shoot();
  G4double norm2 = G4RandGauss::shoot();
  G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
  G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);

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

  if (!alpha)
    return hit;
  if (!beta)
    return hit;
  if (alpha < 0.00001)
    return hit;
  if (beta < 0.00001)
    return hit;

  // spot fluctuations are added to tmax, alpha, beta
  G4double averageTmax = logY - 0.858;
  G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) * logY;
  G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
  G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
  G4double spotBeta = (spotAlpha - 1.0) / spotTmax;

  if (!spotAlpha)
    return hit;
  if (!spotBeta)
    return hit;
  if (spotAlpha < 0.00001)
    return hit;
  if (spotBeta < 0.00001)
    return hit;

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFGflash: Incoming energy = " << energy << " Position (rho,z) = ("
                               << showerStartingPosition.rho() << ", " << showerStartingPosition.z() << ")";

  if (theFillHisto) {
    em_incE->Fill(energy);
    em_ssp_rho->Fill(showerStartingPosition.rho());
    em_ssp_z->Fill(std::abs(showerStartingPosition.z()));
  }
#endif
  //  parameters for lateral distribution and fluctuation
  G4double z1 = 0.0251 + 0.00319 * logEinc;
  G4double z2 = 0.1162 - 0.000381 * tempZCalo;

  G4double k1 = 0.659 - 0.00309 * tempZCalo;
  G4double k2 = 0.645;
  G4double k3 = -2.59;
  G4double k4 = 0.3585 + 0.0421 * logEinc;

  G4double p1 = 2.623 - 0.00094 * tempZCalo;
  G4double p2 = 0.401 + 0.00187 * tempZCalo;
  G4double p3 = 1.313 - 0.0686 * logEinc;

  //   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data
  //   // the width of energy response is not yet tuned.
  const G4double e25Scale = 1.03551;
  z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790e+00 * std::log(energy) - 3.66237e+00);
  p1 *= 0.96;

#ifdef EDM_ML_DEBUG
  G4int nSpots_sd = 0;  // count total number of spots in SD
#endif

  G4double stepLengthLeft = 10000;
  G4double zInX0 = 0.0;               // shower depth in X0 unit
  G4double deltaZInX0 = 0.0;          // segment of depth in X0 unit
  G4double deltaZ = 0.0;              // segment of depth in cm
  G4double stepLengthLeftInX0 = 0.0;  // step length left in X0 unit

  const G4double divisionStepInX0 = 0.1;  //step size in X0 unit

  Genfun::IncompleteGamma gammaDist;

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

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

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

  // Uniqueness of G4Step is important otherwise hits won't be created.
  G4double timeGlobal = preStepPoint->GetGlobalTime();

  theGflashNavigator->SetWorldVolume(
      G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());

  //   // loop for longitudinal integration

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFGflash: Energy = " << energy << " Step Length Left = " << stepLengthLeft;
#endif
  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;

#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFGflash: zInX0 = " << zInX0 << " spotBeta*zInX0 = " << spotBeta * zInX0;
#endif
    if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
        (!(zInX0 * beta != 0)) || (zInX0 * beta < 0.00001))
      return hit;

    G4int nSpotsInStep = 0;

#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFGflash: Energy - Energy Cut off = " << energy - energyCutoff;
#endif

    if (energy > energyCutoff) {
      preEnergyInGamma = energyInGamma;
      gammaDist.a().setValue(alpha);  //alpha

      energyInGamma = gammaDist(beta * zInX0);  //beta
      G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
      deltaEnergy = std::min(energy, energy * 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)
        nSpotsInStep = 1;
    }

    //     // 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.
    G4double tScale = tmax * alpha / (alpha - 1.0) * (1.0 - std::exp(-fluctuatedAlpha));
    G4double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
    G4double rCore = z1 + z2 * tau;
    G4double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));  // @@ check RT3 sign
    G4double p23 = (p2 - tau) / p3;
    G4double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));

    // Deposition of spots according to lateral distr.
    // Apply absolute energy scale
    // Convert into MeV unit
    G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);

#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFGflash: nSpotsInStep = " << nSpotsInStep;
#endif

    for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
      spotCounter++;
      G4double u1 = G4UniformRand();
      G4double u2 = G4UniformRand();
      G4double rInRM = 0.0;

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

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

      //Uniform & random rotation of spot along the azimuthal angle
      G4double azimuthalAngle = twopi * G4UniformRand();

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

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

      G4ThreeVector SpotPosition0 = trajectoryPoint.getPosition() +
                                    rShower * std::cos(azimuthalAngle) * trajectoryPoint.getOrthogonalUnitVector() +
                                    rShower * std::sin(azimuthalAngle) * trajectoryPoint.getCrossUnitVector();

      // Convert into mm unit
      SpotPosition0 *= cm;

      //---------------------------------------------------
      // fill a fake step to send it to hit maker
      //---------------------------------------------------

      // to make a different time for each fake step. (0.03 nsec is corresponding to 1cm step size)
      timeGlobal += 0.0001 * nanosecond;

      //fill equivalent changes to a (fake) step associated with a spot

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

#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HFShower") << "HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot << " , "
                                   << emSpotEnergy / GeV << "emSpotEnergy/GeV =" << emSpotEnergy / GeV;
#endif

      if (zInX0_spot <= 0)
        continue;
      if (emSpotEnergy <= 0)
        continue;
      if (rShower / Gflash::rMoliere[jCalorimeter] <= 0)
        continue;

      oneHit.depth = 1;

#ifdef EDM_ML_DEBUG
      if (theFillHisto) {
        em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy * invgev);
        em_lateral->Fill(rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
        em_2d->Fill(SpotPosition0.z() / cm, rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
      }
#endif

      //if(SpotPosition0 == 0) continue;

      double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));

      if (emSpotEnergy / GeV < 0.0001)
        continue;
      if (energyratio > 80)
        continue;

      double zshift = 0;
      if (SpotPosition0.z() > 0)
        zshift = 18;
      if (SpotPosition0.z() < 0)
        zshift = -18;

      G4ThreeVector gfshift(0, 0, zshift * (pow(100, 0.1) / pow(energy, 0.1)));

      G4ThreeVector SpotPosition = gfshift + SpotPosition0;

      double LengthWeight = std::fabs(std::pow(SpotPosition0.z() / 11370, 1));
      if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
        continue;

      oneHit.position = SpotPosition;
      oneHit.time = timeGlobal;
      oneHit.edep = emSpotEnergy * invgev;
      hit.push_back(oneHit);
#ifdef EDM_ML_DEBUG
      nSpots_sd++;
#endif

    }  // end of for spot iteration

  }  // end of while for longitudinal integration
#ifdef EDM_ML_DEBUG
  if (theFillHisto) {
    em_nSpots_sd->Fill(nSpots_sd);
  }
#endif
  return hit;
}