AliPhysics  df4bbdf (df4bbdf)
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 
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 
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 
591  //..access the standart root output of a bad channel analysis to
592  //..get the necessary histogram
593  TString extRootFileName = Form("./AnalysisOutput/%s/%s/%s",fPeriod.Data(),fTrainNo.Data(),fExternalBadMapName.Data());
594  TFile* outputRoot = TFile::Open(extRootFileName.Data());
595 
596  TH1F* hFlags =(TH1F*)outputRoot->Get("fhCellFlag");
597 
598  //..set info from external source
599  //..Direction of cell ID
600  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
601  {
602  Double_t extFlag = hFlags->GetBinContent(cell+1);
603  //..Cell flagged as dead.
604  //..Flag only if it hasn't been flagged before
605  fFlag[cell] =extFlag;
606  }
607 }
613 //_________________________________________________________________________
615 {
616  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
617  //.. BAD CELLS
618  //.. Flag bad cells with fFlag= 2,3,4,5.. etc
619  //.. this excludes cells from subsequent analysis
620  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
621  TArrayD periodArray;
622  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
623  {
624  periodArray=fAnalysisVector.at(i);
625  PeriodAnalysis(periodArray.At(0),periodArray.At(1),periodArray.At(2),periodArray.At(3));
627  if(fPrint==1)cout<<endl;
628  if(fPrint==1)cout<<endl;
629  }
630  cout<<"o o o End of bad channel analysis o o o"<<endl;
631 }
632 //
633 // Mask an entire SM before doing the BC analysis
634 // This is useful when you get info from QA that there are problems with one SM
635 // and you want to clean up your bad channels beforehand
636 //
637 //________________________________________________________________________
639 {
640  cout<<"o o o Manually mask SM "<<iSM<<" o o o"<<endl;
641  //..Loop over cell ID
642  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
643  {
644  //..check to which SM the cell belongs
645  if(cell>=fStartCellSM[iSM] && cell<fStartCellSM[iSM+1])
646  {
647  fFlag[cell] =1;
648  }
649  }
650 }
659 //________________________________________________________________________
661 {
662  TArrayD periodArray;
663  periodArray.Set(4);
664  periodArray.AddAt((Double_t)criteria,0);
665  periodArray.AddAt(nsigma,1);
666  periodArray.AddAt(emin,2);
667  periodArray.AddAt(emax,3);
668  fAnalysisVector.push_back(periodArray);
669 }
685 //____________________________________________________________________
687 {
688  //.. criterion should be between 1-5
689  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;
690  if(fPrint==1)cout<<"o o o PeriodAnalysis for flag "<<criterion<<" o o o"<<endl;
691  if(fPrint==1)cout<<"o o o This is PeriodAnalysis No. "<<fCriterionCounter<<" in your list o o o"<<endl;
692  if(fPrint==1 && criterion != 3)cout<<"o o o Done in the energy range E "<<emin<<" - "<<emax<<endl;
693  if(fPrint==1 && criterion == 3)cout<<"o o o Done in the time range t "<<emin<<" - "<<emax<<endl;
694 
695  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
696  //.. ANALYSIS OF CELLS WITH ENTRIES
697  //.. Build average distributions and fit them
698  //.. Three tests for bad cells:
699  //.. 1) Average energy per hit
700  //.. 2) Average hit per event
701  //.. 3) Average max of cell time distribution
702  //.. 4) To be implemented: Scaled energy per hit distribution
703  //.. 5) Scaled hit distribution
704  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
705  TH1F* histogram;
706  if(fPrint==1)cout<<"o o o Analyze average cell distributions o o o"<<endl;
707  //..For case 1 or 2
708  if(criterion < 3) histogram = BuildHitAndEnergyMean(criterion, emin, emax);
709  if(criterion > 3) histogram = BuildHitAndEnergyMeanScaled(criterion, emin, emax);
710  if(criterion == 3) histogram = BuildTimeMean(criterion, emin, emax); //.. in case of crit 3 emin is tmin and emax is tmax
711 
712  if(criterion==1 || criterion == 4)
713  {
714 // if(emin>=1.8)FlagAsBad(criterion, histogram, nsigma, -1);//..do not apply a lower boundary
715  if(emin>=fEndLowerBound)FlagAsBad(criterion, histogram, nsigma, -1);//..do not apply a lower boundary
716  else FlagAsBad(criterion, histogram, nsigma, 200);//400
717  }
718  if(criterion==2 || criterion == 5)
719  {
720 // if(emin>=1.8)FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
721  if(emin>=fEndLowerBound)FlagAsBad(criterion, histogram, nsigma, -1);//..do not narrow the integration window
722  else FlagAsBad(criterion, histogram, nsigma, 601);
723  }
724  if(criterion==3) FlagAsBad_Time(criterion, histogram, nsigma);
725 
726 // if(histogram) delete histogram;
727 }
728 
737 //_________________________________________________________________________
739 {
740  if(fPrint==1)cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
741  TH1F *histogram;
742  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);
743  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);
744  histogram->SetXTitle("Abs. Cell Id");
745  if(crit==1)histogram->SetYTitle("Energy per hit");
746  if(crit==2)histogram->SetYTitle("Number of hits in cell");
747  histogram->GetXaxis()->SetNdivisions(510);
748 
749  //..count Events in Energy Range
750  TH1F* pojection = (TH1F*)fCellAmplitude->ProjectionX("Intermediate");
751  fnEventsInRange = pojection->Integral(pojection->GetBinContent(emin),pojection->GetBinContent(emax));
752 
753  //..here the average hit per event and the average energy per hit is caluclated for each cell.
754  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
755  {
756  Double_t Esum = 0;
757  Double_t Nsum = 0;
758 
759  for (Int_t j = 1; j <= fCellAmplitude->GetNbinsX(); j++)
760  {
761  //To Do: I think one should also take the different bin width into account
762  Double_t E = fCellAmplitude->GetXaxis()->GetBinCenter(j);
763  Double_t N = fCellAmplitude->GetBinContent(j, cell+1);
764  //..investigate only cells that were not flagged as dead ore bad
765  if (E < emin || E > emax || fFlag[cell]!=0) continue;
766  Esum += E*N;
767  Nsum += N;
768  }
769  //..Set the values only for cells that are not yet marked as bad
770  if(fFlag[cell]==0)
771  {
772  if(crit==2) histogram->SetBinContent(cell+1, Nsum); //..number of hits (per event -> not anymore)
773  if(Nsum > 0. && crit==1)histogram->SetBinContent(cell+1, Esum/(Nsum)); //..average energy per hit
774  }
775  }
776  return histogram;
777 }
778 
779 
780 //-------------------------------------------------------------------------------------------
781 // This function builds the scaled hit distribution
782 // It normalizes the hits/cell to the mean value of the row and the column of the cell
783 // this is done in an iterative procedure (about 5 iterations are needed)
784 // with this procedure, the effects of the TRD and SM structures on the EMCal can be minimized
785 // The output is a histogram with the hits/cell as a function of cell ID.
786 
787 // ONLY IMPLEMENTED FOR THE HITS PER CELL method
788 
792 // ------------------------------------------------------------------------------------------
794 {
795  if(fPrint==1)cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
796  TH1F *histogram = NULL;
797  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);
798  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);
799  histogram->SetXTitle("Abs. Cell Id");
800  if(crit==4)histogram->SetYTitle("Energy per hit");
801  if(crit==5)histogram->SetYTitle("Number of hits in cell");
802  histogram->GetXaxis()->SetNdivisions(510);
803 
804 
806  TH2F *hEnergyScaled = (TH2F*) fCellAmplitude->Clone("hEnergyScaled");
807  TH1F *hEnergyCol = new TH1F("hEnergyCol", "hEnergyCol", 100,0.,100.);
808  TH1F *hEnergyRow = new TH1F("hEnergyRow", "hEnergyRow", 250,0.,250.);
809 
810 
812  Int_t cellCol, cellRow, trash, col, row;
813  Double_t dCellEnergy;
814  Double_t dNumOfHits = 0.;
815  Double_t *arrHits = new Double_t[fNoOfCells];
816 
818  Double_t NoOfHits = 0.;
819  for(int iCell = 1; iCell <= fNoOfCells; iCell++){
820  NoOfHits = 0.;
821  if(fFlag[iCell -1] == 0){
822  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
823  NoOfHits += hEnergyScaled->GetBinContent(k, iCell);
824  }
825  }
826  arrHits[iCell - 1] = NoOfHits;
827  }
828 
829 
830  // ....................................................................................................
831  // Find cells that are really bad and exclude them from the calculation of the mean column and mean row
832  // ....................................................................................................
833 
834  // Declare Histos for BC finding for scaling
835  TH1D *hHitDistrib_forScaling = new TH1D("hHitDistrib_forScaling","hHitDistrib_forScaling",1000,0.,TMath::Median(fNoOfCells,arrHits)*4);
836  TF1 *fgausForScaling = new TF1("fgausForScaling","gaus", 0.,TMath::Median(fNoOfCells,arrHits)*4);
837  TVectorD fFlagForScaling;
838  fFlagForScaling.ResizeTo(fNoOfCells);
839  Double_t dhits = 0.;
840  for(int i = 1; i <= fNoOfCells; i++){
841  dhits = 0.;
842  dNumOfHits = 0.;
843  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
844  if(crit==5){dhits += hEnergyScaled->GetBinContent(k, i);}
845  if(crit==4){
846  dhits += hEnergyScaled->GetBinContent(k, i + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(k);
847  dNumOfHits += hEnergyScaled->GetBinContent(k, i + 1);
848  }
849  }
850  if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
851  if(dhits > 0.1){
852  hHitDistrib_forScaling->Fill(dhits);
853  }
854  }
855 
857  fgausForScaling->SetParameter(1, hHitDistrib_forScaling->GetBinCenter(hHitDistrib_forScaling-> GetMaximumBin()));
858  hHitDistrib_forScaling->Fit(fgausForScaling,"MQ0");
859  hHitDistrib_forScaling->Fit(fgausForScaling,"MQ0");
860 
861  for(int iCell = 1; iCell <= fNoOfCells; iCell++){
862  dhits = 0.;
863  dNumOfHits = 0.;
864  for(int k = hEnergyScaled->GetXaxis()->FindBin(emin); k <= hEnergyScaled->GetXaxis()->FindBin(emax); k++){
865  if(crit==5){dhits += hEnergyScaled->GetBinContent(k, iCell);}
866  if(crit==4){
867  dhits += hEnergyScaled->GetBinContent(k, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(k);
868  dNumOfHits += hEnergyScaled->GetBinContent(k, iCell + 1);
869  }
870  }
871  if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
872  if(dhits > fgausForScaling->GetParameter(1) + 5* fgausForScaling->GetParameter(2) || dhits < fgausForScaling->GetParameter(1) - 5* fgausForScaling->GetParameter(2)){
873  fFlagForScaling[iCell - 1] = 1;
874  }
875 
876  }
877  delete hHitDistrib_forScaling;
878  delete fgausForScaling;
879  delete [] arrHits;
880 
881  TF1 *fFitCol;
882  TF1 *fFitRow;
883  //...........................................
884  //start iterative process of scaling of cells
885  //...........................................
886  for(int iter = 0; iter < 5; iter++){
887  // array of vectors for calculating the mean hits per col/row
888  std::vector<Double_t> vecCol[100];
889  std::vector<Double_t> vecRow[250];
890 
891 
892  // TH2F *hmap = new TH2F("","",100,0.,100,250,0.,250.);
893  //
894  // // Build the Hitmap
895  // for(int iCell = 1; iCell <= fNoOfCells; iCell++){
896  // if(fFlag[iCell -1 ] == 0 && fFlagForScaling[iCell - 1] == 0){
897  // fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell - 1 ,0,cellCol, cellRow, trash, col, row);
898  // dhits = 0.;
899  // dNumOfHits = 0.;
900  // for(int EBin = hEnergyScaled->GetXaxis()->FindBin(emin); EBin <= hEnergyScaled->GetXaxis()->FindBin(emax); EBin++){
901  // if(crit==5){dhits += hEnergyScaled->GetBinContent(EBin, Cell);}
902  // if(crit==4){
903  // dhits += hEnergyScaled->GetBinContent(EBin, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(EBin);
904  // dNumOfHits += hEnergyScaled->GetBinContent(EBin, iCell + 1);
905  // }
906  // }
907  // if(crit==4 && dNumOfHits > 0){dhits = dhits / dNumOfHits;}
908  // hmap->SetBinContent(col + 1, row + 1, dhits );
909  // }
910  // }
911 
912  // Calculate the mean number of hits of each row and column
913  for(int iCell = 0; iCell < fNoOfCells; iCell++){
914  if(fFlag[iCell] == 0 && fFlagForScaling[iCell] == 0){
915  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
916  dCellEnergy = 0.;
917  dNumOfHits = 0.;
918  for (int EBin = hEnergyScaled->GetXaxis()->FindBin(emin); EBin < hEnergyScaled->GetXaxis()->FindBin(emax); EBin++) {
919  if(crit==5){dCellEnergy += hEnergyScaled->GetBinContent(EBin, iCell + 1);}
920  if(crit==4){
921  dCellEnergy += hEnergyScaled->GetBinContent(EBin, iCell + 1)*hEnergyScaled->GetXaxis()->GetBinCenter(EBin);
922  dNumOfHits += hEnergyScaled->GetBinContent(EBin, iCell + 1);
923  }
924 
925  }
926  if(crit==4 && dNumOfHits > 0){dCellEnergy = dCellEnergy / dNumOfHits;}
927 
928  if( dCellEnergy > 0.){
929  hEnergyCol->Fill(col + 0.5, dCellEnergy);
930  hEnergyRow->Fill(row + 0.5, dCellEnergy);
931  vecCol[col].push_back(dCellEnergy);
932  vecRow[row].push_back(dCellEnergy);
933  }
934  }
935  }
936 
937  // Fill the histogram: mean hit per column
938  for(int iCol = 1; iCol <= hEnergyCol->GetNbinsX() ; iCol++){
939  if(vecCol[iCol -1].size() > 0.){
940  hEnergyCol->SetBinContent(hEnergyCol->FindBin(iCol - 0.5), hEnergyCol->GetBinContent(hEnergyCol->FindBin(iCol - 0.5))/vecCol[iCol-1].size() );
941  }
942  }
943 
944  // Fill the histogram: mean hit per row
945  for(int iRow = 1; iRow <= hEnergyRow->GetNbinsX() ; iRow++){
946  if(vecRow[iRow -1].size() > 0.){
947  hEnergyRow->SetBinContent(hEnergyRow->FindBin(iRow - 0.5), hEnergyRow->GetBinContent(hEnergyRow->FindBin(iRow - 0.5))/vecRow[iRow-1].size() );
948  }
949  }
950 
951  // Global fit to hits per row
952  fFitRow = new TF1("fFitRow","[0]",0.,250.);
953  fFitRow->SetParameter(0, hEnergyRow->GetBinContent(10));
954  Double_t MeanRow = hEnergyRow->GetBinContent(10);
955 
956  // Global fit to hits per column
957  fFitCol = new TF1("fFitCol","[0]",0.,100.);
958  fFitCol->SetParameter(0, hEnergyCol->GetBinContent(10));
959  Double_t MeanCol = hEnergyCol->GetBinContent(10);
960 
961 
962  //Scale each cell by the deviation of the mean of the column and the global mean
963  for(int iCell = 0; iCell < fNoOfCells; iCell++){
964  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
965  if (hEnergyCol->GetBinContent(col + 1) > 0.) {
966  for(int EBin = 1; EBin < hEnergyScaled->GetNbinsX(); EBin++) {
967  hEnergyScaled->SetBinContent(EBin, iCell + 1, hEnergyScaled->GetBinContent(EBin, iCell + 1) * (MeanCol / hEnergyCol->GetBinContent(col + 1)));
968  }
969  }
970  }
971 
972  //Scale each cell by the deviation of the mean of the row and the global mean
973  for(int iCell = 0; iCell < fNoOfCells; iCell++){
974  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(iCell ,0,cellCol, cellRow, trash, col, row);
975  if (hEnergyRow->GetBinContent(row + 1) > 0.) {
976  for(int EBin = 1; EBin < hEnergyScaled->GetNbinsX(); EBin++) {
977  hEnergyScaled->SetBinContent(EBin, iCell + 1, hEnergyScaled->GetBinContent(EBin, iCell + 1) * (MeanRow / hEnergyRow->GetBinContent(row + 1)));
978  }
979  }
980  }
981  //....................
982  // iteration ends here
983  //....................
984  }
985 
986  //............................................................................................
987  //..here the average hit per event and the average energy per hit is caluclated for each cell.
988  //............................................................................................
989  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
990  {
991  Double_t Esum = 0;
992  Double_t Nsum = 0;
993 
994  for (Int_t j = 1; j <= hEnergyScaled->GetNbinsX(); j++)
995  {
996  //To Do: I think one should also take the different bin width into account
997  Double_t E = hEnergyScaled->GetXaxis()->GetBinCenter(j);
998  Double_t N = hEnergyScaled->GetBinContent(j, cell+1);
999  //..investigate only cells that were not flagged as dead ore bad
1000  if (E < emin || E > emax || fFlag[cell]!=0) continue;
1001  Esum += E*N;
1002  Nsum += N;
1003  }
1004  //..Set the values only for cells that are not yet marked as bad
1005  if(fFlag[cell]==0)
1006  {
1007  if(crit==5) histogram->SetBinContent(cell+1, Nsum); //..number of hits (per event -> not anymore)
1008  if(Nsum > 0. && crit==4)histogram->SetBinContent(cell+1, Esum/(Nsum)); //..average energy per hit
1009  }
1010  }
1011  delete hEnergyScaled;
1012  delete hEnergyCol;
1013  delete hEnergyRow;
1014  delete fFitCol;
1015  delete fFitRow;
1016 
1017  return histogram;
1018 }
1019 
1020 
1021 
1026 //_________________________________________________________________________
1028 {
1029  if(fPrint==1)cout<<" o Calculate maximum of cell time distribution "<<endl;
1030  TString name;
1031  TH1F *histogram = NULL;
1032  histogram = new TH1F(Form("timeMax_t%.2f-%.2f",tmin,tmax),Form("Maximum of time distr., %.2f < t < %.2f ns",tmin,tmax), fNoOfCells,0,fNoOfCells);
1033  histogram->SetXTitle("Abs. Cell Id");
1034  histogram->SetYTitle(Form("N_{cells}^{%.0f<t<%.0f}/N_{cells}^{all}",tmin,tmax));
1035  histogram->GetXaxis()->SetNdivisions(510);
1036 
1037  //..Time information
1038  Double_t dHits = 0.;
1039  Double_t dIn = 0.;
1040  Double_t dAll = 0.;
1041  for(int iCell = 1; iCell <= fNoOfCells; iCell ++)
1042  {
1043  if(fFlag[iCell - 1] != 0) { continue; }
1044  dIn = 0.;
1045  dAll = 0.;
1046  for(int iTime = 1; iTime <= fCellTime->GetNbinsX(); iTime++)
1047  {
1048  dHits = fCellTime->GetBinContent(iTime, iCell);
1049  if(fCellTime->GetXaxis()->GetBinCenter(iTime) > tmin && fCellTime->GetXaxis()->GetBinCenter(iTime) < tmax)
1050  {
1051  dIn += dHits;
1052  }
1053  dAll += dHits;
1054  }
1055  if(dIn > 0.)
1056  {
1057  histogram->SetBinContent(iCell, dAll/dIn);
1058  }
1059  else
1060  {
1061  histogram->SetBinContent(iCell, 0.);
1062  }
1063  }
1064  return histogram;
1065 }
1066 
1071 //_________________________________________________________________________
1073 {
1074  Int_t sumOfExcl=0;
1075  Int_t manualMaskCounter=0;
1076  //..sort the elements in manual mask list
1077  std::sort (fManualMask.begin(), fManualMask.end());
1078 
1079  //..Direction of cell ID
1080  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1081  {
1082  Double_t nSum = 0;
1083  //..Direction of amplitude (Checks energies from 0-nBins GeV)
1084  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1085  {
1086  //..cellID+1 = histogram bin
1087  Double_t N = fCellAmplitude->GetBinContent(amp,cell+1);
1088  nSum += N;
1089  }
1090  //..If the amplitude in one cell is basically 0
1091  //..mark the cell as excluded
1092  if(nSum == 0 && fFlag[cell]==0)
1093  {
1094  //..Cell flagged as dead.
1095  //..Flag only if it hasn't been flagged before
1096  fFlag[cell] =1;
1097  sumOfExcl++;
1098  }
1099  //..add here the manual masking
1100  if(manualMaskCounter<(Int_t)fManualMask.size() && fManualMask.at(manualMaskCounter)==cell)
1101  {
1102  fFlag[cell] = 2;
1103  manualMaskCounter++;
1104  }
1105  }
1106  if(fPrint==1)cout<<" o Number of dead cells: "<<sumOfExcl<<endl;
1107  if(fPrint==1)cout<<" ("<<sumOfExcl<<")"<<endl;
1108 }
1109 
1124 //_________________________________________________________________________
1125 void BadChannelAna::FlagAsBad(Int_t crit, TH1F* inhisto, Double_t nsigma, Double_t dnbins)
1126 {
1127  gStyle->SetOptStat(0); //..Do not plot stat boxes
1128  gStyle->SetOptFit(0); //
1129  if(fPrint==1 && crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
1130  if(fPrint==1 && crit==2)cout<<" o Fit average hit per event distribution"<<endl;
1131  if(fPrint==1 && crit==3)cout<<" o Fit average hit maximum distribution"<<endl;
1132 
1133  Int_t cellColumn=0,cellRow=0;
1134  Int_t cellColumnAbs=0,cellRowAbs=0;
1135  Int_t trash;
1136 
1137  TString histoName=inhisto->GetName();
1138  Double_t goodmax = 0. ;
1139  Double_t goodmin = 0. ;
1140  //..set a user range so that the min and max is only found in a certain range
1141  inhisto->GetXaxis()->SetRangeUser(fStartCell,fNoOfCells);//..
1142  Double_t dminVal = inhisto->GetMinimum(0);
1143  Double_t dmaxVal = inhisto->GetMaximum();
1144  Double_t inputBins=dnbins;
1145  Double_t medianOfHisto = 0.;
1146  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1147  //. . .determine settings for the histograms (range and binning)
1148  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1149 
1150  //- - - - - - - - - - - -
1151  //..Determine good min and max range for histograms
1152  //..calculate and print the "median"
1153  //..
1154  Int_t numBins = inhisto->GetXaxis()->GetNbins();
1155 
1156  numBins-=fStartCell;
1157  Double_t *x = new Double_t[numBins];
1158  Double_t* y = new Double_t[numBins];
1159  for (int i = 0; i < numBins; i++)
1160  {
1161  x[i] = inhisto->GetBinContent(i+fStartCell);
1162  y[i] = 1; //each value with weight 1
1163  if(x[i]==0)y[i] = 0;
1164  }
1165  medianOfHisto = TMath::Median(numBins,x,y);
1166 
1167  //..if dmaxVal is too far away from medianOfHisto the histogram
1168  //..range will be too large -> reduce the range
1169  //cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1170  if(medianOfHisto*10<dmaxVal)
1171  {
1172  if(medianOfHisto*100<dmaxVal)
1173  {
1174  dmaxVal=medianOfHisto+0.02*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1175  }
1176  else
1177  {
1178  dmaxVal=medianOfHisto+0.2*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1179  }
1180  //cout<<"- - - median too far away from max range"<<endl;
1181  //cout<<"corrected max value: "<<dmaxVal<<endl;
1182  }
1183 
1184  //- - - - - - - - - - - -
1185  //..Determine number of bins for the histogram
1186  //..find a proper binning automatically (mostly done to avoid steplike structures, 0 entries in bins etc)
1187  //cout<<"max value: "<<dmaxVal<<", min value: "<<dminVal<<endl;
1188  if(crit==2 || crit == 5)//..hits in cell
1189  {
1190  dnbins=dmaxVal-dminVal;
1191  if(dnbins>100)
1192  {
1193  if(dnbins>2000)dnbins=0.01*(dmaxVal-dminVal); //..maximum 5000 bins. changed to 3000 .. lets see..
1194  if(dnbins>2000)dnbins=0.001*(dmaxVal-dminVal);//..maximum 5000 bins.
1195  if(dnbins<100) dnbins=0.02*(dmaxVal-dminVal); //..minimum 100 bins.
1196  }
1197  }
1198  if(crit==1 || crit == 4) //..energy/hit
1199  {
1200  dnbins=(dmaxVal-dminVal)*150; //150 if you have problems with low statistic 500 normal stat
1201  if(dnbins<25)dnbins=50; //..this is the min. binning for E/hit distr.
1202  if(inputBins==-1)
1203  {
1204  //dnbins=200;
1205  dnbins=(dmaxVal-dminVal)*40;
1206  }
1207  }
1208  if(x) delete []x;
1209  if(y) delete []y;
1210 
1211  if(crit==3)
1212  {
1213  //..obtain the binwidth for
1214  Int_t binWidth = fCellTime->GetXaxis()->GetBinWidth(1);
1215  dnbins = fCellTime->GetXaxis()->GetNbins();
1216  dminVal= fCellTime->GetXaxis()->GetBinCenter(1)-(binWidth*1.0)/2.0;
1217  dmaxVal= fCellTime->GetXaxis()->GetBinCenter(dnbins)+(binWidth*1.0)/2.0;
1218  cout<<"set the new histo with following settings- bins: "<<dnbins<<", min val = "<<dminVal<<", max val:"<<dmaxVal<<endl;
1219  }
1220 
1221  if(dmaxVal<dminVal || dnbins<1)
1222  {
1223  cout<<"###############################"<<endl;
1224  cout<<"Big problem: min:"<<dminVal<<", dmaxVal"<<dmaxVal<<endl;
1225  cout<<"dnbins:"<<dnbins<<endl;
1226  cout<<"###############################"<<endl;
1227  }
1228  //..build histos
1229  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1230  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1231  distrib->SetYTitle("Entries");
1232  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistr_wTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1233  TH1F *distrib_woTRDStruc= new TH1F(Form("%sDistr_woTRD",(const char*)histoName), "", dnbins, dminVal, dmaxVal);
1234 
1235  //..build two dimensional histogram with values row vs. column
1236  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);
1237  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1238  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1239 
1240  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1241  //. . .build the distribution of average values
1242  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1243  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1244  {
1245  //..Do that only for cell ids also accepted by the code
1246  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1247 
1248  //..Fill histograms only for cells that are not yet flagged as bad
1249  if(fFlag[cell]==0)
1250  {
1251  //..fill the distribution of avarge cell values
1252  distrib->Fill(inhisto->GetBinContent(cell+1));
1253  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1254  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1255  //..Get Row and Collumn for cell ID
1256  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1257  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1258  {
1259  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1260  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1261  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1262  }
1263  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
1264  //..check TRD support structure
1265  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
1266  {
1267  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1268  }
1269  else
1270  {
1271  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1272  }
1273  }
1274  }
1275 
1276  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1277  //. . .draw histogram + distribution
1278  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1279  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
1280  c1->ToggleEventStatus();
1281  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
1282  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
1283  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
1284  upperPad->Draw();
1285  lowerPadLeft->Draw();
1286  lowerPadRight->Draw();
1287 
1288  upperPad->cd();
1289  upperPad->SetLeftMargin(0.05);
1290  upperPad->SetRightMargin(0.05);
1291  if(crit!=3)upperPad->SetLogy();
1292  inhisto->SetTitleOffset(0.6,"Y");
1293  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
1294  inhisto->GetYaxis()->SetTitleOffset(0.7);
1295  inhisto->SetLineColor(kBlue+1);
1296  inhisto->DrawCopy();
1297 
1298  lowerPadRight->cd();
1299  lowerPadRight->SetLeftMargin(0.09);
1300  lowerPadRight->SetRightMargin(0.12);
1301  plot2D->GetYaxis()->SetTitleOffset(1.3);
1302  plot2D->DrawCopy("colz");
1303 
1304  lowerPadLeft->cd();
1305  lowerPadLeft->SetLeftMargin(0.09);
1306  lowerPadLeft->SetRightMargin(0.07);
1307  lowerPadLeft->SetLogy();
1308  distrib->SetLineColor(kBlue+1);
1309  distrib->GetYaxis()->SetTitleOffset(1.3);
1310  //cout<<medianOfHisto<<endl;
1311  //Double_t drawRangeDown = 0.;//medianOfHisto - 0.5*medianOfHisto;
1312  //Double_t drawRangeUp = 2000.;//medianOfHisto + 0.5*medianOfHisto;
1313  //distrib->GetXaxis()->SetRangeUser(0., 2.);
1314  distrib->DrawCopy("");
1315  distrib_wTRDStruc->SetLineColor(kGreen+1);
1316  distrib_wTRDStruc->DrawCopy("same");
1317  distrib_woTRDStruc->SetLineColor(kMagenta+1);
1318  distrib_woTRDStruc->DrawCopy("same");
1319 
1320 
1321  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1322  //. . .fit histogram
1323  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1324 
1325  //..to exclude first 2 bins from max finding
1326 // distrib->GetXaxis()->SetRangeUser(distrib->GetBinCenter(2),distrib->GetBinCenter(dnbins));
1327  Double_t maxDistr = distrib->GetMaximum();
1328  Int_t maxBin = distrib->GetMaximumBin();
1329  Double_t maxCenter = distrib->GetBinCenter(maxBin);
1330 
1331  //..good range is around the max value as long as the
1332  //..bin content is larger than 2 entries
1333  for(Int_t i = maxBin ; i<=dnbins ; i++)
1334  {
1335  if(distrib->GetBinContent(i)<2) break ;
1336  goodmax = distrib->GetBinCenter(i);
1337  }
1338  for(Int_t i = maxBin ; i>2 ; i--)
1339  {
1340  if(distrib->GetBinContent(i)<2) break ;
1341  goodmin = distrib->GetBinLowEdge(i);
1342  }
1343 
1344  //if(maxBin<2)goodmin = distrib->GetBinLowEdge(2);
1345 
1346  //..Define min/max range of the fit function
1347  Double_t minFitRange=goodmin;
1348  Double_t maxFitRange=goodmax;
1349  //if(crit==2)minFitRange=distrib->GetBinLowEdge(2);//..exclude first 2 bins
1350  //if(crit==2)maxFitRange=dmaxVal;
1351  if(crit==3)minFitRange=-20;
1352  if(crit==3)maxFitRange=20;
1353 
1354  //cout<<"max bin: "<<maxBin<<", max value: "<<maxDistr<<endl;
1355  //cout<<"start mean: "<<maxCenter<<", mean range: 0-"<<dmaxVal<<endl;
1356  //cout<<"fit range: "<<minFitRange<<" - "<<maxFitRange<<endl;
1357  TF1 *fit2 = new TF1("fit2", "gaus",minFitRange,maxFitRange);
1358  //..start the fit with a mean of the highest value
1359  fit2->SetParLimits(0,0,maxDistr); //..do not allow negative ampl values
1360  fit2->SetParameter(1,maxCenter); //..set the start value for the mean
1361  fit2->SetParLimits(1,0,dmaxVal); //..do not allow negative mean values
1362 
1363  //..ELI - TO DO: eventually fit with TRD and without TRD seperatley
1364  distrib->Fit(fit2, "0LQEM", "", minFitRange, maxFitRange);
1365  Double_t sig, mean;// chi2ndf;
1366  mean = fit2->GetParameter(1);
1367  sig = fit2->GetParameter(2);
1368 
1369  //..for case 1 and 2 lower than 0 is an unphysical value
1370 // if(crit=!3 && mean <0.) mean=0.;
1371 
1372  goodmin = mean - nsigma*sig ;
1373  goodmax = mean + nsigma*sig ;
1374  //..for case 1 and 2 lower than 0 is an unphysical value
1375 // if(crit=!3 && goodmin <0.) goodmin=0.;
1376  //..this (below) is a special case for energy ranges where cell do not have many
1377  //..entries typically. One should not apply a lower bound in this case.
1378  if(inputBins==-1) goodmin=-1;
1379  if(goodmin <=0.) goodmin = -100;
1380  if(fPrint==1)cout<<" o Result of fit: "<<endl;
1381  if(fPrint==1)cout<<" o "<<endl;
1382  if(fPrint==1)cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
1383  if(fPrint==1)cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
1384 
1385  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1386  //. . .Add info to histogram
1387  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1388  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1389  lline->SetLineColor(kGreen+2);
1390  lline->SetLineStyle(7);
1391  lline->Draw();
1392 
1393  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1394  rline->SetLineColor(kGreen+2);
1395  rline->SetLineStyle(7);
1396  rline->Draw();
1397 
1398  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1399  leg->AddEntry(lline,"Good region boundary","l");
1400  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1401  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1402  leg->SetBorderSize(0);
1403  leg->Draw("same");
1404 
1405  fit2->SetLineColor(kOrange-3);
1406  fit2->SetLineStyle(1);//7
1407  fit2->Draw("same");
1408 
1409 
1410  TLatex* text = 0x0;
1411  if(crit==1 || crit==4) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1412  if(crit==2 || crit==5) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1413  if(crit==3) text = new TLatex(0.12,0.85,Form("Good range: %.2f-%.2f",goodmin,goodmax));
1414  text->SetTextSize(0.06);
1415  text->SetNDC();
1416  text->SetTextColor(1);
1417  text->Draw();
1418 
1419  upperPad->cd();
1420  TLine *uline = new TLine(0, goodmax,fNoOfCells,goodmax);
1421  uline->SetLineColor(kGreen+2);
1422  uline->SetLineStyle(7);
1423  uline->Draw();
1424  TLine *lowline = new TLine(0, goodmin,fNoOfCells,goodmin);
1425  lowline->SetLineColor(kGreen+2);
1426  lowline->SetLineStyle(7);
1427  lowline->Draw();
1428  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1429  //. . .Save histogram
1430  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1431  c1->Update();
1432  TString name =Form("%s/%s/criteria-_%d.gif",fWorkdir.Data(), fAnalysisOutput.Data(), crit);
1433  if(crit==1)name=Form("%s/%s/AverageEperHit_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1434  if(crit==2)name=Form("%s/%s/AverageHitperEvent_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1435  if(crit==3)name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1436  if(crit==4)name=Form("%s/%s/AverageEperHit_scaled_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1437  if(crit==5)name=Form("%s/%s/AverageHitperEvent_scaled_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1438  c1->SaveAs(name);
1439 
1440  fRootFile->WriteObject(c1,c1->GetName());
1441  fRootFile->WriteObject(plot2D,plot2D->GetName());
1442  fRootFile->WriteObject(distrib,distrib->GetName());
1443  fRootFile->WriteObject(inhisto,inhisto->GetName());
1444  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1445  //. . . Mark the bad cells in the fFlag array
1446  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1447  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1448  //. . .(0 by default - good cell)
1449  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1450  if(fPrint==1)cout<<" o Flag bad cells that are outside the good range "<<endl;
1451  for(Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1452  {
1453  //..cell=0 and bin=1, cell=1 and bin=2
1454  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1455  if(inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]==0)
1456  {
1457  fFlag[cell]=fCriterionCounter;
1458  }
1459  if(inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]==0)
1460  {
1461  fFlag[cell]=fCriterionCounter;
1462  }
1463  }
1464  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;
1465 
1466  /*
1467  delete distrib;
1468  delete plot2D;
1469  delete c1;
1470 
1471 
1472  delete fit2;
1473  delete lline;
1474  delete rline;
1475  delete leg;
1476  delete text;
1477  delete uline;
1478  delete lowline;
1479 */
1480  //delete upperPad;
1481  //delete lowerPadLeft;
1482  //delete lowerPadRight;
1483 }
1484 
1485 
1486 
1487 //.. This function flags bad cells according to their time distribution
1488 //.. The Cut is the upper limit for each cell and is set by hand in the runAnalysisBC macro
1489 //..
1493 
1494 void BadChannelAna::FlagAsBad_Time(Int_t crit, TH1F* inhisto, Double_t nSig)
1495 {
1496  Int_t cellColumn=0,cellRow=0;
1497  Int_t cellColumnAbs=0,cellRowAbs=0;
1498  Int_t trash;
1499 
1500  //..set a user range so that the min and max is only found in a certain range
1501  inhisto->GetXaxis()->SetRangeUser(fStartCell,fNoOfCells);//..
1502  Double_t dminVal = inhisto->GetMinimum(0);
1503  Double_t dmaxVal = inhisto->GetMaximum();
1504  //Double_t dnbins;
1505  Double_t medianOfHisto = 0.;
1506  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1507  //. . .determine settings for the histograms (range and binning)
1508  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1509 
1510  //..Determine good min and max range for histograms
1511  //..A. calculate and print the "median"
1512  Int_t numBins = inhisto->GetXaxis()->GetNbins();
1513 
1514  numBins-=fStartCell;
1515  Double_t *x = new Double_t[numBins];
1516  Double_t* y = new Double_t[numBins];
1517  for (int i = 0; i < numBins; i++)
1518  {
1519  x[i] = inhisto->GetBinContent(i+fStartCell);
1520  y[i] = 1; //each value with weight 1
1521  if(x[i]==0)y[i] = 0;
1522  }
1523  medianOfHisto = TMath::Median(numBins,x,y);
1524 
1525  //..B. if dmaxVal is too far away from medianOfHisto the histogram
1526  //..range will be too large -> reduce the range
1527  cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1528  if(medianOfHisto*10<dmaxVal)
1529  {
1530  if(medianOfHisto*100<dmaxVal)
1531  {
1532  dmaxVal=medianOfHisto+0.002*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1533  }
1534  else
1535  {
1536  dmaxVal=medianOfHisto+0.02*(dmaxVal-medianOfHisto); //..reduce the distance between max and mean drastically to cut out the outliers
1537  }
1538  //cout<<"- - - median too far away from max range"<<endl;
1539  //cout<<"corrected max value: "<<dmaxVal<<endl;
1540  }
1541  cout<<"max value: "<<dmaxVal<<", min: "<<dminVal<<" median of histogram: "<<medianOfHisto<<endl;
1542  if(dmaxVal<dminVal)
1543  {
1544  cout<<"###############################"<<endl;
1545  cout<<"Big problem: min:"<<dminVal<<", dmaxVal"<<dmaxVal<<endl;
1546  cout<<"###############################"<<endl;
1547  }
1548 
1549  TString histoName=inhisto->GetName();
1550  TH1F *distrib = new TH1F(Form("%sDistr",(const char*)histoName), "", 300, dminVal, dmaxVal);
1551  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1552  distrib->SetYTitle("Entries");
1553  TH1F *distrib_wTRDStruc = new TH1F(Form("%sDistrWtrd",(const char*)histoName), "", 300, dminVal, dmaxVal);
1554  TH1F *distrib_woTRDStruc = new TH1F(Form("%sDistrWoTRD",(const char*)histoName), "", 300, dminVal, dmaxVal);
1555  //..build two dimensional histogram with values row vs. column
1556  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);
1557  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1558  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1559 
1560 // cout<<"histoName: "<<histoName<<endl;
1561 // cout<<"distrib: "<<distrib->GetName()<<endl;
1562  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1563  //. . .build the distribution of average values
1564  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1565  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1566  {
1567  //..Do that only for cell ids also accepted by the code
1568  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1569 
1570  //..Fill histograms only for cells that are not yet flagged as bad
1571  if(fFlag[cell]==0)
1572  {
1573  //..fill the distribution of avarge cell values
1574  distrib->Fill(inhisto->GetBinContent(cell+1));
1575  //if(cell<fCellStartDCal)distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1576  //else distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1577  //..Get Row and Collumn for cell ID
1578  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
1579  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
1580  {
1581  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1582  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1583  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1584  }
1585  plot2D->SetBinContent(cellColumnAbs,cellRowAbs,inhisto->GetBinContent(cell+1));
1586  //..check TRD support structure
1587  if(IsCoveredByTRD(cellRowAbs,cellColumnAbs)==1)
1588  {
1589  distrib_wTRDStruc->Fill(inhisto->GetBinContent(cell+1));
1590  }
1591  else
1592  {
1593  distrib_woTRDStruc ->Fill(inhisto->GetBinContent(cell+1));
1594  }
1595  }
1596  }
1597 
1598 
1599  TCanvas *c1 = new TCanvas(histoName,histoName,900,900);
1600  c1->ToggleEventStatus();
1601  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
1602  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
1603  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
1604  upperPad->Draw();
1605  lowerPadLeft->Draw();
1606  lowerPadRight->Draw();
1607 
1608  upperPad->cd();
1609  upperPad->SetLeftMargin(0.05);
1610  upperPad->SetRightMargin(0.05);
1611  upperPad->SetLogy();
1612  inhisto->SetTitleOffset(0.6,"Y");
1613  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
1614  inhisto->GetYaxis()->SetTitleOffset(0.7);
1615  inhisto->SetLineColor(kBlue+1);
1616  inhisto->Draw();
1617 
1618  //- - - - - - - - - - - - - - - - - - - -
1619  lowerPadRight->cd();
1620  lowerPadRight->SetLeftMargin(0.09);
1621  lowerPadRight->SetRightMargin(0.12);
1622  plot2D->GetYaxis()->SetTitleOffset(1.3);
1623  plot2D->Draw("colz");
1624 
1625  //- - - - - - - - - - - - - - - - - - - -
1626  lowerPadLeft->cd();
1627  lowerPadLeft->SetLeftMargin(0.09);
1628  lowerPadLeft->SetRightMargin(0.07);
1629  lowerPadLeft->SetLogy();
1630  distrib->SetLineColor(kBlue+1);
1631  distrib->GetYaxis()->SetTitleOffset(1.3);
1632  distrib->Draw("");
1633  distrib_wTRDStruc->SetLineColor(kGreen+1);
1634  distrib_wTRDStruc->DrawCopy("same");
1635  distrib_woTRDStruc->SetLineColor(kMagenta+1);
1636  distrib_woTRDStruc->DrawCopy("same");
1637 
1638  //- - - - - - - - - - - - - - - - - - - -
1639  //...fit time distr
1640  Double_t maxDistr = distrib->GetMaximum();
1641  Int_t maxBin = distrib->GetMaximumBin();
1642  Double_t maxCenter = distrib->GetBinCenter(maxBin);
1643  Double_t sig=0;
1644  Double_t mean=0;
1645 
1646  Double_t minFitRange=dminVal;
1647  Double_t maxFitRange=dmaxVal;
1648  //cout<<"max bin: "<<maxBin<<", max value: "<<maxDistr<<endl;
1649  //cout<<"start mean: "<<maxCenter<<", mean range: 0-"<<dmaxVal<<endl;
1650  //cout<<"fit range: "<<minFitRange<<" - "<<maxFitRange<<endl;
1651 
1652  TF1 *fitTime = new TF1("fitTime", "gaus",minFitRange,maxFitRange);
1653  //..start the fit with a mean of the highest value
1654  fitTime->SetParLimits(0,maxDistr*0.5,maxDistr*1.2); //..do not allow negative ampl values
1655  fitTime->SetParameter(1,maxCenter); //..set the start value for the mean
1656  fitTime->SetParLimits(1,dminVal,dmaxVal); //..do not allow negative mean values
1657 
1658  distrib->Fit(fitTime, "0LEM", "", minFitRange, maxFitRange);
1659  mean = fitTime->GetParameter(1);
1660  sig = fitTime->GetParameter(2);
1661 
1662  //cout<<"mean:"<<mean<<endl;
1663  //cout<<"sig: "<<sig<<endl;
1664 
1665  Double_t tCutMax=mean+nSig*sig;
1666  Double_t tCutMin=mean-nSig*sig;
1667 
1668 
1669  cout<<"Determined time cut: "<<tCutMin<<" - "<<tCutMax<<endl;
1670  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1671  //. . .Add info to histogram
1672  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1673  TLine *line = new TLine(tCutMax, 0, tCutMax, distrib->GetMaximum());
1674  line->SetLineColor(kGreen+2);
1675  line->SetLineStyle(7);
1676  line->Draw();
1677  TLine *line2 = new TLine(tCutMin, 0, tCutMin, distrib->GetMaximum());
1678  line2->SetLineColor(kGreen+2);
1679  line2->SetLineStyle(7);
1680  line2->Draw();
1681 
1682  fitTime->SetLineColor(kOrange-3);
1683  fitTime->SetLineStyle(1);//7
1684  fitTime->Draw("same");
1685 
1686  TLegend *leg = new TLegend(0.60,0.70,0.9,0.85);
1687  leg->AddEntry(line,"upper limit","l");
1688  leg->AddEntry(distrib_wTRDStruc,"Covered by TRD","l");
1689  leg->AddEntry(distrib_woTRDStruc,"wo TRD structure","l");
1690  leg->SetBorderSize(0);
1691  leg->Draw("same");
1692 
1693  TLatex* text = 0x0;
1694  text = new TLatex(0.12,0.85,Form("Good range: %.3f-%.3f",tCutMin,tCutMax));
1695  text->SetTextSize(0.06);
1696  text->SetNDC();
1697  text->SetTextColor(1);
1698  text->Draw();
1699 
1700  //- - - - -
1701  upperPad->cd();
1702  TLine *uline = new TLine(0, tCutMax,fNoOfCells,tCutMax);
1703  uline->SetLineColor(kGreen+2);
1704  uline->SetLineStyle(7);
1705  uline->Draw();
1706 
1707  TLine *lline = new TLine(0, tCutMin,fNoOfCells,tCutMin);
1708  lline->SetLineColor(kGreen+2);
1709  lline->SetLineStyle(7);
1710  lline->Draw();
1711  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1712  //. . .Save histogram
1713  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1714  c1->Update();
1715  TString name=Form("%s/%s/AverageTimeMax_%s.gif",fWorkdir.Data(), fAnalysisOutput.Data(), (const char*)histoName);
1716  c1->SaveAs(name);
1717 
1718  fRootFile->WriteObject(c1,c1->GetName());
1719  fRootFile->WriteObject(plot2D,plot2D->GetName());
1720  fRootFile->WriteObject(distrib,distrib->GetName());
1721  fRootFile->WriteObject(inhisto,inhisto->GetName());
1722  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1723  //. . . Mark the bad cells in the fFlag array
1724  //. . .(fCriterionCounter= bad because cell average value lower than min allowed)
1725  //. . .(fCriterionCounter= bad because cell average value higher than max allowed)
1726  //. . .(0 by default - good cell)
1727  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1728  if(fPrint==1)cout<<" o Flag bad cells that are outside the good range "<<endl;
1729  for(Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1730  {
1731  //..cell=0 and bin=1, cell=1 and bin=2
1732  //.. <= throws out zeros, might not be a dead cell but still have zero entries in a given energy range
1733  if(inhisto->GetBinContent(cell+1) > tCutMax && fFlag[cell]==0)
1734  {
1735  fFlag[cell]=fCriterionCounter;
1736  }
1737  }
1738  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;
1739 
1740  /*
1741  delete distrib;
1742  delete distrib_wTRDStruc;
1743  delete distrib_woTRDStruc;
1744  delete plot2D;
1745  delete c1;
1746  delete line;
1747  delete leg;
1748  delete text;
1749  delete uline;
1750  */
1751 }
1752 
1753 
1758 //________________________________________________________________________
1760 {
1761  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1762  //.. RESULTS
1763  //.. 1) Print the bad cells
1764  //.. and write the results to a file
1765  //.. for each added period analysis
1766  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1767  TArrayD periodArray;
1768  Double_t emin,emax,sig;
1769  Int_t criterion;
1770  TString output;
1771  Int_t nb1=0;
1772 
1773  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1774  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1775  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1776  TString aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1777  ofstream file2(aliceTwikiTable, ios::out | ios::trunc);
1778  if(file2)
1779  {
1780  file2<<"|*Criterium* | *E_min !GeV* | *E_max !GeV* | *sigma* | *Excluded Cells* |"<<endl;
1781  }
1782  file2.close();
1783 
1784  for(Int_t i=0;i<(Int_t)fAnalysisVector.size();i++)
1785  {
1786  periodArray=fAnalysisVector.at(i);
1787  criterion =periodArray.At(0);
1788  emin =periodArray.At(2);
1789  emax =periodArray.At(3);
1790  sig =periodArray.At(1);
1791 
1792  //..Print the results on the screen and
1793  //..write the results in a file
1794  output.Form("%s/%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt",fWorkdir.Data(), fAnalysisOutput.Data(), criterion,emin,emax);
1795  ofstream file(output, ios::out | ios::trunc);
1796  if(!file)
1797  {
1798  cout<<"#### Major Error. Check the textfile!"<<endl;
1799  }
1800  file<<"fFlag="<<i+2<<"means Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
1801  if(fPrint==1)cout<<" o Criterion : "<<criterion<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<" (Method "<<i<<")"<<endl;
1802 
1803  nb1=0;
1804  for(Int_t cellID=fStartCell ;cellID<fNoOfCells;cellID++)
1805  {
1806  if(fFlag[cellID]==(i+2))
1807  {
1808  nb1++;
1809  file<<cellID<<", ";
1810  }
1811  }
1812  file<<"Total number of bad cells with fFlag=="<<i+2<<endl;
1813  file<<"("<<nb1<<")"<<endl;
1814  file.close();
1815  if(fPrint==1)cout<<" o Total number of bad cells ("<<nb1<<")"<<endl;
1816  if(fPrint==1)cout<<endl;
1817  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1818  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
1819  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1820  ofstream file2(aliceTwikiTable, ios::out | ios::app);
1821  if(file2)
1822  {
1823  file2<<"| "<<criterion<<" | "<<emin<<" | "<<emax<<" | "<<sig<<" | "<<nb1<<" |"<<endl;
1824  }
1825  file2.close();
1826  }
1827 }
1834 //________________________________________________________________________
1836 {
1837  Int_t cellID, nDeadDCalCells = 0, nDeadEMCalCells = 0, nDCalCells = 0, nEMCalCells = 0;
1838  Double_t perDeadEMCal,perDeadDCal,perBadEMCal,perBadDCal,perWarmEMCal,perWarmDCal;
1839  TString aliceTwikiTable, cellSummaryFile, deadPdfName, badPdfName, ratioOfBad,goodCells,goodCellsRatio,cellProp;
1840  TString OADBFile_bad, OADBFile_dead, OADBFile_warm;
1841  TString OADBFile_badRbR, OADBFile_deadRbR, OADBFile_warmRbR;
1842  TH2F* cellAmp_masked = (TH2F*)fCellAmplitude->Clone("hcellAmp_masked");
1843  TH2F* cellTime_masked= (TH2F*)fCellTime->Clone("fCellTime");
1844 
1845  deadPdfName = Form("%s/%s/%s_Dead_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1846  badPdfName = Form("%s/%s/%s_Bad_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1847  ratioOfBad = Form("%s/%s/%s_Bad_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1848  goodCells = Form("%s/%s/%s_Good_Ampl_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1849  goodCellsRatio = Form("%s/%s/%s_Good_Ampl_Ratio_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1850  cellSummaryFile = Form("%s/%s/%s_%s_Bad_Ampl_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1851  OADBFile_bad = Form("%s/%s/%s_%s_OADBFile_Bad_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1852  OADBFile_dead = Form("%s/%s/%s_%s_OADBFile_Dead_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1853  OADBFile_warm = Form("%s/%s/%s_%s_OADBFile_Warm_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1854  if(fRunBRunMap)OADBFile_badRbR = Form("%s/%s/%s_%s_OADBFile_Bad_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1855  if(fRunBRunMap)OADBFile_deadRbR = Form("%s/%s/%s_%s_OADBFile_Dead_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1856  if(fRunBRunMap)OADBFile_warmRbR = Form("%s/%s/%s_%s_OADBFile_Warm_V%i.txt",fWorkdir.Data(), fAnalysisOutputRbR.Data(),fPeriod.Data(), fTrigger.Data() ,fTrial); ;
1857  aliceTwikiTable = Form("%s/%s/%s_TwikiTable_V%i.txt",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial); ;
1858  cellProp = Form("%s/%s/%s_CellProp_V%i.pdf",fWorkdir.Data(), fAnalysisOutput.Data(), fTrigger.Data() ,fTrial);
1859 
1860  cout<<" o Final results o "<<endl;
1861  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1862  //..Create a masked version of the Amp vs. ID and Time vs. ID histograms
1863  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1864  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
1865  {
1866  //..Direction of amplitude (Checks energies from 0-nBins GeV)
1867  for (Int_t amp = 1; amp <= fCellAmplitude->GetNbinsX(); amp++)
1868  {
1869  if(fFlag[cell]!=0)
1870  {
1871  //..cellID+1 = histogram bin
1872  cellAmp_masked->SetBinContent(amp,cell+1,0);
1873  }
1874  }
1875  //..Direction of time (Checks times from -275-975 GeV)
1876  for (Int_t time = 1; time <= fCellTime->GetNbinsX(); time++)
1877  {
1878  if(fFlag[cell]!=0)
1879  {
1880  //..cellID+1 = histogram bin
1881  cellTime_masked->SetBinContent(time,cell+1,0);
1882  }
1883  }
1884  }
1885  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1886  //..Plot some summary canvases
1887  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1888  TCanvas *c1 = new TCanvas("CellProp","I summary of cell properties",1000,1000);
1889  c1->ToggleEventStatus();
1890  c1->Divide(2,2);
1891  c1->cd(1)->SetLogz();
1892  //lowerPadRight->SetLeftMargin(0.09);
1893  //lowerPadRight->SetRightMargin(0.06);
1894  fCellAmplitude->SetXTitle("Cell Energy [GeV]");
1895  fCellAmplitude->SetYTitle("Abs. Cell Id");
1896  fCellAmplitude->GetYaxis()->SetTitleOffset(1.6);
1897  fCellAmplitude->Draw("colz");
1898  c1->cd(2)->SetLogz();
1899  fCellTime->SetXTitle("Cell Time [ns]");
1900  fCellTime->SetYTitle("Abs. Cell Id");
1901  fCellTime->GetYaxis()->SetTitleOffset(1.6);
1902  fCellTime->Draw("colz");
1903  c1->cd(3)->SetLogz();
1904  //lowerPadRight->SetLeftMargin(0.09);
1905  //lowerPadRight->SetRightMargin(0.06);
1906  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1907  cellAmp_masked->SetXTitle("Cell Energy [GeV]");
1908  cellAmp_masked->SetYTitle("Abs. Cell Id");
1909  cellAmp_masked->GetYaxis()->SetTitleOffset(1.6);
1910  cellAmp_masked->Draw("colz");
1911  c1->cd(4)->SetLogz();
1912  cellTime_masked->SetTitle("Masked Cell Time");
1913  cellTime_masked->SetXTitle("Cell Time [ns]");
1914  cellTime_masked->SetYTitle("Abs. Cell Id");
1915  cellTime_masked->GetYaxis()->SetTitleOffset(1.6);
1916  cellTime_masked->Draw("colz");
1917  c1->Update();
1918 
1919  TCanvas *c1_ratio = new TCanvas("CellPropRatio","II summary of cell properties ratio",1000,500);
1920  c1_ratio->ToggleEventStatus();
1921  c1_ratio->Divide(2);
1922  c1_ratio->cd(1)->SetLogz();
1923  cellAmp_masked->SetTitle("Masked Cell Amplitude");
1924  cellAmp_masked->GetZaxis()->SetRangeUser(0.0001,10e7);
1925  cellAmp_masked->Draw("colz");
1926  c1_ratio->cd(2)->SetLogz();
1927 
1928  TH1D *hRefDistr = BuildMeanFromGood();
1929  TH2F* ratio2DAmp =(TH2F*)cellAmp_masked->Clone("ratio2DAmp");
1930  TH2F* Sum2DIdeal =(TH2F*)cellAmp_masked->Clone("Sum2DIdeal");
1931  Sum2DIdeal->Reset();
1932  //..Create an ideal 2D energy distribution for division.
1933  //..Helps to identify whether there are some cells that still
1934  //..need to be masked by hand
1935  for(Int_t eBin=0;eBin<Sum2DIdeal->GetNbinsX();eBin++)
1936  {
1937  Double_t binVal=hRefDistr->GetBinContent(eBin+1);
1938  for(Int_t icell=0;icell<Sum2DIdeal->GetNbinsY();icell++)
1939  {
1940  Sum2DIdeal->SetBinContent(eBin+1,icell+1,binVal);
1941  }
1942  }
1943  ratio2DAmp->SetTitle("Ratio of cell Amplitude to mean cell ampl.");
1944  ratio2DAmp->Divide(Sum2DIdeal);
1945  ratio2DAmp->GetZaxis()->UnZoom();
1946  ratio2DAmp->DrawCopy("colz");
1947  c1_ratio->Update();
1948 
1949  TLatex* textSM = new TLatex(0.1,0.1,"*test*");
1950  textSM->SetTextSize(0.06);
1951  textSM->SetTextColor(1);
1952  textSM->SetNDC();
1953 
1954  TCanvas *c1_proj = new TCanvas("CellPropPProj","III summary of cell properties",1000,500);
1955  c1_proj->ToggleEventStatus();
1956  c1_proj->Divide(2);
1957  c1_proj->cd(1)->SetLogy();
1958  TH1D* projEnergyMask = cellAmp_masked->ProjectionX(Form("%sMask_Proj",cellAmp_masked->GetName()),fStartCell,fNoOfCells);
1959  projEnergyMask->SetXTitle("Cell Energy [GeV]");
1960  projEnergyMask->GetYaxis()->SetTitleOffset(1.6);
1961  projEnergyMask->SetLineColor(kGreen+1);
1962  projEnergyMask->DrawCopy(" hist");
1963  TH1D* projEnergy = fCellAmplitude->ProjectionX(Form("%s_Proj",fCellAmplitude->GetName()),fStartCell,fNoOfCells);
1964  projEnergy->DrawCopy("same hist");
1965  TLegend *leg = new TLegend(0.50,0.75,0.7,0.87);
1966  leg->AddEntry(projEnergy,"all cells","l");
1967  leg->AddEntry(projEnergyMask,"good cells","l");
1968  leg->SetTextSize(0.05);
1969  leg->SetBorderSize(0);
1970  leg->SetFillColorAlpha(10, 0);
1971  leg->Draw("same");
1972  TLegend *legBig = (TLegend*)leg->Clone("legBig");
1973  legBig->SetTextSize(0.08);
1974  legBig->SetX1NDC(0.2);
1975 
1976  c1_proj->cd(2)->SetLogy();
1977  TH1* projTimeMask = cellTime_masked->ProjectionX(Form("%s_Proj",cellTime_masked->GetName()),fStartCell,fNoOfCells);
1978  projTimeMask->SetXTitle("Cell Time [ns]");
1979  projTimeMask->GetYaxis()->SetTitleOffset(1.6);
1980  projTimeMask->GetYaxis()->SetRangeUser(1,projTimeMask->GetMaximum()*20);
1981  projTimeMask->SetLineColor(kGreen+1);
1982  projTimeMask->DrawCopy("hist");
1983  TH1* projTime = fCellTime->ProjectionX(Form("%s_Proj",fCellTime->GetName()),fStartCell,fNoOfCells);
1984  projTime->DrawCopy("same hist");
1985  leg->Draw("same");
1986  c1_proj->Update();
1987 
1988  TCanvas *c1_projSM = new TCanvas("CellPropPProjSM","III summary of cell Energy per SM",1200,900);
1989  c1_projSM->Divide(5,4,0.001,0.001);
1990  TH1* projEnergyMaskSM[20];
1991  TH1* projEnergySM[20];
1992  for(Int_t iSM=0;iSM<20;iSM++)
1993  {
1994  c1_projSM->cd(iSM+1)->SetLogy();
1995  gPad->SetTopMargin(0.03);
1996  gPad->SetBottomMargin(0.11);
1997  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
1998  projEnergyMaskSM[iSM]->SetTitle("");
1999  projEnergyMaskSM[iSM]->SetXTitle(Form("Cell Energy [GeV], SM%i",iSM));
2000  projEnergyMaskSM[iSM]->GetYaxis()->SetTitleOffset(1.6);
2001  projEnergyMaskSM[iSM]->GetYaxis()->SetLabelSize(0.06);
2002  projEnergyMaskSM[iSM]->GetXaxis()->SetLabelSize(0.06);
2003  projEnergyMaskSM[iSM]->GetXaxis()->SetRangeUser(0,20);
2004  projEnergyMaskSM[iSM]->GetXaxis()->SetTitleSize(0.06);
2005  projEnergyMaskSM[iSM]->SetLineColor(kGreen+1);
2006  projEnergyMaskSM[iSM]->DrawCopy(" hist");
2007 
2008  projEnergySM[iSM] = fCellAmplitude->ProjectionX(Form("%s_ProjSM%i",fCellAmplitude->GetName(),iSM),fStartCellSM[iSM],fStartCellSM[iSM+1]-1);
2009  projEnergySM[iSM]->DrawCopy("same hist");
2010  if(iSM==0)legBig->Draw("same");
2011  //textSM->Draw();
2012  textSM->SetTitle(Form("Includes cell IDs %d-%d",fStartCellSM[iSM],fStartCellSM[iSM+1]-1));
2013  textSM->DrawLatex(0.2,0.9,Form("Includes cell IDs %d-%d",fStartCellSM[iSM],fStartCellSM[iSM+1]-1));
2014  }
2015  c1_projSM->Update();
2016 
2017  TCanvas *c1_projRSM = new TCanvas("CellPropPProjRSM","III summary of cell Energy Ratio per SM",1200,900);
2018  c1_projRSM->Divide(5,4,0.001,0.001);
2019  for(Int_t iSM=0;iSM<20;iSM++)
2020  {
2021  c1_projRSM->cd(iSM+1)->SetLogy();
2022  gPad->SetTopMargin(0.03);
2023  gPad->SetBottomMargin(0.11);
2024  //projEnergyMaskSM[iSM]->GetXaxis()->SetRangeUser(0,10);
2025  projEnergyMaskSM[iSM]->SetLineColor(kGray+1);
2026  projEnergyMaskSM[iSM]->Divide(hRefDistr);
2027  projEnergyMaskSM[iSM]->DrawCopy("hist");
2028  }
2029  c1_projRSM->Update();
2030 
2031  TCanvas *c1_projTimeSM = new TCanvas("CellPropPProjTimeSM","III summary of cell Time per SM",1200,900);
2032  c1_projTimeSM->Divide(5,4,0.001,0.001);
2033  TH1* projTimeMaskSM[20];
2034  TH1* projTimeSM[20];
2035  for(Int_t iSM=0;iSM<20;iSM++)
2036  {
2037  c1_projTimeSM->cd(iSM+1)->SetLogy();
2038  gPad->SetTopMargin(0.03);
2039  gPad->SetBottomMargin(0.11);
2040  projTimeMaskSM[iSM] = cellTime_masked->ProjectionX(Form("%sMask_ProjSMTime%i",cellAmp_masked->GetName(),iSM),fStartCellSM[iSM]+1,fStartCellSM[iSM+1]);
2041  projTimeMaskSM[iSM]->SetTitle("");
2042  projTimeMaskSM[iSM]->SetXTitle(Form("Cell Time [ns], SM%i",iSM));
2043  projTimeMaskSM[iSM]->GetYaxis()->SetTitleOffset(1.6);
2044  projTimeMaskSM[iSM]->GetYaxis()->SetLabelSize(0.06);
2045  projTimeMaskSM[iSM]->GetXaxis()->SetLabelSize(0.06);
2046  //projTimeMaskSM[iSM]->GetXaxis()->SetRangeUser(0,20);
2047  projTimeMaskSM[iSM]->GetXaxis()->SetTitleSize(0.06);
2048  projTimeMaskSM[iSM]->SetLineColor(kGreen+1);
2049  projTimeMaskSM[iSM]->DrawCopy(" hist");
2050 
2051  if(iSM==0)legBig->Draw("same");
2052  projTimeSM[iSM] = fCellTime->ProjectionX(Form("%s_ProjSMTime%i",fCellAmplitude->GetName(),iSM),fStartCellSM[iSM],fStartCellSM[iSM+1]-1);
2053  projTimeSM[iSM]->DrawCopy("same hist");
2054  }
2055 
2056  /*
2057  //..This part here eats up too much memory so it was commented out
2058  //..The Canvases are anyway later saved as individual gifs.
2059  //..save to a PDF
2060  c1 ->Print(Form("%s(",cellProp.Data()));
2061  c1_ratio ->Print(Form("%s",cellProp.Data()));
2062  c1_proj ->Print(Form("%s",cellProp.Data()));
2063  c1_projSM ->Print(Form("%s",cellProp.Data()));
2064  c1_projRSM ->Print(Form("%s",cellProp.Data()));
2065  c1_projTimeSM->Print(Form("%s)",cellProp.Data()));
2066  */
2067  //..Scale the histograms by the number of events
2068  //..so that they are more comparable for a run-by-run
2069  //..analysis
2070  Double_t totalevents = fProcessedEvents->Integral();
2071  fCellAmplitude ->Scale(1.0/totalevents);
2072  cellAmp_masked ->Scale(1.0/totalevents);
2073  fCellTime ->Scale(1.0/totalevents);
2074 
2075  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2076  //..Write the final results of dead and bad cells in a file and on screen
2077  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2078  ofstream file(cellSummaryFile, ios::out | ios::trunc);
2079  ofstream fileBad(OADBFile_bad, ios::out | ios::trunc);
2080  ofstream fileDead(OADBFile_dead, ios::out | ios::trunc);
2081  ofstream fileWarm(OADBFile_warm, ios::out | ios::trunc);
2082  /*ofstream fileBadRbR;
2083  ofstream fileDeadRbR;
2084  ofstream fileWarmRbR;
2085  if(fRunBRunMap)
2086  {*/
2087  ofstream fileBadRbR(OADBFile_badRbR, ios::out | ios::trunc);
2088  ofstream fileDeadRbR(OADBFile_deadRbR, ios::out | ios::trunc);
2089  ofstream fileWarmRbR(OADBFile_warmRbR, ios::out | ios::trunc);
2090  //}
2091  if(file)
2092  {
2093  file<<"Dead cells : "<<endl;
2094  cout<<" o Dead cells : "<<endl;
2095  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
2096  {
2097  if(cellID==0)
2098  {
2099  file<<"In EMCal: "<<endl;
2100  }
2101  if(cellID==fCellStartDCal)
2102  {
2103  file<<"\n"<<endl;
2104  file<<"In DCal: "<<endl;
2105  }
2106  if(fFlag[cellID]==1)
2107  {
2108  file<<cellID<<", ";
2109  if(cellID<fCellStartDCal)nDeadEMCalCells++;
2110  else nDeadDCalCells++;
2111  }
2112  if(fFlag[cellID]==1)fileDead<<cellID<<endl;
2113  if(fFlag[cellID]>1)fileBad<<cellID<<endl;
2114  if(fRunBRunMap==1 && fFlag[cellID]==1)fileDeadRbR<<cellID<<endl;
2115  if(fRunBRunMap==1 && fFlag[cellID]>1)fileBadRbR<<cellID<<endl;
2116  }
2117  file<<"\n"<<endl;
2118  perDeadEMCal=100*nDeadEMCalCells/(1.0*fCellStartDCal);
2119  perDeadDCal=100*nDeadDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2120  file<<"EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
2121  cout<<" o EMCal ("<<nDeadEMCalCells<<" ="<<perDeadEMCal<<"%), DCal ("<<nDeadDCalCells<<" ="<<perDeadDCal<<"%)"<<endl;
2122 
2123  file<<"Bad cells: "<<endl;
2124  cout<<" o Bad cells: "<<endl;
2125  for(cellID=fStartCell;cellID<fNoOfCells;cellID++)
2126  {
2127  if(cellID==0)
2128  {
2129  file<<"In EMCal: "<<endl;
2130  }
2131  if(cellID==fCellStartDCal)
2132  {
2133  file<<"\n"<<endl;
2134  file<<"In DCal: "<<endl;
2135  }
2136  if(fFlag[cellID]>1)
2137  {
2138  file<<cellID<<", ";
2139  if(cellID<fCellStartDCal)nEMCalCells++;
2140  else nDCalCells++;
2141  }
2142  }
2143  file<<"\n"<<endl;
2144  perBadEMCal=100*nEMCalCells/(1.0*fCellStartDCal);
2145  perBadDCal =100*nDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2146  file<<"EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
2147  cout<<" o EMCal ("<<nEMCalCells<<" ="<<perBadEMCal<<"%), DCal ("<<nDCalCells<<" ="<<perBadDCal<<"%)"<<endl;
2148  }
2149  file.close();
2150  cout<<" o Total: "<<endl;
2151  cout<<" o Bad+Dead cells: "<<nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells<<", this is "<<(nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells)*100/(1.0*fNoOfCells)<<"% of the whole detector"<<endl;
2152 
2153  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2154  //..Determine the number of warm cells and save the flagged cells to .pdf files
2155  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2156  if(fPrint==1)cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
2157  SaveBadCellsToPDF(1,badPdfName) ;
2158  SaveBadCellsToPDF(10,ratioOfBad) ; //..Special case
2159  if(fTestRoutine==1)SaveBadCellsToPDF(2,goodCells) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
2160  if(fTestRoutine==1)SaveBadCellsToPDF(20,goodCellsRatio) ; //..Special case all good cells to check, should all have a flag naming them *Candidate*
2161 
2162  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2163  //..Fill the histograms with the flag information
2164  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2165  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2166  {
2167  fhCellFlag->SetBinContent(cell+1,fFlag[cell]);
2168  fhCellWarm->SetBinContent(cell+1,fWarmCell[cell]);
2169  if(fWarmCell[cell]==1)fileWarm<<cell<<endl;
2170  if(fRunBRunMap==1 && fWarmCell[cell]==1)fileWarmRbR<<cell<<endl;
2171  }
2172  TCanvas *c2 = new TCanvas("CellFlag","summary of cell flags",1200,800);
2173  c2->ToggleEventStatus();
2174  c2->Divide(1,2);
2175  c2->cd(1);
2176  fhCellFlag->SetTitle("cell flag");
2177  fhCellFlag->SetXTitle("Abs. Cell Id");
2178  fhCellFlag->SetYTitle("flag by which cell was excluded");
2179  fhCellFlag->DrawCopy("hist");
2180  c2->cd(2);
2181  fhCellWarm->SetTitle("Warm cells");
2182  fhCellWarm->SetXTitle("Abs. Cell Id");
2183  fhCellWarm->SetYTitle("warm=1");
2184  fhCellWarm->DrawCopy("hist");
2185  c2->Update();
2186 
2187  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2188  //..Plot the 2D distribution of cells by flag
2189  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2190  PlotFlaggedCells2D(0); //..all good cells
2191  PlotFlaggedCells2D(1); //..all dead cells
2192  PlotFlaggedCells2D(2,fCriterionCounter); //..all bad cells
2193  PlotFlaggedCells2D(0,0); //..Special case - Warm cells
2194 
2195  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2196  //..Add different histograms/canvases to the output root file
2197  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2198  TString name1,name2,name3,name4,name5,name6;
2199  name1 = Form("%s/%s/CellPropertiesRatio.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2200  c1_ratio->SaveAs(name1);
2201  name2 = Form("%s/%s/CellProperties.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2202  c1->SaveAs(name2);
2203  name3 = Form("%s/%s/CellPropertiesProj.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2204  c1_proj->SaveAs(name3);
2205  name4 = Form("%s/%s/CellEnergySM.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2206  c1_projSM->SaveAs(name4);
2207  name5 = Form("%s/%s/CellEnergySMratio.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2208  c1_projRSM->SaveAs(name5);
2209  name6 = Form("%s/%s/CellTime.gif", fWorkdir.Data(),fAnalysisOutput.Data());
2210  c1_projTimeSM->SaveAs(name6);
2211 
2212  fRootFile->WriteObject(c1_ratio,c1_ratio->GetName());
2213  fRootFile->WriteObject(c1,c1->GetName());
2214  fRootFile->WriteObject(c1_proj,c1_proj->GetName());
2215  fRootFile->WriteObject(c1_projSM,c1_projSM->GetName());
2216  fRootFile->WriteObject(c1_projRSM,c1_projRSM->GetName());
2217  fRootFile->WriteObject(c1_projTimeSM,c1_projTimeSM->GetName());
2218 
2219  fRootFile->WriteObject(c2,c2->GetName());
2220  fRootFile->WriteObject(fCellAmplitude,fCellAmplitude->GetName());
2221  fRootFile->WriteObject(cellAmp_masked,cellAmp_masked->GetName());
2222  fRootFile->WriteObject(ratio2DAmp,ratio2DAmp->GetName());
2223  fRootFile->WriteObject(fCellTime,fCellTime->GetName());
2224  fRootFile->WriteObject(fProcessedEvents,fProcessedEvents->GetName());
2225  fRootFile->WriteObject(fhCellFlag,fhCellFlag->GetName());
2226  fRootFile->WriteObject(fhCellWarm,fhCellWarm->GetName());
2227  fRootFile->WriteObject(projEnergyMask,projEnergyMask->GetName());
2228  fRootFile->WriteObject(projEnergy,projEnergy->GetName());
2229  //..Save all amplitudes to the root file
2230  SaveHistoToFile();
2231 
2232  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2233  //..Save also the identified warm channels into a text file.
2234  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2235  Int_t nWEMCalCells =0;
2236  Int_t nWDCalCells =0;
2237  file.open(cellSummaryFile, ios::out | ios::app);
2238  if(file)
2239  {
2240  file<<"Warm cells : "<<endl;
2241  if(fPrint==1)cout<<" o Warm cells : "<<endl;
2242  for(cellID=fStartCell; cellID<fNoOfCells; cellID++)
2243  {
2244  if(cellID==0)
2245  {
2246  file<<"In EMCal: "<<endl;
2247  }
2248  if(cellID==fCellStartDCal)
2249  {
2250  file<<"\n"<<endl;
2251  file<<"In DCal: "<<endl;
2252  }
2253  if(fWarmCell[cellID]==1)
2254  {
2255  file<<cellID<<", ";
2256  if(cellID<fCellStartDCal)nWEMCalCells++;
2257  else nWDCalCells++;
2258  }
2259  }
2260  file<<"\n"<<endl;
2261  perWarmEMCal= 100*nWEMCalCells/(1.0*fCellStartDCal);
2262  perWarmDCal = 100*nWDCalCells/(1.0*fNoOfCells-fCellStartDCal);
2263  file<<"EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
2264  if(fPrint==1)cout<<" o EMCal ("<<nWEMCalCells<<" ="<<perWarmEMCal<<"%), DCal ("<<nWDCalCells<<" ="<<perWarmDCal<<"%)"<<endl;
2265  }
2266  file.close();
2267  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2268  //..Save the results in a tWiki format for the webpage (https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalQABadChannels2)
2269  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2270  ofstream file2(aliceTwikiTable, ios::out | ios::app);
2271  if(file2)
2272  {
2273  file2<<"1=energy/hit, 2= hit/event"<<endl;
2274  file2<<"\n"<<endl;
2275  file2<<"| *Detector* | *No of cells* | *percentage* |"<<endl;
2276  file2<<"| Dead EMCal | "<<nDeadEMCalCells<<" | "<<perDeadEMCal<<"% |"<<endl;
2277  file2<<"| Bad EMCal | "<<nEMCalCells<<" | "<<perBadEMCal<<"% |"<<endl;
2278  file2<<"| - Warm EMCal | "<<nWEMCalCells<<" | "<<perWarmEMCal<<"% |"<<endl;
2279  file2<<"| Dead DCal | "<<nDeadDCalCells<<" | "<<perDeadDCal<<"% |"<<endl;
2280  file2<<"| Bad DCal | "<<nDCalCells<<" | "<<perBadDCal<<"% |"<<endl;
2281  file2<<"| - Warm DCal | "<<nWDCalCells<<" | "<<perWarmDCal<<"% |"<<endl;
2282  file2<<"| Summ D+B | "<<nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells<<" | "<<(nDeadEMCalCells+nEMCalCells+nDeadDCalCells+nDCalCells)*100/(1.0*fNoOfCells)<<"% |"<<endl;
2283  file2<<"\n"<<endl;
2284  }
2285  file2.close();
2286 
2287  /*
2288  if(c1) delete c1;
2289  if(c1_ratio) delete c1_ratio;
2290  if(c1_projSM) delete c1_projSM;
2291  if(c1_projRSM) delete c1_projRSM;
2292  if(c1_projTimeSM) delete c1_projTimeSM;
2293  if(c2) delete c2;
2294  if(textSM) delete textSM;
2295  if(c1_proj) delete c1_proj;
2296  if(leg) delete leg;
2297  if(legBig) delete legBig;
2298  */
2299 }
2300 
2316 //________________________________________________________________________
2318 {
2319  gROOT->SetStyle("Plain");
2320  gStyle->SetOptStat(0);
2321  gStyle->SetFillColor(kWhite);
2322  gStyle->SetTitleFillColor(kWhite);
2323  gStyle->SetPalette(1);
2324 
2325  char title[100];
2326  char name[100];
2327 
2328  TH1D *hRefDistr = BuildMeanFromGood();
2329  fRootFile->WriteObject(hRefDistr,hRefDistr->GetName());
2330  Int_t firstCanvas=0;
2331  Bool_t candidate;
2332  TLatex* text = new TLatex(0.2,0.8,"*Candidate*");
2333  text->SetTextSize(0.07);
2334  text->SetTextColor(kOrange-3);
2335  text->SetNDC();
2336 
2337  TLatex* text2 = new TLatex(0.2,0.8,"*Not a Candidate*");
2338  text2->SetTextSize(0.08);
2339  text2->SetTextColor(8);
2340  text2->SetNDC();
2341 
2342  TLatex* textA = new TLatex(0.65,0.62,"*test*");
2343  textA->SetTextSize(0.04);
2344  textA->SetTextColor(1);
2345  textA->SetNDC();
2346 
2347  //..collect cells in an internal vector.
2348  //..when the vector is of size=9 or at the end of being filled
2349  //..plot the channels into a canvas
2350  std::vector<Int_t> channelVector;
2351  channelVector.clear();
2352  cout<<" o Start printing into .pdf for version: "<<version<<endl;
2353  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
2354  {
2355  if(fFlag[cell]==1 && version==0)channelVector.push_back(cell);
2356  if(fFlag[cell]>1 && version==1)channelVector.push_back(cell);
2357  if(fFlag[cell]==0 && version==2)channelVector.push_back(cell);
2358  if(fFlag[cell]>1 && version==10)channelVector.push_back(cell);
2359  if(fFlag[cell]==0 && version==20)channelVector.push_back(cell);
2360 
2361  if(cell%2000==0)cout<<"...cell No."<<cell<<endl;
2362  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
2363  if(channelVector.size()==9 || cell == fNoOfCells-1)
2364  {
2365  cout<<"."<<flush;
2366 
2367  TString internal_pdfName=pdfName;
2368  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
2369  if(channelVector.size() > 6) c1->Divide(3,3);
2370  else if (channelVector.size() > 3) c1->Divide(3,2);
2371  else c1->Divide(3,1);
2372 
2373  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
2374  for(Int_t i=0; i< (Int_t)channelVector.size() ; i++)
2375  {
2376  if(channelVector.size() >=fNoOfCells) cout<<"Massive problem"<<endl;
2377  sprintf(name, "Cell %d",channelVector.at(i)) ;
2378  TH1 *hCell = fCellAmplitude->ProjectionX(name,channelVector.at(i)+1,channelVector.at(i)+1);
2379  sprintf(title,"Cell No: %d Entries: %d",channelVector.at(i), (Int_t)hCell->GetEntries()) ;
2380 
2381  c1->cd(i%9 + 1);
2382  c1->cd(i%9 + 1)->SetLogy();
2383  if(i==0)
2384  {
2385  leg->AddEntry(hRefDistr,"mean of good","l");
2386  leg->AddEntry(hCell,"current","l");
2387  }
2388  //..Check the distribution whether it looks like a *Candidate* for a miscalibrated warm cell
2389  candidate = CheckDistribution(hCell,hRefDistr);
2390  if(candidate==1)fWarmCell[channelVector.at(i)]=candidate;
2391  if(version>2)//..These are ratio plots of energy distr. of cell and mean of all good cells
2392  {
2393  hCell->Divide(hRefDistr);
2394  }
2395  //.. save histograms to file if you want to double check the output
2396  if(fTestRoutine ==1)
2397  {
2398  if(version==1) fOutputListBad->Add(hCell);
2399  if(version==10)fOutputListBadRatio->Add(hCell);
2400  if(version==2) fOutputListGood->Add(hCell);
2401  if(version==20)fOutputListGoodRatio->Add(hCell);
2402  }
2403  hCell->SetLineColor(kBlue+1);
2404  hCell->GetXaxis()->SetTitle("E (GeV)");
2405  hCell->GetYaxis()->SetTitle("N Entries");
2406  hCell->GetXaxis()->SetRangeUser(0.,10.);
2407  hCell->SetLineWidth(1) ;
2408  hCell->SetTitle(title);
2409  hRefDistr->SetLineColor(kGray+2);
2410  hRefDistr->SetLineWidth(1);
2411 
2412  hCell->Draw("hist");
2413 
2414  if(version==1 || version==2)hRefDistr->Draw("same hist") ;
2415 
2416  //..Mark the histogram that could be miscalibrated and labelled as warm
2417  if(candidate==1 && (version==1 || version==10))
2418  {
2419  gPad->SetFrameFillColor(kYellow-10);
2420  text->Draw();
2421  }
2422  if(version==1)
2423  {
2424  textA->SetTitle(Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
2425  textA->DrawLatex(0.65,0.62,Form("Excluded by No. %d",fFlag[channelVector.at(i)]));
2426  }
2427  if(candidate==0 && (version==2 || version==20))
2428  {
2429  gPad->SetFrameFillColor(kYellow-10);
2430  text2->Draw(); //..Draw a marker in the histogram that was falsley missed as a good candidate
2431  leg->Draw();
2432  }
2433  if(version<2)leg->Draw();
2434  }
2435 
2436  if(channelVector.size()<9 || cell == fNoOfCells-1)
2437  {
2438  internal_pdfName +=")";
2439  c1->Print(internal_pdfName.Data());
2440  }
2441  else if(firstCanvas==0)
2442  {
2443  internal_pdfName +="(";
2444  c1->Print(internal_pdfName.Data());
2445  firstCanvas=1;
2446  }
2447  else
2448  {
2449  c1->Print(internal_pdfName.Data());
2450  }
2451  delete c1;
2452  delete leg;
2453  channelVector.clear();
2454  }
2455  }
2456  cout<<endl;
2457  delete hRefDistr;
2458  delete text;
2459  delete text2;
2460  delete textA;
2461  //..Add the subdirectories to the file
2462  if(version==1) fRootFile->WriteObject(fOutputListBad,fOutputListBad->GetName());
2463  if(version==10)fRootFile->WriteObject(fOutputListBadRatio,fOutputListBadRatio->GetName());
2464  if(version==2) fRootFile->WriteObject(fOutputListGoodRatio,fOutputListGoodRatio->GetName());
2465  if(version==20)fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
2466 
2467  if(fPrint==1)cout<<endl;
2468 }
2472 //_________________________________________________________________________
2474 {
2475  TH1D* hGoodAmp;
2476  TH1D* hgoodMean = (TH1D*)fCellAmplitude->ProjectionX("hgoodMean");
2477  hgoodMean->Reset();
2478  Int_t NrGood=0;
2479 
2480  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2481  {
2482  if(warmIn==0 && fFlag[cell]!=0 )continue;
2483  if(warmIn==1 && fFlag[cell]!=0 && fWarmCell[cell]==0)continue;
2484  if(warmIn==2 && fWarmCell[cell]==0)continue;
2485  NrGood++;
2486 
2487  hGoodAmp = (TH1D*)fCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
2488  hgoodMean->Add(hGoodAmp);
2489  }
2490  hgoodMean->Scale(1.0/NrGood);
2491 
2492  return hgoodMean;
2493 }
2504 //_________________________________________________________________________
2506 {
2507  TString Name = Form("%s_ratio",histogram->GetName());
2508  TH1* ratio=(TH1*)histogram->Clone(Name);
2509  ratio->Divide(reference);
2510 
2511  Double_t percentageOfLast=0.6;
2512  Double_t higherThanMean=5;
2513  Double_t highestRatio=1000;
2514  Double_t maxMean=10;
2515  Double_t minMean=0.1;
2516  Double_t cliffSize=2; //height before cliff shouldn't be double as high than after cliff
2517 
2518  //..By default each cell is a candidate for recalibration
2519  Bool_t candidate=1;
2520  //..Find bin where reference has value 1, and the corresponding x-value
2521  Double_t totalevents = fProcessedEvents->Integral();
2522  Int_t binHeightOne = reference->FindLastBinAbove(1.0/totalevents);//..look at the spectrum until there is 1hit/event GeV
2523  Double_t binCentreHeightOne = reference->GetBinCenter(binHeightOne);
2524  Double_t thirdBinCentre = reference->GetBinCenter(3);
2525 
2526  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2527  //..Check the histogram
2528  //..Different checks to see whether the
2529  //..cell is really bad. Set candidate to 0.
2530 
2531  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2532  //..check end of spectrum, should be larger than "percentageOfLast"% of the end of the mean histogram
2533  if(ratio->FindLastBinAbove(0)<ratio->FindBin(binCentreHeightOne*percentageOfLast))
2534  {
2535  candidate=0;
2536  //cout<<"1"<<endl;
2537  }
2538 
2539  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2540  //..Check maximum of ratio. Cell should not have "highestRatio" times more entries than reference in any bin
2541  //ELI - TO DO: check that crieteria carfully - seems to work but not shure about it
2542  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,10);//..zoom in to find the maximum between "not first 2 bins" - 10 GeV
2543  if(ratio->GetMaximum()>highestRatio)//
2544  {
2545  candidate=0;
2546  //cout<<"2"<<endl;
2547  }
2548 
2549  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2550  //..check whether the ratio is much larger than 1
2551  //..calculate the mean in the relevant energy range
2552  Double_t mean=0;
2553  Int_t nullEntries=0;
2554  for(Int_t i=2;i<binHeightOne;i++)
2555  {
2556  if(ratio->GetBinContent(i)!=0)mean+=ratio->GetBinContent(i);
2557  else nullEntries++;
2558  }
2559  mean*=1.0/(binHeightOne-1-nullEntries);//..divide by number of bins (excluding bins without entries)
2560  if(mean>maxMean || mean<minMean)
2561  {
2562  candidate=0;
2563  //cout<<"3"<<endl;
2564  }
2565 
2566  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2567  //..check whether there are large spikes in the histogram
2568  //..compare bin values to mean of the ratio. If there is a bin value with
2569  //..content "higherThanMean" times lareger than mean it's losing it candidate status
2570  mean=0;
2571  //..Find the maximum in the mean range (0-binHeightOne)
2572  ratio->GetXaxis()->SetRangeUser(0,binCentreHeightOne);
2573  Double_t localMaxBin=ratio->GetMaximumBin();
2574 
2575  for(Int_t i=2;i<binHeightOne;i++)
2576  {
2577  //..Exclude 0 bins and exclude bins near the maximum
2578  if(ratio->GetBinContent(i)<=0) continue;
2579  if(i>localMaxBin-3 && i<localMaxBin+3)continue;
2580  mean+=ratio->GetBinContent(i);
2581  }
2582  mean*=1.0/(binHeightOne-1);//..divide by number of bins
2583  ratio->GetXaxis()->SetRangeUser(thirdBinCentre,binCentreHeightOne);//..zoom in to find the maximum between 0-BinOne
2584  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
2585  if(ratio->GetMaximum()>mean*higherThanMean)
2586  {
2587  candidate=0;
2588  //cout<<"4"<<endl;
2589  }
2590 
2591  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2592  //..Check for a cliff at 4GeV, happens in some cases
2593  Double_t beforeCliff=0;
2594  Double_t afterCliff=0;
2595  Int_t binsBefore=0;
2596  Int_t binsAfter=0;
2597  Int_t cliffBin = ratio->FindBin(4);
2598  for(Int_t i=cliffBin-10;i<cliffBin+11;i++)
2599  {
2600  //..Exclude 0 bins
2601  if(ratio->GetBinContent(i)<=0)continue;
2602  if(i<=cliffBin) binsBefore++;
2603  if(i>cliffBin) binsAfter++;
2604  if(i<=cliffBin) beforeCliff+=ratio->GetBinContent(i);
2605  if(i>cliffBin) afterCliff+=ratio->GetBinContent(i);
2606  beforeCliff*=1.0/binsBefore;
2607  afterCliff*=1.0/binsAfter;
2608  }
2609  if((beforeCliff-afterCliff)>cliffSize*afterCliff)
2610  {
2611  if(beforeCliff!=0 && afterCliff!=0)candidate=0;
2612  if(beforeCliff!=0 && afterCliff!=0)cout<<"5"<<endl;
2613  }
2614 
2615  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2616  //..Employ peak finder
2617 /* if(candidate==1)
2618  {
2619  TSpectrum *spec = new TSpectrum(2,1); //..Nr peaks, resol. 1=3sigma distance
2620  Int_t nfound = spec->Search(ratio,4.3,"nobackground new",0.70);
2621  //GetNPeaks();
2622  //cout<<"found N peaks: "<<nfound<<endl;
2623  }
2624 */
2625  return candidate;
2626 }
2627 
2636 //_________________________________________________________________________
2638 {
2639  //..TRD support structure
2640  //..(determined by eye, could be improved, but is already very acurate):
2641  //..collumn 4,5,6,7,8 33,34,35,36 58,59,60 85,86,87,88,89
2642  //..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
2643  Bool_t coveredByTRDSupportStruc=0;
2644 
2645  if((collumn>3 && collumn<9) || (collumn>32 && collumn<37) || (collumn>57 && collumn<61) || (collumn>84 && collumn<90) ||
2646  (row==1) ||(row>20 && row<25) || (row>44 && row<49) || (row>68 && row<73) || (row>91 && row<96) ||
2647  (row>116 && row<120)|| row==127 || (row>148 && row<152) || (row>172 && row<177) || (row>197 && row<201)
2648  )
2649  {
2650  coveredByTRDSupportStruc=1;
2651  }
2652  return coveredByTRDSupportStruc;
2653 }
2662 //_________________________________________________________________________
2664 {
2665  gStyle->SetPalette(55); //kRainBow==55
2666  gStyle->SetPalette(91); //kPastel==91
2667  //..build two dimensional histogram with values row vs. column
2668  TString histoName;
2669  histoName = Form("2DChannelMap_Flag%d_V%i",flagBegin,fTrial);
2670  if(flagBegin==0 && flagEnd==0)histoName = Form("2DChannelMap_Flag100_V%i",fTrial);
2671 
2672  TH2F *plot2D = new TH2F(histoName,histoName,fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
2673  plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
2674  plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
2675 
2676  Int_t cellColumn=0,cellRow=0;
2677  Int_t cellColumnAbs=0,cellRowAbs=0;
2678  Int_t trash;
2679 
2680  for (Int_t cell = fStartCell; cell < fNoOfCells; cell++)
2681  {
2682  //..Do that only for cell ids also accepted by the code
2683  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
2684 
2685  //..Get Row and Collumn for cell ID c
2686  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
2687  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
2688  {
2689  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
2690  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
2691  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
2692  }
2693  if(flagEnd==-1 && fFlag[cell]==flagBegin) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
2694  if(flagEnd!=0 && flagEnd!=-1 && fFlag[cell]>=flagBegin && fFlag[cell]<=flagEnd)plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1);
2695  if(flagBegin==0 && flagEnd==0 && fWarmCell[cell]==1) plot2D->SetBinContent(cellColumnAbs,cellRowAbs,1); //warm cells
2696 
2697 
2698  }
2699  TCanvas *c1 = new TCanvas(histoName,histoName,500,500);
2700  c1->ToggleEventStatus();
2701  c1->cd();
2702  //lowerPadRight->SetLeftMargin(0.09);
2703  //lowerPadRight->SetRightMargin(0.06);
2704  plot2D->Draw("colz");
2705 
2706  TLatex* text = 0x0;
2707  if(flagBegin==0 && flagEnd==-1) text = new TLatex(0.2,0.8,"Good Cells");
2708  if(flagBegin==1) text = new TLatex(0.2,0.8,"Dead Cells");
2709  if(flagBegin>1) text = new TLatex(0.2,0.8,"Bad Cells");
2710  if(flagBegin==0 && flagEnd==0) text = new TLatex(0.2,0.8,"Warm Cells");
2711  text->SetTextSize(0.06);
2712  text->SetNDC();
2713  text->SetTextColor(1);
2714  text->Draw();
2715 
2716  c1->Update();
2717  TString name =Form("%s/%s/%s_%s.gif", fWorkdir.Data(),fAnalysisOutput.Data(),fPeriod.Data() , histoName.Data());
2718  c1->SaveAs(name);
2720  fRootFile->WriteObject(plot2D,plot2D->GetName());
2721 
2722  /*
2723  delete plot2D;
2724  delete c1;
2725  delete text;
2726  */
2727 }
2731 //_________________________________________________________________________
2733 {
2734  char name[100];
2735  for(Int_t cell=fStartCell;cell<fNoOfCells;cell++)
2736  {
2737  sprintf(name, "Cell %d",cell) ;
2738  TH1 *hCell = fCellAmplitude->ProjectionX(name,cell+1,cell+1);
2739  if(fFlag[cell]==0)fOutputListGood->Add(hCell);
2740  }
2741  fRootFile->WriteObject(fOutputListGood,fOutputListGood->GetName());
2742 }
2743 
2748 //_________________________________________________________________________
2750 {
2751  if(fTrackCellRecord==1)
2752  {
2753  Int_t zeroFlag = std::count(fFlag.begin(), fFlag.end(), 0);
2754  Int_t zeroWarm = std::count(fWarmCell.begin(), fWarmCell.end(), 0);
2755 
2756  cout<<"******* Debug No "<<number<<" ********"<<endl;
2757  cout<<"*** fCriterionCounter: "<<fCriterionCounter<<endl;
2758  cout<<"*** number of period analyses: "<<fAnalysisVector.size()<<endl;
2759  cout<<"*** number of man mask cells: "<<fManualMask.size()<<endl;
2760  cout<<"*** Bad+Dead(!0) fFlag elements: "<<fFlag.size()-zeroFlag<<endl;
2761  cout<<"*** Warm(!0) fWarmCell elements: "<<fWarmCell.size()-zeroWarm<<endl;
2762  cout<<"*** "<<endl;
2763  }
2764 }
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:97
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:98
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
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:96
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.
Definition: BadChannelAna.h:99
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