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