CMS 3D CMS Logo
HcalForwardAnalysis.cc
Go to the documentation of this file.
1 // system include files
2 #include <cmath>
3 #include <iostream>
4 #include <iomanip>
5 
6 // user include files
7 #include "SimG4CMS/ShowerLibraryProducer/interface/HcalForwardAnalysis.h"
8 
9 #include "SimG4Core/Notification/interface/BeginOfRun.h"
10 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
11 #include "SimG4Core/Notification/interface/EndOfEvent.h"
12 
13 #include "DataFormats/Math/interface/Point3D.h"
14 
15 #include "FWCore/Framework/interface/Event.h"
16 #include "FWCore/Framework/interface/EventSetup.h"
17 #include "FWCore/MessageLogger/interface/MessageLogger.h"
18 #include "FWCore/ParameterSet/interface/ParameterSet.h"
19 
20 #include "G4SDManager.hh"
21 #include "G4VProcess.hh"
22 #include "G4HCofThisEvent.hh"
23 #include "CLHEP/Units/GlobalSystemOfUnits.h"
24 #include "CLHEP/Units/GlobalPhysicalConstants.h"
25 
26 HcalForwardAnalysis::HcalForwardAnalysis(const edm::ParameterSet& p) {
27  edm::ParameterSet m_SLP = p.getParameter<edm::ParameterSet>("HFShowerLibraryProducer");
28  theNames = m_SLP.getParameter<std::vector<std::string> >("Names");
29  //LibVer = m_HS.getParameter<std::string> ("LibVer");
30  //produces<HFShowerPhotonCollection> ();
31  init();
32  theEventCounter = 0;
33  nphot = 0;
34  for (int i = 0; i < 10000; ++i) {
35  x[i] = 0.;
36  y[i] = 0.;
37  z[i] = 0.;
38  t[i] = 0.;
39  lambda[i] = 0.;
40  fiberId[i] = 0;
41  }
42  primX = primY = primZ = primT = 0.;
43  primMomX = primMomY = primMomZ = 0.;
44 }
45 
46 HcalForwardAnalysis::~HcalForwardAnalysis() {}
47 
48 //
49 // member functions
50 //
51 
52 void HcalForwardAnalysis::produce(edm::Event& iEvent, const edm::EventSetup&) {
53  if (fillt)
54  fillEvent();
55 }
56 
57 void HcalForwardAnalysis::init() {
58  theTree = theFile->make<TTree>("CherenkovPhotons", "Cherenkov Photons");
59  theTree->Branch("nphot", &nphot, "nphot/I");
60  theTree->Branch("x", &x, "x[nphot]/F");
61  theTree->Branch("y", &y, "y[nphot]/F");
62  theTree->Branch("z", &z, "z[nphot]/F");
63  theTree->Branch("t", &t, "t[nphot]/F");
64  theTree->Branch("lambda", &lambda, "lambda[nphot]/F");
65  theTree->Branch("fiberId", &fiberId, "fiberId[nphot]/I");
66  theTree->Branch("primX", &primX, "primX/F");
67  theTree->Branch("primY", &primY, "primY/F");
68  theTree->Branch("primZ", &primZ, "primZ/F");
69  theTree->Branch("primMomX", &primMomX, "primMomX/F");
70  theTree->Branch("primMomY", &primMomY, "primMomY/F");
71  theTree->Branch("primMomZ", &primMomZ, "primMomZ/F");
72  theTree->Branch("primT", &primT, "primT/F");
73 
74  // counter
75  count = 0;
76  evNum = 0;
77  clear();
78 }
79 
80 void HcalForwardAnalysis::update(const BeginOfRun* run) {
81  int irun = (*run)()->GetRunID();
82  edm::LogVerbatim("HcalForwardLib") << " =====> Begin of Run = " << irun;
83 }
84 
85 void HcalForwardAnalysis::update(const BeginOfEvent* evt) {
86  evNum = (*evt)()->GetEventID();
87  clear();
88  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> Begin of event = " << evNum;
89 }
90 
91 void HcalForwardAnalysis::update(const G4Step* aStep) {}
92 
93 void HcalForwardAnalysis::update(const EndOfEvent* evt) {
94  count++;
95 
96  //fill the buffer
97  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::Fill event " << (*evt)()->GetEventID();
98  setPhotons(evt);
99 
100  int iEvt = (*evt)()->GetEventID();
101  if (iEvt < 10)
102  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
103  else if ((iEvt < 100) && (iEvt % 10 == 0))
104  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
105  else if ((iEvt < 1000) && (iEvt % 100 == 0))
106  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
107  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
108  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Event " << iEvt;
109 }
110 
111 void HcalForwardAnalysis::setPhotons(const EndOfEvent* evt) {
112  fillt = true;
113  int idHC, j;
114  FiberG4HitsCollection* theHC;
115  // Look for the Hit Collection of HCal
116  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
117  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis:: Has " << allHC->GetNumberOfCollections()
118  << " collections";
119  for (int k = 0; k < allHC->GetNumberOfCollections(); ++k) {
120  G4String name = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetName();
121  G4String nameSD = (allHC->GetHC(k) == nullptr) ? "Unknown" : allHC->GetHC(k)->GetSDname();
122  edm::LogVerbatim("HcalForwardLib") << "Collecttion[" << k << "] " << allHC->GetHC(k) << " " << name << ":"
123  << nameSD;
124  }
125  std::string sdName = theNames[0]; //name for fiber hits
126  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
127  theHC = (FiberG4HitsCollection*)allHC->GetHC(idHC);
128  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Hit Collection for " << sdName << " of ID "
129  << idHC << " is obtained at " << theHC;
130  std::vector<HFShowerPhoton> ShortFiberPhotons;
131  std::vector<HFShowerPhoton> LongFiberPhotons;
132  LongFiberPhotons.clear();
133  ShortFiberPhotons.clear();
134  if (idHC >= 0 && theHC != nullptr) {
135  int thehc_entries = theHC->entries();
136  edm::LogVerbatim("HcalForwardLib") << "FiberhitSize " << thehc_entries;
137  for (j = 0; j < thehc_entries; j++) {
138  FiberG4Hit* aHit = (*theHC)[j];
139  std::vector<HFShowerPhoton> thePhotonsFromHit = aHit->photon();
140  edm::LogVerbatim("HcalForwardLib") << "Fiberhit " << j << " has " << thePhotonsFromHit.size() << " photons.";
141  int fTowerId = -1;
142  int fCellId = -1;
143  int fFiberId = -1;
144  parseDetId(aHit->towerId(), fTowerId, fCellId, fFiberId);
145  for (unsigned int iph = 0; iph < thePhotonsFromHit.size(); ++iph) {
146  if (aHit->depth() == 1)
147  LongFiberPhotons.push_back(thePhotonsFromHit[iph]);
148  if (aHit->depth() == 2)
149  ShortFiberPhotons.push_back(thePhotonsFromHit[iph]);
150  }
151  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() NbPhotons " << thePhotonsFromHit.size()
152  << " towerId " << fTowerId << " cellId " << fCellId << " fiberId " << fFiberId
153  << " depth " << aHit->depth();
154  }
155  } else {
156  fillt = false;
157  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Photons!";
158  return;
159  }
160  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
161  << " ShortFibPhotons: " << ShortFiberPhotons.size();
162  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() LongFibPhotons: " << LongFiberPhotons.size()
163  << " ShortFibPhotons: " << ShortFiberPhotons.size();
164 
165  //Chamber hits to find information about primary particle on surface
166  HFShowerG4HitsCollection* theChamberHC;
167  G4HCofThisEvent* allChamberHC = (*evt)()->GetHCofThisEvent();
168  sdName = theNames[1];
169  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
170  theChamberHC = (HFShowerG4HitsCollection*)allChamberHC->GetHC(idHC);
171  math::XYZPoint primPosOnSurf(0, 0, 0);
172  math::XYZPoint primMomDirOnSurf(0, 0, 0);
173  float primTimeOnSurf = 0;
174  // the chamber hit is for primary particle, but step size can be small
175  // (in newer Geant4 versions) and as a result primary particle may have
176  // multiple hits. We want to take last one which is close the HF absorber
177  // if (idHC >= 0 && theChamberHC != nullptr && theChamberHC->entries()>0) {
178  if (idHC >= 0 && theChamberHC != nullptr) {
179  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hits size: "
180  << theChamberHC->entries();
181  int thec_hc_entries = theChamberHC->entries();
182  for (j = 0; j < thec_hc_entries; ++j) {
183  HFShowerG4Hit* aHit = (*theChamberHC)[j];
184  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons() Chamber Hit id " << aHit->hitId()
185  << " track id " << aHit->trackId() << " prim. pos. " << aHit->globalPosition()
186  << " prom mom. dir. " << aHit->primaryMomDir() << " time " << aHit->time();
187  primPosOnSurf.SetXYZ(aHit->globalPosition().x(), aHit->globalPosition().y(), aHit->globalPosition().z());
188  primMomDirOnSurf.SetXYZ(aHit->primaryMomDir().x(), aHit->primaryMomDir().y(), aHit->primaryMomDir().z());
189  primTimeOnSurf = aHit->time();
190  }
191  } else {
192  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): No Chamber hits are stored";
193  fillt = false;
194  return;
195  }
196  primX = primPosOnSurf.x();
197  primY = primPosOnSurf.y();
198  primZ = primPosOnSurf.z();
199  if (primZ < 990) { // there were interactions before HF
200  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis::setPhotons(): First interaction before HF";
201  fillt = false;
202  return;
203  }
204  primT = primTimeOnSurf;
205  primMomX = primMomDirOnSurf.x();
206  primMomY = primMomDirOnSurf.y();
207  primMomZ = primMomDirOnSurf.z();
208  //angles for rotation matrices
209  double theta = primMomDirOnSurf.theta();
210  double phi = primMomDirOnSurf.phi();
211 
212  // my insert ----------------------------------------------------------------
213  double sphi = sin(phi);
214  double cphi = cos(phi);
215  double ctheta = cos(theta);
216  double stheta = sin(theta);
217 
218  double pex = 0, pey = 0, zv = 0;
219  double xx, yy, zz;
220 
221  for (unsigned int k = 0; k < LongFiberPhotons.size(); ++k) {
222  HFShowerPhoton aPhoton = LongFiberPhotons[k];
223  // global coordinates
224  xx = aPhoton.x();
225  yy = aPhoton.y();
226  zz = aPhoton.z();
227 
228  // local coordinates in rotated to shower axis system and vs shower origin
229  pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
230  pey = -xx * sphi + yy * cphi;
231  zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
232 
233  double photonProdTime = aPhoton.t() - primTimeOnSurf;
234  thePhotons.push_back(Photon(1, pex, pey, zv, photonProdTime, aPhoton.lambda()));
235  }
236  for (unsigned int k = 0; k < ShortFiberPhotons.size(); ++k) {
237  HFShowerPhoton aPhoton = ShortFiberPhotons[k];
238  // global coordinates
239  xx = aPhoton.x();
240  yy = aPhoton.y();
241  zz = aPhoton.z();
242 
243  // local coordinates in rotated to shower axis system and vs shower origin
244  pex = xx * ctheta * cphi + yy * ctheta * sphi - zz * stheta;
245  pey = -xx * sphi + yy * cphi;
246  zv = xx * stheta * cphi + yy * stheta * sphi + zz * ctheta - primZ / ctheta;
247 
248  double photonProdTime = aPhoton.t() - primTimeOnSurf;
249  thePhotons.push_back(Photon(2, pex, pey, zv, photonProdTime, aPhoton.lambda()));
250  }
251 }
252 
253 void HcalForwardAnalysis::fillEvent() {
254  /*
255  edm::LogVerbatim("HcalForwardLib") << "HcalForwardAnalysis: =====> filledEvent";
256  */
257  nphot = int(thePhotons.size());
258  for (int i = 0; i < nphot; ++i) {
259  x[i] = thePhotons[i].x;
260  y[i] = thePhotons[i].y;
261  z[i] = thePhotons[i].z;
262  t[i] = thePhotons[i].t;
263  lambda[i] = thePhotons[i].lambda;
264  fiberId[i] = thePhotons[i].fiberId;
265  }
266  theTree->Fill();
267 }
268 
269 void HcalForwardAnalysis::parseDetId(int id, int& tower, int& cell, int& fiber) {
270  tower = id / 10000;
271  cell = id / 10 - tower * 10;
272  fiber = id - tower * 10000 - cell * 10;
273 }
274 
275 void HcalForwardAnalysis::clear() {
276  nphot = 0;
277  for (int i = 0; i < 10000; ++i) {
278  x[i] = 0.;
279  y[i] = 0.;
280  z[i] = 0.;
281  t[i] = 0.;
282  lambda[i] = 0.;
283  fiberId[i] = 0;
284  }
285  primX = primY = primZ = primT = 0.;
286  primMomX = primMomY = primMomZ = 0.;
287 
288  thePhotons.clear();
289 }
HcalForwardAnalysis::fillt
bool fillt
Definition: HcalForwardAnalysis.h:87
HcalForwardAnalysis::y
float y[10000]
Definition: HcalForwardAnalysis.h:80
HFShowerPhoton::y
float y() const
Definition: HFShowerPhoton.h:25
HcalForwardAnalysis::primY
float primY
Definition: HcalForwardAnalysis.h:81
mps_fire.i
i
Definition: mps_fire.py:428
geometryCSVtoXML.zz
zz
Definition: geometryCSVtoXML.py:19
MessageLogger.h
HFShowerG4Hit::primaryMomDir
G4ThreeVector primaryMomDir() const
Definition: HFShowerG4Hit.h:52
HFShowerG4HitsCollection
G4THitsCollection< HFShowerG4Hit > HFShowerG4HitsCollection
Definition: HFShowerG4Hit.h:55
HcalForwardAnalysis::primMomZ
float primMomZ
Definition: HcalForwardAnalysis.h:82
HcalForwardAnalysis::t
float t[10000]
Definition: HcalForwardAnalysis.h:80
HcalForwardAnalysis::thePhotons
std::vector< Photon > thePhotons
Definition: HcalForwardAnalysis.h:85
HcalForwardAnalysis.h
hgcalTowerProducer_cfi.tower
tower
Definition: hgcalTowerProducer_cfi.py:4
HcalForwardAnalysis::primX
float primX
Definition: HcalForwardAnalysis.h:81
HFShowerPhoton::t
float t() const
Definition: HFShowerPhoton.h:28
gpuVertexFinder::zv
float *__restrict__ zv
Definition: gpuFitVertices.h:26
HcalForwardAnalysis::setPhotons
void setPhotons(const EndOfEvent *evt)
Definition: HcalForwardAnalysis.cc:111
HcalForwardAnalysis::init
void init()
Definition: HcalForwardAnalysis.cc:57
EndOfEvent.h
HcalForwardAnalysis::parseDetId
void parseDetId(int id, int &tower, int &cell, int &fiber)
Definition: HcalForwardAnalysis.cc:269
HcalForwardAnalysis::x
float x[10000]
Definition: HcalForwardAnalysis.h:80
HcalForwardAnalysis::theEventCounter
int theEventCounter
Definition: HcalForwardAnalysis.h:77
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
FiberG4Hit
Definition: FiberG4Hit.h:14
HFShowerG4Hit
Definition: HFShowerG4Hit.h:15
BeginOfRun.h
HcalForwardAnalysis::evNum
int evNum
Definition: HcalForwardAnalysis.h:79
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HFShowerG4Hit::hitId
G4int hitId() const
Definition: HFShowerG4Hit.h:46
HFShowerG4Hit::time
G4double time() const
Definition: HFShowerG4Hit.h:49
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
HcalForwardAnalysis::lambda
float lambda[10000]
Definition: HcalForwardAnalysis.h:80
dqmdumpme.k
k
Definition: dqmdumpme.py:60
EndOfEvent
Definition: EndOfEvent.h:6
HFShowerPhoton::x
float x() const
Definition: HFShowerPhoton.h:24
HcalForwardAnalysis::primZ
float primZ
Definition: HcalForwardAnalysis.h:81
HcalForwardAnalysis::nphot
int nphot
Definition: HcalForwardAnalysis.h:83
HcalForwardAnalysis::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HcalForwardAnalysis.cc:52
HcalForwardAnalysis::z
float z[10000]
Definition: HcalForwardAnalysis.h:80
HcalForwardAnalysis::theFile
edm::Service< TFileService > theFile
Definition: HcalForwardAnalysis.h:75
HcalForwardAnalysis::count
int count
Definition: HcalForwardAnalysis.h:78
edm::ParameterSet
Definition: ParameterSet.h:47
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
HcalForwardAnalysis::theTree
TTree * theTree
Definition: HcalForwardAnalysis.h:76
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
FiberG4Hit::depth
G4int depth() const
Definition: FiberG4Hit.h:45
Event.h
geometryCSVtoXML.yy
yy
Definition: geometryCSVtoXML.py:19
BeginOfEvent.h
FiberG4Hit::towerId
G4int towerId() const
Definition: FiberG4Hit.h:44
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
edm::EventSetup
Definition: EventSetup.h:58
HFShowerPhoton::lambda
float lambda() const
Definition: HFShowerPhoton.h:27
nanoDQM_cff.Photon
Photon
Definition: nanoDQM_cff.py:96
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
HcalForwardAnalysis::primMomY
float primMomY
Definition: HcalForwardAnalysis.h:82
HFShowerPhoton::z
float z() const
Definition: HFShowerPhoton.h:26
DDAxes::phi
writedatasetfile.run
run
Definition: writedatasetfile.py:27
HFShowerG4Hit::globalPosition
G4ThreeVector globalPosition() const
Definition: HFShowerG4Hit.h:51
HcalForwardAnalysis::fillEvent
void fillEvent()
Definition: HcalForwardAnalysis.cc:253
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
HcalForwardAnalysis::primT
float primT
Definition: HcalForwardAnalysis.h:81
HcalForwardAnalysis::~HcalForwardAnalysis
~HcalForwardAnalysis() override
Definition: HcalForwardAnalysis.cc:46
HcalForwardAnalysis::primMomX
float primMomX
Definition: HcalForwardAnalysis.h:82
HcalForwardAnalysis::fiberId
int fiberId[10000]
Definition: HcalForwardAnalysis.h:84
HFShowerPhoton
Definition: HFShowerPhoton.h:13
Point3D.h
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
HcalForwardAnalysis::theNames
std::vector< std::string > theNames
Definition: HcalForwardAnalysis.h:86
HcalForwardAnalysis::clear
void clear()
Definition: HcalForwardAnalysis.cc:275
HFShowerG4Hit::trackId
G4int trackId() const
Definition: HFShowerG4Hit.h:47
ParameterSet.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
edm::Event
Definition: Event.h:73
HcalForwardAnalysis::update
void update(const BeginOfRun *run) override
This routine will be called when the appropriate signal arrives.
Definition: HcalForwardAnalysis.cc:80
FiberG4HitsCollection
G4THitsCollection< FiberG4Hit > FiberG4HitsCollection
Definition: FiberG4Hit.h:54
FiberG4Hit::photon
std::vector< HFShowerPhoton > photon() const
Definition: FiberG4Hit.h:50
geometryCSVtoXML.xx
xx
Definition: geometryCSVtoXML.py:19
HcalForwardAnalysis::HcalForwardAnalysis
HcalForwardAnalysis(const edm::ParameterSet &p)
Definition: HcalForwardAnalysis.cc:26
Generated for CMSSW Reference Manual by doxygen 1.8.16