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