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