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