AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 <TList.h>
24 #include <TCanvas.h>
25 #include <TGeoManager.h>
26 #include <TRefArray.h>
27 
28 #include "AliLog.h"
29 #include "AliAnalysisTask.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDEvent.h"
32 #include "AliAODEvent.h"
33 #include "AliVEvent.h"
34 #include "AliESDInputHandler.h"
35 #include "AliAODInputHandler.h"
36 #include "AliESDpid.h"
37 #include "AliTOFcalib.h"
38 #include "AliCDBManager.h"
39 #include "AliRunTag.h"
40 
41 #include "AliTOFT0maker.h"
42 #include "AliVCluster.h"
43 #include "AliESDCaloCluster.h"
44 #include "AliVCaloCells.h"
45 #include "AliESDCaloCells.h"
46 #include "AliAODCaloCluster.h"
47 #include "AliAODCaloCells.h"
48 #include "AliEMCALGeometry.h"
49 
51 
55 
56 using std::cout;
57 using std::endl;
58 
59 //________________________________________________________________________
62 : AliAnalysisTaskSE(name),
63  fRunNumber(-1),
64  fTOFmaker(0),
65  fOutputList(0x0),
66  fgeom(0),
67  fGeometryName(),
68  fMinClusterEnergy(0),
69  fMaxClusterEnergy(0),
70  fMinNcells(0),
71  fMaxNcells(0),
72  fMinLambda0(0),
73  fMaxLambda0(0),
74  fMinLambda0LG(0),
75  fMaxLambda0LG(0),
76  fMaxRtrack(0),
77  fMinCellEnergy(0),
78  fReferenceFileName(),
79  fPileupFromSPD(kFALSE),
80  fMinTime(0),
81  fMaxTime(0),
82  fhcalcEvtTime(0),
83  fhEvtTimeHeader(0),
84  fhEvtTimeDiff(0),
85  fhEventType(0),
86  fhTOFT0vsEventNumber(0),
87  fhTcellvsTOFT0(0),
88  fhTcellvsTOFT0HD(0),
89  fhTcellvsSM(0),
90  fhEneVsAbsIdHG(0),
91  fhEneVsAbsIdLG(0),
92  fhTimeVsBC(0),
93  fhTimeSumSq(),
94  fhTimeEnt(),
95  fhTimeSum(),
96  fhTimeLGSumSq(),
97  fhTimeLGEnt(),
98  fhTimeLGSum(),
99  fhAllAverageBC(),
100  fhAllAverageLGBC(),
101  fhTimeDsup(),
102  fhTimeDsupBC(),
103  fhRawTimeVsIdBC(),
104  fhRawTimeSumBC(),
105  fhRawTimeEntriesBC(),
106  fhRawTimeSumSqBC(),
107  fhRawTimeVsIdLGBC(),
108  fhRawTimeSumLGBC(),
109  fhRawTimeEntriesLGBC(),
110  fhRawTimeSumSqLGBC()
111 
112 {
113  for(Int_t i = 0; i < kNBCmask; i++)
114  {
115  fhAllAverageBC[i]=0;
116  fhAllAverageLGBC[i]=0;
117 
118  fhTimeSumSq[i]=0;
119  fhTimeEnt[i]=0;
120  fhTimeSum[i]=0;
121 
122  fhTimeLGSumSq[i]=0;
123  fhTimeLGEnt[i]=0;
124  fhTimeLGSum[i]=0;
125 
126  fhRawTimeVsIdBC[i]=0;
127  fhRawTimeSumBC[i]=0;
128  fhRawTimeEntriesBC[i]=0;
129  fhRawTimeSumSqBC[i]=0;
130 
131  fhRawTimeVsIdLGBC[i]=0;
132  fhRawTimeSumLGBC[i]=0;
134  fhRawTimeSumSqLGBC[i]=0;
135  }
136 
137  //set default cuts for calibration and geometry name
138  SetDefaultCuts();
139 
140  //T0 TOF time
142 
143  // Define input and output slots here
144  // Input slot #0 works with a TChain
145  DefineInput(0, TChain::Class());
146 
147  // Output slot #0 id reserved by the base class for AOD
148  // Output slot #1 writes into a TH1 container
149  DefineOutput(1, TList::Class());
150 
151 } // End ctor
152 
153 //_____________________________________________________________________
159 //void AliAnalysisTaskEMCALTimeCalib::LocalInit()
160 //{
161 // AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::LocalInit()");
162 //}
163 
165 //_____________________________________________________________________
167 {
168  if(fReferenceFileName.Length()!=0){
169  TFile *myFile = TFile::Open(fReferenceFileName.Data());
170  AliInfo(Form("Reference file: %s, pointer %p",fReferenceFileName.Data(),myFile));
171  if(myFile==0x0) {
172  AliFatal("*** NO REFERENCE FILE");
173  } else {
174  AliDebug(1,"*** OK TFILE");
175  // connect ref run here
176  for(Int_t i = 0; i < kNBCmask; i++)
177  {
178  fhAllAverageBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvBC%d",i));
179  if(fhAllAverageBC[i]==0x0) AliFatal(Form("Reference histogram for BC%d does not exist",i));
180  if(fhAllAverageBC[i]->GetEntries()==0)AliWarning(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
181  fhAllAverageLGBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvLGBC%d",i));
182  if(fhAllAverageLGBC[i]==0x0) AliFatal(Form("Reference LG histogram for BC%d does not exist",i));
183  if(fhAllAverageLGBC[i]->GetEntries()==0)AliFatal(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
184  }
185 
186  AliDebug(1,Form("hAllAverage entries BC0 %d", (Int_t)fhAllAverageBC[0]->GetEntries() ));
187  AliDebug(1,Form("hAllAverage entries BC2 %d",(Int_t)fhAllAverageBC[2]->GetEntries() ));
188  AliDebug(1,Form("hAllAverageLG entries BC0 %d", (Int_t)fhAllAverageLGBC[0]->GetEntries() ));
189  AliDebug(1,Form("hAllAverageLG entries BC2 %d",(Int_t)fhAllAverageLGBC[2]->GetEntries() ));
190 
191  }
192  } else {
193 //end of reference file is provided
194  AliFatal("You require to load reference histos from file but FILENAME is not provided");
195  }
196 } // End of AliAnalysisTaskEMCALTimeCalib::LoadReferenceHistos()
197 
198 //_____________________________________________________________________
202 {
203  AliDebug(1,"AnalysisTaskEMCalTimeCalib::Notify()");
204  AliDebug(2,Form("Notify(): EMCal geometry: fgeom = %p, fGeometryName=%s\n ",fgeom,fGeometryName.Data()));
205 
206  if (!InputEvent())
207  {
208  AliFatal("ERROR: InputEvent not set");
209  return kFALSE;
210  }
211  else AliDebug(1,"Good, InputEvent set");
212 
213  fRunNumber = InputEvent()->GetRunNumber();
214  AliDebug(1,Form("RunNumber %d", fRunNumber));
215 
216  // Init EMCAL geometry
217  if (!fgeom) SetEMCalGeometry();
218  //Init EMCAL geometry done
219 
220  return kTRUE;
221 }
222 
223 //_____________________________________________________________________
226 {
227  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::SetEMCalGeometry()");
228  if(fGeometryName.Length()==0){
229  fgeom=AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
230  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fgeom->GetName(),fRunNumber));
231  } else {
232  fgeom = AliEMCALGeometry::GetInstance(fGeometryName.Data());
233  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeometryName.Data()));
234  }
235 
236  if (!fgeom){
237  AliWarning("Make sure the EMCal geometry is set properly !");
238  } else {
239  AliDebug(1,Form("EMCal geometry properly set: fGeom = %p, fGeometryName=%s",fgeom,fGeometryName.Data()));
240  }
241 
242  return kTRUE;
243 }
244 
245 //_____________________________________________________________________
248 {
249  //method under development
250  AliInfo(Form("<D> -- Run # = %d", fRunNumber));
251  AliInfo("prepare TOFT0maker!!");
252  //cout<<"Run "<< fRunNumber<<" in TOFT0maker"<<endl;
253 
254 
255  AliCDBManager * cdb = AliCDBManager::Instance();
256  cdb->SetDefaultStorage("raw://");
257  cdb->SetRun(fRunNumber);
258 
259  AliESDpid *extPID=new AliESDpid();
260 
261  // Wonder if some have to be declared as private variables??
262  // AliESDpid *extPID = new AliESDpid();
263  // AliTOFcalib * tofCalib = new AliTOFcalib();
264  // tofCalib->SetCalibrateTOFsignal(kTRUE);
265  // tofCalib->Init();
266 
267  fTOFmaker = new AliTOFT0maker(extPID);
268  fTOFmaker->SetTimeResolution(115.0); // if you want set the TOF res
269  // fTOFmaker = new AliTOFT0maker(extPID,tofCalib);
270  // fTOFmaker->SetTimeResolution(130.0);
271 
272  //cout<<"extPID "<<extPID<<" fTOFmaker "<<fTOFmaker<<endl;
273 
274 }// End PrepareTOFT0maker
275 
276 //________________________________________________________________________
280 {
281  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::UserCreateOutputObjects()");
282 
283  Double_t fineTmin = -500;
284  Double_t fineTmax = 400;
285  Double_t fineInterval = 0.20;
286  Int_t nfinebin = (fineTmax-fineTmin)/fineInterval;
287  //cout << "<D> nfinebin = " << nfinebin << endl;
288 
289  Int_t nChannels = 17664;
290  //book histograms
291  fhcalcEvtTime = new TH1F("fhcalcEvtTime","calculated event time from T0",nfinebin, fineTmin,fineTmax);
292  fhcalcEvtTime->GetXaxis()->SetTitle("T ");
293  fhcalcEvtTime->GetYaxis()->SetTitle("Counts (a.u.)");
294 
295  fhEvtTimeHeader = new TH1F("fhEvtTimeHeader","event time from header",nfinebin, fineTmin,fineTmax);
296  fhEvtTimeHeader->GetXaxis()->SetTitle("T ");
297  fhEvtTimeHeader->GetYaxis()->SetTitle("Counts (a.u.)");
298 
299  fhEvtTimeDiff = new TH1F("fhEvtTimeDiff","event time difference",nfinebin, fineTmin,fineTmax);
300  fhEvtTimeDiff->GetXaxis()->SetTitle("#Delta T ");
301  fhEvtTimeDiff->GetYaxis()->SetTitle("Counts (a.u.)");
302 
303  fhEventType = new TH1F("fhEventType","event type",10, 0.,10.);
304  fhEventType ->GetXaxis()->SetTitle("Type ");
305  fhEventType ->GetYaxis()->SetTitle("Counts (a.u.)");
306  fhTcellvsTOFT0 = new TH2F("hTcellvsTOFT0", " T_cell vs TOFT0", 500,-600.0,+400.0,1200,300.0,900.0);
307  fhTcellvsTOFT0HD = new TH2F("hTcellvsTOFT0HD", " T_cell vs TOFT0,HighEnergy", 500,-600.0,+400.0,4000,500.0,700.0);
308  fhTcellvsSM = new TH2F("hTcellvsSM", " T_cell vs SM", 20,0,20,300,300,900);
309  fhEneVsAbsIdHG = new TH2F("fhEneVsAbsIdHG", "energy vs ID for HG",1000,0,18000,200,0,10);
310  fhEneVsAbsIdLG = new TH2F("fhEneVsAbsIdLG", "energy vs ID for LG",1000,0,18000,100,0,20);
311 
312  for (Int_t i = 0; i < kNBCmask ; i++)
313  {
314  //already after correction
315  //high gain
316  fhTimeSumSq[i] = new TH1F(Form("hTimeSumSq%d", i),
317  Form("cell Sum Square time HG, BC %d ", i),
318  nChannels,0.,(Double_t)nChannels);
319  fhTimeSumSq[i]->SetYTitle("Sum Sq Time ");
320  fhTimeSumSq[i]->SetXTitle("AbsId");
321 
322  fhTimeSum[i] = new TH1F(Form("hTimeSum%d", i),
323  Form("cell Sum time HG, BC %d ", i),
324  nChannels,0.,(Double_t)nChannels);
325  fhTimeSum[i]->SetYTitle("Sum Time ");
326  fhTimeSum[i]->SetXTitle("AbsId");
327 
328  fhTimeEnt[i] = new TH1F(Form("hTimeEnt%d", i),
329  Form("cell Entries HG, BC %d ", i),
330  nChannels,0.,(Double_t)nChannels);
331  fhTimeEnt[i]->SetYTitle("Entries for Time ");
332  fhTimeEnt[i]->SetXTitle("AbsId");
333 
334  //low gain
335  fhTimeLGSumSq[i] = new TH1F(Form("hTimeLGSumSq%d", i),
336  Form("cell Sum Square time LG, BC %d ", i),
337  nChannels,0.,(Double_t)nChannels);
338  fhTimeLGSumSq[i]->SetYTitle("Sum Sq Time ");
339  fhTimeLGSumSq[i]->SetXTitle("AbsId");
340 
341  fhTimeLGSum[i] = new TH1F(Form("hTimeLGSum%d", i),
342  Form("cell Sum time LG, BC %d ", i),
343  nChannels,0.,(Double_t)nChannels);
344  fhTimeLGSum[i]->SetYTitle("Sum Time ");
345  fhTimeLGSum[i]->SetXTitle("AbsId");
346 
347  fhTimeLGEnt[i] = new TH1F(Form("hTimeLGEnt%d", i),
348  Form("cell Entries LG, BC %d ", i),
349  nChannels,0.,(Double_t)nChannels);
350  fhTimeLGEnt[i]->SetYTitle("Entries for Time ");
351  fhTimeLGEnt[i]->SetXTitle("AbsId");
352 
353  //raw time histograms
354  //high gain
355  fhRawTimeVsIdBC[i] = new TH2F(Form("RawTimeVsIdBC%d", i),
356  Form("cell raw time vs ID for high gain BC %d ", i),
357  nChannels,0.,(Double_t)nChannels,600,300,900);
358  fhRawTimeVsIdBC[i]->SetXTitle("AbsId");
359  fhRawTimeVsIdBC[i]->SetYTitle("Time");
360 
361  fhRawTimeSumBC[i] = new TH1F(Form("RawTimeSumBC%d", i),
362  Form("sum of cell raw time for high gain BC %d ", i),
363  nChannels,0.,(Double_t)nChannels);
364  fhRawTimeSumBC[i]->SetXTitle("AbsId");
365  fhRawTimeSumBC[i]->SetYTitle("Sum Time");
366 
367  fhRawTimeEntriesBC[i] = new TH1F(Form("RawTimeEntriesBC%d", i),
368  Form("No. entries of cells raw time for high gain BC %d ", i),
369  nChannels,0.,(Double_t)nChannels);
370  fhRawTimeEntriesBC[i]->SetXTitle("AbsId");
371  fhRawTimeEntriesBC[i]->SetYTitle("Entries for Time ");
372 
373  fhRawTimeSumSqBC[i] = new TH1F(Form("RawTimeSumSqBC%d", i),
374  Form("sum of (cell raw time)^2 for high gain BC %d ", i),
375  nChannels,0.,(Double_t)nChannels);
376  fhRawTimeSumSqBC[i]->SetXTitle("AbsId");
377  fhRawTimeSumSqBC[i]->SetYTitle("Sum Sq Time");
378 
379  //low gain
380  fhRawTimeVsIdLGBC[i] = new TH2F(Form("RawTimeVsIdLGBC%d", i),
381  Form("cell raw time vs ID for low gain BC %d ", i),
382  nChannels,0.,(Double_t)nChannels,600,300,900);
383  fhRawTimeVsIdLGBC[i]->SetXTitle("AbsId");
384  fhRawTimeVsIdLGBC[i]->SetYTitle("Time");
385 
386  fhRawTimeSumLGBC[i] = new TH1F(Form("RawTimeSumLGBC%d", i),
387  Form("sum of cell raw time for low gain BC %d ", i),
388  nChannels,0.,(Double_t)nChannels);
389  fhRawTimeSumLGBC[i]->SetXTitle("AbsId");
390  fhRawTimeSumLGBC[i]->SetYTitle("Sum Time");
391 
392  fhRawTimeEntriesLGBC[i] = new TH1F(Form("RawTimeEntriesLGBC%d", i),
393  Form("No. entries of cells raw time for low gain BC %d ", i),
394  nChannels,0.,(Double_t)nChannels);
395  fhRawTimeEntriesLGBC[i]->SetXTitle("AbsId");
396  fhRawTimeEntriesLGBC[i]->SetYTitle("Entries for Time ");
397 
398  fhRawTimeSumSqLGBC[i] = new TH1F(Form("RawTimeSumSqLGBC%d", i),
399  Form("sum of (cell raw time)^2 for low gain BC %d ", i),
400  nChannels,0.,(Double_t)nChannels);
401  fhRawTimeSumSqLGBC[i]->SetXTitle("AbsId");
402  fhRawTimeSumSqLGBC[i]->SetYTitle("Sum Sq Time");
403 
404 
405  for (Int_t j = 0; j < kNSM ; j++)
406  {
407  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,1400,-350.0,350.0);
408  fhTimeDsupBC[j][i]->SetYTitle(" Time (ns) ");
409  fhTimeDsupBC[j][i]->SetXTitle(" E (GeV) ");
410  }
411  }
412 
413  for (Int_t jj = 0; jj < kNSM ; jj++)
414  {
415  fhTimeDsup[jj] = new TH2F(Form("SupMod%d",jj), Form("SupMod %d time_vs_E ",jj),500,0.0,20.0,1400,-350.0,350.0);
416  fhTimeDsup[jj]->SetYTitle(" Time (ns) ");
417  fhTimeDsup[jj]->SetXTitle(" E (GeV) ");
418  }
419 
420  fhTimeVsBC = new TH2F("TimeVsBC"," SupMod time_vs_BC ", 4001,-0.5,4000.5,400,200.0,1000.0);
421 
422 
423  //add histos to list
424  fOutputList = new TList();
425 
429  fOutputList->Add(fhEventType);
432  fOutputList->Add(fhTcellvsSM);
435 
436  for (Int_t i = 0; i < kNBCmask ; i++)
437  {
438  fOutputList->Add(fhTimeSumSq[i]);
439  fOutputList->Add(fhTimeEnt[i]);
440  fOutputList->Add(fhTimeSum[i]);
441 
442  fOutputList->Add(fhTimeLGSumSq[i]);
443  fOutputList->Add(fhTimeLGEnt[i]);
444  fOutputList->Add(fhTimeLGSum[i]);
445 
446  fOutputList->Add(fhRawTimeVsIdBC[i]);
447  fOutputList->Add(fhRawTimeSumBC[i]);
449  fOutputList->Add(fhRawTimeSumSqBC[i]);
450 
452  fOutputList->Add(fhRawTimeSumLGBC[i]);
455 
456  for (Int_t j = 0; j < kNSM ; j++){
457  fOutputList->Add(fhTimeDsupBC[j][i]);
458  }
459  }
460 
461  for (Int_t j = 0; j < kNSM ; j++)
462  {
463  fOutputList->Add(fhTimeDsup[j]);
464  }
465 
466  fOutputList->Add(fhTimeVsBC);
467 
468  fOutputList->SetOwner(kTRUE);
469  PostData(1,fOutputList);
470 
471 
472 } // End of AliAnalysisTaskEMCALTimeCalib::UserCreateOuputObjects()
473 
474 //________________________________________________________________________
477 {
478  // Called for each event
479  AliDebug(2,Form("UserExec: EMCal geometry: fgeom = %p fGeometryName %s",fgeom,fGeometryName.Data()));
480  AliVEvent *event = InputEvent();
481  //cout<<"T0TOF "<<event->GetT0TOF()<<endl;//bad idea
482  //cout<< fEvent->GetTOFHeader()->GetDefaultEventTimeVal()<<endl;
483  AliDebug(2,Form("TOF time from header %f ps",event->GetTOFHeader()->GetDefaultEventTimeVal()));
484  fhEvtTimeHeader->Fill(event->GetTOFHeader()->GetDefaultEventTimeVal());
485 
486  //fEvent = dynamic_cast<AliESDEvent*>(event);
487  if (!event) {
488  AliError("ESD not available, exit");
489  fhEventType->Fill(0.5);
490  return;
491  }
492 
493  if(fPileupFromSPD==kTRUE){
494  if(event->IsPileupFromSPD(3,0.8,3.,2.,5.)){
495  AliDebug(1,"Event: PileUp skip.");
496  fhEventType->Fill(1.5);
497  return;
498  }
499  }
500 
501  TString triggerclasses = event->GetFiredTriggerClasses();
502  if(triggerclasses=="") {
503  fhEventType->Fill(2.5);
504  return;
505  }
506 
507  Int_t eventType = ((AliVHeader*)event->GetHeader())->GetEventType();
508  // physics events eventType=7, select only those
509  AliDebug(1,Form("Triggerclasses %s, eventType %d",triggerclasses.Data(),eventType));
510  if(eventType != 7) {
511  fhEventType->Fill(3.5);
512  return;
513  }
514 
515  // Check trigger
516  Bool_t bMB = kFALSE;
517  Bool_t bL0 = kFALSE;
518  Bool_t bL1G = kFALSE;
519  Bool_t bL1J = kFALSE;
520 
521  if(triggerclasses.Contains("CINT7-B-NOPF-ALLNOTRD") ||
522  triggerclasses.Contains("CINT7-I-NOPF-ALLNOTRD") ||
523  triggerclasses.Contains("CINT1-I-NOPF-ALLNOTRD") ||
524  triggerclasses.Contains("CINT1-B-NOPF-ALLNOTRD") ||
525  triggerclasses.Contains("CINT8") ||
526  triggerclasses.Contains("CINT7") ||
527  triggerclasses.Contains("CPBI2_B1-B-NOPF-ALLNOTRD") ) bMB = kTRUE;
528 
529  if(triggerclasses.Contains("CEMC7-B-NOPF-CENTNOTRD") ||
530  triggerclasses.Contains("CEMC1-B-NOPF-CENTNOTRD") ||
531  triggerclasses.Contains("CEMC7") ||
532  triggerclasses.Contains("CEMC8") ||
533  triggerclasses.Contains("CEMC8-B-NOPF-CENTNOTRD") ) bL0 = kTRUE;
534 
535  if(triggerclasses.Contains("CEMC7EG1-B-NOPF-CENTNOTRD") ||
536  triggerclasses.Contains("CEMC7EG2-B-NOPF-CENTNOTRD") ||
537  triggerclasses.Contains("CEMC8EG1-B-NOPF-CENTNOTRD") ||
538  triggerclasses.Contains("CEMC8EGA") ||
539  triggerclasses.Contains("CEMC7EGA") ||
540  triggerclasses.Contains("CPBI2EGA") ) bL1G = kTRUE;
541 
542 
543  if(triggerclasses.Contains("CEMC7EJ1-B-NOPF-CENTNOTRD") ||
544  triggerclasses.Contains("CEMC7EJ2-B-NOPF-CENTNOTRD") ||
545  triggerclasses.Contains("CEMC8EJ1-B-NOPF-CENTNOTRD") ||
546  triggerclasses.Contains("CEMC7EJE") ||
547  triggerclasses.Contains("CEMC8EJE") ||
548  triggerclasses.Contains("CPBI2EJE") ) bL1J = kTRUE;
549 
550  if( bL1G || bL1J || bL0 ){ fhEventType->Fill(4.5);}
551  if( bMB ){ fhEventType->Fill(5.5);}
552 
553 
554  // if(bL1G || bL1J || bL0){
555 
556  // Prepare TOFT0 maker at the beginning of a run
557 // if (event->GetRunNumber() != fRunNumber)
558 // {
559 // fRunNumber = event->GetRunNumber();
560 // // PrepareTOFT0maker();
561 // // cout<<"tofT0maker per run"<<fRunNumber<<endl;
562 // }// fi Check if run number has changed
563 
564  // --- Use of AliTOFT0maker
565  Double_t calcolot0=0.0;
566  if(!AODEvent()){
567  Double_t* timeTOFtable;
568  timeTOFtable=fTOFmaker->ComputeT0TOF(dynamic_cast<AliESDEvent*>(event));
569  AliDebug(2,Form("TOF time %f ps, resolution %f ps, tracks at TOF %f/used %f",timeTOFtable[0],timeTOFtable[1],timeTOFtable[3],timeTOFtable[7]));
570  //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;
571  calcolot0=timeTOFtable[0];
572  }
573 
574  if (!fhcalcEvtTime) {
575  AliWarning("<E> fhcalcEvtTime not available");
576  return;
577  }// fi no simple histo present
578 
579  fhcalcEvtTime->Fill(calcolot0);
580  if(calcolot0 != 0 && event->GetTOFHeader()->GetDefaultEventTimeVal() != 0 )
581  fhEvtTimeDiff->Fill(calcolot0-event->GetTOFHeader()->GetDefaultEventTimeVal());
582 
583  TRefArray* caloClusters = new TRefArray();
584  event->GetEMCALClusters(caloClusters);
585  // cout << " ###########Bunch Cross nb = " << event->GetBunchCrossNumber() << endl;
586 
587  Int_t BunchCrossNumber =event->GetBunchCrossNumber();
588 
589  Float_t offset=0.;
590 
591  Int_t nBC = 0;
592  nBC = BunchCrossNumber%4;
593  //Int_t nTriggerMask =event->GetTriggerMask();
594  // cout << " nBC " << nBC << " nTriggerMask " << nTriggerMask<< endl;
595  Float_t timeBCoffset = 0.; //manual offest
596  // if( nBC%4 ==0 || nBC%4==1) timeBCoffset = 100.; // correction was for LHC11 when BC was not corrected
597 
598  Int_t nclus = caloClusters->GetEntries();
599  AliDebug(1,Form("###########Bunch Cross nb = %d nclus = %d",nBC,nclus ));
600  //cout << " ###########Bunch Cross nb = " << nBC <<" nclus= "<< nclus<< endl;
601  //Int_t ntracks = event-> GetNumberOfTracks() ;
602 
603  AliVCaloCells &cells= *(event->GetEMCALCells());//it is cluster independent
604  //Variables used plenty in loops
605  Int_t nSupMod=-1, nModule=-1;
606  Int_t iphi=-1, ieta=-1, nIphi=-1, nIeta=-1;
607  Int_t absId=-1;
608  Float_t hkdtime=0.0;
609  Float_t amp=0.0;
610  Bool_t isHighGain=kTRUE;
611 
612  for (Int_t icl = 0; icl < nclus; icl++) {
613  //ESD and AOD CaloCells carries the same information
614  AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
615  if(!AcceptCluster(clus)) continue;
616 
617  //cout<<"nCells="<< clus->GetNCells();<<endl;
618 
619  UShort_t * index = clus->GetCellsAbsId() ;
620 
621  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
622  absId = index[i]; // or clus->GetCellNumber(i) ;
623  hkdtime = cells.GetCellTime(absId) * 1.0e09; // to get ns
624  amp = cells.GetCellAmplitude(absId) ;
625  isHighGain = cells.GetCellHighGain(absId);
626  //cout<<"cell absID: "<<absId<<" cellTime: "<<hkdtime<<" cellaplit: "<< amp<<endl;
627 
628  //main histograms with raw time information
629  if(amp>fMinCellEnergy){
630  if(isHighGain){
631  fhRawTimeVsIdBC[nBC]->Fill(absId,hkdtime);
632  fhRawTimeSumBC[nBC]->Fill(absId,hkdtime);
633  fhRawTimeEntriesBC[nBC]->Fill(absId,1.);
634  fhRawTimeSumSqBC[nBC]->Fill(absId,hkdtime*hkdtime);
635  }else{
636  fhRawTimeVsIdLGBC[nBC]->Fill(absId,hkdtime);
637  fhRawTimeSumLGBC[nBC]->Fill(absId,hkdtime);
638  fhRawTimeEntriesLGBC[nBC]->Fill(absId,1.);
639  fhRawTimeSumSqLGBC[nBC]->Fill(absId,hkdtime*hkdtime);
640  }
641  }
642  //fgeom->PrintCellIndexes(absId);
643  //fgeom->PrintCellIndexes(absId,1);
644 
645  // GEOMETRY tranformations
646  fgeom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
647  fgeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
648 
649  // other histograms for cross-check
650  CheckCellRCU(nSupMod,ieta,iphi);//SM, column, row
651 
652  fhTcellvsSM->Fill(nSupMod,hkdtime);
653  if(isHighGain==kTRUE) {fhEneVsAbsIdHG->Fill(absId,amp);}
654  else {fhEneVsAbsIdLG->Fill(absId,amp);}
655 
656  fhTimeVsBC->Fill(1.*BunchCrossNumber,hkdtime-timeBCoffset);
657  //important remark: We use 'Underflow bin' for absid=0 in OADB for time calibration
658  if(isHighGain==kTRUE){
659  if(fhAllAverageBC[nBC]!=0) {//comming from file after the first iteration
660  offset = (Float_t)(fhAllAverageBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
661  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
662  AliFatal(Form("Reference histogram for BC%d not properly loaded",nBC));
663  }
664  } else {
665  if(fhAllAverageLGBC[nBC]!=0) {//comming from file after the first iteration
666  offset = (Float_t)(fhAllAverageLGBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
667  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
668  AliFatal(Form("Reference LG histogram for BC%d not properly loaded",nBC));
669  }
670  }
671  //if(offset==0)cout<<"offset 0 in SM "<<nSupMod<<endl;
672 
673 
674  if(amp>0.5) {
675  fhTimeDsup[nSupMod]->Fill(amp,hkdtime-offset);
676  fhTimeDsupBC[nSupMod][nBC]->Fill(amp,hkdtime-offset);
677  }
678 
679  if(amp>0.9) {
680  fhTcellvsTOFT0HD->Fill(calcolot0, hkdtime);
681  }
682 
683  fhTcellvsTOFT0->Fill(calcolot0, hkdtime-offset);
684 
685  hkdtime = hkdtime-timeBCoffset;//time corrected by manual offset (default=0)
686  Float_t hkdtimecorr;
687  hkdtimecorr= hkdtime-offset;//time after first iteration
688 
689  //main histograms after the first itereation for calibration constants
690  //if(hkdtimecorr>=-20. && hkdtimecorr<=20. && amp>0.9 ) {
691  if(hkdtimecorr>=fMinTime && hkdtimecorr<=fMaxTime && amp>fMinCellEnergy ) {
692  // per cell
693 // Float_t entriesTime=fhTimeEnt[nBC]->GetBinContent(absId)+1;
694 // Float_t sumTimeSq=(fhTimeSumSq[nBC]->GetBinContent(absId)+(hkdtime*hkdtime));
695 // Float_t sumTime=(fhTimeSum[nBC]->GetBinContent(absId)+hkdtime);
696 //
697 // fhTimeEnt[nBC]->SetBinContent(absId,entriesTime);
698 // fhTimeSumSq[nBC]->SetBinContent(absId,sumTimeSq);
699 // fhTimeSum[nBC]->SetBinContent(absId,sumTime);
700 
701  if(isHighGain){
702  fhTimeEnt[nBC]->Fill(absId,1.);
703  fhTimeSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
704  fhTimeSum[nBC]->Fill(absId,hkdtime);
705  }else{
706  fhTimeLGEnt[nBC]->Fill(absId,1.);
707  fhTimeLGSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
708  fhTimeLGSum[nBC]->Fill(absId,hkdtime);
709  }
710 
711 
712  } // hkdtime:[-20;20]
713  } // end icell
714  } //end cluster
715 
716 
717  // Post output data.
718  //cout<<"Post data and delete caloClusters"<<endl;
719  caloClusters->Delete();
720  delete caloClusters;
721 // } // end if trigger type
722 
723  PostData(1, fOutputList);
724 } // End of AliAnalysisTaskEMCALTimeCalib::UserExec()
725 
726 //________________________________________________________________________
730 {
731  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
732 
733  if(fTOFmaker) delete fTOFmaker;
734 
735 
736  if (!fOutputList)
737  {
738  AliDebug(1,"ERROR: Output list not available");
739  return;
740  }
741 } // End of AliAnalysisTaskEMCALTimeCalib::Terminate
742 
743 //________________________________________________________________________
746 {
747  //fix with noisy EMCAL fee card
748  Int_t nCells = clus->GetNCells();
749 
750  if(clus->IsEMCAL())
751  {
752  if ((clus->E() > fMaxClusterEnergy && nCells > fMaxNcells ) || nCells > fMaxNcells){
753  AliDebug(1,"very big cluster with enormous energy - cluster rejected");
754  return kFALSE;
755  }
756  }
757 
758  // remove other than photonlike
759  Double_t lambda0=clus->GetM02();
760  if (lambda0>fMaxLambda0LG || lambda0<fMinLambda0LG){
761  AliDebug(1,"lambda0 loose cut failed - cluster rejected");
762  return kFALSE;
763  }
764 
765  // remove matched clusters
766  Double_t Dx=clus->GetTrackDx();
767  Double_t Dz=clus->GetTrackDz();
768  Double_t Rtrack = TMath::Sqrt(Dx*Dx+Dz*Dz);
769  if (Rtrack <fMaxRtrack)
770  {
771  AliDebug(1,"track matched - cluster rejected");
772  return kFALSE;
773  }
774 
775  if (nCells<fMinNcells)
776  {
777  AliDebug(1,"single cell cluster - cluster rejected");
778  return kFALSE;
779  }
780 
781  if(clus->E()<fMinClusterEnergy)
782  {
783  AliDebug(1,"cluster energy < 1 GeV- cluster rejected");
784  return kFALSE;
785  }
786 
787 
788  if(!IsLowGainCellInCluster(clus)) {//no low gain cell in cluster
789  //apply more strict lambda0^2 cut
790  if (lambda0>fMaxLambda0 || lambda0<fMinLambda0){
791  AliDebug(1,"lambda0 strict cut failed - cluster rejected");
792  return kFALSE;
793  }
794  }
795 
796 
797 
798 
799  return kTRUE;
800 }//End AliAnalysisTaskEMCALTimeCalib::AcceptCluster
801 
802 //________________________________________________________________________
805  UShort_t * index = clus->GetCellsAbsId() ;
806  AliVCaloCells &cells= *(InputEvent()->GetEMCALCells());
807  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
808  if(cells.GetCellHighGain(index[i])==kFALSE) return kTRUE;//low gain cell
809  }
810  return kFALSE;
811 
812 }
813 
814 //________________________________________________________________________
816 Bool_t AliAnalysisTaskEMCALTimeCalib::CheckCellRCU(Int_t nSupMod,Int_t icol,Int_t irow)
817 {
818  Int_t iRCU;
819  if(nSupMod < 10 || (nSupMod >= 12 && nSupMod <18) )
820  {
821  if (0<=irow&&irow<8) iRCU=0; // first cable row
822  else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
823  //second cable row
824  //RCU1
825  else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
826  //second cable row
827  else if (16<=irow&&irow<24) iRCU=1; // third cable row
828 
829  if (nSupMod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
830  }
831  else
832  {
833  // Last 2 SM have one single SRU, just assign RCU 0
834  iRCU = 0 ;
835  }
836 
837  //cout<<"RCU:"<<iRCU<<endl;
838  if (iRCU<0)
839  AliFatal(Form("Wrong EMCAL/DCAL RCU number = %d\n", iRCU));
840 
841  return kTRUE;
842 }//End AliAnalysisTaskEMCALTimeCalib::CheckCellRCU
843 
844 //________________________________________________________________________
847 {
848  fMinClusterEnergy=1.0;//0.5//0.7
849  fMaxClusterEnergy=500;
850  fMinNcells=2;
851  fMaxNcells=200;
852  fMinLambda0=0.1;
853  fMaxLambda0=0.4;
854  fMinLambda0LG=0.1;
855  fMaxLambda0LG=4.0;
856  fMaxRtrack=0.025;
857  fMinCellEnergy=0.4;//0.1//0.4
858  fReferenceFileName="";//Reference.root
859  fGeometryName="";//EMCAL_COMPLETE12SMV1_DCAL_8SM
860  fPileupFromSPD=kFALSE;
861  fMinTime=-20.;
862  fMaxTime=20.;
863 }
864 
865 //________________________________________________________________________
870 void AliAnalysisTaskEMCALTimeCalib::ProduceCalibConsts(TString inputFile,TString outputFile,Bool_t isFinal)
871 {
872  TFile *file =new TFile(inputFile.Data());
873  if(file==0x0) {
874  //AliWarning("Input file does not exist!");
875  return;
876  }
877 
878  TList *list=(TList*)file->Get("chistolist");
879  if(list==0x0)
880  {
881  //AliWarning("List chistolist does not exist in file!");
882  return;
883  }
884 
885  //high gain
886  TH1F *h1[4];
887  TH1F *h2[4];
888  TH1F *h3[4];
889  TH1F *hAllTimeAvBC[4];
890  TH1F *hAllTimeRMSBC[4];
891 
892  //low gain
893  TH1F *h4[4];
894  TH1F *h5[4];
895  TH1F *h6[4];
896  TH1F *hAllTimeAvLGBC[4];
897  TH1F *hAllTimeRMSLGBC[4];
898 
899  if(isFinal==kFALSE){//first itereation
900  for(Int_t i=0;i<4;i++){
901  h1[i]=(TH1F *)list->FindObject(Form("RawTimeSumBC%d",i));
902  h2[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesBC%d",i));
903  h3[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqBC%d",i));
904 
905  h4[i]=(TH1F *)list->FindObject(Form("RawTimeSumLGBC%d",i));
906  h5[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesLGBC%d",i));
907  h6[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqLGBC%d",i));
908  }
909  } else {//final iteration
910  for(Int_t i=0;i<4;i++){
911  h1[i]=(TH1F *)list->FindObject(Form("hTimeSum%d",i));
912  h2[i]=(TH1F *)list->FindObject(Form("hTimeEnt%d",i));
913  h3[i]=(TH1F *)list->FindObject(Form("hTimeSumSq%d",i));
914 
915  h4[i]=(TH1F *)list->FindObject(Form("hTimeLGSum%d",i));
916  h5[i]=(TH1F *)list->FindObject(Form("hTimeLGEnt%d",i));
917  h6[i]=(TH1F *)list->FindObject(Form("hTimeLGSumSq%d",i));
918  }
919  }
920  //AliWarning("Input histograms read.");
921 
922  for(Int_t i=0;i<4;i++){
923  hAllTimeAvBC[i]=new TH1F(Form("hAllTimeAvBC%d",i),Form("hAllTimeAvBC%d",i),h1[i]->GetNbinsX(),h1[i]->GetXaxis()->GetXmin(),h1[i]->GetXaxis()->GetXmax());
924  hAllTimeRMSBC[i]=new TH1F(Form("hAllTimeRMSBC%d",i),Form("hAllTimeRMSBC%d",i),h3[i]->GetNbinsX(),h3[i]->GetXaxis()->GetXmin(),h3[i]->GetXaxis()->GetXmax());
925 
926  hAllTimeAvLGBC[i]=new TH1F(Form("hAllTimeAvLGBC%d",i),Form("hAllTimeAvLGBC%d",i),h4[i]->GetNbinsX(),h4[i]->GetXaxis()->GetXmin(),h4[i]->GetXaxis()->GetXmax());
927  hAllTimeRMSLGBC[i]=new TH1F(Form("hAllTimeRMSLGBC%d",i),Form("hAllTimeRMSLGBC%d",i),h6[i]->GetNbinsX(),h6[i]->GetXaxis()->GetXmin(),h6[i]->GetXaxis()->GetXmax());
928  }
929 
930  //AliWarning("New histograms booked.");
931 
932  //important remark: we use 'underflow bin' for absid=0 in OADB . That's why there is j-1 below.
933  for(Int_t i=0;i<4;i++){
934  for(Int_t j=1;j<=h1[i]->GetNbinsX();j++){
935  //high gain
936  if(h2[i]->GetBinContent(j)!=0){
937  hAllTimeAvBC[i]->SetBinContent(j-1,h1[i]->GetBinContent(j)/h2[i]->GetBinContent(j));
938  hAllTimeRMSBC[i]->SetBinContent(j-1,TMath::Sqrt(h3[i]->GetBinContent(j)/h2[i]->GetBinContent(j)) );
939  } else {
940  hAllTimeAvBC[i]->SetBinContent(j-1,0.);
941  hAllTimeRMSBC[i]->SetBinContent(j-1,0.);
942  }
943  //low gain
944  if(h5[i]->GetBinContent(j)!=0){
945  hAllTimeAvLGBC[i]->SetBinContent(j-1,h4[i]->GetBinContent(j)/h5[i]->GetBinContent(j));
946  hAllTimeRMSLGBC[i]->SetBinContent(j-1,TMath::Sqrt(h6[i]->GetBinContent(j)/h5[i]->GetBinContent(j)) );
947  } else {
948  hAllTimeAvLGBC[i]->SetBinContent(j-1,0.);
949  hAllTimeRMSLGBC[i]->SetBinContent(j-1,0.);
950  }
951 
952  }
953  }
954 
955  //AliWarning("Average and rms calculated.");
956  TFile *fileNew=new TFile(outputFile.Data(),"recreate");
957  for(Int_t i=0;i<4;i++){
958  hAllTimeAvBC[i]->Write();
959  hAllTimeRMSBC[i]->Write();
960  hAllTimeAvLGBC[i]->Write();
961  hAllTimeRMSLGBC[i]->Write();
962  }
963 
964  //AliWarning(Form("Histograms saved in %s file.",outputFile.Data()));
965 
966  fileNew->Close();
967  delete fileNew;
968  file->Close();
969  delete file;
970 
971  //AliWarning("Pointers deleted. Memory cleaned.");
972 
973 }
Bool_t SetEMCalGeometry()
Set the EMCal Geometry.
TH1F * fhRawTimeEntriesLGBC[kNBCmask]
! 4 BCmask LG
TH2F * fhEneVsAbsIdHG
! energy of each cell for high gain cells with strange time
TH1F * fhRawTimeSumLGBC[kNBCmask]
! 4 BCmask LG
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Bool_t AcceptCluster(AliVCluster *clus)
Selection criteria of good cluster are set here.
TH2F * fhRawTimeVsIdBC[kNBCmask]
! 4 BCmask HG
TH1F * fhRawTimeSumSqBC[kNBCmask]
! 4 BCmask HG
Double_t fMinLambda0
minimum cluster lambda0
TH1F * fhcalcEvtTime
! spectrum calcolot0[0]
Double_t fMinTime
minimum cluster time after correction
AliTOFT0maker * fTOFmaker
pointer to get T0 from TOF
TList * list
void SetDefaultCuts()
Set default cuts for calibration.
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
TH2F * fhTcellvsTOFT0HD
! time of cell vs TOF T0 time for higher energy threshold
Bool_t IsLowGainCellInCluster(AliVCluster *clus)
Check if low gain cell is in a cluster.
Double_t fMaxTime
maximum cluster time after correction
Double_t fMinClusterEnergy
minimum cluster energy
TString fReferenceFileName
! name of reference file
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
Task to work on Time Calibration for EMCal/DCal.
Double_t fMinLambda0LG
minimum cluster lambda0 Low Gain
Double_t fMaxClusterEnergy
maximum cluster energy
TH1F * fhEvtTimeHeader
! spectrum time from header
TH1F * fhEvtTimeDiff
! spectrum time difference
TH2F * fhTcellvsTOFT0
! time of cell vs TOF T0 time
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
TList * fOutputList
pointer to output list
TH1F * fhRawTimeEntriesBC[kNBCmask]
! 4 BCmask HG
TFile * file
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 from file.
TH1F * fhAllAverageLGBC[kNBCmask]
4 BCmask High gain
TH2F * fhTimeDsup[kNSM]
4 BCmask Low gain
TH2F * fhTimeDsupBC[kNSM][kNBCmask]
! 20 x 4