AliPhysics  c923f52 (c923f52)
 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 
32 #include <Riostream.h>
33 #include <stdio.h>
34 #include <fstream>
35 #include <iostream>
36 #include <alloca.h>
37 #include <string>
38 #include <cstring>
39 
40 // --- ANALYSIS system ---
41 #include <AliCalorimeterUtils.h>
42 #include <AliEMCALGeometry.h>
44 #include <AliAODEvent.h>
45 
46 using std::vector;
47 using std::cout;
48 using std::endl;
49 using std::ofstream;
50 using std::flush;
51 using std::ios;
52 
56 
60 //________________________________________________________________________
62 TObject(),
63  fCurrentRunNumber(-1),
64  fPeriod(),
65  fPass(),
66  fTrigger(),
67  fNoOfCells(),
68  fCellStartDCal(12288),
69  fAnalysisOutput(),
70  fAnalysisInput(),
71  fRunList(),
72  fRunListFileName(),
73  fWorkdir(),
74  fQADirect(),
75  fMergedFileName(),
76  fAnalysisVector(),
77  fTrial(),
78  fExternalFileName(),
79  fTestRoutine(),
80  fNMaxCols(48),
81  fNMaxRows(24),
82  fNMaxColsAbs(),
83  fNMaxRowsAbs(),
84  fFlag(),
85  fCriterionCounter(),
86  fWarmCell(),
87  fCaloUtils(),
88  fRootFile(),
89  fCellAmplitude(),
90  fCellTime(),
91  fProcessedEvents(),
92  fhCellFlag(),
93  fhCellWarm()
94 {
95  fCurrentRunNumber = 254381;
96  fPeriod = "LHC16h";
97  fPass = "muon_caloLego";
98  fTrigger = "AnyINT";
99  fWorkdir = ".";
100  fRunListFileName = "runList.txt";
101  fTrial = 0;
102 
103  Init();
104 }
105 
109 //________________________________________________________________________
110 AliAnaCaloChannelAnalysis::AliAnaCaloChannelAnalysis(TString period, TString pass, TString trigger, Int_t runNumber,Int_t trial, TString workDir, TString listName):
111  TObject(),
112  fCurrentRunNumber(-1),
113  fPeriod(),
114  fPass(),
115  fTrigger(),
116  fNoOfCells(),
117  fCellStartDCal(12288),
118  fAnalysisOutput(),
119  fAnalysisInput(),
120  fRunList(),
121  fRunListFileName(),
122  fWorkdir(),
123  fQADirect(),
124  fMergedFileName(),
125  fAnalysisVector(),
126  fTrial(),
127  fExternalFileName(),
128  fTestRoutine(),
129  fNMaxCols(48),
130  fNMaxRows(24),
131  fNMaxColsAbs(),
132  fNMaxRowsAbs(),
133  fFlag(),
134  fCriterionCounter(),
135  fWarmCell(),
136  fCaloUtils(),
137  fRootFile(),
138  fCellAmplitude(),
139  fCellTime(),
140  fProcessedEvents(),
141  fhCellFlag(),
142  fhCellWarm()
143 {
144  fCurrentRunNumber = runNumber;
145  fPeriod = period;
146  fPass = pass; //ELI only for folder structure
147  fTrigger = trigger; //ELI important to select trigger in output file == different wagons in lego train
148  fWorkdir = workDir;
149  fRunListFileName = listName;
150  fTrial = trial;
151 
152  Init();
153 }
154 
158 //________________________________________________________________________
160 {
161  //......................................................
162  //..Default values - can be set by functions
164  fTestRoutine=0;
165 
166  //..Settings for the input/output structure (hard coded)
167  // TO BE CHANGED
168  fAnalysisInput =Form("AnalysisInput/%s",fPeriod.Data());
169  fAnalysisOutput =Form("AnalysisOutput/%s/Version%i",fPeriod.Data(),fTrial);
170 
171  //..Make output directory if it doesn't exist
172  //..first the period folder
173  gSystem->mkdir(Form("%s/AnalysisOutput/%s",fWorkdir.Data(),fPeriod.Data()));
174  //..then the version folder
175  gSystem->mkdir(Form("%s/%s",fWorkdir.Data(),fAnalysisOutput.Data()));
176 
177  fMergedFileName= Form("%s/%s/%s/MergedRuns_%s.root",fWorkdir.Data(),fAnalysisInput.Data(),fPass.Data(),fTrigger.Data());
178  fRunList = Form("%s/%s/%s/%s",fWorkdir.Data(), fAnalysisInput.Data(), fPass.Data(), fRunListFileName.Data());
179  fQADirect = Form("CaloQA_%s",fTrigger.Data());
180 
181  TString fileName = Form("%s/%s/%s%s_Histograms_V%i.root",fWorkdir.Data(),fAnalysisOutput.Data(), fPass.Data(),fTrigger.Data() ,fTrial);
182  fRootFile = new TFile(fileName,"recreate");
183  //.. make sure the vector is empty
184  fAnalysisVector.clear();
185 
186  //......................................................
187  //..Initialize EMCal/DCal geometry
189  //..Create a dummy event for the CaloUtils
190  AliAODEvent* aod = new AliAODEvent();
193 
194  fNoOfCells =fCaloUtils->GetEMCALGeometry()->GetNCells(); //..Very important number, never change after that point!
195  Int_t nModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
196  //..This is how the calorimeter looks like in the current period (defined by example run ID fCurrentRunNumber)
197  cout<<"Called geometry for run number: "<<fCurrentRunNumber<<endl;
198  cout<<"Number of supermod: "<<nModules<<endl;
199  cout<<"Number of cells: "<<fNoOfCells<<endl;
200  cout<<"Cell ID of first DCal cell: "<<fCellStartDCal<<endl;
201  //cout<<"Number of supermod utils: "<<fCaloUtils->GetNumberOfSuperModulesUsed()<<endl; //..will always be 22 unless set by hand
202 
203  //......................................................
204  //..Initialize flag array to store how the cell is categorized
205  //..In the histogram: bin 1= cellID 0, bin 2= cellID 1 etc
206  //..In the array: fFlag[cellID]= some information
207  fFlag = new Int_t[fNoOfCells];
208  fWarmCell = new Bool_t[fNoOfCells];
209  fFlag[fNoOfCells] = {0}; //..flagged as good by default
210  fWarmCell[fNoOfCells] = {0}; //..flagged as not warm by default
211  fCriterionCounter=2; //This value will be written in fflag and updates after each PeriodAnalysis
212  //......................................................
213  //..setings for the 2D histogram
214  //fNMaxCols = 48; //eta direction
215  //fNMaxRows = 24; //phi direction
217  fNMaxRowsAbs = Int_t (nModules/2)*fNMaxRows; //multiply by number of supermodules
218 
219  //......................................................
220  //..Create TLists to organize the outputfile
221  fOutputListBad = new TList();
222  fOutputListGood = new TList();
223  fOutputListBadRatio = new TList();
224  fOutputListGoodRatio = new TList();
225 
226  fOutputListBad ->SetName("BadCell_Amplitudes");
227  fOutputListGood ->SetName("GoodCell_Amplitudes");
228  fOutputListBadRatio ->SetName("BadCell_AmplitudeRatios");
229  fOutputListGoodRatio->SetName("GoodCell_AmplitudeRatios");
230 
231  //fOutputListGood ->SetOwner();//ELI instead of delete in destructor??
232  //fOutputListBadRatio ->SetOwner();
233  //fOutputListGoodRatio ->SetOwner();
234 
235  //......................................................
236  //..Create Histograms to store the flag in a root file
237  fhCellFlag = new TH1F("fhCellFlag","fhCellFlag",fNoOfCells+10,0,fNoOfCells+10); //..cellID+1 = histogram bin
238  fhCellWarm = new TH1F("fhCellWarm","fhCellWarm",fNoOfCells,0,fNoOfCells); //..cellID+1 = histogram bin
239 }
240 
254 //________________________________________________________________________
256 {
257 // cout<<"fired trigger class"<<AliAODEvent::GetFiredTriggerClasses()<<endl;
258 
259  if(fExternalFileName=="")
260  {
261  //..If no extrenal file is provided merge different runs together
262  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
263  cout<<". . .Start process by converting files. . . . . . . . . . . ."<<endl;
264  cout<<endl;
266  if(fMergedFileName.IsNull())
267  {
268  Printf("File not produced, exit");
269  return;
270  }
271  cout<<endl;
272  }
273  else
274  {
275  //..If extrenal file is provided load it
276  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
277  cout<<". . .Start process by loading external file. . . . . . . . . . ."<<endl;
278  fMergedFileName= Form("%s/%s/%s/%s",fWorkdir.Data(),fAnalysisInput.Data(),fPass.Data(),fExternalFileName.Data());
279  }
280  cout<<". . .Load inputfile with name: "<<fMergedFileName<<" . . . . . . . ."<<endl;
281 
282  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
283  //.. Read all the needed input for the Bad/Dead Channel analysis
284  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
285  TFile *mergedFileInput = new TFile(fMergedFileName);
286  if(!mergedFileInput->IsOpen())
287  {
288  Printf("Error! Input file not found, abort");
289  return;
290  }
291  fCellAmplitude = (TH2F*) mergedFileInput->Get("hCellAmplitude");
292  fCellTime = (TH2F*) mergedFileInput->Get("hCellTime");
293  fProcessedEvents = (TH1F*) mergedFileInput->Get("hNEventsProcessedPerRun");
294 
295  cout<<". . .Continue process by . . . . . . . . . . . ."<<endl;
296  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
297  //.. DEAD CELLS
298  //.. Flag dead cells with fFlag=1
299  //.. this excludes cells from analysis (will not appear in results)
300  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
301  cout<<"o o o Flag dead cells o o o"<<endl;
302  FlagAsDead();
303  cout<<endl;
304  cout<<endl;
305 
306  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
307  //.. BAD CELLS
308  //.. Flag dead cells with fFlag=2 and bigger
309  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
310  cout<<"o o o Flag bad cells o o o"<<endl;
311  BCAnalysis();
312 
313  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
314  //..In the end summarize results
315  //..in a .pdf and a .txt file
316  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
317  cout<<"o o o Write .txt for each period analyis with bad cells o o o"<<endl;
319 
320  cout<<"o o o Create summary documents for the entire analysis o o o"<<endl;
322  fRootFile->WriteObject(fFlag,"FlagArray");
323  fRootFile->Close();
324  cout<<endl;
325 
326  //..make a reccomendation about the used energy range to be investigated
327  //..and the binning
328  TH1D *hRefDistr = BuildMeanFromGood();
329  //..Find bin where reference has value 1, and the corresponding x-value
330  Int_t binHeihgtOne = hRefDistr->FindLastBinAbove(1);
331  Double_t binCentreHeightOne = hRefDistr->GetBinCenter(binHeihgtOne);
332  cout<<". . .Recomendation:"<<endl;
333  cout<<". . .With the current statistic on average a cell has 1 hit at "<<binCentreHeightOne<<" GeV"<<endl;
334  cout<<". . .so it makes no sense to select energy ranges >"<<binCentreHeightOne<<" as cells will be maked bad just due to the lack of statistic"<<endl;
335  cout<<". . .End of process . . . . . . . . . . . . . . . . . . . . ."<<endl;
336  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
337 }
338 
345 //________________________________________________________________________
347 {
348  cout<<"o o o Start conversion process o o o"<<endl;
349  cout<<"o o o period: " << fPeriod << ", pass: " << fPass << ", trigger: "<<fTrigger<< endl;
350 
351  //..Create histograms needed for adding all the files together
352  //ELI a little problematic to hard code properties of histograms??
353  TH2F *hCellAmplitude = new TH2F("hCellAmplitude","Cell Amplitude",200,0,10,23040,0,23040);
354  TH2F *hCellTime = new TH2F("hCellTime","Cell Time",250,-275,975,23040,0,23040);
355  TH1F *hNEventsProcessedPerRun = new TH1F("hNEventsProcessedPerRun","Number of processed events in analyzed runs",1,0,1);
356 
357  //..Open the text file with the run list numbers and run index
358  cout<<"o o o Open .txt file with run indices. Name = " << fRunList << endl;
359  FILE *pFile = fopen(fRunList.Data(), "r");
360  if(!pFile)
361  {
362  cout<<"couldn't open file!"<<endl;
363  return "";
364  }
365  Int_t nEntr;
366  Int_t nEntrTot=0;
367  Int_t q;
368  Int_t ncols;
369  Int_t nlines = 0 ;
370  Int_t runId[500] ;
371  while (1)
372  {
373  ncols = fscanf(pFile," %d ",&q);
374  if (ncols< 0) break;
375  runId[nlines]=q;
376  nlines++;
377  }
378  fclose(pFile);
379 
380 
381  //..Open the different .root files with help of the run numbers from the text file
382  const Int_t nRun = nlines ;
383  TString base;
384  TString infile;
385  TString singleRunFileName;
386 
387  cout<<"o o o Start merging process of " << nRun <<" files"<< endl;
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(), fPass.Data(), runId[i]);
392  /*if ((fPass=="cpass1_pass2")||(fPass=="cfPass1-2"))
393  {
394  if (fTrigger=="default")
395  {
396  infile = Form("%s_barrel.root",base.Data());
397  }
398  else
399  {
400  infile = Form("%s_outer.root",base.Data());
401  }
402  }
403  else
404  {*/
405  //..This is a run2 case
406  infile = Form("%s.root",base.Data()) ;
407  //}
408 
409  cout<<" o Open .root file with name: "<<infile<<endl;
410  TFile *f = TFile::Open(infile);
411 
412  //..Do some basic checks
413  if(!f)
414  {
415  cout<<"Couldn't open/find .root file: "<<infile<<endl;
416  continue;
417  }
418  TDirectoryFile *dir = (TDirectoryFile *)f->Get(fQADirect);
419  if(!dir)
420  {
421  cout<<"Couln't open directory file in .root file: "<<infile<<", directory: "<<fQADirect<<endl;
422  continue;
423  }
424  TList *outputList = (TList*)dir->Get(fQADirect);
425  if(!outputList)
426  {
427  cout << "Couln't get TList from directory file: "<<fQADirect<<endl;
428  continue;
429  }
430 
431  //ELI should one maybe clone the hAmpId histos eg to hCellAmplitude, then one does't need to hard code them.
432  TH2F *hAmpId;
433  TH2F *hTimeId;
434  TH1F *hNEvents;
435 
436  hAmpId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
437  if(!hAmpId)
438  {
439  Printf("hAmpId not found");
440  outputList->ls();
441  continue;
442  }
443  hTimeId =(TH2F *)outputList->FindObject("EMCAL_hTimeId");
444  if(!hTimeId)
445  {
446  Printf("hTimeId not found");
447  outputList->ls();
448  continue;
449  }
450  hNEvents =(TH1F *)outputList->FindObject("hNEvents");
451  if(!hNEvents)
452  {
453  Printf("hNEvents not found");
454  outputList->ls();
455  continue;
456  }
457  nEntr = (Int_t)hNEvents->GetEntries();
458 
459  //..does that mean do not merge small files?
460  if (nEntr<100)
461  {
462  cout <<" o File to small to be merged. Only N entries " << nEntr << endl;
463  continue ;
464  }
465  cout <<" o File with N entries " << nEntr<<" will be merged"<< endl;
466  nEntrTot+=nEntr;
467  hCellAmplitude->Add(hAmpId);
468  hCellTime->Add(hTimeId);
469  hNEventsProcessedPerRun->Add(hNEvents);
470 
471  //..Create copies of the original root files just with the bad channel QA
472  //..So that a run by run bad channel analysis can be performed more easily
473  singleRunFileName= Form("%s/%s/%s/%d_Filtered.root",fWorkdir.Data(),fAnalysisInput.Data(),fPass.Data(),runId[i]);
474  TFile *singleRunFile = TFile::Open(singleRunFileName,"recreate");
475  hNEventsProcessedPerRun->Write();
476  hCellAmplitude->Write();
477  hCellTime->Write();
478  singleRunFile->Close();
479 
480  outputList->Delete();
481  dir->Delete();
482  f->Close();
483  delete f;
484  }
485 
486  //.. Save the merged histograms
487  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
488  cout<<"o o o "<<nEntrTot<<" events were merged"<<endl;
489  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
490  hNEventsProcessedPerRun->Write();
491  hCellAmplitude->Write();
492  hCellTime->Write();
493  BCF->Close();
494  cout<<"o o o End conversion process o o o"<<endl;
495  return fMergedFileName;
496 }
497 
498 
504 //_________________________________________________________________________
506 {
507  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
508  //.. BAD CELLS
509  //.. Flag bad cells with fFlag= 2,3,4,5.. etc
510  //.. this excludes cells from subsequent analysis
511  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
512  TArrayD periodArray;
513  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
514  {
515  periodArray=fAnalysisVector.at(i);
516  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
518  cout<<""<<endl;
519  cout<<""<<endl;
520  }
521  cout<<"o o o End of bad channel analysis o o o"<<endl;
522 }
523 
532 //________________________________________________________________________
534 {
535  TArrayD periodArray;
536  periodArray.Set(4);
537  periodArray.AddAt((Double_t)criteria,0);
538  periodArray.AddAt(nsigma,1);
539  periodArray.AddAt(emin,2);
540  periodArray.AddAt(emax,3);
541  fAnalysisVector.push_back(periodArray);
542 }
558 //____________________________________________________________________
560 {
561  //ELI criterion should be between 1-4
562  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;
563  cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
564  cout<<"o o o Done in the energy range E "<<emin<<"-"<<emax<<endl;
565 
566  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
567  //.. ANALYSIS OF CELLS WITH ENTRIES
568  //.. Build average distributions and fit them
569  //.. Three tests for bad cells:
570  //.. 1) Average energy per hit
571  //.. 2) Average hit per event
572  //.. 3) ...
573  //.. 4) ...
574  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
575  TH1F* histogram;
576  if(criterion < 6)cout<<"o o o Analyze average cell distributions o o o"<<endl;
577  //..For case 1 or 2
578  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax,nsigma);
579  //..For case 3, 4 or 5
580  //else if (criterion < 6)
581 
582  //..light systems p+p p+Pb
583  Double_t range=0.001;
584  if(emin>0.49)range=0.0005;
585  if(emin>0.99)range=0.0001;
586  if(emin>1.99)range=0.00005;
587 
588  //..for Pb+Pb
589 /* Double_t range=0.1;
590  if(emin>0.49)range=0.05;
591  if(emin>0.99)range=0.01;
592  if(emin>1.99)range=0.0002;
593 */
594 
595  if(criterion==1) FlagAsBad(criterion, histogram, nsigma, 200,-1);
596  if(criterion==2) FlagAsBad(criterion, histogram, nsigma, 600,range);
597 
598  /*
599  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval3);
600  if(criterion==4) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval1);
601  if(criterion==5) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval2);
602  */
603 }
604 
614 //_________________________________________________________________________
616 {
617  cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
618  TH1F *histogram;
619  if(crit==1)histogram = new TH1F(Form("hCellEtoNtotal_E%.2f-%.2f",emin,emax),Form("Energy per hit, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
620  if(crit==2)histogram = new TH1F(Form("hCellNtotal_E%.2f-%.2f",emin,emax),Form("Number of hits per event, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
621  histogram->SetXTitle("Abs. Cell Id");
622  if(crit==1)histogram->SetYTitle("Energy per hit");
623  if(crit==2)histogram->SetYTitle("Number of hits per event");
624  histogram->GetXaxis()->SetNdivisions(510);
625  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
626 
627  //..here the average hit per event and the average energy per hit is caluclated for each cell.
628  for (Int_t cell = 0; cell < fNoOfCells; cell++)
629  {
630  Double_t Esum = 0;
631  Double_t Nsum = 0;
632 
633  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
634  {
635  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
636  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
637  //..investigate only cells that were not flagged as dead ore bad
638  if (E < emin || E > emax || fFlag[cell]!=0) continue;
639  Esum += E*N;
640  Nsum += N;
641  }
642  //..Set the values only for cells that are not yet marked as bad
643  if(fFlag[cell]==0)
644  {
645  if(totalevents> 0. && crit==2)histogram->SetBinContent(cell+1, Nsum/totalevents); //..number of hits per event
646  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/Nsum); //..average energy per hit
647  }
648  }
649  return histogram;
650 }
655 //_________________________________________________________________________
657 {
658  TH1F *Histogram = new TH1F(Form("hSomethingWithTime_E%.2f-%.2f",emin,emax),Form("I don't know, %.2f < E < %.2f GeV",emin,emax), 2,0,2);
659 
660  return Histogram;
661 }
662 
667 //_________________________________________________________________________
669 {
670  Int_t sumOfExcl=0;
671  //..Direction of cell ID
672  for (Int_t cell = 0; cell < fNoOfCells; cell++)
673  {
674  Double_t nSum = 0;
675  //..Direction of amplitude (Checks energies from 0-10 GeV)
676  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
677  {
678  //..cellID+1 = histogram bin
679  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
680  nSum += N;
681  }
682  //..If the amplitude in one cell is basically 0
683  //..mark the cell as excluded
684  if(nSum < 0.5 && fFlag[cell]==0)
685  {
686  //..Cell flagged as dead.
687  //..Flag only if it hasn't been flagged before
688  fFlag[cell] =1;
689  sumOfExcl++;
690  }
691  }
692  cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
693  cout<<" ("<<sumOfExcl<<")"<<endl;
694 }
695 
710 //_________________________________________________________________________
711 void AliAnaCaloChannelAnalysis::FlagAsBad(Int_t crit, TH1F* inhisto, Double_t nsigma, Int_t dnbins, Double_t dmaxval)
712 {
713  gStyle->SetOptStat(0); // MG modif
714  gStyle->SetOptFit(0); // MG modif
715 
716  if(crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
717  if(crit==2)cout<<" o Fit average hit per event distribution"<<endl;
718 
719  Int_t cellColumn=0,cellRow=0;
720  Int_t cellColumnAbs=0,cellRowAbs=0;
721  Int_t trash;
722 
723  TString histoName=inhisto->GetName();
724  Double_t goodmax= 0. ;
725  Double_t goodmin= 0. ;
726  if (dmaxval < 0.)
727  {
728  dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
729  if(crit==2 && dmaxval > 1) dmaxval =1. ;
730  }
731 
732  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
733  //. . .build the distribution of average values
734  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
735  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
736  //..Be aware that the bin width should never be larger than 1/totalevents because this
737  //..is the minimum differce between the cells. One hit more/event otherwise you will see ugly empty bins in the histogram
738  if(((dmaxval-inhisto->GetMinimum())/(dnbins*1.0))<1.0/totalevents)
739  {
740  dnbins=(dmaxval-inhisto->GetMinimum())*totalevents;
741  cout<<"Problem - Reset dnbins to new value:"<<dnbins<<endl;
742  }
743 
744  //..build histos
745  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
746  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
747  distrib->SetYTitle("Entries");
748  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
749  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
750 
751 
752  //..build two dimensional histogram with values row vs. column
753  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);
754  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
755  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
756 
757  for (Int_t cell = 0; cell < fNoOfCells; cell++)
758  {
759  //..Do that only for cell ids also accepted by the code
760  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
761 
762  //..Fill histograms only for cells that are not yet flagged as bad
763  if(fFlag[cell]==0)
764  {
765  //..fill the distribution of avarge cell values
766  distrib->Fill(inhisto->GetBinContent(cell+1));
767  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
768  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
769  //..Get Row and Collumn for cell ID
770  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
771  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
772  {
773  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
774  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
775  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
776  }
777  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
778  //..check TRD support structure
779  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
780  {
781  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
782  }
783  else
784  {
785  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
786  }
787  }
788  }
789 
790  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
791  //. . .draw histogram + distribution
792  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
793 
794  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
795  c1->ToggleEventStatus();
796  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
797  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
798  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
799  upperPad->Draw();
800  lowerPadLeft->Draw();
801  lowerPadRight->Draw();
802 
803  upperPad->cd();
804  upperPad->SetLeftMargin(0.05);
805  upperPad->SetRightMargin(0.05);
806  upperPad->SetLogy();
807  inhisto->SetTitleOffset(0.6,"Y");
808  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
809  inhisto->GetYaxis()->SetTitleOffset(0.7);
810  inhisto->SetLineColor(kBlue+1);
811  inhisto->Draw();
812 
813  lowerPadRight->cd();
814  lowerPadRight->SetLeftMargin(0.09);
815  lowerPadRight->SetRightMargin(0.12);
816  plot2D->GetYaxis()->SetTitleOffset(1.3);
817  plot2D->Draw("colz");
818 
819  lowerPadLeft->cd();
820  lowerPadLeft->SetLeftMargin(0.09);
821  lowerPadLeft->SetRightMargin(0.07);
822  lowerPadLeft->SetLogy();
823  distrib->SetLineColor(kBlue+1);
824  distrib->GetYaxis()->SetTitleOffset(1.3);
825  distrib->Draw();
826  distrib_wTRDStruc->SetLineColor(kGreen+1);
827  distrib_wTRDStruc->DrawCopy("same");
828  distrib_woTRDStruc->SetLineColor(kMagenta+1);
829  distrib_woTRDStruc->DrawCopy("same");
830 
831  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
832  //. . .fit histogram
833  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
834  Int_t higherbin=0,i;
835  for(i = 2; i <= dnbins; i++)
836  {
837  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i)) higherbin = i ;
838  }
839  //..good range is around the max value as long as the
840  //..bin content is larger than 2 entries
841  for(i = higherbin ; i<=dnbins ; i++)
842  {
843  if(distrib->GetBinContent(i)<2) break ;
844  goodmax = distrib->GetBinCenter(i);
845  }
846  for(i = higherbin ; i>1 ; i--)
847  {
848  if(distrib->GetBinContent(i)<2) break ;
849  goodmin = distrib->GetBinLowEdge(i);
850  }
851  //cout<<"higherbin : "<<higherbin<<endl;
852  //cout<<"good range : "<<goodmin<<" - "<<goodmax<<endl;
853 
854  TF1 *fit2 = new TF1("fit2", "gaus",0,10);
855  //..start the fit with a mean of the highest value
856  fit2->SetParameter(1,higherbin);
857 
858  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
859  Double_t sig, mean;// chi2ndf;
860  // Marie midif to take into account very non gaussian distrig
861  mean = fit2->GetParameter(1);
862  sig = fit2->GetParameter(2);
863  //chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
864 
865  if (mean <0.) mean=0.; //ELI is this not a highly problematic case??
866 
867  goodmin = mean - nsigma*sig ;
868  goodmax = mean + nsigma*sig ;
869 
870  if (goodmin<0) goodmin=0.;
871 
872  cout<<" o Result of fit: "<<endl;
873  cout<<" o "<<endl;
874  cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
875  cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
876 
877  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
878  //. . .Add info to histogram
879  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
880  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
881  lline->SetLineColor(kGreen+2);
882  lline->SetLineStyle(7);
883  lline->Draw();
884 
885  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
886  rline->SetLineColor(kGreen+2);
887  rline->SetLineStyle(7);
888  rline->Draw();
889 
890  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
891  leg->AddEntry(lline,"Good region boundary","l");
892  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
893  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
894  leg->SetBorderSize(0);
895  leg->Draw("same");
896 
897  fit2->SetLineColor(kOrange-3);
898  fit2->SetLineStyle(1);//7
899  fit2->Draw("same");
900 
901  TLatex* text = 0x0;
902  if(crit==1) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
903  if(crit==2) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2fx10^-5",goodmin*100000,goodmax*100000));
904  text->SetTextSize(0.06);
905  text->SetNDC();
906  text->SetTextColor(1);
907  text->Draw();
908 
909 
910  upperPad->cd();
911  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
912  uline->SetLineColor(kGreen+2);
913  uline->SetLineStyle(7);
914  uline->Draw();
915  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
916  lowline->SetLineColor(kGreen+2);
917  lowline->SetLineStyle(7);
918  lowline->Draw();
919  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
920  //. . .Save histogram
921  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
922  c1->Update();
923  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
924  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
925  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
926  c1->SaveAs(name);
927 
928  fRootFile->WriteObject(c1,c1->GetName());
929  fRootFile->WriteObject(plot2D,plot2D->GetName());
930  fRootFile->WriteObject(distrib,distrib->GetName());
931  fRootFile->WriteObject(inhisto,inhisto->GetName());
932  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
933  //. . . Mark the bad cells in the fFlag array
934  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
935  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
936  //. . .(0 by default - good cell)
937  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
938  cout<<" o Flag bad cells that are outside the good range "<<endl;
939  for(Int_t cell = 0; cell < fNoOfCells; cell++)
940  {
941  //..cell=0 and bin=1, cell=1 and bin=2
942  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
943  if (inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)//ELI
944  {
945  fFlag[cell]=fCriterionCounter;
946  }
947  if (inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
948  {
949  fFlag[cell]=fCriterionCounter;
950  }
951  }
952  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;
953 }
954 
955 
956 
957 
962 //________________________________________________________________________
964 {
965  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
966  //.. RESULTS
967  //.. 1) Print the bad cells
968  //.. and write the results to a file
969  //.. for each added period analysis
970  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
971  TArrayD periodArray;
972  Double_t emin,emax,sig;
973  Int_t criterion;
974  TString output;
975  Int_t nb1=0;
976  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
977  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
978  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
979  TString aliceTwikiTable = Form("%s/%s/%s%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
980  ofstream file2(aliceTwikiTable, ios::out | ios::trunc);
981  if(file2)
982  {
983  file2<<"|*Criterium* | *E_min !GeV* | *E_max !GeV* | *sigma* | *Excluded Cells* |"<<endl;
984  }
985  file2.close();
986 
987  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
988  {
989  periodArray=fAnalysisVector.at(i);
990  criterion =periodArray.At(0);
991  emin =periodArray.At(2);
992  emax =periodArray.At(3);
993  sig =periodArray.At(1);
994 
995  //..Print the results on the screen and
996  //..write the results in a file
997  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
998  ofstream file(output, ios::out | ios::trunc);
999  if(!file)
1000  {
1001  cout<<"#### Major Error. Check the textfile!"<<endl;
1002  }
1003  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1004  cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1005 
1006  nb1=0;
1007  for(Int_t cellID=0;cellID<fNoOfCells;cellID++)
1008  {
1009  if(fFlag[cellID]==(i+2))
1010  {
1011  nb1++;
1012  file<<cellID<<", ";
1013  }
1014  }
1015  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1016  file<<"("<<nb1<<")"<<endl;
1017  file.close();
1018  cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1019  cout<<endl;
1020  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1021  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1022  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1023  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1024  if(file2)
1025  {
1026  file2<<"| "<<criterion<<" | "<<emin<<" | "<<emax<<" | "<<sig<<" | "<<nb1<<" |"<<endl;
1027  }
1028  file2.close();
1029  }
1030 }
1037 //________________________________________________________________________
1039 {
1040  Int_t cellID, nDeadDCalCells = 0, nDeadEMCalCells = 0, nDCalCells = 0, nEMCalCells = 0;
1041  Double_t perDeadEMCal,perDeadDCal,perBadEMCal,perBadDCal,perWarmEMCal,perWarmDCal;
1042  TString aliceTwikiTable, cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells;
1043  TH2F* cellAmp_masked= (TH2F*)fCellAmplitude->Clone("cellAmp_masked");
1044  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1045 
1046  deadPdfName = Form("%s/%s/%s%s_Dead_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1047  badPdfName = Form("%s/%s/%s%s_Bad_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1048  cellSummaryFile = Form("%s/%s/%s%s_Bad_Amplitudes_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1049  aliceTwikiTable = Form("%s/%s/%s%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1050  ratioOfBad = Form("%s/%s/%s%s_BCRatio_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1051  goodCells = Form("%s/%s/%s%s_Good_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1052 
1053  cout<<" o Final results o "<<endl;
1054 
1055  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1056  //..Write the final results of dead and bad cells in a file and on screen
1057  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1058  ofstream file(cellSummaryFile, ios::out | ios::trunc);
1059  if(file)
1060  {
1061  file<<"Dead cells : "<<endl;
1062  cout<<" o Dead cells : "<<endl;
1063  for(cellID=0; cellID<fNoOfCells; cellID++)
1064  {
1065  if(cellID==0)
1066  {
1067  file<<"In EMCal: "<<endl;
1068  //cout<<" o In EMCal : "<<endl;
1069  }
1070  if(cellID==fCellStartDCal)
1071  {
1072  file<<"\n"<<endl;
1073  file<<"In DCal: "<<endl;
1074  //cout<<endl;
1075  //cout<<" o In DCal : "<<endl;
1076  }
1077  if(fFlag[cellID]==1)
1078  {
1079  //file<<cellID<<"\n" ;
1080  file<<cellID<<", ";
1081  //cout<<cellID<<"," ;
1082  if(cellID<fCellStartDCal)nDeadEMCalCells++;
1083  else nDeadDCalCells++;
1084  }
1085  }
1086  //cout<<endl;
1087  file<<"\n"<<endl;
1088  perDeadEMCal=100*nDeadEMCalCells/(1.0*fCellStartDCal);
1089  perDeadDCal=100*nDeadDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1090  file<<"EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1091  cout<<" o EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1092 
1093  file<<"Bad cells: "<<endl;
1094  cout<<" o Bad cells: "<<endl;
1095  for(cellID=0;cellID<fNoOfCells;cellID++)
1096  {
1097  if(cellID==0)
1098  {
1099  file<<"In EMCal: "<<endl;
1100  //cout<<" o In EMCal : "<<endl;
1101  }
1102  if(cellID==fCellStartDCal)
1103  {
1104  file<<"\n"<<endl;
1105  file<<"In DCal: "<<endl;
1106  //cout<<endl;
1107  //cout<<" o In DCal : "<<endl;
1108  }
1109  if(fFlag[cellID]>1)
1110  {
1111  //file<<cellID<<"\n" ;
1112  file<<cellID<<", ";
1113  //cout<<cellID<<"," ;
1114  if(cellID<fCellStartDCal)nEMCalCells++;
1115  else nDCalCells++;
1116  }
1117  }
1118  //cout<<endl;
1119  file<<"\n"<<endl;
1120  perBadEMCal=100*nEMCalCells/(1.0*fCellStartDCal);
1121  perBadDCal =100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1122  file<<"EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1123  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1124  }
1125  file.close();
1126  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1127  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1128  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1129  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1130  if(file2)
1131  {
1132  file2<<"1=energy/hit, 2= hit/event"<<endl;
1133  file2<<"\n"<<endl;
1134  file2<<"| *Detector* | *No of cells* | *percentage* |"<<endl;
1135  file2<<"| Dead EMCal | "<<nDeadEMCalCells<<" | "<<perDeadEMCal<<"% |"<<endl;
1136  file2<<"| Bad EMCal | "<<nEMCalCells<<" | "<<perBadEMCal<<"% |"<<endl;
1137  file2<<"EMCal Warm cell row"<<endl;
1138  file2<<"| Dead DCal | "<<nDeadDCalCells<<" | "<<perDeadDCal<<"% |"<<endl;
1139  file2<<"| Bad DCal | "<<nDCalCells<<" | "<<perBadDCal<<"% |"<<endl;
1140  }
1141  file2.close();
1142  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1143  //..Save the flagged cells to .pdf files
1144  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1145  cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1146  SaveBadCellsToPDF(1,badPdfName) ;
1147  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
1148  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1149 
1150 
1151  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1152  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1153  //..And Fill the histograms with the flag information
1154  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1155  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1156  {
1157  //..Direction of amplitude (Checks energies from 0-10 GeV)
1158  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1159  {
1160  if(fFlag[cell]!=0)
1161  {
1162  //..cellID+1 = histogram bin
1163  cellAmp_masked->SetBinContent(amp,cell+1,0);
1164  cellTime_masked->SetBinContent(amp,cell+1,0);
1165  }
1166  }
1167  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
1168  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
1169  }
1170 
1171  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1172  //..Plot the 2D distribution of cells by flag
1173  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1174  PlotFlaggedCells2D(0); //..all good cells
1175  PlotFlaggedCells2D(1); //..all dead cells
1176  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
1177  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
1178 
1179 
1180  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1181  //..Plot some summary canvases
1182  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1183  TCanvas *c1 = new TCanvas("CellProp","summary of cell properties",1000,1000);
1184  c1->ToggleEventStatus();
1185  c1->Divide(2,2);
1186  c1->cd(1)->SetLogz();
1187  //lowerPadRight->SetLeftMargin(0.09);
1188  //lowerPadRight->SetRightMargin(0.06);
1189  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1190  fCellAmplitude->SetYTitle("Abs. Cell Id");
1191  fCellAmplitude->GetYaxis()->SetTitleOffset(1.6);
1192  fCellAmplitude->Draw("colz");
1193  c1->cd(2)->SetLogz();
1194  fCellTime->SetXTitle("Cell Time [ns]");
1195  fCellTime->SetYTitle("Abs. Cell Id");
1196  fCellTime->GetYaxis()->SetTitleOffset(1.6);
1197  fCellTime->Draw("colz");
1198  c1->cd(3)->SetLogz();
1199  //lowerPadRight->SetLeftMargin(0.09);
1200  //lowerPadRight->SetRightMargin(0.06);
1201  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1202  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1203  cellAmp_masked->SetYTitle("Abs. Cell Id");
1204  cellAmp_masked->GetYaxis()->SetTitleOffset(1.6);
1205  cellAmp_masked->Draw("colz");
1206  c1->cd(4)->SetLogz();
1207  cellTime_masked->SetTitle("Masked Cell Time");
1208  cellTime_masked->SetXTitle("Cell Time [ns]");
1209  cellTime_masked->SetYTitle("Abs. Cell Id");
1210  cellTime_masked->GetYaxis()->SetTitleOffset(1.6);
1211  cellTime_masked->Draw("colz");
1212  c1->Update();
1213 
1214  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
1215  c2->ToggleEventStatus();
1216  c2->Divide(1,2);
1217  c2->cd(1);
1218  fhCellFlag->SetTitle("cell flag");
1219  fhCellFlag->SetXTitle("Abs. Cell Id");
1220  fhCellFlag->SetYTitle("flag");
1221  fhCellFlag->DrawCopy("hist");
1222  c2->cd(2);
1223  fhCellWarm->SetTitle("Warm cells");
1224  fhCellWarm->SetXTitle("Abs. Cell Id");
1225  fhCellWarm->SetYTitle("warm=1");
1226  fhCellWarm->DrawCopy("hist");
1227  c2->Update();
1228  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1229  //..Add different histograms/canvases to the output root file
1230  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1231  TString name =Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1232  c1->SaveAs(name);
1233  fRootFile->WriteObject(c1,c1->GetName());
1234  fRootFile->WriteObject(c2,c2->GetName());
1235  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
1236  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
1237  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
1238  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
1239  //..Save all amplitudes to the root file
1240  SaveHistoToFile();
1241 
1242  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1243  //..Save also the identified warm channels into a text file.
1244  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1245  file.open(cellSummaryFile, ios::out | ios::app);
1246  if(file)
1247  {
1248  file<<"Warm cells : "<<endl;
1249  cout<<" o Warm cells : "<<endl;
1250  nEMCalCells =0;
1251  nDCalCells =0;
1252  for(cellID=0; cellID<fNoOfCells; cellID++)
1253  {
1254  if(cellID==0)
1255  {
1256  file<<"In EMCal: "<<endl;
1257  }
1258  if(cellID==fCellStartDCal)
1259  {
1260  file<<"\n"<<endl;
1261  file<<"In DCal: "<<endl;
1262  }
1263  if(fWarmCell[cellID]==1)
1264  {
1265  file<<cellID<<", ";
1266  if(cellID<fCellStartDCal)nEMCalCells++;
1267  else nDCalCells++;
1268  }
1269  }
1270  file<<"\n"<<endl;
1271  perWarmEMCal= 100*nEMCalCells/(1.0*fCellStartDCal);
1272  perWarmDCal = 100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1273  file<<"EMCal ("<<nEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1274  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1275  }
1276  file.close();
1277  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1278  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1279  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1280  ofstream file3(aliceTwikiTable, ios::out | ios::app);
1281  if(file3)
1282  {
1283  file3<<"| Warm DCal | "<<nDCalCells<<" | "<<perWarmDCal<<"% |"<<endl;
1284  file3<<"\n"<<endl;
1285  file3<<"Insert above:"<<endl;
1286  file3<<"| Warm EMCal | "<<nEMCalCells<<" | "<<perWarmEMCal<<"% |"<<endl;
1287  }
1288  file3.close();
1289  //cout<<" o Results can be found in : "<<endl;
1290  //cout<<" o "<<cellSummaryFile<<endl;
1291  //cout<<" o "<<badPdfName<<endl;
1292  //cout<<" o "<<badPdfName<<endl;
1293 }
1294 
1295 
1310 //________________________________________________________________________
1312 {
1313  gROOT->SetStyle("Plain");
1314  gStyle->SetOptStat(0);
1315  gStyle->SetFillColor(kWhite);
1316  gStyle->SetTitleFillColor(kWhite);
1317  gStyle->SetPalette(1);
1318 
1319  char title[100];
1320  char name[100];
1321 
1322  TH1D *hRefDistr = BuildMeanFromGood();
1323  fRootFile->WriteObject(hRefDistr,hRefDistr->GetName());
1324  Int_t firstCanvas=0;
1325  Bool_t candidate;
1326  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
1327  text->SetTextSize(0.07);
1328  text->SetTextColor(kOrange-3);
1329  text->SetNDC();
1330 
1331  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
1332  text2->SetTextSize(0.08);
1333  text2->SetTextColor(8);
1334 
1335  TLatex* textA = new TLatex(0.65,0.62,"*test*");
1336  textA->SetTextSize(0.04);
1337  textA->SetTextColor(1);
1338  textA->SetNDC();
1339 
1340  //..collect cells in an internal vector.
1341  //..when the vector is of size=9 or at the end of being filled
1342  //..plot the channels into a canvas
1343  std::vector<Int_t> channelVector;
1344  channelVector.clear();
1345  cout<<"Start printing into .pdf for version: "<<version<<endl;
1346  for(Int_t cell=0;cell<fNoOfCells;cell++)
1347  {
1348  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
1349  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
1350  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
1351  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
1352 
1353  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
1354  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1355  if(channelVector.size()==9 || cell == fNoOfCells-1)
1356  {
1357  cout<<"."<<flush;
1358 
1359  TString internal_pdfName=pdfName;
1360  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1361  if(channelVector.size() > 6) c1->Divide(3,3);
1362  else if (channelVector.size() > 3) c1->Divide(3,2);
1363  else c1->Divide(3,1);
1364 
1365  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1366  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
1367  {
1368  sprintf(name, "Cell %d",channelVector.at(i)) ;
1369  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
1370  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
1371 
1372  c1->cd(i%9 + 1);
1373  c1->cd(i%9 + 1)->SetLogy();
1374  if(i==0)
1375  {
1376  leg->AddEntry(hRefDistr,"mean of good","l");
1377  leg->AddEntry(hCell,"current","l");
1378  }
1379  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
1380  candidate = CheckDistribution(hCell,hRefDistr);
1381  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
1382  if(version>1)//..These are ratio plots of energy distr. of cell and mean of all good cells
1383  {
1384  hCell->Divide(hRefDistr);
1385  }
1386  //.. save histograms to file
1387  if(version==1) fOutputListBad->Add(hCell);
1388  if(version==10)fOutputListBadRatio->Add(hCell);
1389  if(version==2) fOutputListGoodRatio->Add(hCell);
1390 
1391  hCell->SetLineColor(kBlue+1);
1392  hCell->GetXaxis()->SetTitle("E (GeV)");
1393  hCell->GetYaxis()->SetTitle("N Entries");
1394  hCell->GetXaxis()->SetRangeUser(0.,10.);
1395  hCell->SetLineWidth(1) ;
1396  hCell->SetTitle(title);
1397  hRefDistr->SetLineColor(kGray+2);
1398  hRefDistr->SetLineWidth(1);
1399 
1400  hCell->Draw("hist");
1401 
1402  if(version==1)hRefDistr->Draw("same") ;
1403 
1404  //..Mark the histogram that could be miscalibrated and labelled as warm
1405  if(candidate==1 && (version==1 || version==10))
1406  {
1407  gPad->SetFrameFillColor(kYellow-10);
1408  text->Draw();
1409  }
1410  if(version==1)
1411  {
1412  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
1413  textA->Draw();
1414  }
1415  if(version==2 && candidate==0)
1416  {
1417  gPad->SetFrameFillColor(kYellow-10);
1418  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
1419  leg->Draw();
1420  }
1421  if(version<2)leg->Draw();
1422  }
1423 
1424  if(channelVector.size()<9 || cell == fNoOfCells-1)
1425  {
1426  internal_pdfName +=")";
1427  c1->Print(internal_pdfName.Data());
1428  }
1429  else if(firstCanvas==0)
1430  {
1431  internal_pdfName +="(";
1432  c1->Print(internal_pdfName.Data());
1433  firstCanvas=1;
1434  }
1435  else
1436  {
1437  c1->Print(internal_pdfName.Data());
1438  }
1439  delete c1;
1440  delete leg;
1441  channelVector.clear();
1442  }
1443  }
1444  delete hRefDistr;
1445  //..Add the subdirectories to the file
1446  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
1447  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
1448  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
1449 
1450  cout<<endl;
1451 }
1455 //_________________________________________________________________________
1457 {
1458  TH1D* hGoodAmp;
1459  TH1D* hgoodMean;
1460  Int_t NrGood=0;
1461  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1462  {
1463  if(fFlag[cell]!=0)continue;
1464  NrGood++;
1465  if(NrGood==1)hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean",cell+1,cell+1);
1466  else
1467  {
1468  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1469  hgoodMean->Add(hGoodAmp);
1470  }
1471  }
1472  hgoodMean->Scale(1.0/NrGood);
1473 
1474  return hgoodMean;
1475 }
1486 //_________________________________________________________________________
1488 {
1489  TString Name = Form("%s_ratio",histogram->GetName());
1490  TH1* ratio=(TH1*)histogram->Clone(Name);
1491  ratio->Divide(reference);
1492 
1493  Double_t percentageOfLast=0.6;
1494  Double_t higherThanMean=5;
1495  Double_t highestRatio=1000;
1496  Double_t maxMean=10;
1497  Double_t minMean=0.1;
1498  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
1499 
1500  //..By default each cell is a candidate for recalibration
1501  Bool_t candidate=1;
1502  //..Find bin where reference has value 1, and the corresponding x-value
1503  Int_t binHeihgtOne = reference->FindLastBinAbove(1);
1504  Double_t binCentreHeightOne = reference->GetBinCenter(binHeihgtOne);
1505  Double_t thirdBinCentre = reference->GetBinCenter(3);
1506 
1507  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1508  //..Check the histogram
1509  //..Different checks to see whether the
1510  //..cell is really bad. Set candidate to 0.
1511 
1512  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1513  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
1514  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
1515  {
1516  candidate=0;
1517  }
1518 
1519  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1520  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
1521  //ELI check that crieteria carfully - seems to work but not shure about it
1522  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
1523  if(ratio->GetMaximum()>highestRatio)//
1524  {
1525  candidate=0;
1526  }
1527 
1528  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1529  //..check whether the ratio is much larger than 1
1530  //..calculate the mean in the relevant energy range
1531  Double_t mean=0;
1532  for(Int_t i=2;i<binHeihgtOne;i++)
1533  {
1534  mean+=ratio->GetBinContent(i);
1535  }
1536  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1537  if(mean>maxMean || mean<minMean)
1538  {
1539  candidate=0;
1540  }
1541 
1542  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1543  //..check whether there are large spikes in the histogram
1544  //..compare bin values to mean of the ratio. If there is a bin value with
1545  //..content "higherThanMean" times lareger than mean it's losing it candidate status
1546  mean=0;
1547  //..Find the maximum in the mean range (0-binHeihgtOne)
1548  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
1549  Double_t localMaxBin=ratio->GetMaximumBin();
1550 
1551  for(Int_t i=2;i<binHeihgtOne;i++)
1552  {
1553  //..Exclude 0 bins and exclude bins near the maximum
1554  if(ratio->GetBinContent(i)<=0) continue;
1555  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
1556  mean+=ratio->GetBinContent(i);
1557  }
1558  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1559  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
1560  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1561  if(ratio->GetMaximum()>mean*higherThanMean)
1562  {
1563  candidate=0;
1564  }
1565 
1566  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1567  //..Check for a cliff at 4GeV, happens in some cases
1568  Double_t beforeCliff=0;
1569  Double_t afterCliff=0;
1570  Int_t binsBefore=0;
1571  Int_t binsAfter=0;
1572  Int_t cliffBin = ratio->FindBin(4);
1573  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
1574  {
1575  //..Exclude 0 bins
1576  if(ratio->GetBinContent(i)<=0)continue;
1577  if(i<=cliffBin) binsBefore++;
1578  if(i>cliffBin) binsAfter++;
1579  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
1580  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
1581  beforeCliff*=1.0/binsBefore;
1582  afterCliff*=1.0/binsAfter;
1583  }
1584  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
1585  {
1586  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
1587  }
1588 
1589  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1590  //..Employ peak finder
1591 /* if(candidate==1)
1592  {
1593  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
1594  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
1595  //GetNPeaks();
1596  //cout<<"found N peaks: "<<nfound<<endl;
1597  }
1598 */
1599  return candidate;
1600 }
1601 
1610 //_________________________________________________________________________
1612 {
1613  //..TRD support structure
1614  //..(determined by eye, could be improved, but is already very acurate):
1615  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
1616  //..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
1617  Bool_t coveredByTRDSupportStruc=0;
1618 
1619  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
1620  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
1621  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
1622  )
1623  {
1624  coveredByTRDSupportStruc=1;
1625  }
1626  return coveredByTRDSupportStruc;
1627 }
1636 //_________________________________________________________________________
1638 {
1639  //..build two dimensional histogram with values row vs. column
1640  TString histoName;
1641  histoName = Form("2DChannelMap_Flag%d",flagBegin);
1642  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100");
1643 
1644  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
1645  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1646  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1647 
1648  Int_t cellColumn=0,cellRow=0;
1649  Int_t cellColumnAbs=0,cellRowAbs=0;
1650  Int_t trash;
1651 
1652  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1653  {
1654  //..Do that only for cell ids also accepted by the code
1655  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1656 
1657  //..Get Row and Collumn for cell ID c
1658  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1659  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1660  {
1661  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1662  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1663  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1664  }
1665  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1666  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1667  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
1668 
1669 
1670  }
1671  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
1672  c1->ToggleEventStatus();
1673  c1->cd();
1674  //lowerPadRight->SetLeftMargin(0.09);
1675  //lowerPadRight->SetRightMargin(0.06);
1676  plot2D->Draw("colz");
1677 
1678  TLatex* text = 0x0;
1679  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
1680  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
1681  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
1682  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
1683  text->SetTextSize(0.06);
1684  text->SetNDC();
1685  text->SetTextColor(1);
1686  text->Draw();
1687 
1688  c1->Update();
1689  TString name =Form("%s/%s/%s_%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data() , histoName.Data());
1690  c1->SaveAs(name);
1691 
1692  fRootFile->WriteObject(plot2D,plot2D->GetName());
1693 
1694 }
1698 //_________________________________________________________________________
1700 {
1701  char name[100];
1702  for(Int_t cell=0;cell<fNoOfCells;cell++)
1703  {
1704  sprintf(name, "Cell %d",cell) ;
1705  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
1706  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
1707  }
1708  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1709 }
TH1F * BuildTimeMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
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.
TString fPass
Pass of the analyzed data.
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
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 './'.
void FlagAsBad(Int_t crit, TH1F *inhisto, Double_t nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1.)
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.
const char * pdfName
Definition: DrawAnaELoss.C:30
void PlotFlaggedCells2D(Int_t flagBegin, Int_t flagEnd=-1)
TH1F * BuildHitAndEnergyMean(Int_t crit, Double_t emin=0.1, Double_t emax=2., Double_t nsigma=4.)
TString fTrigger
Selected trigger for the analysis.
Bool_t CheckDistribution(TH1 *ratio, TH1 *reference)
Int_t fNMaxRowsAbs
Maximum No of rows in Calorimeter.
Int_t fNMaxColsAbs
Maximum No of colums in Calorimeter.
void PeriodAnalysis(Int_t criterum=7, Double_t nsigma=4.0, Double_t emin=0.1, Double_t emax=2.0)
Bool_t fTestRoutine
This is a flag, if set true will produce some extra quality check histograms.
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.
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 ...
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
ClassImp(AliAnalysisTaskCascadeTester) AliAnalysisTaskCascadeTester
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