CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HcalForwardAnalysis.cc
Go to the documentation of this file.
1 // system include files
2 #include <cmath>
3 #include <iostream>
4 #include <iomanip>
5 #include <string>
6 #include <vector>
7 
8 // user include files
9 #include "DataFormats/Math/interface/Point3D.h"
10 
11 #include "FWCore/Framework/interface/Event.h"
12 #include "FWCore/Framework/interface/EventSetup.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 #include "FWCore/ParameterSet/interface/ParameterSet.h"
15 #include "SimG4Core/Notification/interface/BeginOfRun.h"
16 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
17 #include "SimG4Core/Notification/interface/EndOfEvent.h"
18 #include "SimG4Core/Watcher/interface/SimProducer.h"
19 #include "SimG4Core/Notification/interface/Observer.h"
20 
21 #include "SimG4CMS/ShowerLibraryProducer/interface/FiberG4Hit.h"
22 #include "SimG4CMS/ShowerLibraryProducer/interface/HFShowerG4Hit.h"
23 
24 #include "G4HCofThisEvent.hh"
25 #include "G4SDManager.hh"
26 #include "G4Step.hh"
27 #include "G4Track.hh"
28 #include "G4ThreeVector.hh"
29 #include "G4VProcess.hh"
30 
31 #include <CLHEP/Units/GlobalSystemOfUnits.h>
32 #include <CLHEP/Units/GlobalPhysicalConstants.h>
33 
34 #include "TFile.h"
35 #include "TTree.h"
36 
37 //#define EDM_ML_DEBUG
38 
39 class HcalForwardAnalysis : public SimProducer,
40  public Observer<const BeginOfRun*>,
41  public Observer<const BeginOfEvent*>,
42  public Observer<const EndOfEvent*>,
43  public Observer<const G4Step*> {
44 public:
45  struct Photon {
46  Photon(int id, float X, float Y, float Z, float T, float Lambda)
47  : fiberId(id), x(X), y(Y), z(Z), t(T), lambda(Lambda) {}
48  int fiberId;
49  float x;
50  float y;
51  float z;
52  float t;
53  float lambda;
54  };
55 
56  HcalForwardAnalysis(const edm::ParameterSet& p);
57  HcalForwardAnalysis(const HcalForwardAnalysis&) = delete; // stop default
58  const HcalForwardAnalysis& operator=(const HcalForwardAnalysis&) = delete;
59  ~HcalForwardAnalysis() override;
60 
61  void produce(edm::Event&, const edm::EventSetup&) override;
62 
63 private:
64  void init();
65 
66  // observer methods
67  void update(const BeginOfRun* run) override;
68  void update(const BeginOfEvent* evt) override;
69  void update(const G4Step* step) override;
70  void update(const EndOfEvent* evt) override;
71  // void write(const EndOfRun * run);
72 
73  //User methods
74  void setPhotons(const EndOfEvent* evt);
75  //void fillEvent(PHcalForwardLibInfo&);
76  void fillEvent();
77  void parseDetId(int id, int& tower, int& cell, int& fiber);
78  void clear();
79 
80  edm::Service<TFileService> theFile;
81  TTree* theTree;
82  int theEventCounter;
83  int count;
84  int evNum;
85  float x[10000], y[10000], z[10000], t[10000], lambda[10000];
86  float primX, primY, primZ, primT;
87  float primMomX, primMomY, primMomZ;
88  int nphot;
89  int fiberId[10000];
90  std::vector<Photon> thePhotons;
91  std::vector<std::string> theNames;
92  bool fillt;
93 };
94 
95 HcalForwardAnalysis::HcalForwardAnalysis(const edm::ParameterSet& p) {
96  edm::ParameterSet m_SLP = p.getParameter<edm::ParameterSet>("HFShowerLibraryProducer");
97  theNames = m_SLP.getParameter<std::vector<std::string> >("Names");
98  //LibVer = m_HS.getParameter<std::string> ("LibVer");
99  //produces<HFShowerPhotonCollection> ();
100  init();
101  theEventCounter = 0;
102  nphot = 0;
103  for (int i = 0; i < 10000; ++i) {
104  x[i] = 0.;
105  y[i] = 0.;
106  z[i] = 0.;
107  t[i] = 0.;
108  lambda[i] = 0.;
109  fiberId[i] = 0;
110  }
111  primX = primY = primZ = primT = 0.;
112  primMomX = primMomY = primMomZ = 0.;
113 }
114 
115 HcalForwardAnalysis::~HcalForwardAnalysis() {}
116 
117 //
118 // member functions
119 //
120 
121 void HcalForwardAnalysis::produce(edm::Event& iEvent, const edm::EventSetup&) {
122  if (fillt)
123  fillEvent();
124 }
125 
126 void HcalForwardAnalysis::init() {
127  theTree = theFile->make<TTree>("CherenkovPhotons", "Cherenkov Photons");
128  theTree->Branch("nphot", &nphot, "nphot/I");
129  theTree->Branch("x", &x, "x[nphot]/F");
130  theTree->Branch("y", &y, "y[nphot]/F");
131  theTree->Branch("z", &z, "z[nphot]/F");
132  theTree->Branch("t", &t, "t[nphot]/F");
133  theTree->Branch("lambda", &lambda, "lambda[nphot]/F");
134  theTree->Branch("fiberId", &fiberId, "fiberId[nphot]/I");
135  theTree->Branch("primX", &primX, "primX/F");
136  theTree->Branch("primY", &primY, "primY/F");
137  theTree->Branch("primZ", &primZ, "primZ/F");
138  theTree->Branch("primMomX", &primMomX, "primMomX/F");
139  theTree->Branch("primMomY", &primMomY, "primMomY/F");
140  theTree->Branch("primMomZ", &primMomZ, "primMomZ/F");
141  theTree->Branch("primT", &primT, "primT/F");
142 
143  // counter
144  count = 0;
145  evNum = 0;
146  clear();
147 }
148 
149 void HcalForwardAnalysis::update(const BeginOfRun* run) {
150  int irun = (*run)()->GetRunID();
151  edm::LogVerbatim("HcalForwardLib") << " =====> Begin of Run = " << irun;
152 }
153 
154 void HcalForwardAnalysis::update(const BeginOfEvent* evt) {
155  evNum = (*evt)()->GetEventID();
156  clear();
157 #ifdef EDM_ML_DEBUG
158  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> Begin of event = " << evNum;
159 #endif
160 }
161 
162 void HcalForwardAnalysis::update(const G4Step* aStep) {}
163 
164 void HcalForwardAnalysis::update(const EndOfEvent* evt) {
165  count++;
166 
167  //fill the buffer
168 #ifdef EDM_ML_DEBUG
169  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::Fill event " << (*evt)()->GetEventID();
170 #endif
171  setPhotons(evt);
172 
173  int iEvt = (*evt)()->GetEventID();
174  if (iEvt < 10)
175  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
176  else if ((iEvt < 100) && (iEvt % 10 == 0))
177  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
178  else if ((iEvt < 1000) && (iEvt % 100 == 0))
179  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
180  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
181  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
182 }
183 
184 void HcalForwardAnalysis::setPhotons(const EndOfEvent* evt) {
185  fillt = true;
186  int idHC, j;
187  FiberG4HitsCollection* theHC;
188  // Look for the Hit Collection of HCal
189  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
190 #ifdef EDM_ML_DEBUG
191  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Has " << allHC->GetNumberOfCollections()
192  << " collections";
193  for (int k = 0; k < allHC->GetNumberOfCollections(); ++k) {
194  G4String name = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetName();
195  G4String nameSD = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetSDname();
196  edm::LogVerbatim("HcalForwardLib") << "Collecttion[" << k << "] " << allHC->GetHC(k) << " " << name << ":"
197  << nameSD;
198  }
199 #endif
200  std::string sdName = theNames[0]; //name for fiber hits
201  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
202  theHC = (FiberG4HitsCollection*)allHC->GetHC(idHC);
203 #ifdef EDM_ML_DEBUG
204  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Hit Collection for " << sdName << " of ID "
205  << idHC << " is obtained at " << theHC;
206 #endif
207  std::vector<HFShowerPhoton> ShortFiberPhotons;
208  std::vector<HFShowerPhoton> LongFiberPhotons;
209  LongFiberPhotons.clear();
210  ShortFiberPhotons.clear();
211  if (idHC >= 0 && theHC != nullptr) {
212  int thehc_entries = theHC->entries();
213 #ifdef EDM_ML_DEBUG
214  edm::LogVerbatim("HcalForwardLib") << "FiberhitSize " << thehc_entries;
215 #endif
216  for (j = 0; j < thehc_entries; j++) {
217  FiberG4Hit* aHit = (*theHC)[j];
218  std::vector<HFShowerPhoton> thePhotonsFromHit = aHit->photon();
219 #ifdef EDM_ML_DEBUG
220  edm::LogVerbatim("HcalForwardLib") << "Fiberhit " << j << " has " << thePhotonsFromHit.size() << " photons.";
221 #endif
222  int fTowerId = -1;
223  int fCellId = -1;
224  int fFiberId = -1;
225  parseDetId(aHit->towerId(), fTowerId, fCellId, fFiberId);
226  for (unsigned int iph = 0; iph < thePhotonsFromHit.size(); ++iph) {
227  if (aHit->depth() == 1)
228  LongFiberPhotons.push_back(thePhotonsFromHit[iph]);
229  if (aHit->depth() == 2)
230  ShortFiberPhotons.push_back(thePhotonsFromHit[iph]);
231  }
232 #ifdef EDM_ML_DEBUG
233  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() NbPhotons " << thePhotonsFromHit.size()
234  << " towerId " << fTowerId << " cellId " << fCellId << " fiberId " << fFiberId
235  << " depth " << aHit->depth();
236 #endif
237  }
238  } else {
239  fillt = false;
240 #ifdef EDM_ML_DEBUG
241  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Photons!";
242 #endif
243  return;
244  }
245  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
246  << " ShortFibPhotons: " << ShortFiberPhotons.size();
247  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
248  << " ShortFibPhotons: " << ShortFiberPhotons.size();
249 
250  //Chamber hits to find information about primary particle on surface
251  HFShowerG4HitsCollection* theChamberHC;
252  G4HCofThisEvent* allChamberHC = (*evt)()->GetHCofThisEvent();
253  sdName = theNames[1];
254  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
255  theChamberHC = (HFShowerG4HitsCollection*)allChamberHC->GetHC(idHC);
256  math::XYZPoint primPosOnSurf(0, 0, 0);
257  math::XYZPoint primMomDirOnSurf(0, 0, 0);
258  float primTimeOnSurf = 0;
259  // the chamber hit is for primary particle, but step size can be small
260  // (in newer Geant4 versions) and as a result primary particle may have
261  // multiple hits. We want to take last one which is close the HF absorber
262  // if (idHC >= 0 && theChamberHC != nullptr && theChamberHC->entries()>0) {
263  if (idHC >= 0 && theChamberHC != nullptr) {
264  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hits size: "
265  << theChamberHC->entries();
266  int thec_hc_entries = theChamberHC->entries();
267  for (j = 0; j < thec_hc_entries; ++j) {
268  HFShowerG4Hit* aHit = (*theChamberHC)[j];
269 #ifdef EDM_ML_DEBUG
270  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hit id " << aHit->hitId()
271  << " track id " << aHit->trackId() << " prim. pos. " << aHit->globalPosition()
272  << " prom mom. dir. " << aHit->primaryMomDir() << " time " << aHit->time();
273 #endif
274  primPosOnSurf.SetXYZ(aHit->globalPosition().x(), aHit->globalPosition().y(), aHit->globalPosition().z());
275  primMomDirOnSurf.SetXYZ(aHit->primaryMomDir().x(), aHit->primaryMomDir().y(), aHit->primaryMomDir().z());
276  primTimeOnSurf = aHit->time();
277  }
278  } else {
279  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Chamber hits are stored";
280  fillt = false;
281  return;
282  }
283  primX = primPosOnSurf.x();
284  primY = primPosOnSurf.y();
285  primZ = primPosOnSurf.z();
286  if (primZ < 990) { // there were interactions before HF
287  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): First interaction before HF";
288  fillt = false;
289  return;
290  }
291  primT = primTimeOnSurf;
292  primMomX = primMomDirOnSurf.x();
293  primMomY = primMomDirOnSurf.y();
294  primMomZ = primMomDirOnSurf.z();
295  //angles for rotation matrices
296  double theta = primMomDirOnSurf.theta();
297  double phi = primMomDirOnSurf.phi();
298 
299  // my insert ----------------------------------------------------------------
300  double sphi = sin(phi);
301  double cphi = cos(phi);
302  double ctheta = cos(theta);
303  double stheta = sin(theta);
304 
305  double pex = 0, pey = 0, zv = 0;
306  double xx, yy, zz;
307 
308  for (unsigned int k = 0; k < LongFiberPhotons.size(); ++k) {
309  HFShowerPhoton aPhoton = LongFiberPhotons[k];
310  // global coordinates
311  xx = aPhoton.x();
312  yy = aPhoton.y();
313  zz = aPhoton.z();
314 
315  // local coordinates in rotated to shower axis system and vs shower origin
316  pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
317  pey = -xx * sphi + yy * cphi;
318  zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
319 
320  double photonProdTime = aPhoton.t() - primTimeOnSurf;
321  thePhotons.push_back(Photon(1, pex, pey, zv, photonProdTime, aPhoton.lambda()));
322  }
323  for (unsigned int k = 0; k < ShortFiberPhotons.size(); ++k) {
324  HFShowerPhoton aPhoton = ShortFiberPhotons[k];
325  // global coordinates
326  xx = aPhoton.x();
327  yy = aPhoton.y();
328  zz = aPhoton.z();
329 
330  // local coordinates in rotated to shower axis system and vs shower origin
331  pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
332  pey = -xx * sphi + yy * cphi;
333  zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
334 
335  double photonProdTime = aPhoton.t() - primTimeOnSurf;
336  thePhotons.push_back(Photon(2, pex, pey, zv, photonProdTime, aPhoton.lambda()));
337  }
338 }
339 
340 void HcalForwardAnalysis::fillEvent() {
341 #ifdef EDM_ML_DEBUG
342  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> filledEvent";
343 #endif
344  nphot = int(thePhotons.size());
345  for (int i = 0; i < nphot; ++i) {
346  x[i] = thePhotons[i].x;
347  y[i] = thePhotons[i].y;
348  z[i] = thePhotons[i].z;
349  t[i] = thePhotons[i].t;
350  lambda[i] = thePhotons[i].lambda;
351  fiberId[i] = thePhotons[i].fiberId;
352  }
353  theTree->Fill();
354 }
355 
356 void HcalForwardAnalysis::parseDetId(int id, int& tower, int& cell, int& fiber) {
357  tower = id / 10000;
358  cell = id / 10 - tower * 10;
359  fiber = id - tower * 10000 - cell * 10;
360 }
361 
362 void HcalForwardAnalysis::clear() {
363  nphot = 0;
364  for (int i = 0; i < 10000; ++i) {
365  x[i] = 0.;
366  y[i] = 0.;
367  z[i] = 0.;
368  t[i] = 0.;
369  lambda[i] = 0.;
370  fiberId[i] = 0;
371  }
372  primX = primY = primZ = primT = 0.;
373  primMomX = primMomY = primMomZ = 0.;
374 
375  thePhotons.clear();
376 }
377 
378 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
379 #include "FWCore/PluginManager/interface/ModuleDef.h"
380 
381 DEFINE_SIMWATCHER(HcalForwardAnalysis);
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
HFShowerG4Hit::primaryMomDir
G4ThreeVector primaryMomDir() const
Definition: HFShowerG4Hit.h:52
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
mps_fire.i
i
Definition: mps_fire.py:428
configurableAnalysis::Photon
constexpr char Photon[]
Definition: modules.cc:14
FiberG4Hit::depth
G4int depth() const
Definition: FiberG4Hit.h:45
HFShowerG4Hit::hitId
G4int hitId() const
Definition: HFShowerG4Hit.h:46
gpuVertexFinder::zv
float *__restrict__ zv
Definition: gpuFitVertices.h:26
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
X
#define X(str)
Definition: MuonsGrabber.cc:38
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
hgcalTowerProducer_cfi.tower
tuple tower
Definition: hgcalTowerProducer_cfi.py:4
HcalForwardAnalysis::setPhotons
void setPhotons(const EndOfEvent *evt)
Definition: HcalForwardAnalysis.cc:184
iEvent
int iEvent
Definition: GenABIO.cc:224
HFShowerPhoton::t
float t() const
Definition: HFShowerPhoton.h:28
HcalForwardAnalysis::Photon::Photon
Photon(int id, float X, float Y, float Z, float T, float Lambda)
Definition: HcalForwardAnalysis.cc:46
HcalForwardAnalysis::thePhotons
std::vector< Photon > thePhotons
Definition: HcalForwardAnalysis.cc:90
HcalForwardAnalysis::operator=
const HcalForwardAnalysis & operator=(const HcalForwardAnalysis &)=delete
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
FiberG4Hit::towerId
G4int towerId() const
Definition: FiberG4Hit.h:44
HFShowerG4Hit::trackId
G4int trackId() const
Definition: HFShowerG4Hit.h:47
FiberG4Hit::photon
std::vector< HFShowerPhoton > photon() const
Definition: FiberG4Hit.h:50
HFShowerPhoton::z
float z() const
Definition: HFShowerPhoton.h:26
HFShowerPhoton::lambda
float lambda() const
Definition: HFShowerPhoton.h:27
HcalForwardAnalysis::update
void update(const BeginOfRun *run) override
This routine will be called when the appropriate signal arrives.
Definition: HcalForwardAnalysis.cc:149
HcalForwardAnalysis::theFile
edm::Service< TFileService > theFile
Definition: HcalForwardAnalysis.cc:80
HFShowerPhoton::x
float x() const
Definition: HFShowerPhoton.h:24
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
HcalForwardAnalysis::theNames
std::vector< std::string > theNames
Definition: HcalForwardAnalysis.cc:91
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
FiberG4HitsCollection
G4THitsCollection< FiberG4Hit > FiberG4HitsCollection
Definition: FiberG4Hit.h:54
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HFShowerG4Hit::globalPosition
G4ThreeVector globalPosition() const
Definition: HFShowerG4Hit.h:51
step
step
Definition: StallMonitor.cc:94
HcalForwardAnalysis::HcalForwardAnalysis
HcalForwardAnalysis(const edm::ParameterSet &p)
Definition: HcalForwardAnalysis.cc:95
T
long double T
Definition: Basic3DVectorLD.h:48
HFShowerG4Hit::time
G4double time() const
Definition: HFShowerG4Hit.h:49
HFShowerG4HitsCollection
G4THitsCollection< HFShowerG4Hit > HFShowerG4HitsCollection
Definition: HFShowerG4Hit.h:55
HcalForwardAnalysis::parseDetId
void parseDetId(int id, int &tower, int &cell, int &fiber)
Definition: HcalForwardAnalysis.cc:356
HFShowerPhoton::y
float y() const
Definition: HFShowerPhoton.h:25
HcalForwardAnalysis::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HcalForwardAnalysis.cc:121
Generated for CMSSW Reference Manual by doxygen 1.8.5