ATLAS Offline Software
Loading...
Searching...
No Matches
MonitoringFile_TGCPostProcess.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4
6//Methods to Process TGC Histograms after merging.
7//Author: Akimasa Ishikawa (akimasa.ishikawa@cern.ch)
8//Date: Aug 2009
10
13
14#include <iostream>
15#include <iomanip>
16#include <algorithm>
17#include <fstream>
18#include <cmath>
19#include <cstdlib>
20#include <sstream>
21#include <vector>
22#include <utility>
23#include <map>
24#include <string>
25
26#include "TH1F.h"
27#include "TH2F.h"
28#include "TFile.h"
29#include "TClass.h"
30#include "TKey.h"
31#include "TMath.h"
32#include "TTree.h"
33#include "TBranch.h"
34#include "TGraph.h"
35
36namespace dqutils {
37 static const float TGCChamberLowOccupancyCut = 1.0e-6;
38 static const float TGCChamberHighOccupancyCut = 0.005;
39 static const float TGCChannelOccupancyCut = 0.01;
40 static const float TGCChamberEfficiencyCut = 0.7;
41 static const float TGCChamberTimingCut = 0.95;
42
43 static const bool tgc_debug = false;
44
45 void
46 MonitoringFile::TGCPostProcess(const std::string& inFilename, bool /* isIncremental */) {
47 //start postprocessing
48 std::vector< std::pair< std::string, float > > NoisyChambers;
49 std::vector< std::pair< std::string, float > > DeadChambers;
50
51 std::vector< std::pair< std::string, float > > NoisyChannels;
52
53 std::vector< std::pair< std::string, float > > LowEfficiencyChambers;
54
55 std::vector< std::pair< std::string, float > > ReadoutTimingChambers;
56 std::vector< std::pair< std::string, float > > TriggerTimingChambers;
57
58 TGCLV1HistogramDivision(inFilename);
59 TGCRawHistogramDivision(inFilename);
60 TGCChamberOccupancy(inFilename, NoisyChambers, DeadChambers);
61 // TGCChannelOccupancy( inFilename, run_dir, NoisyChannels );
62
63 TGCChamberEfficiency(inFilename, LowEfficiencyChambers);
64 TGCChamberTiming(inFilename, ReadoutTimingChambers, TriggerTimingChambers);
65 //DQA results to COOL
66
67 return;
68 }//MonitoringFile::TGCPostProcess
69
70 bool
71 MonitoringFile::TGCCheckHistogram(TFile* f, TString& hname) {
72 if (!(f->Get(hname))) {
73 //std::cerr << "TGC PostProcessing: no such histogram!! : "<< hname << std::endl;
74 gDirectory->pwd();
75 gDirectory->ls();
76 return false;
77 }
78 return true;
79 }//MonitoringFile::TGCCheckHistogram
80
81 void
82 MonitoringFile::TGCChamberOccupancy(const std::string& inFilename,
83 std::vector< std::pair< std::string, float > >& noisychambers,
84 std::vector< std::pair< std::string, float > >& deadchambers) {
85 PostProcessorFileWrapper mf(inFilename, "TGC ChamberOcc Calculations");
86
87 if (!mf.IsOpen()) {
88 //std::cerr << "TGCPostProcess(): "
89 // << "Input file not opened \n";
90 return;
91 }
92 if (mf.GetSize() < 1000.) {
93 //std::cerr << "TGCPostProcess(): "
94 // << "Input file empty \n";
95 return;
96 }
97 /*
98 // get run directory name
99 //Seemingly unnecessary lines are necessary
100 TIter nextcd0(gDirectory->GetListOfKeys());
101 TKey *key0 = (TKey*)nextcd0();
102 TDirectory *dir0= dynamic_cast<TDirectory*> (key0->ReadObj());
103 dir0->cd();
105 */
106 TIter next_run(mf.GetListOfKeys());
107 TKey* key_run(0);
108 while ((key_run = dynamic_cast<TKey*> (next_run())) != 0) {
109 if (!key_run->IsFolder()) continue;
110 TString run_dir = key_run->GetName();
111 if (!run_dir.Contains("run")) continue;
112
113 std::string run_dir2 = run_dir.Data();
114 //int run_number = atoi( (run_dir2.substr(4, run_dir2.size()-4 )).c_str() );
115 //run_number=run_number;
116
117 TString tgc_dir = run_dir + "/Muon/MuonRawDataMonitoring/TGC/";
118
119 TString tgc_global_dir = tgc_dir + "Global/";
120
121 TString tgc_sub_dir[2] = {
122 tgc_dir + "TGCEA/", tgc_dir + "TGCEC/"
123 };
124
125 TString sac[2] = {
126 "_A", "_C"
127 };
128 TString side[2] = {
129 "A", "C"
130 };
131 TString sws[2] = {
132 "Wire_", "Strip_"
133 };
134
135 std::stringstream ss;
136
137 //Summary histograms
138 TString schambertypesummary[2][4][6];//[ws][station][eta]
139 TH1F* chambertypesummary[2][4][6];//
140
141 TString schambersummary[2][2];//[ac][ws]
142 TH1F* chambersummary[2][2];//
143
144 std::string type[17] = {
145 "_F_T1", "_E4_T3", "_E3_T6", "_E2_T7", "_E1_T8",
146 "_F_T2", "_E5_T4", "_E4_T6", "_E3_T7", "_E2_T8", "_E1_T9",
147 "_F_T2", "_E5_T5", "_E4_T6", "_E3_T7", "_E2_T8", "_E1_T9"
148 };
149
150 //get summary histograms
151 for (int ws = 0; ws < 2; ws++) {
152 int ntype = 0;
153 for (int station = 0; station < 3; station++) {
154 for (int eta = 0; eta < 6; eta++) {
155 if (eta == 5 && station == 0) continue;
156 if (eta < 2 && station == 3) continue;
157
158 ss.str("");
159 ss << "Station" << station + 1 << type[ntype];
160 schambertypesummary[ws][station][eta] = tgc_dir + "Global/Summary/"
161 + "Summary_Of_Log10_" + sws[ws] + "Occupancy_Per_Chamber_Type_" +
162 ss.str();
163 if (tgc_debug) std::cout << schambertypesummary[ws][station][eta] << std::endl;
164
165 chambertypesummary[ws][station][eta] = 0;
166 mf.get(schambertypesummary[ws][station][eta], chambertypesummary[ws][station][eta]);
167 if (!chambertypesummary[ws][station][eta]) {
168 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< schambertypesummary[ws][station][eta] <<
169 // std::endl;
170 continue;
171 }
172 TGCResetContents(chambertypesummary[ws][station][eta]);
173
174 ntype++;
175 }//eta
176 }//station
177
178 for (int ac = 0; ac < 2; ac++) {
179 schambersummary[ac][ws] = tgc_sub_dir[ac] + "Summary/"
180 + "Summary_Of_Log10_" + sws[ws] + "Occupancy_Per_GasGap" + sac[ac];
181 if (tgc_debug) std::cout << schambersummary[ac][ws] << std::endl;
182
183 chambersummary[ac][ws] = 0;
184 mf.get(schambersummary[ac][ws], chambersummary[ac][ws]);
185 if (!chambersummary[ac][ws]) {
186 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< schambersummary[ac][ws] << std::endl;
187 continue;
188 }
189 TGCResetContents(chambersummary[ac][ws]);
190 }// ac
191 }// ws
192
193 //get number of events processed
194 TString sentries = tgc_global_dir + "Event_Counter";
195
196 TH1F* hentries = 0;
197 mf.get(sentries, hentries);
198 if (!hentries) {
199 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< sentries << std::endl;
200 continue;
201 }
202
203 double nentries = hentries->GetEntries();
204
205 //calculate occupancies
206 for (int ac = 0; ac < 2; ac++) {
207 for (int ws = 0; ws < 2; ws++) {
208 double min = 1.;
209 double max = 0.;
210
211 TString sprof = tgc_sub_dir[ac] + "Profile/" + sws[ws] + "Profile_Map" + sac[ac];
212 TString soccu = tgc_sub_dir[ac] + "Occupancy/" + sws[ws] + "Occupancy_Map" + sac[ac];
213
214 TH2F* hprof = 0;
215 mf.get(sprof, hprof);
216 if (!hprof) {
217 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< sprof << std::endl;
218 continue;
219 }
220
221 TH2F* hoccu = 0;
222 mf.get(soccu, hoccu);
223 if (!hoccu) {
224 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< soccu << std::endl;
225 continue;
226 }
227 TGCResetContents(hoccu);
228
229 int nbinx = hprof->GetXaxis()->GetNbins();
230 int nbiny = hprof->GetYaxis()->GetNbins();
231
232 for (int binx = 1; binx <= nbinx; binx++) {
233 for (int biny = 1; biny <= nbiny; biny++) {
234 if (ws == 1 && (binx == 2 || binx == 9 || binx == 16 || binx == 23 || binx == 34)) continue; //skip strip
235 // for layer2
236 if ((binx >= 33 && binx <= 43) &&
237 biny % 2 == 0) continue; // skip phi1,3 for forward region
238 if ((binx == 40 || binx == 41) &&
239 (biny == 19 || biny == 35 || biny == 43)) continue; //skip phi9,17,21 for EI
240
241 double num = hprof->GetBinContent(binx, biny);
242 double den = nentries * nTGCWireStripMap(ws, binx - 1, biny);
243 double occu = 0.;
244 double eoccu = 0.;
245 if (num > 0 && nentries > 0) {
246 occu = (double) num / den;
247 eoccu = sqrt(occu * (1. - occu) / den);
248 }
249 hoccu->SetBinContent(binx, biny, occu);
250 hoccu->SetBinError(binx, biny, eoccu);
251
252 // to set minimum and maximum
253 if (occu < min && occu > 1e-8) min = occu;
254 if (occu > max) max = occu;
255
256 if (occu < 1.e-8) {
257 occu = 1.e-8;
258 eoccu = 1.e-8;
259 }
260 if (chambersummary[ac][ws]) chambersummary[ac][ws]->Fill(log10(occu));
261
262 int station = 0;
263 int eta = -1;
264 if (binx < 16) {
265 station = 0;
266 eta = (binx - 1) % 5;
267 } else if (binx < 28) {
268 station = 1;
269 eta = (binx - 16) % 6;
270 } else if (binx < 40) {
271 station = 2;
272 eta = (binx - 28) % 6;
273 } else {
274 station = 3;
275 eta = (binx - 40) % 2;
276 }
277 if (station < 3 && eta >= 0)
278 if (chambertypesummary[ws][station][eta]) chambertypesummary[ws][station][eta]->Fill(log10(occu));
279
280 //high occupancy chamber
281 if (occu > TGCChamberHighOccupancyCut) {
282 ss.str("");
283 TString schamber = hprof->GetXaxis()->GetBinLabel(binx);
284 int sector = (biny - 1) / 4 + 1;
285 int phi = (biny - 1) % 4;
286 ss << side[ac];
287 if (sector < 10) ss << "0";
288 ss << sector << "phi" << phi;
289 if (tgc_debug) std::cout << "TGC noisy channel : " << ss.str() << " " << occu << std::endl;
290
291 std::pair<std::string, float> p(ss.str(), occu);
292 noisychambers.push_back(std::move(p));
293 } else if (occu < TGCChamberLowOccupancyCut) {//too low occupancy
294 ss.str("");
295 TString schamber = hprof->GetXaxis()->GetBinLabel(binx);
296 int sector = (biny - 1) / 4 + 1;
297 int phi = (biny - 1) % 4;
298 ss << side[ac];
299 if (sector < 10) ss << "0";
300 ss << sector << "phi" << phi;
301 if (tgc_debug) std::cout << "TGC low occupancy chamber : " << ss.str() << " " << occu << std::endl;
302
303 std::pair<std::string, float> p(ss.str(), occu);
304 deadchambers.push_back(std::move(p));
305 }
306 }//biny
307 }//binx
308
309 hoccu->SetMinimum(min * 0.95);
310 hoccu->SetMaximum(max * 1.05);
311
312 TString occu_dir = tgc_sub_dir[ac] + "Occupancy/";
313 TDirectory* dir = mf.GetDirectory(occu_dir);
314
315 if (dir) {
316 dir->cd();
317 hoccu->Write("", TObject::kOverwrite);
318 } else {
319 //std::cerr<<"TGCChamberOccupancy: directory "<<occu_dir<<" not found"<<std::endl;
320 }
321 }//ws
322 }//ac
323 mf.Write();
324
325 // if directory is found, save the summary histogram
326 TString sum_dir = tgc_dir + "Global/Summary";
327 TDirectory* dir = mf.GetDirectory(sum_dir);
328
329 if (dir) {
330 dir->cd();
331 for (int ws = 0; ws < 2; ws++) {
332 for (int station = 0; station < 3; station++) {
333 for (int eta = 0; eta < 6; eta++) {
334 if (eta == 5 && station == 0) continue;
335 if (eta < 2 && station == 3) continue;
336 if (chambertypesummary[ws][station][eta]) chambertypesummary[ws][station][eta]->Write("",
337 TObject::kOverwrite);
338
339 }//eta
340 }//station
341 }//ws
342 } else {
343 //std::cerr<<"TGCChamberOccupancy: directory "<<sum_dir<<" not found"<<std::endl;
344 }//if
345 mf.Write();
346
347 for (int ac = 0; ac < 2; ac++) {
348 // if directory is found, save the summary histogram
349 sum_dir = tgc_sub_dir[ac] + "Summary/";
350 dir = mf.GetDirectory(sum_dir);
351 if (dir) {
352 dir->cd();
353 for (int ws = 0; ws < 2; ws++) {
354 if (chambersummary[ac][ws]) chambersummary[ac][ws]->Write("", TObject::kOverwrite);
355 }//ws
356 } else {
357 //std::cerr<<"TGCChamberOccupancy: directory "<<sum_dir<<" not found"<<std::endl;
358 }//if
359 }//ac
360 mf.Write();
361 }//while
362
363 mf.Write();
364 }//MonitoringFile::TGCChamberOccupancy
365
366 void
367 MonitoringFile::TGCChamberEfficiency(const std::string& inFilename,
368 std::vector< std::pair< std::string,
369 float > >& loweffchambers) {
370 PostProcessorFileWrapper mf(inFilename, "TGC ChamberEff Calculations");
371
372 if (!mf.IsOpen()) {
373 //std::cerr << "TGCPostProcess(): "
374 // << "Input file not opened \n";
375 return;
376 }
377 if (mf.GetSize() < 1000.) {
378 //std::cerr << "TGCPostProcess(): "
379 // << "Input file empty \n";
380 return;
381 }
382 /*
383 // get run directory name
384 //Seemingly unnecessary lines are necessary
385 TIter nextcd0(gDirectory->GetListOfKeys());
386 TKey *key0 = (TKey*)nextcd0();
387 TDirectory *dir0= dynamic_cast<TDirectory*> (key0->ReadObj());
388 dir0->cd();
390 */
391 TIter next_run(mf.GetListOfKeys());
392 TKey* key_run(0);
393 while ((key_run = dynamic_cast<TKey*> (next_run())) != 0) {
394 if (!key_run->IsFolder()) continue;
395 TString run_dir = key_run->GetName();
396 if (!run_dir.Contains("run")) continue;
397
398 std::string run_dir2 = run_dir.Data();
399 //int run_number = atoi( (run_dir2.substr(4, run_dir2.size()-4 )).c_str() );
400 //run_number=run_number;
401
402 TString tgc_dir = run_dir + "/Muon/MuonRawDataMonitoring/TGC/";
403
404 TString tgc_sub_dir[2] = {
405 tgc_dir + "TGCEA/", tgc_dir + "TGCEC/"
406 };
407
408 TString sac[2] = {
409 "_A", "_C"
410 };
411 TString sws[2] = {
412 "Wire_", "Strip_"
413 };
414
415 std::stringstream ss;
416
417 //Summary histograms
418 TString schambertypesummary[2][3][6];//[ws][station][eta]
419 TH1F* chambertypesummary[2][3][6];//
420
421 TString schambersummary[2][2];//[ac][ws]
422 TH1F* chambersummary[2][2];//
423
424 std::string type[17] = {
425 "_F_T1", "_E4_T3", "_E3_T6", "_E2_T7", "_E1_T8",
426 "_F_T2", "_E5_T4", "_E4_T6", "_E3_T7", "_E2_T8", "_E1_T9",
427 "_F_T2", "_E5_T5", "_E4_T6", "_E3_T7", "_E2_T8", "_E1_T9"
428 };
429
430 //get summary histograms
431 for (int ws = 0; ws < 2; ws++) {
432 int ntype = 0;
433 for (int station = 0; station < 3; station++) {
434 for (int eta = 0; eta < 6; eta++) {
435 if (eta == 5 && station == 0) continue;
436 if (eta < 2 && station == 3) continue;
437
438 ss.str("");
439 ss << "Station" << station + 1 << type[ntype];
440 schambertypesummary[ws][station][eta] = tgc_dir + "Global/Summary/"
441 + "Summary_Of_" + sws[ws] + "Efficiency_Per_Chamber_Type_" +
442 ss.str();
443 if (tgc_debug) std::cout << schambertypesummary[ws][station][eta] << std::endl;
444
445 chambertypesummary[ws][station][eta] = 0;
446 mf.get(schambertypesummary[ws][station][eta], chambertypesummary[ws][station][eta]);
447 if (!chambertypesummary[ws][station][eta]) {
448 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< schambertypesummary[ws][station][eta] <<
449 // std::endl;
450 continue;
451 }
452 TGCResetContents(chambertypesummary[ws][station][eta]);
453
454 ntype++;
455 }//eta
456 }//station
457
458 for (int ac = 0; ac < 2; ac++) {
459 schambersummary[ac][ws] = tgc_sub_dir[ac] + "Summary/"
460 + "Summary_Of_" + sws[ws] + "Efficiency_Per_GasGap" + sac[ac];
461 if (tgc_debug) std::cout << schambersummary[ac][ws] << std::endl;
462
463 chambersummary[ac][ws] = 0;
464 mf.get(schambersummary[ac][ws], chambersummary[ac][ws]);
465 if (!chambersummary[ac][ws]) {
466 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< schambersummary[ac][ws] << std::endl;
467 continue;
468 }
469 TGCResetContents(chambersummary[ac][ws]);
470 }//ac
471 }//ws
472
473 for (int ac = 0; ac < 2; ac++) {
474 for (int ws = 0; ws < 2; ws++) {
475 TString seff = tgc_sub_dir[ac] + "Efficiency/" + sws[ws] + "Efficiency_Map" + sac[ac];
476 TString sden = tgc_sub_dir[ac] + "Efficiency/NumDenom/" + sws[ws] + "Efficiency_Map" + sac[ac] +
477 "_Denominator";
478
479 TH2F* heff = 0;
480 mf.get(seff, heff);
481 if (!heff) {
482 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< seff << std::endl;
483 continue;
484 }
485
486 TH2F* hden = 0;
487 mf.get(sden, hden);
488 if (!hden) {
489 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< sden << std::endl;
490 }
491
492 int nbinx = heff->GetXaxis()->GetNbins();
493 int nbiny = heff->GetYaxis()->GetNbins();
494
495 for (int binx = 1; binx <= nbinx; binx++) {
496 for (int biny = 1; biny <= nbiny; biny++) {
497 //if( ( binx==5 || binx==10 || binx==15 ||
498 // binx==21 || binx==27 || binx==33 || binx>=39 ) &&
499 // biny%2 == 0 )continue;// skip phi1,3 for forward region
500 //if( ( binx ==40 || binx==42 ) &&
501 // ( biny ==18 || biny==34 || biny==42 ) ) continue;//skip phi9,17,21 for EI
502 if ((binx >= 33 && binx <= 43) &&
503 biny % 2 == 0) continue; // skip phi1,3 for forward region
504 if ((binx == 40 || binx == 41) &&
505 (biny == 19 || biny == 35 || biny == 43)) continue; //skip phi9,17,21 for EI
506
507
508 double eff = heff->GetBinContent(binx, biny);
509 double denom = 1;
510 if (hden) {
511 denom = hden->GetBinContent(binx, biny);
512 }
513
514 if (chambersummary[ac][ws])
515 if (denom > 0) chambersummary[ac][ws]->Fill(eff);
516
517 int station = 0;
518 int eta = -1;
519 if (binx < 16) {
520 station = 0;
521 eta = (binx - 1) % 5;
522 } else if (binx < 28) {
523 station = 1;
524 eta = (binx - 16) % 6;
525 } else if (binx < 40) {
526 station = 2;
527 eta = (binx - 28) % 6;
528 } else {
529 station = 3;
530 eta = (binx - 40) % 2;
531 }
532 if (station < 3 && eta >= 0)
533 if (chambertypesummary[ws][station][eta])
534 if (denom > 0) chambertypesummary[ws][station][eta]->Fill(eff);
535
536 //low efficiency chambers
537 if (eff < TGCChamberEfficiencyCut) {
538 ss.str("");
539 TString schamber = heff->GetXaxis()->GetBinLabel(binx);
540 int sector = (biny - 1) / 4 + 1;
541 int phi = (biny - 1) % 4;
542 ss << sac[ac];
543 if (sector < 10) ss << "0";
544 ss << sector << "phi" << phi;
545
546 if (tgc_debug) std::cout << "TGC low efficiency chamber : " << ss.str() << " " << eff << std::endl;
547
548 std::pair<std::string, float> p(ss.str(), eff);
549 loweffchambers.push_back(std::move(p));
550 }
551 }//biny
552 }//binx
553
554 TString eff_dir = tgc_sub_dir[ac] + "Efficiency/";
555 TDirectory* dir = mf.GetDirectory(eff_dir);
556
557 if (dir) {
558 dir->cd();
559 heff->Write("", TObject::kOverwrite);
560 } else {
561 //std::cerr<<"TGCChamberEfficiency: directory "<<eff_dir<<" not found"<<std::endl;
562 }
563 }//ws
564 }//ac
565
566 mf.Write();
567 // if directory is found, save the summary histogram
568 TString sum_dir = tgc_dir + "Global/Summary";
569 TDirectory* dir = mf.GetDirectory(sum_dir);
570
571
572 if (dir) {
573 dir->cd();
574 for (int ws = 0; ws < 2; ws++) {
575 for (int station = 0; station < 3; station++) {
576 for (int eta = 0; eta < 6; eta++) {
577 if (eta == 5 && station == 0) continue;
578 if (eta < 2 && station == 3) continue;
579 if (chambertypesummary[ws][station][eta]) chambertypesummary[ws][station][eta]->Write("",
580 TObject::kOverwrite);
581
582 }//eta
583 }//station
584 }//ws
585 } else {
586 //std::cerr<<"TGCChamberEfficiency: directory "<<sum_dir<<" not found"<<std::endl;
587 }//if
588 mf.Write();
589
590 for (int ac = 0; ac < 2; ac++) {
591 // if directory is found, save the summary histogram
592 sum_dir = tgc_sub_dir[ac] + "Summary/";
593 dir = mf.GetDirectory(sum_dir);
594 if (dir) {
595 dir->cd();
596 for (int ws = 0; ws < 2; ws++) {
597 if (chambersummary[ac][ws]) chambersummary[ac][ws]->Write("", TObject::kOverwrite);
598 }//ws
599 } else {
600 //std::cerr<<"TGCChamberEfficiency: directory "<<sum_dir<<" not found"<<std::endl;
601 }//if
602 }//ac
603 mf.Write();
604 }//while
605 mf.Write();
606 }//MonitoringFile::TGCChamberEfficiency
607
608 void
609 MonitoringFile::TGCChamberTiming(const std::string& inFilename,
610 std::vector< std::pair< std::string, float > >& badrotimingchambers,
611 std::vector< std::pair< std::string, float > >& badtrgtimingchambers) {
612 PostProcessorFileWrapper mf(inFilename, "TGC Chamber Timing Calculations");
613
614 if (!mf.IsOpen()) {
615 //std::cerr << "TGCPostProcess(): "
616 // << "Input file not opened \n";
617 return;
618 }
619 if (mf.GetSize() < 1000.) {
620 //std::cerr << "TGCPostProcess(): "
621 // << "Input file empty \n";
622 return;
623 }
624 /*
625 // get run directory name
626 //Seemingly unnecessary lines are necessary
627 TIter nextcd0(gDirectory->GetListOfKeys());
628 TKey *key0 = (TKey*)nextcd0();
629 TDirectory *dir0= dynamic_cast<TDirectory*> (key0->ReadObj());
630 dir0->cd();
632 */
633 TIter next_run(mf.GetListOfKeys());
634 TKey* key_run(0);
635 while ((key_run = dynamic_cast<TKey*> (next_run())) != 0) {
636 if (!key_run->IsFolder()) continue;
637 TString run_dir = key_run->GetName();
638 if (!run_dir.Contains("run")) continue;
639
640 std::string run_dir2 = run_dir.Data();
641 //int run_number = atoi( (run_dir2.substr(4, run_dir2.size()-4 )).c_str() );
642 //run_number=run_number;
643
644 TString tgc_dir = run_dir + "/Muon/MuonRawDataMonitoring/TGCLV1/";
645
646 TString tgc_sub_dir[2] = {
647 tgc_dir + "TGCEA/", tgc_dir + "TGCEC/"
648 };
649
650 TString sac[2] = {
651 "_A", "_C"
652 };
653
654 TString ssllpt[2] = {
655 "SL_Timing", "Low_Pt_Timing"
656 };
657
658 std::stringstream ss;
659 //Summary histograms
660 TString schambertypesummary[6];//[station][eta]
661 TH1F* chambertypesummary[6];//
662
663 std::string type[6] = {
664 "_F_T2", "_E5_T5", "_E4_T6", "_E3_T7", "_E2_T8", "_E1_T9"
665 };
666
667 for (int sllpt = 0; sllpt < 2; sllpt++) {//loop over SL/LPT
668 for (int eta = 0; eta < 6; eta++) {
669 //get summary histograms
670 ss.str("");
671 ss << type[eta];
672 schambertypesummary[eta] = tgc_dir + "Global/Summary/"
673 + "Summary_Of_" + ssllpt[sllpt] + "_Per_Chamber_Type" + ss.str();
674
675 chambertypesummary[eta] = 0;
676 mf.get(schambertypesummary[eta], chambertypesummary[eta]);
677 if (!chambertypesummary[eta]) {
678 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< schambertypesummary[eta] << std::endl;
679 continue;
680 }
681 TGCResetContents(chambertypesummary[eta]);
682 }//eta
683
684 for (int ac = 0; ac < 2; ac++) {
685 //get summary histogram
686 TString ssum = tgc_sub_dir[ac] + "Summary/Summary_Of_" + ssllpt[sllpt] + sac[ac];
687 TString stmap = tgc_sub_dir[ac] + ssllpt[sllpt] + "_Map" + sac[ac];
688 TString stfrac = tgc_sub_dir[ac] + ssllpt[sllpt] + "_Fraction_Map" + sac[ac];
689
690 TH1F* hsum = 0;
691 mf.get(ssum, hsum);
692 if (!hsum) {
693 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< ssum << std::endl;
694 continue;
695 }
696
697 TH2F* htmap = 0;
698 mf.get(stmap, htmap);
699 if (!htmap) {
700 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< stmap << std::endl;
701 continue;
702 }
703
704 TH2F* htfrac = 0;
705 mf.get(stfrac, htfrac);
706 if (!htfrac) {
707 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< stfrac << std::endl;
708 continue;
709 }
710
712 TGCResetContents(htfrac);
713
714 int nbinx = htmap->GetXaxis()->GetNbins();
715 int nbiny = htmap->GetYaxis()->GetNbins();
716
717 if (nbinx != 18 || nbiny != 48) continue;
718
719 for (int eta = 0; eta < 6; eta++) {
720 for (int phi48 = 0; phi48 < 48; phi48++) {
721 if (eta == 5 && phi48 % 2 == 1) continue; //Forward
722
723 float p = htmap->GetBinContent(eta + 1, phi48 + 1);
724 float c = htmap->GetBinContent(eta + 7, phi48 + 1);
725 float n = htmap->GetBinContent(eta + 13, phi48 + 1);
726 float tot = p + n + c;
727
728 float fp = 0;
729 float fc = 0;
730 float fn = 0;
731
732 float efp = 0;
733 float efc = 0;
734 float efn = 0;
735
736 if (tot != 0) {
737 fp = p / tot;
738 efp = sqrt(fp * (1. - fp) / tot);
739 fc = c / tot;
740 efc = sqrt(fc * (1. - fc) / tot);
741 fn = n / tot;
742 efn = sqrt(fn * (1. - fn) / tot);
743 }
744
745 htfrac->SetBinContent(eta + 1, phi48 + 1, fp);
746 htfrac->SetBinContent(eta + 7, phi48 + 1, fc);
747 htfrac->SetBinContent(eta + 13, phi48 + 1, fn);
748 htfrac->SetBinError(eta + 1, phi48 + 1, efp);
749 htfrac->SetBinError(eta + 7, phi48 + 1, efc);
750 htfrac->SetBinError(eta + 13, phi48 + 1, efn);
751
752 if (chambertypesummary[eta]) chambertypesummary[eta]->Fill(fc);
753 hsum->Fill(fc);
754
755 //low current fraction chambers
756 if (fc < TGCChamberTimingCut) {
757 ss.str("");
758 //std::string chamber = htiming->GetXaxis()->GetBinLabel(bin);
759 ss << sac[ac] << "phi" << phi48 + 1 << "Eta" << eta;
760
761 if (tgc_debug) std::cout << "TGC " << ssllpt[sllpt] << " low current BCID fraction chamber : " <<
762 ss.str() << " " << fc << std::endl;
763
764 std::pair<std::string, float> p(ss.str(), fc);
765 if (sllpt == 0) {
766 badtrgtimingchambers.push_back(std::move(p));
767 } else {
768 badrotimingchambers.push_back(std::move(p));
769 }
770 }
771 }//phi48
772 }//eta
773
774 //write timing fraction
775 TString timing_dir = tgc_sub_dir[ac];
776 TDirectory* dir = mf.GetDirectory(timing_dir);
777
778 if (dir) {
779 dir->cd();
780 htfrac->Write("", TObject::kOverwrite);
781 } else {
782 //std::cerr<<"TGCChamberTiming : directory "<<timing_dir<<" not found"<<std::endl;
783 }
784
785 //write summary of timing for each side
786 TString sum_dir = tgc_sub_dir[ac] + "Summary/";
787 dir = mf.GetDirectory(sum_dir);
788
789 if (dir) {
790 dir->cd();
791 hsum->Write("", TObject::kOverwrite);
792 } else {
793 //std::cerr<<"TGCChamberTiming : directory "<<sum_dir<<" not found"<<std::endl;
794 }
795 }//ac
796 mf.Write();
797
798 //write summary of timing for each eta
799 TString sum_dir = tgc_dir + "Global/Summary/";
800 TDirectory* dir = mf.GetDirectory(sum_dir);
801
802 if (dir) {
803 dir->cd();
804 for (int eta = 0; eta < 6; eta++) {
805 if (chambertypesummary[eta]) chambertypesummary[eta]->Write("", TObject::kOverwrite);
806 }
807 } else {
808 //std::cerr<<"TGCChamberTiming : directory "<<sum_dir<<" not found"<<std::endl;
809 }
810 mf.Write();
811 }//sllpt
812 }//while
813 mf.Write();
814 }//MonitoringFile::TGCChamberTiming
815
816 void
817 MonitoringFile::TGCChannelOccupancy(const std::string& inFilename,
818 std::vector< std::pair< std::string,
819 float > >& noisychannels) {
820 PostProcessorFileWrapper mf(inFilename, "TGC ChannelOcc Calculations");
821
822 if (!mf.IsOpen()) {
823 //std::cerr << "TGCPostProcess(): "
824 // << "Input file not opened \n";
825 return;
826 }
827 if (mf.GetSize() < 1000.) {
828 //std::cerr << "TGCPostProcess(): "
829 // << "Input file empty \n";
830 return;
831 }
832 /*
833 // get run directory name
834 //Seemingly unnecessary lines are necessary
835 TIter nextcd0(gDirectory->GetListOfKeys());
836 TKey *key0 = (TKey*)nextcd0();
837 TDirectory *dir0= dynamic_cast<TDirectory*> (key0->ReadObj());
838 dir0->cd();
840 */
841 TIter next_run(mf.GetListOfKeys());
842 TKey* key_run(0);
843 while ((key_run = dynamic_cast<TKey*> (next_run())) != 0) {
844 if (!key_run->IsFolder()) continue;
845 TString run_dir = key_run->GetName();
846 if (!run_dir.Contains("run")) continue;
847
848 std::string run_dir2 = run_dir.Data();
849 //int run_number = atoi( (run_dir2.substr(4, run_dir2.size()-4 )).c_str() );
850 //run_number=run_number;
851
852 TString tgc_dir = run_dir + "/Muon/MuonRawDataMonitoring/TGC/";
853
854 TString tgc_global_dir = tgc_dir + "Global/";
855 TString tgc_sub_dir[2] = {
856 tgc_dir + "TGCEA/", tgc_dir + "TGCEC/"
857 };
858
859 TString sac[2] = {
860 "_A", "_C"
861 };
862 TString sws[2] = {
863 "Wire_", "Strip_"
864 };
865 TString slay[7] = {
866 "1", "2", "3", "4", "5", "6", "7"
867 };
868 std::stringstream ss;
869
870 //get number of events processed
871 TString sentries = tgc_global_dir + "Event_Counter";
872
873 TH1F* hentries = 0;
874 mf.get(sentries, hentries);
875 if (!hentries) {
876 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< sentries << std::endl;
877 continue;
878 }
879
880 double nentries = hentries->GetEntries();
881
882 for (int ac = 0; ac < 2; ac++) {
883 for (int ws = 0; ws < 2; ws++) {
884 for (int lay = 0; lay < 7; lay++) {
885 if (ws == 1 && lay == 1) continue;
886 for (int subsect = 1; subsect <= 48; subsect++) {
887 //get each subsector histogram
888 int sector, phi4;
889 TGCsubsect2sectorphi(subsect, sector, phi4);
890
891 ss.str("");
892 if (sector < 10) ss << "0";
893 ss << sector << "_Layer" << slay[lay] << "_Phi" << phi4;
894 TString sprof = tgc_sub_dir[ac] + "Profile/" + sws[ws] + "Hit_Profile" + sac[ac] + ss.str();
895
896 TH1F* hprof = 0;
897 mf.get(sprof, hprof);
898 if (!hprof) {
899 //std::cerr <<"TGC PostProcessing: no such histogram!! "<< sprof << std::endl;
900 continue;
901 }
902
903 int nbin = hprof->GetXaxis()->GetNbins();
904
905 for (int bin = 1; bin <= nbin; bin++) {
906 double num = hprof->GetBinContent(bin);
907 float occu = 0.;
908 if (nentries > 0) occu = (float) num / nentries;
909 if (occu > TGCChannelOccupancyCut) {
910 ss.str("");
911 ss << hprof->GetName() << " " << bin;
912 if (tgc_debug) std::cout << "TGC noisy channel : " << ss.str() << " " << occu << std::endl;
913 noisychannels.emplace_back(ss.str(), occu);
914 }
915 }//bin
916 }//subsect
917 }//lay
918 }//loop over ws
919 }//loop over ac
920 }//while
921 }//MonitoringFile::TGCChannelOccupancy
922
923 void
924 MonitoringFile::TGCsubsectbin2stationeta(int subsect, int bin, int& station, int& eta) {
925 if (subsect % 2 == 1) {//including forward chamber
926 if (bin <= 15) {
927 station = 1;
928 eta = (bin - 1) % 5;
929 } else if (bin <= 27) {
930 station = 2;
931 eta = (bin - 16) % 6;
932 } else {
933 station = 3;
934 eta = (bin - 28) % 6;
935 }
936 } else {
937 if (bin <= 12) {
938 station = 1;
939 eta = (bin - 1) % 4 + 1;
940 } else if (bin <= 22) {
941 station = 2;
942 eta = (bin - 13) % 5 + 1;
943 } else {
944 station = 3;
945 eta = (bin - 23) % 5 + 1;
946 }
947 }
948 }//MonitoringFile::TGCsubsectbin2stationeta
949
950 void
951 MonitoringFile::TGCsubsect2sectorphi(int subsect, int& sector, int& phi4) {
952 sector = (subsect + 1) / 4 + 1;
953 if (sector > 12) sector = 1;
954 phi4 = (subsect + 1) % 4;
955 }//MonitoringFile::TGCsubsect2sectorphi
956
957 int
958 MonitoringFile::nTGCWireStripMap(int ws, int etac, int phi48) {//[0:1][0:42][1:48]
959 int layer = 0;
960 int eta = etac;
961 int st = 0;
962// if ( etac <= 4 ){ layer = 0; eta = 4 - etac; st = 42;}
963// else if ( etac <= 9 ){ layer = 1; eta = 9 - etac; st = 42;}
964// else if ( etac <= 14 ){ layer = 2; eta = 14 - etac; st = 42;}
965// else if ( etac <= 20 ){ layer = 3; eta = 20 - etac; st = 44;}
966// else if ( etac <= 26 ){ layer = 4; eta = 26 - etac; st = 44;}
967// else if ( etac <= 32 ){ layer = 5; eta = 32 - etac; st = 46;}
968// else if ( etac <= 38 ){ layer = 6; eta = 38 - etac; st = 46;}
969// else if ( etac <= 40 ){ layer = 7; eta = 40 - etac; st = 48;}
970// else if ( etac <= 42 ){ layer = 8; eta = 42 - etac; st = 48;}
971
972 int binx = etac + 1;
973
974 if (binx % 7 >= 1 && binx % 7 <= 3 && binx <= 28) {
975 layer = (binx - 1) % 7;
976 eta = 4 - (binx - 1) / 7;
977 st = 42;
978 } else if (binx % 7 >= 4 && binx % 7 <= 5 && binx <= 28) {
979 layer = (binx - 1) % 7;
980 eta = 5 - (binx - 1) / 7;
981 st = 44;
982 } else if ((binx - 1) % 7 >= 5 && binx <= 28) {
983 layer = (binx - 1) % 7;
984 eta = 5 - (binx - 1) / 7;
985 st = 46;
986 } else if (binx > 28 && binx <= 30) {
987 layer = binx - 28 + 3 - 1;
988 eta = 1;
989 st = 44;
990 } else if (binx > 30 && binx <= 32) {
991 layer = binx - 28 + 3 - 1;
992 eta = 1;
993 st = 46;
994 } else if (binx > 32 && binx <= 35) {
995 layer = binx - 32 - 1;
996 eta = 0;
997 st = 42;
998 } else if (binx > 35 && binx <= 37) {
999 layer = binx - 32 - 1;
1000 eta = 0;
1001 st = 44;
1002 } else if (binx > 37 && binx <= 39) {
1003 layer = binx - 32 - 1;
1004 eta = 0;
1005 st = 46;
1006 } else if (binx > 39 && binx <= 41) {
1007 layer = binx - 32 - 1;
1008 eta = 1;
1009 st = 48;
1010 } else if (binx > 41 && binx <= 43) {
1011 layer = binx - 34 - 1;
1012 eta = 0;
1013 st = 48;
1014 }
1015
1016
1017 if (eta == 0) st -= 1;
1018
1019
1020 int nwire = getTGCNumberOfWires(st, layer, eta, phi48);
1021
1022
1023 //number of strips is always 32 except no chamber region which should be nwire==0.
1024 if (ws == 1) {
1025 if (nwire == 0) return 0;
1026
1027 return 32;
1028 }
1029
1030 return nwire;
1031 }
1032
1033 int
1034 MonitoringFile::getTGCNumberOfWires(const int istationName, const int layer, const int istationEta,
1035 const int istationPhi) {
1036 int nWire = 0;
1037
1038 if (istationName == 42) { //TGC1 Endcap //41:T1F 42:T1E 43:T2F 44:T2E 45:T3F 46:T3E 47:T3F 48:T4E (T4 means inner
1039 // small wheel TGCs, EI/FI)
1040 if (layer == 0) {
1041 if (istationEta == 1) { // T3 (E4) wire 1ch-91ch (total 91ch)
1042 nWire = 92;
1043 } else if (istationEta == 2) { //T6 (E3) 91ch-152ch (62ch)
1044 nWire = 61;
1045 } else if (istationEta == 3) { //T7 (E2) 152ch-174ch (23ch)
1046 nWire = 23;
1047 } else if (istationEta == 4) { //T8 (E1) 173ch-196ch(24ch)
1048 nWire = 24;
1049 }
1050 } else if (layer == 1 || layer == 2) {
1051 if (istationEta == 1) { // T3 (E4) wire 1ch-91ch (total 91ch)
1052 nWire = 91;
1053 } else if (istationEta == 2) { //T6 (E3) 91ch-152ch (62ch)
1054 nWire = 62;
1055 } else if (istationEta == 3) { //T7 (E2) 152ch-174ch (23ch)
1056 nWire = 23;
1057 } else if (istationEta == 4) { //T8 (E1) 173ch-196ch(24ch)
1058 nWire = 24;
1059 }
1060 }
1061 } else if (istationName == 41) { // TGC1 Forward
1062 if (layer == 0 || layer == 2) nWire = 105;
1063 else if (layer == 1) nWire = 104;
1064 } else if (istationName == 44) { // TGC2 Endcap
1065 if (istationEta == 1) { // T4 (E5) wire 1ch-91ch (total 91ch)
1066 nWire = 110;
1067 } else if (istationEta == 2) { //T6 (E4) 109ch-211ch (103ch)
1068 nWire = 103;
1069 } else if (istationEta == 3) { //T7 (E3) 211ch-242ch (23ch)
1070 nWire = 32;
1071 } else if (istationEta == 4) { //T8 (E2) 241ch-272ch (32ch)
1072 nWire = 32;
1073 } else if (istationEta == 5) { //T9 (E1) 271ch-302ch(32ch)
1074 nWire = 32;
1075 }
1076 } else if (istationName == 43) { // TGC2 Forward
1077 nWire = 125;
1078 } else if (istationName == 46) { // TGC3 Endcap
1079 if (istationEta == 1) { // T5 (E5) wire 1ch-110ch (total 110ch)
1080 nWire = 96;
1081 } else if (istationEta == 2) { //T6 (E4) 95ch-200ch (106ch)
1082 nWire = 106;
1083 } else if (istationEta == 3) { //T7 (E3) 200ch-231ch (32ch)
1084 nWire = 32;
1085 } else if (istationEta == 4) { //T8 (E2) 231ch-2
1086 nWire = 30;
1087 } else if (istationEta == 5) { //T9 (E1) 260ch-290ch(31ch)
1088 nWire = 31;
1089 }
1090 } else if (istationName == 45) { // TGC3 Forward
1091 nWire = 122;
1092 } else if (istationName == 48) { // EI
1093 //stationPhi 1-8 -> 1-8
1094 // 9-15 -> 10-16
1095 // 16-18 -> 18-20
1096 // 19-21 -> 22-24
1097 if (istationPhi == 2 || istationPhi == 11 || istationPhi == 13 || istationPhi == 14 || istationPhi == 15 ||
1098 istationPhi == 16 || istationPhi == 20 || istationPhi == 21) {
1099 //short 16
1100 nWire = 16;
1101 } else {
1102 //long 24
1103 nWire = 24;
1104 }
1105 } else if (istationName == 47) { // FI
1106 if (istationPhi == 2 || istationPhi == 5 || istationPhi == 8 || istationPhi == 11 || istationPhi == 14 ||
1107 istationPhi == 17 || istationPhi == 20 || istationPhi == 23) {
1108 //short 30
1109 nWire = 30;
1110 } else {
1111 //long 32
1112 nWire = 32;
1113 }
1114 }
1115 return nWire;
1116 }
1117
1118 void
1120 // copy all functions of histogram before resetting
1121 TList* h_funcs;
1122
1123 h_funcs = dynamic_cast<TList*> (h->GetListOfFunctions()->Clone());
1124 h->Reset();
1125 //now insert these functions back into the hist
1126 TIter iterE(h_funcs);
1127 while ((iterE())) {
1128 h->GetListOfFunctions()->Add(*iterE);
1129 }
1130 delete h_funcs;
1131 // are the functions still valid after list deletion?
1132 //http://root.cern.ch/root/html/TList#TList:_TList should be!
1133 // TIter iterTest(EffBCap->GetListOfFunctions());
1134 // while( iterTest() ) std::cout << (*iterTest)->GetName() << std::endl;
1135 }
1136}//namespace
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
static Double_t ss
#define min(a, b)
Definition cfImp.cxx:40
#define max(a, b)
Definition cfImp.cxx:41
double hsum(TH1 *h)
sum the bin contents including the over and underflow bins
Definition chains.cxx:101
Header file for AthHistogramAlgorithm.
static void TGCLV1HistogramDivision(const std::string &inFilename)
static int getTGCNumberOfWires(const int istationName, const int layer, const int istationEta, const int istationPhi)
static void TGCPostProcess(const std::string &inFilename, bool isIncremental=false)
static bool TGCCheckHistogram(TFile *f, TString &hname)
static void TGCsubsectbin2stationeta(int subsect, int bin, int &station, int &eta)
static void TGCChannelOccupancy(const std::string &inFilename, std::vector< std::pair< std::string, float > > &p)
static void TGCChamberOccupancy(const std::string &inFilename, std::vector< std::pair< std::string, float > > &phigh, std::vector< std::pair< std::string, float > > &plow)
static void TGCChamberTiming(const std::string &inFilename, std::vector< std::pair< std::string, float > > &pro, std::vector< std::pair< std::string, float > > &ptrg)
static void TGCsubsect2sectorphi(int subsect, int &sector, int &phi4)
static void TGCChamberEfficiency(const std::string &inFilename, std::vector< std::pair< std::string, float > > &p)
static void TGCRawHistogramDivision(const std::string &inFilename)
static int nTGCWireStripMap(int ws, int etac, int phi48)
static const float TGCChamberHighOccupancyCut
static const bool tgc_debug
static const float TGCChamberTimingCut
static const float TGCChannelOccupancyCut
static const float TGCChamberEfficiencyCut
static const float TGCChamberLowOccupancyCut