File indexing completed on 2024-05-10 02:21:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include "DataFormats/Math/interface/Point3D.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0019 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0020 #include "SimG4CMS/Forward/interface/CastorNumberingScheme.h"
0021
0022 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0023 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0024 #include "SimG4Core/Notification/interface/EndOfRun.h"
0025 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0026 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0027 #include "SimG4Core/Notification/interface/Observer.h"
0028 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0029
0030 #include "G4SDManager.hh"
0031 #include "G4Step.hh"
0032 #include "G4Track.hh"
0033 #include "G4Event.hh"
0034 #include "G4PrimaryVertex.hh"
0035 #include "G4VProcess.hh"
0036 #include "G4HCofThisEvent.hh"
0037 #include "G4UserEventAction.hh"
0038
0039 #include <CLHEP/Units/SystemOfUnits.h>
0040 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0041 #include <CLHEP/Random/Randomize.h>
0042
0043 #include "TROOT.h"
0044 #include "TFile.h"
0045 #include "TH1.h"
0046 #include "TH2.h"
0047 #include "TProfile.h"
0048 #include "TNtuple.h"
0049 #include "TRandom.h"
0050 #include "TLorentzVector.h"
0051 #include "TUnixSystem.h"
0052 #include "TSystem.h"
0053 #include "TMath.h"
0054 #include "TF1.h"
0055
0056 #include <cassert>
0057 #include <cmath>
0058 #include <iostream>
0059 #include <iomanip>
0060 #include <map>
0061 #include <memory>
0062 #include <string>
0063 #include <vector>
0064
0065
0066
0067 class DoCastorAnalysis : public SimWatcher,
0068 public Observer<const BeginOfJob *>,
0069 public Observer<const BeginOfRun *>,
0070 public Observer<const EndOfRun *>,
0071 public Observer<const BeginOfEvent *>,
0072 public Observer<const EndOfEvent *>,
0073 public Observer<const G4Step *> {
0074 public:
0075 DoCastorAnalysis(const edm::ParameterSet &p);
0076 ~DoCastorAnalysis() override;
0077
0078 private:
0079
0080 void update(const BeginOfJob *run) override;
0081 void update(const BeginOfRun *run) override;
0082 void update(const EndOfRun *run) override;
0083 void update(const BeginOfEvent *evt) override;
0084 void update(const EndOfEvent *evt) override;
0085 void update(const G4Step *step) override;
0086
0087 private:
0088 int verbosity;
0089
0090 std::string TreeFileName;
0091
0092 TFile *CastorOutputEventFile;
0093 TTree *CastorTree;
0094
0095 int eventIndex;
0096
0097 std::vector<double> simhit_x, simhit_y, simhit_z;
0098 std::vector<double> simhit_eta, simhit_phi, simhit_energy;
0099 std::vector<int> simhit_sector, simhit_module;
0100
0101 std::vector<double> *psimhit_x, *psimhit_y, *psimhit_z;
0102 std::vector<double> *psimhit_eta, *psimhit_phi, *psimhit_energy;
0103 std::vector<int> *psimhit_sector, *psimhit_module;
0104
0105 double simhit_etot;
0106 };
0107
0108 DoCastorAnalysis::DoCastorAnalysis(const edm::ParameterSet &p) {
0109 edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("DoCastorAnalysis");
0110 verbosity = m_Anal.getParameter<int>("Verbosity");
0111
0112 TreeFileName = m_Anal.getParameter<std::string>("CastorTreeFileName");
0113
0114 if (verbosity > 0) {
0115 edm::LogVerbatim("ForwardSim") << std::endl;
0116 edm::LogVerbatim("ForwardSim") << "============================================================================";
0117 edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis:: Initialized as observer";
0118
0119 edm::LogVerbatim("ForwardSim") << " Castor Tree will be created";
0120 edm::LogVerbatim("ForwardSim") << " Castor Tree will be in file: " << TreeFileName;
0121 #ifdef EDM_ML_DEBUG
0122 getchar();
0123 #endif
0124 edm::LogVerbatim("ForwardSim") << "============================================================================";
0125 edm::LogVerbatim("ForwardSim") << std::endl;
0126 }
0127
0128 edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis: output event root file created";
0129 TString treefilename = TreeFileName;
0130 CastorOutputEventFile = new TFile(treefilename, "RECREATE");
0131
0132 CastorTree = new TTree("Sim", "Sim");
0133
0134 CastorTree->Branch("simhit_x", "std::vector<double>", &psimhit_x);
0135 CastorTree->Branch("simhit_y", "std::vector<double>", &psimhit_y);
0136 CastorTree->Branch("simhit_z", "std::vector<double>", &psimhit_z);
0137
0138 CastorTree->Branch("simhit_eta", "std::vector<double>", &psimhit_eta);
0139 CastorTree->Branch("simhit_phi", "std::vector<double>", &psimhit_phi);
0140 CastorTree->Branch("simhit_energy", "std::vector<double>", &psimhit_energy);
0141
0142 CastorTree->Branch("simhit_sector", "std::vector<int>", &psimhit_sector);
0143 CastorTree->Branch("simhit_module", "std::vector<int>", &psimhit_module);
0144
0145 CastorTree->Branch("simhit_etot", &simhit_etot, "simhit_etot/D");
0146 }
0147
0148 DoCastorAnalysis::~DoCastorAnalysis() {
0149
0150
0151 CastorOutputEventFile->cd();
0152
0153 CastorTree->Write("", TObject::kOverwrite);
0154 edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis: Ntuple event written";
0155 #ifdef EDM_ML_DEBUG
0156 getchar();
0157 #endif
0158 CastorOutputEventFile->Close();
0159 edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis: Event file closed";
0160 #ifdef EDM_ML_DEBUG
0161 getchar();
0162 #endif
0163
0164 if (verbosity > 0) {
0165 edm::LogVerbatim("ForwardSim") << std::endl << "DoCastorAnalysis: end of process";
0166 #ifdef EDM_ML_DEBUG
0167 getchar();
0168 #endif
0169 }
0170 }
0171
0172
0173
0174 void DoCastorAnalysis::update(const BeginOfJob *job) { edm::LogVerbatim("ForwardSim") << " Starting new job "; }
0175
0176
0177
0178 void DoCastorAnalysis::update(const BeginOfRun *run) {
0179 edm::LogVerbatim("ForwardSim") << std::endl << "DoCastorAnalysis: Starting Run";
0180
0181
0182
0183
0184
0185 eventIndex = 1;
0186 }
0187
0188 void DoCastorAnalysis::update(const BeginOfEvent *evt) {
0189 edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis: Processing Event Number: " << eventIndex;
0190 eventIndex++;
0191 }
0192
0193
0194
0195 void DoCastorAnalysis::update(const EndOfEvent *evt) {
0196
0197
0198
0199 G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
0200
0201 int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
0202 CaloG4HitCollection *theCAFI = (CaloG4HitCollection *)allHC->GetHC(CAFIid);
0203
0204 CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
0205
0206 unsigned int volumeID = 0;
0207
0208
0209 int nentries = theCAFI->entries();
0210 #ifdef EDM_ML_DEBUG
0211 edm::LogVerbatim("ForwardSim") << "nentries in CAFI: " << nentries;
0212 getchar();
0213 #endif
0214
0215 psimhit_x = &simhit_x;
0216 psimhit_x->clear();
0217 psimhit_x->reserve(nentries);
0218
0219 psimhit_y = &simhit_y;
0220 psimhit_y->clear();
0221 psimhit_y->reserve(nentries);
0222
0223 psimhit_z = &simhit_z;
0224 psimhit_z->clear();
0225 psimhit_z->reserve(nentries);
0226
0227 psimhit_eta = &simhit_eta;
0228 psimhit_eta->clear();
0229 psimhit_eta->reserve(nentries);
0230
0231 psimhit_phi = &simhit_phi;
0232 psimhit_phi->clear();
0233 psimhit_phi->reserve(nentries);
0234
0235 psimhit_energy = &simhit_energy;
0236 psimhit_energy->clear();
0237 psimhit_energy->reserve(nentries);
0238
0239 psimhit_sector = &simhit_sector;
0240 psimhit_sector->clear();
0241 psimhit_sector->reserve(nentries);
0242
0243 psimhit_module = &simhit_module;
0244 psimhit_module->clear();
0245 psimhit_module->reserve(nentries);
0246
0247 simhit_etot = 0;
0248
0249 if (nentries > 0) {
0250 for (int ihit = 0; ihit < nentries; ihit++) {
0251 CaloG4Hit *aHit = (*theCAFI)[ihit];
0252 volumeID = aHit->getUnitID();
0253
0254
0255 int zside, sector, zmodule;
0256
0257 theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
0258
0259 double energy = aHit->getEnergyDeposit() / CLHEP::GeV;
0260
0261
0262 math::XYZPoint pos = aHit->getPosition();
0263 double theta = pos.theta();
0264 double eta = -log(tan(theta / 2.));
0265 double phi = pos.phi();
0266
0267 psimhit_x->push_back(pos.x());
0268 psimhit_y->push_back(pos.y());
0269 psimhit_z->push_back(pos.z());
0270
0271 psimhit_eta->push_back(eta);
0272 psimhit_phi->push_back(phi);
0273 psimhit_energy->push_back(energy);
0274
0275 psimhit_sector->push_back(sector);
0276 psimhit_module->push_back(zmodule);
0277
0278 simhit_etot += energy;
0279
0280 #ifdef EDM_ML_DEBUG
0281 edm::LogVerbatim("ForwardSim") << "hit " << ihit + 1 << " : x = " << (*psimhit_x)[ihit]
0282 << " , eta = " << (*psimhit_eta)[ihit] << " , phi = " << (*psimhit_phi)[ihit]
0283 << " , energy = " << (*psimhit_energy)[ihit];
0284 #endif
0285 }
0286
0287
0288 #ifdef EDM_ML_DEBUG
0289 getchar();
0290 #endif
0291 CastorTree->Fill();
0292
0293 }
0294 delete theCastorNumScheme;
0295 }
0296
0297 void DoCastorAnalysis::update(const EndOfRun *run) { ; }
0298
0299 void DoCastorAnalysis::update(const G4Step *aStep) { ; }
0300
0301 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0302 #include "FWCore/PluginManager/interface/ModuleDef.h"
0303
0304 DEFINE_SIMWATCHER(DoCastorAnalysis);