AliPhysics  9b6b435 (9b6b435)
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 #include <TVectorD.h>
19 
20 #include <Riostream.h>
21 #include <stdio.h>
22 #include <fstream>
23 #include <iostream>
24 #include <alloca.h>
25 #include <string>
26 #include <cstring>
27 
28 // --- ANALYSIS system ---
29 #include <AliCalorimeterUtils.h>
30 #include <AliEMCALGeometry.h>
31 #include <BadChannelAna.h>
32 #include <AliAODEvent.h>
33 
34 using std::vector;
35 using std::cout;
36 using std::endl;
37 using std::ofstream;
38 using std::flush;
39 using std::ios;
40 
42 ClassImp(BadChannelAna);
44 
48 //________________________________________________________________________
50  TObject(),
51  fCurrentRunNumber(-1),
52  fPeriod(),
53  fTrainNo(),
54  fTrigger(),
55  fNoOfCells(0),
56  fCellStartDCal(12288),
57  fStartCell(0),
58  fEndLowerBound(1),
59  fAnalysisOutput(),
60  fAnalysisInput(),
61  fRunList(),
62  fWorkdir(),
63  fQADirect(),
64  fMergedFileName(),
65  fAnalysisVector(0x0),
66  fTrial(),
67  fExternalFileName(""),
68  fExternalBadMapName(""),
69  fTestRoutine(0),
70  fPrint(0),
71  fTrackCellRecord(0),
72  fNMaxCols(48),
73  fNMaxRows(24),
74  fNMaxColsAbs(),
75  fNMaxRowsAbs(),
76  fFlag(0x0),
77  fCriterionCounter(0x0),
78  fWarmCell(0x0),
79  fCaloUtils(),
80  fRootFile(0x0),
81  fCellAmplitude(),
82  fCellTime(),
83  fProcessedEvents(),
84  fhCellFlag(0x0),
85  fhCellWarm(0x0)
86 {
87  fCurrentRunNumber = 254381;
88  fPeriod = "LHC16h";
89  fTrainNo = "muon_caloLego";
90  fTrigger = "AnyINT";
91  fWorkdir = ".";
92  fRunListFileName = "runList.txt";
93  fTrial = 0;
94  fRunBRunMap = 0;
95 
96  Init();
97  PrintCellInfo(0);
98 }
99 
103 //________________________________________________________________________
104 BadChannelAna::BadChannelAna(TString period, TString train, TString trigger, Int_t runNumber,Int_t trial, TString workDir, TString listName, Bool_t runByRun):
105  TObject(),
106  fCurrentRunNumber(-1),
107  fPeriod(),
108  fTrainNo(),
109  fTrigger(),
110  fNoOfCells(0),
111  fCellStartDCal(12288),
112  fStartCell(0),
113  fEndLowerBound(1),
114  fAnalysisOutput(),
115  fAnalysisInput(),
116  fRunList(),
118  fWorkdir(),
119  fQADirect(),
120  fMergedFileName(),
121  fAnalysisVector(0x0),
122  fTrial(),
123  fExternalFileName(""),
125  fTestRoutine(0),
126  fPrint(0),
127  fTrackCellRecord(0),
128  fNMaxCols(48),
129  fNMaxRows(24),
130  fNMaxColsAbs(),
131  fNMaxRowsAbs(),
132  fFlag(0x0),
133  fCriterionCounter(0x0),
134  fWarmCell(0x0),
135  fCaloUtils(),
136  fRootFile(0x0),
137  fCellAmplitude(),
138  fCellTime(),
140  fhCellFlag(0x0),
141  fhCellWarm(0x0)
142 {
143  fCurrentRunNumber = runNumber;
144  fPeriod = period;
145  fTrainNo = train; //only for folder structure
146  fTrigger = trigger; //important to select trigger in output file == different wagons in lego train
147  fWorkdir = workDir;
148  fRunListFileName = listName;
149  fTrial = trial;
150  fRunBRunMap = runByRun; //..will create a folder with a compact storage of run-by-run output
151 
152  Init();
153  PrintCellInfo(0);
154 }
155 
159 //________________________________________________________________________
161 {
162  gROOT->ProcessLine("gErrorIgnoreLevel = kWarning;"); //Does not work -
163  // fPrint = 1;
164  //......................................................
165  //..Default values - can be set by functions
166  fTestRoutine=0;
167 
168  //..Settings for the input/output structure (hard coded)
169  // TO BE CHANGED
170  fAnalysisInput =Form("AnalysisInput/%s",fPeriod.Data());
171  fAnalysisOutput =Form("AnalysisOutput/%s/%s/Version%i",fPeriod.Data(),fTrainNo.Data(),fTrial);
172  if(fRunBRunMap)fAnalysisOutputRbR =Form("AnalysisOutput/%s/%sRbR/Version%i",fPeriod.Data(),fTrainNo.Data(),fTrial);
173 
174  //..Make output directory if it doesn't exist
175  //..first the general output folder, in case you run this analysis for the very first time
176  gSystem->mkdir(Form("%s/AnalysisOutput",fWorkdir.Data()));
177  //..first the period folder
178  gSystem->mkdir(Form("%s/AnalysisOutput/%s",fWorkdir.Data(),fPeriod.Data()));
179  //..then the Train folder
180  if(fRunBRunMap)gSystem->mkdir(Form("%s/AnalysisOutput/%s/%sRbR",fWorkdir.Data(),fPeriod.Data(),fTrainNo.Data()));
181  gSystem->mkdir(Form("%s/AnalysisOutput/%s/%s",fWorkdir.Data(),fPeriod.Data(),fTrainNo.Data()));
182  //..then the version folder
183  if(fRunBRunMap)gSystem->mkdir(Form("%s/%s",fWorkdir.Data(),fAnalysisOutputRbR.Data()));
184  gSystem->mkdir(Form("%s/%s",fWorkdir.Data(),fAnalysisOutput.Data()));
185 
186  fMergedFileName= Form("%s/%s/%s/MergedRuns_%s.root",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),fTrigger.Data());
187  fRunList = Form("%s/%s/%s/%s",fWorkdir.Data(), fAnalysisInput.Data(), fTrainNo.Data(), fRunListFileName.Data());
188  fQADirect = Form("CaloQA_%s",fTrigger.Data());
189 
190  TString fileName = Form("%s/%s/%s_%s_Histograms_V%i.root",fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data(),fTrigger.Data() ,fTrial);
191  fRootFile = new TFile(fileName,"recreate");
192  //.. make sure the vector is empty
193  fAnalysisVector.clear();
194  fSMMask.clear();
195  //......................................................
196  //..Initialize EMCal/DCal geometry
198  //..Create a dummy event for the CaloUtils
199  AliAODEvent* aod = new AliAODEvent();
202 
203  fNoOfCells =fCaloUtils->GetEMCALGeometry()->GetNCells(); //..Very important number, never change after that point!
204  Int_t nModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
205  //..These are the first cell IDs of each SM and a cell ID of a nonExsting SM20 to mark the end (17664)
206  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};
207  memcpy (fStartCellSM, array_StartCellSM_Value, sizeof (fStartCellSM));
208 
209  //..This is how the calorimeter looks like in the current period (defined by example run ID fCurrentRunNumber)
210  cout<<"Called geometry for run number: "<<fCurrentRunNumber<<endl;
211  cout<<"Number of supermod: "<<nModules<<endl;
212  cout<<"Number of cells: "<<fNoOfCells<<endl;
213  cout<<"Cell ID of first DCal cell: "<<fCellStartDCal<<endl;
214  //cout<<"Number of supermod utils: "<<fCaloUtils->GetNumberOfSuperModulesUsed()<<endl; //..will always be 22 unless set by hand
215 
216  //......................................................
217  //..Initialize flag array to store how the cell is categorized
218  //..In the histogram: bin 1= cellID 0, bin 2= cellID 1 etc
219  //..In the array: fFlag[cellID]= some information
220  fFlag.clear();
221  fFlag.reserve(fNoOfCells);
222  fFlag.resize(fNoOfCells);
223  std::fill(fFlag.begin(), fFlag.end(), 0); //..flagged as good by default
224  fWarmCell.clear();
225  fWarmCell.reserve(fNoOfCells);
226  fWarmCell.resize(fNoOfCells);
227  std::fill(fWarmCell.begin(), fWarmCell.end(), 0);//..flagged as not warm by default
228 
229  fCriterionCounter=2; //This value will be written in fflag and updates after each PeriodAnalysis
230  //......................................................
231  //..setings for the 2D histogram
232  //fNMaxCols = 48; //eta direction
233  //fNMaxRows = 24; //phi direction
235  fNMaxRowsAbs = Int_t (nModules/2)*fNMaxRows; //multiply by number of supermodules
236 
237  //......................................................
238  //..Create TLists to organize the outputfile
239  fOutputListBad = new TList();
240  fOutputListGood = new TList();
241  fOutputListBadRatio = new TList();
242  fOutputListGoodRatio = new TList();
243 
244  fOutputListBad ->SetName("BadCell_Amplitudes");
245  fOutputListGood ->SetName("GoodCell_Amplitudes");
246  fOutputListBadRatio ->SetName("BadCell_AmplitudeRatios");
247  fOutputListGoodRatio->SetName("GoodCell_AmplitudeRatios");
248 
249  fOutputListBad ->SetOwner();
250  fOutputListGood ->SetOwner();
251  fOutputListBadRatio ->SetOwner();
252  fOutputListGoodRatio ->SetOwner();
253 
254  //......................................................
255  //..Create Histograms to store the flag in a root file
256  fhCellFlag = new TH1F("fhCellFlag","fhCellFlag",fNoOfCells+10,0,fNoOfCells+10); //..cellID+1 = histogram bin
257  fhCellWarm = new TH1F("fhCellWarm","fhCellWarm",fNoOfCells+10,0,fNoOfCells+10); //..cellID+1 = histogram bin
258 }
263 {
264 
265  fAnalysisVector.clear();
266  fFlag.clear();
267  fWarmCell.clear();
268  if(fCaloUtils) delete fCaloUtils;
269  if(fOutputListBad) delete fOutputListBad;
270  if(fOutputListGood) delete fOutputListGood;
273  if(fhCellFlag) delete fhCellFlag;
274  if(fhCellWarm) delete fhCellWarm;
275 
276 }
290 //________________________________________________________________________
291 void BadChannelAna::Run(Bool_t mergeOnly)
292 {
293  gStyle->SetPalette(55); //kRainBow==55
294  // cout<<"fired trigger class"<<AliAODEvent::GetFiredTriggerClasses()<<endl;
295  PrintCellInfo(1);
296  if(fExternalFileName=="")
297  {
298  //..If no extrenal file is provided merge different runs together
299  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
300  cout<<". . .Start process by converting files. . . . . . . . . . . ."<<endl;
301  cout<<endl;
303  if(fMergedFileName.IsNull())
304  {
305  Printf("File not produced, exit");
306  return;
307  }
308  cout<<endl;
309  }
310  else
311  {
312  //..If extrenal file is provided load it
313  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
314  cout<<". . .Start process by loading external file. . . . . . . . . . ."<<endl;
315  fMergedFileName= Form("%s/%s/%s/%s",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),fExternalFileName.Data());
316  }
317  PrintCellInfo(2);
318  //..if ==1 only produce filtered and merged files and do not perform a BC analysis
319  if(mergeOnly==0)
320  {
321  cout<<". . .Load inputfile with name: "<<fMergedFileName<<" . . . . . . . ."<<endl;
322  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
323  //.. Read all the needed input for the Bad/Dead Channel analysis
324  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
325  TFile *mergedFileInput = new TFile(fMergedFileName);
326  if(!mergedFileInput->IsOpen())
327  {
328  Printf("Error! Input file not found, abort");
329  return;
330  }
331  fCellAmplitude = (TH2F*) mergedFileInput->Get("hCellAmplitude");
332  fCellTime = (TH2F*) mergedFileInput->Get("hCellTime");
333  fProcessedEvents = (TH1F*) mergedFileInput->Get("hNEvents");
334 
335  //..if you have no external bad map provided that you want
336  //..to test then determine the bad maps here
337  if(fExternalBadMapName=="")
338  {
339  cout<<". . .Continue process by . . . . . . . . . . . ."<<endl;
340  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
341  //.. DEAD CELLS
342  //.. Flag dead cells with fFlag=1
343  //.. this excludes cells from analysis (will not appear in results)
344  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
345  cout<<"o o o Flag dead cells o o o"<<endl;
346  PrintCellInfo(3);
347  FlagAsDead();
348  if(fPrint==1)cout<<endl;
349  if(fPrint==1)cout<<endl;
350  PrintCellInfo(4);
351  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
352  //.. BAD CELLS
353  //.. Flag dead cells with fFlag=2 and bigger
354  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
355  cout<<"o o o Flag bad cells o o o"<<endl;
356  BCAnalysis();
357  PrintCellInfo(5);
358  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
359  //..In the end summarize results
360  //..in a .pdf and a .txt file
361  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
362  if(fPrint==1)cout<<"o o o Write .txt for each period analyis with bad cells o o o"<<endl;
364  PrintCellInfo(6);
365  }
366  //..if you have an external bad map provided load the flags and histograms
367  else
368  {
370  }
371  }
372  if(fSMMask.size()>0)RunMaskSM();
373 
374  if(fPrint==1)cout<<"o o o Create summary documents for the entire analysis o o o"<<endl;
376  PrintCellInfo(7);
377  //..can not save the array directly into the root file. Try with a TTree if needed. fRootFile->WriteObject(fFlag,"FlagArray");
378  fRootFile->Close();
379  cout<<endl;
380 
381  //..make a reccomendation about the used energy range to be investigated
382  //..and the binning
383  TH1D *hRefDistr = BuildMeanFromGood();
384  Double_t totalevents = fProcessedEvents->Integral();
385  //..Find bin where reference has value "totalevents" (means 1 hit/event), and the corresponding x-value
386  Int_t binHeightOne = hRefDistr->FindLastBinAbove(1.0/totalevents);
387  Double_t binCentreHeightOne = hRefDistr->GetBinCenter(binHeightOne);
388  cout<<". . .Recomendation:"<<endl;
389  cout<<". . .With the current statistic on average a cell has 1 hit at "<<binCentreHeightOne<<" GeV"<<endl;
390  cout<<". . .so it makes no sense to select energy ranges >"<<binCentreHeightOne<<"GeV as cells will be"<<endl;
391  cout<<". . .marked bad just due to the lack of statistic"<<endl;
392  cout<<". . .your selected lower bond is "<<fEndLowerBound<<" GeV"<<endl;
393  if(binCentreHeightOne>=fEndLowerBound) cout<<". . .This means you are OK!"<<endl;
394  if(binCentreHeightOne< fEndLowerBound) cout<<". . .#!#!#!#! CAREFUL THIS COULD CAUSE TROUBLE AND THROW OUT MORE CELLS THAN NECESSARY #!#!#!#! "<<endl;
395  cout<<". . .End of process . . . . . . . . . . . . . . . . . . . . ."<<endl;
396  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
397 }
404 //________________________________________________________________________
406 {
407  cout<<"o o o Start conversion process o o o"<<endl;
408  cout<<"o o o period: " << fPeriod << ", pass: " << fTrainNo << ", trigger: "<<fTrigger<< endl;
409 
410  //..Create histograms needed for adding all the files together
411  TH1F *hNEventsProcessedPerRun= new TH1F("hNEvents","Number of processed events in analyzed runs",1,0,1);
412  TH2F *hCellAmplitude;
413  TH2F *hCellTime;
414 
415  //..Open the text file with the run list numbers and run index
416  cout<<"o o o Open .txt file with run indices. Name = " << fRunList << endl;
417  FILE *pFile = fopen(fRunList.Data(), "r");
418  if(!pFile)
419  {
420  cout<<"couldn't open file!"<<endl;
421  return "";
422  }
423  Int_t nEntr;
424  Int_t nEntrTot=0;
425  Int_t q;
426  Int_t ncols;
427  Int_t nlines = 0 ;
428  Int_t runId[500] ;
429  while (1)
430  {
431  ncols = fscanf(pFile," %d ",&q);
432  if (ncols< 0) break;
433  runId[nlines]=q;
434  nlines++;
435  }
436  fclose(pFile);
437 
438 
439  //..Open the different .root files with help of the run numbers from the text file
440  const Int_t nRun = nlines ;
441  TString base;
442  TString infile;
443  TString singleRunFileName;
444 
445  cout<<"o o o Start merging process of " << nRun <<" files"<< endl;
446  Int_t usedFiles=0;
447  //..loop over the amount of run numbers found in the previous text file.
448  for(Int_t i = 0 ; i < nRun ; i++)
449  {
450  base = Form("%s/%s/%s/%d", fWorkdir.Data(), fAnalysisInput.Data(), fTrainNo.Data(), runId[i]);
451  //..This is a run2 case
452  infile = Form("%s.root",base.Data()) ;
453 
454  cout<<" o Open .root file with name: "<<infile<<endl;
455  TFile *f = TFile::Open(infile);
456 
457  //..Do some basic checks
458  if(!f)
459  {
460  cout<<"Couldn't open/find .root file: "<<infile<<endl;
461  continue;
462  }
463  TDirectoryFile *dir = (TDirectoryFile *)f->Get(fQADirect);
464  if(!dir)
465  {
466  cout<<"Couln't open directory file in .root file: "<<infile<<", directory: "<<fQADirect<<endl;
467  continue;
468  }
469  TList *outputList = (TList*)dir->Get(fQADirect);
470  if(!outputList)
471  {
472  cout << "Couln't get TList from directory file: "<<fQADirect<<endl;
473  continue;
474  }
475  usedFiles++;
476  TH2F *hAmpId;
477  TH2F *hTimeId;
478  TH1F *hNEvents;
479 
480  hAmpId =(TH2F*)outputList->FindObject("EMCAL_hAmpId");
481  if(!hAmpId)
482  {
483  Printf("hAmpId not found");
484  outputList->ls();
485  continue;
486  }
487  else
488  {
489  hAmpId->SetName("hCellAmplitude");
490  hAmpId->SetTitle("Cell Amplitude");
491  }
492 
493  hTimeId =(TH2F*)outputList->FindObject("EMCAL_hTimeId");
494  if(!hTimeId)
495  {
496  Printf("hTimeId not found");
497  outputList->ls();
498  continue;
499  }
500  else
501  {
502  hTimeId->SetName("hCellTime");
503  hTimeId->SetTitle("Cell Time");
504  }
505 
506  if(i==0)
507  {
508  //..copy the properties of the mother histogram for adding them all up
509  hCellAmplitude=(TH2F*)hAmpId->Clone("DummyName1");
510  hCellAmplitude->Reset();
511  hCellAmplitude->SetDirectory(0); //otherwise histogram will dissapear when file is closed
512  //..copy the properties of the mother histogram for adding them all up
513  hCellTime=(TH2F*)hTimeId->Clone("DummyName2");
514  hCellTime->Reset();
515  hCellTime->SetDirectory(0); //otherwise histogram will dissapear when file is closed
516  }
517 
518  hNEvents =(TH1F *)outputList->FindObject("hNEvents");
519  if(!hNEvents)
520  {
521  Printf("hNEvents not found");
522  outputList->ls();
523  continue;
524  }
525  nEntr = (Int_t)hNEvents->GetEntries();
526 
527  //..does that mean do not merge small files?
528  if (nEntr<100)
529  {
530  cout <<" o File to small to be merged. Only N entries " << nEntr << endl;
531  continue ;
532  }
533  cout <<" o File with N entries " << nEntr<<" will be merged"<< endl;
534  nEntrTot+=nEntr;
535  hCellAmplitude->Add(hAmpId);
536  hCellTime->Add(hTimeId);
537  hNEventsProcessedPerRun->Add(hNEvents);
538 
539  //..Create copies of the original root files just with the bad channel QA
540  //..So that a run by run bad channel analysis can be performed more easily
541  singleRunFileName= Form("%s/%s/%s/%d_%sFiltered.root",fWorkdir.Data(),fAnalysisInput.Data(),fTrainNo.Data(),runId[i],fTrigger.Data());
542  TFile *singleRunFile = TFile::Open(singleRunFileName,"recreate");
543  //..do not yet normalize by number of events
544  //..due to binning issues in later histograms
545  //..but rather do it at the very end of the analysis
546  hAmpId ->Write();
547  hTimeId->Write();
548  hNEvents->Write();
549  singleRunFile->Close();
550 
551  outputList->Delete();
552  dir->Delete();
553  f->Close();
554  delete f;
555  }
556  //..avoid creating empty files
557  if(usedFiles>0)
558  {
559  //.. Save the merged histograms
560  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
561  cout<<"o o o "<<nEntrTot<<" events were merged"<<endl;
562  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
563  //hNEventsProcessedPerRun->SetName("hNEvents");
564  //hNEventsProcessedPerRun->SetTitle("Number of processed events");
565  hNEventsProcessedPerRun->Write();
566  hCellAmplitude->SetName("hCellAmplitude");
567  hCellAmplitude->SetTitle("Cell Amplitude");
568  hCellAmplitude->Write();
569  hCellTime->SetName("hCellTime");
570  hCellTime->SetTitle("Cell Time");
571  hCellTime->Write();
572  BCF->Close();
573  cout<<"o o o End conversion process o o o"<<endl;
574  if(hNEventsProcessedPerRun) delete hNEventsProcessedPerRun;
575  return fMergedFileName;
576  }
577  else
578  {
579  if(hNEventsProcessedPerRun) delete hNEventsProcessedPerRun;
580  return "";
581  }
582 }
586 //_________________________________________________________________________
588 {
589  if(fExternalBadMapName=="")cout<<"Error - no external Bad Map provided"<<endl;
590  else cout<<"Load external map: "<<fExternalBadMapName<<endl;
591 
592  //..access the standart root output of a bad channel analysis to
593  //..get the necessary histogram
594  TString extRootFileName = Form("./AnalysisOutput/%s/%s/%s",fPeriod.Data(),fTrainNo.Data(),fExternalBadMapName.Data());
595  TFile* outputRoot = TFile::Open(extRootFileName.Data());
596 
597  TH1F* hFlags =(TH1F*)outputRoot->Get("fhCellFlag");
598 
599  //..set info from external source
600  //..Direction of cell ID
601  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
602  {
603  Double_t extFlag = hFlags->GetBinContent(cell+1);
604  //..Cell flagged as dead.
605  //..Flag only if it hasn't been flagged before
606  fFlag[cell] =extFlag;
607  if(extFlag>fCriterionCounter)fCriterionCounter=extFlag;
608  }
609 }
615 //_________________________________________________________________________
617 {
618  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
619  //.. BAD CELLS
620  //.. Flag bad cells with fFlag= 2,3,4,5.. etc
621  //.. this excludes cells from subsequent analysis
622  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
623  TArrayD periodArray;
624  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
625  {
626  periodArray=fAnalysisVector.at(i);
627  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
629  if(fPrint==1)cout<<endl;
630  if(fPrint==1)cout<<endl;
631  }
632  cout<<"o o o End of bad channel analysis o o o"<<endl;
633 }
634 //
635 // Mask an entire SM before doing the BC analysis
636 // This is useful when you get info from QA that there are problems with one SM
637 // and you want to clean up your bad channels beforehand
638 //
639 //________________________________________________________________________
641 {
642  cout<<"o o o Manually mask SM "<<iSM<<" o o o"<<endl;
643  fSMMask.push_back(iSM);
644 }
645 
646 
647 //
648 // Mask an entire SM before doing the BC analysis
649 // This is useful when you get info from QA that there are problems with one SM
650 // and you want to clean up your bad channels beforehand
651 //
652 //________________________________________________________________________
654 {
655  Int_t NoSMmask = fSMMask.size();
656  //..Loop over cell ID
657  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
658  {
659  //..check to which SM the cell belongs
660  for(Int_t iSM=0; iSM<NoSMmask; iSM++)
661  {
662  if(cell>=fStartCellSM[fSMMask.at(iSM)] && cell<fStartCellSM[fSMMask.at(iSM)+1])
663  {
664  fFlag[cell] =1;
665  }
666  }
667  }
668 }
677 //________________________________________________________________________
679 {
680  TArrayD periodArray;
681  periodArray.Set(4);
682  periodArray.AddAt((Double_t)criteria,0);
683  periodArray.AddAt(nsigma,1);
684  periodArray.AddAt(emin,2);
685  periodArray.AddAt(emax,3);
686  fAnalysisVector.push_back(periodArray);
687 }
703 //____________________________________________________________________
705 {
706  //.. criterion should be between 1-5
707  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;
708  if(fPrint==1)cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
709  if(fPrint==1)cout<<"o o o This is PeriodAnalysis No. "<<fCriterionCounter<<" in your list o o o"<<endl;
710  if(fPrint==1 && criterion != 3)cout<<"o o o Done in the energy range E "<<emin<<" - "<<emax<<endl;
711  if(fPrint==1 && criterion == 3)cout<<"o o o Done in the time range t "<<emin<<" - "<<emax<<endl;
712 
713  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
714  //.. ANALYSIS OF CELLS WITH ENTRIES
715  //.. Build average distributions and fit them
716  //.. Three tests for bad cells:
717  //.. 1) Average energy per hit
718  //.. 2) Average hit per event
719  //.. 3) Average max of cell time distribution
720  //.. 4) To be implemented: Scaled energy per hit distribution
721  //.. 5) Scaled hit distribution
722  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
723  TH1F* histogram;
724  if(fPrint==1)cout<<"o o o Analyze average cell distributions o o o"<<endl;
725  //..For case 1 or 2
726  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax);
727  if(criterion > 3) histogram = BuildHitAndEnergyMeanScaled(criterion, emin, emax);
728  if(criterion == 3) histogram = BuildTimeMean(criterion, emin, emax); //.. in case of crit 3 emin is tmin and emax is tmax
729 
730  if(criterion==1 || criterion == 4)
731  {
732  // if(emin>=1.8)FlagAsBad(criterion, histogram, nsigma, -1);//..do not apply a lower boundary
733  if(emin>=fEndLowerBound)FlagAsBad(criterion, histogram, nsigma, -1);//..do not apply a lower boundary
734  else FlagAsBad(criterion, histogram, nsigma, 200);//400
735  }
736  if(criterion==2 || criterion == 5)
737  {
738  // if(emin>=1.8)FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
739  if(emin>=fEndLowerBound)FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
740  else FlagAsBad(criterion, histogram, nsigma, 601);
741  }
742  if(criterion==3) FlagAsBad_Time(criterion, histogram, nsigma);
743 
744  // if(histogram) delete histogram;
745 }
746 
755 //_________________________________________________________________________
757 {
758  if(fPrint==1)cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
759  TH1F *histogram;
760  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);
761  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);
762  histogram->SetXTitle("Abs. Cell Id");
763  if(crit==1)histogram->SetYTitle("Energy per hit");
764  if(crit==2)histogram->SetYTitle("Number of hits in cell");
765  histogram->GetXaxis()->SetNdivisions(510);
766 
767  //..count Events in Energy Range
768  TH1F* pojection = (TH1F*)fCellAmplitude->ProjectionX("Intermediate");
769  fnEventsInRange = pojection->Integral(pojection->GetBinContent(emin),pojection->GetBinContent(emax));
770 
771  //..here the average hit per event and the average energy per hit is caluclated for each cell.
772  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
773  {
774  Double_t Esum = 0;
775  Double_t Nsum = 0;
776 
777  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
778  {
779  //To Do: I think one should also take the different bin width into account
780  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
781  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
782  //..investigate only cells that were not flagged as dead ore bad
783  if (E < emin || E > emax || fFlag[cell]!=0) continue;
784  Esum += E*N;
785  Nsum += N;
786  }
787  //..Set the values only for cells that are not yet marked as bad
788  if(fFlag[cell]==0)
789  {
790  if(crit==2) histogram->SetBinContent(cell+1, Nsum); //..number of hits (per event -> not anymore)
791  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/(Nsum)); //..average energy per hit
792  }
793  }
794  return histogram;
795 }
796 
797 
798 //-------------------------------------------------------------------------------------------
799 // This function builds the scaled hit distribution
800 // It normalizes the hits/cell to the mean value of the row and the column of the cell
801 // this is done in an iterative procedure (about 5 iterations are needed)
802 // with this procedure, the effects of the TRD and SM structures on the EMCal can be minimized
803 // The output is a histogram with the hits/cell as a function of cell ID.
804 
805 // ONLY IMPLEMENTED FOR THE HITS PER CELL method
806 
810 // ------------------------------------------------------------------------------------------
812 {
813  if(fPrint==1)cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
814  TH1F *histogram = NULL;
815  if(crit==4)histogram = new TH1F(Form("hCellEtoN_E%.2f-%.2f",emin,emax),Form("Energy per hit, %.2f < E < %.2f GeV",emin,emax), fNoOfCells,0,fNoOfCells);
816  if(crit==5)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);
817  histogram->SetXTitle("Abs. Cell Id");
818  if(crit==4)histogram->SetYTitle("Energy per hit");
819  if(crit==5)histogram->SetYTitle("Number of hits in cell");
820  histogram->GetXaxis()->SetNdivisions(510);
821 
822 
824  TH2F *hEnergyScaled = (TH2F*) fCellAmplitude->Clone("hEnergyScaled");
825  TH1F *hEnergyCol = new TH1F("hEnergyCol", "hEnergyCol", 100,0.,100.);
826  TH1F *hEnergyRow = new TH1F("hEnergyRow", "hEnergyRow", 250,0.,250.);
827 
828 
830  Int_t cellCol, cellRow, trash, col, row;
831  Double_t dCellEnergy;
832  Double_t dNumOfHits = 0.;
833  Double_t *arrHits = new Double_t[fNoOfCells];
834 
836  Double_t NoOfHits = 0.;
837  for(int iCell = 1; iCell <= fNoOfCells; iCell++){
838  NoOfHits = 0.;
839  if(fFlag[iCell -1] == 0){
840  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
841  NoOfHits += hEnergyScaled->GetBinContent(k, iCell);
842  }
843  }
844  arrHits[iCell - 1] = NoOfHits;
845  }
846 
847 
848  // ....................................................................................................
849  // Find cells that are really bad and exclude them from the calculation of the mean column and mean row
850  // ....................................................................................................
851 
852  // Declare Histos for BC finding for scaling
853  TH1D *hHitDistrib_forScaling = new TH1D("hHitDistrib_forScaling","hHitDistrib_forScaling",1000,0.,TMath::Median(fNoOfCells,arrHits)*4);
854  TF1 *fgausForScaling = new TF1("fgausForScaling","gaus", 0.,TMath::Median(fNoOfCells,arrHits)*4);
855  TVectorD fFlagForScaling;
856  fFlagForScaling.ResizeTo(fNoOfCells);
857  Double_t dhits = 0.;
858  for(int i = 1; i <= fNoOfCells; i++){
859  dhits = 0.;
860  dNumOfHits = 0.;
861  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
862  if(crit==5){dhits += hEnergyScaled->GetBinContent(k, i);}
863  if(crit==4){
864  dhits += hEnergyScaled->GetBinContent(k, i + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(k);
865  dNumOfHits += hEnergyScaled->GetBinContent(k, i + 1);
866  }
867  }
868  if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
869  if(dhits > 0.1){
870  hHitDistrib_forScaling->Fill(dhits);
871  }
872  }
873 
875  fgausForScaling->SetParameter(1, hHitDistrib_forScaling->GetBinCenter(hHitDistrib_forScaling-> GetMaximumBin()));
876  hHitDistrib_forScaling->Fit(fgausForScaling,"MQ0");
877  hHitDistrib_forScaling->Fit(fgausForScaling,"MQ0");
878 
879  for(int iCell = 1; iCell <= fNoOfCells; iCell++){
880  dhits = 0.;
881  dNumOfHits = 0.;
882  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
883  if(crit==5){dhits += hEnergyScaled->GetBinContent(k, iCell);}
884  if(crit==4){
885  dhits += hEnergyScaled->GetBinContent(k, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(k);
886  dNumOfHits += hEnergyScaled->GetBinContent(k, iCell + 1);
887  }
888  }
889  if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
890  if(dhits > fgausForScaling->GetParameter(1) + 5* fgausForScaling->GetParameter(2) || dhits < fgausForScaling->GetParameter(1) - 5* fgausForScaling->GetParameter(2)){
891  fFlagForScaling[iCell - 1] = 1;
892  }
893 
894  }
895  delete hHitDistrib_forScaling;
896  delete fgausForScaling;
897  delete [] arrHits;
898 
899  TF1 *fFitCol;
900  TF1 *fFitRow;
901  //...........................................
902  //start iterative process of scaling of cells
903  //...........................................
904  for(int iter = 0; iter < 5; iter++){
905  // array of vectors for calculating the mean hits per col/row
906  std::vector<Double_t> vecCol[100];
907  std::vector<Double_t> vecRow[250];
908 
909 
910  // TH2F *hmap = new TH2F("","",100,0.,100,250,0.,250.);
911  //
912  // // Build the Hitmap
913  // for(int iCell = 1; iCell <= fNoOfCells; iCell++){
914  // if(fFlag[iCell -1 ] == 0 && fFlagForScaling[iCell - 1] == 0){
915  // fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell - 1 ,0,cellCol, cellRow, trash, col, row);
916  // dhits = 0.;
917  // dNumOfHits = 0.;
918  // for(int EBin = hEnergyScaled->GetXaxis()->FindBin(emin); EBin <= hEnergyScaled->GetXaxis()->FindBin(emax); EBin++){
919  // if(crit==5){dhits += hEnergyScaled->GetBinContent(EBin, Cell);}
920  // if(crit==4){
921  // dhits += hEnergyScaled->GetBinContent(EBin, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(EBin);
922  // dNumOfHits += hEnergyScaled->GetBinContent(EBin, iCell + 1);
923  // }
924  // }
925  // if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
926  // hmap->SetBinContent(col + 1, row + 1, dhits );
927  // }
928  // }
929 
930  // Calculate the mean number of hits of each row and column
931  for(int iCell = 0; iCell < fNoOfCells; iCell++){
932  if(fFlag[iCell] == 0 && fFlagForScaling[iCell] == 0){
933  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
934  dCellEnergy = 0.;
935  dNumOfHits = 0.;
936  for (int EBin = hEnergyScaled->GetXaxis()->FindBin(emin); EBin < hEnergyScaled->GetXaxis()->FindBin(emax); EBin++) {
937  if(crit==5){dCellEnergy += hEnergyScaled->GetBinContent(EBin, iCell + 1);}
938  if(crit==4){
939  dCellEnergy += hEnergyScaled->GetBinContent(EBin, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(EBin);
940  dNumOfHits += hEnergyScaled->GetBinContent(EBin, iCell + 1);
941  }
942 
943  }
944  if(crit==4 && dNumOfHits > 0){dCellEnergy = dCellEnergy / dNumOfHits;}
945 
946  if( dCellEnergy > 0.){
947  hEnergyCol->Fill(col + 0.5, dCellEnergy);
948  hEnergyRow->Fill(row + 0.5, dCellEnergy);
949  vecCol[col].push_back(dCellEnergy);
950  vecRow[row].push_back(dCellEnergy);
951  }
952  }
953  }
954 
955  // Fill the histogram: mean hit per column
956  for(int iCol = 1; iCol <= hEnergyCol->GetNbinsX() ; iCol++){
957  if(vecCol[iCol -1].size() > 0.){
958  hEnergyCol->SetBinContent(hEnergyCol->FindBin(iCol - 0.5), hEnergyCol->GetBinContent(hEnergyCol->FindBin(iCol - 0.5))/vecCol[iCol-1].size() );
959  }
960  }
961 
962  // Fill the histogram: mean hit per row
963  for(int iRow = 1; iRow <= hEnergyRow->GetNbinsX() ; iRow++){
964  if(vecRow[iRow -1].size() > 0.){
965  hEnergyRow->SetBinContent(hEnergyRow->FindBin(iRow - 0.5), hEnergyRow->GetBinContent(hEnergyRow->FindBin(iRow - 0.5))/vecRow[iRow-1].size() );
966  }
967  }
968 
969  // Global fit to hits per row
970  fFitRow = new TF1("fFitRow","[0]",0.,250.);
971  fFitRow->SetParameter(0, hEnergyRow->GetBinContent(10));
972  Double_t MeanRow = hEnergyRow->GetBinContent(10);
973 
974  // Global fit to hits per column
975  fFitCol = new TF1("fFitCol","[0]",0.,100.);
976  fFitCol->SetParameter(0, hEnergyCol->GetBinContent(10));
977  Double_t MeanCol = hEnergyCol->GetBinContent(10);
978 
979 
980  //Scale each cell by the deviation of the mean of the column and the global mean
981  for(int iCell = 0; iCell < fNoOfCells; iCell++){
982  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
983  if (hEnergyCol->GetBinContent(col + 1) > 0.) {
984  for(int EBin = 1; EBin < hEnergyScaled->GetNbinsX(); EBin++) {
985  hEnergyScaled->SetBinContent(EBin, iCell + 1, hEnergyScaled->GetBinContent(EBin, iCell + 1) * (MeanCol / hEnergyCol->GetBinContent(col + 1)));
986  }
987  }
988  }
989 
990  //Scale each cell by the deviation of the mean of the row and the global mean
991  for(int iCell = 0; iCell < fNoOfCells; iCell++){
992  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
993  if (hEnergyRow->GetBinContent(row + 1) > 0.) {
994  for(int EBin = 1; EBin < hEnergyScaled->GetNbinsX(); EBin++) {
995  hEnergyScaled->SetBinContent(EBin, iCell + 1, hEnergyScaled->GetBinContent(EBin, iCell + 1) * (MeanRow / hEnergyRow->GetBinContent(row + 1)));
996  }
997  }
998  }
999  //....................
1000  // iteration ends here
1001  //....................
1002  }
1003 
1004  //............................................................................................
1005  //..here the average hit per event and the average energy per hit is caluclated for each cell.
1006  //............................................................................................
1007  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1008  {
1009  Double_t Esum = 0;
1010  Double_t Nsum = 0;
1011 
1012  for (Int_t j = 1; j <= hEnergyScaled->GetNbinsX(); j++)
1013  {
1014  //To Do: I think one should also take the different bin width into account
1015  Double_t E = hEnergyScaled->GetXaxis()->GetBinCenter(j);
1016  Double_t N = hEnergyScaled->GetBinContent(j, cell+1);
1017  //..investigate only cells that were not flagged as dead ore bad
1018  if (E < emin || E > emax || fFlag[cell]!=0) continue;
1019  Esum += E*N;
1020  Nsum += N;
1021  }
1022  //..Set the values only for cells that are not yet marked as bad
1023  if(fFlag[cell]==0)
1024  {
1025  if(crit==5) histogram->SetBinContent(cell+1, Nsum); //..number of hits (per event -> not anymore)
1026  if(Nsum > 0. && crit==4)histogram->SetBinContent(cell+1, Esum/(Nsum)); //..average energy per hit
1027  }
1028  }
1029  delete hEnergyScaled;
1030  delete hEnergyCol;
1031  delete hEnergyRow;
1032  delete fFitCol;
1033  delete fFitRow;
1034 
1035  return histogram;
1036 }
1037 
1038 
1039 
1044 //_________________________________________________________________________
1046 {
1047  if(fPrint==1)cout<<" o Calculate maximum of cell time distribution "<<endl;
1048  TString name;
1049  TH1F *histogram = NULL;
1050  histogram = new TH1F(Form("timeMax_t%.2f-%.2f",tmin,tmax),Form("Maximum of time distr., %.2f < t < %.2f ns",tmin,tmax), fNoOfCells,0,fNoOfCells);
1051  histogram->SetXTitle("Abs. Cell Id");
1052  histogram->SetYTitle(Form("N_{cells}^{%.0f<t<%.0f}/N_{cells}^{all}",tmin,tmax));
1053  histogram->GetXaxis()->SetNdivisions(510);
1054 
1055  //..Time information
1056  Double_t dHits = 0.;
1057  Double_t dIn = 0.;
1058  Double_t dAll = 0.;
1059  for(int iCell = 1; iCell <= fNoOfCells; iCell ++)
1060  {
1061  if(fFlag[iCell - 1] != 0) { continue; }
1062  dIn = 0.;
1063  dAll = 0.;
1064  for(int iTime = 1; iTime <= fCellTime->GetNbinsX(); iTime++)
1065  {
1066  dHits = fCellTime->GetBinContent(iTime, iCell);
1067  if(fCellTime->GetXaxis()->GetBinCenter(iTime) > tmin && fCellTime->GetXaxis()->GetBinCenter(iTime) < tmax)
1068  {
1069  dIn += dHits;
1070  }
1071  dAll += dHits;
1072  }
1073  if(dIn > 0.)
1074  {
1075  histogram->SetBinContent(iCell, dAll/dIn);
1076  }
1077  else
1078  {
1079  histogram->SetBinContent(iCell, 0.);
1080  }
1081  }
1082  return histogram;
1083 }
1084 
1089 //_________________________________________________________________________
1091 {
1092  Int_t sumOfExcl=0;
1093  Int_t manualMaskCounter=0;
1094  //..sort the elements in manual mask list
1095  std::sort (fManualMask.begin(), fManualMask.end());
1096 
1097  //..Direction of cell ID
1098  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1099  {
1100  Double_t nSum = 0;
1101  //..Direction of amplitude (Checks energies from 0-nBins GeV)
1102  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1103  {
1104  //..cellID+1 = histogram bin
1105  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
1106  nSum += N;
1107  }
1108  //..If the amplitude in one cell is basically 0
1109  //..mark the cell as excluded
1110  if(nSum == 0 && fFlag[cell]==0)
1111  {
1112  //..Cell flagged as dead.
1113  //..Flag only if it hasn't been flagged before
1114  fFlag[cell] =1;
1115  sumOfExcl++;
1116  }
1117  //..add here the manual masking
1118  if(manualMaskCounter<(Int_t)fManualMask.size() && fManualMask.at(manualMaskCounter)==cell)
1119  {
1120  fFlag[cell] = 2;
1121  manualMaskCounter++;
1122  }
1123  }
1124  if(fPrint==1)cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
1125  if(fPrint==1)cout<<" ("<<sumOfExcl<<")"<<endl;
1126 }
1127 
1142 //_________________________________________________________________________
1143 void BadChannelAna::FlagAsBad(Int_t crit, TH1F* inhisto, Double_t nsigma, Double_t dnbins)
1144 {
1145  gStyle->SetOptStat(0); //..Do not plot stat boxes
1146  gStyle->SetOptFit(0); //
1147  if(fPrint==1 && crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
1148  if(fPrint==1 && crit==2)cout<<" o Fit average hit per event distribution"<<endl;
1149  if(fPrint==1 && crit==3)cout<<" o Fit average hit maximum distribution"<<endl;
1150 
1151  Int_t cellColumn=0,cellRow=0;
1152  Int_t cellColumnAbs=0,cellRowAbs=0;
1153  Int_t trash;
1154 
1155  TString histoName=inhisto->GetName();
1156  Double_t goodmax = 0. ;
1157  Double_t goodmin = 0. ;
1158  //..set a user range so that the min and max is only found in a certain range
1159  inhisto->GetXaxis()->SetRangeUser(fStartCell,fNoOfCells);//..
1160  Double_t dminVal = inhisto->GetMinimum(0);
1161  Double_t dmaxVal = inhisto->GetMaximum();
1162  Double_t inputBins=dnbins;
1163  Double_t medianOfHisto = 0.;
1164  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1165  //. . .determine settings for the histograms (range and binning)
1166  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1167 
1168  //- - - - - - - - - - - -
1169  //..Determine good min and max range for histograms
1170  //..calculate and print the "median"
1171  //..
1172  Int_t numBins = inhisto->GetXaxis()->GetNbins();
1173 
1174  numBins-=fStartCell;
1175  Double_t *x = new Double_t[numBins];
1176  Double_t* y = new Double_t[numBins];
1177  for (int i = 0; i < numBins; i++)
1178  {
1179  x[i] = inhisto->GetBinContent(i+fStartCell);
1180  y[i] = 1; //each value with weight 1
1181  if(x[i]==0)y[i] = 0;
1182  }
1183  medianOfHisto = TMath::Median(numBins,x,y);
1184 
1185  //..if dmaxVal is too far away from medianOfHisto the histogram
1186  //..range will be too large -> reduce the range
1187  //cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1188  if(medianOfHisto*10<dmaxVal)
1189  {
1190  if(medianOfHisto*100<dmaxVal)
1191  {
1192  dmaxVal=medianOfHisto+0.02*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1193  }
1194  else
1195  {
1196  dmaxVal=medianOfHisto+0.2*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1197  }
1198  //cout<<"- - - median too far away from max range"<<endl;
1199  //cout<<"corrected max value: "<<dmaxVal<<endl;
1200  }
1201 
1202  //- - - - - - - - - - - -
1203  //..Determine number of bins for the histogram
1204  //..find a proper binning automatically (mostly done to avoid steplike structures, 0 entries in bins etc)
1205  //cout<<"max value: "<<dmaxVal<<", min value: "<<dminVal<<endl;
1206  if(crit==2 || crit == 5)//..hits in cell
1207  {
1208  dnbins=dmaxVal-dminVal;
1209  if(dnbins>100)
1210  {
1211  if(dnbins>2000)dnbins=0.01*(dmaxVal-dminVal); //..maximum 5000 bins. changed to 3000 .. lets see..
1212  if(dnbins>2000)dnbins=0.001*(dmaxVal-dminVal);//..maximum 5000 bins.
1213  if(dnbins<100) dnbins=0.02*(dmaxVal-dminVal); //..minimum 100 bins.
1214  }
1215  }
1216  if(crit==1 || crit == 4) //..energy/hit
1217  {
1218  //..Finding a good binning automatically is a bit messy. IF you know a better
1219  //..way please go ahead and implenent it.
1220  dnbins=(dmaxVal-dminVal)*150; //150 if you have problems with low statistic 500 normal stat
1221  if(dnbins>500)dnbins=(dmaxVal-dminVal)*50; //150 if you have problems with low statistic 500 normal stat
1222 
1223  if(dnbins<25)dnbins=dnbins*3; //..this is the min. binning for E/hit distr.
1224  if(dnbins<25)dnbins=dnbins*2; //..this is the min. binning for E/hit distr.
1225 
1226  if(inputBins==-1)
1227  {
1228  //dnbins=200;
1229  dnbins=(dmaxVal-dminVal)*40;
1230  if(dnbins>500)dnbins=(dmaxVal-dminVal)*15;
1231  if(dnbins<25)dnbins=dnbins*3; //..this is the min. binning for E/hit distr.
1232  if(dnbins<25)dnbins=dnbins*2; //..this is the min. binning for E/hit distr.
1233  }
1234  }
1235  if(x) delete []x;
1236  if(y) delete []y;
1237 
1238  if(crit==3)
1239  {
1240  //..obtain the binwidth for
1241  Int_t binWidth = fCellTime->GetXaxis()->GetBinWidth(1);
1242  dnbins = fCellTime->GetXaxis()->GetNbins();
1243  dminVal= fCellTime->GetXaxis()->GetBinCenter(1)-(binWidth*1.0)/2.0;
1244  dmaxVal= fCellTime->GetXaxis()->GetBinCenter(dnbins)+(binWidth*1.0)/2.0;
1245  cout<<"set the new histo with following settings- bins: "<<dnbins<<", min val = "<<dminVal<<", max val:"<<dmaxVal<<endl;
1246  }
1247 
1248  if(dmaxVal<dminVal || dnbins<1)
1249  {
1250  cout<<"###############################"<<endl;
1251  cout<<"Big problem: min:"<<dminVal<<", dmaxVal"<<dmaxVal<<endl;
1252  cout<<"dnbins:"<<dnbins<<endl;
1253  cout<<"###############################"<<endl;
1254  }
1255  //..build histos
1256  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1257  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1258  distrib->SetYTitle("Entries");
1259  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1260  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1261 
1262  //..build two dimensional histogram with values row vs. column
1263  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);
1264  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1265  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1266 
1267  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1268  //. . .build the distribution of average values
1269  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1270  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1271  {
1272  //..Do that only for cell ids also accepted by the code
1273  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1274 
1275  //..Fill histograms only for cells that are not yet flagged as bad
1276  if(fFlag[cell]==0)
1277  {
1278  //..fill the distribution of avarge cell values
1279  distrib->Fill(inhisto->GetBinContent(cell+1));
1280  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1281  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1282  //..Get Row and Collumn for cell ID
1283  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1284  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1285  {
1286  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1287  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1288  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1289  }
1290  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
1291  //..check TRD support structure
1292  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
1293  {
1294  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1295  }
1296  else
1297  {
1298  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1299  }
1300  }
1301  }
1302 
1303  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1304  //. . .draw histogram + distribution
1305  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1306  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
1307  c1->ToggleEventStatus();
1308  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
1309  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
1310  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
1311  upperPad->Draw();
1312  lowerPadLeft->Draw();
1313  lowerPadRight->Draw();
1314 
1315  upperPad->cd();
1316  upperPad->SetLeftMargin(0.05);
1317  upperPad->SetRightMargin(0.05);
1318  if(crit!=3)upperPad->SetLogy();
1319  inhisto->SetTitleOffset(0.6,"Y");
1320  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
1321  inhisto->GetYaxis()->SetTitleOffset(0.7);
1322  inhisto->SetLineColor(kBlue+1);
1323  inhisto->DrawCopy();
1324 
1325  lowerPadRight->cd();
1326  lowerPadRight->SetLeftMargin(0.09);
1327  lowerPadRight->SetRightMargin(0.12);
1328  plot2D->GetYaxis()->SetTitleOffset(1.3);
1329  plot2D->DrawCopy("colz");
1330 
1331  lowerPadLeft->cd();
1332  lowerPadLeft->SetLeftMargin(0.09);
1333  lowerPadLeft->SetRightMargin(0.07);
1334  lowerPadLeft->SetLogy();
1335  distrib->SetLineColor(kBlue+1);
1336  distrib->GetYaxis()->SetTitleOffset(1.3);
1337  //cout<<medianOfHisto<<endl;
1338  //Double_t drawRangeDown = 0.;//medianOfHisto - 0.5*medianOfHisto;
1339  //Double_t drawRangeUp = 2000.;//medianOfHisto + 0.5*medianOfHisto;
1340  //distrib->GetXaxis()->SetRangeUser(0., 2.);
1341  distrib->DrawCopy("");
1342  distrib_wTRDStruc->SetLineColor(kGreen+1);
1343  distrib_wTRDStruc->DrawCopy("same");
1344  distrib_woTRDStruc->SetLineColor(kMagenta+1);
1345  distrib_woTRDStruc->DrawCopy("same");
1346 
1347 
1348  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1349  //. . .fit histogram
1350  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1351 
1352  //..to exclude first 2 bins from max finding
1353  // distrib->GetXaxis()->SetRangeUser(distrib->GetBinCenter(2),distrib->GetBinCenter(dnbins));
1354  Double_t maxDistr = distrib->GetMaximum();
1355  Int_t maxBin = distrib->GetMaximumBin();
1356  Double_t maxCenter = distrib->GetBinCenter(maxBin);
1357 
1358  //..good range is around the max value as long as the
1359  //..bin content is larger than 2 entries
1360  for(Int_t i = maxBin ; i<=dnbins ; i++)
1361  {
1362  if(distrib->GetBinContent(i)<2) break ;
1363  goodmax = distrib->GetBinCenter(i);
1364  }
1365  for(Int_t i = maxBin ; i>2 ; i--)
1366  {
1367  if(distrib->GetBinContent(i)<2) break ;
1368  goodmin = distrib->GetBinLowEdge(i);
1369  }
1370 
1371  //if(maxBin<2)goodmin = distrib->GetBinLowEdge(2);
1372 
1373  //..Define min/max range of the fit function
1374  Double_t minFitRange=goodmin;
1375  Double_t maxFitRange=goodmax;
1376  //if(crit==2)minFitRange=distrib->GetBinLowEdge(2);//..exclude first 2 bins
1377  //if(crit==2)maxFitRange=dmaxVal;
1378  if(crit==3)minFitRange=-20;
1379  if(crit==3)maxFitRange=20;
1380 
1381  //cout<<"max bin: "<<maxBin<<", max value: "<<maxDistr<<endl;
1382  //cout<<"start mean: "<<maxCenter<<", mean range: 0-"<<dmaxVal<<endl;
1383  //cout<<"fit range: "<<minFitRange<<" - "<<maxFitRange<<endl;
1384  TF1 *fit2 = new TF1("fit2", "gaus",minFitRange,maxFitRange);
1385  //..start the fit with a mean of the highest value
1386  fit2->SetParLimits(0,0,maxDistr); //..do not allow negative ampl values
1387  fit2->SetParameter(1,maxCenter); //..set the start value for the mean
1388  fit2->SetParLimits(1,0,dmaxVal); //..do not allow negative mean values
1389 
1390  //..ELI - TO DO: eventually fit with TRD and without TRD seperatley
1391  distrib->Fit(fit2, "0LQEM", "", minFitRange, maxFitRange);
1392  Double_t sig, mean;// chi2ndf;
1393  mean = fit2->GetParameter(1);
1394  sig = fit2->GetParameter(2);
1395 
1396  //..for case 1 and 2 lower than 0 is an unphysical value
1397  // if(crit=!3 && mean <0.) mean=0.;
1398 
1399  goodmin = mean - nsigma*sig ;
1400  goodmax = mean + nsigma*sig ;
1401  //..for case 1 and 2 lower than 0 is an unphysical value
1402  // if(crit=!3 && goodmin <0.) goodmin=0.;
1403  //..this (below) is a special case for energy ranges where cell do not have many
1404  //..entries typically. One should not apply a lower bound in this case.
1405  if(inputBins==-1) goodmin=-1;
1406  if(goodmin <=0.) goodmin = -100;
1407  if(fPrint==1)cout<<" o Result of fit: "<<endl;
1408  if(fPrint==1)cout<<" o "<<endl;
1409  if(fPrint==1)cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
1410  if(fPrint==1)cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
1411 
1412  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1413  //. . .Add info to histogram
1414  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1415  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1416  lline->SetLineColor(kGreen+2);
1417  lline->SetLineStyle(7);
1418  lline->Draw();
1419 
1420  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1421  rline->SetLineColor(kGreen+2);
1422  rline->SetLineStyle(7);
1423  rline->Draw();
1424 
1425  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1426  leg->AddEntry(lline,"Good region boundary","l");
1427  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1428  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1429  leg->SetBorderSize(0);
1430  leg->Draw("same");
1431 
1432  fit2->SetLineColor(kOrange-3);
1433  fit2->SetLineStyle(1);//7
1434  fit2->Draw("same");
1435 
1436 
1437  TLatex* text = 0x0;
1438  if(crit==1 || crit==4) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1439  if(crit==2 || crit==5) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1440  if(crit==3) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1441  text->SetTextSize(0.06);
1442  text->SetNDC();
1443  text->SetTextColor(1);
1444  text->Draw();
1445 
1446  upperPad->cd();
1447  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
1448  uline->SetLineColor(kGreen+2);
1449  uline->SetLineStyle(7);
1450  uline->Draw();
1451  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
1452  lowline->SetLineColor(kGreen+2);
1453  lowline->SetLineStyle(7);
1454  lowline->Draw();
1455  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1456  //. . .Save histogram
1457  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1458  c1->Update();
1459  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
1460  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1461  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1462  if(crit==3)name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1463  if(crit==4)name=Form("%s/%s/AverageEperHit_scaled_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1464  if(crit==5)name=Form("%s/%s/AverageHitperEvent_scaled_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1465  c1->SaveAs(name);
1466 
1467  fRootFile->WriteObject(c1,c1->GetName());
1468  fRootFile->WriteObject(plot2D,plot2D->GetName());
1469  fRootFile->WriteObject(distrib,distrib->GetName());
1470  fRootFile->WriteObject(inhisto,inhisto->GetName());
1471  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1472  //. . . Mark the bad cells in the fFlag array
1473  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1474  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1475  //. . .(0 by default - good cell)
1476  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1477  if(fPrint==1)cout<<" o Flag bad cells that are outside the good range "<<endl;
1478  for(Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1479  {
1480  //..cell=0 and bin=1, cell=1 and bin=2
1481  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1482  if(inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)
1483  {
1484  fFlag[cell]=fCriterionCounter;
1485  }
1486  if(inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
1487  {
1488  fFlag[cell]=fCriterionCounter;
1489  }
1490  }
1491  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;
1492 
1493  /*
1494  delete distrib;
1495  delete plot2D;
1496  delete c1;
1497 
1498 
1499  delete fit2;
1500  delete lline;
1501  delete rline;
1502  delete leg;
1503  delete text;
1504  delete uline;
1505  delete lowline;
1506  */
1507  //delete upperPad;
1508  //delete lowerPadLeft;
1509  //delete lowerPadRight;
1510 }
1511 
1512 
1513 
1514 //.. This function flags bad cells according to their time distribution
1515 //.. The Cut is the upper limit for each cell and is set by hand in the runAnalysisBC macro
1516 //..
1520 
1521 void BadChannelAna::FlagAsBad_Time(Int_t crit, TH1F* inhisto, Double_t nSig)
1522 {
1523  Int_t cellColumn=0,cellRow=0;
1524  Int_t cellColumnAbs=0,cellRowAbs=0;
1525  Int_t trash;
1526 
1527  //..set a user range so that the min and max is only found in a certain range
1528  inhisto->GetXaxis()->SetRangeUser(fStartCell,fNoOfCells);//..
1529  Double_t dminVal = inhisto->GetMinimum(0);
1530  Double_t dmaxVal = inhisto->GetMaximum();
1531  //Double_t dnbins;
1532  Double_t medianOfHisto = 0.;
1533  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1534  //. . .determine settings for the histograms (range and binning)
1535  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1536 
1537  //..Determine good min and max range for histograms
1538  //..A. calculate and print the "median"
1539  Int_t numBins = inhisto->GetXaxis()->GetNbins();
1540 
1541  numBins-=fStartCell;
1542  Double_t *x = new Double_t[numBins];
1543  Double_t* y = new Double_t[numBins];
1544  for (int i = 0; i < numBins; i++)
1545  {
1546  x[i] = inhisto->GetBinContent(i+fStartCell);
1547  y[i] = 1; //each value with weight 1
1548  if(x[i]==0)y[i] = 0;
1549  }
1550  medianOfHisto = TMath::Median(numBins,x,y);
1551 
1552  //..B. if dmaxVal is too far away from medianOfHisto the histogram
1553  //..range will be too large -> reduce the range
1554  cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1555  if(medianOfHisto*10<dmaxVal)
1556  {
1557  if(medianOfHisto*100<dmaxVal)
1558  {
1559  dmaxVal=medianOfHisto+0.002*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1560  }
1561  else if(medianOfHisto*20<dmaxVal)
1562  {
1563  dmaxVal=medianOfHisto+0.02*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1564  }
1565  else
1566  {
1567  dmaxVal=medianOfHisto+0.05*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1568  }
1569  //cout<<"- - - median too far away from max range"<<endl;
1570  //cout<<"corrected max value: "<<dmaxVal<<endl;
1571  }
1572  cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1573  if(dmaxVal<dminVal)
1574  {
1575  cout<<"###############################"<<endl;
1576  cout<<"Big problem: min:"<<dminVal<<", dmaxVal"<<dmaxVal<<endl;
1577  cout<<"###############################"<<endl;
1578  }
1579 
1580  TString histoName=inhisto->GetName();
1581  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", 80, dminVal, dmaxVal);
1582  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1583  distrib->SetYTitle("Entries");
1584  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistrWtrd",(const char*)histoName), "", 80, dminVal, dmaxVal);
1585  TH1F *distrib_woTRDStruc = new TH1F(Form("%sDistrWoTRD",(const char*)histoName), "", 80, dminVal, dmaxVal);
1586  //..build two dimensional histogram with values row vs. column
1587  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);
1588  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1589  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1590 
1591  // cout<<"histoName: "<<histoName<<endl;
1592  // cout<<"distrib: "<<distrib->GetName()<<endl;
1593  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1594  //. . .build the distribution of average values
1595  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1596  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1597  {
1598  //..Do that only for cell ids also accepted by the code
1599  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1600 
1601  //..Fill histograms only for cells that are not yet flagged as bad
1602  if(fFlag[cell]==0)
1603  {
1604  //..fill the distribution of avarge cell values
1605  distrib->Fill(inhisto->GetBinContent(cell+1));
1606  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1607  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1608  //..Get Row and Collumn for cell ID
1609  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1610  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1611  {
1612  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1613  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1614  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1615  }
1616  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
1617  //..check TRD support structure
1618  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
1619  {
1620  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1621  }
1622  else
1623  {
1624  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1625  }
1626  }
1627  }
1628 
1629 
1630  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
1631  c1->ToggleEventStatus();
1632  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
1633  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
1634  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
1635  upperPad->Draw();
1636  lowerPadLeft->Draw();
1637  lowerPadRight->Draw();
1638 
1639  upperPad->cd();
1640  upperPad->SetLeftMargin(0.05);
1641  upperPad->SetRightMargin(0.05);
1642  upperPad->SetLogy();
1643  inhisto->SetTitleOffset(0.6,"Y");
1644  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
1645  inhisto->GetYaxis()->SetTitleOffset(0.7);
1646  inhisto->SetLineColor(kBlue+1);
1647  inhisto->Draw();
1648 
1649  //- - - - - - - - - - - - - - - - - - - -
1650  lowerPadRight->cd();
1651  lowerPadRight->SetLeftMargin(0.09);
1652  lowerPadRight->SetRightMargin(0.12);
1653  plot2D->GetYaxis()->SetTitleOffset(1.3);
1654  plot2D->Draw("colz");
1655 
1656  //- - - - - - - - - - - - - - - - - - - -
1657  lowerPadLeft->cd();
1658  lowerPadLeft->SetLeftMargin(0.09);
1659  lowerPadLeft->SetRightMargin(0.07);
1660  lowerPadLeft->SetLogy();
1661  distrib->SetLineColor(kBlue+1);
1662  distrib->GetYaxis()->SetTitleOffset(1.3);
1663  distrib->Draw("");
1664  distrib_wTRDStruc->SetLineColor(kGreen+1);
1665  distrib_wTRDStruc->DrawCopy("same");
1666  distrib_woTRDStruc->SetLineColor(kMagenta+1);
1667  distrib_woTRDStruc->DrawCopy("same");
1668 
1669  //- - - - - - - - - - - - - - - - - - - -
1670  //...fit time distr
1671  Double_t maxDistr = distrib->GetMaximum();
1672  Int_t maxBin = distrib->GetMaximumBin();
1673  Double_t maxCenter = distrib->GetBinCenter(maxBin);
1674  Double_t sig=0;
1675  Double_t mean=0;
1676 
1677  Double_t minFitRange=dminVal;
1678  Double_t maxFitRange=dmaxVal;
1679  //cout<<"max bin: "<<maxBin<<", max value: "<<maxDistr<<endl;
1680  //cout<<"start mean: "<<maxCenter<<", mean range: 0-"<<dmaxVal<<endl;
1681  //cout<<"fit range: "<<minFitRange<<" - "<<maxFitRange<<endl;
1682 
1683  TF1 *fitTime = new TF1("fitTime", "gaus",minFitRange,maxFitRange);
1684  //..start the fit with a mean of the highest value
1685  fitTime->SetParLimits(0,maxDistr*0.5,maxDistr*1.2); //..do not allow negative ampl values
1686  fitTime->SetParameter(1,maxCenter); //..set the start value for the mean
1687  fitTime->SetParLimits(1,dminVal,dmaxVal); //..do not allow negative mean values
1688 
1689  distrib->Fit(fitTime, "0LEM", "", minFitRange, maxFitRange);
1690  mean = fitTime->GetParameter(1);
1691  sig = fitTime->GetParameter(2);
1692 
1693  //cout<<"mean:"<<mean<<endl;
1694  //cout<<"sig: "<<sig<<endl;
1695 
1696  Double_t tCutMax=mean+nSig*sig;
1697  Double_t tCutMin=mean-nSig*sig;
1698 
1699  //tCutMax=1.9; //..Hard coded value for difficult cases
1700 
1701  cout<<"Determined time cut: "<<tCutMin<<" - "<<tCutMax<<endl;
1702  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1703  //. . .Add info to histogram
1704  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1705  TLine *line = new TLine(tCutMax, 0, tCutMax, distrib->GetMaximum());
1706  line->SetLineColor(kGreen+2);
1707  line->SetLineStyle(7);
1708  line->Draw();
1709  TLine *line2 = new TLine(tCutMin, 0, tCutMin, distrib->GetMaximum());
1710  line2->SetLineColor(kGreen+2);
1711  line2->SetLineStyle(7);
1712  line2->Draw();
1713 
1714  fitTime->SetLineColor(kOrange-3);
1715  fitTime->SetLineStyle(1);//7
1716  fitTime->Draw("same");
1717 
1718  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1719  leg->AddEntry(line,"upper limit","l");
1720  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1721  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1722  leg->SetBorderSize(0);
1723  leg->Draw("same");
1724 
1725  TLatex* text = 0x0;
1726  text = new TLatex(0.12,0.85,Form("Good range: %.3f-%.3f",tCutMin,tCutMax));
1727  text->SetTextSize(0.06);
1728  text->SetNDC();
1729  text->SetTextColor(1);
1730  text->Draw();
1731 
1732  //- - - - -
1733  upperPad->cd();
1734  TLine *uline = new TLine(0, tCutMax,fNoOfCells,tCutMax);
1735  uline->SetLineColor(kGreen+2);
1736  uline->SetLineStyle(7);
1737  uline->Draw();
1738 
1739  TLine *lline = new TLine(0, tCutMin,fNoOfCells,tCutMin);
1740  lline->SetLineColor(kGreen+2);
1741  lline->SetLineStyle(7);
1742  lline->Draw();
1743  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1744  //. . .Save histogram
1745  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1746  c1->Update();
1747  TString name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1748  c1->SaveAs(name);
1749 
1750  fRootFile->WriteObject(c1,c1->GetName());
1751  fRootFile->WriteObject(plot2D,plot2D->GetName());
1752  fRootFile->WriteObject(distrib,distrib->GetName());
1753  fRootFile->WriteObject(inhisto,inhisto->GetName());
1754  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1755  //. . . Mark the bad cells in the fFlag array
1756  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1757  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1758  //. . .(0 by default - good cell)
1759  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1760  if(fPrint==1)cout<<" o Flag bad cells that are outside the good range "<<endl;
1761  for(Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1762  {
1763  //..cell=0 and bin=1, cell=1 and bin=2
1764  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1765  if(inhisto->GetBinContent(cell+1) > tCutMax && fFlag[cell]==0)
1766  {
1767  fFlag[cell]=fCriterionCounter;
1768  }
1769  }
1770  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;
1771 
1772  /*
1773  delete distrib;
1774  delete distrib_wTRDStruc;
1775  delete distrib_woTRDStruc;
1776  delete plot2D;
1777  delete c1;
1778  delete line;
1779  delete leg;
1780  delete text;
1781  delete uline;
1782  */
1783 }
1784 
1785 
1790 //________________________________________________________________________
1792 {
1793  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1794  //.. RESULTS
1795  //.. 1) Print the bad cells
1796  //.. and write the results to a file
1797  //.. for each added period analysis
1798  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1799  TArrayD periodArray;
1800  Double_t emin,emax,sig;
1801  Int_t criterion;
1802  TString output;
1803  Int_t nb1=0;
1804 
1805  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1806  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1807  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1808  TString aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1809  ofstream file2(aliceTwikiTable, ios::out | ios::trunc);
1810  if(file2)
1811  {
1812  file2<<"|*Criterium* | *E_min !GeV* | *E_max !GeV* | *sigma* | *Excluded Cells* |"<<endl;
1813  }
1814  file2.close();
1815 
1816  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
1817  {
1818  periodArray=fAnalysisVector.at(i);
1819  criterion =periodArray.At(0);
1820  emin =periodArray.At(2);
1821  emax =periodArray.At(3);
1822  sig =periodArray.At(1);
1823 
1824  //..Print the results on the screen and
1825  //..write the results in a file
1826  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
1827  ofstream file(output, ios::out | ios::trunc);
1828  if(!file)
1829  {
1830  cout<<"#### Major Error. Check the textfile!"<<endl;
1831  }
1832  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1833  if(fPrint==1)cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1834 
1835  nb1=0;
1836  for(Int_t cellID=fStartCell ;cellID<fNoOfCells;cellID++)
1837  {
1838  if(fFlag[cellID]==(i+2))
1839  {
1840  nb1++;
1841  file<<cellID<<", ";
1842  }
1843  }
1844  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1845  file<<"("<<nb1<<")"<<endl;
1846  file.close();
1847  if(fPrint==1)cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1848  if(fPrint==1)cout<<endl;
1849  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1850  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1851  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1852  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1853  if(file2)
1854  {
1855  file2<<"| "<<criterion<<" | "<<emin<<" | "<<emax<<" | "<<sig<<" | "<<nb1<<" |"<<endl;
1856  }
1857  file2.close();
1858  }
1859 }
1866 //________________________________________________________________________
1868 {
1869  Int_t cellID, nDeadDCalCells = 0, nDeadEMCalCells = 0, nDCalCells = 0, nEMCalCells = 0;
1870  Double_t perDeadEMCal,perDeadDCal,perBadEMCal,perBadDCal,perWarmEMCal,perWarmDCal;
1871  TString aliceTwikiTable, cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells,goodCellsRatio,cellProp;
1872  TString OADBFile_bad, OADBFile_dead, OADBFile_warm;
1873  TString OADBFile_badRbR, OADBFile_deadRbR, OADBFile_warmRbR;
1874  TH2F* cellAmp_masked = (TH2F*)fCellAmplitude->Clone("hcellAmp_masked");
1875  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1876 
1877  deadPdfName = Form("%s/%s/%s_Dead_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1878  badPdfName = Form("%s/%s/%s_Bad_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1879  ratioOfBad = Form("%s/%s/%s_Bad_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1880  goodCells = Form("%s/%s/%s_Good_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1881  goodCellsRatio = Form("%s/%s/%s_Good_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1882  cellSummaryFile = Form("%s/%s/%s_%s_Bad_Ampl_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1883  OADBFile_bad = Form("%s/%s/%s_%s_OADBFile_Bad_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1884  OADBFile_dead = Form("%s/%s/%s_%s_OADBFile_Dead_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1885  OADBFile_warm = Form("%s/%s/%s_%s_OADBFile_Warm_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1886  if(fRunBRunMap)OADBFile_badRbR = Form("%s/%s/%s_%s_OADBFile_Bad_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1887  if(fRunBRunMap)OADBFile_deadRbR = Form("%s/%s/%s_%s_OADBFile_Dead_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1888  if(fRunBRunMap)OADBFile_warmRbR = Form("%s/%s/%s_%s_OADBFile_Warm_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1889  aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1890  cellProp = Form("%s/%s/%s_CellProp_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1891 
1892  cout<<" o Final results o "<<endl;
1893  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1894  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1895  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1896  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1897  {
1898  //..Direction of amplitude (Checks energies from 0-nBins GeV)
1899  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1900  {
1901  if(fFlag[cell]!=0)
1902  {
1903  //..cellID+1 = histogram bin
1904  cellAmp_masked->SetBinContent(amp,cell+1,0);
1905  }
1906  }
1907  //..Direction of time (Checks times from -275-975 GeV)
1908  for (Int_t time = 1; time <= fCellTime->GetNbinsX(); time++)
1909  {
1910  if(fFlag[cell]!=0)
1911  {
1912  //..cellID+1 = histogram bin
1913  cellTime_masked->SetBinContent(time,cell+1,0);
1914  }
1915  }
1916  }
1917  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1918  //..Plot some summary canvases
1919  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1920  TCanvas *c1 = new TCanvas("CellProp","I summary of cell properties",1000,1000);
1921  c1->ToggleEventStatus();
1922  c1->Divide(2,2);
1923  c1->cd(1)->SetLogz();
1924  //lowerPadRight->SetLeftMargin(0.09);
1925  //lowerPadRight->SetRightMargin(0.06);
1926  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1927  fCellAmplitude->SetYTitle("Abs. Cell Id");
1928  fCellAmplitude->GetYaxis()->SetTitleOffset(1.6);
1929  fCellAmplitude->Draw("colz");
1930  c1->cd(2)->SetLogz();
1931  fCellTime->SetXTitle("Cell Time [ns]");
1932  fCellTime->SetYTitle("Abs. Cell Id");
1933  fCellTime->GetYaxis()->SetTitleOffset(1.6);
1934  fCellTime->Draw("colz");
1935  c1->cd(3)->SetLogz();
1936  //lowerPadRight->SetLeftMargin(0.09);
1937  //lowerPadRight->SetRightMargin(0.06);
1938  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1939  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1940  cellAmp_masked->SetYTitle("Abs. Cell Id");
1941  cellAmp_masked->GetYaxis()->SetTitleOffset(1.6);
1942  cellAmp_masked->Draw("colz");
1943  c1->cd(4)->SetLogz();
1944  cellTime_masked->SetTitle("Masked Cell Time");
1945  cellTime_masked->SetXTitle("Cell Time [ns]");
1946  cellTime_masked->SetYTitle("Abs. Cell Id");
1947  cellTime_masked->GetYaxis()->SetTitleOffset(1.6);
1948  cellTime_masked->Draw("colz");
1949  c1->Update();
1950 
1951  TCanvas *c1_ratio = new TCanvas("CellPropRatio","II summary of cell properties ratio",1000,500);
1952  c1_ratio->ToggleEventStatus();
1953  c1_ratio->Divide(2);
1954  c1_ratio->cd(1)->SetLogz();
1955  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1956  cellAmp_masked->GetZaxis()->SetRangeUser(0.0001,10e7);
1957  cellAmp_masked->Draw("colz");
1958  c1_ratio->cd(2)->SetLogz();
1959 
1960  TH1D *hRefDistr = BuildMeanFromGood();
1961  TH2F* ratio2DAmp =(TH2F*)cellAmp_masked->Clone("ratio2DAmp");
1962  TH2F* Sum2DIdeal =(TH2F*)cellAmp_masked->Clone("Sum2DIdeal");
1963  Sum2DIdeal->Reset();
1964  //..Create an ideal 2D energy distribution for division.
1965  //..Helps to identify whether there are some cells that still
1966  //..need to be masked by hand
1967  for(Int_t eBin=0;eBin<Sum2DIdeal->GetNbinsX();eBin++)
1968  {
1969  Double_t binVal=hRefDistr->GetBinContent(eBin+1);
1970  for(Int_t icell=0;icell<Sum2DIdeal->GetNbinsY();icell++)
1971  {
1972  Sum2DIdeal->SetBinContent(eBin+1,icell+1,binVal);
1973  }
1974  }
1975  ratio2DAmp->SetTitle("Ratio of cell Amplitude to mean cell ampl.");
1976  ratio2DAmp->Divide(Sum2DIdeal);
1977  ratio2DAmp->GetZaxis()->UnZoom();
1978  ratio2DAmp->DrawCopy("colz");
1979  c1_ratio->Update();
1980 
1981  TLatex* textSM = new TLatex(0.1,0.1,"*test*");
1982  textSM->SetTextSize(0.06);
1983  textSM->SetTextColor(1);
1984  textSM->SetNDC();
1985 
1986  TCanvas *c1_proj = new TCanvas("CellPropPProj","III summary of cell properties",1000,500);
1987  c1_proj->ToggleEventStatus();
1988  c1_proj->Divide(2);
1989  c1_proj->cd(1)->SetLogy();
1990  TH1D* projEnergyMask = cellAmp_masked->ProjectionX(Form("%sMask_Proj",cellAmp_masked->GetName()),fStartCell,fNoOfCells);
1991  projEnergyMask->SetXTitle("Cell Energy [GeV]");
1992  projEnergyMask->GetYaxis()->SetTitleOffset(1.6);
1993  projEnergyMask->SetLineColor(kGreen+1);
1994  projEnergyMask->DrawCopy(" hist");
1995  TH1D* projEnergy = fCellAmplitude->ProjectionX(Form("%s_Proj",fCellAmplitude->GetName()),fStartCell,fNoOfCells);
1996  projEnergy->DrawCopy("same hist");
1997  TLegend *leg = new TLegend(0.50,0.75,0.7,0.87);
1998  leg->AddEntry(projEnergy,"all cells","l");
1999  leg->AddEntry(projEnergyMask,"good cells","l");
2000  leg->SetTextSize(0.05);
2001  leg->SetBorderSize(0);
2002  leg->SetFillColorAlpha(10, 0);
2003  leg->Draw("same");
2004  TLegend *legBig = (TLegend*)leg->Clone("legBig");
2005  legBig->SetTextSize(0.08);
2006  legBig->SetX1NDC(0.2);
2007 
2008  c1_proj->cd(2)->SetLogy();
2009  TH1* projTimeMask = cellTime_masked->ProjectionX(Form("%s_Proj",cellTime_masked->GetName()),fStartCell,fNoOfCells);
2010  projTimeMask->SetXTitle("Cell Time [ns]");
2011  projTimeMask->GetYaxis()->SetTitleOffset(1.6);
2012  projTimeMask->GetYaxis()->SetRangeUser(1,projTimeMask->GetMaximum()*20);
2013  projTimeMask->SetLineColor(kGreen+1);
2014  projTimeMask->DrawCopy("hist");
2015  TH1* projTime = fCellTime->ProjectionX(Form("%s_Proj",fCellTime->GetName()),fStartCell,fNoOfCells);
2016  projTime->DrawCopy("same hist");
2017  leg->Draw("same");
2018  c1_proj->Update();
2019 
2020  TCanvas *c1_projSM = new TCanvas("CellPropPProjSM","III summary of cell Energy per SM",1200,900);
2021  c1_projSM->Divide(5,4,0.001,0.001);
2022  TH1* projEnergyMaskSM[20];
2023  TH1* projEnergySM[20];
2024  for(Int_t iSM=0;iSM<20;iSM++)
2025  {
2026  c1_projSM->cd(iSM+1)->SetLogy();
2027  gPad->SetTopMargin(0.03);
2028  gPad->SetBottomMargin(0.11);
2029  projEnergyMaskSM[iSM] = cellAmp_masked->ProjectionX(Form("%sMask_ProjSM%i",cellAmp_masked->GetName(),iSM),fStartCellSM[iSM]+1,fStartCellSM[iSM+1]); //histogram bin 1 has cell ID0
2030  projEnergyMaskSM[iSM]->SetTitle("");
2031  projEnergyMaskSM[iSM]->SetXTitle(Form("Cell Energy [GeV], SM%i",iSM));
2032  projEnergyMaskSM[iSM]->GetYaxis()->SetTitleOffset(1.6);
2033  projEnergyMaskSM[iSM]->GetYaxis()->SetLabelSize(0.06);
2034  projEnergyMaskSM[iSM]->GetXaxis()->SetLabelSize(0.06);
2035  projEnergyMaskSM[iSM]->GetXaxis()->SetRangeUser(0,20);
2036  projEnergyMaskSM[iSM]->GetXaxis()->SetTitleSize(0.06);
2037  projEnergyMaskSM[iSM]->SetLineColor(kGreen+1);
2038  projEnergyMaskSM[iSM]->DrawCopy(" hist");
2039 
2040  projEnergySM[iSM] = fCellAmplitude->ProjectionX(Form("%s_ProjSM%i",fCellAmplitude->GetName(),iSM),fStartCellSM[iSM],fStartCellSM[iSM+1]-1);
2041  projEnergySM[iSM]->DrawCopy("same hist");
2042  if(iSM==0)legBig->Draw("same");
2043  //textSM->Draw();
2044  textSM->SetTitle(Form("Includes cell IDs %d-%d",fStartCellSM[iSM],fStartCellSM[iSM+1]-1));
2045  textSM->DrawLatex(0.2,0.9,Form("Includes cell IDs %d-%d",fStartCellSM[iSM],fStartCellSM[iSM+1]-1));
2046  }
2047  c1_projSM->Update();
2048 
2049  TCanvas *c1_projRSM = new TCanvas("CellPropPProjRSM","III summary of cell Energy Ratio per SM",1200,900);
2050  c1_projRSM->Divide(5,4,0.001,0.001);
2051  for(Int_t iSM=0;iSM<20;iSM++)
2052  {
2053  c1_projRSM->cd(iSM+1)->SetLogy();
2054  gPad->SetTopMargin(0.03);
2055  gPad->SetBottomMargin(0.11);
2056  //projEnergyMaskSM[iSM]->GetXaxis()->SetRangeUser(0,10);
2057  projEnergyMaskSM[iSM]->SetLineColor(kGray+1);
2058  projEnergyMaskSM[iSM]->Divide(hRefDistr);
2059  projEnergyMaskSM[iSM]->DrawCopy("hist");
2060  }
2061  c1_projRSM->Update();
2062 
2063  TCanvas *c1_projTimeSM = new TCanvas("CellPropPProjTimeSM","III summary of cell Time per SM",1200,900);
2064  c1_projTimeSM->Divide(5,4,0.001,0.001);
2065  TH1* projTimeMaskSM[20];
2066  TH1* projTimeSM[20];
2067  for(Int_t iSM=0;iSM<20;iSM++)
2068  {
2069  c1_projTimeSM->cd(iSM+1)->SetLogy();
2070  gPad->SetTopMargin(0.03);
2071  gPad->SetBottomMargin(0.11);
2072  projTimeMaskSM[iSM] = cellTime_masked->ProjectionX(Form("%sMask_ProjSMTime%i",cellAmp_masked->GetName(),iSM),fStartCellSM[iSM]+1,fStartCellSM[iSM+1]);
2073  projTimeMaskSM[iSM]->SetTitle("");
2074  projTimeMaskSM[iSM]->SetXTitle(Form("Cell Time [ns], SM%i",iSM));
2075  projTimeMaskSM[iSM]->GetYaxis()->SetTitleOffset(1.6);
2076  projTimeMaskSM[iSM]->GetYaxis()->SetLabelSize(0.06);
2077  projTimeMaskSM[iSM]->GetXaxis()->SetLabelSize(0.06);
2078  //projTimeMaskSM[iSM]->GetXaxis()->SetRangeUser(0,20);
2079  projTimeMaskSM[iSM]->GetXaxis()->SetTitleSize(0.06);
2080  projTimeMaskSM[iSM]->SetLineColor(kGreen+1);
2081  projTimeMaskSM[iSM]->DrawCopy(" hist");
2082 
2083  if(iSM==0)legBig->Draw("same");
2084  projTimeSM[iSM] = fCellTime->ProjectionX(Form("%s_ProjSMTime%i",fCellAmplitude->GetName(),iSM),fStartCellSM[iSM],fStartCellSM[iSM+1]-1);
2085  projTimeSM[iSM]->DrawCopy("same hist");
2086  }
2087 
2088  /*
2089  //..This part here eats up too much memory so it was commented out
2090  //..The Canvases are anyway later saved as individual gifs.
2091  //..save to a PDF
2092  c1 ->Print(Form("%s(",cellProp.Data()));
2093  c1_ratio ->Print(Form("%s",cellProp.Data()));
2094  c1_proj ->Print(Form("%s",cellProp.Data()));
2095  c1_projSM ->Print(Form("%s",cellProp.Data()));
2096  c1_projRSM ->Print(Form("%s",cellProp.Data()));
2097  c1_projTimeSM->Print(Form("%s)",cellProp.Data()));
2098  */
2099  //..Scale the histograms by the number of events
2100  //..so that they are more comparable for a run-by-run
2101  //..analysis
2102  Double_t totalevents = fProcessedEvents->Integral();
2103  fCellAmplitude ->Scale(1.0/totalevents);
2104  cellAmp_masked ->Scale(1.0/totalevents);
2105  fCellTime ->Scale(1.0/totalevents);
2106 
2107  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2108  //..Write the final results of dead and bad cells in a file and on screen
2109  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2110  ofstream file(cellSummaryFile, ios::out | ios::trunc);
2111  ofstream fileBad(OADBFile_bad, ios::out | ios::trunc);
2112  ofstream fileDead(OADBFile_dead, ios::out | ios::trunc);
2113  ofstream fileWarm(OADBFile_warm, ios::out | ios::trunc);
2114  /*ofstream fileBadRbR;
2115  ofstream fileDeadRbR;
2116  ofstream fileWarmRbR;
2117  if(fRunBRunMap)
2118  {*/
2119  ofstream fileBadRbR(OADBFile_badRbR, ios::out | ios::trunc);
2120  ofstream fileDeadRbR(OADBFile_deadRbR, ios::out | ios::trunc);
2121  ofstream fileWarmRbR(OADBFile_warmRbR, ios::out | ios::trunc);
2122  //}
2123  if(file)
2124  {
2125  file<<"Dead cells : "<<endl;
2126  cout<<" o Dead cells : "<<endl;
2127  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
2128  {
2129  if(cellID==0)
2130  {
2131  file<<"In EMCal: "<<endl;
2132  }
2133  if(cellID==fCellStartDCal)
2134  {
2135  file<<"\n"<<endl;
2136  file<<"In DCal: "<<endl;
2137  }
2138  if(fFlag[cellID]==1)
2139  {
2140  file<<cellID<<", ";
2141  if(cellID<fCellStartDCal)nDeadEMCalCells++;
2142  else nDeadDCalCells++;
2143  }
2144  if(fFlag[cellID]==1)fileDead<<cellID<<endl;
2145  if(fFlag[cellID]>1)fileBad<<cellID<<endl;
2146  if(fRunBRunMap==1 && fFlag[cellID]==1)fileDeadRbR<<cellID<<endl;
2147  if(fRunBRunMap==1 && fFlag[cellID]>1)fileBadRbR<<cellID<<endl;
2148  }
2149  file<<"\n"<<endl;
2150  perDeadEMCal=100*nDeadEMCalCells/(1.0*fCellStartDCal);
2151  perDeadDCal=100*nDeadDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2152  file<<"EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
2153  cout<<" o EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
2154 
2155  file<<"Bad cells: "<<endl;
2156  cout<<" o Bad cells: "<<endl;
2157  for(cellID=fStartCell;cellID<fNoOfCells;cellID++)
2158  {
2159  if(cellID==0)
2160  {
2161  file<<"In EMCal: "<<endl;
2162  }
2163  if(cellID==fCellStartDCal)
2164  {
2165  file<<"\n"<<endl;
2166  file<<"In DCal: "<<endl;
2167  }
2168  if(fFlag[cellID]>1)
2169  {
2170  file<<cellID<<", ";
2171  if(cellID<fCellStartDCal)nEMCalCells++;
2172  else nDCalCells++;
2173  }
2174  }
2175  file<<"\n"<<endl;
2176  perBadEMCal=100*nEMCalCells/(1.0*fCellStartDCal);
2177  perBadDCal =100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2178  file<<"EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
2179  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
2180  }
2181  file.close();
2182  cout<<" o Total: "<<endl;
2183  cout<<" o Bad+Dead cells: "<<nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells<<", this is "<<(nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells)*100/(1.0*fNoOfCells)<<"% of the whole detector"<<endl;
2184 
2185  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2186  //..Determine the number of warm cells and save the flagged cells to .pdf files
2187  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2188  if(fPrint==1)cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
2189  SaveBadCellsToPDF(1,badPdfName) ;
2190  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
2191  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
2192  if(fTestRoutine==1)SaveBadCellsToPDF(20,goodCellsRatio) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
2193 
2194  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2195  //..Fill the histograms with the flag information
2196  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2197  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2198  {
2199  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
2200  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
2201  if(fWarmCell[cell]==1)fileWarm<<cell<<endl;
2202  if(fRunBRunMap==1 && fWarmCell[cell]==1)fileWarmRbR<<cell<<endl;
2203  }
2204  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
2205  c2->ToggleEventStatus();
2206  c2->Divide(1,2);
2207  c2->cd(1);
2208  fhCellFlag->SetTitle("cell flag");
2209  fhCellFlag->SetXTitle("Abs. Cell Id");
2210  fhCellFlag->SetYTitle("flag by which cell was excluded");
2211  fhCellFlag->DrawCopy("hist");
2212  c2->cd(2);
2213  fhCellWarm->SetTitle("Warm cells");
2214  fhCellWarm->SetXTitle("Abs. Cell Id");
2215  fhCellWarm->SetYTitle("warm=1");
2216  fhCellWarm->DrawCopy("hist");
2217  c2->Update();
2218 
2219  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2220  //..Plot the 2D distribution of cells by flag
2221  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2222  PlotFlaggedCells2D(0); //..all good cells
2223  PlotFlaggedCells2D(1); //..all dead cells
2224  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
2225  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
2226 
2227  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2228  //..Add different histograms/canvases to the output root file
2229  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2230  TString name1,name2,name3,name4,name5,name6;
2231  name1 = Form("%s/%s/CellPropertiesRatio.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2232  c1_ratio->SaveAs(name1);
2233  name2 = Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2234  c1->SaveAs(name2);
2235  name3 = Form("%s/%s/CellPropertiesProj.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2236  c1_proj->SaveAs(name3);
2237  name4 = Form("%s/%s/CellEnergySM.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2238  c1_projSM->SaveAs(name4);
2239  name5 = Form("%s/%s/CellEnergySMratio.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2240  c1_projRSM->SaveAs(name5);
2241  name6 = Form("%s/%s/CellTimeSM.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2242  c1_projTimeSM->SaveAs(name6);
2243 
2244  fRootFile->WriteObject(c1_ratio,c1_ratio->GetName());
2245  fRootFile->WriteObject(c1,c1->GetName());
2246  fRootFile->WriteObject(c1_proj,c1_proj->GetName());
2247  fRootFile->WriteObject(c1_projSM,c1_projSM->GetName());
2248  fRootFile->WriteObject(c1_projRSM,c1_projRSM->GetName());
2249  fRootFile->WriteObject(c1_projTimeSM,c1_projTimeSM->GetName());
2250 
2251  fRootFile->WriteObject(c2,c2->GetName());
2252  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
2253  fRootFile->WriteObject(cellAmp_masked,cellAmp_masked->GetName());
2254  fRootFile->WriteObject(ratio2DAmp,ratio2DAmp->GetName());
2255  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
2256  fRootFile->WriteObject(fProcessedEvents,fProcessedEvents->GetName());
2257  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
2258  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
2259  fRootFile->WriteObject(projEnergyMask,projEnergyMask->GetName());
2260  fRootFile->WriteObject(projEnergy,projEnergy->GetName());
2261  //..Save all amplitudes to the root file
2262  SaveHistoToFile();
2263 
2264  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2265  //..Save also the identified warm channels into a text file.
2266  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2267  Int_t nWEMCalCells =0;
2268  Int_t nWDCalCells =0;
2269  file.open(cellSummaryFile, ios::out | ios::app);
2270  if(file)
2271  {
2272  file<<"Warm cells : "<<endl;
2273  if(fPrint==1)cout<<" o Warm cells : "<<endl;
2274  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
2275  {
2276  if(cellID==0)
2277  {
2278  file<<"In EMCal: "<<endl;
2279  }
2280  if(cellID==fCellStartDCal)
2281  {
2282  file<<"\n"<<endl;
2283  file<<"In DCal: "<<endl;
2284  }
2285  if(fWarmCell[cellID]==1)
2286  {
2287  file<<cellID<<", ";
2288  if(cellID<fCellStartDCal)nWEMCalCells++;
2289  else nWDCalCells++;
2290  }
2291  }
2292  file<<"\n"<<endl;
2293  perWarmEMCal= 100*nWEMCalCells/(1.0*fCellStartDCal);
2294  perWarmDCal = 100*nWDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2295  file<<"EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
2296  if(fPrint==1)cout<<" o EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
2297  }
2298  file.close();
2299  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2300  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
2301  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2302  ofstream file2(aliceTwikiTable, ios::out | ios::app);
2303  if(file2)
2304  {
2305  file2<<"1=energy/hit, 2= hit/event, 3=cell time"<<endl;
2306  file2<<"\n"<<endl;
2307  file2<<"| *Detector* | *No of cells* | *percentage* |"<<endl;
2308  file2<<"| Dead EMCal | "<<nDeadEMCalCells<<" | "<<perDeadEMCal<<"% |"<<endl;
2309  file2<<"| Bad EMCal | "<<nEMCalCells<<" | "<<perBadEMCal<<"% |"<<endl;
2310  file2<<"| - Warm EMCal | "<<nWEMCalCells<<" | "<<perWarmEMCal<<"% |"<<endl;
2311  file2<<"| Dead DCal | "<<nDeadDCalCells<<" | "<<perDeadDCal<<"% |"<<endl;
2312  file2<<"| Bad DCal | "<<nDCalCells<<" | "<<perBadDCal<<"% |"<<endl;
2313  file2<<"| - Warm DCal | "<<nWDCalCells<<" | "<<perWarmDCal<<"% |"<<endl;
2314  file2<<"| Summ D+B | "<<nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells<<" | "<<(nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells)*100/(1.0*fNoOfCells)<<"% |"<<endl;
2315  file2<<"\n"<<endl;
2316  }
2317  file2.close();
2318 
2319  /*
2320  if(c1) delete c1;
2321  if(c1_ratio) delete c1_ratio;
2322  if(c1_projSM) delete c1_projSM;
2323  if(c1_projRSM) delete c1_projRSM;
2324  if(c1_projTimeSM) delete c1_projTimeSM;
2325  if(c2) delete c2;
2326  if(textSM) delete textSM;
2327  if(c1_proj) delete c1_proj;
2328  if(leg) delete leg;
2329  if(legBig) delete legBig;
2330  */
2331 }
2332 
2348 //________________________________________________________________________
2350 {
2351  gROOT->SetStyle("Plain");
2352  gStyle->SetOptStat(0);
2353  gStyle->SetFillColor(kWhite);
2354  gStyle->SetTitleFillColor(kWhite);
2355  gStyle->SetPalette(91);
2356 
2357  char title[100];
2358  char name[100];
2359 
2360  TH1D *hRefDistr = BuildMeanFromGood();
2361  fRootFile->WriteObject(hRefDistr,hRefDistr->GetName());
2362  Int_t firstCanvas=0;
2363  Bool_t candidate;
2364  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
2365  text->SetTextSize(0.07);
2366  text->SetTextColor(kOrange-3);
2367  text->SetNDC();
2368 
2369  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
2370  text2->SetTextSize(0.08);
2371  text2->SetTextColor(8);
2372  text2->SetNDC();
2373 
2374  TLatex* textA = new TLatex(0.65,0.62,"*test*");
2375  textA->SetTextSize(0.04);
2376  textA->SetTextColor(1);
2377  textA->SetNDC();
2378 
2379  //..collect cells in an internal vector.
2380  //..when the vector is of size=9 or at the end of being filled
2381  //..plot the channels into a canvas
2382  std::vector<Int_t> channelVector;
2383  channelVector.clear();
2384  cout<<" o Start printing into .pdf for version: "<<version<<endl;
2385  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
2386  {
2387  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
2388  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
2389  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
2390  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
2391  if(fFlag[cell]==0 && version==20)channelVector.push_back(cell);
2392 
2393  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
2394  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
2395  if(channelVector.size()==9 || cell == fNoOfCells-1)
2396  {
2397  cout<<"."<<flush;
2398 
2399  TString internal_pdfName=pdfName;
2400  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
2401  if(channelVector.size() > 6) c1->Divide(3,3);
2402  else if (channelVector.size() > 3) c1->Divide(3,2);
2403  else c1->Divide(3,1);
2404 
2405  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
2406  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
2407  {
2408  if(channelVector.size() >=fNoOfCells) cout<<"Massive problem"<<endl;
2409  sprintf(name, "Cell %d",channelVector.at(i)) ;
2410  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
2411  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
2412 
2413  c1->cd(i%9 + 1);
2414  c1->cd(i%9 + 1)->SetLogy();
2415  if(i==0)
2416  {
2417  leg->AddEntry(hRefDistr,"mean of good","l");
2418  leg->AddEntry(hCell,"current","l");
2419  }
2420  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
2421  candidate = CheckDistribution(hCell,hRefDistr);
2422  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
2423  if(version>2)//..These are ratio plots of energy distr. of cell and mean of all good cells
2424  {
2425  hCell->Divide(hRefDistr);
2426  }
2427  //.. save histograms to file if you want to double check the output
2428  if(fTestRoutine ==1)
2429  {
2430  if(version==1) fOutputListBad->Add(hCell);
2431  if(version==10)fOutputListBadRatio->Add(hCell);
2432  if(version==2) fOutputListGood->Add(hCell);
2433  if(version==20)fOutputListGoodRatio->Add(hCell);
2434  }
2435  hCell->SetLineColor(kBlue+1);
2436  hCell->GetXaxis()->SetTitle("E (GeV)");
2437  hCell->GetYaxis()->SetTitle("N Entries");
2438  hCell->GetXaxis()->SetRangeUser(0.,10.);
2439  hCell->SetLineWidth(1) ;
2440  hCell->SetTitle(title);
2441  hRefDistr->SetLineColor(kGray+2);
2442  hRefDistr->SetLineWidth(1);
2443 
2444  hCell->Draw("hist");
2445 
2446  if(version==1 || version==2)hRefDistr->Draw("same hist") ;
2447 
2448  //..Mark the histogram that could be miscalibrated and labelled as warm
2449  if(candidate==1 && (version==1 || version==10))
2450  {
2451  gPad->SetFrameFillColor(kYellow-10);
2452  text->Draw();
2453  }
2454  if(version==1)
2455  {
2456  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
2457  textA->DrawLatex(0.65,0.62,Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
2458  }
2459  if(candidate==0 && (version==2 || version==20))
2460  {
2461  gPad->SetFrameFillColor(kYellow-10);
2462  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
2463  leg->Draw();
2464  }
2465  if(version<2)leg->Draw();
2466  }
2467 
2468  if(channelVector.size()<9 || cell == fNoOfCells-1)
2469  {
2470  internal_pdfName +=")";
2471  c1->Print(internal_pdfName.Data());
2472  }
2473  else if(firstCanvas==0)
2474  {
2475  internal_pdfName +="(";
2476  c1->Print(internal_pdfName.Data());
2477  firstCanvas=1;
2478  }
2479  else
2480  {
2481  c1->Print(internal_pdfName.Data());
2482  }
2483  delete c1;
2484  delete leg;
2485  channelVector.clear();
2486  }
2487  }
2488  cout<<endl;
2489  delete hRefDistr;
2490  delete text;
2491  delete text2;
2492  delete textA;
2493  //..Add the subdirectories to the file
2494  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
2495  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
2496  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
2497  if(version==20)fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
2498 
2499  if(fPrint==1)cout<<endl;
2500 }
2504 //_________________________________________________________________________
2506 {
2507  TH1D* hGoodAmp;
2508  TH1D* hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean");
2509  hgoodMean->Reset();
2510  Int_t NrGood=0;
2511 
2512  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2513  {
2514  if(warmIn==0 && fFlag[cell]!=0 )continue;
2515  if(warmIn==1 && fFlag[cell]!=0 && fWarmCell[cell]==0)continue;
2516  if(warmIn==2 && fWarmCell[cell]==0)continue;
2517  NrGood++;
2518 
2519  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
2520  hgoodMean->Add(hGoodAmp);
2521  }
2522  hgoodMean->Scale(1.0/NrGood);
2523 
2524  return hgoodMean;
2525 }
2536 //_________________________________________________________________________
2538 {
2539  TString Name = Form("%s_ratio",histogram->GetName());
2540  TH1* ratio=(TH1*)histogram->Clone(Name);
2541  ratio->Divide(reference);
2542 
2543  Double_t percentageOfLast=0.6;
2544  Double_t higherThanMean=5;
2545  Double_t highestRatio=1000;
2546  Double_t maxMean=10;
2547  Double_t minMean=0.1;
2548  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
2549 
2550  //..By default each cell is a candidate for recalibration
2551  Bool_t candidate=1;
2552  //..Find bin where reference has value 1, and the corresponding x-value
2553  Double_t totalevents = fProcessedEvents->Integral();
2554  Int_t binHeightOne = reference->FindLastBinAbove(1.0/totalevents);//..look at the spectrum until there is 1hit/event GeV
2555  Double_t binCentreHeightOne = reference->GetBinCenter(binHeightOne);
2556  Double_t thirdBinCentre = reference->GetBinCenter(3);
2557 
2558  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2559  //..Check the histogram
2560  //..Different checks to see whether the
2561  //..cell is really bad. Set candidate to 0.
2562 
2563  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2564  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
2565  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
2566  {
2567  candidate=0;
2568  //cout<<"1"<<endl;
2569  }
2570 
2571  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2572  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
2573  //ELI - TO DO: check that crieteria carfully - seems to work but not shure about it
2574  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
2575  if(ratio->GetMaximum()>highestRatio)//
2576  {
2577  candidate=0;
2578  //cout<<"2"<<endl;
2579  }
2580 
2581  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2582  //..check whether the ratio is much larger than 1
2583  //..calculate the mean in the relevant energy range
2584  Double_t mean=0;
2585  Int_t nullEntries=0;
2586  for(Int_t i=2;i<binHeightOne;i++)
2587  {
2588  if(ratio->GetBinContent(i)!=0)mean+=ratio->GetBinContent(i);
2589  else nullEntries++;
2590  }
2591  mean*=1.0/(binHeightOne-1-nullEntries);//..divide by number of bins (excluding bins without entries)
2592  if(mean>maxMean || mean<minMean)
2593  {
2594  candidate=0;
2595  //cout<<"3"<<endl;
2596  }
2597 
2598  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2599  //..check whether there are large spikes in the histogram
2600  //..compare bin values to mean of the ratio. If there is a bin value with
2601  //..content "higherThanMean" times lareger than mean it's losing it candidate status
2602  mean=0;
2603  //..Find the maximum in the mean range (0-binHeightOne)
2604  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
2605  Double_t localMaxBin=ratio->GetMaximumBin();
2606 
2607  for(Int_t i=2;i<binHeightOne;i++)
2608  {
2609  //..Exclude 0 bins and exclude bins near the maximum
2610  if(ratio->GetBinContent(i)<=0) continue;
2611  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
2612  mean+=ratio->GetBinContent(i);
2613  }
2614  mean*=1.0/(binHeightOne-1);//..divide by number of bins
2615  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
2616  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
2617  if(ratio->GetMaximum()>mean*higherThanMean)
2618  {
2619  candidate=0;
2620  //cout<<"4"<<endl;
2621  }
2622 
2623  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2624  //..Check for a cliff at 4GeV, happens in some cases
2625  Double_t beforeCliff=0;
2626  Double_t afterCliff=0;
2627  Int_t binsBefore=0;
2628  Int_t binsAfter=0;
2629  Int_t cliffBin = ratio->FindBin(4);
2630  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
2631  {
2632  //..Exclude 0 bins
2633  if(ratio->GetBinContent(i)<=0)continue;
2634  if(i<=cliffBin) binsBefore++;
2635  if(i>cliffBin) binsAfter++;
2636  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
2637  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
2638  beforeCliff*=1.0/binsBefore;
2639  afterCliff*=1.0/binsAfter;
2640  }
2641  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
2642  {
2643  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
2644  if(beforeCliff!=0 && afterCliff!=0)cout<<"5"<<endl;
2645  }
2646 
2647  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2648  //..Employ peak finder
2649  /* if(candidate==1)
2650  {
2651  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
2652  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
2653  //GetNPeaks();
2654  //cout<<"found N peaks: "<<nfound<<endl;
2655  }
2656  */
2657  return candidate;
2658 }
2659 
2668 //_________________________________________________________________________
2670 {
2671  //..TRD support structure
2672  //..(determined by eye, could be improved, but is already very acurate):
2673  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
2674  //..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
2675  Bool_t coveredByTRDSupportStruc=0;
2676 
2677  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
2678  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
2679  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
2680  )
2681  {
2682  coveredByTRDSupportStruc=1;
2683  }
2684  return coveredByTRDSupportStruc;
2685 }
2694 //_________________________________________________________________________
2696 {
2697  gStyle->SetPalette(55); //kRainBow==55
2698  gStyle->SetPalette(91); //kPastel==91
2699  //..build two dimensional histogram with values row vs. column
2700  TString histoName;
2701  histoName = Form("2DChannelMap_Flag%d_V%i",flagBegin,fTrial);
2702  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100_V%i",fTrial);
2703 
2704  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
2705  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
2706  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
2707 
2708  Int_t cellColumn=0,cellRow=0;
2709  Int_t cellColumnAbs=0,cellRowAbs=0;
2710  Int_t trash;
2711 
2712  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2713  {
2714  //..Do that only for cell ids also accepted by the code
2715  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
2716 
2717  //..Get Row and Collumn for cell ID c
2718  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
2719  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
2720  {
2721  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
2722  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
2723  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
2724  }
2725  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
2726  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
2727  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
2728 
2729 
2730  }
2731  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
2732  c1->ToggleEventStatus();
2733  c1->cd();
2734  //lowerPadRight->SetLeftMargin(0.09);
2735  //lowerPadRight->SetRightMargin(0.06);
2736  plot2D->Draw("colz");
2737 
2738  TLatex* text = 0x0;
2739  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
2740  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
2741  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
2742  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
2743  text->SetTextSize(0.06);
2744  text->SetNDC();
2745  text->SetTextColor(1);
2746  text->Draw();
2747 
2748  c1->Update();
2749  TString name =Form("%s/%s/%s_%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data() , histoName.Data());
2750  c1->SaveAs(name);
2752  fRootFile->WriteObject(plot2D,plot2D->GetName());
2753 
2754  /*
2755  delete plot2D;
2756  delete c1;
2757  delete text;
2758  */
2759 }
2763 //_________________________________________________________________________
2765 {
2766  char name[100];
2767  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
2768  {
2769  sprintf(name, "Cell %d",cell) ;
2770  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
2771  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
2772  }
2773  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
2774 }
2775 
2780 //_________________________________________________________________________
2782 {
2783  if(fTrackCellRecord==1)
2784  {
2785  Int_t zeroFlag = std::count(fFlag.begin(), fFlag.end(), 0);
2786  Int_t zeroWarm = std::count(fWarmCell.begin(), fWarmCell.end(), 0);
2787 
2788  cout<<"******* Debug No "<<number<<" ********"<<endl;
2789  cout<<"*** fCriterionCounter: "<<fCriterionCounter<<endl;
2790  cout<<"*** number of period analyses: "<<fAnalysisVector.size()<<endl;
2791  cout<<"*** number of man mask cells: "<<fManualMask.size()<<endl;
2792  cout<<"*** Bad+Dead(!0) fFlag elements: "<<fFlag.size()-zeroFlag<<endl;
2793  cout<<"*** Warm(!0) fWarmCell elements: "<<fWarmCell.size()-zeroWarm<<endl;
2794  cout<<"*** "<<endl;
2795  }
2796 }
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
void FlagAsBad_Time(Int_t crit, TH1F *inhisto, Double_t nSig=3)
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:98
TString fWorkdir
Directory which contains the folders fMergeOutput, fAnalysisInput and fAnalysisOutput. By default it is &#39;./&#39;.
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.
std::vector< Int_t > fFlag
! fFlag[CellID] = 0 (ok),1 (dead),2 and higher (bad certain criteria) start at 0 (cellID 0 = histobin...
TString fRunListFileName
This is the name of the file with the run numbers to be merged, by default it&#39;s &#39;runList.txt&#39;.
TList * fOutputListBad
! list with bad channel amplitudes, stored in fRootFile
Double_t fEndLowerBound
Lower bound.
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 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.
Bool_t fRunBRunMap
Produce truely run-by-run maps in a separate folder.
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:99
Int_t fNMaxColsAbs
Maximum No of colums in Calorimeter.
Int_t fNoOfCells
Number of cells in EMCal and DCal.
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.
std::vector< Bool_t > fWarmCell
! fWarmCell[CellID] = 0 (really bad), fWarmCell[CellID] = 1 (candidate for warm), ...
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
std::vector< Int_t > fSMMask
! fSMMask is filled with SM numbers that need to be masked by hand.
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:97
void AddMaskSM(Int_t iSM)
TString fAnalysisOutputRbR
For a compact summary of true run-by-run BC maps.
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.)
Bool_t fTrackCellRecord
Track the non-zero elements in the flags throughout the routine.
TFile * file
TList with histograms for a given trigger.
Int_t fTrial
Number of trial that this specific analyis is. By default &#39;0&#39; so one can try different settings witho...
TH1F * BuildHitAndEnergyMeanScaled(Int_t crit, Double_t emin=0.1, Double_t emax=2.)
Bool_t IsCoveredByTRD(Int_t row, Int_t collumn)
Int_t fStartCell
ID of the first cell you want to check.
Int_t fNMaxCols
Maximum No of colums in module (eta direction)
void PrintCellInfo(Int_t number)
bool Bool_t
Definition: External.C:53
TString fTrigger
Selected trigger for the analysis.
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.
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