Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-07 02:12:54

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/GlobalPhysicalConstants.h"
0021 #include "CLHEP/Units/GlobalSystemOfUnits.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     : hcalConstants_(hcons), 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_ = hcalConstants_->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, hcalConstants_, 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   fibre_ = std::make_unique<HFFibre>(name, hcalConstants_, hps, p);
0101   attLMeanInv_ = fibre_->attLength(lambdaMean);
0102   edm::LogVerbatim("HFShower") << "att. length used for (lambda=" << lambdaMean
0103                                << ") = " << 1 / (attLMeanInv_ * CLHEP::cm) << " cm";
0104 }
0105 
0106 HFShowerParam::~HFShowerParam() {}
0107 
0108 std::vector<HFShowerParam::Hit> HFShowerParam::getHits(const G4Step* aStep, double weight, bool& isKilled) {
0109   auto const preStepPoint = aStep->GetPreStepPoint();
0110   auto const track = aStep->GetTrack();
0111   bool isEM = G4TrackToParticleID::isGammaElectronPositron(track);
0112   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0113   double zv = std::abs(hitPoint.z()) - gpar_[4] - 0.5 * gpar_[1];
0114   G4ThreeVector localPoint = G4ThreeVector(hitPoint.x(), hitPoint.y(), zv);
0115 
0116   double pin = (preStepPoint->GetTotalEnergy()) / CLHEP::GeV;
0117   double zint = hitPoint.z();
0118   double zz = std::abs(zint) - gpar_[4];
0119 
0120 #ifdef EDM_ML_DEBUG
0121   edm::LogVerbatim("HFShower") << "HFShowerParam: getHits " << track->GetDefinition()->GetParticleName()
0122                                << " of energy " << pin << " GeV Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y()
0123                                << "," << zint << " (" << zz << "," << localPoint.z() << ", "
0124                                << (localPoint.z() + 0.5 * gpar_[1]) << ") Local " << localPoint;
0125 #endif
0126   std::vector<HFShowerParam::Hit> hits;
0127   HFShowerParam::Hit hit;
0128   hit.position = hitPoint;
0129 
0130   // look for other charged particles
0131   bool other = false;
0132   double pBeta = track->GetDynamicParticle()->GetTotalMomentum() / track->GetDynamicParticle()->GetTotalEnergy();
0133   double dirz = (track->GetDynamicParticle()->GetMomentumDirection()).z();
0134   if (hitPoint.z() < 0)
0135     dirz *= -1.;
0136 #ifdef EDM_ML_DEBUG
0137   edm::LogVerbatim("HFShower") << "HFShowerParam: getHits Momentum "
0138                                << track->GetDynamicParticle()->GetMomentumDirection() << " HitPoint " << hitPoint
0139                                << " dirz " << dirz;
0140 #endif
0141   if (!isEM && track->GetDefinition()->GetPDGCharge() != 0 && pBeta > (1 / ref_index_) &&
0142       aStep->GetTotalEnergyDeposit() > 0.) {
0143     other = true;
0144   }
0145 
0146   // take only e+-/gamma/or special particles
0147   if (isEM || other) {
0148     // Leave out the last part
0149     double edep = 0.;
0150     if ((!trackEM_) && ((zz < (gpar_[1] - gpar_[2])) || parametrizeLast_) && (!other)) {
0151       edep = pin;
0152       isKilled = true;
0153     } else if ((track->GetDefinition()->GetPDGCharge() != 0) && (pBeta > (1 / ref_index_)) && (dirz > aperture_)) {
0154       edep = (aStep->GetTotalEnergyDeposit()) / GeV;
0155     }
0156     std::string path = "ShowerLibrary";
0157 #ifdef EDM_ML_DEBUG
0158     edm::LogVerbatim("HFShower") << "HFShowerParam: getHits edep = " << edep << " weight " << weight << " final "
0159                                  << edep * weight << ", Kill = " << isKilled << ", pin = " << pin
0160                                  << ", edMin = " << edMin_ << " Other " << other;
0161 #endif
0162     edep *= weight;
0163     if (edep > 0) {
0164       if ((showerLibrary_.get() || gflash_.get()) && isKilled && pin > edMin_ && (!other)) {
0165         if (showerLibrary_.get()) {
0166           std::vector<HFShowerLibrary::Hit> hitSL = showerLibrary_->getHits(aStep, isKilled, weight, onlyLong_);
0167           for (unsigned int i = 0; i < hitSL.size(); i++) {
0168             bool ok = true;
0169 #ifdef EDM_ML_DEBUG
0170             edm::LogVerbatim("HFShower") << "HFShowerParam: getHits applyFidCut = " << applyFidCut_;
0171 #endif
0172             if (applyFidCut_) {  // @@ For showerlibrary no z-cut for Short (no z)
0173               int npmt = HFFibreFiducial::PMTNumber(hitSL[i].position);
0174               if (npmt <= 0)
0175                 ok = false;
0176             }
0177             if (ok) {
0178               hit.position = hitSL[i].position;
0179               hit.depth = hitSL[i].depth;
0180               hit.time = hitSL[i].time;
0181               hit.edep = 1;
0182               hits.push_back(hit);
0183 #ifdef plotDebug
0184               if (fillHisto_) {
0185                 double zv = std::abs(hit.position.z()) - gpar_[4];
0186                 hzvem_->Fill(zv);
0187                 em_long_sl_->Fill(hit.position.z() / CLHEP::cm);
0188                 double sq = sqrt(pow(hit.position.x() / CLHEP::cm, 2) + pow(hit.position.y() / CLHEP::cm, 2));
0189                 double zp = hit.position.z() / CLHEP::cm;
0190                 if (hit.depth == 1) {
0191                   em_2d_1_->Fill(zp, sq);
0192                   em_lateral_1_->Fill(sq);
0193                   em_long_1_->Fill(zp);
0194                 } else if (hit.depth == 2) {
0195                   em_2d_2_->Fill(zp, sq);
0196                   em_lateral_2_->Fill(sq);
0197                   em_long_2_->Fill(zp);
0198                 }
0199               }
0200 #endif
0201 #ifdef EDM_ML_DEBUG
0202               edm::LogVerbatim("HFShower")
0203                   << "HFShowerParam: Hit at depth " << hit.depth << " with edep " << hit.edep << " Time " << hit.time;
0204 #endif
0205             }
0206           }
0207         } else {  // GFlash clusters with known z
0208           std::vector<HFGflash::Hit> hitSL = gflash_->gfParameterization(aStep, onlyLong_);
0209           for (unsigned int i = 0; i < hitSL.size(); ++i) {
0210             bool ok = true;
0211             G4ThreeVector pe_effect(hitSL[i].position.x(), hitSL[i].position.y(), hitSL[i].position.z());
0212             double zv = std::abs(pe_effect.z()) - gpar_[4];
0213             //depth
0214             int depth = 1;
0215             int npmt = 0;
0216             if (zv < 0. || zv > gpar_[1]) {
0217 #ifdef mkdebug
0218               edm::LogVerbatim("HFShower") << "-#Zcut-HFShowerParam::getHits:z=" << zv << ",m=" << gpar_[1];
0219 #endif
0220               ok = false;
0221             }
0222             if (ok && applyFidCut_) {
0223               npmt = HFFibreFiducial::PMTNumber(pe_effect);
0224 #ifdef EDM_ML_DEBUG
0225               edm::LogVerbatim("HFShower") << "HFShowerParam::getHits:#PMT= " << npmt << ",z = " << zv;
0226 #endif
0227               if (npmt <= 0) {
0228 #ifdef EDM_ML_DEBUG
0229                 edm::LogVerbatim("HFShower") << "-#PMT=0 cut-HFShowerParam::getHits: npmt = " << npmt;
0230 #endif
0231                 ok = false;
0232               } else if (npmt > 24) {  // a short fibre
0233                 if (zv > gpar_[0]) {
0234                   depth = 2;
0235                 } else {
0236 #ifdef EDM_ML_DEBUG
0237                   edm::LogVerbatim("HFShower") << "-SHORT cut-HFShowerParam::getHits:zMin=" << gpar_[0];
0238 #endif
0239                   ok = false;
0240                 }
0241               }
0242 #ifdef EDM_ML_DEBUG
0243               edm::LogVerbatim("HFShower")
0244                   << "HFShowerParam: npmt " << npmt << " zv " << std::abs(pe_effect.z()) << ":" << gpar_[4] << ":" << zv
0245                   << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0246 #endif
0247             } else {
0248               if (G4UniformRand() > 0.5)
0249                 depth = 2;
0250               if (depth == 2 && zv < gpar_[0])
0251                 ok = false;
0252             }
0253             //attenuation
0254             double dist = fibre_->zShift(localPoint, depth, 0);  // distance to PMT
0255             double r1 = G4UniformRand();
0256 #ifdef EDM_ML_DEBUG
0257             edm::LogVerbatim("HFShower") << "HFShowerParam:Distance to PMT (" << npmt << ") " << dist
0258                                          << ", exclusion flag " << (r1 > exp(-attLMeanInv_ * zv));
0259 #endif
0260             if (r1 > exp(-attLMeanInv_ * dist))
0261               ok = false;
0262             if (ok) {
0263               double r2 = G4UniformRand();
0264 #ifdef EDM_ML_DEBUG
0265               edm::LogVerbatim("HFShower")
0266                   << "HFShowerParam:Extra exclusion " << r2 << ">" << weight << " " << (r2 > weight);
0267 #endif
0268               if (r2 < weight) {
0269                 double time = (equalizeTimeShift_) ? (fibre_->tShift(localPoint, depth, -1))
0270                                                    : (fibre_->tShift(localPoint, depth, 0));
0271 
0272                 hit.position = hitSL[i].position;
0273                 hit.depth = depth;
0274                 hit.time = time + hitSL[i].time;
0275                 hit.edep = 1;
0276                 hits.push_back(hit);
0277 #ifdef plotDebug
0278                 if (fillHisto_) {
0279                   em_long_gflash_->Fill(pe_effect.z() / CLHEP::cm, hitSL[i].edep);
0280                   hzvem_->Fill(zv);
0281                   double sq = sqrt(pow(hit.position.x() / CLHEP::cm, 2) + pow(hit.position.y() / CLHEP::cm, 2));
0282                   double zp = hit.position.z() / CLHEP::cm;
0283                   if (hit.depth == 1) {
0284                     em_2d_1_->Fill(zp, sq);
0285                     em_lateral_1_->Fill(s);
0286                     em_long_1_->Fill(zp);
0287                   } else if (hit.depth == 2) {
0288                     em_2d_2_->Fill(zp, sq);
0289                     em_lateral_2_->Fill(sq);
0290                     em_long_2_->Fill(zp);
0291                   }
0292                 }
0293 #endif
0294 #ifdef EDM_ML_DEBUG
0295                 edm::LogVerbatim("HFShower")
0296                     << "HFShowerParam: Hit at depth " << hit.depth << " with edep " << hit.edep << " Time " << hit.time;
0297 #endif
0298               }
0299             }
0300           }
0301         }
0302       } else {
0303         path = "Rest";
0304         edep *= pePerGeV_;
0305         double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0306         double time = (equalizeTimeShift_) ? (fibre_->tShift(localPoint, 1, -1))
0307                                            : (fibre_->tShift(localPoint, 1, 0));  // remaining part
0308         bool ok = true;
0309         if (applyFidCut_) {  // @@ For showerlibrary no z-cut for Short (no z)
0310           int npmt = HFFibreFiducial::PMTNumber(hitPoint);
0311           if (npmt <= 0)
0312             ok = false;
0313         }
0314 #ifdef EDM_ML_DEBUG
0315         edm::LogVerbatim("HFShower") << "HFShowerParam: getHits hitPoint " << hitPoint << " flag " << ok;
0316 #endif
0317         if (ok) {
0318           hit.depth = 1;
0319           hit.time = tSlice + time;
0320           hit.edep = edep;
0321           hits.push_back(hit);
0322 #ifdef EDM_ML_DEBUG
0323           edm::LogVerbatim("HFShower") << "HFShowerParam: Hit at depth 1 with edep " << edep << " Time " << tSlice
0324                                        << ":" << time << ":" << hit.time;
0325 #endif
0326 #ifdef plotDebug
0327           double zv = std::abs(hitPoint.z()) - gpar_[4];
0328           if (fillHisto_) {
0329             hzvhad_->Fill(zv);
0330           }
0331 #endif
0332           if (zz >= gpar_[0]) {
0333             time = (equalizeTimeShift_) ? (fibre_->tShift(localPoint, 2, -1)) : (fibre_->tShift(localPoint, 2, 0));
0334             hit.depth = 2;
0335             hit.time = tSlice + time;
0336             hits.push_back(hit);
0337 #ifdef EDM_ML_DEBUG
0338             edm::LogVerbatim("HFShower") << "HFShowerParam: Hit at depth 2 with edep " << edep << " Time " << tSlice
0339                                          << ":" << time << hit.time;
0340 #endif
0341 #ifdef plotDebug
0342             if (fillHisto_) {
0343               hzvhad_->Fill(zv);
0344             }
0345 #endif
0346           }
0347         }
0348       }
0349 #ifdef EDM_ML_DEBUG
0350       for (unsigned int ii = 0; ii < hits.size(); ++ii) {
0351         double zv = std::abs(hits[ii].position.z());
0352         if (zv > 12790)
0353           edm::LogVerbatim("HFShower") << "HFShowerParam: Abnormal hit along " << path << " in "
0354                                        << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName() << " at "
0355                                        << hits[ii].position << " zz " << zv << " Edep " << edep << " due to "
0356                                        << track->GetDefinition()->GetParticleName() << " time " << hit.time;
0357       }
0358       edm::LogVerbatim("HFShower") << "HFShowerParam: getHits kill (" << isKilled << ") track " << track->GetTrackID()
0359                                    << " at " << hitPoint << " and deposit " << edep << " " << hits.size() << " times"
0360                                    << " ZZ " << zz << " " << gpar_[0];
0361 #endif
0362     }
0363   }
0364   return hits;
0365 }