CMS 3D CMS Logo
MaterialBudgetForward.cc
Go to the documentation of this file.
1 #include "Validation/Geometry/interface/MaterialBudgetForward.h"
2 
3 #include "FWCore/Utilities/interface/Exception.h"
4 #include "FWCore/MessageLogger/interface/MessageLogger.h"
5 #include "FWCore/ServiceRegistry/interface/Service.h"
6 #include "CommonTools/UtilAlgos/interface/TFileService.h"
7 
8 #include "SimG4Core/Notification/interface/BeginOfRun.h"
9 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
10 #include "SimG4Core/Notification/interface/EndOfTrack.h"
11 
12 #include "G4LogicalVolumeStore.hh"
13 #include "G4Step.hh"
14 #include "G4Track.hh"
15 
16 #include <iostream>
17 //#define DebugLog
18 const int MaterialBudgetForward::maxSet;
19 
20 MaterialBudgetForward::MaterialBudgetForward(const edm::ParameterSet& p) {
21  edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetForward");
22  detTypes = m_p.getParameter<std::vector<std::string> >("DetectorTypes");
23  constituents = m_p.getParameter<std::vector<int> >("Constituents");
24  stackOrder = m_p.getParameter<std::vector<int> >("StackOrder");
25  detNames = m_p.getParameter<std::vector<std::string> >("DetectorNames");
26  detLevels = m_p.getParameter<std::vector<int> >("DetectorLevels");
27  etaRegions = m_p.getParameter<std::vector<double> >("EtaBoundaries");
28  regionTypes = m_p.getParameter<std::vector<int> >("RegionTypes");
29  boundaries = m_p.getParameter<std::vector<double> >("Boundaries");
30  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward initialized for " << detTypes.size() << " detector types";
31  unsigned int nc = 0;
32  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
33  edm::LogInfo("MaterialBudget") << "Type[" << ii << "] : " << detTypes[ii] << " with " << constituents[ii]
34  << ", order " << stackOrder[ii] << " constituents --> ";
35  for (int kk = 0; kk < constituents[ii]; ++kk) {
36  std::string name = "Unknown";
37  int level = -1;
38  if (nc < detNames.size()) {
39  name = detNames[nc];
40  level = detLevels[nc];
41  ++nc;
42  }
43  edm::LogInfo("MaterialBudget") << " Constituent[" << kk << "]: " << name << " at level " << level;
44  }
45  }
46  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward Stop condition for " << etaRegions.size() << " eta regions";
47  for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
48  edm::LogInfo("MaterialBudget") << "Region[" << ii << "] : Eta < " << etaRegions[ii] << " boundary type "
49  << regionTypes[ii] << " limit " << boundaries[ii];
50  }
51  book(m_p);
52 }
53 
54 MaterialBudgetForward::~MaterialBudgetForward() {}
55 
56 void MaterialBudgetForward::update(const BeginOfRun*) {
57  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
58  std::vector<G4LogicalVolume*>::const_iterator lvcite;
59 
60  unsigned int kount = detNames.size();
61  for (unsigned int ii = 0; ii < kount; ++ii)
62  logVolumes.push_back(nullptr);
63 
64  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
65  for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
66  if (strcmp(detNames[ii].c_str(), (*lvcite)->GetName().c_str()) == 0) {
67  logVolumes[ii] = (*lvcite);
68  kount--;
69  break;
70  }
71  }
72  if (kount <= 0)
73  break;
74  }
75  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward::Finds " << detNames.size() - kount << " out of "
76  << detNames.size() << " LV addresses";
77  for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
78  std::string name("Unknown");
79  if (logVolumes[ii] != nullptr)
80  name = logVolumes[ii]->GetName();
81  edm::LogInfo("MaterialBudget") << "LV[" << ii << "] : " << detNames[ii] << " Address " << logVolumes[ii] << " | "
82  << name;
83  }
84 
85  for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
86  stepLen.push_back(0);
87  radLen.push_back(0);
88  intLen.push_back(0);
89  }
90  stackOrder.push_back(0);
91 }
92 
93 void MaterialBudgetForward::update(const BeginOfTrack* trk) {
94  for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
95  stepLen[ii] = 0;
96  radLen[ii] = 0;
97  intLen[ii] = 0;
98  }
99 
100  const G4Track* aTrack = (*trk)(); // recover G4 pointer if wanted
101  const G4ThreeVector& mom = aTrack->GetMomentum();
102  if (mom.theta() != 0) {
103  eta_ = mom.eta();
104  } else {
105  eta_ = -99;
106  }
107  phi_ = mom.phi();
108  stepT = 0;
109 
110 #ifdef DebugLog
111  double theEnergy = aTrack->GetTotalEnergy();
112  int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
113  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code " << theID
114  << " Energy " << theEnergy / CLHEP::GeV << " GeV; Eta " << eta_ << " Phi "
115  << phi_ / CLHEP::deg << " PT " << mom.perp() / CLHEP::GeV << " GeV *****";
116 #endif
117 }
118 
119 void MaterialBudgetForward::update(const G4Step* aStep) {
120  //---------- each step
121  G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
122  double step = aStep->GetStepLength();
123  double radl = material->GetRadlen();
124  double intl = material->GetNuclearInterLength();
125 
126  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
127  unsigned int indx = detTypes.size();
128  unsigned int nc = 0;
129  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
130  for (int kk = 0; kk < constituents[ii]; ++kk) {
131  if (detLevels[nc + kk] <= (int)((touch->GetHistoryDepth()) + 1)) {
132  int jj = (int)((touch->GetHistoryDepth()) + 1) - detLevels[nc + kk];
133  if ((touch->GetVolume(jj)->GetLogicalVolume()) == logVolumes[nc + kk]) {
134  indx = ii;
135  break;
136  }
137  }
138  }
139  nc += (unsigned int)(constituents[ii]);
140  if (indx == ii)
141  break;
142  }
143 
144  stepT += step;
145  stepLen[indx] = stepT;
146  radLen[indx] += (step / radl);
147  intLen[indx] += (step / intl);
148 #ifdef DebugLog
149  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward::Step in "
150  << touch->GetVolume(0)->GetLogicalVolume()->GetName() << " Index " << indx << " Step "
151  << step << " RadL " << step / radl << " IntL " << step / intl;
152 #endif
153  //----- Stop tracking after selected position
154  if (stopAfter(aStep)) {
155  G4Track* track = aStep->GetTrack();
156  track->SetTrackStatus(fStopAndKill);
157  }
158 }
159 
160 void MaterialBudgetForward::update(const EndOfTrack* trk) {
161  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
162  for (unsigned int jj = 0; jj <= detTypes.size(); ++jj) {
163  if (stackOrder[jj] == (int)(ii + 1)) {
164  for (unsigned int kk = 0; kk <= detTypes.size(); ++kk) {
165  if (stackOrder[kk] == (int)(ii)) {
166  radLen[jj] += radLen[kk];
167  intLen[jj] += intLen[kk];
168 #ifdef DebugLog
169  edm::LogInfo("MaterialBudget")
170  << "MaterialBudgetForward::Add " << kk << ":" << stackOrder[kk] << " to " << jj << ":" << stackOrder[jj]
171  << " RadL " << radLen[kk] << " : " << radLen[jj] << " IntL " << intLen[kk] << " : " << intLen[jj]
172  << " StepL " << stepLen[kk] << " : " << stepLen[jj];
173 #endif
174  break;
175  }
176  }
177  break;
178  }
179  }
180  }
181 
182  for (unsigned int ii = 0; ii <= detTypes.size(); ++ii) {
183  me100[ii]->Fill(eta_, radLen[ii]);
184  me200[ii]->Fill(eta_, intLen[ii]);
185  me300[ii]->Fill(eta_, stepLen[ii]);
186  me400[ii]->Fill(eta_);
187  me500[ii]->Fill(eta_, phi_, radLen[ii]);
188  me600[ii]->Fill(eta_, phi_, intLen[ii]);
189  me700[ii]->Fill(eta_, phi_, stepLen[ii]);
190  me800[ii]->Fill(eta_, phi_);
191 
192  std::string name("Unknown");
193  if (ii < detTypes.size())
194  name = detTypes[ii];
195 #ifdef DebugLog
196  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward::Volume[" << ii << "]: " << name << " eta " << eta_
197  << " == Step " << stepLen[ii] << " RadL " << radLen[ii] << " IntL " << intLen[ii];
198 #endif
199  }
200 }
201 
202 void MaterialBudgetForward::book(const edm::ParameterSet& m_p) {
203  // Book histograms
204  edm::Service<TFileService> tfile;
205 
206  if (!tfile.isAvailable())
207  throw cms::Exception("BadConfig") << "TFileService unavailable: "
208  << "please add it to config file";
209 
210  int binEta = m_p.getUntrackedParameter<int>("NBinEta", 320);
211  int binPhi = m_p.getUntrackedParameter<int>("NBinPhi", 180);
212  double minEta = m_p.getUntrackedParameter<double>("MinEta", -8.0);
213  double maxEta = m_p.getUntrackedParameter<double>("MaxEta", 8.0);
214  double maxPhi = CLHEP::pi;
215  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward: Booking user "
216  << "histos === with " << binEta << " bins "
217  << "in eta from " << minEta << " to " << maxEta << " and " << binPhi << " bins "
218  << "in phi from " << -maxPhi << " to " << maxPhi;
219 
220  char name[20], title[80];
221  std::string named;
222  for (int i = 0; i < std::min((int)(detTypes.size() + 1), maxSet); i++) {
223  if (i >= (int)(detTypes.size()))
224  named = "Unknown";
225  else
226  named = detTypes[i];
227  sprintf(name, "%d", i + 100);
228  sprintf(title, "MB(X0) prof Eta in %s", named.c_str());
229  me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
230  sprintf(name, "%d", i + 200);
231  sprintf(title, "MB(L0) prof Eta in %s", named.c_str());
232  me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
233  sprintf(name, "%d", i + 300);
234  sprintf(title, "MB(Step) prof Eta in %s", named.c_str());
235  me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
236  sprintf(name, "%d", i + 400);
237  sprintf(title, "Eta in %s", named.c_str());
238  me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
239  sprintf(name, "%d", i + 500);
240  sprintf(title, "MB(X0) prof Eta Phi in %s", named.c_str());
241  me500[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
242  sprintf(name, "%d", i + 600);
243  sprintf(title, "MB(L0) prof Eta Phi in %s", named.c_str());
244  me600[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
245  sprintf(name, "%d", i + 700);
246  sprintf(title, "MB(Step) prof Eta Phi in %s", named.c_str());
247  me700[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
248  sprintf(name, "%d", i + 800);
249  sprintf(title, "Eta vs Phi in %s", named.c_str());
250  me800[i] = tfile->make<TH2F>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
251  }
252 
253  edm::LogInfo("MaterialBudget") << "MaterialBudgetForward: Booking user "
254  << "histos done ===";
255 }
256 
257 bool MaterialBudgetForward::stopAfter(const G4Step* aStep) {
258  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
259  if (aStep->GetPostStepPoint() != nullptr)
260  hitPoint = aStep->GetPostStepPoint()->GetPosition();
261  double rr = hitPoint.perp();
262  double zz = std::abs(hitPoint.z());
263 
264  bool flag(false);
265  for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
266 #ifdef DebugLog
267  edm::LogInfo("MaterialBudget") << " MaterialBudgetForward::Eta " << eta_ << " in Region[" << ii << "] with "
268  << etaRegions[ii] << " type " << regionTypes[ii] << "|" << boundaries[ii];
269 #endif
270  if (fabs(eta_) < etaRegions[ii]) {
271  if (regionTypes[ii] == 0) {
272  if (rr >= boundaries[ii] - 0.001)
273  flag = true;
274  } else {
275  if (zz >= boundaries[ii] - 0.001)
276  flag = true;
277  }
278 #ifdef DebugLog
279  if (flag)
280  edm::LogInfo("MaterialBudget") << " MaterialBudgetForward::Stop after R = " << rr << " and Z = " << zz;
281 #endif
282  break;
283  }
284  }
285 #ifdef DebugLog
286  edm::LogInfo("MaterialBudget") << " MaterialBudgetForward:: R = " << rr << " and Z = " << zz << " Flag " << flag;
287 #endif
288  return flag;
289 }
personalPlayback.level
level
Definition: personalPlayback.py:22
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
mps_fire.i
i
Definition: mps_fire.py:428
geometryCSVtoXML.zz
zz
Definition: geometryCSVtoXML.py:19
HLT_FULL_cff.track
track
Definition: HLT_FULL_cff.py:11713
MessageLogger.h
MaterialBudgetForward::detNames
std::vector< std::string > detNames
Definition: MaterialBudgetForward.h:45
step
step
Definition: StallMonitor.cc:94
MaterialBudgetForward::me200
TProfile * me200[maxSet]
Definition: MaterialBudgetForward.h:52
min
T min(T a, T b)
Definition: MathUtil.h:58
MaterialBudgetForward::logVolumes
std::vector< G4LogicalVolume * > logVolumes
Definition: MaterialBudgetForward.h:48
MaterialBudgetForward::constituents
std::vector< int > constituents
Definition: MaterialBudgetForward.h:46
findQualityFiles.rr
string rr
Definition: findQualityFiles.py:185
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
MaterialBudgetForward::maxSet
static const int maxSet
Definition: MaterialBudgetForward.h:49
MaterialBudgetForward::stepT
double stepT
Definition: MaterialBudgetForward.h:55
EndOfTrack.h
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
MaterialBudgetForward.h
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
EndOfTrack
Definition: EndOfTrack.h:6
BeginOfTrack.h
HLT_FULL_cff.maxPhi
maxPhi
Definition: HLT_FULL_cff.py:52987
BeginOfRun.h
MaterialBudgetForward::me300
TProfile * me300[maxSet]
Definition: MaterialBudgetForward.h:52
MaterialBudgetForward::boundaries
std::vector< double > boundaries
Definition: MaterialBudgetForward.h:47
MaterialBudgetForward::me100
TProfile * me100[maxSet]
Definition: MaterialBudgetForward.h:52
Service.h
tfile
Definition: tfile.py:1
maxEta
double maxEta
Definition: PFJetBenchmarkAnalyzer.cc:76
MaterialBudgetForward::me400
TH1F * me400[maxSet]
Definition: MaterialBudgetForward.h:50
BeginOfTrack
Definition: BeginOfTrack.h:6
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
MaterialBudgetForward::regionTypes
std::vector< int > regionTypes
Definition: MaterialBudgetForward.h:46
MaterialBudgetForward::update
void update(const BeginOfRun *) override
This routine will be called when the appropriate signal arrives.
Definition: MaterialBudgetForward.cc:56
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TFileService.h
edm::ParameterSet
Definition: ParameterSet.h:47
MaterialBudgetForward::eta_
double eta_
Definition: MaterialBudgetForward.h:55
MaterialBudgetForward::phi_
double phi_
Definition: MaterialBudgetForward.h:55
MaterialBudgetForward::stopAfter
bool stopAfter(const G4Step *)
Definition: MaterialBudgetForward.cc:257
GeV
const double GeV
Definition: MathUtil.h:16
MaterialBudgetForward::radLen
std::vector< double > radLen
Definition: MaterialBudgetForward.h:54
edm::Service< TFileService >
createfilelist.int
int
Definition: createfilelist.py:10
MaterialBudgetForward::stackOrder
std::vector< int > stackOrder
Definition: MaterialBudgetForward.h:46
BeginOfRun
Definition: BeginOfRun.h:6
MaterialBudgetForward::me800
TH2F * me800[maxSet]
Definition: MaterialBudgetForward.h:51
compare.tfile
tfile
Definition: compare.py:325
MaterialBudgetForward::~MaterialBudgetForward
~MaterialBudgetForward() override
Definition: MaterialBudgetForward.cc:54
MaterialBudgetForward::me500
TProfile2D * me500[maxSet]
Definition: MaterialBudgetForward.h:53
MaterialBudgetForward::me600
TProfile2D * me600[maxSet]
Definition: MaterialBudgetForward.h:53
MaterialBudgetForward::detLevels
std::vector< int > detLevels
Definition: MaterialBudgetForward.h:46
findQualityFiles.jj
string jj
Definition: findQualityFiles.py:188
MaterialBudgetForward::detTypes
std::vector< std::string > detTypes
Definition: MaterialBudgetForward.h:45
MaterialBudgetForward::stepLen
std::vector< double > stepLen
Definition: MaterialBudgetForward.h:54
MaterialBudgetForward::MaterialBudgetForward
MaterialBudgetForward(const edm::ParameterSet &)
Definition: MaterialBudgetForward.cc:20
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
Exception.h
pi
const Double_t pi
Definition: trackSplitPlot.h:36
cms::Exception
Definition: Exception.h:70
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
MaterialBudgetForward::book
void book(const edm::ParameterSet &)
Definition: MaterialBudgetForward.cc:202
EgHLTOffEleSelection_cfi.minEta
minEta
Definition: EgHLTOffEleSelection_cfi.py:11
cuy.ii
ii
Definition: cuy.py:590
MaterialBudgetForward::etaRegions
std::vector< double > etaRegions
Definition: MaterialBudgetForward.h:47
MaterialBudgetForward::intLen
std::vector< double > intLen
Definition: MaterialBudgetForward.h:54
MaterialBudgetForward::me700
TProfile2D * me700[maxSet]
Definition: MaterialBudgetForward.h:53
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:116
Generated for CMSSW Reference Manual by doxygen 1.8.16