DAQ/macros/Efficiency_ModuleDisplay.C

0001 ///////////////////////////////////////////////////////////
0002 //         Jennifer James, Charles Hughes                //
0003 //          derived from files created by                //
0004 //         Aditya Dash and Thomas Marshall               //
0005 //                  May 19, 2023                         //
0006 ///////////////////////////////////////////////////////////
0007
0008 // includes
0009 #ifndef __CINT__
0010 #include <algorithm>
0011 #include <cmath>
0012 #include "TCanvas.h"
0013 #include "TApplication.h"
0014 #include "TH1D.h"
0015 #include "TH1F.h"
0016 #include "TH2D.h"
0017 #include "TH2F.h"
0018 #include "TH3D.h"
0019 #include "TH3F.h"
0020 #include "TString.h"
0021 #include "TStyle.h"
0022 #include "TFile.h"
0023 #include "TGraphPolar.h"
0024 #include "TGraphPolargram.h"
0025 #include "TAxis.h"
0026 #include "TLatex.h"
0027 #include "TLegend.h"
0028 #include "TApplication.h"
0029 #include "TCanvas.h"
0030 #include "TMath.h"
0031 #include "TVirtualFitter.h"
0032 #include "Math/MinimizerOptions.h"
0033 #include "TFitResultPtr.h"
0034 #include "TFitResult.h"
0035 #include "TF1.h"
0036 #include "TPaveLabel.h"
0037 #include <string.h>
0038 #include <vector>
0039 #include <iostream>
0040 #include <RtypesCore.h>
0041 #endif
0042
0043 using namespace std;
0044
0045 void Locate(Int_t id, Bool_t is_ASIDE, Double_t *rbin, Double_t *thbin);
0046
0047
0048 void Efficiency_ModuleDisplay(){
0049
0050   bool skip_empty_fees =true;  // true will skip FEEs with 0 channel output and false will keep them in the plot/analysis
0051
0052   std::vector<pair<int,int>> vec;
0053
0054   const TString filename3( Form( "./pedestal-26540-outfile.root") ); // change the run number as needed
0055
0056   TFile *infile3 = TFile::Open(filename3);
0057
0058   if(!infile3) return;
0059
0060   TNtuple * liveTntuple ;
0061   liveTntuple = (TNtuple*) infile3->Get("h_Alive");
0062
0063   TNtuple * totTntuple ;
0064   totTntuple = (TNtuple*) infile3->Get("h_AliveTot");
0065
0066   TH3F * dm2 = new TH3F("dm2","TPC map",26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
0067   liveTntuple->Draw("module_id:sec_id:fee_id>>dm2","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
0068
0069   TH3F * h_AlivePedestalNoise = new TH3F("h_AlivePedestalNoise","TPC Alive Pedestal Noise of Std. Dev,",26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
0070   liveTntuple->Draw("module_id:sec_id:fee_id>>h_AlivePedestalNoise","pedStdi",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
0071
0072   TH3F * tot = new TH3F("tot"," TPC Tot Map", 26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
0073   totTntuple->Draw("module_id:sec_id:fee_id>>tot","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
0074
0075   dm2->Print();
0076   tot->Print();
0077
0078   std::vector<Float_t> sub_arrA;
0079   std::vector<Float_t> sub_arrC;
0080
0081   Float_t frac[dm2->GetNbinsX()][dm2->GetNbinsY()]; // array to store fractions
0082   std::vector<Float_t> module_std_dev;
0083
0084   for (Int_t i = 0; i < dm2->GetNbinsY(); i++) { // i is looping over sec ID
0085     for (Int_t j = 0; j < dm2->GetNbinsZ(); j++) { // j is looping over Module ID
0086
0087       Float_t num=0; // numerator (live)
0088       Float_t denom=0; // denominator (total)
0089       Float_t sum_noise=0; // getting noise if needed -- it is not plotted
0090
0091       for (Int_t k = 0; k < dm2->GetNbinsX(); k++) { // k is looping over FEE ID
0092
0093         if( skip_empty_fees){ //if you want to skip empty fees
0094           if( dm2->GetBinContent(k+1,i+1,j+1) >= 1){ //only add to numerator and denominator if FEE has at least 1 live channel
0095             num+=dm2->GetBinContent(k+1,i+1,j+1);
0096             denom+=tot->GetBinContent(k+1,i+1,j+1);
0097             sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
0098           }
0099         }
0100         else { //if you don't want to skip empty fees
0101       num+=dm2->GetBinContent(k+1,i+1,j+1);
0102       denom+=tot->GetBinContent(k+1,i+1,j+1);
0103       sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
0104     }
0105       }
0106
0107       Float_t frac_val= ( num / denom) * 100.0; // calculate the live channel fraction
0108       Float_t noise_value = sum_noise/num;
0109
0110       frac[i+1][j+1] = frac_val; // store fraction in array
0111       std::cout << "Sec ID = " << i << ", Module ID = " << j+1 << ", Live fraction = " << frac_val << "%" << std::endl;
0112       if (i < 12) {
0113     sub_arrC.push_back(frac_val); //North Side
0114       }else {
0115     sub_arrA.push_back(frac_val); //South Side
0116       }
0117     }
0118   }
0119   // std::cout << "Live Fraction =" << frac_val << std::endl;
0120   // std::cout <<  " Size A = " << sub_arrA.size() << std::endl;
0121   // std::cout <<  " Size C = " << sub_arrC.size() << std::endl;
0122
0123   ////////////////////////////////////////////////////////////////////////
0124   //       A Side   South Label Conventions                             //
0125   ////////////////////////////////////////////////////////////////////////
0126   //   1 - 12_R1   16 - 17_R1   31 - 22_R1
0127   //   2 - 12_R2   17 - 17_R2   32 - 22_R2
0128   //   3 - 12_R3   18 - 17_R3   33 - 22_R3
0129   //   4 - 13_R1   19 - 18_R1   34 - 23_R1
0130   //   5 - 13_R2   20 - 18_R2   35 - 23_R2
0131   //   6 - 13_R3   21 - 18_R3   36 - 23_R3
0132   //   7 - 14_R1   22 - 19_R1
0133   //   8 - 14_R2   23 - 19_R2
0134   //   9 - 14_R3   24 - 19_R3
0135   //  10 - 15_R1   25 - 20_R1
0136   //  11 - 15_R2   26 - 20_R2
0137   //  12 - 15_R3   27 - 20_R3
0138   //  13 - 16_R1   28 - 21_R1
0139   //  14 - 16_R2   29 - 21_R2
0140   //  15 - 16_R3   30 - 21_R3
0141
0142   ////////////////////////////////////////////////////////////////////////
0143   //       C Side North Label Conventions                               //
0144   ////////////////////////////////////////////////////////////////////////
0145   //   1 - 00_R1   16 - 05_R1   31 - 10_R1
0146   //   2 - 00_R2   17 - 05_R2   32 - 10_R2
0147   //   3 - 00_R3   18 - 05_R3   33 - 10_R3
0148   //   4 - 01_R1   19 - 06_R1   34 - 11_R1
0149   //   5 - 01_R2   20 - 06_R2   35 - 11_R2
0150   //   6 - 01_R3   21 - 06_R3   36 - 11_R3
0151   //   7 - 02_R1   22 - 07_R1
0152   //   8 - 02_R2   23 - 07_R2
0153   //   9 - 02_R3   24 - 07_R3
0154   //  10 - 03_R1   25 - 08_R1
0155   //  11 - 03_R2   26 - 08_R2
0156   //  12 - 03_R3   27 - 08_R3
0157   //  13 - 04_R1   28 - 09_R1
0158   //  14 - 04_R2   29 - 09_R2
0159   //  15 - 04_R3   30 - 09_R3
0160
0161
0162   gStyle->SetOptStat(0);
0163
0164   const Int_t N_rBins = 4; //(one inner bin NOT to be filled, 2nd bin is R1, 3rd bin is R2, 4th bin is R3)
0165   const Int_t N_thBins = 12; //(12 theta bins of uniform angle (360/12 = 30 degrees = TMath::Pi()/6 ~= 0.523 rad) )
0166
0167   Double_t rBin_edges[N_rBins + 1] = {0.0, 0.256, 0.504, 0.752, 1.00}; //variable edges for the radial dimensions
0168
0169   TGraphPolargram* polardummy1 = new TGraphPolargram("polardummy1",0,1,0,2.*TMath::Pi()); //dummy plots to get the canvas right (not to be filled)
0170   polardummy1->SetToGrad();
0171   polardummy1->SetNdivPolar(N_thBins);
0172   polardummy1->SetLineColor(0);
0173   polardummy1->SetRadialLabelSize(0);
0174
0175   TGraphPolargram* polardummy2 = new TGraphPolargram("polardummy2",0,1,0,2.*TMath::Pi());
0176   polardummy2->SetToGrad();
0177   polardummy2->SetNdivPolar( N_thBins);
0178   polardummy2->SetLineColor(0);
0179   polardummy2->SetRadialLabelSize(0);
0180
0181   for(Int_t i = 0 ; i < 12 ; i++){ //setting the axis label (CCW from horizontal right axis)
0182     char labelstr1[128];
0183     char labelstr2[128];
0184     if(i<=9){ // i -> [0:9]
0185       sprintf(labelstr2,"C0%d",i);
0186
0187       if(i<=6){ // i -> [0:6] (halfway)
0188     sprintf(labelstr1,"A%d",18-i);
0189       }
0190       else if(i>6){ // i -> [7:9]
0191     sprintf(labelstr1,"A%d",30-i);
0192       }
0193
0194     }
0195     else { // i -> [10:11]
0196       sprintf(labelstr2,"C%d",i);
0197       sprintf(labelstr1,"A%d",30-i);
0198     }
0199     polardummy1->SetPolarLabel(i,labelstr1);
0200     polardummy2->SetPolarLabel(i,labelstr2);
0201   }
0202
0203   TH2D* ErrASide = new TH2D( "ASide" , "ADC Counts North Side" , N_thBins, -TMath::Pi()/12. , 23.*TMath::Pi()/12. , N_rBins , rBin_edges ); // X maps to theta, Y maps to R
0204
0205   TH2D* ErrCSide = new TH2D( "CSide" , "ADC Counts South Side" , N_thBins, -TMath::Pi()/12. , 23.*TMath::Pi()/12. , N_rBins , rBin_edges ); // X maps to theta, Y maps to R
0206
0207
0208   Double_t r, theta;
0209   Int_t trip_count_total = 0;
0210   Bool_t is_ASIDE = true;
0211
0212   for (Int_t i = 0; i < 36; i++) {
0213     Locate(i + 1, true, &r, &theta);
0214     ErrASide->Fill(theta, r, sub_arrA[i]);
0215     // cout<<"Region # A "<<(i)<<" Alive Fraction = "<<sub_arrA[i]<<endl;
0216   }
0217
0218   for (Int_t i = 0; i < 36; i++) {
0219     Locate(i + 1,false, &r, &theta);
0220     ErrCSide->Fill(theta, r, sub_arrC[i]);
0221     // cout<<"Region # C "<<(i)<<" Alive Fraction = "<<sub_arrC[i]<<endl;
0222   }
0223
0224   // change the titles as needed
0225   TH2D* dummy_his1 = new TH2D("dummy1", "26540-Alive Channel Fraction North Side (%)", 100, -1.5, 1.5, 100, -1.5, 1.5); //dummy histos for titles
0226   TH2D* dummy_his2 = new TH2D("dummy2", "26540-Alive Channel Fraction South Side (%)", 100, -1.5, 1.5, 100, -1.5, 1.5);
0227   //TPaveLabels for sector labels
0228   TPaveLabel* A00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "18" );
0229   TPaveLabel* A01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "17" );
0230   TPaveLabel* A02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "16" );
0231   TPaveLabel* A03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "15" );
0232   TPaveLabel* A04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "14" );
0233   TPaveLabel* A05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "13" );
0234   TPaveLabel* A06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "12" );
0235   TPaveLabel* A07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "23" );
0236   TPaveLabel* A08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "22" );
0237   TPaveLabel* A09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "21" );
0238   TPaveLabel* A10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "20" );
0239   TPaveLabel* A11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "19" );
0240
0241   TPaveLabel* C00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "00" );
0242   TPaveLabel* C01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "01" );
0243   TPaveLabel* C02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "02" );
0244   TPaveLabel* C03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "03" );
0245   TPaveLabel* C04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "04" );
0246   TPaveLabel* C05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "05" );
0247   TPaveLabel* C06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "06" );
0248   TPaveLabel* C07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "07" );
0249   TPaveLabel* C08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "08" );
0250   TPaveLabel* C09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "09" );
0251   TPaveLabel* C10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "10" );
0252   TPaveLabel* C11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "11" );
0253
0254   A00->SetFillColor(0);
0255   A01->SetFillColor(0);
0256   A02->SetFillColor(0);
0257   A03->SetFillColor(0);
0258   A04->SetFillColor(0);
0259   A05->SetFillColor(0);
0260   A06->SetFillColor(0);
0261   A07->SetFillColor(0);
0262   A08->SetFillColor(0);
0263   A09->SetFillColor(0);
0264   A10->SetFillColor(0);
0265   A11->SetFillColor(0);
0266
0267   C00->SetFillColor(0);
0268   C01->SetFillColor(0);
0269   C02->SetFillColor(0);
0270   C03->SetFillColor(0);
0271   C04->SetFillColor(0);
0272   C05->SetFillColor(0);
0273   C06->SetFillColor(0);
0274   C07->SetFillColor(0);
0275   C08->SetFillColor(0);
0276   C09->SetFillColor(0);
0277   C10->SetFillColor(0);
0278   C11->SetFillColor(0);
0279
0280   gStyle->SetPalette(kBird);
0281
0282   TCanvas *Error_Viz = new TCanvas("Error_Viz", "Error_Viz", 1248, 598);
0283   Error_Viz->Divide(2,1);
0284   Error_Viz->cd(1);
0285   dummy_his1->Draw("");
0286   //polardummy1->Draw("same");
0287   ErrCSide->Draw("colpolzsame0");
0288   C00->Draw("same");
0289   C01->Draw("same");
0290   C02->Draw("same");
0291   C03->Draw("same");
0292   C04->Draw("same");
0293   C05->Draw("same");
0294   C06->Draw("same");
0295   C07->Draw("same");
0296   C08->Draw("same");
0297   C09->Draw("same");
0298   C10->Draw("same");
0299   C11->Draw("same");
0300   Error_Viz->cd(2);
0301   dummy_his2->Draw("");
0302   //polardummy2->Draw("");
0303   ErrASide->Draw("colpolzsame0");
0304   A00->Draw("same");
0305   A01->Draw("same");
0306   A02->Draw("same");
0307   A03->Draw("same");
0308   A04->Draw("same");
0309   A05->Draw("same");
0310   A06->Draw("same");
0311   A07->Draw("same");
0312   A08->Draw("same");
0313   A09->Draw("same");
0314   A10->Draw("same");
0315   A11->Draw("same");
0316
0317   ErrCSide->SetMaximum(100);
0318   ErrASide->SetMaximum(100);
0319
0320   //change minimum live fraction value as needed
0321   ErrCSide->SetMinimum(80);
0322   ErrASide->SetMinimum(80);
0323
0324
0325   //Set Same Scale for A and C side displays
0326
0327   Double_t Maxval = TMath::Max(ErrASide->GetBinContent(ErrASide->GetMaximumBin()),ErrCSide->GetBinContent(ErrCSide->GetMaximumBin()));
0328
0329   TFile *outf = new TFile("Trip_Histos.root","RECREATE");
0330   Error_Viz->Write();
0331
0332   outf->Write();
0333 }
0334
0335 void Locate(Int_t id, Bool_t is_ASIDE, Double_t *rbin, Double_t *thbin) {
0336
0337   Double_t ASIDE_angle_bins[12] = { 0.1*2.*TMath::Pi()/12 , 1.1*2.*TMath::Pi()/12 , 2.1*2.*TMath::Pi()/12 , 3.1*2.*TMath::Pi()/12 , 4.1*2.*TMath::Pi()/12 , 5.1*2.*TMath::Pi()/12 , 6.1*2.*TMath::Pi()/12 , 7.1*2.*TMath::Pi()/12 , 8.1*2.*TMath::Pi()/12 , 9.1*2.*TMath::Pi()/12 , 10.1*2.*TMath::Pi()/12 , 11.1*2.*TMath::Pi()/12 }; //CCW from x = 0 (RHS horizontal)
0338
0339   Double_t CSIDE_angle_bins[12] = { 6.1*2.*TMath::Pi()/12 , 5.1*2.*TMath::Pi()/12 , 4.1*2.*TMath::Pi()/12 , 3.1*2.*TMath::Pi()/12 , 2.1*2.*TMath::Pi()/12 , 1.1*2.*TMath::Pi()/12 , 0.1*2.*TMath::Pi()/12 , 11.1*2.*TMath::Pi()/12 , 10.1*2.*TMath::Pi()/12 , 9.1*2.*TMath::Pi()/12 , 8.1*2.*TMath::Pi()/12 , 7.1*2.*TMath::Pi()/12  }; //CCW from x = 0 (RHS horizontal)
0340
0341   Int_t modid3 = id % 3;
0342
0343   switch(modid3) {
0344   case 1:
0345     *rbin = 0.4; //R1
0346     break;
0347   case 2:
0348     *rbin = 0.6; //R2
0349     break;
0350   case 0:
0351     *rbin = 0.8; //R3
0352     break;
0353   }
0354
0355
0356   if( is_ASIDE ){
0357     *thbin = CSIDE_angle_bins[TMath::FloorNint((id-1)/3)];
0358   }
0359   else{
0360     *thbin = ASIDE_angle_bins[TMath::FloorNint((id-1)/3)];
0361   }
0362
0363 }