CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MaterialBudgetCastorHistos.cc
Go to the documentation of this file.
1 #include "Validation/Geometry/interface/MaterialBudgetCastorHistos.h"
2 
3 #include "DetectorDescription/Core/interface/DDFilter.h"
4 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
5 #include "DetectorDescription/Core/interface/DDSplit.h"
6 #include "DetectorDescription/Core/interface/DDValue.h"
7 #include "FWCore/Utilities/interface/Exception.h"
8 #include "FWCore/MessageLogger/interface/MessageLogger.h"
9 #include "FWCore/ServiceRegistry/interface/Service.h"
10 #include "CommonTools/UtilAlgos/interface/TFileService.h"
11 
12 #include "CLHEP/Units/GlobalPhysicalConstants.h"
13 #include "CLHEP/Units/GlobalSystemOfUnits.h"
14 
15 MaterialBudgetCastorHistos::MaterialBudgetCastorHistos(const edm::ParameterSet &p){
16 
17  binEta = p.getUntrackedParameter<int>("NBinEta", 100);
18  binPhi = p.getUntrackedParameter<int>("NBinPhi", 180);
19  etaLow = p.getUntrackedParameter<double>("EtaLow", 5.0);
20  etaHigh = p.getUntrackedParameter<double>("EtaHigh", 7.0);
21  fillHistos = p.getUntrackedParameter<bool>("FillHisto", true);
22  printSum = p.getUntrackedParameter<bool>("PrintSummary", false);
23  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: FillHisto : "
24  << fillHistos << " PrintSummary " << printSum
25  << " == Eta plot: NX " << binEta << " Range "
26  << etaLow << " (" << -etaHigh << ") : "
27  << etaHigh << " (" << -etaLow <<") Phi plot: "
28  << "NX " << binPhi << " Range " << -pi << ":"
29  << pi << " (Eta limit " << etaLow << ":"
30  << etaHigh <<")";
31  if (fillHistos) book();
32 
33 }
34 
35 MaterialBudgetCastorHistos::~MaterialBudgetCastorHistos() {
36  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Save user "
37  << "histos ===";
38 }
39 
40 void MaterialBudgetCastorHistos::fillStartTrack(const G4Track* aTrack) {
41 
42  id1 = id2 = steps = 0;
43  radLen = intLen = stepLen = 0;
44 
45  const G4ThreeVector& dir = aTrack->GetMomentum() ;
46  if (dir.theta() != 0 ) {
47  eta = dir.eta();
48  } else {
49  eta = -99;
50  }
51  phi = dir.phi();
52  double theEnergy = aTrack->GetTotalEnergy();
53  int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
54 
55  if (printSum) {
56  matList.clear();
57  stepLength.clear();
58  radLength.clear();
59  intLength.clear();
60  }
61 
62  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Track "
63  << aTrack->GetTrackID() << " Code " << theID
64  << " Energy " << theEnergy/GeV << " GeV; Eta "
65  << eta << " Phi " << phi/deg << " PT "
66  << dir.perp()/GeV << " GeV *****";
67 }
68 
69 
70 void MaterialBudgetCastorHistos::fillPerStep(const G4Step* aStep) {
71 
72  G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
73  double step = aStep->GetStepLength();
74  double radl = material->GetRadlen();
75  double intl = material->GetNuclearInterLength();
76  double density = material->GetDensity() / (g/cm3);
77 
78  int id1Old = id1;
79  int id2Old = id2;
80  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
81  std::string name = touch->GetVolume(0)->GetName();
82  std::string matName = material->GetName();
83  if (printSum) {
84  bool found = false;
85  for (unsigned int ii=0; ii<matList.size(); ii++) {
86  if (matList[ii] == matName) {
87  stepLength[ii] += step;
88  radLength[ii] += (step/radl);
89  intLength[ii] += (step/intl);
90  found = true;
91  break;
92  }
93  }
94  if (!found) {
95  matList.push_back(matName);
96  stepLength.push_back(step);
97  radLength.push_back(step/radl);
98  intLength.push_back(step/intl);
99  }
100  edm::LogInfo("MaterialBudget") << name << " " << step << " " << matName
101  << " " << stepLen << " " << step/radl << " "
102  << radLen << " " <<step/intl << " " <<intLen;
103  } else {
104  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Step at "
105  << name << " Length " << step << " in "
106  << matName << " of density " << density
107  << " g/cc; Radiation Length " <<radl <<" mm;"
108  << " Interaction Length " << intl << " mm\n"
109  << " Position "
110  << aStep->GetPreStepPoint()->GetPosition()
111  << " Cylindrical R "
112  <<aStep->GetPreStepPoint()->GetPosition().perp()
113  << " Length (so far) " << stepLen << " L/X0 "
114  << step/radl << "/" << radLen << " L/Lambda "
115  << step/intl << "/" << intLen;
116  }
117 
118  int level = ((touch->GetHistoryDepth())+1);
119  std::string name1="XXXX", name2="XXXX";
120  if (level>3) name1 = touch->GetVolume(level-4)->GetName();
121  if (level>4) name2 = touch->GetVolume(level-5)->GetName();
122  if (name1 == "CAST") {
123  id1 = 1;
124  if (name2 == "CAEC") id2 = 2;
125  else if (name2 == "CAHC") id2 = 3;
126  else if (name2 == "CEDC") id2 = 4;
127  else if (name2 == "CHDC") id2 = 5;
128  else id2 = 0;
129  } else {
130  id1 = id2 = 0;
131  }
132  LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos: Level " << level
133  << " Volume " << name1 << " and " << name2
134  << " ID1 " << id1 << " (" << id1Old << ") ID2 "
135  << id2 << " (" << id2Old << ")";
136 
137  if (fillHistos) {
138  if (id1 != id1Old) {
139  if (id1 == 0) fillHisto(id1Old, 1);
140  else fillHisto(id1, 0);
141  }
142  if (id2 != id2Old) {
143  if (id2 == 0) fillHisto(id2Old, 1);
144  else fillHisto(id2, 0);
145  }
146  }
147 
148  stepLen += step;
149  radLen += step/radl;
150  intLen += step/intl;
151 }
152 
153 
154 void MaterialBudgetCastorHistos::fillEndTrack() {
155 
156  if (fillHistos) {
157  if (id1 != 0) fillHisto(id1, 1);
158  if (id2 != 0) fillHisto(id2, 1);
159  }
160  if (printSum) {
161  for (unsigned int ii=0; ii<matList.size(); ii++) {
162  edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
163  << "\t" << radLength[ii] << "\t"
164  << intLength[ii];
165  }
166  }
167 }
168 
169 void MaterialBudgetCastorHistos::book() {
170 
171  // Book histograms
172  edm::Service<TFileService> tfile;
173 
174  if ( !tfile.isAvailable() )
175  throw cms::Exception("BadConfig") << "TFileService unavailable: "
176  << "please add it to config file";
177 
178  double maxPhi=pi;
179  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
180  << "histos === with " << binEta << " bins "
181  << "in eta from " << etaLow << " to "
182  << etaHigh << " and " << binPhi << " bins "
183  << "in phi from " << -maxPhi << " to "
184  << maxPhi;
185 
186  char name[10], title[80], tag1[10], tag2[15];
187  // total X0
188  for (int i=0; i<maxSet; i++) {
189  double minEta=etaLow;
190  double maxEta=etaHigh;
191  int ireg = i;
192  if (i > 9) {
193  minEta = -etaHigh;
194  maxEta = -etaLow;
195  ireg -= 10;
196  sprintf (tag2, " (-ve Eta Side)");
197  } else {
198  sprintf (tag2, " (+ve Eta Side)");
199  }
200  if ((i%2) == 0) {
201  ireg /= 2;
202  sprintf (tag1, " == Start");
203  } else {
204  ireg = (ireg-1)/2;
205  sprintf (tag1, " == End");
206  }
207  sprintf(name, "%d", i+100);
208  sprintf(title, "MB(X0) prof Eta in region %d%s%s", ireg, tag1, tag2);
209  me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
210  sprintf(name, "%d", i+200);
211  sprintf(title, "MB(L0) prof Eta in region %d%s%s", ireg, tag1, tag2);
212  me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
213  sprintf(name, "%d", i+300);
214  sprintf(title, "MB(Step) prof Eta in region %d%s%s", ireg, tag1, tag2);
215  me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
216  sprintf(name, "%d", i+400);
217  sprintf(title, "Eta in region %d%s%s", ireg, tag1, tag2);
218  me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
219  sprintf(name, "%d", i+500);
220  sprintf(title, "MB(X0) prof Ph in region %d%s%s", ireg, tag1, tag2);
221  me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
222  sprintf(name, "%d", i+600);
223  sprintf(title, "MB(L0) prof Ph in region %d%s%s", ireg, tag1, tag2);
224  me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
225  sprintf(name, "%d", i+700);
226  sprintf(title, "MB(Step) prof Ph in region %d%s%s", ireg, tag1, tag2);
227  me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
228  sprintf(name, "%d", i+800);
229  sprintf(title, "Phi in region %d%s%s", ireg, tag1, tag2);
230  me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
231  sprintf(name, "%d", i+900);
232  sprintf(title, "MB(X0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
233  me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
234  binPhi/2, -maxPhi, maxPhi);
235  sprintf(name, "%d", i+1000);
236  sprintf(title, "MB(L0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
237  me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
238  binPhi/2, -maxPhi, maxPhi);
239  sprintf(name, "%d", i+1100);
240  sprintf(title, "MB(Step) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
241  me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
242  binPhi/2, -maxPhi, maxPhi);
243  sprintf(name, "%d", i+1200);
244  sprintf(title, "Eta vs Phi in region %d%s%s", ireg, tag1, tag2);
245  me1200[i]= tfile->make<TH2F>(name, title, binEta/2, minEta, maxEta,
246  binPhi/2, -maxPhi, maxPhi);
247  }
248 
249  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
250  << "histos done ===";
251 
252 }
253 
254 void MaterialBudgetCastorHistos::fillHisto(int id, int ix) {
255 
256  int ii = 2*(id-1) + ix;
257  double etaAbs = eta;
258  if (eta < 0) {
259  etaAbs = -eta;
260  ii += 10;
261  }
262  LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos:FillHisto "
263  << "called with index " << id << ":" << ix
264  << ":" << ii << " eta " << etaAbs << " ("
265  << eta << ") integrated step " << stepLen
266  << " X0 " << radLen << " Lamda " << intLen;
267 
268  me100[ii]->Fill(eta, radLen);
269  me200[ii]->Fill(eta, intLen);
270  me300[ii]->Fill(eta, stepLen);
271  me400[ii]->Fill(eta);
272 
273  if (etaAbs >= etaLow && etaAbs <= etaHigh) {
274  me500[ii]->Fill(phi, radLen);
275  me600[ii]->Fill(phi, intLen);
276  me700[ii]->Fill(phi, stepLen);
277  me800[ii]->Fill(phi);
278  }
279 
280  me900[ii]->Fill(eta, phi, radLen);
281  me1000[ii]->Fill(eta, phi, intLen);
282  me1100[ii]->Fill(eta, phi, stepLen);
283  me1200[ii]->Fill(eta, phi);
284 
285 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
GeV
const double GeV
Definition: MathUtil.h:16
MaterialBudgetCastorHistos::matList
std::vector< std::string > matList
Definition: MaterialBudgetCastorHistos.h:39
benchmark_cfg.minEta
tuple minEta
Definition: benchmark_cfg.py:51
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
MaterialBudgetCastorHistos::MaterialBudgetCastorHistos
MaterialBudgetCastorHistos(const edm::ParameterSet &p)
Definition: MaterialBudgetCastorHistos.cc:15
maxEta
double maxEta
Definition: PFJetBenchmarkAnalyzer.cc:78
cuy.ii
int ii
Definition: cuy.py:588
g
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
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
pi
const Double_t pi
Definition: trackSplitPlot.h:36
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:46
relval_parameters_module.step
string step
Definition: relval_parameters_module.py:36
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
testEve_cfg.level
tuple level
Definition: testEve_cfg.py:34
dir
dbl *** dir
Definition: mlp_gen.cc:35
Generated for CMSSW Reference Manual by doxygen 1.8.5