AliPhysics  9fe175b (9fe175b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaCaloChannelAnalysis.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TFile.h>
18 #include <TH1F.h>
19 #include <TH2F.h>
20 #include <TF1.h>
21 #include <TROOT.h>
22 #include <TLine.h>
23 #include <TCanvas.h>
24 #include <TStyle.h>
25 #include <TSystem.h>
26 #include <TLegend.h>
27 #include <TString.h>
28 #include <TLatex.h>
29 
30 #include <Riostream.h>
31 #include <stdio.h>
32 #include <fstream>
33 #include <iostream>
34 #include <alloca.h>
35 #include <string>
36 #include <cstring>
37 
38 // --- ANALYSIS system ---
39 #include <AliCalorimeterUtils.h>
40 #include <AliEMCALGeometry.h>
42 #include <AliAODEvent.h>
43 
44 using namespace std;
45 
49 
53 //________________________________________________________________________
55  TObject(),
56  fCurrentRunNumber(-1),fPeriod(),fPass(),fTrigger(),fNoOfCells(),
57  fMergeOutput(),fAnalysisOutput(),fAnalysisInput(),fRunList(),
58  fQADirect(), fMergedFileName(), fAnalysisVector(),
59  fRunListFileName(),fWorkdir(),fTrial(),fExternalFileName(),
60  fCaloUtils()
61 {
62  fCurrentRunNumber = 254381;
63  fPeriod = "LHC16h";
64  fPass = "muon_caloLego";
65  fTrigger = "AnyINT";
66 
67  Init();
68 }
69 
73 //________________________________________________________________________
74 AliAnaCaloChannelAnalysis::AliAnaCaloChannelAnalysis(TString period, TString pass, TString trigger, Int_t RunNumber):
75  TObject(),
76  fCurrentRunNumber(-1),fPeriod(),fPass(),fTrigger(),fNoOfCells(),
77  fMergeOutput(),fAnalysisOutput(),fAnalysisInput(),fRunList(),
78  fQADirect(), fMergedFileName(), fAnalysisVector(),
79  fRunListFileName(),fWorkdir(),fTrial(),fExternalFileName(),
80  fCaloUtils()
81 {
82  fCurrentRunNumber = RunNumber;
83  fPeriod = period;
84  fPass = pass;
85  fTrigger = trigger;
86 
87  Init();
88 }
89 
93 //________________________________________________________________________
95 {
96  //......................................................
97  //..Default values - can be set by functions
98  fWorkdir="./";
99  fRunListFileName="runList.txt";
100  fTrial = 0;
102 
103  //..Settings for the input/output structure (hard coded)
104  fAnalysisInput ="AnalysisInput";
105  fMergeOutput ="ConvertOutput";
106  fAnalysisOutput ="AnalysisOutput";
107  //..Stuff for the convert function
108  gSystem->mkdir(fMergeOutput);
109  gSystem->mkdir(fAnalysisOutput);
110  fMergedFileName= Form("%s/%s_%s_Merged.root", fMergeOutput.Data(), fPeriod.Data(),fPass.Data());
111  fQADirect = Form("CaloQA_%s",fTrigger.Data());
112  fRunList = Form("%s/%s/%s/%s", fAnalysisInput.Data(), fPeriod.Data(), fPass.Data(), fRunListFileName.Data());
113 
114  //.. make sure the vector is empty
115  fAnalysisVector.clear();
116 
117  //......................................................
118  //..Initialize EMCal/DCal geometry
120  //..Create a dummy event for the CaloUtils
121  AliAODEvent* aod = new AliAODEvent();
124  //fCaloUtils->AccessOADB(aod);//ELI do we need that. Works not anyway
125 
126  fNoOfCells =fCaloUtils->GetEMCALGeometry()->GetNCells(); //..Very important number, never change after that point!
127  Int_t NModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
128  fCellStartDCal=12288; //..ELI this should be automatized from the geometry information!!
129 
130  //..This is how the calorimeter looks like in the current period (defined by example run ID fCurrentRunNumber)
131  cout<<"Called geometry for run number: "<<fCurrentRunNumber<<endl;
132  cout<<"Number of supermod: "<<NModules<<endl;
133  cout<<"Number of cells: "<<fNoOfCells<<endl;
134  cout<<"Cell ID of first DCal cell: "<<fCellStartDCal<<endl;
135  //cout<<"Number of supermod utils: "<<fCaloUtils->GetNumberOfSuperModulesUsed()<<endl; //..will always be 22 unless set by hand
136 
137  //......................................................
138  //..Initialize flag array to store how the cell is categorized
139  //..In the histogram: bin 1= cellID 0, bin 2= cellID 1 etc
140  //..In the array: fFlag[cellID]= some information
141  fFlag = new Int_t[fNoOfCells];
142  fFlag[fNoOfCells]={0}; //..flagged as good by default
143 
144  //......................................................
145  //..setings for the 2D histogram
146  fNMaxCols = 48; //eta direction
147  fNMaxRows = 24; //phi direction
149  fNMaxRowsAbs = Int_t (NModules/2)*fNMaxRows; //multiply by number of supermodules
150 }
151 
160 //________________________________________________________________________
162 {
163  TString inputfile = "";
164 
165  if(fExternalFileName=="")
166  {
167  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
168  cout<<". . .Start process by converting files. . . . . . . . . . . ."<<endl;
169  cout<<endl;
171  cout<<endl;
172  }
173  else
174  {
175  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
176  cout<<". . .Start process by loading external file. . . . . . . . . . ."<<endl;
177  fMergedFileName = Form("%s/%s", fMergeOutput.Data(), fExternalFileName.Data());
178  }
179 
180  cout<<". . .Load inputfile with name: "<<inputfile<<" . . . . . . . ."<<endl;
181  cout<<". . .Continue process by . . . . . . . . . . . ."<<endl;
182  cout<<endl;
183  cout<<"o o o Bad channel analysis o o o"<<endl;
184 
185 
186  BCAnalysis();
187 
188 
189  cout<<endl;
190  cout<<". . .End of process . . . . . . . . . . . . . . . . . . . . ."<<endl;
191  cout<<". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."<<endl;
192 
193 }
194 
198 //________________________________________________________________________
200 {
201  cout<<"o o o Start conversion process o o o"<<endl;
202  cout<<"o o o period: " << fPeriod << ", pass: " << fPass << ", trigger: "<<fTrigger<< endl;
203 
204  //..Create histograms needed for adding all the files together
205  TH1D *hNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun","Number of processed events vs run number",200000,100000,300000);
206  //ELI a little problematic to hard code properties of histograms??
207  TH2F *hCellAmplitude = new TH2F("hCellAmplitude","Cell Amplitude",200,0,10,23040,0,23040);
208  TH2F *hCellTime = new TH2F("hCellTime","Cell Time",250,-275,975,23040,0,23040);
209 
210  //..Open the text file with the run list numbers and run index
211  cout<<"o o o Open .txt file with run indices. Name = " << fRunList << endl;
212  FILE *pFile = fopen(fRunList.Data(), "r");
213  if(!pFile)cout<<"count't open file!"<<endl;
214  Int_t Nentr;
215  Int_t q;
216  Int_t ncols;
217  Int_t nlines = 0 ;
218  Int_t RunId[500] ;
219  while (1)
220  {
221  ncols = fscanf(pFile," %d ",&q);
222  if (ncols< 0) break;
223  RunId[nlines]=q;
224  nlines++;
225  }
226  fclose(pFile);
227 
228 
229  //..Open the different .root files with help of the run numbers from the text file
230  const Int_t nRun = nlines ;
231  TString base;
232  TString infile;
233 
234  cout<<"o o o Start merging process of " << nRun <<" files"<< endl;
235  //..loop over the amount of run numbers found in the previous text file.
236  for(Int_t i = 0 ; i < nRun ; i++)
237  {
238  base = Form("%s/%s/%s/%d", fAnalysisInput.Data(), fPeriod.Data(), fPass.Data(), RunId[i]);
239  if ((fPass=="cpass1_pass2")||(fPass=="cfPass1-2"))
240  {
241  if (fTrigger=="default")
242  {
243  infile = Form("%s_barrel.root",base.Data());
244  }
245  else
246  {
247  infile = Form("%s_outer.root",base.Data());
248  }
249  }
250  else
251  { //..This is a run2 case
252  infile = Form("%s.root",base.Data()) ;
253  }
254 
255  cout<<" o Open .root file with name: "<<infile<<endl;
256  TFile *f = TFile::Open(infile);
257 
258  //..Do some basic checks
259  if(!f)
260  {
261  cout<<"Couldn't open/find .root file: "<<infile<<endl;
262  continue;
263  }
264  TDirectoryFile *dir = (TDirectoryFile *)f->Get(fQADirect);
265  if(!dir)
266  {
267  cout<<"Couln't open directory file in .root file: "<<infile<<", directory: "<<fQADirect<<endl;
268  continue;
269  }
270  TList *outputList = (TList*)dir->Get(fQADirect);
271  if(!outputList)
272  {
273  cout << "Couln't get TList from directory file: "<<fQADirect<<endl;
274  continue;
275  }
276 
277  //ELI should one maybe clone the hAmpId histos eg to hCellAmplitude, then one does't need to hard code them.
278  TH2F *hAmpId;
279  TH2F *hTimeId;
280  TH2F *hNEvents;
281 
282  hAmpId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
283  if(!hAmpId)
284  {
285  Printf("hAmpId not found");
286  outputList->ls();
287  continue;
288  }
289  hTimeId =(TH2F *)outputList->FindObject("EMCAL_hTimeId");
290  if(!hTimeId)
291  {
292  Printf("hTimeId not found");
293  outputList->ls();
294  continue;
295  }
296  hNEvents =(TH2F *)outputList->FindObject("hNEvents");
297  if(!hNEvents)
298  {
299  Printf("hNEvents not found");
300  outputList->ls();
301  continue;
302  }
303  Nentr = (Int_t)hNEvents->GetEntries();
304 
305  //..does that mean do not merge small files?
306  if (Nentr<100)
307  {
308  cout <<" o File to small to be merged. Only N entries " << Nentr << endl;
309  continue ;
310  }
311  cout <<" o File with N entries " << Nentr<<" will be merged"<< endl;
312 
313  hNEventsProcessedPerRun->SetBinContent(RunId[i]-100000,(Double_t)Nentr);
314  hCellAmplitude->Add(hAmpId);
315  hCellTime->Add(hTimeId);
316 
317  outputList->Delete();
318  dir->Delete();
319  f->Close();
320  delete f;
321  }
322 
323  //.. Save the merged histograms
324  cout<<"o o o Save the merged histogramms to .root file with name: "<<fMergedFileName<<endl;
325  TFile *BCF = TFile::Open(fMergedFileName,"recreate");
326  hNEventsProcessedPerRun->Write();
327  hCellAmplitude->Write();
328  hCellTime->Write();
329  BCF->Close();
330  cout<<"o o o End conversion process o o o"<<endl;
331  return fMergedFileName;
332 }
333 
334 
338 //..Run over the list of stored period analysed settings (in fAnalysisVector) and
339 //..pass them to the function PeriodAnalysis()
340 //..The last PeriodAnalysis() is always executed with criteria 7
341 //_________________________________________________________________________
343 {
344  TFile::Open(fMergedFileName);
345  cout<<"o o o Bad channel analysis o o o"<<endl;
346  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
347  //.. DEAD CELLS
348  //.. Flag dead cells with fFlag=1
349  //.. this excludes cells from analysis (will not appear in results)
350  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
351  cout<<"o o o Flag Dead Cells o o o"<<endl;
352  FlagAsDead();
353 
354  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
355  //.. BAD CELLS
356  //.. Flag bad cells with fFlag= 2 or 3
357  //.. this excludes cells from subsequent analysis
358  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
359  TArrayD PeriodArray;
360  for(Int_t i=0;i<fAnalysisVector.size();i++)
361  {
362  PeriodArray=fAnalysisVector.at(i);
363  PeriodAnalysis(PeriodArray.At(0),PeriodArray.At(1),PeriodArray.At(2),PeriodArray.At(3));
364  }
365 
366  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
367  //..In the end summarize results
368  //..in a .pdf and a .txt file
369  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
371  cout<<"o o o End of bad channel analysis o o o"<<endl;
372 }
373 
381 //________________________________________________________________________
382 void AliAnaCaloChannelAnalysis::AddPeriodAnalysis(Int_t criteria, Double_t nsigma, Double_t emin, Double_t emax)
383 {
384  TArrayD periodArray;
385  periodArray.Set(4);
386  periodArray.AddAt((Double_t)criteria,0);
387  periodArray.AddAt(nsigma,1);
388  periodArray.AddAt(emin,2);
389  periodArray.AddAt(emax,3);
390  fAnalysisVector.push_back(periodArray);
391 }
392 
393 //..This function does perform different checks depending on the given criterium variable
394 //..diffrent possibilities for criterium are:
395 // 1 : average E for E>Emin and E<Emax
396 // 2 : entries for E>Emin and E<Emax
397 
398 // 3 : kiĀ²/ndf (from fit of each cell Amplitude between Emin and Emax)
399 // 4 : A parameter (from fit of each cell Amplitude between Emin and Emax)
400 // 5 : B parameter (from fit of each cell Amplitude between Emin and Emax)
401 // 6 :
402 // 7 : give bad + dead list
403 //_________________________________________________________________________
404 void AliAnaCaloChannelAnalysis::PeriodAnalysis(Int_t criterum, Double_t nsigma, Double_t emin, Double_t emax)
405 {
406  //ELI criterum should be between 1-4
407 
408  cout<<""<<endl;
409  cout<<""<<endl;
410  cout<<""<<endl;
411  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;
412  cout<<"o o o PeriodAnalysis for flag "<<criterum<<" o o o"<<endl;
413  cout<<"o o o Done in the energy range E "<<emin<<"-"<<emax<<endl;
414 
415  Int_t CellID, nb1=0, nb2=0;
416  TString output;
417 
418  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
419  //.. ANALYSIS OF CELLS WITH ENTRIES
420  //.. Build average distributions and fit them
421  //.. Three tests for bad cells:
422  //.. 1) Average energy per hit
423  //.. 2) Average hit per event
424  //.. 3) ...
425  //.. 4) ...
426  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
427  TH1F* hisogram;
428  if(criterum < 6)cout<<"o o o Analyze average cell distributions o o o"<<endl;
429  //..For case 1 or 2
430  if(criterum < 3) hisogram = BuildHitAndEnergyMean(criterum, emin, emax,nsigma);
431  //..For case 3, 4 or 5
432  else if (criterum < 6) TestCellShapes(criterum, emin, emax, nsigma);
433 
434  Int_t dnbins = 200;
435  if(criterum==1) FlagAsBad(criterum, hisogram, nsigma, dnbins,-1);
436  if(criterum==2 && emin==0.5) FlagAsBad(criterum, hisogram, nsigma, dnbins*9000,-1); //ELI I did massivley increase the binning now but it helps a lot
437  if(criterum==2 && emin>0.5) FlagAsBad(criterum, hisogram, nsigma, dnbins*17,-1);
438 
439  /*
440  if(criterum==3) FlagAsBad(criterum, hisogram, nsigma, dnbins, maxval3);
441  if(criterum==4) FlagAsBad(criterum, hisogram, nsigma, dnbins, maxval1);
442  if(criterum==5) FlagAsBad(criterum, hisogram, nsigma, dnbins, maxval2);
443  */
444 
445  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
446  //.. RESULTS
447  //.. 1) Print the bad cells
448  //.. and write the results to a file
449  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
450 
451  //..Print the results on the screen and
452  //..write the results in a file
453  output.Form("%s/Criterion%d_Emin-%.2f_Emax-%.2f.txt", fAnalysisOutput.Data(), criterum,emin,emax);
454  ofstream file(output, ios::out | ios::trunc);
455  if(!file)
456  {
457  cout<<"#### Major Error. Check the textfile!"<<endl;
458  }
459  file<<"Criterion : "<<criterum<<", emin = "<<emin<<" GeV"<<", emax = "<<emax<<" GeV"<<endl;
460  file<<"Bad by lower value : "<<endl;
461  cout<<" o bad cells by lower value (for cell E between "<<emin<<"-"<<emax<<")"<<endl;
462  cout<<" ";
463  nb1=0;
464  for(CellID=0;CellID<fNoOfCells;CellID++)
465  {
466  if(fFlag[CellID]==2)
467  {
468  nb1++;
469  file<<CellID<<", ";
470  cout<<CellID<<",";
471  }
472  }
473  file<<"("<<nb1<<")"<<endl;
474  cout<<"("<<nb1<<")"<<endl;
475  file<<"Bad by higher value : "<<endl;
476  cout<<" o bad cells by higher value (for cell E between "<<emin<<"-"<<emax<<")"<<endl;
477  cout<<" ";
478  nb2=0;
479  for(CellID=0;CellID<fNoOfCells;CellID++)
480  {
481  if(fFlag[CellID]==3)
482  {
483  nb2++;
484  file<<CellID<<", ";
485  cout<<CellID<<",";
486  }
487  }
488  file<<"("<<nb2<<")"<<endl;
489  cout<<"("<<nb2<<")"<<endl;
490 
491  file<<"Total number of bad cells"<<endl;
492  file<<"("<<nb1+nb2<<")"<<endl;
493  file.close();
494  cout<<" o Total number of bad cells "<<endl;
495  cout<<" ("<<nb1+nb2<<")"<<endl;
496 
497 }
498 
502 //_________________________________________________________________________
503 TH1F* AliAnaCaloChannelAnalysis::BuildHitAndEnergyMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
504 {
505  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
506 
507  cout<<" o Calculate average cell hit per event and average cell energy per hit "<<endl;
508  //..binning parameters
509  Int_t ncells = hCellAmplitude->GetNbinsY();
510  Double_t amin = hCellAmplitude->GetYaxis()->GetXmin();
511  Double_t amax = hCellAmplitude->GetYaxis()->GetXmax();
512  TH1F *Histogram;
513  if(crit==1)Histogram = new TH1F(Form("hCellEtoNtotal_E%.2f-%.2f",emin,emax),Form("Average energy per hit, %.2f < E < %.2f GeV",emin,emax), ncells,amin,amax);
514  if(crit==2)Histogram = new TH1F(Form("hCellNtotal_E%.2f-%.2f",emin,emax),Form("Number of hits per events, %.2f < E < %.2f GeV",emin,emax), ncells,amin,amax);
515  Histogram->SetXTitle("AbsId");
516  Histogram->SetYTitle("Av. hits per events");
517  Histogram->GetXaxis()->SetNdivisions(505);
518 
519  TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
520  Double_t totalevents = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
521 
522  //..here the average hit per event and the average energy per hit is caluclated for each cell.
523  for (Int_t cell = 0; cell < fNoOfCells; cell++)
524  {
525  Double_t Esum = 0;
526  Double_t Nsum = 0;
527 
528  for (Int_t j = 1; j <= hCellAmplitude->GetNbinsX(); j++)
529  {
530  Double_t E = hCellAmplitude->GetXaxis()->GetBinCenter(j);
531  Double_t N = hCellAmplitude->GetBinContent(j, cell+1);
532  //..investigate only cells that were not flagged as dead ore bad
533  if (E < emin || E > emax || fFlag[cell]!=0) continue;
534  Esum += E*N;
535  Nsum += N;
536  }
537  //..Set the values only for cells that are not yet marked as bad
538  if(fFlag[cell]==0)
539  {
540  if(totalevents> 0. && crit==2)Histogram->SetBinContent(cell+1, Nsum/totalevents); //..number of hits per event
541  if(Nsum > 0. && crit==1)Histogram->SetBinContent(cell+1, Esum/Nsum); //..average energy per hit
542  }
543  }
544  delete hCellAmplitude;
545 
546  return Histogram;
547 }
549 //_________________________________________________________________________
550 TH1F* AliAnaCaloChannelAnalysis::BuildTimeMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
551 {
552  TH2 *hCellTime = (TH2*) gFile->Get("hCellTime");
553  TH1F *Histogram = new TH1F(Form("hSomethingWithTime_E%.2f-%.2f",emin,emax),Form("I don't know, %.2f < E < %.2f GeV",emin,emax), 2,0,2);
554 
555  return Histogram;
556 }
560 //_________________________________________________________________________
561 void AliAnaCaloChannelAnalysis::TestCellShapes(Int_t crit, Double_t fitemin, Double_t fitemax, Double_t nsigma)
562 {
563  TString hname= "hCellAmplitude";
564  Int_t dnbins = 1000;
565  TH2 *hCellAmplitude = (TH2*) gFile->Get(Form("%s",(const char*)hname));
566 
567  // binning parameters
568  Int_t ncells = hCellAmplitude->GetNbinsY();
569  Double_t amin = hCellAmplitude->GetYaxis()->GetXmin();
570  Double_t amax = hCellAmplitude->GetYaxis()->GetXmax();
571  cout << "ncells " << ncells << " amin = " << amin << "amax = " << amax<< endl;
572 
573  // initialize histograms
574  TH1 *hFitA = new TH1F(Form("hFitA_%s",(const char*)hname),"Fit A value", ncells,amin,amax);
575  hFitA->SetXTitle("AbsId");
576  hFitA->SetYTitle("A");
577 
578  TH1 *hFitB = new TH1F(Form("hFitB_%s",(const char*)hname),"Fit B value", ncells,amin,amax);
579  hFitB->SetXTitle("AbsId");
580  hFitB->SetYTitle("B");
581 
582  TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",(const char*)hname),"Fit #chi^{2}/ndf value", ncells,amin,amax);
583  hFitChi2Ndf->SetXTitle("AbsId");
584  hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
585 
586  Double_t maxval1=0., maxval2=0., maxval3=0.;
587  Double_t prev=0., MSA=0., AvA = 0. ; //those param are used to automaticaly determined a reasonable maxval1
588  Double_t prev2=0., MSB=0., AvB = 0. ; //those param are used to automaticaly determined a reasonable maxval2
589  Double_t prev3=0., MSki2=0., Avki2 = 0. ; //those param are used to automaticaly determined a reasonable maxval3
590  Double_t ki2=0.0 ;
591  for (Int_t k = 1; k <= ncells; k++)
592  {
593  TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
594  TH1 *hCell = hCellAmplitude->ProjectionX("",k,k);
595  if (hCell->GetEntries() == 0) continue;
596  // hCell->Rebin(3);
597  hCell->Fit(fit, "0QEM", "", fitemin, fitemax);
598  delete hCell;
599 
600  if(fit->GetParameter(0) < 5000.)
601  {
602  hFitA->SetBinContent(k, fit->GetParameter(0));
603  if(k<3000)
604  {
605  AvA += fit->GetParameter(0);
606  if(k==2999) maxval1 = AvA/3000. ;
607  if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
608  else MSA -= (fit->GetParameter(0) - prev) ;
609  prev = fit->GetParameter(0);
610  }
611  else
612  {
613  if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
614  {
615  maxval1 = fit->GetParameter(0);
616  }
617  }
618  }
619  else hFitA->SetBinContent(k, 5000.);
620 
621  if(fit->GetParameter(1) < 5000.)
622  {
623  hFitB->SetBinContent(k, fit->GetParameter(1));
624  if(k<3000)
625  {
626  AvB += fit->GetParameter(1);
627  if(k==2999) maxval2 = AvB/3000. ;
628  if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
629  else MSB -= (fit->GetParameter(1) - prev2) ;
630  prev2 = fit->GetParameter(1);
631  }
632  else
633  {
634  if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
635  {
636  maxval2 = fit->GetParameter(1);
637  }
638  }
639  }
640  else hFitB->SetBinContent(k, 5000.);
641 
642 
643  if (fit->GetNDF() != 0 ) ki2 = fit->GetChisquare()/fit->GetNDF();
644  else ki2 = 1000.;
645 
646  if(ki2 < 1000.)
647  {
648  hFitChi2Ndf->SetBinContent(k, ki2);
649  if(k<3000)
650  {
651  Avki2 += ki2;
652  if(k==2999) maxval3 = Avki2/3000. ;
653  if (prev3 < ki2) MSki2 += ki2 - prev3;
654  else MSki2 -= (ki2 - prev3) ;
655  prev3 = ki2;
656  }
657  else
658  {
659  if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
660  {
661  maxval3 = ki2;
662  }
663  }
664  }
665  else hFitChi2Ndf->SetBinContent(k, 1000.);
666 
667  delete fit ;
668  }
669 
670  delete hCellAmplitude;
671 
672  // if you have problem with automatic parameter :
673  // maxval1 =
674  // maxval2 =
675  // maxval3 =
676  /*
677  if(crit==3)
678  Process(crit, hFitChi2Ndf, nsigma, dnbins, maxval3);
679  if(crit==4)
680  Process(crit, hFitA, nsigma, dnbins, maxval1);
681  if(crit==5)
682  Process(crit, hFitB, nsigma, dnbins, maxval2);
683  */
684 
685  //ELI something like thie in the future:
686  //return histogram;
687 }
688 
689 
694 //_________________________________________________________________________
696 {
697  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
698  Int_t SumOfExcl=0;
699  //..Direction of cell ID
700  for (Int_t cell = 0; cell < fNoOfCells; cell++)
701  {
702  Double_t Nsum = 0;
703  //..Direction of amplitude
704  for (Int_t amp = 1; amp <= hCellAmplitude->GetNbinsX(); amp++)
705  {
706  //..cellID+1 = histogram bin
707  Double_t N = hCellAmplitude->GetBinContent(amp,cell+1);
708  Nsum += N;
709  }
710  //..If the amplitude in one cell is basically 0
711  //..mark the cell as excluded
712  //ELI I just wonder how you can have less than one but more than 0.5
713  //shouldnt everything below 1 be excluded?
714  if(Nsum >= 0.5 && Nsum < 1)cout<<"-----------------------small but non zero!!!!"<<endl;
715  if(Nsum < 0.5 && Nsum != 0)cout<<"-----------------------non zero!!!!"<<endl;
716 
717  if(Nsum < 0.5 && fFlag[cell]==0)
718  {
719  //..Cell flagged as dead.
720  //..Flag only if it hasn't been flagged before
721  fFlag[cell] =1;
722  SumOfExcl++;
723  }
724  }
725  delete hCellAmplitude;
726  cout<<" o Number of dead cells: "<<SumOfExcl<<endl;
727  cout<<" ("<<SumOfExcl<<")"<<endl;
728 }
729 
739 //_________________________________________________________________________
740 void AliAnaCaloChannelAnalysis::FlagAsBad(Int_t crit, TH1* inhisto, Double_t nsigma, Int_t dnbins, Double_t dmaxval)
741 {
742  gStyle->SetOptStat(1); // MG modif
743  gStyle->SetOptFit(1); // MG modif
744 
745  if(crit==1)cout<<" o Fit average energy per hit distribution"<<endl;
746  if(crit==2)cout<<" o Fit average hit per event distribution"<<endl;
747 
748  Int_t CellColumn=0,CellRow=0;
749  Int_t CellColumnAbs=0,CellRowAbs=0;
750  Int_t Trash;
751 
752  TString HistoName=inhisto->GetName();
753  Double_t goodmax= 0. ;
754  Double_t goodmin= 0. ;
755  if (dmaxval < 0.)
756  {
757  dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
758  if(crit==2 && dmaxval > 1) dmaxval =1. ;
759  }
760 
761  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
762  //. . .build the distribution of average values
763  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
764  TH1 *distrib = new TH1F(Form("%sDistr",(const char*)HistoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
765  distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
766  distrib->SetYTitle("Entries");
767  TH1 *distrib_EMCal = new TH1F(Form("%sDistr_EMCal",(const char*)HistoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
768  TH1 *distrib_DCal = new TH1F(Form("%sDistr_DCal",(const char*)HistoName), "", dnbins, inhisto->GetMinimum(), dmaxval);
769 
770  //..build two dimensional histogram with values row vs. column
771  TH2F *Plot2D = new TH2F(Form("%s_HitRowColumn",(const char*)HistoName),Form("%s_HitRowColumn",(const char*)HistoName),fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
772  Plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
773  Plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
774 
775  for (Int_t cell = 0; cell < fNoOfCells; cell++)
776  {
777  //..Do that only for cell ids also accepted by the code
778  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
779 
780  //..Fill histograms only for cells that are not yet flagged as bad
781  if(fFlag[cell]==0)
782  {
783  //..fill the distribution of avarge cell values
784  distrib->Fill(inhisto->GetBinContent(cell+1));
785  if(cell<fCellStartDCal)distrib_EMCal->Fill(inhisto->GetBinContent(cell+1));
786  else distrib_DCal ->Fill(inhisto->GetBinContent(cell+1));
787  //..Get Row and Collumn for cell ID
788  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,CellColumn,CellRow,Trash,CellColumnAbs,CellRowAbs);
789  if(CellColumnAbs> fNMaxColsAbs || CellRowAbs>fNMaxRowsAbs)
790  {
791  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
792  cout<<"current col: "<<CellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
793  cout<<"current row: "<<CellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
794  }
795  Plot2D->SetBinContent(CellColumnAbs,CellRowAbs,inhisto->GetBinContent(cell+1));
796  }
797  }
798 
799  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
800  //. . .draw histogram + distribution
801  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
802 
803  TCanvas *c1 = new TCanvas(HistoName,HistoName,900,900);
804  c1->ToggleEventStatus();
805  TPad* upperPad = new TPad("upperPad", "upperPad",.005, .5, .995, .995);
806  TPad* lowerPadLeft = new TPad("lowerPadL", "lowerPadL",.005, .005, .5, .5);
807  TPad* lowerPadRight = new TPad("lowerPadR", "lowerPadR",.5, .005, .995, .5);
808  upperPad->Draw();
809  lowerPadLeft->Draw();
810  lowerPadRight->Draw();
811 
812  upperPad->cd();
813  upperPad->SetLeftMargin(0.045);
814  upperPad->SetRightMargin(0.05);
815  upperPad->SetLogy();
816  inhisto->SetTitleOffset(0.6,"Y");
817  inhisto->GetXaxis()->SetRangeUser(0,fNoOfCells+1);
818 
819  inhisto->SetLineColor(kBlue+1);
820  inhisto->Draw();
821 
822  lowerPadRight->cd();
823  lowerPadRight->SetLeftMargin(0.09);
824  lowerPadRight->SetRightMargin(0.1);
825  Plot2D->Draw("colz");
826 
827  lowerPadLeft->cd();
828  lowerPadLeft->SetLeftMargin(0.09);
829  lowerPadLeft->SetRightMargin(0.07);
830  lowerPadLeft->SetLogy();
831  distrib->SetLineColor(kBlue+1);
832  distrib->Draw();
833  distrib_EMCal->SetLineColor(kGreen+1);
834  distrib_EMCal->DrawCopy("same");
835  distrib_DCal->SetLineColor(kMagenta+1);
836  distrib_DCal->DrawCopy("same");
837 
838  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
839  //. . .fit histogram
840  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
841  Int_t higherbin=0,i;
842  for(i = 2; i <= dnbins; i++)
843  {
844  if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i)) higherbin = i ;
845  }
846  //..good range is around the max value as long as the
847  //..bin content is larger than 2 entries
848  for(i = higherbin ; i<=dnbins ; i++)
849  {
850  if(distrib->GetBinContent(i)<2) break ;
851  goodmax = distrib->GetBinCenter(i);
852  }
853  for(i = higherbin ; i>1 ; i--)
854  {
855  if(distrib->GetBinContent(i)<2) break ;
856  goodmin = distrib->GetBinLowEdge(i);
857  }
858  //cout<<"higherbin : "<<higherbin<<endl;
859  //cout<<"good range : "<<goodmin<<" - "<<goodmax<<endl;
860 
861  TF1 *fit2 = new TF1("fit2", "gaus");
862  //..start the fit with a mean of the highest value
863  fit2->SetParameter(1,higherbin);
864 
865  distrib->Fit(fit2, "0LQEM", "", goodmin, goodmax);
866  Double_t sig, mean, chi2ndf;
867  // Marie midif to take into account very non gaussian distrig
868  mean = fit2->GetParameter(1);
869  sig = fit2->GetParameter(2);
870  chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
871 
872  if (mean <0.) mean=0.; //ELI is this not a highly problematic case??
873 
874  goodmin = mean - nsigma*sig ;
875  goodmax = mean + nsigma*sig ;
876 
877  if (goodmin<0) goodmin=0.;
878 
879  cout<<" o Result of fit: "<<endl;
880  cout<<" o "<<endl;
881  cout<<" o Mean: "<<mean <<" sigma: "<<sig<<endl;
882  cout<<" o good range : "<<goodmin <<" - "<<goodmax<<endl;
883 
884  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
885  //. . .Add info to histogram
886  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
887  TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
888  lline->SetLineColor(kGreen+2);
889  lline->SetLineStyle(7);
890  lline->Draw();
891 
892  TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
893  rline->SetLineColor(kGreen+2);
894  rline->SetLineStyle(7);
895  rline->Draw();
896 
897  TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
898  leg->AddEntry(lline, "Good region boundary","l");
899  leg->Draw("same");
900 
901  fit2->SetLineColor(kOrange-3);
902  fit2->SetLineStyle(1);//7
903  fit2->Draw("same");
904 
905  TLatex* text = 0x0;
906  if(crit==1) text = new TLatex(0.2,0.8,Form("Good range: %.2f-%.2f",goodmin,goodmax));
907  if(crit==2) text = new TLatex(0.2,0.8,Form("Good range: %.2f-%.2fx10^-5",goodmin*100000,goodmax*100000));
908  text->SetTextSize(0.06);
909  text->SetNDC();
910  text->SetTextColor(1);
911  text->Draw();
912  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
913  //. . .Save histogram
914  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
915  c1->Update();
916  TString name =Form("%s/criteria-_%d.gif", fAnalysisOutput.Data(), crit);
917  if(crit==1)name=Form("%s/AverageEperHit_%s.gif", fAnalysisOutput.Data(), (const char*)HistoName);
918  if(crit==2)name=Form("%s/AverageHitperEvent_%s.gif", fAnalysisOutput.Data(), (const char*)HistoName);
919  c1->SaveAs(name);
920 
921 
922  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
923  //. . . Mark the bad cells in the pflag array
924  //. . .(0= bad because cell average value lower than min allowed)
925  //. . .(2= bad because cell average value higher than max allowed)
926  //. . .(1 by default)
927  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
928  cout<<" o Flag bad cells that are outside the good range "<<endl;
929  for(Int_t cell = 0; cell < fNoOfCells; cell++)
930  {
931  //cel=0 and bin=1, cel=1 and bin=2
932  if (inhisto->GetBinContent(cell+1) <= goodmin && fFlag[cell]!=1) //ELI <= throws out zeros, might not be dead but still have zero entries in a given energy range
933  {
934  fFlag[cell]=2;
935  if(inhisto->GetBinContent(cell+1)==goodmin)cout<<"cell number with zero= "<<cell<<endl;
936  }
937  if (inhisto->GetBinContent(cell+1) > goodmax && fFlag[cell]!=1)
938  {
939  fFlag[cell]=3;
940  }
941  }
942  cout<<" o "<<endl;
943 
944 }
945 
950 //________________________________________________________________________
952 {
953  Int_t CellID, DCalCells=0, EMCalCells=0;
954  TString CellSummaryFile, DeadPdfName, BadPdfName, RatioOfBad,GoodCells;
955 
956  DeadPdfName = Form("%s/%s%sDC_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
957  BadPdfName = Form("%s/%s%sBC_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
958  CellSummaryFile = Form("%s/%s%sBC_SummaryResults_V%i.txt", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial); ;
959  RatioOfBad = Form("%s/%s%sBCRatio_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
960  GoodCells = Form("%s/%s%sGoodCells_SummaryResults_V%i.pdf", fAnalysisOutput.Data(), fPeriod.Data(), fPass.Data(), fTrial);
961  cout<<" o Final results o "<<endl;
962  cout<<" o write results into .txt file: "<<CellSummaryFile<<endl;
963  cout<<" o write results into .pdf file: "<<BadPdfName<<endl;
964  ofstream file(CellSummaryFile, ios::out | ios::trunc);
965  if(file)
966  {
967  file<<"Dead cells : "<<endl;
968  cout<<" o Dead cells : "<<endl;
969  EMCalCells =0;
970  DCalCells =0;
971  for(CellID=0; CellID<fNoOfCells; CellID++)
972  {
973  if(CellID==0)
974  {
975  file<<"In EMCal : "<<endl;
976  cout<<" o In EMCal : "<<endl;
977  }
978  if(CellID==fCellStartDCal)
979  {
980  file<<"In DCal : "<<endl;
981  cout<<endl;
982  cout<<" o In DCal : "<<endl;
983  }
984  if(fFlag[CellID]==1)
985  {
986  file<<CellID<<"\n" ;
987  cout<<CellID<<"," ;
988  if(CellID<fCellStartDCal)EMCalCells++;
989  else DCalCells++;
990  }
991  }
992  cout<<endl;
993  file<<"EMCal ("<<EMCalCells<<" ="<<100*EMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<DCalCells<<" ="<<100*DCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
994  cout<<" o EMCal ("<<EMCalCells<<" ="<<100*EMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<DCalCells<<" ="<<100*DCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
995 
996  file<<"Bad cells: "<<endl;
997  cout<<" o Bad cells: "<<endl;
998  EMCalCells =0;
999  DCalCells =0;
1000  for(CellID=0;CellID<fNoOfCells;CellID++)
1001  {
1002  if(CellID==0)
1003  {
1004  file<<"In EMCal : "<<endl;
1005  cout<<" o In EMCal : "<<endl;
1006  }
1007  if(CellID==fCellStartDCal)
1008  {
1009  file<<"In DCal : "<<endl;
1010  cout<<endl;
1011  cout<<" o In DCal : "<<endl;
1012  }
1013  if(fFlag[CellID]>1)
1014  {
1015  file<<CellID<<"\n" ;
1016  cout<<CellID<<"," ;
1017  if(CellID<fCellStartDCal)EMCalCells++;
1018  else DCalCells++;
1019  }
1020  }
1021  cout<<endl;
1022  file<<"EMCal ("<<EMCalCells<<" ="<<100*EMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<DCalCells<<" ="<<100*DCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1023  cout<<" o EMCal ("<<EMCalCells<<" ="<<100*EMCalCells/(1.0*fCellStartDCal)<<"%), DCal ("<<DCalCells<<" ="<<100*DCalCells/(1.0*fNoOfCells-fCellStartDCal)<<"%)"<<endl;
1024  }
1025  file.close();
1026 
1027  TFile::Open(fMergedFileName);
1028  //cout<<" o Save the Dead channel spectra to a .pdf file"<<endl;
1029  //SaveBadCellsToPDF(0,DeadPdfName);
1030  cout<<" o Save the bad channel spectra to a .pdf file"<<endl;
1031  SaveBadCellsToPDF(1,BadPdfName) ;
1032  SaveBadCellsToPDF(10,RatioOfBad) ; //..Special case
1033  //SaveBadCellsToPDF(2,GoodCells) ; //..Special case all good cells to check
1034 
1035 
1036  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1037  //..Plot some summary canvases
1038  //. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1039  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
1040  TH2 *hCellTime = (TH2*) gFile->Get("hCellTime");
1041 
1042  TCanvas *c1 = new TCanvas("CellProp","summary of cell properties",1000,500);
1043  c1->ToggleEventStatus();
1044  c1->Divide(2);
1045  c1->cd(1);
1046  //lowerPadRight->SetLeftMargin(0.09);
1047  //lowerPadRight->SetRightMargin(0.06);
1048  hCellAmplitude->Draw("colz");
1049  c1->cd(2);
1050  hCellTime->Draw("colz");
1051  c1->Update();
1052  TString name =Form("%s/CellProperties.gif", fAnalysisOutput.Data());
1053  c1->SaveAs(name);
1054 
1055 
1056  PlotFlaggedCells2D(0); //..all good cells
1057  PlotFlaggedCells2D(1); //..all dead cells
1058  PlotFlaggedCells2D(2,3); //..all bad cells
1059 }
1060 
1061 
1062 
1066 //________________________________________________________________________
1067 void AliAnaCaloChannelAnalysis::SaveBadCellsToPDF(Int_t version, TString pdfName)
1068 {
1069  //..version=0 ->Print dead cells
1070  //..version=1 ->Print bad cells
1071  //..ELI can be refined with more versions!
1072  gROOT->SetStyle("Plain");
1073  gStyle->SetOptStat(0);
1074  gStyle->SetFillColor(kWhite);
1075  gStyle->SetTitleFillColor(kWhite);
1076  gStyle->SetPalette(1);
1077 
1078  char title[100];
1079  char name[100];
1080 
1081  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
1082  TH1 *hRefDistr = BuildMeanFromGood();
1083  Int_t FirstCanvas=0;
1084  Bool_t suspicious;
1085  TLatex* text = 0x0;
1086  text = new TLatex(0.2,0.8,"*Candidate*");
1087  text->SetTextSize(0.06);
1088  text->SetTextColor(8);
1089  text->SetNDC();
1090  text->SetTextColor(1);
1091  //..collect cells in an internal vector.
1092  //..when the vector is of size=9 or at the end of being filled
1093  //..plot the channels into a canvas
1094  std::vector<Int_t> ChannelVector;
1095  ChannelVector.clear();
1096  for(Int_t cell=0;cell<fNoOfCells;cell++)
1097  {
1098  if(fFlag[cell]==1 && version==0)ChannelVector.push_back(cell);
1099  if(fFlag[cell]>1 && version==1)ChannelVector.push_back(cell);
1100  if(fFlag[cell]==0 && version==2)ChannelVector.push_back(cell);
1101  if(fFlag[cell]>1 && version==10)ChannelVector.push_back(cell);
1102 
1103  //..when 9 bad cells are collected or we are at the end of the list, fill the canvas
1104  if(ChannelVector.size()==9 || cell == fNoOfCells-1)
1105  {
1106  TString Internal_pdfName=pdfName;
1107  TCanvas *c1 = new TCanvas("badcells","badcells",1000,750);
1108  if(ChannelVector.size() > 6) c1->Divide(3,3);
1109  else if (ChannelVector.size() > 3) c1->Divide(3,2);
1110  else c1->Divide(3,1);
1111 
1112  TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
1113  for(Int_t i=0; i<ChannelVector.size() ; i++)
1114  {
1115  sprintf(name, "Cell %d",ChannelVector.at(i)) ;
1116  TH1 *hCell = hCellAmplitude->ProjectionX(name,ChannelVector.at(i)+1,ChannelVector.at(i)+1);
1117  sprintf(title,"Cell No: %d Entries: %d",ChannelVector.at(i), (Int_t)hCell->GetEntries()) ;
1118 
1119  c1->cd(i%9 + 1);
1120  c1->cd(i%9 + 1)->SetLogy();
1121  if(i==0)
1122  {
1123  leg->AddEntry(hRefDistr,"mean of good","l");
1124  leg->AddEntry(hCell,"current","l");
1125  }
1126  if(version>1)
1127  {
1128  //cout<<"bin nr: "<<ChannelVector.at(i)<<endl;
1129  hCell->Divide(hRefDistr);
1130  suspicious = CheckDistribution(hCell,hRefDistr);
1131  }
1132 
1133  hCell->SetLineColor(2);
1134  hCell->GetXaxis()->SetTitle("E (GeV)");
1135  hCell->GetYaxis()->SetTitle("N Entries");
1136  //hCell->GetYaxis()->SetRangeUser(0,1e6);
1137  hCell->GetXaxis()->SetRangeUser(0.,10.);
1138  hCell->SetLineWidth(1) ;
1139  hCell->SetTitle(title);
1140  hRefDistr->SetLineWidth(1);
1141 
1142  hCell->Draw();
1143  if(version<2)hRefDistr->Draw("same") ;
1144  if(version>1 && suspicious==0)text->Draw(); //..Draw a marker in the histogram that could be miscalibrated and labelled as warm
1145  leg->Draw();
1146  }
1147 
1148  if(ChannelVector.size()<9 || cell == fNoOfCells-1)
1149  {
1150  Internal_pdfName +=")";
1151  c1->Print(Internal_pdfName.Data());
1152  }
1153  else if(FirstCanvas==0)
1154  {
1155  Internal_pdfName +="(";
1156  c1->Print(Internal_pdfName.Data());
1157  FirstCanvas=1;
1158  }
1159  else
1160  {
1161  c1->Print(Internal_pdfName.Data());
1162  }
1163  delete c1;
1164  delete leg;
1165  ChannelVector.clear();
1166  }
1167  }
1168  delete hRefDistr;
1169 }
1173 //_________________________________________________________________________
1175 {
1176  TH2 *hCellAmplitude = (TH2*) gFile->Get("hCellAmplitude");
1177  TH1* hGoodAmp;
1178  TH1* hgoodMean;
1179  Int_t NrGood=0;
1180  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1181  {
1182  if(fFlag[cell]!=0)continue;
1183  NrGood++;
1184  if(NrGood==1)hgoodMean = hCellAmplitude->ProjectionX("hgoodMean",cell+1,cell+1);
1185  else
1186  {
1187  hGoodAmp = hCellAmplitude->ProjectionX("hGoodCells",cell+1,cell+1);
1188  hgoodMean->Add(hGoodAmp);
1189  }
1190  }
1191  hgoodMean->Scale(1.0/NrGood);
1192  return hgoodMean;
1193 }
1197 //_________________________________________________________________________
1198 Bool_t AliAnaCaloChannelAnalysis::CheckDistribution(TH1* ratio, TH1* reference)
1199 {
1200  //histo
1201  Bool_t suspicious=0;
1202  //..Find bin where reference has value 1, and the corresponding x-value
1203  Int_t BinHeihgtOne = reference->FindLastBinAbove(1);
1204  Double_t X_of_BinHeightOne = reference->GetBinCenter(BinHeihgtOne);
1205  Double_t X_of_SecondBin = reference->GetBinCenter(2);
1206  //cout<<"specrta is one at: "<<BinOne<<", this is at: "<<reference->GetBinCenter(BinOne)<<endl;
1207  //..Check the histogram
1208  //..Different checks to see whether the
1209  //..cell is really bad. Set suspicious to 1.
1210 
1211  //..check end of spectrum, should be bigger than 60% of the mean histogram
1212  if(ratio->FindLastBinAbove(0)<ratio->FindBin(X_of_BinHeightOne*0.60))
1213  {
1214  suspicious=1;
1215  }
1216  //..check whether there are large spikes in the histogram
1217  //..compare bin values to mean. If there is a bin value with
1218  //..content 5 times lareger than mean it's suspicious
1219  Double_t mean=0;
1220  for(Int_t i=1;i<BinHeihgtOne;i++)
1221  {
1222  if(ratio->GetBinContent(i)>0)mean+=ratio->GetBinContent(i);
1223  }
1224  mean*=1.0/(BinHeihgtOne-1);//..divide by number of bins
1225  ratio->GetXaxis()->SetRangeUser(X_of_SecondBin,X_of_BinHeightOne);//..zoom in to find the maximum between 0-BinOne
1226  //cout<<"mean: "<<mean<<", max: "<<ratio->GetMaximum()<<endl;
1227  if(ratio->GetMaximum()>mean*5)
1228  {
1229  suspicious=1;
1230  }
1231 
1232 
1233  return suspicious;
1234 }
1236 //_________________________________________________________________________
1237 void AliAnaCaloChannelAnalysis::PlotFlaggedCells2D(Int_t flag1,Int_t flag2,Int_t flag3)
1238 {
1239  //..build two dimensional histogram with values row vs. column
1240  TString HistoName;
1241  HistoName=Form("HitRowColumn_Flag%d",flag1);
1242  TH2F *Plot2D = new TH2F(HistoName,HistoName,fNMaxColsAbs+2,-1.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+2,-1.5,fNMaxRowsAbs+0.5);
1243  Plot2D->GetXaxis()->SetTitle("cell column (#eta direction)");
1244  Plot2D->GetYaxis()->SetTitle("cell row (#phi direction)");
1245 
1246  Int_t CellColumn=0,CellRow=0;
1247  Int_t CellColumnAbs=0,CellRowAbs=0;
1248  Int_t Trash;
1249 
1250  for (Int_t cell = 0; cell < fNoOfCells; cell++)
1251  {
1252  //..Do that only for cell ids also accepted by the code
1253  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
1254 
1255  //..Get Row and Collumn for cell ID c
1256  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,CellColumn,CellRow,Trash,CellColumnAbs,CellRowAbs);
1257  if(CellColumnAbs> fNMaxColsAbs || CellRowAbs>fNMaxRowsAbs)
1258  {
1259  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
1260  cout<<"current col: "<<CellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
1261  cout<<"current row: "<<CellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
1262  }
1263  if(fFlag[cell]==flag1)Plot2D->SetBinContent(CellColumnAbs,CellRowAbs,1);
1264  if(flag2!=-1 && fFlag[cell]==flag2)Plot2D->SetBinContent(CellColumnAbs,CellRowAbs,1);
1265  if(flag3!=-1 && fFlag[cell]==flag3)Plot2D->SetBinContent(CellColumnAbs,CellRowAbs,1);
1266  }
1267  TCanvas *c1 = new TCanvas(HistoName,HistoName,500,500);
1268  c1->ToggleEventStatus();
1269  c1->cd();
1270  //lowerPadRight->SetLeftMargin(0.09);
1271  //lowerPadRight->SetRightMargin(0.06);
1272  Plot2D->Draw("colz");
1273 
1274  TLatex* text = 0x0;
1275  if(flag1==0) text = new TLatex(0.2,0.8,"Good Cells");
1276  if(flag1==1) text = new TLatex(0.2,0.8,"Dead Cells");
1277  if(flag1>1) text = new TLatex(0.2,0.8,"Bad Cells");
1278  text->SetTextSize(0.06);
1279  text->SetNDC();
1280  text->SetTextColor(1);
1281  text->Draw();
1282 
1283  c1->Update();
1284  TString name =Form("%s/2DChannelMap_Flag%d.gif", fAnalysisOutput.Data(), flag1);
1285  c1->SaveAs(name);
1286 
1287 }
TH1F * BuildTimeMean(Int_t crit, Double_t emin, Double_t emax, Double_t nsigma)
Possibility to add this check too, if the time if calibrated.
void SetRunNumber(Int_t run)
const char * title
Definition: MakeQAPdf.C:26
void AddPeriodAnalysis(Int_t criteria, Double_t nsigma, Double_t emin, Double_t emax)
Int_t fNMaxRows
Maximum No of rows in module (phi direction)
Int_t fTrial
Number of trial that this specific analyis is. By default '0'.
TSystem * gSystem
std::vector< TArrayD > fAnalysisVector
Vector of analysis information. Each place is filled with 4 doubles: version, sigma, lower, and upper energy range.
Int_t fNoOfCells
Number of cells in EMCal and DCal.
AliCalorimeterUtils * fCaloUtils
! Calorimeter information for the investigated runs
Int_t fCellStartDCal
ID of the first cell in the DCal.
TString fPass
Pass of the analyzed data.
Int_t * fFlag
! fFlag[CellID] = 0 (ok),1 (dead),2 (bad by lower),3 (bad by upper) start at 0 (cellID 0 = histobin 1...
Int_t fNMaxCols
Maximum No of colums in module (eta direction)
TString fMergeOutput
Here the merged files of a period are saved for a later analysis.
TString fMergedFileName
Filename of the .root file containing the merged runs.
AliEMCALGeometry * GetEMCALGeometry() const
TString fAnalysisInput
Here the .root files of each run of the period are saved.
void TestCellShapes(Int_t crit, Double_t fitemin, Double_t fitemax, Double_t nsigma=4.)
TString fQADirect
Dierctory in the QA.root files where the input histograms are stored.
TString fWorkdir
Directory which contains the folders fMergeOutput, fAnalysisInput and fAnalysisOutput. By default it is './'.
void FlagAsBad(Int_t crit, TH1 *inhisto, Double_t nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1.)
Double_t nsigma
TString fExternalFileName
If you have already a file that contains many runs merged together you can place it in fMergeOutput a...
TString fPeriod
The name of the analyzed period.
TH1F * BuildHitAndEnergyMean(Int_t crit, Double_t emin=0.1, Double_t emax=2., Double_t nsigma=4.)
TString fTrigger
Selected trigger for the analysis.
Bool_t CheckDistribution(TH1 *ratio, TH1 *reference)
Int_t fNMaxRowsAbs
Maximum No of rows in Calorimeter.
Int_t fNMaxColsAbs
Maximum No of colums in Calorimeter.
void PeriodAnalysis(Int_t criterum=7, Double_t nsigma=4.0, Double_t emin=0.1, Double_t emax=2.0)
Analyses cell properties and identifies bad cells.
TString fAnalysisOutput
The list with bad channels and histograms are saved in this folder.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TFile * file
Class with utils specific to calorimeter clusters/cells.
void AccessGeometry(AliVEvent *inputEvent)
void PlotFlaggedCells2D(Int_t flag1, Int_t flag2=-1, Int_t flag3=-1)
Plots a 2D map of flagged cells, Dead, Bad, Good.
void SaveBadCellsToPDF(Int_t version, TString pdfName)
TString fRunListFileName
This is the name of the file with the run numbers to be merged, by default it's 'runList.txt'.
Int_t fCurrentRunNumber
A run number of an analyzed period. This is important for the AliCalorimeterUtils initialization...
TString fRunList
Thats the full path and name of the file which contains a list of all runs to be merged together...
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