AliPhysics  2797316 (2797316)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaCaloChannelAnalysis.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TFile.h>
18 #include <TH1F.h>
19 #include <TH2F.h>
20 #include <TF1.h>
21 #include <TSpectrum.h>
22 #include <TROOT.h>
23 #include <TLine.h>
24 #include <TCanvas.h>
25 #include <TStyle.h>
26 #include <TSystem.h>
27 #include <TLegend.h>
28 #include <TString.h>
29 #include <TLatex.h>
30 
31 #include <Riostream.h>
32 #include <stdio.h>
33 #include <fstream>
34 #include <iostream>
35 #include <alloca.h>
36 #include <string>
37 #include <cstring>
38 
39 // --- ANALYSIS system ---
40 #include <AliCalorimeterUtils.h>
41 #include <AliEMCALGeometry.h>
43 #include <AliAODEvent.h>
44 
45 using namespace std;
46 
50 
54 //________________________________________________________________________
56  TObject(),
57  fCurrentRunNumber(-1),fPeriod(),fPass(),fTrigger(),fNoOfCells(),fCellStartDCal(),
58  fMergeOutput(),fAnalysisOutput(),fAnalysisInput(),fRunList(),
59  fQADirect(), fMergedFileName(), fAnalysisVector(),
60  fRunListFileName(),fWorkdir(),fTrial(),fExternalFileName(),fTestRoutine(),
61  fNMaxCols(),fNMaxRows(),fNMaxColsAbs(),fNMaxRowsAbs(),
62  fFlag(),fCaloUtils()
63 {
64  fCurrentRunNumber = 254381;
65  fPeriod = "LHC16h";
66  fPass = "muon_caloLego";
67  fTrigger = "AnyINT";
68  fWorkdir = "./";
69  fRunListFileName = "runList.txt";
70 
71  Init();
72 }
73 
77 //________________________________________________________________________
78 AliAnaCaloChannelAnalysis::AliAnaCaloChannelAnalysis(TString period, TString pass, TString trigger, Int_t runNumber, TString workDir, TString listName):
79  TObject(),
80  fCurrentRunNumber(-1),fPeriod(),fPass(),fTrigger(),fNoOfCells(),fCellStartDCal(),
81  fMergeOutput(),fAnalysisOutput(),fAnalysisInput(),fRunList(),
82  fQADirect(), fMergedFileName(), fAnalysisVector(),
83  fRunListFileName(),fWorkdir(),fTrial(),fExternalFileName(),fTestRoutine(),
84  fNMaxCols(),fNMaxRows(),fNMaxColsAbs(),fNMaxRowsAbs(),
85  fFlag(),fCaloUtils()
86 {
87  fCurrentRunNumber = runNumber;
88  fPeriod = period;
89  fPass = pass;
90  fTrigger = trigger;
91  fWorkdir = workDir;
92  fRunListFileName = listName;
93 
94  Init();
95 }
96 
100 //________________________________________________________________________
102 {
103  //......................................................
104  //..Default values - can be set by functions
105  fTrial = 0;
107  fTestRoutine=0;
108 
109  //..Settings for the input/output structure (hard coded)
110  fAnalysisInput ="AnalysisInput";
111  fMergeOutput ="ConvertOutput";
112  fAnalysisOutput ="AnalysisOutput";
113  //..Stuff for the convert function
114  gSystem->mkdir(fMergeOutput);
115  gSystem->mkdir(fAnalysisOutput);
116  fMergedFileName= Form("%s/%s_%s_Merged.root", fMergeOutput.Data(), fPeriod.Data(),fPass.Data());
117  fQADirect = Form("CaloQA_%s",fTrigger.Data());
118  fRunList = Form("%s/%s/%s/%s", fAnalysisInput.Data(), fPeriod.Data(), fPass.Data(), fRunListFileName.Data());
119 
120  //.. make sure the vector is empty
121  fAnalysisVector.clear();
122 
123  //......................................................
124  //..Initialize EMCal/DCal geometry
126  //..Create a dummy event for the CaloUtils
127  AliAODEvent* aod = new AliAODEvent();
130 
131  fNoOfCells =fCaloUtils->GetEMCALGeometry()->GetNCells(); //..Very important number, never change after that point!
132  Int_t nModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
133  fCellStartDCal=12288; //..ELI this should be automatized from the geometry information!!
134 
135  //..This is how the calorimeter looks like in the current period (defined by example run ID fCurrentRunNumber)
136  cout<<"Called geometry for run number: "<<fCurrentRunNumber<<endl;
137  cout<<"Number of supermod: "<<nModules<<endl;
138  cout<<"Number of cells: "<<fNoOfCells<<endl;
139  cout<<"Cell ID of first DCal cell: "<<fCellStartDCal<<endl;
140  //cout<<"Number of supermod utils: "<<fCaloUtils->GetNumberOfSuperModulesUsed()<<endl; //..will always be 22 unless set by hand
141 
142  //......................................................
143  //..Initialize flag array to store how the cell is categorized
144  //..In the histogram: bin 1= cellID 0, bin 2= cellID 1 etc
145  //..In the array: fFlag[cellID]= some information
146  fFlag = new Int_t[fNoOfCells];
147  fFlag[fNoOfCells] = {0}; //..flagged as good by default
148 
149  //......................................................
150  //..setings for the 2D histogram
151  fNMaxCols = 48; //eta direction
152  fNMaxRows = 24; //phi direction
154  fNMaxRowsAbs = Int_t (nModules/2)*fNMaxRows; //multiply by number of supermodules
155 }
156 
170 //________________________________________________________________________
172 {
173  if(fExternalFileName=="")
174  {
175  //..If no extrenal file is provided merge different runs together
176  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
177  cout<<". . .Start process by converting files. . . . . . . . . . . ."<<endl;
178  cout<<endl;
180  cout<<endl;
181  }
182  else
183  {
184  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
185  cout<<". . .Start process by loading external file. . . . . . . . . . ."<<endl;
186  fMergedFileName = Form("%s/%s", fMergeOutput.Data(), fExternalFileName.Data());
187  }
188  cout<<". . .Load inputfile with name: "<<fMergedFileName<<" . . . . . . . ."<<endl;
189 
190  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
191  //.. Read all the needed input for the Bad Channel analysis
192  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
193  TFile *mergedFileInput = new TFile(fMergedFileName);
194  if(!mergedFileInput->IsOpen()){
195  Printf("Error! Input file not found, abort");
196  return;
197  }
198  fCellAmplitude = (TH2F*) mergedFileInput->Get("hCellAmplitude");
199  fCellTime = (TH2F*) mergedFileInput->Get("hCellTime");
200  fProcessedEvents = (TH1F*) mergedFileInput->Get("hNEventsProcessedPerRun");
201 
202  cout<<". . .Continue process by . . . . . . . . . . . ."<<endl;
203  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
204  //.. DEAD CELLS
205  //.. Flag dead cells with fFlag=1
206  //.. this excludes cells from analysis (will not appear in results)
207  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
208  cout<<"o o o Flag Dead Cells o o o"<<endl;
209  FlagAsDead();
210  cout<<endl;
211  cout<<endl;
212 
213  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
214  //.. BAD CELLS
215  //.. Flag dead cells with fFlag=2,3
216  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
217  cout<<"o o o Bad channel analysis o o o"<<endl;
218  BCAnalysis();
219 
220  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
221  //..In the end summarize results
222  //..in a .pdf and a .txt file
223  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
225 
226  cout<<endl;
227  cout<<". . .End of process . . . . . . . . . . . . . . . . . . . . ."<<endl;
228  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
229 }
230 
237 //________________________________________________________________________
239 {
240  cout<<"o o o Start conversion process o o o"<<endl;
241  cout<<"o o o period: " << fPeriod << ", pass: " << fPass << ", trigger: "<<fTrigger<< endl;
242 
243  //..Create histograms needed for adding all the files together
244  TH1D *hNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun","Number of processed events vs run number",200000,100000,300000);
245  //ELI a little problematic to hard code properties of histograms??
246  TH2F *hCellAmplitude = new TH2F("hCellAmplitude","Cell Amplitude",200,0,10,23040,0,23040);
247  TH2F *hCellTime = new TH2F("hCellTime","Cell Time",250,-275,975,23040,0,23040);
248 
249  //..Open the text file with the run list numbers and run index
250  cout<<"o o o Open .txt file with run indices. Name = " << fRunList << endl;
251  FILE *pFile = fopen(fRunList.Data(), "r");
252  if(!pFile)cout<<"count't open file!"<<endl;
253  Int_t Nentr;
254  Int_t q;
255  Int_t ncols;
256  Int_t nlines = 0 ;
257  Int_t RunId[500] ;
258  while (1)
259  {
260  ncols = fscanf(pFile," %d ",&q);
261  if (ncols< 0) break;
262  RunId[nlines]=q;
263  nlines++;
264  }
265  fclose(pFile);
266 
267 
268  //..Open the different .root files with help of the run numbers from the text file
269  const Int_t nRun = nlines ;
270  TString base;
271  TString infile;
272 
273  cout<<"o o o Start merging process of " << nRun <<" files"<< endl;
274  //..loop over the amount of run numbers found in the previous text file.
275  for(Int_t i = 0 ; i < nRun ; i++)
276  {
277  base = Form("%s/%s/%s/%d", fAnalysisInput.Data(), fPeriod.Data(), fPass.Data(), RunId[i]);
278  if ((fPass=="cpass1_pass2")||(fPass=="cfPass1-2"))
279  {
280  if (fTrigger=="default")
281  {
282  infile = Form("%s_barrel.root",base.Data());
283  }
284  else
285  {
286  infile = Form("%s_outer.root",base.Data());
287  }
288  }
289  else
290  { //..This is a run2 case
291  infile = Form("%s.root",base.Data()) ;
292  }
293 
294  cout<<" o Open .root file with name: "<<infile<<endl;
295  TFile *f = TFile::Open(infile);
296 
297  //..Do some basic checks
298  if(!f)
299  {
300  cout<<"Couldn't open/find .root file: "<<infile<<endl;
301  continue;
302  }
303  TDirectoryFile *dir = (TDirectoryFile *)f->Get(fQADirect);
304  if(!dir)
305  {
306  cout<<"Couln't open directory file in .root file: "<<infile<<", directory: "<<fQADirect<<endl;
307  continue;
308  }
309  TList *outputList = (TList*)dir->Get(fQADirect);
310  if(!outputList)
311  {
312  cout << "Couln't get TList from directory file: "<<fQADirect<<endl;
313  continue;
314  }
315 
316  //ELI should one maybe clone the hAmpId histos eg to hCellAmplitude, then one does't need to hard code them.
317  TH2F *hAmpId;
318  TH2F *hTimeId;
319  TH2F *hNEvents;
320 
321  hAmpId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
322  if(!hAmpId)
323  {
324  Printf("hAmpId not found");
325  outputList->ls();
326  continue;
327  }
328  hTimeId =(TH2F *)outputList->FindObject("EMCAL_hTimeId");
329  if(!hTimeId)
330  {
331  Printf("hTimeId not found");
332  outputList->ls();
333  continue;
334  }
335  hNEvents =(TH2F *)outputList->FindObject("hNEvents");
336  if(!hNEvents)
337  {
338  Printf("hNEvents not found");
339  outputList->ls();
340  continue;
341  }
342  Nentr = (Int_t)hNEvents->GetEntries();
343 
344  //..does that mean do not merge small files?
345  if (Nentr<100)
346  {
347  cout <<" o File to small to be merged. Only N entries " << Nentr << endl;
348  continue ;
349  }
350  cout <<" o File with N entries " << Nentr<<" will be merged"<< endl;
351 
352  hNEventsProcessedPerRun->SetBinContent(RunId[i]-100000,(Double_t)Nentr);
353  hCellAmplitude->Add(hAmpId);
354  hCellTime->Add(hTimeId);
355 
356  outputList->Delete();
357  dir->Delete();
358  f->Close();
359  delete f;
360  }
361 
362  //.. Save the merged histograms
363  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
364  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
365  hNEventsProcessedPerRun->Write();
366  hCellAmplitude->Write();
367  hCellTime->Write();
368  BCF->Close();
369  cout<<"o o o End conversion process o o o"<<endl;
370  return fMergedFileName;
371 }
372 
373 
379 //_________________________________________________________________________
381 {
382  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
383  //.. BAD CELLS
384  //.. Flag bad cells with fFlag= 2 or 3
385  //.. this excludes cells from subsequent analysis
386  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
387  TArrayD periodArray;
388  for(Int_t i=0;i<fAnalysisVector.size();i++)
389  {
390  periodArray=fAnalysisVector.at(i);
391  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
392  cout<<""<<endl;
393  cout<<""<<endl;
394  }
395  cout<<"o o o End of bad channel analysis o o o"<<endl;
396 }
397 
406 //________________________________________________________________________
407 void AliAnaCaloChannelAnalysis::AddPeriodAnalysis(Int_t criteria, Double_t nsigma, Double_t emin, Double_t emax)
408 {
409  TArrayD periodArray;
410  periodArray.Set(4);
411  periodArray.AddAt((Double_t)criteria,0);
412  periodArray.AddAt(nsigma,1);
413  periodArray.AddAt(emin,2);
414  periodArray.AddAt(emax,3);
415  fAnalysisVector.push_back(periodArray);
416 }
426 //____________________________________________________________________
427 void AliAnaCaloChannelAnalysis::PeriodAnalysis(Int_t criterion, Double_t nsigma, Double_t emin, Double_t emax)
428 {
429  //ELI criterion should be between 1-4
430 
431  cout<<"o o o o o o o o o o o o o o o o o o o o o o o o o"<<endl;
432  cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
433  cout<<"o o o Done in the energy range E "<<emin<<"-"<<emax<<endl;
434 
435  Int_t cellID, nb1=0, nb2=0;
436  TString output;
437 
438  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
439  //.. ANALYSIS OF CELLS WITH ENTRIES
440  //.. Build average distributions and fit them
441  //.. Three tests for bad cells:
442  //.. 1) Average energy per hit
443  //.. 2) Average hit per event
444  //.. 3) ...
445  //.. 4) ...
446  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
447  TH1F* histogram;
448  if(criterion < 6)cout<<"o o o Analyze average cell distributions o o o"<<endl;
449  //..For case 1 or 2
450  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax,nsigma);
451  //..For case 3, 4 or 5
452  else if (criterion < 6) TestCellShapes(criterion, emin, emax, nsigma);
453 
454  Int_t dnbins = 200;
455  if(criterion==1) FlagAsBad(criterion, histogram, nsigma, dnbins,-1);
456  if(criterion==2 && emin==0.5) FlagAsBad(criterion, histogram, nsigma, dnbins*9000,-1);
457  if(criterion==2 && emin>0.5) FlagAsBad(criterion, histogram, nsigma, dnbins*17,-1);
458 
459  /*
460  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval3);
461  if(criterion==4) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval1);
462  if(criterion==5) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval2);
463  */
464 
465  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
466  //.. RESULTS
467  //.. 1) Print the bad cells
468  //.. and write the results to a file
469  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
470 
471  //..Print the results on the screen and
472  //..write the results in a file
473  output.Form("%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt", fAnalysisOutput.Data(), criterion,emin,emax);
474  ofstream file(output, ios::out | ios::trunc);
475  if(!file)
476  {
477  cout<<"#### Major Error. Check the textfile!"<<endl;
478  }
479  file<<"Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
480  file<<"Bad by lower value : "<<endl;
481  cout<<" o bad cells by lower value (for cell E between "<<emin<<"-"<<emax<<")"<<endl;
482  cout<<" ";
483  nb1=0;
484  for(cellID=0;cellID<fNoOfCells;cellID++)
485  {
486  if(fFlag[cellID]==2)
487  {
488  nb1++;
489  file<<cellID<<", ";
490  }
491  }
492  file<<"("<<nb1<<")"<<endl;
493  cout<<"("<<nb1<<")"<<endl;
494  file<<"Bad by higher value : "<<endl;
495  cout<<" o bad cells by higher value (for cell E between "<<emin<<"-"<<emax<<")"<<endl;
496  cout<<" ";
497  nb2=0;
498  for(cellID=0;cellID<fNoOfCells;cellID++)
499  {
500  if(fFlag[cellID]==3)
501  {
502  nb2++;
503  file<<cellID<<", ";
504  }
505  }
506  file<<"("<<nb2<<")"<<endl;
507  cout<<"("<<nb2<<")"<<endl;
508 
509  file<<"Total number of bad cells"<<endl;
510  file<<"("<<nb1+nb2<<")"<<endl;
511  file.close();
512  cout<<" o Total number of bad cells "<<endl;
513  cout<<" ("<<nb1+nb2<<")"<<endl;
514 }
515 
520 //_________________________________________________________________________
521 TH1F* AliAnaCaloChannelAnalysis::BuildHitAndEnergyMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
522 {
523  cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
524  TH1F *histogram;
525  if(crit==1)histogram = new TH1F(Form("hCellEtoNtotal_E%.2f-%.2f",emin,emax),Form("Energy per hit, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
526  if(crit==2)histogram = new TH1F(Form("hCellNtotal_E%.2f-%.2f",emin,emax),Form("Number of hits per event, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
527  histogram->SetXTitle("Abs. Cell Id");
528  if(crit==1)histogram->SetYTitle("Energy per hit");
529  if(crit==2)histogram->SetYTitle("Number of hits per event");
530  histogram->GetXaxis()->SetNdivisions(505);
531 
532  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
533 
534  //..here the average hit per event and the average energy per hit is caluclated for each cell.
535  for (Int_t cell = 0; cell < fNoOfCells; cell++)
536  {
537  Double_t Esum = 0;
538  Double_t Nsum = 0;
539 
540  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
541  {
542  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
543  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
544  //..investigate only cells that were not flagged as dead ore bad
545  if (E < emin || E > emax || fFlag[cell]!=0) continue;
546  Esum += E*N;
547  Nsum += N;
548  }
549  //..Set the values only for cells that are not yet marked as bad
550  if(fFlag[cell]==0)
551  {
552  if(totalevents> 0. && crit==2)histogram->SetBinContent(cell+1, Nsum/totalevents); //..number of hits per event
553  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/Nsum); //..average energy per hit
554  }
555  }
556  return histogram;
557 }
562 //_________________________________________________________________________
563 TH1F* AliAnaCaloChannelAnalysis::BuildTimeMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
564 {
565  TH1F *Histogram = new TH1F(Form("hSomethingWithTime_E%.2f-%.2f",emin,emax),Form("I don't know, %.2f < E < %.2f GeV",emin,emax), 2,0,2);
566 
567  return Histogram;
568 }
574 //_________________________________________________________________________
575 void AliAnaCaloChannelAnalysis::TestCellShapes(Int_t crit, Double_t fitemin, Double_t fitemax, Double_t nsigma)
576 {
577  Int_t dnbins = 1000;
578  // binning parameters
579  Int_t ncells = fCellAmplitude->GetNbinsY();
580  Double_t amin = fCellAmplitude->GetYaxis()->GetXmin();
581  Double_t amax = fCellAmplitude->GetYaxis()->GetXmax();
582  cout << "ncells " << ncells << " amin = " << amin << "amax = " << amax<< endl;
583 
584  // initialize histograms
585  TH1 *hFitA = new TH1F("hFitA_hCellAmplitude","Fit A value", ncells,amin,amax);
586  hFitA->SetXTitle("AbsId");
587  hFitA->SetYTitle("A");
588 
589  TH1 *hFitB = new TH1F("hFitB_hCellAmplitude","Fit B value", ncells,amin,amax);
590  hFitB->SetXTitle("AbsIdhname");
591  hFitB->SetYTitle("B");
592 
593  TH1 *hFitChi2Ndf = new TH1F("hFitChi2Ndf_hCellAmplitude","Fit #chi^{2}/ndf value", ncells,amin,amax);
594  hFitChi2Ndf->SetXTitle("AbsId");
595  hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
596 
597  Double_t maxval1=0., maxval2=0., maxval3=0.;
598  Double_t prev=0., MSA=0., AvA = 0. ; //those param are used to automaticaly determined a reasonable maxval1
599  Double_t prev2=0., MSB=0., AvB = 0. ; //those param are used to automaticaly determined a reasonable maxval2
600  Double_t prev3=0., MSki2=0., Avki2 = 0. ; //those param are used to automaticaly determined a reasonable maxval3
601  Double_t ki2=0.0 ;
602  for (Int_t k = 1; k <= fNoOfCells; k++)
603  {
604  TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
605  TH1 *hCell = fCellAmplitude->ProjectionX("",k,k);
606  if (hCell->GetEntries() == 0) continue;
607  // hCell->Rebin(3);
608  hCell->Fit(fit, "0QEM", "", fitemin, fitemax);
609  delete hCell;
610 
611  if(fit->GetParameter(0) < 5000.)
612  {
613  hFitA->SetBinContent(k, fit->GetParameter(0));
614  if(k<3000)
615  {
616  AvA += fit->GetParameter(0);
617  if(k==2999) maxval1 = AvA/3000. ;
618  if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
619  else MSA -= (fit->GetParameter(0) - prev) ;
620  prev = fit->GetParameter(0);
621  }
622  else
623  {
624  if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
625  {
626  maxval1 = fit->GetParameter(0);
627  }
628  }
629  }
630  else hFitA->SetBinContent(k, 5000.);
631 
632  if(fit->GetParameter(1) < 5000.)
633  {
634  hFitB->SetBinContent(k, fit->GetParameter(1));
635  if(k<3000)
636  {
637  AvB += fit->GetParameter(1);
638  if(k==2999) maxval2 = AvB/3000. ;
639  if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
640  else MSB -= (fit->GetParameter(1) - prev2) ;
641  prev2 = fit->GetParameter(1);
642  }
643  else
644  {
645  if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
646  {
647  maxval2 = fit->GetParameter(1);
648  }
649  }
650  }
651  else hFitB->SetBinContent(k, 5000.);
652 
653 
654  if (fit->GetNDF() != 0 ) ki2 = fit->GetChisquare()/fit->GetNDF();
655  else ki2 = 1000.;
656 
657  if(ki2 < 1000.)
658  {
659  hFitChi2Ndf->SetBinContent(k, ki2);
660  if(k<3000)
661  {
662  Avki2 += ki2;
663  if(k==2999) maxval3 = Avki2/3000. ;
664  if (prev3 < ki2) MSki2 += ki2 - prev3;
665  else MSki2 -= (ki2 - prev3) ;
666  prev3 = ki2;
667  }
668  else
669  {
670  if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
671  {
672  maxval3 = ki2;
673  }
674  }
675  }
676  else hFitChi2Ndf->SetBinContent(k, 1000.);
677 
678  delete fit ;
679  }
680 
681  // if you have problem with automatic parameter :
682  // maxval1 =
683  // maxval2 =
684  // maxval3 =
685  /*
686  if(crit==3)
687  Process(crit, hFitChi2Ndf, nsigma, dnbins, maxval3);
688  if(crit==4)
689  Process(crit, hFitA, nsigma, dnbins, maxval1);
690  if(crit==5)
691  Process(crit, hFitB, nsigma, dnbins, maxval2);
692  */
693 
694  //ELI something like thie in the future:
695  //return histogram;
696 }
697 
698 
703 //_________________________________________________________________________
705 {
706  Int_t sumOfExcl=0;
707  //..Direction of cell ID
708  for (Int_t cell = 0; cell < fNoOfCells; cell++)
709  {
710  Double_t nSum = 0;
711  //..Direction of amplitude (Checks energies from 0-10 GeV)
712  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
713  {
714  //..cellID+1 = histogram bin
715  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
716  nSum += N;
717  }
718  //..If the amplitude in one cell is basically 0
719  //..mark the cell as excluded
720  if(nSum < 0.5 && fFlag[cell]==0)
721  {
722  //..Cell flagged as dead.
723  //..Flag only if it hasn't been flagged before
724  fFlag[cell] =1;
725  sumOfExcl++;
726  }
727  }
728  cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
729  cout<<" ("<<sumOfExcl<<")"<<endl;
730 }
731 
746 //_________________________________________________________________________
747 void AliAnaCaloChannelAnalysis::FlagAsBad(Int_t crit, TH1* inhisto, Double_t nsigma, Int_t dnbins, Double_t dmaxval)
748 {
749  gStyle->SetOptStat(1); // MG modif
750  gStyle->SetOptFit(1); // MG modif
751 
752  if(crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
753  if(crit==2)cout<<" o Fit average hit per event distribution"<<endl;
754 
755  Int_t cellColumn=0,cellRow=0;
756  Int_t cellColumnAbs=0,cellRowAbs=0;
757  Int_t trash;
758 
759  TString histoName=inhisto->GetName();
760  Double_t goodmax= 0. ;
761  Double_t goodmin= 0. ;
762  if (dmaxval < 0.)
763  {
764  dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
765  if(crit==2 && dmaxval > 1) dmaxval =1. ;
766  }
767 
768  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
769  //. . .build the distribution of average values
770  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
771  TH1 *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
772  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
773  distrib->SetYTitle("Entries");
774  TH1 *distrib_EMCal = new TH1F(Form("%sDistr_EMCal",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
775  TH1 *distrib_DCal = new TH1F(Form("%sDistr_DCal",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
776 
777  //EMCAL and DCAL are similar - ELI separate TRD supportstructer areas into two different histograms
778  //..TRD support structure:
779  //..collumn 3,4,5,6,7 32,33,34,35 57,58,59 84,85,86,87,88
780  //..row 21,22,23 44,45,46 68,69,70 92,93,94 116,117 126 148,149,150 173,174 197
781 
782  //..build two dimensional histogram with values row vs. column
783  TH2F *plot2D = new TH2F(Form("%s_HitRowColumn",(const char*)histoName),Form("%s_HitRowColumn",(const char*)histoName),fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
784  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
785  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
786 
787  for (Int_t cell = 0; cell < fNoOfCells; cell++)
788  {
789  //..Do that only for cell ids also accepted by the code
790  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
791 
792  //..Fill histograms only for cells that are not yet flagged as bad
793  if(fFlag[cell]==0)
794  {
795  //..fill the distribution of avarge cell values
796  distrib->Fill(inhisto->GetBinContent(cell+1));
797  //if(cell<fCellStartDCal)distrib_EMCal->Fill(inhisto->GetBinContent(cell+1));
798  //else distrib_DCal ->Fill(inhisto->GetBinContent(cell+1));
799  //..Get Row and Collumn for cell ID
800  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
801  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
802  {
803  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
804  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
805  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
806  }
807  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
808  // check TRD support structure
809  if((cellColumnAbs>2 && cellColumnAbs<8) || (cellColumnAbs>31 && cellColumnAbs<36) || (cellColumnAbs>56 && cellColumnAbs<60) || (cellColumnAbs>83 && cellColumnAbs<89) ||
810  (cellRowAbs>20 && cellRowAbs<24) || (cellRowAbs>43 && cellRowAbs<47) || (cellRowAbs>67 && cellRowAbs<71) || (cellRowAbs>91 && cellRowAbs<95) ||
811  (cellRowAbs>115 && cellRowAbs<118)|| cellRowAbs==126 || (cellRowAbs>147 && cellRowAbs<151) || (cellRowAbs>172 && cellRowAbs<175) || cellRowAbs==197
812  )
813  {
814  distrib_EMCal->Fill(inhisto->GetBinContent(cell+1));
815  }
816  else
817  {
818  distrib_DCal ->Fill(inhisto->GetBinContent(cell+1));
819  }
820  }
821  }
822 
823  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
824  //. . .draw histogram + distribution
825  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
826 
827  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
828  c1->ToggleEventStatus();
829  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
830  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
831  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
832  upperPad->Draw();
833  lowerPadLeft->Draw();
834  lowerPadRight->Draw();
835 
836  upperPad->cd();
837  upperPad->SetLeftMargin(0.045);
838  upperPad->SetRightMargin(0.05);
839  upperPad->SetLogy();
840  inhisto->SetTitleOffset(0.6,"Y");
841  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
842 
843  inhisto->SetLineColor(kBlue+1);
844  inhisto->Draw();
845 
846  lowerPadRight->cd();
847  lowerPadRight->SetLeftMargin(0.09);
848  lowerPadRight->SetRightMargin(0.1);
849  plot2D->Draw("colz");
850 
851  lowerPadLeft->cd();
852  lowerPadLeft->SetLeftMargin(0.09);
853  lowerPadLeft->SetRightMargin(0.07);
854  lowerPadLeft->SetLogy();
855  distrib->SetLineColor(kBlue+1);
856  distrib->Draw();
857  distrib_EMCal->SetLineColor(kGreen+1);
858  distrib_EMCal->DrawCopy("same");
859  distrib_DCal->SetLineColor(kMagenta+1);
860  distrib_DCal->DrawCopy("same");
861 
862  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
863  //. . .fit histogram
864  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
865  Int_t higherbin=0,i;
866  for(i = 2; i <= dnbins; i++)
867  {
868  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i)) higherbin = i ;
869  }
870  //..good range is around the max value as long as the
871  //..bin content is larger than 2 entries
872  for(i = higherbin ; i<=dnbins ; i++)
873  {
874  if(distrib->GetBinContent(i)<2) break ;
875  goodmax = distrib->GetBinCenter(i);
876  }
877  for(i = higherbin ; i>1 ; i--)
878  {
879  if(distrib->GetBinContent(i)<2) break ;
880  goodmin = distrib->GetBinLowEdge(i);
881  }
882  //cout<<"higherbin : "<<higherbin<<endl;
883  //cout<<"good range : "<<goodmin<<" - "<<goodmax<<endl;
884 
885  TF1 *fit2 = new TF1("fit2", "gaus");
886  //..start the fit with a mean of the highest value
887  fit2->SetParameter(1,higherbin);
888 
889  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
890  Double_t sig, mean, chi2ndf;
891  // Marie midif to take into account very non gaussian distrig
892  mean = fit2->GetParameter(1);
893  sig = fit2->GetParameter(2);
894  chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
895 
896  if (mean <0.) mean=0.; //ELI is this not a highly problematic case??
897 
898  goodmin = mean - nsigma*sig ;
899  goodmax = mean + nsigma*sig ;
900 
901  if (goodmin<0) goodmin=0.;
902 
903  cout<<" o Result of fit: "<<endl;
904  cout<<" o "<<endl;
905  cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
906  cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
907 
908  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
909  //. . .Add info to histogram
910  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
911  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
912  lline->SetLineColor(kGreen+2);
913  lline->SetLineStyle(7);
914  lline->Draw();
915 
916  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
917  rline->SetLineColor(kGreen+2);
918  rline->SetLineStyle(7);
919  rline->Draw();
920 
921  TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
922  leg->AddEntry(lline, "Good region boundary","l");
923  leg->Draw("same");
924 
925  fit2->SetLineColor(kOrange-3);
926  fit2->SetLineStyle(1);//7
927  fit2->Draw("same");
928 
929  TLatex* text = 0x0;
930  if(crit==1) text = new TLatex(0.2,0.8,Form("Good range: %.2f-%.2f",goodmin,goodmax));
931  if(crit==2) text = new TLatex(0.2,0.8,Form("Good range: %.2f-%.2fx10^-5",goodmin*100000,goodmax*100000));
932  text->SetTextSize(0.06);
933  text->SetNDC();
934  text->SetTextColor(1);
935  text->Draw();
936  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
937  //. . .Save histogram
938  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
939  c1->Update();
940  TString name =Form("%s/criteria-_%d.gif", fAnalysisOutput.Data(), crit);
941  if(crit==1)name=Form("%s/AverageEperHit_%s.gif", fAnalysisOutput.Data(), (const char*)histoName);
942  if(crit==2)name=Form("%s/AverageHitperEvent_%s.gif", fAnalysisOutput.Data(), (const char*)histoName);
943  c1->SaveAs(name);
944 
945  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
946  //. . . Mark the bad cells in the fFlag array
947  //. . .(2= bad because cell average value lower than min allowed)
948  //. . .(3= bad because cell average value higher than max allowed)
949  //. . .(0 by default - good cell)
950  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
951  cout<<" o Flag bad cells that are outside the good range "<<endl;
952  for(Int_t cell = 0; cell < fNoOfCells; cell++)
953  {
954  //cel=0 and bin=1, cel=1 and bin=2
955  // <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
956  if (inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]!=1)
957  {
958  fFlag[cell]=2;
959  }
960  if (inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]!=1)
961  {
962  fFlag[cell]=3;
963  }
964  }
965  cout<<" o "<<endl;
966 }
967 
974 //________________________________________________________________________
976 {
977  Int_t cellID, nDCalCells = 0, nEMCalCells = 0;
978  TString cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells;
979 
980  deadPdfName = Form("%s/%s%sDC_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
981  badPdfName = Form("%s/%s%sBC_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
982  cellSummaryFile = Form("%s/%s%sBC_SummaryResults_V%i.txt", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial); ;
983  ratioOfBad = Form("%s/%s%sBCRatio_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
984  goodCells = Form("%s/%s%sGoodCells_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
985 
986  cout<<" o Final results o "<<endl;
987  cout<<" o write results into .txt file: "<<cellSummaryFile<<endl;
988  cout<<" o write results into .pdf file: "<<badPdfName<<endl;
989  ofstream file(cellSummaryFile, ios::out | ios::trunc);
990  if(file)
991  {
992  file<<"Dead cells : "<<endl;
993  cout<<" o Dead cells : "<<endl;
994  nEMCalCells =0;
995  nDCalCells =0;
996  for(cellID=0; cellID<fNoOfCells; cellID++)
997  {
998  if(cellID==0)
999  {
1000  file<<"In EMCal : "<<endl;
1001  cout<<" o In EMCal : "<<endl;
1002  }
1003  if(cellID==fCellStartDCal)
1004  {
1005  file<<"In DCal : "<<endl;
1006  cout<<endl;
1007  cout<<" o In DCal : "<<endl;
1008  }
1009  if(fFlag[cellID]==1)
1010  {
1011  file<<cellID<<"\n" ;
1012  cout<<cellID<<"," ;
1013  if(cellID<fCellStartDCal)nEMCalCells++;
1014  else nDCalCells++;
1015  }
1016  }
1017  cout<<endl;
1018  file<<"EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1019  cout<<" o EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1020 
1021  file<<"Bad cells: "<<endl;
1022  cout<<" o Bad cells: "<<endl;
1023  nEMCalCells =0;
1024  nDCalCells =0;
1025  for(cellID=0;cellID<fNoOfCells;cellID++)
1026  {
1027  if(cellID==0)
1028  {
1029  file<<"In EMCal : "<<endl;
1030  cout<<" o In EMCal : "<<endl;
1031  }
1032  if(cellID==fCellStartDCal)
1033  {
1034  file<<"In DCal : "<<endl;
1035  cout<<endl;
1036  cout<<" o In DCal : "<<endl;
1037  }
1038  if(fFlag[cellID]>1)
1039  {
1040  file<<cellID<<"\n" ;
1041  cout<<cellID<<"," ;
1042  if(cellID<fCellStartDCal)nEMCalCells++;
1043  else nDCalCells++;
1044  }
1045  }
1046  cout<<endl;
1047  file<<"EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1048  cout<<" o EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1049  }
1050  file.close();
1051 
1052  //cout<<" o Save the Dead channel spectra to a .pdf file"<<endl;
1053  //SaveBadCellsToPDF(0,deadPdfName);
1054  cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1055  SaveBadCellsToPDF(1,badPdfName) ;
1056  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
1057  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1058 
1059  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1060  //..Plot some summary canvases
1061  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1062 
1063  TCanvas *c1 = new TCanvas("CellProp","summary of cell properties",1000,500);
1064  c1->ToggleEventStatus();
1065  c1->Divide(2);
1066  c1->cd(1);
1067  //lowerPadRight->SetLeftMargin(0.09);
1068  //lowerPadRight->SetRightMargin(0.06);
1069  fCellAmplitude->Draw("colz");
1070  c1->cd(2);
1071  fCellTime->Draw("colz");
1072  c1->Update();
1073  TString name =Form("%s/CellProperties.gif", fAnalysisOutput.Data());
1074  c1->SaveAs(name);
1075 
1076 
1077  PlotFlaggedCells2D(0); //..all good cells
1078  PlotFlaggedCells2D(1); //..all dead cells
1079  PlotFlaggedCells2D(2,3); //..all bad cells
1080 }
1081 
1082 
1083 
1098 //________________________________________________________________________
1099 void AliAnaCaloChannelAnalysis::SaveBadCellsToPDF(Int_t version, TString pdfName)
1100 {
1101  gROOT->SetStyle("Plain");
1102  gStyle->SetOptStat(0);
1103  gStyle->SetFillColor(kWhite);
1104  gStyle->SetTitleFillColor(kWhite);
1105  gStyle->SetPalette(1);
1106 
1107  char title[100];
1108  char name[100];
1109 
1110  TH1 *hRefDistr = BuildMeanFromGood();
1111  Int_t firstCanvas=0;
1112  Bool_t suspicious;
1113  Bool_t candidate;
1114  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
1115  text->SetTextSize(0.06);
1116  text->SetTextColor(8);
1117  text->SetNDC();
1118  text->SetTextColor(1);
1119 
1120  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
1121  text2->SetTextSize(0.08);
1122  text2->SetTextColor(8);
1123 
1124  //..collect cells in an internal vector.
1125  //..when the vector is of size=9 or at the end of being filled
1126  //..plot the channels into a canvas
1127  std::vector<Int_t> channelVector;
1128  channelVector.clear();
1129  cout<<"Start printing into .pdf for version: "<<version<<endl;
1130  for(Int_t cell=0;cell<fNoOfCells;cell++)
1131  {
1132  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
1133  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
1134  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
1135  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
1136 
1137  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
1138  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1139  if(channelVector.size()==9 || cell == fNoOfCells-1)
1140  {
1141  cout<<"."<<flush;
1142 
1143  TString internal_pdfName=pdfName;
1144  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1145  if(channelVector.size() > 6) c1->Divide(3,3);
1146  else if (channelVector.size() > 3) c1->Divide(3,2);
1147  else c1->Divide(3,1);
1148 
1149  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1150  for(Int_t i=0; i<channelVector.size() ; i++)
1151  {
1152  sprintf(name, "Cell %d",channelVector.at(i)) ;
1153  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
1154  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
1155 
1156  c1->cd(i%9 + 1);
1157  c1->cd(i%9 + 1)->SetLogy();
1158  if(i==0)
1159  {
1160  leg->AddEntry(hRefDistr,"mean of good","l");
1161  leg->AddEntry(hCell,"current","l");
1162  }
1163  if(version>1)//..These are ratio plots of energy distr. of cell and mean of all good cells
1164  {
1165  hCell->Divide(hRefDistr);
1166  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
1167  candidate = CheckDistribution(hCell,hRefDistr);
1168  }
1169 
1170  hCell->SetLineColor(kBlue+1);
1171  hCell->GetXaxis()->SetTitle("E (GeV)");
1172  hCell->GetYaxis()->SetTitle("N Entries");
1173  hCell->GetXaxis()->SetRangeUser(0.,10.);
1174  hCell->SetLineWidth(1) ;
1175  hCell->SetTitle(title);
1176  hRefDistr->SetLineColor(kGray+2);
1177  hRefDistr->SetLineWidth(1);
1178 
1179  if(version!=10)hCell->Draw("hist");
1180  if(version==10)
1181  {
1182  if(fTestRoutine==0 && candidate==1)hCell->Draw("hist"); //..Draw those cells that are labled as a candidate for wam cells
1183  if(fTestRoutine==1) hCell->Draw("hist"); //..In case we are running in QA mode plot all Bad Cell distributions
1184  }
1185  if(version==1)hRefDistr->Draw("same") ;
1186 
1187  if(version==10 && candidate==1)text->Draw(); //..Draw a marker in the histogram that could be miscalibrated and labelled as warm
1188  if(version==2 && candidate==0)text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
1189 
1190  if(version==2 && candidate==0)leg->Draw();
1191  if(version<2)leg->Draw();
1192  }
1193 
1194  if(channelVector.size()<9 || cell == fNoOfCells-1)
1195  {
1196  internal_pdfName +=")";
1197  c1->Print(internal_pdfName.Data());
1198  }
1199  else if(firstCanvas==0)
1200  {
1201  internal_pdfName +="(";
1202  c1->Print(internal_pdfName.Data());
1203  firstCanvas=1;
1204  }
1205  else
1206  {
1207  c1->Print(internal_pdfName.Data());
1208  }
1209  delete c1;
1210  delete leg;
1211  channelVector.clear();
1212  }
1213  }
1214  delete hRefDistr;
1215  cout<<endl;
1216 }
1220 //_________________________________________________________________________
1222 {
1223  TH1* hGoodAmp;
1224  TH1* hgoodMean;
1225  Int_t NrGood=0;
1226  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1227  {
1228  if(fFlag[cell]!=0)continue;
1229  NrGood++;
1230  if(NrGood==1)hgoodMean = fCellAmplitude->ProjectionX("hgoodMean",cell+1,cell+1);
1231  else
1232  {
1233  hGoodAmp = fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1234  hgoodMean->Add(hGoodAmp);
1235  }
1236  }
1237  hgoodMean->Scale(1.0/NrGood);
1238  return hgoodMean;
1239 }
1247 //_________________________________________________________________________
1248 Bool_t AliAnaCaloChannelAnalysis::CheckDistribution(TH1* ratio, TH1* reference)
1249 {
1250  Double_t percentageOfLast=0.6;
1251  Double_t higherThanMean=5;
1252  Double_t highestRatio=1000;
1253  Double_t cliffSize=2; //before cliff shouldn't be double as high than after cliff
1254 
1255  //..By default each cell is a candidate for recalibration
1256  Bool_t candidate=1;
1257  //..Find bin where reference has value 1, and the corresponding x-value
1258  Int_t binHeihgtOne = reference->FindLastBinAbove(1);
1259  Double_t binCentreHeightOne = reference->GetBinCenter(binHeihgtOne);
1260  Double_t thirdBinCentre = reference->GetBinCenter(3);
1261 
1262  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1263  //..Check the histogram
1264  //..Different checks to see whether the
1265  //..cell is really bad. Set suspicious to 1.
1266  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1267 
1268  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1269  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
1270  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
1271  {
1272  candidate=0;
1273  }
1274 
1275  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
1276  //ELI check that crieteria carfully - seems to work but not shure about it
1277  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
1278  if(ratio->GetMaximum()>highestRatio)//
1279  {
1280  candidate=0;
1281  }
1282 
1283 
1284  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1285  //..check whether there are large spikes in the histogram
1286  //..compare bin values to mean of the ratio. If there is a bin value with
1287  //..content "higherThanMean" times lareger than mean it's suspicious
1288  Double_t mean=0;
1289  //..Find the maximum in the mean range (0-binHeihgtOne)
1290  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
1291  Double_t localMaxBin=ratio->GetMaximumBin();
1292 
1293  for(Int_t i=2;i<binHeihgtOne;i++)
1294  {
1295  //..Exclude 0 bins and exclude bins near the maximum
1296  if(ratio->GetBinContent(i)<=0) continue;
1297  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
1298  mean+=ratio->GetBinContent(i);
1299  }
1300  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1301  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
1302  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1303  if(ratio->GetMaximum()>mean*higherThanMean)
1304  {
1305  candidate=0;
1306  }
1307 
1308  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1309  //..Check for a cliff at 4GeV, happens in some cases
1310  Double_t beforeCliff=0;
1311  Double_t afterCliff=0;
1312  Int_t binsBefore=0;
1313  Int_t binsAfter=0;
1314  Int_t cliffBin = ratio->FindBin(4);
1315  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
1316  {
1317  //..Exclude 0 bins
1318  if(ratio->GetBinContent(i)<=0)continue;
1319  if(i<=cliffBin) binsBefore++;
1320  if(i>cliffBin) binsAfter++;
1321  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
1322  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
1323  beforeCliff*=1.0/binsBefore;
1324  afterCliff*=1.0/binsAfter;
1325  }
1326  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
1327  {
1328  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
1329  }
1330 
1331  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1332  //..Employ peak finder
1333 /* if(candidate==1)
1334  {
1335  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
1336  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
1337  //GetNPeaks();
1338  //cout<<"found N peaks: "<<nfound<<endl;
1339  }
1340 */
1341  return candidate;
1342 }
1347 //_________________________________________________________________________
1348 void AliAnaCaloChannelAnalysis::PlotFlaggedCells2D(Int_t flag1,Int_t flag2,Int_t flag3)
1349 {
1350  //..build two dimensional histogram with values row vs. column
1351  TString histoName = Form("HitRowColumn_Flag%d",flag1);
1352  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
1353  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1354  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1355 
1356  Int_t cellColumn=0,cellRow=0;
1357  Int_t cellColumnAbs=0,cellRowAbs=0;
1358  Int_t trash;
1359 
1360  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1361  {
1362  //..Do that only for cell ids also accepted by the code
1363  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1364 
1365  //..Get Row and Collumn for cell ID c
1366  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1367  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1368  {
1369  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1370  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1371  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1372  }
1373  if(fFlag[cell]==flag1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1374  if(flag2!=-1 && fFlag[cell]==flag2)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1375  if(flag3!=-1 && fFlag[cell]==flag3)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1376  }
1377  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
1378  c1->ToggleEventStatus();
1379  c1->cd();
1380  //lowerPadRight->SetLeftMargin(0.09);
1381  //lowerPadRight->SetRightMargin(0.06);
1382  plot2D->Draw("colz");
1383 
1384  TLatex* text = 0x0;
1385  if(flag1==0) text = new TLatex(0.2,0.8,"Good Cells");
1386  if(flag1==1) text = new TLatex(0.2,0.8,"Dead Cells");
1387  if(flag1>1) text = new TLatex(0.2,0.8,"Bad Cells");
1388  text->SetTextSize(0.06);
1389  text->SetNDC();
1390  text->SetTextColor(1);
1391  text->Draw();
1392 
1393  c1->Update();
1394  TString name =Form("%s/2DChannelMap_Flag%d.gif", fAnalysisOutput.Data(), flag1);
1395  c1->SaveAs(name);
1396 
1397 }
TH1F * BuildTimeMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
void SetRunNumber(Int_t run)
const char * title
Definition: MakeQAPdf.C:26
void AddPeriodAnalysis(Int_t criteria, Double_t nsigma, Double_t emin, Double_t emax)
TH2F * fCellTime
! possible histogram for the analysis. Cell ID vs. time, read from the input merged file ...
Int_t fNMaxRows
Maximum No of rows in module (phi direction)
Int_t fTrial
Number of trial that this specific analyis is. By default '0' so one can try different settings witho...
TSystem * gSystem
std::vector< TArrayD > fAnalysisVector
Vector of analysis information. Each place is filled with 4 doubles: version, sigma, lower, and upper energy range.
TH1F * fProcessedEvents
! Stores the number of events in the run
Int_t fNoOfCells
Number of cells in EMCal and DCal.
AliCalorimeterUtils * fCaloUtils
! Calorimeter information for the investigated runs
Int_t fCellStartDCal
ID of the first cell in the DCal.
TString fPass
Pass of the analyzed data.
Int_t * fFlag
! fFlag[CellID] = 0 (ok),1 (dead),2 (bad by lower),3 (bad by upper) start at 0 (cellID 0 = histobin 1...
Int_t fNMaxCols
Maximum No of colums in module (eta direction)
TString fMergeOutput
Here the merged files of a period are saved for a later analysis.
TString fMergedFileName
Filename of the .root file containing the merged runs.
AliEMCALGeometry * GetEMCALGeometry() const
TString fAnalysisInput
Here the .root files of each run of the period are saved.
void TestCellShapes(Int_t crit, Double_t fitemin, Double_t fitemax, Double_t nsigma=4.)
TString fQADirect
Dierctory in the QA.root files where the input histograms are stored.
TString fWorkdir
Directory which contains the folders fMergeOutput, fAnalysisInput and fAnalysisOutput. By default it is './'.
void FlagAsBad(Int_t crit, TH1 *inhisto, Double_t nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1.)
Double_t nsigma
TString fExternalFileName
If you have already a file that contains many runs merged together you can place it in fMergeOutput a...
TString fPeriod
The name of the analyzed period.
TH1F * BuildHitAndEnergyMean(Int_t crit, Double_t emin=0.1, Double_t emax=2., Double_t nsigma=4.)
TString fTrigger
Selected trigger for the analysis.
Bool_t CheckDistribution(TH1 *ratio, TH1 *reference)
Int_t fNMaxRowsAbs
Maximum No of rows in Calorimeter.
Int_t fNMaxColsAbs
Maximum No of colums in Calorimeter.
void PeriodAnalysis(Int_t criterum=7, Double_t nsigma=4.0, Double_t emin=0.1, Double_t emax=2.0)
Bool_t fTestRoutine
This is a flag, if set true will produce some extra quality check histograms.
Analyses cell properties and identifies bad cells.
TString fAnalysisOutput
The list with bad channels and histograms are saved in this folder.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TFile * file
Class with utils specific to calorimeter clusters/cells.
void AccessGeometry(AliVEvent *inputEvent)
TH2F * fCellAmplitude
! main histogram for the analysis. Cell ID vs. amplitude, read from the input merged file ...
void PlotFlaggedCells2D(Int_t flag1, Int_t flag2=-1, Int_t flag3=-1)
void SaveBadCellsToPDF(Int_t version, TString pdfName)
TString fRunListFileName
This is the name of the file with the run numbers to be merged, by default it's 'runList.txt'.
Int_t fCurrentRunNumber
A run number of an analyzed period. This is important for the AliCalorimeterUtils initialization...
TString fRunList
Thats the full path and name of the file which contains a list of all runs to be merged together...
Int_t GetModuleNumberCellIndexesAbsCaloMap(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU, Int_t &icolAbs, Int_t &irowAbs) const