CMS 3D CMS Logo
MaterialBudget.cc
Go to the documentation of this file.
1 #include "Validation/Geometry/interface/MaterialBudget.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 #include "DD4hep/Filter.h"
16 
17 #include <iostream>
18 
19 MaterialBudget::MaterialBudget(const edm::ParameterSet& p) {
20  edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudget");
21  detTypes = m_p.getParameter<std::vector<std::string> >("DetectorTypes");
22  constituents = m_p.getParameter<std::vector<int> >("Constituents");
23  stackOrder = m_p.getParameter<std::vector<int> >("StackOrder");
24  detNames = m_p.getParameter<std::vector<std::string> >("DetectorNames");
25  detLevels = m_p.getParameter<std::vector<int> >("DetectorLevels");
26  etaRegions = m_p.getParameter<std::vector<double> >("EtaBoundaries");
27  regionTypes = m_p.getParameter<std::vector<int> >("RegionTypes");
28  boundaries = m_p.getParameter<std::vector<double> >("Boundaries");
29  edm::LogInfo("MaterialBudget") << "MaterialBudget initialized for " << detTypes.size() << " detector types";
30  unsigned int nc = 0;
31  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
32  edm::LogInfo("MaterialBudget") << "Type[" << ii << "] : " << detTypes[ii] << " with " << constituents[ii]
33  << ", order " << stackOrder[ii] << " constituents --> ";
34  for (int kk = 0; kk < constituents[ii]; ++kk) {
35  std::string name = "Unknown";
36  int level = -1;
37  if (nc < detNames.size()) {
38  name = detNames[nc];
39  level = detLevels[nc];
40  ++nc;
41  }
42  edm::LogInfo("MaterialBudget") << " Constituent[" << kk << "]: " << name << " at level " << level;
43  }
44  }
45  edm::LogInfo("MaterialBudget") << "MaterialBudget Stop condition for " << etaRegions.size() << " eta regions";
46  for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
47  edm::LogInfo("MaterialBudget") << "Region[" << ii << "] : Eta < " << etaRegions[ii] << " boundary type "
48  << regionTypes[ii] << " limit " << boundaries[ii];
49  }
50  book(m_p);
51 }
52 
53 MaterialBudget::~MaterialBudget() {}
54 
55 void MaterialBudget::update(const BeginOfRun*) {
56  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
57  std::vector<G4LogicalVolume*>::const_iterator lvcite;
58 
59  unsigned int kount = detNames.size();
60  for (unsigned int ii = 0; ii < kount; ++ii)
61  logVolumes.push_back(nullptr);
62 
63  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
64  for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
65  if (strcmp(detNames[ii].c_str(),
66  static_cast<std::string>(dd4hep::dd::noNamespace((*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") << "MaterialBudget::Finds " << detNames.size() - kount << " out of " << detNames.size()
76  << " LV addresses";
77  for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
78  std::string name("Unknown");
79  if (logVolumes[ii] != nullptr)
80  name = dd4hep::dd::noNamespace(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 MaterialBudget::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 EDM_ML_DEBUG
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 MaterialBudget::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 EDM_ML_DEBUG
149  edm::LogInfo("MaterialBudget") << "MaterialBudget::Step in "
150  << dd4hep::dd::noNamespace(touch->GetVolume(0)->GetLogicalVolume()->GetName())
151  << " Index " << indx << " Step " << step << " RadL " << step / radl << " IntL "
152  << step / intl;
153 #endif
154  //----- Stop tracking after selected position
155  if (stopAfter(aStep)) {
156  G4Track* track = aStep->GetTrack();
157  track->SetTrackStatus(fStopAndKill);
158  }
159 }
160 
161 void MaterialBudget::update(const EndOfTrack* trk) {
162  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
163  for (unsigned int jj = 0; jj <= detTypes.size(); ++jj) {
164  if (stackOrder[jj] == (int)(ii + 1)) {
165  for (unsigned int kk = 0; kk <= detTypes.size(); ++kk) {
166  if (stackOrder[kk] == (int)(ii)) {
167  radLen[jj] += radLen[kk];
168  intLen[jj] += intLen[kk];
169 #ifdef EDM_ML_DEBUG
170  edm::LogInfo("MaterialBudget")
171  << "MaterialBudget::Add " << kk << ":" << stackOrder[kk] << " to " << jj << ":" << stackOrder[jj]
172  << " RadL " << radLen[kk] << " : " << radLen[jj] << " IntL " << intLen[kk] << " : " << intLen[jj]
173  << " StepL " << stepLen[kk] << " : " << stepLen[jj];
174 #endif
175  break;
176  }
177  }
178  break;
179  }
180  }
181  }
182 
183  for (unsigned int ii = 0; ii <= detTypes.size(); ++ii) {
184  me100[ii]->Fill(eta_, radLen[ii]);
185  me200[ii]->Fill(eta_, intLen[ii]);
186  me300[ii]->Fill(eta_, stepLen[ii]);
187  me400[ii]->Fill(phi_, radLen[ii]);
188  me500[ii]->Fill(phi_, intLen[ii]);
189  me600[ii]->Fill(phi_, stepLen[ii]);
190 #ifdef EDM_ML_DEBUG
191  std::string name("Unknown");
192  if (ii < detTypes.size())
193  name = detTypes[ii];
194  edm::LogInfo("MaterialBudget") << "MaterialBudget::Volume[" << ii << "]: " << name << " eta " << eta_ << " == Step "
195  << stepLen[ii] << " RadL " << radLen[ii] << " IntL " << intLen[ii];
196 #endif
197  }
198 }
199 
200 void MaterialBudget::book(const edm::ParameterSet& m_p) {
201  // Book histograms
202  edm::Service<TFileService> tfile;
203 
204  if (!tfile.isAvailable())
205  throw cms::Exception("BadConfig") << "TFileService unavailable: "
206  << "please add it to config file";
207 
208  int binEta = m_p.getUntrackedParameter<int>("NBinEta", 320);
209  int binPhi = m_p.getUntrackedParameter<int>("NBinPhi", 180);
210  double minEta = m_p.getUntrackedParameter<double>("MinEta", -8.0);
211  double maxEta = m_p.getUntrackedParameter<double>("MaxEta", 8.0);
212  double maxPhi = CLHEP::pi;
213  edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
214  << "histos === with " << binEta << " bins "
215  << "in eta from " << minEta << " to " << maxEta << " and " << binPhi << " bins "
216  << "in phi from " << -maxPhi << " to " << maxPhi;
217 
218  char name[20], title[80];
219  std::string named;
220  for (unsigned int i = 0; i <= detTypes.size(); i++) {
221  if (i >= detTypes.size())
222  named = "Unknown";
223  else
224  named = detTypes[i];
225  sprintf(name, "%d", i + 100);
226  sprintf(title, "MB(X0) prof Eta in %s", named.c_str());
227  me100.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
228  sprintf(name, "%d", i + 200);
229  sprintf(title, "MB(L0) prof Eta in %s", named.c_str());
230  me200.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
231  sprintf(name, "%d", i + 300);
232  sprintf(title, "MB(Step) prof Eta in %s", named.c_str());
233  me300.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
234  sprintf(name, "%d", i + 400);
235  sprintf(title, "MB(X0) prof Phi in %s", named.c_str());
236  me400.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
237  sprintf(name, "%d", i + 500);
238  sprintf(title, "MB(L0) prof Phi in %s", named.c_str());
239  me500.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
240  sprintf(name, "%d", i + 600);
241  sprintf(title, "MB(Step) prof Phi in %s", named.c_str());
242  me600.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
243  }
244 
245  edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
246  << "histos done ===";
247 }
248 
249 bool MaterialBudget::stopAfter(const G4Step* aStep) {
250  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
251  if (aStep->GetPostStepPoint() != nullptr)
252  hitPoint = aStep->GetPostStepPoint()->GetPosition();
253  double rr = hitPoint.perp();
254  double zz = std::abs(hitPoint.z());
255 
256  bool flag(false);
257  for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
258 #ifdef EDM_ML_DEBUG
259  edm::LogInfo("MaterialBudget") << " MaterialBudget::Eta " << eta_ << " in Region[" << ii << "] with "
260  << etaRegions[ii] << " type " << regionTypes[ii] << "|" << boundaries[ii];
261 #endif
262  if (fabs(eta_) < etaRegions[ii]) {
263  if (regionTypes[ii] == 0) {
264  if (rr >= boundaries[ii] - 0.001)
265  flag = true;
266  } else {
267  if (zz >= boundaries[ii] - 0.001)
268  flag = true;
269  }
270 #ifdef EDM_ML_DEBUG
271  if (flag)
272  edm::LogInfo("MaterialBudget") << " MaterialBudget::Stop after R = " << rr << " and Z = " << zz;
273 #endif
274  break;
275  }
276  }
277 #ifdef EDM_ML_DEBUG
278  edm::LogInfo("MaterialBudget") << " MaterialBudget:: R = " << rr << " and Z = " << zz << " Flag " << flag;
279 #endif
280  return flag;
281 }
personalPlayback.level
level
Definition: personalPlayback.py:22
MaterialBudget.h
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
MaterialBudget::me300
std::vector< TProfile * > me300
Definition: MaterialBudget.h:49
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:11724
MessageLogger.h
step
step
Definition: StallMonitor.cc:94
MaterialBudget::book
void book(const edm::ParameterSet &)
Definition: MaterialBudget.cc:200
findQualityFiles.rr
string rr
Definition: findQualityFiles.py:185
MaterialBudget::~MaterialBudget
~MaterialBudget() override
Definition: MaterialBudget.cc:53
MaterialBudget::radLen
std::vector< double > radLen
Definition: MaterialBudget.h:50
MaterialBudget::constituents
std::vector< int > constituents
Definition: MaterialBudget.h:46
MaterialBudget::me100
std::vector< TProfile * > me100
Definition: MaterialBudget.h:49
MaterialBudget::boundaries
std::vector< double > boundaries
Definition: MaterialBudget.h:47
EndOfTrack.h
MaterialBudget::logVolumes
std::vector< G4LogicalVolume * > logVolumes
Definition: MaterialBudget.h:48
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
MaterialBudget::stepLen
std::vector< double > stepLen
Definition: MaterialBudget.h:50
EndOfTrack
Definition: EndOfTrack.h:6
BeginOfTrack.h
HLT_FULL_cff.maxPhi
maxPhi
Definition: HLT_FULL_cff.py:52995
BeginOfRun.h
MaterialBudget::me200
std::vector< TProfile * > me200
Definition: MaterialBudget.h:49
Service.h
tfile
Definition: tfile.py:1
maxEta
double maxEta
Definition: PFJetBenchmarkAnalyzer.cc:76
BeginOfTrack
Definition: BeginOfTrack.h:6
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
MaterialBudget::me600
std::vector< TProfile * > me600
Definition: MaterialBudget.h:49
TFileService.h
MaterialBudget::regionTypes
std::vector< int > regionTypes
Definition: MaterialBudget.h:46
MaterialBudget::detLevels
std::vector< int > detLevels
Definition: MaterialBudget.h:46
edm::ParameterSet
Definition: ParameterSet.h:47
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
GeV
const double GeV
Definition: MathUtil.h:16
MaterialBudget::me400
std::vector< TProfile * > me400
Definition: MaterialBudget.h:49
MaterialBudget::MaterialBudget
MaterialBudget(const edm::ParameterSet &)
Definition: MaterialBudget.cc:19
edm::Service< TFileService >
createfilelist.int
int
Definition: createfilelist.py:10
MaterialBudget::intLen
std::vector< double > intLen
Definition: MaterialBudget.h:50
MaterialBudget::stackOrder
std::vector< int > stackOrder
Definition: MaterialBudget.h:46
BeginOfRun
Definition: BeginOfRun.h:6
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
compare.tfile
tfile
Definition: compare.py:324
MaterialBudget::stopAfter
bool stopAfter(const G4Step *)
Definition: MaterialBudget.cc:249
MaterialBudget::detTypes
std::vector< std::string > detTypes
Definition: MaterialBudget.h:45
MaterialBudget::eta_
double eta_
Definition: MaterialBudget.h:51
findQualityFiles.jj
string jj
Definition: findQualityFiles.py:188
MaterialBudget::phi_
double phi_
Definition: MaterialBudget.h:51
MaterialBudget::me500
std::vector< TProfile * > me500
Definition: MaterialBudget.h:49
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
MaterialBudget::detNames
std::vector< std::string > detNames
Definition: MaterialBudget.h:45
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
MaterialBudget::etaRegions
std::vector< double > etaRegions
Definition: MaterialBudget.h:47
MaterialBudget::stepT
double stepT
Definition: MaterialBudget.h:51
EgHLTOffEleSelection_cfi.minEta
minEta
Definition: EgHLTOffEleSelection_cfi.py:11
cuy.ii
ii
Definition: cuy.py:589
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:117
MaterialBudget::update
void update(const BeginOfRun *) override
This routine will be called when the appropriate signal arrives.
Definition: MaterialBudget.cc:55
Generated for CMSSW Reference Manual by doxygen 1.8.16