Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-06 04:01:07 

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   G4double stepLengthLeft = 10000;
0237   G4int nSpots_sd = 0;                // count total number of spots in SD
0238   G4double zInX0 = 0.0;               // shower depth in X0 unit
0239   G4double deltaZInX0 = 0.0;          // segment of depth in X0 unit
0240   G4double deltaZ = 0.0;              // segment of depth in cm
0241   G4double stepLengthLeftInX0 = 0.0;  // step length left in X0 unit
0242 
0243   const G4double divisionStepInX0 = 0.1;  //step size in X0 unit
0244 
0245   Genfun::IncompleteGamma gammaDist;
0246 
0247   G4double energyInGamma = 0.0;     // energy in a specific depth(z) according to Gamma distribution
0248   G4double preEnergyInGamma = 0.0;  // energy calculated in a previous depth
0249   G4double sigmaInGamma = 0.;       // sigma of energy in a specific depth(z) according to Gamma distribution
0250   G4double preSigmaInGamma = 0.0;   // sigma of energy in a previous depth
0251 
0252   //energy segment in Gamma distribution of shower in each step
0253   G4double deltaEnergy = 0.0;  // energy in deltaZ
0254   G4int spotCounter = 0;       // keep track of number of spots generated
0255 
0256   //step increment along the shower direction
0257   G4double deltaStep = 0.0;
0258 
0259   // Uniqueness of G4Step is important otherwise hits won't be created.
0260   G4double timeGlobal = preStepPoint->GetGlobalTime();
0261 
0262   theGflashNavigator->SetWorldVolume(
0263       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0264 
0265   //   // loop for longitudinal integration
0266 
0267 #ifdef EDM_ML_DEBUG
0268   edm::LogVerbatim("HFShower") << "HFGflash: Energy = " << energy << " Step Length Left = " << stepLengthLeft;
0269 #endif
0270   while (energy > 0.0 && stepLengthLeft > 0.0) {
0271     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
0272 
0273     if (stepLengthLeftInX0 < divisionStepInX0) {
0274       deltaZInX0 = stepLengthLeftInX0;
0275       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0276       stepLengthLeft = 0.0;
0277     } else {
0278       deltaZInX0 = divisionStepInX0;
0279       deltaZ = deltaZInX0 * Gflash::radLength[jCalorimeter];
0280       stepLengthLeft -= deltaZ;
0281     }
0282 
0283     zInX0 += deltaZInX0;
0284 
0285 #ifdef EDM_ML_DEBUG
0286     edm::LogVerbatim("HFShower") << "HFGflash: zInX0 = " << zInX0 << " spotBeta*zInX0 = " << spotBeta * zInX0;
0287 #endif
0288     if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
0289         (!(zInX0 * beta != 0)) || (zInX0 * beta < 0.00001))
0290       return hit;
0291 
0292     G4int nSpotsInStep = 0;
0293 
0294 #ifdef EDM_ML_DEBUG
0295     edm::LogVerbatim("HFShower") << "HFGflash: Energy - Energy Cut off = " << energy - energyCutoff;
0296 #endif
0297 
0298     if (energy > energyCutoff) {
0299       preEnergyInGamma = energyInGamma;
0300       gammaDist.a().setValue(alpha);  //alpha
0301 
0302       energyInGamma = gammaDist(beta * zInX0);  //beta
0303       G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
0304       deltaEnergy = std::min(energy, energy * energyInDeltaZ);
0305 
0306       preSigmaInGamma = sigmaInGamma;
0307       gammaDist.a().setValue(spotAlpha);           //alpha spot
0308       sigmaInGamma = gammaDist(spotBeta * zInX0);  //beta spot
0309       nSpotsInStep = std::max(1, int(nSpots * (sigmaInGamma - preSigmaInGamma)));
0310     } else {
0311       deltaEnergy = energy;
0312       preSigmaInGamma = sigmaInGamma;
0313       nSpotsInStep = std::max(1, int(nSpots * (1.0 - preSigmaInGamma)));
0314     }
0315 
0316     if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
0317       deltaEnergy = energy;
0318 
0319     energy -= deltaEnergy;
0320 
0321     if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
0322       nSpotsInStep = maxNumberOfSpots - spotCounter;
0323       if (nSpotsInStep < 1)
0324         nSpotsInStep = 1;
0325     }
0326 
0327     //     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
0328     deltaStep += 0.5 * deltaZ;
0329     pathLength += deltaStep;
0330     deltaStep = 0.5 * deltaZ;
0331 
0332     //lateral shape and fluctuations for  homogenous calo.
0333     G4double tScale = tmax * alpha / (alpha - 1.0) * (1.0 - std::exp(-fluctuatedAlpha));
0334     G4double tau = std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
0335     G4double rCore = z1 + z2 * tau;
0336     G4double rTail = k1 * (std::exp(k3 * (tau - k2)) + std::exp(k4 * (tau - k2)));  // @@ check RT3 sign
0337     G4double p23 = (p2 - tau) / p3;
0338     G4double probabilityWeight = p1 * std::exp(p23 - std::exp(p23));
0339 
0340     // Deposition of spots according to lateral distr.
0341     // Apply absolute energy scale
0342     // Convert into MeV unit
0343     G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
0344 
0345 #ifdef EDM_ML_DEBUG
0346     edm::LogVerbatim("HFShower") << "HFGflash: nSpotsInStep = " << nSpotsInStep;
0347 #endif
0348 
0349     for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
0350       spotCounter++;
0351       G4double u1 = G4UniformRand();
0352       G4double u2 = G4UniformRand();
0353       G4double rInRM = 0.0;
0354 
0355       if (u1 < probabilityWeight) {
0356         rInRM = rCore * std::sqrt(u2 / (1.0 - u2));
0357       } else {
0358         rInRM = rTail * std::sqrt(u2 / (1.0 - u2));
0359       }
0360 
0361       G4double rShower = rInRM * Gflash::rMoliere[jCalorimeter];
0362 
0363       //Uniform & random rotation of spot along the azimuthal angle
0364       G4double azimuthalAngle = twopi * G4UniformRand();
0365 
0366       //Compute global position of generated spots with taking into account magnetic field
0367       //Divide deltaZ into nSpotsInStep and give a spot a global position
0368       G4double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
0369 
0370       // trajectoryPoint give a spot an imaginary point along the shower development
0371       GflashTrajectoryPoint trajectoryPoint;
0372       theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
0373 
0374       G4ThreeVector SpotPosition0 = trajectoryPoint.getPosition() +
0375                                     rShower * std::cos(azimuthalAngle) * trajectoryPoint.getOrthogonalUnitVector() +
0376                                     rShower * std::sin(azimuthalAngle) * trajectoryPoint.getCrossUnitVector();
0377 
0378       //!V.Ivanchenko - not clear if it is correct
0379       // Convert into mm unit
0380       SpotPosition0 *= cm;
0381 
0382       //---------------------------------------------------
0383       // fill a fake step to send it to hit maker
0384       //---------------------------------------------------
0385 
0386       // to make a different time for each fake step. (0.03 nsec is corresponding to 1cm step size)
0387       timeGlobal += 0.0001 * nanosecond;
0388 
0389       //fill equivalent changes to a (fake) step associated with a spot
0390 
0391       G4double zInX0_spot = std::abs(pathLength + incrementPath - pathLength0) / Gflash::radLength[jCalorimeter];
0392 
0393 #ifdef EDM_ML_DEBUG
0394       edm::LogVerbatim("HFShower") << "HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot << " , "
0395                                    << emSpotEnergy / GeV << "emSpotEnergy/GeV =" << emSpotEnergy / GeV;
0396 #endif
0397 
0398       if (zInX0_spot <= 0)
0399         continue;
0400       if (emSpotEnergy <= 0)
0401         continue;
0402       if (rShower / Gflash::rMoliere[jCalorimeter] <= 0)
0403         continue;
0404 
0405       oneHit.depth = 1;
0406 
0407 #ifdef EDM_ML_DEBUG
0408       if (theFillHisto) {
0409         em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy * invgev);
0410         em_lateral->Fill(rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0411         em_2d->Fill(SpotPosition0.z() / cm, rShower / Gflash::rMoliere[jCalorimeter], emSpotEnergy * invgev);
0412       }
0413 #endif
0414 
0415       //!V.Ivanchenko what this means??
0416       //if(SpotPosition0 == 0) continue;
0417 
0418       double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
0419 
0420       if (emSpotEnergy / GeV < 0.0001)
0421         continue;
0422       if (energyratio > 80)
0423         continue;
0424 
0425       double zshift = 0;
0426       if (SpotPosition0.z() > 0)
0427         zshift = 18;
0428       if (SpotPosition0.z() < 0)
0429         zshift = -18;
0430 
0431       G4ThreeVector gfshift(0, 0, zshift * (pow(100, 0.1) / pow(energy, 0.1)));
0432 
0433       G4ThreeVector SpotPosition = gfshift + SpotPosition0;
0434 
0435       double LengthWeight = std::fabs(std::pow(SpotPosition0.z() / 11370, 1));
0436       if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
0437         continue;
0438 
0439       oneHit.position = SpotPosition;
0440       oneHit.time = timeGlobal;
0441       oneHit.edep = emSpotEnergy * invgev;
0442       hit.push_back(oneHit);
0443       nSpots_sd++;
0444 
0445     }  // end of for spot iteration
0446 
0447   }  // end of while for longitudinal integration
0448 #ifdef EDM_ML_DEBUG
0449   if (theFillHisto) {
0450     em_nSpots_sd->Fill(nSpots_sd);
0451   }
0452 #endif
0453   return hit;
0454 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!