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