File indexing completed on 2023-03-17 11:24:56
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
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;
0140
0141
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) {
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 {
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
0174 #include "FWCore/PluginManager/interface/ModuleDef.h"
0175 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0176
0177 DEFINE_SIMWATCHER(GflashG4Watcher);