Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Efficiency_ModuleDisplay.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Efficiency_ModuleDisplay.C
1 
2 // Jennifer James, Charles Hughes //
3 // derived from files created by //
4 // Aditya Dash and Thomas Marshall //
5 // May 19, 2023 //
7 
8 // includes
9 #ifndef __CINT__
10 #include <algorithm>
11 #include <cmath>
12 #include "TCanvas.h"
13 #include "TApplication.h"
14 #include "TH1D.h"
15 #include "TH1F.h"
16 #include "TH2D.h"
17 #include "TH2F.h"
18 #include "TH3D.h"
19 #include "TH3F.h"
20 #include "TString.h"
21 #include "TStyle.h"
22 #include "TFile.h"
23 #include "TGraphPolar.h"
24 #include "TGraphPolargram.h"
25 #include "TAxis.h"
26 #include "TLatex.h"
27 #include "TLegend.h"
28 #include "TApplication.h"
29 #include "TCanvas.h"
30 #include "TMath.h"
31 #include "TVirtualFitter.h"
32 #include "Math/MinimizerOptions.h"
33 #include "TFitResultPtr.h"
34 #include "TFitResult.h"
35 #include "TF1.h"
36 #include "TPaveLabel.h"
37 #include <string.h>
38 #include <vector>
39 #include <iostream>
40 #include <RtypesCore.h>
41 #endif
42 
43 using namespace std;
44 
45 void Locate(Int_t id, Bool_t is_ASIDE, Double_t *rbin, Double_t *thbin);
46 
47 
48 void Efficiency_ModuleDisplay(){
49 
50  bool skip_empty_fees =true; // true will skip FEEs with 0 channel output and false will keep them in the plot/analysis
51 
52  std::vector<pair<int,int>> vec;
53 
54  const TString filename3( Form( "./pedestal-26540-outfile.root") ); // change the run number as needed
55 
56  TFile *infile3 = TFile::Open(filename3);
57 
58  if(!infile3) return;
59 
60  TNtuple * liveTntuple ;
61  liveTntuple = (TNtuple*) infile3->Get("h_Alive");
62 
63  TNtuple * totTntuple ;
64  totTntuple = (TNtuple*) infile3->Get("h_AliveTot");
65 
66  TH3F * dm2 = new TH3F("dm2","TPC map",26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
67  liveTntuple->Draw("module_id:sec_id:fee_id>>dm2","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
68 
69  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 );
70  liveTntuple->Draw("module_id:sec_id:fee_id>>h_AlivePedestalNoise","pedStdi",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
71 
72  TH3F * tot = new TH3F("tot"," TPC Tot Map", 26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
73  totTntuple->Draw("module_id:sec_id:fee_id>>tot","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
74 
75  dm2->Print();
76  tot->Print();
77 
78  std::vector<Float_t> sub_arrA;
79  std::vector<Float_t> sub_arrC;
80 
81  Float_t frac[dm2->GetNbinsX()][dm2->GetNbinsY()]; // array to store fractions
82  std::vector<Float_t> module_std_dev;
83 
84  for (Int_t i = 0; i < dm2->GetNbinsY(); i++) { // i is looping over sec ID
85  for (Int_t j = 0; j < dm2->GetNbinsZ(); j++) { // j is looping over Module ID
86 
87  Float_t num=0; // numerator (live)
88  Float_t denom=0; // denominator (total)
89  Float_t sum_noise=0; // getting noise if needed -- it is not plotted
90 
91  for (Int_t k = 0; k < dm2->GetNbinsX(); k++) { // k is looping over FEE ID
92 
93  if( skip_empty_fees){ //if you want to skip empty fees
94  if( dm2->GetBinContent(k+1,i+1,j+1) >= 1){ //only add to numerator and denominator if FEE has at least 1 live channel
95  num+=dm2->GetBinContent(k+1,i+1,j+1);
96  denom+=tot->GetBinContent(k+1,i+1,j+1);
97  sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
98  }
99  }
100  else { //if you don't want to skip empty fees
101  num+=dm2->GetBinContent(k+1,i+1,j+1);
102  denom+=tot->GetBinContent(k+1,i+1,j+1);
103  sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
104  }
105  }
106 
107  Float_t frac_val= ( num / denom) * 100.0; // calculate the live channel fraction
108  Float_t noise_value = sum_noise/num;
109 
110  frac[i+1][j+1] = frac_val; // store fraction in array
111  std::cout << "Sec ID = " << i << ", Module ID = " << j+1 << ", Live fraction = " << frac_val << "%" << std::endl;
112  if (i < 12) {
113  sub_arrC.push_back(frac_val); //North Side
114  }else {
115  sub_arrA.push_back(frac_val); //South Side
116  }
117  }
118  }
119  // std::cout << "Live Fraction =" << frac_val << std::endl;
120  // std::cout << " Size A = " << sub_arrA.size() << std::endl;
121  // std::cout << " Size C = " << sub_arrC.size() << std::endl;
122 
124  // A Side South Label Conventions //
126  // 1 - 12_R1 16 - 17_R1 31 - 22_R1
127  // 2 - 12_R2 17 - 17_R2 32 - 22_R2
128  // 3 - 12_R3 18 - 17_R3 33 - 22_R3
129  // 4 - 13_R1 19 - 18_R1 34 - 23_R1
130  // 5 - 13_R2 20 - 18_R2 35 - 23_R2
131  // 6 - 13_R3 21 - 18_R3 36 - 23_R3
132  // 7 - 14_R1 22 - 19_R1
133  // 8 - 14_R2 23 - 19_R2
134  // 9 - 14_R3 24 - 19_R3
135  // 10 - 15_R1 25 - 20_R1
136  // 11 - 15_R2 26 - 20_R2
137  // 12 - 15_R3 27 - 20_R3
138  // 13 - 16_R1 28 - 21_R1
139  // 14 - 16_R2 29 - 21_R2
140  // 15 - 16_R3 30 - 21_R3
141 
143  // C Side North Label Conventions //
145  // 1 - 00_R1 16 - 05_R1 31 - 10_R1
146  // 2 - 00_R2 17 - 05_R2 32 - 10_R2
147  // 3 - 00_R3 18 - 05_R3 33 - 10_R3
148  // 4 - 01_R1 19 - 06_R1 34 - 11_R1
149  // 5 - 01_R2 20 - 06_R2 35 - 11_R2
150  // 6 - 01_R3 21 - 06_R3 36 - 11_R3
151  // 7 - 02_R1 22 - 07_R1
152  // 8 - 02_R2 23 - 07_R2
153  // 9 - 02_R3 24 - 07_R3
154  // 10 - 03_R1 25 - 08_R1
155  // 11 - 03_R2 26 - 08_R2
156  // 12 - 03_R3 27 - 08_R3
157  // 13 - 04_R1 28 - 09_R1
158  // 14 - 04_R2 29 - 09_R2
159  // 15 - 04_R3 30 - 09_R3
160 
161 
162  gStyle->SetOptStat(0);
163 
164  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)
165  const Int_t N_thBins = 12; //(12 theta bins of uniform angle (360/12 = 30 degrees = TMath::Pi()/6 ~= 0.523 rad) )
166 
167  Double_t rBin_edges[N_rBins + 1] = {0.0, 0.256, 0.504, 0.752, 1.00}; //variable edges for the radial dimensions
168 
169  TGraphPolargram* polardummy1 = new TGraphPolargram("polardummy1",0,1,0,2.*TMath::Pi()); //dummy plots to get the canvas right (not to be filled)
170  polardummy1->SetToGrad();
171  polardummy1->SetNdivPolar(N_thBins);
172  polardummy1->SetLineColor(0);
173  polardummy1->SetRadialLabelSize(0);
174 
175  TGraphPolargram* polardummy2 = new TGraphPolargram("polardummy2",0,1,0,2.*TMath::Pi());
176  polardummy2->SetToGrad();
177  polardummy2->SetNdivPolar( N_thBins);
178  polardummy2->SetLineColor(0);
179  polardummy2->SetRadialLabelSize(0);
180 
181  for(Int_t i = 0 ; i < 12 ; i++){ //setting the axis label (CCW from horizontal right axis)
182  char labelstr1[128];
183  char labelstr2[128];
184  if(i<=9){ // i -> [0:9]
185  sprintf(labelstr2,"C0%d",i);
186 
187  if(i<=6){ // i -> [0:6] (halfway)
188  sprintf(labelstr1,"A%d",18-i);
189  }
190  else if(i>6){ // i -> [7:9]
191  sprintf(labelstr1,"A%d",30-i);
192  }
193 
194  }
195  else { // i -> [10:11]
196  sprintf(labelstr2,"C%d",i);
197  sprintf(labelstr1,"A%d",30-i);
198  }
199  polardummy1->SetPolarLabel(i,labelstr1);
200  polardummy2->SetPolarLabel(i,labelstr2);
201  }
202 
203  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
204 
205  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
206 
207 
208  Double_t r, theta;
209  Int_t trip_count_total = 0;
210  Bool_t is_ASIDE = true;
211 
212  for (Int_t i = 0; i < 36; i++) {
213  Locate(i + 1, true, &r, &theta);
214  ErrASide->Fill(theta, r, sub_arrA[i]);
215  // cout<<"Region # A "<<(i)<<" Alive Fraction = "<<sub_arrA[i]<<endl;
216  }
217 
218  for (Int_t i = 0; i < 36; i++) {
219  Locate(i + 1,false, &r, &theta);
220  ErrCSide->Fill(theta, r, sub_arrC[i]);
221  // cout<<"Region # C "<<(i)<<" Alive Fraction = "<<sub_arrC[i]<<endl;
222  }
223 
224  // change the titles as needed
225  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
226  TH2D* dummy_his2 = new TH2D("dummy2", "26540-Alive Channel Fraction South Side (%)", 100, -1.5, 1.5, 100, -1.5, 1.5);
227  //TPaveLabels for sector labels
228  TPaveLabel* A00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "18" );
229  TPaveLabel* A01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "17" );
230  TPaveLabel* A02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "16" );
231  TPaveLabel* A03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "15" );
232  TPaveLabel* A04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "14" );
233  TPaveLabel* A05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "13" );
234  TPaveLabel* A06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "12" );
235  TPaveLabel* A07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "23" );
236  TPaveLabel* A08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "22" );
237  TPaveLabel* A09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "21" );
238  TPaveLabel* A10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "20" );
239  TPaveLabel* A11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "19" );
240 
241  TPaveLabel* C00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "00" );
242  TPaveLabel* C01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "01" );
243  TPaveLabel* C02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "02" );
244  TPaveLabel* C03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "03" );
245  TPaveLabel* C04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "04" );
246  TPaveLabel* C05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "05" );
247  TPaveLabel* C06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "06" );
248  TPaveLabel* C07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "07" );
249  TPaveLabel* C08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "08" );
250  TPaveLabel* C09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "09" );
251  TPaveLabel* C10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "10" );
252  TPaveLabel* C11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "11" );
253 
254  A00->SetFillColor(0);
255  A01->SetFillColor(0);
256  A02->SetFillColor(0);
257  A03->SetFillColor(0);
258  A04->SetFillColor(0);
259  A05->SetFillColor(0);
260  A06->SetFillColor(0);
261  A07->SetFillColor(0);
262  A08->SetFillColor(0);
263  A09->SetFillColor(0);
264  A10->SetFillColor(0);
265  A11->SetFillColor(0);
266 
267  C00->SetFillColor(0);
268  C01->SetFillColor(0);
269  C02->SetFillColor(0);
270  C03->SetFillColor(0);
271  C04->SetFillColor(0);
272  C05->SetFillColor(0);
273  C06->SetFillColor(0);
274  C07->SetFillColor(0);
275  C08->SetFillColor(0);
276  C09->SetFillColor(0);
277  C10->SetFillColor(0);
278  C11->SetFillColor(0);
279 
280  gStyle->SetPalette(kBird);
281 
282  TCanvas *Error_Viz = new TCanvas("Error_Viz", "Error_Viz", 1248, 598);
283  Error_Viz->Divide(2,1);
284  Error_Viz->cd(1);
285  dummy_his1->Draw("");
286  //polardummy1->Draw("same");
287  ErrCSide->Draw("colpolzsame0");
288  C00->Draw("same");
289  C01->Draw("same");
290  C02->Draw("same");
291  C03->Draw("same");
292  C04->Draw("same");
293  C05->Draw("same");
294  C06->Draw("same");
295  C07->Draw("same");
296  C08->Draw("same");
297  C09->Draw("same");
298  C10->Draw("same");
299  C11->Draw("same");
300  Error_Viz->cd(2);
301  dummy_his2->Draw("");
302  //polardummy2->Draw("");
303  ErrASide->Draw("colpolzsame0");
304  A00->Draw("same");
305  A01->Draw("same");
306  A02->Draw("same");
307  A03->Draw("same");
308  A04->Draw("same");
309  A05->Draw("same");
310  A06->Draw("same");
311  A07->Draw("same");
312  A08->Draw("same");
313  A09->Draw("same");
314  A10->Draw("same");
315  A11->Draw("same");
316 
317  ErrCSide->SetMaximum(100);
318  ErrASide->SetMaximum(100);
319 
320  //change minimum live fraction value as needed
321  ErrCSide->SetMinimum(80);
322  ErrASide->SetMinimum(80);
323 
324 
325  //Set Same Scale for A and C side displays
326 
327  Double_t Maxval = TMath::Max(ErrASide->GetBinContent(ErrASide->GetMaximumBin()),ErrCSide->GetBinContent(ErrCSide->GetMaximumBin()));
328 
329  TFile *outf = new TFile("Trip_Histos.root","RECREATE");
330  Error_Viz->Write();
331 
332  outf->Write();
333 }
334 
335 void Locate(Int_t id, Bool_t is_ASIDE, Double_t *rbin, Double_t *thbin) {
336 
337  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)
338 
339  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)
340 
341  Int_t modid3 = id % 3;
342 
343  switch(modid3) {
344  case 1:
345  *rbin = 0.4; //R1
346  break;
347  case 2:
348  *rbin = 0.6; //R2
349  break;
350  case 0:
351  *rbin = 0.8; //R3
352  break;
353  }
354 
355 
356  if( is_ASIDE ){
357  *thbin = CSIDE_angle_bins[TMath::FloorNint((id-1)/3)];
358  }
359  else{
360  *thbin = ASIDE_angle_bins[TMath::FloorNint((id-1)/3)];
361  }
362 
363 }