Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFShowerParam.cc
0003 // Description: Parametrized version of HF hits
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HFShowerParam.h"
0007 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
0008 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0009 
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0013 
0014 #include "G4VPhysicalVolume.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "G4NavigationHistory.hh"
0018 #include "Randomize.hh"
0019 
0020 #include <CLHEP/Units/PhysicalConstants.h>
0021 #include <CLHEP/Units/SystemOfUnits.h>
0022 
0023 #include <iostream>
0024 
0025 //#define EDM_ML_DEBUG
0026 //#define plotDebug
0027 //#define mkdebug
0028 
0029 HFShowerParam::HFShowerParam(const std::string& name,
0030                              const HcalDDDSimConstants* hcons,
0031                              const HcalSimulationParameters* hps,
0032                              edm::ParameterSet const& p)
0033     : fibre_(hcons, hps, p), fillHisto_(false) {
0034   edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
0035   edm::ParameterSet m_HF2 = m_HF.getParameter<edm::ParameterSet>("HFShowerBlock");
0036   pePerGeV_ = m_HF.getParameter<double>("PEPerGeV");
0037   trackEM_ = m_HF.getParameter<bool>("TrackEM");
0038   bool useShowerLibrary = m_HF.getParameter<bool>("UseShowerLibrary");
0039   bool useGflash = m_HF.getParameter<bool>("UseHFGflash");
0040   edMin_ = m_HF.getParameter<double>("EminLibrary");
0041   equalizeTimeShift_ = m_HF2.getParameter<bool>("EqualizeTimeShift");
0042   onlyLong_ = m_HF2.getParameter<bool>("OnlyLong");
0043   ref_index_ = m_HF.getParameter<double>("RefIndex");
0044   double lambdaMean = m_HF.getParameter<double>("LambdaMean");
0045   aperture_ = cos(asin(m_HF.getParameter<double>("Aperture")));
0046   applyFidCut_ = m_HF.getParameter<bool>("ApplyFiducialCut");
0047   parametrizeLast_ = m_HF.getUntrackedParameter<bool>("ParametrizeLast", false);
0048   gpar_ = hcons->getGparHF();
0049 
0050   edm::LogVerbatim("HFShower") << "HFShowerParam::Use of shower library is set to " << useShowerLibrary
0051                                << " Use of Gflash is set to " << useGflash << " P.E. per GeV " << pePerGeV_
0052                                << ", ref. index of fibre " << ref_index_ << ", Track EM Flag " << trackEM_ << ", edMin "
0053                                << edMin_ << " GeV, use of Short fibre info in"
0054                                << " shower library set to " << !(onlyLong_)
0055                                << " equalize flag for time shift in fibre is set to " << equalizeTimeShift_
0056                                << ", use of parametrization for last part set to " << parametrizeLast_
0057                                << ", Mean lambda " << lambdaMean << ", aperture (cutoff) " << aperture_
0058                                << ", Application of Fiducial Cut " << applyFidCut_;
0059 
0060 #ifdef plotDebug
0061   edm::Service<TFileService> tfile;
0062   if (tfile.isAvailable()) {
0063     fillHisto_ = true;
0064     edm::LogVerbatim("HFShower") << "HFShowerParam::Save histos in directory "
0065                                  << "ProfileFromParam";
0066     TFileDirectory showerDir = tfile->mkdir("ProfileFromParam");
0067     hzvem_ = showerDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part);Number of PE", 330, 0.0, 1650.0);
0068     hzvhad_ = showerDir.make<TH1F>("hzvhad", "Longitudinal Profile (Had Part);Number of PE", 330, 0.0, 1650.0);
0069     em_2d_1_ = showerDir.make<TH2F>(
0070         "em_2d_1", "Lateral Profile vs. Shower Depth;cm;Events", 800, 800.0, 1600.0, 100, 50.0, 150.0);
0071     em_long_1_ =
0072         showerDir.make<TH1F>("em_long_1", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
0073     em_long_1_tuned_ = showerDir.make<TH1F>(
0074         "em_long_1_tuned", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
0075     em_lateral_1_ = showerDir.make<TH1F>("em_lateral_1", "Lateral Profile;cm;Events", 100, 50.0, 150.0);
0076     em_2d_2_ = showerDir.make<TH2F>(
0077         "em_2d_2", "Lateral Profile vs. Shower Depth;cm;Events", 800, 800.0, 1600.0, 100, 50.0, 150.0);
0078     em_long_2_ =
0079         showerDir.make<TH1F>("em_long_2", "Longitudinal Profile;Radiation Length;Number of Spots", 800, 800.0, 1600.0);
0080     em_lateral_2_ = showerDir.make<TH1F>("em_lateral_2", "Lateral Profile;cm;Events", 100, 50.0, 150.0);
0081     em_long_gflash_ = showerDir.make<TH1F>(
0082         "em_long_gflash", "Longitudinal Profile From GFlash;cm;Number of Spots", 800, 800.0, 1600.0);
0083     em_long_sl_ = showerDir.make<TH1F>(
0084         "em_long_sl", "Longitudinal Profile From Shower Library;cm;Number of Spots", 800, 800.0, 1600.0);
0085   } else {
0086     fillHisto_ = false;
0087     edm::LogVerbatim("HFShower") << "HFShowerParam::No file is available for saving histos so the "
0088                                  << "flag is set to false";
0089   }
0090 #endif
0091 
0092   if (useShowerLibrary)
0093     showerLibrary_ = std::make_unique<HFShowerLibrary>(name, hcons, hps, p);
0094   else
0095     showerLibrary_.reset(nullptr);
0096   if (useGflash)
0097     gflash_ = std::make_unique<HFGflash>(p);
0098   else
0099     gflash_.reset(nullptr);
0100   attLMeanInv_ = fibre_.attLength(lambdaMean);
0101   edm::LogVerbatim("HFShower") << "att. length used for (lambda=" << lambdaMean
0102                                << ") = " << 1 / (attLMeanInv_ * CLHEP::cm) << " cm";
0103 }
0104 
0105 std::vector<HFShowerParam::Hit> HFShowerParam::getHits(const G4Step* aStep, double weight, bool& isKilled) {
0106   auto const preStepPoint = aStep->GetPreStepPoint();
0107   auto const track = aStep->GetTrack();
0108   bool isEM = G4TrackToParticleID::isGammaElectronPositron(track);
0109   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0110   double zv = std::abs(hitPoint.z()) - gpar_[4] - 0.5 * gpar_[1];
0111   G4ThreeVector localPoint = G4ThreeVector(hitPoint.x(), hitPoint.y(), zv);
0112 
0113   double pin = (preStepPoint->GetTotalEnergy()) / CLHEP::GeV;
0114   double zint = hitPoint.z();
0115   double zz = std::abs(zint) - gpar_[4];
0116 
0117 #ifdef EDM_ML_DEBUG
0118   edm::LogVerbatim("HFShower") << "HFShowerParam: getHits " << track->GetDefinition()->GetParticleName()
0119                                << " of energy " << pin << " GeV Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y()
0120                                << "," << zint << " (" << zz << "," << localPoint.z() << ", "
0121                                << (localPoint.z() + 0.5 * gpar_[1]) << ") Local " << localPoint;
0122 #endif
0123   std::vector<HFShowerParam::Hit> hits;
0124   HFShowerParam::Hit hit;
0125   hit.position = hitPoint;
0126 
0127   // look for other charged particles
0128   bool other = false;
0129   double pBeta = track->GetDynamicParticle()->GetTotalMomentum() / track->GetDynamicParticle()->GetTotalEnergy();
0130   double dirz = (track->GetDynamicParticle()->GetMomentumDirection()).z();
0131   if (hitPoint.z() < 0)
0132     dirz *= -1.;
0133 #ifdef EDM_ML_DEBUG
0134   edm::LogVerbatim("HFShower") << "HFShowerParam: getHits Momentum "
0135                                << track->GetDynamicParticle()->GetMomentumDirection() << " HitPoint " << hitPoint
0136                                << " dirz " << dirz;
0137 #endif
0138   if (!isEM && track->GetDefinition()->GetPDGCharge() != 0 && pBeta > (1 / ref_index_) &&
0139       aStep->GetTotalEnergyDeposit() > 0.) {
0140     other = true;
0141   }
0142 
0143   // take only e+-/gamma/or special particles
0144   if (isEM || other) {
0145     // Leave out the last part
0146     double edep = 0.;
0147     if ((!trackEM_) && ((zz < (gpar_[1] - gpar_[2])) || parametrizeLast_) && (!other)) {
0148       edep = pin;
0149       isKilled = true;
0150     } else if ((track->GetDefinition()->GetPDGCharge() != 0) && (pBeta > (1 / ref_index_)) && (dirz > aperture_)) {
0151       edep = (aStep->GetTotalEnergyDeposit()) / CLHEP::GeV;
0152     }
0153     std::string path = "ShowerLibrary";
0154 #ifdef EDM_ML_DEBUG
0155     edm::LogVerbatim("HFShower") << "HFShowerParam: getHits edep = " << edep << " weight " << weight << " final "
0156                                  << edep * weight << ", Kill = " << isKilled << ", pin = " << pin
0157                                  << ", edMin = " << edMin_ << " Other " << other;
0158 #endif
0159     edep *= weight;
0160     if (edep > 0) {
0161       if ((showerLibrary_.get() || gflash_.get()) && isKilled && pin > edMin_ && (!other)) {
0162         if (showerLibrary_.get()) {
0163           std::vector<HFShowerLibrary::Hit> hitSL = showerLibrary_->getHits(aStep, isKilled, weight, onlyLong_);
0164           for (unsigned int i = 0; i < hitSL.size(); i++) {
0165             bool ok = true;
0166 #ifdef EDM_ML_DEBUG
0167             edm::LogVerbatim("HFShower") << "HFShowerParam: getHits applyFidCut = " << applyFidCut_;
0168 #endif
0169             if (applyFidCut_) {  // @@ For showerlibrary no z-cut for Short (no z)
0170               int npmt = HFFibreFiducial::PMTNumber(hitSL[i].position);
0171               if (npmt <= 0)
0172                 ok = false;
0173             }
0174             if (ok) {
0175               hit.position = hitSL[i].position;
0176               hit.depth = hitSL[i].depth;
0177               hit.time = hitSL[i].time;
0178               hit.edep = 1;
0179               hits.push_back(hit);
0180 #ifdef plotDebug
0181               if (fillHisto_) {
0182                 double zv = std::abs(hit.position.z()) - gpar_[4];
0183                 hzvem_->Fill(zv);
0184                 em_long_sl_->Fill(hit.position.z() / CLHEP::cm);
0185                 double sq = sqrt(pow(hit.position.x() / CLHEP::cm, 2) + pow(hit.position.y() / CLHEP::cm, 2));
0186                 double zp = hit.position.z() / CLHEP::cm;
0187                 if (hit.depth == 1) {
0188                   em_2d_1_->Fill(zp, sq);
0189                   em_lateral_1_->Fill(sq);
0190                   em_long_1_->Fill(zp);
0191                 } else if (hit.depth == 2) {
0192                   em_2d_2_->Fill(zp, sq);
0193                   em_lateral_2_->Fill(sq);
0194                   em_long_2_->Fill(zp);
0195                 }
0196               }
0197 #endif
0198 #ifdef EDM_ML_DEBUG
0199               edm::LogVerbatim("HFShower")
0200                   << "HFShowerParam: Hit at depth " << hit.depth << " with edep " << hit.edep << " Time " << hit.time;
0201 #endif
0202             }
0203           }
0204         } else {  // GFlash clusters with known z
0205           std::vector<HFGflash::Hit> hitSL = gflash_->gfParameterization(aStep, onlyLong_);
0206           for (unsigned int i = 0; i < hitSL.size(); ++i) {
0207             bool ok = true;
0208             G4ThreeVector pe_effect(hitSL[i].position.x(), hitSL[i].position.y(), hitSL[i].position.z());
0209             double zv = std::abs(pe_effect.z()) - gpar_[4];
0210             //depth
0211             int depth = 1;
0212             int npmt = 0;
0213             if (zv < 0. || zv > gpar_[1]) {
0214 #ifdef mkdebug
0215               edm::LogVerbatim("HFShower") << "-#Zcut-HFShowerParam::getHits:z=" << zv << ",m=" << gpar_[1];
0216 #endif
0217               ok = false;
0218             }
0219             if (ok && applyFidCut_) {
0220               npmt = HFFibreFiducial::PMTNumber(pe_effect);
0221 #ifdef EDM_ML_DEBUG
0222               edm::LogVerbatim("HFShower") << "HFShowerParam::getHits:#PMT= " << npmt << ",z = " << zv;
0223 #endif
0224               if (npmt <= 0) {
0225 #ifdef EDM_ML_DEBUG
0226                 edm::LogVerbatim("HFShower") << "-#PMT=0 cut-HFShowerParam::getHits: npmt = " << npmt;
0227 #endif
0228                 ok = false;
0229               } else if (npmt > 24) {  // a short fibre
0230                 if (zv > gpar_[0]) {
0231                   depth = 2;
0232                 } else {
0233 #ifdef EDM_ML_DEBUG
0234                   edm::LogVerbatim("HFShower") << "-SHORT cut-HFShowerParam::getHits:zMin=" << gpar_[0];
0235 #endif
0236                   ok = false;
0237                 }
0238               }
0239 #ifdef EDM_ML_DEBUG
0240               edm::LogVerbatim("HFShower")
0241                   << "HFShowerParam: npmt " << npmt << " zv " << std::abs(pe_effect.z()) << ":" << gpar_[4] << ":" << zv
0242                   << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0243 #endif
0244             } else {
0245               if (G4UniformRand() > 0.5)
0246                 depth = 2;
0247               if (depth == 2 && zv < gpar_[0])
0248                 ok = false;
0249             }
0250             //attenuation
0251             double dist = fibre_.zShift(localPoint, depth, 0);  // distance to PMT
0252             double r1 = G4UniformRand();
0253 #ifdef EDM_ML_DEBUG
0254             edm::LogVerbatim("HFShower") << "HFShowerParam:Distance to PMT (" << npmt << ") " << dist
0255                                          << ", exclusion flag " << (r1 > exp(-attLMeanInv_ * zv));
0256 #endif
0257             if (r1 > exp(-attLMeanInv_ * dist))
0258               ok = false;
0259             if (ok) {
0260               double r2 = G4UniformRand();
0261 #ifdef EDM_ML_DEBUG
0262               edm::LogVerbatim("HFShower")
0263                   << "HFShowerParam:Extra exclusion " << r2 << ">" << weight << " " << (r2 > weight);
0264 #endif
0265               if (r2 < weight) {
0266                 double time = (equalizeTimeShift_) ? (fibre_.tShift(localPoint, depth, -1))
0267                                                    : (fibre_.tShift(localPoint, depth, 0));
0268 
0269                 hit.position = hitSL[i].position;
0270                 hit.depth = depth;
0271                 hit.time = time + hitSL[i].time;
0272                 hit.edep = 1;
0273                 hits.push_back(hit);
0274 #ifdef plotDebug
0275                 if (fillHisto_) {
0276                   em_long_gflash_->Fill(pe_effect.z() / CLHEP::cm, hitSL[i].edep);
0277                   hzvem_->Fill(zv);
0278                   double sq = sqrt(pow(hit.position.x() / CLHEP::cm, 2) + pow(hit.position.y() / CLHEP::cm, 2));
0279                   double zp = hit.position.z() / CLHEP::cm;
0280                   if (hit.depth == 1) {
0281                     em_2d_1_->Fill(zp, sq);
0282                     em_lateral_1_->Fill(s);
0283                     em_long_1_->Fill(zp);
0284                   } else if (hit.depth == 2) {
0285                     em_2d_2_->Fill(zp, sq);
0286                     em_lateral_2_->Fill(sq);
0287                     em_long_2_->Fill(zp);
0288                   }
0289                 }
0290 #endif
0291 #ifdef EDM_ML_DEBUG
0292                 edm::LogVerbatim("HFShower")
0293                     << "HFShowerParam: Hit at depth " << hit.depth << " with edep " << hit.edep << " Time " << hit.time;
0294 #endif
0295               }
0296             }
0297           }
0298         }
0299       } else {
0300         path = "Rest";
0301         edep *= pePerGeV_;
0302         double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0303         double time = (equalizeTimeShift_) ? (fibre_.tShift(localPoint, 1, -1))
0304                                            : (fibre_.tShift(localPoint, 1, 0));  // remaining part
0305         bool ok = true;
0306         if (applyFidCut_) {  // @@ For showerlibrary no z-cut for Short (no z)
0307           int npmt = HFFibreFiducial::PMTNumber(hitPoint);
0308           if (npmt <= 0)
0309             ok = false;
0310         }
0311 #ifdef EDM_ML_DEBUG
0312         edm::LogVerbatim("HFShower") << "HFShowerParam: getHits hitPoint " << hitPoint << " flag " << ok;
0313 #endif
0314         if (ok) {
0315           hit.depth = 1;
0316           hit.time = tSlice + time;
0317           hit.edep = edep;
0318           hits.push_back(hit);
0319 #ifdef EDM_ML_DEBUG
0320           edm::LogVerbatim("HFShower") << "HFShowerParam: Hit at depth 1 with edep " << edep << " Time " << tSlice
0321                                        << ":" << time << ":" << hit.time;
0322 #endif
0323 #ifdef plotDebug
0324           double zv = std::abs(hitPoint.z()) - gpar_[4];
0325           if (fillHisto_) {
0326             hzvhad_->Fill(zv);
0327           }
0328 #endif
0329           if (zz >= gpar_[0]) {
0330             time = (equalizeTimeShift_) ? (fibre_.tShift(localPoint, 2, -1)) : (fibre_.tShift(localPoint, 2, 0));
0331             hit.depth = 2;
0332             hit.time = tSlice + time;
0333             hits.push_back(hit);
0334 #ifdef EDM_ML_DEBUG
0335             edm::LogVerbatim("HFShower") << "HFShowerParam: Hit at depth 2 with edep " << edep << " Time " << tSlice
0336                                          << ":" << time << hit.time;
0337 #endif
0338 #ifdef plotDebug
0339             if (fillHisto_) {
0340               hzvhad_->Fill(zv);
0341             }
0342 #endif
0343           }
0344         }
0345       }
0346 #ifdef EDM_ML_DEBUG
0347       for (unsigned int ii = 0; ii < hits.size(); ++ii) {
0348         double zv = std::abs(hits[ii].position.z());
0349         if (zv > 12790)
0350           edm::LogVerbatim("HFShower") << "HFShowerParam: Abnormal hit along " << path << " in "
0351                                        << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName() << " at "
0352                                        << hits[ii].position << " zz " << zv << " Edep " << edep << " due to "
0353                                        << track->GetDefinition()->GetParticleName() << " time " << hit.time;
0354       }
0355       edm::LogVerbatim("HFShower") << "HFShowerParam: getHits kill (" << isKilled << ") track " << track->GetTrackID()
0356                                    << " at " << hitPoint << " and deposit " << edep << " " << hits.size() << " times"
0357                                    << " ZZ " << zz << " " << gpar_[0];
0358 #endif
0359     }
0360   }
0361   return hits;
0362 }
[ 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!