Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:04:01

0001 #include "SimG4CMS/Calo/interface/HFGflash.h"
0002 
0003 #include "G4VPhysicalVolume.hh"
0004 #include "G4Step.hh"
0005 #include "G4Track.hh"
0006 #include "G4Navigator.hh"
0007 #include "G4NavigationHistory.hh"
0008 #include "CLHEP/Units/PhysicalConstants.h"
0009 #include "CLHEP/Units/SystemOfUnits.h"
0010 
0011 #include "Randomize.hh"
0012 #include "G4TransportationManager.hh"
0013 #include "G4VPhysicalVolume.hh"
0014 #include "G4LogicalVolume.hh"
0015 #include "G4VSensitiveDetector.hh"
0016 #include "G4EventManager.hh"
0017 #include "G4SteppingManager.hh"
0018 #include "G4FastTrack.hh"
0019 #include "G4ParticleTable.hh"
0020 
0021 #include "CLHEP/GenericFunctions/IncompleteGamma.hh"
0022 
0023 #include "SimG4Core/Application/interface/SteppingAction.h"
0024 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
0025 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
0026 
0027 #include "FWCore/ServiceRegistry/interface/Service.h"
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029 
0030 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0031 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0032 
0033 #include <cmath>
0034 
0035 //#define EDM_ML_DEBUG
0036 
0037 HFGflash::HFGflash(edm::ParameterSet const& p) {
0038   edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFGflash");
0039   theBField = m_HF.getUntrackedParameter<double>("BField", 3.8);
0040   theWatcherOn = m_HF.getUntrackedParameter<bool>("WatcherOn", true);
0041   theFillHisto = m_HF.getUntrackedParameter<bool>("FillHisto", true);
0042   edm::LogVerbatim("HFShower") << "HFGFlash:: Set B-Field to " << theBField << ", WatcherOn to " << theWatcherOn
0043                                << " and FillHisto to " << theFillHisto;
0044 
0045   theHelix = std::make_unique<GflashTrajectory>();
0046   theGflashStep = std::make_unique<G4Step>();
0047   theGflashNavigator = std::make_unique<G4Navigator>();
0048   theGflashTouchableHandle = std::make_unique<G4TouchableHistory>();
0049 
0050 #ifdef EDM_ML_DEBUG
0051   if (theFillHisto) {
0052     edm::Service<TFileService> tfile;
0053     if (tfile.isAvailable()) {
0054       TFileDirectory showerDir = tfile->mkdir("GflashEMShowerProfile");
0055 
0056       em_incE = showerDir.make<TH1F>("em_incE", "Incoming energy (GeV)", 500, 0, 500.);
0057       em_ssp_rho =
0058           showerDir.make<TH1F>("em_ssp_rho", "Shower starting position;#rho (cm);Number of Events", 100, 100.0, 200.0);
0059       em_ssp_z =
0060           showerDir.make<TH1F>("em_ssp_z", "Shower starting position;z (cm);Number of Events", 2000, 0.0, 2000.0);
0061       em_long =
0062           showerDir.make<TH1F>("em_long", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
0063       em_lateral = showerDir.make<TH1F>("em_lateral", "Lateral Profile;Radiation Length;Moliere Radius", 100, 0.0, 5.0);
0064       em_2d = showerDir.make<TH2F>("em_2d",
0065                                    "Lateral Profile vs. Shower Depth;Radiation Length;Moliere Radius",
0066                                    800,
0067                                    800.0,
0068                                    1600.0,
0069                                    100,
0070                                    0.0,
0071                                    5.0);
0072       em_long_sd = showerDir.make<TH1F>("em_long_sd",
0073                                         "Longitudinal Profile in Sensitive Detector;Radiation Length;Number of Spots",
0074                                         800,
0075                                         800.0,
0076                                         1600.0);
0077       em_lateral_sd =
0078           showerDir.make<TH1F>("em_lateral_sd",
0079                                "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
0080                                100,
0081                                0.0,
0082                                5.0);
0083       em_2d_sd =
0084           showerDir.make<TH2F>("em_2d_sd",
0085                                "Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",
0086                                800,
0087                                800.0,
0088                                1600.0,
0089                                100,
0090                                0.0,
0091                                5.0);
0092       em_ze = showerDir.make<TH2F>(
0093           "em_ze", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 1.0);
0094       em_ratio = showerDir.make<TH2F>(
0095           "em_ratio", "Profile vs. Energy;Radiation Length;Moliere Radius", 800, 800.0, 1600.0, 1000, 0.0, 100.0);
0096       em_ratio_selected = showerDir.make<TH2F>("em_ratio_selected",
0097                                                "Profile vs. Energy;Radiation Length;Moliere Radius",
0098                                                800,
0099                                                800.0,
0100                                                1600.0,
0101                                                1000,
0102                                                0.0,
0103                                                100.0);
0104       em_nSpots_sd =
0105           showerDir.make<TH1F>("em_nSpots_sd",
0106                                "Number of Gflash Spots in Sensitive Detector;Number of Spots;Number of Events",
0107                                100,
0108                                0.0,
0109                                100);
0110       em_ze_ratio = showerDir.make<TH1F>("em_ze_ratio", "Ratio of Energy and Z Position", 1000, 0.0, 0.001);
0111     } else {
0112       theFillHisto = false;
0113       edm::LogVerbatim("HFShower") << "HFGFlash::No file is available for saving"
0114                                    << " histos so the flag is set to false";
0115     }
0116   }
0117 #endif
0118   jCalorimeter = Gflash::kHF;
0119 }
0120 
0121 HFGflash::~HFGflash() {}
0122 
0123 std::vector<HFGflash::Hit> HFGflash::gfParameterization(const G4Step* aStep, bool onlyLong) {
0124   double tempZCalo = 26;            // Gflash::Z[jCalorimeter]
0125   double hfcriticalEnergy = 0.021;  // Gflash::criticalEnergy
0126 
0127   std::vector<HFGflash::Hit> hit;
0128   HFGflash::Hit oneHit;
0129 
0130   auto const preStepPoint = aStep->GetPreStepPoint();
0131   auto const track = aStep->GetTrack();
0132 
0133   // This part of code is copied from the original GFlash Fortran code.
0134   // reference : hep-ex/0001020v1
0135 
0136   const G4double energyCutoff = 1;
0137   const G4int maxNumberOfSpots = 10000000;
0138 
0139   G4ThreeVector showerStartingPosition = track->GetPosition() / cm;
0140   G4ThreeVector showerMomentum = preStepPoint->GetMomentum() / GeV;
0141   jCalorimeter = Gflash::kHF;
0142 
0143   const G4double invgev = 1.0 / CLHEP::GeV;
0144   G4double energy = preStepPoint->GetTotalEnergy() * invgev;  // energy left in GeV
0145   G4double logEinc = std::log(energy);
0146 
0147   G4double y = energy / hfcriticalEnergy;  // y = E/Ec, criticalEnergy is in GeV
0148   G4double logY = std::log(y);
0149 
0150   G4double nSpots = 93.0 * std::log(tempZCalo) * energy;  // total number of spot due linearization
0151   if (energy < 1.6)
0152     nSpots = 140.4 * std::log(tempZCalo) * std::pow(energy, 0.876);
0153 
0154   //   // implementing magnetic field effects
0155   double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
0156   theHelix->initializeTrajectory(showerMomentum, showerStartingPosition, charge, theBField);
0157 
0158   G4double pathLength0 = theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
0159   G4double pathLength = pathLength0;  // this will grow along the shower development
0160 
0161   //--- 2.2  Fix intrinsic properties of em. showers.
0162 
0163   G4double fluctuatedTmax = std::log(logY - 0.7157);
0164   G4double fluctuatedAlpha = std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) * logY);
0165 
0166   G4double sigmaTmax = 1.0 / (-1.4 + 1.26 * logY);
0167   G4double sigmaAlpha = 1.0 / (-0.58 + 0.86 * logY);
0168   G4double rho = 0.705 - 0.023 * logY;
0169   G4double sqrtPL = std::sqrt((1.0 + rho) / 2.0);
0170   G4double sqrtLE = std::sqrt((1.0 - rho) / 2.0);
0171 
0172   G4double norm1 = G4RandGauss::shoot();
0173   G4double norm2 = G4RandGauss::shoot();
0174   G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
0175   G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
0176 
0177   // tmax, alpha, beta : parameters of gamma distribution
0178   G4double tmax = std::exp(tempTmax);
0179   G4double alpha = std::exp(tempAlpha);
0180   G4double beta = (alpha - 1.0) / tmax;
0181 
0182   if (!alpha)
0183     return hit;
0184   if (!beta)
0185     return hit;
0186   if (alpha < 0.00001)
0187     return hit;
0188   if (beta < 0.00001)
0189     return hit;
0190 
0191   // spot fluctuations are added to tmax, alpha, beta
0192   G4double averageTmax = logY - 0.858;
0193   G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) * logY;
0194   G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
0195   G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
0196   G4double spotBeta = (spotAlpha - 1.0) / spotTmax;
0197 
0198   if (!spotAlpha)
0199     return hit;
0200   if (!spotBeta)
0201     return hit;
0202   if (spotAlpha < 0.00001)
0203     return hit;
0204   if (spotBeta < 0.00001)
0205     return hit;
0206 
0207 #ifdef EDM_ML_DEBUG
0208   edm::LogVerbatim("HFShower") << "HFGflash: Incoming energy = " << energy << " Position (rho,z) = ("
0209                                << showerStartingPosition.rho() << ", " << showerStartingPosition.z() << ")";
0210 
0211   if (theFillHisto) {
0212     em_incE->Fill(energy);
0213     em_ssp_rho->Fill(showerStartingPosition.rho());
0214     em_ssp_z->Fill(std::abs(showerStartingPosition.z()));
0215   }
0216 #endif
0217   //  parameters for lateral distribution and fluctuation
0218   G4double z1 = 0.0251 + 0.00319 * logEinc;
0219   G4double z2 = 0.1162 - 0.000381 * tempZCalo;
0220 
0221   G4double k1 = 0.659 - 0.00309 * tempZCalo;
0222   G4double k2 = 0.645;
0223   G4double k3 = -2.59;
0224   G4double k4 = 0.3585 + 0.0421 * logEinc;
0225 
0226   G4double p1 = 2.623 - 0.00094 * tempZCalo;
0227   G4double p2 = 0.401 + 0.00187 * tempZCalo;
0228   G4double p3 = 1.313 - 0.0686 * logEinc;
0229 
0230   //   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data
0231   //   // the width of energy response is not yet tuned.
0232   const G4double e25Scale = 1.03551;
0233   z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790e+00 * std::log(energy) - 3.66237e+00);
0234   p1 *= 0.96;
0235 
0236 #ifdef EDM_ML_DEBUG
0237   G4int nSpots_sd = 0;  // count total number of spots in SD
0238 #endif
0239 
0240   G4double stepLengthLeft = 10000;
0241   G4double zInX0 = 0.0;               // shower depth in X0 unit
0242   G4double deltaZInX0 = 0.0;          // segment of depth in X0 unit
0243   G4double deltaZ = 0.0;              // segment of depth in cm
0244   G4double stepLengthLeftInX0 = 0.0;  // step length left in X0 unit
0245 
0246   const G4double divisionStepInX0 = 0.1;  //step size in X0 unit
0247 
0248   Genfun::IncompleteGamma gammaDist;
0249 
0250   G4double energyInGamma = 0.0;     // energy in a specific depth(z) according to Gamma distribution
0251   G4double preEnergyInGamma = 0.0;  // energy calculated in a previous depth
0252   G4double sigmaInGamma = 0.;       // sigma of energy in a specific depth(z) according to Gamma distribution
0253   G4double preSigmaInGamma = 0.0;   // sigma of energy in a previous depth
0254 
0255   //energy segment in Gamma distribution of shower in each step
0256   G4double deltaEnergy = 0.0;  // energy in deltaZ
0257   G4int spotCounter = 0;       // keep track of number of spots generated
0258 
0259   //step increment along the shower direction
0260   G4double deltaStep = 0.0;
0261 
0262   // Uniqueness of G4Step is important otherwise hits won't be created.
0263   G4double timeGlobal = preStepPoint->GetGlobalTime();
0264 
0265   theGflashNavigator->SetWorldVolume(
0266       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0267 
0268   //   // loop for longitudinal integration
0269 
0270 #ifdef EDM_ML_DEBUG
0271   edm::LogVerbatim("HFShower") << "HFGflash: Energy = " << energy << " Step Length Left = " << stepLengthLeft;
0272 #endif
0273   while (energy > 0.0 && stepLengthLeft > 0.0) {
0274     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
0275 
0276     if (stepLengthLeftInX0 < divisionStepInX0) {
0277       deltaZInX0 = stepLengthLeftInX0;
0278       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0279       stepLengthLeft = 0.0;
0280     } else {
0281       deltaZInX0 = divisionStepInX0;
0282       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0283       stepLengthLeft -= deltaZ;
0284     }
0285 
0286     zInX0 += deltaZInX0;
0287 
0288 #ifdef EDM_ML_DEBUG
0289     edm::LogVerbatim("HFShower") << "HFGflash: zInX0 = " << zInX0 << " spotBeta*zInX0 = " << spotBeta * zInX0;
0290 #endif
0291     if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
0292         (!(zInX0 * beta != 0)) || (zInX0 * beta < 0.00001))
0293       return hit;
0294 
0295     G4int nSpotsInStep = 0;
0296 
0297 #ifdef EDM_ML_DEBUG
0298     edm::LogVerbatim("HFShower") << "HFGflash: Energy - Energy Cut off = " << energy - energyCutoff;
0299 #endif
0300 
0301     if (energy > energyCutoff) {
0302       preEnergyInGamma = energyInGamma;
0303       gammaDist.a().setValue(alpha);  //alpha
0304 
0305       energyInGamma = gammaDist(beta * zInX0);  //beta
0306       G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
0307       deltaEnergy = std::min(energy, energy * energyInDeltaZ);
0308 
0309       preSigmaInGamma = sigmaInGamma;
0310       gammaDist.a().setValue(spotAlpha);           //alpha spot
0311       sigmaInGamma = gammaDist(spotBeta * zInX0);  //beta spot
0312       nSpotsInStep = std::max(1, int(nSpots * (sigmaInGamma - preSigmaInGamma)));
0313     } else {
0314       deltaEnergy = energy;
0315       preSigmaInGamma = sigmaInGamma;
0316       nSpotsInStep = std::max(1, int(nSpots * (1.0 - preSigmaInGamma)));
0317     }
0318 
0319     if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
0320       deltaEnergy = energy;
0321 
0322     energy -= deltaEnergy;
0323 
0324     if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
0325       nSpotsInStep = maxNumberOfSpots - spotCounter;
0326       if (nSpotsInStep < 1)
0327         nSpotsInStep = 1;
0328     }
0329 
0330     //     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
0331     deltaStep += 0.5 * deltaZ;
0332     pathLength += deltaStep;
0333     deltaStep = 0.5 * deltaZ;
0334 
0335     //lateral shape and fluctuations for  homogenous calo.
0336     G4double tScale = tmax * alpha / (alpha - 1.0) * (1.0 - std::exp(-fluctuatedAlpha));
0337     G4double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
0338     G4double rCore = z1 + z2 * tau;
0339     G4double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));  // @@ check RT3 sign
0340     G4double p23 = (p2 - tau) / p3;
0341     G4double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));
0342 
0343     // Deposition of spots according to lateral distr.
0344     // Apply absolute energy scale
0345     // Convert into MeV unit
0346     G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
0347 
0348 #ifdef EDM_ML_DEBUG
0349     edm::LogVerbatim("HFShower") << "HFGflash: nSpotsInStep = " << nSpotsInStep;
0350 #endif
0351 
0352     for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
0353       spotCounter++;
0354       G4double u1 = G4UniformRand();
0355       G4double u2 = G4UniformRand();
0356       G4double rInRM = 0.0;
0357 
0358       if (u1 < probabilityWeight) {
0359         rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
0360       } else {
0361         rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
0362       }
0363 
0364       G4double rShower = rInRM * Gflash::rMoliere[jCalorimeter];
0365 
0366       //Uniform & random rotation of spot along the azimuthal angle
0367       G4double azimuthalAngle = twopi * G4UniformRand();
0368 
0369       //Compute global position of generated spots with taking into account magnetic field
0370       //Divide deltaZ into nSpotsInStep and give a spot a global position
0371       G4double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0372 
0373       // trajectoryPoint give a spot an imaginary point along the shower development
0374       GflashTrajectoryPoint trajectoryPoint;
0375       theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
0376 
0377       G4ThreeVector SpotPosition0 = trajectoryPoint.getPosition() +
0378                                     rShower * std::cos(azimuthalAngle) * trajectoryPoint.getOrthogonalUnitVector() +
0379                                     rShower * std::sin(azimuthalAngle) * trajectoryPoint.getCrossUnitVector();
0380 
0381       //!V.Ivanchenko - not clear if it is correct
0382       // Convert into mm unit
0383       SpotPosition0 *= cm;
0384 
0385       //---------------------------------------------------
0386       // fill a fake step to send it to hit maker
0387       //---------------------------------------------------
0388 
0389       // to make a different time for each fake step. (0.03 nsec is corresponding to 1cm step size)
0390       timeGlobal += 0.0001 * nanosecond;
0391 
0392       //fill equivalent changes to a (fake) step associated with a spot
0393 
0394       G4double zInX0_spot = std::abs(pathLength + incrementPath - pathLength0) / Gflash::radLength[jCalorimeter];
0395 
0396 #ifdef EDM_ML_DEBUG
0397       edm::LogVerbatim("HFShower") << "HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot << " , "
0398                                    << emSpotEnergy / GeV << "emSpotEnergy/GeV =" << emSpotEnergy / GeV;
0399 #endif
0400 
0401       if (zInX0_spot <= 0)
0402         continue;
0403       if (emSpotEnergy <= 0)
0404         continue;
0405       if (rShower / Gflash::rMoliere[jCalorimeter] <= 0)
0406         continue;
0407 
0408       oneHit.depth = 1;
0409 
0410 #ifdef EDM_ML_DEBUG
0411       if (theFillHisto) {
0412         em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy * invgev);
0413         em_lateral->Fill(rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0414         em_2d->Fill(SpotPosition0.z() / cm, rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0415       }
0416 #endif
0417 
0418       //!V.Ivanchenko what this means??
0419       //if(SpotPosition0 == 0) continue;
0420 
0421       double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
0422 
0423       if (emSpotEnergy / GeV < 0.0001)
0424         continue;
0425       if (energyratio > 80)
0426         continue;
0427 
0428       double zshift = 0;
0429       if (SpotPosition0.z() > 0)
0430         zshift = 18;
0431       if (SpotPosition0.z() < 0)
0432         zshift = -18;
0433 
0434       G4ThreeVector gfshift(0, 0, zshift * (pow(100, 0.1) / pow(energy, 0.1)));
0435 
0436       G4ThreeVector SpotPosition = gfshift + SpotPosition0;
0437 
0438       double LengthWeight = std::fabs(std::pow(SpotPosition0.z() / 11370, 1));
0439       if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
0440         continue;
0441 
0442       oneHit.position = SpotPosition;
0443       oneHit.time = timeGlobal;
0444       oneHit.edep = emSpotEnergy * invgev;
0445       hit.push_back(oneHit);
0446 #ifdef EDM_ML_DEBUG
0447       nSpots_sd++;
0448 #endif
0449 
0450     }  // end of for spot iteration
0451 
0452   }  // end of while for longitudinal integration
0453 #ifdef EDM_ML_DEBUG
0454   if (theFillHisto) {
0455     em_nSpots_sd->Fill(nSpots_sd);
0456   }
0457 #endif
0458   return hit;
0459 }