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