HGCalTestBeam/plugins/HGCalTBMB.cc

0001 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ServiceRegistry/interface/Service.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0008 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0009 #include "SimG4Core/Notification/interface/Observer.h"
0010 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0011
0012 #include "G4LogicalVolumeStore.hh"
0013 #include "G4Step.hh"
0014 #include "G4Track.hh"
0015
0016 #include <TH1F.h>
0017
0018 #include <iostream>
0019 #include <string>
0020 #include <vector>
0021
0022 //#define EDM_ML_DEBUG
0023
0024 class HGCalTBMB : public SimWatcher,
0025                   public Observer<const BeginOfTrack*>,
0026                   public Observer<const G4Step*>,
0027                   public Observer<const EndOfTrack*> {
0028 public:
0029   HGCalTBMB(const edm::ParameterSet&);
0030   HGCalTBMB(const HGCalTBMB&) = delete;                   // stop default
0031   const HGCalTBMB& operator=(const HGCalTBMB&) = delete;  // ...
0032   ~HGCalTBMB() override;
0033
0034 private:
0035   void update(const BeginOfTrack*) override;
0036   void update(const G4Step*) override;
0037   void update(const EndOfTrack*) override;
0038
0039   bool stopAfter(const G4Step*);
0040   int findVolume(const G4VTouchable* touch, bool stop) const;
0041
0042   std::vector<std::string> listNames_;
0043   std::string stopName_;
0044   double stopZ_;
0045   unsigned int nList_;
0046   std::vector<double> radLen_, intLen_, stepLen_;
0047   std::vector<TH1D*> me100_, me200_, me300_;
0048 };
0049
0050 HGCalTBMB::HGCalTBMB(const edm::ParameterSet& p) {
0051   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("HGCalTBMB");
0052   listNames_ = m_p.getParameter<std::vector<std::string> >("DetectorNames");
0053   stopName_ = m_p.getParameter<std::string>("StopName");
0054   stopZ_ = m_p.getParameter<double>("MaximumZ");
0055   nList_ = listNames_.size();
0056   edm::LogVerbatim("HGCSim") << "HGCalTBMB initialized for " << nList_ << " volumes";
0057   for (unsigned int k = 0; k < nList_; ++k)
0058     edm::LogVerbatim("HGCSim") << " [" << k << "] " << listNames_[k];
0059   edm::LogVerbatim("HGCSim") << "Stop after " << stopZ_ << " or reaching volume " << stopName_;
0060
0061   edm::Service<TFileService> tfile;
0062   if (!tfile.isAvailable())
0063     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0064                                       << "please add it to config file";
0065   char name[20], title[80];
0066   TH1D* hist;
0067   for (unsigned int i = 0; i <= nList_; i++) {
0068     std::string named = (i == nList_) ? "Total" : listNames_[i];
0069     sprintf(name, "RadL%d", i);
0070     sprintf(title, "MB(X0) for (%s)", named.c_str());
0071     hist = tfile->make<TH1D>(name, title, 100000, 0.0, 100.0);
0072     hist->Sumw2(true);
0073     me100_.push_back(hist);
0074     sprintf(name, "IntL%d", i);
0075     sprintf(title, "MB(L0) for (%s)", named.c_str());
0076     hist = tfile->make<TH1D>(name, title, 100000, 0.0, 10.0);
0077     hist->Sumw2(true);
0078     me200_.push_back(hist);
0079     sprintf(name, "StepL%d", i);
0080     sprintf(title, "MB(Step) for (%s)", named.c_str());
0081     hist = tfile->make<TH1D>(name, title, 100000, 0.0, 50000.0);
0082     hist->Sumw2(true);
0083     me300_.push_back(hist);
0084   }
0085   edm::LogVerbatim("HGCSim") << "HGCalTBMB: Booking user histos done ===";
0086 }
0087
0088 HGCalTBMB::~HGCalTBMB() {}
0089
0090 void HGCalTBMB::update(const BeginOfTrack* trk) {
0091   radLen_ = std::vector<double>(nList_ + 1, 0);
0092   intLen_ = std::vector<double>(nList_ + 1, 0);
0093   stepLen_ = std::vector<double>(nList_ + 1, 0);
0094
0095 #ifdef EDM_ML_DEBUG
0096   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0097   const G4ThreeVector& mom = aTrack->GetMomentum();
0098   double theEnergy = aTrack->GetTotalEnergy();
0099   int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0100   edm::LogVerbatim("HGCSim") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code " << theID
0101                              << " Energy " << theEnergy / CLHEP::GeV << " GeV; Momentum " << mom;
0102 #endif
0103 }
0104
0105 void HGCalTBMB::update(const G4Step* aStep) {
0106   G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0107   double step = aStep->GetStepLength();
0108   double radl = material->GetRadlen();
0109   double intl = material->GetNuclearInterLength();
0110
0111   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0112   int indx = findVolume(touch, false);
0113
0114   if (indx >= 0) {
0115     stepLen_[indx] += step;
0116     radLen_[indx] += (step / radl);
0117     intLen_[indx] += (step / intl);
0118   }
0119   stepLen_[nList_] += step;
0120   radLen_[nList_] += (step / radl);
0121   intLen_[nList_] += (step / intl);
0122 #ifdef EDM_ML_DEBUG
0123   edm::LogVerbatim("HGCSim") << "HGCalTBMB::Step in " << touch->GetVolume(0)->GetLogicalVolume()->GetName() << " Index "
0124                              << indx << " Step " << step << " RadL " << step / radl << " IntL " << step / intl;
0125 #endif
0126
0127   if (stopAfter(aStep)) {
0128     G4Track* track = aStep->GetTrack();
0129     track->SetTrackStatus(fStopAndKill);
0130   }
0131 }
0132
0133 void HGCalTBMB::update(const EndOfTrack* trk) {
0134   for (unsigned int ii = 0; ii <= nList_; ++ii) {
0135     me100_[ii]->Fill(radLen_[ii]);
0136     me200_[ii]->Fill(intLen_[ii]);
0137     me300_[ii]->Fill(stepLen_[ii]);
0138 #ifdef EDM_ML_DEBUG
0139     std::string name("Total");
0140     if (ii < nList_)
0141       name = listNames_[ii];
0142     edm::LogVerbatim("HGCSim") << "HGCalTBMB::Volume[" << ii << "]: " << name << " == Step " << stepLen_[ii] << " RadL "
0143                                << radLen_[ii] << " IntL " << intLen_[ii];
0144 #endif
0145   }
0146 }
0147
0148 bool HGCalTBMB::stopAfter(const G4Step* aStep) {
0149   bool flag(false);
0150   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0151   G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0152   if (aStep->GetPostStepPoint() != nullptr)
0153     hitPoint = aStep->GetPostStepPoint()->GetPosition();
0154   double zz = hitPoint.z();
0155
0156   if ((findVolume(touch, true) == 0) || (zz > stopZ_))
0157     flag = true;
0158 #ifdef EDM_ML_DEBUG
0159   edm::LogVerbatim("HGCSim") << " HGCalTBMB::Name " << touch->GetVolume(0)->GetName() << " z " << zz << " Flag" << flag;
0160 #endif
0161   return flag;
0162 }
0163
0164 int HGCalTBMB::findVolume(const G4VTouchable* touch, bool stop) const {
0165   int ivol = -1;
0166   int level = (touch->GetHistoryDepth()) + 1;
0167   for (int ii = 0; ii < level; ii++) {
0168     std::string name = touch->GetVolume(ii)->GetName();
0169     if (stop) {
0170       if (strcmp(name.c_str(), stopName_.c_str()) == 0)
0171         ivol = 0;
0172     } else {
0173       for (unsigned int k = 0; k < nList_; ++k) {
0174         if (strcmp(name.c_str(), listNames_[k].c_str()) == 0) {
0175           ivol = k;
0176           break;
0177         }
0178       }
0179     }
0180     if (ivol >= 0)
0181       break;
0182   }
0183   return ivol;
0184 }
0185
0186 #include "FWCore/Framework/interface/MakerMacros.h"
0187 #include "FWCore/PluginManager/interface/ModuleDef.h"
0188 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0189
0190 DEFINE_SIMWATCHER(HGCalTBMB);