Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:24:40 

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 
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 //#define EDM_ML_DEBUG

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;  // stop default

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   // observer methods

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   //  void write(const EndOfRun * run);

0072 
0073   //User methods

0074   void setPhotons(const EndOfEvent* evt);
0075   //void fillEvent(PHcalForwardLibInfo&);

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   //LibVer = m_HS.getParameter<std::string> ("LibVer");

0098   //produces<HFShowerPhotonCollection> ();

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 // member functions

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   // counter

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   //fill the buffer

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   // Look for the Hit Collection of HCal

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];  //name for fiber hits

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   //Chamber hits to find information about primary particle on surface

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   //    the chamber hit is for primary particle, but step size can be small

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

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

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

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) {  // there were interactions before HF

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   //angles for rotation matrices

0296   double theta = primMomDirOnSurf.theta();
0297   double phi = primMomDirOnSurf.phi();
0298 
0299   // my insert ----------------------------------------------------------------

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     // global coordinates

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

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     // global coordinates

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

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);
[ 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!