File indexing completed on 2023-03-17 11:24:40
0001
0002 #include <cmath>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <string>
0006 #include <vector>
0007
0008
0009 #include "DataFormats/Math/interface/Point3D.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0016 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0017 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0018 #include "SimG4Core/Watcher/interface/SimProducer.h"
0019 #include "SimG4Core/Notification/interface/Observer.h"
0020
0021 #include "SimG4CMS/ShowerLibraryProducer/interface/FiberG4Hit.h"
0022 #include "SimG4CMS/ShowerLibraryProducer/interface/HFShowerG4Hit.h"
0023
0024 #include "G4HCofThisEvent.hh"
0025 #include "G4SDManager.hh"
0026 #include "G4Step.hh"
0027 #include "G4Track.hh"
0028 #include "G4ThreeVector.hh"
0029 #include "G4VProcess.hh"
0030
0031 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0032 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0033
0034 #include "TFile.h"
0035 #include "TTree.h"
0036
0037
0038
0039 class HcalForwardAnalysis : public SimProducer,
0040 public Observer<const BeginOfRun*>,
0041 public Observer<const BeginOfEvent*>,
0042 public Observer<const EndOfEvent*>,
0043 public Observer<const G4Step*> {
0044 public:
0045 struct Photon {
0046 Photon(int id, float X, float Y, float Z, float T, float Lambda)
0047 : fiberId(id), x(X), y(Y), z(Z), t(T), lambda(Lambda) {}
0048 int fiberId;
0049 float x;
0050 float y;
0051 float z;
0052 float t;
0053 float lambda;
0054 };
0055
0056 HcalForwardAnalysis(const edm::ParameterSet& p);
0057 HcalForwardAnalysis(const HcalForwardAnalysis&) = delete;
0058 const HcalForwardAnalysis& operator=(const HcalForwardAnalysis&) = delete;
0059 ~HcalForwardAnalysis() override;
0060
0061 void produce(edm::Event&, const edm::EventSetup&) override;
0062
0063 private:
0064 void init();
0065
0066
0067 void update(const BeginOfRun* run) override;
0068 void update(const BeginOfEvent* evt) override;
0069 void update(const G4Step* step) override;
0070 void update(const EndOfEvent* evt) override;
0071
0072
0073
0074 void setPhotons(const EndOfEvent* evt);
0075
0076 void fillEvent();
0077 void parseDetId(int id, int& tower, int& cell, int& fiber);
0078 void clear();
0079
0080 TTree* theTree;
0081 int theEventCounter;
0082 int count;
0083 int evNum;
0084 float x[10000], y[10000], z[10000], t[10000], lambda[10000];
0085 float primX, primY, primZ, primT;
0086 float primMomX, primMomY, primMomZ;
0087 int nphot;
0088 int fiberId[10000];
0089 std::vector<Photon> thePhotons;
0090 std::vector<std::string> theNames;
0091 bool fillt;
0092 };
0093
0094 HcalForwardAnalysis::HcalForwardAnalysis(const edm::ParameterSet& p) {
0095 edm::ParameterSet m_SLP = p.getParameter<edm::ParameterSet>("HFShowerLibraryProducer");
0096 theNames = m_SLP.getParameter<std::vector<std::string> >("Names");
0097
0098
0099 init();
0100 theEventCounter = 0;
0101 nphot = 0;
0102 for (int i = 0; i < 10000; ++i) {
0103 x[i] = 0.;
0104 y[i] = 0.;
0105 z[i] = 0.;
0106 t[i] = 0.;
0107 lambda[i] = 0.;
0108 fiberId[i] = 0;
0109 }
0110 primX = primY = primZ = primT = 0.;
0111 primMomX = primMomY = primMomZ = 0.;
0112 }
0113
0114 HcalForwardAnalysis::~HcalForwardAnalysis() {}
0115
0116
0117
0118
0119
0120 void HcalForwardAnalysis::produce(edm::Event& iEvent, const edm::EventSetup&) {
0121 if (fillt)
0122 fillEvent();
0123 }
0124
0125 void HcalForwardAnalysis::init() {
0126 edm::Service<TFileService> theFile;
0127 theTree = theFile->make<TTree>("CherenkovPhotons", "Cherenkov Photons");
0128 theTree->Branch("nphot", &nphot, "nphot/I");
0129 theTree->Branch("x", &x, "x[nphot]/F");
0130 theTree->Branch("y", &y, "y[nphot]/F");
0131 theTree->Branch("z", &z, "z[nphot]/F");
0132 theTree->Branch("t", &t, "t[nphot]/F");
0133 theTree->Branch("lambda", &lambda, "lambda[nphot]/F");
0134 theTree->Branch("fiberId", &fiberId, "fiberId[nphot]/I");
0135 theTree->Branch("primX", &primX, "primX/F");
0136 theTree->Branch("primY", &primY, "primY/F");
0137 theTree->Branch("primZ", &primZ, "primZ/F");
0138 theTree->Branch("primMomX", &primMomX, "primMomX/F");
0139 theTree->Branch("primMomY", &primMomY, "primMomY/F");
0140 theTree->Branch("primMomZ", &primMomZ, "primMomZ/F");
0141 theTree->Branch("primT", &primT, "primT/F");
0142
0143
0144 count = 0;
0145 evNum = 0;
0146 clear();
0147 }
0148
0149 void HcalForwardAnalysis::update(const BeginOfRun* run) {
0150 int irun = (*run)()->GetRunID();
0151 edm::LogVerbatim("HcalForwardLib") << " =====> Begin of Run = " << irun;
0152 }
0153
0154 void HcalForwardAnalysis::update(const BeginOfEvent* evt) {
0155 evNum = (*evt)()->GetEventID();
0156 clear();
0157 #ifdef EDM_ML_DEBUG
0158 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> Begin of event = " << evNum;
0159 #endif
0160 }
0161
0162 void HcalForwardAnalysis::update(const G4Step* aStep) {}
0163
0164 void HcalForwardAnalysis::update(const EndOfEvent* evt) {
0165 count++;
0166
0167
0168 #ifdef EDM_ML_DEBUG
0169 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::Fill event " << (*evt)()->GetEventID();
0170 #endif
0171 setPhotons(evt);
0172
0173 int iEvt = (*evt)()->GetEventID();
0174 if (iEvt < 10)
0175 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0176 else if ((iEvt < 100) && (iEvt % 10 == 0))
0177 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0178 else if ((iEvt < 1000) && (iEvt % 100 == 0))
0179 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0180 else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0181 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0182 }
0183
0184 void HcalForwardAnalysis::setPhotons(const EndOfEvent* evt) {
0185 fillt = true;
0186 int idHC, j;
0187 FiberG4HitsCollection* theHC;
0188
0189 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0190 #ifdef EDM_ML_DEBUG
0191 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Has " << allHC->GetNumberOfCollections()
0192 << " collections";
0193 for (int k = 0; k < allHC->GetNumberOfCollections(); ++k) {
0194 G4String name = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetName();
0195 G4String nameSD = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetSDname();
0196 edm::LogVerbatim("HcalForwardLib") << "Collecttion[" << k << "] " << allHC->GetHC(k) << " " << name << ":"
0197 << nameSD;
0198 }
0199 #endif
0200 std::string sdName = theNames[0];
0201 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0202 theHC = (FiberG4HitsCollection*)allHC->GetHC(idHC);
0203 #ifdef EDM_ML_DEBUG
0204 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Hit Collection for " << sdName << " of ID "
0205 << idHC << " is obtained at " << theHC;
0206 #endif
0207 std::vector<HFShowerPhoton> ShortFiberPhotons;
0208 std::vector<HFShowerPhoton> LongFiberPhotons;
0209 LongFiberPhotons.clear();
0210 ShortFiberPhotons.clear();
0211 if (idHC >= 0 && theHC != nullptr) {
0212 int thehc_entries = theHC->entries();
0213 #ifdef EDM_ML_DEBUG
0214 edm::LogVerbatim("HcalForwardLib") << "FiberhitSize " << thehc_entries;
0215 #endif
0216 for (j = 0; j < thehc_entries; j++) {
0217 FiberG4Hit* aHit = (*theHC)[j];
0218 std::vector<HFShowerPhoton> thePhotonsFromHit = aHit->photon();
0219 #ifdef EDM_ML_DEBUG
0220 edm::LogVerbatim("HcalForwardLib") << "Fiberhit " << j << " has " << thePhotonsFromHit.size() << " photons.";
0221 #endif
0222 int fTowerId = -1;
0223 int fCellId = -1;
0224 int fFiberId = -1;
0225 parseDetId(aHit->towerId(), fTowerId, fCellId, fFiberId);
0226 for (unsigned int iph = 0; iph < thePhotonsFromHit.size(); ++iph) {
0227 if (aHit->depth() == 1)
0228 LongFiberPhotons.push_back(thePhotonsFromHit[iph]);
0229 if (aHit->depth() == 2)
0230 ShortFiberPhotons.push_back(thePhotonsFromHit[iph]);
0231 }
0232 #ifdef EDM_ML_DEBUG
0233 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() NbPhotons " << thePhotonsFromHit.size()
0234 << " towerId " << fTowerId << " cellId " << fCellId << " fiberId " << fFiberId
0235 << " depth " << aHit->depth();
0236 #endif
0237 }
0238 } else {
0239 fillt = false;
0240 #ifdef EDM_ML_DEBUG
0241 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Photons!";
0242 #endif
0243 return;
0244 }
0245 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
0246 << " ShortFibPhotons: " << ShortFiberPhotons.size();
0247 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
0248 << " ShortFibPhotons: " << ShortFiberPhotons.size();
0249
0250
0251 HFShowerG4HitsCollection* theChamberHC;
0252 G4HCofThisEvent* allChamberHC = (*evt)()->GetHCofThisEvent();
0253 sdName = theNames[1];
0254 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0255 theChamberHC = (HFShowerG4HitsCollection*)allChamberHC->GetHC(idHC);
0256 math::XYZPoint primPosOnSurf(0, 0, 0);
0257 math::XYZPoint primMomDirOnSurf(0, 0, 0);
0258 float primTimeOnSurf = 0;
0259
0260
0261
0262
0263 if (idHC >= 0 && theChamberHC != nullptr) {
0264 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hits size: "
0265 << theChamberHC->entries();
0266 int thec_hc_entries = theChamberHC->entries();
0267 for (j = 0; j < thec_hc_entries; ++j) {
0268 HFShowerG4Hit* aHit = (*theChamberHC)[j];
0269 #ifdef EDM_ML_DEBUG
0270 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hit id " << aHit->hitId()
0271 << " track id " << aHit->trackId() << " prim. pos. " << aHit->globalPosition()
0272 << " prom mom. dir. " << aHit->primaryMomDir() << " time " << aHit->time();
0273 #endif
0274 primPosOnSurf.SetXYZ(aHit->globalPosition().x(), aHit->globalPosition().y(), aHit->globalPosition().z());
0275 primMomDirOnSurf.SetXYZ(aHit->primaryMomDir().x(), aHit->primaryMomDir().y(), aHit->primaryMomDir().z());
0276 primTimeOnSurf = aHit->time();
0277 }
0278 } else {
0279 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Chamber hits are stored";
0280 fillt = false;
0281 return;
0282 }
0283 primX = primPosOnSurf.x();
0284 primY = primPosOnSurf.y();
0285 primZ = primPosOnSurf.z();
0286 if (primZ < 990) {
0287 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): First interaction before HF";
0288 fillt = false;
0289 return;
0290 }
0291 primT = primTimeOnSurf;
0292 primMomX = primMomDirOnSurf.x();
0293 primMomY = primMomDirOnSurf.y();
0294 primMomZ = primMomDirOnSurf.z();
0295
0296 double theta = primMomDirOnSurf.theta();
0297 double phi = primMomDirOnSurf.phi();
0298
0299
0300 double sphi = sin(phi);
0301 double cphi = cos(phi);
0302 double ctheta = cos(theta);
0303 double stheta = sin(theta);
0304
0305 double pex = 0, pey = 0, zv = 0;
0306 double xx, yy, zz;
0307
0308 for (unsigned int k = 0; k < LongFiberPhotons.size(); ++k) {
0309 HFShowerPhoton aPhoton = LongFiberPhotons[k];
0310
0311 xx = aPhoton.x();
0312 yy = aPhoton.y();
0313 zz = aPhoton.z();
0314
0315
0316 pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
0317 pey = -xx * sphi + yy * cphi;
0318 zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
0319
0320 double photonProdTime = aPhoton.t() - primTimeOnSurf;
0321 thePhotons.push_back(Photon(1, pex, pey, zv, photonProdTime, aPhoton.lambda()));
0322 }
0323 for (unsigned int k = 0; k < ShortFiberPhotons.size(); ++k) {
0324 HFShowerPhoton aPhoton = ShortFiberPhotons[k];
0325
0326 xx = aPhoton.x();
0327 yy = aPhoton.y();
0328 zz = aPhoton.z();
0329
0330
0331 pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
0332 pey = -xx * sphi + yy * cphi;
0333 zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
0334
0335 double photonProdTime = aPhoton.t() - primTimeOnSurf;
0336 thePhotons.push_back(Photon(2, pex, pey, zv, photonProdTime, aPhoton.lambda()));
0337 }
0338 }
0339
0340 void HcalForwardAnalysis::fillEvent() {
0341 #ifdef EDM_ML_DEBUG
0342 edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> filledEvent";
0343 #endif
0344 nphot = int(thePhotons.size());
0345 for (int i = 0; i < nphot; ++i) {
0346 x[i] = thePhotons[i].x;
0347 y[i] = thePhotons[i].y;
0348 z[i] = thePhotons[i].z;
0349 t[i] = thePhotons[i].t;
0350 lambda[i] = thePhotons[i].lambda;
0351 fiberId[i] = thePhotons[i].fiberId;
0352 }
0353 theTree->Fill();
0354 }
0355
0356 void HcalForwardAnalysis::parseDetId(int id, int& tower, int& cell, int& fiber) {
0357 tower = id / 10000;
0358 cell = id / 10 - tower * 10;
0359 fiber = id - tower * 10000 - cell * 10;
0360 }
0361
0362 void HcalForwardAnalysis::clear() {
0363 nphot = 0;
0364 for (int i = 0; i < 10000; ++i) {
0365 x[i] = 0.;
0366 y[i] = 0.;
0367 z[i] = 0.;
0368 t[i] = 0.;
0369 lambda[i] = 0.;
0370 fiberId[i] = 0;
0371 }
0372 primX = primY = primZ = primT = 0.;
0373 primMomX = primMomY = primMomZ = 0.;
0374
0375 thePhotons.clear();
0376 }
0377
0378 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0379 #include "FWCore/PluginManager/interface/ModuleDef.h"
0380
0381 DEFINE_SIMWATCHER(HcalForwardAnalysis);