PWGGA/CaloTrackCorrelations/macros/QA/BadChannelAnalysis.C

   1 // This macro has been developed to find badcells candidates in EMCal based on cells amplitude distributions
   2 //
   3 // Input needed :
   4 // can be either outputs QA from AliAnaCalorimeterQA task (BadChannelAnalysis() function)
   5 // Or from merged output of AliAnalysisTaskCaloCellsQA (use BCAnalysis() function)
   6 //
   7 // Ouput:
   8 // in Results.txt file and cell control plots for bad channel candidates will be created in BadCandidate.pdf
   9 //
  10 // To classify between bad/warm/good (this  last case can also occur in bad candidates)
  11 // the user will have to check manually (by eye)control plots. In those in red the bad cell candidate energy distrib in black a reference choosen
  12 // just to guide the eye
  13 //
  14 //
  15 // Author : Alexis Mas (SUBATECH), based on getCellsRunQA.C from Olga Driga (SUBATECH)
  16
  17 #if !defined(__CINT__) || defined(__MAKECINT__)
  18 #include <TFile.h>
  19 #include <TH1F.h>
  20 #include <TH2F.h>
  21 #include <TH3D.h>
  22 #include <TLine.h>
  23 #include <Riostream.h>
  24 #include <TCanvas.h>
  25 #include <TGraphErrors.h>
  26 #include <TGrid.h>
  27 #include <TStyle.h>
  28 #include <TFileMerger.h>
  29 #include <TMultiGraph.h>
  30 #include <TROOT.h>
  31 #include <TLegend.h>
  32 #include <TString.h>
  33 #include <TGridCollection.h>
  34 #include <TGridResult.h>
  35 #include <TClonesArray.h>
  36 #include <TObjString.h>
  37 #include <stdio.h>
  38 #include <fstream>
  39 #include <iostream>
  40 #endif
  41 using namespace std;
  42
  43
  44 void Draw(Int_t cell[], Int_t iBC, Int_t nBC, const Int_t cellref=151){
  45   //Allow to produce a pdf file with badcells candidates (red) compared to a refence cell (black)
  46
  47   gROOT ->SetStyle("Plain");
  48   gStyle->SetOptStat(0);
  49   gStyle->SetFillColor(kWhite);
  50   gStyle->SetTitleFillColor(kWhite);
  51   gStyle->SetPalette(1);
  52
  53   char out[120]; char title[100]; char name[100];
  54   TString slide = "Cells ";
  55   slide+=cell[iBC];
  56   slide+="-";
  57   slide+=cell[iBC+nBC-1];
  58   sprintf(out,"%d-%d.gif",cell[iBC],cell[iBC+nBC-1]);
  59
  60   TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
  61   TH1 *hCellref       = hCellAmplitude->ProjectionY("badcells",cellref+1,cellref+1);
  62
  63   Int_t i;
  64   TCanvas *c1 = new TCanvas("badcells","badcells",1000,750) ;
  65   c1->SetLogy();
  66   if      (nBC > 6) c1->Divide(3,3) ;
  67   else if (nBC > 3) c1->Divide(3,2) ;
  68   else              c1->Divide(3,1);
  69   // hCellref->Rebin(3);
  70
  71   TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
  72   for(i=0; i<nBC ; i++){
  73     sprintf(name,"Cell %d",cell[iBC+i]) ;
  74     TH1 *hCell = hCellAmplitude->ProjectionY(name,cell[iBC+i]+1,cell[iBC+i]+1);
  75     c1->cd(i%9 + 1) ;
  76     c1->cd(i%9 + 1)->SetLogy();
  77     sprintf(title,"Cell %d      Entries : %d",cell[iBC+i], (Int_t)hCell->GetEntries()) ;
  78     hCell->SetLineColor(2)  ;
  79     // cout<<title<<endl ;
  80     hCell->SetMaximum(1e5);
  81     // hCell->Rebin(3);
  82     hCell->SetAxisRange(0.,4.);
  83     hCell->GetXaxis()->SetTitle("E (GeV)");
  84     hCell->GetYaxis()->SetTitle("N Entries");
  85     hCellref->SetAxisRange(0.,4.);
  86     hCell->SetLineWidth(1) ;
  87     hCellref->SetLineWidth(1) ;
  88     hCell->SetTitle(title);
  89     hCellref->SetLineColor(1)  ;
  90     if(i==0){
  91       leg->AddEntry(hCellref,"reference","l");
  92       leg->AddEntry(hCell,"current","l");;
  93     }
  94     hCell->Draw() ;
  95     hCellref->Draw("same") ;
  96     leg->Draw();
  97   }
  98
  99   //CREATE A PDF FILE
 100   TString PdfName = "BadCellsCandidate.pdf";
 101   if(nBC<9) {
 102     PdfName +=")";
 103     c1->Print(PdfName.Data());}
 104   else if(iBC==0) {
 105     PdfName +="(";
 106     c1->Print(PdfName.Data());}
 107   else  c1->Print(PdfName.Data());
 108
 109
 110   //  c1->SaveAs(out);
 111   delete hCellref ;
 112   delete c1 ;
 113   delete leg ;
 114
 115 }
 116
 117 //_________________________________________________________________________
 118 //_________________________________________________________________________
 119
 120 void Convert(TString fCalorimeter = "EMCAL", TString period = "LHC11h", TString pass = "pass1_HLT", TString trigger= "default"){
 121
 122   // Create one file for the analysis from several outputs QA files listed in runlist.txt
 123   // You need :
 124   //   runlist.txt with runs listed
 125   //   outputsQA  e.g  period/pass/123456.root
 126
 127   TH2F * hCellAmplitude = new TH2F("hCellAmplitude","Cell Amplitude",11520,0,11520,200,0,10);
 128   TH1D * hNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun","Number of processed events vs run number",200000,100000,300000);
 129   FILE * pFile;
 130   TString file = "/scratch/alicehp2/mas/analyse/QA/"+period+"/"+ pass + "/runlistMB.txt" ;
 131   cout<<file<<endl;
 132   pFile = fopen(file.Data(), "r"); //open the text file where include the run list and correct run index
 133
 134   cout<<file<<endl;
 135   cout << " fcalo: " << fCalorimeter << "; period: " << period << "; pass: " << pass << "  trigger "<<trigger<<  endl;
 136
 137   Int_t ix,iy;
 138   Int_t Nentr;
 139   Int_t p;
 140   Int_t q;
 141   Int_t ncols;
 142   Int_t nlines = 0 ;
 143   Int_t RunId[500] ;
 144
 145   Double_t x[500] ;
 146   Double_t xrun[500] ;
 147
 148   while (1){
 149     ncols = fscanf(pFile,"%d  %d ",&p,&q);
 150     if (ncols< 0) break;
 151     x[nlines]=p;
 152     RunId[nlines]=q;
 153     xrun[nlines]=1.*q;
 154     nlines++;
 155   }
 156
 157   fclose(pFile);
 158
 159   const Int_t nRun = nlines ;
 160   Double_t content;
 161   TString base ;
 162   TString BCfile ;
 163   TString direct = "CaloQA_";
 164   direct += trigger;
 165
 166   for(Int_t i = 0 ; i < nRun ; i++) {
 167     base = "/scratch/alicehp2/germain/QA/";
 168     BCfile = base + period ;
 169     BCfile += trigger ;
 170     BCfile += ".root";
 171     base += period ;
 172     base += "/";
 173     base += pass ;
 174     base += "/";
 175     base += RunId[i] ;
 176     TString infile ;
 177     infile = base + ".root" ;
 178     TFile *f = TFile::Open(infile);
 179     base += "/" ;
 180     base += trigger ;
 181     TDirectoryFile *dir = (TDirectoryFile *)f->Get(direct);
 182     TList *outputList = (TList*)dir->Get(direct);
 183     TH2F *hAmpId;
 184     TH2F *hNEvents;
 185     hAmpId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
 186     hNEvents =(TH2F *)outputList->FindObject("hNEvents");
 187     Nentr =  (Int_t)hNEvents->GetEntries();
 188     if (Nentr<100) continue ;
 189     hNEventsProcessedPerRun->SetBinContent(RunId[i]-100000,(Double_t)Nentr);
 190
 191     for(ix=1;ix<=200;ix++){
 192       for(iy=1; iy<=11520; iy++){
 193         content = 0.0 ;
 194         content +=  hAmpId->GetBinContent(ix,iy);
 195         content +=  hCellAmplitude->GetBinContent(iy,ix);
 196         if(content > 0.5)
 197           hCellAmplitude->SetBinContent(iy,ix,content);
 198       }
 199     }
 200
 201     //cout<<i<<endl;
 202     if(i==0){ cout<<"Merging/Converting procedure ..." ; cout.flush();}
 203     else { cout<<"..." ; cout.flush();}
 204     outputList->Delete();
 205     dir->Delete();
 206     f->Close();
 207     delete f;
 208   }
 209
 210
 211   TFile *BCF = TFile::Open(BCfile,"recreate");
 212   hNEventsProcessedPerRun->Write();
 213   hCellAmplitude->Write();
 214   BCF->Close();
 215   cout<<"DONE !"<<endl;
 216
 217 }
 218
 219 //_________________________________________________________________________
 220 //_________________________________________________________________________
 221
 222 void Process(Int_t *pflag[11520][7], TH1* inhisto, Double_t Nsigma = 4., Int_t dnbins = 200, Double_t dmaxval = -1.)
 223 {
 224   //  1) create a distribution for the input histogram;
 225   //  2) fit the distribution with a gaussian
 226   //  3) define good area within +-Nsigma to identfy badcells
 227   //
 228   // inhisto -- input histogram;
 229   // dnbins -- number of bins in distribution;
 230   // dmaxval -- maximum value on distribution histogram.
 231
 232   Int_t crit = *pflag[0][0] ; //identify the criterum processed
 233   Double_t goodmax= 0. ;
 234   Double_t goodmin= 0. ;
 235   *pflag[0][0] =1;
 236   if (dmaxval < 0.) {
 237     dmaxval = inhisto->GetMaximum()*1.01;  // 1.01 - to see the last bin
 238     if(crit==2 && dmaxval > 1) dmaxval =1. ;
 239   }
 240
 241   TH1 *distrib = new TH1F(Form("%sDistr",inhisto->GetName()), "", dnbins, inhisto->GetMinimum(), dmaxval);
 242   distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
 243   distrib->SetYTitle("Entries");
 244
 245   // fill distribution
 246   for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
 247     distrib->Fill(inhisto->GetBinContent(c));
 248
 249   // draw histogram + distribution
 250   TCanvas *c1 = new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400);
 251   c1->Divide(2,1);
 252
 253   c1->cd(1);
 254   gPad->SetLeftMargin(0.14);
 255   gPad->SetRightMargin(0.06);
 256   gPad->SetLogy();
 257   inhisto->SetTitleOffset(1.7,"Y");
 258   inhisto->Draw();
 259
 260   c1->cd(2);
 261   gPad->SetLeftMargin(0.14);
 262   gPad->SetRightMargin(0.10);
 263   gPad->SetLogy();
 264   distrib->Draw();
 265
 266   Int_t higherbin=0,i;
 267   for (i = 2; i <= distrib->GetNbinsX(); i++){
 268     if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i))
 269       higherbin = i ;
 270   }
 271
 272   for(i=higherbin ; i<=dnbins ; i++)
 273     if(distrib->GetBinContent(i)<2) break ;
 274   goodmax = distrib->GetBinCenter(i);
 275
 276   for(i=higherbin ; i>0 ; i--)
 277     if(distrib->GetBinContent(i)<2) break ;
 278   goodmin = distrib->GetBinLowEdge(i);
 279
 280   TF1 *fit2 = new TF1("fit2", "gaus");
 281
 282   distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
 283   Double_t sig, mean, chi2ndf;
 284   mean = fit2->GetParameter(1);
 285   sig = fit2->GetParameter(2);
 286   chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
 287   goodmin = mean - Nsigma*sig ;
 288   goodmax = mean + Nsigma*sig ;
 289
 290   // lines
 291   TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
 292   lline->SetLineColor(kOrange);
 293   lline->SetLineStyle(7);
 294   lline->Draw();
 295
 296   TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
 297   rline->SetLineColor(kOrange);
 298   rline->SetLineStyle(7);
 299   rline->Draw();
 300
 301   // legend
 302   TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
 303   leg->AddEntry(lline, "Good region boundary","l");
 304   leg->Draw("same");
 305   fit2->Draw("same");
 306
 307   c1->Update();
 308   TString name = "criteria-" ;
 309   name+= crit;
 310   name+= ".gif";
 311
 312   c1->SaveAs(name);
 313
 314   Int_t ntot = 0, cel;
 315
 316   for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++) {
 317     cel=(Int_t)(inhisto->GetBinLowEdge(c)+0.1);
 318     if (inhisto->GetBinContent(c) < goodmin) {
 319       ntot++;
 320       *pflag[cel][crit]=0;
 321     }
 322     else if (inhisto->GetBinContent(c) > goodmax) {
 323       ntot++;
 324       *pflag[cel][crit]=2;
 325     }
 326   }
 327 }
 328
 329 //_________________________________________________________________________
 330 //_________________________________________________________________________
 331
 332 void TestCellEandN(Int_t *pflag[11520][7], Double_t Emin = 0.1, Double_t Nsigma = 4., char* hname = "hCellAmplitude", Int_t dnbins = 200)
 333 {
 334
 335
 336   // Three more tests for bad cells:
 337   //  1) total deposited energy;
 338   //  2) total number of entries;
 339   //  3) average energy = [total deposited energy]/[total number of entries].
 340   //
 341
 342
 343   // input; X axis -- absId numbers
 344   TH2 *hCellAmplitude = (TH2*) gFile->Get(hname);
 345
 346   // binning parameters
 347   Int_t ncells = hCellAmplitude->GetNbinsX();
 348   Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
 349   Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
 350
 351   TH1* hCellEtotal = new TH1F(Form("%s_hCellEtotal_E%.2f",hname,Emin),
 352                               Form("Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax);
 353   hCellEtotal->SetXTitle("AbsId");
 354   hCellEtotal->SetYTitle("Energy, GeV");
 355
 356   TH1F *hCellNtotal = new TH1F(Form("%s_hCellNtotal_E%.2f",hname,Emin),
 357                                Form("Number of entries per events, E > %.2f GeV",Emin), ncells,amin,amax);
 358   hCellNtotal->SetXTitle("AbsId");
 359   hCellNtotal->SetYTitle("Entries");
 360
 361   TH1F *hCellEtoNtotal = new TH1F(Form("%s_hCellEtoNtotal_E%.2f",hname,Emin),
 362                                   Form("Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax);
 363   hCellEtoNtotal->SetXTitle("AbsId");
 364   hCellEtoNtotal->SetYTitle("Energy, GeV");
 365
 366   TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
 367   Double_t totalevents = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
 368
 369   // fill cells
 370   for (Int_t c = 1; c <= ncells; c++) {
 371     Double_t Esum = 0;
 372     Double_t Nsum = 0;
 373
 374
 375     for (Int_t j = 1; j <= hCellAmplitude->GetNbinsY(); j++) {
 376       Double_t E = hCellAmplitude->GetYaxis()->GetBinCenter(j);
 377       Double_t N = hCellAmplitude->GetBinContent(c, j);
 378       if (E < Emin) continue;
 379       Esum += E*N;
 380       Nsum += N;
 381     }
 382
 383     hCellEtotal->SetBinContent(c, Esum);
 384     hCellNtotal->SetBinContent(c, Nsum/totalevents);
 385
 386     if (Nsum > 0.5)  // number of entries >= 1
 387       hCellEtoNtotal->SetBinContent(c, Esum/Nsum);
 388
 389   }
 390
 391   delete hCellAmplitude;
 392
 393   // Process(hCellEtotal,   dnbins );
 394   if(*pflag[0][0]==1)
 395     Process(pflag, hCellEtoNtotal, Nsigma, dnbins );
 396   if(*pflag[0][0]==2)
 397     Process(pflag, hCellNtotal, Nsigma,  dnbins );
 398 }
 399
 400 //_________________________________________________________________________
 401 //_________________________________________________________________________
 402
 403 void TestCellShapes(Int_t *pflag[11520][7], Double_t fitEmin, Double_t fitEmax, Double_t Nsigma =4., char* hname = "hCellAmplitude", Int_t dnbins = 1000)
 404 {
 405   // Test cells shape using fit function f(x)=A*exp(-B*x)/x^2.
 406   // Produce values per cell + distributions for A,B and chi2/ndf parameters.
 407
 408   TH2 *hCellAmplitude = (TH2*) gFile->Get(Form("%s",hname));
 409
 410   // binning parameters
 411   Int_t  ncells = hCellAmplitude->GetNbinsX();
 412   Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
 413   Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
 414
 415   // initialize histograms
 416   TH1 *hFitA = new TH1F(Form("hFitA_%s",hname),"Fit A value", ncells,amin,amax);
 417   hFitA->SetXTitle("AbsId");
 418   hFitA->SetYTitle("A");
 419
 420   TH1 *hFitB = new TH1F(Form("hFitB_%s",hname),"Fit B value", ncells,amin,amax);
 421   hFitB->SetXTitle("AbsId");
 422   hFitB->SetYTitle("B");
 423
 424   TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",hname),"Fit #chi^{2}/ndf value", ncells,amin,amax);
 425   hFitChi2Ndf->SetXTitle("AbsId");
 426   hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
 427
 428   Double_t maxval1=0., maxval2=0., maxval3=0.;
 429   Double_t prev=0., MSA=0., AvA = 0. ; //those param are used to automaticaly determined a reasonable maxval1
 430   Double_t prev2=0., MSB=0., AvB = 0.  ; //those param are used to automaticaly determined a reasonable maxval2
 431   Double_t prev3=0., MSki2=0., Avki2 = 0. ; //those param are used to automaticaly determined a reasonable maxval3
 432   Double_t ki2=0.0 ;
 433
 434   for (Int_t k = 1; k <= ncells; k++) {
 435     TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
 436     TH1 *hCell = hCellAmplitude->ProjectionY("",k,k);
 437     if (hCell->GetEntries() == 0) continue;
 438     // hCell->Rebin(3);
 439     hCell->Fit(fit, "0QEM", "", fitEmin, fitEmax);
 440     delete hCell;
 441
 442     if(fit->GetParameter(0) < 5000.){
 443       hFitA->SetBinContent(k, fit->GetParameter(0));
 444       if(k<3000) {
 445         AvA +=  fit->GetParameter(0);
 446         if(k==2999)  maxval1  = AvA/3000. ;
 447         if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
 448         else MSA -= (fit->GetParameter(0) - prev) ;
 449         prev = fit->GetParameter(0);
 450       }
 451
 452       else
 453       {
 454
 455         if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
 456         {
 457           maxval1 = fit->GetParameter(0);
 458         }
 459       }
 460     }
 461     else hFitA->SetBinContent(k, 5000.);
 462
 463
 464
 465     if(fit->GetParameter(1) < 5000.){
 466       hFitB->SetBinContent(k, fit->GetParameter(1));
 467       if(k<3000) {
 468         AvB +=  fit->GetParameter(1);
 469         if(k==2999)  maxval2  = AvB/3000. ;
 470         if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
 471         else MSB -= (fit->GetParameter(1) - prev2) ;
 472         prev2 = fit->GetParameter(1);
 473       }
 474
 475       else
 476       {
 477
 478         if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
 479         {
 480           maxval2 = fit->GetParameter(1);
 481         }
 482       }
 483     }
 484     else hFitB->SetBinContent(k, 5000.);
 485
 486
 487     if (fit->GetNDF() != 0 ) ki2 =  fit->GetChisquare()/fit->GetNDF();
 488     else ki2 = 1000.;
 489
 490     if(ki2 < 1000.){
 491       hFitChi2Ndf->SetBinContent(k, ki2);
 492       if(k<3000) {
 493         Avki2 +=  ki2;
 494         if(k==2999)  maxval3  = Avki2/3000. ;
 495         if (prev3 < ki2) MSki2 += ki2 - prev3;
 496         else MSki2 -= (ki2 - prev3) ;
 497         prev3 = ki2;
 498       }
 499
 500       else
 501       {
 502
 503         if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
 504         {
 505           maxval3 = ki2;
 506         }
 507       }
 508     }
 509     else hFitChi2Ndf->SetBinContent(k, 1000.);
 510
 511
 512     delete fit ;
 513   }
 514
 515   delete hCellAmplitude;
 516
 517   // if you have problem with automatic parameter :
 518   //  maxval1 =
 519   //  maxval2 =
 520   //  maxval3 =
 521
 522
 523   if(*pflag[0][0]==3)
 524     Process(pflag, hFitChi2Ndf, Nsigma, dnbins, maxval3);
 525
 526
 527   if(*pflag[0][0]==4)
 528     Process(pflag, hFitA, Nsigma, dnbins,  maxval1);
 529
 530
 531   if(*pflag[0][0]==5)
 532     Process(pflag, hFitB, Nsigma, dnbins, maxval2);
 533
 534 }
 535
 536 //_________________________________________________________________________
 537 //_________________________________________________________________________
 538
 539 void ExcludeCells(Int_t *pexclu[11520]) {
 540   //find the cell with 0 entrie for excluding
 541   TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
 542
 543
 544   for (Int_t c = 1; c <= 11520; c++) {
 545     Double_t Nsum = 0;
 546
 547     for (Int_t l = 1; l <= hCellAmplitude->GetNbinsY(); l++) {
 548       Double_t N = hCellAmplitude->GetBinContent(c, l);
 549       Nsum += N;
 550     }
 551     if(Nsum < 0.5 && *pexclu[c-1]!=1) *pexclu[c-1]=1; //trick for criterum 7
 552     //if(Nsum < 0.5 ) *pexclu[c-1]=1;
 553     else *pexclu[c-1]=0;
 554   }
 555   delete hCellAmplitude;
 556 }
 557
 558 //_________________________________________________________________________
 559 //_________________________________________________________________________
 560
 561 void KillCells(Int_t filter[], Int_t nbc) {
 562   // kill a cell : put it to 0 entrie
 563   TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
 564
 565   for(Int_t i =0; i<nbc; i++){
 566     for(Int_t j=0; j<= hCellAmplitude->GetNbinsY() ;j++){
 567       hCellAmplitude->SetBinContent(filter[i]+1,j,0) ; }}
 568
 569   TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
 570
 571   TFile *tf = new TFile("filter.root","recreate");
 572   hCellAmplitude->Write();
 573   hNEventsProcessedPerRun->Write();
 574   tf->Write();
 575   tf->Close();
 576   delete hCellAmplitude; delete hNEventsProcessedPerRun;
 577 }
 578
 579 //_________________________________________________________________________
 580 //_________________________________________________________________________
 581
 582 void PeriodAnalysis(Int_t criterum=7, Double_t Nsigma = 4.0, Double_t Emin=0.1, Double_t Emax=1.0, TString file ="none") {
 583
 584   // what it does in function of criterum value
 585
 586   // 1 : average E for E>Emin
 587   // 2 : entries for E>Emin
 588   // 3 : ki²/ndf  (from fit of each cell Amplitude between Emin and Emax)
 589   // 4 : A parameter (from fit of each cell Amplitude between Emin and Emax)
 590   // 5 : B parameter (from fit of each cell Amplitude between Emin and Emax)
 591   // 6 :
 592   // 7 : give bad + dead list
 593
 594   Int_t newBC[11520]; // newBC[0] donne l'id de la premiere BC trouvée
 595   Int_t *pexclu[11520] ;
 596   Int_t exclu[11520];
 597   Int_t *pflag[11520][7] ;
 598   Int_t flag[11520][7];
 599   Int_t bad[1000] ;
 600   Int_t i, j, nb=0;
 601
 602   //INIT
 603   TString output, bilan;
 604   if(criterum == 7) bilan = "Results.txt" ;
 605   output.Form("Criterum-%d_Emin-%.2f.txt",criterum,Emin);
 606   for(i=0;i<11520;i++) { exclu[i]=0; pexclu[i] =&exclu[i];
 607     for(j=0;j<7;j++) { flag[i][j] =1 ; pflag[i][j] = &flag[i][j];}}
 608   flag[0][0]=criterum ; //to identify the criterum tested
 609
 610
 611   //CELLS EXCLUDED
 612   ExcludeCells(pexclu); //exclude cells from analysis (will not appear in results)
 613   if(criterum < 7){
 614     cout<<"Excluded/dead cells : "<<endl;
 615     for(i=0;i<11520;i++) {if(exclu[i]!=0) {cout<<i<<", " ; nb++;}}
 616     cout<<"("<<nb<<")"<<endl; nb=0;}
 617
 618
 619   //CRITERUM 7 : FINAL RESULT
 620   if(criterum ==7) {
 621     cout<<"FINAL RESULTS"<<endl;
 622     ofstream fichier(bilan, ios::out | ios::trunc);
 623     if(fichier){
 624       fichier<<"Dead cells : "<<endl;
 625       cout<<"Dead cells : "<<endl;
 626       for(i=0;i<11520;i++) {
 627         if(exclu[i]!=0) {fichier<<i<<", " ; cout<<i<<", " ; nb++;}}
 628       fichier<<"("<<nb<<")"<<endl; cout<<"("<<nb<<")"<<endl; nb=0;
 629
 630       TFile::Open("filter.root");
 631       ExcludeCells(pexclu);
 632       fichier<<"Bad cells candidates : "<<endl; cout<<"Bad cells candidates : "<<endl;
 633       for(i=0;i<11520;i++) {
 634         if(exclu[i]!=0) {bad[nb]=i; fichier<<i<<", " ; cout<<i<<", " ;
 635           nb++; if(nb==999){ cout<<"TO MUCH BAD CELLS"<<endl ; break;}}}
 636       fichier<<"("<<nb<<")"<<endl; cout<<"("<<nb<<")"<<endl;}
 637     fichier.close();
 638
 639     if(file!="none"){
 640       TFile::Open(file);
 641       Int_t w=0 ;
 642       Int_t c;
 643       for(w=0; (w*9)<=nb; w++) {
 644         if(9<=(nb-w*9)) c = 9 ;
 645         else c = nb-9*w ;
 646         Draw(bad, w*9, c) ;
 647       }}
 648
 649     nb=0;
 650   }
 651
 652
 653   //ANALYSIS
 654   if (criterum < 3)  TestCellEandN(pflag, Emin, Nsigma);
 655   else if (criterum < 6)
 656     TestCellShapes(pflag, Emin, Emax, Nsigma);
 657
 658
 659   //RESULTS
 660   if(criterum < 6) { nb=0;
 661     cout<<"bad by lower value : "<<endl;
 662     for(i=0;i<11520;i++) {
 663       if(flag[i][criterum]==0 && exclu[i]==0){nb++;
 664         cout<<i<<", " ;}} cout<<"("<<nb<<")"<<endl; nb=0;
 665
 666     cout<<"bad by higher value : "<<endl;
 667     for(i=0;i<11520;i++) {
 668       if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
 669         cout<<i<<", " ;}} cout<<"("<<nb<<")"<<endl; nb=0;
 670
 671     cout<<"total bad "<<endl;
 672     for(i=0;i<11520;i++) {
 673       if(flag[i][criterum]!=1 && exclu[i]==0) {
 674         newBC[nb]=i;
 675         nb++;
 676         cout<<i<<", " ; }} cout<<"("<<nb<<")"<<endl;
 677
 678
 679     //create a filtered file
 680     KillCells(newBC,nb) ; nb=0;
 681
 682     //write in a file the results
 683     ofstream fichier(output, ios::out | ios::trunc);
 684     if(fichier)
 685     {
 686       fichier <<"criterum : "<<criterum<<", Emin = "<<Emin<<" GeV"<<", Emax = "<<Emax<<" GeV"<<endl;
 687       fichier<<"bad by lower value : "<<endl;
 688       for(i=0;i<11520;i++) {
 689         if(flag[i][criterum]==0 && exclu[i]==0){nb++;
 690           fichier<<i<<", " ;}} fichier<<"("<<nb<<")"<<endl; nb=0;
 691
 692       fichier<<"bad by higher value : "<<endl;
 693       for(i=0;i<11520;i++) {
 694         if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
 695           fichier<<i<<", " ;}} fichier<<"("<<nb<<")"<<endl; nb=0;
 696
 697       fichier<<"total bad "<<endl;
 698       for(i=0;i<11520;i++) {
 699         if(flag[i][criterum]!=1 && exclu[i]==0) {
 700           newBC[nb]=i;
 701           nb++;
 702           fichier<<i<<", " ; }} fichier<<"("<<nb<<")"<<endl;
 703       fichier.close();
 704     }
 705     else
 706       cerr << "opening error" << endl;
 707
 708   }
 709
 710 }
 711
 712
 713 //_________________________________________________________________________
 714 //_________________________________________________________________________
 715
 716 void BCAnalysis(TString file, TString trigger = "default"){
 717
 718   //Configure a complete analysis with different criteria, it provides bad+dead cells lists
 719   //You can manage criteria used and their order, the first criteria will use the original
 720   //output file from AliAnalysisTaskCaloCellsQA task, then after each criteria it will use a filtered file without the badchannel previously identified
 721
 722   if(trigger=="default"){
 723
 724     TFile::Open(file);
 725     PeriodAnalysis(3, 8., 0.1, 2.5);
 726     TFile::Open("filter.root");
 727     PeriodAnalysis(2, 4., 0.1, 2.5);
 728     TFile::Open("filter.root");
 729     PeriodAnalysis(2, 4., 0.5, 2.5);
 730     TFile::Open("filter.root");
 731     PeriodAnalysis(1, 4., 0.1, 2.5);
 732     TFile::Open("filter.root");
 733     PeriodAnalysis(4, 4., 0.1,2.5);
 734
 735   }
 736
 737   else { //you have the possibility to change analysis configuration  in function of trigger type
 738
 739     TFile::Open(file);
 740     PeriodAnalysis(3, 8., 0.1, 2.5);
 741     TFile::Open("filter.root");
 742     PeriodAnalysis(2, 4., 0.1, 2.5);
 743     TFile::Open("filter.root");
 744     PeriodAnalysis(2, 4., 0.5, 2.5);
 745     TFile::Open("filter.root");
 746     PeriodAnalysis(1, 4., 0.1, 2.5);
 747     TFile::Open("filter.root");
 748     PeriodAnalysis(4, 4., 0.1, 2.5);
 749
 750   }
 751
 752   TFile::Open(file);
 753   PeriodAnalysis(7,0.,0.,0.,file); //provide dead cells list from original file and draw bad cells candidate from indicated file
 754
 755 }
 756
 757 //_________________________________________________________________________
 758 //_________________________________________________________________________
 759
 760 void BadChannelAnalysis(TString fCalorimeter = "EMCAL", TString period = "LHC11h", TString pass = "pass1_HLT", TString trigger= "default"){
 761   Convert(fCalorimeter, period, pass, trigger);
 762   TString inputfile =  "/scratch/alicehp2/germain/QA/" + period;
 763   inputfile += trigger;
 764   inputfile += ".root";
 765   BCAnalysis(inputfile, trigger);
 766 }