CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MaterialBudgetHcalHistos.cc
Go to the documentation of this file.
2 
11 
12 #include "CLHEP/Units/GlobalPhysicalConstants.h"
13 #include "CLHEP/Units/GlobalSystemOfUnits.h"
14 
16 
17  binEta = p.getUntrackedParameter<int>("NBinEta", 260);
18  binPhi = p.getUntrackedParameter<int>("NBinPhi", 180);
19  maxEta = p.getUntrackedParameter<double>("MaxEta", 5.2);
20  etaLow = p.getUntrackedParameter<double>("EtaLow", -5.2);
21  etaHigh = p.getUntrackedParameter<double>("EtaHigh", 5.2);
22  fillHistos = p.getUntrackedParameter<bool>("FillHisto", true);
23  printSum = p.getUntrackedParameter<bool>("PrintSummary", false);
24  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: FillHisto : "
25  << fillHistos << " PrintSummary " << printSum
26  << " == Eta plot: NX " << binEta << " Range "
27  << -maxEta << ":" << maxEta <<" Phi plot: NX "
28  << binPhi << " Range " << -pi << ":" << pi
29  << " (Eta limit " << etaLow << ":" << etaHigh
30  <<")";
31  if (fillHistos) book();
32 
33 }
34 
36 
37  if (fillHistos) {
38  std::string attribute = "ReadOutName";
39  std::string value = "HcalHits";
40  DDSpecificsFilter filter1;
41  DDValue ddv1(attribute,value,0);
42  filter1.setCriteria(ddv1,DDCompOp::equals);
43  DDFilteredView fv1(cpv);
44  fv1.addFilter(filter1);
45  sensitives = getNames(fv1);
46  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
47  << "tested for " << attribute << " = "
48  << value << " has " << sensitives.size()
49  << " elements";
50  for (unsigned int i=0; i<sensitives.size(); i++)
51  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: sensitives["
52  << i << "] = " << sensitives[i];
53 
54  attribute = "Volume";
55  value = "HF";
56  DDSpecificsFilter filter2;
57  DDValue ddv2(attribute,value,0);
58  filter2.setCriteria(ddv2,DDCompOp::equals);
59  DDFilteredView fv2(cpv);
60  fv2.addFilter(filter2);
61  hfNames = getNames(fv2);
62  fv2.firstChild();
64  std::vector<double> temp = getDDDArray("Levels",sv);
65  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
66  << "tested for " << attribute << " = "
67  << value << " has " << hfNames.size()
68  << " elements";
69  for (unsigned int i=0; i<hfNames.size(); i++) {
70  int level = static_cast<int>(temp[i]);
71  hfLevels.push_back(level);
72  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: HF[" << i
73  << "] = " << hfNames[i] << " at level "
74  << hfLevels[i];
75  }
76 
77  std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
78  attribute = "ReadOutName";
79  for (int k=0; k<2; k++) {
80  value = ecalRO[k];
81  DDSpecificsFilter filter3;
82  DDValue ddv3(attribute,value,0);
83  filter3.setCriteria(ddv3,DDCompOp::equals);
84  DDFilteredView fv3(cpv);
85  fv3.addFilter(filter3);
86  std::vector<std::string> senstmp = getNames(fv3);
87  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be"
88  << " tested for " << attribute << " = "
89  << value << " has " << senstmp.size()
90  << " elements";
91  for (unsigned int i=0; i<senstmp.size(); i++)
92  sensitiveEC.push_back(senstmp[i]);
93  }
94  for (unsigned int i=0; i<sensitiveEC.size(); i++)
95  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos:sensitiveEC["
96  << i << "] = " << sensitiveEC[i];
97  }
98 }
99 
100 void MaterialBudgetHcalHistos::fillStartTrack(const G4Track* aTrack) {
101 
102  id = layer = steps = 0;
103  radLen = intLen = stepLen = 0;
104  nlayHB = nlayHE = nlayHF = nlayHO = 0;
105 
106  const G4ThreeVector& dir = aTrack->GetMomentum() ;
107  if (dir.theta() != 0 ) {
108  eta = dir.eta();
109  } else {
110  eta = -99;
111  }
112  phi = dir.phi();
113  double theEnergy = aTrack->GetTotalEnergy();
114  int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
115 
116  if (printSum) {
117  matList.clear();
118  stepLength.clear();
119  radLength.clear();
120  intLength.clear();
121  }
122 
123  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Track "
124  << aTrack->GetTrackID() << " Code " << theID
125  << " Energy " << theEnergy/GeV << " GeV; Eta "
126  << eta << " Phi " << phi/deg << " PT "
127  << dir.perp()/GeV << " GeV *****";
128 }
129 
130 
131 void MaterialBudgetHcalHistos::fillPerStep(const G4Step* aStep) {
132 
133  G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
134  double step = aStep->GetStepLength();
135  double radl = material->GetRadlen();
136  double intl = material->GetNuclearInterLength();
137  double density = material->GetDensity() / (g/cm3);
138 
139  int idOld = id;
140  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
141  std::string name = touch->GetVolume(0)->GetName();
142  std::string matName = material->GetName();
143  if (printSum) {
144  bool found = false;
145  for (unsigned int ii=0; ii<matList.size(); ii++) {
146  if (matList[ii] == matName) {
147  stepLength[ii] += step;
148  radLength[ii] += (step/radl);
149  intLength[ii] += (step/intl);
150  found = true;
151  break;
152  }
153  }
154  if (!found) {
155  matList.push_back(matName);
156  stepLength.push_back(step);
157  radLength.push_back(step/radl);
158  intLength.push_back(step/intl);
159  }
160  edm::LogInfo("MaterialBudget") << name << " " << step << " " << matName
161  << " " << stepLen << " " << step/radl << " "
162  << radLen << " " <<step/intl << " " <<intLen;
163  } else {
164  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Step at "
165  << name << " Length " << step << " in "
166  << matName << " of density " << density
167  << " g/cc; Radiation Length " <<radl <<" mm;"
168  << " Interaction Length " << intl << " mm\n"
169  << " Position "
170  << aStep->GetPreStepPoint()->GetPosition()
171  << " Cylindrical R "
172  <<aStep->GetPreStepPoint()->GetPosition().perp()
173  << " Length (so far) " << stepLen << " L/X0 "
174  << step/radl << "/" << radLen << " L/Lambda "
175  << step/intl << "/" << intLen;
176  }
177 
178  int det=0, lay=0;
179  if (fillHistos) {
180  if (isItEC(name)) {
181  det = 1;
182  lay = 1;
183  } else {
184  if (isSensitive(name)) {
185  if (isItHF(touch)) {
186  det = 5;
187  lay = 21;
188  if (lay != layer) nlayHF++;
189  } else {
190  det = (touch->GetReplicaNumber(1))/1000;
191  lay = (touch->GetReplicaNumber(0)/10)%100 + 3;
192  if (det == 4) {
193  double abeta = std::abs(eta);
194  if (abeta < 1.479) lay = layer + 1;
195  else lay--;
196  if (lay < 3) lay = 3;
197  if (lay == layer) lay++;
198  if (lay > 20) lay = 20;
199  if (lay != layer) nlayHE++;
200  } else if (lay != layer) {
201  if (lay >= 20) nlayHO++;
202  else nlayHB++;
203  }
204  }
205  LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: Det " << det
206  << " Layer " << lay << " Eta " << eta
207  << " Phi " << phi/deg;
208  } else if (layer == 1) {
209  det = -1;
210  lay = 2;
211  }
212  }
213  if (det != 0) {
214  if (lay != layer) {
215  id = lay;
216  layer = lay;
217  }
218  }
219 
220  if (id > idOld) {
221  // edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Step at " << name;
222  fillHisto(id-1);
223  }
224  }
225 
226  stepLen += step;
227  radLen += step/radl;
228  intLen += step/intl;
229  if (fillHistos) {
230  if (layer == 21 && det == 5) {
231  if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
232  LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: After HF in "
233  << aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName();
234  fillHisto(id);
235  id++;
236  layer = 0;
237  }
238  }
239  }
240 }
241 
242 
244  edm::LogInfo("MaterialBudget") << "Number of layers hit in HB:" << nlayHB
245  << " HE:" << nlayHE << " HO:" << nlayHO
246  << " HF:" << nlayHF;
247  if (fillHistos) {
248  fillHisto(maxSet-1);
249  fillLayer();
250  }
251  if (printSum) {
252  for (unsigned int ii=0; ii<matList.size(); ii++) {
253  edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
254  << "\t" << radLength[ii] << "\t"
255  << intLength[ii];
256  }
257  }
258 }
259 
261 
262  // Book histograms
264 
265  if ( !tfile.isAvailable() )
266  throw cms::Exception("BadConfig") << "TFileService unavailable: "
267  << "please add it to config file";
268 
269  double maxPhi=pi;
270  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
271  << "histos === with " << binEta << " bins "
272  << "in eta from " << -maxEta << " to "
273  << maxEta << " and " << binPhi << " bins "
274  << "in phi from " << -maxPhi << " to "
275  << maxPhi;
276 
277  char name[10], title[40];
278  // total X0
279  for (int i=0; i<maxSet; i++) {
280  sprintf(name, "%d", i+100);
281  sprintf(title, "MB(X0) prof Eta in region %d", i);
282  me100[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
283  sprintf(name, "%d", i+200);
284  sprintf(title, "MB(L0) prof Eta in region %d", i);
285  me200[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
286  sprintf(name, "%d", i+300);
287  sprintf(title, "MB(Step) prof Eta in region %d", i);
288  me300[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
289  sprintf(name, "%d", i+400);
290  sprintf(title, "Eta in region %d", i);
291  me400[i] = tfile->make<TH1F>(name, title, binEta, -maxEta, maxEta);
292  sprintf(name, "%d", i+500);
293  sprintf(title, "MB(X0) prof Ph in region %d", i);
294  me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
295  sprintf(name, "%d", i+600);
296  sprintf(title, "MB(L0) prof Ph in region %d", i);
297  me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
298  sprintf(name, "%d", i+700);
299  sprintf(title, "MB(Step) prof Ph in region %d", i);
300  me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
301  sprintf(name, "%d", i+800);
302  sprintf(title, "Phi in region %d", i);
303  me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
304  sprintf(name, "%d", i+900);
305  sprintf(title, "MB(X0) prof Eta Phi in region %d", i);
306  me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
307  binPhi/2, -maxPhi, maxPhi);
308  sprintf(name, "%d", i+1000);
309  sprintf(title, "MB(L0) prof Eta Phi in region %d", i);
310  me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
311  binPhi/2, -maxPhi, maxPhi);
312  sprintf(name, "%d", i+1100);
313  sprintf(title, "MB(Step) prof Eta Phi in region %d", i);
314  me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
315  binPhi/2, -maxPhi, maxPhi);
316  sprintf(name, "%d", i+1200);
317  sprintf(title, "Eta vs Phi in region %d", i);
318  me1200[i]= tfile->make<TH2F>(name, title, binEta/2, -maxEta, maxEta,
319  binPhi/2, -maxPhi, maxPhi);
320  }
321  for (int i=0; i<maxSet2; i++) {
322  sprintf(name, "%d", i+1300);
323  sprintf(title, "Events with layers Hit (0 all, 1 HB, ..) for %d", i);
324  me1300[i]= tfile->make<TH1F>(name, title, binEta, -maxEta, maxEta);
325  sprintf(name, "%d", i+1400);
326  sprintf(title, "Eta vs Phi for layers hit in %d", i);
327  me1400[i]= tfile->make<TH2F>(name, title, binEta/2, -maxEta, maxEta,
328  binPhi/2, -maxPhi, maxPhi);
329  sprintf(name, "%d", i+1500);
330  sprintf(title, "Number of layers crossed (0 all, 1 HB, ..) for %d", i);
331  me1500[i]= tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
332  }
333 
334  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
335  << "histos done ===";
336 
337 }
338 
340 
341  LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:FillHisto called "
342  << "with index " << ii << " integrated step "
343  << stepLen << " X0 " << radLen << " Lamda "
344  << intLen;
345 
346  if (ii >=0 && ii < maxSet) {
347  me100[ii]->Fill(eta, radLen);
348  me200[ii]->Fill(eta, intLen);
349  me300[ii]->Fill(eta, stepLen);
350  me400[ii]->Fill(eta);
351 
352  if (eta >= etaLow && eta <= etaHigh) {
353  me500[ii]->Fill(phi, radLen);
354  me600[ii]->Fill(phi, intLen);
355  me700[ii]->Fill(phi, stepLen);
356  me800[ii]->Fill(phi);
357  }
358 
359  me900[ii]->Fill(eta, phi, radLen);
360  me1000[ii]->Fill(eta, phi, intLen);
361  me1100[ii]->Fill(eta, phi, stepLen);
362  me1200[ii]->Fill(eta, phi);
363 
364  }
365 }
366 
368 
369  me1300[0]->Fill(eta);
370  me1400[0]->Fill(eta,phi);
371  if (nlayHB > 0) {
372  me1300[1]->Fill(eta);
373  me1400[1]->Fill(eta,phi);
374  }
375  if (nlayHB >= 16) {
376  me1300[2]->Fill(eta);
377  me1400[2]->Fill(eta,phi);
378  }
379  if (nlayHE > 0) {
380  me1300[3]->Fill(eta);
381  me1400[3]->Fill(eta,phi);
382  }
383  if (nlayHE >= 16) {
384  me1300[4]->Fill(eta);
385  me1400[4]->Fill(eta,phi);
386  }
387  if (nlayHO > 0) {
388  me1300[5]->Fill(eta);
389  me1400[5]->Fill(eta,phi);
390  }
391  if (nlayHO >= 2) {
392  me1300[6]->Fill(eta);
393  me1400[6]->Fill(eta,phi);
394  }
395  if (nlayHF > 0) {
396  me1300[7]->Fill(eta);
397  me1400[7]->Fill(eta,phi);
398  }
399  if (nlayHB > 0 || nlayHE > 0 || (nlayHF > 0 && std::abs(eta) > 3.0)) {
400  me1300[8]->Fill(eta);
401  me1400[8]->Fill(eta,phi);
402  }
403  me1500[0]->Fill(eta,(double)(nlayHB+nlayHO+nlayHE+nlayHF));
404  me1500[1]->Fill(eta,(double)(nlayHB));
405  me1500[2]->Fill(eta,(double)(nlayHE));
406  me1500[4]->Fill(eta,(double)(nlayHF));
407 }
408 
410  edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Save user "
411  << "histos ===";
412 }
413 
414 std::vector<std::string> MaterialBudgetHcalHistos::getNames(DDFilteredView& fv) {
415 
416  std::vector<std::string> tmp;
417  bool dodet = fv.firstChild();
418  while (dodet) {
419  const DDLogicalPart & log = fv.logicalPart();
420  std::string namx = log.name().name();
421  bool ok = true;
422  for (unsigned int i=0; i<tmp.size(); i++)
423  if (namx == tmp[i]) ok = false;
424  if (ok) tmp.push_back(namx);
425  dodet = fv.next();
426  }
427  return tmp;
428 }
429 
430 std::vector<double> MaterialBudgetHcalHistos::getDDDArray(const std::string & str,
431  const DDsvalues_type & sv) {
432 
433  LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:getDDDArray called "
434  << "for " << str;
435  DDValue value(str);
436  if (DDfetch(&sv,value)) {
437  LogDebug("MaterialBudget") << value;
438  const std::vector<double> & fvec = value.doubles();
439  int nval = fvec.size();
440  if (nval < 1) {
441  edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : # of "
442  << str << " bins " << nval
443  << " < 1 ==> illegal";
444  throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
445  << "nval < 1 for array " << str << "\n";
446  }
447 
448  return fvec;
449  } else {
450  edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : cannot get "
451  << "array " << str;
452  throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
453  << "cannot get array " << str << "\n";
454  }
455 }
456 
458 
459  std::vector<std::string>::const_iterator it = sensitives.begin();
460  std::vector<std::string>::const_iterator itEnd = sensitives.end();
461  for (; it != itEnd; ++it)
462  if (name == *it) return true;
463  return false;
464 }
465 
466 bool MaterialBudgetHcalHistos::isItHF (const G4VTouchable* touch) {
467 
468  // std::vector<std::string>::const_iterator it = hfNames.begin();
469  int levels = ((touch->GetHistoryDepth())+1);
470  for (unsigned int it=0; it < hfNames.size(); it++) {
471  if (levels >= hfLevels[it]) {
472  std::string name = touch->GetVolume(levels-hfLevels[it])->GetName();
473  if (name == hfNames[it]) return true;
474  }
475  }
476  return false;
477 }
478 
480 
481  std::vector<std::string>::const_iterator it = sensitiveEC.begin();
482  std::vector<std::string>::const_iterator itEnd = sensitiveEC.end();
483  for (; it != itEnd; ++it)
484  if (name == *it) return true;
485  return false;
486 }
487 
#define LogDebug(id)
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:139
const double GeV
Definition: MathUtil.h:16
const N & name() const
Definition: DDBase.h:78
std::vector< double > stepLength
void addFilter(const DDFilter &, DDLogOp op=DDLogOp::AND)
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
int ii
Definition: cuy.py:588
type of data representation of DDCompactView
Definition: DDCompactView.h:77
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:80
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
const Double_t pi
std::vector< std::string > getNames(DDFilteredView &fv)
std::vector< std::string > hfNames
void fillBeginJob(const DDCompactView &)
std::vector< std::string > matList
void fillStartTrack(const G4Track *)
bool next()
set current node to the next node in the filtered tree
std::vector< double > radLength
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
Definition: DDsvalues.h:19
bool isAvailable() const
Definition: Service.h:46
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
A DDLogicalPart aggregates information concerning material, solid and sensitveness ...
Definition: DDLogicalPart.h:88
std::vector< std::string > sensitiveEC
DDsvalues_type mergedSpecifics() const
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
bool firstChild()
set the current node to the first child ...
MaterialBudgetHcalHistos(const edm::ParameterSet &p)
std::vector< std::string > sensitives
tuple level
Definition: testEve_cfg.py:34
dbl *** dir
Definition: mlp_gen.cc:35
std::vector< double > getDDDArray(const std::string &str, const DDsvalues_type &sv)
const std::string & name() const
Returns the name.
Definition: DDName.cc:87
void setCriteria(const DDValue &nameVal, DDCompOp, DDLogOp l=DDLogOp::AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:245
std::vector< double > intLength
bool isItHF(const G4VTouchable *)
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:32