Additional protections, corrected ownership of the reco. parameters (needed to run...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom.h>
25 #include <TString.h>
26 #include <TSystem.h>
27 #include <TInterpreter.h>
28 #include <TChain.h>
29 #include <TFile.h>
30 #include <TKey.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3F.h>
34 #include <TProfile.h>
35 #include <TList.h>
36 #include <TLorentzVector.h>
37 #include <TRandom3.h>
38 #include <TClonesArray.h>
39 #include "TDatabasePDG.h"
40
41 #include "AliAnalysisTaskJetServices.h"
42 #include "AliCentrality.h"
43 #include "AliAnalysisDataContainer.h"
44 #include "AliAnalysisDataSlot.h"
45 #include "AliAnalysisManager.h"
46 #include "AliJetFinder.h"
47 #include "AliJetHeader.h"
48 #include "AliJetReader.h"
49 #include "AliJetReaderHeader.h"
50 #include "AliUA1JetHeaderV1.h"
51 #include "AliESDEvent.h"
52 #include "AliAODEvent.h"
53 #include "AliAODHandler.h"
54 #include "AliAODTrack.h"
55 #include "AliAODJet.h"
56 #include "AliAODMCParticle.h"
57 #include "AliMCEventHandler.h"
58 #include "AliMCEvent.h"
59 #include "AliStack.h"
60 #include "AliGenPythiaEventHeader.h"
61 #include "AliJetKineReaderHeader.h"
62 #include "AliGenCocktailEventHeader.h"
63 #include "AliInputEventHandler.h"
64 #include "AliPhysicsSelection.h"
65 #include "AliTriggerAnalysis.h"
66
67 #include "AliAnalysisHelperJetTasks.h"
68
69 ClassImp(AliAnalysisTaskJetServices)
70
71 AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
72 TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
73
74 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): 
75   AliAnalysisTaskSE(),
76   fUseAODInput(kFALSE),
77   fUsePhysicsSelection(kFALSE),
78   fMC(kFALSE),
79   fFilterAODCollisions(kFALSE),
80   fPhysicsSelectionFlag(AliVEvent::kMB),
81   fSelectionInfoESD(0),
82   fEventCutInfoESD(0),
83   fFilterMask(0),
84   fRPSubeventMethod(0),
85   fCollisionType(kPbPb),
86   fAvgTrials(1),
87   fVtxXMean(0),
88   fVtxYMean(0),
89   fVtxZMean(0),
90   fVtxRCut(1.),
91   fVtxZCut(8.),
92   fPtMinCosmic(5.),
93   fRIsolMinCosmic(3.),
94   fMaxCosmicAngle(0.01),
95   fCentrality(101),
96   fTrackRecEtaWindow(0.9),
97   fMinTrackPt(0.15),
98   fRPAngle(0),
99   fRandomizer(0),
100   fNonStdFile(""),
101   fh1Xsec(0x0),
102   fh1Trials(0x0),
103   fh1PtHard(0x0),
104   fh1PtHardTrials(0x0),
105   fh1SelectionInfoESD(0x0),
106   fh1EventCutInfoESD(0),
107   fh1CentralityESD(0),
108   fh1Centrality(0),
109   fh1RP(0),
110   fh2TriggerCount(0x0),
111   fh2ESDTriggerCount(0x0),
112   fh2TriggerVtx(0x0),
113   fh2ESDTriggerVtx(0x0),
114   fh2ESDTriggerRun(0x0),
115   fh2VtxXY(0x0),
116   fh1NCosmicsPerEvent(0x0),
117   fh2RPSubevents(0x0),
118   fh2RPCentrality(0x0),
119   fh2RPDeltaRP(0x0),
120   fh2RPQxQy(0x0),
121   fh2RPCosDeltaRP(0x0),
122   fh3PhiWeights(0x0),
123   fh3RPPhiTracks(0x0),
124   fTriggerAnalysis(0x0),
125   fHistList(0x0)  
126 {
127   fRunRange[0] = fRunRange[1] = 0; 
128   fFlatA[0] =   fFlatA[1] = 0;
129   fFlatB[0] =   fFlatB[1] = 0;
130   fDeltaQxy[0] =   fDeltaQxy[1] = 0; 
131
132 }
133
134 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
135   AliAnalysisTaskSE(name),
136   fUseAODInput(kFALSE),
137   fUsePhysicsSelection(kFALSE),
138   fMC(kFALSE),
139   fFilterAODCollisions(kFALSE),
140   fPhysicsSelectionFlag(AliVEvent::kMB),
141   fSelectionInfoESD(0),
142   fEventCutInfoESD(0),
143   fFilterMask(0),
144   fRPSubeventMethod(0),
145   fCollisionType(kPbPb),
146   fAvgTrials(1),
147   fVtxXMean(0),
148   fVtxYMean(0),
149   fVtxZMean(0),
150   fVtxRCut(1.),
151   fVtxZCut(8.),
152   fPtMinCosmic(5.),
153   fRIsolMinCosmic(3.),
154   fMaxCosmicAngle(0.01),
155   fCentrality(101),
156   fTrackRecEtaWindow(0.9),
157   fMinTrackPt(0.15),
158   fRPAngle(0),
159   fRandomizer(0),
160   fNonStdFile(""),
161   fh1Xsec(0x0),
162   fh1Trials(0x0),
163   fh1PtHard(0x0),
164   fh1PtHardTrials(0x0),
165   fh1SelectionInfoESD(0x0),
166   fh1EventCutInfoESD(0),
167   fh1CentralityESD(0),
168   fh1Centrality(0),
169   fh1RP(0),
170   fh2TriggerCount(0x0),
171   fh2ESDTriggerCount(0x0),
172   fh2TriggerVtx(0x0),
173   fh2ESDTriggerVtx(0x0),
174   fh2ESDTriggerRun(0x0),
175   fh2VtxXY(0x0),
176   fh1NCosmicsPerEvent(0x0),
177   fh2RPSubevents(0x0),
178   fh2RPCentrality(0x0),
179   fh2RPDeltaRP(0x0),
180   fh2RPQxQy(0x0),
181   fh2RPCosDeltaRP(0x0),
182   fh3PhiWeights(0x0),
183   fh3RPPhiTracks(0x0),
184
185   fTriggerAnalysis(0x0),
186   fHistList(0x0)  
187 {
188   fRunRange[0] = fRunRange[1] = 0; 
189   fFlatA[0] =   fFlatA[1] = 0;
190   fFlatB[0] =   fFlatB[1] = 0;
191   fDeltaQxy[0] =   fDeltaQxy[1] = 0; 
192   DefineOutput(1,TList::Class());
193 }
194
195
196
197 Bool_t AliAnalysisTaskJetServices::Notify()
198 {
199   //
200   // Implemented Notify() to read the cross sections
201   // and number of trials from pyxsec.root
202   // 
203
204   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
205   Float_t xsection = 0;
206   Float_t ftrials  = 1;
207
208   fAvgTrials = 1;
209   if(tree){
210     TFile *curfile = tree->GetCurrentFile();
211     if (!curfile) {
212       Error("Notify","No current file");
213       return kFALSE;
214     }
215     if(!fh1Xsec||!fh1Trials){
216       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
217       return kFALSE;
218     }
219     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
220     fh1Xsec->Fill("<#sigma>",xsection);
221     // construct a poor man average trials 
222     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
223     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
224   }  
225   return kTRUE;
226 }
227
228 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
229 {
230
231   //
232   // Create the output container
233   //
234
235   fRandomizer = new TRandom3(0);
236   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
237
238   OpenFile(1);
239   if(!fHistList)fHistList = new TList();
240   fHistList->SetOwner();
241   PostData(1, fHistList); // post data in any case once
242
243   Bool_t oldStatus = TH1::AddDirectoryStatus();
244   TH1::AddDirectory(kFALSE);
245   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
246   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
247   fHistList->Add(fh1Xsec);
248
249   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
250   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
251   fHistList->Add(fh1Trials);
252
253   const Int_t nBinPt = 100;
254   Double_t binLimitsPt[nBinPt+1];
255   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
256     if(iPt == 0){
257       binLimitsPt[iPt] = 0.0;
258     }
259     else {// 1.0
260       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
261     }
262   }
263   
264
265   fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
266   fHistList->Add(fh1CentralityESD);
267   
268   fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
269   fHistList->Add(fh1Centrality);
270
271   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
272   fHistList->Add(fh1RP);
273
274   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5); 
275   fHistList->Add(fh2TriggerCount);
276
277   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5); 
278   fHistList->Add(fh2ESDTriggerCount);
279   const Int_t nBins = 6*kConstraints;
280   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
281   fHistList->Add(fh2TriggerVtx);
282
283   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
284   fHistList->Add(fh2ESDTriggerVtx);
285   
286
287   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
288   fHistList->Add(fh1PtHard);
289   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
290   fHistList->Add(fh1PtHardTrials);
291   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
292                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
293   fHistList->Add(fh1SelectionInfoESD);
294
295   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
296                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
297   fHistList->Add(fh1EventCutInfoESD);
298
299   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
300   // 3 triggers BB BE/EB EE
301
302   fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Eventclass vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
303   fHistList->Add(fh2ESDTriggerRun);
304
305   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
306   fHistList->Add(fh2VtxXY);
307   // =========== Switch on Sumw2 for all histos ===========
308   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
309     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
310     if (h1){
311       h1->Sumw2();
312       continue;
313     }
314     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
315     if(hn)hn->Sumw2();
316   }
317
318   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
319
320   
321   fh2RPSubevents = new TH2F("fh2RPSubevents" ,"Reaction Plane Angle" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
322   fHistList->Add( fh2RPSubevents);
323
324   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
325   fHistList->Add(fh2RPCentrality);
326
327   fh2RPDeltaRP   = new TH2F("fh2DeltaRP" ,"Delta Reaction Plane Angle" , 100, -TMath::Pi()/2, TMath::Pi()/2,20,0.,100.0);
328   fHistList->Add(fh2RPDeltaRP);
329
330   fh2RPQxQy      = new TH2F("fh2RPQxQy" ,"" , 100, -100,100,100,-100,100);
331   fHistList->Add(fh2RPQxQy);
332
333   fh2RPCosDeltaRP = new TH2F("fh2RPCosDeltaRP" ,"" , 20, 0.001,100.001,100,-1,1);
334   fHistList->Add(fh2RPCosDeltaRP);
335
336   fh3RPPhiTracks = new TH3F("fh3RPPhiTracks","Phi Tracks Pt Centrality", 10, 0.,100.,20,-5,5,180, 0, 2*TMath::Pi());
337   fHistList->Add(fh3RPPhiTracks);
338   
339
340   fHistList->Add(fh1NCosmicsPerEvent),
341
342
343   TH1::AddDirectory(oldStatus);
344
345   // Add an AOD branch for replication
346   if(fNonStdFile.Length()){
347      if (fDebug > 1) AliInfo("Replicating header");
348      fgAODHeader = new AliAODHeader;
349      AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
350      if (fDebug > 1) AliInfo("Replicating primary vertices");
351      fgAODVertices = new TClonesArray("AliAODVertex",3);
352      fgAODVertices->SetName("vertices");
353      AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
354   }
355 }
356
357 void AliAnalysisTaskJetServices::Init()
358 {
359   //
360   // Initialization
361   //
362   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
363
364 }
365
366 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
367 {
368
369   //
370   // Execute analysis for current event
371   //
372  
373   AliAODEvent *aod = 0;
374   AliESDEvent *esd = 0;
375
376
377   
378   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
379   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
380   AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
381   fSelectionInfoESD = 0; // reset
382   fEventCutInfoESD = 0; // reset
383   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
384
385
386   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
387    
388
389
390   if(fUseAODInput){    
391     aod = dynamic_cast<AliAODEvent*>(InputEvent());
392     if(!aod){
393       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
394       return;
395     }
396     // fethc the header
397   }
398   else{
399     //  assume that the AOD is in the general output...
400     aod  = AODEvent();
401     if(!aod){
402       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
403       return;
404     }
405     esd = dynamic_cast<AliESDEvent*>(InputEvent());
406   }
407   if(aod&&fDebug>2){
408     aod->Print();
409     Printf("Vertices %d",aod->GetNumberOfVertices());
410     Printf("tracks %d",aod->GetNumberOfTracks());
411     Printf("jets %d",aod->GetNJets());
412   }
413   fSelectionInfoESD |= kNoEventCut;
414   fEventCutInfoESD |= kNoEventCut;
415
416   Bool_t esdVtxValid = false;
417   Bool_t esdVtxIn = false;
418   Bool_t aodVtxValid = false;
419   Bool_t aodVtxIn = false;
420
421   if(esd){
422     // trigger analyisis
423     if(!fTriggerAnalysis){
424       fTriggerAnalysis = new AliTriggerAnalysis;
425       fTriggerAnalysis->SetAnalyzeMC(fMC);
426       fTriggerAnalysis->SetSPDGFOThreshhold(1);
427     }
428     //    fTriggerAnalysis->FillTriggerClasses(esd);
429     Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
430     Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
431     Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
432     Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
433     Bool_t spdFO      = fTriggerAnalysis->SPDFiredChips(esd, 0);;
434     if(v0A)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0A;
435     if(v0C)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0C;
436     if(!(v0ABG||v0CBG))fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kNoV0BG;
437     if(spdFO)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kSPDFO;
438   }
439  
440   // Apply additional constraints
441   Bool_t esdEventSelected = IsEventSelected(esd);
442   Bool_t esdEventPileUp = IsEventPileUp(esd);
443   Bool_t esdEventCosmic = IsEventCosmic(esd);
444
445   Bool_t aodEventSelected = IsEventSelected(aod);
446
447   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
448
449
450   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
451   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
452
453   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
454   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
455   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
456   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
457
458
459   // here we have all selection information, fill histogram
460   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
461     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
462   }
463   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
464
465   if(esd&&aod&&fDebug){
466     if(esdEventSelected&&!aodEventSelected){
467       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
468       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
469       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
470       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
471       vtxESD->Print();
472       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
473       vtxAOD->Print();
474     }
475   }
476
477   // loop over all possible triggers for esd 
478
479   Float_t cent = 0;
480   if(fCollisionType==kPbPb){
481     if(aod)cent = aod->GetHeader()->GetCentrality();
482     if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
483     if(cent<0)cent = 101;
484   }
485   fCentrality = cent;
486   fRPAngle = 0;
487
488   if(esd){
489     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
490     esdVtxValid = IsVertexValid(vtxESD);
491     esdVtxIn = IsVertexIn(vtxESD);
492     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
493       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
494       if(cent<=80&&esdVtxIn){
495         aodH->SetFillAOD(kTRUE);
496         aodH->SetFillExtension(kTRUE);
497       }
498     }
499
500
501     Float_t zvtx = vtxESD->GetZ();
502     Int_t  iCl = GetEventClass(esd);
503     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
504     Bool_t cand = physicsSelection;
505
506     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
507     fh2ESDTriggerCount->Fill(0.,kAllTriggered); 
508     fh2ESDTriggerCount->Fill(iCl,kAllTriggered); 
509     if(cand){
510       fh2ESDTriggerCount->Fill(0.,kSelectedALICE); 
511       fh2ESDTriggerCount->Fill(iCl,kSelectedALICE); 
512       fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
513     }
514     //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
515     if(esdVtxValid){
516       fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
517       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
518       fh2ESDTriggerVtx->Fill(iCl,zvtx);
519       if(esdVtxIn){
520         fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
521         fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
522         fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
523       }
524       if(cand){
525         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
526         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
527         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
528       }
529     }
530
531     if(cand&&esdVtxIn&&iCl<5){
532       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
533       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
534       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
535       fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
536       fh2ESDTriggerCount->Fill(iCl,kSelected);
537       fh2ESDTriggerCount->Fill(0.,kSelected);
538       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
539       if(esd->GetCentrality()){
540         Float_t tmpCent = 100;
541         tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
542         if(tmpCent<0)tmpCent = 101;
543         fh1CentralityESD->Fill(tmpCent);
544       }
545     }
546   }
547
548
549
550   if(aod){
551     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
552     aodVtxValid = IsVertexValid(vtxAOD);
553     aodVtxIn = IsVertexIn(vtxAOD);
554     Float_t zvtx = vtxAOD->GetZ();
555     Int_t  iCl = GetEventClass(aod);
556     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
557     Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
558     if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
559     fh2TriggerCount->Fill(0.,kAllTriggered); 
560     fh2TriggerCount->Fill(iCl,kAllTriggered); 
561     if(cand){
562       fh2TriggerCount->Fill(0.,kSelectedALICE); 
563       fh2TriggerCount->Fill(iCl,kSelectedALICE); 
564       fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
565     }
566     if(aodVtxValid){
567       fh2TriggerCount->Fill(0.,kTriggeredVertex);
568       fh2TriggerCount->Fill(iCl,kTriggeredVertex);
569       fh2TriggerVtx->Fill(iCl,zvtx);
570       if(aodVtxIn){
571         fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
572         fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
573         fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
574       }
575       if(cand){
576         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
577         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
578         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
579       }
580     }
581
582     if(cand&&aodVtxIn&&iCl<5){
583       fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
584       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
585       fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
586       fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
587       fh2TriggerCount->Fill(iCl,kSelected);
588       fh2TriggerCount->Fill(0.,kSelected);
589       fh1Centrality->Fill(cent);
590       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
591       if(aodH&&cand&&fFilterAODCollisions&&!esd){
592         if(fCentrality<=80&&aodVtxIn){
593           aodH->SetFillAOD(kTRUE);
594           aodH->SetFillExtension(kTRUE);
595         }
596       }
597
598       TList recTracks;
599       GetListOfTracks(&recTracks);
600       CalculateReactionPlaneAngle(&recTracks);
601       AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
602       if(fUseAODInput&&fCentrality<=80){
603         if(fFilterAODCollisions&&aod){
604           aodH->SetFillAOD(kTRUE);
605         }
606       }
607     }    
608   }
609
610   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
611
612   
613   Double_t ptHard = 0; 
614   Double_t nTrials = 1; // Trials for MC trigger 
615
616   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
617   AliMCEvent* mcEvent = MCEvent();
618   //    AliStack *pStack = 0; 
619   if(mcEvent){
620     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
621     if(pythiaGenHeader){
622       nTrials = pythiaGenHeader->Trials();
623       ptHard  = pythiaGenHeader->GetPtHard();
624       int iProcessType = pythiaGenHeader->ProcessType();
625       // 11 f+f -> f+f
626       // 12 f+barf -> f+barf
627       // 13 f+barf -> g+g
628       // 28 f+g -> f+g
629       // 53 g+g -> f+barf
630       // 68 g+g -> g+g
631       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
632       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
633       fh1PtHard->Fill(ptHard);
634       fh1PtHardTrials->Fill(ptHard,nTrials);
635
636     }// if pythia gen header
637   }
638
639   // trigger selection
640   
641   // replication of 
642   if(fNonStdFile.Length()&&aod){
643     if (fgAODHeader){
644       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
645       Double_t q[2] = {fRPAngle,fRPAngle};
646       fgAODHeader->SetQTheta(q,2);
647     }
648     if(fgAODVertices){
649       fgAODVertices->Delete();
650       TClonesArray &vertices = *fgAODVertices;
651       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
652       // we only use some basic information, 
653       
654
655       Double_t pos[3];
656       Double_t covVtx[6];
657
658       vtxAOD->GetXYZ(pos); // position                                                                
659       vtxAOD->GetCovMatrix(covVtx); //covariance matrix                                                      
660       Int_t jVertices = 0;
661       AliAODVertex * vtx = new(vertices[jVertices++])
662         AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
663       vtx->SetName(vtxAOD->GetName());
664       vtx->SetTitle(vtxAOD->GetTitle());
665       TString vtitle = vtxAOD->GetTitle();
666       vtx->SetNContributors(vtxAOD->GetNContributors());
667
668       // Add SPD "main" vertex                                                                    
669       const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
670       vtxS->GetXYZ(pos); // position
671       vtxS->GetCovMatrix(covVtx); //covariance matrix                                                      
672       AliAODVertex * mVSPD = new(vertices[jVertices++])
673         AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
674       mVSPD->SetName(vtxS->GetName());
675       mVSPD->SetTitle(vtxS->GetTitle());
676       mVSPD->SetNContributors(vtxS->GetNContributors());
677       
678       // Add tpc only vertex
679       if(esd){
680         const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
681         vtxT->GetXYZ(pos); // position                                                    
682         vtxT->GetCovMatrix(covVtx); //covariance matrix                                               
683         AliAODVertex * mVTPC = new(vertices[jVertices++])
684           AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
685         mVTPC->SetName(vtxT->GetName());
686         mVTPC->SetTitle(vtxT->GetTitle());
687         mVTPC->SetNContributors(vtxT->GetNContributors());
688       }
689     }
690   }
691   
692   PostData(1, fHistList);
693 }
694
695 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
696   if(!esd)return kFALSE;
697   const AliESDVertex *vtx = esd->GetPrimaryVertex();
698   return IsVertexIn(vtx); // vertex in calls vertex valid
699 }
700
701 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
702   if(fgAODVertices){
703     fgAODVertices->Delete();
704     delete fgAODVertices;
705   }
706   delete fRandomizer;
707   if(fgAODHeader)delete fgAODHeader;
708 }
709
710
711 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
712   if(!aod)return kFALSE;
713   const AliAODVertex *vtx = aod->GetPrimaryVertex();
714   return IsVertexIn(vtx); // VertexIn calls VertexValid
715 }
716
717 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
718
719   // We can check the type of the vertex by its name and title
720   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
721   // vtx       name                  title
722   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
723   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
724   // TPC       TPCVertex             VertexerTracksNoConstraint
725   // TPC       TPCVertex             VertexerTracksWithConstraint
726   // SPD       SPDVertex             vertexer: 3D
727   // SPD       SPDVertex             vertexer: Z
728   
729   Int_t nCont = vtx->GetNContributors();
730   if(nCont>=1){
731     fEventCutInfoESD |= kContributorsCut1;    
732     if(nCont>=2){
733       fEventCutInfoESD |= kContributorsCut2;    
734       if(nCont>=3){
735         fEventCutInfoESD |= kContributorsCut3;    
736       }
737     }
738   }
739   
740   if(nCont<3)return kFALSE;
741
742   // do not want tpc only primary vertex
743   TString vtxName(vtx->GetName());
744   if(vtxName.Contains("TPCVertex")){
745     fEventCutInfoESD |= kVertexTPC;
746     return kFALSE;
747   }
748   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
749   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
750
751
752   TString vtxTitle(vtx->GetTitle());
753   if(vtxTitle.Contains("vertexer: Z")){
754     if(vtx->GetDispersion()>0.02)return kFALSE;   
755   }
756   fEventCutInfoESD |= kSPDDispersionCut;
757   return kTRUE;
758 }
759
760
761 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
762
763   // We can check the type of the vertex by its name and title
764   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
765   // vtx       name                  title
766   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
767   // TPC       TPCVertex             VertexerTracksNoConstraint
768   // SPD       SPDVertex             vertexer: 3D
769   // SPD       SPDVertex             vertexer: Z
770
771   if(fDebug){
772     Printf(" n contrib %d",vtx->GetNContributors());
773     vtx->Print();
774   }
775   
776   //  if(vtx->GetNContributors()<3)return kFALSE;
777   // do not want tpc only primary vertex
778   TString vtxName(vtx->GetName());
779   if(vtxName.Contains("TPCVertex"))return kFALSE;
780
781   // no dispersion yet...
782   /* 
783   TString vtxTitle(vtx->GetTitle());
784   if(vtxTitle.Contains("vertexer: Z")){
785     if(vtx->GetDispersion()>0.02)return kFALSE;
786   }
787   */
788   return kTRUE;
789 }
790
791
792 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
793
794   if(!IsVertexValid(vtx))return kFALSE;
795
796   Float_t zvtx = vtx->GetZ();
797   Float_t yvtx = vtx->GetY();
798   Float_t xvtx = vtx->GetX();
799
800   xvtx -= fVtxXMean;
801   yvtx -= fVtxYMean;
802   zvtx -= fVtxZMean;
803
804
805
806   if(TMath::Abs(zvtx)>fVtxZCut){
807     return kFALSE;
808   }
809   fEventCutInfoESD |= kVertexZCut;  
810   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
811   if(r2>(fVtxRCut*fVtxRCut)){
812     return kFALSE;
813   }
814   fEventCutInfoESD |= kVertexRCut;  
815   return kTRUE;
816 }
817
818
819 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
820
821   if(!IsVertexValid(vtx))return kFALSE;
822
823   Float_t zvtx = vtx->GetZ();
824   Float_t yvtx = vtx->GetY();
825   Float_t xvtx = vtx->GetX();
826
827   xvtx -= fVtxXMean;
828   yvtx -= fVtxYMean;
829   zvtx -= fVtxZMean;
830
831   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
832
833   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
834   return vertexIn;
835 }
836
837 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
838   if(!esd)return kFALSE;
839   return esd->IsPileupFromSPD();
840 }
841
842 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
843   if(!esd)return kFALSE;
844   // add track cuts for which we look for cosmics...
845
846   Bool_t isCosmic = kFALSE;
847   Int_t nTracks = esd->GetNumberOfTracks();
848   Int_t nCosmicCandidates = 0;
849
850   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
851     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
852     if (!track1)  continue;
853     UInt_t status1 = track1->GetStatus();
854     //If track is ITS stand alone track, skip the track
855     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
856     if(track1->Pt()<fPtMinCosmic) continue;
857     //Start 2nd track loop to look for correlations
858     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
859       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
860       if(!track2) continue;
861       UInt_t status2 = track2->GetStatus();
862       //If track is ITS stand alone track, skip the track
863       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
864       if(track2->Pt()<fPtMinCosmic) continue;
865       //Check if back-to-back
866       Double_t mom1[3],mom2[3];
867       track1->GetPxPyPz(mom1);
868       track2->GetPxPyPz(mom2);
869       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
870       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
871       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
872       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
873
874       Float_t deltaPhi = track1->Phi()-track2->Phi();
875       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
876
877       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
878       if(rIsol<fRIsolMinCosmic) continue;
879
880       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
881         nCosmicCandidates+=1;
882         isCosmic = kTRUE;
883       }
884       
885     }
886   }
887
888   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
889
890   return isCosmic;
891 }
892
893
894 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
895
896   Float_t cent = 0;
897   if(fCollisionType==kPbPb){
898     if(esd->GetCentrality()){
899       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
900     }
901     if(cent<0)cent = 101;
902     if(cent>80||cent<0)return 5;
903     if(cent>50)return 4;
904     if(cent>30)return 3;
905     if(cent>10)return 2;
906     return 1;
907   }
908   return 1;
909 }
910
911
912 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
913
914   if(fCollisionType==kPbPb){
915     Float_t cent = aod->GetHeader()->GetCentrality();
916     if(cent>80||cent<0)return 5;
917     if(cent>50)return 4;
918     if(cent>30)return 3;
919     if(cent>10)return 2;
920     return 1;
921   }
922   return 1;
923
924 }
925
926
927 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
928 {
929   // Terminate analysis
930   //
931 }
932
933 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngle(const TList *trackList)
934 {
935
936   if(!trackList)return kFALSE;
937   fRPAngle=0;
938
939   // need to get this info from elsewhere??
940
941   Double_t fPsiRP =0,fDeltaPsiRP = 0;
942    
943    
944     
945   TVector2 mQ,mQ1,mQ2;
946   Float_t mQx= fDeltaQxy[0], mQy=fDeltaQxy[1];
947   
948   Float_t mQx1=fDeltaQxy[0], mQy1=fDeltaQxy[1];
949   Float_t mQx2=fDeltaQxy[0], mQy2=fDeltaQxy[1];
950   
951   AliVParticle *track=0x0;
952   Int_t count[3]={0,0,0};
953   
954
955   for (Int_t iter=0;iter<trackList->GetEntries();iter++){
956
957     track=(AliVParticle*)trackList->At(iter);
958     
959     //cuts already applied before
960     // Comment DCA not correctly implemented yet for AOD tracks
961     
962     Double_t momentum;
963     if(track->Charge()>0){momentum=track->Pt();}
964     else{momentum=-track->Pt();}
965
966        
967
968     // For Weighting
969     fh3RPPhiTracks->Fill(fCentrality,momentum,track->Phi());
970     count[0]++;
971
972     Double_t phiweight=GetPhiWeight(track->Phi(),momentum);
973     //    Double_t phiweight=1; 
974     Double_t weight=2;
975     if(track->Pt()<2){weight=track->Pt();}
976     
977
978     mQx += (cos(2*track->Phi()))*weight*phiweight;
979     mQy += (sin(2*track->Phi()))*weight*phiweight;
980
981     // Make random Subevents
982
983     if(fRPSubeventMethod==0){
984       if(fRandomizer->Binomial(1,0.5)){
985         mQx1 += (cos(2*track->Phi()))*weight*phiweight;
986         mQy1 += (sin(2*track->Phi()))*weight*phiweight;
987         count[1]++;}
988       else{
989         mQx2 += (cos(2*track->Phi()))*weight*phiweight;
990         mQy2 += (sin(2*track->Phi()))*weight*phiweight;
991         count[2]++;}
992     }
993     else if(fRPSubeventMethod==1){
994       // Make eta dependent subevents
995       if(track->Eta()>0){
996         mQx1 += (cos(2*track->Phi()))*weight*phiweight;
997         mQy1 += (sin(2*track->Phi()))*weight*phiweight;
998         count[1]++;}
999       else{
1000         mQx2 += (cos(2*track->Phi()))*weight*phiweight;
1001         mQy2 += (sin(2*track->Phi()))*weight*phiweight;
1002         count[2]++;}
1003     }
1004
1005   }
1006
1007
1008
1009   //If no track passes the cuts, the ,Q.Phi() will return Pi and a peak at Pi/2 in the RP Angular Distribution will appear
1010   if(count[0]==0||count[1]==0||count[2]==0){
1011     return kFALSE;
1012   }
1013
1014   mQ.Set(mQx,mQy);
1015   mQ1.Set(mQx1,mQy1);
1016   mQ2.Set(mQx2,mQy2);
1017
1018   // cout<<"MQ"<<mQx<<" " <<mQy<<" psi"<<endl;
1019
1020   fPsiRP=mQ.Phi()/2;
1021     
1022   //Correction
1023   fPsiRP+=fFlatA[0]*TMath::Cos(2*fPsiRP)+fFlatB[0]*TMath::Sin(2*fPsiRP)+fFlatA[1]*TMath::Cos(4*fPsiRP)+fFlatB[1]*TMath::Sin(4*fPsiRP);
1024
1025   Double_t fPsiRP1=mQ1.Phi()/2;
1026   fPsiRP1+=fFlatA[0]*TMath::Cos(2*fPsiRP1)+fFlatB[0]*TMath::Sin(2*fPsiRP1)+fFlatA[1]*TMath::Cos(4*fPsiRP1)+fFlatB[1]*TMath::Sin(4*fPsiRP1);
1027   Double_t fPsiRP2=mQ2.Phi()/2;
1028   fPsiRP2+=fFlatA[0]*TMath::Cos(2*fPsiRP2)+fFlatB[0]*TMath::Sin(2*fPsiRP2)+fFlatA[1]*TMath::Cos(4*fPsiRP2)+fFlatB[1]*TMath::Sin(4*fPsiRP2);
1029   fDeltaPsiRP=fPsiRP1-fPsiRP2;
1030   
1031   if(fPsiRP>TMath::Pi()){fPsiRP-=TMath::Pi();}
1032   if(fPsiRP<0){fPsiRP+=TMath::Pi();}
1033   
1034   // reactionplaneangle + Pi() is the same angle
1035   if(TMath::Abs(fDeltaPsiRP)>TMath::Pi()/2){
1036     if(fDeltaPsiRP>0)fDeltaPsiRP-=TMath::Pi();
1037     else fDeltaPsiRP+=TMath::Pi();
1038   }
1039   
1040   Double_t cos2deltaRP=TMath::Cos(2*fDeltaPsiRP);
1041   
1042   // FillHistograms
1043   fh2RPSubevents->Fill(fPsiRP1,fPsiRP2);
1044   fh1RP->Fill(fPsiRP);
1045   fh2RPCentrality->Fill(fCentrality,fPsiRP);
1046   fh2RPDeltaRP->Fill(fDeltaPsiRP,fCentrality);
1047   fh2RPQxQy->Fill(mQx,mQy);
1048   fh2RPCosDeltaRP->Fill(fCentrality,cos2deltaRP);
1049   
1050   fRPAngle=fPsiRP;  
1051   return kTRUE;
1052 }
1053
1054 Double_t AliAnalysisTaskJetServices::GetPhiWeight(Double_t phi,Double_t signedpt){
1055   if(!fh3PhiWeights)return 1;
1056   else return fh3PhiWeights->GetBinContent(fh3PhiWeights->GetXaxis()->FindBin(fCentrality),fh3PhiWeights->GetYaxis()->FindBin(signedpt),fh3PhiWeights->GetZaxis()->FindBin(phi));
1057 }
1058
1059  //________________________________________________________________________
1060
1061 Int_t  AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1062   Int_t iCount = 0;
1063   AliAODEvent *aod = 0;
1064   if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1065   else aod = AODEvent();
1066   if(!aod){
1067     return iCount;
1068   }
1069   for(int it = 0;it < aod->GetNumberOfTracks();++it){
1070     AliAODTrack *tr = aod->GetTrack(it);
1071     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1072     if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1073     if(tr->Pt()<fMinTrackPt)continue;
1074     list->Add(tr);
1075     iCount++;
1076   }
1077   list->Sort();
1078   return iCount;
1079
1080 }
1081