AliPhysics  master (3d17d9d)
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 <vector>
17 #include <TChain.h>
18 #include <TTree.h>
19 #include <TFile.h>
20 #include <TF1.h>
21 #include <TH1F.h>
22 #include <TH1D.h>
23 #include <TH2D.h>
24 #include <TH2F.h>
25 #include <TH1C.h>
26 #include <TList.h>
27 #include <TCanvas.h>
28 #include <TGeoManager.h>
29 #include <TRefArray.h>
30 #include <TKey.h>
31 
32 #include "AliLog.h"
33 #include "AliAnalysisTask.h"
34 #include "AliAnalysisManager.h"
35 #include "AliESDEvent.h"
36 #include "AliAODEvent.h"
37 #include "AliVEvent.h"
38 #include "AliESDInputHandler.h"
39 #include "AliAODInputHandler.h"
40 //#include "AliESDpid.h"
41 //#include "AliTOFcalib.h"
42 #include "AliCDBManager.h"
43 #include "AliRunTag.h"
44 
45 //#include "AliTOFT0maker.h"
46 #include "AliVCluster.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliVCaloCells.h"
49 #include "AliESDCaloCells.h"
50 #include "AliAODCaloCluster.h"
51 #include "AliAODCaloCells.h"
52 #include "AliEMCALGeometry.h"
53 #include "AliOADBContainer.h"
54 #include "AliDataFile.h"
55 
57 
61 
62 //using std::cout;
63 //using std::endl;
64 
65 //________________________________________________________________________
68 : AliAnalysisTaskSE(name),
69  fPARvec(),
70  fCurrentPARs(),
71  fCurrentPARIndex(0),
72  fIsPARRun(0),
73  fRunNumber(-1),
74  fOutputList(0x0),
75  fgeom(0),
76  fGeometryName(0),
77  fMinClusterEnergy(0),
78  fMaxClusterEnergy(0),
79  fMinNcells(0),
80  fMaxNcells(0),
81  fMinLambda0(0),
82  fMaxLambda0(0),
83  fMinLambda0LG(0),
84  fMaxLambda0LG(0),
85  fMaxRtrack(0),
86  fMinCellEnergy(0),
87  fReferenceFileName(0),
88  fReferenceRunByRunFileName(0),
89  fPileupFromSPD(kFALSE),
90  fMinTime(0),
91  fMaxTime(0),
92  fMostEneCellOnly(kFALSE),
93  fRawTimeNbins (0),
94  fRawTimeMin (0),
95  fRawTimeMax (0),
96  fPassTimeNbins(0),
97  fPassTimeMin (0),
98  fPassTimeMax (0),
99  fEnergyNbins (0),
100  fEnergyMin(0),
101  fEnergyMax(0),
102  fEnergyLGNbins (0),
103  fEnergyLGMin(0),
104  fEnergyLGMax(0),
105  fFineNbins(0),
106  fFineTmin(0),
107  fFineTmax(0),
108  fL1PhaseList(0),
109  fBadReco(kFALSE),
110  fFillHeavyHisto(kFALSE),
111  fOneHistAllBCs(kFALSE),
112  fBadChannelMapArray(0),
113  fBadChannelMapSet(kFALSE),
114  fSetBadChannelMapSource(0),
115  fBadChannelFileName(0),
116  fhcalcEvtTime(0),
117  fhEvtTimeHeader(0),
118  fhEvtTimeDiff(0),
119  fhEventType(0),
120  fhTOFT0vsEventNumber(0),
121  fhTcellvsTOFT0(0),
122  fhTcellvsTOFT0HD(0),
123  fhTcellvsSM(0),
124  fhEneVsAbsIdHG(0),
125  fhEneVsAbsIdLG(0),
126  fhTimeVsBC(0),
127  fhTimeSumSq(),
128  fhTimeEnt(),
129  fhTimeSum(),
130  fhTimeSumSqAllBCs(0x0),
131  fhTimeEntAllBCs(0x0),
132  fhTimeSumAllBCs(0x0),
133  fhTimeLGSumSq(),
134  fhTimeLGEnt(),
135  fhTimeLGSum(),
136  fhTimeLGSumSqAllBCs(0x0),
137  fhTimeLGEntAllBCs(0x0),
138  fhTimeLGSumAllBCs(0x0),
139  fhAllAverageBC(),
140  fhAllAverageLGBC(),
141  fhAllAverageAllBCs(0x0),
142  fhAllAverageLGAllBCs(0x0),
143  fhRefRuns(0),
144  fhTimeDsup(),
145  fhTimeDsupBC(),
146  fhTimeDsupLG(),
147  fhTimeDsupLGBC(),
148  fhRawTimeVsIdBC(),
149  fhRawTimeSumBC(),
150  fhRawTimeEntriesBC(),
151  fhRawTimeSumSqBC(),
152  fhRawTimeVsIdLGBC(),
153  fhRawTimeSumLGBC(),
154  fhRawTimeEntriesLGBC(),
155  fhRawTimeSumSqLGBC(),
156  fhRawTimePARs(),
157  fhRawTimeLGPARs(),
158  fhRawCorrTimeVsIdBC(),
159  fhRawCorrTimeVsIdLGBC(),
160  fhTimeVsIdBC(),
161  fhTimeVsIdLGBC(),
162  fhTimeVsIdAllBCs(0x0),
163  fhTimeVsIdLGAllBCs(0x0)
164 {
165  for(Int_t i = 0; i < kNBCmask; i++)
166  {
167  fhAllAverageBC[i]=0;
168  fhAllAverageLGBC[i]=0;
169 
170  fhTimeSumSq[i]=0;
171  fhTimeEnt[i]=0;
172  fhTimeSum[i]=0;
173 
174  fhTimeLGSumSq[i]=0;
175  fhTimeLGEnt[i]=0;
176  fhTimeLGSum[i]=0;
177 
178  fhRawTimeVsIdBC[i]=0;
179  fhRawTimeSumBC[i]=0;
180  fhRawTimeEntriesBC[i]=0;
181  fhRawTimeSumSqBC[i]=0;
182 
183  fhRawTimeVsIdLGBC[i]=0;
184  fhRawTimeSumLGBC[i]=0;
186  fhRawTimeSumSqLGBC[i]=0;
187  fhRawCorrTimeVsIdBC[i]=0;
189  fhTimeVsIdBC[i]=0;
190  fhTimeVsIdLGBC[i]=0;
191  }
192 
193  //set default cuts for calibration and geometry name
194  SetDefaultCuts();
195 
196  //T0 TOF time
198 
199  // Define input and output slots here
200  // Input slot #0 works with a TChain
201  DefineInput(0, TChain::Class());
202 
203  // Output slot #0 id reserved by the base class for AOD
204  // Output slot #1 writes into a TH1 container
205  DefineOutput(1, TList::Class());
206 
207 } // End ctor
208 
209 //_____________________________________________________________________
215 //void AliAnalysisTaskEMCALTimeCalib::LocalInit()
216 //{
217 // AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::LocalInit()");
218 //}
219 
221 //_____________________________________________________________________
223 {
224  if(fReferenceFileName.Length()!=0){
225  TFile *myFile = TFile::Open(fReferenceFileName.Data());
226  AliInfo(Form("Reference file: %s, pointer %p",fReferenceFileName.Data(),myFile));
227  if(myFile==0x0) {
228  AliFatal("*** NO REFERENCE FILE");
229  } else {
230  AliDebug(1,"*** OK TFILE");
231  // connect ref run here
232  if(!fOneHistAllBCs){
233  for(Int_t i = 0; i < kNBCmask; i++)
234  {
235  fhAllAverageBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvBC%d",i));
236  if(fhAllAverageBC[i]==0x0) AliFatal(Form("Reference histogram for BC%d does not exist",i));
237  if(fhAllAverageBC[i]->GetEntries()==0)AliWarning(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
238  fhAllAverageLGBC[i]=(TH1F*) myFile->Get(Form("hAllTimeAvLGBC%d",i));
239  if(fhAllAverageLGBC[i]==0x0) AliFatal(Form("Reference LG histogram for BC%d does not exist",i));
240  if(fhAllAverageLGBC[i]->GetEntries()==0)AliFatal(Form("fhAllAverageLGBC[%d]->GetEntries() = 0",i));
241  }
242 
243  AliDebug(1,Form("hAllAverage entries BC0 %d", (Int_t)fhAllAverageBC[0]->GetEntries() ));
244  AliDebug(1,Form("hAllAverage entries BC2 %d",(Int_t)fhAllAverageBC[2]->GetEntries() ));
245  AliDebug(1,Form("hAllAverageLG entries BC0 %d", (Int_t)fhAllAverageLGBC[0]->GetEntries() ));
246  AliDebug(1,Form("hAllAverageLG entries BC2 %d",(Int_t)fhAllAverageLGBC[2]->GetEntries() ));
247  }else{
248  fhAllAverageAllBCs=(TH1S*) myFile->Get("hAllTimeAv");
249  if(fhAllAverageAllBCs==0x0) AliFatal("Reference histogram for All BCs does not exist");
250  if(fhAllAverageAllBCs->GetEntries()==0)AliWarning("fhAllAverageLGAllBCs->GetEntries() = 0");
251  fhAllAverageLGAllBCs=(TH1S*) myFile->Get("hAllTimeAvLG");
252  if(fhAllAverageLGAllBCs==0x0) AliFatal("Reference LG histogram for all BCs does not exist");
253  if(fhAllAverageLGAllBCs->GetEntries()==0)AliFatal("fhAllAverageLGAllBCs->GetEntries() = 0");
254 
255  AliDebug(1,Form("fhAllAverageAllBCs entries %d", (Int_t)fhAllAverageAllBCs->GetEntries() ));
256  AliDebug(1,Form("fhAllAverageLGAllBCs entries %d", (Int_t)fhAllAverageLGAllBCs->GetEntries() ));
257  }
258 
259  }
260  } else { //end of reference file is provided
261  AliFatal("You require to load reference histos from file but FILENAME is not provided");
262  }
263 } // End of AliAnalysisTaskEMCALTimeCalib::LoadReferenceHistos()
264 
267 //_____________________________________________________________________
269 {
270  // connect ref run here
271  if(fReferenceRunByRunFileName.Length()!=0){
272  TFile *referenceFile = TFile::Open(fReferenceRunByRunFileName.Data());
273  if(referenceFile==0x0) {
274  AliFatal("*** NO REFERENCE R-B-R FILE");
275  return;
276  } else {
277  AliInfo(Form("Reference R-b-R file: %s, pointer %p",fReferenceRunByRunFileName.Data(),referenceFile));
278 
279  //load L1 phases to memory
280  fL1PhaseList=new TObjArray(referenceFile->GetNkeys());
281  TIter next(referenceFile->GetListOfKeys());
282  TKey *key;
283  while ((key=(TKey*)next())) {
284  fL1PhaseList->AddLast((TH1F*)referenceFile->Get(key->GetName()) );
285  //printf("key: %s points to an object of class: %s at %dn",key->GetName(),key->GetClassName(),key->GetSeekKey());
286  }
287  }
288  } else { //reference file is not provided
289  AliFatal("You require to load reference run-by-run histos from file but FILENAME is not provided");
290  return;
291  }
292 } // End of AliAnalysisTaskEMCALTimeCalib::LoadReferenceRunByRunHistos()
293 
298 {
299  fhRefRuns=NULL;
300  if(!fL1PhaseList) {
301  AliFatal("Array with reference L1 phase histograms do not exist in memory");
302  return;
303  }
304  if(fRunNumber<0) {
305  AliFatal("Negative run number");
306  return;
307  }
308 
309  fhRefRuns=(TH1C*)fL1PhaseList->FindObject(Form("h%d",fRunNumber));
310  if(fhRefRuns==0x0){
311  AliError(Form("Reference histogram for run %d does not exist. Use Default",fRunNumber));
312  fhRefRuns=(TH1C*)fL1PhaseList->FindObject("h0");
313  }
314  if(fhRefRuns==0x0) {
315  AliFatal(Form("No default histogram with L1 phases! Add default histogram to file %s!!!",fReferenceRunByRunFileName.Data()));
316  return;
317  }
318 
319  AliDebug(1,Form("Reference R-b-R histo %p, list %p, run number %d",fhRefRuns,fL1PhaseList,fRunNumber));
320  if(fhRefRuns->GetEntries()==0)AliWarning("fhRefRuns->GetEntries() = 0");
321  AliDebug(1,Form("hRefRuns entries %d", (Int_t)fhRefRuns->GetEntries() ));
322 }
323 
324 //_____________________________________________________________________
325 // Function for Setting L1 Phase reference in case of multiple phases for PAR
326 
328  if(fCurrentPARs.runNumber == 0){
329  AliFatal("fCurrentPARs not properly set! Unable to get PAR information.");
330  return;
331  }
332 
333  //if Reference is set, check if it is for correct PAR region
334  if(fhRefRuns!=0x0){
335  TString refName(fhRefRuns->GetName());
336  TString correctName;
337  if(fCurrentPARIndex==0){//before any PAR
338  correctName = Form("h%d", fRunNumber);
339  //if(fCurrentPARIndex < fCurrentPARs.numPARs){
340  }else{
341  correctName = Form("h%d_%llu", fRunNumber, fCurrentPARs.PARGlobalBCs[fCurrentPARIndex-1]);
342  }
343  if(refName.CompareTo(correctName)==0) return;
344  }
345 
346  fhRefRuns=NULL;
347  if(!fL1PhaseList) {
348  AliFatal("Array with reference L1 phase histograms do not exist in memory");
349  return;
350  }
351  if(fRunNumber<0) {
352  AliFatal("Negative run number");
353  return;
354  }
355  if(fCurrentPARIndex == 0){
356  fhRefRuns=(TH1C*)fL1PhaseList->FindObject(Form("h%d", fRunNumber));
357  }else{
358  fhRefRuns=(TH1C*)fL1PhaseList->FindObject(Form("h%d_%llu", fRunNumber, fCurrentPARs.PARGlobalBCs[fCurrentPARIndex-1]));
359  }
360 
361  if(fhRefRuns==0x0){
362  AliFatal(Form("No Reference R-b-R histo found for run %d PAR %d!", fRunNumber, fCurrentPARIndex));
363  return;
364  }
365  if(fhRefRuns->GetEntries()==0)AliWarning("fhRefRuns->GetEntries() = 0");
366  AliDebug(1,Form("hRefRuns entries %d", (Int_t)fhRefRuns->GetEntries() ));
367 }
368 
369 //_____________________________________________________________________
373 {
374  AliDebug(1,"AnalysisTaskEMCalTimeCalib::NotifyRun()");
375  AliDebug(2,Form("Notify(): EMCal geometry: fgeom = %p, fGeometryName=%s\n ",fgeom,fGeometryName.Data()));
376 
377  if (!InputEvent())
378  {
379  AliFatal("ERROR: InputEvent not set");
380  return;
381  }
382  else AliDebug(1,"Good, InputEvent set");
383 
384  // AliInfo(Form("NotifyRun, fCurrentRunnumber %d",fCurrentRunNumber));
385  fRunNumber = InputEvent()->GetRunNumber();
386  AliDebug(1,Form("RunNumber %d", fRunNumber));
387 
388  // Init EMCAL geometry
389  if (!fgeom) SetEMCalGeometry();
390  //Init EMCAL geometry done
391 
392  AliInfo(Form("Run number in NotifyRun %d",fRunNumber));
394 
395  //set L1 phases for current run
396  if(fReferenceRunByRunFileName.Length()!=0){
397  if(fIsPARRun){
399  }else{
401  }
402  }
403 
404  // set bad channel map
406 
407  return;
408 }
409 
410 //_____________________________________________________________________
413 {
414  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::SetEMCalGeometry()");
415  if(fGeometryName.Length()==0){
416  fgeom=AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
417  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fgeom->GetName(),fRunNumber));
418  } else {
419  fgeom = AliEMCALGeometry::GetInstance(fGeometryName.Data());
420  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeometryName.Data()));
421  }
422 
423  if (!fgeom){
424  AliWarning("Make sure the EMCal geometry is set properly !");
425  } else {
426  AliDebug(1,Form("EMCal geometry properly set: fGeom = %p, fGeometryName=%s",fgeom,fGeometryName.Data()));
427  }
428 
429  return kTRUE;
430 }
431 
432 //_____________________________________________________________________
435 {
436  //method under development
437  AliInfo(Form("<D> -- Run # = %d", fRunNumber));
438  AliInfo("prepare TOFT0maker!!");
439  //cout<<"Run "<< fRunNumber<<" in TOFT0maker"<<endl;
440 
441 
442  AliCDBManager * cdb = AliCDBManager::Instance();
443  cdb->SetDefaultStorage("raw://");
444  cdb->SetRun(fRunNumber);
445 
446 // AliESDpid *extPID=new AliESDpid();
447 //
448 // // Wonder if some have to be declared as private variables??
449 // // AliESDpid *extPID = new AliESDpid();
450 // // AliTOFcalib * tofCalib = new AliTOFcalib();
451 // // tofCalib->SetCalibrateTOFsignal(kTRUE);
452 // // tofCalib->Init();
453 //
454 // fTOFmaker = new AliTOFT0maker(extPID);
455 // fTOFmaker->SetTimeResolution(115.0); // if you want set the TOF res
456 // // fTOFmaker = new AliTOFT0maker(extPID,tofCalib);
457 // // fTOFmaker->SetTimeResolution(130.0);
458 //
459 // //cout<<"extPID "<<extPID<<" fTOFmaker "<<fTOFmaker<<endl;
460 
461 }// End PrepareTOFT0maker
462 
463 //________________________________________________________________________
467 {
468  AliDebug(1,"AliAnalysisTaskEMCALTimeCalib::UserCreateOutputObjects()");
469 
470  const Int_t nChannels = 17664;
471  //book histograms
472  if(fFillHeavyHisto){
473  fhcalcEvtTime = new TH1F("fhcalcEvtTime","calculated event time from T0",fFineNbins, fFineNbins,fFineTmax);
474  fhcalcEvtTime->GetXaxis()->SetTitle("T ");
475  fhcalcEvtTime->GetYaxis()->SetTitle("Counts (a.u.)");
476 
477  fhEvtTimeHeader = new TH1F("fhEvtTimeHeader","event time from header",fFineNbins, fFineNbins,fFineTmax);
478  fhEvtTimeHeader->GetXaxis()->SetTitle("T ");
479  fhEvtTimeHeader->GetYaxis()->SetTitle("Counts (a.u.)");
480 
481  fhEvtTimeDiff = new TH1F("fhEvtTimeDiff","event time difference",fFineNbins, fFineNbins,fFineTmax);
482  fhEvtTimeDiff->GetXaxis()->SetTitle("#Delta T ");
483  fhEvtTimeDiff->GetYaxis()->SetTitle("Counts (a.u.)");
484  }
485 
486  fhEventType = new TH1F("fhEventType","event type",10, 0.,10.);
487  //fhEventType ->GetXaxis()->SetTitle("Type ");
488  fhEventType->GetXaxis()->SetBinLabel(1 ,"1=No ESD");
489  fhEventType->GetXaxis()->SetBinLabel(2 ,"2=Pileup");
490  fhEventType->GetXaxis()->SetBinLabel(3 ,"3=No Trigger");
491  fhEventType->GetXaxis()->SetBinLabel(4 ,"4=Evt Type != 7");
492  fhEventType->GetXaxis()->SetBinLabel(5 ,"5=INT7,8");
493  fhEventType->GetXaxis()->SetBinLabel(6 ,"6=EMC7,8");
494  fhEventType->GetXaxis()->SetBinLabel(7 ,"7=L1 EMCal");
495  fhEventType->GetXaxis()->SetBinLabel(8 ,"8=DMC7,8");
496  fhEventType->GetXaxis()->SetBinLabel(9 ,"9=L1 DCal");
497 
498 
499  fhEventType ->GetYaxis()->SetTitle("Counts (a.u.)");
500  if(fFillHeavyHisto){
501  fhTcellvsTOFT0 = new TH2F("hTcellvsTOFT0", " T_cell vs TOFT0", 500,-600.0,+400.0,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
502  fhTcellvsTOFT0HD = new TH2F("hTcellvsTOFT0HD", " T_cell vs TOFT0,HighEnergy", 500,-600.0,+400.0,4*fRawTimeNbins,fRawTimeMin,fRawTimeMax);
503  }
504  fhTcellvsSM = new TH2F("hTcellvsSM", " T_cell vs SM", (Int_t)kNSM,0,(Double_t)kNSM,(Int_t)(fRawTimeNbins/2),fRawTimeMin,fRawTimeMax);
505 
506  if(fFillHeavyHisto){
507  fhEneVsAbsIdHG = new TH2F("fhEneVsAbsIdHG", "energy vs ID for HG",1000,0,18000,200,0,10);
508  fhEneVsAbsIdLG = new TH2F("fhEneVsAbsIdLG", "energy vs ID for LG",1000,0,18000,200,0,40);
509  }
510 
511  //Set-up Info for PAR histograms
512  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
513  if(!mgr) AliFatal("No Analysis Manager available...\n");
514  Int_t runNum = mgr->GetRunFromPath();
515  AliInfo(Form("Run number from path %d",runNum));
516  if(runNum == 0){
517  runNum = TString(gSystem->Getenv("RUNNO")).Atoi();
518  AliInfo(Form("Run number from RUNNO variable %d",runNum));
519  if(runNum < 200000){
520  AliFatal("Run Number not correctly set in UserCreateOutputObjects()!");
521  }
522  }
523  GetPARInfoForRunNumber(runNum);
524 
525  if(fIsPARRun){
526  for (Int_t iPAR = 0; iPAR <= fCurrentPARs.numPARs; iPAR++){
527  TH2F* fRawTimeSinglePAR;
528  TH2F* fRawTimeLGSinglePAR;
529  std::vector<TH2F*> vecRawTimePAR;
530  std::vector<TH2F*> vecRawTimeLGPAR;
531  for(Int_t iBC = 0; iBC < kNBCmask; iBC++){
532  fRawTimeSinglePAR = new TH2F(Form("RawTimeBeforePAR%dBC%d",iPAR+1, iBC),
533  Form("cell raw time vs ID for high gain BC %d ", iBC),
534  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
535  fRawTimeLGSinglePAR = new TH2F(Form("RawTimeLGBeforePAR%dBC%d",iPAR+1, iBC),
536  Form("cell raw time vs ID for low gain BC %d ", iBC),
537  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
538  vecRawTimePAR.push_back(fRawTimeSinglePAR);
539  vecRawTimeLGPAR.push_back(fRawTimeLGSinglePAR);
540  }
541  fhRawTimePARs.push_back(vecRawTimePAR);
542  fhRawTimeLGPARs.push_back(vecRawTimeLGPAR);
543  }
544  }
545 
546 
547  if(fOneHistAllBCs){
548  //high gain
549  fhTimeSumSqAllBCs = new TH1F("hTimeSumSqAllBCs",
550  "cell Sum Square time HG, All BCs",
551  nChannels,0.,(Double_t)nChannels);
552  fhTimeSumSqAllBCs->SetYTitle("Sum Sq Time ");
553  fhTimeSumSqAllBCs->SetXTitle("AbsId");
554 
555  fhTimeSumAllBCs = new TH1F("hTimeSumAllBCs",
556  "cell Sum time HG, All BCs",
557  nChannels,0.,(Double_t)nChannels);
558  fhTimeSumAllBCs->SetYTitle("Sum Time ");
559  fhTimeSumAllBCs->SetXTitle("AbsId");
560 
561  fhTimeEntAllBCs = new TH1F("hTimeEntAllBCs",
562  "cell Entries HG, All BCs",
563  nChannels,0.,(Double_t)nChannels);
564  fhTimeEntAllBCs->SetYTitle("Entries for Time ");
565  fhTimeEntAllBCs->SetXTitle("AbsId");
566 
567  //low gain
568  fhTimeLGSumSqAllBCs = new TH1F("hTimeLGSumSqAllBCs",
569  "cell Sum Square time LG, All BCs",
570  nChannels,0.,(Double_t)nChannels);
571  fhTimeLGSumSqAllBCs->SetYTitle("Sum Sq Time ");
572  fhTimeLGSumSqAllBCs->SetXTitle("AbsId");
573 
574  fhTimeLGSumAllBCs = new TH1F("hTimeLGSumAllBCs",
575  "cell Sum time LG, All BCs",
576  nChannels,0.,(Double_t)nChannels);
577  fhTimeLGSumAllBCs->SetYTitle("Sum Time ");
578  fhTimeLGSumAllBCs->SetXTitle("AbsId");
579 
580  fhTimeLGEntAllBCs = new TH1F("hTimeLGEntAllBCs",
581  "cell Entries LG, All BCs",
582  nChannels,0.,(Double_t)nChannels);
583  fhTimeLGEntAllBCs->SetYTitle("Entries for Time ");
584  fhTimeLGEntAllBCs->SetXTitle("AbsId");
585 
586  //histograms with corrected raw time for L1 shift and 100ns + new L1 phase
587  if(fReferenceRunByRunFileName.Length()!=0 && fFillHeavyHisto){
588  fhTimeVsIdAllBCs = new TH2F("TimeVsIdAllBCs",
589  "cell time corrected for L1 shift, 100ns and L1 phase vs ID for high gain All BCs",
590  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
591  fhTimeVsIdAllBCs->SetXTitle("AbsId");
592  fhTimeVsIdAllBCs->SetYTitle("Time");
593 
594  fhTimeVsIdLGAllBCs = new TH2F("TimeVsIdLGAllBCs",
595  "cell time corrected for L1 shift, 100ns and L1 phase vs ID for low gain All BCs",
596  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
597  fhTimeVsIdLGAllBCs->SetXTitle("AbsId");
598  fhTimeVsIdLGAllBCs->SetYTitle("Time");
599  }
600 
601  }
602 
603  for (Int_t i = 0; i < kNBCmask ; i++)
604  {
605  //already after correction
606  if(!fOneHistAllBCs){
607  //high gain
608  fhTimeSumSq[i] = new TH1F(Form("hTimeSumSq%d", i),
609  Form("cell Sum Square time HG, BC %d ", i),
610  nChannels,0.,(Double_t)nChannels);
611  fhTimeSumSq[i]->SetYTitle("Sum Sq Time ");
612  fhTimeSumSq[i]->SetXTitle("AbsId");
613 
614  fhTimeSum[i] = new TH1F(Form("hTimeSum%d", i),
615  Form("cell Sum time HG, BC %d ", i),
616  nChannels,0.,(Double_t)nChannels);
617  fhTimeSum[i]->SetYTitle("Sum Time ");
618  fhTimeSum[i]->SetXTitle("AbsId");
619 
620  fhTimeEnt[i] = new TH1F(Form("hTimeEnt%d", i),
621  Form("cell Entries HG, BC %d ", i),
622  nChannels,0.,(Double_t)nChannels);
623  fhTimeEnt[i]->SetYTitle("Entries for Time ");
624  fhTimeEnt[i]->SetXTitle("AbsId");
625 
626  //low gain
627  fhTimeLGSumSq[i] = new TH1F(Form("hTimeLGSumSq%d", i),
628  Form("cell Sum Square time LG, BC %d ", i),
629  nChannels,0.,(Double_t)nChannels);
630  fhTimeLGSumSq[i]->SetYTitle("Sum Sq Time ");
631  fhTimeLGSumSq[i]->SetXTitle("AbsId");
632 
633  fhTimeLGSum[i] = new TH1F(Form("hTimeLGSum%d", i),
634  Form("cell Sum time LG, BC %d ", i),
635  nChannels,0.,(Double_t)nChannels);
636  fhTimeLGSum[i]->SetYTitle("Sum Time ");
637  fhTimeLGSum[i]->SetXTitle("AbsId");
638 
639  fhTimeLGEnt[i] = new TH1F(Form("hTimeLGEnt%d", i),
640  Form("cell Entries LG, BC %d ", i),
641  nChannels,0.,(Double_t)nChannels);
642  fhTimeLGEnt[i]->SetYTitle("Entries for Time ");
643  fhTimeLGEnt[i]->SetXTitle("AbsId");
644  }
645 
646  //raw time histograms
647  //high gain
648  if(fFillHeavyHisto){
649  fhRawTimeVsIdBC[i] = new TH2F(Form("RawTimeVsIdBC%d", i),
650  Form("cell raw time vs ID for high gain BC %d ", i),
651  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
652  fhRawTimeVsIdBC[i]->SetXTitle("AbsId");
653  fhRawTimeVsIdBC[i]->SetYTitle("Time");
654  }
655 
656  fhRawTimeSumBC[i] = new TH1F(Form("RawTimeSumBC%d", i),
657  Form("sum of cell raw time for high gain BC %d ", i),
658  nChannels,0.,(Double_t)nChannels);
659  fhRawTimeSumBC[i]->SetXTitle("AbsId");
660  fhRawTimeSumBC[i]->SetYTitle("Sum Time");
661 
662  fhRawTimeEntriesBC[i] = new TH1F(Form("RawTimeEntriesBC%d", i),
663  Form("No. entries of cells raw time for high gain BC %d ", i),
664  nChannels,0.,(Double_t)nChannels);
665  fhRawTimeEntriesBC[i]->SetXTitle("AbsId");
666  fhRawTimeEntriesBC[i]->SetYTitle("Entries for Time ");
667 
668  fhRawTimeSumSqBC[i] = new TH1F(Form("RawTimeSumSqBC%d", i),
669  Form("sum of (cell raw time)^2 for high gain BC %d ", i),
670  nChannels,0.,(Double_t)nChannels);
671  fhRawTimeSumSqBC[i]->SetXTitle("AbsId");
672  fhRawTimeSumSqBC[i]->SetYTitle("Sum Sq Time");
673 
674  //low gain
675  if(fFillHeavyHisto){
676  fhRawTimeVsIdLGBC[i] = new TH2F(Form("RawTimeVsIdLGBC%d", i),
677  Form("cell raw time vs ID for low gain BC %d ", i),
678  nChannels,0.,(Double_t)nChannels,fRawTimeNbins,fRawTimeMin,fRawTimeMax);
679  fhRawTimeVsIdLGBC[i]->SetXTitle("AbsId");
680  fhRawTimeVsIdLGBC[i]->SetYTitle("Time");
681  }
682 
683  fhRawTimeSumLGBC[i] = new TH1F(Form("RawTimeSumLGBC%d", i),
684  Form("sum of cell raw time for low gain BC %d ", i),
685  nChannels,0.,(Double_t)nChannels);
686  fhRawTimeSumLGBC[i]->SetXTitle("AbsId");
687  fhRawTimeSumLGBC[i]->SetYTitle("Sum Time");
688 
689  fhRawTimeEntriesLGBC[i] = new TH1F(Form("RawTimeEntriesLGBC%d", i),
690  Form("No. entries of cells raw time for low gain BC %d ", i),
691  nChannels,0.,(Double_t)nChannels);
692  fhRawTimeEntriesLGBC[i]->SetXTitle("AbsId");
693  fhRawTimeEntriesLGBC[i]->SetYTitle("Entries for Time ");
694 
695  fhRawTimeSumSqLGBC[i] = new TH1F(Form("RawTimeSumSqLGBC%d", i),
696  Form("sum of (cell raw time)^2 for low gain BC %d ", i),
697  nChannels,0.,(Double_t)nChannels);
698  fhRawTimeSumSqLGBC[i]->SetXTitle("AbsId");
699  fhRawTimeSumSqLGBC[i]->SetYTitle("Sum Sq Time");
700 
701  //histograms with corrected raw time for L1 shift and 100ns
702  if(fBadReco && fFillHeavyHisto){
703  fhRawCorrTimeVsIdBC[i] = new TH2F(Form("RawCorrTimeVsIdBC%d", i),
704  Form("cell L1 shift and 100ns corrected raw time vs ID for high gain BC %d ", i),
705  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
706  fhRawCorrTimeVsIdBC[i]->SetXTitle("AbsId");
707  fhRawCorrTimeVsIdBC[i]->SetYTitle("Time");
708 
709  fhRawCorrTimeVsIdLGBC[i] = new TH2F(Form("RawCorrTimeVsIdLGBC%d", i),
710  Form("cell L1 shift and 100ns corrected raw time vs ID for low gain BC %d ", i),
711  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
712  fhRawCorrTimeVsIdLGBC[i]->SetXTitle("AbsId");
713  fhRawCorrTimeVsIdLGBC[i]->SetYTitle("Time");
714  }
715 
716  //histograms with corrected raw time for L1 shift and 100ns + new L1 phase
718  fhTimeVsIdBC[i] = new TH2F(Form("TimeVsIdBC%d", i),
719  Form("cell time corrected for L1 shift, 100ns and L1 phase vs ID for high gain BC %d ", i),
720  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
721  fhTimeVsIdBC[i]->SetXTitle("AbsId");
722  fhTimeVsIdBC[i]->SetYTitle("Time");
723 
724  fhTimeVsIdLGBC[i] = new TH2F(Form("TimeVsIdLGBC%d", i),
725  Form("cell time corrected for L1 shift, 100ns and L1 phase vs ID for low gain BC %d ", i),
726  nChannels,0.,(Double_t)nChannels,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
727  fhTimeVsIdLGBC[i]->SetXTitle("AbsId");
728  fhTimeVsIdLGBC[i]->SetYTitle("Time");
729  }
730 
731  if(!fOneHistAllBCs){
732  for (Int_t j = 0; j < kNSM ; j++)
733  {
734  //High gain
735  //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);
736  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);
737  fhTimeDsupBC[j][i]->SetYTitle(" Time (ns) ");
738  fhTimeDsupBC[j][i]->SetXTitle(" E (GeV) ");
739 
740  //low gain
741  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);
742  fhTimeDsupLGBC[j][i]->SetYTitle(" Time (ns) ");
743  fhTimeDsupLGBC[j][i]->SetXTitle(" E (GeV) ");
744  }
745  }
746 
747  }
748 
749  for (Int_t jj = 0; jj < kNSM ; jj++)
750  {
751  //high gain
752  fhTimeDsup[jj] = new TH2F(Form("SupMod%d",jj), Form("SupMod %d time_vs_E, high gain",jj),fEnergyNbins,fEnergyMin,fEnergyMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
753  fhTimeDsup[jj]->SetYTitle(" Time (ns) ");
754  fhTimeDsup[jj]->SetXTitle(" E (GeV) ");
755 
756  //low gain
757  fhTimeDsupLG[jj] = new TH2F(Form("SupMod%dLG",jj), Form("SupMod %d time_vs_E, low gain ",jj),fEnergyLGNbins,fEnergyLGMin,fEnergyLGMax,fPassTimeNbins,fPassTimeMin,fPassTimeMax);
758  fhTimeDsupLG[jj]->SetYTitle(" Time (ns) ");
759  fhTimeDsupLG[jj]->SetXTitle(" E (GeV) ");
760  }
761 
762  fhTimeVsBC = new TH2F("TimeVsBC"," SupMod time_vs_BC ", 4001,-0.5,4000.5,(Int_t)(fRawTimeNbins/2.),fRawTimeMin,fRawTimeMax);
763 
764 
765  //add histos to list
766  fOutputList = new TList();
767  fOutputList->Add(fhEventType);
768  if(fFillHeavyHisto){
772 
777  }
778  fOutputList->Add(fhTcellvsSM);
779 
780  if(fIsPARRun && fFillHeavyHisto){
781  for (Int_t iPAR = 0; iPAR <= fCurrentPARs.numPARs; iPAR++){
782  for(Int_t iBC = 0; iBC < kNBCmask; iBC++){
783  fOutputList->Add(fhRawTimePARs[iPAR][iBC]);
784  fOutputList->Add(fhRawTimeLGPARs[iPAR][iBC]);
785  }
786  }
787  }
788 
789  if(fOneHistAllBCs){
793 
797 
798  if(fReferenceRunByRunFileName.Length()!=0 && fFillHeavyHisto) {
801  }
802  }
803 
804  for (Int_t i = 0; i < kNBCmask ; i++)
805  {
806 
807  if(!fOneHistAllBCs){
808  fOutputList->Add(fhTimeSumSq[i]);
809  fOutputList->Add(fhTimeEnt[i]);
810  fOutputList->Add(fhTimeSum[i]);
811 
812  fOutputList->Add(fhTimeLGSumSq[i]);
813  fOutputList->Add(fhTimeLGEnt[i]);
814  fOutputList->Add(fhTimeLGSum[i]);
815  }
816 
817  if(fFillHeavyHisto) {
818  fOutputList->Add(fhRawTimeVsIdBC[i]);
820  }
821  fOutputList->Add(fhRawTimeSumBC[i]);
823  fOutputList->Add(fhRawTimeSumSqBC[i]);
824 
825  fOutputList->Add(fhRawTimeSumLGBC[i]);
828 
829  if(fBadReco && fFillHeavyHisto) {
832  }
834  fOutputList->Add(fhTimeVsIdBC[i]);
835  fOutputList->Add(fhTimeVsIdLGBC[i]);
836  }
837 
838  if(!fOneHistAllBCs) {
839  for (Int_t j = 0; j < kNSM ; j++){
840  fOutputList->Add(fhTimeDsupBC[j][i]);
841  fOutputList->Add(fhTimeDsupLGBC[j][i]);
842  }
843  }
844  }
845 
846  for (Int_t j = 0; j < kNSM ; j++)
847  {
848  fOutputList->Add(fhTimeDsup[j]);
849  fOutputList->Add(fhTimeDsupLG[j]);
850  }
851 
852  fOutputList->Add(fhTimeVsBC);
853 
854  fOutputList->SetOwner(kTRUE);
855  PostData(1,fOutputList);
856 
857 
858 } // End of AliAnalysisTaskEMCALTimeCalib::UserCreateOuputObjects()
859 
860 //________________________________________________________________________
863 {
864  // Called for each event
865  AliDebug(2,Form("UserExec: EMCal geometry: fgeom = %p fGeometryName %s",fgeom,fGeometryName.Data()));
866  AliVEvent *event = InputEvent();
867  //cout<<"T0TOF "<<event->GetT0TOF()<<endl;//bad idea
868  //cout<< fEvent->GetTOFHeader()->GetDefaultEventTimeVal()<<endl;
869  AliDebug(2,Form("TOF time from header %f ps",event->GetTOFHeader()->GetDefaultEventTimeVal()));
870  if(fFillHeavyHisto) fhEvtTimeHeader->Fill(event->GetTOFHeader()->GetDefaultEventTimeVal());
871 
872  //fEvent = dynamic_cast<AliESDEvent*>(event);
873  if (!event) {
874  AliError("ESD not available, exit");
875  fhEventType->Fill(0.5);
876  return;
877  }
878 
879  if(fPileupFromSPD==kTRUE){
880  if(event->IsPileupFromSPD(3,0.8,3.,2.,5.)){
881  AliDebug(1,"Event: PileUp skip.");
882  fhEventType->Fill(1.5);
883  return;
884  }
885  }
886 
887  TString triggerclasses = event->GetFiredTriggerClasses();
888  if(triggerclasses=="") {
889  fhEventType->Fill(2.5);
890  return;
891  }
892 
893  Int_t eventType = ((AliVHeader*)event->GetHeader())->GetEventType();
894  // physics events eventType=7, select only those
895  AliDebug(1,Form("Triggerclasses %s, eventType %d",triggerclasses.Data(),eventType));
896  if(eventType != 7) {
897  fhEventType->Fill(3.5);
898  return;
899  }
900 
901  // Check trigger
902  Bool_t bMB = kFALSE;
903  Bool_t bL0 = kFALSE;
904  Bool_t bL1G = kFALSE;
905  Bool_t bL1J = kFALSE;
906  Bool_t bDL0 = kFALSE;
907  Bool_t bDL1G = kFALSE;
908  Bool_t bDL1J = kFALSE;
909 
910  if(triggerclasses.Contains("CINT7-B-NOPF-ALLNOTRD") ||
911  triggerclasses.Contains("CINT7-I-NOPF-ALLNOTRD") ||
912  triggerclasses.Contains("CINT1-I-NOPF-ALLNOTRD") ||
913  triggerclasses.Contains("CINT1-B-NOPF-ALLNOTRD") ||
914  triggerclasses.Contains("CINT8") ||
915  triggerclasses.Contains("CINT7") ||
916  triggerclasses.Contains("CPBI2_B1-B-NOPF-ALLNOTRD") ) bMB = kTRUE;
917 
918  if(triggerclasses.Contains("CEMC7-B-NOPF-CENTNOTRD") ||
919  triggerclasses.Contains("CEMC1-B-NOPF-CENTNOTRD") ||
920  triggerclasses.Contains("CEMC7") ||
921  triggerclasses.Contains("CEMC8") ||
922  triggerclasses.Contains("CEMC8-B-NOPF-CENTNOTRD") ) bL0 = kTRUE;
923 
924  if(triggerclasses.Contains("CDMC7-B-NOPF-CENTNOTRD") ||
925  triggerclasses.Contains("CDMC1-B-NOPF-CENTNOTRD") ||
926  triggerclasses.Contains("CDMC7") ||
927  triggerclasses.Contains("CDMC8") ||
928  triggerclasses.Contains("CDMC8-B-NOPF-CENTNOTRD") ) bDL0 = kTRUE;
929 
930  if(triggerclasses.Contains("CEMC7EG1-B-NOPF-CENTNOTRD") ||
931  triggerclasses.Contains("CEMC7EG2-B-NOPF-CENTNOTRD") ||
932  triggerclasses.Contains("CEMC8EG1-B-NOPF-CENTNOTRD") ||
933  triggerclasses.Contains("CEMC8EGA") ||
934  triggerclasses.Contains("CEMC7EGA") ||
935  triggerclasses.Contains("CEMC7EG1-B") ||
936  triggerclasses.Contains("CEMC7EG2-B") ||
937  triggerclasses.Contains("CPBI2EGA") ) bL1G = kTRUE;
938 
939  if(triggerclasses.Contains("CDMC7DG1-B-NOPF-CENTNOTRD") ||
940  triggerclasses.Contains("CDMC7DG2-B-NOPF-CENTNOTRD") ||
941  triggerclasses.Contains("CDMC8DG1-B-NOPF-CENTNOTRD") ||
942  triggerclasses.Contains("CDMC8DGA") ||
943  triggerclasses.Contains("CDMC7DGA") ||
944  triggerclasses.Contains("CDMC7DG1-B") ||
945  triggerclasses.Contains("CDMC7DG2-B") ||
946  triggerclasses.Contains("CPBI2DGA") ) bDL1G = kTRUE;
947 
948  if(triggerclasses.Contains("CEMC7EJ1-B-NOPF-CENTNOTRD") ||
949  triggerclasses.Contains("CEMC7EJ2-B-NOPF-CENTNOTRD") ||
950  triggerclasses.Contains("CEMC8EJ1-B-NOPF-CENTNOTRD") ||
951  triggerclasses.Contains("CEMC7EJE") ||
952  triggerclasses.Contains("CEMC8EJE") ||
953  triggerclasses.Contains("CEMC7EJ1-B") ||
954  triggerclasses.Contains("CEMC7EJ2-B") ||
955  triggerclasses.Contains("CPBI2EJE") ) bL1J = kTRUE;
956 
957  if(triggerclasses.Contains("CDMC7DJ1-B-NOPF-CENTNOTRD") ||
958  triggerclasses.Contains("CDMC7DJ2-B-NOPF-CENTNOTRD") ||
959  triggerclasses.Contains("CDMC8DJ1-B-NOPF-CENTNOTRD") ||
960  triggerclasses.Contains("CDMC7DJE") ||
961  triggerclasses.Contains("CDMC8DJE") ||
962  triggerclasses.Contains("CDMC7DJ1-B") ||
963  triggerclasses.Contains("CDMC7DJ2-B") ||
964  triggerclasses.Contains("CPBI2DJE") ) bDL1J = kTRUE;
965 
966  if( bMB ){ fhEventType->Fill(4.5);}//INT7,8
967  if( bL0 ){ fhEventType->Fill(5.5);}//EMC7,EMC8
968  if( bL1G || bL1J ){ fhEventType->Fill(6.5);}//L1 EMCal
969  if( bDL0 ){ fhEventType->Fill(7.5);}//DMC7,DMC8
970  if( bDL1G || bDL1J ){ fhEventType->Fill(8.5);}//L1 DCal
971 
972  // if(bL1G || bL1J || bL0){
973 
974 // Prepare TOFT0 maker at the beginning of a run
975 // if (event->GetRunNumber() != fRunNumber){
976 // AliInfo(Form("Runno per event %d",event->GetRunNumber()));
977 // fRunNumber = event->GetRunNumber();
978 // // PrepareTOFT0maker();
979 // // cout<<"tofT0maker per run"<<fRunNumber<<endl;
980 // }// fi Check if run number has changed
981 
982 // // --- Use of AliTOFT0maker
983 // Double_t calcolot0=0.0;
984 // if(!AODEvent()){
985 // Double_t* timeTOFtable;
986 // timeTOFtable=fTOFmaker->ComputeT0TOF(dynamic_cast<AliESDEvent*>(event));
987 // AliDebug(2,Form("TOF time %f ps, resolution %f ps, tracks at TOF %f/used %f",timeTOFtable[0],timeTOFtable[1],timeTOFtable[3],timeTOFtable[7]));
988 // //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;
989 // calcolot0=timeTOFtable[0];
990 // }
991 
992 // if(fFillHeavyHisto) {
993 // fhcalcEvtTime->Fill(calcolot0);
994 // if(calcolot0 != 0 && event->GetTOFHeader()->GetDefaultEventTimeVal() != 0 )
995 // fhEvtTimeDiff->Fill(calcolot0-event->GetTOFHeader()->GetDefaultEventTimeVal());
996 // }
997 
998  TRefArray* caloClusters = new TRefArray();
999  event->GetEMCALClusters(caloClusters);
1000  // cout << " ###########Bunch Cross nb = " << event->GetBunchCrossNumber() << endl;
1001 
1002  Int_t BunchCrossNumber =event->GetBunchCrossNumber();
1003 
1004  Float_t offset=0.;
1005  Float_t offsetPerSM=0.;
1006  Int_t L1phaseshift=0;
1007  Int_t L1phase=0;
1008  Int_t L1shiftOffset=0;
1009 
1010  Int_t nBC = 0;
1011  nBC = BunchCrossNumber%4;
1012  //Int_t nTriggerMask =event->GetTriggerMask();
1013  // cout << " nBC " << nBC << " nTriggerMask " << nTriggerMask<< endl;
1014  Float_t timeBCoffset = 0.; //manual offset
1015  // if( nBC%4 ==0 || nBC%4==1) timeBCoffset = 100.; // correction was for LHC11 when BC was not corrected
1016 
1017  Int_t nclus = caloClusters->GetEntries();
1018  AliDebug(1,Form("###########Bunch Cross nb = %d nclus = %d",nBC,nclus ));
1019  //cout << " ###########Bunch Cross nb = " << nBC <<" nclus= "<< nclus<< endl;
1020  //Int_t ntracks = event-> GetNumberOfTracks() ;
1021 
1022  AliVCaloCells &cells= *(event->GetEMCALCells());//it is cluster independent
1023  //Variables used plenty in loops
1024  Int_t nSupMod=-1, nModule=-1;
1025  Int_t iphi=-1, ieta=-1, nIphi=-1, nIeta=-1;
1026  Int_t absId=-1;
1027  Float_t hkdtime=0.0;
1028  Float_t amp=0.0;
1029  Bool_t isHighGain=kTRUE;
1030  Int_t mostEneId=-1;
1031  Float_t mostEneEn=0.;
1032 
1033  fCurrentPARIndex = 0;//before any PAR
1034  if(fIsPARRun){
1035  ULong64_t eventBC = (ULong64_t)event->GetBunchCrossNumber();
1036  ULong64_t eventOrbit = ((ULong64_t)(3564))*((ULong64_t)event->GetOrbitNumber());
1037  ULong64_t eventPeriod = ((ULong64_t)(59793994260))*((ULong64_t)(event->GetPeriodNumber()));
1038  //ULong64_t globalBC = event->GetBunchCrossNumber() + 3564*event->GetOrbitNumber() + 59793994260*event->GetPeriodNumber();
1039  ULong64_t globalBC = eventBC + eventOrbit + eventPeriod;
1040  for(int ipar = 0; ipar < fCurrentPARs.numPARs; ipar++){
1041  if(globalBC >= fCurrentPARs.PARGlobalBCs[ipar]){
1042  fCurrentPARIndex ++;
1043  }
1044  }
1045  }
1046  if(fReferenceRunByRunFileName.Length()!=0 && fIsPARRun){
1048  }
1049  for (Int_t icl = 0; icl < nclus; icl++) {
1050  //ESD and AOD CaloCells carries the same information
1051  AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
1052  if(!AcceptCluster(clus)) continue;
1053 
1054  //cout<<"nCells="<< clus->GetNCells();<<endl;
1055 
1056  UShort_t * index = clus->GetCellsAbsId() ;
1057 
1058  // find index of the most energetic cell in cluster
1059  mostEneEn=0.;
1060  mostEneId=-1;
1061  if(fMostEneCellOnly) {
1062  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
1063  absId = index[i];
1064  amp = cells.GetCellAmplitude(absId) ;
1065  if(amp > mostEneEn){
1066  mostEneEn = amp;
1067  mostEneId = absId;
1068  }
1069  }
1070  }//works only for fMostEneCellOnly=kTRUE
1071 
1072  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
1073  absId = index[i]; // or clus->GetCellNumber(i) ;
1074  if(fMostEneCellOnly && absId != mostEneId) {
1075  //printf("tr.%s.cl.%d.cell.%d.rejected\n",triggerclasses.Data(),icl,i);
1076  continue;
1077  }
1078  //printf("tr.%s.cl.%d.cell.%d.accepted\n",triggerclasses.Data(),icl,i);
1079  hkdtime = cells.GetCellTime(absId) * 1.0e09; // to get ns
1080  amp = cells.GetCellAmplitude(absId) ;
1081  isHighGain = cells.GetCellHighGain(absId);
1082  //cout<<"cell absID: "<<absId<<" cellTime: "<<hkdtime<<" cellaplit: "<< amp<<endl;
1083  // GEOMETRY tranformations
1084  fgeom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1085  fgeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1086 
1087  //bad channel check. 0: good channel, 1-5: bad channel
1088  if(fSetBadChannelMapSource==1){
1089  if(GetEMCALChannelStatus(nSupMod,ieta,iphi)) continue;//printf("bad\n");
1090  } else if(fSetBadChannelMapSource==2){
1091  if(GetEMCALChannelStatus(absId)) continue;//printf("bad\n");
1092  }
1093 
1094  //main histograms with raw time information
1095  if(amp>fMinCellEnergy){
1096 
1097  if(isHighGain){
1098  if(fFillHeavyHisto){
1099  fhRawTimeVsIdBC[nBC]->Fill(absId,hkdtime);
1100  if(fIsPARRun){
1101  if(fhRawTimePARs[fCurrentPARIndex][nBC]==0x0)AliFatal(Form("No Histogram for PAR number %d! Problem with Create Output Objects", fCurrentPARIndex));
1102  fhRawTimePARs[fCurrentPARIndex][nBC]->Fill(absId, hkdtime);
1103  }
1104  }
1105  fhRawTimeSumBC[nBC]->Fill(absId,hkdtime);
1106  fhRawTimeEntriesBC[nBC]->Fill(absId,1.);
1107  fhRawTimeSumSqBC[nBC]->Fill(absId,hkdtime*hkdtime);
1108  }else{
1109  if(fFillHeavyHisto){
1110  fhRawTimeVsIdLGBC[nBC]->Fill(absId,hkdtime);
1111  if(fIsPARRun){
1112  fhRawTimeLGPARs[fCurrentPARIndex][nBC]->Fill(absId, hkdtime);
1113  }
1114  }
1115  fhRawTimeSumLGBC[nBC]->Fill(absId,hkdtime);
1116  fhRawTimeEntriesLGBC[nBC]->Fill(absId,1.);
1117  fhRawTimeSumSqLGBC[nBC]->Fill(absId,hkdtime*hkdtime);
1118  }
1119  }
1120  //fgeom->PrintCellIndexes(absId);
1121  //fgeom->PrintCellIndexes(absId,1);
1122 
1123  // other histograms for cross-check
1124  CheckCellRCU(nSupMod,ieta,iphi);//SM, column, row
1125 
1126  fhTcellvsSM->Fill(nSupMod,hkdtime);
1127  if(fFillHeavyHisto) {
1128  if(isHighGain==kTRUE) {fhEneVsAbsIdHG->Fill(absId,amp);}
1129  else {fhEneVsAbsIdLG->Fill(absId,amp);}
1130  }
1131  fhTimeVsBC->Fill(1.*BunchCrossNumber,hkdtime-timeBCoffset);
1132  //important remark: We use 'Underflow bin' for absid=0 in OADB for time calibration
1133  if(!fOneHistAllBCs){
1134  if(isHighGain==kTRUE){
1135  if(fhAllAverageBC[nBC]!=0) {//comming from file after the first iteration
1136  offset = (Float_t)(fhAllAverageBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
1137  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
1138  AliFatal(Form("Reference histogram for BC%d not properly loaded",nBC));
1139  }
1140  } else {
1141  if(fhAllAverageLGBC[nBC]!=0) {//comming from file after the first iteration
1142  offset = (Float_t)(fhAllAverageLGBC[nBC]->GetBinContent(absId));//channel absId=0 has histogram bin=0
1143  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
1144  AliFatal(Form("Reference LG histogram for BC%d not properly loaded",nBC));
1145  }
1146  }
1147  }else{
1148  if(isHighGain==kTRUE){
1149  if(fhAllAverageAllBCs!=0) {//comming from file after the first iteration
1150  offset = (Float_t)(fhAllAverageAllBCs->GetBinContent(absId));//channel absId=0 has histogram bin=0
1151  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
1152  AliFatal("Reference histogram for all BCs not properly loaded");
1153  }
1154  } else {
1155  if(fhAllAverageLGAllBCs!=0) {//comming from file after the first iteration
1156  offset = (Float_t)(fhAllAverageLGAllBCs->GetBinContent(absId));//channel absId=0 has histogram bin=0
1157  } else if(fReferenceFileName.Length()!=0){//protection against missing reference histogram
1158  AliFatal("Reference LG histogram for all BCs not properly loaded");
1159  }
1160  }
1161  }
1162  //if(offset==0)cout<<"offset 0 in SM "<<nSupMod<<endl;
1163 
1164  // Solution for 2015 data where L1 phase and L1 shift is not correct in data. We need to calibrate run by run.
1165  // The shift and phase are done per SM (0-19).
1166  // L1 phase is necessary in run 2.
1167  // L1 shift is necessary only for bad reconstructed runs (muon_calo_pass1 lhc15f-m)
1168  if(fhRefRuns!=0) {//comming from file after the first iteration
1169  L1phaseshift = (Int_t)(fhRefRuns->GetBinContent(nSupMod));//SM0 = bin0
1170 
1171  // to correct for L1 phase
1172  // this part works for both: muon_calo_pass1 of LHC15n (pp@2.76) and later reconstructions
1173  // wrong reconstruction done before in run2
1174  L1phase = L1phaseshift & 3; //bit operation
1175  if(nBC >= L1phase)
1176  offsetPerSM = (nBC - L1phase)*25;
1177  else
1178  offsetPerSM = (nBC - L1phase + 4)*25;
1179 
1180  // to correct for L1 shift
1181  // this part is only for wrongly reconstructed runs before LHC15n in run2
1182  if(fBadReco){
1183  L1shiftOffset = L1phaseshift>>2; //bit operation
1184  L1shiftOffset*=25;
1185  //(we subtract it here because we subtract the whole wrong offset later --=+)
1186  if(nBC==0 || nBC==1) L1shiftOffset-=100.;//additional shift for muon_calo_pass1 up to lhc15f-m
1187  }
1188  } else if(fReferenceRunByRunFileName.Length()!=0){//protection against missing reference histogram
1189  AliFatal("Reference histogram run-by-run not properly loaded");
1190  }
1191  //end of load additional offset
1192 
1193  //fill the raw time with L1 shift correction and 100ns
1194  if(fBadReco && fFillHeavyHisto && amp>fMinCellEnergy){
1195  if(isHighGain){
1196  fhRawCorrTimeVsIdBC[nBC]->Fill(absId,hkdtime-L1shiftOffset);
1197  }else{
1198  fhRawCorrTimeVsIdLGBC[nBC]->Fill(absId,hkdtime-L1shiftOffset);
1199  }
1200  }
1201 
1202  //fill time after L1 shift correction and 100ns and new L1 phase
1204  if(isHighGain){
1205  fhTimeVsIdBC[nBC]->Fill(absId,hkdtime-offset-offsetPerSM-L1shiftOffset);
1206  }else{
1207  fhTimeVsIdLGBC[nBC]->Fill(absId,hkdtime-offset-offsetPerSM-L1shiftOffset);
1208  }
1209  }
1210 
1212  if(isHighGain){
1213  fhTimeVsIdAllBCs->Fill(absId,hkdtime-offset-offsetPerSM-L1shiftOffset);
1214  }else{
1215  fhTimeVsIdLGAllBCs->Fill(absId,hkdtime-offset-offsetPerSM-L1shiftOffset);
1216  }
1217  }
1218 
1219  //other control histograms
1220  if(amp>0.5) {
1221  if(isHighGain){
1222  fhTimeDsup[nSupMod]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
1223  if(!fOneHistAllBCs) fhTimeDsupBC[nSupMod][nBC]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
1224  }else{
1225  fhTimeDsupLG[nSupMod]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
1226  if(!fOneHistAllBCs) fhTimeDsupLGBC[nSupMod][nBC]->Fill(amp,hkdtime-offset-offsetPerSM-L1shiftOffset);
1227  }
1228  }
1229 
1230 // if(fFillHeavyHisto) {
1231 // if(amp>0.9) {
1232 // fhTcellvsTOFT0HD->Fill(calcolot0, hkdtime);
1233 // }
1234 // fhTcellvsTOFT0->Fill(calcolot0, hkdtime-offset-offsetPerSM-L1shiftOffset);
1235 // }
1236 
1237  hkdtime = hkdtime-timeBCoffset;//time corrected by manual offset (default=0)
1238  Float_t hkdtimecorr;
1239  hkdtimecorr= hkdtime-offset-offsetPerSM-L1shiftOffset;//time after first iteration
1240 
1241  //main histograms after the first itereation for calibration constants
1242  //if(hkdtimecorr>=-20. && hkdtimecorr<=20. && amp>0.9 ) {
1243  if(hkdtimecorr>=fMinTime && hkdtimecorr<=fMaxTime && amp>fMinCellEnergy ) {
1244  // per cell
1245 // Float_t entriesTime=fhTimeEnt[nBC]->GetBinContent(absId)+1;
1246 // Float_t sumTimeSq=(fhTimeSumSq[nBC]->GetBinContent(absId)+(hkdtime*hkdtime));
1247 // Float_t sumTime=(fhTimeSum[nBC]->GetBinContent(absId)+hkdtime);
1248 //
1249 // fhTimeEnt[nBC]->SetBinContent(absId,entriesTime);
1250 // fhTimeSumSq[nBC]->SetBinContent(absId,sumTimeSq);
1251 // fhTimeSum[nBC]->SetBinContent(absId,sumTime);
1252 
1253  //correction in 2015 for wrong L1 phase and L1 shift
1254  hkdtime = hkdtime - offsetPerSM - L1shiftOffset;
1255 
1256  if(!fOneHistAllBCs){
1257  if(isHighGain){
1258  fhTimeEnt[nBC]->Fill(absId,1.);
1259  fhTimeSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
1260  fhTimeSum[nBC]->Fill(absId,hkdtime);
1261  }else{
1262  fhTimeLGEnt[nBC]->Fill(absId,1.);
1263  fhTimeLGSumSq[nBC]->Fill(absId,hkdtime*hkdtime);
1264  fhTimeLGSum[nBC]->Fill(absId,hkdtime);
1265  }
1266  }else{
1267  if(isHighGain){
1268  fhTimeEntAllBCs->Fill(absId,1.);
1269  fhTimeSumSqAllBCs->Fill(absId,hkdtime*hkdtime);
1270  fhTimeSumAllBCs->Fill(absId,hkdtime);
1271  }else{
1272  fhTimeLGEntAllBCs->Fill(absId,1.);
1273  fhTimeLGSumSqAllBCs->Fill(absId,hkdtime*hkdtime);
1274  fhTimeLGSumAllBCs->Fill(absId,hkdtime);
1275  }
1276  }
1277 
1278 
1279  } // hkdtime:[-20;20]
1280  } // end icell
1281  } //end cluster
1282 
1283 
1284  // Post output data.
1285  //cout<<"Post data and delete caloClusters"<<endl;
1286  caloClusters->Delete();
1287  delete caloClusters;
1288 // } // end if trigger type
1289 
1290  PostData(1, fOutputList);
1291 } // End of AliAnalysisTaskEMCALTimeCalib::UserExec()
1292 
1293 //________________________________________________________________________
1297 {
1298  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
1299 
1300  // if(fTOFmaker) delete fTOFmaker;
1301 
1302  if(fL1PhaseList) {
1303  fL1PhaseList->SetOwner();
1304  fL1PhaseList->Clear();
1305  delete fL1PhaseList;
1306  }
1307 
1308  if (fBadChannelMapArray) {
1309  fBadChannelMapArray->Clear();
1310  delete fBadChannelMapArray;
1311  }
1312 
1313  if (!fOutputList)
1314  {
1315  AliDebug(1,"ERROR: Output list not available");
1316  return;
1317  }
1318 } // End of AliAnalysisTaskEMCALTimeCalib::Terminate
1319 
1320 //________________________________________________________________________
1323 {
1324  //fix with noisy EMCAL fee card
1325  Int_t nCells = clus->GetNCells();
1326 
1327  if(clus->IsEMCAL())
1328  {
1329  if ((clus->E() > fMaxClusterEnergy && nCells > fMaxNcells ) || nCells > fMaxNcells){
1330  AliDebug(1,"very big cluster with enormous energy - cluster rejected");
1331  return kFALSE;
1332  }
1333  }
1334 
1335  // remove other than photonlike
1336  Double_t lambda0=clus->GetM02();
1337  if (lambda0>fMaxLambda0LG || lambda0<fMinLambda0LG){
1338  AliDebug(1,"lambda0 loose cut failed - cluster rejected");
1339  return kFALSE;
1340  }
1341 
1342  // remove matched clusters
1343  Double_t Dx=clus->GetTrackDx();
1344  Double_t Dz=clus->GetTrackDz();
1345  Double_t Rtrack = TMath::Sqrt(Dx*Dx+Dz*Dz);
1346  if (Rtrack <fMaxRtrack)
1347  {
1348  AliDebug(1,"track matched - cluster rejected");
1349  return kFALSE;
1350  }
1351 
1352  if (nCells<fMinNcells)
1353  {
1354  AliDebug(1,"single cell cluster - cluster rejected");
1355  return kFALSE;
1356  }
1357 
1358  if(clus->E()<fMinClusterEnergy)
1359  {
1360  AliDebug(1,"cluster energy < 1 GeV- cluster rejected");
1361  return kFALSE;
1362  }
1363 
1364 
1365  if(!IsLowGainCellInCluster(clus)) {//no low gain cell in cluster
1366  //apply more strict lambda0^2 cut
1367  if (lambda0>fMaxLambda0 || lambda0<fMinLambda0){
1368  AliDebug(1,"lambda0 strict cut failed - cluster rejected");
1369  return kFALSE;
1370  }
1371  }
1372 
1373 
1374 
1375 
1376  return kTRUE;
1377 }//End AliAnalysisTaskEMCALTimeCalib::AcceptCluster
1378 
1379 //________________________________________________________________________
1382  UShort_t * index = clus->GetCellsAbsId() ;
1383  AliVCaloCells &cells= *(InputEvent()->GetEMCALCells());
1384  for(Int_t i = 0; i < clus->GetNCells() ; i++) {
1385  if(cells.GetCellHighGain(index[i])==kFALSE) return kTRUE;//low gain cell
1386  }
1387  return kFALSE;
1388 
1389 }
1390 
1391 //________________________________________________________________________
1394 {
1395  Int_t iRCU;
1396  if(nSupMod < 10 || (nSupMod >= 12 && nSupMod <18) )
1397  {
1398  if (0<=irow&&irow<8) iRCU=0; // first cable row
1399  else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
1400  //second cable row
1401  //RCU1
1402  else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
1403  //second cable row
1404  else if (16<=irow&&irow<24) iRCU=1; // third cable row
1405 
1406  if (nSupMod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1407  }
1408  else
1409  {
1410  // Last 2 SM have one single SRU, just assign RCU 0
1411  iRCU = 0 ;
1412  }
1413 
1414  //cout<<"RCU:"<<iRCU<<endl;
1415  if (iRCU<0)
1416  AliFatal(Form("Wrong EMCAL/DCAL RCU number = %d\n", iRCU));
1417 
1418  return kTRUE;
1419 }//End AliAnalysisTaskEMCALTimeCalib::CheckCellRCU
1420 
1421 //________________________________________________________________________
1424 {
1425  fMinClusterEnergy=1.0;//0.5//0.7
1426  fMaxClusterEnergy=500;
1427  fMinNcells=2;
1428  fMaxNcells=200;
1429  fMinLambda0=0.1;
1430  fMaxLambda0=0.4;
1431  fMinLambda0LG=0.1;
1432  fMaxLambda0LG=4.0;
1433  fMaxRtrack=0.025;
1434  fMinCellEnergy=0.4;//0.1//0.4
1435  fReferenceFileName="";//Reference.root
1437  fBadReco=kFALSE;
1438  fFillHeavyHisto=kFALSE;
1439  fGeometryName="";//EMCAL_COMPLETE12SMV1_DCAL_8SM
1440  fPileupFromSPD=kFALSE;
1441  fMinTime=-20.;
1442  fMaxTime=20.;
1443  fMostEneCellOnly=kFALSE;
1444 
1445  fBadChannelMapSet=kFALSE;
1448 
1449  //histograms
1450  fRawTimeNbins = 400; // Raw time settings should be like that all the time
1451  fRawTimeMin = 400.; // importent in pass1
1452  fRawTimeMax = 800.;
1453  fPassTimeNbins = 1000; // in pass2 should be (400,400,800)
1454  fPassTimeMin = -250.;// in pass3 should be (1000,-250,250)
1455  fPassTimeMax = 250.;
1456  fEnergyNbins = 100; // default settings was 500
1457  fEnergyMin = 0.;
1458  fEnergyMax = 20.;
1459  fEnergyLGNbins = 200; // default settings
1460  fEnergyLGMin = 0.;
1461  fEnergyLGMax = 100.;
1462  fFineNbins = 90; //was 4500 for T0 time studies
1463  fFineTmin = -500;
1464  fFineTmax = 400;
1465 }
1466 
1467 //________________________________________________________________________
1472 void AliAnalysisTaskEMCALTimeCalib::ProduceCalibConsts(TString inputFile,TString outputFile,Bool_t isFinal, Bool_t oneHistoAllBCs, Bool_t isPAR)
1473 {
1474  TFile *file =new TFile(inputFile.Data());
1475  if(file==0x0) {
1476  printf("Input file does not exist!\n");
1477  return;
1478  }
1479 
1480  TList *list=(TList*)file->Get("chistolist");
1481  if(list==0x0)
1482  {
1483  printf("List chistolist does not exist in file, trying chistosingle!\n");
1484  list=(TList*)file->Get("chistosingle");
1485  if(list==0x0)
1486  {
1487  printf("List chistosingle does not exist either in file, returning!\n");
1488  return;
1489  }
1490  }
1491  Int_t numPARs = 0;
1492  Int_t counter = 0;
1493  if(isPAR){
1494  TIter next(list);
1495  TObject* obj;
1496  while((obj = next())){
1497  TString name(obj->GetName());
1498  if(name.BeginsWith("RawTimeBeforePAR")) counter++;
1499  }
1500  }
1501  numPARs = Int_t(counter/4) - 1;
1502  printf("number of PARs found to be %d!\n",numPARs);
1503 
1504  if(numPARs == -1) isPAR = kFALSE;
1505 
1506  if(!oneHistoAllBCs){
1507  //high gain
1508  TH1F *h1[4];
1509  TH1F *h2[4];
1510  TH1F *h3[4];
1511  TH1F *hAllTimeAvBC[4];
1512  TH1F *hAllTimeRMSBC[4];
1513 
1514  //low gain
1515  TH1F *h4[4];
1516  TH1F *h5[4];
1517  TH1F *h6[4];
1518  TH1F *hAllTimeAvLGBC[4];
1519  TH1F *hAllTimeRMSLGBC[4];
1520 
1521  //PAR histos
1522  TH1F *h1PAR[numPARs+1][4];
1523  TH1F *h2PAR[numPARs+1][4];
1524  //TH1F *h3PAR[numPARs+1][4];
1525  TH1F *hAllTimeAvBCPAR[numPARs+1][4];
1526  //TH1F *hAllTimeRMSBCPAR[numPARs+1][4];
1527 
1528  TH1F *h4PAR[numPARs+1][4];
1529  TH1F *h5PAR[numPARs+1][4];
1530  //TH1F *h6PAR[numPARs+1][4];
1531  TH1F *hAllTimeAvLGBCPAR[numPARs+1][4];
1532  //TH1F *hAllTimeRMSLGBCPAR[numPARs+1][4];
1533 
1534  TH2D* raw2D[4];
1535  TH2D* rawLG2D[4];
1536 
1537  if(isFinal==kFALSE){//first itereation
1538  for(Int_t i=0;i<4;i++){
1539  h1[i]=(TH1F *)list->FindObject(Form("RawTimeSumBC%d",i));
1540  h2[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesBC%d",i));
1541  h3[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqBC%d",i));
1542 
1543  h4[i]=(TH1F *)list->FindObject(Form("RawTimeSumLGBC%d",i));
1544  h5[i]=(TH1F *)list->FindObject(Form("RawTimeEntriesLGBC%d",i));
1545  h6[i]=(TH1F *)list->FindObject(Form("RawTimeSumSqLGBC%d",i));
1546 
1547  if(isPAR){ //set-up histograms for different PAR time regions
1548  for(Int_t iPAR = 0; iPAR <= numPARs; iPAR++){
1549  raw2D[i] = (TH2D*)list->FindObject(Form("RawTimeBeforePAR%dBC%d", iPAR+1, i));
1550  rawLG2D[i] = (TH2D*)list->FindObject(Form("RawTimeLGBeforePAR%dBC%d", iPAR+1, i));
1551  h1PAR[iPAR][i] = new TH1F(Form("hAllTimeSumPAR%dBC%d",iPAR, i), Form("hAlltimeSumPAR%dBC%d",iPAR, i), raw2D[i]->GetXaxis()->GetNbins(), raw2D[i]->GetXaxis()->GetXmin(), raw2D[i]->GetXaxis()->GetXmax());
1552  hAllTimeAvBCPAR[iPAR][i] = new TH1F(Form("hAllTimeAvPAR%dBC%d",iPAR, i), Form("hAlltimeAvPAR%dBC%d",iPAR, i), raw2D[i]->GetXaxis()->GetNbins(), raw2D[i]->GetXaxis()->GetXmin(), raw2D[i]->GetXaxis()->GetXmax());
1553  h2PAR[iPAR][i] = (TH1F*)raw2D[i]->ProjectionX(Form("hAllTimeEntriesPAR%dBC%d",iPAR, i), 0, raw2D[i]->GetYaxis()->GetNbins());
1554 
1555  h4PAR[iPAR][i] = new TH1F(Form("hAllTimeSumLGPAR%dBC%d",iPAR, i), Form("hAllTimeSumLGPAR%dBC%d",iPAR, i), raw2D[i]->GetXaxis()->GetNbins(), raw2D[i]->GetXaxis()->GetXmin(), raw2D[i]->GetXaxis()->GetXmax());
1556  hAllTimeAvLGBCPAR[iPAR][i] = new TH1F(Form("hAllTimeAvLGPAR%dBC%d",iPAR, i), Form("hAlltimeAvLGPAR%dBC%d",iPAR, i), raw2D[i]->GetXaxis()->GetNbins(), raw2D[i]->GetXaxis()->GetXmin(), raw2D[i]->GetXaxis()->GetXmax());
1557  h5PAR[iPAR][i] = (TH1F*)raw2D[i]->ProjectionX(Form("hAllTimeEntriesPAR%dLGBC%d",iPAR, i), 0, raw2D[i]->GetYaxis()->GetNbins());
1558  for(int ixbin = 0; ixbin < raw2D[i]->GetXaxis()->GetNbins(); ixbin++){
1559  float sumtime = 0.0;
1560  float sumLGtime = 0.0;
1561  for(int iybin = 0; iybin < raw2D[i]->GetYaxis()->GetNbins(); iybin++){
1562  sumtime += raw2D[i]->GetBinContent(ixbin, iybin)*raw2D[i]->GetYaxis()->GetBinCenter(iybin);
1563  sumLGtime += rawLG2D[i]->GetBinContent(ixbin, iybin)*rawLG2D[i]->GetYaxis()->GetBinCenter(iybin);
1564  }//end of loop over y-bins
1565  h1PAR[iPAR][i]->SetBinContent(ixbin, sumtime);
1566  h4PAR[iPAR][i]->SetBinContent(ixbin, sumLGtime);
1567  if(h2PAR[iPAR][i]->GetBinContent(ixbin) ==0){
1568  hAllTimeAvBCPAR[iPAR][i]->SetBinContent(ixbin, 0);
1569  }else{
1570  hAllTimeAvBCPAR[iPAR][i]->SetBinContent(ixbin, h1PAR[iPAR][i]->GetBinContent(ixbin)/h2PAR[iPAR][i]->GetBinContent(ixbin));
1571  }
1572 
1573  if(h5PAR[iPAR][i]->GetBinContent(ixbin) ==0){
1574  hAllTimeAvLGBCPAR[iPAR][i]->SetBinContent(ixbin, 0);
1575  }else{
1576  hAllTimeAvLGBCPAR[iPAR][i]->SetBinContent(ixbin, h4PAR[iPAR][i]->GetBinContent(ixbin)/h5PAR[iPAR][i]->GetBinContent(ixbin));
1577  }
1578  }//end of loop over x-bins
1579 
1580  }//end of loop over PARs
1581  }//end of if(isPAR)
1582  }//end of loop over BC
1583  } else {//final iteration
1584  for(Int_t i=0;i<4;i++){
1585  h1[i]=(TH1F *)list->FindObject(Form("hTimeSum%d",i));
1586  h2[i]=(TH1F *)list->FindObject(Form("hTimeEnt%d",i));
1587  h3[i]=(TH1F *)list->FindObject(Form("hTimeSumSq%d",i));
1588 
1589  h4[i]=(TH1F *)list->FindObject(Form("hTimeLGSum%d",i));
1590  h5[i]=(TH1F *)list->FindObject(Form("hTimeLGEnt%d",i));
1591  h6[i]=(TH1F *)list->FindObject(Form("hTimeLGSumSq%d",i));
1592  }
1593  }
1594  //AliWarning("Input histograms read.");
1595 
1596  for(Int_t i=0;i<4;i++){
1597  hAllTimeAvBC[i]=new TH1F(Form("hAllTimeAvBC%d",i),Form("hAllTimeAvBC%d",i),h1[i]->GetNbinsX(),h1[i]->GetXaxis()->GetXmin(),h1[i]->GetXaxis()->GetXmax());
1598  hAllTimeRMSBC[i]=new TH1F(Form("hAllTimeRMSBC%d",i),Form("hAllTimeRMSBC%d",i),h3[i]->GetNbinsX(),h3[i]->GetXaxis()->GetXmin(),h3[i]->GetXaxis()->GetXmax());
1599 
1600  hAllTimeAvLGBC[i]=new TH1F(Form("hAllTimeAvLGBC%d",i),Form("hAllTimeAvLGBC%d",i),h4[i]->GetNbinsX(),h4[i]->GetXaxis()->GetXmin(),h4[i]->GetXaxis()->GetXmax());
1601  hAllTimeRMSLGBC[i]=new TH1F(Form("hAllTimeRMSLGBC%d",i),Form("hAllTimeRMSLGBC%d",i),h6[i]->GetNbinsX(),h6[i]->GetXaxis()->GetXmin(),h6[i]->GetXaxis()->GetXmax());
1602  }
1603 
1604  //AliWarning("New histograms booked.");
1605 
1606  //important remark: we use 'underflow bin' for absid=0 in OADB . That's why there is j-1 below.
1607  for(Int_t i=0;i<4;i++){
1608  for(Int_t j=1;j<=h1[i]->GetNbinsX();j++){
1609  //high gain
1610  if(h2[i]->GetBinContent(j)!=0){
1611  hAllTimeAvBC[i]->SetBinContent(j-1,h1[i]->GetBinContent(j)/h2[i]->GetBinContent(j));
1612  hAllTimeRMSBC[i]->SetBinContent(j-1,TMath::Sqrt(h3[i]->GetBinContent(j)/h2[i]->GetBinContent(j)) );
1613  } else {
1614  hAllTimeAvBC[i]->SetBinContent(j-1,0.);
1615  hAllTimeRMSBC[i]->SetBinContent(j-1,0.);
1616  }
1617  //low gain
1618  if(h5[i]->GetBinContent(j)!=0){
1619  hAllTimeAvLGBC[i]->SetBinContent(j-1,h4[i]->GetBinContent(j)/h5[i]->GetBinContent(j));
1620  hAllTimeRMSLGBC[i]->SetBinContent(j-1,TMath::Sqrt(h6[i]->GetBinContent(j)/h5[i]->GetBinContent(j)) );
1621  } else {
1622  hAllTimeAvLGBC[i]->SetBinContent(j-1,0.);
1623  hAllTimeRMSLGBC[i]->SetBinContent(j-1,0.);
1624  }
1625 
1626  }
1627  }
1628 
1629  //AliWarning("Average and rms calculated.");
1630  TFile *fileNew=new TFile(outputFile.Data(),"recreate");
1631  for(Int_t i=0;i<4;i++){
1632  if(isPAR){
1633  for(Int_t iPAR = 0; iPAR <= numPARs; iPAR++){
1634  hAllTimeAvBCPAR[iPAR][i]->Write();
1635  //hAllTimeRMSBCPAR[iPAR][i]->Write();
1636  hAllTimeAvLGBCPAR[iPAR][i]->Write();
1637  //hAllTimeRMSLGBCPAR[iPAR][i]->Write();
1638  }
1639  }else{
1640  hAllTimeAvBC[i]->Write();
1641  hAllTimeRMSBC[i]->Write();
1642  hAllTimeAvLGBC[i]->Write();
1643  hAllTimeRMSLGBC[i]->Write();
1644  }
1645  }
1646 
1647  //AliWarning(Form("Histograms saved in %s file.",outputFile.Data()));
1648 
1649  fileNew->Close();
1650  delete fileNew;
1651 
1652  for(Int_t i=0;i<4;i++){
1653  delete hAllTimeAvBC[i];
1654  delete hAllTimeRMSBC[i];
1655  delete hAllTimeAvLGBC[i];
1656  delete hAllTimeRMSLGBC[i];
1657 
1658  if(isPAR){ //set-up histograms for different PAR time regions
1659  for(Int_t iPAR = 0; iPAR <= numPARs; iPAR++){
1660  delete h1PAR[iPAR][i];
1661  delete hAllTimeAvBCPAR[iPAR][i];
1662  delete h2PAR[iPAR][i];
1663  delete h4PAR[iPAR][i];
1664  delete hAllTimeAvLGBCPAR[iPAR][i];
1665  delete h5PAR[iPAR][i];
1666  }
1667  }
1668  }
1669  list->SetOwner(1);
1670  delete list;
1671  file->Close();
1672  delete file;
1673  }else{
1674 
1675  //high gain
1676  TH1F *h1;
1677  TH1F *h2;
1678  TH1F *h3;
1679  TH1S *hAllTimeAvBC;
1680  TH1S *hAllTimeRMSBC;
1681 
1682  //low gain
1683  TH1F *h4;
1684  TH1F *h5;
1685  TH1F *h6;
1686  TH1S *hAllTimeAvLGBC;
1687  TH1S *hAllTimeRMSLGBC;
1688 
1689  h1=(TH1F *)list->FindObject("hTimeSumAllBCs");
1690  h2=(TH1F *)list->FindObject("hTimeEntAllBCs");
1691  h3=(TH1F *)list->FindObject("hTimeSumSqAllBCs");
1692 
1693  h4=(TH1F *)list->FindObject("hTimeLGSumAllBCs");
1694  h5=(TH1F *)list->FindObject("hTimeLGEntAllBCs");
1695  h6=(TH1F *)list->FindObject("hTimeLGSumSqAllBCs");
1696  //AliWarning("Input histograms read.");
1697 
1698  hAllTimeAvBC=new TH1S("hAllTimeAv","hAllTimeAv",h1->GetNbinsX(),h1->GetXaxis()->GetXmin(),h1->GetXaxis()->GetXmax());
1699  hAllTimeRMSBC=new TH1S("hAllTimeRMS","hAllTimeRMS",h3->GetNbinsX(),h3->GetXaxis()->GetXmin(),h3->GetXaxis()->GetXmax());
1700 
1701  hAllTimeAvLGBC=new TH1S("hAllTimeAvLG","hAllTimeAvLG",h4->GetNbinsX(),h4->GetXaxis()->GetXmin(),h4->GetXaxis()->GetXmax());
1702  hAllTimeRMSLGBC=new TH1S("hAllTimeRMSLG","hAllTimeRMSLG",h6->GetNbinsX(),h6->GetXaxis()->GetXmin(),h6->GetXaxis()->GetXmax());
1703 
1704  //AliWarning("New histograms booked.");
1705 
1706  //important remark: we use 'underflow bin' for absid=0 in OADB . That's why there is j-1 below.
1707  for(Int_t j=1;j<=h1->GetNbinsX();j++){
1708  //high gain
1709  if(h2->GetBinContent(j)!=0){
1710  hAllTimeAvBC->SetBinContent(j-1,h1->GetBinContent(j)/h2->GetBinContent(j));
1711  hAllTimeRMSBC->SetBinContent(j-1,TMath::Sqrt(h3->GetBinContent(j)/h2->GetBinContent(j)) );
1712  } else {
1713  hAllTimeAvBC->SetBinContent(j-1,0.);
1714  hAllTimeRMSBC->SetBinContent(j-1,0.);
1715  }
1716  //low gain
1717  if(h5->GetBinContent(j)!=0){
1718  hAllTimeAvLGBC->SetBinContent(j-1,h4->GetBinContent(j)/h5->GetBinContent(j));
1719  hAllTimeRMSLGBC->SetBinContent(j-1,TMath::Sqrt(h6->GetBinContent(j)/h5->GetBinContent(j)) );
1720  } else {
1721  hAllTimeAvLGBC->SetBinContent(j-1,0.);
1722  hAllTimeRMSLGBC->SetBinContent(j-1,0.);
1723  }
1724 
1725  }
1726 
1727  //AliWarning("Average and rms calculated.");
1728  TFile *fileNew=new TFile(outputFile.Data(),"recreate");
1729 
1730  hAllTimeAvBC->Write();
1731  hAllTimeRMSBC->Write();
1732  hAllTimeAvLGBC->Write();
1733  hAllTimeRMSLGBC->Write();
1734 
1735  //AliWarning(Form("Histograms saved in %s file.",outputFile.Data()));
1736 
1737  fileNew->Close();
1738  delete fileNew;
1739 
1740  delete hAllTimeAvBC;
1741  delete hAllTimeRMSBC;
1742  delete hAllTimeAvLGBC;
1743  delete hAllTimeRMSLGBC;
1744 
1745  list->SetOwner(1);
1746  delete list;
1747  file->Close();
1748  delete file;
1749 
1750  }
1751 
1752  //AliWarning("Pointers deleted. Memory cleaned.");
1753 }
1754 
1755 //________________________________________________________________________
1759 void AliAnalysisTaskEMCALTimeCalib::ProduceOffsetForSMsV2(Int_t runNumber,TString inputFile,TString outputFile, Bool_t offset100, Bool_t justL1phase, TString PARFilename){
1760 
1761  const Double_t lowerLimit[]={
1762  0,
1763  1152,
1764  2304,
1765  3456,
1766  4608,
1767  5760,
1768  6912,
1769  8064,
1770  9216,
1771  10368,
1772  11520,
1773  11904,
1774  12288,
1775  13056,
1776  13824,
1777  14592,
1778  15360,
1779  16128,
1780  16896,
1781  17280};
1782 
1783  const Double_t upperLimit[]={
1784  1151 ,
1785  2303 ,
1786  3455 ,
1787  4607 ,
1788  5759 ,
1789  6911 ,
1790  8063 ,
1791  9215 ,
1792  10367,
1793  11519,
1794  11903,
1795  12287,
1796  13055,
1797  13823,
1798  14591,
1799  15359,
1800  16127,
1801  16895,
1802  17279,
1803  17663};
1804 
1805  PARInfo info;
1806  info.numPARs = 0;
1807  Bool_t isPAR = kFALSE;
1808  if(PARFilename.Length() != 0){
1809  std::ifstream input;
1810  int inputrunnumber = 0, numPARs = 0;
1811  ULong64_t PAR = 0;
1812  input.open(PARFilename.Data());
1813  if(!input.good()){
1814  printf("PAR info file not accessable: %s\n", PARFilename.Data());
1815  return;
1816  }
1817  while(input.good()){
1818  input >> inputrunnumber >> numPARs;
1819  if(!input.good()) break;
1820  info.runNumber = inputrunnumber;
1821  info.numPARs = numPARs;
1822  //printf("\n\n!!!!\n\n from file: runnumber = %d, numPars = %d\n\n", info.runNumber, info.numPARs);
1823  if(numPARs <= 0 || numPARs > 10){
1824  printf("Number of PARS incorrectly found to be %d!\n", numPARs);
1825  return;
1826  }
1827  for(int iPAR = 0; iPAR < numPARs; iPAR++){
1828  input >> PAR;
1829  if(info.runNumber == runNumber){
1830  info.PARGlobalBCs.push_back(PAR);
1831  }
1832  }
1833  if(info.runNumber == runNumber) break;
1834  }
1835  input.close();
1836 
1837  if(info.runNumber != runNumber){
1838  isPAR = kFALSE;
1839  info.numPARs = 0;
1840  }else{
1841  isPAR = kTRUE;
1842  printf("---- NEW RUN NUMBER ----\n");
1843  printf("info.runNumber = %d\n", info.runNumber);
1844  printf("info.numPARs = %d\n", info.numPARs);
1845  for(int i = 0; i < info.numPARs; i++){
1846  printf("info.PARGlobalBCs[%d] = %llu\n", i, info.PARGlobalBCs[i]);
1847  }
1848  }
1849  }
1850 
1851  TFile *file =new TFile(inputFile.Data());
1852  if(file==0x0) return;
1853 
1854  TH1F *ccBC[4];
1855  Bool_t shouldBeEmpty[4];
1856  TH1F *ccBCPAR[info.numPARs+1][4];
1857  Int_t emptyCounter;
1858  Bool_t shouldBeEmptyPAR[info.numPARs+1][4];
1859 
1860  for(Int_t i = 0; i < kNBCmask; i++){
1861  if(isPAR){
1862  for(Int_t iPAR = 0; iPAR <= info.numPARs; iPAR++){
1863  ccBCPAR[iPAR][i] = (TH1F*)file->Get(Form("hAllTimeAvPAR%dBC%d", iPAR, i));
1864  shouldBeEmptyPAR[iPAR][i]=kFALSE;
1865  emptyCounter=0;
1866  for(Int_t j=0;j<upperLimit[19];j++){
1867  if(ccBCPAR[iPAR][i]->GetBinContent(j)>0.) emptyCounter++;
1868  }
1869  if(emptyCounter<400) shouldBeEmptyPAR[iPAR][i]=kTRUE;
1870  printf("Non-zero channels %d BC %d PAR %d should be empty: %d \n",emptyCounter,i,iPAR,shouldBeEmptyPAR[iPAR][i]);
1871  }
1872  //it cannot be empty after par when befor was not empty
1873  //need to correct for this for events after PAR(s)
1874  for(Int_t iPAR = 1; iPAR <= info.numPARs; iPAR++){
1875  if(shouldBeEmptyPAR[iPAR][i] && !shouldBeEmptyPAR[0][i] ) {
1876  shouldBeEmptyPAR[iPAR][i] = kFALSE;
1877  printf("BC %d PAR %d can NOT be empty because before any PAR was filled. Correct to %d \n",i,iPAR,shouldBeEmptyPAR[iPAR][i]);
1878  }
1879  }
1880  }else{
1881  ccBC[i]=(TH1F*) file->Get(Form("hAllTimeAvBC%d",i));
1882  shouldBeEmpty[i]=kFALSE;
1883  emptyCounter=0;
1884  for(Int_t j=0;j<upperLimit[19];j++){
1885  if(ccBC[i]->GetBinContent(j)>0.) emptyCounter++;
1886  }
1887  if(emptyCounter<1500) shouldBeEmpty[i]=kTRUE;
1888  printf("Non-zero channels %d BC %d should be empty: %d \n",emptyCounter,i,shouldBeEmpty[i]);
1889  }
1890  }
1891 
1892  TH1C *hRun=new TH1C(Form("h%d",runNumber),Form("h%d",runNumber),19,0,19);
1893  TH1C *hPARRun[info.numPARs+1];
1894  Int_t fitResult=0;
1895  Double_t minimumValue=10000.;
1896  Int_t minimumIndex=-1;
1897  Double_t meanBC[4];
1898 
1899  Double_t fitParameter=0;
1900  TF1 *f1=new TF1("f1","pol0",0,17664);
1901  Bool_t orderTest=kTRUE;
1902  Int_t iorder=0;//order index
1903  Int_t j=0;//BC index
1904  Int_t L1shift=0;
1905  Int_t totalValue=0;
1906 
1907  for(Int_t iPAR = 0; iPAR <= info.numPARs; iPAR++){
1908  if(iPAR ==0){//iPAR=0 means before any PAR
1909  hPARRun[iPAR] =new TH1C(Form("h%d", runNumber), Form("h%d", runNumber),19,0,19);
1910  }else{
1911  hPARRun[iPAR] =new TH1C(Form("h%d_%llu", runNumber, (ULong64_t)info.PARGlobalBCs[iPAR-1]), Form("h%d_%llu", runNumber, (ULong64_t)info.PARGlobalBCs[iPAR-1]),19,0,19);
1912  }
1913  for(Int_t i=0;i<20;i++){
1914  minimumValue=10000;
1915  for(j=0;j<kNBCmask;j++){
1916  if(isPAR){
1917  if(shouldBeEmptyPAR[iPAR][j]){
1918  meanBC[j]=-1;
1919  continue;
1920  }
1921  }else{
1922  if(shouldBeEmpty[j]) {
1923  meanBC[j]=-1;
1924  continue;
1925  }
1926  }
1927  if(isPAR){
1928  fitResult=ccBCPAR[iPAR][j]->Fit("f1", "CQN", "", lowerLimit[i],upperLimit[i]);
1929  }else{
1930  fitResult=ccBC[j]->Fit("f1","CQN","",lowerLimit[i],upperLimit[i]);
1931  }
1932  if(fitResult<0){
1933  //hRun->SetBinContent(i,0);//correct it please
1934  meanBC[j]=-1;
1935  if(isPAR){
1936  printf("Fit failed for SM %d BC%d PAR %d, integral %f\n",i,j,iPAR,ccBCPAR[iPAR][j]->Integral(lowerLimit[i],upperLimit[i]));
1937  }else{
1938  printf("Fit failed for SM %d BC%d, integral %f\n",i,j,ccBC[j]->Integral(lowerLimit[i],upperLimit[i]));
1939  }
1940  continue;
1941  } else {
1942  fitParameter = f1->GetParameter(0);
1943  }
1944  if(offset100 && (j==0 || j==1)) {
1945  //the 100 ns offset was removed in LHC15n muon_calo_pass1 and further reconstructions
1946  fitParameter+=100;
1947  }
1948  meanBC[j]=fitParameter;
1949 
1950  if(fitParameter>0 && fitParameter<minimumValue){
1951  minimumValue = fitParameter;
1952  minimumIndex = j;
1953  }
1954  }//end of loop over BCs
1955 
1956  if( minimumValue/25-(Int_t)(minimumValue/25)>0.5 ) {
1957  L1shift=(Int_t)(minimumValue/25.)+1;
1958  } else {
1959  L1shift=(Int_t)(minimumValue/25.);
1960  }
1961 
1962  if(TMath::Abs(minimumValue/25-(Int_t)(minimumValue/25)-0.5)<0.05)
1963  printf("Run %d, PAR %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,iPAR,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);
1964 
1965  if(justL1phase) totalValue = minimumIndex;
1966  else totalValue = L1shift<<2 | minimumIndex ;
1967  //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);
1968 
1969  if(isPAR){
1970  hPARRun[iPAR]->SetBinContent(i,totalValue);
1971  }else{
1972  hRun->SetBinContent(i,totalValue);
1973  }
1974  orderTest=kTRUE;
1975  for(iorder=minimumIndex;iorder<minimumIndex+4-1;iorder++){
1976  if( meanBC[(iorder+1)%4] <= meanBC[iorder%4] ) orderTest=kFALSE;
1977  }
1978 
1979  if(!orderTest){
1980  if(isPAR)
1981  printf("run %d, PAR %d, SM %d, min index %d meanBC %f %f %f %f, order ok? %d\n",runNumber,iPAR,i,minimumIndex,meanBC[0],meanBC[1],meanBC[2],meanBC[3],orderTest);
1982  else
1983  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);
1984  }
1985 
1986  //patch for runs with not filled one, two or three BCs
1987  //manual patch for LHC16q - pPb@5TeV - only BC0 is filled and phase rotate
1988  if(isPAR){// PAR case
1989  if(shouldBeEmptyPAR[iPAR][0] || shouldBeEmptyPAR[iPAR][1] || shouldBeEmptyPAR[iPAR][3] || shouldBeEmptyPAR[iPAR][3]){
1990  Double_t newMean = meanBC[minimumIndex]-600;
1991  if(newMean<=12.5){
1992  hPARRun[iPAR]->SetBinContent(i,minimumIndex);
1993  } else {
1994  Int_t minIndexTmp=-1;
1995  if(newMean/25. - (Int_t)(newMean/25.) <0.5)
1996  minIndexTmp = (Int_t)(newMean/25.);
1997  else
1998  minIndexTmp = 1+(Int_t)(newMean/25.);
1999 
2000  hPARRun[iPAR]->SetBinContent(i,(4-minIndexTmp+minimumIndex)%4);
2001  }
2002  printf("run with missing BC; PAR %d; new L1 phase in SM%d set to %d\n",iPAR,i,(Int_t)hPARRun[iPAR]->GetBinContent(i));
2003  }
2004  } else {//regular case
2005  if(shouldBeEmpty[0] || shouldBeEmpty[1] || shouldBeEmpty[2] || shouldBeEmpty[3]){
2006  Double_t newMean = meanBC[minimumIndex]-600;
2007  if(newMean<=12.5){
2008  hRun->SetBinContent(i,minimumIndex);
2009  } else {
2010  Int_t minIndexTmp=-1;
2011  if(newMean/25. - (Int_t)(newMean/25.) <0.5)
2012  minIndexTmp = (Int_t)(newMean/25.);
2013  else
2014  minIndexTmp = 1+(Int_t)(newMean/25.);
2015 
2016  hRun->SetBinContent(i,(4-minIndexTmp+minimumIndex)%4);
2017  //cout<<newMean/25.<<" int "<<(Int_t)(newMean/25.)<<" dif "<< newMean/25.-(Int_t)(newMean/25.)<<endl;
2018  }
2019  printf("run with missing BC; new L1 phase in SM%d set to %d\n",i,(Int_t)hRun->GetBinContent(i));
2020  }
2021  }//end of patch for LHC16q and other runs with not filled BCs
2022 
2023  }//end of loop over SM
2024  }//end of loop over PARs
2025 
2026  delete f1;
2027  TFile *fileNew=new TFile(outputFile.Data(),"update");
2028  if(isPAR){
2029  for(Int_t iPAR = 0; iPAR <= info.numPARs; iPAR++){
2030  hPARRun[iPAR]->Write();
2031  delete hPARRun[iPAR];
2032  }
2033  // create tree for PAR global IDs
2034  ULong64_t ParGlobalBCs;
2035  TTree *treePAR=new TTree(Form("t%d_GID",runNumber),"Tree with Global ID");
2036  treePAR->Branch("GID",&ParGlobalBCs,"GID/l");
2037 
2038  for(Int_t iPAR = 0; iPAR < info.numPARs; iPAR++){
2039  ParGlobalBCs=(ULong64_t)(info.PARGlobalBCs[iPAR]);
2040  //printf("infoPAR %llu in tree %llu",info.PARGlobalBCs[iPAR],ParGlobalBCs);
2041  treePAR->Fill();
2042  }
2043  treePAR->Write();
2044  }else{
2045  hRun->Write();
2046  delete hRun;
2047  }
2048  fileNew->Close();
2049  delete fileNew;
2050 
2051  file->Close();
2052  delete file;
2053 }
2054 
2055 //____________________________________________________
2057 {
2058  if(fBadChannelMapSet) return;
2059  AliOADBContainer *contBC=new AliOADBContainer("");
2060  contBC->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALBadChannels.root").data(),"AliEMCALBadChannels");
2061  printf("contBC %p, ent %d\n",contBC,contBC->GetNumberOfEntries());
2062  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRunNumber);
2063  if(arrayBC) {
2064  AliInfo("Remove EMCAL bad cells");
2066  for (Int_t i=0; i<kNSM; ++i) {
2067  TH2I *hbm = (TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
2068  if (!hbm) {
2069  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
2070  continue;
2071  }
2072  hbm->SetDirectory(0);
2073  fBadChannelMapArray->AddAt(hbm,i);
2074 
2075  } // loop over SMs
2076  } else AliInfo("Do NOT remove EMCAL bad channels\n"); // run array
2077 
2078  delete contBC;
2079  fBadChannelMapSet = kFALSE;//BC map is not fixed at the beginning but can change r-by-r
2080 } // Bad channel map loaded
2081 
2082 //____________________________________________________
2084 {
2085  if(fBadChannelMapSet) return;
2086 
2087  TFile *referenceFile = TFile::Open(fBadChannelFileName.Data());
2088  if(referenceFile==0x0) {
2089  AliFatal("*** NO bad channel map FILE");
2090  }
2091 
2092  TH1F *hbm = (TH1F*)referenceFile->Get("h1");
2093  if (!hbm) {
2094  AliError("Can not get EMCALBadChannelMap");
2095  }
2096  fBadChannelMapArray = new TObjArray(1);
2097  fBadChannelMapArray->AddAt(hbm,0);
2098  fBadChannelMapSet=kTRUE;//BC map is fixed at the beginning for whole dataset
2099 } // Bad channel map loaded
2100 
2101 
2102 //_____________________________________________________________________
2107 }
2108 
2109 //_____________________________________________________________________
2110 // Load PAR info from text file
2112  std::ifstream input;
2113  int runnumber = 0, numPARs = 0, numRuns=0;
2114  ULong64_t PAR = 0;
2115  gSystem->ExpandPathName(PARFileName);
2116  //handle case of PAR file in Alien location, needs to be copied to working directory before ifstream can open.
2117  if(PARFileName.Contains("alien://")){
2118  TString localFileName(gSystem->BaseName(PARFileName.Data()));
2119  TFile::Cp(PARFileName.Data(), localFileName.Data());
2120  PARFileName = localFileName;
2121  }
2122  input.open(PARFileName.Data());
2123  if(!input.good()){
2124  AliFatal(Form("PAR info file not accessable: %s", PARFileName.Data()));
2125  }
2126  while(input.good()){
2127  input >> runnumber >> numPARs;
2128  if(!input.good()) break;
2129  PARInfo info;
2130  info.runNumber = runnumber;
2131  info.numPARs = numPARs;
2132  //printf("\n\n!!!!\n\n from file: runnumber = %d, numPars = %d\n\n", info.runNumber, info.numPARs);
2133  if(numPARs <= 0 || numPARs > 10){
2134  AliFatal(Form("Number of PARS incorrectly found to be %d!", numPARs));
2135  }
2136  for(int iPAR = 0; iPAR < numPARs; iPAR++){
2137  input >> PAR;
2138  info.PARGlobalBCs.push_back(PAR);
2139  }
2140  fPARvec.push_back(info);
2141  numRuns++;
2142  }
2143  printf("number of runs processed in PAR file: %d\n", numRuns);
2144  input.close();
2145 }
2146 
2147 //_______________________________________________________________________
2148 // Get Par info for the current run number, set-up PAR info variables
2150  if(runnum < 200000) AliFatal(Form("Bad Run Number %d passed to GetPARInfo!", runnum));
2151  if(fRunNumber!=runnum) fRunNumber = runnum;
2152  fIsPARRun = kFALSE;
2154  for(unsigned int iPARrun = 0; iPARrun < fPARvec.size(); iPARrun++){
2155  //printf("from stored vectors: %d %d", fPARvec[iPARrun].runNumber, fPARvec[iPARrun].numPARs);
2156  //for(int i = 0; i < fPARvec[iPARrun].numPARs; i++){
2157  // printf(" %llu", fPARvec[iPARrun].PARGlobalBCs[i]);
2158  //}
2159  //printf("\n");
2160  if (fRunNumber==fPARvec[iPARrun].runNumber){
2161  //set PAR flag & setup copy of specific PAR info
2162  fIsPARRun = kTRUE;
2164  fCurrentPARs.numPARs = fPARvec[iPARrun].numPARs;
2165  for(int ipar = 0; ipar < fPARvec[iPARrun].numPARs; ipar++){
2166  fCurrentPARs.PARGlobalBCs.push_back(fPARvec[iPARrun].PARGlobalBCs[ipar]);
2167  }
2168  }
2169  }
2170 }
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 fCurrentPARIndex
Par Info for current Run Number.
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
std::vector< std::vector< TH2F * > > fhRawTimeLGPARs
!
TSystem * gSystem
Double_t fMinTime
minimum cluster time after correction
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
Bool_t fMostEneCellOnly
flag to use calibration on most energetic cell in cluster only
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
std::vector< std::vector< TH2F * > > fhRawTimePARs
!
Double_t fEnergyMax
upper range of histo with energy HG
Bool_t fOneHistAllBCs
flag to use one histogram for all the BCs instead of four
static void ProduceCalibConsts(TString inputFile="time186319testWOL0.root", TString outputFile="Reference.root", Bool_t isFinal=kFALSE, Bool_t oneHistoAllBCs=kFALSE, Bool_t isPAR=kFALSE)
Double_t fMinClusterEnergy
minimum cluster energy
float Float_t
Definition: External.C:68
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
std::vector< PARInfo > fPARvec
vector of PAR info for all runs
Definition: External.C:228
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 ProduceOffsetForSMsV2(Int_t runNumber, TString inputFile="Reference.root", TString outputFile="ReferenceSM.root", Bool_t offset100=kTRUE, Bool_t justL1phase=kTRUE, TString PARfilename="")
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
Bool_t fIsPARRun
Which PAR the currnt event is after.
TH2F * fhTimeVsIdBC[kNBCmask]
! 4 BCmask HG
TFile * file
TList with histograms for a given trigger.
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
void GetPARInfoForRunNumber(Int_t runnum)
Does current run have PAR info?
TH2F * fhTimeDsupBC[kNSM][kNBCmask]
! 20 x 4 high gain
TH2F * fhTimeDsupLGBC[kNSM][kNBCmask]
! 20 x 4 low gain