Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-17 23:03:07 

0001 // system include files

0002 #include <cmath>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <string>
0006 #include <vector>
0007 
0008 // user include files

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 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 
0023 #include "SimG4CMS/ShowerLibraryProducer/interface/FiberG4Hit.h"
0024 #include "SimG4CMS/ShowerLibraryProducer/interface/HFShowerG4Hit.h"
0025 
0026 #include "G4HCofThisEvent.hh"
0027 #include "G4SDManager.hh"
0028 #include "G4Step.hh"
0029 #include "G4Track.hh"
0030 #include "G4ThreeVector.hh"
0031 #include "G4VProcess.hh"
0032 
0033 #include <CLHEP/Units/SystemOfUnits.h>
0034 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0035 
0036 #include "TFile.h"
0037 #include "TTree.h"
0038 
0039 //#define EDM_ML_DEBUG

0040 
0041 class HcalForwardAnalysis : public SimProducer,
0042                             public Observer<const BeginOfRun*>,
0043                             public Observer<const BeginOfEvent*>,
0044                             public Observer<const EndOfEvent*>,
0045                             public Observer<const G4Step*> {
0046 public:
0047   struct Photon {
0048     Photon(int id, float X, float Y, float Z, float T, float Lambda)
0049         : fiberId(id), x(X), y(Y), z(Z), t(T), lambda(Lambda) {}
0050     int fiberId;
0051     float x;
0052     float y;
0053     float z;
0054     float t;
0055     float lambda;
0056   };
0057 
0058   HcalForwardAnalysis(const edm::ParameterSet& p);
0059   HcalForwardAnalysis(const HcalForwardAnalysis&) = delete;  // stop default

0060   const HcalForwardAnalysis& operator=(const HcalForwardAnalysis&) = delete;
0061   ~HcalForwardAnalysis() override;
0062 
0063   void produce(edm::Event&, const edm::EventSetup&) override;
0064 
0065 private:
0066   void init();
0067 
0068   // observer methods

0069   void update(const BeginOfRun* run) override;
0070   void update(const BeginOfEvent* evt) override;
0071   void update(const G4Step* step) override;
0072   void update(const EndOfEvent* evt) override;
0073   //  void write(const EndOfRun * run);

0074 
0075   //User methods

0076   void setPhotons(const EndOfEvent* evt);
0077   //void fillEvent(PHcalForwardLibInfo&);

0078   void fillEvent();
0079   void parseDetId(int id, int& tower, int& cell, int& fiber);
0080   void clear();
0081 
0082   TTree* theTree;
0083   int theEventCounter;
0084   int count;
0085   int evNum;
0086   float x[10000], y[10000], z[10000], t[10000], lambda[10000];
0087   float primX, primY, primZ, primT;
0088   float primMomX, primMomY, primMomZ;
0089   int nphot;
0090   int fiberId[10000];
0091   std::vector<Photon> thePhotons;
0092   std::vector<std::string> theNames;
0093   bool fillt;
0094 };
0095 
0096 HcalForwardAnalysis::HcalForwardAnalysis(const edm::ParameterSet& p) {
0097   edm::ParameterSet m_SLP = p.getParameter<edm::ParameterSet>("HFShowerLibraryProducer");
0098   theNames = m_SLP.getParameter<std::vector<std::string> >("Names");
0099   //LibVer = m_HS.getParameter<std::string> ("LibVer");

0100   //produces<HFShowerPhotonCollection> ();

0101   init();
0102   theEventCounter = 0;
0103   nphot = 0;
0104   for (int i = 0; i < 10000; ++i) {
0105     x[i] = 0.;
0106     y[i] = 0.;
0107     z[i] = 0.;
0108     t[i] = 0.;
0109     lambda[i] = 0.;
0110     fiberId[i] = 0;
0111   }
0112   primX = primY = primZ = primT = 0.;
0113   primMomX = primMomY = primMomZ = 0.;
0114 }
0115 
0116 HcalForwardAnalysis::~HcalForwardAnalysis() {}
0117 
0118 //

0119 // member functions

0120 //

0121 
0122 void HcalForwardAnalysis::produce(edm::Event& iEvent, const edm::EventSetup&) {
0123   if (fillt)
0124     fillEvent();
0125 }
0126 
0127 void HcalForwardAnalysis::init() {
0128   edm::Service<TFileService> theFile;
0129   theTree = theFile->make<TTree>("CherenkovPhotons", "Cherenkov Photons");
0130   theTree->Branch("nphot", &nphot, "nphot/I");
0131   theTree->Branch("x", &x, "x[nphot]/F");
0132   theTree->Branch("y", &y, "y[nphot]/F");
0133   theTree->Branch("z", &z, "z[nphot]/F");
0134   theTree->Branch("t", &t, "t[nphot]/F");
0135   theTree->Branch("lambda", &lambda, "lambda[nphot]/F");
0136   theTree->Branch("fiberId", &fiberId, "fiberId[nphot]/I");
0137   theTree->Branch("primX", &primX, "primX/F");
0138   theTree->Branch("primY", &primY, "primY/F");
0139   theTree->Branch("primZ", &primZ, "primZ/F");
0140   theTree->Branch("primMomX", &primMomX, "primMomX/F");
0141   theTree->Branch("primMomY", &primMomY, "primMomY/F");
0142   theTree->Branch("primMomZ", &primMomZ, "primMomZ/F");
0143   theTree->Branch("primT", &primT, "primT/F");
0144 
0145   // counter

0146   count = 0;
0147   evNum = 0;
0148   clear();
0149 }
0150 
0151 void HcalForwardAnalysis::update(const BeginOfRun* run) {
0152   int irun = (*run)()->GetRunID();
0153   edm::LogVerbatim("HcalForwardLib") << " =====> Begin of Run = " << irun;
0154 }
0155 
0156 void HcalForwardAnalysis::update(const BeginOfEvent* evt) {
0157   evNum = (*evt)()->GetEventID();
0158   clear();
0159 #ifdef EDM_ML_DEBUG
0160   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> Begin of event = " << evNum;
0161 #endif
0162 }
0163 
0164 void HcalForwardAnalysis::update(const G4Step* aStep) {}
0165 
0166 void HcalForwardAnalysis::update(const EndOfEvent* evt) {
0167   count++;
0168 
0169   //fill the buffer

0170 #ifdef EDM_ML_DEBUG
0171   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::Fill event " << (*evt)()->GetEventID();
0172 #endif
0173   setPhotons(evt);
0174 
0175   int iEvt = (*evt)()->GetEventID();
0176   if (iEvt < 10)
0177     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0178   else if ((iEvt < 100) && (iEvt % 10 == 0))
0179     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0180   else if ((iEvt < 1000) && (iEvt % 100 == 0))
0181     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0182   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0183     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
0184 }
0185 
0186 void HcalForwardAnalysis::setPhotons(const EndOfEvent* evt) {
0187   fillt = true;
0188   int idHC, j;
0189   FiberG4HitsCollection* theHC;
0190   // Look for the Hit Collection of HCal

0191   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0192 #ifdef EDM_ML_DEBUG
0193   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Has " << allHC->GetNumberOfCollections()
0194                                      << " collections";
0195   for (int k = 0; k < allHC->GetNumberOfCollections(); ++k) {
0196     G4String name = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetName();
0197     G4String nameSD = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetSDname();
0198     edm::LogVerbatim("HcalForwardLib") << "Collecttion[" << k << "] " << allHC->GetHC(k) << "  " << name << ":"
0199                                        << nameSD;
0200   }
0201 #endif
0202   std::string sdName = theNames[0];  //name for fiber hits

0203   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0204   theHC = (FiberG4HitsCollection*)allHC->GetHC(idHC);
0205 #ifdef EDM_ML_DEBUG
0206   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Hit Collection for " << sdName << " of ID "
0207                                      << idHC << " is obtained at " << theHC;
0208 #endif
0209   std::vector<HFShowerPhoton> ShortFiberPhotons;
0210   std::vector<HFShowerPhoton> LongFiberPhotons;
0211   LongFiberPhotons.clear();
0212   ShortFiberPhotons.clear();
0213   if (idHC >= 0 && theHC != nullptr) {
0214     int thehc_entries = theHC->entries();
0215 #ifdef EDM_ML_DEBUG
0216     edm::LogVerbatim("HcalForwardLib") << "FiberhitSize " << thehc_entries;
0217 #endif
0218     for (j = 0; j < thehc_entries; j++) {
0219       FiberG4Hit* aHit = (*theHC)[j];
0220       std::vector<HFShowerPhoton> thePhotonsFromHit = aHit->photon();
0221 #ifdef EDM_ML_DEBUG
0222       edm::LogVerbatim("HcalForwardLib") << "Fiberhit " << j << " has " << thePhotonsFromHit.size() << " photons.";
0223 #endif
0224       int fTowerId = -1;
0225       int fCellId = -1;
0226       int fFiberId = -1;
0227       parseDetId(aHit->towerId(), fTowerId, fCellId, fFiberId);
0228       for (unsigned int iph = 0; iph < thePhotonsFromHit.size(); ++iph) {
0229         if (aHit->depth() == 1)
0230           LongFiberPhotons.push_back(thePhotonsFromHit[iph]);
0231         if (aHit->depth() == 2)
0232           ShortFiberPhotons.push_back(thePhotonsFromHit[iph]);
0233       }
0234 #ifdef EDM_ML_DEBUG
0235       edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() NbPhotons " << thePhotonsFromHit.size()
0236                                          << " towerId " << fTowerId << " cellId " << fCellId << " fiberId " << fFiberId
0237                                          << " depth " << aHit->depth();
0238 #endif
0239     }
0240   } else {
0241     fillt = false;
0242 #ifdef EDM_ML_DEBUG
0243     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Photons!";
0244 #endif
0245     return;
0246   }
0247   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
0248                                      << " ShortFibPhotons: " << ShortFiberPhotons.size();
0249   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
0250                                      << " ShortFibPhotons: " << ShortFiberPhotons.size();
0251 
0252   //Chamber hits to find information about primary particle on surface

0253   HFShowerG4HitsCollection* theChamberHC;
0254   G4HCofThisEvent* allChamberHC = (*evt)()->GetHCofThisEvent();
0255   sdName = theNames[1];
0256   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0257   theChamberHC = (HFShowerG4HitsCollection*)allChamberHC->GetHC(idHC);
0258   math::XYZPoint primPosOnSurf(0, 0, 0);
0259   math::XYZPoint primMomDirOnSurf(0, 0, 0);
0260   float primTimeOnSurf = 0;
0261   //    the chamber hit is for primary particle, but step size can be small

0262   //    (in newer Geant4 versions) and as a result primary particle may have

0263   //    multiple hits. We want to take last one which is close the HF absorber

0264   //  if (idHC >= 0 && theChamberHC != nullptr && theChamberHC->entries()>0) {

0265   if (idHC >= 0 && theChamberHC != nullptr) {
0266     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hits size: "
0267                                        << theChamberHC->entries();
0268     int thec_hc_entries = theChamberHC->entries();
0269     for (j = 0; j < thec_hc_entries; ++j) {
0270       HFShowerG4Hit* aHit = (*theChamberHC)[j];
0271 #ifdef EDM_ML_DEBUG
0272       edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hit id " << aHit->hitId()
0273                                          << " track id " << aHit->trackId() << " prim. pos. " << aHit->globalPosition()
0274                                          << " prom mom. dir. " << aHit->primaryMomDir() << " time " << aHit->time();
0275 #endif
0276       primPosOnSurf.SetXYZ(aHit->globalPosition().x(), aHit->globalPosition().y(), aHit->globalPosition().z());
0277       primMomDirOnSurf.SetXYZ(aHit->primaryMomDir().x(), aHit->primaryMomDir().y(), aHit->primaryMomDir().z());
0278       primTimeOnSurf = aHit->time();
0279     }
0280   } else {
0281     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Chamber hits are stored";
0282     fillt = false;
0283     return;
0284   }
0285   primX = primPosOnSurf.x();
0286   primY = primPosOnSurf.y();
0287   primZ = primPosOnSurf.z();
0288   if (primZ < 990) {  // there were interactions before HF

0289     edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): First interaction before HF";
0290     fillt = false;
0291     return;
0292   }
0293   primT = primTimeOnSurf;
0294   primMomX = primMomDirOnSurf.x();
0295   primMomY = primMomDirOnSurf.y();
0296   primMomZ = primMomDirOnSurf.z();
0297   //angles for rotation matrices

0298   double theta = primMomDirOnSurf.theta();
0299   double phi = primMomDirOnSurf.phi();
0300 
0301   // my insert ----------------------------------------------------------------

0302   double sphi = sin(phi);
0303   double cphi = cos(phi);
0304   double ctheta = cos(theta);
0305   double stheta = sin(theta);
0306 
0307   double pex = 0, pey = 0, zv = 0;
0308   double xx, yy, zz;
0309 
0310   for (unsigned int k = 0; k < LongFiberPhotons.size(); ++k) {
0311     HFShowerPhoton aPhoton = LongFiberPhotons[k];
0312     // global coordinates

0313     xx = aPhoton.x();
0314     yy = aPhoton.y();
0315     zz = aPhoton.z();
0316 
0317     // local coordinates in rotated to shower axis system and vs shower origin

0318     pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
0319     pey = -xx * sphi + yy * cphi;
0320     zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
0321 
0322     double photonProdTime = aPhoton.t() - primTimeOnSurf;
0323     thePhotons.push_back(Photon(1, pex, pey, zv, photonProdTime, aPhoton.lambda()));
0324   }
0325   for (unsigned int k = 0; k < ShortFiberPhotons.size(); ++k) {
0326     HFShowerPhoton aPhoton = ShortFiberPhotons[k];
0327     // global coordinates

0328     xx = aPhoton.x();
0329     yy = aPhoton.y();
0330     zz = aPhoton.z();
0331 
0332     // local coordinates in rotated to shower axis system and vs shower origin

0333     pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
0334     pey = -xx * sphi + yy * cphi;
0335     zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
0336 
0337     double photonProdTime = aPhoton.t() - primTimeOnSurf;
0338     thePhotons.push_back(Photon(2, pex, pey, zv, photonProdTime, aPhoton.lambda()));
0339   }
0340 }
0341 
0342 void HcalForwardAnalysis::fillEvent() {
0343 #ifdef EDM_ML_DEBUG
0344   edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> filledEvent";
0345 #endif
0346   nphot = int(thePhotons.size());
0347   for (int i = 0; i < nphot; ++i) {
0348     x[i] = thePhotons[i].x;
0349     y[i] = thePhotons[i].y;
0350     z[i] = thePhotons[i].z;
0351     t[i] = thePhotons[i].t;
0352     lambda[i] = thePhotons[i].lambda;
0353     fiberId[i] = thePhotons[i].fiberId;
0354   }
0355   theTree->Fill();
0356 }
0357 
0358 void HcalForwardAnalysis::parseDetId(int id, int& tower, int& cell, int& fiber) {
0359   tower = id / 10000;
0360   cell = id / 10 - tower * 10;
0361   fiber = id - tower * 10000 - cell * 10;
0362 }
0363 
0364 void HcalForwardAnalysis::clear() {
0365   nphot = 0;
0366   for (int i = 0; i < 10000; ++i) {
0367     x[i] = 0.;
0368     y[i] = 0.;
0369     z[i] = 0.;
0370     t[i] = 0.;
0371     lambda[i] = 0.;
0372     fiberId[i] = 0;
0373   }
0374   primX = primY = primZ = primT = 0.;
0375   primMomX = primMomY = primMomZ = 0.;
0376 
0377   thePhotons.clear();
0378 }
0379 
0380 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0381 #include "FWCore/PluginManager/interface/ModuleDef.h"
0382 
0383 DEFINE_SIMWATCHER(HcalForwardAnalysis);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!