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  LogDebug("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  LogDebug("MaterialBudget") << name << " " << step << " " << matName
101  << " " << stepLen << " " << step/radl << " "
102  << radLen << " " << step/intl << " " << intLen;
103  } else {
104  LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos: Step at " << name
105  << " Length " << step << " in " << matName
106  << " of density " << density << " g/cc;"
107  << " 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  if (fillHistos) {
156  if (id1 != 0) fillHisto(id1, 1);
157  if (id2 != 0) fillHisto(id2, 1);
158  }
159  if (printSum) {
160  for (unsigned int ii=0; ii<matList.size(); ii++) {
161  edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
162  << "\t" << radLength[ii] << "\t"
163  << intLength[ii];
164  }
165  }
166 }
167 
168 void MaterialBudgetCastorHistos::book() {
169 
170  // Book histograms
171  edm::Service<TFileService> tfile;
172 
173  if ( !tfile.isAvailable() )
174  throw cms::Exception("BadConfig") << "TFileService unavailable: "
175  << "please add it to config file";
176 
177  double maxPhi=pi;
178  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
179  << "histos === with " << binEta << " bins "
180  << "in eta from " << etaLow << " to "
181  << etaHigh << " and " << binPhi << " bins "
182  << "in phi from " << -maxPhi << " to "
183  << maxPhi;
184 
185  char name[10], title[80], tag1[10], tag2[15];
186  // total X0
187  for (int i=0; i<maxSet; i++) {
188  double minEta=etaLow;
189  double maxEta=etaHigh;
190  int ireg = i;
191  if (i > 9) {
192  minEta = -etaHigh;
193  maxEta = -etaLow;
194  ireg -= 10;
195  sprintf (tag2, " (-ve Eta Side)");
196  } else {
197  sprintf (tag2, " (+ve Eta Side)");
198  }
199  if ((i%2) == 0) {
200  ireg /= 2;
201  sprintf (tag1, " == Start");
202  } else {
203  ireg = (ireg-1)/2;
204  sprintf (tag1, " == End");
205  }
206  sprintf(name, "%d", i+100);
207  sprintf(title, "MB(X0) prof Eta in region %d%s%s", ireg, tag1, tag2);
208  me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
209  sprintf(name, "%d", i+200);
210  sprintf(title, "MB(L0) prof Eta in region %d%s%s", ireg, tag1, tag2);
211  me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
212  sprintf(name, "%d", i+300);
213  sprintf(title, "MB(Step) prof Eta in region %d%s%s", ireg, tag1, tag2);
214  me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
215  sprintf(name, "%d", i+400);
216  sprintf(title, "Eta in region %d%s%s", ireg, tag1, tag2);
217  me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
218  sprintf(name, "%d", i+500);
219  sprintf(title, "MB(X0) prof Ph in region %d%s%s", ireg, tag1, tag2);
220  me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
221  sprintf(name, "%d", i+600);
222  sprintf(title, "MB(L0) prof Ph in region %d%s%s", ireg, tag1, tag2);
223  me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
224  sprintf(name, "%d", i+700);
225  sprintf(title, "MB(Step) prof Ph in region %d%s%s", ireg, tag1, tag2);
226  me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
227  sprintf(name, "%d", i+800);
228  sprintf(title, "Phi in region %d%s%s", ireg, tag1, tag2);
229  me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
230  sprintf(name, "%d", i+900);
231  sprintf(title, "MB(X0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
232  me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
233  binPhi/2, -maxPhi, maxPhi);
234  sprintf(name, "%d", i+1000);
235  sprintf(title, "MB(L0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
236  me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
237  binPhi/2, -maxPhi, maxPhi);
238  sprintf(name, "%d", i+1100);
239  sprintf(title, "MB(Step) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
240  me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
241  binPhi/2, -maxPhi, maxPhi);
242  sprintf(name, "%d", i+1200);
243  sprintf(title, "Eta vs Phi in region %d%s%s", ireg, tag1, tag2);
244  me1200[i]= tfile->make<TH2F>(name, title, binEta/2, minEta, maxEta,
245  binPhi/2, -maxPhi, maxPhi);
246  }
247 
248  edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
249  << "histos done ===";
250 
251 }
252 
253 void MaterialBudgetCastorHistos::fillHisto(int id, int ix) {
254 
255  int ii = 2*(id-1) + ix;
256  double etaAbs = eta;
257  if (eta < 0) {
258  etaAbs = -eta;
259  ii += 10;
260  }
261  LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos:FillHisto called "
262  << "with index " << id << ":" << ix << ":" << ii
263  << " eta " << etaAbs << " (" << eta << ")"
264  << " integrated step " << stepLen << " X0 "
265  << radLen << " Lamda " << intLen;
266 
267  me100[ii]->Fill(eta, radLen);
268  me200[ii]->Fill(eta, intLen);
269  me300[ii]->Fill(eta, stepLen);
270  me400[ii]->Fill(eta);
271 
272  if (etaAbs >= etaLow && etaAbs <= etaHigh) {
273  me500[ii]->Fill(phi, radLen);
274  me600[ii]->Fill(phi, intLen);
275  me700[ii]->Fill(phi, stepLen);
276  me800[ii]->Fill(phi);
277  }
278 
279  me900[ii]->Fill(eta, phi, radLen);
280  me1000[ii]->Fill(eta, phi, intLen);
281  me1100[ii]->Fill(eta, phi, stepLen);
282  me1200[ii]->Fill(eta, phi);
283 
284 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:510
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
MaterialBudgetCastorHistos::matList
std::vector< std::string > matList
Definition: MaterialBudgetCastorHistos.h:39
benchmark_cfg.minEta
tuple minEta
Definition: benchmark_cfg.py:51
MaterialBudgetCastorHistos::MaterialBudgetCastorHistos
MaterialBudgetCastorHistos(const edm::ParameterSet &p)
Definition: MaterialBudgetCastorHistos.cc:15
maxEta
double maxEta
Definition: PFJetBenchmarkAnalyzer.cc:77
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
AlCaRecoCosmics_cfg.name
tuple name
Definition: AlCaRecoCosmics_cfg.py:24
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:47
ExpressReco_HICollisions_FallBack.step
tuple step
Definition: ExpressReco_HICollisions_FallBack.py:1252
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:59
testEve_cfg.level
tuple level
Definition: testEve_cfg.py:34
dir
dbl *** dir
Definition: mlp_gen.cc:35
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
Generated for CMSSW Reference Manual by doxygen 1.8.5