AliPhysics  32b88a8 (32b88a8)
 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) TestCellShapes(criterion, emin, emax, nsigma);
581 
582  Double_t range=0.001;
583  if(emin>0.49)range=0.0005;
584  if(emin>0.99)range=0.0001;
585  if(emin>1.99)range=0.00005;
586 
587  if(criterion==1) FlagAsBad(criterion, histogram, nsigma, 200,-1);
588  if(criterion==2) FlagAsBad(criterion, histogram, nsigma, 600,range);
589 
590  /*
591  if(criterion==3) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval3);
592  if(criterion==4) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval1);
593  if(criterion==5) FlagAsBad(criterion, histogram, nsigma, dnbins, maxval2);
594  */
595 }
596 
606 //_________________________________________________________________________
608 {
609  cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
610  TH1F *histogram;
611  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);
612  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);
613  histogram->SetXTitle("Abs. Cell Id");
614  if(crit==1)histogram->SetYTitle("Energy per hit");
615  if(crit==2)histogram->SetYTitle("Number of hits per event");
616  histogram->GetXaxis()->SetNdivisions(510);
617  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
618 
619  //..here the average hit per event and the average energy per hit is caluclated for each cell.
620  for (Int_t cell = 0; cell < fNoOfCells; cell++)
621  {
622  Double_t Esum = 0;
623  Double_t Nsum = 0;
624 
625  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
626  {
627  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
628  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
629  //..investigate only cells that were not flagged as dead ore bad
630  if (E < emin || E > emax || fFlag[cell]!=0) continue;
631  Esum += E*N;
632  Nsum += N;
633  }
634  //..Set the values only for cells that are not yet marked as bad
635  if(fFlag[cell]==0)
636  {
637  if(totalevents> 0. && crit==2)histogram->SetBinContent(cell+1, Nsum/totalevents); //..number of hits per event
638  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/Nsum); //..average energy per hit
639  }
640  }
641  return histogram;
642 }
647 //_________________________________________________________________________
649 {
650  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);
651 
652  return Histogram;
653 }
664 //_________________________________________________________________________
666 {
667  //..binning parameters
668  Int_t ncells = fCellAmplitude->GetNbinsY();
669  Double_t amin = fCellAmplitude->GetYaxis()->GetXmin();
670  Double_t amax = fCellAmplitude->GetYaxis()->GetXmax();
671  cout << "ncells " << ncells << " amin = " << amin << "amax = " << amax<< endl;
672 
673  //..initialize histograms
674  TH1 *hFitA = new TH1F("hFitA_hCellAmplitude","Fit A value", ncells,amin,amax);
675  hFitA->SetXTitle("AbsId");
676  hFitA->SetYTitle("A");
677 
678  TH1 *hFitB = new TH1F("hFitB_hCellAmplitude","Fit B value", ncells,amin,amax);
679  hFitB->SetXTitle("AbsIdhname");
680  hFitB->SetYTitle("B");
681 
682  TH1 *hFitChi2Ndf = new TH1F("hFitChi2Ndf_hCellAmplitude","Fit #chi^{2}/ndf value", ncells,amin,amax);
683  hFitChi2Ndf->SetXTitle("AbsId");
684  hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
685 
686  Double_t maxval1=0., maxval2=0., maxval3=0.;
687  Double_t prev=0., MSA=0., AvA = 0. ; //those param are used to automaticaly determined a reasonable maxval1
688  Double_t prev2=0., MSB=0., AvB = 0. ; //those param are used to automaticaly determined a reasonable maxval2
689  Double_t prev3=0., MSki2=0., Avki2 = 0. ; //those param are used to automaticaly determined a reasonable maxval3
690  Double_t ki2=0.0 ;
691  for (Int_t k = 1; k <= fNoOfCells; k++)
692  {
693  TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
694  TH1 *hCell = fCellAmplitude->ProjectionX("",k,k);
695  if (hCell->GetEntries() == 0) continue;
696  // hCell->Rebin(3);
697  hCell->Fit(fit, "0QEM", "", fitemin, fitemax);
698  delete hCell;
699 
700  if(fit->GetParameter(0) < 5000.)
701  {
702  hFitA->SetBinContent(k, fit->GetParameter(0));
703  if(k<3000)
704  {
705  AvA += fit->GetParameter(0);
706  if(k==2999) maxval1 = AvA/3000. ;
707  if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
708  else MSA -= (fit->GetParameter(0) - prev) ;
709  prev = fit->GetParameter(0);
710  }
711  else
712  {
713  if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
714  {
715  maxval1 = fit->GetParameter(0);
716  }
717  }
718  }
719  else hFitA->SetBinContent(k, 5000.);
720 
721  if(fit->GetParameter(1) < 5000.)
722  {
723  hFitB->SetBinContent(k, fit->GetParameter(1));
724  if(k<3000)
725  {
726  AvB += fit->GetParameter(1);
727  if(k==2999) maxval2 = AvB/3000. ;
728  if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
729  else MSB -= (fit->GetParameter(1) - prev2) ;
730  prev2 = fit->GetParameter(1);
731  }
732  else
733  {
734  if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
735  {
736  maxval2 = fit->GetParameter(1);
737  }
738  }
739  }
740  else hFitB->SetBinContent(k, 5000.);
741 
742 
743  if (fit->GetNDF() != 0 ) ki2 = fit->GetChisquare()/fit->GetNDF();
744  else ki2 = 1000.;
745 
746  if(ki2 < 1000.)
747  {
748  hFitChi2Ndf->SetBinContent(k, ki2);
749  if(k<3000)
750  {
751  Avki2 += ki2;
752  if(k==2999) maxval3 = Avki2/3000. ;
753  if (prev3 < ki2) MSki2 += ki2 - prev3;
754  else MSki2 -= (ki2 - prev3) ;
755  prev3 = ki2;
756  }
757  else
758  {
759  if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
760  {
761  maxval3 = ki2;
762  }
763  }
764  }
765  else hFitChi2Ndf->SetBinContent(k, 1000.);
766 
767  delete fit ;
768  }
769 
770  // if you have problem with automatic parameter :
771  // maxval1 =
772  // maxval2 =
773  // maxval3 =
774  /*
775  if(crit==3)
776  Process(crit, hFitChi2Ndf, nsigma, dnbins, maxval3);
777  if(crit==4)
778  Process(crit, hFitA, nsigma, dnbins, maxval1);
779  if(crit==5)
780  Process(crit, hFitB, nsigma, dnbins, maxval2);
781  */
782 
783  //ELI something like thie in the future:
784  //return histogram;
785 }
786 
787 
792 //_________________________________________________________________________
794 {
795  Int_t sumOfExcl=0;
796  //..Direction of cell ID
797  for (Int_t cell = 0; cell < fNoOfCells; cell++)
798  {
799  Double_t nSum = 0;
800  //..Direction of amplitude (Checks energies from 0-10 GeV)
801  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
802  {
803  //..cellID+1 = histogram bin
804  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
805  nSum += N;
806  }
807  //..If the amplitude in one cell is basically 0
808  //..mark the cell as excluded
809  if(nSum < 0.5 && fFlag[cell]==0)
810  {
811  //..Cell flagged as dead.
812  //..Flag only if it hasn't been flagged before
813  fFlag[cell] =1;
814  sumOfExcl++;
815  }
816  }
817  cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
818  cout<<" ("<<sumOfExcl<<")"<<endl;
819 }
820 
835 //_________________________________________________________________________
836 void AliAnaCaloChannelAnalysis::FlagAsBad(Int_t crit, TH1F* inhisto, Double_t nsigma, Int_t dnbins, Double_t dmaxval)
837 {
838  gStyle->SetOptStat(0); // MG modif
839  gStyle->SetOptFit(0); // MG modif
840 
841  if(crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
842  if(crit==2)cout<<" o Fit average hit per event distribution"<<endl;
843 
844  Int_t cellColumn=0,cellRow=0;
845  Int_t cellColumnAbs=0,cellRowAbs=0;
846  Int_t trash;
847 
848  TString histoName=inhisto->GetName();
849  Double_t goodmax= 0. ;
850  Double_t goodmin= 0. ;
851  if (dmaxval < 0.)
852  {
853  dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
854  if(crit==2 && dmaxval > 1) dmaxval =1. ;
855  }
856 
857  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
858  //. . .build the distribution of average values
859  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
860  Double_t totalevents = fProcessedEvents->Integral(1, fProcessedEvents->GetNbinsX());
861  //..Be aware that the bin width should never be larger than 1/totalevents because this
862  //..is the minimum differce between the cells. One hit more/event otherwise you will see ugly empty bins in the histogram
863  if(((dmaxval-inhisto->GetMinimum())/(dnbins*1.0))<1.0/totalevents)
864  {
865  dnbins=(dmaxval-inhisto->GetMinimum())*totalevents;
866  cout<<"Problem - Reset dnbins to new value:"<<dnbins<<endl;
867  }
868 
869  //..build histos
870  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
871  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
872  distrib->SetYTitle("Entries");
873  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
874  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
875 
876 
877  //..build two dimensional histogram with values row vs. column
878  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);
879  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
880  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
881 
882  for (Int_t cell = 0; cell < fNoOfCells; cell++)
883  {
884  //..Do that only for cell ids also accepted by the code
885  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
886 
887  //..Fill histograms only for cells that are not yet flagged as bad
888  if(fFlag[cell]==0)
889  {
890  //..fill the distribution of avarge cell values
891  distrib->Fill(inhisto->GetBinContent(cell+1));
892  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
893  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
894  //..Get Row and Collumn for cell ID
895  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
896  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
897  {
898  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
899  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
900  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
901  }
902  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
903  //..check TRD support structure
904  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
905  {
906  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
907  }
908  else
909  {
910  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
911  }
912  }
913  }
914 
915  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
916  //. . .draw histogram + distribution
917  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
918 
919  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
920  c1->ToggleEventStatus();
921  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
922  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
923  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
924  upperPad->Draw();
925  lowerPadLeft->Draw();
926  lowerPadRight->Draw();
927 
928  upperPad->cd();
929  upperPad->SetLeftMargin(0.045);
930  upperPad->SetRightMargin(0.05);
931  upperPad->SetLogy();
932  inhisto->SetTitleOffset(0.6,"Y");
933  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
934 
935  inhisto->SetLineColor(kBlue+1);
936  inhisto->Draw();
937 
938  lowerPadRight->cd();
939  lowerPadRight->SetLeftMargin(0.09);
940  lowerPadRight->SetRightMargin(0.1);
941  plot2D->Draw("colz");
942 
943  lowerPadLeft->cd();
944  lowerPadLeft->SetLeftMargin(0.09);
945  lowerPadLeft->SetRightMargin(0.07);
946  lowerPadLeft->SetLogy();
947  distrib->SetLineColor(kBlue+1);
948  distrib->Draw();
949  distrib_wTRDStruc->SetLineColor(kGreen+1);
950  distrib_wTRDStruc->DrawCopy("same");
951  distrib_woTRDStruc->SetLineColor(kMagenta+1);
952  distrib_woTRDStruc->DrawCopy("same");
953 
954  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
955  //. . .fit histogram
956  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
957  Int_t higherbin=0,i;
958  for(i = 2; i <= dnbins; i++)
959  {
960  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i)) higherbin = i ;
961  }
962  //..good range is around the max value as long as the
963  //..bin content is larger than 2 entries
964  for(i = higherbin ; i<=dnbins ; i++)
965  {
966  if(distrib->GetBinContent(i)<2) break ;
967  goodmax = distrib->GetBinCenter(i);
968  }
969  for(i = higherbin ; i>1 ; i--)
970  {
971  if(distrib->GetBinContent(i)<2) break ;
972  goodmin = distrib->GetBinLowEdge(i);
973  }
974  //cout<<"higherbin : "<<higherbin<<endl;
975  //cout<<"good range : "<<goodmin<<" - "<<goodmax<<endl;
976 
977  TF1 *fit2 = new TF1("fit2", "gaus",0,10);
978  //..start the fit with a mean of the highest value
979  fit2->SetParameter(1,higherbin);
980 
981  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
982  Double_t sig, mean;// chi2ndf;
983  // Marie midif to take into account very non gaussian distrig
984  mean = fit2->GetParameter(1);
985  sig = fit2->GetParameter(2);
986  //chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
987 
988  if (mean <0.) mean=0.; //ELI is this not a highly problematic case??
989 
990  goodmin = mean - nsigma*sig ;
991  goodmax = mean + nsigma*sig ;
992 
993  if (goodmin<0) goodmin=0.;
994 
995  cout<<" o Result of fit: "<<endl;
996  cout<<" o "<<endl;
997  cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
998  cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
999 
1000  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1001  //. . .Add info to histogram
1002  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1003  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1004  lline->SetLineColor(kGreen+2);
1005  lline->SetLineStyle(7);
1006  lline->Draw();
1007 
1008  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1009  rline->SetLineColor(kGreen+2);
1010  rline->SetLineStyle(7);
1011  rline->Draw();
1012 
1013  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1014  leg->AddEntry(lline,"Good region boundary","l");
1015  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1016  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1017  leg->SetBorderSize(0);
1018  leg->Draw("same");
1019 
1020  fit2->SetLineColor(kOrange-3);
1021  fit2->SetLineStyle(1);//7
1022  fit2->Draw("same");
1023 
1024  TLatex* text = 0x0;
1025  if(crit==1) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1026  if(crit==2) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2fx10^-5",goodmin*100000,goodmax*100000));
1027  text->SetTextSize(0.06);
1028  text->SetNDC();
1029  text->SetTextColor(1);
1030  text->Draw();
1031 
1032 
1033  upperPad->cd();
1034  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
1035  uline->SetLineColor(kGreen+2);
1036  uline->SetLineStyle(7);
1037  uline->Draw();
1038  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
1039  lowline->SetLineColor(kGreen+2);
1040  lowline->SetLineStyle(7);
1041  lowline->Draw();
1042  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1043  //. . .Save histogram
1044  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1045  c1->Update();
1046  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
1047  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1048  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1049  c1->SaveAs(name);
1050 
1051  fRootFile->WriteObject(c1,c1->GetName());
1052  fRootFile->WriteObject(plot2D,plot2D->GetName());
1053  fRootFile->WriteObject(distrib,distrib->GetName());
1054  fRootFile->WriteObject(inhisto,inhisto->GetName());
1055  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1056  //. . . Mark the bad cells in the fFlag array
1057  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1058  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1059  //. . .(0 by default - good cell)
1060  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1061  cout<<" o Flag bad cells that are outside the good range "<<endl;
1062  for(Int_t cell = 0; cell < fNoOfCells; cell++)
1063  {
1064  //..cell=0 and bin=1, cell=1 and bin=2
1065  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1066  if (inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)//ELI
1067  {
1068  fFlag[cell]=fCriterionCounter;
1069  }
1070  if (inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
1071  {
1072  fFlag[cell]=fCriterionCounter;
1073  }
1074  }
1075  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;
1076 }
1077 
1078 
1079 
1080 
1085 //________________________________________________________________________
1087 {
1088  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1089  //.. RESULTS
1090  //.. 1) Print the bad cells
1091  //.. and write the results to a file
1092  //.. for each added period analysis
1093  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1094  TArrayD periodArray;
1095  Double_t emin,emax;
1096  Int_t criterion;
1097  TString output;
1098  Int_t nb1=0;
1099 
1100  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
1101  {
1102  periodArray=fAnalysisVector.at(i);
1103  criterion =periodArray.At(0);
1104  emin =periodArray.At(2);
1105  emax =periodArray.At(3);
1106 
1107  //..Print the results on the screen and
1108  //..write the results in a file
1109  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
1110  ofstream file(output, ios::out | ios::trunc);
1111  if(!file)
1112  {
1113  cout<<"#### Major Error. Check the textfile!"<<endl;
1114  }
1115  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1116  cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1117 
1118  nb1=0;
1119  for(Int_t cellID=0;cellID<fNoOfCells;cellID++)
1120  {
1121  if(fFlag[cellID]==(i+2))
1122  {
1123  nb1++;
1124  file<<cellID<<", ";
1125  }
1126  }
1127  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1128  file<<"("<<nb1<<")"<<endl;
1129  file.close();
1130  cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1131  cout<<endl;
1132  }
1133 }
1140 //________________________________________________________________________
1142 {
1143  Int_t cellID, nDCalCells = 0, nEMCalCells = 0;
1144  TString cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells;
1145  TH2F* cellAmp_masked= (TH2F*)fCellAmplitude->Clone("cellAmp_masked");
1146  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1147 
1148  deadPdfName = Form("%s/%s/%s%s_Dead_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1149  badPdfName = Form("%s/%s/%s%s_Bad_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1150  cellSummaryFile = Form("%s/%s/%s%s_Bad_Amplitudes_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial); ;
1151  ratioOfBad = Form("%s/%s/%s%s_BCRatio_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1152  goodCells = Form("%s/%s/%s%s_Good_Amplitudes_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fPass.Data(), fTrigger.Data() ,fTrial);
1153 
1154  cout<<" o Final results o "<<endl;
1155 
1156  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1157  //..Write the final results of dead and bad cells in a file and on screen
1158  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1159  ofstream file(cellSummaryFile, ios::out | ios::trunc);
1160  if(file)
1161  {
1162  file<<"Dead cells : "<<endl;
1163  cout<<" o Dead cells : "<<endl;
1164  nEMCalCells =0;
1165  nDCalCells =0;
1166  for(cellID=0; cellID<fNoOfCells; cellID++)
1167  {
1168  if(cellID==0)
1169  {
1170  file<<"In EMCal: "<<endl;
1171  //cout<<" o In EMCal : "<<endl;
1172  }
1173  if(cellID==fCellStartDCal)
1174  {
1175  file<<"\n"<<endl;
1176  file<<"In DCal: "<<endl;
1177  //cout<<endl;
1178  //cout<<" o In DCal : "<<endl;
1179  }
1180  if(fFlag[cellID]==1)
1181  {
1182  //file<<cellID<<"\n" ;
1183  file<<cellID<<", ";
1184  //cout<<cellID<<"," ;
1185  if(cellID<fCellStartDCal)nEMCalCells++;
1186  else nDCalCells++;
1187  }
1188  }
1189  //cout<<endl;
1190  file<<"\n"<<endl;
1191  file<<"EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1192  cout<<" o EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1193 
1194  file<<"Bad cells: "<<endl;
1195  cout<<" o Bad cells: "<<endl;
1196  nEMCalCells =0;
1197  nDCalCells =0;
1198  for(cellID=0;cellID<fNoOfCells;cellID++)
1199  {
1200  if(cellID==0)
1201  {
1202  file<<"In EMCal: "<<endl;
1203  //cout<<" o In EMCal : "<<endl;
1204  }
1205  if(cellID==fCellStartDCal)
1206  {
1207  file<<"\n"<<endl;
1208  file<<"In DCal: "<<endl;
1209  //cout<<endl;
1210  //cout<<" o In DCal : "<<endl;
1211  }
1212  if(fFlag[cellID]>1)
1213  {
1214  //file<<cellID<<"\n" ;
1215  file<<cellID<<", ";
1216  //cout<<cellID<<"," ;
1217  if(cellID<fCellStartDCal)nEMCalCells++;
1218  else nDCalCells++;
1219  }
1220  }
1221  //cout<<endl;
1222  file<<"\n"<<endl;
1223  file<<"EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1224  cout<<" o EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1225  }
1226  file.close();
1227 
1228  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1229  //..Save the flagged cells to .pdf files
1230  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1231  cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1232  SaveBadCellsToPDF(1,badPdfName) ;
1233  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
1234  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
1235 
1236 
1237  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1238  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1239  //..And Fill the histograms with the flag information
1240  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1241  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1242  {
1243  //..Direction of amplitude (Checks energies from 0-10 GeV)
1244  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1245  {
1246  if(fFlag[cell]!=0)
1247  {
1248  //..cellID+1 = histogram bin
1249  cellAmp_masked->SetBinContent(amp,cell+1,0);
1250  cellTime_masked->SetBinContent(amp,cell+1,0);
1251  }
1252  }
1253  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
1254  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
1255  }
1256 
1257  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1258  //..Plot the 2D distribution of cells by flag
1259  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1260  PlotFlaggedCells2D(0); //..all good cells
1261  PlotFlaggedCells2D(1); //..all dead cells
1262  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
1263  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
1264 
1265 
1266  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1267  //..Plot some summary canvases
1268  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1269  TCanvas *c1 = new TCanvas("CellProp","summary of cell properties",1000,1000);
1270  c1->ToggleEventStatus();
1271  c1->Divide(2,2);
1272  c1->cd(1)->SetLogz();
1273  //lowerPadRight->SetLeftMargin(0.09);
1274  //lowerPadRight->SetRightMargin(0.06);
1275  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1276  fCellAmplitude->SetYTitle("Abs. Cell Id");
1277  fCellAmplitude->Draw("colz");
1278  c1->cd(2)->SetLogz();
1279  fCellTime->SetXTitle("Cell Time [ns]");
1280  fCellTime->SetYTitle("Abs. Cell Id");
1281  fCellTime->Draw("colz");
1282  c1->cd(3)->SetLogz();
1283  //lowerPadRight->SetLeftMargin(0.09);
1284  //lowerPadRight->SetRightMargin(0.06);
1285  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1286  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1287  cellAmp_masked->SetYTitle("Abs. Cell Id");
1288  cellAmp_masked->Draw("colz");
1289  c1->cd(4)->SetLogz();
1290  cellTime_masked->SetTitle("Masked Cell Time");
1291  cellTime_masked->SetXTitle("Cell Time [ns]");
1292  cellTime_masked->SetYTitle("Abs. Cell Id");
1293  cellTime_masked->Draw("colz");
1294  c1->Update();
1295 
1296  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
1297  c2->ToggleEventStatus();
1298  c2->Divide(1,2);
1299  c2->cd(1);
1300  fhCellFlag->SetTitle("cell flag");
1301  fhCellFlag->SetXTitle("Abs. Cell Id");
1302  fhCellFlag->SetYTitle("flag");
1303  fhCellFlag->DrawCopy("hist");
1304  c2->cd(2);
1305  fhCellWarm->SetTitle("Warm cells");
1306  fhCellWarm->SetXTitle("Abs. Cell Id");
1307  fhCellWarm->SetYTitle("warm=1");
1308  fhCellWarm->DrawCopy("hist");
1309  c2->Update();
1310  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1311  //..Add different histograms/canvases to the output root file
1312  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1313  TString name =Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
1314  c1->SaveAs(name);
1315  fRootFile->WriteObject(c1,c1->GetName());
1316  fRootFile->WriteObject(c2,c2->GetName());
1317  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
1318  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
1319  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
1320  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
1321  //..Save all amplitudes to the root file
1322  SaveHistoToFile();
1323 
1324  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1325  //..Save also the identified warm channels into a text file.
1326  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1327  file.open(cellSummaryFile, ios::out | ios::app);
1328  if(file)
1329  {
1330  file<<"Warm cells : "<<endl;
1331  cout<<" o Warm cells : "<<endl;
1332  nEMCalCells =0;
1333  nDCalCells =0;
1334  for(cellID=0; cellID<fNoOfCells; cellID++)
1335  {
1336  if(cellID==0)
1337  {
1338  file<<"In EMCal: "<<endl;
1339  }
1340  if(cellID==fCellStartDCal)
1341  {
1342  file<<"\n"<<endl;
1343  file<<"In DCal: "<<endl;
1344  }
1345  if(fWarmCell[cellID]==1)
1346  {
1347  file<<cellID<<", ";
1348  if(cellID<fCellStartDCal)nEMCalCells++;
1349  else nDCalCells++;
1350  }
1351  }
1352  file<<"\n"<<endl;
1353  file<<"EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1354  cout<<" o EMCal ("<<nEMCalCells<<" ="<<100*nEMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<nDCalCells<<" ="<<100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1355  }
1356  file.close();
1357 
1358  //cout<<" o Results can be found in : "<<endl;
1359  //cout<<" o "<<cellSummaryFile<<endl;
1360  //cout<<" o "<<badPdfName<<endl;
1361  //cout<<" o "<<badPdfName<<endl;
1362 }
1363 
1364 
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 
1404  TLatex* textA = new TLatex(0.65,0.62,"*test*");
1405  textA->SetTextSize(0.04);
1406  textA->SetTextColor(1);
1407  textA->SetNDC();
1408 
1409  //..collect cells in an internal vector.
1410  //..when the vector is of size=9 or at the end of being filled
1411  //..plot the channels into a canvas
1412  std::vector<Int_t> channelVector;
1413  channelVector.clear();
1414  cout<<"Start printing into .pdf for version: "<<version<<endl;
1415  for(Int_t cell=0;cell<fNoOfCells;cell++)
1416  {
1417  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
1418  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
1419  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
1420  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
1421 
1422  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
1423  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1424  if(channelVector.size()==9 || cell == fNoOfCells-1)
1425  {
1426  cout<<"."<<flush;
1427 
1428  TString internal_pdfName=pdfName;
1429  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1430  if(channelVector.size() > 6) c1->Divide(3,3);
1431  else if (channelVector.size() > 3) c1->Divide(3,2);
1432  else c1->Divide(3,1);
1433 
1434  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1435  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
1436  {
1437  sprintf(name, "Cell %d",channelVector.at(i)) ;
1438  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
1439  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
1440 
1441  c1->cd(i%9 + 1);
1442  c1->cd(i%9 + 1)->SetLogy();
1443  if(i==0)
1444  {
1445  leg->AddEntry(hRefDistr,"mean of good","l");
1446  leg->AddEntry(hCell,"current","l");
1447  }
1448  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
1449  candidate = CheckDistribution(hCell,hRefDistr);
1450  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
1451  if(version>1)//..These are ratio plots of energy distr. of cell and mean of all good cells
1452  {
1453  hCell->Divide(hRefDistr);
1454  }
1455  //.. save histograms to file
1456  if(version==1) fOutputListBad->Add(hCell);
1457  if(version==10)fOutputListBadRatio->Add(hCell);
1458  if(version==2) fOutputListGoodRatio->Add(hCell);
1459 
1460  hCell->SetLineColor(kBlue+1);
1461  hCell->GetXaxis()->SetTitle("E (GeV)");
1462  hCell->GetYaxis()->SetTitle("N Entries");
1463  hCell->GetXaxis()->SetRangeUser(0.,10.);
1464  hCell->SetLineWidth(1) ;
1465  hCell->SetTitle(title);
1466  hRefDistr->SetLineColor(kGray+2);
1467  hRefDistr->SetLineWidth(1);
1468 
1469  hCell->Draw("hist");
1470 
1471  if(version==1)hRefDistr->Draw("same") ;
1472 
1473  //..Mark the histogram that could be miscalibrated and labelled as warm
1474  if(candidate==1 && (version==1 || version==10))
1475  {
1476  gPad->SetFrameFillColor(kYellow-10);
1477  text->Draw();
1478  }
1479  if(version==1)
1480  {
1481  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
1482  textA->Draw();
1483  }
1484  if(version==2 && candidate==0)
1485  {
1486  gPad->SetFrameFillColor(kYellow-10);
1487  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
1488  leg->Draw();
1489  }
1490  if(version<2)leg->Draw();
1491  }
1492 
1493  if(channelVector.size()<9 || cell == fNoOfCells-1)
1494  {
1495  internal_pdfName +=")";
1496  c1->Print(internal_pdfName.Data());
1497  }
1498  else if(firstCanvas==0)
1499  {
1500  internal_pdfName +="(";
1501  c1->Print(internal_pdfName.Data());
1502  firstCanvas=1;
1503  }
1504  else
1505  {
1506  c1->Print(internal_pdfName.Data());
1507  }
1508  delete c1;
1509  delete leg;
1510  channelVector.clear();
1511  }
1512  }
1513  delete hRefDistr;
1514  //..Add the subdirectories to the file
1515  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
1516  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
1517  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
1518 
1519  cout<<endl;
1520 }
1524 //_________________________________________________________________________
1526 {
1527  TH1D* hGoodAmp;
1528  TH1D* hgoodMean;
1529  Int_t NrGood=0;
1530  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1531  {
1532  if(fFlag[cell]!=0)continue;
1533  NrGood++;
1534  if(NrGood==1)hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean",cell+1,cell+1);
1535  else
1536  {
1537  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1538  hgoodMean->Add(hGoodAmp);
1539  }
1540  }
1541  hgoodMean->Scale(1.0/NrGood);
1542 
1543  return hgoodMean;
1544 }
1555 //_________________________________________________________________________
1557 {
1558  TString Name = Form("%s_ratio",histogram->GetName());
1559  TH1* ratio=(TH1*)histogram->Clone(Name);
1560  ratio->Divide(reference);
1561 
1562  Double_t percentageOfLast=0.6;
1563  Double_t higherThanMean=5;
1564  Double_t highestRatio=1000;
1565  Double_t maxMean=10;
1566  Double_t minMean=0.1;
1567  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
1568 
1569  //..By default each cell is a candidate for recalibration
1570  Bool_t candidate=1;
1571  //..Find bin where reference has value 1, and the corresponding x-value
1572  Int_t binHeihgtOne = reference->FindLastBinAbove(1);
1573  Double_t binCentreHeightOne = reference->GetBinCenter(binHeihgtOne);
1574  Double_t thirdBinCentre = reference->GetBinCenter(3);
1575 
1576  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1577  //..Check the histogram
1578  //..Different checks to see whether the
1579  //..cell is really bad. Set candidate to 0.
1580 
1581  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1582  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
1583  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
1584  {
1585  candidate=0;
1586  }
1587 
1588  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1589  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
1590  //ELI check that crieteria carfully - seems to work but not shure about it
1591  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
1592  if(ratio->GetMaximum()>highestRatio)//
1593  {
1594  candidate=0;
1595  }
1596 
1597  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1598  //..check whether the ratio is much larger than 1
1599  //..calculate the mean in the relevant energy range
1600  Double_t mean=0;
1601  for(Int_t i=2;i<binHeihgtOne;i++)
1602  {
1603  mean+=ratio->GetBinContent(i);
1604  }
1605  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1606  if(mean>maxMean || mean<minMean)
1607  {
1608  candidate=0;
1609  }
1610 
1611  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1612  //..check whether there are large spikes in the histogram
1613  //..compare bin values to mean of the ratio. If there is a bin value with
1614  //..content "higherThanMean" times lareger than mean it's losing it candidate status
1615  mean=0;
1616  //..Find the maximum in the mean range (0-binHeihgtOne)
1617  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
1618  Double_t localMaxBin=ratio->GetMaximumBin();
1619 
1620  for(Int_t i=2;i<binHeihgtOne;i++)
1621  {
1622  //..Exclude 0 bins and exclude bins near the maximum
1623  if(ratio->GetBinContent(i)<=0) continue;
1624  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
1625  mean+=ratio->GetBinContent(i);
1626  }
1627  mean*=1.0/(binHeihgtOne-1);//..divide by number of bins
1628  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
1629  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1630  if(ratio->GetMaximum()>mean*higherThanMean)
1631  {
1632  candidate=0;
1633  }
1634 
1635  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1636  //..Check for a cliff at 4GeV, happens in some cases
1637  Double_t beforeCliff=0;
1638  Double_t afterCliff=0;
1639  Int_t binsBefore=0;
1640  Int_t binsAfter=0;
1641  Int_t cliffBin = ratio->FindBin(4);
1642  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
1643  {
1644  //..Exclude 0 bins
1645  if(ratio->GetBinContent(i)<=0)continue;
1646  if(i<=cliffBin) binsBefore++;
1647  if(i>cliffBin) binsAfter++;
1648  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
1649  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
1650  beforeCliff*=1.0/binsBefore;
1651  afterCliff*=1.0/binsAfter;
1652  }
1653  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
1654  {
1655  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
1656  }
1657 
1658  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1659  //..Employ peak finder
1660 /* if(candidate==1)
1661  {
1662  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
1663  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
1664  //GetNPeaks();
1665  //cout<<"found N peaks: "<<nfound<<endl;
1666  }
1667 */
1668  return candidate;
1669 }
1670 
1679 //_________________________________________________________________________
1681 {
1682  //..TRD support structure
1683  //..(determined by eye, could be improved, but is already very acurate):
1684  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
1685  //..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
1686  Bool_t coveredByTRDSupportStruc=0;
1687 
1688  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
1689  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
1690  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
1691  )
1692  {
1693  coveredByTRDSupportStruc=1;
1694  }
1695  return coveredByTRDSupportStruc;
1696 }
1705 //_________________________________________________________________________
1707 {
1708  //..build two dimensional histogram with values row vs. column
1709  TString histoName;
1710  histoName = Form("2DChannelMap_Flag%d",flagBegin);
1711  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100");
1712 
1713  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
1714  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1715  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1716 
1717  Int_t cellColumn=0,cellRow=0;
1718  Int_t cellColumnAbs=0,cellRowAbs=0;
1719  Int_t trash;
1720 
1721  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1722  {
1723  //..Do that only for cell ids also accepted by the code
1724  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1725 
1726  //..Get Row and Collumn for cell ID c
1727  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1728  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1729  {
1730  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1731  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1732  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1733  }
1734  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1735  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
1736  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
1737 
1738 
1739  }
1740  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
1741  c1->ToggleEventStatus();
1742  c1->cd();
1743  //lowerPadRight->SetLeftMargin(0.09);
1744  //lowerPadRight->SetRightMargin(0.06);
1745  plot2D->Draw("colz");
1746 
1747  TLatex* text = 0x0;
1748  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
1749  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
1750  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
1751  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
1752  text->SetTextSize(0.06);
1753  text->SetNDC();
1754  text->SetTextColor(1);
1755  text->Draw();
1756 
1757  c1->Update();
1758  TString name =Form("%s/%s/%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(), histoName.Data());
1759  c1->SaveAs(name);
1760 
1761  fRootFile->WriteObject(plot2D,plot2D->GetName());
1762 
1763 }
1767 //_________________________________________________________________________
1769 {
1770  char name[100];
1771  for(Int_t cell=0;cell<fNoOfCells;cell++)
1772  {
1773  sprintf(name, "Cell %d",cell) ;
1774  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
1775  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
1776  }
1777  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
1778 }
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.
void TestCellShapes(Int_t crit, Double_t fitemin, Double_t fitemax, Double_t nsigma=4.)
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.
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