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