Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:16 

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