AliPhysics  239578e (239578e)
CalculateParamChi2MCxTalkDataPerSM.C
Go to the documentation of this file.
1 
14 #if !defined(__CINT__) || defined(__MAKECINT__)
15 
16 #include <TFile.h>
17 #include <TDirectoryFile.h>
18 #include <TList.h>
19 #include <TString.h>
20 #include <TROOT.h>
21 #include <TStyle.h>
22 #include <TH1D.h>
23 #include <TH2F.h>
24 #include <TH3F.h>
25 #include <TCanvas.h>
26 #include <TPad.h>
27 #include <TLegend.h>
28 #include <TObject.h>
29 #include <TAxis.h>
30 #include <TGaxis.h>
31 #include <TLine.h>
32 #include <TF1.h>
33 #include <TMath.h>
34 #include <TArrayD.h>
35 #include <TGraphErrors.h>
36 #endif
37 
38 void Chi2( TH1D* h, TH1D* hD,TH1D*& hChi2, TH1D*& hDiff, Bool_t debug = kFALSE );
39 
40 TString fileFormat = ".eps";
41 
42 Int_t firstEbin = 0;
43 
44 //------------------------------------------------------------------------------
57 //------------------------------------------------------------------------------
58 void CalculateParamChi2MCxTalkDataPerSM
59 (
60  TString histoName = "SMM02NoCut",
61  TString inputFileName = "",
62  Int_t dataProd = 0,
63  Int_t firstSM = 0,
64  Int_t lastSM = 9,
65  TArrayD binE = 0,
66  TArrayD prodParam = 0,
67  Int_t rebin = 1,
68  Bool_t debug = kFALSE
69  )
70 {
71  Int_t color[] = {1,4,2,kYellow-2,8,kCyan,kOrange+2,kViolet,kOrange-2,4,2,6,8,9,kYellow-6,kCyan,kOrange+2,kViolet,kOrange-2};
72 
73  // Define Chi2 minimum finding window
74  Float_t xmin = 0.1;
75  Float_t xmax = 1.0;
76  if ( histoName.Contains("SMM02") )
77  {
78  xmin = 0.1;
79  xmax = 1.0;
80  }
81  if ( histoName == "SMNCell" && !histoName.Contains("Module") )
82  {
83  xmin = 2 ;
84  xmax = 12 ;
85  }
86  if ( histoName.Contains("SMM20Low") )
87  {
88  xmin = 0.05; xmax = 0.25;
89  }
90  if ( histoName.Contains("SMM20Hig") )
91  {
92  xmin = 0.05; xmax = 0.30;
93  }
94  if ( histoName.Contains("Module") )
95  {
96  xmin = 0.0; xmax = 100000;
97  }
98 
99  // Open the file with projected distributions
100  TFile * file = TFile::Open(Form("figures/%s/Projections_%s.root",histoName.Data(),inputFileName.Data()));;
101  if(!file) return;
102 
103  // Energy bins
104  const Int_t nEBins = binE.GetSize()-1;
105  Double_t binEErr[nEBins-1];
106  for(Int_t ie = 0; ie < nEBins; ie++)
107  binEErr[ie] = (binE.At(ie+1)-binE.At(ie))/2;
108 
109  // MC production parameter
110  const Int_t nProd = prodParam.GetSize();
111  Double_t prodParamE[nProd];
112  for(Int_t iprod = 0; iprod < nProd-1; iprod++)
113  {
114  prodParamE[iprod] = (prodParam.At(iprod+1)-prodParam.At(iprod))/2;
115  if ( debug ) printf("mc iprod %d: param %2.5f, next %2.5f, err %2.5f\n",
116  iprod,prodParam.At(iprod),prodParam.At(iprod+1),prodParamE[iprod]);
117  }
118  prodParamE[nProd] = prodParamE[nProd-1];
119 
120  // Total number of SM inspected
121  const Int_t nSM = lastSM-firstSM+1;
122  if ( debug )
123  printf("N MC prod %d N e bins %d, nSM %d\n", nProd, nEBins, nSM);
124 
125  //
126  // Recover the histograms from file
127  //
128  TH1D* hData [nEBins][nSM];
129  TH1D* hSimu[nProd][nEBins][nSM];
130 
131  for(Int_t iebin = 0; iebin < nEBins; iebin++)
132  {
133  for(Int_t ism = firstSM; ism <= lastSM; ism++)
134  {
135  hData[iebin][ism-firstSM] = (TH1D *) file->Get(Form("Prod%d_Ebin%d_Param%d",dataProd,iebin,ism));
136  if ( rebin > 1 )
137  hData[iebin][ism-firstSM]->Rebin(rebin);
138  if ( debug )
139  printf("data ie %d, ism %d, %p\n",iebin,ism,hData[iebin][ism-firstSM]);
140 
141  hSimu[dataProd][iebin][ism-firstSM] = 0;
142  for(Int_t iprod = 0; iprod < nProd; iprod++)
143  {
144  hSimu[iprod][iebin][ism-firstSM] = (TH1D *) file->Get(Form("Prod%d_Ebin%d_Param%d",iprod+dataProd+1,iebin,ism));
145  if ( rebin > 1 )
146  hSimu[iprod][iebin][ism-firstSM]->Rebin(rebin);
147 
148  if ( debug )
149  printf("mc %d ie %d, ism %d, %p\n",iprod,iebin,ism,hSimu[iprod][iebin][ism-firstSM]);
150 
151  } // prod
152  } // sm
153  } // ebin
154 
156  // Make Chi2 and difference distributions
158 
159  TH1D* hChi2 [nProd][nEBins][nSM];
160  TH1D* hDiff [nProd][nEBins][nSM];
161 
162  if ( debug )
163  printf("Make Chi2 distributions");
164 
165  for(Int_t iebin = 0; iebin < nEBins; iebin++)
166  {
167  for(Int_t iprod = 0; iprod < nProd; iprod++)
168  {
169  for(Int_t ism = firstSM; ism <= lastSM; ism++)
170  {
171  hChi2[iprod][iebin][ism-firstSM] =0;
172  hDiff[iprod][iebin][ism-firstSM] = 0;
173 
174  if ( debug )
175  printf("prod %d, ebin %d, sm %d\n",iprod,iebin,ism);
176 
177  Chi2(hSimu[iprod][iebin][ism-firstSM],hData[iebin][ism-firstSM],
178  hChi2[iprod][iebin][ism-firstSM],hDiff[iprod][iebin][ism-firstSM], debug);
179 
180  if ( debug )
181  printf("\t %p %p\n",hChi2[iprod][iebin][ism-firstSM],hDiff[iprod][iebin][ism-firstSM]);
182 
183  if(hChi2[iprod][iebin][ism-firstSM])hChi2[iprod][iebin][ism-firstSM]->SetTitle(Form("SM = %d",ism));
184  if(hDiff[iprod][iebin][ism-firstSM])hDiff[iprod][iebin][ism-firstSM]->SetTitle(Form("SM = %d",ism));
185  }// sm
186 
187  }// prod
188 
189  } // energy bin
190 
192  // Make Chi2 and Difference distributions
193  // Each file an energy bin,
194  // each file contains pads per SM
196  if ( debug ) printf("Plot Chi2-Diff\n");
197 
198  Int_t ncol = 4;
199  Int_t nrow = 3;
200  GetCanvasColRowNumber(nSM, ncol, nrow); // PlotUtils.C
201 
202  TString fileName = "";
203 
204  for(Int_t iebin = 0; iebin < nEBins; iebin++)
205  {
206  if ( debug ) printf("iebin %d",iebin);
207 
208  gStyle->SetOptTitle(1);
209  gStyle->SetPadTopMargin(0.1);
210 
211  // Chi2
212  //
213  TCanvas * cchi2 = new TCanvas
214  (Form("cchi2_ebin%d_%s",iebin,histoName.Data()),
215  Form("Chi2 %2.1f < E < %2.1f, %s",binE[iebin],binE[iebin+1],histoName.Data()),
216  ncol*2000,nrow*2000);
217 
218  cchi2->Divide(ncol,nrow);
219 
220  TLegend *l = new TLegend(-0.04,0,1,1);
221  l->SetFillColor(0);
222  l->SetFillStyle(0);
223  l->SetLineColor(0);
224  l->SetBorderSize(0);
225  l->SetTextSize(0.07);
226  if((histoName.Contains("MM02") || histoName.Contains("MM20")) && !histoName.Contains("NoCut") )
227  l->SetHeader(Form(" %2.1f < #it{E} < %2.1f GeV, #it{n}^{w}_{cells} > 4",binE[iebin],binE[iebin+1]));
228  else
229  l->SetHeader(Form(" %2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
230 
231  for(Int_t ism = firstSM; ism <= lastSM; ism++)
232  {
233  cchi2->cd(ism+1);
234 
235  gPad->SetLogy();
236  gPad->SetGridy();
237 
238  if ( debug ) printf("iE %d ism %d\n",iebin,ism);
239  if ( !hData[iebin][ism-firstSM] ) continue;
240 
241  for(Int_t iprod = 0; iprod < nProd; iprod++)
242  {
243  if ( debug )
244  printf("\t Chi2 prod %d %p %p\n",iprod,hChi2[iprod][iebin][ism-firstSM],hDiff[iprod][iebin][ism-firstSM]);
245 
246  if(!hChi2[iprod][iebin][ism-firstSM]) continue;
247 
248  hChi2[iprod][iebin][ism-firstSM]->SetTitleOffset(1.8,"Y");
249  hChi2[iprod][iebin][ism-firstSM]->SetYTitle("(x_{simu}-x_{Data})^{2}/(#sigma_{x,simu}^{2}+#sigma_{x,data}^{2})");
250  // hChi2[iprod][iebin][ism-firstSM]->SetAxisRange(0,0.8,"X");
251 
252  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[iebin][ism-firstSM]);
253 
254  if(iprod==1) hChi2[iprod][iebin][ism-firstSM]->Draw("H");
255  else hChi2[iprod][iebin][ism-firstSM]->Draw("H same");
256 
257  //hChi2[iprod][iebin][ism-firstSM]->SetMaximum(100);
258  hChi2[iprod][iebin][ism-firstSM]->SetMinimum(1e-2);
259  if(hChi2[1][iebin][ism-firstSM] && hChi2[iprod][iebin][ism-firstSM]->GetMaximum() > hChi2[1][iebin][ism-firstSM]->GetMaximum())
260  hChi2[1][iebin][ism-firstSM]->SetMaximum(hChi2[iprod][iebin][ism-firstSM]->GetMaximum()*1.2);
261 
262  if(ism==0)
263  l->AddEntry(hChi2[iprod][iebin][ism-firstSM],Form("%f",prodParam[iprod]),"PL");
264  } // iprod
265  } // param
266 
267  cchi2->cd(ncol*nrow);
268  l->Draw();
269 
270  fileName = Form("figures/%s/Comparison_Chi2_Ebin%d_%s",
271  histoName.Data(),iebin,inputFileName.Data());
272  fileName+=fileFormat;
273  cchi2->Print(fileName);
274 
275  // Difference
276  //
277  TCanvas * cdiff = new TCanvas(Form("cdiff_ebin%d_%s",
278  iebin,histoName.Data()),
279  Form("MC-Data %2.1f < E < %2.1f, %s",
280  binE[iebin],binE[iebin+1],
281  histoName.Data()),
282  ncol*2000,nrow*2000);
283  cdiff->Divide(ncol,nrow);
284 
285  // TLegend *l = new TLegend(-0.04,0,1,1);
286  // l->SetFillColor(0);
287  // l->SetFillStyle(0);
288  // l->SetLineColor(0);
289  // l->SetBorderSize(0);
290  // l->SetTextSize(0.07);
291  // l->SetHeader(Form(" %2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
292  for(Int_t ism = firstSM; ism <= lastSM; ism++)
293  {
294  cdiff->cd(ism+1);
295 
296  //gPad->SetLogy();
297  gPad->SetGridy();
298 
299  if( debug ) printf("iE %d ism %d\n",iebin,ism);
300  if ( !hData[iebin][ism-firstSM] ) continue;
301 
302  for(Int_t iprod = 0; iprod < nProd; iprod++)
303  {
304  if( debug )
305  printf("\t Diff A prod %d %p\n",iprod,hDiff[iprod][iebin][ism-firstSM]);
306 
307  if(!hDiff[iprod][iebin][ism-firstSM]) continue;
308 
309  hDiff[iprod][iebin][ism-firstSM]->SetTitleOffset(1.8,"Y");
310  //h[iebin][0]->SetLineColor(1);
311  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
312 
313  hDiff[iprod][iebin][ism-firstSM]->SetYTitle("x_{simu}-x_{Data}");
314 
315  // hDiff[iprod][iebin][ism-firstSM]->SetAxisRange(0,0.8,"X");
316 
317  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[iebin][ism-firstSM]);
318 
319  if(iprod==1) hDiff[iprod][iebin][ism-firstSM]->Draw("H");
320  else hDiff[iprod][iebin][ism-firstSM]->Draw("H same");
321 
322  //hDiff[iprod][iebin][ism-firstSM]->SetMaximum(100);
323  //hDiff[iprod][iebin][ism-firstSM]->SetMinimum(1e-2);
324  if(hDiff[1][iebin][ism-firstSM] && hDiff[iprod][iebin][ism-firstSM]->GetMaximum() > hDiff[1][iebin][ism-firstSM]->GetMaximum())
325  hDiff[1][iebin][ism-firstSM]->SetMaximum(hDiff[iprod][iebin][ism-firstSM]->GetMaximum()*1.2);
326 
327  // if(ism==0)
328  // l->AddEntry(hDiff[iprod][iebin][ism-firstSM],Form("%s",prodLeg[iprod].Data()),"PL");
329  } // iprod
330  } // sm
331 
332  cdiff->cd(ncol*nrow);
333  l->Draw();
334 
335  fileName = Form("figures/%s/Comparison_DiffMCData_Ebin%d_%s",
336  histoName.Data(),iebin,inputFileName.Data());
337  fileName+=fileFormat;
338  cdiff->Print(fileName);
339  } // ebin
340 
342  // For each MC different parameter
343  // assing a Chi2 value, do plots with
344  // this distribution and the mininimum
345  // of this distribution
347  if ( debug ) printf ("Compare param\n");
348 
349  const Int_t nX = hData[0][0]->GetNbinsX();
350  const Int_t nXbins = nX;
351  //Double_t * xaxis = hData[0][0]->GetXaxis()->GetXbins()->GetArray();
352  Double_t xaxis[nXbins];
353  Double_t xaxisErr[nXbins];
354  Double_t yaxisErr[nXbins];
355  Double_t yaxis[nXbins];
356  Double_t yaxisA[nXbins];
357 
358  TGraphErrors* gChi2 [nEBins][nSM];
359  TGraphErrors* gChi2Min [nSM];
360  TGraphErrors* gDiffMin [nEBins][nSM];
361 
362  Double_t chi2axis [nProd];
363  Double_t chi2axisErr [nProd];
364 
365  Double_t chi2axisMin [nEBins];
366  Double_t chi2axisMinErr [nEBins];
367 
368  // Get x array
369  for(Int_t ix = 0; ix < nXbins; ix++)
370  {
371  yaxisErr[ix] = 0.05;
372  xaxisErr[ix] = hData[0][0]->GetBinWidth(ix+1)/2;
373  xaxis [ix] = hData[0][0]->GetBinCenter(ix+1);
374  //printf("bin %d, x %2.4f err %2.4f\n",ix,xaxis[ix],xaxisErr[ix]);
375  }
376 
377  //
378  // Fill y array, get min MC-Data of all cases or Chi2 in selected range
379  //
381  for(Int_t ism = firstSM; ism <= lastSM; ism++)
382  {
383  for(Int_t iebin = 0; iebin < nEBins; iebin++)
384  {
385  for(Int_t ix = 0; ix < nXbins; ix++)
386  {
387  yaxis[ix] = -1;
388  Double_t min = 1e6;
389 
390  for(Int_t iprod = 0; iprod < nProd; iprod++)
391  {
392  //if( debug )
393  // printf("\t Diff B prod %d %p\n",iprod,hDiff[iprod][iebin][ism-firstSM]);
394 
395  if ( !hDiff[iprod][iebin][ism-firstSM] ) continue;
396 
397  Double_t content = TMath::Abs(hDiff[iprod][iebin][ism-firstSM]->GetBinContent(ix+1));
398  if(min > content && iprod > 0)
399  {
400  min = content;
401  yaxis[ix] = prodParam[iprod];
402  }
403  } // prod
404  //if(ism==3 && iebin==3) printf("ie %d, ip %d, ix %d, x %2.2f, y %2.2f\n",iebin,ism,ix,xaxis[ix],yaxis[ix]);
405  } // x
406 
407  gDiffMin[iebin][ism-firstSM] = new TGraphErrors(nXbins,xaxis,yaxis,xaxisErr,yaxisErr);
408 
409  //gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetAxisRange(xmin,xmax,"X");
410  if ( hDiff[1][iebin][ism-firstSM] )
411  gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetXTitle(hDiff[1][iebin][ism-firstSM]->GetXaxis()->GetTitle());
412 
413  //
414  // Chi2
415  //
416  Double_t minChi2SM = 10000000;
417  for(Int_t iprod = 1; iprod < nProd-1; iprod++)
418  {
419  if ( !hChi2[iprod][iebin][ism-firstSM] ) continue;
420 
421  Int_t minChi2Bin = hChi2[iprod][iebin][ism-firstSM]->FindBin(xmin);
422  Int_t maxChi2Bin = hChi2[iprod][iebin][ism-firstSM]->FindBin(xmax);
423  Double_t error = 0;
424  if ( maxChi2Bin-minChi2Bin > 0 )
425  {
426  chi2axis [iprod-1] = hChi2[iprod][iebin][ism-firstSM]->
427  IntegralAndError(minChi2Bin,maxChi2Bin,error)/(maxChi2Bin-minChi2Bin);
428  chi2axisErr[iprod-1] = error / (maxChi2Bin-minChi2Bin);
429  // if(ism==3 && iebin > 2 && iebin < 5)
430  // printf("sm %d, ie %d, iprod %d, mu %1.2f, chi2 %2.2f, chierr %2.2f\n",
431  // ism,iebin,iprod,prodParam[iprod-1],chi2axis[iprod],chi2axisErr[iprod]);
432  }
433  else
434  {
435  chi2axis [iprod-1] = 0;
436  chi2axisErr[iprod-1] = 0;
437  }
438 
439  if ( minChi2SM > chi2axis[iprod-1] )
440  {
441  minChi2SM = chi2axis[iprod-1];
442  chi2axisMin [iebin] = prodParam[iprod-1];
443 
444  chi2axisMinErr[iebin] = 0.1;//0.05;
445  }
446  } // prod
447 
448  gChi2[iebin][ism-firstSM] = new TGraphErrors
449  (nProd,prodParam.GetArray(),chi2axis,prodParamE,chi2axisErr);
450  gChi2[iebin][ism-firstSM]->GetHistogram()->SetXTitle("#mu_{1} (%)");
451 
452  } // energy bin
453 
454  gChi2Min[ism-firstSM] = new TGraphErrors(nEBins,binE.GetArray(),chi2axisMin,binEErr,chi2axisMinErr);
455  gChi2Min[ism-firstSM]->GetHistogram()->SetXTitle("E_{cluster} (GeV)");
456 
457  }// sm
458 
460  // Plot difference graphs
462  if ( debug )
463  printf("Plot Diff Graphs!\n");
464 
465  // Diff
466  //
467  TCanvas * cDiffgraph = new TCanvas(Form("cDiff_graph_%s",
468  histoName.Data()),
469  Form("Diff, %s",
470  histoName.Data()),
471  ncol*2000,nrow*2000);
472  cDiffgraph->Divide(ncol,nrow);
473 
474  TLegend *lg = new TLegend(-0.04,0,1,1);
475  lg->SetFillColor(0);
476  lg->SetFillStyle(0);
477  lg->SetLineColor(0);
478  lg->SetBorderSize(0);
479  lg->SetTextSize(0.07);
480  //l->SetHeader(Form(" %2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
481  for(Int_t ism = firstSM; ism <= lastSM; ism++)
482  {
483  cDiffgraph->cd(ism+1);
484 
485  //gPad->SetLogy();
486  gPad->SetGridy();
487  //gPad->SetGridx();
488 
489  //printf("iE %d ism %d\n",iebin,ism);
490  if(!gDiffMin[firstEbin][ism-firstSM]) continue;
491 
492  for(Int_t iebin = firstEbin; iebin < nEBins; iebin++)
493  {
494  //printf("\t prod %d %p\n",iprod,gDiffMin[iebin][ism-firstSM]);
495  if(!gDiffMin[iebin][ism-firstSM]) continue;
496 
497  gDiffMin[iebin][ism-firstSM]->SetMarkerColor(color[iebin]);
498  gDiffMin[iebin][ism-firstSM]->SetLineColor (color[iebin]);
499  gDiffMin[iebin][ism-firstSM]->SetMarkerStyle(24);
500  gDiffMin[iebin][ism-firstSM]->SetMarkerSize(3);
501 
502  gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetTitleOffset(1.8,"Y");
503  //h[iebin][0]->SetLineColor(1);
504  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
505 
506  //hDiff[iprod][iebin][ism-firstSM]->SetYTitle("(x_{simu}-x_{Data})^{2}/(#sigma_{x,simu}^{2}+#sigma_{x,data}^{2})");
507  gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetYTitle("#mu_{1} %");
508  gDiffMin[iebin][ism-firstSM]->SetTitle(Form("SM %d",ism));
509 
510  // gDiffMin[iebin][ism-firstSM]->SetAxisRange(0,0.8,"X");
511 
512  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[iebin][ism-firstSM]);
513 
514  if(iebin==firstEbin)
515  gDiffMin[iebin][ism-firstSM]->Draw("APL");
516  else
517  gDiffMin[iebin][ism-firstSM]->Draw("PL");
518 
519  gDiffMin[iebin][ism-firstSM]->SetMaximum(2);
520  gDiffMin[iebin][ism-firstSM]->SetMinimum(-0.1);
521  // if(hDiff[iprod][iebin][ism-firstSM]->GetMaximum() > hData[iebin][ism-firstSM]->GetMaximum())
522  // hData[iebin][ism-firstSM]->SetMaximum(hDiff[iprod][iebin][ism-firstSM]->GetMaximum()*1.2);
523 
524  if(ism==0)
525  {
526  lg->AddEntry(gDiffMin[iebin][ism-firstSM],
527  Form("%2.1f < #it{E} < %2.1f GeV", binE[iebin], binE[iebin+1]), "PL");
528  }
529 
530  } // iebin
531 
532  //printf("\t end\n");
533  } // param
534 
535  cDiffgraph->cd(ncol*nrow);
536  lg->Draw();
537 
538  fileName = Form("figures/%s/Comparison_MinDiff_PerSM_Ebin_%s",
539  histoName.Data(),inputFileName.Data());
540  fileName+=fileFormat;
541  cDiffgraph->Print(fileName);
542 
543  // Min diff
544  //
545  TCanvas * cDiffgraph2 = new TCanvas(Form("cMinDiff_graph2_%s",
546  histoName.Data()),
547  Form("MinDiff, %s",
548  histoName.Data()),
549  3*2000,2*2000);
550  cDiffgraph2->Divide(3,2);
551 
552  TLegend *lg2 = new TLegend(-0.04,0,1,1);
553  lg2->SetFillColor(0);
554  lg2->SetFillStyle(0);
555  lg2->SetLineColor(0);
556  lg2->SetBorderSize(0);
557  lg2->SetTextSize(0.07);
558  for(Int_t iebin = firstEbin; iebin < nEBins; iebin++)
559  {
560  cDiffgraph2->cd(iebin-2);
561 
562  //gPad->SetLogy();
563  gPad->SetGridy();
564 
565  //printf("iE %d ism %d\n",iebin,ism);
566  if(!gDiffMin[iebin][0]) continue;
567 
568  for(Int_t ism = firstSM; ism <= lastSM; ism++)
569  {
570  //printf("\t prod %d %p\n",iprod,gDiffMin[iebin][ism-firstSM]);
571  if(!gDiffMin[iebin][ism-firstSM]) continue;
572 
573  gDiffMin[iebin][ism-firstSM]->SetMarkerColor(color[ism-firstSM]);
574  gDiffMin[iebin][ism-firstSM]->SetLineColor (color[ism-firstSM]);
575  gDiffMin[iebin][ism-firstSM]->SetMarkerStyle(24);
576  gDiffMin[iebin][ism-firstSM]->SetMarkerSize(3);
577 
578  gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetTitleOffset(1.8,"Y");
579  //h[iebin][0]->SetLineColor(1);
580  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
581 
582  //hDiff[iprod][iebin][ism-firstSM]->SetYTitle("(x_{simu}-x_{Data})^{2}/(#sigma_{x,simu}^{2}+#sigma_{x,data}^{2})");
583  gDiffMin[iebin][ism-firstSM]->GetHistogram()->SetYTitle("#mu_{1} %");
584 
585  gDiffMin[iebin][ism-firstSM]->SetTitle(Form("%2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
586  // gDiffMin[iebin][ism-firstSM]->SetAxisRange(0,0.8,"X");
587 
588  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[iebin][ism-firstSM]);
589 
590  if(ism==0)
591  gDiffMin[iebin][ism-firstSM]->Draw("APL");
592  else
593  gDiffMin[iebin][ism-firstSM]->Draw("PL");
594 
595  gDiffMin[iebin][ism-firstSM]->SetMaximum(2);
596  gDiffMin[iebin][ism-firstSM]->SetMinimum(-0.1);
597  // if(hDiff[iprod][iebin][ism-firstSM]->GetMaximum() > hData[iebin][ism-firstSM]->GetMaximum())
598  // hData[iebin][ism-firstSM]->SetMaximum(hDiff[iprod][iebin][ism-firstSM]->GetMaximum()*1.2);
599 
600  if(iebin==firstEbin)
601  {
602  lg2->AddEntry(gDiffMin[iebin][ism-firstSM],(Form("SM %d",ism)),"PL");
603  }
604 
605  } // iebin
606 
607  //printf("\t end\n");
608  } // param
609 
610  cDiffgraph2->cd(6);
611  lg2->Draw();
612 
613  fileName = Form("figures/%s/Comparison_MinDiff_PerEBin_SM_%s",
614  histoName.Data(),inputFileName.Data());
615  fileName+=fileFormat;
616  cDiffgraph2->Print(fileName);
617 
619  // Plot Chi2 graphs
621  if ( debug )
622  printf("Plot Chi2 Graphs!\n");
623 
624  // Chi2
625  //
626  TCanvas * gChi2graph = new TCanvas(Form("cChi2_graph_%s",
627  histoName.Data()),
628  Form("Chi2, %s",
629  histoName.Data()),
630  ncol*2000,nrow*2000);
631  gChi2graph->Divide(ncol,nrow);
632 
633  TLegend *lgchi2 = new TLegend(-0.04,0,1,1);
634  lgchi2->SetFillColor(0);
635  lgchi2->SetFillStyle(0);
636  lgchi2->SetLineColor(0);
637  lgchi2->SetBorderSize(0);
638  lgchi2->SetTextSize(0.07);
639  //l->SetHeader(Form(" %2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
640  for(Int_t ism = firstSM; ism <= lastSM; ism++)
641  {
642  gChi2graph->cd(ism+1);
643 
644  gPad->SetLogy();
645  gPad->SetGridy();
646  gPad->SetGridx();
647 
648  //printf("iE %d ism %d\n",iebin,ism);
649  if(!gChi2[firstEbin][ism-firstSM]) continue;
650 
651  for(Int_t iebin = firstEbin; iebin < nEBins; iebin++)
652  {
653  //printf("\t prod %d %p\n",iprod,gChi2[iebin][ism-firstSM]);
654  if(!gChi2[iebin][ism-firstSM]) continue;
655 
656  gChi2[iebin][ism-firstSM]->SetMarkerColor(color[iebin]);
657  gChi2[iebin][ism-firstSM]->SetLineColor (color[iebin]);
658  gChi2[iebin][ism-firstSM]->SetMarkerStyle(24);
659  gChi2[iebin][ism-firstSM]->SetMarkerSize(3);
660 
661  gChi2[iebin][ism-firstSM]->GetHistogram()->SetTitleOffset(1.6,"Y");
662  //h[iebin][0]->SetLineColor(1);
663  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
664 
665  //hDiff[iprod][iebin][ism-firstSM]->SetYTitle("(x_{simu}-x_{Data})^{2}/(#sigma_{x,simu}^{2}+#sigma_{x,data}^{2})");
666  gChi2[iebin][ism-firstSM]->GetHistogram()->SetYTitle("#chi^{2} / #nu");
667  gChi2[iebin][ism-firstSM]->SetTitle(Form("SM %d",ism));
668 
669  // gChi2[iebin][ism-firstSM]->SetAxisRange(0,0.8,"X");
670 
671  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[iebin][ism-firstSM]);
672 
673  if(iebin==firstEbin)
674  gChi2[iebin][ism-firstSM]->Draw("APL");
675  else
676  gChi2[iebin][ism-firstSM]->Draw("PL");
677 
678  gChi2[iebin][ism-firstSM]->SetMaximum(200);
679  gChi2[iebin][ism-firstSM]->SetMinimum(0.25);
680  // if(hDiff[iprod][iebin][ism-firstSM]->GetMaximum() > hData[iebin][ism-firstSM]->GetMaximum())
681  // hData[iebin][ism-firstSM]->SetMaximum(hDiff[iprod][iebin][ism-firstSM]->GetMaximum()*1.2);
682 
683  if(ism==0)
684  {
685  lgchi2->AddEntry(gChi2[iebin][ism-firstSM],
686  Form("%2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]),"PL");
687  }
688 
689  } // iebin
690 
691  //printf("\t end\n");
692  } // param
693 
694  gChi2graph->cd(ncol*nrow);
695  lgchi2->Draw();
696 
697  fileName = Form("figures/%s/Comparison_Chi2_PerSM_Ebin_%s",
698  histoName.Data(),inputFileName.Data());
699  fileName+=fileFormat;
700  gChi2graph->Print(fileName);
701 
702  // Min Chi2
703  //
704  TCanvas * gChi2Mingraph = new TCanvas(Form("cChi2Min_graph_%s",
705  histoName.Data()),
706  Form("Diff, %s",
707  histoName.Data()),
708  ncol*2000,nrow*2000);
709  gChi2Mingraph->Divide(ncol,nrow);
710 
711  TLegend *lgChi2Min = new TLegend(-0.04,0,1,1);
712  lgChi2Min->SetFillColor(0);
713  lgChi2Min->SetFillStyle(0);
714  lgChi2Min->SetLineColor(0);
715  lgChi2Min->SetBorderSize(0);
716  lgChi2Min->SetTextSize(0.07);
717 
718  for(Int_t ism = firstSM; ism <= lastSM; ism++)
719  {
720  gChi2Mingraph->cd(ism+1);
721 
722  //gPad->SetLogy();
723  //gPad->SetGridy();
724 
725  //printf("iE %d ism %d\n",iebin,ism);
726  if(!gChi2Min[ism-firstSM]) continue;
727 
728  //printf("\t prod %d %p\n",iprod,gChi2Min[ism-firstSM]);
729  if(!gChi2Min[ism-firstSM]) continue;
730 
731  gChi2Min[ism-firstSM]->SetMarkerColor(color[ism-firstSM]);
732  gChi2Min[ism-firstSM]->SetLineColor (color[ism-firstSM]);
733  gChi2Min[ism-firstSM]->SetMarkerStyle(24);
734  gChi2Min[ism-firstSM]->SetMarkerSize(3);
735 
736  gChi2Min[ism-firstSM]->GetHistogram()->SetTitleOffset(1.8,"Y");
737  //hData->SetLineColor(1);
738  //hData->SetAxisRange(0.1,2.5,"X");
739 
740  //hDiff[iprod][ism-firstSM]->SetYTitle("(x_{simu}-x_{Data})^{2}/(#sigma_{x,simu}^{2}+#sigma_{x,data}^{2})");
741  gChi2Min[ism-firstSM]->GetHistogram()->SetYTitle("#mu_{1} (%)");
742  gChi2Min[ism-firstSM]->SetTitle(Form("SM %d",ism));
743 
744  // if(firstEbin == 3) gChi2Min[ism-firstSM]->GetHistogram()->SetAxisRange(7,20,"X");
745  // else gChi2Min[ism-firstSM]->GetHistogram()->SetAxisRange(15,60,"X");
746 
747  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,hData[ism-firstSM]);
748 
749  gChi2Min[ism-firstSM]->Draw("APL");
750 
751  gChi2Min[ism-firstSM]->SetMaximum(1.7);
752  gChi2Min[ism-firstSM]->SetMinimum(-0.1);
753 
754  } // param
755 
756  gChi2Mingraph->cd(ncol*nrow);
757  lgChi2Min->Draw();
758 
759  fileName = Form("figures/%s/Comparison_Chi2Min_PerSM_Ebin_%s",
760  histoName.Data(),inputFileName.Data());
761  fileName+=fileFormat;
762  gChi2Mingraph->Print(fileName);
763 
764  // Write to file
765  TString inputFileName2 = inputFileName;
766  inputFileName2.ReplaceAll("ForChi2","");
767 
768  TFile* fout = new TFile(Form("figures/%s/Chi2Distributions_%s_PerSM_Ebin.root",
769  histoName.Data(),inputFileName2.Data()),"recreate");
770  for(Int_t ism = firstSM; ism <= lastSM; ism++)
771  {
772  gChi2Min[ism-firstSM]->SetName(Form("Chi2Min_SM%d",ism));
773  gChi2Min[ism-firstSM]->Write();
774  for(Int_t iebin = firstEbin; iebin < nEBins; iebin++)
775  {
776  gChi2[iebin][ism-firstSM]->SetName(Form("Chi2_SM%d_EBin%d",ism,iebin));
777  gChi2[iebin][ism-firstSM]->SetTitle(Form("%2.1f < #it{E} < %2.1f GeV",binE[iebin],binE[iebin+1]));
778  gChi2[iebin][ism-firstSM]->Write();
779  }
780  }
781  fout->Close();
782 
783 }
784 
785 //------------------------------------------------------------------------------
793 //------------------------------------------------------------------------------
794 void Chi2( TH1D* h, TH1D* hD, TH1D*& hChi2, TH1D*& hDiff, Bool_t debug)
795 {
796  if ( !h || !hD ) return ;
797 
798  hChi2 = (TH1D*) h->Clone(Form("%s_Chi2",h->GetName()));
799  hDiff = (TH1D*) h->Clone(Form("%s_Diff",h->GetName()));
800 
801  //printf("%p %p, N bins %d\n",h,hD,h->GetNbinsX());
802 
803  for(Int_t ibin = 1; ibin < h->GetNbinsX(); ibin++)
804  {
805  //printf("bin %d\n",ibin);
806  Double_t content = h ->GetBinContent(ibin);
807  Double_t contentD = hD->GetBinContent(ibin);
808  Double_t value2 = content-contentD;
809  value2*=value2;
810 
811  Double_t econtent2 = h->GetBinError(ibin);
812  Double_t econtent2D = hD->GetBinError(ibin);
813 
814  Double_t esum = econtent2D+econtent2;
815  esum *=esum;
816 
817  econtent2 *=econtent2;
818  econtent2D*=econtent2D;
819 
820  Double_t esum2 = econtent2D+econtent2;
821 
822  if ( esum2 > 1e-19 )
823  {
824  if(debug) printf("ibin %d, value2 %e, esum %e\n",ibin,value2,esum2);
825 
826  hDiff->SetBinContent(ibin,content-contentD);
827  hDiff->SetBinError (ibin,h->GetBinError(ibin)+hD->GetBinError(ibin));
828 
829  value2/=esum2;
830  esum/=esum2;
831 
832  if ( debug )
833  printf("\t ibin %d, x %2.2f, vS %2.4f, vD %2.4f, |vS-vD| %2.4e |vS-vD|^2 %2.4e, Err|vS-vD|^2 %2.4e\n",
834  ibin, hD->GetBinCenter(ibin), content, contentD, TMath::Abs(content-contentD),value2,esum);
835 
836  hChi2->SetBinContent(ibin,value2);
837  hChi2->SetBinError(ibin,esum);
838  }
839  else
840  {
841  hChi2->SetBinContent(ibin,0);
842  hChi2->SetBinError(ibin,0);
843  hDiff->SetBinContent(ibin,0);
844  hDiff->SetBinError(ibin,0);
845  }
846 
847  }
848  if ( debug )
849  printf("\t func %p %p\n",hChi2,hDiff);
850 }
851 
852 
color
Int_t color[]
print message on plot with ok/not ok
Definition: DrawAnaCaloTrackQA.C:99
Double_t
double Double_t
Definition: External.C:58
fileName
TString fileName
Definition: ProjectCombinHFAndFit.C:25
lastSM
int lastSM
Definition: HVrecalculation.C:79
CalculateParamChi2MCxTalkDataPerSM
void CalculateParamChi2MCxTalkDataPerSM(TString histoName="SMM02NoCut", TString inputFileName="", Int_t dataProd=0, Int_t firstSM=0, Int_t lastSM=9, TArrayD binE=0, TArrayD prodParam=0, Int_t rebin=1, Bool_t debug=kFALSE)
Definition: CalculateParamChi2MCxTalkDataPerSM.C:59
Chi2
void Chi2(TH1D *h, TH1D *hD, TH1D *&hChi2, TH1D *&hDiff, Bool_t debug=kFALSE)
Definition: CalculateParamChi2MCxTalkDataPerSM.C:794
fileFormat
TString fileFormat
Definition: CalculateParamChi2MCxTalkDataPerSM.C:40
Int_t
int Int_t
Definition: External.C:63
nProd
const Int_t nProd
Definition: DrawProductionComparison.C:102
Float_t
float Float_t
Definition: External.C:68
TH1D
Definition: External.C:212
nEBins
const Int_t nEBins
Definition: CalculateShowerShapeFractionPerNCellCut.C:52
firstEbin
Int_t firstEbin
Definition: CalculateParamChi2MCxTalkDataPerSM.C:42
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:89
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
binE
Double_t binE[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:53
GetCanvasColRowNumber
void GetCanvasColRowNumber(Int_t npad, Int_t &ncol, Int_t &nrow, Bool_t square=kFALSE)
Definition: PlotUtils.C:304
Bool_t
bool Bool_t
Definition: External.C:53
fout
TFile * fout
input train file
Definition: DrawAnaCaloTrackQA.C:90
binEErr
Double_t binEErr[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:54