GFlash/plugins/GflashG4Watcher.cc

0001 #include "G4Electron.hh"
0002 #include "G4Positron.hh"
0003 #include "G4PrimaryParticle.hh"
0004 #include "G4PrimaryVertex.hh"
0005 #include "G4SDManager.hh"
0006 #include "G4Step.hh"
0007 #include "G4VProcess.hh"
0008
0009 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0010 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0011 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0012
0013 #include "SimG4Core/GFlash/interface/GflashG4Watcher.h"
0014
0015 #include <TVector2.h>
0016
0017 using namespace CLHEP;
0018
0019 // constructors and destructor
0020 //
0021 GflashG4Watcher::GflashG4Watcher(const edm::ParameterSet &p) {
0022   edm::ParameterSet myP = p.getParameter<edm::ParameterSet>("GflashG4Watcher");
0023   histFileName_ = myP.getParameter<std::string>("histFileName");
0024   recoEnergyScaleEB_ = myP.getParameter<double>("recoEnergyScaleEB");
0025   recoEnergyScaleEE_ = myP.getParameter<double>("recoEnergyScaleEE");
0026
0027   histFile_ = new TFile(histFileName_.c_str(), "RECREATE");
0028
0029   TH1::AddDirectory(kTRUE);
0030
0031   em_incE = new TH1F("em_incE", "Incoming energy at Ecal;E (GeV);Number of Events", 500, 0.0, 500.0);
0032   em_vtx_rho = new TH1F("em_vtx_rho", "vertex position;#rho (cm);Number of Events", 100, 0.0, 10.0);
0033   em_vtx_z = new TH1F("em_vtx_z", "vertex position;z (cm);Number of Events", 100, -10.0, 10.0);
0034
0035   eb_ssp_rho = new TH1F("eb_ssp_rho", "Shower starting position;#rho (cm);Number of Events", 200, 0.0, 200.0);
0036   eb_hit_long = new TH1F("eb_hit_long",
0037                          "longitudinal hit position;shower depth (cm);Number "
0038                          "of energy weighted hits",
0039                          400,
0040                          0.0,
0041                          200.0);
0042   eb_hit_lat = new TH1F("eb_hit_lat", "lateral hit position;arm (cm);Number of energy weighted hits", 100, 0.0, 5.0);
0043   eb_hit_rz = new TH2F(
0044       "eb_hit_rz", "hit position along the shower direction;shower depth (cm);arm (cm)", 400, 0.0, 200.0, 100, 0.0, 5.0);
0045   eb_hit_long_sd = new TH1F("eb_hit_long_sd",
0046                             "longitudinal hit position in Sensitive Detector;shower depth "
0047                             "(cm);Number of energy weighted hits",
0048                             400,
0049                             0.0,
0050                             200.0);
0051   eb_hit_lat_sd = new TH1F("eb_hit_lat_sd",
0052                            "lateral hit position in Sensitive Detector;arm "
0053                            "(cm);Number of energy weighted hits",
0054                            100,
0055                            0.0,
0056                            5.0);
0057   eb_hit_rz_sd = new TH2F("eb_hit_rz_sd",
0058                           "hit position along the shower direction in "
0059                           "Sensitive Detector;shower depth (cm);arm (cm)",
0060                           400,
0061                           0.0,
0062                           200.0,
0063                           100,
0064                           0.0,
0065                           5.0);
0066
0067   ee_ssp_z = new TH1F("ee_ssp_z", "Shower starting position;z (cm);Number of Events", 800, -400.0, 400.0);
0068   ee_hit_long = new TH1F("ee_hit_long",
0069                          "longitudinal hit position;shower depth (cm);Number "
0070                          "of energy weighted hits",
0071                          800,
0072                          0.0,
0073                          400.0);
0074   ee_hit_lat = new TH1F("ee_hit_lat", "lateral hit position;arm (cm);Number of energy weighted hits", 100, 0.0, 5.0);
0075   ee_hit_rz = new TH2F(
0076       "ee_hit_rz", "hit position along the shower direction;shower depth (cm);arm (cm)", 800, 0.0, 400.0, 100, 0.0, 5.0);
0077   ee_hit_long_sd = new TH1F("ee_hit_long_sd",
0078                             "longitudinal hit position in Sensitive Detector;shower depth "
0079                             "(cm);Number of energy weighted hits",
0080                             800,
0081                             0.0,
0082                             400.0);
0083   ee_hit_lat_sd = new TH1F("ee_hit_lat_sd",
0084                            "lateral hit position in Sensitive Detector;arm "
0085                            "(cm);Number of energy weighted hits",
0086                            100,
0087                            0.0,
0088                            5.0);
0089   ee_hit_rz_sd = new TH2F("ee_hit_rz_sd",
0090                           "hit position along the shower direction in "
0091                           "Sensitive Detector;shower depth (cm);arm (cm)",
0092                           800,
0093                           0.0,
0094                           400.0,
0095                           100,
0096                           0.0,
0097                           5.0);
0098 }
0099
0100 GflashG4Watcher::~GflashG4Watcher() {
0101   histFile_->cd();
0102   histFile_->Write();
0103   histFile_->Close();
0104 }
0105
0106 void GflashG4Watcher::update(const BeginOfEvent *g4Event) {
0107   inc_flag = false;
0108
0109   const G4Event *evt = (*g4Event)();
0110   inc_vertex = evt->GetPrimaryVertex(0)->GetPosition();
0111   inc_position = inc_vertex;
0112   inc_direction = evt->GetPrimaryVertex(0)->GetPrimary(0)->GetMomentum().unit();
0113   inc_energy = evt->GetPrimaryVertex(0)->GetPrimary(0)->GetMomentum().mag();
0114   out_energy = 0;
0115
0116   em_incE->Fill(inc_energy / GeV);
0117   em_vtx_rho->Fill(inc_vertex.rho() / cm);
0118   em_vtx_z->Fill(inc_vertex.z() / cm);
0119
0120   if (std::abs(inc_direction.eta()) < 1.5)
0121     eb_ssp_rho->Fill(inc_position.rho() / cm);
0122   else
0123     ee_ssp_z->Fill(inc_position.z() / cm);
0124 }
0125
0126 void GflashG4Watcher::update(const EndOfEvent *g4Event) {}
0127
0128 void GflashG4Watcher::update(const G4Step *aStep) {
0129   if (aStep == nullptr)
0130     return;
0131
0132   double hitEnergy = aStep->GetTotalEnergyDeposit();
0133
0134   if (hitEnergy < 1.0e-6)
0135     return;
0136
0137   bool inEB = std::abs(inc_direction.eta()) < 1.5;
0138
0139   out_energy += hitEnergy;  // to check outgoing energy
0140
0141   // This is to calculate shower depth and arm of hits from the shower direction
0142   G4ThreeVector hitPosition = aStep->GetPreStepPoint()->GetPosition();
0143   G4ThreeVector diff = hitPosition - inc_position;
0144   double angle = diff.angle(inc_direction);
0145   double diff_z = std::abs(diff.mag() * std::cos(angle));
0146   double diff_r = std::abs(diff.mag() * std::sin(angle));
0147
0148   G4VSensitiveDetector *aSensitive = aStep->GetPreStepPoint()->GetSensitiveDetector();
0149
0150   if (inEB) {  // showers in barrel crystals
0151     hitEnergy *= recoEnergyScaleEB_;
0152     eb_hit_long->Fill(diff_z / cm, hitEnergy / GeV);
0153     eb_hit_lat->Fill(diff_r / cm, hitEnergy / GeV);
0154     eb_hit_rz->Fill(diff_z / cm, diff_r / cm, hitEnergy / GeV);
0155     if (aSensitive) {
0156       eb_hit_long_sd->Fill(diff_z / cm, hitEnergy / GeV);
0157       eb_hit_lat_sd->Fill(diff_r / cm, hitEnergy / GeV);
0158       eb_hit_rz_sd->Fill(diff_z / cm, diff_r / cm, hitEnergy / GeV);
0159     }
0160   } else {  // showers in endcap crystals
0161     hitEnergy *= recoEnergyScaleEE_;
0162     ee_hit_long->Fill(diff_z / cm, hitEnergy / GeV);
0163     ee_hit_lat->Fill(diff_r / cm, hitEnergy / GeV);
0164     ee_hit_rz->Fill(diff_z / cm, diff_r / cm, hitEnergy / GeV);
0165     if (aSensitive) {
0166       ee_hit_long_sd->Fill(diff_z / cm, hitEnergy / GeV);
0167       ee_hit_lat_sd->Fill(diff_r / cm, hitEnergy / GeV);
0168       ee_hit_rz_sd->Fill(diff_z / cm, diff_r / cm, hitEnergy / GeV);
0169     }
0170   }
0171 }
0172
0173 // define this as a plug-in
0174 #include "FWCore/PluginManager/interface/ModuleDef.h"
0175 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0176
0177 DEFINE_SIMWATCHER(GflashG4Watcher);