AliPhysics  b1a33bd (b1a33bd)
 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 <TH2D.h>
21 #include <TList.h>
22 #include <TCanvas.h>
23 #include <TGeoManager.h>
24 #include <TRefArray.h>
25 
26 #include "AliLog.h"
27 #include "AliAnalysisTask.h"
28 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliAODEvent.h"
31 #include "AliVEvent.h"
32 #include "AliESDInputHandler.h"
33 #include "AliAODInputHandler.h"
34 #include "AliESDpid.h"
35 #include "AliTOFcalib.h"
36 #include "AliCDBManager.h"
37 #include "AliRunTag.h"
38 
39 #include "AliTOFT0maker.h"
40 #include "AliESDCaloCluster.h"
41 #include "AliESDCaloCells.h"
42 #include "AliAODCaloCluster.h"
43 #include "AliAODCaloCells.h"
44 #include "AliEMCALGeometry.h"
45 
47 
51 
52 using std::cout;
53 using std::endl;
54 
55 //________________________________________________________________________
58 : AliAnalysisTaskSE(name),
59  fEvent(0),
60  fRunNumber(-1),
61  fTOFmaker(0),
62  fOutputList(0x0),
63  fgeom(0),
64  fGeometryName(),
65  fMinClusterEnergy(0),
66  fMaxClusterEnergy(0),
67  fMinNcells(0),
68  fMaxNcells(0),
69  fMinLambda0(0),
70  fMaxLambda0(0),
71  fMaxRtrack(0),
72  fMinCellEnergy(0),
73  fReferenceFileName(),
74  fhcalcEvtTime(0),
75  fhEventType(0),
76  fhTOFT0vsEventNumber(0),
77  fhTcellvsTOFT0(0),
78  fhTcellvsTOFT0HD(0),
79  fhTcellvsSM(0),
80  fhEneVsAbsIdHG(0),
81  fhEneVsAbsIdLG(0),
82  fhTimeVsBC(0),
83  fhTimeSumSq(),
84  fhTimeEnt(),
85  fhTimeSum(),
86  fhAllAverageBC(),
87  fhTimeDsup(),
88  fhTimeDsupBC(),
89  fhRawTimeVsIdBC(),
90  fhRawTimeSumBC(),
91  fhRawTimeEntriesBC(),
92  fhRawTimeSumSqBC(),
93  fhRawTimeVsIdLGBC(),
94  fhRawTimeSumLGBC(),
95  fhRawTimeEntriesLGBC(),
96  fhRawTimeSumSqLGBC()
97 
98 {
99  for(Int_t i = 0; i < kNBCmask; i++)
100  {
101  fhAllAverageBC[i]=0;
102  fhTimeSumSq[i]=0;
103  fhTimeEnt[i]=0;
104  fhTimeSum[i]=0;
105 
106  fhRawTimeVsIdBC[i]=0;
107  fhRawTimeSumBC[i]=0;
108  fhRawTimeEntriesBC[i]=0;
109  fhRawTimeSumSqBC[i]=0;
110 
111  fhRawTimeVsIdLGBC[i]=0;
112  fhRawTimeSumLGBC[i]=0;
114  fhRawTimeSumSqLGBC[i]=0;
115  }
116 
117  //set default cuts for calibration
118  SetDefaultCuts();
119 
120  //T0 TOF time
121  //PrepareTOFT0maker();
122 
123  // Define input and output slots here
124  // Input slot #0 works with a TChain
125  DefineInput(0, TChain::Class());
126 
127  // Output slot #0 id reserved by the base class for AOD
128  // Output slot #1 writes into a TH1 container
129  DefineOutput(1, TList::Class());
130 
131 } // End ctor
132 
133 //_____________________________________________________________________
140 {
141  //AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
142  //(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
143 
144  AliDebug(1,"*** the code version is: AliAnalysisTaskEMCALTimeCalib ***");
145  //TFile *myFile = TFile::Open("RefLHC13bcCorrected.root");
146  //TFile *myFile = TFile::Open("Reference.root");
147 
148  TFile *myFile = TFile::Open(fReferenceFileName.Data());
149 
150  if(myFile==0x0)
151  AliDebug(1,"*** NO REFERENCE FILE");
152  else
153  {
154  AliDebug(1,"*** OK TFILE");
155  // printf("*** OK TFILE \n");
156 
157  // connect ref run here
158 
159  for(Int_t i = 0; i < kNBCmask; i++)
160  {
161  fhAllAverageBC[i]=(TH1D*) myFile->Get(Form("hAllTimeAvBC%d",i));
162  }
163 
164  AliDebug(1,Form("hAllAverage entries %d", (Int_t)fhAllAverageBC[0]->GetEntries() ));
165  AliDebug(1,Form("hAllAverage entries2 %d",(Int_t)fhAllAverageBC[2]->GetEntries() ));
166 
167  }
168 } // End of AliAnalysisTaskTimeTaskMB2::Init()
169 
170 //_____________________________________________________________________
174 {
175  AliDebug(1,"AnalysisTaskEMCalGeom::Notify()");
176 
177  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
178  if (!tree) {
179  AliFatal("ERROR: Could not read chain from input slot 0");
180  } else if(TString(tree->GetName()).Contains("esdTree")) {
181  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
182  if (!esdH) {
183  AliFatal("ERROR: Could not get ESDInputHandler");
184  } else {
185  fEvent = dynamic_cast<AliESDEvent*> (esdH->GetEvent());
186  }
187  } else if(TString(tree->GetName()).Contains("aodTree")) {
188  AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
189  if (!aodH) {
190  AliFatal("ERROR: Could not get AODInputHandler");
191  } else {
192  fEvent = dynamic_cast<AliAODEvent*> (aodH->GetEvent());
193  }
194  } else {
195  AliFatal("FATAL: Only ESD or AOD. You have used other tree.");
196  }
197 
198  if (!fEvent) AliFatal("ERROR: fEvent not set");
199  else AliDebug(1,"Good, fEvent set");
200 
201 // // connect ref run here
202 // TH2F *hAllAverage = (TH2F*)hAllAverage;
203 // cout << " hAllAverage entries " << hAllAverage->GetEntries() << endl;
204 
205  // Init EMCAL geometry
206  Bool_t ok = SetEMCalGeometry();
207  if(ok) {
208  AliDebug(1,"EMCal geometry properly set...");
209  return kTRUE;
210  } else {
211  AliDebug(1,"Make sure the EMCal geometry is set properly !");
212  return kFALSE;
213  }
214 }
215 
216 //_____________________________________________________________________
219 {
220  AliDebug(1,"Read and set MisalMatrix");
221 
222  if (!fEvent)
223  {
224  AliDebug(1,"ERROR: fEvent not available");
225  return kFALSE;
226  }
227 
228  // Define EMCAL geometry to be able to read ESDs
229  // fgeom = new AliEMCALGeoUtils("EMCAL_COMPLETEv1","EMCAL");
230  // fgeom = new AliEMCALGeometry("EMCAL_COMPLETE12SMv1","EMCAL");
231  // fgeom = new AliEMCALGeometry("EMCAL_COMPLETE12SMV1_DCAL_8SM","EMCAL");
232  fgeom = new AliEMCALGeometry(fGeometryName.Data(),"EMCAL");
233  AliDebug(1,Form("Geometry defined as %s",fGeometryName.Data()));
234 
235  fRunNumber = fEvent->GetRunNumber();
236  AliDebug(1,Form("RunNumber %d", fRunNumber));
237 
238  return kTRUE;
239 }
240 
241 //_____________________________________________________________________
244 {
245  //method under development
246  AliInfo(Form("<D> -- Run # = %d", fRunNumber));
247  AliInfo("prepare TOFT0maker!!");
248 
249  AliCDBManager * cdb = AliCDBManager::Instance();
250  cdb->SetDefaultStorage("raw://");
251 
252  cdb->SetRun(fRunNumber);
253 
254  // Wonder if some have to be declared as private variables??
255  // AliESDpid *extPID = new AliESDpid();
256  // AliTOFcalib * tofCalib = new AliTOFcalib();
257  // tofCalib->SetCalibrateTOFsignal(kTRUE);
258  // tofCalib->Init();
259  // fTOFmaker = new AliTOFT0maker(extPID,tofCalib);
260  // fTOFmaker->SetTimeResolution(130.0);
261 
262  // read also the geometry..
263 
264 }// End
265 
266 //________________________________________________________________________
270 {
271  Double_t fineTmin = -500;
272  Double_t fineTmax = 400;
273  Double_t fineInterval = 0.20;
274  Int_t nfinebin = (fineTmax-fineTmin)/fineInterval;
275  //cout << "<D> nfinebin = " << nfinebin << endl;
276 
277  //book histograms
278  fhcalcEvtTime = new TH1F("fhcalcEvtTime","calculated event time from T0",nfinebin, fineTmin,fineTmax);
279  fhcalcEvtTime->GetXaxis()->SetTitle("T ");
280  fhcalcEvtTime->GetYaxis()->SetTitle("Counts (a.u.)");
281 
282  fhEventType = new TH1F("fhEventType","event type",10, 0.,10.);
283  fhEventType ->GetXaxis()->SetTitle("Type ");
284  fhEventType ->GetYaxis()->SetTitle("Counts (a.u.)");
285  fhTcellvsTOFT0 = new TH2F("hTcellvsTOFT0", " T_cell vs TOFT0", 500,-600.0,+400.0,1200,300.0,900.0);
286  fhTcellvsTOFT0HD = new TH2F("hTcellvsTOFT0HD", " T_cell vs TOFT0,HighEnergy", 500,-600.0,+400.0,4000,500.0,700.0);
287  fhTcellvsSM = new TH2F("hTcellvsSM", " T_cell vs SM", 20,0,20,300,300,900);
288  fhEneVsAbsIdHG = new TH2F("fhEneVsAbsIdHG", "energy vs ID for HG",1000,0,18000,200,0,10);
289  fhEneVsAbsIdLG = new TH2F("fhEneVsAbsIdLG", "energy vs ID for LG",1000,0,18000,100,0,20);
290 
291  for (Int_t i = 0; i < kNBCmask ; i++)
292  {
293  fhTimeSumSq[i] = new TH1F(Form("hTimeSumSq%d", i),
294  Form("cell Sum Square time BC %d ", i),
295  12289,0.,12289.);
296  fhTimeSumSq[i]->SetYTitle("Sum Sq Time ");
297  fhTimeSumSq[i]->SetXTitle("AbsId");
298 
299  fhTimeSum[i] = new TH1F(Form("hTimeSum%d", i),
300  Form("cell Sum time BC %d ", i),
301  12289,0.,12289.);
302  fhTimeSum[i]->SetYTitle("Sum Time ");
303  fhTimeSum[i]->SetXTitle("AbsId");
304 
305  fhTimeEnt[i] = new TH1F(Form("hTimeEnt%d", i),
306  Form("cell Entries BC %d ", i),
307  12289,0.,12289.);
308  fhTimeEnt[i]->SetYTitle("Entries for Time ");
309  fhTimeEnt[i]->SetXTitle("AbsId");
310 
311  //by adam
312  //high gain
313  fhRawTimeVsIdBC[i] = new TH2F(Form("RawTimeVsIdBC%d", i),
314  Form("cell raw time vs ID for high gain BC %d ", i),
315  17664,0.,17664.,600,300,900);
316  fhRawTimeVsIdBC[i]->SetXTitle("AbsId");
317  fhRawTimeVsIdBC[i]->SetYTitle("Time");
318 
319  fhRawTimeSumBC[i] = new TH1F(Form("RawTimeSumBC%d", i),
320  Form("sum of cell raw time for high gain BC %d ", i),
321  17664,0.,17664.);
322  fhRawTimeSumBC[i]->SetXTitle("AbsId");
323  fhRawTimeSumBC[i]->SetYTitle("Sum Time");
324 
325  fhRawTimeEntriesBC[i] = new TH1F(Form("RawTimeEntriesBC%d", i),
326  Form("No. entries of cells raw time for high gain BC %d ", i),
327  17664,0.,17664.);
328  fhRawTimeEntriesBC[i]->SetXTitle("AbsId");
329  fhRawTimeEntriesBC[i]->SetYTitle("Entries for Time ");
330 
331  fhRawTimeSumSqBC[i] = new TH1F(Form("RawTimeSumSqBC%d", i),
332  Form("sum of (cell raw time)^2 for high gain BC %d ", i),
333  17664,0.,17664.);
334  fhRawTimeSumSqBC[i]->SetXTitle("AbsId");
335  fhRawTimeSumSqBC[i]->SetYTitle("Sum Sq Time");
336 
337  //low gain
338  fhRawTimeVsIdLGBC[i] = new TH2F(Form("RawTimeVsIdLGBC%d", i),
339  Form("cell raw time vs ID for low gain BC %d ", i),
340  17664,0.,17664.,600,300,900);
341  fhRawTimeVsIdLGBC[i]->SetXTitle("AbsId");
342  fhRawTimeVsIdLGBC[i]->SetYTitle("Time");
343 
344  fhRawTimeSumLGBC[i] = new TH1F(Form("RawTimeSumLGBC%d", i),
345  Form("sum of cell raw time for low gain BC %d ", i),
346  17664,0.,17664.);
347  fhRawTimeSumLGBC[i]->SetXTitle("AbsId");
348  fhRawTimeSumLGBC[i]->SetYTitle("Sum Time");
349 
350  fhRawTimeEntriesLGBC[i] = new TH1F(Form("RawTimeEntriesLGBC%d", i),
351  Form("No. entries of cells raw time for low gain BC %d ", i),
352  17664,0.,17664.);
353  fhRawTimeEntriesLGBC[i]->SetXTitle("AbsId");
354  fhRawTimeEntriesLGBC[i]->SetYTitle("Entries for Time ");
355 
356  fhRawTimeSumSqLGBC[i] = new TH1F(Form("RawTimeSumSqLGBC%d", i),
357  Form("sum of (cell raw time)^2 for low gain BC %d ", i),
358  17664,0.,17664.);
359  fhRawTimeSumSqLGBC[i]->SetXTitle("AbsId");
360  fhRawTimeSumSqLGBC[i]->SetYTitle("Sum Sq Time");
361 
362 
363  for (Int_t j = 0; j < kNSM ; j++)
364  {
365  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);
366  fhTimeDsupBC[j][i]->SetYTitle(" Time (ns) ");
367  fhTimeDsupBC[j][i]->SetXTitle(" E (GeV) ");
368  }
369  }
370 
371  for (Int_t jj = 0; jj < kNSM ; jj++)
372  {
373  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);
374  fhTimeDsup[jj]->SetYTitle(" Time (ns) ");
375  fhTimeDsup[jj]->SetXTitle(" E (GeV) ");
376  }
377 
378  fhTimeVsBC = new TH2F("TimeVsBC"," SupMod time_vs_BC ", 4001,-0.5,4000.5,400,200.0,1000.0);
379 
380 
381  //add histos to list
382  fOutputList = new TList();
383 
385  fOutputList->Add(fhEventType);
388  fOutputList->Add(fhTcellvsSM);
391 
392  for (Int_t i = 0; i < kNBCmask ; i++)
393  {
394  fOutputList->Add(fhTimeSumSq[i]);
395  fOutputList->Add(fhTimeEnt[i]);
396  fOutputList->Add(fhTimeSum[i]);
397 
398  fOutputList->Add(fhRawTimeVsIdBC[i]);
399  fOutputList->Add(fhRawTimeSumBC[i]);
401  fOutputList->Add(fhRawTimeSumSqBC[i]);
402 
404  fOutputList->Add(fhRawTimeSumLGBC[i]);
407 
408  for (Int_t j = 0; j < kNSM ; j++){
409  fOutputList->Add(fhTimeDsupBC[j][i]);
410  }
411  }
412 
413  for (Int_t j = 0; j < kNSM ; j++)
414  {
415  fOutputList->Add(fhTimeDsup[j]);
416  }
417 
418  fOutputList->Add(fhTimeVsBC);
419 
420  fOutputList->SetOwner(kTRUE);
421  PostData(1,fOutputList);
422 
423 
424 } // End of AliAnalysisTaskTimeTaskMB2::UserCreateOuputObjects()
425 
426 //________________________________________________________________________
429 {
430  // Called for each event
431 
432  // Check whether ESD or AOD
433  if(AODEvent()!=0) {
434  fEvent = dynamic_cast<AliAODEvent*>(AODEvent());
435  if (!fEvent) {
436  AliError("AOD not available, exit");
437  fhEventType->Fill(1.5);
438  return;
439  }
440  } else {
441  AliVEvent* event = InputEvent();
442  //cout<<"T0TOF "<<event->GetT0TOF()<<endl;//bad idea
443  //cout<< event->GetEventNumberInFile()<<endl;
444  //cout<< event->GetTOFHeader()->GetDefaultEventTimeVal()<<endl;
445  fEvent = dynamic_cast<AliESDEvent*>(event);
446  if (!fEvent) {
447  AliError("ESD not available, exit");
448  fhEventType->Fill(0.5);
449  return;
450  }
451  }
452 
453  if(fEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
454  AliDebug(1,"Event: PileUp skip.");
455  fhEventType->Fill(2.5);
456  return;
457  }
458 
459  TString triggerclasses = fEvent->GetFiredTriggerClasses();
460  if(triggerclasses=="") {
461  fhEventType->Fill(3.5);
462  return;
463  }
464 
465  Int_t eventType = ((AliVHeader*)InputEvent()->GetHeader())->GetEventType();
466  // physics events eventType=7, select only those
467  AliDebug(1,Form("Triggerclasses %s, eventType %d",triggerclasses.Data(),eventType));
468  if(eventType != 7) {
469  fhEventType->Fill(4.5);
470  return;
471  }
472 
473  // Check trigger
474  Bool_t bMB = kFALSE;
475  Bool_t bL0 = kFALSE;
476  Bool_t bL1G = kFALSE;
477  Bool_t bL1J = kFALSE;
478 
479  if(triggerclasses.Contains("CINT7-B-NOPF-ALLNOTRD") ||
480  triggerclasses.Contains("CINT7-I-NOPF-ALLNOTRD") ||
481  triggerclasses.Contains("CINT1-I-NOPF-ALLNOTRD") ||
482  triggerclasses.Contains("CINT1-B-NOPF-ALLNOTRD") ||
483  triggerclasses.Contains("CINT8") ||
484  triggerclasses.Contains("CINT7") ||
485  triggerclasses.Contains("CPBI2_B1-B-NOPF-ALLNOTRD") ) bMB = kTRUE;
486 
487  if(triggerclasses.Contains("CEMC7-B-NOPF-CENTNOTRD") ||
488  triggerclasses.Contains("CEMC1-B-NOPF-CENTNOTRD") ||
489  triggerclasses.Contains("CEMC7") ||
490  triggerclasses.Contains("CEMC8") ||
491  triggerclasses.Contains("CEMC8-B-NOPF-CENTNOTRD") ) bL0 = kTRUE;
492 
493  if(triggerclasses.Contains("CEMC7EG1-B-NOPF-CENTNOTRD") ||
494  triggerclasses.Contains("CEMC7EG2-B-NOPF-CENTNOTRD") ||
495  triggerclasses.Contains("CEMC8EG1-B-NOPF-CENTNOTRD") ||
496  triggerclasses.Contains("CEMC8EGA") ||
497  triggerclasses.Contains("CEMC7EGA") ||
498  triggerclasses.Contains("CPBI2EGA") ) bL1G = kTRUE;
499 
500 
501  if(triggerclasses.Contains("CEMC7EJ1-B-NOPF-CENTNOTRD") ||
502  triggerclasses.Contains("CEMC7EJ2-B-NOPF-CENTNOTRD") ||
503  triggerclasses.Contains("CEMC8EJ1-B-NOPF-CENTNOTRD") ||
504  triggerclasses.Contains("CEMC7EJE") ||
505  triggerclasses.Contains("CEMC8EJE") ||
506  triggerclasses.Contains("CPBI2EJE") ) bL1J = kTRUE;
507 
508  if( bL1G || bL1J || bL0 ){ fhEventType->Fill(5.5);}
509  if( bMB ){ fhEventType->Fill(6.5);}
510 
511 
512  // if(bL1G || bL1J || bL0){
513 
514  // Prepare TOFT0 maker at the beginning of a run
515  if (fEvent->GetRunNumber() != fRunNumber)
516  {
517  fRunNumber = fEvent->GetRunNumber();
518  // PrepareTOFT0maker();
519  // cout<<"tofT0maker per run"<<fRunNumber<<endl;
520  }// fi Check if run number has changed
521 
522  // --- Use of AliTOFT0maker
523  //Double_t *calcolot;//=0.0;
524  // calcolot = fTOFmaker->ComputeT0TOF(fEvent);
525 
526  Double_t calcolot0=0.0;
527  //calcolot0=calcolot[0];
528  // cout<<"tof time "<<calcolot0<<endl;
529 
530  if (!fhcalcEvtTime)
531  {
532  AliWarning("<E> fhcalcEvtTime not available");
533  return;
534  }// fi no simple histo present
535 
536  fhcalcEvtTime->Fill(calcolot0);
537 
538  TRefArray* caloClusters = new TRefArray();
539  fEvent->GetEMCALClusters(caloClusters);
540  // cout << " ###########Bunch Cross nb = " << fEvent->GetBunchCrossNumber() << endl;
541 
542  Int_t BunchCrossNumber =fEvent->GetBunchCrossNumber();
543 
544  Float_t offset=0;
545 
546  Int_t nBC = 0;
547  nBC = BunchCrossNumber%4;
548  //Int_t nTriggerMask =fEvent->GetTriggerMask();
549  // cout << " nBC " << nBC << " nTriggerMask " << nTriggerMask<< endl;
550  Float_t timeBCoffset = 0.; //manual offest
551  // if( nBC%4 ==0 || nBC%4==1) timeBCoffset = 100.; // correction was for LHC11 when BC was not corrected
552 
553  Int_t nclus = caloClusters->GetEntries();
554  AliDebug(1,Form("###########Bunch Cross nb = %d nclus = %d",nBC,nclus ));
555  //cout << " ###########Bunch Cross nb = " << nBC <<" nclus= "<< nclus<< endl;
556  //Int_t ntracks = fEvent-> GetNumberOfTracks() ;
557 
558  for (Int_t icl = 0; icl < nclus; icl++) {
559  //ESD and AOD CaloCells carries the same information
560  AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
561  if(!AcceptCluster(clus)) continue;
562 
563  AliESDCaloCells &cells= *(dynamic_cast<AliESDCaloCells*>(fEvent->GetEMCALCells()));
564  //cout<<"nCells="<< clus->GetNCells();<<endl;
565 
566  UShort_t * index = clus->GetCellsAbsId() ;
567 
568  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
569  Int_t absId = index[i]; // or clus->GetCellNumber(i) ;
570  Float_t hkdtime = cells.GetCellTime(absId) * 1.0e09; // to get ns
571  Float_t amp = cells.GetCellAmplitude(absId) ;
572  Bool_t isHighGain = cells.GetCellHighGain(absId);
573  //cout<<"cell absID: "<<absId<<" cellTime: "<<hkdtime<<" cellaplit: "<< amp<<endl;
574  Int_t nSupMod, nModule;//,iRCU;
575  Int_t iphi, ieta, nIphi, nIeta;
576 
577 
578  //main histograms with time information
579  if(amp>fMinCellEnergy){
580  if(isHighGain){
581  fhRawTimeVsIdBC[nBC]->Fill(absId,hkdtime);
582  fhRawTimeSumBC[nBC]->Fill(absId,hkdtime);
583  fhRawTimeEntriesBC[nBC]->Fill(absId,1.);
584  fhRawTimeSumSqBC[nBC]->Fill(absId,hkdtime*hkdtime);
585  }else{
586  fhRawTimeVsIdLGBC[nBC]->Fill(absId,hkdtime);
587  fhRawTimeSumLGBC[nBC]->Fill(absId,hkdtime);
588  fhRawTimeEntriesLGBC[nBC]->Fill(absId,1.);
589  fhRawTimeSumSqLGBC[nBC]->Fill(absId,hkdtime*hkdtime);
590  }
591  }
592  //fgeom->PrintCellIndexes(absId);
593  //fgeom->PrintCellIndexes(absId,1);
594 
595  // 2015 GEOMETRY
596  fgeom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
597  fgeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
598 
599  // other histograms for cross-check
600  CheckCellRCU(nSupMod,ieta,iphi);//SM, column, row
601 
602  fhTcellvsSM->Fill(nSupMod,hkdtime);
603  if(isHighGain==kTRUE) {fhEneVsAbsIdHG->Fill(absId,amp);}
604  else {fhEneVsAbsIdLG->Fill(absId,amp);}
605 
606  fhTimeVsBC->Fill(1.*BunchCrossNumber,hkdtime-timeBCoffset);
607  if(fhAllAverageBC[nBC]!=0){
608  offset = fhAllAverageBC[nBC]->GetBinContent(absId);//comming from file after first iteration
609  //if(offset==0)cout<<"offset 0 in SM "<<nSupMod<<endl;
610  }
611  if(amp>0.5) {
612  fhTimeDsup[nSupMod]->Fill(amp,hkdtime-offset);
613  fhTimeDsupBC[nSupMod][nBC]->Fill(amp,hkdtime-offset);
614  }
615 
616  if(amp>0.9) {
617  fhTcellvsTOFT0HD->Fill(calcolot0, hkdtime);
618  }
619 
620  fhTcellvsTOFT0->Fill(calcolot0, hkdtime-offset);
621 
622  hkdtime = hkdtime-timeBCoffset;
623  Float_t hkdtimecorr;
624  hkdtimecorr= hkdtime-offset;
625 
626  if(hkdtimecorr>=-20. && hkdtimecorr<=20. && amp>0.9 ) {
627  // per cell
628  Float_t entriesTime=fhTimeEnt[nBC]->GetBinContent(absId)+1;
629  Float_t sumTimeSq=(fhTimeSumSq[nBC]->GetBinContent(absId)+(hkdtime*hkdtime));
630  Float_t sumTime=(fhTimeSum[nBC]->GetBinContent(absId)+hkdtime);
631 
632  fhTimeEnt[nBC]->SetBinContent(absId,entriesTime);
633  fhTimeSumSq[nBC]->SetBinContent(absId,sumTimeSq);
634  fhTimeSum[nBC]->SetBinContent(absId,sumTime);
635  } // hkdtime:[-20;20]
636  } // end icell
637  } //end cluster
638 
639 
640  // Post output data.
641  //cout<<"Post data and delete caloClusters"<<endl;
642  caloClusters->Delete();
643 // } // end if trigger type
644 
645  PostData(1, fOutputList);
646 } // End of AliAnalysisTaskEMCALTimeCalib::UserExec()
647 
648 //________________________________________________________________________
652 {
653  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
654  fgeom->Delete();
655 
656  if (!fOutputList)
657  {
658  AliDebug(1,"ERROR: Output list not available");
659  return;
660  }
661 } // Enf of AliAnalysisTaskEMCALTimeCalib::Terminate
662 
663 //________________________________________________________________________
666 {
667  //fix with noisy EMCAL fee card
668  Int_t nCells = clus->GetNCells();
669 
670  if(clus->IsEMCAL())
671  {
672  if ((clus->E() > fMaxClusterEnergy && nCells > fMaxNcells ) || nCells > fMaxNcells){
673  AliDebug(1,"very big cluster with enormous energy - cluster rejected");
674  return kFALSE;
675  }
676  }
677 
678  // remove other than photonlike
679  Double_t lambda0=clus->GetM02();
680  if (lambda0>fMaxLambda0 || lambda0<fMinLambda0){
681  AliDebug(1,"lambda0 failed - cluster rejected");
682  return kFALSE;
683  }
684 
685  // remove matched clusters
686  Double_t Dx=clus->GetTrackDx();
687  Double_t Dz=clus->GetTrackDz();
688  Double_t Rtrack = TMath::Sqrt(Dx*Dx+Dz*Dz);
689  if (Rtrack <fMaxRtrack)
690  {
691  AliDebug(1,"track matched - cluster rejected");
692  return kFALSE;
693  }
694 
695  if (nCells<fMinNcells)
696  {
697  AliDebug(1,"single cell cluster - cluster rejected");
698  return kFALSE;
699  }
700 
701  if(clus->E()<fMinClusterEnergy)
702  {
703  AliDebug(1,"cluster energy < 1 GeV- cluster rejected");
704  return kFALSE;
705  }
706 
707  return kTRUE;
708 }
709 
710 //________________________________________________________________________
712 Bool_t AliAnalysisTaskEMCALTimeCalib::CheckCellRCU(Int_t nSupMod,Int_t icol,Int_t irow)
713 {
714  Int_t iRCU;
715  if(nSupMod < 10 || (nSupMod >= 12 && nSupMod <18) )
716  {
717  if (0<=irow&&irow<8) iRCU=0; // first cable row
718  else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
719  //second cable row
720  //RCU1
721  else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
722  //second cable row
723  else if (16<=irow&&irow<24) iRCU=1; // third cable row
724 
725  if (nSupMod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
726  }
727  else
728  {
729  // Last 2 SM have one single SRU, just assign RCU 0
730  iRCU = 0 ;
731  }
732 
733  //cout<<"RCU:"<<iRCU<<endl;
734  if (iRCU<0)
735  AliFatal(Form("Wrong EMCAL/DCAL RCU number = %d\n", iRCU));
736 
737  return kTRUE;
738 }
739 
740 //________________________________________________________________________
743 {
744  fMinClusterEnergy=1.0;//0.5//0.7
745  fMaxClusterEnergy=500;
746  fMinNcells=2;
747  fMaxNcells=200;
748  fMinLambda0=0.1;
749  fMaxLambda0=0.4;
750  fMaxRtrack=0.025;
751  fMinCellEnergy=0.4;//0.1//0.4
752  fReferenceFileName="Reference.root";
753  fGeometryName="EMCAL_COMPLETE12SMV1_DCAL_8SM";
754 }
755 
756 //________________________________________________________________________
760 //void AliAnalysisTaskEMCALTimeCalib::ProduceCalibConsts(TString inputFile="time186319testWOL0.root",TString outputFile="Reference.root"){
761 void AliAnalysisTaskEMCALTimeCalib::ProduceCalibConsts(TString inputFile,TString outputFile)
762 {
763  TFile *file =new TFile(inputFile.Data());
764  if(file==0x0) {
765  AliWarning("Input file does not exist!");
766  return;
767  }
768 
769  TList *list=(TList*)file->Get("chistolist");
770  if(list==0x0)
771  {
772  AliWarning("List chistolist does not exist in file!");
773  return;
774  }
775 
776  //high gain
777  TH1F *h1[4];
778  TH1F *h2[4];
779  TH1F *h3[4];
780  TH1F *hAllTimeAvBC[4];
781  TH1F *hAllTimeRMSBC[4];
782 
783  //low gain
784  TH1F *h4[4];
785  TH1F *h5[4];
786  TH1F *h6[4];
787  TH1F *hAllTimeAvLGBC[4];
788  TH1F *hAllTimeRMSLGBC[4];
789 
790 
791  for(Int_t i=0;i<4;i++){
792  h1[i]=(TH1F *)list->FindObject(Form("RawTimeSumBC%d",i));
793  h2[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesBC%d",i));
794  h3[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqBC%d",i));
795 
796  h4[i]=(TH1F *)list->FindObject(Form("RawTimeSumLGBC%d",i));
797  h5[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesLGBC%d",i));
798  h6[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqLGBC%d",i));
799  }
800 
801 
802  for(Int_t i=0;i<4;i++){
803  hAllTimeAvBC[i]=new TH1F(Form("hAllTimeAvBC%d",i),Form("hAllTimeAvBC%d",i),h1[i]->GetNbinsX(),h1[i]->GetXaxis()->GetXmin(),h1[i]->GetXaxis()->GetXmax());
804  hAllTimeRMSBC[i]=new TH1F(Form("hAllTimeRMSBC%d",i),Form("hAllTimeRMSBC%d",i),h3[i]->GetNbinsX(),h3[i]->GetXaxis()->GetXmin(),h3[i]->GetXaxis()->GetXmax());
805 
806  hAllTimeAvLGBC[i]=new TH1F(Form("hAllTimeAvLGBC%d",i),Form("hAllTimeAvLGBC%d",i),h4[i]->GetNbinsX(),h4[i]->GetXaxis()->GetXmin(),h4[i]->GetXaxis()->GetXmax());
807  hAllTimeRMSLGBC[i]=new TH1F(Form("hAllTimeRMSLGBC%d",i),Form("hAllTimeRMSLGBC%d",i),h6[i]->GetNbinsX(),h6[i]->GetXaxis()->GetXmin(),h6[i]->GetXaxis()->GetXmax());
808  }
809 
810  AliWarning("Histograms booked.");
811 
812  for(Int_t i=0;i<4;i++){
813  for(Int_t j=1;j<=h1[i]->GetNbinsX();j++){
814  //high gain
815  if(h2[i]->GetBinContent(j)!=0){
816  hAllTimeAvBC[i]->SetBinContent(j,h1[i]->GetBinContent(j)/h2[i]->GetBinContent(j));
817  hAllTimeRMSBC[i]->SetBinContent(j,TMath::Sqrt(h3[i]->GetBinContent(j)/h2[i]->GetBinContent(j)) );
818  } else {
819  hAllTimeAvBC[i]->SetBinContent(j,0.);
820  hAllTimeRMSBC[i]->SetBinContent(j,0.);
821  }
822  //low gain
823  if(h5[i]->GetBinContent(j)!=0){
824  hAllTimeAvLGBC[i]->SetBinContent(j,h4[i]->GetBinContent(j)/h5[i]->GetBinContent(j));
825  hAllTimeRMSLGBC[i]->SetBinContent(j,TMath::Sqrt(h6[i]->GetBinContent(j)/h5[i]->GetBinContent(j)) );
826  } else {
827  hAllTimeAvLGBC[i]->SetBinContent(j,0.);
828  hAllTimeRMSLGBC[i]->SetBinContent(j,0.);
829  }
830 
831  }
832  }
833 
834  AliWarning("Average and rms calculated.");
835  TFile *fileNew=new TFile(outputFile.Data(),"recreate");
836  for(Int_t i=0;i<4;i++){
837  hAllTimeAvBC[i]->Write();
838  hAllTimeRMSBC[i]->Write();
839  hAllTimeAvLGBC[i]->Write();
840  hAllTimeRMSLGBC[i]->Write();
841  }
842 
843  AliWarning(Form("Histograms saved in %s file.",outputFile.Data()));
844 
845  fileNew->Close();
846  delete fileNew;
847  file->Close();
848  delete file;
849 
850  AliWarning("Pointers deleted. Memory cleaned.");
851 
852 }
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]
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
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.
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.
AliVEvent * fEvent
pointer to ESD or AOD object
Double_t fMaxClusterEnergy
maximum cluster energy
TH2F * fhTcellvsTOFT0
! time of cell vs TOF T0 time
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
TH1F * fhEventType
! spectrum calcolot0[0]
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 ProduceCalibConsts(TString inputFile="time186319testWOL0.root", TString outputFile="Reference.root")
TH1D * fhAllAverageBC[kNBCmask]
! 4 BCmask
TH2F * fhTimeDsupBC[kNSM][kNBCmask]
! 20 x 4