28 #include "TObjString.h"
35 #include "TProfile2D.h"
37 #include "TDirectory.h"
53 LArCellsEmptyMonitoring::LArCellsEmptyMonitoring()
88 TString htitle,
hname,textfilename,rootfilename;
90 double Ecum_cut_PS = 20.;
91 double Ecum_cut_calo = 10.;
94 const double GeV = 1/1000.;
96 double nsigmaHits = 10.;
102 unsigned int onlid = 0;
113 double MeanEventEnergy_max = -999;
114 double MeanEventEnergy_min = 999;
121 vector<int, std::allocator<int> > DQLBList;
123 vector<int, std::allocator<int> > BadLBList;
126 printf(
"Getting histogram limits... ");
135 printf(
"Done.\n LBmin: %d LBmax: %d Emin: %d Emax: %d Qmin: %d Qmax: %d\n",
lbmin,
lbmax,emin,emax,qmin,qmax);
139 FILE * pFile=
nullptr;
142 textfilename.Form(
"Output/BadCellList_run%d.txt",
runNumber);
143 pFile = fopen (textfilename ,
"w");
145 printf(
"Cannot open output text file\n");
151 pFile = fopen (textfilename ,
"w");
153 printf(
"Cannot open output text file\n");
160 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);
163 int EventCount=0., qcount=0.;
169 rootfilename.Form(
"Output/BadCells_run%d.root",
runNumber);
172 TFile*
fout =
new TFile(rootfilename,
"RECREATE");
175 TH1F* h1_lb =
nullptr;
176 TH1F* h1_elb =
nullptr;
177 TH1F* h1_qlb =
nullptr;
178 TH1F* h1_e =
nullptr;
179 TH1F* h1_q =
nullptr;
180 TH2F* h2_elb =
nullptr;
186 printf(
"Creating histograms (nEbins = %d, nQbins = %d, nLBbins = %d)... ",ne,nq,nlb);
187 TH1F hene =
TH1F(
"hEne",
"",ne,emin,emax); hene.SetXTitle(
"Energy [GeV]"); hene.SetYTitle(
"Events");
188 TH1F hqua =
TH1F(
"hQua",
"",nq/100.,qmin,qmax); hqua.SetXTitle(
"Quality"); hqua.SetYTitle(
"Events");
189 TH1F hlb =
TH1F(
"hLB",
"",nlb,
lbmin,
lbmax); hlb.SetXTitle(
"LB"); hlb.SetYTitle(
"Events");
190 TH2F henelb =
TH2F(
"hEnelb",
"",nlb,
lbmin,
lbmax,ne,emin,emax); henelb.SetXTitle(
"LB"); henelb.SetYTitle(
"Energy [GeV]");
197 printf(
"Creating cut test histograms...");
198 static constexpr
int npl=4;
199 TH2F tempHist =
TH2F(
"tempNEvLB",
"",1000,0,1000.,1000,0.,1000.);
200 std::array<TH1F*, npl> hNoise{};
201 std::array<TH1F*, npl> hCellsPerLB{};
202 std::array<TH2F*, npl> hNEvVsEMean{};
203 std::array<TH2F*, npl> hNEvVsECum{};
204 std::array<TH2F*, npl> hECumVsEMean{};
206 for (
int i=0;
i<npl;
i++){
207 hname.Form(
"hNEvVsEMean_Layer%i",
i);
208 hNEvVsEMean[
i] = (
TH2F*)tempHist.Clone(
hname);
209 hNEvVsEMean[
i]->GetXaxis()->SetTitle(
"Number of Events (E > 10 #sigma) / LB"); hNEvVsEMean[
i]->GetYaxis()->SetTitle(
"E_{mean} per Event / LB [GeV]");
211 hname.Form(
"hNEvVsECum_Layer%i",
i);
213 hNEvVsECum[
i]->GetXaxis()->SetTitle(
"Number of Events (E > 10 #sigma) / LB"); hNEvVsECum[
i]->GetYaxis()->SetTitle(
"E_{cumulative} per Cell / LB [GeV]");
215 hname.Form(
"hECumVsEMean_Layer%i",
i);
216 hECumVsEMean[
i] = (
TH2F*)tempHist.Clone(
hname);
217 hECumVsEMean[
i]->GetXaxis()->SetTitle(
"E_{cumulative} per Cell / LB [GeV]"); hECumVsEMean[
i]->GetYaxis()->SetTitle(
"E_{mean} per Event / LB [GeV]");
219 hname.Form(
"hNoise_Layer%i",
i);
221 hNoise[
i]->GetXaxis()->SetTitle(
"Noise [GeV]"); hNoise[
i]->GetYaxis()->SetTitle(
"Number Of Cells");
223 hname.Form(
"hCellsPerLB_Layer%i",
i);
225 hCellsPerLB[
i]->GetXaxis()->SetTitle(
"LB"); hNoise[
i]->GetYaxis()->SetTitle(
"Number Of Cells");
234 unsigned int nchannels = tuple->
nChannels();
236 printf(
"Begin looping over cells...\n");
248 onlid = cellInfo->
onlid();
251 unsigned int ndigits =
hist->nData();
252 if(ndigits<(
unsigned int)
threshold){ nskipped++;
continue; }
256 hname.Form(
"h0x%x_Energy",onlid);
258 hname.Form(
"h0x%x_Quality",onlid);
260 hname.Form(
"h0x%x_LB",onlid);
262 hname.Form(
"h0x%x_hNEvVsEMean",onlid);
264 hname.Form(
"h0x%x_Energy_LB",onlid);
266 hname.Form(
"h0x%x_Quality_LB",onlid);
271 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
275 if(!Evdata){ noEvdata++;
continue; }
284 h1_q->Fill(
data->quality());
292 if (
data->quality() > qthreshold){
306 if ( CellCounter%1000 == 0 ) {
307 printf(
"Scanning Cell No. = %d\n",CellCounter);
309 MeanEventEnergy_max = -999;
310 MeanEventEnergy_min = 999;
317 if (h1_lb->GetBinContent(
ilb) > 0.){
318 double EtotInLB = h1_elb->GetBinContent(
ilb);
319 double NEventsInLB = h1_lb->GetBinContent(
ilb);
321 int LB_number = h1_lb->GetBinLowEdge(
ilb);
323 double MeanEventEnergy = EtotInLB / NEventsInLB;
326 if (cellInfo->
layer()==0){
327 if (EtotInLB>Ecum_cut_PS){
328 hNEvVsEMean[cellInfo->
layer()]->Fill(NEventsInLB,MeanEventEnergy);
329 hNEvVsECum[cellInfo->
layer()]->Fill(NEventsInLB,EtotInLB);
330 hECumVsEMean[cellInfo->
layer()]->Fill(EtotInLB,MeanEventEnergy);
331 hCellsPerLB[cellInfo->
layer()]->Fill(LB_number);
334 else if (EtotInLB>Ecum_cut_calo){
335 hNEvVsEMean[cellInfo->
layer()]->Fill(NEventsInLB,MeanEventEnergy);
336 hNEvVsECum[cellInfo->
layer()]->Fill(NEventsInLB,EtotInLB);
337 hECumVsEMean[cellInfo->
layer()]->Fill(EtotInLB,MeanEventEnergy);
338 hCellsPerLB[cellInfo->
layer()]->Fill(LB_number);
341 if (EtotInLB>Ecum_cut_PS && cellInfo->
layer()==0){
342 FlagCell=
true; FlagCount++;
343 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
344 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
346 else if (EtotInLB>Ecum_cut_calo && cellInfo->
layer()==1){
347 FlagCell=
true; FlagCount++;
348 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
349 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
351 else if (EtotInLB>Ecum_cut_calo && cellInfo->
layer()==2){
352 FlagCell=
true; FlagCount++;
353 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
354 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
356 else if (EtotInLB>Ecum_cut_calo && cellInfo->
layer()==3){
357 FlagCell=
true; FlagCount++;
358 if (MeanEventEnergy_max<MeanEventEnergy){ MeanEventEnergy_max = MeanEventEnergy; }
359 if (MeanEventEnergy_min>MeanEventEnergy){ MeanEventEnergy_min = MeanEventEnergy; }
364 fr_q4k = (
double)qcount/(
double)EventCount;
371 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 printf(
"Skipped %d cells. Cells scanned = %d\n",nskipped,nbad);
388 printf(
"Cells put in Bad Cell List = %d\n",nBadCells);
389 if (noEvdata>0) printf(
"No Evdata pointer: %d\n",noEvdata);
394 for(
int i=0;
i<npl;
i++){
395 hNEvVsEMean[
i]->Write();
396 hNEvVsECum[
i]->Write();
397 hECumVsEMean[
i]->Write();
399 hCellsPerLB[
i]->Write();
416 TString htitle,
hname,textfilename,rootfilename;
428 printf(
"WARNING: Received instruction to select recurring bad cells only in SetSelectRecurringBadCells(bool flag)\n");
429 printf(
"-------> Only one algorithm selected for cell selection criteria in SetAlgo(int algoindex)\n");
430 printf(
"-------> Recurring cell selection overides algorithm specification\n");
435 printf(
"ERROR: Partition not specified in SetPartition(bool set_part, std::string partname)\n");
442 int selectionFlag = -1;
443 const double GeV = 1/1000.;
453 unsigned int onlid = 0;
465 int nEvents_E_gt_ecut=0;
466 double EventEnergySum=0.;
467 double EventEnergyMean=0.;
470 double MeanHits=0, rmsHits=0;
471 int iMeanCellHitsSelection=0,iEnergyThresholdSelection=0,iOR_selection=0,iAND_selection=0;
477 vector<int, std::allocator<int> > DQLBList;
478 vector<int, std::allocator<int> > BadLBList;
480 printf(
"Reading in DQ Defect LB List...\n");
485 printf(
"Getting bad LB list...\n");
487 printf(
"...Done.\n");
488 nBadLBs = BadLBList.size();
491 BadLBList = std::move(DQLBList);
492 nBadLBs = BadLBList.size();
496 printf(
"Getting histogram limits... ");
507 printf(
"Done.\n LBmin: %d LBmax: %d Emin: %d Emax: %d Qmin: %d Qmax: %d\n",
lbmin,
lbmax,emin,emax,qmin,qmax);
512 textfilename.Form(
"Output/BadCellSelection_run%d.txt",
runNumber);
515 std::cout<<(
"Cannot open output selection summary file " + textfilename)<<
"\n";
521 std::cout<<(
"Cannot open output selection summary file " + textfilename)<<
"\n";
527 FILE * pFile=
nullptr;
530 textfilename.Form(
"Output/BadCellList_run%d.txt",
runNumber);
531 pFile = fopen (textfilename ,
"w");
533 printf(
"Cannot open output text file\n");
538 pFile = fopen (textfilename ,
"w");
540 printf(
"Cannot open output text file\n");
549 printf(
"Determining mean cell hits...\n ");
552 printf(
"Set threshold at %4.3f counts per cell for LB range. \n",(MeanHits+(nsigmaHits*rmsHits))*((
double)nlb_corr-(
double)nBadLBs));
556 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);
559 int EventCount=0., qcount=0.;
566 rootfilename.Form(
"Output/BadCells_run%d.root",
runNumber);
569 TFile*
fout =
new TFile(rootfilename,
"RECREATE");
577 Float_t meanECell = 0;
579 TH1F* h1_lb =
nullptr;
580 TH1F* h1_e =
nullptr;
581 TH1F* h1_q =
nullptr;
582 TH2F* h2_elb =
nullptr;
583 TH2D* h2_pulse =
nullptr;
584 TH2D* h2_t_LB =
nullptr;
586 TH1F* h1_ADCmax =
nullptr;
587 TTree* tree_cells =
nullptr;
589 tree_cells =
new TTree(
"BadTree",
"LAr Tree ordered by OnlID");
590 tree_cells->Branch(
"Energy",&h1_e);
591 tree_cells->Branch(
"Quality",&h1_q);
592 tree_cells->Branch(
"LB",&h1_lb);
593 tree_cells->Branch(
"larid",&larid);
594 tree_cells->Branch(
"onlid",&onlid);
595 tree_cells->Branch(
"Run",&t_run);
596 tree_cells->Branch(
"selectionFlag",&selectionFlag);
597 tree_cells->Branch(
"nhits_ensig",&n_ensig);
598 tree_cells->Branch(
"nhits_e1GeV",&n_ecut);
599 tree_cells->Branch(
"Quality_4k",&fr_q4k);
600 tree_cells->Branch(
"LBfraction",&fr_LB);
601 tree_cells->Branch(
"EnergyVsLB",&h2_elb);
602 tree_cells->Branch(
"PulseShape",&h2_pulse);
603 tree_cells->Branch(
"PulsePeakTimeVsLB",&h2_t_LB);
604 tree_cells->Branch(
"noise",&
noise);
605 tree_cells->Branch(
"MeanCellEnergy",&meanECell);
606 tree_cells->Branch(
"PulseShapeTProf",&TProf_pulse);
607 tree_cells->Branch(
"ADCmax",&h1_ADCmax);
612 TH2D* NormPulse =
new TH2D(
"Normalised_pulse",
"",100,-200,200,400,-10,10); NormPulse->SetXTitle(
"Time [ns]"); NormPulse->SetYTitle(
"Value [ADC counts] / ADCmax");
613 TH2D** Pulsemaps =
new TH2D*[npl];
616 TH2D** t_LBmaps =
new TH2D*[npl];
617 TH1F** CellsFlagged_LB_part =
new TH1F*[npl];
618 TH1F** Emean_NbrEvents_part =
new TH1F*[npl];
619 TH2F* E_LB =
new TH2F(
"E_LB",
"",nlb,
lbmin,
lbmax,ne,emin,emax); E_LB->SetXTitle(
"LB"); E_LB->SetYTitle(
"Energy [GeV]");
620 TH2D* t_LB=
new TH2D(
"t_LB",
"",nlb,
lbmin,
lbmax,400,-200,200); t_LB->GetXaxis()->SetTitle(
"LB"); t_LB->GetYaxis()->SetTitle(
"Time(maxSample) + ofcTime [ns]");
621 TH1F* CF_LB =
new TH1F(
"CF_LB",
"",nlb,
lbmin,
lbmax); CF_LB->GetXaxis()->SetTitle(
"LB"); CF_LB->GetYaxis()->SetTitle(
"Number of cells flagged");
622 TH1F* ME_EV =
new TH1F(
"ME_EV",
"",1000,0,1000); ME_EV->GetXaxis()->SetTitle(
"Nbr events >1 GeV"); ME_EV->GetYaxis()->SetTitle(
"Mean Energy [GeV]");
627 for (
int ii=0;ii<npl;ii++){
629 Cellmaps[ii]->GetXaxis()->SetTitle(
"#eta"); Cellmaps[ii]->GetYaxis()->SetTitle(
"#Phi");
631 Pulsemaps[ii] = (TH2D*)NormPulse->Clone(
hname);
633 E_LBmaps[ii] = (
TH2F*)E_LB->Clone(
hname);
635 t_LBmaps[ii] = (TH2D*)t_LB->Clone(
hname);
637 CellsFlagged_LB_part[ii] = (
TH1F*)CF_LB->Clone(
hname);
640 Emean_NbrEvents_part[ii] = (
TH1F*)ME_EV->Clone(
hname);
645 Eeta->GetXaxis()->SetTitle(
"#eta"); Eeta->GetYaxis()->SetTitle(
"Energy [GeV]");
646 TH1F* CellEnergy =
new TH1F(
"CellEnergy",
"",ne,emin,emax);
647 CellEnergy->GetXaxis()->SetTitle(
"Cell Mean Energy [GeV]"); CellEnergy->GetYaxis()->SetTitle(
"Cells");
648 TH1F* LBfrac =
new TH1F(
"LBfrac",
"",100,0.,1.);
649 LBfrac->GetXaxis()->SetTitle(
"LB fraction"); LBfrac->GetYaxis()->SetTitle(
"Cells");
650 TH1F* Qfrac =
new TH1F(
"Qfrac",
"",100.,0.,1.);
651 Qfrac->GetXaxis()->SetTitle(
"Fraction of Events Q>4000"); Qfrac->GetYaxis()->SetTitle(
"Cells");
653 CellsFlagged_LB->GetXaxis()->SetTitle(
"LB"); CellsFlagged_LB->GetYaxis()->SetTitle(
"Number of cells flagged");
654 TH1F* Emean_NbrEvents =
new TH1F(
"Emean_NEvents_1GeV",
"",1000,0,1000);
655 Emean_NbrEvents->GetXaxis()->SetTitle(
"N events > 1 GeV"); Emean_NbrEvents->GetYaxis()->SetTitle(
"Mean Energy [GeV]");
660 printf(
"Creating histograms (nEbins = %d, nQbins = %d, nLBbins = %d)... ",ne,nq,nlb);
661 TH1F hene =
TH1F(
"hEne",
"",ne,emin,emax); hene.SetXTitle(
"Energy [GeV]"); hene.SetYTitle(
"Events");
662 TH1F hqua =
TH1F(
"hQua",
"",nq/100.,qmin,qmax); hqua.SetXTitle(
"Quality"); hqua.SetYTitle(
"Events");
663 TH1F hlb =
TH1F(
"hLB",
"",nlb,
lbmin,
lbmax); hlb.SetXTitle(
"LB"); hlb.SetYTitle(
"Events");
664 TH2F henelb =
TH2F(
"hEnelb",
"",nlb,
lbmin,
lbmax,ne,emin,emax); henelb.SetXTitle(
"LB"); henelb.SetYTitle(
"Energy [GeV]");
665 TH2D hpulse = TH2D(
"hPulse",
"",100,-200,200,400, -10,10); hpulse.SetXTitle(
"Time [ns]"); hpulse.SetYTitle(
"Value [ADC counts] / ADCmax");
666 TProfile TProfpulse =
TProfile(
"",
"",5, 0, 5,
"s"); TProfpulse.SetXTitle(
"Sample Number"); TProfpulse.SetYTitle(
"Value [ADC counts]");
667 TH2D ht_LB= TH2D(
"ht_LB",
"",nlb,
lbmin,
lbmax,400,-200,200); ht_LB.GetXaxis()->SetTitle(
"LB"); ht_LB.GetYaxis()->SetTitle(
"Time(maxSample) + ofcTime [ns]");
668 TH1F hADCmax =
TH1F(
"",
"",110,-200,2000); hADCmax.SetXTitle(
"ADCmax [ADC counts]"); hADCmax.SetYTitle(
"Events");
675 unsigned int nchannels = tuple->
nChannels();
677 fprintf(
outputFile,
"Onlid \t MeanCellHitBased \t EnergyThresholdBased \t OR \t AND\n");
680 printf(
"Begin looping over cells...\n");
689 double eta = cellInfo->
eta();
690 double phi = cellInfo->
phi();
695 onlid = cellInfo->
onlid();
703 unsigned int ndigits =
hist->nData();
704 if(ndigits<(
unsigned int)
threshold){ nskipped++;
continue; }
710 hname.Form(
"h0x%x_Energy",onlid);
712 hname.Form(
"h0x%x_Quality",onlid);
714 hname.Form(
"h0x%x_LB",onlid);
716 hname.Form(
"h0x%x_Energy_LB",onlid);
718 hname.Form(
"h0x%x_Pulse",onlid);
719 h2_pulse = (TH2D*)hpulse.Clone(
hname);
720 hname.Form(
"TProf0x%x_Pulse",onlid);
722 hname.Form(
"h0x%x_ADCmax",onlid);
724 hname.Form(
"h0x%x_t_LB",onlid);
725 h2_t_LB = (TH2D*)ht_LB.Clone(
hname);
728 bool CellInList =
false;
729 int LBFlaggedIn = -1;
732 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
736 if(!Evdata){ noEvdata++;
continue; }
745 if (isBadLB)
continue;
752 h1_e->Fill(
data->energy()*
GeV);
753 h1_q->Fill(
data->quality());
756 h1_ADCmax->Fill(
data->adcMax());
759 for(
unsigned int isample=0;isample<
data->nSamples();isample++) {
760 TProf_pulse->Fill(isample,
data->value(isample));
761 h2_pulse->Fill(
data->time(isample)-
time,
data->value(isample)/
data->adcMax());
766 if (
data->quality() > qthreshold){
774 if (
data->energy()*
GeV > ecut) {
776 EventEnergySum+=(
data->energy()*
GeV);
781 selectionFlag=
CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
782 if(selectionFlag==2){
786 else if(selectionFlag==1){
796 else if(selectionFlag==0){
813 if ( CellCounter%1000 == 0 ) {
814 printf(
"Scanning Cell No. = %d\n",CellCounter);
816 selectionFlag=
CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
819 if (selectionFlag >= 0){
820 bool ListOutput =
true;
821 int MeanCellHitsSelection = 0;
822 int EnergyThresholdSelection = 0;
823 int OR_selection = 0;
824 int AND_selection = 0;
826 if (selectionFlag == 2){
827 MeanCellHitsSelection = 1;
828 EnergyThresholdSelection = 1;
834 else if (selectionFlag == 1){
835 if (
m_Algo == 1) ListOutput =
false;
836 MeanCellHitsSelection = 1;
837 EnergyThresholdSelection = 0;
842 else if (selectionFlag == 0){
843 if (
m_Algo == 0) ListOutput =
false;
844 MeanCellHitsSelection = 0;
845 EnergyThresholdSelection = 1;
851 iMeanCellHitsSelection += MeanCellHitsSelection;
852 iEnergyThresholdSelection += EnergyThresholdSelection;
853 iOR_selection += OR_selection;
854 iAND_selection += AND_selection;
861 fprintf(
outputFile,
"0x%x \t %d \t \t \t %d \t \t %d \t %d\n",onlid,MeanCellHitsSelection,EnergyThresholdSelection,OR_selection,AND_selection);
869 if (h1_lb->GetBinContent(
ilb) > 0.)
nLB++;
872 fr_LB = (
double)
nLB/(
double)(nlb_corr-nBadLBs);
874 meanECell = h1_e->GetMean();
875 n_ensig = EventCount;
877 fr_q4k = (
double)qcount/(
double)EventCount;
879 n_ecut = nEvents_E_gt_ecut;
885 if (h1_lb->GetBinContent(
ilb) > 0.)
nLB++;
887 fr_LB = (
double)
nLB/(
double)(nlb_corr-nBadLBs);
889 meanECell = h1_e->GetMean();
890 n_ensig = EventCount;
892 fr_q4k = (
double)qcount/(
double)EventCount;
894 n_ecut = nEvents_E_gt_ecut;
896 Eeta->Fill(
eta,meanECell);
897 CellEnergy->Fill(meanECell);
900 Emean_NbrEvents->Fill(n_ecut,EventEnergySum/n_ecut);
901 Emean_NbrEvents_part[
index]->Fill(n_ecut,EventEnergySum/n_ecut);
905 Pulsemaps[
index]->Add(Pulsemaps[
index],h2_pulse);
907 t_LBmaps[
index]->Add(t_LBmaps[
index],h2_t_LB);
908 CellsFlagged_LB_part[
index]->Fill(LBFlaggedIn);
909 CellsFlagged_LB->Fill(LBFlaggedIn);
918 if (nEvents_E_gt_ecut>0) { EventEnergyMean = EventEnergySum/nEvents_E_gt_ecut; }
919 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);
926 if(h1_e)
delete h1_e;
927 if(h1_q)
delete h1_q;
928 if(h1_lb)
delete h1_lb;
929 if(h2_elb)
delete h2_elb;
930 if(h2_pulse)
delete h2_pulse;
931 if(TProf_pulse)
delete TProf_pulse;
932 if(h1_ADCmax)
delete h1_ADCmax;
933 if(h2_t_LB)
delete h2_t_LB;
936 nEvents_E_gt_ecut = 0.;
938 EventEnergyMean = 0.;
940 printf(
"Cells selected: MeanCellHitsCut \t EnergyThreshodCut \t Either List \t Both Lists\n");
941 printf(
"\t \t %d \t \t \t \t %d \t \t %d \t %d\n",iMeanCellHitsSelection,iEnergyThresholdSelection,iOR_selection,iAND_selection);
948 printf(
"Skipped %d cells. Cells scanned = %d\n",nskipped,nbad);
949 printf(
"Cells put in Bad Cell List = %d\n",nBadCells);
950 if (noEvdata>0) printf(
"No Evdata pointer: %d\n",noEvdata);
955 for(
int i=0;
i<npl;
i++){
960 if(CellsFlagged_LB_part[
i]->
GetEntries()>0) CellsFlagged_LB_part[
i]->Write();
961 if( Emean_NbrEvents_part[
i]->
GetEntries()>0) Emean_NbrEvents_part[
i]->Write();
963 CellsFlagged_LB->Write();
967 Emean_NbrEvents->Write();
983 unsigned int nchannels = tuple->
nChannels();
986 qmin = emin =
lbmin = 99999;
987 qmax = emax =
lbmax = -1;
997 unsigned int ndigits =
hist->nData();
998 if(ndigits==0){
continue; }
1001 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
1004 if(!Evdata)
continue;
1008 if (
data->energy() != 0. &&
data->noise() != 0.){
1011 if (isBadLB)
continue;
1017 double energyGeV =
data->energy()/1000.;
1018 if(energyGeV > emax) emax = energyGeV;
1019 if(energyGeV < emin) emin = energyGeV;
1020 if(
data->quality() > qmax) qmax =
data->quality();
1021 if(
data->quality() < qmin) qmin =
data->quality();
1031 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventCells_tot;
1032 std::map< std::pair<unsigned int, unsigned int>,
double > eventEnergy_tot;
1033 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventLumiBlock;
1034 std::map< std::pair<std::string, unsigned int>,
unsigned int > eventCells;
1036 double NSIG = nsigma;
1040 unsigned int nchannels = tuple->
nChannels();
1051 int calo = cellInfo->
calo();
1060 unsigned int ndigits =
hist->nData();
1061 if(ndigits==0){
continue; }
1064 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
1068 if(!Evdata)
continue;
1073 double energyGeV =
data->energy()/1000.;
1076 if (lumiBlock <= lbmax && lumiBlock >=
lbmin ) {
1077 if (
data->energy() > sqrt((NSIG*
data->noise())*(NSIG*
data->noise())) &&
data->energy() != 0. &&
data->noise() != 0.){
1079 std::pair<unsigned int, unsigned int>
ev(
data->run(),
data->event());
1080 std::pair<std::string, unsigned int> ev_cryo(
GetCryostat(calo),
data->event());
1081 eventCells_tot[
ev]++;
1082 eventEnergy_tot[
ev] += energyGeV;
1083 eventLumiBlock[
ev] =
lb;
1084 eventCells[ev_cryo]++;
1093 std::map< std::string, TProfile*> hCellsEvLB;
1094 std::map< std::string, TH1F*> hp_cryo;
1099 std::string Cryo[8] = {
"EMBA",
"EMBC",
"EMECA",
"EMECC",
"HECA",
"HECC",
"FCALA",
"FCALC"};
1100 for(
int k = 0;
k< nCryo;
k++){
1101 hCellsEvLB[Cryo[
k]] = (
TProfile*)hCells_Ev_LB->Clone();
1103 for(
unsigned int ievent = 0; ievent < tuple->
nEvents(); ievent++) {
1105 std::pair<unsigned int, unsigned int>
id(evtData->
run(), evtData->
event());
1106 if(eventCells_tot[
id]==0)
continue;
1108 for(
int j=0; j<nCryo ;j++){
1109 std::pair<std::string, unsigned int> id_cryo(Cryo[j], evtData->
event());
1110 hCellsEvLB[Cryo[j]]->Fill(eventLumiBlock[
id], eventCells[id_cryo]);
1113 hCells_Ev_LB->Fill(eventLumiBlock[
id], eventCells_tot[
id]);
1116 TH1F* hp =
new TH1F(
"",
"",(hCells_Ev_LB->GetBinContent(hCells_Ev_LB->GetMaximumBin())),0,hCells_Ev_LB->GetBinContent(hCells_Ev_LB->GetMaximumBin()));
1117 for(
int k = 0;
k< nCryo;
k++){
1118 TH1F* hp_temp =
new TH1F(
"",
"",(hCellsEvLB[Cryo[
k]]->GetBinContent(hCellsEvLB[Cryo[
k]]->GetMaximumBin())*10),0,hCellsEvLB[Cryo[
k]]->GetBinContent(hCellsEvLB[Cryo[
k]]->GetMaximumBin()));
1119 std::string
name1 =
"hp_" + Cryo[
k];
1120 hp_cryo[Cryo[
k]] = (
TH1F*)hp_temp->Clone(
name1.c_str());
1123 for (
int p=1;
p<=nlb;
p++){
1124 for(
int j=0; j<nCryo; j++){
1125 double bincont = hCellsEvLB[Cryo[j]]->GetBinContent(
p);
1126 hp_cryo[Cryo[j]]->Fill(bincont);
1128 double bincont1 = hCells_Ev_LB->GetBinContent(
p);
1132 std::map< std::string, double>
Mean;
1133 std::map< std::string, double>
RMS;
1134 std::map< int, vector<std::string> > BadLB_cryo;
1136 for(
int k=0;
k<nCryo;
k++)
1138 Mean[Cryo[
k]] = hp_cryo[Cryo[
k]]->GetMean();
1139 RMS[Cryo[
k]] = hp_cryo[Cryo[
k]]->GetRMS();
1142 double MEAN2 = hp->GetMean();
1143 double RMS2 = hp->GetRMS();
1148 vector<int> BadLB_tot=DQLBList;
1150 double value1 = (MEAN2+(RMS2*3.));
1151 int numberFlagged=0;
1152 for (
int l=1;
l<=nlb;
l++){
1153 for(
int j=0; j<nCryo; j++){
1154 double value = MEAN2 + (
RMS[Cryo[j]]*3);
1155 if(hCellsEvLB[Cryo[j]]->GetBinContent(
l) >
value){
1156 int badLB = hCellsEvLB[Cryo[j]]->GetBinLowEdge(
l);
1157 BadLB_cryo[badLB].push_back(Cryo[j]);
1159 for(
unsigned int p=0;
p< BadLB_tot.size();
p++) {
1160 if(badLB == BadLB_tot[
p]) isadd =
false;
1163 BadLB_tot.push_back(badLB);
1165 printf(
"Removing bad LB %d in %s\n",badLB,Cryo[j].c_str());
1170 if (hCells_Ev_LB->GetBinContent(
l) > value1){
1171 BadLB.push_back(hCells_Ev_LB->GetBinLowEdge(
l));
1175 printf(
"Number of Bad LBs found: %d\n",numberFlagged);
1176 printf(
"Number of LBs to be removed: %zu", BadLB_tot.size());
1186 for (
unsigned int yy=0;
yy<BadLBList.size();
yy++){
1202 unsigned int nchannels = tuple->
nChannels();
1205 double TotalRecordedHits=0;
1206 std::vector<int, std::allocator<int> > HitsPerLB;
1220 unsigned int ndigits =
hist->nData();
1223 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
1226 if(!Evdata)
continue;
1232 if (isBadLB)
continue;
1243 if (nHits > 0.)
nCells++;
1256 TH1F* tp_ev =
new TH1F(
"",
"",th1_Hits->GetBinContent(th1_Hits->GetMaximumBin())*100,0,th1_Hits->GetBinContent(th1_Hits->GetMaximumBin()));
1258 for (
int i=1;
i<=nlb;
i++){
1259 if(hNLB->GetBinContent(
i)>0){
1261 tp_ev->Fill(th1_Hits->GetBinContent(
i));
1262 TotalRecordedHits+=th1_Hits->GetBinContent(
i);
1263 HitsPerLB.push_back(th1_Hits->GetBinContent(
i));
1267 MeanHits = ((
double)TotalRecordedHits/(
double)nlb_corr)/(
double)
nCells;
1271 for (
unsigned int j=0;j<HitsPerLB.size();j++){
1275 rmsHits =
var/nlb_corr;
1277 printf(
"Mean number of hits/cell for 1 LB = %4.3f, RMS = %4.3f\n",MeanHits,rmsHits);
1278 printf(
"Number of cells firing at E > %d sigma = %d\n",nsigma,
nCells);
1279 printf(
"Total number of LBs included = %d\n",nlb_corr);
1282 if (th1_Hits)
delete th1_Hits;
1292 }
else if (lumiBlock <= lbmax && lumiBlock >=
lbmin &&
energy > sqrt((nsigma*
noise)*(nsigma*
noise)) ) {
1315 if ((EventCount > (MeanHits+(nsigmaHits*rmsHits))*((
double)nlb-(
double)nBadLBs)))
a=
true;
1330 if ((
a &&
b) || (
a &&
c)){
1332 }
else if (
a && !
b && !
c){
1334 }
else if ((!
a &&
b) || (!
a &&
c)){
1366 bool isFoundPartition =
false;
1367 for(
int i=0;
i< npl;
i++){
1371 isFoundPartition =
true;
1374 if(!isFoundPartition){
1375 printf(
"ERROR: Partition %s does not exist! Invalid input name.\n",
partname.c_str());
1376 printf(
"Possible candidates are:\n");
1377 for (
int j=0;j<npl;j++){ printf(
"%s\n",
m_LarIdTranslator->GetPartitonLayerName(j)); }
1385 std::string cryostat =
"NotGiven";
1386 if(calo == 1) cryostat =
"EMBA";
1387 else if(calo == -1) cryostat =
"EMBC";
1388 else if(calo == 5) cryostat =
"FCALA";
1389 else if(calo == -5) cryostat =
"FCALC";
1390 else if(calo == 4) cryostat =
"HECA";
1391 else if(calo == -4) cryostat =
"HECC";
1392 else if(calo == 2 || calo == 3) cryostat =
"EMECA";
1393 else if(calo == -2 || calo == -3) cryostat =
"EMECC";
1416 std::vector<int> LBList;
1418 std::ifstream
infile(LBfile.Data());
1425 if(
list->GetEntries() == 0){
1426 printf(
"No LB filtering specified, or bad format. Exiting.\n");
1428 LBList.push_back(0);
1432 for(
int k = 0;
k <
list->GetEntries();
k++){
1433 TObjString* tobs = (TObjString*)(
list->At(
k));
1434 LBList.push_back((
int)(tobs->String()).Atoi());
1437 printf(
"LB List: %d\n",(
int)LBList.size());
1476 ULong64_t onlid = 0;
1477 std::map<ULong64_t,unsigned int>* idmap =
new std::map<ULong64_t,unsigned int>;
1479 int nskipped=0,nrepeated=0;
1484 unsigned int nchannels = tuple->
nChannels();
1497 onlid = cellInfo->
onlid();
1503 if(onlid<=0) printf(
"%u: Bad Cell Onlid = 0x%x (%+.2f,%+.2f)\n",
ichan,(
unsigned int)onlid,cellInfo->
eta(),cellInfo->
phi());
1505 idmap_itr = idmap->find(onlid);
1508 if(idmap_itr == idmap->end()){
1509 (*idmap)[onlid] =
ichan;
1517 printf(
"Skipped %d cells.\n",nskipped);
1518 printf(
"Number of onlids: Unique=%d, Repeated=%d\n",(
int)idmap->size(),nrepeated);
1537 printf(
"LArCellsEmptyMonitoring::DoEtaPhiMonitoring: Reading options: %s %s.\n",optionplot,optionsave);
1539 bool knormplot =
false;
1540 if(strstr(optionplot,
"Precentage")) knormplot =
true;
1543 double nsigmacut = 0.;
1544 double qualitycut = 0.;
1545 if(strstr(optionplot,
"Sigma")){
1547 if(strstr(optionplot,
"4Sigma")) nsigmacut = 4.;
1548 if(strstr(optionplot,
"5Sigma")) nsigmacut = 5.;
1549 }
else if(strstr(optionplot,
"Quality")){
1557 ULong64_t onlid = 0;
1562 unsigned int nchannels = tuple->
nChannels();
1569 for(
int j=0;j<
nhists;j++){
1572 hmap_counts_all[j]->SetName(
hname);
1576 hmap_energy_cut[j]->SetName(
hname);
1580 hmap_quality_cut[j]->SetName(
hname);
1591 unsigned int ndigits =
hist->nData();
1592 if(ndigits==0){ nskipped++;
continue; }
1597 onlid = cellInfo->
onlid();
1598 if(onlid<=0) printf(
"%u: Bad Cell Onlid = 0x%x (%+.2f,%+.2f)\n",
ichan,(
unsigned int)onlid,cellInfo->
eta(),cellInfo->
phi());
1601 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
1609 if(kcuttype==1 &&
data->noise()>0.){
1619 printf(
"Skipped %d cells. Bad Cells = %d\n",nskipped,nbad);
1623 for(
int j=0;j<
nhists;j++){
1624 hmap_energy_cut[j]->Divide(hmap_counts_all[j]);
1625 hmap_quality_cut[j]->Divide(hmap_counts_all[j]);
1629 TFile* tout =
nullptr;
1630 TCanvas*
c0 =
nullptr, *
c1 =
nullptr;
1632 if(!strcmp(optionsave,
"root")){
1633 tout =
new TFile(
"EtaPhiMonitoring.root",
"recreate");
1635 c0->SetName(
"Normalization");
1637 for(
int j=0;j<
nhists;j++) hmap_counts_all[j]->Write();
1640 c0->SetName(
"EnergyCut");
1642 for(
int j=0;j<
nhists;j++) hmap_energy_cut[j]->Write();
1645 c0 =
m_LarIdTranslator->CaloPartitionLayerDisplay((TH1**)hmap_quality_cut,
"QualityCut",1);
1646 c0->SetName(
"QualityCut");
1648 for(
int j=0;j<
nhists;j++) hmap_quality_cut[j]->Write();
1650 tout->Close();
delete tout;
1653 c0->SaveAs(
"Normalization.png");
1654 c1 =
new TCanvas(
"c1",
"");
1655 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); }
1658 c0->SaveAs(
"EnergyCut.png");
1660 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); }
1663 c0 =
m_LarIdTranslator->CaloPartitionLayerDisplay((TH1**)hmap_quality_cut,
"QualityCut",1);
1664 c0->SaveAs(
"QualityCut.png");
1666 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); }
1687 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventCells;
1688 std::map< std::pair<unsigned int, unsigned int>,
double > eventEnergy_PS;
1689 std::map< std::pair<unsigned int, unsigned int>,
double > eventEnergy_LAr;
1690 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventCells_PS;
1691 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventLayer;
1695 unsigned int nchannels = tuple->
nChannels();
1707 unsigned int ndigits =
hist->nData();
1708 if(ndigits==0){
continue; }
1711 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
1716 double energyGeV =
data->energy()/1000.;
1719 std::pair<unsigned int, unsigned int>
ev(
data->run(),
data->event());
1721 eventEnergy_LAr[
ev] += energyGeV;
1724 eventEnergy_PS[
ev] += energyGeV;
1725 eventCells_PS[
ev]++;
1733 const int nbins = 42;
1734 double* binE =
new double[
nbins+1];
1737 double de = 0.,ene = 0.;
1765 TH1F* E_LAr_EM5_pass =
new TH1F(
"E_LAr_EM5_pass",
"",
nbins,binE);
1766 TH1F* E_LAr_EM5_tot =
new TH1F(
"E_LAr_EM5_tot",
"",
nbins,binE);
1767 TH1F* E_LAr_TAU8_pass =
new TH1F(
"E_LAr_TAU8_pass",
"",
nbins,binE);
1768 TH1F* E_LAr_TAU8_tot =
new TH1F(
"E_LAr_TAU8_tot",
"",
nbins,binE);
1769 TH1F* E_LAr_J10_pass =
new TH1F(
"E_LAr_J10_pass",
"",
nbins,binE);
1770 TH1F* E_LAr_J10_tot =
new TH1F(
"E_LAr_J10_tot",
"",
nbins,binE);
1771 TH1F* nCellsPS_pass =
new TH1F(
"nCellsPS_pass",
"",100, 0, 100);
1772 TH1F* nCellsPS_tot =
new TH1F(
"nCellsPS_tot",
"",100, 0, 100);
1773 TH1F* nCellsPS_EM5_pass =
new TH1F(
"nCellsPS_EM5_pass",
"",100, 0, 100);
1774 TH1F* nCellsPS_EM5_tot =
new TH1F(
"nCellsPS_EM5_tot",
"",100, 0, 100);
1775 TH1F* nCellsPS_TAU8_pass =
new TH1F(
"nCellsPS_TAU8_pass",
"",100, 0, 100);
1776 TH1F* nCellsPS_TAU8_tot =
new TH1F(
"nCellsPS_TAU8_tot",
"",100, 0, 100);
1777 TH1F* nCellsPS_J10_tot =
new TH1F(
"nCellsPS_J10_tot",
"",100, 0, 100);
1778 TH1F* nCellsPS_J10_pass =
new TH1F(
"nCellsPS_J10_pass",
"",100, 0, 100);
1783 for(
unsigned int ievent = 0; ievent < tuple->
nEvents(); ievent++) {
1785 std::pair<unsigned int, unsigned int>
id(evtData->
run(), evtData->
event());
1786 if(eventCells[
id]==0)
continue;
1787 if(eventEnergy_LAr[
id]<=0)
continue;
1788 bool isEM5 = evtData->
isPassed(
"L1_EM5_EMPTY");
1789 bool isTAU8 = evtData->
isPassed(
"L1_TAU8_EMPTY");
1790 bool isJ10 = evtData->
isPassed(
"L1_J10_EMPTY");
1792 E_LAr_tot->Fill(eventEnergy_LAr[
id]);
1793 nCellsPS_tot->Fill(eventCells_PS[
id]);
1795 E_LAr_EM5_tot->Fill(eventEnergy_LAr[
id]);
1796 nCellsPS_EM5_tot->Fill(eventCells_PS[
id]);
1799 E_LAr_TAU8_tot->Fill(eventEnergy_LAr[
id]);
1800 nCellsPS_TAU8_tot->Fill(eventCells_PS[
id]);
1803 E_LAr_J10_tot->Fill(eventEnergy_LAr[
id]);
1804 nCellsPS_J10_tot->Fill(eventCells_PS[
id]);
1808 if((eventEnergy_LAr[
id])!=0)
1809 ratio = eventEnergy_PS[
id]/(eventEnergy_LAr[
id]);
1810 bool isL2_PreS =
false;
1811 if(
ratio > fractionInPS) isL2_PreS =
true;
1814 E_LAr_pass->Fill(eventEnergy_LAr[
id]);
1815 nCellsPS_pass->Fill(eventCells_PS[
id]);
1817 E_LAr_EM5_pass->Fill(eventEnergy_LAr[
id]);
1818 nCellsPS_EM5_pass->Fill(eventCells_PS[
id]);
1821 E_LAr_TAU8_pass->Fill(eventEnergy_LAr[
id]);
1822 nCellsPS_TAU8_pass->Fill(eventCells_PS[
id]);
1825 E_LAr_J10_pass->Fill(eventEnergy_LAr[
id]);
1826 nCellsPS_J10_pass->Fill(eventCells_PS[
id]);
1832 sprintf(
fname,
"TriggerEfficiency_%d.root",
run);
1833 TFile* mfile =
new TFile(
fname,
"recreate");
1834 E_LAr_pass->Write();
1836 E_LAr_EM5_pass->Write();
1837 E_LAr_EM5_tot->Write();
1838 E_LAr_TAU8_pass->Write();
1839 E_LAr_TAU8_tot->Write();
1840 E_LAr_J10_pass->Write();
1841 E_LAr_J10_tot->Write();
1842 nCellsPS_pass->Write();
1843 nCellsPS_tot->Write();
1844 nCellsPS_EM5_pass->Write();
1845 nCellsPS_EM5_tot->Write();
1846 nCellsPS_TAU8_pass->Write();
1847 nCellsPS_TAU8_tot->Write();
1848 nCellsPS_J10_pass->Write();
1849 nCellsPS_J10_tot->Write();
1851 mfile->Close();
delete mfile;