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