AliPhysics  1168478 (1168478)
 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 <TLegend.h>
27 #include <TString.h>
28 #include <TSystem.h>
29 #include <TList.h>
30 #include <TLatex.h>
31 #include <TError.h>
32 
33 #include <Riostream.h>
34 #include <stdio.h>
35 #include <fstream>
36 #include <iostream>
37 #include <alloca.h>
38 #include <string>
39 #include <cstring>
40 
41 // --- ANALYSIS system ---
42 #include <AliCalorimeterUtils.h>
43 #include <AliEMCALGeometry.h>
45 #include <AliAODEvent.h>
46 
47 using std::vector;
48 using std::cout;
49 using std::endl;
50 using std::ofstream;
51 using std::flush;
52 using std::ios;
53 
57 
61 //________________________________________________________________________
63 TObject(),
64  fCurrentRunNumber(-1),
65  fPeriod(),
66  fTrainNo(),
67  fTrigger(),
68  fNoOfCells(),
69  fCellStartDCal(12288),
70  fStartCell(0),
71  fAnalysisOutput(),
72  fAnalysisInput(),
73  fRunList(),
74  fRunListFileName(),
75  fWorkdir(),
76  fQADirect(),
77  fMergedFileName(),
78  fAnalysisVector(),
79  fTrial(),
80  fExternalFileName(),
81  fTestRoutine(),
82  fPrint(0),
83  fNMaxCols(48),
84  fNMaxRows(24),
85  fNMaxColsAbs(),
86  fNMaxRowsAbs(),
87  fFlag(),
88  fCriterionCounter(),
89  fWarmCell(),
90  fCaloUtils(),
91  fRootFile(),
92  fCellAmplitude(),
93  fCellTime(),
94  fProcessedEvents(),
95  fhCellFlag(),
96  fhCellWarm()
97 {
98  fCurrentRunNumber = 254381;
99  fPeriod = "LHC16h";
100  fTrainNo = "muon_caloLego";
101  fTrigger = "AnyINT";
102  fWorkdir = ".";
103  fRunListFileName = "runList.txt";
104  fTrial = 0;
105 
106  Init();
107 }
108 
112 //________________________________________________________________________
113 AliAnaCaloChannelAnalysis::AliAnaCaloChannelAnalysis(TString period, TString train, TString trigger, Int_t runNumber,Int_t trial, TString workDir, TString listName):
114  TObject(),
115  fCurrentRunNumber(-1),
116  fPeriod(),
117  fTrainNo(),
118  fTrigger(),
119  fNoOfCells(),
120  fCellStartDCal(12288),
121  fStartCell(0),
122  fAnalysisOutput(),
123  fAnalysisInput(),
124  fRunList(),
125  fRunListFileName(),
126  fWorkdir(),
127  fQADirect(),
128  fMergedFileName(),
129  fAnalysisVector(),
130  fTrial(),
131  fExternalFileName(),
132  fTestRoutine(),
133  fPrint(0),
134  fNMaxCols(48),
135  fNMaxRows(24),
136  fNMaxColsAbs(),
137  fNMaxRowsAbs(),
138  fFlag(),
139  fCriterionCounter(),
140  fWarmCell(),
141  fCaloUtils(),
142  fRootFile(),
143  fCellAmplitude(),
144  fCellTime(),
145  fProcessedEvents(),
146  fhCellFlag(),
147  fhCellWarm()
148 {
149  fCurrentRunNumber = runNumber;
150  fPeriod = period;
151  fTrainNo = train; //only for folder structure
152  fTrigger = trigger; //important to select trigger in output file == different wagons in lego train
153  fWorkdir = workDir;
154  fRunListFileName = listName;
155  fTrial = trial;
156 
157  Init();
158 }
159 
163 //________________________________________________________________________
165 {
166  gROOT->ProcessLine("gErrorIgnoreLevel = kWarning;"); //Does not work -
167  //......................................................
168  //..Default values - can be set by functions
170  fTestRoutine=0;
171 
172  //..Settings for the input/output structure (hard coded)
173  // TO BE CHANGED
174  fAnalysisInput =Form("AnalysisInput/%s",fPeriod.Data());
175  fAnalysisOutput =Form("AnalysisOutput/%s/%s/Version%i",fPeriod.Data(),fTrainNo.Data(),fTrial);
176 
177  //..Make output directory if it doesn't exist
178  //..first the period folder
179  gSystem->mkdir(Form("%s/AnalysisOutput/%s",fWorkdir.Data(),fPeriod.Data()));
180  //..then the Train folder
181  gSystem->mkdir(Form("%s/AnalysisOutput/%s/%s",fWorkdir.Data(),fPeriod.Data(),fTrainNo.Data()));
182  //..then the version folder
183  gSystem->mkdir(Form("%s/%s",fWorkdir.Data(),fAnalysisOutput.Data()));
184 
185  fMergedFileName= Form("%s/%s/%s/MergedRuns_%s.root",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),fTrigger.Data());
186  fRunList = Form("%s/%s/%s/%s",fWorkdir.Data(), fAnalysisInput.Data(), fTrainNo.Data(), fRunListFileName.Data());
187  fQADirect = Form("CaloQA_%s",fTrigger.Data());
188 
189  TString fileName = Form("%s/%s/%s_%s_Histograms_V%i.root",fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data(),fTrigger.Data() ,fTrial);
190  fRootFile = new TFile(fileName,"recreate");
191  //.. make sure the vector is empty
192  fAnalysisVector.clear();
193 
194  //......................................................
195  //..Initialize EMCal/DCal geometry
197  //..Create a dummy event for the CaloUtils
198  AliAODEvent* aod = new AliAODEvent();
201 
202  fNoOfCells =fCaloUtils->GetEMCALGeometry()->GetNCells(); //..Very important number, never change after that point!
203  Int_t nModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
204 
205  //..This is how the calorimeter looks like in the current period (defined by example run ID fCurrentRunNumber)
206  cout<<"Called geometry for run number: "<<fCurrentRunNumber<<endl;
207  cout<<"Number of supermod: "<<nModules<<endl;
208  cout<<"Number of cells: "<<fNoOfCells<<endl;
209  cout<<"Cell ID of first DCal cell: "<<fCellStartDCal<<endl;
210  //cout<<"Number of supermod utils: "<<fCaloUtils->GetNumberOfSuperModulesUsed()<<endl; //..will always be 22 unless set by hand
211 
212  //......................................................
213  //..Initialize flag array to store how the cell is categorized
214  //..In the histogram: bin 1= cellID 0, bin 2= cellID 1 etc
215  //..In the array: fFlag[cellID]= some information
216  fFlag = new Int_t[fNoOfCells];
217  fWarmCell = new Bool_t[fNoOfCells];
218  fFlag[fNoOfCells] = {0}; //..flagged as good by default
219  fWarmCell[fNoOfCells] = {0}; //..flagged as not warm by default
220  fCriterionCounter=2; //This value will be written in fflag and updates after each PeriodAnalysis
221  //......................................................
222  //..setings for the 2D histogram
223  //fNMaxCols = 48; //eta direction
224  //fNMaxRows = 24; //phi direction
226  fNMaxRowsAbs = Int_t (nModules/2)*fNMaxRows; //multiply by number of supermodules
227 
228  //......................................................
229  //..Create TLists to organize the outputfile
230  fOutputListBad = new TList();
231  fOutputListGood = new TList();
232  fOutputListBadRatio = new TList();
233  fOutputListGoodRatio = new TList();
234 
235  fOutputListBad ->SetName("BadCell_Amplitudes");
236  fOutputListGood ->SetName("GoodCell_Amplitudes");
237  fOutputListBadRatio ->SetName("BadCell_AmplitudeRatios");
238  fOutputListGoodRatio->SetName("GoodCell_AmplitudeRatios");
239 
240  //fOutputListGood ->SetOwner();//ELI instead of delete in destructor??
241  //fOutputListBadRatio ->SetOwner();
242  //fOutputListGoodRatio ->SetOwner();
243 
244  //......................................................
245  //..Create Histograms to store the flag in a root file
246  fhCellFlag = new TH1F("fhCellFlag","fhCellFlag",fNoOfCells+10,0,fNoOfCells+10); //..cellID+1 = histogram bin
247  fhCellWarm = new TH1F("fhCellWarm","fhCellWarm",fNoOfCells+10,0,fNoOfCells+10); //..cellID+1 = histogram bin
248 }
249 
263 //________________________________________________________________________
265 {
266 // cout<<"fired trigger class"<<AliAODEvent::GetFiredTriggerClasses()<<endl;
267 
268  if(fExternalFileName=="")
269  {
270  //..If no extrenal file is provided merge different runs together
271  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
272  cout<<". . .Start process by converting files. . . . . . . . . . . ."<<endl;
273  cout<<endl;
275  if(fMergedFileName.IsNull())
276  {
277  Printf("File not produced, exit");
278  return;
279  }
280  cout<<endl;
281  }
282  else
283  {
284  //..If extrenal file is provided load it
285  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
286  cout<<". . .Start process by loading external file. . . . . . . . . . ."<<endl;
287  fMergedFileName= Form("%s/%s/%s/%s",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),fExternalFileName.Data());
288  }
289  //..if ==1 only produce filtered and merged files and do not perform a BC analysis
290  if(mergeOnly==0)
291  {
292  cout<<". . .Load inputfile with name: "<<fMergedFileName<<" . . . . . . . ."<<endl;
293  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
294  //.. Read all the needed input for the Bad/Dead Channel analysis
295  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
296  TFile *mergedFileInput = new TFile(fMergedFileName);
297  if(!mergedFileInput->IsOpen())
298  {
299  Printf("Error! Input file not found, abort");
300  return;
301  }
302  fCellAmplitude = (TH2F*) mergedFileInput->Get("hCellAmplitude");
303  fCellTime = (TH2F*) mergedFileInput->Get("hCellTime");
304  fProcessedEvents = (TH1F*) mergedFileInput->Get("hNEvents");
305 
306  cout<<". . .Continue process by . . . . . . . . . . . ."<<endl;
307  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
308  //.. DEAD CELLS
309  //.. Flag dead cells with fFlag=1
310  //.. this excludes cells from analysis (will not appear in results)
311  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
312  cout<<"o o o Flag dead cells o o o"<<endl;
313  FlagAsDead();
314  if(fPrint==1)cout<<endl;
315  if(fPrint==1)cout<<endl;
316 
317  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
318  //.. BAD CELLS
319  //.. Flag dead cells with fFlag=2 and bigger
320  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
321  cout<<"o o o Flag bad cells o o o"<<endl;
322  BCAnalysis();
323 
324  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
325  //..In the end summarize results
326  //..in a .pdf and a .txt file
327  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
328  if(fPrint==1)cout<<"o o o Write .txt for each period analyis with bad cells o o o"<<endl;
330 
331  if(fPrint==1)cout<<"o o o Create summary documents for the entire analysis o o o"<<endl;
333  fRootFile->WriteObject(fFlag,"FlagArray");
334  fRootFile->Close();
335  cout<<endl;
336 
337  //..make a reccomendation about the used energy range to be investigated
338  //..and the binning
339  TH1D *hRefDistr = BuildMeanFromGood();
340  Double_t totalevents = fProcessedEvents->Integral();
341  //..Find bin where reference has value "totalevents" (means 1 hit/event), and the corresponding x-value
342  Int_t binHeightOne = hRefDistr->FindLastBinAbove(1.0/totalevents);
343  Double_t binCentreHeightOne = hRefDistr->GetBinCenter(binHeightOne);
344  cout<<". . .Recomendation:"<<endl;
345  cout<<". . .With the current statistic on average a cell has 1 hit at "<<binCentreHeightOne<<" GeV"<<endl;
346  cout<<". . .so it makes no sense to select energy ranges >"<<binCentreHeightOne<<" as cells will be"<<endl;
347  cout<<". . .marked bad just due to the lack of statistic"<<endl;
348  cout<<". . .End of process . . . . . . . . . . . . . . . . . . . . ."<<endl;
349  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
350  }
351 }
352 
359 //________________________________________________________________________
361 {
362  cout<<"o o o Start conversion process o o o"<<endl;
363  cout<<"o o o period: " << fPeriod << ", pass: " << fTrainNo << ", trigger: "<<fTrigger<< endl;
364 
365  //..Create histograms needed for adding all the files together
366  TH1F *hNEventsProcessedPerRun= new TH1F("hNEvents","Number of processed events in analyzed runs",1,0,1);
367  TH2F *hCellAmplitude;
368  TH2F *hCellTime;
369 
370  //..Open the text file with the run list numbers and run index
371  cout<<"o o o Open .txt file with run indices. Name = " << fRunList << endl;
372  FILE *pFile = fopen(fRunList.Data(), "r");
373  if(!pFile)
374  {
375  cout<<"couldn't open file!"<<endl;
376  return "";
377  }
378  Int_t nEntr;
379  Int_t nEntrTot=0;
380  Int_t q;
381  Int_t ncols;
382  Int_t nlines = 0 ;
383  Int_t runId[500] ;
384  while (1)
385  {
386  ncols = fscanf(pFile," %d ",&q);
387  if (ncols< 0) break;
388  runId[nlines]=q;
389  nlines++;
390  }
391  fclose(pFile);
392 
393 
394  //..Open the different .root files with help of the run numbers from the text file
395  const Int_t nRun = nlines ;
396  TString base;
397  TString infile;
398  TString singleRunFileName;
399 
400  cout<<"o o o Start merging process of " << nRun <<" files"<< endl;
401  Int_t usedFiles=0;
402  //..loop over the amount of run numbers found in the previous text file.
403  for(Int_t i = 0 ; i < nRun ; i++)
404  {
405  base = Form("%s/%s/%s/%d", fWorkdir.Data(), fAnalysisInput.Data(), fTrainNo.Data(), runId[i]);
406  //..This is a run2 case
407  infile = Form("%s.root",base.Data()) ;
408 
409  cout<<" o Open .root file with name: "<<infile<<endl;
410  TFile *f = TFile::Open(infile);
411 
412  //..Do some basic checks
413  if(!f)
414  {
415  cout<<"Couldn't open/find .root file: "<<infile<<endl;
416  continue;
417  }
418  TDirectoryFile *dir = (TDirectoryFile *)f->Get(fQADirect);
419  if(!dir)
420  {
421  cout<<"Couln't open directory file in .root file: "<<infile<<", directory: "<<fQADirect<<endl;
422  continue;
423  }
424  TList *outputList = (TList*)dir->Get(fQADirect);
425  if(!outputList)
426  {
427  cout << "Couln't get TList from directory file: "<<fQADirect<<endl;
428  continue;
429  }
430  usedFiles++;
431  TH2F *hAmpId;
432  TH2F *hTimeId;
433  TH1F *hNEvents;
434 
435  hAmpId =(TH2F*)outputList->FindObject("EMCAL_hAmpId");
436  if(!hAmpId)
437  {
438  Printf("hAmpId not found");
439  outputList->ls();
440  continue;
441  }
442  else
443  {
444  hAmpId->SetName("hCellAmplitude");
445  hAmpId->SetTitle("Cell Amplitude");
446  }
447 
448  hTimeId =(TH2F*)outputList->FindObject("EMCAL_hTimeId");
449  if(!hTimeId)
450  {
451  Printf("hTimeId not found");
452  outputList->ls();
453  continue;
454  }
455  else
456  {
457  hTimeId->SetName("hCellTime");
458  hTimeId->SetTitle("Cell Time");
459  }
460 
461  if(i==0)
462  {
463  //..copy the properties of the mother histogram for adding them all up
464  hCellAmplitude=(TH2F*)hAmpId->Clone("DummyName1");
465  hCellAmplitude->Reset();
466  hCellAmplitude->SetDirectory(0); //otherwise histogram will dissapear when file is closed
467  //..copy the properties of the mother histogram for adding them all up
468  hCellTime=(TH2F*)hTimeId->Clone("DummyName2");
469  hCellTime->Reset();
470  hCellTime->SetDirectory(0); //otherwise histogram will dissapear when file is closed
471  }
472 
473  hNEvents =(TH1F *)outputList->FindObject("hNEvents");
474  if(!hNEvents)
475  {
476  Printf("hNEvents not found");
477  outputList->ls();
478  continue;
479  }
480  nEntr = (Int_t)hNEvents->GetEntries();
481 
482  //..does that mean do not merge small files?
483  if (nEntr<100)
484  {
485  cout <<" o File to small to be merged. Only N entries " << nEntr << endl;
486  continue ;
487  }
488  cout <<" o File with N entries " << nEntr<<" will be merged"<< endl;
489  nEntrTot+=nEntr;
490  hCellAmplitude->Add(hAmpId);
491  hCellTime->Add(hTimeId);
492  hNEventsProcessedPerRun->Add(hNEvents);
493 
494  //..Create copies of the original root files just with the bad channel QA
495  //..So that a run by run bad channel analysis can be performed more easily
496  singleRunFileName= Form("%s/%s/%s/%d_%sFiltered.root",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),runId[i],fTrigger.Data());
497  TFile *singleRunFile = TFile::Open(singleRunFileName,"recreate");
498  //..do not yet normalize by number of events
499  //..due to binning issues in later histograms
500  //..but rather do it at the very end of the analysis
501  hAmpId ->Write();
502  hTimeId->Write();
503  hNEvents->Write();
504  singleRunFile->Close();
505 
506  outputList->Delete();
507  dir->Delete();
508  f->Close();
509  delete f;
510  }
511  //..avoid creating empty files
512  if(usedFiles>0)
513  {
514  //.. Save the merged histograms
515  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
516  cout<<"o o o "<<nEntrTot<<" events were merged"<<endl;
517  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
518  //hNEventsProcessedPerRun->SetName("hNEvents");
519  //hNEventsProcessedPerRun->SetTitle("Number of processed events");
520  hNEventsProcessedPerRun->Write();
521  hCellAmplitude->SetName("hCellAmplitude");
522  hCellAmplitude->SetTitle("Cell Amplitude");
523  hCellAmplitude->Write();
524  hCellTime->SetName("hCellTime");
525  hCellTime->SetTitle("Cell Time");
526  hCellTime->Write();
527  BCF->Close();
528  cout<<"o o o End conversion process o o o"<<endl;
529  return fMergedFileName;
530  }
531  else
532  {
533  return "";
534  }
535 }
536 
537 
543 //_________________________________________________________________________
545 {
546  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
547  //.. BAD CELLS
548  //.. Flag bad cells with fFlag= 2,3,4,5.. etc
549  //.. this excludes cells from subsequent analysis
550  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
551  TArrayD periodArray;
552  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
553  {
554  periodArray=fAnalysisVector.at(i);
555  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
557  if(fPrint==1)cout<<endl;
558  if(fPrint==1)cout<<endl;
559  }
560  cout<<"o o o End of bad channel analysis o o o"<<endl;
561 }
562 
571 //________________________________________________________________________
573 {
574  TArrayD periodArray;
575  periodArray.Set(4);
576  periodArray.AddAt((Double_t)criteria,0);
577  periodArray.AddAt(nsigma,1);
578  periodArray.AddAt(emin,2);
579  periodArray.AddAt(emax,3);
580  fAnalysisVector.push_back(periodArray);
581 }
597 //____________________________________________________________________
599 {
600  //.. criterion should be between 1-4
601  if(fPrint==1)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;
602  if(fPrint==1)cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
603  if(fPrint==1 && criterion < 3)cout<<"o o o Done in the energy range E "<<emin<<" - "<<emax<<endl;
604  if(fPrint==1 && criterion == 3)cout<<"o o o Done in the time range t "<<emin<<" - "<<emax<<endl;
605 
606  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
607  //.. ANALYSIS OF CELLS WITH ENTRIES
608  //.. Build average distributions and fit them
609  //.. Three tests for bad cells:
610  //.. 1) Average energy per hit
611  //.. 2) Average hit per event
612  //.. 3) Average max of cell time distribution
613  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
614  TH1F* histogram;
615  if(fPrint==1)cout<<"o o o Analyze average cell distributions o o o"<<endl;
616  //..For case 1 or 2
617  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax);
618  if(criterion == 3) histogram = BuildTimeMean(criterion, emin, emax); //.. in case of crit 3 emin is tmin and emax is tmax
619 
620  if(criterion==1)
621  {
622  if(emin>2)FlagAsBad(criterion, histogram, nsigma, -1);//..do not apply a lower boundary
623  else FlagAsBad(criterion, histogram, nsigma, 400);
624  }
625  if(criterion==2)
626  {
627  if(emin>2)FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
628  else FlagAsBad(criterion, histogram, nsigma, 601);
629  //FlagAsBad(criterion, histogram, nsigma,601);// ELIANE ELIANE
630  //FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
631  }
632  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, 602);
633 }
634 
643 //_________________________________________________________________________
645 {
646  if(fPrint==1)cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
647  TH1F *histogram;
648  if(crit==1)histogram = new TH1F(Form("hCellEtoN_E%.2f-%.2f",emin,emax),Form("Energy per hit, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
649  if(crit==2)histogram = new TH1F(Form("hCellNtoEvt_E%.2f-%.2f",emin,emax),Form("Number of hits in cell, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
650  histogram->SetXTitle("Abs. Cell Id");
651  if(crit==1)histogram->SetYTitle("Energy per hit");
652  if(crit==2)histogram->SetYTitle("Number of hits in cell");
653  histogram->GetXaxis()->SetNdivisions(510);
654 
655  //..count Events in Energy Range
656  TH1F* pojection = (TH1F*)fCellAmplitude->ProjectionX("Intermediate");
657  fnEventsInRange = pojection->Integral(pojection->GetBinContent(emin),pojection->GetBinContent(emax));
658 
659  //..here the average hit per event and the average energy per hit is caluclated for each cell.
660  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
661  {
662  Double_t Esum = 0;
663  Double_t Nsum = 0;
664 
665  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
666  {
667  //To Do: I think one should also take the different bin width into account
668  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
669  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
670  //..investigate only cells that were not flagged as dead ore bad
671  if (E < emin || E > emax || fFlag[cell]!=0) continue;
672  Esum += E*N;
673  Nsum += N;
674  }
675  //..Set the values only for cells that are not yet marked as bad
676  if(fFlag[cell]==0)
677  {
678  if(crit==2) histogram->SetBinContent(cell+1, Nsum); //..number of hits (per event -> not anymore)
679  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/(Nsum)); //..average energy per hit
680  }
681  }
682  return histogram;
683 }
688 //_________________________________________________________________________
690 {
691  if(fPrint==1)cout<<" o Calculate maximum of cell time distribution "<<endl;
692  TString name;
693  TH1F *histogram;
694  histogram = new TH1F(Form("timeMax_t%.2f-%.2f",tmin,tmax),Form("Maximum of time distr., %.2f < t < %.2f GeV",tmin,tmax), fNoOfCells,0,fNoOfCells);
695  histogram->SetXTitle("Abs. Cell Id");
696  histogram->SetYTitle("time max");
697  histogram->GetXaxis()->SetNdivisions(510);
698 
699  //..Time information
700  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
701  {
702  Double_t timeMax = -100;
703  //..search for the maximum bin between tmin and tmax
704  name=Form("Cell %d",cell);
705  //sprintf(name, "Cell %d",channelVector.at(i)) ;
706  TH1 *hCell = fCellTime->ProjectionX(name,cell+1,cell+1);
707  hCell->GetXaxis()->SetRangeUser(tmin,tmax);
708 
709  timeMax = hCell->GetXaxis()->GetBinCenter(hCell->GetMaximumBin());
710  if(cell>55 &&cell<65)cout<<"Integral: "<<hCell->Integral()<<endl;
711  if(hCell->Integral()<=0)timeMax=-40;
712  //..Set the values only for cells that are not yet marked as bad
713  if(fFlag[cell]==0)
714  {
715  histogram->SetBinContent(cell+1, timeMax); //..number of hits per event
716  }
717  else
718  {
719  histogram->SetBinContent(cell+1, -50); //..number of hits per event
720  }
721 
722  }
723  return histogram;
724 }
725 
730 //_________________________________________________________________________
732 {
733  Int_t sumOfExcl=0;
734  Int_t manualMaskCounter=0;
735  //..sort the elements in manual mask list
736  std::sort (fManualMask.begin(), fManualMask.end());
737 
738  //..Direction of cell ID
739  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
740  {
741  Double_t nSum = 0;
742  //..Direction of amplitude (Checks energies from 0-nBins GeV)
743  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
744  {
745  //..cellID+1 = histogram bin
746  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
747  nSum += N;
748  }
749  //..If the amplitude in one cell is basically 0
750  //..mark the cell as excluded
751  if(nSum == 0 && fFlag[cell]==0)
752  {
753  //..Cell flagged as dead.
754  //..Flag only if it hasn't been flagged before
755  fFlag[cell] =1;
756  sumOfExcl++;
757  }
758  //..add here the manual masking
759  if(manualMaskCounter<(Int_t)fManualMask.size() && fManualMask.at(manualMaskCounter)==cell)
760  {
761  fFlag[cell] = 2;
762  manualMaskCounter++;
763  }
764  }
765  if(fPrint==1)cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
766  if(fPrint==1)cout<<" ("<<sumOfExcl<<")"<<endl;
767 }
768 
783 //_________________________________________________________________________
785 {
786  gStyle->SetOptStat(0); //..Do not plot stat boxes
787  gStyle->SetOptFit(0); //
788  if(fPrint==1 && crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
789  if(fPrint==1 && crit==2)cout<<" o Fit average hit per event distribution"<<endl;
790  if(fPrint==1 && crit==3)cout<<" o Fit average hit maximum distribution"<<endl;
791 
792  Int_t cellColumn=0,cellRow=0;
793  Int_t cellColumnAbs=0,cellRowAbs=0;
794  Int_t trash;
795 
796  TString histoName=inhisto->GetName();
797  Double_t goodmax = 0. ;
798  Double_t goodmin = 0. ;
799  //..set a user range so that the min and max is only found in a certain range
800  inhisto->GetXaxis()->SetRangeUser(fStartCell,fNoOfCells);//..
801  Double_t dminVal = inhisto->GetMinimum(0);
802  Double_t dmaxVal = inhisto->GetMaximum();
803  Double_t inputBins=dnbins;
804  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
805  //. . .determine settings for the histograms (range and binning)
806  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
807  //cout<<"max value: "<<dmaxVal<<", min value: "<<dminVal<<endl;
808  if(crit==2 && inputBins==-1) dnbins=dmaxVal-dminVal;
809  if(crit==1 && inputBins==-1) dnbins=400;
810 
811  if(crit==2 && inputBins!=-1)
812  {
813  //..calculate and print the "median"
814  Int_t numBins = inhisto->GetXaxis()->GetNbins();
815 
816  numBins-=fStartCell;
817  Double_t *x = new Double_t[numBins];
818  Double_t* y = new Double_t[numBins];
819  for (int i = 0; i < numBins; i++)
820  {
821  x[i] = inhisto->GetBinContent(i+fStartCell);
822  y[i] = 1; //each value with weight 1
823  if(x[i]==0)y[i] = 0;
824  }
825  Double_t medianOfHisto = TMath::Median(numBins,x,y);
826 
827  //..if dmaxVal is too far away from medianOfHisto the histogram
828  //..range will be too large -> reduce the range
829  //cout<<"max value: "<<dmaxVal<<" median of histogram: "<<medianOfHisto<<endl;
830  if(medianOfHisto*10<dmaxVal)
831  {
832  //cout<<"- - - median too far away from max range"<<endl;
833  dmaxVal=medianOfHisto+0.2*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
834  }
835  dnbins=dmaxVal-dminVal;
836 
837  if(dmaxVal-dminVal>100)
838  {
839  if(dnbins>2000)dnbins=0.01*(dmaxVal-dminVal); //..maximum 5000 bins. changed to 3000 .. lets see..
840  if(dnbins>2000)dnbins=0.001*(dmaxVal-dminVal);//..maximum 5000 bins.
841  if(dnbins<100) dnbins=0.02*(dmaxVal-dminVal); //..minimum 100 bins.
842  }
843  }
844 
845  if(crit==3)
846  {
847  //..obtain the binwidth for
848  Int_t binWidth = fCellTime->GetXaxis()->GetBinWidth(1);
849  dnbins = fCellTime->GetXaxis()->GetNbins();
850  dminVal= fCellTime->GetXaxis()->GetBinCenter(1)-(binWidth*1.0)/2.0;
851  dmaxVal= fCellTime->GetXaxis()->GetBinCenter(dnbins)+(binWidth*1.0)/2.0;
852  cout<<"set the new histo with following settings- bins: "<<dnbins<<", min val = "<<dminVal<<", max val:"<<dmaxVal<<endl;
853  }
854  //..build histos
855  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
856  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
857  distrib->SetYTitle("Entries");
858  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
859  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
860 
861  //..build two dimensional histogram with values row vs. column
862  TH2F *plot2D = new TH2F(Form("%s_HitRowColumn",(const char*)histoName),Form("%s_HitRowColumn",(const char*)histoName),fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
863  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
864  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
865 
866  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
867  //. . .build the distribution of average values
868  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
869  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
870  {
871  //..Do that only for cell ids also accepted by the code
872  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
873 
874  //..Fill histograms only for cells that are not yet flagged as bad
875  if(fFlag[cell]==0)
876  {
877  //..fill the distribution of avarge cell values
878  distrib->Fill(inhisto->GetBinContent(cell+1));
879  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
880  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
881  //..Get Row and Collumn for cell ID
882  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
883  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
884  {
885  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
886  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
887  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
888  }
889  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
890  //..check TRD support structure
891  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
892  {
893  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
894  }
895  else
896  {
897  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
898  }
899  }
900  }
901 
902  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
903  //. . .draw histogram + distribution
904  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
905  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
906  c1->ToggleEventStatus();
907  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
908  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
909  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
910  upperPad->Draw();
911  lowerPadLeft->Draw();
912  lowerPadRight->Draw();
913 
914  upperPad->cd();
915  upperPad->SetLeftMargin(0.05);
916  upperPad->SetRightMargin(0.05);
917  if(crit<3)upperPad->SetLogy();
918  inhisto->SetTitleOffset(0.6,"Y");
919  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
920  inhisto->GetYaxis()->SetTitleOffset(0.7);
921  inhisto->SetLineColor(kBlue+1);
922  inhisto->Draw();
923 
924  lowerPadRight->cd();
925  lowerPadRight->SetLeftMargin(0.09);
926  lowerPadRight->SetRightMargin(0.12);
927  plot2D->GetYaxis()->SetTitleOffset(1.3);
928  plot2D->Draw("colz");
929 
930  lowerPadLeft->cd();
931  lowerPadLeft->SetLeftMargin(0.09);
932  lowerPadLeft->SetRightMargin(0.07);
933  lowerPadLeft->SetLogy();
934  distrib->SetLineColor(kBlue+1);
935  distrib->GetYaxis()->SetTitleOffset(1.3);
936  distrib->Draw();
937  distrib_wTRDStruc->SetLineColor(kGreen+1);
938  distrib_wTRDStruc->DrawCopy("same");
939  distrib_woTRDStruc->SetLineColor(kMagenta+1);
940  distrib_woTRDStruc->DrawCopy("same");
941 
942  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
943  //. . .fit histogram
944  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
945 
946  //..to exclude first 2 bins from max finding
947  distrib->GetXaxis()->SetRangeUser(distrib->GetBinCenter(2),distrib->GetBinCenter(dnbins));
948  Double_t maxDistr = distrib->GetMaximum();
949  Int_t maxBin = distrib->GetMaximumBin();
950  Double_t maxCenter = distrib->GetBinCenter(maxBin);
951 
952  //..good range is around the max value as long as the
953  //..bin content is larger than 2 entries
954  for(Int_t i = maxBin ; i<=dnbins ; i++)
955  {
956  if(distrib->GetBinContent(i)<2) break ;
957  goodmax = distrib->GetBinCenter(i);
958  }
959  for(Int_t i = maxBin ; i>2 ; i--)
960  {
961  if(distrib->GetBinContent(i)<2) break ;
962  goodmin = distrib->GetBinLowEdge(i);
963  }
964  //if(maxBin<2)goodmin = distrib->GetBinLowEdge(2);
965 
966  //..Define min/max range of the fit function
967  Double_t minFitRange=goodmin;
968  Double_t maxFitRange=goodmax;
969  //if(crit==2)minFitRange=distrib->GetBinLowEdge(2);//..exclude first 2 bins
970  //if(crit==2)maxFitRange=dmaxVal;
971  if(crit==3)minFitRange=-20;
972  if(crit==3)maxFitRange=20;
973 
974  //cout<<"max bin: "<<maxBin<<", max value: "<<maxDistr<<endl;
975  //cout<<"start mean: "<<maxCenter<<", mean range: 0-"<<dmaxVal<<endl;
976  //cout<<"fit range: "<<minFitRange<<" - "<<maxFitRange<<endl;
977  TF1 *fit2 = new TF1("fit2", "gaus",minFitRange,maxFitRange);
978  //..start the fit with a mean of the highest value
979  fit2->SetParLimits(0,0,maxDistr); //..do not allow negative ampl values
980  fit2->SetParameter(1,maxCenter); //..set the start value for the mean
981  fit2->SetParLimits(1,0,dmaxVal); //..do not allow negative mean values
982 
983  //..ELI - TO DO: eventually fit with TRD and without TRD seperatley
984  distrib->Fit(fit2, "0LQEM", "", minFitRange, maxFitRange);
985  Double_t sig, mean;// chi2ndf;
986  mean = fit2->GetParameter(1);
987  sig = fit2->GetParameter(2);
988 
989  //..for case 1 and 2 lower than 0 is an unphysical value
990  if(crit<3 && mean <0.) mean=0.;
991 
992  goodmin = mean - nsigma*sig ;
993  goodmax = mean + nsigma*sig ;
994  //..for case 1 and 2 lower than 0 is an unphysical value
995  if(crit<3 && goodmin <0.) goodmin=0.;
996  if(inputBins==-1) goodmin=-1; //..this is a special case for the very last histogram 3-40 GeV
997 
998  if(fPrint==1)cout<<" o Result of fit: "<<endl;
999  if(fPrint==1)cout<<" o "<<endl;
1000  if(fPrint==1)cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
1001  if(fPrint==1)cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
1002 
1003  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1004  //. . .Add info to histogram
1005  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1006  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1007  lline->SetLineColor(kGreen+2);
1008  lline->SetLineStyle(7);
1009  lline->Draw();
1010 
1011  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1012  rline->SetLineColor(kGreen+2);
1013  rline->SetLineStyle(7);
1014  rline->Draw();
1015 
1016  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1017  leg->AddEntry(lline,"Good region boundary","l");
1018  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1019  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1020  leg->SetBorderSize(0);
1021  leg->Draw("same");
1022 
1023  fit2->SetLineColor(kOrange-3);
1024  fit2->SetLineStyle(1);//7
1025  fit2->Draw("same");
1026 
1027  TLatex* text = 0x0;
1028  if(crit==1) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1029  if(crit==2) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1030  if(crit==3) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1031  text->SetTextSize(0.06);
1032  text->SetNDC();
1033  text->SetTextColor(1);
1034  text->Draw();
1035 
1036  upperPad->cd();
1037  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
1038  uline->SetLineColor(kGreen+2);
1039  uline->SetLineStyle(7);
1040  uline->Draw();
1041  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
1042  lowline->SetLineColor(kGreen+2);
1043  lowline->SetLineStyle(7);
1044  lowline->Draw();
1045  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1046  //. . .Save histogram
1047  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1048  c1->Update();
1049  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
1050  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1051  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1052  if(crit==3)name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1053  c1->SaveAs(name);
1054 
1055  fRootFile->WriteObject(c1,c1->GetName());
1056  fRootFile->WriteObject(plot2D,plot2D->GetName());
1057  fRootFile->WriteObject(distrib,distrib->GetName());
1058  fRootFile->WriteObject(inhisto,inhisto->GetName());
1059  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1060  //. . . Mark the bad cells in the fFlag array
1061  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1062  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1063  //. . .(0 by default - good cell)
1064  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1065  if(fPrint==1)cout<<" o Flag bad cells that are outside the good range "<<endl;
1066  for(Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1067  {
1068  //..cell=0 and bin=1, cell=1 and bin=2
1069  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1070  if(inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)
1071  {
1072  //cout<<"smaller than min value: "<<inhisto->GetBinContent(cell+1)<<endl;
1073  fFlag[cell]=fCriterionCounter;
1074  }
1075  if(inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
1076  {
1077  //cout<<"bigger than max value: "<<inhisto->GetBinContent(cell+1)<<endl;
1078  fFlag[cell]=fCriterionCounter;
1079  }
1080  }
1081  if(fPrint==1)cout<<" o o o o o o o o o o o o o o o o o o o o o o o"<<endl;
1082 }
1083 
1088 //________________________________________________________________________
1090 {
1091  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1092  //.. RESULTS
1093  //.. 1) Print the bad cells
1094  //.. and write the results to a file
1095  //.. for each added period analysis
1096  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1097  TArrayD periodArray;
1098  Double_t emin,emax,sig;
1099  Int_t criterion;
1100  TString output;
1101  Int_t nb1=0;
1102 
1103  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1104  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1105  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1106  TString aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1107  ofstream file2(aliceTwikiTable, ios::out | ios::trunc);
1108  if(file2)
1109  {
1110  file2<<"|*Criterium* | *E_min !GeV* | *E_max !GeV* | *sigma* | *Excluded Cells* |"<<endl;
1111  }
1112  file2.close();
1113 
1114  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
1115  {
1116  periodArray=fAnalysisVector.at(i);
1117  criterion =periodArray.At(0);
1118  emin =periodArray.At(2);
1119  emax =periodArray.At(3);
1120  sig =periodArray.At(1);
1121 
1122  //..Print the results on the screen and
1123  //..write the results in a file
1124  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
1125  ofstream file(output, ios::out | ios::trunc);
1126  if(!file)
1127  {
1128  cout<<"#### Major Error. Check the textfile!"<<endl;
1129  }
1130  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1131  if(fPrint==1)cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1132 
1133  nb1=0;
1134  for(Int_t cellID=fStartCell ;cellID<fNoOfCells;cellID++)
1135  {
1136  if(fFlag[cellID]==(i+2))
1137  {
1138  nb1++;
1139  file<<cellID<<", ";
1140  }
1141  }
1142  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1143  file<<"("<<nb1<<")"<<endl;
1144  file.close();
1145  if(fPrint==1)cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1146  if(fPrint==1)cout<<endl;
1147  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1148  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1149  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1150  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1151  if(file2)
1152  {
1153  file2<<"| "<<criterion<<" | "<<emin<<" | "<<emax<<" | "<<sig<<" | "<<nb1<<" |"<<endl;
1154  }
1155  file2.close();
1156  }
1157 }
1164 //________________________________________________________________________
1166 {
1167  Int_t cellID, nDeadDCalCells = 0, nDeadEMCalCells = 0, nDCalCells = 0, nEMCalCells = 0;
1168  Double_t perDeadEMCal,perDeadDCal,perBadEMCal,perBadDCal,perWarmEMCal,perWarmDCal;
1169  TString aliceTwikiTable, cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells,goodCellsRatio,cellProp;
1170  TH2F* cellAmp_masked = (TH2F*)fCellAmplitude->Clone("hcellAmp_masked");
1171  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1172 
1173  deadPdfName = Form("%s/%s/%s_Dead_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1174  badPdfName = Form("%s/%s/%s_Bad_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1175  ratioOfBad = Form("%s/%s/%s_Bad_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1176  goodCells = Form("%s/%s/%s_Good_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1177  goodCellsRatio = Form("%s/%s/%s_Good_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1178  cellSummaryFile = Form("%s/%s/%s_%s_Bad_Ampl_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1179  aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1180  cellProp = Form("%s/%s/%s_CellProp_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1181 
1182  cout<<" o Final results o "<<endl;
1183  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1184  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1185  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1186  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1187  {
1188  //..Direction of amplitude (Checks energies from 0-nBins GeV)
1189  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1190  {
1191  if(fFlag[cell]!=0)
1192  {
1193  //..cellID+1 = histogram bin
1194  cellAmp_masked->SetBinContent(amp,cell+1,0);
1195  }
1196  }
1197  //..Direction of time (Checks times from -275-975 GeV)
1198  for (Int_t time = 1; time <= fCellTime->GetNbinsX(); time++)
1199  {
1200  if(fFlag[cell]!=0)
1201  {
1202  //..cellID+1 = histogram bin
1203  cellTime_masked->SetBinContent(time,cell+1,0);
1204  }
1205  }
1206  }
1207  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1208  //..Plot some summary canvases
1209  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1210  TCanvas *c1 = new TCanvas("CellProp","I summary of cell properties",1000,1000);
1211  c1->ToggleEventStatus();
1212  c1->Divide(2,2);
1213  c1->cd(1)->SetLogz();
1214  //lowerPadRight->SetLeftMargin(0.09);
1215  //lowerPadRight->SetRightMargin(0.06);
1216  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1217  fCellAmplitude->SetYTitle("Abs. Cell Id");
1218  fCellAmplitude->GetYaxis()->SetTitleOffset(1.6);
1219  fCellAmplitude->Draw("colz");
1220  c1->cd(2)->SetLogz();
1221  fCellTime->SetXTitle("Cell Time [ns]");
1222  fCellTime->SetYTitle("Abs. Cell Id");
1223  fCellTime->GetYaxis()->SetTitleOffset(1.6);
1224  fCellTime->Draw("colz");
1225  c1->cd(3)->SetLogz();
1226  //lowerPadRight->SetLeftMargin(0.09);
1227  //lowerPadRight->SetRightMargin(0.06);
1228  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1229  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1230  cellAmp_masked->SetYTitle("Abs. Cell Id");
1231  cellAmp_masked->GetYaxis()->SetTitleOffset(1.6);
1232  cellAmp_masked->Draw("colz");
1233  c1->cd(4)->SetLogz();
1234  cellTime_masked->SetTitle("Masked Cell Time");
1235  cellTime_masked->SetXTitle("Cell Time [ns]");
1236  cellTime_masked->SetYTitle("Abs. Cell Id");
1237  cellTime_masked->GetYaxis()->SetTitleOffset(1.6);
1238  cellTime_masked->Draw("colz");
1239  c1->Update();
1240 
1241  TCanvas *c1_ratio = new TCanvas("CellPropRatio","II summary of cell properties ratio",1000,500);
1242  c1_ratio->ToggleEventStatus();
1243  c1_ratio->Divide(2);
1244  c1_ratio->cd(1)->SetLogz();
1245  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1246  cellAmp_masked->GetZaxis()->SetRangeUser(0.0001,10e7);
1247  cellAmp_masked->Draw("colz");
1248  c1_ratio->cd(2)->SetLogz();
1249 
1250  TH1D *hRefDistr = BuildMeanFromGood();
1251  TH2F* ratio2DAmp =(TH2F*)cellAmp_masked->Clone("ratio2DAmp");
1252  TH2F* Sum2DIdeal =(TH2F*)cellAmp_masked->Clone("Sum2DIdeal");
1253  Sum2DIdeal->Reset();
1254  //..Create an ideal 2D energy distribution for division.
1255  //..Helps to identify whether there are some cells that still
1256  //..need to be masked by hand
1257  for(Int_t eBin=0;eBin<Sum2DIdeal->GetNbinsX();eBin++)
1258  {
1259  Double_t binVal=hRefDistr->GetBinContent(eBin+1);
1260  for(Int_t icell=0;icell<Sum2DIdeal->GetNbinsY();icell++)
1261  {
1262  Sum2DIdeal->SetBinContent(eBin+1,icell+1,binVal);
1263  }
1264  }
1265  ratio2DAmp->SetTitle("Ratio of cell Amplitude to mean cell ampl.");
1266  ratio2DAmp->Divide(Sum2DIdeal);
1267  ratio2DAmp->GetZaxis()->UnZoom();
1268  ratio2DAmp->Draw("colz");
1269 
1270  TCanvas *c1_proj = new TCanvas("CellPropPProj","III summary of cell properties",1000,500);
1271  c1_proj->ToggleEventStatus();
1272  c1_proj->Divide(2);
1273  c1_proj->cd(1)->SetLogy();
1274  TH1* projEnergyMask = cellAmp_masked->ProjectionX(Form("%sMask_Proj",cellAmp_masked->GetName()),fStartCell,fNoOfCells);
1275  projEnergyMask->SetXTitle("Cell Energy [GeV]");
1276  projEnergyMask->GetYaxis()->SetTitleOffset(1.6);
1277  projEnergyMask->SetLineColor(kTeal+3);
1278  projEnergyMask->DrawCopy(" hist");
1279 
1280  TH1* projEnergy = fCellAmplitude->ProjectionX(Form("%s_Proj",fCellAmplitude->GetName()),fStartCell,fNoOfCells);
1281  projEnergy->DrawCopy("same hist");
1282 
1283  c1_proj->cd(2)->SetLogy();
1284  TH1* projTimeMask = cellTime_masked->ProjectionX(Form("%s_Proj",cellTime_masked->GetName()),fStartCell,fNoOfCells);
1285  projTimeMask->SetXTitle("Cell Time [ns]");
1286  projTimeMask->GetYaxis()->SetTitleOffset(1.6);
1287  projTimeMask->SetLineColor(kGreen+3);
1288  projTimeMask->DrawCopy("hist");
1289  TH1* projTime = fCellTime->ProjectionX(Form("%s_Proj",fCellTime->GetName()),fStartCell,fNoOfCells);
1290  projTime->DrawCopy("same hist");
1291  c1_proj->Update();
1292 
1293  //..save to a PDF
1294  c1 ->Print(Form("%s(",cellProp.Data()));
1295  c1_ratio ->Print(Form("%s",cellProp.Data()));
1296  c1_proj ->Print(Form("%s)",cellProp.Data()));
1297  //..Scale the histogtams by the number of events
1298  //..so that they are more comparable for a run-by-run
1299  //..analysis
1300  Double_t totalevents = fProcessedEvents->Integral();
1301  fCellAmplitude ->Scale(1.0/totalevents);
1302  cellAmp_masked ->Scale(1.0/totalevents);
1303  fCellTime ->Scale(1.0/totalevents);
1304 
1305  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1306  //..Write the final results of dead and bad cells in a file and on screen
1307  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1308  ofstream file(cellSummaryFile, ios::out | ios::trunc);
1309  if(file)
1310  {
1311  file<<"Dead cells : "<<endl;
1312  cout<<" o Dead cells : "<<endl;
1313  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
1314  {
1315  if(cellID==0)
1316  {
1317  file<<"In EMCal: "<<endl;
1318  }
1319  if(cellID==fCellStartDCal)
1320  {
1321  file<<"\n"<<endl;
1322  file<<"In DCal: "<<endl;
1323  }
1324  if(fFlag[cellID]==1)
1325  {
1326  file<<cellID<<", ";
1327  if(cellID<fCellStartDCal)nDeadEMCalCells++;
1328  else nDeadDCalCells++;
1329  }
1330  }
1331  file<<"\n"<<endl;
1332  perDeadEMCal=100*nDeadEMCalCells/(1.0*fCellStartDCal);
1333  perDeadDCal=100*nDeadDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1334  file<<"EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1335  cout<<" o EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1336 
1337  file<<"Bad cells: "<<endl;
1338  cout<<" o Bad cells: "<<endl;
1339  for(cellID=fStartCell;cellID<fNoOfCells;cellID++)
1340  {
1341  if(cellID==0)
1342  {
1343  file<<"In EMCal: "<<endl;
1344  }
1345  if(cellID==fCellStartDCal)
1346  {
1347  file<<"\n"<<endl;
1348  file<<"In DCal: "<<endl;
1349  }
1350  if(fFlag[cellID]>1)
1351  {
1352  file<<cellID<<", ";
1353  if(cellID<fCellStartDCal)nEMCalCells++;
1354  else nDCalCells++;
1355  }
1356  }
1357  file<<"\n"<<endl;
1358  perBadEMCal=100*nEMCalCells/(1.0*fCellStartDCal);
1359  perBadDCal =100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1360  file<<"EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1361  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1362  }
1363  file.close();
1364 
1365  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1366  //..Determine the number of warm cells and save the flagged cells to .pdf files
1367  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1368  if(fPrint==1)cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1369  SaveBadCellsToPDF(1,badPdfName) ;
1370  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
1371  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1372  if(fTestRoutine==1)SaveBadCellsToPDF(20,goodCellsRatio) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1373 
1374  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1375  //..Fill the histograms with the flag information
1376  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1377  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1378  {
1379  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
1380  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
1381  }
1382  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
1383  c2->ToggleEventStatus();
1384  c2->Divide(1,2);
1385  c2->cd(1);
1386  fhCellFlag->SetTitle("cell flag");
1387  fhCellFlag->SetXTitle("Abs. Cell Id");
1388  fhCellFlag->SetYTitle("flag by which cell was excluded");
1389  fhCellFlag->DrawCopy("hist");
1390  c2->cd(2);
1391  fhCellWarm->SetTitle("Warm cells");
1392  fhCellWarm->SetXTitle("Abs. Cell Id");
1393  fhCellWarm->SetYTitle("warm=1");
1394  fhCellWarm->DrawCopy("hist");
1395  c2->Update();
1396 
1397  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1398  //..Plot the 2D distribution of cells by flag
1399  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1400  PlotFlaggedCells2D(0); //..all good cells
1401  PlotFlaggedCells2D(1); //..all dead cells
1402  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
1403  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
1404 
1405  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1406  //..Add different histograms/canvases to the output root file
1407  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1408  TString name1,name2,name3;
1409  name1 = Form("%s/%s/CellPropertiesRatio.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1410  c1_ratio->SaveAs(name1);
1411  name2 = Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1412  c1->SaveAs(name2);
1413  name3 = Form("%s/%s/CellPropertiesProj.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1414  c1_proj->SaveAs(name3);
1415 
1416  fRootFile->WriteObject(c1,c1->GetName());
1417  fRootFile->WriteObject(c1_ratio,c1_ratio->GetName());
1418  fRootFile->WriteObject(c1_proj,c1_proj->GetName());
1419  fRootFile->WriteObject(c2,c2->GetName());
1420  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
1421  fRootFile->WriteObject(cellAmp_masked,cellAmp_masked->GetName());
1422  fRootFile->WriteObject(ratio2DAmp,ratio2DAmp->GetName());
1423  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
1424  fRootFile->WriteObject(fProcessedEvents,fProcessedEvents->GetName());
1425  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
1426  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
1427  fRootFile->WriteObject(projEnergyMask,projEnergyMask->GetName());
1428  fRootFile->WriteObject(projEnergy,projEnergy->GetName());
1429  //..Save all amplitudes to the root file
1430  SaveHistoToFile();
1431 
1432  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1433  //..Save also the identified warm channels into a text file.
1434  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1435  Int_t nWEMCalCells =0;
1436  Int_t nWDCalCells =0;
1437  file.open(cellSummaryFile, ios::out | ios::app);
1438  if(file)
1439  {
1440  file<<"Warm cells : "<<endl;
1441  if(fPrint==1)cout<<" o Warm cells : "<<endl;
1442  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
1443  {
1444  if(cellID==0)
1445  {
1446  file<<"In EMCal: "<<endl;
1447  }
1448  if(cellID==fCellStartDCal)
1449  {
1450  file<<"\n"<<endl;
1451  file<<"In DCal: "<<endl;
1452  }
1453  if(fWarmCell[cellID]==1)
1454  {
1455  file<<cellID<<", ";
1456  if(cellID<fCellStartDCal)nWEMCalCells++;
1457  else nWDCalCells++;
1458  }
1459  }
1460  file<<"\n"<<endl;
1461  perWarmEMCal= 100*nWEMCalCells/(1.0*fCellStartDCal);
1462  perWarmDCal = 100*nWDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1463  file<<"EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1464  if(fPrint==1)cout<<" o EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1465  }
1466  file.close();
1467  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1468  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1469  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1470  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1471  if(file2)
1472  {
1473  file2<<"1=energy/hit, 2= hit/event"<<endl;
1474  file2<<"\n"<<endl;
1475  file2<<"| *Detector* | *No of cells* | *percentage* |"<<endl;
1476  file2<<"| Dead EMCal | "<<nDeadEMCalCells<<" | "<<perDeadEMCal<<"% |"<<endl;
1477  file2<<"| Bad EMCal | "<<nEMCalCells<<" | "<<perBadEMCal<<"% |"<<endl;
1478  file2<<"| - Warm EMCal | "<<nWEMCalCells<<" | "<<perWarmEMCal<<"% |"<<endl;
1479  file2<<"| Dead DCal | "<<nDeadDCalCells<<" | "<<perDeadDCal<<"% |"<<endl;
1480  file2<<"| Bad DCal | "<<nDCalCells<<" | "<<perBadDCal<<"% |"<<endl;
1481  file2<<"| - Warm DCal | "<<nWDCalCells<<" | "<<perWarmDCal<<"% |"<<endl;
1482  file2<<"\n"<<endl;
1483  }
1484  file2.close();
1485 
1486 }
1487 
1503 //________________________________________________________________________
1505 {
1506  gROOT->SetStyle("Plain");
1507  gStyle->SetOptStat(0);
1508  gStyle->SetFillColor(kWhite);
1509  gStyle->SetTitleFillColor(kWhite);
1510  gStyle->SetPalette(1);
1511 
1512  char title[100];
1513  char name[100];
1514 
1515  TH1D *hRefDistr = BuildMeanFromGood();
1516  fRootFile->WriteObject(hRefDistr,hRefDistr->GetName());
1517  Int_t firstCanvas=0;
1518  Bool_t candidate;
1519  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
1520  text->SetTextSize(0.07);
1521  text->SetTextColor(kOrange-3);
1522  text->SetNDC();
1523 
1524  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
1525  text2->SetTextSize(0.08);
1526  text2->SetTextColor(8);
1527  text2->SetNDC();
1528 
1529  TLatex* textA = new TLatex(0.65,0.62,"*test*");
1530  textA->SetTextSize(0.04);
1531  textA->SetTextColor(1);
1532  textA->SetNDC();
1533 
1534  //..collect cells in an internal vector.
1535  //..when the vector is of size=9 or at the end of being filled
1536  //..plot the channels into a canvas
1537  std::vector<Int_t> channelVector;
1538  channelVector.clear();
1539  cout<<" o Start printing into .pdf for version: "<<version<<endl;
1540  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
1541  {
1542  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
1543  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
1544  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
1545  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
1546  if(fFlag[cell]==0 && version==20)channelVector.push_back(cell);
1547 
1548  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
1549  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1550  if(channelVector.size()==9 || cell == fNoOfCells-1)
1551  {
1552  cout<<"."<<flush;
1553 
1554  TString internal_pdfName=pdfName;
1555  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1556  if(channelVector.size() > 6) c1->Divide(3,3);
1557  else if (channelVector.size() > 3) c1->Divide(3,2);
1558  else c1->Divide(3,1);
1559 
1560  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1561  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
1562  {
1563  sprintf(name, "Cell %d",channelVector.at(i)) ;
1564  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
1565  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
1566 
1567  c1->cd(i%9 + 1);
1568  c1->cd(i%9 + 1)->SetLogy();
1569  if(i==0)
1570  {
1571  leg->AddEntry(hRefDistr,"mean of good","l");
1572  leg->AddEntry(hCell,"current","l");
1573  }
1574  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
1575  candidate = CheckDistribution(hCell,hRefDistr);
1576  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
1577  if(version>2)//..These are ratio plots of energy distr. of cell and mean of all good cells
1578  {
1579  hCell->Divide(hRefDistr);
1580  }
1581  //.. save histograms to file
1582  if(version==1) fOutputListBad->Add(hCell);
1583  if(version==10)fOutputListBadRatio->Add(hCell);
1584  if(version==2) fOutputListGood->Add(hCell);
1585  if(version==20)fOutputListGoodRatio->Add(hCell);
1586 
1587  hCell->SetLineColor(kBlue+1);
1588  hCell->GetXaxis()->SetTitle("E (GeV)");
1589  hCell->GetYaxis()->SetTitle("N Entries");
1590  hCell->GetXaxis()->SetRangeUser(0.,10.);
1591  hCell->SetLineWidth(1) ;
1592  hCell->SetTitle(title);
1593  hRefDistr->SetLineColor(kGray+2);
1594  hRefDistr->SetLineWidth(1);
1595 
1596  hCell->Draw("hist");
1597 
1598  if(version==1 || version==2)hRefDistr->Draw("same") ;
1599 
1600  //..Mark the histogram that could be miscalibrated and labelled as warm
1601  if(candidate==1 && (version==1 || version==10))
1602  {
1603  gPad->SetFrameFillColor(kYellow-10);
1604  text->Draw();
1605  }
1606  if(version==1)
1607  {
1608  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
1609  textA->DrawLatex(0.65,0.62,Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
1610  }
1611  if(candidate==0 && (version==2 || version==20))
1612  {
1613  gPad->SetFrameFillColor(kYellow-10);
1614  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
1615  leg->Draw();
1616  }
1617  if(version<2)leg->Draw();
1618  }
1619 
1620  if(channelVector.size()<9 || cell == fNoOfCells-1)
1621  {
1622  internal_pdfName +=")";
1623  c1->Print(internal_pdfName.Data());
1624  }
1625  else if(firstCanvas==0)
1626  {
1627  internal_pdfName +="(";
1628  c1->Print(internal_pdfName.Data());
1629  firstCanvas=1;
1630  }
1631  else
1632  {
1633  c1->Print(internal_pdfName.Data());
1634  }
1635  delete c1;
1636  delete leg;
1637  channelVector.clear();
1638  }
1639  }
1640  cout<<endl;
1641  delete hRefDistr;
1642  //..Add the subdirectories to the file
1643  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
1644  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
1645  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
1646  if(version==20)fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1647 
1648  if(fPrint==1)cout<<endl;
1649 }
1653 //_________________________________________________________________________
1655 {
1656  TH1D* hGoodAmp;
1657  TH1D* hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean");
1658  hgoodMean->Reset();
1659  Int_t NrGood=0;
1660 
1661  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1662  {
1663  if(warmIn==0 && fFlag[cell]!=0 )continue;
1664  if(warmIn==1 && fFlag[cell]!=0 && fWarmCell[cell]==0)continue;
1665  if(warmIn==2 && fWarmCell[cell]==0)continue;
1666  NrGood++;
1667 
1668  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1669  hgoodMean->Add(hGoodAmp);
1670  }
1671  hgoodMean->Scale(1.0/NrGood);
1672 
1673  return hgoodMean;
1674 }
1685 //_________________________________________________________________________
1687 {
1688  TString Name = Form("%s_ratio",histogram->GetName());
1689  TH1* ratio=(TH1*)histogram->Clone(Name);
1690  ratio->Divide(reference);
1691 
1692  Double_t percentageOfLast=0.6;
1693  Double_t higherThanMean=5;
1694  Double_t highestRatio=1000;
1695  Double_t maxMean=10;
1696  Double_t minMean=0.1;
1697  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
1698 
1699  //..By default each cell is a candidate for recalibration
1700  Bool_t candidate=1;
1701  //..Find bin where reference has value 1, and the corresponding x-value
1702  Double_t totalevents = fProcessedEvents->Integral();
1703  Int_t binHeightOne = reference->FindLastBinAbove(1.0/totalevents);//..look at the spectrum until there is 1hit/event GeV
1704  Double_t binCentreHeightOne = reference->GetBinCenter(binHeightOne);
1705  Double_t thirdBinCentre = reference->GetBinCenter(3);
1706 
1707  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1708  //..Check the histogram
1709  //..Different checks to see whether the
1710  //..cell is really bad. Set candidate to 0.
1711 
1712  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1713  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
1714  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
1715  {
1716  candidate=0;
1717  //cout<<"1"<<endl;
1718  }
1719 
1720  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1721  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
1722  //ELI - TO DO: check that crieteria carfully - seems to work but not shure about it
1723  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
1724  if(ratio->GetMaximum()>highestRatio)//
1725  {
1726  candidate=0;
1727  //cout<<"2"<<endl;
1728  }
1729 
1730  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1731  //..check whether the ratio is much larger than 1
1732  //..calculate the mean in the relevant energy range
1733  Double_t mean=0;
1734  Int_t nullEntries=0;
1735  for(Int_t i=2;i<binHeightOne;i++)
1736  {
1737  if(ratio->GetBinContent(i)!=0)mean+=ratio->GetBinContent(i);
1738  else nullEntries++;
1739  }
1740  mean*=1.0/(binHeightOne-1-nullEntries);//..divide by number of bins (excluding bins without entries)
1741  if(mean>maxMean || mean<minMean)
1742  {
1743  candidate=0;
1744  //cout<<"3"<<endl;
1745  }
1746 
1747  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1748  //..check whether there are large spikes in the histogram
1749  //..compare bin values to mean of the ratio. If there is a bin value with
1750  //..content "higherThanMean" times lareger than mean it's losing it candidate status
1751  mean=0;
1752  //..Find the maximum in the mean range (0-binHeightOne)
1753  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
1754  Double_t localMaxBin=ratio->GetMaximumBin();
1755 
1756  for(Int_t i=2;i<binHeightOne;i++)
1757  {
1758  //..Exclude 0 bins and exclude bins near the maximum
1759  if(ratio->GetBinContent(i)<=0) continue;
1760  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
1761  mean+=ratio->GetBinContent(i);
1762  }
1763  mean*=1.0/(binHeightOne-1);//..divide by number of bins
1764  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
1765  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1766  if(ratio->GetMaximum()>mean*higherThanMean)
1767  {
1768  candidate=0;
1769  //cout<<"4"<<endl;
1770  }
1771 
1772  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1773  //..Check for a cliff at 4GeV, happens in some cases
1774  Double_t beforeCliff=0;
1775  Double_t afterCliff=0;
1776  Int_t binsBefore=0;
1777  Int_t binsAfter=0;
1778  Int_t cliffBin = ratio->FindBin(4);
1779  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
1780  {
1781  //..Exclude 0 bins
1782  if(ratio->GetBinContent(i)<=0)continue;
1783  if(i<=cliffBin) binsBefore++;
1784  if(i>cliffBin) binsAfter++;
1785  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
1786  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
1787  beforeCliff*=1.0/binsBefore;
1788  afterCliff*=1.0/binsAfter;
1789  }
1790  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
1791  {
1792  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
1793  if(beforeCliff!=0 && afterCliff!=0)cout<<"5"<<endl;
1794  }
1795 
1796  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1797  //..Employ peak finder
1798 /* if(candidate==1)
1799  {
1800  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
1801  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
1802  //GetNPeaks();
1803  //cout<<"found N peaks: "<<nfound<<endl;
1804  }
1805 */
1806  return candidate;
1807 }
1808 
1817 //_________________________________________________________________________
1819 {
1820  //..TRD support structure
1821  //..(determined by eye, could be improved, but is already very acurate):
1822  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
1823  //..row 1 (21),22,23,24 45,46,47,(48) 69,70,71,(72) (92),93,94,95 117,118,(119) 127 149,150,151 (173),174,175,(176) 198,199,200
1824  Bool_t coveredByTRDSupportStruc=0;
1825 
1826  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
1827  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
1828  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
1829  )
1830  {
1831  coveredByTRDSupportStruc=1;
1832  }
1833  return coveredByTRDSupportStruc;
1834 }
1843 //_________________________________________________________________________
1845 {
1846  //..build two dimensional histogram with values row vs. column
1847  TString histoName;
1848  histoName = Form("2DChannelMap_Flag%d",flagBegin);
1849  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100");
1850 
1851  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
1852  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1853  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1854 
1855  Int_t cellColumn=0,cellRow=0;
1856  Int_t cellColumnAbs=0,cellRowAbs=0;
1857  Int_t trash;
1858 
1859  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1860  {
1861  //..Do that only for cell ids also accepted by the code
1862  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1863 
1864  //..Get Row and Collumn for cell ID c
1865  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1866  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1867  {
1868  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1869  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1870  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1871  }
1872  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1873  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1874  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
1875 
1876 
1877  }
1878  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
1879  c1->ToggleEventStatus();
1880  c1->cd();
1881  //lowerPadRight->SetLeftMargin(0.09);
1882  //lowerPadRight->SetRightMargin(0.06);
1883  plot2D->Draw("colz");
1884 
1885  TLatex* text = 0x0;
1886  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
1887  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
1888  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
1889  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
1890  text->SetTextSize(0.06);
1891  text->SetNDC();
1892  text->SetTextColor(1);
1893  text->Draw();
1894 
1895  c1->Update();
1896  TString name =Form("%s/%s/%s_%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data() , histoName.Data());
1897  c1->SaveAs(name);
1899  fRootFile->WriteObject(plot2D,plot2D->GetName());
1900 
1901 }
1905 //_________________________________________________________________________
1907 {
1908  char name[100];
1909  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
1910  {
1911  sprintf(name, "Cell %d",cell) ;
1912  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
1913  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
1914  }
1915  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1916 }
TH1F * BuildTimeMean(Int_t crit, Double_t tmin, Double_t tmax)
double Double_t
Definition: External.C:58
void SetRunNumber(Int_t run)
Definition: External.C:236
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...
TString fileName
TFile * fRootFile
! root file with all histograms from this analysis
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.
Int_t * fFlag
! fFlag[CellID] = 0 (ok),1 (dead),2 and higher (bad certain criteria) start at 0 (cellID 0 = histobin...
Int_t fNMaxCols
Maximum No of colums in module (eta direction)
Int_t fCriterionCounter
! This value will be written in fflag and updates after each PeriodAnalysis, to distinguish the steps...
TString fMergedFileName
Filename of the .root file containing the merged runs.
AliEMCALGeometry * GetEMCALGeometry() const
int Int_t
Definition: External.C:63
TList * fOutputListBadRatio
! list with bad channel amplitude ratios, stored in fRootFile
TH1F * BuildHitAndEnergyMean(Int_t crit, Double_t emin=0.1, Double_t emax=2.)
TString fAnalysisInput
Here the .root files of each run of the period are saved.
std::vector< Int_t > fManualMask
! Is a list of cells that should be addidionally masked by hand.
TString fQADirect
Dierctory in the QA.root files where the input histograms are stored.
Bool_t * fWarmCell
! fWarmCell[CellID] = 0 (really bad), fWarmCell[CellID] = 1 (candidate for warm), ...
Definition: External.C:212
TString fWorkdir
Directory which contains the folders fMergeOutput, fAnalysisInput and fAnalysisOutput. By default it is './'.
TList * fOutputListBad
! list with bad channel amplitudes, stored in fRootFile
Bool_t IsCoveredByTRD(Int_t row, Int_t collumn)
TH1F * fhCellFlag
! histogram that stores by which flag the cell has been excluded
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.
void FlagAsBad(Int_t crit, TH1F *inhisto, Double_t nsigma=4., Double_t dnbins=200)
const char * pdfName
Definition: DrawAnaELoss.C:30
void PlotFlaggedCells2D(Int_t flagBegin, Int_t flagEnd=-1)
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.
Bool_t fPrint
If set true more couts with information of the excluded cells will be printed.
TList * fOutputListGood
! list with good channel amplitudes, stored in fRootFile
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
TList with histograms for a given trigger.
TH1F * fhCellWarm
! histogram that stores whether the cell was marked as warm
Int_t fStartCell
ID of the first cell you want to check.
TList * fOutputListGoodRatio
! list with good channel amplitude ratios, stored in fRootFile
bool Bool_t
Definition: External.C:53
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 ...
TString fTrainNo
Train number of the analyszed data (can deduce pass & trigger from that etc.)
void SaveBadCellsToPDF(Int_t version, TString pdfName)
TH1D * BuildMeanFromGood(Int_t warmIn=0)
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...
Definition: External.C:196
TString fRunList
Thats the full path and name of the file which contains a list of all runs to be merged together...
TDirectoryFile * dir
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