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