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