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