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 
16 #include <iostream>
17 //#define DebugLog
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(), (*lvcite)->GetName().c_str()) == 0) {
66  logVolumes[ii] = (*lvcite);
67  kount--;
68  break;
69  }
70  }
71  if (kount <= 0)
72  break;
73  }
74  edm::LogInfo("MaterialBudget") << "MaterialBudget::Finds " << detNames.size() - kount << " out of " << detNames.size()
75  << " LV addresses";
76  for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
77  std::string name("Unknown");
78  if (logVolumes[ii] != nullptr)
79  name = logVolumes[ii]->GetName();
80  edm::LogInfo("MaterialBudget") << "LV[" << ii << "] : " << detNames[ii] << " Address " << logVolumes[ii] << " | "
81  << name;
82  }
83 
84  for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
85  stepLen.push_back(0);
86  radLen.push_back(0);
87  intLen.push_back(0);
88  }
89  stackOrder.push_back(0);
90 }
91 
92 void MaterialBudget::update(const BeginOfTrack* trk) {
93  for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
94  stepLen[ii] = 0;
95  radLen[ii] = 0;
96  intLen[ii] = 0;
97  }
98 
99  const G4Track* aTrack = (*trk)(); // recover G4 pointer if wanted
100  const G4ThreeVector& mom = aTrack->GetMomentum();
101  if (mom.theta() != 0) {
102  eta_ = mom.eta();
103  } else {
104  eta_ = -99;
105  }
106  phi_ = mom.phi();
107  stepT = 0;
108 
109 #ifdef DebugLog
110  double theEnergy = aTrack->GetTotalEnergy();
111  int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
112  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code " << theID
113  << " Energy " << theEnergy / CLHEP::GeV << " GeV; Eta " << eta_ << " Phi "
114  << phi_ / CLHEP::deg << " PT " << mom.perp() / CLHEP::GeV << " GeV *****";
115 #endif
116 }
117 
118 void MaterialBudget::update(const G4Step* aStep) {
119  //---------- each step
120  G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
121  double step = aStep->GetStepLength();
122  double radl = material->GetRadlen();
123  double intl = material->GetNuclearInterLength();
124 
125  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
126  unsigned int indx = detTypes.size();
127  unsigned int nc = 0;
128  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
129  for (int kk = 0; kk < constituents[ii]; ++kk) {
130  if (detLevels[nc + kk] <= (int)((touch->GetHistoryDepth()) + 1)) {
131  int jj = (int)((touch->GetHistoryDepth()) + 1) - detLevels[nc + kk];
132  if ((touch->GetVolume(jj)->GetLogicalVolume()) == logVolumes[nc + kk]) {
133  indx = ii;
134  break;
135  }
136  }
137  }
138  nc += (unsigned int)(constituents[ii]);
139  if (indx == ii)
140  break;
141  }
142 
143  stepT += step;
144  stepLen[indx] = stepT;
145  radLen[indx] += (step / radl);
146  intLen[indx] += (step / intl);
147 #ifdef DebugLog
148  edm::LogInfo("MaterialBudget") << "MaterialBudget::Step in " << touch->GetVolume(0)->GetLogicalVolume()->GetName()
149  << " Index " << indx << " Step " << step << " RadL " << step / radl << " IntL "
150  << step / intl;
151 #endif
152  //----- Stop tracking after selected position
153  if (stopAfter(aStep)) {
154  G4Track* track = aStep->GetTrack();
155  track->SetTrackStatus(fStopAndKill);
156  }
157 }
158 
159 void MaterialBudget::update(const EndOfTrack* trk) {
160  for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
161  for (unsigned int jj = 0; jj <= detTypes.size(); ++jj) {
162  if (stackOrder[jj] == (int)(ii + 1)) {
163  for (unsigned int kk = 0; kk <= detTypes.size(); ++kk) {
164  if (stackOrder[kk] == (int)(ii)) {
165  radLen[jj] += radLen[kk];
166  intLen[jj] += intLen[kk];
167 #ifdef DebugLog
168  edm::LogInfo("MaterialBudget")
169  << "MaterialBudget::Add " << kk << ":" << stackOrder[kk] << " to " << jj << ":" << stackOrder[jj]
170  << " RadL " << radLen[kk] << " : " << radLen[jj] << " IntL " << intLen[kk] << " : " << intLen[jj]
171  << " StepL " << stepLen[kk] << " : " << stepLen[jj];
172 #endif
173  break;
174  }
175  }
176  break;
177  }
178  }
179  }
180 
181  for (unsigned int ii = 0; ii <= detTypes.size(); ++ii) {
182  me100[ii]->Fill(eta_, radLen[ii]);
183  me200[ii]->Fill(eta_, intLen[ii]);
184  me300[ii]->Fill(eta_, stepLen[ii]);
185  me400[ii]->Fill(phi_, radLen[ii]);
186  me500[ii]->Fill(phi_, intLen[ii]);
187  me600[ii]->Fill(phi_, stepLen[ii]);
188 #ifdef DebugLog
189  std::string name("Unknown");
190  if (ii < detTypes.size())
191  name = detTypes[ii];
192  edm::LogInfo("MaterialBudget") << "MaterialBudget::Volume[" << ii << "]: " << name << " eta " << eta_ << " == Step "
193  << stepLen[ii] << " RadL " << radLen[ii] << " IntL " << intLen[ii];
194 #endif
195  }
196 }
197 
198 void MaterialBudget::book(const edm::ParameterSet& m_p) {
199  // Book histograms
200  edm::Service<TFileService> tfile;
201 
202  if (!tfile.isAvailable())
203  throw cms::Exception("BadConfig") << "TFileService unavailable: "
204  << "please add it to config file";
205 
206  int binEta = m_p.getUntrackedParameter<int>("NBinEta", 320);
207  int binPhi = m_p.getUntrackedParameter<int>("NBinPhi", 180);
208  double minEta = m_p.getUntrackedParameter<double>("MinEta", -8.0);
209  double maxEta = m_p.getUntrackedParameter<double>("MaxEta", 8.0);
210  double maxPhi = CLHEP::pi;
211  edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
212  << "histos === with " << binEta << " bins "
213  << "in eta from " << minEta << " to " << maxEta << " and " << binPhi << " bins "
214  << "in phi from " << -maxPhi << " to " << maxPhi;
215 
216  char name[20], title[80];
217  std::string named;
218  for (unsigned int i = 0; i <= detTypes.size(); i++) {
219  if (i >= detTypes.size())
220  named = "Unknown";
221  else
222  named = detTypes[i];
223  sprintf(name, "%d", i + 100);
224  sprintf(title, "MB(X0) prof Eta in %s", named.c_str());
225  me100.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
226  sprintf(name, "%d", i + 200);
227  sprintf(title, "MB(L0) prof Eta in %s", named.c_str());
228  me200.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
229  sprintf(name, "%d", i + 300);
230  sprintf(title, "MB(Step) prof Eta in %s", named.c_str());
231  me300.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
232  sprintf(name, "%d", i + 400);
233  sprintf(title, "MB(X0) prof Phi in %s", named.c_str());
234  me400.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
235  sprintf(name, "%d", i + 500);
236  sprintf(title, "MB(L0) prof Phi in %s", named.c_str());
237  me500.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
238  sprintf(name, "%d", i + 600);
239  sprintf(title, "MB(Step) prof Phi in %s", named.c_str());
240  me600.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
241  }
242 
243  edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
244  << "histos done ===";
245 }
246 
247 bool MaterialBudget::stopAfter(const G4Step* aStep) {
248  G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
249  if (aStep->GetPostStepPoint() != nullptr)
250  hitPoint = aStep->GetPostStepPoint()->GetPosition();
251  double rr = hitPoint.perp();
252  double zz = std::abs(hitPoint.z());
253 
254  bool flag(false);
255  for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
256 #ifdef DebugLog
257  edm::LogInfo("MaterialBudget") << " MaterialBudget::Eta " << eta_ << " in Region[" << ii << "] with "
258  << etaRegions[ii] << " type " << regionTypes[ii] << "|" << boundaries[ii];
259 #endif
260  if (fabs(eta_) < etaRegions[ii]) {
261  if (regionTypes[ii] == 0) {
262  if (rr >= boundaries[ii] - 0.001)
263  flag = true;
264  } else {
265  if (zz >= boundaries[ii] - 0.001)
266  flag = true;
267  }
268 #ifdef DebugLog
269  if (flag)
270  edm::LogInfo("MaterialBudget") << " MaterialBudget::Stop after R = " << rr << " and Z = " << zz;
271 #endif
272  break;
273  }
274  }
275 #ifdef DebugLog
276  edm::LogInfo("MaterialBudget") << " MaterialBudget:: R = " << rr << " and Z = " << zz << " Flag " << flag;
277 #endif
278  return flag;
279 }
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:11713
MessageLogger.h
step
step
Definition: StallMonitor.cc:94
MaterialBudget::book
void book(const edm::ParameterSet &)
Definition: MaterialBudget.cc:198
findQualityFiles.rr
string rr
Definition: findQualityFiles.py:185
MaterialBudget::~MaterialBudget
~MaterialBudget() override
Definition: MaterialBudget.cc:53
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
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:52983
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
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
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
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
compare.tfile
tfile
Definition: compare.py:325
MaterialBudget::stopAfter
bool stopAfter(const G4Step *)
Definition: MaterialBudget.cc:247
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:590
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:116
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