AliPhysics  97344c9 (97344c9)
 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  else
451  {
452  //correct time in case there is nothing in AODB!
453  }
454  hNEvents =(TH1F *)outputList->FindObject("hNEvents");
455  if(!hNEvents)
456  {
457  Printf("hNEvents not found");
458  outputList->ls();
459  continue;
460  }
461  nEntr = (Int_t)hNEvents->GetEntries();
462 
463  //..does that mean do not merge small files?
464  if (nEntr<100)
465  {
466  cout <<" o File to small to be merged. Only N entries " << nEntr << endl;
467  continue ;
468  }
469  cout <<" o File with N entries " << nEntr<<" will be merged"<< endl;
470  nEntrTot+=nEntr;
471  hCellAmplitude->Add(hAmpId);
472  hCellTime->Add(hTimeId);
473  hNEventsProcessedPerRun->Add(hNEvents);
474 
475  //..Create copies of the original root files just with the bad channel QA
476  //..So that a run by run bad channel analysis can be performed more easily
477  singleRunFileName= Form("%s/%s/%s/%d_Filtered.root",fWorkdir.Data(),fAnalysisInput.Data(),fPass.Data(),runId[i]);
478  TFile *singleRunFile = TFile::Open(singleRunFileName,"recreate");
479  hNEventsProcessedPerRun->Write();
480  hCellAmplitude->Write();
481  hCellTime->Write();
482  singleRunFile->Close();
483 
484  outputList->Delete();
485  dir->Delete();
486  f->Close();
487  delete f;
488  }
489 
490  //.. Save the merged histograms
491  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
492  cout<<"o o o "<<nEntrTot<<" events were merged"<<endl;
493  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
494  hNEventsProcessedPerRun->Write();
495  hCellAmplitude->Write();
496  hCellTime->Write();
497  BCF->Close();
498  cout<<"o o o End conversion process o o o"<<endl;
499  return fMergedFileName;
500 }
501 
502 
508 //_________________________________________________________________________
510 {
511  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
512  //.. BAD CELLS
513  //.. Flag bad cells with fFlag= 2,3,4,5.. etc
514  //.. this excludes cells from subsequent analysis
515  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
516  TArrayD periodArray;
517  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
518  {
519  periodArray=fAnalysisVector.at(i);
520  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
522  cout<<""<<endl;
523  cout<<""<<endl;
524  }
525  cout<<"o o o End of bad channel analysis o o o"<<endl;
526 }
527 
536 //________________________________________________________________________
538 {
539  TArrayD periodArray;
540  periodArray.Set(4);
541  periodArray.AddAt((Double_t)criteria,0);
542  periodArray.AddAt(nsigma,1);
543  periodArray.AddAt(emin,2);
544  periodArray.AddAt(emax,3);
545  fAnalysisVector.push_back(periodArray);
546 }
562 //____________________________________________________________________
564 {
565  //ELI criterion should be between 1-4
566  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;
567  cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
568  if(criterion < 3)cout<<"o o o Done in the energy range E "<<emin<<" - "<<emax<<endl;
569  if(criterion == 3)cout<<"o o o Done in the time range t "<<emin<<" - "<<emax<<endl;
570 
571  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
572  //.. ANALYSIS OF CELLS WITH ENTRIES
573  //.. Build average distributions and fit them
574  //.. Three tests for bad cells:
575  //.. 1) Average energy per hit
576  //.. 2) Average hit per event
577  //.. 3) Average max of cell time distribution
578  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
579  TH1F* histogram;
580  cout<<"o o o Analyze average cell distributions o o o"<<endl;
581  //..For case 1 or 2
582  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax);
583  if(criterion == 3) histogram = BuildTimeMean(criterion, emin, emax); //.. in case of crit 3 emin is tmin and emax is tmax
584  //..For case 3, 4 or 5
585  //else if (criterion < 6)
586 
587  //..light systems p+p p+Pb
588 /* Double_t range=0.001;
589  if(emin>0.49)range=0.0005;
590  if(emin>0.99)range=0.0001;
591  if(emin>1.99)range=0.00005;
592 */
593  //..for Pb+Pb
594  Double_t range=0.1;
595  if(emin>0.49)range=0.05;
596  if(emin>0.99)range=0.01;
597  if(emin>1.99)range=0.0002;
598 
599 
600  if(criterion==1) FlagAsBad(criterion, histogram, nsigma, 400,-1);
601  if(criterion==2) FlagAsBad(criterion, histogram, nsigma, 600,range);
602  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, 600,30);
603 
604  /*
605  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval3);
606  if(criterion==4) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval1);
607  if(criterion==5) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval2);
608  */
609 }
610 
619 //_________________________________________________________________________
621 {
622  cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
623  TH1F *histogram;
624  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);
625  if(crit==2)histogram = new TH1F(Form("hCellNtoEvt_E%.2f-%.2f",emin,emax),Form("Number of hits per event, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
626  histogram->SetXTitle("Abs. Cell Id");
627  if(crit==1)histogram->SetYTitle("Energy per hit");
628  if(crit==2)histogram->SetYTitle("Number of hits per event");
629  histogram->GetXaxis()->SetNdivisions(510);
630  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
631 
632  //..here the average hit per event and the average energy per hit is caluclated for each cell.
633  for (Int_t cell = 0; cell < fNoOfCells; cell++)
634  {
635  Double_t Esum = 0;
636  Double_t Nsum = 0;
637 
638  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
639  {
640  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
641  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
642  //..investigate only cells that were not flagged as dead ore bad
643  if (E < emin || E > emax || fFlag[cell]!=0) continue;
644  Esum += E*N;
645  Nsum += N;
646  }
647  //..Set the values only for cells that are not yet marked as bad
648  if(fFlag[cell]==0)
649  {
650  if(totalevents> 0. && crit==2)histogram->SetBinContent(cell+1, Nsum/totalevents); //..number of hits per event
651  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/Nsum); //..average energy per hit
652  }
653  }
654  return histogram;
655 }
660 //_________________________________________________________________________
662 {
663  cout<<" o Calculate maximum of cell time distribution "<<endl;
664  TString name;
665  TH1F *histogram;
666  histogram = new TH1F(Form("timeMax_t%.2f-%.2f",tmin,tmax),Form("Maximum of time distr., %.2f < t < %.2f GeV",tmin,tmax), fNoOfCells,0,fNoOfCells);
667  histogram->SetXTitle("Abs. Cell Id");
668  histogram->SetYTitle("time max");
669  histogram->GetXaxis()->SetNdivisions(510);
670  //Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
671 
672  //..Time information
673  for (Int_t cell = 0; cell < fNoOfCells; cell++)
674  {
675  Double_t timeMax = -100;
676  //..search for the maximum bin between tmin and tmax
677  name=Form("Cell %d",cell);
678  //sprintf(name, "Cell %d",channelVector.at(i)) ;
679  TH1 *hCell = fCellTime->ProjectionX(name,cell+1,cell+1);
680  hCell->GetXaxis()->SetRangeUser(tmin,tmax);
681 
682  timeMax = hCell->GetXaxis()->GetBinCenter(hCell->GetMaximumBin());
683  if(cell>55 &&cell<65)cout<<"Integral: "<<hCell->Integral()<<endl;
684  if(hCell->Integral()<=0)timeMax=-40;
685  //..Set the values only for cells that are not yet marked as bad
686  if(fFlag[cell]==0)
687  {
688  histogram->SetBinContent(cell+1, timeMax); //..number of hits per event
689  }
690  else
691  {
692  histogram->SetBinContent(cell+1, -50); //..number of hits per event
693  }
694 
695  }
696  return histogram;
697 }
698 
703 //_________________________________________________________________________
705 {
706  Int_t sumOfExcl=0;
707  //..Direction of cell ID
708  for (Int_t cell = 0; cell < fNoOfCells; cell++)
709  {
710  Double_t nSum = 0;
711  //..Direction of amplitude (Checks energies from 0-10 GeV)
712  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
713  {
714  //..cellID+1 = histogram bin
715  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
716  nSum += N;
717  }
718  //..If the amplitude in one cell is basically 0
719  //..mark the cell as excluded
720  if(nSum < 0.5 && fFlag[cell]==0)
721  {
722  //..Cell flagged as dead.
723  //..Flag only if it hasn't been flagged before
724  fFlag[cell] =1;
725  sumOfExcl++;
726  }
727  }
728  cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
729  cout<<" ("<<sumOfExcl<<")"<<endl;
730 }
731 
746 //_________________________________________________________________________
747 void AliAnaCaloChannelAnalysis::FlagAsBad(Int_t crit, TH1F* inhisto, Double_t nsigma, Int_t dnbins, Double_t dmaxVal)
748 {
749  gStyle->SetOptStat(0); // MG modif
750  gStyle->SetOptFit(0); // MG modif
751 
752  if(crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
753  if(crit==2)cout<<" o Fit average hit per event distribution"<<endl;
754  if(crit==3)cout<<" o Fit average hit maximum distribution"<<endl;
755 
756  Int_t cellColumn=0,cellRow=0;
757  Int_t cellColumnAbs=0,cellRowAbs=0;
758  Int_t trash;
759 
760  TString histoName=inhisto->GetName();
761  Double_t goodmax= 0. ;
762  Double_t goodmin= 0. ;
763  Double_t dminVal=inhisto->GetMinimum();
764  if (dmaxVal < 0.)
765  {
766  dmaxVal = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
767  if(crit==2 && dmaxVal > 1) dmaxVal =1. ;
768  }
769 
770  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
771  //. . .build the distribution of average values
772  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
773  if(crit<3)
774  {
775  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
776  //..Be aware that the bin width should never be larger than 1/totalevents because this
777  //..is the minimum differce between the cells. One hit more/event otherwise you will see ugly empty bins in the histogram
778  if(((dmaxVal-inhisto->GetMinimum())/(dnbins*1.0))<1.0/totalevents)
779  {
780  dnbins=(dmaxVal-inhisto->GetMinimum())*totalevents;
781  cout<<"Problem - Reset dnbins to new value:"<<dnbins<<endl;
782  }
783  }
784  if(crit==3)
785  {
786  //..obtain the binwidth for
787  Int_t binWidth = fCellTime->GetXaxis()->GetBinWidth(1);
788  dnbins = fCellTime->GetXaxis()->GetNbins();
789  dminVal= fCellTime->GetXaxis()->GetBinCenter(1)-(binWidth*1.0)/2.0;
790  dmaxVal= fCellTime->GetXaxis()->GetBinCenter(dnbins)+(binWidth*1.0)/2.0;
791 
792  cout<<"set the new histo with following settings- bins: "<<dnbins<<", min val = "<<dminVal<<", max val:"<<dmaxVal<<endl;
793  }
794  //..build histos
795  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
796  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
797  distrib->SetYTitle("Entries");
798  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
799  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
800 
801 
802  //..build two dimensional histogram with values row vs. column
803  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);
804  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
805  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
806 
807  for (Int_t cell = 0; cell < fNoOfCells; cell++)
808  {
809  //..Do that only for cell ids also accepted by the code
810  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
811 
812  //..Fill histograms only for cells that are not yet flagged as bad
813  if(fFlag[cell]==0)
814  {
815  //..fill the distribution of avarge cell values
816  distrib->Fill(inhisto->GetBinContent(cell+1));
817  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
818  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
819  //..Get Row and Collumn for cell ID
820  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
821  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
822  {
823  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
824  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
825  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
826  }
827  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
828  //..check TRD support structure
829  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
830  {
831  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
832  }
833  else
834  {
835  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
836  }
837  }
838  }
839 
840  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
841  //. . .draw histogram + distribution
842  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
843 
844  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
845  c1->ToggleEventStatus();
846  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
847  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
848  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
849  upperPad->Draw();
850  lowerPadLeft->Draw();
851  lowerPadRight->Draw();
852 
853  upperPad->cd();
854  upperPad->SetLeftMargin(0.05);
855  upperPad->SetRightMargin(0.05);
856  if(crit<3)upperPad->SetLogy();
857  inhisto->SetTitleOffset(0.6,"Y");
858  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
859  inhisto->GetYaxis()->SetTitleOffset(0.7);
860  inhisto->SetLineColor(kBlue+1);
861  inhisto->Draw();
862 
863  lowerPadRight->cd();
864  lowerPadRight->SetLeftMargin(0.09);
865  lowerPadRight->SetRightMargin(0.12);
866  plot2D->GetYaxis()->SetTitleOffset(1.3);
867  plot2D->Draw("colz");
868 
869  lowerPadLeft->cd();
870  lowerPadLeft->SetLeftMargin(0.09);
871  lowerPadLeft->SetRightMargin(0.07);
872  lowerPadLeft->SetLogy();
873  distrib->SetLineColor(kBlue+1);
874  distrib->GetYaxis()->SetTitleOffset(1.3);
875  distrib->Draw();
876  distrib_wTRDStruc->SetLineColor(kGreen+1);
877  distrib_wTRDStruc->DrawCopy("same");
878  distrib_woTRDStruc->SetLineColor(kMagenta+1);
879  distrib_woTRDStruc->DrawCopy("same");
880 
881  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
882  //. . .fit histogram
883  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
884  Int_t higherbin=0,i;
885  for(i = 2; i <= dnbins; i++)
886  {
887  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i)) higherbin = i ;
888  }
889  //..good range is around the max value as long as the
890  //..bin content is larger than 2 entries
891  for(i = higherbin ; i<=dnbins ; i++)
892  {
893  if(distrib->GetBinContent(i)<2) break ;
894  goodmax = distrib->GetBinCenter(i);
895  }
896  for(i = higherbin ; i>1 ; i--)
897  {
898  if(distrib->GetBinContent(i)<2) break ;
899  goodmin = distrib->GetBinLowEdge(i);
900  }
901 
902  //..Define min/max range of the fit function
903  Double_t minFitRange=0;
904  Int_t maxFitRange=10;
905  if(crit==3)minFitRange=-20;
906  if(crit==3)maxFitRange=20;
907 
908  TF1 *fit2 = new TF1("fit2", "gaus",minFitRange,maxFitRange);
909  //..start the fit with a mean of the highest value
910  fit2->SetParameter(1,higherbin);
911 
912  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
913  Double_t sig, mean;// chi2ndf;
914  // Marie midif to take into account very non gaussian distrig
915  mean = fit2->GetParameter(1);
916  sig = fit2->GetParameter(2);
917  //chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
918 
919  //..for case 1 and 2 lower than 0 is an unphysical value
920  if(crit<3 && mean <0.) mean=0.;
921 
922  goodmin = mean - nsigma*sig ;
923  goodmax = mean + nsigma*sig ;
924 
925  //..for case 1 and 2 lower than 0 is an unphysical value
926  if(crit<3 && goodmin <0.) goodmin=0.;
927 
928  cout<<" o Result of fit: "<<endl;
929  cout<<" o "<<endl;
930  cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
931  cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
932 
933  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
934  //. . .Add info to histogram
935  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
936  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
937  lline->SetLineColor(kGreen+2);
938  lline->SetLineStyle(7);
939  lline->Draw();
940 
941  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
942  rline->SetLineColor(kGreen+2);
943  rline->SetLineStyle(7);
944  rline->Draw();
945 
946  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
947  leg->AddEntry(lline,"Good region boundary","l");
948  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
949  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
950  leg->SetBorderSize(0);
951  leg->Draw("same");
952 
953  fit2->SetLineColor(kOrange-3);
954  fit2->SetLineStyle(1);//7
955  fit2->Draw("same");
956 
957  TLatex* text = 0x0;
958  if(crit==1) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
959  if(crit==2) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2fx10^-5",goodmin*100000,goodmax*100000));
960  if(crit==3) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
961  text->SetTextSize(0.06);
962  text->SetNDC();
963  text->SetTextColor(1);
964  text->Draw();
965 
966 
967  upperPad->cd();
968  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
969  uline->SetLineColor(kGreen+2);
970  uline->SetLineStyle(7);
971  uline->Draw();
972  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
973  lowline->SetLineColor(kGreen+2);
974  lowline->SetLineStyle(7);
975  lowline->Draw();
976  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
977  //. . .Save histogram
978  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
979  c1->Update();
980  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
981  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
982  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
983  if(crit==3)name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
984  c1->SaveAs(name);
985 
986  fRootFile->WriteObject(c1,c1->GetName());
987  fRootFile->WriteObject(plot2D,plot2D->GetName());
988  fRootFile->WriteObject(distrib,distrib->GetName());
989  fRootFile->WriteObject(inhisto,inhisto->GetName());
990  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
991  //. . . Mark the bad cells in the fFlag array
992  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
993  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
994  //. . .(0 by default - good cell)
995  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
996  cout<<" o Flag bad cells that are outside the good range "<<endl;
997  for(Int_t cell = 0; cell < fNoOfCells; cell++)
998  {
999  //..cell=0 and bin=1, cell=1 and bin=2
1000  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1001  if (inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)//ELI
1002  {
1003  fFlag[cell]=fCriterionCounter;
1004  }
1005  if (inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
1006  {
1007  fFlag[cell]=fCriterionCounter;
1008  }
1009  }
1010  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;
1011 }
1012 
1013 
1014 
1015 
1020 //________________________________________________________________________
1022 {
1023  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1024  //.. RESULTS
1025  //.. 1) Print the bad cells
1026  //.. and write the results to a file
1027  //.. for each added period analysis
1028  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1029  TArrayD periodArray;
1030  Double_t emin,emax,sig;
1031  Int_t criterion;
1032  TString output;
1033  Int_t nb1=0;
1034  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1035  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1036  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1037  TString aliceTwikiTable = Form("%s/%s/%s%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1038  ofstream file2(aliceTwikiTable, ios::out | ios::trunc);
1039  if(file2)
1040  {
1041  file2<<"|*Criterium* | *E_min !GeV* | *E_max !GeV* | *sigma* | *Excluded Cells* |"<<endl;
1042  }
1043  file2.close();
1044 
1045  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
1046  {
1047  periodArray=fAnalysisVector.at(i);
1048  criterion =periodArray.At(0);
1049  emin =periodArray.At(2);
1050  emax =periodArray.At(3);
1051  sig =periodArray.At(1);
1052 
1053  //..Print the results on the screen and
1054  //..write the results in a file
1055  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
1056  ofstream file(output, ios::out | ios::trunc);
1057  if(!file)
1058  {
1059  cout<<"#### Major Error. Check the textfile!"<<endl;
1060  }
1061  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1062  cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1063 
1064  nb1=0;
1065  for(Int_t cellID=0;cellID<fNoOfCells;cellID++)
1066  {
1067  if(fFlag[cellID]==(i+2))
1068  {
1069  nb1++;
1070  file<<cellID<<", ";
1071  }
1072  }
1073  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1074  file<<"("<<nb1<<")"<<endl;
1075  file.close();
1076  cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1077  cout<<endl;
1078  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1079  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1080  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1081  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1082  if(file2)
1083  {
1084  file2<<"| "<<criterion<<" | "<<emin<<" | "<<emax<<" | "<<sig<<" | "<<nb1<<" |"<<endl;
1085  }
1086  file2.close();
1087  }
1088 }
1095 //________________________________________________________________________
1097 {
1098  Int_t cellID, nDeadDCalCells = 0, nDeadEMCalCells = 0, nDCalCells = 0, nEMCalCells = 0;
1099  Double_t perDeadEMCal,perDeadDCal,perBadEMCal,perBadDCal,perWarmEMCal,perWarmDCal;
1100  TString aliceTwikiTable, cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells,goodCellsRatio;
1101  TH2F* cellAmp_masked= (TH2F*)fCellAmplitude->Clone("cellAmp_masked");
1102  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1103 
1104  deadPdfName = Form("%s/%s/%s%s_Dead_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1105  badPdfName = Form("%s/%s/%s%s_Bad_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1106  ratioOfBad = Form("%s/%s/%s%s_Bad_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1107  goodCells = Form("%s/%s/%s%s_Good_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1108  goodCellsRatio = Form("%s/%s/%s%s_Good_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1109  cellSummaryFile = Form("%s/%s/%s%s_Bad_Ampl_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1110  aliceTwikiTable = Form("%s/%s/%s%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1111 
1112  cout<<" o Final results o "<<endl;
1113 
1114  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1115  //..Write the final results of dead and bad cells in a file and on screen
1116  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1117  ofstream file(cellSummaryFile, ios::out | ios::trunc);
1118  if(file)
1119  {
1120  file<<"Dead cells : "<<endl;
1121  cout<<" o Dead cells : "<<endl;
1122  for(cellID=0; cellID<fNoOfCells; cellID++)
1123  {
1124  if(cellID==0)
1125  {
1126  file<<"In EMCal: "<<endl;
1127  //cout<<" o In EMCal : "<<endl;
1128  }
1129  if(cellID==fCellStartDCal)
1130  {
1131  file<<"\n"<<endl;
1132  file<<"In DCal: "<<endl;
1133  //cout<<endl;
1134  //cout<<" o In DCal : "<<endl;
1135  }
1136  if(fFlag[cellID]==1)
1137  {
1138  //file<<cellID<<"\n" ;
1139  file<<cellID<<", ";
1140  //cout<<cellID<<"," ;
1141  if(cellID<fCellStartDCal)nDeadEMCalCells++;
1142  else nDeadDCalCells++;
1143  }
1144  }
1145  //cout<<endl;
1146  file<<"\n"<<endl;
1147  perDeadEMCal=100*nDeadEMCalCells/(1.0*fCellStartDCal);
1148  perDeadDCal=100*nDeadDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1149  file<<"EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1150  cout<<" o EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
1151 
1152  file<<"Bad cells: "<<endl;
1153  cout<<" o Bad cells: "<<endl;
1154  for(cellID=0;cellID<fNoOfCells;cellID++)
1155  {
1156  if(cellID==0)
1157  {
1158  file<<"In EMCal: "<<endl;
1159  //cout<<" o In EMCal : "<<endl;
1160  }
1161  if(cellID==fCellStartDCal)
1162  {
1163  file<<"\n"<<endl;
1164  file<<"In DCal: "<<endl;
1165  //cout<<endl;
1166  //cout<<" o In DCal : "<<endl;
1167  }
1168  if(fFlag[cellID]>1)
1169  {
1170  //file<<cellID<<"\n" ;
1171  file<<cellID<<", ";
1172  //cout<<cellID<<"," ;
1173  if(cellID<fCellStartDCal)nEMCalCells++;
1174  else nDCalCells++;
1175  }
1176  }
1177  //cout<<endl;
1178  file<<"\n"<<endl;
1179  perBadEMCal=100*nEMCalCells/(1.0*fCellStartDCal);
1180  perBadDCal =100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1181  file<<"EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1182  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
1183  }
1184  file.close();
1185  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1186  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1187  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1188  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1189  if(file2)
1190  {
1191  file2<<"1=energy/hit, 2= hit/event"<<endl;
1192  file2<<"\n"<<endl;
1193  file2<<"| *Detector* | *No of cells* | *percentage* |"<<endl;
1194  file2<<"| Dead EMCal | "<<nDeadEMCalCells<<" | "<<perDeadEMCal<<"% |"<<endl;
1195  file2<<"| Bad EMCal | "<<nEMCalCells<<" | "<<perBadEMCal<<"% |"<<endl;
1196  file2<<"EMCal Warm cell row"<<endl;
1197  file2<<"| Dead DCal | "<<nDeadDCalCells<<" | "<<perDeadDCal<<"% |"<<endl;
1198  file2<<"| Bad DCal | "<<nDCalCells<<" | "<<perBadDCal<<"% |"<<endl;
1199  }
1200  file2.close();
1201  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1202  //..Save the flagged cells to .pdf files
1203  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1204  cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1205  SaveBadCellsToPDF(1,badPdfName) ;
1206  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
1207  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1208  if(fTestRoutine==1)SaveBadCellsToPDF(20,goodCellsRatio) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1209 
1210 
1211  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1212  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1213  //..And Fill the histograms with the flag information
1214  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1215  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1216  {
1217  //..Direction of amplitude (Checks energies from 0-10 GeV)
1218  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1219  {
1220  if(fFlag[cell]!=0)
1221  {
1222  //..cellID+1 = histogram bin
1223  cellAmp_masked->SetBinContent(amp,cell+1,0);
1224  }
1225  }
1226  //..Direction of time (Checks times from -275-975 GeV)
1227  for (Int_t time = 1; time <= fCellTime->GetNbinsX(); time++)
1228  {
1229  if(fFlag[cell]!=0)
1230  {
1231  //..cellID+1 = histogram bin
1232  cellTime_masked->SetBinContent(time,cell+1,0);
1233  }
1234  }
1235  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
1236  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
1237  }
1238 
1239  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1240  //..Plot the 2D distribution of cells by flag
1241  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1242  PlotFlaggedCells2D(0); //..all good cells
1243  PlotFlaggedCells2D(1); //..all dead cells
1244  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
1245  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
1246 
1247 
1248  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1249  //..Plot some summary canvases
1250  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1251  TCanvas *c1 = new TCanvas("CellProp","summary of cell properties",1000,1000);
1252  c1->ToggleEventStatus();
1253  c1->Divide(2,2);
1254  c1->cd(1)->SetLogz();
1255  //lowerPadRight->SetLeftMargin(0.09);
1256  //lowerPadRight->SetRightMargin(0.06);
1257  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1258  fCellAmplitude->SetYTitle("Abs. Cell Id");
1259  fCellAmplitude->GetYaxis()->SetTitleOffset(1.6);
1260  fCellAmplitude->Draw("colz");
1261  c1->cd(2)->SetLogz();
1262  fCellTime->SetXTitle("Cell Time [ns]");
1263  fCellTime->SetYTitle("Abs. Cell Id");
1264  fCellTime->GetYaxis()->SetTitleOffset(1.6);
1265  fCellTime->Draw("colz");
1266  c1->cd(3)->SetLogz();
1267  //lowerPadRight->SetLeftMargin(0.09);
1268  //lowerPadRight->SetRightMargin(0.06);
1269  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1270  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1271  cellAmp_masked->SetYTitle("Abs. Cell Id");
1272  cellAmp_masked->GetYaxis()->SetTitleOffset(1.6);
1273  cellAmp_masked->Draw("colz");
1274  c1->cd(4)->SetLogz();
1275  cellTime_masked->SetTitle("Masked Cell Time");
1276  cellTime_masked->SetXTitle("Cell Time [ns]");
1277  cellTime_masked->SetYTitle("Abs. Cell Id");
1278  cellTime_masked->GetYaxis()->SetTitleOffset(1.6);
1279  cellTime_masked->Draw("colz");
1280  c1->Update();
1281 
1282  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
1283  c2->ToggleEventStatus();
1284  c2->Divide(1,2);
1285  c2->cd(1);
1286  fhCellFlag->SetTitle("cell flag");
1287  fhCellFlag->SetXTitle("Abs. Cell Id");
1288  fhCellFlag->SetYTitle("flag by which cell was excluded");
1289  fhCellFlag->DrawCopy("hist");
1290  c2->cd(2);
1291  fhCellWarm->SetTitle("Warm cells");
1292  fhCellWarm->SetXTitle("Abs. Cell Id");
1293  fhCellWarm->SetYTitle("warm=1");
1294  fhCellWarm->DrawCopy("hist");
1295  c2->Update();
1296  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1297  //..Add different histograms/canvases to the output root file
1298  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1299  TString name =Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1300  c1->SaveAs(name);
1301  fRootFile->WriteObject(c1,c1->GetName());
1302  fRootFile->WriteObject(c2,c2->GetName());
1303  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
1304  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
1305  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
1306  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
1307  //..Save all amplitudes to the root file
1308  SaveHistoToFile();
1309 
1310  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1311  //..Save also the identified warm channels into a text file.
1312  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1313  file.open(cellSummaryFile, ios::out | ios::app);
1314  if(file)
1315  {
1316  file<<"Warm cells : "<<endl;
1317  cout<<" o Warm cells : "<<endl;
1318  nEMCalCells =0;
1319  nDCalCells =0;
1320  for(cellID=0; cellID<fNoOfCells; cellID++)
1321  {
1322  if(cellID==0)
1323  {
1324  file<<"In EMCal: "<<endl;
1325  }
1326  if(cellID==fCellStartDCal)
1327  {
1328  file<<"\n"<<endl;
1329  file<<"In DCal: "<<endl;
1330  }
1331  if(fWarmCell[cellID]==1)
1332  {
1333  file<<cellID<<", ";
1334  if(cellID<fCellStartDCal)nEMCalCells++;
1335  else nDCalCells++;
1336  }
1337  }
1338  file<<"\n"<<endl;
1339  perWarmEMCal= 100*nEMCalCells/(1.0*fCellStartDCal);
1340  perWarmDCal = 100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
1341  file<<"EMCal ("<<nEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1342  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
1343  }
1344  file.close();
1345  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1346  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1347  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1348  ofstream file3(aliceTwikiTable, ios::out | ios::app);
1349  if(file3)
1350  {
1351  file3<<"| Warm DCal | "<<nDCalCells<<" | "<<perWarmDCal<<"% |"<<endl;
1352  file3<<"\n"<<endl;
1353  file3<<"Insert above:"<<endl;
1354  file3<<"| Warm EMCal | "<<nEMCalCells<<" | "<<perWarmEMCal<<"% |"<<endl;
1355  }
1356  file3.close();
1357  //cout<<" o Results can be found in : "<<endl;
1358  //cout<<" o "<<cellSummaryFile<<endl;
1359  //cout<<" o "<<badPdfName<<endl;
1360  //cout<<" o "<<badPdfName<<endl;
1361 }
1362 
1363 
1379 //________________________________________________________________________
1381 {
1382  gROOT->SetStyle("Plain");
1383  gStyle->SetOptStat(0);
1384  gStyle->SetFillColor(kWhite);
1385  gStyle->SetTitleFillColor(kWhite);
1386  gStyle->SetPalette(1);
1387 
1388  char title[100];
1389  char name[100];
1390 
1391  TH1D *hRefDistr = BuildMeanFromGood();
1392  fRootFile->WriteObject(hRefDistr,hRefDistr->GetName());
1393  Int_t firstCanvas=0;
1394  Bool_t candidate;
1395  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
1396  text->SetTextSize(0.07);
1397  text->SetTextColor(kOrange-3);
1398  text->SetNDC();
1399 
1400  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
1401  text2->SetTextSize(0.08);
1402  text2->SetTextColor(8);
1403  text2->SetNDC();
1404 
1405  TLatex* textA = new TLatex(0.65,0.62,"*test*");
1406  textA->SetTextSize(0.04);
1407  textA->SetTextColor(1);
1408  textA->SetNDC();
1409 
1410  //..collect cells in an internal vector.
1411  //..when the vector is of size=9 or at the end of being filled
1412  //..plot the channels into a canvas
1413  std::vector<Int_t> channelVector;
1414  channelVector.clear();
1415  cout<<"Start printing into .pdf for version: "<<version<<endl;
1416  for(Int_t cell=0;cell<fNoOfCells;cell++)
1417  {
1418  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
1419  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
1420  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
1421  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
1422  if(fFlag[cell]==0 && version==20)channelVector.push_back(cell);
1423 
1424  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
1425  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1426  if(channelVector.size()==9 || cell == fNoOfCells-1)
1427  {
1428  cout<<"."<<flush;
1429 
1430  TString internal_pdfName=pdfName;
1431  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1432  if(channelVector.size() > 6) c1->Divide(3,3);
1433  else if (channelVector.size() > 3) c1->Divide(3,2);
1434  else c1->Divide(3,1);
1435 
1436  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1437  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
1438  {
1439  sprintf(name, "Cell %d",channelVector.at(i)) ;
1440  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
1441  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
1442 
1443  c1->cd(i%9 + 1);
1444  c1->cd(i%9 + 1)->SetLogy();
1445  if(i==0)
1446  {
1447  leg->AddEntry(hRefDistr,"mean of good","l");
1448  leg->AddEntry(hCell,"current","l");
1449  }
1450  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
1451  candidate = CheckDistribution(hCell,hRefDistr);
1452  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
1453  if(version>2)//..These are ratio plots of energy distr. of cell and mean of all good cells
1454  {
1455  hCell->Divide(hRefDistr);
1456  }
1457  //.. save histograms to file
1458  if(version==1) fOutputListBad->Add(hCell);
1459  if(version==10)fOutputListBadRatio->Add(hCell);
1460  if(version==2) fOutputListGood->Add(hCell);
1461  if(version==20)fOutputListGoodRatio->Add(hCell);
1462 
1463  hCell->SetLineColor(kBlue+1);
1464  hCell->GetXaxis()->SetTitle("E (GeV)");
1465  hCell->GetYaxis()->SetTitle("N Entries");
1466  hCell->GetXaxis()->SetRangeUser(0.,10.);
1467  hCell->SetLineWidth(1) ;
1468  hCell->SetTitle(title);
1469  hRefDistr->SetLineColor(kGray+2);
1470  hRefDistr->SetLineWidth(1);
1471 
1472  hCell->Draw("hist");
1473 
1474  if(version==1 || version==10)hRefDistr->Draw("same") ;
1475 
1476  //..Mark the histogram that could be miscalibrated and labelled as warm
1477  if(candidate==1 && (version==1 || version==10))
1478  {
1479  gPad->SetFrameFillColor(kYellow-10);
1480  text->Draw();
1481  }
1482  if(version==1)
1483  {
1484  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
1485  textA->Draw();
1486  }
1487  if(candidate==0 && (version==2 || version==20))
1488  {
1489  gPad->SetFrameFillColor(kYellow-10);
1490  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
1491  leg->Draw();
1492  }
1493  if(version<2)leg->Draw();
1494  }
1495 
1496  if(channelVector.size()<9 || cell == fNoOfCells-1)
1497  {
1498  internal_pdfName +=")";
1499  c1->Print(internal_pdfName.Data());
1500  }
1501  else if(firstCanvas==0)
1502  {
1503  internal_pdfName +="(";
1504  c1->Print(internal_pdfName.Data());
1505  firstCanvas=1;
1506  }
1507  else
1508  {
1509  c1->Print(internal_pdfName.Data());
1510  }
1511  delete c1;
1512  delete leg;
1513  channelVector.clear();
1514  }
1515  }
1516  delete hRefDistr;
1517  //..Add the subdirectories to the file
1518  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
1519  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
1520  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
1521  if(version==20)fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1522 
1523  cout<<endl;
1524 }
1528 //_________________________________________________________________________
1530 {
1531  TH1D* hGoodAmp;
1532  TH1D* hgoodMean;
1533  Int_t NrGood=0;
1534  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1535  {
1536  if(fFlag[cell]!=0)continue;
1537  NrGood++;
1538  if(NrGood==1)hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean",cell+1,cell+1);
1539  else
1540  {
1541  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1542  hgoodMean->Add(hGoodAmp);
1543  }
1544  }
1545  hgoodMean->Scale(1.0/NrGood);
1546 
1547  return hgoodMean;
1548 }
1559 //_________________________________________________________________________
1561 {
1562  TString Name = Form("%s_ratio",histogram->GetName());
1563  TH1* ratio=(TH1*)histogram->Clone(Name);
1564  ratio->Divide(reference);
1565 
1566  Double_t percentageOfLast=0.6;
1567  Double_t higherThanMean=5;
1568  Double_t highestRatio=1000;
1569  Double_t maxMean=10;
1570  Double_t minMean=0.1;
1571  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
1572 
1573  //..By default each cell is a candidate for recalibration
1574  Bool_t candidate=1;
1575  //..Find bin where reference has value 1, and the corresponding x-value
1576  Int_t binHeihgtOne = reference->FindLastBinAbove(1);
1577  Double_t binCentreHeightOne = reference->GetBinCenter(binHeihgtOne);
1578  Double_t thirdBinCentre = reference->GetBinCenter(3);
1579 
1580  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1581  //..Check the histogram
1582  //..Different checks to see whether the
1583  //..cell is really bad. Set candidate to 0.
1584 
1585  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1586  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
1587  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
1588  {
1589  candidate=0;
1590  }
1591 
1592  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1593  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
1594  //ELI check that crieteria carfully - seems to work but not shure about it
1595  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
1596  if(ratio->GetMaximum()>highestRatio)//
1597  {
1598  candidate=0;
1599  }
1600 
1601  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1602  //..check whether the ratio is much larger than 1
1603  //..calculate the mean in the relevant energy range
1604  Double_t mean=0;
1605  for(Int_t i=2;i<binHeihgtOne;i++)
1606  {
1607  mean+=ratio->GetBinContent(i);
1608  }
1609  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1610  if(mean>maxMean || mean<minMean)
1611  {
1612  candidate=0;
1613  }
1614 
1615  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1616  //..check whether there are large spikes in the histogram
1617  //..compare bin values to mean of the ratio. If there is a bin value with
1618  //..content "higherThanMean" times lareger than mean it's losing it candidate status
1619  mean=0;
1620  //..Find the maximum in the mean range (0-binHeihgtOne)
1621  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
1622  Double_t localMaxBin=ratio->GetMaximumBin();
1623 
1624  for(Int_t i=2;i<binHeihgtOne;i++)
1625  {
1626  //..Exclude 0 bins and exclude bins near the maximum
1627  if(ratio->GetBinContent(i)<=0) continue;
1628  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
1629  mean+=ratio->GetBinContent(i);
1630  }
1631  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1632  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
1633  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1634  if(ratio->GetMaximum()>mean*higherThanMean)
1635  {
1636  candidate=0;
1637  }
1638 
1639  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1640  //..Check for a cliff at 4GeV, happens in some cases
1641  Double_t beforeCliff=0;
1642  Double_t afterCliff=0;
1643  Int_t binsBefore=0;
1644  Int_t binsAfter=0;
1645  Int_t cliffBin = ratio->FindBin(4);
1646  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
1647  {
1648  //..Exclude 0 bins
1649  if(ratio->GetBinContent(i)<=0)continue;
1650  if(i<=cliffBin) binsBefore++;
1651  if(i>cliffBin) binsAfter++;
1652  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
1653  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
1654  beforeCliff*=1.0/binsBefore;
1655  afterCliff*=1.0/binsAfter;
1656  }
1657  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
1658  {
1659  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
1660  }
1661 
1662  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1663  //..Employ peak finder
1664 /* if(candidate==1)
1665  {
1666  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
1667  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
1668  //GetNPeaks();
1669  //cout<<"found N peaks: "<<nfound<<endl;
1670  }
1671 */
1672  return candidate;
1673 }
1674 
1683 //_________________________________________________________________________
1685 {
1686  //..TRD support structure
1687  //..(determined by eye, could be improved, but is already very acurate):
1688  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
1689  //..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
1690  Bool_t coveredByTRDSupportStruc=0;
1691 
1692  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
1693  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
1694  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
1695  )
1696  {
1697  coveredByTRDSupportStruc=1;
1698  }
1699  return coveredByTRDSupportStruc;
1700 }
1709 //_________________________________________________________________________
1711 {
1712  //..build two dimensional histogram with values row vs. column
1713  TString histoName;
1714  histoName = Form("2DChannelMap_Flag%d",flagBegin);
1715  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100");
1716 
1717  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
1718  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1719  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1720 
1721  Int_t cellColumn=0,cellRow=0;
1722  Int_t cellColumnAbs=0,cellRowAbs=0;
1723  Int_t trash;
1724 
1725  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1726  {
1727  //..Do that only for cell ids also accepted by the code
1728  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1729 
1730  //..Get Row and Collumn for cell ID c
1731  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1732  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1733  {
1734  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1735  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1736  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1737  }
1738  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1739  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1740  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
1741 
1742 
1743  }
1744  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
1745  c1->ToggleEventStatus();
1746  c1->cd();
1747  //lowerPadRight->SetLeftMargin(0.09);
1748  //lowerPadRight->SetRightMargin(0.06);
1749  plot2D->Draw("colz");
1750 
1751  TLatex* text = 0x0;
1752  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
1753  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
1754  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
1755  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
1756  text->SetTextSize(0.06);
1757  text->SetNDC();
1758  text->SetTextColor(1);
1759  text->Draw();
1760 
1761  c1->Update();
1762  TString name =Form("%s/%s/%s_%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data() , histoName.Data());
1763  c1->SaveAs(name);
1764 
1765  fRootFile->WriteObject(plot2D,plot2D->GetName());
1766 
1767 }
1771 //_________________________________________________________________________
1773 {
1774  char name[100];
1775  for(Int_t cell=0;cell<fNoOfCells;cell++)
1776  {
1777  sprintf(name, "Cell %d",cell) ;
1778  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
1779  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
1780  }
1781  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1782 }
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.
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
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 './'.
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)
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.
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 ...
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