fix indentations, open trigger cluster at BC=0 time window definition
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliCaloTrackReader.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and 
18 // Central Barrel Tracking detectors (CTS).
19 // Not all MC particles/tracks/clusters are kept, some kinematical/fiducial restrictions are done.
20 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
21 //                 : AliCaloTrackMCReader : Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
22 //                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS) 
23 //-- Author: Gustavo Conesa (LNF-INFN) 
24 //////////////////////////////////////////////////////////////////////////////
25
26
27 // --- ROOT system ---
28 #include <TFile.h>
29 #include <TGeoManager.h>
30
31 // ---- ANALYSIS system ----
32 #include "AliMCEvent.h"
33 #include "AliAODMCHeader.h"
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliESDEvent.h"
36 #include "AliAODEvent.h"
37 #include "AliVTrack.h"
38 #include "AliVParticle.h"
39 #include "AliMixedEvent.h"
40 #include "AliESDtrack.h"
41 #include "AliESDtrackCuts.h"
42 #include "AliTriggerAnalysis.h"
43 #include "AliESDVZERO.h"
44 #include "AliVCaloCells.h"
45 #include "AliAnalysisManager.h"
46 #include "AliInputEventHandler.h"
47
48 // ---- Detectors ----
49 #include "AliPHOSGeoUtils.h"
50 #include "AliEMCALGeometry.h"
51 #include "AliEMCALRecoUtils.h"
52
53 // ---- CaloTrackCorr ---
54 #include "AliCalorimeterUtils.h"
55 #include "AliCaloTrackReader.h"
56
57 ClassImp(AliCaloTrackReader)
58   
59   
60 //________________________________________
61 AliCaloTrackReader::AliCaloTrackReader() : 
62 TObject(),                   fEventNumber(-1), //fCurrentFileName(""),
63 fDataType(0),                fDebug(0), 
64 fFiducialCut(0x0),           fCheckFidCut(kFALSE), 
65 fComparePtHardAndJetPt(0),   fPtHardAndJetPtFactor(0),
66 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
67 fCTSPtMin(0),                fEMCALPtMin(0),                  fPHOSPtMin(0), 
68 fCTSPtMax(0),                fEMCALPtMax(0),                  fPHOSPtMax(0),
69 fUseEMCALTimeCut(1),         fUseParamTimeCut(0),             fUseTrackTimeCut(0),
70 fEMCALTimeCutMin(-10000),    fEMCALTimeCutMax(10000),
71 fEMCALParamTimeCutMin(),     fEMCALParamTimeCutMax(),
72 fTrackTimeCutMin(-10000),    fTrackTimeCutMax(10000),
73 fUseTrackDCACut(0),
74 fAODBranchList(0x0),
75 fCTSTracks(0x0),             fEMCALClusters(0x0),             fPHOSClusters(0x0),
76 fEMCALCells(0x0),            fPHOSCells(0x0),
77 fInputEvent(0x0),            fOutputEvent(0x0),fMC(0x0),
78 fFillCTS(0),                 fFillEMCAL(0),                   fFillPHOS(0),
79 fFillEMCALCells(0),          fFillPHOSCells(0), 
80 fRecalculateClusters(kFALSE),fSelectEmbeddedClusters(kFALSE),
81 fTrackStatus(0),             fTrackFilterMask(0),             
82 fESDtrackCuts(0),            fESDtrackComplementaryCuts(0),   fConstrainTrack(kFALSE),
83 fSelectHybridTracks(0),      fSelectSPDHitTracks(kFALSE),
84 fTrackMult(0),               fTrackMultEtaCut(0.9),
85 fReadStack(kFALSE),          fReadAODMCParticles(kFALSE), 
86 fDeltaAODFileName(""),       fFiredTriggerClassName(""),      
87 fEventTriggerMask(0),        fMixEventTriggerMask(0),         fEventTriggerAtSE(0), 
88 fAnaLED(kFALSE),
89 fTaskName(""),               fCaloUtils(0x0), 
90 fMixedEvent(NULL),           fNMixedEvent(0),                 fVertex(NULL), 
91 fListMixedTracksEvents(),    fListMixedCaloEvents(),
92 fLastMixedTracksEvent(-1),   fLastMixedCaloEvent(-1),
93 fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
94 fEMCALClustersListName(""),  fZvtxCut(0.),                    
95 fAcceptFastCluster(kFALSE),  fRemoveLEDEvents(kTRUE),
96 //Trigger rejection
97 fRemoveBadTriggerEvents(0),  fTriggerPatchClusterMatch(0),
98 fTriggerPatchTimeWindow(),   fTriggerEventThreshold(0),
99 fTriggerClusterBC(0),        fTriggerClusterIndex(0),         fTriggerClusterId(0),
100 fIsExoticEvent(0),           fIsBadCellEvent(0),              fIsBadMaxCellEvent(0),
101 fIsTriggerMatch(0),
102
103 fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),
104 fDoVertexBCEventSelection(kFALSE),
105 fDoRejectNoTrackEvents(kFALSE),
106 fUseEventsWithPrimaryVertex(kFALSE),
107 fTriggerAnalysis (0x0),      fTimeStampEventSelect(0),
108 fTimeStampEventFracMin(0),   fTimeStampEventFracMax(0),
109 fTimeStampRunMin(0),         fTimeStampRunMax(0),
110 fNPileUpClusters(-1),        fNNonPileUpClusters(-1),         fNPileUpClustersCut(3),
111 fVertexBC(-200),             fRecalculateVertexBC(0),
112 fCentralityClass(""),        fCentralityOpt(0),
113 fEventPlaneMethod(""),       fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
114 {
115   //Ctor
116
117   //Initialize parameters
118   InitParameters();
119 }
120
121 //_______________________________________
122 AliCaloTrackReader::~AliCaloTrackReader() 
123 {
124   //Dtor
125   
126   delete fFiducialCut ;
127         
128   if(fAODBranchList)
129   {
130     fAODBranchList->Delete();
131     delete fAODBranchList ;
132   }  
133   
134   if(fCTSTracks)
135   {
136     if(fDataType!=kMC)fCTSTracks->Clear() ; 
137     else              fCTSTracks->Delete() ; 
138     delete fCTSTracks ;
139   }
140   
141   if(fEMCALClusters)
142   {
143     if(fDataType!=kMC)fEMCALClusters->Clear("C") ; 
144     else              fEMCALClusters->Delete() ; 
145     delete fEMCALClusters ;
146   }
147   
148   if(fPHOSClusters)
149   {
150     if(fDataType!=kMC)fPHOSClusters->Clear("C") ; 
151     else              fPHOSClusters->Delete() ; 
152     delete fPHOSClusters ;
153   }
154   
155   if(fVertex)
156   {
157     for (Int_t i = 0; i < fNMixedEvent; i++) 
158     {
159       delete [] fVertex[i] ;
160       
161     }
162     delete [] fVertex ;
163         }
164   
165   delete fESDtrackCuts;
166   delete fESDtrackComplementaryCuts;
167   delete fTriggerAnalysis;
168   
169   //  Pointers not owned, done by the analysis frame
170   //  if(fInputEvent)  delete fInputEvent ;
171   //  if(fOutputEvent) delete fOutputEvent ;
172   //  if(fMC)          delete fMC ;  
173   //  Pointer not owned, deleted by maker
174   //  if (fCaloUtils) delete fCaloUtils ;
175   
176 }
177
178 //________________________________________________________________________
179 Bool_t  AliCaloTrackReader::AcceptDCA(const Float_t pt, const Float_t dca)
180 {
181   // Accept track if DCA is smaller than function
182   
183   Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
184   
185   if(TMath::Abs(dca) < cut)
186     return kTRUE;
187   else
188     return kFALSE;
189   
190 }
191
192 //________________________________________________
193 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt()
194 {
195   // Check the event, if the requested ptHard is much smaller than the jet pT, then there is a problem.
196   // Only for PYTHIA.
197   
198   //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
199   
200   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
201   {
202     TParticle * jet =  0;
203     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
204     Int_t nTriggerJets =  pygeh->NTriggerJets();
205     Float_t ptHard = pygeh->GetPtHard();
206     
207     if(fDebug > 1) 
208       printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
209     
210     Float_t tmpjet[]={0,0,0,0};
211     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
212     {
213       pygeh->TriggerJet(ijet, tmpjet);
214       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
215       
216       if(fDebug > 1) 
217         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - jet %d; pycell jet pT %f\n",ijet, jet->Pt());
218       
219       //Compare jet pT and pt Hard
220       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
221       {
222         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
223                 ptHard, jet->Pt(), fPtHardAndJetPtFactor);
224         return kFALSE;
225       }
226     }
227     
228     if(jet) delete jet; 
229   }
230   
231   return kTRUE ;
232   
233 }
234
235 //____________________________________________________________________
236 Bool_t AliCaloTrackReader::ComparePtHardAndClusterPt()
237 {
238   // Check the event, if the requested ptHard is smaller than the calorimeter cluster E, then there is a problem.
239   // Only for PYTHIA.
240   
241   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
242   {
243     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
244     Float_t ptHard = pygeh->GetPtHard();
245     
246     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
247     for (Int_t iclus =  0; iclus <  nclusters; iclus++) 
248     {
249       AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ; 
250       Float_t ecluster = clus->E();
251       
252       if(ecluster > fPtHardAndClusterPtFactor*ptHard) 
253       {
254         printf("AliCaloTrackReader::ComparePtHardAndClusterPt() - Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f\n",ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard);
255
256         return kFALSE;
257       }
258     }
259     
260   }
261   
262   return kTRUE ;
263   
264 }
265
266 //____________________________________________
267 AliStack* AliCaloTrackReader::GetStack() const 
268 {
269   //Return pointer to stack
270   if(fMC)
271     return fMC->Stack();
272   else
273   {
274     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n"); 
275     return 0x0 ;
276   }
277 }
278
279 //__________________________________________________
280 TString AliCaloTrackReader::GetFiredTriggerClasses() 
281
282   // List of triggered classes in a TString
283
284   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetInputEvent());
285   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetInputEvent());
286     
287   if     (esdevent) return esdevent->GetFiredTriggerClasses();
288   else if(aodevent) return aodevent->GetFiredTriggerClasses();
289   else              return ""; // Mixed Event, MC event, does not have this trigger info
290
291 }
292
293 //______________________________________________
294 AliHeader* AliCaloTrackReader::GetHeader() const 
295 {
296   //Return pointer to header
297   if(fMC)
298   {
299     return fMC->Header();
300   }
301   else
302   {
303     printf("AliCaloTrackReader::Header is not available\n"); 
304     return 0x0 ;
305   }
306 }
307
308 //______________________________________________________________
309 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const 
310 {
311   //Return pointer to Generated event header
312   if     (ReadStack() && fMC)
313   {
314     return fMC->GenEventHeader();
315   }
316   else if(ReadAODMCParticles() && GetAODMCHeader())
317   {
318     //printf("AliCaloTrackReader::GetGenEventHeader() - N headers %d\n",GetAODMCHeader()->GetNCocktailHeaders());
319     if( GetAODMCHeader()->GetNCocktailHeaders() > 0)
320       return GetAODMCHeader()->GetCocktailHeader(0) ;
321     else 
322       return 0x0;
323   }
324   else 
325   {
326     //printf("AliCaloTrackReader::GetGenEventHeader() - MC header not available! \n");
327     return 0;
328   }
329 }
330
331 //____________________________________________________________________
332 TClonesArray* AliCaloTrackReader::GetAODMCParticles() const 
333 {
334   //Return list of particles in AOD. Do it for the corresponding input event.
335   
336   TClonesArray * rv = NULL ; 
337   if(fDataType == kAOD)
338   {
339     //Normal input AOD
340     AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
341     if(evt)
342       rv = (TClonesArray*)evt->FindListObject("mcparticles");
343     else  
344       printf("AliCaloTrackReader::GetAODMCParticles() - Null AOD event \n"); 
345   } 
346   else 
347   {
348     printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n"); 
349   }
350   
351   return rv ; 
352 }
353
354 //________________________________________________________
355 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const 
356 {
357   //Return MC header in AOD. Do it for the corresponding input event.
358   
359   AliAODMCHeader *mch = NULL;
360   
361   if(fDataType == kAOD)
362   {
363     AliAODEvent * aod = dynamic_cast<AliAODEvent*> (fInputEvent);
364     if(aod) mch = dynamic_cast<AliAODMCHeader*>(aod->FindListObject("mcHeader"));
365   }
366   else 
367   {
368     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
369   }
370   
371   return mch;
372 }
373
374 //___________________________________________________________
375 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
376 {
377   // Get the vertex BC
378     
379   Int_t vertexBC=vtx->GetBC();
380   if(!fRecalculateVertexBC) return vertexBC;
381   
382   // In old AODs BC not stored, recalculate it
383   // loop over the global track and select those which have small DCA to primary vertex (e.g. primary).
384   // If at least one of these primaries has valid BC != 0, then this vertex is a pile-up candidate.
385   // Execute after CTS
386   Double_t bz  = fInputEvent->GetMagneticField();
387   Bool_t   bc0 = kFALSE;
388   Int_t    ntr = GetCTSTracks()->GetEntriesFast();
389   //printf("N Tracks %d\n",ntr);
390   
391   for(Int_t i = 0 ; i < ntr ; i++)
392   {
393     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
394     
395     //Check if has TOF info, if not skip
396     ULong_t status  = track->GetStatus();
397     Bool_t  okTOF   = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
398     vertexBC        = track->GetTOFBunchCrossing(bz);
399     Float_t pt      = track->Pt();
400     
401     if(!okTOF) continue;
402     
403     // Get DCA x, y
404     Double_t dca[2]   = {1e6,1e6};
405     Double_t covar[3] = {1e6,1e6,1e6};
406     track->PropagateToDCA(vtx,bz,100.,dca,covar);
407     
408     if(AcceptDCA(pt,dca[0]))
409     {
410       if     (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
411       else if(vertexBC == 0)                                    bc0 = kTRUE;
412     }
413   }
414   
415   if( bc0 ) vertexBC = 0 ;
416   else      vertexBC = AliVTrack::kTOFBCNA ;
417   
418   return vertexBC;
419   
420 }
421
422 //_____________________________
423 void AliCaloTrackReader::Init()
424 {
425   //Init reader. Method to be called in AliAnaPartCorrMaker
426
427   //printf(" AliCaloTrackReader::Init() %p \n",gGeoManager);
428
429   if(fReadStack && fReadAODMCParticles)
430   {
431     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
432     fReadStack          = kFALSE;
433     fReadAODMCParticles = kFALSE;
434   }
435   
436   // Init geometry, I do not like much to do it like this ...
437   if(fImportGeometryFromFile && !gGeoManager) 
438   {
439     if(fImportGeometryFilePath=="") // If not specified, set a default location
440     fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root"; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
441
442     printf("AliCaloTrackReader::Init() - Import %s\n",fImportGeometryFilePath.Data());
443     TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
444   }
445
446   if(!fESDtrackCuts)
447     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); //initialize with TPC only tracks
448
449 }
450
451 //_______________________________________
452 void AliCaloTrackReader::InitParameters()
453 {
454   //Initialize the parameters of the analysis.
455   fDataType   = kESD ;
456   fCTSPtMin   = 0.1 ;
457   fEMCALPtMin = 0.1 ;
458   fPHOSPtMin  = 0.1 ;
459   fCTSPtMax   = 1000. ;
460   fEMCALPtMax = 1000. ;
461   fPHOSPtMax  = 1000. ;
462   
463   //Track DCA cuts
464   // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
465   fTrackDCACut[0] = 0.0105; 
466   fTrackDCACut[1] = 0.0350; 
467   fTrackDCACut[2] = 1.1;
468   
469   //Do not filter the detectors input by default.
470   fFillEMCAL      = kFALSE;
471   fFillPHOS       = kFALSE;
472   fFillCTS        = kFALSE;
473   fFillEMCALCells = kFALSE;
474   fFillPHOSCells  = kFALSE;
475   
476   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
477   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
478   fDeltaAODFileName      = "deltaAODPartCorr.root";
479   fFiredTriggerClassName = "";
480   fEventTriggerMask      = AliVEvent::kAny;
481   fMixEventTriggerMask   = AliVEvent::kAnyINT;
482   fEventTriggerAtSE      = kTRUE; // Use only events that pass event selection at SE base class
483   
484   fAcceptFastCluster = kTRUE;
485   fAnaLED            = kFALSE;
486   
487   //We want tracks fitted in the detectors:
488   //fTrackStatus=AliESDtrack::kTPCrefit;
489   //fTrackStatus|=AliESDtrack::kITSrefit; 
490   fTrackStatus     = 0;
491   fTrackFilterMask = 128; //For AODs, but what is the difference between fTrackStatus and fTrackFilterMask?
492   
493   fESDtrackCuts = 0;
494   fESDtrackComplementaryCuts = 0;
495   
496   fConstrainTrack = kFALSE ; // constrain tracks to vertex
497   
498   fV0ADC[0] = 0;   fV0ADC[1] = 0; 
499   fV0Mul[0] = 0;   fV0Mul[1] = 0; 
500   
501   fZvtxCut   = 10.;
502   
503   fNMixedEvent = 1;
504   
505   fPtHardAndJetPtFactor     = 7.;
506   fPtHardAndClusterPtFactor = 1.;
507
508   //Centrality
509   fCentralityClass  = "V0M";
510   fCentralityOpt    = 10;
511   fCentralityBin[0] = fCentralityBin[1]=-1;
512   
513   fEventPlaneMethod = "V0";
514
515   // Allocate memory (not sure this is the right place)
516   fCTSTracks       = new TObjArray();
517   fEMCALClusters   = new TObjArray();
518   fPHOSClusters    = new TObjArray(); 
519   fTriggerAnalysis = new AliTriggerAnalysis;
520   fAODBranchList   = new TList ;
521
522   fImportGeometryFromFile = kFALSE;
523   
524   fPileUpParamSPD[0] = 3   ; fPileUpParamSPD[1] = 0.8 ;
525   fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
526   
527   // Parametrized time cut (LHC11d)
528   fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 3.5 ; fEMCALParamTimeCutMin[3] = 1.  ;   
529   fEMCALParamTimeCutMax[0] = 5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.45; fEMCALParamTimeCutMax[3] = 1.25;   
530
531   // Parametrized time cut (LHC11c)
532   //fEMCALParamTimeCutMin[0] =-5;   fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;   
533   //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
534   
535   fTimeStampRunMin = -1;
536   fTimeStampRunMax = 1e12;
537   fTimeStampEventFracMin = -1;
538   fTimeStampEventFracMax = 2;
539
540   for(Int_t i = 0; i < 19; i++)
541   {
542     fEMCalBCEvent   [i] = 0;
543     fEMCalBCEventCut[i] = 0;
544     fTrackBCEvent   [i] = 0;
545     fTrackBCEventCut[i] = 0;
546   }
547   
548   // Trigger match-rejection
549   fTriggerPatchTimeWindow[0] = 8;
550   fTriggerPatchTimeWindow[1] = 9;
551   
552   fTriggerClusterBC      = -10000 ;
553   fTriggerEventThreshold =  2.;
554   fTriggerClusterIndex   = -1;
555   fTriggerClusterId      = -1;
556
557 }
558
559 //___________________________________________________________
560 Bool_t AliCaloTrackReader::IsInTimeWindow(const Double_t tof, const Float_t energy) const
561 {
562   // Cluster time selection window
563   
564   // Parametrized cut depending on E
565   if(fUseParamTimeCut)
566   {
567     Float_t minCut= fEMCALParamTimeCutMin[0]+fEMCALParamTimeCutMin[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMin[2])/fEMCALParamTimeCutMin[3]);
568     Float_t maxCut= fEMCALParamTimeCutMax[0]+fEMCALParamTimeCutMax[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMax[2])/fEMCALParamTimeCutMax[3]);
569     //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
570     if( tof < minCut || tof > maxCut )  return kFALSE ;
571   }
572   
573   //In any case, the time should to be larger than the fixed window ...
574   if( tof < fEMCALTimeCutMin  || tof > fEMCALTimeCutMax )  return kFALSE ;
575   
576   return kTRUE ;
577 }
578
579 //________________________________________________
580 Bool_t AliCaloTrackReader::IsPileUpFromSPD() const
581 {
582   // Check if event is from pile-up determined by SPD
583   // Default values: (3, 0.8, 3., 2., 5.)
584   return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] , 
585                                               fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] ); 
586   //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
587
588 }
589
590 //__________________________________________________
591 Bool_t AliCaloTrackReader::IsPileUpFromEMCal() const
592 {
593   // Check if event is from pile-up determined by EMCal
594   if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
595   else                                       return kFALSE;
596 }
597
598 //________________________________________________________
599 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndEMCal() const
600 {
601   // Check if event is from pile-up determined by SPD and EMCal
602   if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
603   else                                          return kFALSE;
604 }
605
606 //_______________________________________________________
607 Bool_t AliCaloTrackReader::IsPileUpFromSPDOrEMCal() const
608 {
609   // Check if event is from pile-up determined by SPD or EMCal
610   if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
611   else                                          return kFALSE;
612 }
613
614 //___________________________________________________________
615 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndNotEMCal() const
616 {
617   // Check if event is from pile-up determined by SPD and not by EMCal
618   if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
619   else                                          return kFALSE;
620 }
621
622 //___________________________________________________________
623 Bool_t AliCaloTrackReader::IsPileUpFromEMCalAndNotSPD() const
624 {
625   // Check if event is from pile-up determined by EMCal, not by SPD
626   if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
627   else                                           return kFALSE;
628 }
629
630 //______________________________________________________________
631 Bool_t AliCaloTrackReader::IsPileUpFromNotSPDAndNotEMCal() const
632 {
633   // Check if event not from pile-up determined neither by SPD nor by EMCal
634   if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
635   else                                            return kFALSE;
636 }
637
638 //________________________________________________________
639 void AliCaloTrackReader::Print(const Option_t * opt) const
640 {
641   
642   //Print some relevant parameters set for the analysis
643   if(! opt)
644     return;
645   
646   printf("***** Print: %s %s ******\n",    GetName(), GetTitle() ) ;
647   printf("Task name      : %s\n",          fTaskName.Data()) ;
648   printf("Data type      : %d\n",          fDataType) ;
649   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
650   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
651   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
652   printf("CTS Max pT     : %2.1f GeV/c\n", fCTSPtMax) ;
653   printf("EMCAL Max pT   : %2.1f GeV/c\n", fEMCALPtMax) ;
654   printf("PHOS Max pT    : %2.1f GeV/c\n", fPHOSPtMax) ;
655   printf("EMCAL Time Cut: %3.1f < TOF  < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
656   printf("Use CTS         =     %d\n",     fFillCTS) ;
657   printf("Use EMCAL       =     %d\n",     fFillEMCAL) ;
658   printf("Use PHOS        =     %d\n",     fFillPHOS) ;
659   printf("Use EMCAL Cells =     %d\n",     fFillEMCALCells) ;
660   printf("Use PHOS  Cells =     %d\n",     fFillPHOSCells) ;
661   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
662   printf("AODs Track filter mask  =  %d or hybrid %d, SPD hit %d\n", (Int_t) fTrackFilterMask,fSelectHybridTracks,fSelectSPDHitTracks) ;
663   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
664   printf("Write delta AOD =     %d\n",     fWriteOutputDeltaAOD) ;
665   printf("Recalculate Clusters = %d\n",    fRecalculateClusters) ;
666   
667   printf("Use Triggers selected in SE base class %d; If not what trigger Mask? %d; Trigger max for mixed %d \n", 
668          fEventTriggerAtSE, fEventTriggerMask,fMixEventTriggerMask);
669   
670   if(fComparePtHardAndClusterPt)
671     printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
672   
673   if(fComparePtHardAndClusterPt)
674     printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
675   
676   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
677   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
678   printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
679   
680   printf("    \n") ;
681   
682
683
684 //_________________________________________________________________________
685 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, 
686                                           const char * /*currentFileName*/) 
687 {
688   //Fill the event counter and input lists that are needed, called by the analysis maker.
689   
690   fEventNumber         = iEntry;
691   fTriggerClusterIndex = -1;
692   fTriggerClusterId    = -1;
693   fIsTriggerMatch      = kFALSE;
694   fTriggerClusterBC    = -10000;
695   fIsExoticEvent       = kFALSE;
696   fIsBadCellEvent      = kFALSE;
697   fIsBadMaxCellEvent   = kFALSE;
698   
699   //fCurrentFileName = TString(currentFileName);
700   if(!fInputEvent)
701   {
702     if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
703     return kFALSE;
704   }
705   
706   //Select events only fired by a certain trigger configuration if it is provided
707   Int_t eventType = 0;
708   if(fInputEvent->GetHeader())
709     eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
710   
711   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster) 
712   {
713     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
714     return kFALSE;
715   }
716   
717   
718   //-------------------------------------------------------------------------------------
719   // Reject event if large clusters with large energy
720   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
721   // If clusterzer NxN or V2 it does not help
722   //-------------------------------------------------------------------------------------
723   Int_t run = fInputEvent->GetRunNumber();
724   if( fRemoveLEDEvents && run > 146857  && run < 146861 )
725   {
726     Bool_t reject = RejectLEDEvents();
727     if(reject) return kFALSE;
728   }// Remove LED events
729   
730   //------------------------
731   // Reject pure LED events?
732   //-------------------------
733   if( fFiredTriggerClassName  !="" && !fAnaLED)
734   {
735     //printf("Event type %d\n",eventType);
736     if(eventType!=7)
737       return kFALSE; //Only physics event, do not use for simulated events!!!
738     
739     if(fDebug > 0) 
740       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
741              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
742     
743     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
744     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
745   }
746   else if(fAnaLED)
747   {
748     //    kStartOfRun =       1,    // START_OF_RUN
749     //    kEndOfRun =         2,    // END_OF_RUN
750     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
751     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
752     //    kStartOfBurst =     5,    // START_OF_BURST
753     //    kEndOfBurst =       6,    // END_OF_BURST
754     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
755     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
756     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
757     //    kStartOfData =      10,   // START_OF_DATA
758     //    kEndOfData =        11,   // END_OF_DATA
759     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
760     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
761     
762           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
763           if(eventType!=8)return kFALSE;
764   }
765   
766   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard 
767   if(fComparePtHardAndJetPt) 
768   {
769     if(!ComparePtHardAndJetPt()) return kFALSE ;
770   }
771   
772   if(fComparePtHardAndClusterPt) 
773   {
774     if(!ComparePtHardAndClusterPt()) return kFALSE ;
775   }
776   
777   //Fill Vertex array
778   FillVertexArray();
779   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
780   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;  
781   
782   //------------------------------------------------------
783   //Event rejection depending on vertex, pileup, v0and
784   //------------------------------------------------------
785   if(fDataType==kESD && fTimeStampEventSelect)
786   {
787     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
788     if(esd)
789     {
790       Int_t timeStamp = esd->GetTimeStamp();
791       Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
792       
793       //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
794       
795       if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
796     }
797     //printf("\t accept time stamp\n");
798   }
799
800   
801   //------------------------------------------------------
802   //Event rejection depending on vertex, pileup, v0and
803   //------------------------------------------------------
804
805   if(fUseEventsWithPrimaryVertex)
806   {
807     if( !CheckForPrimaryVertex() )              return kFALSE;
808     if( TMath::Abs(fVertex[0][0] ) < 1.e-6 && 
809         TMath::Abs(fVertex[0][1] ) < 1.e-6 && 
810         TMath::Abs(fVertex[0][2] ) < 1.e-6    ) return kFALSE;
811   }
812   
813   //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
814   
815   if(fDoEventSelection)
816   {
817     if(!fCaloFilterPatch)
818     {
819       // Do not analyze events with pileup
820       Bool_t bPileup = IsPileUpFromSPD();
821       //IsPileupFromSPDInMultBins() // method to try
822       //printf("pile-up %d, %d, %2.2f, %2.2f, %2.2f, %2.2f\n",bPileup, (Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
823       if(bPileup) return kFALSE;
824       
825       if(fDoV0ANDEventSelection)
826       {
827         Bool_t bV0AND = kTRUE; 
828         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
829         if(esd) 
830           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
831         //else bV0AND = //FIXME FOR AODs
832         if(!bV0AND) return kFALSE;
833       }
834     }//CaloFilter patch
835     else
836     { 
837       if(fInputEvent->GetNumberOfCaloClusters() > 0) 
838       {
839         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
840         if(calo->GetNLabels() == 4)
841         {
842           Int_t * selection = calo->GetLabels();
843           Bool_t bPileup = selection[0];
844           if(bPileup) return kFALSE;
845           
846           Bool_t bGoodV = selection[1]; 
847           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
848           
849           if(fDoV0ANDEventSelection)
850           {
851             Bool_t bV0AND = selection[2]; 
852             if(!bV0AND) return kFALSE;
853           }
854           
855           fTrackMult = selection[3];
856           if(fTrackMult == 0) return kFALSE;
857         } 
858         else 
859         {
860           //First filtered AODs, track multiplicity stored there.  
861           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
862           if(fTrackMult == 0) return kFALSE;          
863         }
864       }//at least one cluster
865       else 
866       {
867         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
868         //Remove events with  vertex (0,0,0), bad vertex reconstruction
869         if(fUseEventsWithPrimaryVertex && TMath::Abs(fVertex[0][0]) < 1.e-6 && TMath::Abs(fVertex[0][1]) < 1.e-6 && TMath::Abs(fVertex[0][2]) < 1.e-6) return kFALSE;
870         
871         //First filtered AODs, track multiplicity stored there.  
872         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
873         if(fTrackMult == 0) return kFALSE;
874       }// no cluster
875     }// CaloFileter patch
876   }// Event selection/AliceSoft/AliRoot/trunk/PWG/CaloTrackCorrBase/AliCaloTrackReader.h
877   
878   //------------------------------------------------------
879   
880   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
881   //If we need a centrality bin, we select only those events in the corresponding bin.
882   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100)
883   {
884     Int_t cen = GetEventCentrality();
885     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
886   }
887     
888   //Fill the arrays with cluster/tracks/cells data
889   
890   if(!fEventTriggerAtSE)
891   {
892     // In case of mixing analysis, accept MB events, not only Trigger
893     // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
894     // via de method in the base class FillMixedEventPool()
895     
896     AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
897     AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
898     
899     if(!inputHandler) return kFALSE ;  // to content coverity
900     
901     UInt_t isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
902     UInt_t isMB      = inputHandler->IsEventSelected() & fMixEventTriggerMask;
903     
904     if(!isTrigger && !isMB) return kFALSE;
905     
906     //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
907   }
908   
909   //----------------------------------------------------------------------
910   // Do not count events that where likely triggered by an exotic cluster
911   // or out BC cluster
912   //----------------------------------------------------------------------
913
914   //Get Patches that triggered
915   TArrayI patches = GetL0TriggerPatches();
916 /*
917   if(fRemoveExoticEvents)
918   {
919     RejectExoticEvents(patches);
920     if(fIsExoticEvent)
921     {
922       //printf("AliCaloTrackReader::FillInputEvent() - REJECT exotic triggered event \n");
923       return kFALSE;
924     }
925   }
926   
927   RejectTriggeredEventsByPileUp(patches);
928   //printf("AliCaloTrackReader::FillInputEvent(), Trigger BC = %d\n",fTriggerClusterBC);
929   
930   if(fRemoveTriggerOutBCEvents)
931   {
932     if(fTriggerClusterBC != 0 && fTriggerClusterBC != 6)
933     {
934       //printf("\t REJECT, bad trigger cluster BC\n");
935       return kFALSE;
936     }
937   }
938 */
939  
940   MatchTriggerCluster(patches);
941   
942   if(fRemoveBadTriggerEvents)
943   {
944     //printf("ACCEPT triggered event? - exotic? %d - bad cell %d - bad Max cell %d - BC %d  - Matched %d\n",
945     //       fIsExoticEvent,fIsBadCellEvent, fIsBadMaxCellEvent, fTriggerClusterBC,fIsTriggerMatch);
946     if     (fIsExoticEvent)         return kFALSE;
947     else if(fIsBadCellEvent)        return kFALSE;
948     else if(fTriggerClusterBC == 0) return kFALSE;
949     //printf("\t *** YES\n");
950   }
951   
952   patches.Reset();
953   
954   // Get the main vertex BC, in case not available
955   // it is calculated in FillCTS checking the BC of tracks
956   // with DCA small (if cut applied, if open)
957   fVertexBC=fInputEvent->GetPrimaryVertex()->GetBC();
958   
959   if(fFillCTS)
960   {
961     FillInputCTS();
962     //Accept events with at least one track
963     if(fTrackMult == 0 && fDoRejectNoTrackEvents) return kFALSE ;
964   }
965   
966   if(fDoVertexBCEventSelection)
967   {
968     if(fVertexBC!=0 && fVertexBC!=AliVTrack::kTOFBCNA) return kFALSE ;
969   }
970   
971   if(fFillEMCALCells)
972     FillInputEMCALCells();
973   
974   if(fFillPHOSCells)
975     FillInputPHOSCells();
976   
977   if(fFillEMCAL)
978     FillInputEMCAL();
979   
980   if(fFillPHOS)
981     FillInputPHOS();
982   
983   FillInputVZERO();
984
985   
986   return kTRUE ;
987 }
988
989 //___________________________________
990 void AliCaloTrackReader::ResetLists() 
991 {
992   //  Reset lists, called by the analysis maker 
993   
994   if(fCTSTracks)       fCTSTracks     -> Clear();
995   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
996   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
997   
998   fV0ADC[0] = 0;   fV0ADC[1] = 0; 
999   fV0Mul[0] = 0;   fV0Mul[1] = 0;
1000   
1001 }
1002
1003 //____________________________________________________________
1004 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)  
1005 {
1006   fInputEvent  = input;
1007   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ; 
1008   if (fMixedEvent) 
1009     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ; 
1010   
1011   //Delete previous vertex
1012   if(fVertex)
1013   {
1014     for (Int_t i = 0; i < fNMixedEvent; i++) 
1015     {
1016       delete [] fVertex[i] ; 
1017     }
1018     delete [] fVertex ;
1019   }
1020   
1021   fVertex = new Double_t*[fNMixedEvent] ; 
1022   for (Int_t i = 0; i < fNMixedEvent; i++) 
1023   {
1024     fVertex[i] = new Double_t[3] ; 
1025     fVertex[i][0] = 0.0 ; 
1026     fVertex[i][1] = 0.0 ; 
1027     fVertex[i][2] = 0.0 ; 
1028   }
1029 }
1030
1031 //__________________________________________________
1032 Int_t AliCaloTrackReader::GetEventCentrality() const 
1033 {
1034   //Return current event centrality
1035   
1036   if(GetCentrality())
1037   {
1038     if     (fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1039     else if(fCentralityOpt==10)  return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1040     else if(fCentralityOpt==20)  return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1041     else 
1042     {
1043       printf("AliCaloTrackReader::GetEventCentrality() - Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
1044       return -1;
1045     } 
1046   }
1047   else return -1;
1048   
1049 }
1050
1051 //_____________________________________________________
1052 Double_t AliCaloTrackReader::GetEventPlaneAngle() const 
1053 {
1054   //Return current event centrality
1055   
1056   if(GetEventPlane())
1057   {
1058     Float_t ep =  GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1059     
1060     if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi())) 
1061     {
1062       if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <Q> method : %f\n",ep);
1063       return -1000;
1064     }
1065     else if(GetEventPlaneMethod().Contains("V0")  ) 
1066     {
1067       if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1068       {
1069         if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep);
1070         return -1000;
1071       }
1072       
1073       ep+=TMath::Pi()/2; // put same range as for <Q> method
1074       
1075     }
1076   
1077     //printf("AliCaloTrackReader::GetEventPlaneAngle() = %f\n",ep);
1078     if(fDebug > 0 )
1079     {
1080       if     (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1081       else if(ep < 0          ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1082     }
1083     
1084     return ep;
1085   }
1086   else
1087   {
1088     if(fDataType!=kMC && fDebug > 0) printf("AliCaloTrackReader::GetEventPlaneAngle() -  No EP pointer\n");
1089     return -1000;
1090   } 
1091   
1092 }
1093
1094 //__________________________________________________________
1095 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const 
1096 {
1097   //Return vertex position to be used for single event analysis
1098   vertex[0]=fVertex[0][0];  
1099   vertex[1]=fVertex[0][1];  
1100   vertex[2]=fVertex[0][2];
1101 }
1102
1103 //____________________________________________________________
1104 void AliCaloTrackReader::GetVertex(Double_t vertex[3], 
1105                                    const Int_t evtIndex) const 
1106 {
1107   //Return vertex position for mixed event, recover the vertex in a particular event.
1108   
1109   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
1110   
1111 }
1112
1113 //________________________________________
1114 void AliCaloTrackReader::FillVertexArray() 
1115 {
1116   
1117   //Fill data member with vertex
1118   //In case of Mixed event, multiple vertices
1119   
1120   //Delete previous vertex
1121   if(fVertex)
1122   {
1123     for (Int_t i = 0; i < fNMixedEvent; i++) 
1124     {
1125       delete [] fVertex[i] ; 
1126     }
1127     delete [] fVertex ;  
1128   }
1129   
1130   fVertex = new Double_t*[fNMixedEvent] ; 
1131   for (Int_t i = 0; i < fNMixedEvent; i++) 
1132   {
1133     fVertex[i] = new Double_t[3] ; 
1134     fVertex[i][0] = 0.0 ; 
1135     fVertex[i][1] = 0.0 ; 
1136     fVertex[i][2] = 0.0 ; 
1137   }          
1138   
1139   if (!fMixedEvent) 
1140   { //Single event analysis
1141     if(fDataType!=kMC)
1142     {
1143       
1144       if(fInputEvent->GetPrimaryVertex())
1145       {
1146         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]); 
1147       }
1148       else 
1149       {
1150         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
1151         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1152       }//Primary vertex pointer do not exist
1153       
1154     } else
1155     {//MC read event 
1156       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1157     }
1158     
1159     if(fDebug > 1)
1160       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
1161     
1162   } else 
1163   { // MultiEvent analysis
1164     for (Int_t iev = 0; iev < fNMixedEvent; iev++) 
1165     {
1166       if (fMixedEvent->GetVertexOfEvent(iev))
1167         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1168       else
1169       { // no vertex found !!!!
1170         AliWarning("No vertex found");
1171       }
1172       
1173       if(fDebug > 1)
1174         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
1175       
1176     }
1177   }
1178   
1179 }
1180
1181 //_____________________________________
1182 void AliCaloTrackReader::FillInputCTS() 
1183 {
1184   //Return array with Central Tracking System (CTS) tracks
1185   
1186   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
1187   
1188   Double_t pTrack[3] = {0,0,0};
1189   
1190   Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1191   fTrackMult    = 0;
1192   Int_t nstatus = 0;
1193   Double_t bz   = GetInputEvent()->GetMagneticField();
1194
1195   for(Int_t i = 0; i < 19; i++)
1196   {
1197     fTrackBCEvent   [i] = 0;
1198     fTrackBCEventCut[i] = 0;
1199   }
1200   
1201   Bool_t   bc0  = kFALSE;
1202   if(fRecalculateVertexBC) fVertexBC=AliVTrack::kTOFBCNA;
1203   
1204   for (Int_t itrack =  0; itrack <  nTracks; itrack++)
1205   {////////////// track loop
1206     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1207     
1208     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1209     ULong_t status = track->GetStatus();
1210
1211     if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1212       continue ;
1213     
1214     nstatus++;
1215     
1216     Float_t dcaTPC =-999;
1217     
1218     if     (fDataType==kESD)
1219     {
1220       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (track);
1221       
1222       if(esdTrack)
1223       {
1224         if(fESDtrackCuts->AcceptTrack(esdTrack))
1225         {
1226           track->GetPxPyPz(pTrack) ;
1227
1228           if(fConstrainTrack)
1229           {
1230             if(esdTrack->GetConstrainedParam())
1231             {
1232               const AliExternalTrackParam* constrainParam = esdTrack->GetConstrainedParam();
1233               esdTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
1234               esdTrack->GetConstrainedPxPyPz(pTrack);
1235             }
1236             else continue;
1237             
1238           } // use constrained tracks
1239           
1240           if(fSelectSPDHitTracks)
1241           {//Not much sense to use with TPC only or Hybrid tracks
1242             if(!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1)) continue ;
1243           }
1244         }
1245         // Complementary track to global : Hybrids (make sure that the previous selection is for Global)
1246         else  if(fESDtrackComplementaryCuts && fESDtrackComplementaryCuts->AcceptTrack(esdTrack))
1247         {
1248           // constrain the track
1249           if(esdTrack->GetConstrainedParam())
1250           {
1251             esdTrack->Set(esdTrack->GetConstrainedParam()->GetX(),esdTrack->GetConstrainedParam()->GetAlpha(),esdTrack->GetConstrainedParam()->GetParameter(),esdTrack->GetConstrainedParam()->GetCovariance());
1252           
1253             track->GetPxPyPz(pTrack) ;
1254
1255           }
1256           else continue;
1257         }
1258         else continue;
1259       }
1260     } // ESD
1261     else if(fDataType==kAOD)
1262     {
1263       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
1264       
1265       if(aodtrack)
1266       {
1267        if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS():AOD track type: %d (primary %d), hybrid? %d \n",
1268                               aodtrack->GetType(),AliAODTrack::kPrimary,
1269                               aodtrack->IsHybridGlobalConstrainedGlobal());
1270         
1271         
1272         if (fSelectHybridTracks)
1273         {
1274           if (!aodtrack->IsHybridGlobalConstrainedGlobal())       continue ;
1275         }
1276         else 
1277         {
1278           if ( aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue ;
1279         }
1280         
1281         if(fSelectSPDHitTracks)
1282         {//Not much sense to use with TPC only or Hybrid tracks
1283           if(!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1)) continue ;
1284         }
1285         
1286         if (aodtrack->GetType()!= AliAODTrack::kPrimary)          continue ;
1287         
1288         if (fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS(): \t accepted track! \n");
1289         
1290         //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1291         // info stored here
1292         dcaTPC = aodtrack->DCA();
1293         
1294         track->GetPxPyPz(pTrack) ;
1295         
1296       } // aod track exists
1297       else continue ;
1298       
1299     } // AOD
1300     
1301     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
1302     
1303     Bool_t okTOF  = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1304     Double_t tof  = -1000;
1305     Int_t trackBC = -1000 ;
1306     
1307     if(okTOF)
1308     {
1309       trackBC = track->GetTOFBunchCrossing(bz);
1310       SetTrackEventBC(trackBC+9);
1311
1312       tof = track->GetTOFsignal()*1e-3;
1313     }
1314                 
1315     if(fUseTrackDCACut)
1316     {
1317       //normal way to get the dca, cut on dca_xy
1318       if(dcaTPC==-999)
1319       {
1320         Double_t dca[2]   = {1e6,1e6};
1321         Double_t covar[3] = {1e6,1e6,1e6};
1322         Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1323         if( okDCA) okDCA = AcceptDCA(momentum.Pt(),dca[0]);
1324         if(!okDCA)
1325         {
1326           //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f, BC %d\n",momentum.Pt(),dca[0],trackBC);
1327           continue ;
1328         }
1329       }
1330     }// DCA cuts
1331     
1332     if(okTOF)
1333     {
1334       //SetTrackEventBCcut(bc);
1335       SetTrackEventBCcut(trackBC+9);
1336       
1337       //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1338       if(fRecalculateVertexBC)
1339       {
1340         if     (trackBC !=0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1341         else if(trackBC == 0)                                  bc0 = kTRUE;
1342       }
1343
1344       //In any case, the time should to be larger than the fixed window ...
1345       if( fUseTrackTimeCut && (trackBC!=0 || tof < fTrackTimeCutMin  || tof > fTrackTimeCutMax) )
1346       {
1347         //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1348         continue ;
1349       }
1350       //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1351       
1352     }
1353     
1354     //Count the tracks in eta < 0.9
1355     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
1356     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
1357     
1358     if(fCTSPtMin > momentum.Pt() || fCTSPtMax < momentum.Pt()) continue ;
1359     
1360     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
1361     
1362     if(fDebug > 2 && momentum.Pt() > 0.1)
1363       printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1364              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1365     
1366     if (fMixedEvent)  track->SetID(itrack);
1367
1368     fCTSTracks->Add(track);
1369     
1370   }// track loop
1371         
1372   if(fVertexBC ==0 || fVertexBC == AliVTrack::kTOFBCNA)
1373   {
1374     if( bc0 ) fVertexBC = 0 ;
1375     else      fVertexBC = AliVTrack::kTOFBCNA ;
1376   }
1377   
1378   if(fDebug > 1)
1379     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
1380   
1381 }
1382
1383 //__________________________________________________________________
1384 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, 
1385                                                  const Int_t iclus) 
1386 {
1387   //Fill the EMCAL data in the array, do it
1388     
1389   Int_t vindex = 0 ;  
1390   if (fMixedEvent) 
1391     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1392   
1393   if(fRecalculateClusters)
1394   {
1395     //Recalibrate the cluster energy
1396     if(GetCaloUtils()->IsRecalibrationOn())
1397     {
1398       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
1399       
1400       clus->SetE(energy);
1401       //printf("Recalibrated Energy %f\n",clus->E());
1402       
1403       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
1404       GetCaloUtils()->RecalculateClusterPID(clus);
1405       
1406     } // recalculate E
1407     
1408     //Recalculate distance to bad channels, if new list of bad channels provided
1409     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
1410     
1411     //Recalculate cluster position
1412     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1413     {
1414       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
1415       //clus->GetPosition(pos);
1416       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1417     }
1418     
1419     // Recalculate TOF
1420     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1421     {
1422       Double_t tof      = clus->GetTOF();
1423       Float_t  frac     =-1;
1424       Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1425       
1426       if(fDataType==AliCaloTrackReader::kESD)
1427       {
1428         tof = fEMCALCells->GetCellTime(absIdMax);
1429       }
1430       
1431       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1432       
1433       clus->SetTOF(tof);
1434       
1435     }// Time recalibration
1436   }
1437   
1438   //Reject clusters with bad channels, close to borders and exotic;
1439   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
1440   
1441   //Mask all cells in collumns facing ALICE thick material if requested
1442   if(GetCaloUtils()->GetNMaskCellColumns())
1443   {
1444     Int_t absId   = -1;
1445     Int_t iSupMod = -1;
1446     Int_t iphi    = -1;
1447     Int_t ieta    = -1;
1448     Bool_t shared = kFALSE;
1449     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
1450     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
1451   }
1452   
1453   if(fSelectEmbeddedClusters)
1454   {
1455     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1456     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
1457   }
1458   
1459   //Float_t pos[3];
1460   //clus->GetPosition(pos);
1461   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1462   
1463   //Correct non linearity
1464   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn())
1465   {
1466     GetCaloUtils()->CorrectClusterEnergy(clus) ;
1467     //printf("Linearity Corrected Energy %f\n",clus->E());  
1468     
1469     //In case of MC analysis, to match resolution/calibration in real data
1470     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1471     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
1472     clus->SetE(rdmEnergy);
1473   }
1474   
1475   Double_t tof = clus->GetTOF()*1e9;
1476
1477   Int_t bc = TMath::Nint(tof/50) + 9;
1478   //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1479   
1480   SetEMCalEventBC(bc);
1481   
1482   if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) return ;
1483
1484   TLorentzVector momentum ;
1485   
1486   clus->GetMomentum(momentum, fVertex[vindex]);
1487   
1488   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1489   
1490   SetEMCalEventBCcut(bc);
1491
1492   if(!IsInTimeWindow(tof,clus->E()))
1493   {
1494     fNPileUpClusters++ ;
1495     if(fUseEMCALTimeCut) return ;
1496   }
1497   else
1498     fNNonPileUpClusters++;
1499   
1500   if(fDebug > 2 && momentum.E() > 0.1) 
1501     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1502            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1503   
1504   if (fMixedEvent) 
1505     clus->SetID(iclus) ; 
1506   
1507   //Correct MC label for AODs
1508   
1509   fEMCALClusters->Add(clus);
1510   
1511 }
1512
1513 //_______________________________________
1514 void AliCaloTrackReader::FillInputEMCAL() 
1515 {
1516   //Return array with EMCAL clusters in aod format
1517   
1518   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
1519   
1520   // First recalibrate cells, time or energy
1521   //  if(GetCaloUtils()->IsRecalibrationOn())
1522   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(), 
1523   //                                                          GetEMCALCells(), 
1524   //                                                          fInputEvent->GetBunchCrossNumber());
1525   
1526   fNPileUpClusters    = 0; // Init counter
1527   fNNonPileUpClusters = 0; // Init counter
1528   for(Int_t i = 0; i < 19; i++)
1529   {
1530     fEMCalBCEvent   [i] = 0;
1531     fEMCalBCEventCut[i] = 0;
1532   }
1533   
1534   //Loop to select clusters in fiducial cut and fill container with aodClusters
1535   if(fEMCALClustersListName=="")
1536   {
1537     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1538     for (Int_t iclus =  0; iclus <  nclusters; iclus++) 
1539     {
1540       AliVCluster * clus = 0;
1541       if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) 
1542       {
1543         if (IsEMCALCluster(clus))
1544         {          
1545           FillInputEMCALAlgorithm(clus, iclus);
1546         }//EMCAL cluster
1547       }// cluster exists
1548     }// cluster loop
1549     
1550     //Recalculate track matching
1551     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
1552     
1553   }//Get the clusters from the input event
1554   else
1555   {
1556     TClonesArray * clusterList = 0x0; 
1557     
1558     if      (fInputEvent->FindListObject(fEMCALClustersListName))
1559     {
1560       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1561     }
1562     else if(fOutputEvent)
1563     {
1564       clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1565     }
1566     
1567     if(!clusterList)
1568     {
1569         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
1570         return;
1571     }
1572     
1573     Int_t nclusters = clusterList->GetEntriesFast();
1574     for (Int_t iclus =  0; iclus <  nclusters; iclus++) 
1575     {
1576       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1577       //printf("E %f\n",clus->E());
1578       if (clus) FillInputEMCALAlgorithm(clus, iclus);
1579       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
1580     }// cluster loop
1581     
1582     // Recalculate the pile-up time, in case long time clusters removed during clusterization
1583     //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1584
1585     fNPileUpClusters    = 0; // Init counter
1586     fNNonPileUpClusters = 0; // Init counter
1587     for(Int_t i = 0; i < 19; i++)
1588     {
1589       fEMCalBCEvent   [i] = 0;
1590       fEMCalBCEventCut[i] = 0;
1591     }
1592     
1593     for (Int_t iclus =  0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
1594     {
1595       AliVCluster * clus = 0;
1596       
1597       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1598       {
1599         if (IsEMCALCluster(clus))
1600         {
1601           
1602           Float_t  frac     =-1;
1603           Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1604           Double_t tof = clus->GetTOF();
1605           GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1606           tof*=1e9;
1607
1608           //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
1609
1610           //Reject clusters with bad channels, close to borders and exotic;
1611           if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber()))  continue;
1612
1613           Int_t bc = TMath::Nint(tof/50) + 9;
1614           SetEMCalEventBC(bc);
1615           
1616           if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
1617
1618           TLorentzVector momentum ;
1619           
1620           clus->GetMomentum(momentum, fVertex[0]);
1621           
1622           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1623
1624           SetEMCalEventBCcut(bc);
1625
1626           if(!IsInTimeWindow(tof,clus->E()))
1627             fNPileUpClusters++ ;
1628           else
1629             fNNonPileUpClusters++;
1630           
1631         }
1632       }
1633     }
1634     
1635     //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1636     
1637     // Recalculate track matching, not necessary if already done in the reclusterization task.
1638     // in case it was not done ...
1639     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
1640     
1641   }
1642     
1643   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d, n pile-up clusters %d, n non pile-up %d \n",  fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters);
1644   
1645 }
1646
1647 //______________________________________
1648 void AliCaloTrackReader::FillInputPHOS() 
1649 {
1650   //Return array with PHOS clusters in aod format
1651   
1652   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
1653   
1654   //Loop to select clusters in fiducial cut and fill container with aodClusters
1655   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1656   for (Int_t iclus = 0; iclus < nclusters; iclus++) 
1657   {
1658     AliVCluster * clus = 0;
1659     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) 
1660     {
1661       if (IsPHOSCluster(clus))
1662       {
1663         //Check if the cluster contains any bad channel and if close to calorimeter borders
1664         Int_t vindex = 0 ;  
1665         if (fMixedEvent) 
1666           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1667         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) 
1668           continue;
1669         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex)) 
1670           continue;
1671         
1672         if(fRecalculateClusters)
1673         {
1674           //Recalibrate the cluster energy 
1675           if(GetCaloUtils()->IsRecalibrationOn()) 
1676           {
1677             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1678             clus->SetE(energy);
1679           }
1680         }
1681         
1682         TLorentzVector momentum ;
1683         
1684         clus->GetMomentum(momentum, fVertex[vindex]);      
1685         
1686         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
1687         
1688         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
1689         
1690         if(fDebug > 2 && momentum.E() > 0.1) 
1691           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1692                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());        
1693         
1694         
1695         if (fMixedEvent) 
1696         {
1697           clus->SetID(iclus) ; 
1698         }              
1699         
1700         fPHOSClusters->Add(clus);       
1701         
1702       }//PHOS cluster
1703     }//cluster exists
1704   }//esd cluster loop
1705   
1706   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1707   
1708 }
1709
1710 //____________________________________________
1711 void AliCaloTrackReader::FillInputEMCALCells() 
1712 {
1713   //Return array with EMCAL cells in aod format
1714   
1715   fEMCALCells = fInputEvent->GetEMCALCells(); 
1716   
1717 }
1718
1719 //___________________________________________
1720 void AliCaloTrackReader::FillInputPHOSCells() 
1721 {
1722   //Return array with PHOS cells in aod format
1723   
1724   fPHOSCells = fInputEvent->GetPHOSCells(); 
1725   
1726 }
1727
1728 //_______________________________________
1729 void AliCaloTrackReader::FillInputVZERO()
1730 {
1731   //Fill VZERO information in data member, add all the channels information.
1732   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1733   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1734   
1735   if (v0) 
1736   {
1737     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1738     for (Int_t i = 0; i < 32; i++)
1739     {
1740       if(esdV0)
1741       {//Only available in ESDs
1742         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1743         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1744       }
1745       
1746       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1747       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1748     }
1749     if(fDebug > 0)
1750       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1751   }
1752   else
1753   {
1754     if(fDebug > 0)
1755       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1756   }
1757 }
1758
1759
1760 //___________________________________________________________________
1761 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const 
1762 {
1763   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1764   // different number and need to patch here
1765   
1766   if(fDataType==kAOD && fOldAOD)
1767   {
1768     if (cluster->GetType() == 2) return kTRUE;
1769     else                         return kFALSE;
1770   }
1771   else 
1772   {
1773     return cluster->IsEMCAL();
1774   }
1775   
1776 }
1777
1778 //___________________________________________________________________
1779 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const 
1780 {
1781   //Check if it is a cluster from PHOS.For old AODs cluster type has
1782   // different number and need to patch here
1783   
1784   if(fDataType==kAOD && fOldAOD)
1785   {
1786     Int_t type = cluster->GetType();
1787     if (type == 0 || type == 1) return kTRUE;
1788     else                        return kFALSE;
1789   }
1790   else 
1791   {
1792     return cluster->IsPHOS();
1793   }
1794   
1795 }
1796
1797 //________________________________________________
1798 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1799 {
1800   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1801   //Only for ESDs ...
1802   
1803   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1804   if(!event) return kTRUE;
1805   
1806   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
1807   {
1808     return kTRUE;
1809   }
1810   
1811   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) 
1812   {
1813     // SPD vertex
1814     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) 
1815     {
1816       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1817       return kTRUE;
1818       
1819     }
1820     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
1821     {
1822       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1823       return kFALSE;
1824     }
1825   }
1826   
1827   return kFALSE;  
1828   
1829 }
1830
1831 //_______________________________________________
1832 TArrayI AliCaloTrackReader::GetL0TriggerPatches()
1833 {
1834   //Select the patches that triggered
1835   
1836   // some variables
1837   Int_t  trigtimes[30], globCol, globRow,ntimes, i;
1838   Int_t  absId  = -1; //[100];
1839   Int_t  nPatch = 0;
1840   
1841   TArrayI patches(0);
1842   
1843   // get object pointer
1844   AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
1845   
1846   // class is not empty
1847   if( caloTrigger->GetEntries() > 0 )
1848   {
1849     // must reset before usage, or the class will fail
1850     caloTrigger->Reset();
1851     
1852     // go throuth the trigger channels
1853     while( caloTrigger->Next() )
1854     {
1855       // get position in global 2x2 tower coordinates
1856       caloTrigger->GetPosition( globCol, globRow );
1857       
1858       // get dimension of time arrays
1859       caloTrigger->GetNL0Times( ntimes );
1860       
1861       // no L0s in this channel
1862       // presence of the channel in the iterator still does not guarantee that L0 was produced!!
1863       if( ntimes < 1 )
1864         continue;
1865       
1866       // get timing array
1867       caloTrigger->GetL0Times( trigtimes );
1868       
1869       //printf("trigger time window %d - %d\n",fTriggerPatchTimeWindow[0],fTriggerPatchTimeWindow[1]);
1870       // go through the array
1871       for( i = 0; i < ntimes; i++ )
1872       {
1873         // check if in cut - 8,9 shall be accepted in 2011
1874         if( trigtimes[i] >= fTriggerPatchTimeWindow[0] && trigtimes[i] <= fTriggerPatchTimeWindow[1] )
1875         {
1876           //printf("Accepted trigger time %d \n",trigtimes[i]);
1877           //if(nTrig > 99) continue;
1878           GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);//absIDTrig[nTrig]) ;
1879           //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absIDTrig[nTrig]);
1880           patches.Set(nPatch+1);
1881           patches.AddAt(absId,nPatch++);
1882         }
1883       } // trigger time array
1884       
1885     } // trigger iterator
1886   } // go thorough triggers
1887
1888   //printf("N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nTrig, patches.At(0), patches.At(patches.GetSize()-1));
1889   
1890   return patches;
1891 }
1892
1893
1894 //______________________________________________________________________
1895 void  AliCaloTrackReader::MatchTriggerCluster(TArrayI patches)
1896 {
1897   // Finds the cluster that triggered
1898   
1899   // Init info from previous event
1900   fTriggerClusterIndex = -1;
1901   fTriggerClusterId    = -1;
1902   fIsTriggerMatch      = kFALSE;
1903   fTriggerClusterBC    = -10000;
1904   fIsExoticEvent       = kFALSE;
1905   fIsBadCellEvent      = kFALSE;
1906   fIsBadMaxCellEvent   = kFALSE;
1907   
1908   // Do only analysis for triggered events
1909   if(!GetFiredTriggerClasses().Contains("EMC"))
1910   {
1911     fTriggerClusterBC = 0;
1912     return;
1913   }
1914   
1915   //Recover the list of clusters
1916   TClonesArray * clusterList = 0;
1917   if      (fInputEvent->FindListObject(fEMCALClustersListName))
1918   {
1919     clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1920   }
1921   else if(fOutputEvent)
1922   {
1923     clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1924   }
1925   
1926   // Get number of clusters and of trigger patches
1927   Int_t   nclusters   = fInputEvent->GetNumberOfCaloClusters();
1928   if(clusterList)
1929          nclusters    = clusterList->GetEntriesFast();
1930
1931   Int_t   nPatch      = patches.GetSize();
1932   Float_t exoDiffTime = GetCaloUtils()->GetEMCALRecoUtils()->GetExoticCellDiffTimeCut();
1933
1934   //Init some variables used in the cluster loop
1935   Float_t tofPatchMax = 100000;
1936   Float_t ePatchMax   =-1;
1937   
1938   Float_t tofMax      = 100000;
1939   Float_t eMax        =-1;
1940   
1941   Int_t   clusMax     =-1;
1942   Int_t   idclusMax   =-1;
1943   Bool_t  badClMax    = kFALSE;
1944   Bool_t  badCeMax    = kFALSE;
1945   Bool_t  exoMax      = kFALSE;
1946   
1947   Int_t   nOfHighECl  = 0 ;
1948
1949   // Loop on the clusters, check if there is any that falls into one of the patches
1950   for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1951   {
1952     AliVCluster * clus = 0;
1953     if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
1954     else            clus = fInputEvent->GetCaloCluster(iclus);
1955     
1956     if ( !clus )                continue ;
1957     
1958     if ( !IsEMCALCluster(clus)) continue ;
1959       
1960     //Skip clusters with too low energy to be triggering
1961     if ( clus->E() < fTriggerEventThreshold / 2. ) continue ;
1962     
1963     Float_t  frac       = -1;
1964     Int_t    absIdMax   = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
1965     
1966     Bool_t   badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
1967                                                                                    clus->GetCellsAbsId(),clus->GetNCells());
1968     UShort_t cellMax[]  = {absIdMax};
1969     Bool_t   badCell    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
1970     
1971     // if cell is bad, it can happen that time calibration is not available,
1972     // when calculating if it is exotic, this can make it to be exotic by default
1973     // open it temporarily for this cluster
1974     if(badCell)
1975      GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(10000000);
1976     
1977     Bool_t   exotic     = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
1978     
1979     // Set back the time cut on exotics
1980     if(badCell)
1981       GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(exoDiffTime);
1982     
1983     // Energy threshold for exotic Ecross typically at 4 GeV,
1984     // for lower energy, check that there are more than 1 cell in the cluster
1985     if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
1986     
1987     Float_t  energy     = clus->E();
1988     Int_t    idclus     = clus->GetID();
1989     
1990     Double_t tof        = clus->GetTOF();
1991     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1992       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1993     tof *=1.e9;
1994     
1995 //    printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
1996 //           iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
1997
1998     // Find the highest energy cluster, avobe trigger threshold
1999     // in the event in case no match to trigger is found
2000     if( energy > eMax )
2001     {
2002       tofMax    = tof;
2003       eMax      = energy;
2004       badClMax  = badCluster;
2005       badCeMax  = badCell;
2006       exoMax    = exotic;
2007       clusMax   = iclus;
2008       idclusMax = idclus;
2009     }
2010     
2011     // count the good clusters in the event avobe the trigger threshold
2012     // to check the exotic events 
2013     if(!badCluster && !exotic)
2014       nOfHighECl++;
2015     
2016     // Find match to trigger
2017     if(fTriggerPatchClusterMatch)
2018     {
2019       for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2020       {
2021         Int_t absIDCell[4];
2022         GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2023         //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2024         //                     clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2025         
2026         for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2027         {
2028           if(absIdMax == absIDCell[ipatch])
2029           {
2030             //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2031             if(energy > ePatchMax)
2032             {
2033               tofPatchMax          = tof;
2034               ePatchMax            = energy;
2035               fIsBadCellEvent      = badCluster;
2036               fIsBadMaxCellEvent   = badCell;
2037               fIsExoticEvent       = exotic;
2038               fTriggerClusterIndex = iclus;
2039               fTriggerClusterId    = idclus;
2040               fIsTriggerMatch      = kTRUE;
2041             }
2042           }
2043         }// cell patch loop
2044       }// trigger patch loop
2045     } // Do trigger patch matching
2046     
2047   }// Cluster loop
2048
2049   // If there was no match, assign as trigger
2050   // the highest energy cluster in the event
2051   if(!fIsTriggerMatch)
2052   {
2053     tofPatchMax          = tofMax;
2054     ePatchMax            = eMax;
2055     fIsBadCellEvent      = badClMax;
2056     fIsBadMaxCellEvent   = badCeMax;
2057     fIsExoticEvent       = exoMax;
2058     fTriggerClusterIndex = clusMax;
2059     fTriggerClusterId    = idclusMax;
2060   }
2061   
2062   Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2063     
2064   if     (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2065   else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2066   else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2067   else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2068   else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2069   else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2070   else
2071   {
2072     //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2073     if(fTriggerClusterIndex >= 0) fTriggerClusterBC = 6 ;
2074     else
2075     {
2076       fTriggerClusterIndex = -2;
2077       fTriggerClusterId    = -2;
2078     }
2079   }
2080   
2081   if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2082   
2083 //  printf("AliCaloTrackReader::MatchTriggerCluster(TArrayI patches) - Trigger cluster: index %d, ID %d, E = %2.2f, tof = %2.2f (BC = %d), bad cluster? %d, bad cell? %d, exotic? %d, patch match? %d, n High E cluster %d\n",
2084 //         fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2085 //         fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl);
2086 //  
2087 //  if(!fIsTriggerMatch)  printf("\t highest energy cluster:  index %d, ID %d, E = %2.2f, tof = %2.2f, bad cluster? %d, bad cell? %d, exotic? %d\n",clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2088 }
2089
2090 //__________________________________________
2091 Bool_t  AliCaloTrackReader::RejectLEDEvents()
2092 {
2093   // LED Events in period LHC11a contaminated sample, simple method
2094   // to reject such events
2095   
2096   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
2097   Int_t ncellsSM3 = 0;
2098   for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
2099   {
2100     Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
2101     Int_t sm    = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
2102     if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
2103   }
2104   
2105   Int_t ncellcut = 21;
2106   if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
2107     
2108   if(ncellsSM3 >= ncellcut)
2109   {
2110     if(fDebug > 0)
2111       printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d, cut %d, trig %s\n",
2112              ncellsSM3,ncellcut,GetFiredTriggerClasses().Data());
2113     return kTRUE;
2114   }
2115   
2116   return kFALSE;
2117   
2118 }
2119
2120 //____________________________________________________________
2121 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
2122
2123   // Set Track cuts
2124   
2125   if(fESDtrackCuts) delete fESDtrackCuts ;
2126   
2127   fESDtrackCuts = cuts ; 
2128   
2129 }                 
2130
2131 //_________________________________________________________________________
2132 void  AliCaloTrackReader::SetTrackComplementaryCuts(AliESDtrackCuts * cuts)
2133 {
2134   // Set Track cuts for complementary tracks (hybrids)
2135   
2136   if(fESDtrackComplementaryCuts) delete fESDtrackComplementaryCuts ;
2137   
2138   fESDtrackComplementaryCuts = cuts ;
2139   
2140 }
2141
2142