AliPhysics  cda3415 (cda3415)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BadChannelAnalysis.C
Go to the documentation of this file.
1 // This macro has been developed to find badcells candidates in EMCal based on cells amplitude distributions
2 // Input needed can be either outputs QA from AliAnaCalorimeterQA task (BadChannelAnalysis() function)
3 // Or from merged output of AliAnalysisTaskCaloCellsQA (use BCAnalysis() function)
4 //
5 // Author : Alexis Mas (SUBATECH) & M. Germain, based on getCellsRunQA.C from Olga Driga (SUBATECH)
6 //
7 //-----------------
8 // Main method:
9 //---------------
10 // BadChannelAnalysis:
11 //
12 // step 1 : Convert()
13 // read list of mergeable runs in your working directory (in example below the $workdir is "/scratch/alicehp2/germain/QANew2/"
14 // The QAresults.root files should be aleady copied from alien and be in $workdir/<period>/<pass>/runnb.root
15 // read/merge the histos"EMCAL_hAmpId" and "EMCAL_hTimeId" from QAresults.root file and write them in
16 // $workdir/period/pass/<period><pass>Runlist0New.root" !!! this is hardcoded !!!!!
17 // step 1 has to be called only the first time runing on a new list
18 //
19 // step 2 BCanalysis() main method to analyse previously created file (hardcoded)
20 //
21 // call of different Periodanalysis(criterium,..) functions according to the wanted tests (critreria)
22 // 1 : average E for E>Emin
23 // 2 : entries for E>Emin
24 // 3 : ki²/ndf (from fit of each cell Amplitude between Emin and Emax)
25 // 4 : A parameter (from fit of each cell Amplitude between Emin and Emax)
26 // 5 : B parameter (from fit of each cell Amplitude between Emin and Emax)
27 // 6 :
28 // 7 : give bad + dead list
29 //
30 // Mainly used: 1 and 2 (with different settings (chi2, intervals of energy : this is quite dependent of the stat you may have )
31 // Further improvement: implement tests 1 and 2 on time distribustion histogram
32 //
33 // ---------------------
34 // Running the macro
35 // ---------------------
36 // root [2] .L BadChannelAnalysis.C++
37 // root [2] BadChannelAnalysis("EMCAL","LHC15o","muon_caloLego","AnyINTnoBC",trial=0)
38 //
39 // !!! pay attention the trigger name depends on the caloQA_triggername you want to analyse check first in QAresults.root what is abvailable
40 // !!! it is generally not good to run it on triggered data for the following reasons:
41 //
42 // ------------------
43 // outputs
44 //-------------------
45 // the output of this analysis povides you:
46 // intermediate steps files: (those will be recreated each time you rerun so pa attention to save the different files when changing period/listof runs....
47 // - Criterum-xx_Emin-xx_Emax-xx.txt : list of identified bad for the different test/Emin/Emax
48 // - <period><pass>.txt file with list of dead/bad cells identified
49 // - a pdf file with all energy distributions plots of all bad cells candidates (compared to a reference one (hard coded see Draw function to change)
50 //
51 // Further improvement: implement tests 1 and 2 on time distribution histogram
52 //
54 
55 
56 
57 
58 #if !defined(__CINT__) || defined(__MAKECINT__)
59 #include <TFile.h>
60 #include <TH1F.h>
61 #include <TH2F.h>
62 #include <TH3D.h>
63 #include <TLine.h>
64 #include <Riostream.h>
65 #include <TCanvas.h>
66 #include <TGraphErrors.h>
67 #include <TGrid.h>
68 #include <TStyle.h>
69 #include <TFileMerger.h>
70 #include <TMultiGraph.h>
71 #include <TROOT.h>
72 #include <TLegend.h>
73 #include <TString.h>
74 #include <TGridCollection.h>
75 #include <TGridResult.h>
76 #include <TClonesArray.h>
77 #include <TObjString.h>
78 #include <stdio.h>
79 #include <fstream>
80 #include <iostream>
81 #include <alloca.h>
82 #include <string>
83 #include <cstring>
84 #endif
85 using namespace std;
86 
87 void Draw2(Int_t cell, Int_t cellref=400) {
88  gROOT->SetStyle("Plain");
89  gStyle->SetOptStat(0);
90  gStyle->SetFillColor(kWhite);
91  gStyle->SetTitleFillColor(kWhite);
92  gStyle->SetPalette(1);
93  char out[120]; char title[100]; char name[100];char name2[100];
94  TString slide(Form("Cells %d-%d",cell,cell));
95 
96  sprintf(out,"%d.gif",cell);
97  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
98  TH1 *hCellref = hCellAmplitude->ProjectionX("badcells",cellref+1,cellref+1);
99 
100  TCanvas *c1 = new TCanvas("badcells","badcells",600,600) ;
101  c1->SetLogy();
102 
103  // hCellref->Rebin(3);
104  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
105 
106  sprintf(name,"Cell %d",cell) ;
107  TH1 *hCell = hCellAmplitude->ProjectionX(name,cell+1,cell+1);
108 
109  sprintf(title,"Cell %d Entries : %d Enties ref: %d",cell, (Int_t)hCell->GetEntries(),(Int_t)hCellref->GetEntries()) ;
110  hCell->SetLineColor(2) ;
111  // cout<<title<<endl ;
112  hCell->SetMaximum(1e5);
113  // hCell->Rebin(3);
114  hCell->SetAxisRange(0.,10.);
115  hCell->GetXaxis()->SetTitle("E (GeV)");
116  hCell->GetYaxis()->SetTitle("N Entries");
117  hCellref->SetAxisRange(0.,10.);
118  hCell->SetLineWidth(1) ;
119  hCellref->SetLineWidth(1) ;
120  hCell->SetTitle(title);
121  hCellref->SetLineColor(1) ;
122  leg->AddEntry(hCellref,"reference","l");
123  leg->AddEntry(hCell,"current","l");
124  hCell->Draw() ;
125  hCellref->Draw("same") ;
126  leg->Draw();
127  sprintf(name2,"Cell%dLHC13MB.gif",cell) ;
128  c1->SaveAs(name2);
129 
130 }
131 
132 void Draw(Int_t cell[], Int_t iBC, Int_t nBC,TString datapath="/scratch/alicehp2/germain/QANew2/", TString period="LHC15f", TString pass="pass2", Int_t trial=0,const Int_t cellref=2377){
133  //Allow to produce a pdf file with badcells candidates (red) compared to a refence cell (black)
134 
135  gROOT->SetStyle("Plain");
136  gStyle->SetOptStat(0);
137  gStyle->SetFillColor(kWhite);
138  gStyle->SetTitleFillColor(kWhite);
139  gStyle->SetPalette(1);
140  char out[120]; char title[100]; char name[100];
141  TString slide(Form("Cells %d-%d",cell[iBC],cell[iBC+nBC-1]));
142 
143  TString reflegend = Form("reference Cell %i",cellref);
144  sprintf(out,"%d-%d.gif",cell[iBC],cell[iBC+nBC-1]);
145  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
146  TH1 *hCellref = hCellAmplitude->ProjectionX("badcells",cellref+1,cellref+1);
147  Int_t i;
148  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750) ;
149  c1->SetLogy();
150  if(nBC > 6) c1->Divide(3,3) ;
151  else if (nBC > 3) c1->Divide(3,2) ;
152  else c1->Divide(3,1);
153  // hCellref->Rebin(3);
154  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
155  for(i=0; i<nBC ; i++){
156  sprintf(name,"Cell %d",cell[iBC+i]) ;
157  TH1 *hCell = hCellAmplitude->ProjectionX(name,cell[iBC+i]+1,cell[iBC+i]+1);
158 
159  c1->cd(i%9 + 1) ;
160  c1->cd(i%9 + 1)->SetLogy();
161  sprintf(title,"Cell %d Entries : %d Ref : %d",cell[iBC+i], (Int_t)hCell->GetEntries(), (Int_t)hCellref->GetEntries() ) ;
162  hCell->SetLineColor(2) ;
163  // cout<<title<<endl ;
164  hCell->SetMaximum(1e6);
165  // hCell->Rebin(3);
166  hCell->SetAxisRange(0.,10.);
167  hCell->GetXaxis()->SetTitle("E (GeV)");
168  hCell->GetYaxis()->SetTitle("N Entries");
169  hCellref->SetAxisRange(0.,8.);
170  hCell->SetLineWidth(1) ;
171  hCellref->SetLineWidth(1) ;
172  hCell->SetTitle(title);
173  hCellref->SetLineColor(1) ;
174 
175 
176 
177  if(i==0){
178 
179  leg->AddEntry(hCellref,reflegend,"l");
180  leg->AddEntry(hCell,"current","l");
181  }
182  hCell->Draw() ;
183  hCellref->Draw("same") ;
184  leg->Draw();
185  }
186 
187  //CREATE A PDF FILE
188  TString PdfName = Form("%s/BadChannelNew/2016/%s%sList1Test%i.pdf",datapath.Data(),period.Data(),pass.Data(),trial);
189 
190  // TString PdfName = Form("/scratch/alicehp2/germain/QANew2/BadChannelNew/2016/%s%sList0Test%i.pdf",period.Data(),pass.Data(),trial);
191  if(nBC<9) {
192  PdfName +=")";
193  c1->Print(PdfName.Data());}
194  else if(iBC==0) {
195  PdfName +="(";
196  c1->Print(PdfName.Data());}
197  else c1->Print(PdfName.Data());
198 
199 
200  // c1->SaveAs(out);
201  delete hCellref ;
202  delete c1 ;
203  delete leg ;
204 
205 }
206 
207 //_________________________________________________________________________
208 //_________________________________________________________________________
209 
210 void Convert(TString datapath= "/scratch/alicehp2/germain/QANew2",TString fCalorimeter = "EMCAL", TString period = "LHC11h", TString pass = "pass1_HLT", TString trigger= "default"){
211 
212  //Create one file for the analysis from several outputs QA files listed in runlist.txt
213  //You need :
214  // runlist.txt with runs listed
215  // outputsQA e.g period/pass/123456.root
216 
217  TH2F *hCellAmplitude = new TH2F("hCellAmplitude","Cell Amplitude",200,0,10,23040,0,23040);
218  TH2F *hCellTime = new TH2F("hCellTime","Cell Time",250,-275,975,23040,0,23040);
219 
220  TH1D *hNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun","Number of processed events vs run number",200000,100000,300000);
221  FILE * pFile;
222 
223  TString file = Form("%s/%s%sBC1.txt",datapath.Data(),period.Data(),pass.Data());
224 
225  // TString file = Form("/scratch/alicehp2/germain/QANew2/%s%sBC0.txt",period.Data(),pass.Data());
226 cout << " file = " << file << endl;;
227  pFile = fopen(file.Data(), "r"); //open the text file where include the run list and correct run index
228 
229 
230  cout << " fcalo: " << fCalorimeter << "; period: " << period << "; pass: " << pass << " trigger "<<trigger<< endl;
231 
232  Int_t ix,iy;
233  Int_t Nentr;
234  Int_t p;
235  Int_t q;
236  Int_t ncols;
237  Int_t nlines = 0 ;
238  Int_t RunId[500] ;
239  Double_t x[500] ;
240  Double_t xrun[500] ;
241  while (1){
242  ncols = fscanf(pFile," %d ",&q);
243  if (ncols< 0) break;
244  // x[nlines]=p;
245  RunId[nlines]=q;
246  //xrun[nlines]=1.*q;
247  nlines++;
248  }
249  fclose(pFile);
250  const Int_t nRun = nlines ;
251  Double_t content;
252  TString base ;
253  TString base2 ;
254  TString BCfile ;
255 
256 
257  TString direct(Form("CaloQA_%s",trigger.Data()));
258  // TString direct ="EMCAL_TrigAnyINT_Cl" ;
259 
260  for(Int_t i = 0 ; i < nRun ; i++) {
261  // base2 = Form("/scratch/alicehp2/germain/QANew2/%s/%s/",period.Data(),pass.Data());
262  // base = Form("/scratch/alicehp2/germain/QANew2/%s/%s/%d",period.Data(),pass.Data(),RunId[i]);
263  base2 = Form("%s/%s/%s/",datapath.Data(),period.Data(),pass.Data());
264  base = Form("%s/%s/%s/%d",datapath.Data(),period.Data(),pass.Data(),RunId[i]);
265  BCfile=Form("%s%s%sRunlist1New.root",base2.Data(),period.Data(),pass.Data());
266 
267  TString infile ;
268 
269  // ICI on met le nom period/pass/runblabla.root
270  if ((pass=="cpass1_pass2")||(pass=="cpass1-2")){
271  if (trigger=="default"){
272  infile = Form("%s_barrel.root",base.Data());}
273  else {infile = Form("%s_outer.root",base.Data());}
274  }
275  else
276  infile = Form("%s.root",base.Data()) ;
277  cout<<"file : "<<infile<<endl;
278  TFile *f = TFile::Open(infile);
279  // base += "/" ;
280  // base += trigger ;
281  base=Form("%s/%s",base.Data(),trigger.Data());
282  if (!f) continue;
283  cout << " jusqu'ici ca va "<< endl;
284  TDirectoryFile *dir = (TDirectoryFile *)f->Get(direct);
285  if (!dir) continue;
286  cout << " jusqu'ici ca va dir OK "<< endl;
287  TList *outputList = (TList*)dir->Get(direct);
288  if (!outputList) continue;
289  // TList *outputList = (TList *)f->Get(direct.Data());
290  cout << " jusqu'ici ca va List OK "<< endl;
291 
292  TH2F *hAmpId;
293  TH2F *hNEvents;
294  TH2F * hTimeId;
295 
296  hAmpId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
297  hAmpId->Draw();
298 
299  hTimeId =(TH2F *)outputList->FindObject("EMCAL_hTimeId");
300  hTimeId->Draw();
301 
302  cout<<"file : "<<infile<<endl;
303  hNEvents =(TH2F *)outputList->FindObject("hNEvents");
304  if (!hNEvents)continue;
305  Nentr = (Int_t)hNEvents->GetEntries();
306  cout << " N entries " << Nentr << endl;
307  if (Nentr<100) continue ;
308  hNEventsProcessedPerRun->SetBinContent(RunId[i]-100000,(Double_t)Nentr);
309 
310 
311 
312  hCellAmplitude->Add(hAmpId);
313  hCellTime->Add(hTimeId);
314 
315 
316  //cout<<i<<endl;
317  if(i==0){ cout<<"Merging/Converting procedure ..." ; cout.flush();}
318  else { cout<<"..." ; cout.flush();}
319  outputList->Delete();
320  dir->Delete();
321  f->Close();
322  delete f;
323  }
324 
325 
326  TFile *BCF = TFile::Open(BCfile,"recreate");
327  hNEventsProcessedPerRun->Write();
328  hCellAmplitude->Write();
329  hCellTime->Write();
330  BCF->Close();
331  cout<<"DONE !"<<endl;
332 
333 }
334 
335 //_________________________________________________________________________
336 //_________________________________________________________________________
337 
338 void Process(Int_t *pflag[23040][7], TH1* inhisto, Double_t Nsigma = 4., Int_t dnbins = 200, Double_t dmaxval = -1., Int_t compteur = 1)
339 {
340  // 1) create a distribution for the input histogram;
341  // 2) fit the distribution with a gaussian
342  // 3) define good area within +-Nsigma to identfy badcells
343  //
344  // inhisto -- input histogram;
345  // dnbins -- number of bins in distribution;
346  // dmaxval -- maximum value on distribution histogram.
347 
348  gStyle->SetOptStat(1); // MG modif
349  gStyle->SetOptFit(1); // MG modif
350  Int_t crit = *pflag[0][0] ; //identify the criterum processed
351  Double_t goodmax= 0. ;
352  Double_t goodmin= 0. ;
353  *pflag[0][0] =1;
354  if (dmaxval < 0.) {
355  dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
356  if(crit==2 && dmaxval > 1) dmaxval =1. ;
357  }
358 
359  TH1 *distrib = new TH1F(Form("%sDistr",inhisto->GetName()), "", dnbins, inhisto->GetMinimum(), dmaxval);
360  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
361  distrib->SetYTitle("Entries");
362 
363  // fill distribution
364  for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
365  distrib->Fill(inhisto->GetBinContent(c));
366 
367  // draw histogram + distribution
368  TCanvas *c1 = new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400);
369  c1->Divide(2,1);
370 
371  c1->cd(1);
372  gPad->SetLeftMargin(0.14);
373  gPad->SetRightMargin(0.06);
374  gPad->SetLogy();
375  inhisto->SetTitleOffset(1.7,"Y");
376  inhisto->Draw();
377 
378  c1->cd(2);
379  gPad->SetLeftMargin(0.14);
380  gPad->SetRightMargin(0.10);
381  gPad->SetLogy();
382  distrib->Draw();
383 
384  Int_t higherbin=0,i;
385  for (i = 2; i <= distrib->GetNbinsX(); i++){
386  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i))
387  higherbin = i ;
388  }
389 
390  for(i=higherbin ; i<=dnbins ; i++)
391  if(distrib->GetBinContent(i)<2) break ;
392  goodmax = distrib->GetBinCenter(i);
393 
394  for(i=higherbin ; i>0 ; i--)
395  if(distrib->GetBinContent(i)<2) break ;
396  goodmin = distrib->GetBinLowEdge(i);
397 
398  TF1 *fit2 = new TF1("fit2", "gaus");
399 
400  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
401  Double_t sig, mean, chi2ndf;
402  // Marie midif to take into account very non gaussian distrig
403  mean = fit2->GetParameter(1);
404  sig = fit2->GetParameter(2);
405  chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
406 
407 
408 // mean = distrib->GetMean(); // MG
409 // sig = distrib->GetRMS(); // MG
410 
411  // cout << "----------------------------------------------"<< endl;
412 // cout <<" pass " << compteur << " mean " << mean << " rms" << sig << endl;
413 // cout << "----------------------------------------------"<< endl;
414  if (mean <0.) mean=0.;
415 
416  goodmin = mean - Nsigma*sig ;
417  goodmax = mean + Nsigma*sig ;
418 
419  if (goodmin<0) goodmin=0.;
420  // if (compteur==0){
421 // goodmin = 0.;
422 // goodmax = 2.*mean;
423 // }
424 
425  cout << "-----------------------------------------------------------------"<< endl;
426  cout << " pass " << compteur << " mean " << mean << " rms" << sig << " goodmin " << goodmin<< " goodmax" << goodmax << endl;
427  cout << "-----------------------------------------------------------------"<< endl;
428 
429  // lines
430  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
431  lline->SetLineColor(kOrange);
432  lline->SetLineStyle(7);
433  lline->Draw();
434 
435  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
436  rline->SetLineColor(kOrange);
437  rline->SetLineStyle(7);
438  rline->Draw();
439 
440  // legend
441  TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
442  leg->AddEntry(lline, "Good region boundary","l");
443  leg->Draw("same");
444  // fit2->Draw("same");
445 
446  c1->Update();
447  TString name = "criteria-" ;
448  name+= crit;
449  name+= ".gif";
450 
451  c1->SaveAs(name);
452 
453  Int_t ntot = 0, cel;
454 
455  for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++) {
456  cel=(Int_t)(inhisto->GetBinLowEdge(c)+0.1);
457  if (inhisto->GetBinContent(c) <= goodmin) {
458  ntot++;
459  *pflag[cel][crit]=0;
460  }
461  else if (inhisto->GetBinContent(c) > goodmax) {
462  ntot++;
463  *pflag[cel][crit]=2;
464  }
465  }
466 }
467 
468 //_________________________________________________________________________
469 //_________________________________________________________________________
470 
471 void TestCellEandN(Int_t *pflag[23040][7], Double_t Emin = 0.1, Double_t Emax=2., Double_t Nsigma = 4., Int_t compteur = 1, char const * hname = "hCellAmplitude", Int_t dnbins = 200)
472 {
473 
474 
475  // Three more tests for bad cells:
476  // 1) total deposited energy;
477  // 2) total number of entries;
478  // 3) average energy = [total deposited energy]/[total number of entries].
479  //
480 
481  // Int_t count;
482 // count = compteur;
483 
484  // input; X axis -- absId numbers
485 
486  TH2 *hCellAmplitude = (TH2*) gFile->Get(hname);
487 
488  // binning parameters
489  Int_t ncells = hCellAmplitude->GetNbinsY();
490  Double_t amin = hCellAmplitude->GetYaxis()->GetXmin();
491  Double_t amax = hCellAmplitude->GetYaxis()->GetXmax();
492 
493  cout << "ncells " << ncells << " amin = " << amin << "amax = " << amax<< endl;
494 
495 
496  TH1* hCellEtotal = new TH1F(Form("%s_hCellEtotal_E%.2f",hname,Emin),
497  Form("Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax);
498  hCellEtotal->SetXTitle("AbsId");
499  hCellEtotal->SetYTitle("Energy, GeV");
500 
501  TH1F *hCellNtotal = new TH1F(Form("%s_hCellNtotal_E%.2f",hname,Emin),
502  Form("Number of entries per events, E > %.2f GeV",Emin), ncells,amin,amax);
503  hCellNtotal->SetXTitle("AbsId");
504  hCellNtotal->SetYTitle("Entries");
505 
506  TH1F *hCellEtoNtotal = new TH1F(Form("%s_hCellEtoNtotal_E%.2f",hname,Emin),
507  Form("Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax);
508  hCellEtoNtotal->SetXTitle("AbsId");
509  hCellEtoNtotal->SetYTitle("Energy, GeV");
510 
511  TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
512  Double_t totalevents = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
513 
514  // fill cells
515  for (Int_t c = 1; c <= ncells; c++) {
516  Double_t Esum = 0;
517  Double_t Nsum = 0;
518 
519 
520  for (Int_t j = 1; j <= hCellAmplitude->GetNbinsX(); j++) {
521  // for (Int_t j = hCellAmplitude->GetXaxis()->FindBin(Emin); j <= hCellAmplitude->GetXaxis()->FindBin(Emax); j++) {
522  Double_t E = hCellAmplitude->GetXaxis()->GetBinCenter(j);
523  Double_t N = hCellAmplitude->GetBinContent(j, c);
524  if (E < Emin || E>Emax) continue;
525  // if (E > Emin && E< Emax) {
526  Esum += E*N;
527  Nsum += N;
528  //}
529  }
530 
531  hCellEtotal->SetBinContent(c, Esum);
532  hCellNtotal->SetBinContent(c, Nsum/totalevents);
533  // hCellNtotal->SetBinContent(c, Nsum);
534 
535  if (Nsum > 0.) // number of entries >= 1
536  hCellEtoNtotal->SetBinContent(c, Esum/Nsum);
537 
538  }
539 
540  delete hCellAmplitude;
541 
542  // Process(hCellEtotal, dnbins );
543  if(*pflag[0][0]==1)
544  Process(pflag, hCellEtoNtotal, Nsigma, dnbins, -1,compteur);
545  if(*pflag[0][0]==2)
546  Process(pflag, hCellNtotal, Nsigma, dnbins,-1, compteur);
547 }
548 
549 //_________________________________________________________________________
550 //_________________________________________________________________________
551 
552 void TestCellShapes(Int_t *pflag[23040][7], Double_t fitEmin, Double_t fitEmax, Double_t Nsigma =4., Int_t compteur= 1, char const * hname = "hCellAmplitude", Int_t dnbins = 1000)
553 {
554  // Test cells shape using fit function f(x)=A*exp(-B*x)/x^2.
555  // Produce values per cell + distributions for A,B and chi2/ndf parameters.
556 // Int_t count;
557 // count = compteur;
558 
559  TH2 *hCellAmplitude = (TH2*) gFile->Get(Form("%s",hname));
560 
561  // binning parameters
562  Int_t ncells = hCellAmplitude->GetNbinsY();
563  Double_t amin = hCellAmplitude->GetYaxis()->GetXmin();
564  Double_t amax = hCellAmplitude->GetYaxis()->GetXmax();
565  cout << "ncells " << ncells << " amin = " << amin << "amax = " << amax<< endl;
566 
567  // initialize histograms
568  TH1 *hFitA = new TH1F(Form("hFitA_%s",hname),"Fit A value", ncells,amin,amax);
569  hFitA->SetXTitle("AbsId");
570  hFitA->SetYTitle("A");
571 
572  TH1 *hFitB = new TH1F(Form("hFitB_%s",hname),"Fit B value", ncells,amin,amax);
573  hFitB->SetXTitle("AbsId");
574  hFitB->SetYTitle("B");
575 
576  TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",hname),"Fit #chi^{2}/ndf value", ncells,amin,amax);
577  hFitChi2Ndf->SetXTitle("AbsId");
578  hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
579 
580  Double_t maxval1=0., maxval2=0., maxval3=0.;
581  Double_t prev=0., MSA=0., AvA = 0. ; //those param are used to automaticaly determined a reasonable maxval1
582  Double_t prev2=0., MSB=0., AvB = 0. ; //those param are used to automaticaly determined a reasonable maxval2
583  Double_t prev3=0., MSki2=0., Avki2 = 0. ; //those param are used to automaticaly determined a reasonable maxval3
584  Double_t ki2=0.0 ;
585  for (Int_t k = 1; k <= ncells; k++) {
586  TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
587  TH1 *hCell = hCellAmplitude->ProjectionX("",k,k);
588  if (hCell->GetEntries() == 0) continue;
589  // hCell->Rebin(3);
590  hCell->Fit(fit, "0QEM", "", fitEmin, fitEmax);
591  delete hCell;
592 
593  if(fit->GetParameter(0) < 5000.){
594  hFitA->SetBinContent(k, fit->GetParameter(0));
595  if(k<3000) {
596  AvA += fit->GetParameter(0);
597  if(k==2999) maxval1 = AvA/3000. ;
598  if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
599  else MSA -= (fit->GetParameter(0) - prev) ;
600  prev = fit->GetParameter(0);
601  }
602 
603  else
604  {
605 
606  if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
607  {
608  maxval1 = fit->GetParameter(0);
609  }
610  }
611  }
612  else hFitA->SetBinContent(k, 5000.);
613 
614 
615 
616  if(fit->GetParameter(1) < 5000.){
617  hFitB->SetBinContent(k, fit->GetParameter(1));
618  if(k<3000) {
619  AvB += fit->GetParameter(1);
620  if(k==2999) maxval2 = AvB/3000. ;
621  if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
622  else MSB -= (fit->GetParameter(1) - prev2) ;
623  prev2 = fit->GetParameter(1);
624  }
625 
626  else
627  {
628 
629  if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
630  {
631  maxval2 = fit->GetParameter(1);
632  }
633  }
634  }
635  else hFitB->SetBinContent(k, 5000.);
636 
637 
638  if (fit->GetNDF() != 0 ) ki2 = fit->GetChisquare()/fit->GetNDF();
639  else ki2 = 1000.;
640 
641  if(ki2 < 1000.){
642  hFitChi2Ndf->SetBinContent(k, ki2);
643  if(k<3000) {
644  Avki2 += ki2;
645  if(k==2999) maxval3 = Avki2/3000. ;
646  if (prev3 < ki2) MSki2 += ki2 - prev3;
647  else MSki2 -= (ki2 - prev3) ;
648  prev3 = ki2;
649  }
650 
651  else
652  {
653 
654  if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
655  {
656  maxval3 = ki2;
657  }
658  }
659  }
660  else hFitChi2Ndf->SetBinContent(k, 1000.);
661 
662 
663  delete fit ;
664  }
665 
666  delete hCellAmplitude;
667 
668  // if you have problem with automatic parameter :
669  // maxval1 =
670  // maxval2 =
671  // maxval3 =
672 
673 
674  if(*pflag[0][0]==3)
675  Process(pflag, hFitChi2Ndf, Nsigma, dnbins, maxval3,compteur);
676 
677 
678  if(*pflag[0][0]==4)
679  Process(pflag, hFitA, Nsigma, dnbins, maxval1,compteur);
680 
681 
682  if(*pflag[0][0]==5)
683  Process(pflag, hFitB, Nsigma, dnbins, maxval2,compteur);
684 
685 }
686 
687 //_________________________________________________________________________
688 //_________________________________________________________________________
689 
690 void ExcludeCells(Int_t *pexclu[23040]) {
691  //find the cell with 0 entrie for excluding
692  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
693 
694 
695  for (Int_t c = 1; c <= 23040; c++) {
696  Double_t Nsum = 0;
697  // cout << "exclude cells ca va "<<endl
698  for (Int_t l = 1; l <= hCellAmplitude->GetNbinsX(); l++) {
699  Double_t N = hCellAmplitude->GetBinContent(l, c);
700  Nsum += N;
701  }
702  // cout << "2emem exclude cells ca va "<< c-1 << endl;
703  if(Nsum < 0.5 && *pexclu[c-1]!=1)
704  { *pexclu[c-1]=1; }//trick for criterum 7
705  //if(Nsum < 0.5 ) *pexclu[c-1]=1;
706  else *pexclu[c-1]=0;
707  }
708  delete hCellAmplitude;
709 }
710 
711 //_________________________________________________________________________
712 //_________________________________________________________________________
713 
714 void KillCells(Int_t filter[], Int_t nbc) {
715  // kill a cell : put it to 0 entrie
716  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
717 
718  for(Int_t i =0; i<nbc; i++){
719  for(Int_t j=0; j<= hCellAmplitude->GetNbinsX() ;j++){
720  hCellAmplitude->SetBinContent(j,filter[i]+1,0) ; }}
721 
722  TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
723 
724  TFile *tf = new TFile("filter.root","recreate");
725  hCellAmplitude->Write();
726  hNEventsProcessedPerRun->Write();
727  tf->Write();
728  tf->Close();
729  delete hCellAmplitude; delete hNEventsProcessedPerRun;
730 }
731 
732 //_________________________________________________________________________
733 //_________________________________________________________________________
734 
735 void PeriodAnalysis(Int_t criterum=7, Double_t Nsigma = 4.0, Double_t Emin=0.1, Double_t Emax=2.0, Int_t compteur = 1, TString datapath="/scratch/alicehp2/germain/QANew2",TString period = "LHC15f", TString pass = "pass2", Int_t trial=0, TString file ="none"){
736 
737  // what it does in function of criterum value
738 
739  // 1 : average E for E>Emin
740  // 2 : entries for E>Emin
741  // 3 : ki²/ndf (from fit of each cell Amplitude between Emin and Emax)
742  // 4 : A parameter (from fit of each cell Amplitude between Emin and Emax)
743  // 5 : B parameter (from fit of each cell Amplitude between Emin and Emax)
744  // 6 :
745  // 7 : give bad + dead list
746 
747  Int_t newBC[23040]; // newBC[0] donne l'id de la premiere BC trouvée
748  Int_t *pexclu[23040] ;
749  Int_t exclu[23040];
750  Int_t *pflag[23040][7] ;
751  Int_t flag[23040][7];
752  Int_t bad[10000] ;
753  Int_t i, j, nb=0;
754 
755  //INIT
756  TString output, bilan;
757  //bilan = Form("%s%sBC0Test%s.txt",period.Data());
758  if(criterum == 7) bilan = Form("%s%sBC0Test%i.txt",period.Data(),pass.Data(),trial);
759  // if(criterum == 7) bilan = "ResultsLHC15oBCRunList0Test1.txt" ;
760  output.Form("Criterum-%d_Emin-%.2f_Emax-%.2f.txt",criterum,Emin,Emax);
761  for(i=0;i<23040;i++) { exclu[i]=0; pexclu[i] =&exclu[i];
762  for(j=0;j<7;j++) { flag[i][j] =1 ; pflag[i][j] = &flag[i][j];}}
763  flag[0][0]=criterum ; //to identify the criterum tested
764 
765 
766  //CELLS EXCLUDED
767  ExcludeCells(pexclu); //exclude cells from analysis (will not appear in results)
768  // if(criterum < 7){
769  // cout<<"Excluded/dead cells : "<<endl;
770  // for(i=0;i<23040;i++) {if(exclu[i]!=0) {cout<<i<<", " ; nb++;}}
771  // cout<<"("<<nb<<")"<<endl; nb=0;}
772 
773 
774  //CRITERUM 7 : FINAL RESULT
775  if(criterum ==7) {
776  cout<<"FINAL RESULTS"<<endl;
777  ofstream fichier(bilan, ios::out | ios::trunc);
778  if(fichier){
779  fichier<<"Dead cells : "<<endl;
780  cout<<"Dead cells : "<<endl;
781  // MG modif
782  //for(i=0;i<23040;i++) {
783  for(i=0;i<17665;i++) {
784  //if(exclu[i]!=0) {fichier<<i<<", " ; cout<<i<<", " ; nb++;}}
785  if(exclu[i]!=0) {fichier<<i<<"\n" ; cout<<i<<", " ; nb++;}}
786  // fichier<<i<<"\n" ; cout<<i<<", " ; nb++;}
787  fichier<<"("<<nb<<")"<<endl; cout<<"("<<nb<<")"<<endl; nb=0;
788 
789  TFile::Open("filter.root");
790  ExcludeCells(pexclu);
791  fichier<<"Bad cells : "<<endl; cout<<"Bad cells : "<<endl;
792  for(i=0;i<17665;i++) {
793  // if(exclu[i]!=0) {bad[nb]=i; fichier<<i<<", " ; cout<<i<<", " ;
794  if(exclu[i]!=0) {bad[nb]=i; fichier<<i<<"\n" ; cout<<i<<", " ;
795  nb++;
796  //if(nb==999){ cout<<"TO MUCH BAD CELLS"<<endl ; break;}
797  }
798  }
799  fichier<<"("<<nb<<")"<<endl; cout<<"("<<nb<<")"<<endl;}
800  fichier.close();
801 
802  if(file!="none"){
803  TFile::Open(file);
804  Int_t w=0 ;
805  Int_t c;
806  for(w=0; (w*9)<=nb; w++) {
807  if(9<=(nb-w*9)) c = 9 ;
808  else c = nb-9*w ;
809  Draw(bad, w*9, c,datapath,period,pass,trial) ;
810  }}
811 
812  nb=0;
813  }
814 
815 
816  //ANALYSIS
817  if (criterum < 3) TestCellEandN(pflag, Emin, Emax,Nsigma,compteur);
818  else if (criterum < 6)
819  TestCellShapes(pflag, Emin, Emax, Nsigma,compteur);
820 
821 
822  //RESULTS
823  if(criterum < 6) { nb=0;
824  cout<<"bad by lower value Emin : "<< Emin << " Emax = " << Emax << endl;
825  for(i=0;i<17665;i++) {
826  if(flag[i][criterum]==0 && exclu[i]==0){nb++;
827  cout<<i<<", " ;}} cout<<"("<<nb<<")"<<endl; nb=0;
828 
829  cout<<"bad by higher value Emin : "<< Emin << " Emax = " << Emax <<endl;
830  for(i=0;i<17665;i++) {
831  if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
832  cout<<i<<", " ;}} cout<<"("<<nb<<")"<<endl; nb=0;
833 
834  cout<<"total bad "<<endl;
835  for(i=0;i<17665;i++) {
836  if(flag[i][criterum]!=1 && exclu[i]==0) {
837  newBC[nb]=i;
838  nb++;
839  cout<<i<<", " ; }} cout<<"("<<nb<<")"<<endl;
840 
841 
842  //create a filtered file
843  KillCells(newBC,nb) ; nb=0;
844 
845  //write in a file the results
846  ofstream fichier(output, ios::out | ios::trunc);
847  if(fichier)
848  {
849  fichier <<"criterum : "<<criterum<<", Emin = "<<Emin<<" GeV"<<", Emax = "<<Emax<<" GeV"<<endl;
850  fichier<<"bad by lower value : "<<endl;
851  for(i=0;i<17665;i++) {
852  if(flag[i][criterum]==0 && exclu[i]==0){nb++;
853  fichier<<i<<", " ;}} fichier<<"("<<nb<<")"<<endl; nb=0;
854 
855  fichier<<"bad by higher value : "<<endl;
856  for(i=0;i<17665;i++) {
857  if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
858  fichier<<i<<", " ;}} fichier<<"("<<nb<<")"<<endl; nb=0;
859  fichier<<"total bad "<<endl;
860  for(i=0;i<17665;i++) {
861  if(flag[i][criterum]!=1 && exclu[i]==0) {
862  newBC[nb]=i;
863  nb++;
864  fichier<<i<<", " ; }} fichier<<"("<<nb<<")"<<endl;
865  fichier.close();
866  }
867  else
868  cerr << "opening error" << endl;
869 
870 
871 
872 }
873 
874 }
875 //_________________________________________________________________________
876 //_________________________________________________________________________
877 
878  void BCAnalysis(TString file, TString datapath="scratch/alicehp2/germain/QANew2",TString trigger = "default",TString period = "LHC15f", TString pass = "pass2",Int_t trial = 0){
879 
880  //Configure a complete analysis with different criteria, it provides bad+dead cells lists
881  //You can manage criteria used and their order, the first criteria will use the original output file from AliAnalysisTaskCaloCellsQA task, then after each criteria it will use a filtered file without the badchannel previously identified
882 
883 
884 
885  Int_t criter;
886  Double_t Emini, Emaxi, Nsig;
887 
888 
889 
890 
891 // Default Configuration:
892 
893  //if(trigger=="default"){
894  if(trigger=="default"||trigger=="INT7"||trigger=="DMC7"||trigger=="AnyINTnoBC"){
895 
896 
897  TFile::Open(file);
898  PeriodAnalysis(2, 4., 0.2, 0.5,1,datapath,period,pass,trial); // nb ent emin emax
899  TFile::Open("filter.root");
900  PeriodAnalysis(2, 4., 0.5, 1.,1,datapath,period,pass,trial); // nb ent emin emax
901  TFile::Open("filter.root");
902  PeriodAnalysis(1, 6., 0.5, 1.,1,datapath,period,pass,trial); // energy mea emin emax
903  TFile::Open("filter.root");
904  PeriodAnalysis(2, 4., 1., 2.,1,datapath,period,pass,trial); // nb ent emin emax
905  TFile::Open("filter.root");
906  PeriodAnalysis(1, 6., 1., 2.,1,datapath,period,pass,trial); // energy mea emin emax
907  TFile::Open("filter.root");
908  PeriodAnalysis(2, 4., 1., 10.,1,datapath, period,pass,trial); //nb ent emin emax
909  // TFile::Open("filter.root");
910  // PeriodAnalysis(1, 6., 1., 10.,1,datapath,period,pass,trial); //energy mea emin emax
911 
912  TFile::Open("filter.root");
913  PeriodAnalysis(2, 4., 2., 3.,1,datapath,period,pass,trial); //entriesmea emin emax
914  TFile::Open("filter.root");
915  PeriodAnalysis(2, 4., 3., 4.,1,datapath,period,pass,trial); //entries mea emin emax
916  TFile::Open("filter.root");
917  PeriodAnalysis(2, 4., 4., 5.,1,datapath,period,pass,trial); //entriesmea emin emax
918  TFile::Open("filter.root");
919  PeriodAnalysis(2, 4., 5., 10.,1,datapath,period,pass,trial); //entries mea emin emax
920 
921 // TFile::Open("filter.root");
922 // PeriodAnalysis(4, 4., 0.1, 2.,1,period,pass,trial); // param A fit E = A exp(-Bx)/x^2
923 
924  }
925 
926 
927  else { //you have the possibility to change analysis configuration in function of trigger type
928 
929  TFile::Open(file);
930  PeriodAnalysis(2, 6., 0.5, 2.,1,datapath,period,pass,trial); // nb ent emin emax
931  TFile::Open("filter.root");
932  PeriodAnalysis(1, 6., 0.5, 2.,1,datapath,period,pass,trial); // energy mea emin emax
933  TFile::Open("filter.root");
934  PeriodAnalysis(2, 6., 2., 5.,1,datapath,period,pass,trial); // nb ent emin emax
935  TFile::Open("filter.root");
936  PeriodAnalysis(1, 6., 2., 5.,1,datapath,period,pass,trial); // energy mea emin emax
937  TFile::Open("filter.root");
938  PeriodAnalysis(2, 6., 5., 10.,1,datapath,period,pass,trial); // nb ent emin emax
939  TFile::Open("filter.root");
940  PeriodAnalysis(1, 6., 5., 10.,1,datapath,period,pass,trial); // energy mea emin emax
941  TFile::Open("filter.root");
942  PeriodAnalysis(1, 6., 5., 10.,1,datapath,period,pass,trial); // energy mea emin emax
943 
944  // this was old settings for EMC riggers checks
945 // TFile::Open(file);
946 // // PeriodAnalysis(3, 8., 1., 3.,1,period,pass,trial);
947 
948 
949  }
950 
951  TFile::Open(file);
952  PeriodAnalysis(7,0.,0.,0.,1,datapath,period,pass,trial,file); //provide dead cells list from original file and draw bad cells candidate from indicated file
953 
954 }
955 
956 //_________________________________________________________________________
957 //________________________________________________________________________
958 
959 void BadChannelAnalysis(TString datapath= "/scratch/alicehp2,germain/QANew2",TString fCalorimeter = "EMCAL", TString period = "LHC15f", TString pass = "pass2", TString trigger= "default",Int_t trial=0){
960  //Convert(datapath,fCalorimeter, period, pass, trigger);
961  // TString inputfile(Form( "/scratch/alicehp2/germain/QANew2/%s/%s/%s%sRunlist0New.root",period.Data(),pass.Data(),period.Data(),pass.Data(),trigger.Data()));
962 
963  TString inputfile(Form( "%s/%s/%s/%s%sRunlist1New.root",datapath.Data(),period.Data(),pass.Data(),period.Data(),pass.Data(),trigger.Data()));
964 
965  BCAnalysis(inputfile,datapath,trigger,period,pass,trial);
966 }
void TestCellShapes(Int_t *pflag[23040][7], Double_t fitEmin, Double_t fitEmax, Double_t Nsigma=4., Int_t compteur=1, char const *hname="hCellAmplitude", Int_t dnbins=1000)
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
void Draw(Int_t cell[], Int_t iBC, Int_t nBC, TString datapath="/scratch/alicehp2/germain/QANew2/", TString period="LHC15f", TString pass="pass2", Int_t trial=0, const Int_t cellref=2377)
TCanvas * c
Definition: TestFitELoss.C:172
int Int_t
Definition: External.C:63
void BCAnalysis(TString file, TString datapath="scratch/alicehp2/germain/QANew2", TString trigger="default", TString period="LHC15f", TString pass="pass2", Int_t trial=0)
void PeriodAnalysis(Int_t criterum=7, Double_t Nsigma=4.0, Double_t Emin=0.1, Double_t Emax=2.0, Int_t compteur=1, TString datapath="/scratch/alicehp2/germain/QANew2", TString period="LHC15f", TString pass="pass2", Int_t trial=0, TString file="none")
void BadChannelAnalysis(TString datapath="/scratch/alicehp2,germain/QANew2", TString fCalorimeter="EMCAL", TString period="LHC15f", TString pass="pass2", TString trigger="default", Int_t trial=0)
Definition: External.C:212
void Convert(TString datapath="/scratch/alicehp2/germain/QANew2", TString fCalorimeter="EMCAL", TString period="LHC11h", TString pass="pass1_HLT", TString trigger="default")
void KillCells(Int_t filter[], Int_t nbc)
Definition: External.C:220
void ExcludeCells(Int_t *pexclu[23040])
TFile * file
void Process(Int_t *pflag[23040][7], TH1 *inhisto, Double_t Nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1., Int_t compteur=1)
void Draw2(Int_t cell, Int_t cellref=400)
Definition: External.C:196
void TestCellEandN(Int_t *pflag[23040][7], Double_t Emin=0.1, Double_t Emax=2., Double_t Nsigma=4., Int_t compteur=1, char const *hname="hCellAmplitude", Int_t dnbins=200)