AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEMCALTimeCalib.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 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 #include <TChain.h>
17 #include <TTree.h>
18 #include <TFile.h>
19 #include <TH1F.h>
20 #include <TH1D.h>
21 #include <TH2D.h>
22 #include <TH2F.h>
23 #include <TH1C.h>
24 #include <TList.h>
25 #include <TCanvas.h>
26 #include <TGeoManager.h>
27 #include <TRefArray.h>
28 #include <TKey.h>
29 
30 #include "AliLog.h"
31 #include "AliAnalysisTask.h"
32 #include "AliAnalysisManager.h"
33 #include "AliESDEvent.h"
34 #include "AliAODEvent.h"
35 #include "AliVEvent.h"
36 #include "AliESDInputHandler.h"
37 #include "AliAODInputHandler.h"
38 #include "AliESDpid.h"
39 #include "AliTOFcalib.h"
40 #include "AliCDBManager.h"
41 #include "AliRunTag.h"
42 
43 #include "AliTOFT0maker.h"
44 #include "AliVCluster.h"
45 #include "AliESDCaloCluster.h"
46 #include "AliVCaloCells.h"
47 #include "AliESDCaloCells.h"
48 #include "AliAODCaloCluster.h"
49 #include "AliAODCaloCells.h"
50 #include "AliEMCALGeometry.h"
51 #include "AliOADBContainer.h"
52 
54 
58 
59 using std::cout;
60 using std::endl;
61 
62 //________________________________________________________________________
65 : AliAnalysisTaskSE(name),
66  fRunNumber(-1),
67  fTOFmaker(0),
68  fOutputList(0x0),
69  fgeom(0),
70  fGeometryName(),
71  fMinClusterEnergy(0),
72  fMaxClusterEnergy(0),
73  fMinNcells(0),
74  fMaxNcells(0),
75  fMinLambda0(0),
76  fMaxLambda0(0),
77  fMinLambda0LG(0),
78  fMaxLambda0LG(0),
79  fMaxRtrack(0),
80  fMinCellEnergy(0),
81  fReferenceFileName(),
82  fReferenceRunByRunFileName(),
83  fPileupFromSPD(kFALSE),
84  fMinTime(0),
85  fMaxTime(0),
86  fRawTimeNbins (0),
87  fRawTimeMin (0),
88  fRawTimeMax (0),
89  fPassTimeNbins(0),
90  fPassTimeMin (0),
91  fPassTimeMax (0),
92  fEnergyNbins (0),
93  fEnergyMin(0),
94  fEnergyMax(0),
95  fEnergyLGNbins (0),
96  fEnergyLGMin(0),
97  fEnergyLGMax(0),
98  fFineNbins(0),
99  fFineTmin(0),
100  fFineTmax(0),
101  fL1PhaseList(0),
102  fBadReco(kFALSE),
103  fFillHeavyHisto(kFALSE),
104  fBadChannelMapArray(),
105  fBadChannelMapSet(kFALSE),
106  fSetBadChannelMapSource(0),
107  fBadChannelFileName(),
108  fhcalcEvtTime(0),
109  fhEvtTimeHeader(0),
110  fhEvtTimeDiff(0),
111  fhEventType(0),
112  fhTOFT0vsEventNumber(0),
113  fhTcellvsTOFT0(0),
114  fhTcellvsTOFT0HD(0),
115  fhTcellvsSM(0),
116  fhEneVsAbsIdHG(0),
117  fhEneVsAbsIdLG(0),
118  fhTimeVsBC(0),
119  fhTimeSumSq(),
120  fhTimeEnt(),
121  fhTimeSum(),
122  fhTimeLGSumSq(),
123  fhTimeLGEnt(),
124  fhTimeLGSum(),
125  fhAllAverageBC(),
126  fhAllAverageLGBC(),
127  fhRefRuns(0),
128  fhTimeDsup(),
129  fhTimeDsupBC(),
130  fhTimeDsupLG(),
131  fhTimeDsupLGBC(),
132  fhRawTimeVsIdBC(),
133  fhRawTimeSumBC(),
134  fhRawTimeEntriesBC(),
135  fhRawTimeSumSqBC(),
136  fhRawTimeVsIdLGBC(),
137  fhRawTimeSumLGBC(),
138  fhRawTimeEntriesLGBC(),
139  fhRawTimeSumSqLGBC(),
140  fhRawCorrTimeVsIdBC(),
141  fhRawCorrTimeVsIdLGBC(),
142  fhTimeVsIdBC(),
143  fhTimeVsIdLGBC()
144 {
145  for(Int_t i = 0; i < kNBCmask; i++)
146  {
147  fhAllAverageBC[i]=0;
148  fhAllAverageLGBC[i]=0;
149 
150  fhTimeSumSq[i]=0;
151  fhTimeEnt[i]=0;
152  fhTimeSum[i]=0;
153 
154  fhTimeLGSumSq[i]=0;
155  fhTimeLGEnt[i]=0;
156  fhTimeLGSum[i]=0;
157 
158  fhRawTimeVsIdBC[i]=0;
159  fhRawTimeSumBC[i]=0;
160  fhRawTimeEntriesBC[i]=0;
161  fhRawTimeSumSqBC[i]=0;
162 
163  fhRawTimeVsIdLGBC[i]=0;
164  fhRawTimeSumLGBC[i]=0;
166  fhRawTimeSumSqLGBC[i]=0;
167  fhRawCorrTimeVsIdBC[i]=0;
169  fhTimeVsIdBC[i]=0;
170  fhTimeVsIdLGBC[i]=0;
171 
172  }
173 
174  //set default cuts for calibration and geometry name
175  SetDefaultCuts();
176 
177  //T0 TOF time
179 
180  // Define input and output slots here
181  // Input slot #0 works with a TChain
182  DefineInput(0, TChain::Class());
183 
184  // Output slot #0 id reserved by the base class for AOD
185  // Output slot #1 writes into a TH1 container
186  DefineOutput(1, TList::Class());
187 
188 } // End ctor
189 
190 //_____________________________________________________________________
196 //void AliAnalysisTaskEMCALTimeCalib::LocalInit()
197 //{
198 // AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::LocalInit()");
199 //}
200 
202 //_____________________________________________________________________
204 {
205  if(fReferenceFileName.Length()!=0){
206  TFile *myFile = TFile::Open(fReferenceFileName.Data());
207  AliInfo(Form("Reference file: %s, pointer %p",fReferenceFileName.Data(),myFile));
208  if(myFile==0x0) {
209  AliFatal("*** NO REFERENCE FILE");
210  } else {
211  AliDebug(1,"*** OK TFILE");
212  // connect ref run here
213  for(Int_t i = 0; i < kNBCmask; i++)
214  {
215  fhAllAverageBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvBC%d",i));
216  if(fhAllAverageBC[i]==0x0) AliFatal(Form("Reference histogram for BC%d does not exist",i));
217  if(fhAllAverageBC[i]->GetEntries()==0)AliWarning(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
218  fhAllAverageLGBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvLGBC%d",i));
219  if(fhAllAverageLGBC[i]==0x0) AliFatal(Form("Reference LG histogram for BC%d does not exist",i));
220  if(fhAllAverageLGBC[i]->GetEntries()==0)AliFatal(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
221  }
222 
223  AliDebug(1,Form("hAllAverage entries BC0 %d", (Int_t)fhAllAverageBC[0]->GetEntries() ));
224  AliDebug(1,Form("hAllAverage entries BC2 %d",(Int_t)fhAllAverageBC[2]->GetEntries() ));
225  AliDebug(1,Form("hAllAverageLG entries BC0 %d", (Int_t)fhAllAverageLGBC[0]->GetEntries() ));
226  AliDebug(1,Form("hAllAverageLG entries BC2 %d",(Int_t)fhAllAverageLGBC[2]->GetEntries() ));
227 
228  }
229  } else { //end of reference file is provided
230  AliFatal("You require to load reference histos from file but FILENAME is not provided");
231  }
232 } // End of AliAnalysisTaskEMCALTimeCalib::LoadReferenceHistos()
233 
236 //_____________________________________________________________________
238 {
239  // connect ref run here
240  if(fReferenceRunByRunFileName.Length()!=0){
241  TFile *referenceFile = TFile::Open(fReferenceRunByRunFileName.Data());
242  if(referenceFile==0x0) {
243  AliFatal("*** NO REFERENCE R-B-R FILE");
244  return;
245  } else {
246  AliInfo(Form("Reference R-b-R file: %s, pointer %p",fReferenceRunByRunFileName.Data(),referenceFile));
247 
248  //load L1 phases to memory
249  fL1PhaseList=new TObjArray(referenceFile->GetNkeys());
250  TIter next(referenceFile->GetListOfKeys());
251  TKey *key;
252  while ((key=(TKey*)next())) {
253  fL1PhaseList->AddLast((TH1F*)referenceFile->Get(key->GetName()) );
254  //printf("key: %s points to an object of class: %s at %dn",key->GetName(),key->GetClassName(),key->GetSeekKey());
255  }
256  }
257  } else { //reference file is not provided
258  AliFatal("You require to load reference run-by-run histos from file but FILENAME is not provided");
259  return;
260  }
261 } // End of AliAnalysisTaskEMCALTimeCalib::LoadReferenceRunByRunHistos()
262 
267 {
268  fhRefRuns=NULL;
269  if(!fL1PhaseList) {
270  AliFatal("Array with reference L1 phase histograms do not exist in memory");
271  return;
272  }
273  if(fRunNumber<0) {
274  AliFatal("Negative run number");
275  return;
276  }
277 
278  fhRefRuns=(TH1C*)fL1PhaseList->FindObject(Form("h%d",fRunNumber));
279  if(fhRefRuns==0x0){
280  AliError(Form("Reference histogram for run %d does not exist. Use Default",fRunNumber));
281  fhRefRuns=(TH1C*)fL1PhaseList->FindObject("h0");
282  }
283  if(fhRefRuns==0x0) {
284  AliFatal(Form("No default histogram with L1 phases! Add default histogram to file %s!!!",fReferenceRunByRunFileName.Data()));
285  return;
286  }
287 
288  AliDebug(1,Form("Reference R-b-R histo %p, list %p, run number %d",fhRefRuns,fL1PhaseList,fRunNumber));
289  if(fhRefRuns->GetEntries()==0)AliWarning("fhRefRuns->GetEntries() = 0");
290  AliDebug(1,Form("hRefRuns entries %d", (Int_t)fhRefRuns->GetEntries() ));
291 }
292 
293 //_____________________________________________________________________
297 {
298  AliDebug(1,"AnalysisTaskEMCalTimeCalib::NotifyRun()");
299  AliDebug(2,Form("Notify(): EMCal geometry: fgeom = %p, fGeometryName=%s\n ",fgeom,fGeometryName.Data()));
300 
301  if (!InputEvent())
302  {
303  AliFatal("ERROR: InputEvent not set");
304  return;
305  }
306  else AliDebug(1,"Good, InputEvent set");
307 
308  // AliInfo(Form("NotifyRun, fCurrentRunnumber %d",fCurrentRunNumber));
309  fRunNumber = InputEvent()->GetRunNumber();
310  AliDebug(1,Form("RunNumber %d", fRunNumber));
311 
312  // Init EMCAL geometry
313  if (!fgeom) SetEMCalGeometry();
314  //Init EMCAL geometry done
315 
316  //set L1 phases for current run
317  if(fReferenceRunByRunFileName.Length()!=0)
319 
320  // set bad channel map
322 
323  return;
324 }
325 
326 //_____________________________________________________________________
329 {
330  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::SetEMCalGeometry()");
331  if(fGeometryName.Length()==0){
332  fgeom=AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
333  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fgeom->GetName(),fRunNumber));
334  } else {
335  fgeom = AliEMCALGeometry::GetInstance(fGeometryName.Data());
336  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeometryName.Data()));
337  }
338 
339  if (!fgeom){
340  AliWarning("Make sure the EMCal geometry is set properly !");
341  } else {
342  AliDebug(1,Form("EMCal geometry properly set: fGeom = %p, fGeometryName=%s",fgeom,fGeometryName.Data()));
343  }
344 
345  return kTRUE;
346 }
347 
348 //_____________________________________________________________________
351 {
352  //method under development
353  AliInfo(Form("<D> -- Run # = %d", fRunNumber));
354  AliInfo("prepare TOFT0maker!!");
355  //cout<<"Run "<< fRunNumber<<" in TOFT0maker"<<endl;
356 
357 
358  AliCDBManager * cdb = AliCDBManager::Instance();
359  cdb->SetDefaultStorage("raw://");
360  cdb->SetRun(fRunNumber);
361 
362  AliESDpid *extPID=new AliESDpid();
363 
364  // Wonder if some have to be declared as private variables??
365  // AliESDpid *extPID = new AliESDpid();
366  // AliTOFcalib * tofCalib = new AliTOFcalib();
367  // tofCalib->SetCalibrateTOFsignal(kTRUE);
368  // tofCalib->Init();
369 
370  fTOFmaker = new AliTOFT0maker(extPID);
371  fTOFmaker->SetTimeResolution(115.0); // if you want set the TOF res
372  // fTOFmaker = new AliTOFT0maker(extPID,tofCalib);
373  // fTOFmaker->SetTimeResolution(130.0);
374 
375  //cout<<"extPID "<<extPID<<" fTOFmaker "<<fTOFmaker<<endl;
376 
377 }// End PrepareTOFT0maker
378 
379 //________________________________________________________________________
383 {
384  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::UserCreateOutputObjects()");
385 
386  const Int_t nChannels = 17664;
387  //book histograms
388  if(fFillHeavyHisto){
389  fhcalcEvtTime = new TH1F("fhcalcEvtTime","calculated event time from T0",fFineNbins, fFineNbins,fFineTmax);
390  fhcalcEvtTime->GetXaxis()->SetTitle("T ");
391  fhcalcEvtTime->GetYaxis()->SetTitle("Counts (a.u.)");
392 
393  fhEvtTimeHeader = new TH1F("fhEvtTimeHeader","event time from header",fFineNbins, fFineNbins,fFineTmax);
394  fhEvtTimeHeader->GetXaxis()->SetTitle("T ");
395  fhEvtTimeHeader->GetYaxis()->SetTitle("Counts (a.u.)");
396 
397  fhEvtTimeDiff = new TH1F("fhEvtTimeDiff","event time difference",fFineNbins, fFineNbins,fFineTmax);
398  fhEvtTimeDiff->GetXaxis()->SetTitle("#Delta T ");
399  fhEvtTimeDiff->GetYaxis()->SetTitle("Counts (a.u.)");
400  }
401 
402  fhEventType = new TH1F("fhEventType","event type",10, 0.,10.);
403  fhEventType ->GetXaxis()->SetTitle("Type ");
404  fhEventType ->GetYaxis()->SetTitle("Counts (a.u.)");
405  if(fFillHeavyHisto){
406  fhTcellvsTOFT0 = new TH2F("hTcellvsTOFT0", " T_cell vs TOFT0", 500,-600.0,+400.0,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
407  fhTcellvsTOFT0HD = new TH2F("hTcellvsTOFT0HD", " T_cell vs TOFT0,HighEnergy", 500,-600.0,+400.0,4*fRawTimeNbins,fRawTimeMin,fRawTimeMax);
408  }
409  fhTcellvsSM = new TH2F("hTcellvsSM", " T_cell vs SM", (Int_t)kNSM,0,(Double_t)kNSM,(Int_t)(fRawTimeNbins/2),fRawTimeMin,fRawTimeMax);
410 
411  if(fFillHeavyHisto){
412  fhEneVsAbsIdHG = new TH2F("fhEneVsAbsIdHG", "energy vs ID for HG",1000,0,18000,200,0,10);
413  fhEneVsAbsIdLG = new TH2F("fhEneVsAbsIdLG", "energy vs ID for LG",1000,0,18000,200,0,40);
414  }
415 
416  for (Int_t i = 0; i < kNBCmask ; i++)
417  {
418  //already after correction
419  //high gain
420  fhTimeSumSq[i] = new TH1F(Form("hTimeSumSq%d", i),
421  Form("cell Sum Square time HG, BC %d ", i),
422  nChannels,0.,(Double_t)nChannels);
423  fhTimeSumSq[i]->SetYTitle("Sum Sq Time ");
424  fhTimeSumSq[i]->SetXTitle("AbsId");
425 
426  fhTimeSum[i] = new TH1F(Form("hTimeSum%d", i),
427  Form("cell Sum time HG, BC %d ", i),
428  nChannels,0.,(Double_t)nChannels);
429  fhTimeSum[i]->SetYTitle("Sum Time ");
430  fhTimeSum[i]->SetXTitle("AbsId");
431 
432  fhTimeEnt[i] = new TH1F(Form("hTimeEnt%d", i),
433  Form("cell Entries HG, BC %d ", i),
434  nChannels,0.,(Double_t)nChannels);
435  fhTimeEnt[i]->SetYTitle("Entries for Time ");
436  fhTimeEnt[i]->SetXTitle("AbsId");
437 
438  //low gain
439  fhTimeLGSumSq[i] = new TH1F(Form("hTimeLGSumSq%d", i),
440  Form("cell Sum Square time LG, BC %d ", i),
441  nChannels,0.,(Double_t)nChannels);
442  fhTimeLGSumSq[i]->SetYTitle("Sum Sq Time ");
443  fhTimeLGSumSq[i]->SetXTitle("AbsId");
444 
445  fhTimeLGSum[i] = new TH1F(Form("hTimeLGSum%d", i),
446  Form("cell Sum time LG, BC %d ", i),
447  nChannels,0.,(Double_t)nChannels);
448  fhTimeLGSum[i]->SetYTitle("Sum Time ");
449  fhTimeLGSum[i]->SetXTitle("AbsId");
450 
451  fhTimeLGEnt[i] = new TH1F(Form("hTimeLGEnt%d", i),
452  Form("cell Entries LG, BC %d ", i),
453  nChannels,0.,(Double_t)nChannels);
454  fhTimeLGEnt[i]->SetYTitle("Entries for Time ");
455  fhTimeLGEnt[i]->SetXTitle("AbsId");
456 
457  //raw time histograms
458  //high gain
459  if(fFillHeavyHisto){
460  fhRawTimeVsIdBC[i] = new TH2F(Form("RawTimeVsIdBC%d", i),
461  Form("cell raw time vs ID for high gain BC %d ", i),
462  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
463  fhRawTimeVsIdBC[i]->SetXTitle("AbsId");
464  fhRawTimeVsIdBC[i]->SetYTitle("Time");
465  }
466 
467  fhRawTimeSumBC[i] = new TH1F(Form("RawTimeSumBC%d", i),
468  Form("sum of cell raw time for high gain BC %d ", i),
469  nChannels,0.,(Double_t)nChannels);
470  fhRawTimeSumBC[i]->SetXTitle("AbsId");
471  fhRawTimeSumBC[i]->SetYTitle("Sum Time");
472 
473  fhRawTimeEntriesBC[i] = new TH1F(Form("RawTimeEntriesBC%d", i),
474  Form("No. entries of cells raw time for high gain BC %d ", i),
475  nChannels,0.,(Double_t)nChannels);
476  fhRawTimeEntriesBC[i]->SetXTitle("AbsId");
477  fhRawTimeEntriesBC[i]->SetYTitle("Entries for Time ");
478 
479  fhRawTimeSumSqBC[i] = new TH1F(Form("RawTimeSumSqBC%d", i),
480  Form("sum of (cell raw time)^2 for high gain BC %d ", i),
481  nChannels,0.,(Double_t)nChannels);
482  fhRawTimeSumSqBC[i]->SetXTitle("AbsId");
483  fhRawTimeSumSqBC[i]->SetYTitle("Sum Sq Time");
484 
485  //low gain
486  if(fFillHeavyHisto){
487  fhRawTimeVsIdLGBC[i] = new TH2F(Form("RawTimeVsIdLGBC%d", i),
488  Form("cell raw time vs ID for low gain BC %d ", i),
489  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
490  fhRawTimeVsIdLGBC[i]->SetXTitle("AbsId");
491  fhRawTimeVsIdLGBC[i]->SetYTitle("Time");
492  }
493 
494  fhRawTimeSumLGBC[i] = new TH1F(Form("RawTimeSumLGBC%d", i),
495  Form("sum of cell raw time for low gain BC %d ", i),
496  nChannels,0.,(Double_t)nChannels);
497  fhRawTimeSumLGBC[i]->SetXTitle("AbsId");
498  fhRawTimeSumLGBC[i]->SetYTitle("Sum Time");
499 
500  fhRawTimeEntriesLGBC[i] = new TH1F(Form("RawTimeEntriesLGBC%d", i),
501  Form("No. entries of cells raw time for low gain BC %d ", i),
502  nChannels,0.,(Double_t)nChannels);
503  fhRawTimeEntriesLGBC[i]->SetXTitle("AbsId");
504  fhRawTimeEntriesLGBC[i]->SetYTitle("Entries for Time ");
505 
506  fhRawTimeSumSqLGBC[i] = new TH1F(Form("RawTimeSumSqLGBC%d", i),
507  Form("sum of (cell raw time)^2 for low gain BC %d ", i),
508  nChannels,0.,(Double_t)nChannels);
509  fhRawTimeSumSqLGBC[i]->SetXTitle("AbsId");
510  fhRawTimeSumSqLGBC[i]->SetYTitle("Sum Sq Time");
511 
512  //histograms with corrected raw time for L1 shift and 100ns
513  if(fBadReco && fFillHeavyHisto){
514  fhRawCorrTimeVsIdBC[i] = new TH2F(Form("RawCorrTimeVsIdBC%d", i),
515  Form("cell L1 shift and 100ns corrected raw time vs ID for high gain BC %d ", i),
516  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
517  fhRawCorrTimeVsIdBC[i]->SetXTitle("AbsId");
518  fhRawCorrTimeVsIdBC[i]->SetYTitle("Time");
519 
520  fhRawCorrTimeVsIdLGBC[i] = new TH2F(Form("RawCorrTimeVsIdLGBC%d", i),
521  Form("cell L1 shift and 100ns corrected raw time vs ID for low gain BC %d ", i),
522  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
523  fhRawCorrTimeVsIdLGBC[i]->SetXTitle("AbsId");
524  fhRawCorrTimeVsIdLGBC[i]->SetYTitle("Time");
525  }
526 
527  //histograms with corrected raw time for L1 shift and 100ns + new L1 phase
528  if(fReferenceRunByRunFileName.Length()!=0 && fFillHeavyHisto){
529  fhTimeVsIdBC[i] = new TH2F(Form("TimeVsIdBC%d", i),
530  Form("cell time corrected for L1 shift, 100ns and L1 phase vs ID for high gain BC %d ", i),
531  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
532  fhTimeVsIdBC[i]->SetXTitle("AbsId");
533  fhTimeVsIdBC[i]->SetYTitle("Time");
534 
535  fhTimeVsIdLGBC[i] = new TH2F(Form("TimeVsIdLGBC%d", i),
536  Form("cell time corrected for L1 shift, 100ns and L1 phase vs ID for low gain BC %d ", i),
537  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
538  fhTimeVsIdLGBC[i]->SetXTitle("AbsId");
539  fhTimeVsIdLGBC[i]->SetYTitle("Time");
540  }
541 
542  for (Int_t j = 0; j < kNSM ; j++)
543  {
544  //High gain
545  //fhTimeDsupBC[j][i]= new TH2F(Form("SupMod%dBC%d",j,i), Form("SupMod %d time_vs_E BC %d",j,i),500,0.0,20.0,2200,-350.0,750.0);
546  fhTimeDsupBC[j][i]= new TH2F(Form("SupMod%dBC%d",j,i), Form("SupMod %d time_vs_E, high gain, BC %d",j,i),fEnergyNbins,fEnergyMin,fEnergyMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
547  fhTimeDsupBC[j][i]->SetYTitle(" Time (ns) ");
548  fhTimeDsupBC[j][i]->SetXTitle(" E (GeV) ");
549 
550  //low gain
551  fhTimeDsupLGBC[j][i]= new TH2F(Form("SupMod%dBC%dLG",j,i), Form("SupMod %d time_vs_E, low gain, BC %d",j,i),fEnergyLGNbins,fEnergyLGMin,fEnergyLGMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
552  fhTimeDsupLGBC[j][i]->SetYTitle(" Time (ns) ");
553  fhTimeDsupLGBC[j][i]->SetXTitle(" E (GeV) ");
554  }
555  }
556 
557  for (Int_t jj = 0; jj < kNSM ; jj++)
558  {
559  //high gain
560  fhTimeDsup[jj] = new TH2F(Form("SupMod%d",jj), Form("SupMod %d time_vs_E, high gain",jj),fEnergyNbins,fEnergyMin,fEnergyMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
561  fhTimeDsup[jj]->SetYTitle(" Time (ns) ");
562  fhTimeDsup[jj]->SetXTitle(" E (GeV) ");
563 
564  //low gain
565  fhTimeDsupLG[jj] = new TH2F(Form("SupMod%dLG",jj), Form("SupMod %d time_vs_E, low gain ",jj),fEnergyLGNbins,fEnergyLGMin,fEnergyLGMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
566  fhTimeDsupLG[jj]->SetYTitle(" Time (ns) ");
567  fhTimeDsupLG[jj]->SetXTitle(" E (GeV) ");
568  }
569 
570  fhTimeVsBC = new TH2F("TimeVsBC"," SupMod time_vs_BC ", 4001,-0.5,4000.5,(Int_t)(fRawTimeNbins/2.),fRawTimeMin,fRawTimeMax);
571 
572 
573  //add histos to list
574  fOutputList = new TList();
575  fOutputList->Add(fhEventType);
576  if(fFillHeavyHisto){
580 
585  }
586  fOutputList->Add(fhTcellvsSM);
587 
588  for (Int_t i = 0; i < kNBCmask ; i++)
589  {
590  fOutputList->Add(fhTimeSumSq[i]);
591  fOutputList->Add(fhTimeEnt[i]);
592  fOutputList->Add(fhTimeSum[i]);
593 
594  fOutputList->Add(fhTimeLGSumSq[i]);
595  fOutputList->Add(fhTimeLGEnt[i]);
596  fOutputList->Add(fhTimeLGSum[i]);
597 
598  if(fFillHeavyHisto) {
599  fOutputList->Add(fhRawTimeVsIdBC[i]);
601  }
602 
603  fOutputList->Add(fhRawTimeSumBC[i]);
605  fOutputList->Add(fhRawTimeSumSqBC[i]);
606 
607  fOutputList->Add(fhRawTimeSumLGBC[i]);
610 
611  if(fBadReco && fFillHeavyHisto) {
614  }
615  if(fReferenceRunByRunFileName.Length()!=0 && fFillHeavyHisto) {
616  fOutputList->Add(fhTimeVsIdBC[i]);
617  fOutputList->Add(fhTimeVsIdLGBC[i]);
618  }
619 
620  for (Int_t j = 0; j < kNSM ; j++){
621  fOutputList->Add(fhTimeDsupBC[j][i]);
622  fOutputList->Add(fhTimeDsupLGBC[j][i]);
623  }
624  }
625 
626  for (Int_t j = 0; j < kNSM ; j++)
627  {
628  fOutputList->Add(fhTimeDsup[j]);
629  fOutputList->Add(fhTimeDsupLG[j]);
630  }
631 
632  fOutputList->Add(fhTimeVsBC);
633 
634  fOutputList->SetOwner(kTRUE);
635  PostData(1,fOutputList);
636 
637 
638 } // End of AliAnalysisTaskEMCALTimeCalib::UserCreateOuputObjects()
639 
640 //________________________________________________________________________
643 {
644  // Called for each event
645  AliDebug(2,Form("UserExec: EMCal geometry: fgeom = %p fGeometryName %s",fgeom,fGeometryName.Data()));
646  AliVEvent *event = InputEvent();
647  //cout<<"T0TOF "<<event->GetT0TOF()<<endl;//bad idea
648  //cout<< fEvent->GetTOFHeader()->GetDefaultEventTimeVal()<<endl;
649  AliDebug(2,Form("TOF time from header %f ps",event->GetTOFHeader()->GetDefaultEventTimeVal()));
650  if(fFillHeavyHisto) fhEvtTimeHeader->Fill(event->GetTOFHeader()->GetDefaultEventTimeVal());
651 
652  //fEvent = dynamic_cast<AliESDEvent*>(event);
653  if (!event) {
654  AliError("ESD not available, exit");
655  fhEventType->Fill(0.5);
656  return;
657  }
658 
659  if(fPileupFromSPD==kTRUE){
660  if(event->IsPileupFromSPD(3,0.8,3.,2.,5.)){
661  AliDebug(1,"Event: PileUp skip.");
662  fhEventType->Fill(1.5);
663  return;
664  }
665  }
666 
667  TString triggerclasses = event->GetFiredTriggerClasses();
668  if(triggerclasses=="") {
669  fhEventType->Fill(2.5);
670  return;
671  }
672 
673  Int_t eventType = ((AliVHeader*)event->GetHeader())->GetEventType();
674  // physics events eventType=7, select only those
675  AliDebug(1,Form("Triggerclasses %s, eventType %d",triggerclasses.Data(),eventType));
676  if(eventType != 7) {
677  fhEventType->Fill(3.5);
678  return;
679  }
680 
681  // Check trigger
682  Bool_t bMB = kFALSE;
683  Bool_t bL0 = kFALSE;
684  Bool_t bL1G = kFALSE;
685  Bool_t bL1J = kFALSE;
686 
687  if(triggerclasses.Contains("CINT7-B-NOPF-ALLNOTRD") ||
688  triggerclasses.Contains("CINT7-I-NOPF-ALLNOTRD") ||
689  triggerclasses.Contains("CINT1-I-NOPF-ALLNOTRD") ||
690  triggerclasses.Contains("CINT1-B-NOPF-ALLNOTRD") ||
691  triggerclasses.Contains("CINT8") ||
692  triggerclasses.Contains("CINT7") ||
693  triggerclasses.Contains("CPBI2_B1-B-NOPF-ALLNOTRD") ) bMB = kTRUE;
694 
695  if(triggerclasses.Contains("CEMC7-B-NOPF-CENTNOTRD") ||
696  triggerclasses.Contains("CEMC1-B-NOPF-CENTNOTRD") ||
697  triggerclasses.Contains("CEMC7") ||
698  triggerclasses.Contains("CEMC8") ||
699  triggerclasses.Contains("CEMC8-B-NOPF-CENTNOTRD") ) bL0 = kTRUE;
700 
701  if(triggerclasses.Contains("CEMC7EG1-B-NOPF-CENTNOTRD") ||
702  triggerclasses.Contains("CEMC7EG2-B-NOPF-CENTNOTRD") ||
703  triggerclasses.Contains("CEMC8EG1-B-NOPF-CENTNOTRD") ||
704  triggerclasses.Contains("CEMC8EGA") ||
705  triggerclasses.Contains("CEMC7EGA") ||
706  triggerclasses.Contains("CPBI2EGA") ) bL1G = kTRUE;
707 
708 
709  if(triggerclasses.Contains("CEMC7EJ1-B-NOPF-CENTNOTRD") ||
710  triggerclasses.Contains("CEMC7EJ2-B-NOPF-CENTNOTRD") ||
711  triggerclasses.Contains("CEMC8EJ1-B-NOPF-CENTNOTRD") ||
712  triggerclasses.Contains("CEMC7EJE") ||
713  triggerclasses.Contains("CEMC8EJE") ||
714  triggerclasses.Contains("CPBI2EJE") ) bL1J = kTRUE;
715 
716  if( bL1G || bL1J || bL0 ){ fhEventType->Fill(4.5);}
717  if( bMB ){ fhEventType->Fill(5.5);}
718 
719 
720  // if(bL1G || bL1J || bL0){
721 
722 // Prepare TOFT0 maker at the beginning of a run
723 // if (event->GetRunNumber() != fRunNumber){
724 // AliInfo(Form("Runno per event %d",event->GetRunNumber()));
725 // fRunNumber = event->GetRunNumber();
726 // // PrepareTOFT0maker();
727 // // cout<<"tofT0maker per run"<<fRunNumber<<endl;
728 // }// fi Check if run number has changed
729 
730  // --- Use of AliTOFT0maker
731  Double_t calcolot0=0.0;
732  if(!AODEvent()){
733  Double_t* timeTOFtable;
734  timeTOFtable=fTOFmaker->ComputeT0TOF(dynamic_cast<AliESDEvent*>(event));
735  AliDebug(2,Form("TOF time %f ps, resolution %f ps, tracks at TOF %f/used %f",timeTOFtable[0],timeTOFtable[1],timeTOFtable[3],timeTOFtable[7]));
736  //cout<<"event time "<<timeTOFtable[0]<<" resolution "<<timeTOFtable[1]<<"ps av. ev. time "<<timeTOFtable[2]<<" trks at TOF "<<timeTOFtable[3]<<" calc evnt time "<<timeTOFtable[4]<<" resolution "<<timeTOFtable[5]<<" tracks used "<<timeTOFtable[7]<<endl;
737  calcolot0=timeTOFtable[0];
738  }
739 
740  if(fFillHeavyHisto) {
741  fhcalcEvtTime->Fill(calcolot0);
742  if(calcolot0 != 0 && event->GetTOFHeader()->GetDefaultEventTimeVal() != 0 )
743  fhEvtTimeDiff->Fill(calcolot0-event->GetTOFHeader()->GetDefaultEventTimeVal());
744  }
745 
746  TRefArray* caloClusters = new TRefArray();
747  event->GetEMCALClusters(caloClusters);
748  // cout << " ###########Bunch Cross nb = " << event->GetBunchCrossNumber() << endl;
749 
750  Int_t BunchCrossNumber =event->GetBunchCrossNumber();
751 
752  Float_t offset=0.;
753  Float_t offsetPerSM=0.;
754  Int_t L1phaseshift=0;
755  Int_t L1phase=0;
756  Int_t L1shiftOffset=0;
757 
758  Int_t nBC = 0;
759  nBC = BunchCrossNumber%4;
760  //Int_t nTriggerMask =event->GetTriggerMask();
761  // cout << " nBC " << nBC << " nTriggerMask " << nTriggerMask<< endl;
762  Float_t timeBCoffset = 0.; //manual offset
763  // if( nBC%4 ==0 || nBC%4==1) timeBCoffset = 100.; // correction was for LHC11 when BC was not corrected
764 
765  Int_t nclus = caloClusters->GetEntries();
766  AliDebug(1,Form("###########Bunch Cross nb = %d nclus = %d",nBC,nclus ));
767  //cout << " ###########Bunch Cross nb = " << nBC <<" nclus= "<< nclus<< endl;
768  //Int_t ntracks = event-> GetNumberOfTracks() ;
769 
770  AliVCaloCells &cells= *(event->GetEMCALCells());//it is cluster independent
771  //Variables used plenty in loops
772  Int_t nSupMod=-1, nModule=-1;
773  Int_t iphi=-1, ieta=-1, nIphi=-1, nIeta=-1;
774  Int_t absId=-1;
775  Float_t hkdtime=0.0;
776  Float_t amp=0.0;
777  Bool_t isHighGain=kTRUE;
778 
779  for (Int_t icl = 0; icl < nclus; icl++) {
780  //ESD and AOD CaloCells carries the same information
781  AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
782  if(!AcceptCluster(clus)) continue;
783 
784  //cout<<"nCells="<< clus->GetNCells();<<endl;
785 
786  UShort_t * index = clus->GetCellsAbsId() ;
787 
788  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
789  absId = index[i]; // or clus->GetCellNumber(i) ;
790  hkdtime = cells.GetCellTime(absId) * 1.0e09; // to get ns
791  amp = cells.GetCellAmplitude(absId) ;
792  isHighGain = cells.GetCellHighGain(absId);
793  //cout<<"cell absID: "<<absId<<" cellTime: "<<hkdtime<<" cellaplit: "<< amp<<endl;
794  // GEOMETRY tranformations
795  fgeom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
796  fgeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
797 
798  //bad channel check. 0: good channel, 1-5: bad channel
800  if(GetEMCALChannelStatus(nSupMod,ieta,iphi)) continue;//printf("bad\n");
801  } else if(fSetBadChannelMapSource==2){
802  if(GetEMCALChannelStatus(absId)) continue;//printf("bad\n");
803  }
804 
805  //main histograms with raw time information
806  if(amp>fMinCellEnergy){
807  if(isHighGain){
808  if(fFillHeavyHisto) fhRawTimeVsIdBC[nBC]->Fill(absId,hkdtime);
809  fhRawTimeSumBC[nBC]->Fill(absId,hkdtime);
810  fhRawTimeEntriesBC[nBC]->Fill(absId,1.);
811  fhRawTimeSumSqBC[nBC]->Fill(absId,hkdtime*hkdtime);
812  }else{
813  if(fFillHeavyHisto) fhRawTimeVsIdLGBC[nBC]->Fill(absId,hkdtime);
814  fhRawTimeSumLGBC[nBC]->Fill(absId,hkdtime);
815  fhRawTimeEntriesLGBC[nBC]->Fill(absId,1.);
816  fhRawTimeSumSqLGBC[nBC]->Fill(absId,hkdtime*hkdtime);
817  }
818  }
819  //fgeom->PrintCellIndexes(absId);
820  //fgeom->PrintCellIndexes(absId,1);
821 
822  // other histograms for cross-check
823  CheckCellRCU(nSupMod,ieta,iphi);//SM, column, row
824 
825  fhTcellvsSM->Fill(nSupMod,hkdtime);
826  if(fFillHeavyHisto) {
827  if(isHighGain==kTRUE) {fhEneVsAbsIdHG->Fill(absId,amp);}
828  else {fhEneVsAbsIdLG->Fill(absId,amp);}
829  }
830  fhTimeVsBC->Fill(1.*BunchCrossNumber,hkdtime-timeBCoffset);
831  //important remark: We use 'Underflow bin' for absid=0 in OADB for time calibration
832  if(isHighGain==kTRUE){
833  if(fhAllAverageBC[nBC]!=0) {//comming from file after the first iteration
834  offset = (Float_t)(fhAllAverageBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
835  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
836  AliFatal(Form("Reference histogram for BC%d not properly loaded",nBC));
837  }
838  } else {
839  if(fhAllAverageLGBC[nBC]!=0) {//comming from file after the first iteration
840  offset = (Float_t)(fhAllAverageLGBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
841  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
842  AliFatal(Form("Reference LG histogram for BC%d not properly loaded",nBC));
843  }
844  }
845  //if(offset==0)cout<<"offset 0 in SM "<<nSupMod<<endl;
846 
847  // Solution for 2015 data where L1 phase and L1 shift is not correct in data. We need to calibrate run by run.
848  // The shift and phase are done per SM (0-19).
849  // L1 phase is necessary in run 2.
850  // L1 shift is necessary only for bad reconstructed runs (muon_calo_pass1 lhc15f-m)
851  if(fhRefRuns!=0) {//comming from file after the first iteration
852  L1phaseshift = (Int_t)(fhRefRuns->GetBinContent(nSupMod));//SM0 = bin0
853 
854  // to correct for L1 phase
855  // this part works for both: muon_calo_pass1 of LHC15n (pp@2.76) and later reconstructions
856  // wrong reconstruction done before in run2
857  L1phase = L1phaseshift & 3; //bit operation
858  if(nBC >= L1phase)
859  offsetPerSM = (nBC - L1phase)*25;
860  else
861  offsetPerSM = (nBC - L1phase + 4)*25;
862 
863  // to correct for L1 shift
864  // this part is only for wrongly reconstructed runs before LHC15n in run2
865  if(fBadReco){
866  L1shiftOffset = L1phaseshift>>2; //bit operation
867  L1shiftOffset*=25;
868  //(we subtract it here because we subtract the whole wrong offset later --=+)
869  if(nBC==0 || nBC==1) L1shiftOffset-=100.;//additional shift for muon_calo_pass1 up to lhc15f-m
870  }
871  } else if(fReferenceRunByRunFileName.Length()!=0){//protection against missing reference histogram
872  AliFatal("Reference histogram run-by-run not properly loaded");
873  }
874  //end of load additional offset
875 
876  //fill the raw time with L1 shift correction and 100ns
878  if(isHighGain){
879  fhRawCorrTimeVsIdBC[nBC]->Fill(absId,hkdtime-L1shiftOffset);
880  }else{
881  fhRawCorrTimeVsIdLGBC[nBC]->Fill(absId,hkdtime-L1shiftOffset);
882  }
883  }
884 
885  //fill time after L1 shift correction and 100ns and new L1 phase
887  if(isHighGain){
888  fhTimeVsIdBC[nBC]->Fill(absId,hkdtime-L1shiftOffset-offsetPerSM);
889  }else{
890  fhTimeVsIdLGBC[nBC]->Fill(absId,hkdtime-L1shiftOffset-offsetPerSM);
891  }
892  }
893 
894  //other control histograms
895  if(amp>0.5) {
896  if(isHighGain){
897  fhTimeDsup[nSupMod]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
898  fhTimeDsupBC[nSupMod][nBC]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
899  }else{
900  fhTimeDsupLG[nSupMod]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
901  fhTimeDsupLGBC[nSupMod][nBC]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
902  }
903  }
904 
905  if(fFillHeavyHisto) {
906  if(amp>0.9) {
907  fhTcellvsTOFT0HD->Fill(calcolot0, hkdtime);
908  }
909  fhTcellvsTOFT0->Fill(calcolot0, hkdtime-offset-offsetPerSM-L1shiftOffset);
910  }
911 
912  hkdtime = hkdtime-timeBCoffset;//time corrected by manual offset (default=0)
913  Float_t hkdtimecorr;
914  hkdtimecorr= hkdtime-offset-offsetPerSM-L1shiftOffset;//time after first iteration
915 
916  //main histograms after the first itereation for calibration constants
917  //if(hkdtimecorr>=-20. && hkdtimecorr<=20. && amp>0.9 ) {
918  if(hkdtimecorr>=fMinTime && hkdtimecorr<=fMaxTime && amp>fMinCellEnergy ) {
919  // per cell
920 // Float_t entriesTime=fhTimeEnt[nBC]->GetBinContent(absId)+1;
921 // Float_t sumTimeSq=(fhTimeSumSq[nBC]->GetBinContent(absId)+(hkdtime*hkdtime));
922 // Float_t sumTime=(fhTimeSum[nBC]->GetBinContent(absId)+hkdtime);
923 //
924 // fhTimeEnt[nBC]->SetBinContent(absId,entriesTime);
925 // fhTimeSumSq[nBC]->SetBinContent(absId,sumTimeSq);
926 // fhTimeSum[nBC]->SetBinContent(absId,sumTime);
927 
928  //correction in 2015 for wrong L1 phase and L1 shift
929  hkdtime = hkdtime - offsetPerSM - L1shiftOffset;
930 
931  if(isHighGain){
932  fhTimeEnt[nBC]->Fill(absId,1.);
933  fhTimeSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
934  fhTimeSum[nBC]->Fill(absId,hkdtime);
935  }else{
936  fhTimeLGEnt[nBC]->Fill(absId,1.);
937  fhTimeLGSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
938  fhTimeLGSum[nBC]->Fill(absId,hkdtime);
939  }
940 
941 
942  } // hkdtime:[-20;20]
943  } // end icell
944  } //end cluster
945 
946 
947  // Post output data.
948  //cout<<"Post data and delete caloClusters"<<endl;
949  caloClusters->Delete();
950  delete caloClusters;
951 // } // end if trigger type
952 
953  PostData(1, fOutputList);
954 } // End of AliAnalysisTaskEMCALTimeCalib::UserExec()
955 
956 //________________________________________________________________________
960 {
961  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
962 
963  if(fTOFmaker) delete fTOFmaker;
964 
965  if(fL1PhaseList) {
966  fL1PhaseList->SetOwner();
967  fL1PhaseList->Clear();
968  delete fL1PhaseList;
969  }
970 
971  if (fBadChannelMapArray) {
972  fBadChannelMapArray->Clear();
973  delete fBadChannelMapArray;
974  }
975 
976  if (!fOutputList)
977  {
978  AliDebug(1,"ERROR: Output list not available");
979  return;
980  }
981 } // End of AliAnalysisTaskEMCALTimeCalib::Terminate
982 
983 //________________________________________________________________________
986 {
987  //fix with noisy EMCAL fee card
988  Int_t nCells = clus->GetNCells();
989 
990  if(clus->IsEMCAL())
991  {
992  if ((clus->E() > fMaxClusterEnergy && nCells > fMaxNcells ) || nCells > fMaxNcells){
993  AliDebug(1,"very big cluster with enormous energy - cluster rejected");
994  return kFALSE;
995  }
996  }
997 
998  // remove other than photonlike
999  Double_t lambda0=clus->GetM02();
1000  if (lambda0>fMaxLambda0LG || lambda0<fMinLambda0LG){
1001  AliDebug(1,"lambda0 loose cut failed - cluster rejected");
1002  return kFALSE;
1003  }
1004 
1005  // remove matched clusters
1006  Double_t Dx=clus->GetTrackDx();
1007  Double_t Dz=clus->GetTrackDz();
1008  Double_t Rtrack = TMath::Sqrt(Dx*Dx+Dz*Dz);
1009  if (Rtrack <fMaxRtrack)
1010  {
1011  AliDebug(1,"track matched - cluster rejected");
1012  return kFALSE;
1013  }
1014 
1015  if (nCells<fMinNcells)
1016  {
1017  AliDebug(1,"single cell cluster - cluster rejected");
1018  return kFALSE;
1019  }
1020 
1021  if(clus->E()<fMinClusterEnergy)
1022  {
1023  AliDebug(1,"cluster energy < 1 GeV- cluster rejected");
1024  return kFALSE;
1025  }
1026 
1027 
1028  if(!IsLowGainCellInCluster(clus)) {//no low gain cell in cluster
1029  //apply more strict lambda0^2 cut
1030  if (lambda0>fMaxLambda0 || lambda0<fMinLambda0){
1031  AliDebug(1,"lambda0 strict cut failed - cluster rejected");
1032  return kFALSE;
1033  }
1034  }
1035 
1036 
1037 
1038 
1039  return kTRUE;
1040 }//End AliAnalysisTaskEMCALTimeCalib::AcceptCluster
1041 
1042 //________________________________________________________________________
1045  UShort_t * index = clus->GetCellsAbsId() ;
1046  AliVCaloCells &cells= *(InputEvent()->GetEMCALCells());
1047  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
1048  if(cells.GetCellHighGain(index[i])==kFALSE) return kTRUE;//low gain cell
1049  }
1050  return kFALSE;
1051 
1052 }
1053 
1054 //________________________________________________________________________
1057 {
1058  Int_t iRCU;
1059  if(nSupMod < 10 || (nSupMod >= 12 && nSupMod <18) )
1060  {
1061  if (0<=irow&&irow<8) iRCU=0; // first cable row
1062  else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
1063  //second cable row
1064  //RCU1
1065  else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
1066  //second cable row
1067  else if (16<=irow&&irow<24) iRCU=1; // third cable row
1068 
1069  if (nSupMod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1070  }
1071  else
1072  {
1073  // Last 2 SM have one single SRU, just assign RCU 0
1074  iRCU = 0 ;
1075  }
1076 
1077  //cout<<"RCU:"<<iRCU<<endl;
1078  if (iRCU<0)
1079  AliFatal(Form("Wrong EMCAL/DCAL RCU number = %d\n", iRCU));
1080 
1081  return kTRUE;
1082 }//End AliAnalysisTaskEMCALTimeCalib::CheckCellRCU
1083 
1084 //________________________________________________________________________
1087 {
1088  fMinClusterEnergy=1.0;//0.5//0.7
1089  fMaxClusterEnergy=500;
1090  fMinNcells=2;
1091  fMaxNcells=200;
1092  fMinLambda0=0.1;
1093  fMaxLambda0=0.4;
1094  fMinLambda0LG=0.1;
1095  fMaxLambda0LG=4.0;
1096  fMaxRtrack=0.025;
1097  fMinCellEnergy=0.4;//0.1//0.4
1098  fReferenceFileName="";//Reference.root
1100  fBadReco=kFALSE;
1101  fFillHeavyHisto=kFALSE;
1102  fGeometryName="";//EMCAL_COMPLETE12SMV1_DCAL_8SM
1103  fPileupFromSPD=kFALSE;
1104  fMinTime=-20.;
1105  fMaxTime=20.;
1106 
1107  fBadChannelMapSet=kFALSE;
1110 
1111  //histograms
1112  fRawTimeNbins = 400; // Raw time settings should be like that all the time
1113  fRawTimeMin = 400.; // importent in pass1
1114  fRawTimeMax = 800.;
1115  fPassTimeNbins = 1000; // in pass2 should be (400,400,800)
1116  fPassTimeMin = -250.;// in pass3 should be (1000,-250,250)
1117  fPassTimeMax = 250.;
1118  fEnergyNbins = 100; // default settings was 500
1119  fEnergyMin = 0.;
1120  fEnergyMax = 20.;
1121  fEnergyLGNbins = 200; // default settings
1122  fEnergyLGMin = 0.;
1123  fEnergyLGMax = 100.;
1124  fFineNbins = 90; //was 4500 for T0 time studies
1125  fFineTmin = -500;
1126  fFineTmax = 400;
1127 }
1128 
1129 //________________________________________________________________________
1135 {
1136  TFile *file =new TFile(inputFile.Data());
1137  if(file==0x0) {
1138  //AliWarning("Input file does not exist!");
1139  return;
1140  }
1141 
1142  TList *list=(TList*)file->Get("chistolist");
1143  if(list==0x0)
1144  {
1145  //AliWarning("List chistolist does not exist in file!");
1146  return;
1147  }
1148 
1149  //high gain
1150  TH1F *h1[4];
1151  TH1F *h2[4];
1152  TH1F *h3[4];
1153  TH1F *hAllTimeAvBC[4];
1154  TH1F *hAllTimeRMSBC[4];
1155 
1156  //low gain
1157  TH1F *h4[4];
1158  TH1F *h5[4];
1159  TH1F *h6[4];
1160  TH1F *hAllTimeAvLGBC[4];
1161  TH1F *hAllTimeRMSLGBC[4];
1162 
1163  if(isFinal==kFALSE){//first itereation
1164  for(Int_t i=0;i<4;i++){
1165  h1[i]=(TH1F *)list->FindObject(Form("RawTimeSumBC%d",i));
1166  h2[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesBC%d",i));
1167  h3[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqBC%d",i));
1168 
1169  h4[i]=(TH1F *)list->FindObject(Form("RawTimeSumLGBC%d",i));
1170  h5[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesLGBC%d",i));
1171  h6[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqLGBC%d",i));
1172  }
1173  } else {//final iteration
1174  for(Int_t i=0;i<4;i++){
1175  h1[i]=(TH1F *)list->FindObject(Form("hTimeSum%d",i));
1176  h2[i]=(TH1F *)list->FindObject(Form("hTimeEnt%d",i));
1177  h3[i]=(TH1F *)list->FindObject(Form("hTimeSumSq%d",i));
1178 
1179  h4[i]=(TH1F *)list->FindObject(Form("hTimeLGSum%d",i));
1180  h5[i]=(TH1F *)list->FindObject(Form("hTimeLGEnt%d",i));
1181  h6[i]=(TH1F *)list->FindObject(Form("hTimeLGSumSq%d",i));
1182  }
1183  }
1184  //AliWarning("Input histograms read.");
1185 
1186  for(Int_t i=0;i<4;i++){
1187  hAllTimeAvBC[i]=new TH1F(Form("hAllTimeAvBC%d",i),Form("hAllTimeAvBC%d",i),h1[i]->GetNbinsX(),h1[i]->GetXaxis()->GetXmin(),h1[i]->GetXaxis()->GetXmax());
1188  hAllTimeRMSBC[i]=new TH1F(Form("hAllTimeRMSBC%d",i),Form("hAllTimeRMSBC%d",i),h3[i]->GetNbinsX(),h3[i]->GetXaxis()->GetXmin(),h3[i]->GetXaxis()->GetXmax());
1189 
1190  hAllTimeAvLGBC[i]=new TH1F(Form("hAllTimeAvLGBC%d",i),Form("hAllTimeAvLGBC%d",i),h4[i]->GetNbinsX(),h4[i]->GetXaxis()->GetXmin(),h4[i]->GetXaxis()->GetXmax());
1191  hAllTimeRMSLGBC[i]=new TH1F(Form("hAllTimeRMSLGBC%d",i),Form("hAllTimeRMSLGBC%d",i),h6[i]->GetNbinsX(),h6[i]->GetXaxis()->GetXmin(),h6[i]->GetXaxis()->GetXmax());
1192  }
1193 
1194  //AliWarning("New histograms booked.");
1195 
1196  //important remark: we use 'underflow bin' for absid=0 in OADB . That's why there is j-1 below.
1197  for(Int_t i=0;i<4;i++){
1198  for(Int_t j=1;j<=h1[i]->GetNbinsX();j++){
1199  //high gain
1200  if(h2[i]->GetBinContent(j)!=0){
1201  hAllTimeAvBC[i]->SetBinContent(j-1,h1[i]->GetBinContent(j)/h2[i]->GetBinContent(j));
1202  hAllTimeRMSBC[i]->SetBinContent(j-1,TMath::Sqrt(h3[i]->GetBinContent(j)/h2[i]->GetBinContent(j)) );
1203  } else {
1204  hAllTimeAvBC[i]->SetBinContent(j-1,0.);
1205  hAllTimeRMSBC[i]->SetBinContent(j-1,0.);
1206  }
1207  //low gain
1208  if(h5[i]->GetBinContent(j)!=0){
1209  hAllTimeAvLGBC[i]->SetBinContent(j-1,h4[i]->GetBinContent(j)/h5[i]->GetBinContent(j));
1210  hAllTimeRMSLGBC[i]->SetBinContent(j-1,TMath::Sqrt(h6[i]->GetBinContent(j)/h5[i]->GetBinContent(j)) );
1211  } else {
1212  hAllTimeAvLGBC[i]->SetBinContent(j-1,0.);
1213  hAllTimeRMSLGBC[i]->SetBinContent(j-1,0.);
1214  }
1215 
1216  }
1217  }
1218 
1219  //AliWarning("Average and rms calculated.");
1220  TFile *fileNew=new TFile(outputFile.Data(),"recreate");
1221  for(Int_t i=0;i<4;i++){
1222  hAllTimeAvBC[i]->Write();
1223  hAllTimeRMSBC[i]->Write();
1224  hAllTimeAvLGBC[i]->Write();
1225  hAllTimeRMSLGBC[i]->Write();
1226  }
1227 
1228  //AliWarning(Form("Histograms saved in %s file.",outputFile.Data()));
1229 
1230  fileNew->Close();
1231  delete fileNew;
1232 
1233  for(Int_t i=0;i<4;i++){
1234  delete hAllTimeAvBC[i];
1235  delete hAllTimeRMSBC[i];
1236  delete hAllTimeAvLGBC[i];
1237  delete hAllTimeRMSLGBC[i];
1238 
1239  delete h1[i];
1240  delete h2[i];
1241  delete h3[i];
1242  delete h4[i];
1243  delete h5[i];
1244  delete h6[i];
1245  }
1246 
1247  file->Close();
1248  delete file;
1249 
1250  //AliWarning("Pointers deleted. Memory cleaned.");
1251 }
1252 
1253 //________________________________________________________________________
1257 void AliAnalysisTaskEMCALTimeCalib::ProduceOffsetForSMsV2(Int_t runNumber,TString inputFile,TString outputFile, Bool_t offset100, Bool_t justL1phase){
1258 
1259 const Double_t lowerLimit[]={
1260  0,
1261  1152,
1262  2304,
1263  3456,
1264  4608,
1265  5760,
1266  6912,
1267  8064,
1268  9216,
1269  10368,
1270  11520,
1271  11904,
1272  12288,
1273  13056,
1274  13824,
1275  14592,
1276  15360,
1277  16128,
1278  16896,
1279  17280};
1280 
1281 const Double_t upperLimit[]={
1282  1151 ,
1283  2303 ,
1284  3455 ,
1285  4607 ,
1286  5759 ,
1287  6911 ,
1288  8063 ,
1289  9215 ,
1290  10367,
1291  11519,
1292  11903,
1293  12287,
1294  13055,
1295  13823,
1296  14591,
1297  15359,
1298  16127,
1299  16895,
1300  17279,
1301  17663};
1302 
1303  TFile *file =new TFile(inputFile.Data());
1304  if(file==0x0) return;
1305 
1306  TH1F *ccBC[4];
1307  for(Int_t i = 0; i < kNBCmask; i++){
1308  ccBC[i]=(TH1F*) file->Get(Form("hAllTimeAvBC%d",i));
1309  }
1310 
1311  TH1C *hRun=new TH1C(Form("h%d",runNumber),Form("h%d",runNumber),19,0,19);
1312  Int_t fitResult=0;
1313  Double_t minimumValue=10000.;
1314  Int_t minimumIndex=-1;
1315  Double_t meanBC[4];
1316 
1317  Double_t fitParameter=0;
1318  TF1 *f1=new TF1("f1","pol0",0,17664);
1319  Bool_t orderTest=kTRUE;
1320  Int_t iorder=0;//order index
1321  Int_t j=0;//BC index
1322  Int_t L1shift=0;
1323  Int_t totalValue=0;
1324 
1325  for(Int_t i=0;i<20;i++){
1326  minimumValue=10000;
1327  for(j=0;j<kNBCmask;j++){
1328  fitResult=ccBC[j]->Fit("f1","CQ","",lowerLimit[i],upperLimit[i]);
1329  if(fitResult<0){
1330  //hRun->SetBinContent(i,0);//correct it please
1331  meanBC[j]=-1;
1332  printf("Fit failed for SM %d BC%d, integral %f\n",i,j,ccBC[j]->Integral(lowerLimit[i],upperLimit[i]));
1333  continue;
1334  } else {
1335  fitParameter = f1->GetParameter(0);
1336  }
1337  if(offset100 && (j==0 || j==1)) {
1338  //the 100 ns offset was removed in LHC15n muon_calo_pass1 and further reconstructions
1339  fitParameter+=100;
1340  }
1341  meanBC[j]=fitParameter;
1342 
1343  if(fitParameter>0 && fitParameter<minimumValue){
1344  minimumValue = fitParameter;
1345  minimumIndex = j;
1346  }
1347  }
1348 
1349  if( minimumValue/25-(Int_t)(minimumValue/25)>0.5 ) {
1350  L1shift=(Int_t)(minimumValue/25.)+1;
1351  } else {
1352  L1shift=(Int_t)(minimumValue/25.);
1353  }
1354 
1355  if(TMath::Abs(minimumValue/25-(Int_t)(minimumValue/25)-0.5)<0.05)
1356  printf("Run %d, SM %d, min %f, next_min %f, next+1_min %f, next+2_min %f, min/25 %f, min%%25 %d, next_min/25 %f, next+1_min/25 %f, next+2_min/25 %f, SMmin %d\n",runNumber,i,minimumValue,meanBC[(minimumIndex+1)%4],meanBC[(minimumIndex+2)%4],meanBC[(minimumIndex+3)%4],minimumValue/25., (Int_t)((Int_t)minimumValue%25), meanBC[(minimumIndex+1)%4]/25., meanBC[(minimumIndex+2)%4]/25., meanBC[(minimumIndex+3)%4]/25., L1shift*25);
1357 
1358  if(justL1phase) totalValue = minimumIndex;
1359  else totalValue = L1shift<<2 | minimumIndex ;
1360  //printf("L1 phase %d, L1 shift %d *25ns= %d, L1p+L1s %d, total %d, L1pback %d, L1sback %d\n",minimumIndex,L1shift,L1shift*25,minimumIndex+L1shift,totalValue,totalValue&3,totalValue>>2);
1361 
1362  hRun->SetBinContent(i,totalValue);
1363  orderTest=kTRUE;
1364  for(iorder=minimumIndex;iorder<minimumIndex+4-1;iorder++){
1365  if( meanBC[(iorder+1)%4] <= meanBC[iorder%4] ) orderTest=kFALSE;
1366  }
1367  if(!orderTest)
1368  printf("run %d, SM %d, min index %d meanBC %f %f %f %f, order ok? %d\n",runNumber,i,minimumIndex,meanBC[0],meanBC[1],meanBC[2],meanBC[3],orderTest);
1369  }
1370  delete f1;
1371  TFile *fileNew=new TFile(outputFile.Data(),"update");
1372  hRun->Write();
1373  fileNew->Close();
1374  delete fileNew;
1375 
1376  file->Close();
1377  delete file;
1378 }
1379 
1380 //____________________________________________________
1382 {
1383  if(fBadChannelMapSet) return;
1384  AliOADBContainer *contBC=new AliOADBContainer("");
1385  contBC->InitFromFile(Form("%s/EMCALBadChannels.root","$ALICE_PHYSICS/OADB/EMCAL"),"AliEMCALBadChannels");
1386  printf("contBC %p, ent %d\n",contBC,contBC->GetNumberOfEntries());
1387  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRunNumber);
1388  if(arrayBC) {
1389  AliInfo("Remove EMCAL bad cells");
1391  for (Int_t i=0; i<kNSM; ++i) {
1392  TH2I *hbm = (TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
1393  if (!hbm) {
1394  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
1395  continue;
1396  }
1397  hbm->SetDirectory(0);
1398  fBadChannelMapArray->AddAt(hbm,i);
1399 
1400  } // loop over SMs
1401  } else AliInfo("Do NOT remove EMCAL bad channels\n"); // run array
1402 
1403  delete contBC;
1404  fBadChannelMapSet=kTRUE;
1405 } // Bad channel map loaded
1406 
1407 //____________________________________________________
1409 {
1410  if(fBadChannelMapSet) return;
1411 
1412  TFile *referenceFile = TFile::Open(fBadChannelFileName.Data());
1413  if(referenceFile==0x0) {
1414  AliFatal("*** NO bad channel map FILE");
1415  }
1416 
1417  TH1F *hbm = (TH1F*)referenceFile->Get("h1");
1418  if (!hbm) {
1419  AliError("Can not get EMCALBadChannelMap");
1420  }
1421  fBadChannelMapArray = new TObjArray(1);
1422  fBadChannelMapArray->AddAt(hbm,0);
1423  fBadChannelMapSet=kTRUE;
1424 } // Bad channel map loaded
1425 
1426 
1427 //_____________________________________________________________________
1432 }
Bool_t SetEMCalGeometry()
Set the EMCal Geometry.
TH1F * fhRawTimeEntriesLGBC[kNBCmask]
! 4 BCmask LG
Double_t fPassTimeMax
upper range of histo with time in passX
TH2F * fhEneVsAbsIdHG
! energy of each cell for high gain cells with strange time
TH1F * fhRawTimeSumLGBC[kNBCmask]
! 4 BCmask LG
Bool_t AcceptCluster(AliVCluster *clus)
Selection criteria of good cluster are set here.
TH2F * fhRawTimeVsIdBC[kNBCmask]
! 4 BCmask HG
double Double_t
Definition: External.C:58
Int_t fFineNbins
number of bins of histo with T0 time
Definition: External.C:236
TH1F * fhRawTimeSumSqBC[kNBCmask]
! 4 BCmask HG
Double_t fMinLambda0
minimum cluster lambda0
Int_t fEnergyLGNbins
number of bins of histo with energy LG
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow) const
TH1F * fhcalcEvtTime
! spectrum calcolot0[0]
Bool_t fBadReco
flag to apply 100ns shift and L1 shift
Double_t fEnergyMin
lower range of histo with energy HG
TString fReferenceRunByRunFileName
name of reference file (run-by-run)
TObjArray * fL1PhaseList
array with phases for set of runs
Double_t fMinTime
minimum cluster time after correction
AliTOFT0maker * fTOFmaker
pointer to get T0 from TOF
TList * list
Double_t fFineTmax
upper range of histo with T0 time
void SetDefaultCuts()
Set default cuts for calibration.
Int_t fSetBadChannelMapSource
switch to load BC map 0-no BC,1-OADB,2-file
void LoadBadChannelMap()
Load Bad Channel Map from different source.
Double_t fRawTimeMin
lower range of histo with raw time
TH2F * fhTimeDsupLG[kNSM]
! 20 SM low gain
TH1F * fhRawTimeSumSqLGBC[kNBCmask]
! 4 BCmask LG
Double_t fMaxRtrack
maximum cluster track distance
Double_t fMaxLambda0
maximum cluster lambda0
AliEMCALGeometry * fgeom
pointer to EMCal geometry
Double_t fRawTimeMax
upper range of histo with raw time
TH2F * fhTcellvsTOFT0HD
! time of cell vs TOF T0 time for higher energy threshold
TH2F * fhRawCorrTimeVsIdBC[kNBCmask]
! 4 BCmask HG
TH2F * fhRawCorrTimeVsIdLGBC[kNBCmask]
! 4 BCmask LG
Bool_t IsLowGainCellInCluster(AliVCluster *clus)
Check if low gain cell is in a cluster.
Double_t fMaxTime
maximum cluster time after correction
int Int_t
Definition: External.C:63
Double_t fEnergyMax
upper range of histo with energy HG
Double_t fMinClusterEnergy
minimum cluster energy
float Float_t
Definition: External.C:68
static void ProduceOffsetForSMsV2(Int_t runNumber, TString inputFile="Reference.root", TString outputFile="ReferenceSM.root", Bool_t offset100=kTRUE, Bool_t justL1phase=kTRUE)
TObjArray * fBadChannelMapArray
bad channel map array
Double_t fEnergyLGMax
upper range of histo with energy LG
Double_t fFineTmin
lower range of histo with T0 time
TString fReferenceFileName
! name of reference file (for one period)
Int_t fMaxNcells
maximum number of cells in cluster
virtual void UserExec(Option_t *option)
Main loop executed for each event.
virtual void PrepareTOFT0maker()
Get T0 time from TOF.
Bool_t fPileupFromSPD
flag to set PileupFromSPD
TH1F * fhRawTimeSumBC[kNBCmask]
! 4 BCmask HG
Double_t fMinCellEnergy
minimum cell energy
Int_t fMinNcells
minimum number of cells in cluster
Double_t fEnergyLGMin
lower range of histo with energy LG
Task to work on Time Calibration for EMCal/DCal.
Double_t fMinLambda0LG
minimum cluster lambda0 Low Gain
Double_t fMaxClusterEnergy
maximum cluster energy
Bool_t fBadChannelMapSet
flag whether bad channel map is set
TH1F * fhEvtTimeHeader
! spectrum time from header
TH1F * fhEvtTimeDiff
! spectrum time difference
TH2F * fhTcellvsTOFT0
! time of cell vs TOF T0 time
Double_t fPassTimeMin
lower range of histo with time in passX
Double_t fMaxLambda0LG
maximum cluster lambda0 Low Gain
static void ProduceCalibConsts(TString inputFile="time186319testWOL0.root", TString outputFile="Reference.root", Bool_t isFinal=kFALSE)
TH2F * fhEneVsAbsIdLG
! energy of each cell for low gain cells with strange time
Int_t fPassTimeNbins
number of bins of histo with time in passX
TList * fOutputList
pointer to output list
TH1F * fhRawTimeEntriesBC[kNBCmask]
! 4 BCmask HG
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TH2F * fhTimeVsIdBC[kNBCmask]
! 4 BCmask HG
TFile * file
TH2F * fhTimeVsIdLGBC[kNBCmask]
! 4 BCmask LG
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
Bool_t CheckCellRCU(Int_t nSupMod, Int_t icol, Int_t irow)
Check RCU for cell given by Super Module, column index, row index.
TH2F * fhRawTimeVsIdLGBC[kNBCmask]
! 4 BCmask LG
void LoadReferenceHistos()
Load reference Histograms (for one period) from file.
bool Bool_t
Definition: External.C:53
TString fBadChannelFileName
name of file with bad channels
TH1F * fhAllAverageLGBC[kNBCmask]
4 BCmask High gain
Bool_t fFillHeavyHisto
flag to fill heavy histograms
Int_t fEnergyNbins
number of bins of histo with energy HG
Int_t fRawTimeNbins
number of bins of histo with raw time
TH2F * fhTimeDsup[kNSM]
! 20 SM high gain
TH2F * fhTimeDsupBC[kNSM][kNBCmask]
! 20 x 4 high gain
TH2F * fhTimeDsupLGBC[kNSM][kNBCmask]
! 20 x 4 low gain