425 TString htitle,
hname,textfilename,rootfilename;
437 printf(
"WARNING: Received instruction to select recurring bad cells only in SetSelectRecurringBadCells(bool flag)\n");
438 printf(
"-------> Only one algorithm selected for cell selection criteria in SetAlgo(int algoindex)\n");
439 printf(
"-------> Recurring cell selection overides algorithm specification\n");
444 printf(
"ERROR: Partition not specified in SetPartition(bool set_part, std::string partname)\n");
451 int selectionFlag = -1;
452 const double GeV = 1/1000.;
462 unsigned int onlid = 0;
474 int nEvents_E_gt_ecut=0;
475 double EventEnergySum=0.;
476 double EventEnergyMean=0.;
479 double MeanHits=0, rmsHits=0;
480 int iMeanCellHitsSelection=0,iEnergyThresholdSelection=0,iOR_selection=0,iAND_selection=0;
486 vector<int, std::allocator<int> > DQLBList;
487 vector<int, std::allocator<int> > BadLBList;
489 printf(
"Reading in DQ Defect LB List...\n");
494 printf(
"Getting bad LB list...\n");
496 printf(
"...Done.\n");
497 nBadLBs = BadLBList.size();
500 BadLBList = std::move(DQLBList);
501 nBadLBs = BadLBList.size();
505 printf(
"Getting histogram limits... ");
516 printf(
"Done.\n LBmin: %d LBmax: %d Emin: %d Emax: %d Qmin: %d Qmax: %d\n",
lbmin,
lbmax,emin,emax,qmin,qmax);
521 textfilename.Form(
"Output/BadCellSelection_run%d.txt",
runNumber);
524 std::cout<<(
"Cannot open output selection summary file " + textfilename)<<
"\n";
530 std::cout<<(
"Cannot open output selection summary file " + textfilename)<<
"\n";
536 FILE * pFile=
nullptr;
539 textfilename.Form(
"Output/BadCellList_run%d.txt",
runNumber);
540 pFile = fopen (textfilename ,
"w");
542 printf(
"Cannot open output text file\n");
547 pFile = fopen (textfilename ,
"w");
549 printf(
"Cannot open output text file\n");
558 printf(
"Determining mean cell hits...\n ");
561 printf(
"Set threshold at %4.3f counts per cell for LB range. \n",(MeanHits+(nsigmaHits*rmsHits))*((
double)nlb_corr-(
double)nBadLBs));
565 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);
568 int EventCount=0., qcount=0.;
575 rootfilename.Form(
"Output/BadCells_run%d.root",
runNumber);
578 auto fout = std::make_unique<TFile>(rootfilename,
"RECREATE");
586 Float_t meanECell = 0;
588 TH1F* h1_lb =
nullptr;
589 TH1F* h1_e =
nullptr;
590 TH1F* h1_q =
nullptr;
591 TH2F* h2_elb =
nullptr;
592 TH2D* h2_pulse =
nullptr;
593 TH2D* h2_t_LB =
nullptr;
595 TH1F* h1_ADCmax =
nullptr;
596 std::unique_ptr<TTree> tree_cells;
598 tree_cells.reset(
new TTree(
"BadTree",
"LAr Tree ordered by OnlID"));
599 tree_cells->Branch(
"Energy",&h1_e);
600 tree_cells->Branch(
"Quality",&h1_q);
601 tree_cells->Branch(
"LB",&h1_lb);
602 tree_cells->Branch(
"larid",&larid);
603 tree_cells->Branch(
"onlid",&onlid);
604 tree_cells->Branch(
"Run",&t_run);
605 tree_cells->Branch(
"selectionFlag",&selectionFlag);
606 tree_cells->Branch(
"nhits_ensig",&n_ensig);
607 tree_cells->Branch(
"nhits_e1GeV",&n_ecut);
608 tree_cells->Branch(
"Quality_4k",&fr_q4k);
609 tree_cells->Branch(
"LBfraction",&fr_LB);
610 tree_cells->Branch(
"EnergyVsLB",&h2_elb);
611 tree_cells->Branch(
"PulseShape",&h2_pulse);
612 tree_cells->Branch(
"PulsePeakTimeVsLB",&h2_t_LB);
613 tree_cells->Branch(
"noise",&
noise);
614 tree_cells->Branch(
"MeanCellEnergy",&meanECell);
615 tree_cells->Branch(
"PulseShapeTProf",&TProf_pulse);
616 tree_cells->Branch(
"ADCmax",&h1_ADCmax);
621 TH2D NormPulse(
"Normalised_pulse",
"",100,-200,200,400,-10,10);
622 NormPulse.SetXTitle(
"Time [ns]");
623 NormPulse.SetYTitle(
"Value [ADC counts] / ADCmax");
634 E_LB.SetXTitle(
"LB");
635 E_LB.SetYTitle(
"Energy [GeV]");
637 TH2D t_LB(
"t_LB",
"",nlb,
lbmin,
lbmax,400,-200,200);
638 t_LB.GetXaxis()->SetTitle(
"LB");
639 t_LB.GetYaxis()->SetTitle(
"Time(maxSample) + ofcTime [ns]");
642 CF_LB.GetXaxis()->SetTitle(
"LB");
643 CF_LB.GetYaxis()->SetTitle(
"Number of cells flagged");
645 TH1F ME_EV(
"ME_EV",
"",1000,0,1000);
646 ME_EV.GetXaxis()->SetTitle(
"Nbr events >1 GeV");
647 ME_EV.GetYaxis()->SetTitle(
"Mean Energy [GeV]");
652 for (
int ii=0;ii<npl;ii++){
654 Cellmaps[ii]->GetXaxis()->SetTitle(
"#eta"); Cellmaps[ii]->GetYaxis()->SetTitle(
"#Phi");
656 Pulsemaps[ii].reset((TH2D*)NormPulse.Clone(
hname));
658 E_LBmaps[ii].reset((
TH2F*)E_LB.Clone(
hname));
660 t_LBmaps[ii].reset((TH2D*)t_LB.Clone(
hname));
662 CellsFlagged_LB_part[ii].reset((
TH1F*)CF_LB.Clone(
hname));
665 Emean_NbrEvents_part[ii].reset((
TH1F*)ME_EV.Clone(
hname));
670 Eeta.GetXaxis()->SetTitle(
"#eta");
671 Eeta.GetYaxis()->SetTitle(
"Energy [GeV]");
673 TH1F CellEnergy(
"CellEnergy",
"",ne,emin,emax);
674 CellEnergy.GetXaxis()->SetTitle(
"Cell Mean Energy [GeV]");
675 CellEnergy.GetYaxis()->SetTitle(
"Cells");
677 TH1F LBfrac(
"LBfrac",
"",100,0.,1.);
678 LBfrac.GetXaxis()->SetTitle(
"LB fraction");
679 LBfrac.GetYaxis()->SetTitle(
"Cells");
681 TH1F Qfrac(
"Qfrac",
"",100.,0.,1.);
682 Qfrac.GetXaxis()->SetTitle(
"Fraction of Events Q>4000");
683 Qfrac.GetYaxis()->SetTitle(
"Cells");
686 CellsFlagged_LB.GetXaxis()->SetTitle(
"LB");
687 CellsFlagged_LB.GetYaxis()->SetTitle(
"Number of cells flagged");
689 TH1F Emean_NbrEvents(
"Emean_NEvents_1GeV",
"",1000,0,1000);
690 Emean_NbrEvents.GetXaxis()->SetTitle(
"N events > 1 GeV");
691 Emean_NbrEvents.GetYaxis()->SetTitle(
"Mean Energy [GeV]");
696 printf(
"Creating histograms (nEbins = %d, nQbins = %d, nLBbins = %d)... ",ne,nq,nlb);
697 TH1F hene =
TH1F(
"hEne",
"",ne,emin,emax); hene.SetXTitle(
"Energy [GeV]");
698 hene.SetYTitle(
"Events");
699 TH1F hqua =
TH1F(
"hQua",
"",nq/100.,qmin,qmax); hqua.SetXTitle(
"Quality");
700 hqua.SetYTitle(
"Events");
702 hlb.SetYTitle(
"Events");
704 henelb.SetYTitle(
"Energy [GeV]");
705 TH2D hpulse = TH2D(
"hPulse",
"",100,-200,200,400, -10,10); hpulse.SetXTitle(
"Time [ns]");
706 hpulse.SetYTitle(
"Value [ADC counts] / ADCmax");
707 TProfile TProfpulse =
TProfile(
"",
"",5, 0, 5,
"s"); TProfpulse.SetXTitle(
"Sample Number");
708 TProfpulse.SetYTitle(
"Value [ADC counts]");
709 TH2D ht_LB= TH2D(
"ht_LB",
"",nlb,
lbmin,
lbmax,400,-200,200); ht_LB.GetXaxis()->SetTitle(
"LB");
710 ht_LB.GetYaxis()->SetTitle(
"Time(maxSample) + ofcTime [ns]");
711 TH1F hADCmax =
TH1F(
"",
"",110,-200,2000); hADCmax.SetXTitle(
"ADCmax [ADC counts]"); hADCmax.SetYTitle(
"Events");
718 unsigned int nchannels = tuple->
nChannels();
720 fprintf(
outputFile,
"Onlid \t MeanCellHitBased \t EnergyThresholdBased \t OR \t AND\n");
723 printf(
"Begin looping over cells...\n");
732 double eta = cellInfo->
eta();
733 double phi = cellInfo->
phi();
738 onlid = cellInfo->
onlid();
746 unsigned int ndigits =
hist->nData();
747 if(ndigits<(
unsigned int)
threshold){ nskipped++;
continue; }
753 hname.Form(
"h0x%x_Energy",onlid);
755 hname.Form(
"h0x%x_Quality",onlid);
757 hname.Form(
"h0x%x_LB",onlid);
759 hname.Form(
"h0x%x_Energy_LB",onlid);
761 hname.Form(
"h0x%x_Pulse",onlid);
762 h2_pulse = (TH2D*)hpulse.Clone(
hname);
763 hname.Form(
"TProf0x%x_Pulse",onlid);
765 hname.Form(
"h0x%x_ADCmax",onlid);
767 hname.Form(
"h0x%x_t_LB",onlid);
768 h2_t_LB = (TH2D*)ht_LB.Clone(
hname);
771 bool CellInList =
false;
772 int LBFlaggedIn = -1;
775 for(
unsigned int idigit = 0; idigit < ndigits; idigit++){
779 if(!Evdata){ noEvdata++;
continue; }
788 if (isBadLB)
continue;
795 h1_e->Fill(
data->energy()*
GeV);
796 h1_q->Fill(
data->quality());
799 h1_ADCmax->Fill(
data->adcMax());
802 for(
unsigned int isample=0;isample<
data->nSamples();isample++) {
803 TProf_pulse->Fill(isample,
data->value(isample));
804 h2_pulse->Fill(
data->time(isample)-time,
data->value(isample)/
data->adcMax());
809 if (
data->quality() > qthreshold){
817 if (
data->energy()*
GeV > ecut) {
819 EventEnergySum+=(
data->energy()*
GeV);
824 selectionFlag=
CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
825 if(selectionFlag==2){
829 else if(selectionFlag==1){
839 else if(selectionFlag==0){
856 if ( CellCounter%1000 == 0 ) {
857 printf(
"Scanning Cell No. = %d\n",CellCounter);
859 selectionFlag=
CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
862 if (selectionFlag >= 0){
863 bool ListOutput =
true;
864 int MeanCellHitsSelection = 0;
865 int EnergyThresholdSelection = 0;
866 int OR_selection = 0;
867 int AND_selection = 0;
869 if (selectionFlag == 2){
870 MeanCellHitsSelection = 1;
871 EnergyThresholdSelection = 1;
877 else if (selectionFlag == 1){
878 if (
m_Algo == 1) ListOutput =
false;
879 MeanCellHitsSelection = 1;
880 EnergyThresholdSelection = 0;
885 else if (selectionFlag == 0){
886 if (
m_Algo == 0) ListOutput =
false;
887 MeanCellHitsSelection = 0;
888 EnergyThresholdSelection = 1;
894 iMeanCellHitsSelection += MeanCellHitsSelection;
895 iEnergyThresholdSelection += EnergyThresholdSelection;
896 iOR_selection += OR_selection;
897 iAND_selection += AND_selection;
904 fprintf(
outputFile,
"0x%x \t %d \t \t \t %d \t \t %d \t %d\n",onlid,MeanCellHitsSelection,EnergyThresholdSelection,OR_selection,AND_selection);
912 if (h1_lb->GetBinContent(
ilb) > 0.)
nLB++;
915 fr_LB = (
double)
nLB/(
double)(nlb_corr-nBadLBs);
917 meanECell = h1_e->GetMean();
918 n_ensig = EventCount;
920 fr_q4k = (
double)qcount/(
double)EventCount;
922 n_ecut = nEvents_E_gt_ecut;
928 if (h1_lb->GetBinContent(
ilb) > 0.)
nLB++;
930 fr_LB = (
double)
nLB/(
double)(nlb_corr-nBadLBs);
932 meanECell = h1_e->GetMean();
933 n_ensig = EventCount;
935 fr_q4k = (
double)qcount/(
double)EventCount;
937 n_ecut = nEvents_E_gt_ecut;
938 Cellmaps[
index]->Fill(eta,phi,meanECell);
939 Eeta.Fill(eta,meanECell);
940 CellEnergy.Fill(meanECell);
943 Emean_NbrEvents.Fill(n_ecut,EventEnergySum/n_ecut);
944 Emean_NbrEvents_part[
index]->Fill(n_ecut,EventEnergySum/n_ecut);
951 CellsFlagged_LB_part[
index]->Fill(LBFlaggedIn);
952 CellsFlagged_LB.Fill(LBFlaggedIn);
960 if (nEvents_E_gt_ecut>0) { EventEnergyMean = EventEnergySum/nEvents_E_gt_ecut; }
961 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);
970 nEvents_E_gt_ecut = 0.;
972 EventEnergyMean = 0.;
974 printf(
"Cells selected: MeanCellHitsCut \t EnergyThreshodCut \t Either List \t Both Lists\n");
975 printf(
"\t \t %d \t \t \t \t %d \t \t %d \t %d\n",iMeanCellHitsSelection,iEnergyThresholdSelection,iOR_selection,iAND_selection);
982 printf(
"Skipped %d cells. Cells scanned = %d\n",nskipped,nbad);
983 printf(
"Cells put in Bad Cell List = %d\n",nBadCells);
984 if (noEvdata>0) printf(
"No Evdata pointer: %d\n",noEvdata);
988 for(
int i=0;
i<npl;
i++){
993 if(CellsFlagged_LB_part[
i]->
GetEntries()>0) CellsFlagged_LB_part[
i]->Write();
994 if( Emean_NbrEvents_part[
i]->
GetEntries()>0) Emean_NbrEvents_part[
i]->Write();
996 CellsFlagged_LB.Write();
1000 Emean_NbrEvents.Write();