AliPhysics  1811c8f (1811c8f)
MultiplyPi0CalibrationFactors_TextToHisto_Final.C
Go to the documentation of this file.
1 #include <TChain.h>
2 #include <TNtuple.h>
3 #include <TObjArray.h>
4 #include <TSystem.h>
5 #include <TString.h>
6 #include <TH1F.h>
7 #include <TVector.h>
8 #include <TRefArray.h>
9 #include <TArrayS.h>
10 #include "TError.h"
11 #include "TTree.h"
12 #include "TClonesArray.h"
13 #include "TGraphErrors.h"
14 #include "TPostScript.h"
15 #include "TLegend.h"
16 #include "TH2I.h"
17 #include "TF1.h"
18 #include "TStyle.h"
19 #include "TCanvas.h"
20 #include "TPolyLine.h"
21 #include "TLine.h"
22 #include "TFile.h"
23 #include "TMath.h"
24 #include "TLeaf.h"
25 #include "TBranch.h"
26 
27 
28 //#include "/cebaf/cebaf/EMCAL/cosmicsAnalysis/macros/defineMyPalette40All.C"
29 
30 
42 
47 {int a,b,c,d,aTot,bTot,cTot,dTot,i,j,jFile,kNbFiles,iSM;
48  float e,eTot;
49  double minHist,maxHist,minHistProduct;
50 
51 
52  char SMP2Name[][100]={"SMA0","SMC0","SMA1","SMC1","SMA2","SMC2","SMA3","SMC3","SMA4","SMC4","SMA5","SMC5","SMA9","SMC9","SMA10","SMC10","SMA11","SMC11","SMA12","SMC12"};
53  char SMnumber[][100]={"0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19"};
54 
56  int detTypeType[]={kEMCAL,kEMCALthird,kDCAL,kDCALthird};
57  char detTypeString[][100]={"EMCAL","EMCALthird","DCAL","DCALthird"};
59  const int kNbColEMCAL=48;
60  const int kNbRowEMCAL=24;
61  const int kNbSMEMCAL=10;
62  const int kNbColEMCALthird=kNbColEMCAL;
63  const int kNbRowEMCALthird=(int)(kNbRowEMCAL/3);
64  const int kNbSMEMCALthird=2;
65  const int kNbColDCAL=32;
66  const int kNbRowDCAL=kNbRowEMCAL;
67  const int kNbSMDCAL=6;
70  const int kNbSMDCALthird=2;
71  const int kNbSMtot=kNbSMEMCAL+kNbSMEMCALthird+kNbSMDCAL+kNbSMDCALthird;
72  const int kTabNbCol[4]={kNbColEMCAL,kNbColEMCALthird,kNbColDCAL,kNbColDCALthird};
73  const int kTabNbRow[4]={kNbRowEMCAL,kNbRowEMCALthird,kNbRowDCAL,kNbRowDCALthird};
74  const int kTabNbSM[4]={kNbSMEMCAL,kNbSMEMCALthird,kNbSMDCAL,kNbSMDCALthird};
75  const int kNbColMax=kNbColEMCAL;
76  const int kNbRowMax=kNbRowEMCAL;
77  const int kNbColOffsetDCAL=kNbColEMCAL-kNbColDCAL;
78 
79 
80  //CUSTOMIZE customize :
81  kNbFiles=3;
82 
83 
84  char psfile[150];
85  FILE *txtFileOut = NULL;
86  txtFileOut = fopen("multiplyPi0CalibrationFactors_TextToHisto_Final.txt","w");
87  TFile * f = new TFile("multiplyPi0CalibrationFactors_TextToHisto_Final.root","recreate");
88  sprintf(psfile,"multiplyPi0CalibrationFactors_TextToHisto_Final.ps");
89 
90  //defineMyPalette40All(30,5);
91 
92  //Coeff files to be read and multiplied :
93  FILE **calibCoeffsInput;
94  calibCoeffsInput = new FILE*[kNbFiles];
95  //CUSTOMIZE customize :
96  calibCoeffsInput[0] = fopen("/cebaf/cebaf/EMCAL/calibPi0_run2/calibPi0_4_with2015data/output/pass0_mergedDCALandThirdSMs/output_calibPi0_coeffs_clean_veryHighTowers.txt","r"); //Chosen for pass0.
97  calibCoeffsInput[1] = fopen("/cebaf/cebaf/EMCAL/calibPi0_run2/calibPi0_4_with2015data/output/pass1_DCALandThirdSMs/output_calibPi0_coeffs_clean_veryHighTowers.txt","r"); //Chosen for pass1.
98  calibCoeffsInput[2] = fopen("/cebaf/cebaf/EMCAL/calibPi0_run2/calibPi0_4_with2015data/output/pass2_DCALandThirdSMsVeryHighTowers/output_calibPi0_coeffs_clean_finalFile4HVcalculation_setUntrustedToOne.txt","r"); //Chosen for pass2.
99  //Histoes definition :
100  TH2F **h2Coeffs;
101  h2Coeffs = new TH2F*[kNbSMtot];
102  for (i=0;i<kNbSMtot;i++)
103  {h2Coeffs[i] = new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),kNbColMax,0,kNbColMax,kNbRowMax,0,kNbRowMax); //All SMs use the same size.
104  }
105 
106  TH2F **hSpace;
107  TH1F **hDistr; //Coeff for all SMs, per step
108  TH1F **hDistrPerSM; //Coeff for all steps, per SM
109  TH1F **hDistrNoone; //Coeff for all SMs, per step, exclude towers at 1.
110  TH1F **hDistrNoonePerSM; //Coeff for all steps, per SM, exclude towers at 1.
111  hSpace = new TH2F*[(kNbFiles+1)*kNbSMtot];
112  hDistr = new TH1F*[(kNbFiles+1)*2];
113  hDistrPerSM = new TH1F*[kNbSMtot];
114  hDistrNoone = new TH1F*[(kNbFiles+1)*2];
115  hDistrNoonePerSM = new TH1F*[kNbSMtot];
116  for (i=0;i<kNbSMtot;i++)
117  {hDistrPerSM[i] = new TH1F(Form("hDistrPerSM%d",i),Form("hDistrPerSM%d",i),100,0.6,1.4);
118  hDistrPerSM[i]->SetXTitle("Coefficient");
119  hDistrPerSM[i]->SetYTitle("Counts");
120  hDistrPerSM[i]->SetStats(0);
121  hDistrNoonePerSM[i] = new TH1F(Form("hDistrNoonePerSM%d",i),Form("hDistrNoonePerSM%d",i),100,0.6,1.4);
122  hDistrNoonePerSM[i]->SetXTitle("Coefficient");
123  hDistrNoonePerSM[i]->SetYTitle("Counts");
124  hDistrNoonePerSM[i]->SetStats(0);
125  for (j=0;j<kNbFiles+1;j++)
126  {hSpace[(kNbFiles+1)*i+j] = new TH2F(Form("hSpace_%d_%d",i,j),Form("hSpace_%d_%d",i,j),kNbColMax,-0.5,kNbColMax-0.5,kNbRowMax,-0.5,kNbRowMax-0.5);
127  hSpace[(kNbFiles+1)*i+j]->SetXTitle("Column");
128  hSpace[(kNbFiles+1)*i+j]->SetYTitle("Row");
129  hSpace[(kNbFiles+1)*i+j]->SetStats(0);
130  hSpace[(kNbFiles+1)*i+j]->SetContour(30);
131  }
132  }
133  for (j=0;j<kNbFiles+1;j++)
134  {hDistr[j] = new TH1F(Form("hDistr_%d",j),Form("hDistr_%d",j),100,0.0,3.0);
135  hDistr[j]->SetMaximum(20.);
136  hDistr[(kNbFiles+1)+j] = new TH1F(Form("hDistrZm_%d",j),Form("hDistrZm_%d",j),100,0.85,1.15);
137  hDistrNoone[j] = new TH1F(Form("hDistrNoone_%d",j),Form("hDistrNoone_%d",j),100,0.0,3.0);
138  hDistrNoone[j]->SetMaximum(20.);
139  hDistrNoone[(kNbFiles+1)+j] = new TH1F(Form("hDistrNooneZm_%d",j),Form("hDistrNooneZm_%d",j),100,0.85,1.15);
140  for (i=0;i<2;i++)
141  {hDistr[(kNbFiles+1)*i+j]->SetXTitle("Coefficient");
142  hDistr[(kNbFiles+1)*i+j]->SetYTitle("Counts");
143  hDistr[(kNbFiles+1)*i+j]->SetStats(0);
144  hDistr[(kNbFiles+1)*i+j]->SetLineColor(4);
145  hDistrNoone[(kNbFiles+1)*i+j]->SetXTitle("Coefficient");
146  hDistrNoone[(kNbFiles+1)*i+j]->SetYTitle("Counts");
147  hDistrNoone[(kNbFiles+1)*i+j]->SetStats(0);
148  hDistrNoone[(kNbFiles+1)*i+j]->SetLineColor(4);
149  }
150  }
151 
152  //Calculation of coeffs :
153  for (iSM=0;iSM<kNbSMtot;iSM++)
154  {for(i = 0;i<kTabNbCol[SMdetType[iSM]]*kTabNbRow[SMdetType[iSM]];i++)
155  {jFile=0;
156  fscanf(calibCoeffsInput[jFile], " %d %d %d %d %f\n", &a, &b, &c,&d,&e);
157  //printf( " %d %d %d %d %f\n", a, b, c,d,e);
158  if (b != iSM) printf("File reading out of sync : received SM %d while expected %d.\n",b,iSM);
159  hSpace[(kNbFiles+1)*b+jFile]->Fill(c,d,e);
160  hDistr[jFile]->Fill(e);
161  hDistr[(kNbFiles+1)+jFile]->Fill(e);
162  if (e != 1.)
163  {hDistrNoone[jFile]->Fill(e);
164  hDistrNoone[(kNbFiles+1)+jFile]->Fill(e);
165  }
166  aTot=a;
167  bTot=b;
168  cTot=c;
169  dTot=d;
170  eTot=e;
171  for (jFile=1;jFile<kNbFiles;jFile++)
172  {fscanf(calibCoeffsInput[jFile], " %d %d %d %d %f\n", &a, &b, &c,&d,&e);
173  hSpace[(kNbFiles+1)*b+jFile]->Fill(c,d,e);
174  hDistr[jFile]->Fill(e);
175  hDistr[(kNbFiles+1)+jFile]->Fill(e);
176  if (e != 1.)
177  {hDistrNoone[jFile]->Fill(e);
178  hDistrNoone[(kNbFiles+1)+jFile]->Fill(e);
179  }
180  if (bTot != b) printf("Files reading out of sync : (%d,%d,%d) vs (%d,%d,%d)\n",b,c,d,bTot,cTot,dTot);
181  if (cTot != c) printf("Files reading out of sync : (%d,%d,%d) vs (%d,%d,%d)\n",b,c,d,bTot,cTot,dTot);
182  if (dTot != d) printf("Files reading out of sync : (%d,%d,%d) vs (%d,%d,%d)\n",b,c,d,bTot,cTot,dTot);
183  eTot*=e;
184  if (jFile == (kNbFiles-1))
185  {if (e == 1) eTot=1.0; //If tower not trusted at last iteration, don't trust it for the whole process.
186  }
187  }
188  hSpace[(kNbFiles+1)*bTot+kNbFiles]->Fill(cTot,dTot,eTot);
189  hDistr[kNbFiles]->Fill(eTot);
190  hDistr[(kNbFiles+1)+kNbFiles]->Fill(eTot);
191  hDistrPerSM[bTot]->Fill(eTot);
192  if (eTot != 1.)
193  {hDistrNoone[kNbFiles]->Fill(eTot);
194  hDistrNoone[(kNbFiles+1)+kNbFiles]->Fill(eTot);
195  hDistrNoonePerSM[bTot]->Fill(eTot);
196  }
197  fprintf(txtFileOut,"%d %d %d %f\n",bTot,cTot,dTot,eTot);
198  if(eTot==0) printf("\n###### Coeff cannot be 0 ! Please check tower : %d %d %d %f #######\n\n",bTot,cTot,dTot,eTot);
199  else h2Coeffs[bTot]->SetBinContent(cTot,dTot,1./eTot);
200  }
201  }
202 
203 
204  // Draw plots :
205 
206  const int cWidth=500;
207  const int cHeight=(int)(500*(29./21.));
208  TCanvas *c1 = new TCanvas("c1","EMCal cosmics analysis",cWidth,cHeight);
209  TPostScript *ps = new TPostScript(psfile,111);
210 
211  //Spatial distrib of coeffs, per SM, for each pass and integrated over all passes.
212  for (i=0;i<kNbSMtot;i++)
213  {minHist=1.;
214  maxHist=1.;
215  minHistProduct=1.;
216  for (j=0;j<kNbFiles;j++)
217  {//if (hSpace[(kNbFiles+1)*i+j]->GetMinimum() < minHist) minHist=hSpace[(kNbFiles+1)*i+j]->GetMinimum();
218  if (hSpace[(kNbFiles+1)*i+j]->GetMaximum() > maxHist) maxHist=hSpace[(kNbFiles+1)*i+j]->GetMaximum();
219  for (int iCol=0;iCol<kTabNbCol[SMdetType[i]];iCol++)
220  {for (int iRow=0;iRow<kTabNbRow[SMdetType[i]];iRow++)
221  {if (hSpace[(kNbFiles+1)*i+j]->GetBinContent(iCol+1,iRow+1) < minHist) minHist=hSpace[(kNbFiles+1)*i+j]->GetBinContent(iCol+1,iRow+1);
222  }
223  }
224  }
225  for (int iCol=0;iCol<kTabNbCol[SMdetType[i]];iCol++)
226  {for (int iRow=0;iRow<kTabNbRow[SMdetType[i]];iRow++)
227  {if (hSpace[(kNbFiles+1)*i+kNbFiles]->GetBinContent(iCol+1,iRow+1) < minHistProduct) minHistProduct=hSpace[(kNbFiles+1)*i+kNbFiles]->GetBinContent(iCol+1,iRow+1);
228  }
229  }
230  ps->NewPage();
231  c1->Clear();
232  c1->Divide(2,5);
233  for (j=0;j<kNbFiles;j++)
234  {c1->cd(j+1);
235  hSpace[(kNbFiles+1)*i+j]->SetMinimum(minHist);
236  hSpace[(kNbFiles+1)*i+j]->SetMaximum(maxHist);
237  hSpace[(kNbFiles+1)*i+j]->Draw("COLZ");
238  hSpace[(kNbFiles+1)*i+j]->Write();
239  }
240  c1->cd(kNbFiles+1);
241  hSpace[(kNbFiles+1)*i+kNbFiles]->SetMinimum(minHistProduct);
242  hSpace[(kNbFiles+1)*i+kNbFiles]->Draw("COLZ");
243  hSpace[(kNbFiles+1)*i+kNbFiles]->Write();
244  c1->Update();
245  }
246 
247  //1-D distribs of all calib coeffs, per pass and all passes integrated.
248  //Page "1" : zoomed (in y) ; page "0" : not zoomed.
249  for (i=0;i<2;i++)
250  {ps->NewPage();
251  c1->Clear();
252  c1->Divide(2,5);
253  for (j=0;j<kNbFiles;j++)
254  {c1->cd(j+1);
255  /*printf("Histo coeffs pass %d ",j+1);
256  switch (i)
257  {case 0 : printf("(large range) ");
258  break;
259  case 1 : printf("(reduced range) ");
260  break;
261  default : printf("##### unknown histo i=%d ##### \n",i);
262  }
263  printf(": mean = %f, RMS = %f\n",hDistr[(kNbFiles+1)*i+j]->GetMean(),hDistr[(kNbFiles+1)*i+j]->GetRMS());*/
264  hDistr[(kNbFiles+1)*i+j]->Draw();
265  hDistr[(kNbFiles+1)*i+j]->Write();
266  }
267  c1->cd(kNbFiles+1);
268  /*printf("Histo 'product of all coeffs' ");
269  switch (i)
270  {case 0 : printf("(large range) ");
271  break;
272  case 1 : printf("(reduced range) ");
273  break;
274  default : printf("##### unknown histo i=%d ##### \n",i);
275  }
276  printf(": mean = %f, RMS = %f\n",hDistr[(kNbFiles+1)*i+kNbFiles]->GetMean(),hDistr[(kNbFiles+1)*i+kNbFiles]->GetRMS());*/
277  /*if (i == 1)
278  {TF1 *func = new TF1("func","gaus(0)",0.85,1.15);
279  func->SetParLimits(0,10.,500.);
280  func->SetParLimits(1,0.9,1.1);
281  func->SetParLimits(2,0.01,0.2);
282  func->SetParameter(0,hDistr[(kNbFiles+1)+kNbFiles]->GetMaximum());
283  func->SetParameter(1,1.);
284  func->SetParameter(2,0.1);
285  printf("\n");
286  hDistr[(kNbFiles+1)+kNbFiles]->Fit(func,"BRM"); //Option I crashe
287  //printf("Histo moy = %f, RMS = %f\n",hDistr[(kNbFiles+1)+kNbFiles]->GetMean(),hDistr[(kNbFiles+1)+kNbFiles]->GetRMS());
288  printf("\nFit over 'product of all coeffs' (reduced range) : mean = %f, RMS = %f\n",func->GetParameter(1),func->GetParameter(2));
289  }*/
290  hDistr[(kNbFiles+1)*i+kNbFiles]->Draw();
291  hDistr[(kNbFiles+1)*i+kNbFiles]->Write();
292  c1->Update();
293  }
294 
295  //Same thing now without the towers with coeff at 1 :
296  for (i=0;i<2;i++)
297  {ps->NewPage();
298  c1->Clear();
299  c1->Divide(2,5);
300  for (j=0;j<kNbFiles;j++)
301  {c1->cd(j+1);
302  printf("Histo coeffs pass %d ",j+1);
303  switch (i)
304  {case 0 : printf("(large range) ");
305  break;
306  case 1 : printf("(reduced range) ");
307  break;
308  default : printf("##### unknown histo i=%d ##### \n",i);
309  }
310  printf(": mean = %f, RMS = %f\n",hDistrNoone[(kNbFiles+1)*i+j]->GetMean(),hDistrNoone[(kNbFiles+1)*i+j]->GetRMS());
311  hDistrNoone[(kNbFiles+1)*i+j]->Draw();
312  hDistrNoone[(kNbFiles+1)*i+j]->Write();
313  }
314  c1->cd(kNbFiles+1);
315  printf("Histo 'product of all coeffs' ");
316  switch (i)
317  {case 0 : printf("(large range) ");
318  break;
319  case 1 : printf("(reduced range) ");
320  break;
321  default : printf("##### unknown histo i=%d ##### \n",i);
322  }
323  printf(": mean = %f, RMS = %f\n",hDistrNoone[(kNbFiles+1)*i+kNbFiles]->GetMean(),hDistrNoone[(kNbFiles+1)*i+kNbFiles]->GetRMS());
324  if (i == 1)
325  {TF1 *func = new TF1("func","gaus(0)",0.85,1.15);
326  func->SetParLimits(0,10.,500.);
327  func->SetParLimits(1,0.9,1.1);
328  func->SetParLimits(2,0.01,0.2);
329  func->SetParameter(0,hDistrNoone[(kNbFiles+1)+kNbFiles]->GetMaximum());
330  func->SetParameter(1,1.);
331  func->SetParameter(2,0.1);
332  printf("\n");
333  hDistrNoone[(kNbFiles+1)+kNbFiles]->Fit(func,"BRM"); //Option I crashe
334  //printf("Histo moy = %f, RMS = %f\n",hDistrNoone[(kNbFiles+1)+kNbFiles]->GetMean(),hDistrNoone[(kNbFiles+1)+kNbFiles]->GetRMS());
335  printf("\nFit over 'product of all coeffs' (reduced range) : mean = %f, RMS = %f\n",func->GetParameter(1),func->GetParameter(2));
336  }
337  hDistrNoone[(kNbFiles+1)*i+kNbFiles]->Draw();
338  hDistrNoone[(kNbFiles+1)*i+kNbFiles]->Write();
339  c1->Update();
340  }
341 
342  //Coeff (all iterations multiplied) distrib SM per SM
343  printf("\n\nHistoes 'product of all coeffs' for each SM :\n");
344 
345  ps->NewPage();
346  c1->Clear();
347  c1->Divide(2,5);
348  for (j=0;j<kNbSMEMCAL;j++)
349  {c1->cd(j+1);
350  //printf(" SM %d : mean = %f, RMS = %f\n",j,hDistrPerSM[j]->GetMean(),hDistrPerSM[j]->GetRMS());
351  hDistrPerSM[j]->Draw();
352  hDistrPerSM[j]->Write();
353  }
354  c1->Update();
355 
356  ps->NewPage();
357  c1->Clear();
358  c1->Divide(2,5);
359  for (j=kNbSMEMCAL;j<kNbSMtot;j++)
360  {c1->cd(j-kNbSMEMCAL+1);
361  //printf(" SM %d : mean = %f, RMS = %f\n",j,hDistrPerSM[j]->GetMean(),hDistrPerSM[j]->GetRMS());
362  hDistrPerSM[j]->Draw();
363  hDistrPerSM[j]->Write();
364  }
365  c1->Update();
366 
367  //Same thing excluding the towers which have coeff equal to 1 :
368  ps->NewPage();
369  c1->Clear();
370  c1->Divide(2,5);
371  for (j=0;j<kNbSMEMCAL;j++)
372  {c1->cd(j+1);
373  printf(" SM %d : mean = %f, RMS = %f\n",j,hDistrNoonePerSM[j]->GetMean(),hDistrNoonePerSM[j]->GetRMS());
374  hDistrNoonePerSM[j]->Draw();
375  hDistrNoonePerSM[j]->Write();
376  }
377  c1->Update();
378 
379  ps->NewPage();
380  c1->Clear();
381  c1->Divide(2,5);
382  for (j=kNbSMEMCAL;j<kNbSMtot;j++)
383  {c1->cd(j-kNbSMEMCAL+1);
384  printf(" SM %d : mean = %f, RMS = %f\n",j,hDistrNoonePerSM[j]->GetMean(),hDistrNoonePerSM[j]->GetRMS());
385  hDistrNoonePerSM[j]->Draw();
386  hDistrNoonePerSM[j]->Write();
387  }
388  c1->Update();
389 
390 
391  ps->Close();
392 
393  for (i=0;i<kNbSMtot;i++) h2Coeffs[i]->Write();
394 
395  f->Close();
396 
397  fclose(txtFileOut);
398 
399  for (jFile=kNbFiles-1;jFile>-1;jFile--) fclose(calibCoeffsInput[jFile]);
400 
401  return;
402 
403  }
404 
405 
406 
407 
void Draw(const char *filename, const char *title="", const char *others="ALL", const char *options="DEFAULT", const char *outFlg="ALL", UShort_t rebin=5, Float_t eff=0, const char *base="")
Definition: DrawdNdeta.C:3603
const int kNbColEMCALthird
Definition: External.C:236
const int kTabNbCol[4]
char detTypeString[][100]
const int kNbRowMax
TCanvas * c
Definition: TestFitELoss.C:172
const int kNbSMtot
const int kNbSMEMCAL
const int kNbRowEMCAL
const int kNbColOffsetDCAL
int SMdetType[]
const int kNbRowDCAL
const int kNbRowDCALthird
char SMnumber[][100]
const int kNbColMax
const int kNbSMEMCALthird
const int kNbColDCAL
char SMP2Name[][100]
const int kTabNbRow[4]
const int kNbColEMCAL
int detTypeType[]
const int kNbSMDCALthird
const int kNbRowEMCALthird
const int kNbSMDCAL
const int kTabNbSM[4]
const int kNbColDCALthird
void MultiplyPi0CalibrationFactors_TextToHisto_Final()