Forward/plugins/DoCastorAnalysis.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     Forward
0004 // Class  :     DoCastorAnalysis
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author: P. Katsas
0010 //         Created: 02/2007
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 //#define EDM_ML_DEBUG
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   // observer classes
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   //destructor of DoCastorAnalysis
0150
0151   CastorOutputEventFile->cd();
0152   //-- CastorOutputEventFile->Write();
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 //=================================================================== per EVENT
0173
0174 void DoCastorAnalysis::update(const BeginOfJob *job) { edm::LogVerbatim("ForwardSim") << " Starting new job "; }
0175
0176 //==================================================================== per RUN
0177
0178 void DoCastorAnalysis::update(const BeginOfRun *run) {
0179   edm::LogVerbatim("ForwardSim") << std::endl << "DoCastorAnalysis: Starting Run";
0180
0181   // edm::LogVerbatim("ForwardSim") << "DoCastorAnalysis: output event root file created";
0182   // TString treefilename = TreeFileName;
0183   // CastorOutputEventFile = new TFile(treefilename,"RECREATE");
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 //================= End of EVENT ===============
0194
0195 void DoCastorAnalysis::update(const EndOfEvent *evt) {
0196   // Look for the Hit Collection
0197
0198   // access to the G4 hit collections
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   // std::map<int,float,std::less<int> > themap;
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       //themap[volumeID] += aHit->getEnergyDeposit();
0255       int zside, sector, zmodule;
0256
0257       theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
0258
0259       double energy = aHit->getEnergyDeposit() / CLHEP::GeV;
0260       //double time     = aHit->getTimeSlice();
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     //if (debug) edm::LogVerbatim("ForwardSim") << " total energy = " << simhit_etot;
0288 #ifdef EDM_ML_DEBUG
0289     getchar();
0290 #endif
0291     CastorTree->Fill();
0292
0293   }  // nentries > 0
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);