ATLAS Offline Software
Loading...
Searching...
No Matches
LArCellsEmptyMonitoring.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5// LArCellsEmptyMonitoring
6// -- Basic class aimed at running simple LAr monitoring function on ntuples
7// produced from that very same package.
8//
9// -- Author(s): Emma Kuwertz (KTH), Maud Schwoerer (LAPP), Olivier Simard (CEA-Saclay)
10//
11//
12//
13
14// Base includes
15
16// Package includes
17#include "LArSamplesMon/LArCellsEmptyMonitoring.h"
18#include "LArSamplesMon/Data.h"
19#include "LArSamplesMon/Interface.h"
20#include "LArSamplesMon/HistoryIterator.h"
21#include "LArSamplesMon/History.h"
22#include "LArCafJobs/CellInfo.h"
23#include "LArCafJobs/ShapeInfo.h"
24#include "LArCafJobs/EventData.h"
25#include "LArCafJobs/CaloId.h"
26// ROOT includes
27#include "TObject.h"
28#include "TObjString.h"
29#include "TString.h"
30#include "TH1F.h"
31#include "TH2F.h"
32#include "TH2I.h"
33#include "TH3F.h"
34#include "TProfile.h"
35#include "TProfile2D.h"
36#include "TGraph.h"
37#include "TDirectory.h"
38#include "TFile.h"
39#include "TTree.h"
40#include "TCanvas.h"
41#include "TAxis.h"
42#include "TKey.h"
43
44#include <array>
45#include <map>
46#include <cstdio>//fopen, fprintf
47#include <fstream>
48#include <iostream>
49
50using namespace std;
51using namespace LArSamples;
52
53template <class T>
54using THVec=std::vector<std::unique_ptr<T>>;
55
56template <class T, std::size_t n>
57using THArray = std::array<std::unique_ptr<T>, n >;
58
59using TH1Fp = std::unique_ptr<TH1F>;
60using TH2Fp = std::unique_ptr<TH2F>;
61using TProfilep = std::unique_ptr<TProfile>;
62
63LArCellsEmptyMonitoring::LArCellsEmptyMonitoring()
64 : m_LarIdTranslator (std::make_unique<LArIdTranslatorHelper>("LarIdTree.root"))
65{
66 // class constructor - set your class variables to default values here
67
68 m_nexpected = 195072; // entries expected from Calo container
69
70
71 m_SaveRootFile = false;
72 m_SaveTextFile = true;
73 m_Algo = 2;
74 m_SetLumiblockRange = false;
75 m_input_lbmin = 0;
76 m_input_lbmax = 0;
77 m_RemoveBadLB = true;
78 m_Qthreshold = 4000;
79 m_nsigma = 4.;
80 m_Ethreshold = 1.;
81 m_UpperCountThreshold = 250;
82 m_LowerCountThreshold = 25;
83 m_LowerCellEnergyThreshold = 1.5;
84 m_UpperCellEnergyThreshold = 30.;
85 m_SelectRecurringBadCells = false;
86 m_nsigmaHits = 10.;
87 m_SaveRootTree = false;
88 m_SetPartition = false;
89}
90
94//____________________________________________________________________________________________________
95void LArCellsEmptyMonitoring::TestRun(const TString& inputfile)
96{
97
98 TString htitle,hname,textfilename,rootfilename;
99 //define cumulative energy cuts here (in GeV)
100 double Ecum_cut_PS = 20.;
101 double Ecum_cut_calo = 10.;
102
103 // define other constants and counters
104 const double GeV = 1/1000.;
105 int index = 0;
106 double nsigmaHits = 10.;
107 const int threshold = 1.;
108 const int qthreshold = m_Qthreshold;
109 int nskipped=0;
110 int nBadCells=0;
111 int noEvdata=0;
112 unsigned int onlid = 0;
113 int ne=0;
114 int emin=0;
115 int emax=0;
116 int nq=0;
117 int qmin=0;
118 int qmax=0;
119 int nlb=0;
120 int lbmin=0;
121 int lbmax=0;
122 int runNumber=0;
123 double MeanEventEnergy_max = -999;
124 double MeanEventEnergy_min = 999;
125
126 int CellCounter=0;
127
128
129 // Deal with bad LBs
130 // @ TODO: Do we still need this option? Not if we are always reading in the defect list at merging...?
131 vector<int, std::allocator<int> > DQLBList;
132 // placeholder
133 vector<int, std::allocator<int> > BadLBList;
134
135 // find minimum and maximum LB, cell energy, q factor...
136 printf("Getting histogram limits... ");
137 GetLimits_EqLB(inputfile,lbmin,lbmax,emin,emax,qmin,qmax,runNumber,BadLBList);
138 if (m_SetLumiblockRange){
139 lbmin = m_input_lbmin;
140 lbmax = m_input_lbmax;
141 }
142 nlb=lbmax-lbmin;
143 ne=emax-emin;
144 nq=qmax-qmin;
145 printf("Done.\n LBmin: %d LBmax: %d Emin: %d Emax: %d Qmin: %d Qmax: %d\n",lbmin,lbmax,emin,emax,qmin,qmax);
146
147
148 // Declare and open output textfile
149 FILE * pFile=nullptr;
150 if (m_SaveTextFile){
151 if (!m_SetPartition){
152 textfilename.Form("Output/BadCellList_run%d.txt",runNumber);
153 pFile = fopen (textfilename , "w");
154 if (!pFile){
155 printf("Cannot open output text file\n");
156 std::abort();}
157 } else {
158 // placeholder - in case someone wants to look at individual partitions only!
159 if (m_SetPartition){
160 textfilename.Form("Output/BadCellList_%s_run%d.txt",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex),runNumber);
161 pFile = fopen (textfilename , "w");
162 if (!pFile){
163 printf("Cannot open output text file\n");
164 std::abort();}
165 }
166 }
167 }
168 // Open output text file to write column titles
169 if (m_SaveTextFile){
170 fprintf(pFile, "onlid // part // FT // SL // CH // # Events E>%d sigma // Mean Event Energy [GeV] (min // max) // # LBs // # LBs affected // fraction events Q>%d \n",(int)nsigmaHits,qthreshold);
171 }
172
173 int EventCount=0., qcount=0.;
174
175 // Declare root file (LAr investigations - not intended to be produced by default)
176 if (m_SetPartition){
177 rootfilename.Form("Output/BadCells_%s_run%d.root",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex),runNumber);
178 } else {
179 rootfilename.Form("Output/BadCells_run%d.root",runNumber);
180 }
181
182 auto fout = std::make_unique<TFile>(rootfilename,"RECREATE");
183 Float_t noise = 0;
184 Float_t fr_q4k = 0;
185 TH1Fp h1_lb;
186 TH1Fp h1_elb;
187 TH1Fp h1_qlb;
188 TH1Fp h1_e;
189 TH1Fp h1_q;
190 TH2Fp h2_elb;
191
192
193 // -------------------------------------------------------------------------
194 // create individual cell histos
195 // -------------------------------------------------------------------------
196 printf("Creating histograms (nEbins = %d, nQbins = %d, nLBbins = %d)... ",ne,nq,nlb);
197 TH1F hene = TH1F("hEne","",ne,emin,emax); hene.SetXTitle("Energy [GeV]"); hene.SetYTitle("Events");
198 TH1F hqua = TH1F("hQua","",nq/100.,qmin,qmax); hqua.SetXTitle("Quality"); hqua.SetYTitle("Events");
199 TH1F hlb = TH1F("hLB","",nlb,lbmin,lbmax); hlb.SetXTitle("LB"); hlb.SetYTitle("Events");
200 TH2F henelb = TH2F("hEnelb","",nlb,lbmin,lbmax,ne,emin,emax); henelb.SetXTitle("LB"); henelb.SetYTitle("Energy [GeV]");
201 // -------------------------------------------------------------------------
202 printf("Done.\n");
203
204 // -------------------------------------------------------------------------
205 // New test histos:
206 // -------------------------------------------------------------------------
207 printf("Creating cut test histograms...");
208 static constexpr int npl=4;
209 TH2F tempHist = TH2F("tempNEvLB","",1000,0,1000.,1000,0.,1000.);
210 //THArray, defined at file scope, holds unique_ptr
211 THArray<TH1F, npl> hNoise{};
212 THArray<TH1F, npl> hCellsPerLB{};
213 THArray<TH2F, npl> hNEvVsEMean{};
214 THArray<TH2F, npl> hNEvVsECum{};
215 THArray<TH2F, npl> hECumVsEMean{};
216
217 for (int i=0;i<npl;i++){
218 hname.Form("hNEvVsEMean_Layer%i",i);
219 hNEvVsEMean[i].reset((TH2F*)tempHist.Clone(hname));
220 hNEvVsEMean[i]->GetXaxis()->SetTitle("Number of Events (E > 10 #sigma) / LB");
221 hNEvVsEMean[i]->GetYaxis()->SetTitle("E_{mean} per Event / LB [GeV]");
222
223 hname.Form("hNEvVsECum_Layer%i",i);
224 hNEvVsECum[i].reset((TH2F*)tempHist.Clone(hname));
225 hNEvVsECum[i]->GetXaxis()->SetTitle("Number of Events (E > 10 #sigma) / LB");
226 hNEvVsECum[i]->GetYaxis()->SetTitle("E_{cumulative} per Cell / LB [GeV]");
227
228 hname.Form("hECumVsEMean_Layer%i",i);
229 hECumVsEMean[i].reset((TH2F*)tempHist.Clone(hname));
230 hECumVsEMean[i]->GetXaxis()->SetTitle("E_{cumulative} per Cell / LB [GeV]");
231 hECumVsEMean[i]->GetYaxis()->SetTitle("E_{mean} per Event / LB [GeV]");
232
233 hname.Form("hNoise_Layer%i",i);
234 hNoise[i].reset(new TH1F(hname,"",60,0.,30.));
235 hNoise[i]->GetXaxis()->SetTitle("Noise [GeV]");
236 hNoise[i]->GetYaxis()->SetTitle("Number Of Cells");
237
238 hname.Form("hCellsPerLB_Layer%i",i);
239 hCellsPerLB[i].reset((TH1F*)hlb.Clone(hname));
240 hCellsPerLB[i]->GetXaxis()->SetTitle("LB");
241 hNoise[i]->GetYaxis()->SetTitle("Number Of Cells");
242 }
243 printf("Done\n");
244 // -------------------------------------------------------------------------
245
246
247 // Opening file:
248 std::unique_ptr<LArSamples::Interface> tuple (Interface::open(inputfile));
249 printf("Number of events: %u %u\n",tuple->nEvents(),tuple->nChannels());
250 unsigned int nchannels = tuple->nChannels();
251
252 printf("Begin looping over cells...\n");
253 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
254
255 // pass empty cells
256 if(tuple->historySize(ichan)==0){ nskipped++; continue; }
257
258 // process non empty cells
259 const LArSamples::History* hist = tuple->cellHistory(ichan);
260 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
261
262 // index for partition
263 index = m_LarIdTranslator->GetPartitionLayerIndex(cellInfo->calo(),cellInfo->layer());
264 onlid = cellInfo->onlid();
265
266 // if cell contains less hits than the threshold requirement, skip it
267 unsigned int ndigits = hist->nData();
268 if(ndigits<(unsigned int)threshold){ nskipped++; continue; }
269
270 if(m_SaveTextFile || m_SaveRootFile){
271 // set individual cell histos
272 hname.Form("h0x%x_Energy",onlid);
273 h1_e.reset((TH1F*)hene.Clone(hname));
274 hname.Form("h0x%x_Quality",onlid);
275 h1_q.reset((TH1F*)hqua.Clone(hname));
276 hname.Form("h0x%x_LB",onlid);
277 h1_lb.reset((TH1F*)hlb.Clone(hname));
278 hname.Form("h0x%x_hNEvVsEMean",onlid);
279 h1_elb.reset((TH1F*)hlb.Clone(hname));
280 hname.Form("h0x%x_Energy_LB",onlid);
281 h2_elb.reset((TH2F*)henelb.Clone(hname));
282 hname.Form("h0x%x_Quality_LB",onlid);
283 h1_qlb.reset((TH1F*)hlb.Clone(hname));
284 }
285
286 // loop on the events for each cell
287 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
288 const LArSamples::Data* data = hist->data(idigit);
289 double noise = data->noise()*GeV;
290 const LArSamples::EventData& Evdata = data->eventData();
291 int lumiBlock = Evdata.lumiBlock();
292 double energy = data->energy()*GeV;
293
294 if (energy > nsigmaHits*noise){ // E>10sigma
295
296 // Fill individual cell histos
297 if(m_SaveTextFile || m_SaveRootFile){
298 h1_e->Fill(energy);
299 h1_q->Fill(data->quality());
300 h1_lb->Fill(lumiBlock);
301 h1_elb->Fill(lumiBlock, energy);
302 h2_elb->Fill(lumiBlock, energy);
303 h1_qlb->Fill(lumiBlock, data->quality());
304 }
305
306 // count the number of processed events per cell with a bad q-factor
307 if (data->quality() > qthreshold){
308 qcount = qcount+1;
309 }
310
311 // count the number of events processed per cell
312 EventCount++;
313
314 } // end of E>10sigma
315
316 } // end of event loop
317 if(m_SaveTextFile || m_SaveRootFile){
318 hNoise[cellInfo->layer()]->Fill(noise);
319 }
320 CellCounter++;
321 if ( CellCounter%1000 == 0 ) {
322 printf("Scanning Cell No. = %d\n",CellCounter);
323 }
324 MeanEventEnergy_max = -999;
325 MeanEventEnergy_min = 999;
326 bool FlagCell=false;
327 int FlagCount=0;
328 int NFullLB = 0;
329
330 if(m_SaveTextFile || m_SaveRootFile){
331 for (int ilb=1;ilb<=nlb;ilb++){
332 if (h1_lb->GetBinContent(ilb) > 0.){
333 double EtotInLB = h1_elb->GetBinContent(ilb);
334 double NEventsInLB = h1_lb->GetBinContent(ilb);
335 //double MeanQ = h1_qlb->GetBinContent(ilb) / NEventsInLB; // change to fraction..
336 int LB_number = h1_lb->GetBinLowEdge(ilb);
337 NFullLB += 1;
338 double MeanEventEnergy = EtotInLB / NEventsInLB;
339
340
341 if (cellInfo->layer()==0){
342 if (EtotInLB>Ecum_cut_PS){
343 hNEvVsEMean[cellInfo->layer()]->Fill(NEventsInLB,MeanEventEnergy);
344 hNEvVsECum[cellInfo->layer()]->Fill(NEventsInLB,EtotInLB);
345 hECumVsEMean[cellInfo->layer()]->Fill(EtotInLB,MeanEventEnergy);
346 hCellsPerLB[cellInfo->layer()]->Fill(LB_number);
347 }
348 }
349 else if (EtotInLB>Ecum_cut_calo){
350 hNEvVsEMean[cellInfo->layer()]->Fill(NEventsInLB,MeanEventEnergy);
351 hNEvVsECum[cellInfo->layer()]->Fill(NEventsInLB,EtotInLB);
352 hECumVsEMean[cellInfo->layer()]->Fill(EtotInLB,MeanEventEnergy);
353 hCellsPerLB[cellInfo->layer()]->Fill(LB_number);
354 }
355
356 if (EtotInLB>Ecum_cut_PS && cellInfo->layer()==0){
357 FlagCell=true; FlagCount++;
358 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
359 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
360 }
361 else if (EtotInLB>Ecum_cut_calo && cellInfo->layer()==1){
362 FlagCell=true; FlagCount++;
363 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
364 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
365 }
366 else if (EtotInLB>Ecum_cut_calo && cellInfo->layer()==2){
367 FlagCell=true; FlagCount++;
368 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
369 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
370 }
371 else if (EtotInLB>Ecum_cut_calo && cellInfo->layer()==3){
372 FlagCell=true; FlagCount++;
373 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
374 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
375 }
376 }
377 }
378 if (EventCount!=0){
379 fr_q4k = (double)qcount/(double)EventCount;
380 }
381 }
382 // Write bad channel quantities to text file
383 if (m_SaveTextFile){
384 if (FlagCell){
385 nBadCells++;
386 fprintf(pFile, "0x%x \t %7s \t %d \t %d \t %d \t %d \t %4.3f \t %4.3f \t %d \t %d \t %4.3f \n",onlid,m_LarIdTranslator->GetPartitonLayerName(index),cellInfo->feedThrough(),cellInfo->slot(),cellInfo->channel(),EventCount,MeanEventEnergy_min,MeanEventEnergy_max,NFullLB,FlagCount,fr_q4k);
387
388 }
389 }
390
391 qcount = 0.;
392 EventCount = 0.;
393 } // end of cell loop
394
395 int nbad = m_nexpected-nskipped;
396 printf("Skipped %d cells. Cells scanned = %d\n",nskipped,nbad);
397 printf("Cells put in Bad Cell List = %d\n",nBadCells);
398 if (noEvdata>0) printf("No Evdata pointer: %d\n",noEvdata);
399
400 // Write test histos to output rootfile
401 if(m_SaveRootFile){
402 fout->cd();
403 for(int i=0;i<npl;i++){
404 hNEvVsEMean[i]->Write();
405 hNEvVsECum[i]->Write();
406 hECumVsEMean[i]->Write();
407 hNoise[i]->Write();
408 hCellsPerLB[i]->Write();
409 }
410 }
411
412 fout->Close();
413 fclose(pFile);
414 return;
415
416}
417
418
419
420//____________________________________________________________________________________________________
421void LArCellsEmptyMonitoring::Run(const TString& inputfile)
422{
423
424 TString htitle,hname,textfilename,rootfilename;
425
426 // partition/layer indexing entrirely defined by LarIdTranslator
427 const int npl = m_LarIdTranslator->GetNpl();
428
429
430 // ----------------------------------------------------
431 // Setter error messages
432 // ----------------------------------------------------
433
434
435 if (m_SelectRecurringBadCells && m_Algo != 2){
436 printf("WARNING: Received instruction to select recurring bad cells only in SetSelectRecurringBadCells(bool flag)\n");
437 printf("-------> Only one algorithm selected for cell selection criteria in SetAlgo(int algoindex)\n");
438 printf("-------> Recurring cell selection overides algorithm specification\n");
439 std::abort();
440 }
441
442 if (m_SetPartition && m_PartitionIndex == -1){
443 printf("ERROR: Partition not specified in SetPartition(bool set_part, std::string partname)\n");
444 std::abort();
445 } else if (m_SetPartition){
446 printf("Running cell selection on %s only.\n",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex));
447 }
448 // ----------------------------------------------------
449 //
450 int selectionFlag = -1;
451 const double GeV = 1/1000.;
452 const double ecut = m_Ethreshold; // GeV
453 int index = 0;
454 double nsigma = m_nsigma;
455 double nsigmaHits = m_nsigmaHits;
456 const int threshold = 10.;
457 const int qthreshold = m_Qthreshold;
458 int nskipped=0;
459 int nBadCells=0;
460 int noEvdata=0;
461 unsigned int onlid = 0;
462 int ne=0;
463 int emin=0;
464 int emax=0;
465 int nq=0;
466 int qmin=0;
467 int qmax=0;
468 int nlb=0;
469 int nlb_corr=0; // count only non-empty LBs
470 int lbmin=0;
471 int lbmax=0;
472 int runNumber=0;
473 int nEvents_E_gt_ecut=0;
474 double EventEnergySum=0.;
475 double EventEnergyMean=0.;
476
477 int nBadLBs=0;
478 double MeanHits=0, rmsHits=0;
479 int iMeanCellHitsSelection=0,iEnergyThresholdSelection=0,iOR_selection=0,iAND_selection=0;
480 int CellCounter=0;
481
482
483 // Deal with bad LBs
484 // @ TODO: Do we still need this option? Not if we are always reading in the defect list at merging...?
485 vector<int, std::allocator<int> > DQLBList;
486 vector<int, std::allocator<int> > BadLBList;
487 if(m_ReadDefectLBList){
488 printf("Reading in DQ Defect LB List...\n");
489 DQLBList=ReadBadLBList(m_LBfile);
490 printf("...Done\n");
491 }
492 if (m_RemoveBadLB){
493 printf("Getting bad LB list...\n");
494 BadLBList=GetBadLBList(inputfile,lbmin,lbmax,nsigma,nlb,DQLBList);
495 printf("...Done.\n");
496 nBadLBs = BadLBList.size();
497 }
498 else if (m_ReadDefectLBList){
499 BadLBList = std::move(DQLBList);
500 nBadLBs = BadLBList.size();
501 }
502
503 // find minimum and maximum LB, cell energy, q factor...
504 printf("Getting histogram limits... ");
505 GetLimits_EqLB(inputfile,lbmin,lbmax,emin,emax,qmin,qmax,runNumber,BadLBList);
506 if (m_SetLumiblockRange){
507 lbmin = m_input_lbmin;
508 lbmax = m_input_lbmax;
509 }
510
511
512 nlb=lbmax-lbmin;
513 ne=emax-emin;
514 nq=qmax-qmin;
515 printf("Done.\n LBmin: %d LBmax: %d Emin: %d Emax: %d Qmin: %d Qmax: %d\n",lbmin,lbmax,emin,emax,qmin,qmax);
516
517 // Declare and open summary textfile
518 FILE * outputFile=nullptr;
519 if (!m_SetPartition){
520 textfilename.Form("Output/BadCellSelection_run%d.txt",runNumber);
521 outputFile = fopen (textfilename , "w");
522 if (!outputFile){
523 std::cout<<("Cannot open output selection summary file " + textfilename)<< "\n";
524 std::abort();}
525 } else {
526 textfilename.Form("BadCellSelection_%s_run%d.txt",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex),runNumber);
527 outputFile = fopen (textfilename , "w");
528 if (!outputFile){
529 std::cout<<("Cannot open output selection summary file " + textfilename)<< "\n";
530 std::abort();}
531 }
532
533
534 // Declare and open output textfile
535 FILE * pFile=nullptr;
536 if (m_SaveTextFile){
537 if (!m_SetPartition){
538 textfilename.Form("Output/BadCellList_run%d.txt",runNumber);
539 pFile = fopen (textfilename , "w");
540 if (!pFile){
541 printf("Cannot open output text file\n");
542 std::abort();}
543 } else {
544 if (m_SetPartition){
545 textfilename.Form("BadCellList_%s_run%d.txt",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex),runNumber);
546 pFile = fopen (textfilename , "w");
547 if (!pFile){
548 printf("Cannot open output text file\n");
549 std::abort();}
550 }
551 }
552 }
553
554
555
556
557printf("Determining mean cell hits...\n ");
558 GetMeanCellHits(inputfile, nlb, lbmin, lbmax, nsigmaHits, BadLBList, MeanHits, rmsHits, nlb_corr);
559printf("Done.\n");
560printf("Set threshold at %4.3f counts per cell for LB range. \n",(MeanHits+(nsigmaHits*rmsHits))*((double)nlb_corr-(double)nBadLBs));
561
562
563 if (m_SaveTextFile){
564 fprintf(pFile, "onlid // part // FT // SL // CH // Events (E>%d sigma) // Events E>1GeV // Mean Energy [GeV] // # LBs // fraction events Q>%d // Flagged in LB // selectionFlag\n",(int)nsigmaHits,qthreshold);
565 }
566
567 int EventCount=0., qcount=0.;
568 // ----------------------------------------------------
569 // Prepare output tree listing bad channels
570 // ----------------------------------------------------
571 if (m_SetPartition){
572 rootfilename.Form("Output/BadCells_%s_run%d.root",m_LarIdTranslator->GetPartitonLayerName(m_PartitionIndex),runNumber);
573 } else {
574 rootfilename.Form("Output/BadCells_run%d.root",runNumber);
575 }
576
577 auto fout = std::make_unique<TFile>(rootfilename,"RECREATE");
578 Int_t larid = 0;
579 Int_t t_run =0;
580 Int_t n_ensig = 0;
581 Int_t n_ecut = 0;
582 Float_t fr_LB = 0;
583 Float_t nbrLB = 0;
584 Float_t noise = 0;
585 Float_t meanECell = 0;
586 Float_t fr_q4k = 0;
587 TH1F* h1_lb = nullptr;
588 TH1F* h1_e = nullptr;
589 TH1F* h1_q = nullptr;
590 TH2F* h2_elb = nullptr;
591 TH2D* h2_pulse = nullptr;
592 TH2D* h2_t_LB = nullptr;
593 TProfile* TProf_pulse = nullptr;
594 TH1F* h1_ADCmax = nullptr;
595 std::unique_ptr<TTree> tree_cells;
596 if(m_SaveRootFile){
597 tree_cells.reset(new TTree("BadTree","LAr Tree ordered by OnlID"));
598 tree_cells->Branch("Energy",&h1_e);
599 tree_cells->Branch("Quality",&h1_q);
600 tree_cells->Branch("LB",&h1_lb);
601 tree_cells->Branch("larid",&larid);
602 tree_cells->Branch("onlid",&onlid);
603 tree_cells->Branch("Run",&t_run);
604 tree_cells->Branch("selectionFlag",&selectionFlag);
605 tree_cells->Branch("nhits_ensig",&n_ensig);
606 tree_cells->Branch("nhits_e1GeV",&n_ecut);
607 tree_cells->Branch("Quality_4k",&fr_q4k);
608 tree_cells->Branch("LBfraction",&fr_LB);
609 tree_cells->Branch("EnergyVsLB",&h2_elb);
610 tree_cells->Branch("PulseShape",&h2_pulse);
611 tree_cells->Branch("PulsePeakTimeVsLB",&h2_t_LB);
612 tree_cells->Branch("noise",&noise);
613 tree_cells->Branch("MeanCellEnergy",&meanECell);
614 tree_cells->Branch("PulseShapeTProf",&TProf_pulse);
615 tree_cells->Branch("ADCmax",&h1_ADCmax);
616 }
617
618 Double_t xmin=-4.9,xmax=4.9;
619 Double_t xw=0.1;
620 TH2D NormPulse("Normalised_pulse","",100,-200,200,400,-10,10);
621 NormPulse.SetXTitle("Time [ns]");
622 NormPulse.SetYTitle("Value [ADC counts] / ADCmax");
623 //
624 //THVec, defined at file scope, holds unique_ptr, defined a filescope
625 THVec<TH2D> Pulsemaps(npl);
626 THVec<TH2> Cellmaps(npl);
627 THVec<TH2F> E_LBmaps(npl);
628 THVec<TH2D> t_LBmaps(npl);
629 THVec<TH1F> CellsFlagged_LB_part(npl);
630 THVec<TH1F> Emean_NbrEvents_part(npl);
631 //
632 TH2F E_LB("E_LB","",nlb,lbmin,lbmax,ne,emin,emax);
633 E_LB.SetXTitle("LB");
634 E_LB.SetYTitle("Energy [GeV]");
635 //
636 TH2D t_LB("t_LB","",nlb,lbmin,lbmax,400,-200,200);
637 t_LB.GetXaxis()->SetTitle("LB");
638 t_LB.GetYaxis()->SetTitle("Time(maxSample) + ofcTime [ns]");
639 //
640 TH1F CF_LB("CF_LB","",nlb,lbmin,lbmax);
641 CF_LB.GetXaxis()->SetTitle("LB");
642 CF_LB.GetYaxis()->SetTitle("Number of cells flagged");
643 //
644 TH1F ME_EV("ME_EV","",1000,0,1000);
645 ME_EV.GetXaxis()->SetTitle("Nbr events >1 GeV");
646 ME_EV.GetYaxis()->SetTitle("Mean Energy [GeV]");
647
648
649
650 if(m_SaveRootFile){
651 for (int ii=0;ii<npl;ii++){
652 Cellmaps[ii].reset(m_LarIdTranslator->GetCaloPartitionLayerMap(ii));
653 Cellmaps[ii]->GetXaxis()->SetTitle("#eta"); Cellmaps[ii]->GetYaxis()->SetTitle("#Phi");
654 hname.Form("%s_PulseShape_%dsigma",m_LarIdTranslator->GetPartitonLayerName(ii),(int)nsigmaHits);
655 Pulsemaps[ii].reset((TH2D*)NormPulse.Clone(hname));
656 hname.Form("%s_EvergyVsLB",m_LarIdTranslator->GetPartitonLayerName(ii));
657 E_LBmaps[ii].reset((TH2F*)E_LB.Clone(hname));
658 hname.Form("%s_PulseTimeVsLB",m_LarIdTranslator->GetPartitonLayerName(ii));
659 t_LBmaps[ii].reset((TH2D*)t_LB.Clone(hname));
660 hname.Form("%s_CellsFlaggedVsLB",m_LarIdTranslator->GetPartitonLayerName(ii));
661 CellsFlagged_LB_part[ii].reset((TH1F*)CF_LB.Clone(hname));
662 // en plus
663 hname.Form("%s_EmeanVsNbrEvents",m_LarIdTranslator->GetPartitonLayerName(ii));
664 Emean_NbrEvents_part[ii].reset((TH1F*)ME_EV.Clone(hname));
665 //
666 }
667 }
668 TH2F Eeta("Eeta","",(int)floor((xmax-xmin)/xw),xmin,xmax,ne,emin,emax);
669 Eeta.GetXaxis()->SetTitle("#eta");
670 Eeta.GetYaxis()->SetTitle("Energy [GeV]");
671 //
672 TH1F CellEnergy("CellEnergy","",ne,emin,emax);
673 CellEnergy.GetXaxis()->SetTitle("Cell Mean Energy [GeV]");
674 CellEnergy.GetYaxis()->SetTitle("Cells");
675 //
676 TH1F LBfrac("LBfrac","",100,0.,1.);
677 LBfrac.GetXaxis()->SetTitle("LB fraction");
678 LBfrac.GetYaxis()->SetTitle("Cells");
679 //
680 TH1F Qfrac("Qfrac","",100.,0.,1.);
681 Qfrac.GetXaxis()->SetTitle("Fraction of Events Q>4000");
682 Qfrac.GetYaxis()->SetTitle("Cells");
683 //
684 TH1F CellsFlagged_LB("CellsFlagged_LB","",nlb,lbmin,lbmax);
685 CellsFlagged_LB.GetXaxis()->SetTitle("LB");
686 CellsFlagged_LB.GetYaxis()->SetTitle("Number of cells flagged");
687 //
688 TH1F Emean_NbrEvents("Emean_NEvents_1GeV","",1000,0,1000);
689 Emean_NbrEvents.GetXaxis()->SetTitle("N events > 1 GeV");
690 Emean_NbrEvents.GetYaxis()->SetTitle("Mean Energy [GeV]");
691 //
692 // -------------------------------------------------------------------------
693 // create individual cell histos
694 // -------------------------------------------------------------------------
695 printf("Creating histograms (nEbins = %d, nQbins = %d, nLBbins = %d)... ",ne,nq,nlb);
696 TH1F hene = TH1F("hEne","",ne,emin,emax); hene.SetXTitle("Energy [GeV]");
697 hene.SetYTitle("Events");
698 TH1F hqua = TH1F("hQua","",nq/100.,qmin,qmax); hqua.SetXTitle("Quality");
699 hqua.SetYTitle("Events");
700 TH1F hlb = TH1F("hLB","",nlb,lbmin,lbmax); hlb.SetXTitle("LB");
701 hlb.SetYTitle("Events");
702 TH2F henelb = TH2F("hEnelb","",nlb,lbmin,lbmax,ne,emin,emax); henelb.SetXTitle("LB");
703 henelb.SetYTitle("Energy [GeV]");
704 TH2D hpulse = TH2D("hPulse","",100,-200,200,400, -10,10); hpulse.SetXTitle("Time [ns]");
705 hpulse.SetYTitle("Value [ADC counts] / ADCmax");
706 TProfile TProfpulse = TProfile("", "",5, 0, 5, "s"); TProfpulse.SetXTitle("Sample Number");
707 TProfpulse.SetYTitle("Value [ADC counts]");
708 TH2D ht_LB= TH2D("ht_LB","",nlb,lbmin,lbmax,400,-200,200); ht_LB.GetXaxis()->SetTitle("LB");
709 ht_LB.GetYaxis()->SetTitle("Time(maxSample) + ofcTime [ns]");
710 TH1F hADCmax = TH1F("","",110,-200,2000); hADCmax.SetXTitle("ADCmax [ADC counts]"); hADCmax.SetYTitle("Events");
711 // -------------------------------------------------------------------------
712 printf("Done.\n");
713
714 // Opening file:
715 std::unique_ptr<LArSamples::Interface> tuple = Interface::open(inputfile);
716 printf("Number of events: %u %u\n",tuple->nEvents(),tuple->nChannels()); //tuple->ShowEvents("energy>0.");
717 unsigned int nchannels = tuple->nChannels();
718
719 fprintf(outputFile,"Onlid \t MeanCellHitBased \t EnergyThresholdBased \t OR \t AND\n");
720 // -------------------------------------------------------------------------
721
722 printf("Begin looping over cells...\n");
723 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
724 // pass empty cells
725 if(tuple->historySize(ichan)==0){ nskipped++; continue; }
726
727
728 // process non empty cells
729 const LArSamples::History* hist = tuple->cellHistory(ichan);
730 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
731 double eta = cellInfo->eta();
732 double phi = cellInfo->phi();
733
734
735 // index for partition
736 index = m_LarIdTranslator->GetPartitionLayerIndex(cellInfo->calo(),cellInfo->layer());
737 onlid = cellInfo->onlid();
738
739 // if only checking a specific partition, skip any cells not associated with that partition
740 if(m_SetPartition){
741 if (!CheckPartition(index)){ nskipped++; continue; }
742 }
743
744 // if cell contains less hits that the threshold requirement, skip it
745 unsigned int ndigits = hist->nData();
746 if(ndigits<(unsigned int)threshold){ nskipped++; continue; }
747
748
749 // initialise individual cell histograms
750 if(m_SaveTextFile || m_SaveRootFile){
751 // set individual cell histos
752 hname.Form("h0x%x_Energy",onlid);
753 h1_e = (TH1F*)hene.Clone(hname);
754 hname.Form("h0x%x_Quality",onlid);
755 h1_q = (TH1F*)hqua.Clone(hname);
756 hname.Form("h0x%x_LB",onlid);
757 h1_lb = (TH1F*)hlb.Clone(hname);
758 hname.Form("h0x%x_Energy_LB",onlid);
759 h2_elb = (TH2F*)henelb.Clone(hname);
760 hname.Form("h0x%x_Pulse",onlid);
761 h2_pulse = (TH2D*)hpulse.Clone(hname);
762 hname.Form("TProf0x%x_Pulse",onlid);
763 TProf_pulse = (TProfile*)TProfpulse.Clone(hname);
764 hname.Form("h0x%x_ADCmax",onlid);
765 h1_ADCmax = (TH1F*)hADCmax.Clone(hname);
766 hname.Form("h0x%x_t_LB",onlid);
767 h2_t_LB = (TH2D*)ht_LB.Clone(hname);
768 }
769
770 bool CellInList = false;
771 int LBFlaggedIn = -1;
772
773// loop on the events for each cell
774 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
775 const LArSamples::Data* data = hist->data(idigit);
776 noise = data->noise()*GeV;
777 const LArSamples::EventData& Evdata = data->eventData();
778 int lumiBlock = Evdata.lumiBlock();
779 double energy = data->energy()*GeV;
780
781 // checks whether or not an event is in a bad LB when the bad LB list is read in manually
782 // @ TODO: Do we still need this functionality? Aren't we always going to read in the list at the merging stage now?
783 bool isBadLB = CheckBadLBList(lumiBlock,BadLBList);
784
785 // skip bad LBs
786 if (isBadLB) continue;
787
788 // check the event satisfies selection criteria
789 bool ProcessEvent = CheckEventSelectionCriteria(lumiBlock, nsigmaHits, energy, noise, lbmin, lbmax);
790 if (ProcessEvent){
791
792 if(m_SaveTextFile || m_SaveRootFile){
793 h1_e->Fill(data->energy()*GeV);
794 h1_q->Fill(data->quality());
795 h1_lb->Fill(lumiBlock);
796 h2_elb->Fill(lumiBlock, data->energy()*GeV);
797 h1_ADCmax->Fill(data->adcMax());
798 double time = data->maxPosition()*25.0+data->ofcTime();
799 h2_t_LB->Fill(lumiBlock,time);
800 for(unsigned int isample=0;isample<data->nSamples();isample++) {
801 TProf_pulse->Fill(isample, data->value(isample));
802 h2_pulse->Fill(data->time(isample)-time, data->value(isample)/data->adcMax());
803 }
804 }
805
806 // count the number of processed events per cell with a bad q-factor
807 if (data->quality() > qthreshold){
808 qcount = qcount+1;
809 }
810
811 // count the number of events processed per cell
812 EventCount++;
813
814 // count number of processed events with energy > ecut and calculate cumulative energy
815 if (data->energy()*GeV > ecut) {
816 nEvents_E_gt_ecut++;
817 EventEnergySum+=(data->energy()*GeV);
818 }
819
820 // if the cell isn't already in the list check whether or not it should be
821 if(!CellInList){
822 selectionFlag=CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
823 if(selectionFlag==2){
824 LBFlaggedIn = lumiBlock;
825 CellInList = true;
826 }
827 else if(selectionFlag==1){
828 if(m_Algo==2){
829 LBFlaggedIn = lumiBlock;
830 CellInList = true;
831 }
832 else if(m_Algo==0){
833 LBFlaggedIn = lumiBlock;
834 CellInList = true;
835 }
836 }
837 else if(selectionFlag==0){
838 if(m_Algo==2){
839 LBFlaggedIn = lumiBlock;
840 CellInList = true;
841 }
842 else if(m_Algo==1){
843 LBFlaggedIn = lumiBlock;
844 CellInList = true;
845 }
846 }
847 }
848
849 }
850
851 } // end of event loop
852
853 CellCounter++;
854 if ( CellCounter%1000 == 0 ) {
855 printf("Scanning Cell No. = %d\n",CellCounter);
856 }
857 selectionFlag=CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
858
859 // check whether or not cell needs to be put in list
860 if (selectionFlag >= 0){
861 bool ListOutput = true;
862 int MeanCellHitsSelection = 0;
863 int EnergyThresholdSelection = 0;
864 int OR_selection = 0;
865 int AND_selection = 0;
866
867 if (selectionFlag == 2){
868 MeanCellHitsSelection = 1;
869 EnergyThresholdSelection = 1;
870 OR_selection = 1;
871 AND_selection = 1;
872 }
873
874
875 else if (selectionFlag == 1){
876 if (m_Algo == 1) ListOutput = false;
877 MeanCellHitsSelection = 1;
878 EnergyThresholdSelection = 0;
879 OR_selection = 1;
880 AND_selection = 0;
881 }
882
883 else if (selectionFlag == 0){
884 if (m_Algo == 0) ListOutput = false;
885 MeanCellHitsSelection = 0;
886 EnergyThresholdSelection = 1;
887 OR_selection = 1;
888 AND_selection = 0;
889 }
890
891 // record number of cells selected by each algorithm
892 iMeanCellHitsSelection += MeanCellHitsSelection;
893 iEnergyThresholdSelection += EnergyThresholdSelection;
894 iOR_selection += OR_selection;
895 iAND_selection += AND_selection;
896
897 // does the cell occur using both algorithms?
898 if (m_SelectRecurringBadCells){ if(AND_selection == 0) ListOutput = false; }
899 if (m_SelectRecurringBadCells){ if(AND_selection == 1) ListOutput = true; }
900
901 if (ListOutput){
902 fprintf(outputFile,"0x%x \t %d \t \t \t %d \t \t %d \t %d\n",onlid,MeanCellHitsSelection,EnergyThresholdSelection,OR_selection,AND_selection);
903
904
905
906 nBadCells++;
907 if(m_SaveTextFile){
908 int nLB = 0;
909 for (int ilb=1;ilb<=nlb;ilb++){
910 if (h1_lb->GetBinContent(ilb) > 0.) nLB++;
911
912 }
913 fr_LB = (double)nLB/(double)(nlb_corr-nBadLBs);
914 nbrLB = (double)nLB;
915 meanECell = h1_e->GetMean();
916 n_ensig = EventCount;
917 if (EventCount!=0){
918 fr_q4k = (double)qcount/(double)EventCount;
919 }
920 n_ecut = nEvents_E_gt_ecut;
921 }
922
923 if (m_SaveRootFile){
924 int nLB = 0;
925 for (int ilb=1;ilb<=nlb;ilb++){
926 if (h1_lb->GetBinContent(ilb) > 0.) nLB++;
927 }
928 fr_LB = (double)nLB/(double)(nlb_corr-nBadLBs);
929 nbrLB = (double)nLB;
930 meanECell = h1_e->GetMean();
931 n_ensig = EventCount;
932 if (EventCount!=0){
933 fr_q4k = (double)qcount/(double)EventCount;
934 }
935 n_ecut = nEvents_E_gt_ecut;
936 Cellmaps[index]->Fill(eta,phi,meanECell);
937 Eeta.Fill(eta,meanECell);
938 CellEnergy.Fill(meanECell);
939 LBfrac.Fill(fr_LB);
940 if (n_ecut>0) {
941 Emean_NbrEvents.Fill(n_ecut,EventEnergySum/n_ecut);
942 Emean_NbrEvents_part[index]->Fill(n_ecut,EventEnergySum/n_ecut);
943 }
944
945 Qfrac.Fill(fr_q4k);
946 Pulsemaps[index]->Add(Pulsemaps[index].get(),h2_pulse);
947 E_LBmaps[index]->Add(E_LBmaps[index].get(),h2_elb);
948 t_LBmaps[index]->Add(t_LBmaps[index].get(),h2_t_LB);
949 CellsFlagged_LB_part[index]->Fill(LBFlaggedIn);
950 CellsFlagged_LB.Fill(LBFlaggedIn);
951 }
952
953 if (m_SaveRootTree){ tree_cells->Fill(); }
954
955
956
957 if (m_SaveTextFile){
958 if (nEvents_E_gt_ecut>0) { EventEnergyMean = EventEnergySum/nEvents_E_gt_ecut; }
959 fprintf(pFile, "0x%x \t %7s \t %d \t %d \t %d \t %d \t %d \t %f \t %4.0f \t %4.3f \t %d \t %d\n",onlid,m_LarIdTranslator->GetPartitonLayerName(index),cellInfo->feedThrough(),cellInfo->slot(),cellInfo->channel(),EventCount,nEvents_E_gt_ecut,EventEnergyMean,nbrLB,fr_q4k,LBFlaggedIn,selectionFlag);
960
961 }
962
963 }
964 }
965
966 qcount = 0.;
967 EventCount = 0.;
968 nEvents_E_gt_ecut = 0.;
969 EventEnergySum = 0.;
970 EventEnergyMean = 0.;
971 } // end of cell loop
972 printf("Cells selected: MeanCellHitsCut \t EnergyThreshodCut \t Either List \t Both Lists\n");
973 printf("\t \t %d \t \t \t \t %d \t \t %d \t %d\n",iMeanCellHitsSelection,iEnergyThresholdSelection,iOR_selection,iAND_selection);
974 if (m_SaveTextFile){
975 fclose (pFile);
976 }
977 fclose (outputFile);
978
979 int nbad = m_nexpected-nskipped;
980 printf("Skipped %d cells. Cells scanned = %d\n",nskipped,nbad);
981 printf("Cells put in Bad Cell List = %d\n",nBadCells);
982 if (noEvdata>0) printf("No Evdata pointer: %d\n",noEvdata);
983
984 if(m_SaveRootFile){
985 fout->cd();
986 for(int i=0;i<npl;i++){
987 if(Cellmaps[i]->GetEntries()>0) Cellmaps[i]->Write();
988 if(E_LBmaps[i]->GetEntries()>0) E_LBmaps[i]->Write();
989 if(t_LBmaps[i]->GetEntries()>0) t_LBmaps[i]->Write();
990 if(Pulsemaps[i]->GetEntries()>0) Pulsemaps[i]->Write();
991 if(CellsFlagged_LB_part[i]->GetEntries()>0) CellsFlagged_LB_part[i]->Write();
992 if( Emean_NbrEvents_part[i]->GetEntries()>0) Emean_NbrEvents_part[i]->Write();
993 }
994 CellsFlagged_LB.Write();
995 Eeta.Write();
996 CellEnergy.Write();
997 Qfrac.Write();
998 Emean_NbrEvents.Write();
999 LBfrac.Write();
1000 }
1001 fout->Close();
1002 delete h1_lb;
1003 delete h1_e;
1004 delete h1_q;
1005 delete h2_elb;
1006 delete h2_pulse;
1007 delete h2_t_LB;
1008 delete TProf_pulse;
1009 delete h1_ADCmax;
1010 return;
1011
1012}
1013
1014//____________________________________________________________________________________________________
1015void LArCellsEmptyMonitoring::GetLimits_EqLB(const char* inputfile, int& lbmin, int& lbmax, int& emin, int& emax, int& qmin, int& qmax, int& runNumber, const std::vector<int, std::allocator<int> >& BadLBList){
1016 int lb=0.;
1017
1018 // Opening file:
1019 std::unique_ptr<LArSamples::Interface> tuple = Interface::open(inputfile);
1020
1021 unsigned int nchannels = tuple->nChannels();
1022
1023 // -------------------------------------------------------------------------
1024 qmin = emin = lbmin = 99999;
1025 qmax = emax = lbmax = -1;
1026 runNumber = 0.;
1027
1028 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1029 // pass empty cells
1030 if(tuple->historySize(ichan)==0){ continue; }
1031
1032 // process non empty cells
1033 const LArSamples::History* hist = tuple->cellHistory(ichan);
1034
1035 unsigned int ndigits = hist->nData();
1036 if(ndigits==0){ continue; }
1037
1038 // loop on the events for each cells
1039 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
1040 const LArSamples::Data* data = hist->data(idigit);
1041 const LArSamples::EventData& Evdata = data->eventData();
1042 int lumiBlock = Evdata.lumiBlock();
1043
1044
1045 if (data->energy() != 0. && data->noise() != 0.){ // record only events with real energy/noise values
1046 bool isBadLB = CheckBadLBList(lumiBlock,BadLBList);
1047 //skip bad LBs
1048 if (isBadLB) continue;
1049
1050
1051 lb = (int)lumiBlock;
1052 if(lb>lbmax) lbmax = lb;
1053 if(lb<lbmin) lbmin = lb;
1054 double energyGeV = data->energy()/1000.;
1055 if(energyGeV > emax) emax = energyGeV;
1056 if(energyGeV < emin) emin = energyGeV;
1057 if(data->quality() > qmax) qmax = data->quality();
1058 if(data->quality() < qmin) qmin = data->quality();
1059
1060 runNumber = data->run();
1061 }
1062 }
1063 }
1064}
1065
1066//____________________________________________________________________________________________________ USING ALL LBS
1067std::vector<int, std::allocator<int> > LArCellsEmptyMonitoring::GetBadLBList(const char* inputfile, int lbmin, int lbmax, double nsigma, int nlb, const std::vector<int, std::allocator<int> >& DQLBList){
1068 std::map< std::pair<unsigned int, unsigned int>, unsigned int > eventCells_tot;
1069 std::map< std::pair<unsigned int, unsigned int>, double > eventEnergy_tot;
1070 std::map< std::pair<unsigned int, unsigned int>, unsigned int > eventLumiBlock;
1071 std::map< std::pair<std::string, unsigned int>, unsigned int > eventCells;
1072 int lb=0.;
1073 double NSIG = nsigma;
1074
1075 // Opening file:
1076 std::unique_ptr<LArSamples::Interface> tuple = Interface::open(inputfile);
1077 unsigned int nchannels = tuple->nChannels();
1078
1079 // -------------------------------------------------------------------------
1080 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1081 // pass empty cells
1082 if(tuple->historySize(ichan)==0){ continue; }
1083
1084 // process non empty cells
1085 const LArSamples::History* hist = tuple->cellHistory(ichan);
1086 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
1087 int index = m_LarIdTranslator->GetPartitionLayerIndex(cellInfo->calo(),cellInfo->layer());
1088 int calo = cellInfo->calo();
1089
1090 if (m_MaskPresampler){
1091 if (CheckForPresamplerCells(index)){ continue; }
1092 }
1093 if (m_MaskCalo){
1094 if (!CheckForPresamplerCells(index)){ continue; }
1095 }
1096
1097 unsigned int ndigits = hist->nData();
1098 if(ndigits==0){ continue; }
1099
1100 // loop on the events for each cell
1101 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
1102 const LArSamples::Data* data = hist->data(idigit);
1103
1104 const LArSamples::EventData& Evdata = data->eventData();
1105 int lumiBlock = Evdata.lumiBlock();
1106 //int lumiBlock = data->lumiBlock();
1107
1108 lb = (int)lumiBlock;
1109 double energyGeV = data->energy()/1000.;
1110
1111
1112 if (lumiBlock <= lbmax && lumiBlock >= lbmin ) { // use only LBs in range
1113 if (data->energy() > sqrt((NSIG*data->noise())*(NSIG*data->noise())) && data->energy() != 0. && data->noise() != 0.){ // record only events with E>NSIG and real energy/noise values
1114 // cell-event mapping
1115 std::pair<unsigned int, unsigned int> ev(data->run(), data->event());
1116 std::pair<std::string, unsigned int> ev_cryo(GetCryostat(calo), data->event());
1117 eventCells_tot[ev]++;
1118 eventEnergy_tot[ev] += energyGeV;
1119 eventLumiBlock[ev] = lb;
1120 eventCells[ev_cryo]++;
1121 }
1122 }
1123 }
1124 }
1125 // -------------------------------------------------------------------------
1126 // Find Bad LBs -> make list (BadLB vector output)
1127 // -------------------------------------------------------------------------
1128 auto hCells_Ev_LB = std::make_unique<TProfile>("","",nlb,lbmin,lbmax,-100,100000);//averaged number of cells per event per LB
1129 std::map< std::string, TProfilep> hCellsEvLB;
1130 std::map< std::string, TH1Fp> hp_cryo;
1131 // -------------------------------------------------------------------------
1132 // loop over events (based on Interface::ShowEvents).
1133 // -------------------------------------------------------------------------
1134 int nCryo = 8;
1135 std::string Cryo[8] = {"EMBA","EMBC","EMECA","EMECC","HECA", "HECC","FCALA","FCALC"};
1136 for(int k = 0; k< nCryo; k++){
1137 hCellsEvLB[Cryo[k]].reset((TProfile*)hCells_Ev_LB->Clone());
1138 }
1139 for(unsigned int ievent = 0; ievent < tuple->nEvents(); ievent++) {
1140 const LArSamples::EventData* evtData = tuple->eventData(ievent);
1141 std::pair<unsigned int, unsigned int> id(evtData->run(), evtData->event());
1142 if(eventCells_tot[id]==0) continue;
1143
1144 for(int j=0; j<nCryo ;j++){
1145 std::pair<std::string, unsigned int> id_cryo(Cryo[j], evtData->event());
1146 hCellsEvLB[Cryo[j]]->Fill(eventLumiBlock[id], eventCells[id_cryo]);
1147 }
1148
1149 hCells_Ev_LB->Fill(eventLumiBlock[id], eventCells_tot[id]);
1150 }
1151
1152 auto hp = std::make_unique<TH1F>("","",(hCells_Ev_LB->GetBinContent(hCells_Ev_LB->GetMaximumBin())),0,hCells_Ev_LB->GetBinContent(hCells_Ev_LB->GetMaximumBin()));// histo to find average number of cells firing per event per LB
1153 for(int k = 0; k< nCryo; k++){
1154 auto hp_temp = std::make_unique< TH1F>("","",(hCellsEvLB[Cryo[k]]->GetBinContent(hCellsEvLB[Cryo[k]]->GetMaximumBin())*10),0,hCellsEvLB[Cryo[k]]->GetBinContent(hCellsEvLB[Cryo[k]]->GetMaximumBin()));// histo to find average number of cells firing per event per LB
1155 std::string name1 = "hp_" + Cryo[k];
1156 hp_cryo[Cryo[k]].reset((TH1F*)hp_temp->Clone(name1.c_str()));
1157 }
1158
1159 for (int p=1;p<=nlb;p++){
1160 for(int j=0; j<nCryo; j++){
1161 double bincont = hCellsEvLB[Cryo[j]]->GetBinContent(p);
1162 hp_cryo[Cryo[j]]->Fill(bincont);
1163 }
1164 double bincont1 = hCells_Ev_LB->GetBinContent(p);
1165 hp->Fill(bincont1);
1166 }
1167
1168 std::map< std::string, double> Mean;
1169 std::map< std::string, double> RMS;
1170 std::map< int, vector<std::string> > BadLB_cryo;
1171
1172 for(int k=0; k<nCryo; k++)
1173 {
1174 Mean[Cryo[k]] = hp_cryo[Cryo[k]]->GetMean();
1175 RMS[Cryo[k]] = hp_cryo[Cryo[k]]->GetRMS();
1176 }
1177
1178 double MEAN2 = hp->GetMean();
1179 double RMS2 = hp->GetRMS();
1180 m_Mean_checkBadLBs = MEAN2;
1181 m_RMS_checkBadLBs = RMS2;
1182 vector<int> BadLB;
1183 // vector<int> BadLB_tot;
1184 vector<int> BadLB_tot=DQLBList;
1185
1186 double value1 = (MEAN2+(RMS2*3.));
1187 int numberFlagged=0;
1188 for (int l=1;l<=nlb;l++){
1189 for(int j=0; j<nCryo; j++){
1190 double value = MEAN2 + (RMS[Cryo[j]]*3);
1191 if(hCellsEvLB[Cryo[j]]->GetBinContent(l) > value){
1192 int badLB = hCellsEvLB[Cryo[j]]->GetBinLowEdge(l);
1193 BadLB_cryo[badLB].push_back(Cryo[j]);
1194 bool isadd = true;
1195 for(unsigned int p=0; p< BadLB_tot.size(); p++) {
1196 if(badLB == BadLB_tot[p]) isadd = false;
1197 }
1198 if(isadd){
1199 BadLB_tot.push_back(badLB);
1200 numberFlagged++;
1201 printf("Removing bad LB %d in %s\n",badLB,Cryo[j].c_str());
1202 }
1203 }
1204 }
1205
1206 if (hCells_Ev_LB->GetBinContent(l) > value1){
1207 BadLB.push_back(hCells_Ev_LB->GetBinLowEdge(l));
1208 }
1209 }
1210 printf("\n");
1211 printf("Number of Bad LBs found: %d\n",numberFlagged);
1212 printf("Number of LBs to be removed: %zu", BadLB_tot.size());
1213 printf("\n");
1214
1215 return BadLB_tot;
1216}
1217
1218
1219//____________________________________________________________________________________________________
1220bool LArCellsEmptyMonitoring::CheckBadLBList(int lumiBlock, const std::vector<int, std::allocator<int> >& BadLBList){
1221 for (unsigned int yy=0;yy<BadLBList.size();yy++){
1222 if (lumiBlock == BadLBList[yy]){
1223 return true;
1224 }
1225 }
1226 return false;
1227}
1228
1229//____________________________________________________________________________________________________
1230void LArCellsEmptyMonitoring::GetMeanCellHits(const char* inputfile, int nlb, int lbmin, int lbmax, int nsigma, const std::vector<int, std::allocator<int> >& BadLBList, double& MeanHits, double& rmsHits, int& nlb_corr)
1231{
1232
1233
1234 int nHits = 0., lumiBlock = 0.,nCells = 0.;
1235 double energy = 0, noise = 0;
1236 std::unique_ptr<LArSamples::Interface> tuple(Interface::open(inputfile));
1237 unsigned int nchannels = tuple->nChannels();
1238
1239 double TotalRecordedHits=0;
1240 std::vector<int, std::allocator<int> > HitsPerLB;
1241 double var=0;
1242 TH1Fp th1_Hits = std::make_unique<TH1F>("","",nlb,lbmin,lbmax); // number of events (E>nsig) in each LB for each cell
1243 TH1Fp hNLB = std::make_unique<TH1F>("","",nlb,lbmin,lbmax); // number of lumiblocks
1244
1245 // -------------------------------------------------------------------------
1246 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1247 // pass empty cells
1248 if(tuple->historySize(ichan)==0){ continue; }
1249
1250
1251 // process non empty cells
1252 const LArSamples::History* hist = tuple->cellHistory(ichan);
1253
1254 unsigned int ndigits = hist->nData();
1255
1256 // loop on the events for each cells
1257 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
1258 const LArSamples::Data* data = hist->data(idigit);
1259 const LArSamples::EventData& Evdata = data->eventData();
1260 lumiBlock = Evdata.lumiBlock();
1261
1262 energy = data->energy();
1263 noise = data->noise();
1264 bool isBadLB = CheckBadLBList(lumiBlock,BadLBList);
1265 if (isBadLB) continue; //skip bad LBs
1266 hNLB->Fill(lumiBlock);
1267
1268 bool ProcessEvent = CheckEventSelectionCriteria(lumiBlock, nsigma, energy, noise, lbmin, lbmax);
1269 if (ProcessEvent){
1270
1271
1272 th1_Hits->Fill(lumiBlock);
1273 nHits++;
1274 }
1275 }
1276 if (nHits > 0.) nCells++;
1277
1278 nHits = 0.;
1279
1280 }
1281
1282
1283 // -------------------------------------------------------------------------
1284 // determine average #hits per cell per LB
1285 // -------------------------------------------------------------------------
1286
1287
1288
1289 TH1Fp tp_ev = std::make_unique<TH1F>("","",th1_Hits->GetBinContent(th1_Hits->GetMaximumBin())*100,0,th1_Hits->GetBinContent(th1_Hits->GetMaximumBin()));
1290
1291 for (int i=1;i<=nlb;i++){
1292 if(hNLB->GetBinContent(i)>0){
1293 nlb_corr++;
1294 tp_ev->Fill(th1_Hits->GetBinContent(i));
1295 TotalRecordedHits+=th1_Hits->GetBinContent(i);
1296 HitsPerLB.push_back(th1_Hits->GetBinContent(i));
1297 }
1298 }
1299
1300 MeanHits = (nCells==0||nlb_corr==0) ? 0 : static_cast<double>(TotalRecordedHits)/(nlb_corr*nCells);
1301 if (nCells > 0){
1302 for (unsigned int j=0;j<HitsPerLB.size();j++){
1303 var += (((double)HitsPerLB[j]/(double)nCells) - MeanHits)*(((double)HitsPerLB[j]/(double)nCells) - MeanHits);
1304 }
1305 } else {
1306 std::cout<<"nCells is " << nCells << " in LArCellsEmptyMonitoring::GetMeanCellHits"<<std::endl;
1307 return;
1308 }
1309 if (nlb_corr > 0){
1310 rmsHits = var/nlb_corr;
1311 }
1312 printf("Mean number of hits/cell for 1 LB = %4.3f, RMS = %4.3f\n",MeanHits,rmsHits);
1313 printf("Number of cells firing at E > %d sigma = %d\n",nsigma,nCells);
1314 printf("Total number of LBs included = %d\n",nlb_corr);
1315
1316
1317}
1318
1319//____________________________________________________________________________________________________
1320bool LArCellsEmptyMonitoring::CheckEventSelectionCriteria(int lumiBlock, int nsigma, double energy, Float_t noise, int lbmin, int lbmax){
1321
1322
1323 if (energy == 0 || noise == 0){
1324 return false;
1325 } else if (lumiBlock <= lbmax && lumiBlock >= lbmin && energy > sqrt((nsigma*noise)*(nsigma*noise)) ) { // use only LBs in range
1326 return true;
1327 }
1328 return false;
1329
1330}
1331
1332
1333//____________________________________________________________________________________________________
1334int LArCellsEmptyMonitoring::CheckCellSelectionCriteria(int EventCount, int nsigmaHits, double MeanHits, double rmsHits, int nEvents_E_gt_ecut, double EventEnergySum, int nBadLBs, int nlb) const{
1335 //FLAG KEY:
1336 //selectionFlag == 0 -> Used EventEnergyCut
1337 //selectionFlag == 1 -> Used MeanCellHitCut
1338 //selectionFlag == 2 -> Used EventEnergyCut && MeanCellHitCut
1339
1340
1341 bool a=false;
1342 bool b=false;
1343 bool c=false;
1344
1345 //pas bon if ((EventCount > (MeanHits[GetSectionIndex(index)]+(nsigmaHits*rmsHits[GetSectionIndex(index)]))*(nlb-nBadLBs[GetCryoIndex(index)]))) a=true;
1346
1347 //original
1348 if ((EventCount > (MeanHits+(nsigmaHits*rmsHits))*((double)nlb-(double)nBadLBs))) a=true;
1349 // reverse cut
1350 // if ((EventCount < (MeanHits+(nsigmaHits*rmsHits))*((double)nlb-(double)nBadLBs))) a=true;
1351
1352 //original
1353 if (nEvents_E_gt_ecut > m_UpperCountThreshold && (EventEnergySum/(double)nEvents_E_gt_ecut) > m_LowerCellEnergyThreshold) b=true;
1354 // reverse cut
1355 // if (nEvents_E_gt_ecut < m_UpperCountThreshold && (EventEnergySum/(double)nEvents_E_gt_ecut) < m_LowerCellEnergyThreshold) b=true;
1356
1357 //original
1358 if (nEvents_E_gt_ecut > m_LowerCountThreshold && (EventEnergySum/(double)nEvents_E_gt_ecut) > m_UpperCellEnergyThreshold) c=true;
1359 // reverse cut
1360 // if (nEvents_E_gt_ecut < m_LowerCountThreshold && (EventEnergySum/(double)nEvents_E_gt_ecut) < m_UpperCellEnergyThreshold) c=true;
1361
1362
1363 if ((a && b) || (a && c)){
1364 return 2;
1365 } else if (a && !b && !c){
1366 return 1;
1367 } else if ((!a && b) || (!a && c)){
1368 return 0;
1369 }
1370
1371 return -1;
1372}
1373
1374//____________________________________________________________________________________________________
1375bool LArCellsEmptyMonitoring::CheckForPresamplerCells(int index){
1376
1377 return index == 0 || index == 4;
1378}
1379
1380//____________________________________________________________________________________________________
1384
1385//____________________________________________________________________________________________________
1386void LArCellsEmptyMonitoring::SetPartition(bool setpart, const std::string& partname)
1387{
1388 const int npl = m_LarIdTranslator->GetNpl();
1389
1390 m_SetPartition = setpart;
1391
1392 if(!setpart)
1393 {
1394 m_PartitionName = "";
1395 m_PartitionIndex = -1;
1396 }
1397 else if(setpart)
1398 {
1399 bool isFoundPartition = false;
1400 for(int i=0; i< npl; i++){
1401 if(partname == m_LarIdTranslator->GetPartitonLayerName(i)){
1402 m_PartitionName = "_" +partname;
1403 m_PartitionIndex = i;
1404 isFoundPartition = true;
1405 }
1406 }
1407 if(!isFoundPartition){
1408 printf("ERROR: Partition %s does not exist! Invalid input name.\n",partname.c_str());
1409 printf("Possible candidates are:\n");
1410 for (int j=0;j<npl;j++){ printf("%s\n",m_LarIdTranslator->GetPartitonLayerName(j)); }
1411 std::abort();
1412 }
1413 }
1414}
1415//____________________________________________________________________________________________________
1416std::string LArCellsEmptyMonitoring::GetCryostat(int calo)
1417{
1418 std::string cryostat = "NotGiven";
1419 if(calo == 1) cryostat = "EMBA";
1420 else if(calo == -1) cryostat = "EMBC";
1421 else if(calo == 5) cryostat = "FCALA";
1422 else if(calo == -5) cryostat = "FCALC";
1423 else if(calo == 4) cryostat = "HECA";
1424 else if(calo == -4) cryostat = "HECC";
1425 else if(calo == 2 || calo == 3) cryostat = "EMECA";
1426 else if(calo == -2 || calo == -3) cryostat = "EMECC";
1427 return cryostat;
1428}
1429
1430
1431//____________________________________________________________________________________________________
1432void LArCellsEmptyMonitoring::ReadDefectLBList(bool ReadList, const TString& LBfile)
1433{
1434
1435
1436 if (!ReadList){
1437 m_LBfile="";
1438 m_ReadDefectLBList=false;
1439 }
1440 else{
1441 m_LBfile=LBfile;
1442 m_ReadDefectLBList=true;
1443 }
1444}
1445
1446//____________________________________________________________________________________________________
1447std::vector<int> LArCellsEmptyMonitoring::ReadBadLBList(const TString& LBfile)
1448{
1449 std::vector<int> LBList;
1450
1451 std::ifstream infile(LBfile.Data());
1452 std::string line;
1453
1454 // assume single-line format with coma-separated LBs (from python)
1455 std::getline(infile,line,'\n');
1456 TString filter(line.c_str());
1457 std::unique_ptr<TObjArray> list(filter.Tokenize(", ")); // coma\space delimiters
1458 if(list->GetEntries() == 0){
1459 printf("No LB filtering specified, or bad format. Exiting.\n");
1460 LBList.push_back(0);
1461 return LBList;
1462 }
1463
1464 for(int k = 0; k < list->GetEntries(); k++){
1465 TObjString* tobs = (TObjString*)(list->At(k));
1466 LBList.push_back((int)(tobs->String()).Atoi());
1467 }
1468 printf("LB List: %d\n",(int)LBList.size());
1469 return LBList;
1470}
1471/*
1472
1473std::vector<int, std::allocator<int> > LArCellsEmptyMonitoring::ReadBadLBList(const TString LBfile)
1474{
1475 // Reads in defect LB list if necessary
1476 std::vector<int, std::allocator<int> > tlist;
1477 std::ifstream infile(LBfile);
1478 std::string argStr;
1479 std::getline(infile,argStr);
1480 for (size_t i = 0,n; i <= argStr.length(); i = n+1)
1481 {
1482 n = argStr.find_first_of(',',i);
1483 if (n == std::string::npos) n = argStr.length();
1484 std::string tmp = argStr.substr(i,n-i);
1485 char ctmp[6];
1486 for (int ic=0;ic<6;ic++){
1487 ctmp[ic]=tmp[ic];
1488 }
1489 int temp = atoi (ctmp);
1490 // printf("temp: %d\n",temp);
1491 tlist.push_back(temp);
1492 }
1493 return tlist;
1494
1495
1496
1497}
1498
1499*/
1500//____________________________________________________________________________________________________
1501void LArCellsEmptyMonitoring::ScanOnlids(const TString& inputfile)
1502{
1503 // Scans all channel entries in the ntuple and looks for duplicates onlids.
1504
1505 int index = -1;
1506 ULong64_t onlid = 0;
1507 std::map<ULong64_t,unsigned int> idmap;
1508 std::map<ULong64_t,unsigned int>::iterator idmap_itr;
1509 int nskipped=0,nrepeated=0;
1510
1511 // Opening file:
1512 std::unique_ptr<LArSamples::Interface> tuple(Interface::open(inputfile));
1513 printf("Number of events: %u %u\n",tuple->nEvents(),tuple->nChannels()); //tuple->ShowEvents("energy>0.");
1514 unsigned int nchannels = tuple->nChannels();
1515
1516 // -------------------------------------------------------------------------
1517 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1518 // pass empty cells
1519 if(tuple->historySize(ichan)==0){ nskipped++; continue; }
1520
1521 // process non empty cells
1522 const LArSamples::History* hist = tuple->cellHistory(ichan);
1523 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
1524
1525 // index for partition
1526 index = m_LarIdTranslator->GetPartitionLayerIndex(cellInfo->calo(),cellInfo->layer());
1527 onlid = cellInfo->onlid();
1528 /*
1529 if(m_LarIdTranslator->FindChannel(index,cellInfo->eta(),cellInfo->phi())){
1530 onlid = m_LarIdTranslator->onlid;
1531 } else onlid = 0;
1532 */
1533 if(onlid<=0) printf("%u: Bad Cell Onlid = 0x%x (%+.2f,%+.2f)\n",ichan,(unsigned int)onlid,cellInfo->eta(),cellInfo->phi());
1534
1535 idmap_itr = idmap.find(onlid);
1536
1537 // new onlid
1538 if(idmap_itr == idmap.end()){
1539 idmap[onlid] = ichan;
1540 } else {
1541 nrepeated+=1;
1542 printf("Onlid 0x%x (%d,%s,%+.2f,%+.2f)\n",(unsigned int)onlid,index,m_LarIdTranslator->GetPartitonLayerName(index),
1543 m_LarIdTranslator->eta,m_LarIdTranslator->phi);
1544 }
1545 }
1546
1547 printf("Skipped %d cells.\n",nskipped);
1548 printf("Number of onlids: Unique=%d, Repeated=%d\n",(int)idmap.size(),nrepeated);
1549
1550 return;
1551}
1552
1553//____________________________________________________________________________________________________
1554void LArCellsEmptyMonitoring::DoEtaPhiMonitoring(const char* inputfile,const char* optionplot,const char* optionsave)
1555{
1556 // Main monitoring function with primary goal to produce simple DQM-equivalent (eta,phi) maps.
1557 // Depending on the argument, the function will put different quantities in the histograms.
1558 // Currently supported options:
1559 //
1560 // - optionplot = "NAbove4Sigma","PercentageAbove4Sigma",
1561 // "NAbove5Sigma","PercentageAbove5Sigma",
1562 // "NQualityAbove4000","PercentageQualityAbove4000"
1563 // - optionsave = "root","png"
1564 //
1565 //
1566
1567 printf("LArCellsEmptyMonitoring::DoEtaPhiMonitoring: Reading options: %s %s.\n",optionplot,optionsave);
1568
1569 bool knormplot = false;
1570 if(strstr(optionplot,"Precentage")) knormplot = true; // needs normalization plots
1571
1572 int kcuttype = 0; // no cut
1573 double nsigmacut = 0.;
1574 double qualitycut = 0.;
1575 if(strstr(optionplot,"Sigma")){
1576 kcuttype = 1; // energy
1577 if(strstr(optionplot,"4Sigma")) nsigmacut = 4.;
1578 if(strstr(optionplot,"5Sigma")) nsigmacut = 5.;
1579 } else if(strstr(optionplot,"Quality")){
1580 kcuttype = 2; // quality
1581 qualitycut = 4000.; // not much control for now
1582 }
1583
1584 Char_t hname[512];
1585 int index = 0;
1586 int nskipped=0;
1587 ULong64_t onlid = 0;
1588
1589 // Opening file:
1590 std::unique_ptr<LArSamples::Interface> tuple(Interface::open(inputfile));
1591 printf("Number of events: %u %u\n",tuple->nEvents(),tuple->nChannels());
1592 unsigned int nchannels = tuple->nChannels();
1593
1594 // eta-phi maps
1595 const int nhists = m_LarIdTranslator->GetNpl();
1596 std::vector<TH2*> hmap_counts_all (nhists);
1597 std::vector<TH2*> hmap_energy_cut (nhists);
1598 std::vector<TH2*> hmap_quality_cut (nhists);
1599 for(int j=0;j<nhists;j++){
1600 hmap_counts_all[j] = m_LarIdTranslator->GetCaloPartitionLayerMap(j);
1601 sprintf(hname,"counst_all_%s_%d",m_LarIdTranslator->GetPartitonLayerName(j),j);
1602 hmap_counts_all[j]->SetName(hname);
1603 //
1604 hmap_energy_cut[j] = m_LarIdTranslator->GetCaloPartitionLayerMap(j);
1605 sprintf(hname,"energy_cut_%s_%d",m_LarIdTranslator->GetPartitonLayerName(j),j);
1606 hmap_energy_cut[j]->SetName(hname);
1607 //
1608 hmap_quality_cut[j] = m_LarIdTranslator->GetCaloPartitionLayerMap(j);
1609 sprintf(hname,"quality_cut_%s_%d",m_LarIdTranslator->GetPartitonLayerName(j),j);
1610 hmap_quality_cut[j]->SetName(hname);
1611 }
1612
1613 // -------------------------------------------------------------------------
1614 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1615 // pass empty cells
1616 if(tuple->historySize(ichan)==0){ nskipped++; continue; }
1617
1618 // process non empty cells
1619 const LArSamples::History* hist = tuple->cellHistory(ichan);
1620 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
1621 unsigned int ndigits = hist->nData();
1622 if(ndigits==0){ nskipped++; continue; }
1623
1624 // index for partition; this can return a negative value in case of error
1625 index = m_LarIdTranslator->GetPartitionLayerIndex(cellInfo->calo(),cellInfo->layer());
1626 if (index<0){
1627 std::cout<<"GetPartitionLayerIndex returned -1 in LArCellsEmptyMonitoring::DoEtaPhiMonitoring"<<std::endl;
1628 continue;
1629 }
1630 onlid = cellInfo->onlid();
1631 if(onlid<=0) printf("%u: Bad Cell Onlid = 0x%x (%+.2f,%+.2f)\n",ichan,(unsigned int)onlid,cellInfo->eta(),cellInfo->phi());
1632
1633 // loop on the events for each cells
1634 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
1635 const LArSamples::Data* data = hist->data(idigit);
1636
1637 // all stats
1638 hmap_counts_all[index]->Fill(m_LarIdTranslator->eta,m_LarIdTranslator->phi);
1639
1640 // energy-like
1641 if(kcuttype==1 && data->noise()>0.){
1642 if((data->energy()/data->noise())>nsigmacut) hmap_energy_cut[index]->Fill(m_LarIdTranslator->eta,m_LarIdTranslator->phi);
1643 }
1644 if(kcuttype==2){
1645 if(data->quality()>qualitycut) hmap_quality_cut[index]->Fill(m_LarIdTranslator->eta,m_LarIdTranslator->phi);
1646 }
1647 }
1648 }
1649
1650 int nbad = m_nexpected-nskipped;
1651 printf("Skipped %d cells. Bad Cells = %d\n",nskipped,nbad);
1652
1653 // normalize, if needed
1654 if(knormplot){
1655 for(int j=0;j<nhists;j++){
1656 hmap_energy_cut[j]->Divide(hmap_counts_all[j]);
1657 hmap_quality_cut[j]->Divide(hmap_counts_all[j]);
1658 }
1659 }
1660
1661 std::unique_ptr<TFile> tout;
1662 TCanvas* c0 = nullptr, *c1 = nullptr;
1663
1664 if(!strcmp(optionsave,"root")){
1665 tout.reset(new TFile("EtaPhiMonitoring.root","recreate"));
1666 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_counts_all.data(),"Counts",1);
1667 c0->SetName("Normalization");
1668 c0->Write();
1669 for(int j=0;j<nhists;j++) hmap_counts_all[j]->Write();
1670 if(kcuttype==1){
1671 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_energy_cut.data(),"EnergyCut",1);
1672 c0->SetName("EnergyCut");
1673 c0->Write();
1674 for(int j=0;j<nhists;j++) hmap_energy_cut[j]->Write();
1675 }
1676 if(kcuttype==2){
1677 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_quality_cut.data(),"QualityCut",1);
1678 c0->SetName("QualityCut");
1679 c0->Write();
1680 for(int j=0;j<nhists;j++) hmap_quality_cut[j]->Write();
1681 }
1682 tout->Close();
1683 } else {
1684 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_counts_all.data(),"Counts",1);
1685 c0->SaveAs("Normalization.png");
1686 c1 = new TCanvas("c1","");
1687 for(int j=0;j<nhists;j++){ hmap_counts_all[j]->Draw("colz"); sprintf(hname,"%s.png",hmap_counts_all[j]->GetName()); c1->SaveAs(hname); }
1688 if(kcuttype==1){
1689 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_energy_cut.data(),"EnergyCut",1);
1690 c0->SaveAs("EnergyCut.png");
1691 c1->cd();
1692 for(int j=0;j<nhists;j++){ hmap_energy_cut[j]->Draw("colz"); sprintf(hname,"%s.png",hmap_energy_cut[j]->GetName()); c1->SaveAs(hname); }
1693 }
1694 if(kcuttype==2){
1695 c0 = m_LarIdTranslator->CaloPartitionLayerDisplay(hmap_quality_cut.data(),"QualityCut",1);
1696 c0->SaveAs("QualityCut.png");
1697 c1->cd();
1698 for(int j=0;j<nhists;j++){ hmap_quality_cut[j]->Draw("colz"); sprintf(hname,"%s.png",hmap_quality_cut[j]->GetName()); c1->SaveAs(hname); }
1699 }
1700 }
1701 return;
1702}
1703
1704//____________________________________________________________________________________________________
1705void LArCellsEmptyMonitoring::TriggerEfficiency(const char* inputfile,float fractionInPS)
1706{
1707 // Function aimed at extracting the trigger performance of a specific chain
1708 // by emulating its online selection using the calibration stream data.
1709 //
1710 // This is currently used to monitor the rejection performance of the
1711 // PreSampler noise in the physics_CosmicCalo stream (L2_PreS* chains)
1712 // which runs on each family of triggers (EM,TAU,J).
1713 // One needs to set the fraction of energy deposited in the PreSampler with
1714 // fractionInPS.
1715 //
1716
1717
1718 std::map< std::pair<unsigned int, unsigned int>, unsigned int > eventCells;
1719 std::map< std::pair<unsigned int, unsigned int>, double > eventEnergy_PS;
1720 std::map< std::pair<unsigned int, unsigned int>, double > eventEnergy_LAr;
1721 std::map< std::pair<unsigned int, unsigned int>, unsigned int > eventCells_PS;
1722 std::map< std::pair<unsigned int, unsigned int>, unsigned int > eventLayer;
1723 int run=0;
1724 // Opening file:
1725 std::unique_ptr<LArSamples::Interface> tuple(Interface::open(inputfile));
1726 unsigned int nchannels = tuple->nChannels();
1727
1728 // -------------------------------------------------------------------------
1729 for(unsigned int ichan = 0; ichan < nchannels; ichan++){
1730 // pass empty cells
1731 if(tuple->historySize(ichan)==0){ continue; }
1732
1733 // process non empty cells
1734 const LArSamples::History* hist = tuple->cellHistory(ichan);
1735 const LArSamples::CellInfo* cellInfo = hist->cellInfo();
1736 unsigned int layer = cellInfo->layer();
1737
1738 unsigned int ndigits = hist->nData();
1739 if(ndigits==0){ continue; }
1740
1741 // loop on the events for each cells
1742 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
1743 const LArSamples::Data* data = hist->data(idigit);
1744 run = data->run();
1745
1746 // convert to GeV
1747 double energyGeV = data->energy()/1000.;
1748
1749 // cell-event mapping
1750 std::pair<unsigned int, unsigned int> ev(data->run(), data->event());
1751 eventCells[ev]++;
1752 eventEnergy_LAr[ev] += energyGeV;
1753 eventLayer[ev] = layer;
1754 if(layer == 0){
1755 eventEnergy_PS[ev] += energyGeV;
1756 eventCells_PS[ev]++;
1757 }
1758 }
1759 }
1760
1761 // non-constant binning for a wider spetrum
1762 //double binE[48]={0,2,4,6,8,10,13,16,19,22,26,30,34,38,42,47,52,57,62,68,75,82,90,100,110,120,130,140,150,160,170,180,195,210,225,240,255,270,285,300,320,340,360,380,400,430,460,500};
1763
1764 constexpr int nbins = 42;
1765 double binE[nbins+1]{};
1766 int i=0,j=0;
1767 int stop=0;
1768 double de = 0.,ene = 0.;
1769 // steps of 2GeV for the first 10 bins [0-20]
1770 de = 2.;
1771 stop = (int)(20./de);
1772 for(i=0;i<stop;i++,j++){
1773 ene = de*(double)i;
1774 binE[j] = ene;
1775 }
1776 // steps of 5GeV for [20-100]
1777 de = 5.;
1778 stop = (int)((100.-20.)/de);
1779 for(i=0;i<stop;i++,j++){
1780 ene += de;
1781 binE[j] = ene;
1782 }
1783 // steps of 25GeV for [100-500]
1784 de = 25.;
1785 stop = (int)((500.-100.)/de);
1786 for(i=0;i<stop;i++,j++){
1787 ene += de;
1788 binE[j] = ene;
1789 }
1790 // add extra bin to cover the rest
1791 binE[j] = 3500.;
1792
1793
1794 TH1F E_LAr_pass("E_LAr_pass","",nbins,binE);
1795 TH1F E_LAr_tot("E_LAr_tot","",nbins,binE);
1796 TH1F E_LAr_EM5_pass("E_LAr_EM5_pass","",nbins,binE);
1797 TH1F E_LAr_EM5_tot("E_LAr_EM5_tot","",nbins,binE);
1798 TH1F E_LAr_TAU8_pass("E_LAr_TAU8_pass","",nbins,binE);
1799 TH1F E_LAr_TAU8_tot("E_LAr_TAU8_tot","",nbins,binE);
1800 TH1F E_LAr_J10_pass("E_LAr_J10_pass","",nbins,binE);
1801 TH1F E_LAr_J10_tot("E_LAr_J10_tot","",nbins,binE);
1802 TH1F nCellsPS_pass("nCellsPS_pass","",100, 0, 100);
1803 TH1F nCellsPS_tot("nCellsPS_tot","",100, 0, 100);
1804 TH1F nCellsPS_EM5_pass("nCellsPS_EM5_pass","",100, 0, 100);
1805 TH1F nCellsPS_EM5_tot("nCellsPS_EM5_tot","",100, 0, 100);
1806 TH1F nCellsPS_TAU8_pass("nCellsPS_TAU8_pass","",100, 0, 100);
1807 TH1F nCellsPS_TAU8_tot("nCellsPS_TAU8_tot","",100, 0, 100);
1808 TH1F nCellsPS_J10_tot("nCellsPS_J10_tot","",100, 0, 100);
1809 TH1F nCellsPS_J10_pass("nCellsPS_J10_pass","",100, 0, 100);
1810
1811 // -------------------------------------------------------------------------
1812 // loop over events (based on Interface::ShowEvents).
1813 // -------------------------------------------------------------------------
1814 for(unsigned int ievent = 0; ievent < tuple->nEvents(); ievent++) {
1815 const LArSamples::EventData* evtData = tuple->eventData(ievent);
1816 std::pair<unsigned int, unsigned int> id(evtData->run(), evtData->event());
1817 if(eventCells[id]==0) continue;
1818 if(eventEnergy_LAr[id]<=0) continue;
1819 bool isEM5 = evtData->isPassed("L1_EM5_EMPTY");
1820 bool isTAU8 = evtData->isPassed("L1_TAU8_EMPTY");
1821 bool isJ10 = evtData->isPassed("L1_J10_EMPTY");
1822
1823 E_LAr_tot.Fill(eventEnergy_LAr[id]);
1824 nCellsPS_tot.Fill(eventCells_PS[id]);
1825 if(isEM5){
1826 E_LAr_EM5_tot.Fill(eventEnergy_LAr[id]);
1827 nCellsPS_EM5_tot.Fill(eventCells_PS[id]);
1828 }
1829 if(isTAU8){
1830 E_LAr_TAU8_tot.Fill(eventEnergy_LAr[id]);
1831 nCellsPS_TAU8_tot.Fill(eventCells_PS[id]);
1832 }
1833 if(isJ10){
1834 E_LAr_J10_tot.Fill(eventEnergy_LAr[id]);
1835 nCellsPS_J10_tot.Fill(eventCells_PS[id]);
1836 }
1837
1838 double ratio = 0;
1839 if((eventEnergy_LAr[id])!=0)
1840 ratio = eventEnergy_PS[id]/(eventEnergy_LAr[id]);
1841 bool isL2_PreS = false;
1842 if(ratio > fractionInPS) isL2_PreS = true;
1843
1844 if(isL2_PreS){
1845 E_LAr_pass.Fill(eventEnergy_LAr[id]);
1846 nCellsPS_pass.Fill(eventCells_PS[id]);
1847 if(isEM5){
1848 E_LAr_EM5_pass.Fill(eventEnergy_LAr[id]);
1849 nCellsPS_EM5_pass.Fill(eventCells_PS[id]);
1850 }
1851 if(isTAU8){
1852 E_LAr_TAU8_pass.Fill(eventEnergy_LAr[id]);
1853 nCellsPS_TAU8_pass.Fill(eventCells_PS[id]);
1854 }
1855 if(isJ10){
1856 E_LAr_J10_pass.Fill(eventEnergy_LAr[id]);
1857 nCellsPS_J10_pass.Fill(eventCells_PS[id]);
1858 }
1859 }
1860 }
1861
1862 char fname[50]{};
1863 sprintf(fname,"TriggerEfficiency_%d.root",run);
1864 auto mfile = std::make_unique<TFile>(fname,"recreate");
1865 E_LAr_pass.Write();
1866 E_LAr_tot.Write();
1867 E_LAr_EM5_pass.Write();
1868 E_LAr_EM5_tot.Write();
1869 E_LAr_TAU8_pass.Write();
1870 E_LAr_TAU8_tot.Write();
1871 E_LAr_J10_pass.Write();
1872 E_LAr_J10_tot.Write();
1873 nCellsPS_pass.Write();
1874 nCellsPS_tot.Write();
1875 nCellsPS_EM5_pass.Write();
1876 nCellsPS_EM5_tot.Write();
1877 nCellsPS_TAU8_pass.Write();
1878 nCellsPS_TAU8_tot.Write();
1879 nCellsPS_J10_pass.Write();
1880 nCellsPS_J10_tot.Write();
1881
1882 mfile->Close();
1883
1884}
1885
1886
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
data
char data[hepevt_bytes_allocation_ATLAS]
Definition HepEvt.cxx:11
TProfilep
std::unique_ptr< TProfile > TProfilep
Definition LArCellsEmptyMonitoring.cxx:61
THVec
std::vector< std::unique_ptr< T > > THVec
Definition LArCellsEmptyMonitoring.cxx:54
TH2Fp
std::unique_ptr< TH2F > TH2Fp
Definition LArCellsEmptyMonitoring.cxx:60
TH1Fp
std::unique_ptr< TH1F > TH1Fp
Definition LArCellsEmptyMonitoring.cxx:59
THArray
std::array< std::unique_ptr< T >, n > THArray
Definition LArCellsEmptyMonitoring.cxx:57
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
nHits
static const uint32_t nHits
Definition PixelByteStreamModuleMask.h:55
GetEntries
TGraphErrors * GetEntries(TH2F *histo)
Definition TRTCalib_makeplots.cxx:4025
LArSamples::CellInfo::eta
double eta() const
Definition CellInfo.h:93
LArSamples::CellInfo::slot
short slot() const
Definition CellInfo.h:68
LArSamples::CellInfo::layer
short layer() const
Definition CellInfo.h:53
LArSamples::CellInfo::feedThrough
short feedThrough() const
Definition CellInfo.h:65
LArSamples::CellInfo::phi
double phi() const
Definition CellInfo.h:94
LArSamples::CellInfo::onlid
ULong64_t onlid() const
Definition CellInfo.h:90
LArSamples::CellInfo::calo
CaloId calo() const
Definition CellInfo.h:50
LArSamples::CellInfo::channel
short channel() const
Definition CellInfo.h:76
LArSamples::EventData::isPassed
bool isPassed(const TString &bitName) const
Definition EventData.cxx:88
LArSamples::Interface::open
static std::unique_ptr< Interface > open(const TString &fileName)
Definition Interface.cxx:35
LArSamples::LArCellsEmptyMonitoring::SetPartition
void SetPartition(bool setpart, const std::string &partname)
Definition LArCellsEmptyMonitoring.cxx:1386
LArSamples::LArCellsEmptyMonitoring::GetLimits_EqLB
void GetLimits_EqLB(const char *inputfile, int &lbmin, int &lbmax, int &emin, int &emax, int &qmin, int &qmax, int &runNumber, const std::vector< int, std::allocator< int > > &BadLBList)
Definition LArCellsEmptyMonitoring.cxx:1015
LArSamples::LArCellsEmptyMonitoring::TriggerEfficiency
void TriggerEfficiency(const char *inputfile, float fractionInPS=0.8)
Definition LArCellsEmptyMonitoring.cxx:1705
LArSamples::LArCellsEmptyMonitoring::DoEtaPhiMonitoring
void DoEtaPhiMonitoring(const char *inputfile, const char *optionplot, const char *optionsave)
Definition LArCellsEmptyMonitoring.cxx:1554
LArSamples::LArCellsEmptyMonitoring::m_LarIdTranslator
std::unique_ptr< LArIdTranslatorHelper > m_LarIdTranslator
Definition LArCellsEmptyMonitoring.h:84
LArSamples::LArCellsEmptyMonitoring::GetBadLBList
std::vector< int, std::allocator< int > > GetBadLBList(const char *inputfile, int lbmin, int lbmax, double nsigma, int nlb, const std::vector< int, std::allocator< int > > &DQLBList)
Definition LArCellsEmptyMonitoring.cxx:1067
LArSamples::LArCellsEmptyMonitoring::CheckBadLBList
bool CheckBadLBList(int lumiBlock, const std::vector< int, std::allocator< int > > &BadLBList)
Definition LArCellsEmptyMonitoring.cxx:1220
LArSamples::LArCellsEmptyMonitoring::CheckEventSelectionCriteria
bool CheckEventSelectionCriteria(int lumiBlock, int nsigma, double energy, Float_t noise, int lbmin, int lbmax)
Definition LArCellsEmptyMonitoring.cxx:1320
LArSamples::LArCellsEmptyMonitoring::CheckCellSelectionCriteria
int CheckCellSelectionCriteria(int EventCount, int nsigmaHits, double MeanHits, double rmsHits, int nEvents_E_gt_ecut, double EventEnergySum, int nBadLBs, int nlb) const
Definition LArCellsEmptyMonitoring.cxx:1334
LArSamples::LArCellsEmptyMonitoring::GetMeanCellHits
void GetMeanCellHits(const char *inputfile, int nlb, int lbmin, int lbmax, int nsigma, const std::vector< int, std::allocator< int > > &BadLBList, double &MeanHits, double &rmsHits, int &nlb_corr)
Definition LArCellsEmptyMonitoring.cxx:1230
LArSamples::LArCellsEmptyMonitoring::ReadDefectLBList
void ReadDefectLBList(bool ReadList, const TString &LBfile)
Definition LArCellsEmptyMonitoring.cxx:1432
LArSamples::LArCellsEmptyMonitoring::ReadBadLBList
std::vector< int, std::allocator< int > > ReadBadLBList(const TString &LBfile)
Definition LArCellsEmptyMonitoring.cxx:1447
lb
int lb
Definition globals.cxx:23
ev
int ev
Definition globals.cxx:25
get
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
Definition hcg.cxx:130
xmax
double xmax
Definition listroot.cxx:61
fout
static TFile * fout
Definition listroot.cxx:40
xmin
double xmin
Definition listroot.cxx:60
index
Definition index.py:1
std
STL namespace.
run
int run(int argc, char *argv[])
Definition ttree2hdf5.cxx:28
Generated on for ATLAS Offline Software by doxygen 1.14.0