422{
423
424 TString htitle,
hname,textfilename,rootfilename;
425
426
428
429
430
431
432
433
434
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
443 printf("ERROR: Partition not specified in SetPartition(bool set_part, std::string partname)\n");
444 std::abort();
447 }
448
449
450 int selectionFlag = -1;
451 const double GeV = 1/1000.;
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;
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
484
485 vector<int, std::allocator<int> > DQLBList;
486 vector<int, std::allocator<int> > BadLBList;
488 printf("Reading in DQ Defect LB List...\n");
490 printf("...Done\n");
491 }
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 }
499 BadLBList = std::move(DQLBList);
500 nBadLBs = BadLBList.size();
501 }
502
503
504 printf("Getting histogram limits... ");
505 GetLimits_EqLB(inputfile,lbmin,lbmax,emin,emax,qmin,qmax,runNumber,BadLBList);
509 }
510
511
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
520 textfilename.Form("Output/BadCellSelection_run%d.txt",runNumber);
522 if (!outputFile){
523 std::cout<<("Cannot open output selection summary file " + textfilename)<< "\n";
524 std::abort();}
525 } else {
528 if (!outputFile){
529 std::cout<<("Cannot open output selection summary file " + textfilename)<< "\n";
530 std::abort();}
531 }
532
533
534
535 FILE * pFile=nullptr;
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 {
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
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
570
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;
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;
594 TH1F* h1_ADCmax =
nullptr;
595 std::unique_ptr<TTree> tree_cells;
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
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
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
651 for (int ii=0;ii<npl;ii++){
653 Cellmaps[ii]->GetXaxis()->SetTitle("#eta"); Cellmaps[ii]->GetYaxis()->SetTitle("#Phi");
655 Pulsemaps[ii].reset((TH2D*)NormPulse.Clone(hname));
657 E_LBmaps[ii].reset((TH2F*)E_LB.Clone(hname));
659 t_LBmaps[ii].reset((TH2D*)t_LB.Clone(hname));
661 CellsFlagged_LB_part[ii].reset((TH1F*)CF_LB.Clone(hname));
662
664 Emean_NbrEvents_part[ii].reset((TH1F*)ME_EV.Clone(hname));
665
666 }
667 }
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
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
715 std::unique_ptr<LArSamples::Interface> tuple =
Interface::open(inputfile);
716 printf("Number of events: %u %u\n",tuple->nEvents(),tuple->nChannels());
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
725 if(tuple->historySize(ichan)==0){ nskipped++; continue; }
726
727
728
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
737 onlid = cellInfo->
onlid();
738
739
742 }
743
744
745 unsigned int ndigits =
hist->nData();
746 if(ndigits<(unsigned int)threshold){ nskipped++; continue; }
747
748
749
751
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
774 for(unsigned int idigit = 0; idigit < ndigits; idigit++){
775 const LArSamples::Data*
data =
hist->data(idigit);
777 const LArSamples::EventData& Evdata =
data->eventData();
780
781
782
784
785
786 if (isBadLB) continue;
787
788
790 if (ProcessEvent){
791
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());
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
807 if (
data->quality() > qthreshold){
808 qcount = qcount+1;
809 }
810
811
812 EventCount++;
813
814
815 if (
data->energy()*
GeV > ecut) {
816 nEvents_E_gt_ecut++;
817 EventEnergySum+=(
data->energy()*
GeV);
818 }
819
820
821 if(!CellInList){
822 selectionFlag=
CheckCellSelectionCriteria(EventCount,nsigmaHits,MeanHits,rmsHits,nEvents_E_gt_ecut,EventEnergySum,nBadLBs,nlb_corr);
823 if(selectionFlag==2){
825 CellInList = true;
826 }
827 else if(selectionFlag==1){
830 CellInList = true;
831 }
834 CellInList = true;
835 }
836 }
837 else if(selectionFlag==0){
840 CellInList = true;
841 }
844 CellInList = true;
845 }
846 }
847 }
848
849 }
850
851 }
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
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
892 iMeanCellHitsSelection += MeanCellHitsSelection;
893 iEnergyThresholdSelection += EnergyThresholdSelection;
894 iOR_selection += OR_selection;
895 iAND_selection += AND_selection;
896
897
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++;
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);
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
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);
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;
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
954
955
956
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 }
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);
975 fclose (pFile);
976 }
977 fclose (outputFile);
978
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
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 }
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}
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
std::vector< std::unique_ptr< T > > THVec
TGraphErrors * GetEntries(TH2F *histo)
short feedThrough() const
bool CheckPartition(int index) const
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)
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)
int CheckCellSelectionCriteria(int EventCount, int nsigmaHits, double MeanHits, double rmsHits, int nEvents_E_gt_ecut, double EventEnergySum, int nBadLBs, int nlb) const
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)
std::vector< int, std::allocator< int > > ReadBadLBList(const TString &LBfile)
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
time(flags, cells_name, *args, **kw)
TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
TProfile(*args, **kwargs)
TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)