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