Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Noise_ModuleDisplay.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Noise_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 Noise_ModuleDisplay(){
49 
50  // std::vector <pair<int,int>> vec1;
51 
52  bool skip_empty_fees = false;
53 
54  std::vector<pair<int,int>> vec;
55 
56  char name[100];
57  int run_num;
58 
59  cout << "Input run number: ";
60  cin >> run_num;
61 
62  const TString filename3( Form( "./pedestal-%d-outfile.root",run_num,run_num) );
63  //const TString filename3( Form( "/sphenix/u/jamesj3j3/workfest_Charles_mistake/sPHENIXProjects/pedestal-10131-outfile.root") );
64 
65  // std::cout << "Analyze - filename2: " << filename2 << std::endl;
66  //
67  TFile *infile3 = TFile::Open(filename3);
68  // TFile* infile2 = TFile::Open(filename2);
69  //
70  if(!infile3) return;
71 
72  TNtuple * liveTntuple ;
73  liveTntuple = (TNtuple*) infile3->Get("h_Alive");
74 
75  TNtuple * totTntuple ;
76  totTntuple = (TNtuple*) infile3->Get("h_AliveTot");
77 
78  TH3F * dm2 = new TH3F("dm2","TPC map",26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
79  liveTntuple->Draw("module_id:sec_id:fee_id>>dm2","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
80 
81  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 );
82  liveTntuple->Draw("module_id:sec_id:fee_id>>h_AlivePedestalNoise","pedStdi * (pedStdi < 20)",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
83 
84  TH3F * tot = new TH3F("tot"," TPC Tot Map", 26, -0.5, 25.5, 24, -0.5, 23.5, 3, 0.5, 3.5 );
85  totTntuple->Draw("module_id:sec_id:fee_id>>tot","",""); //x axis = FEE ID, y axis = SEC ID, z axis = MODULE ID
86  //h3->Draw("colz");
87 
88  dm2->Print();
89  tot->Print();
90 
91  // std::vector<Float_t> frac_val;
92  std::vector<Float_t> sub_arrA;
93  std::vector<Float_t> sub_arrC;
94 
95  Float_t frac[dm2->GetNbinsX()][dm2->GetNbinsY()]; // array to store fractions
96  std::vector<Float_t> module_std_dev;
97 
98  for (Int_t i = 0; i < dm2->GetNbinsY(); i++) { // i is looping over sec ID
99  // if(fee->GetBinContent(i)<1){
100  // continue;
101  // }
102  for (Int_t j = 0; j < dm2->GetNbinsZ(); j++) { // j is looping over Module ID
103 
104  Float_t num=0; // numerator (live)
105  Float_t denom=0; // denominator (total)
106  Float_t sum_noise=0;
107 
108  for (Int_t k = 0; k < dm2->GetNbinsX(); k++) { // k is looping over FEE ID
109 
110  Float_t FEE_noise_sum = 0.0;
111 
112  if(h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1) > 0.0){
113  if( skip_empty_fees){ //if you want to skip empty fees
114  if( dm2->GetBinContent(k+1,i+1,j+1) >= 1){ //only add to numerator and denominator if FEE has at least 1 live channel
115  num+=dm2->GetBinContent(k+1,i+1,j+1);
116  denom+=tot->GetBinContent(k+1,i+1,j+1);
117  sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
118  }
119  }
120  else { //if you don't want to skip empty fees
121  num+=dm2->GetBinContent(k+1,i+1,j+1);
122  denom+=tot->GetBinContent(k+1,i+1,j+1);
123  sum_noise+=h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
124  }
125  }
126  FEE_noise_sum = h_AlivePedestalNoise->GetBinContent(k+1,i+1,j+1);
127  std::cout << "Sec ID = " << i+1 << ", Module ID = " << j+1 << ", FEE ID = " << k+1 << ", Std. Dev. sum = " << FEE_noise_sum << std::endl;
128  }
129 
130  Float_t frac_val= ( num / denom) * 100.0; // calculate the fraction
131  Float_t noise_value = sum_noise/num;
132 
133  frac[i+1][j+1] = frac_val; // store fraction in array
134  std::cout << "Sec ID = " << i+1 << ", Module ID = " << j+1 << ", Live fraction = " << frac_val << "%" << std::endl;
135  std::cout << "Sec ID = " << i+1 << ", Module ID = " << j+1 << ", Avg. Std. Dev. = " << noise_value << std::endl;
136  if (i < 12) {
137  sub_arrC.push_back(noise_value);
138  //std::cout << " Live fraction C = " << sub_arrC.size() << " %" << std::endl;
139  //std::cout << "sub_arrC[" << i << "] = " << frac_val << std::endl;
140  }else {
141  sub_arrA.push_back(noise_value);
142  //std::cout << "sub_arrA[" << i << "] = " << frac_val << std::endl;
143  //std::cout << " Size = " << sub_arrA.size() <<" %" << std::endl;
144  }
145  }
146  }
147 
149  // A Side South Label Conventions //
151  // 1 - 12_R1 16 - 17_R1 31 - 22_R1
152  // 2 - 12_R2 17 - 17_R2 32 - 22_R2
153  // 3 - 12_R3 18 - 17_R3 33 - 22_R3
154  // 4 - 13_R1 19 - 18_R1 34 - 23_R1
155  // 5 - 13_R2 20 - 18_R2 35 - 23_R2
156  // 6 - 13_R3 21 - 18_R3 36 - 23_R3
157  // 7 - 14_R1 22 - 19_R1
158  // 8 - 14_R2 23 - 19_R2
159  // 9 - 14_R3 24 - 19_R3
160  // 10 - 15_R1 25 - 20_R1
161  // 11 - 15_R2 26 - 20_R2
162  // 12 - 15_R3 27 - 20_R3
163  // 13 - 16_R1 28 - 21_R1
164  // 14 - 16_R2 29 - 21_R2
165  // 15 - 16_R3 30 - 21_R3
166 
168  // C Side North Label Conventions //
170  // 1 - 00_R1 16 - 05_R1 31 - 10_R1
171  // 2 - 00_R2 17 - 05_R2 32 - 10_R2
172  // 3 - 00_R3 18 - 05_R3 33 - 10_R3
173  // 4 - 01_R1 19 - 06_R1 34 - 11_R1
174  // 5 - 01_R2 20 - 06_R2 35 - 11_R2
175  // 6 - 01_R3 21 - 06_R3 36 - 11_R3
176  // 7 - 02_R1 22 - 07_R1
177  // 8 - 02_R2 23 - 07_R2
178  // 9 - 02_R3 24 - 07_R3
179  // 10 - 03_R1 25 - 08_R1
180  // 11 - 03_R2 26 - 08_R2
181  // 12 - 03_R3 27 - 08_R3
182  // 13 - 04_R1 28 - 09_R1
183  // 14 - 04_R2 29 - 09_R2
184  // 15 - 04_R3 30 - 09_R3
185 
186 
187  gStyle->SetOptStat(0);
188 
189  // if have lists inside input file - create pointer to list
190 
191  /* // Open up the histograms from the infile and give new names
192  THnSparseD* sparse = (THnSparseD *)infile1->Get("histsparse"); //get the sparse object and store it in memory
193  if(!sparse){
194  cout<<"sparse does not exist in "<<input1<<endl; //double check that it's in the file it should be !
195  return;
196  }
197  sparse->THnSparse::GetAxis(2)->SetRange(1,15); //limit to ONLY real trips
198  */
199  //____________________________________________________________________________________________________//
200 
201  // Heat map Visualization
202 
203  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)
204  const Int_t N_thBins = 12; //(12 theta bins of uniform angle (360/12 = 30 degrees = TMath::Pi()/6 ~= 0.523 rad) )
205 
206  Double_t rBin_edges[N_rBins + 1] = {0.0, 0.256, 0.504, 0.752, 1.00}; //variable edges for the radial dimensions
207 
208  TGraphPolargram* polardummy1 = new TGraphPolargram("polardummy1",0,1,0,2.*TMath::Pi()); //dummy plots to get the canvas right (not to be filled)
209  polardummy1->SetToGrad();
210  polardummy1->SetNdivPolar(N_thBins);
211  polardummy1->SetLineColor(0);
212  polardummy1->SetRadialLabelSize(0);
213 
214  TGraphPolargram* polardummy2 = new TGraphPolargram("polardummy2",0,1,0,2.*TMath::Pi());
215  polardummy2->SetToGrad();
216  polardummy2->SetNdivPolar( N_thBins);
217  polardummy2->SetLineColor(0);
218  polardummy2->SetRadialLabelSize(0);
219 
220  for(Int_t i = 0 ; i < 12 ; i++){ //setting the axis label (CCW from horizontal right axis)
221  char labelstr1[128];
222  char labelstr2[128];
223  if(i<=9){ // i -> [0:9]
224  sprintf(labelstr2,"C0%d",i);
225 
226  if(i<=6){ // i -> [0:6] (halfway)
227  sprintf(labelstr1,"A%d",18-i);
228  }
229  else if(i>6){ // i -> [7:9]
230  sprintf(labelstr1,"A%d",30-i);
231  }
232 
233  }
234  else { // i -> [10:11]
235  sprintf(labelstr2,"C%d",i);
236  sprintf(labelstr1,"A%d",30-i);
237  }
238  polardummy1->SetPolarLabel(i,labelstr1);
239  polardummy2->SetPolarLabel(i,labelstr2);
240  }
241 
242  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
243 
244  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
245 
246 
247  Double_t r, theta;
248  Int_t trip_count_total = 0;
249  Bool_t is_ASIDE = true;
250 
251  for (Int_t i = 0; i < 36; i++) {
252  if(sub_arrA[i] > 0.0){
253  Locate(i + 1, true, &r, &theta);
254  ErrASide->Fill(theta, r, sub_arrA[i]);
255  // cout<<"Region # A "<<(i)<<" Alive Fraction = "<<sub_arrA[i]<<endl;
256  }
257  }
258 
259  for (Int_t i = 0; i < 36; i++) {
260  if(sub_arrC[i] > 0.0){
261  Locate(i + 1,false, &r, &theta);
262  ErrCSide->Fill(theta, r, sub_arrC[i]);
263  // cout<<"Region # C "<<(i)<<" Alive Fraction = "<<sub_arrC[i]<<endl;
264  }
265  }
266 
267  sprintf(name, "%d-Noise Std. Dev. in ADC (North Side)",run_num);
268  TH2D* dummy_his1 = new TH2D("dummy1", name, 100, -1.5, 1.5, 100, -1.5, 1.5);
269  //TH2D* dummy_his1 = new TH2D("dummy1", "10616-Noise map in ADC unit North Side", 100, -1.5, 1.5, 100, -1.5, 1.5); //dummy histos for titles
270  sprintf(name, "%d-Noise Std. Dev. in ADC (South Side)",run_num);
271  TH2D* dummy_his2 = new TH2D("dummy2", name, 100, -1.5, 1.5, 100, -1.5, 1.5);
272  //TH2D* dummy_his2 = new TH2D("dummy2", "10616-Noise map in ADC unit South Side", 100, -1.5, 1.5, 100, -1.5, 1.5);
273  //TPaveLabels for sector labels
274  TPaveLabel* A00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "18" );
275  TPaveLabel* A01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "17" );
276  TPaveLabel* A02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "16" );
277  TPaveLabel* A03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "15" );
278  TPaveLabel* A04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "14" );
279  TPaveLabel* A05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "13" );
280  TPaveLabel* A06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "12" );
281  TPaveLabel* A07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "23" );
282  TPaveLabel* A08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "22" );
283  TPaveLabel* A09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "21" );
284  TPaveLabel* A10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "20" );
285  TPaveLabel* A11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "19" );
286 
287  TPaveLabel* C00 = new TPaveLabel( 1.046586,-0.1938999,1.407997,0.2144871, "00" );
288  TPaveLabel* C01 = new TPaveLabel( 0.962076,0.4382608,1.323487,0.8466479 , "01" );
289  TPaveLabel* C02 = new TPaveLabel( 0.4801947,0.8802139,0.8416056,1.288601 , "02" );
290  TPaveLabel* C03 = new TPaveLabel( -0.1823921,1.011681,0.1790189,1.425662, "03" );
291  TPaveLabel* C04 = new TPaveLabel( -0.8449788,0.8690253,-0.4835679,1.288601 , "04" );
292  TPaveLabel* C05 = new TPaveLabel( -1.30879,0.441058,-0.9473786,0.8550394 , "05" );
293  TPaveLabel* C06 = new TPaveLabel( -1.411009,-0.2050886,-1.049598,0.2144871, "06" );
294  TPaveLabel* C07 = new TPaveLabel( -1.302585,-0.7757116,-0.9471979,-0.3561359 , "07" );
295  TPaveLabel* C08 = new TPaveLabel( -0.8449788,-1.309971,-0.4835679,-0.8848013 , "08" );
296  TPaveLabel* C09 = new TPaveLabel( -0.1823921,-1.426557,0.1790189,-1.006982 , "09" );
297  TPaveLabel* C10 = new TPaveLabel( 0.4801947,-1.309076,0.8416056,-0.8839062 , "10" );
298  TPaveLabel* C11 = new TPaveLabel( 0.9622567,-0.7785088,1.323668,-0.3533387 , "11" );
299 
300  A00->SetFillColor(0);
301  A01->SetFillColor(0);
302  A02->SetFillColor(0);
303  A03->SetFillColor(0);
304  A04->SetFillColor(0);
305  A05->SetFillColor(0);
306  A06->SetFillColor(0);
307  A07->SetFillColor(0);
308  A08->SetFillColor(0);
309  A09->SetFillColor(0);
310  A10->SetFillColor(0);
311  A11->SetFillColor(0);
312 
313  C00->SetFillColor(0);
314  C01->SetFillColor(0);
315  C02->SetFillColor(0);
316  C03->SetFillColor(0);
317  C04->SetFillColor(0);
318  C05->SetFillColor(0);
319  C06->SetFillColor(0);
320  C07->SetFillColor(0);
321  C08->SetFillColor(0);
322  C09->SetFillColor(0);
323  C10->SetFillColor(0);
324  C11->SetFillColor(0);
325 
326  gStyle->SetPalette(kBird);
327 
328  TCanvas *Error_Viz = new TCanvas("Error_Viz", "Error_Viz", 1248, 598);
329  Error_Viz->Divide(2,1);
330  Error_Viz->cd(1);
331  dummy_his1->Draw("");
332  //polardummy1->Draw("same");
333  ErrCSide->Draw("colpolzsame0");
334  C00->Draw("same");
335  C01->Draw("same");
336  C02->Draw("same");
337  C03->Draw("same");
338  C04->Draw("same");
339  C05->Draw("same");
340  C06->Draw("same");
341  C07->Draw("same");
342  C08->Draw("same");
343  C09->Draw("same");
344  C10->Draw("same");
345  C11->Draw("same");
346  Error_Viz->cd(2);
347  dummy_his2->Draw("");
348  //polardummy2->Draw("");
349  ErrASide->Draw("colpolzsame0");
350  A00->Draw("same");
351  A01->Draw("same");
352  A02->Draw("same");
353  A03->Draw("same");
354  A04->Draw("same");
355  A05->Draw("same");
356  A06->Draw("same");
357  A07->Draw("same");
358  A08->Draw("same");
359  A09->Draw("same");
360  A10->Draw("same");
361  A11->Draw("same");
362 
363  ErrCSide->SetMaximum(4);
364  ErrASide->SetMaximum(4);
365 
366  ErrCSide->SetMinimum(0);
367  ErrASide->SetMinimum(0);
368 
369 
370  //Set Same Scale for A and C side displays
371 
372  Double_t Maxval = TMath::Max(ErrASide->GetBinContent(ErrASide->GetMaximumBin()),ErrCSide->GetBinContent(ErrCSide->GetMaximumBin()));
373  //ErrASide->SetMaximum(Maxval);
374  //ErrCSide->SetMaximum(Maxval);
375 
376  //ErrASide->SetMinimum(0);
377  //ErrCSide->SetMinimum(0);
378 
379  // std::cout <<"A Side Entries=" << ErrASide->GetBinContent() << std::endl;
380  // std::cout <<"C Side Entries=" << ErrCSide->GetBinContent() << std::endl;
381  // std::cout <<" Min A Side Bin" << ErrASide->GetMinimumBin() << std::endl;
382  // std::cout <<" Min C Side Bin" << ErrCSide->GetMinimumBin() << std::endl;
383 
384  //____________________________________________________________________________________________________//
385 
386 
387  TFile *outf = new TFile("Trip_Histos.root","RECREATE");
388  Error_Viz->Write();
389  //Trip_per_stack_dist->Scale(1./trip_count_total);
390  //Trip_per_stack_dist->Fit("PoissonFit","ERS");
391  //cout<< "\nThe total trips = " << trip_count_total << endl;
392  //delete sparsforVIZ;
393 
394  sprintf(name,"/sphenix/user/llegnosky/run_%d/run_%d_Noise.png",run_num,run_num);
395  Error_Viz->Print(name);
396 
397  outf->Write();
398 }
399 
400 void Locate(Int_t id, Bool_t is_ASIDE, Double_t *rbin, Double_t *thbin) {
401 
402  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)
403  /*
404  Double_t CSIDE_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 }; // CW from x = 0 (RHS horizontal)*/
405 
406  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)
407  /*
408  Double_t ASIDE_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)
409  */
410  Int_t modid3 = id % 3;
411 
412  switch(modid3) {
413  case 1:
414  *rbin = 0.4; //R1
415  break;
416  case 2:
417  *rbin = 0.6; //R2
418  break;
419  case 0:
420  *rbin = 0.8; //R3
421  break;
422  }
423 
424 
425  if( is_ASIDE ){
426  *thbin = CSIDE_angle_bins[TMath::FloorNint((id-1)/3)];
427  }
428  else{
429  *thbin = ASIDE_angle_bins[TMath::FloorNint((id-1)/3)];
430  }
431 
432 }