f17cb81dca6b54e9a720743e65dde98ea23f6044
[u/mrichter/AliRoot.git] / JETAN / DEV / AliAnalysisTaskJetCluster.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 <TRandom3.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include <TChain.h>
28 #include <TRefArray.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 <TClonesArray.h>
38 #include  "TDatabasePDG.h"
39
40 #include "AliAnalysisTaskJetCluster.h"
41 #include "AliAnalysisManager.h"
42 #include "AliJetFinder.h"
43 #include "AliJetHeader.h"
44 #include "AliJetReader.h"
45 #include "AliESDEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliAODHandler.h"
48 #include "AliAODExtension.h"
49 #include "AliAODTrack.h"
50 #include "AliAODJet.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
54 #include "AliStack.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
59 #include "AliAODJetEventBackground.h"
60
61 #include "fastjet/PseudoJet.hh"
62 #include "fastjet/ClusterSequenceArea.hh"
63 #include "fastjet/AreaDefinition.hh"
64 #include "fastjet/JetDefinition.hh"
65 // get info on how fastjet was configured
66 #include "fastjet/config.h"
67
68
69 ClassImp(AliAnalysisTaskJetCluster)
70
71 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
72   //
73   // Destructor
74   //
75   delete fRef;
76   delete fRandom;
77
78   if(fTCAJetsOut)fTCAJetsOut->Delete();
79   delete fTCAJetsOut;
80   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
81   delete fTCAJetsOutRan;
82   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
83   delete fTCARandomConesOut;
84   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
85   delete fTCARandomConesOutRan;
86   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
87   delete fAODJetBackgroundOut;
88 }
89
90 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): 
91   AliAnalysisTaskSE(),
92   fAOD(0x0),
93   fAODExtension(0x0),
94   fRef(new TRefArray),
95   fUseAODTrackInput(kFALSE),
96   fUseAODMCInput(kFALSE),
97   fUseBackgroundCalc(kFALSE),
98   fEventSelection(kFALSE),     
99   fFilterMask(0),
100   fFilterType(0),
101   fJetTypes(kJet),
102   fTrackTypeRec(kTrackUndef),
103   fTrackTypeGen(kTrackUndef),  
104   fNSkipLeadingRan(0),
105   fNSkipLeadingCone(0),
106   fNRandomCones(0),
107   fAvgTrials(1),
108   fExternalWeight(1),
109   fTrackEtaWindow(0.9),    
110   fRecEtaWindow(0.5),
111   fTrackPtCut(0.),                                                      
112   fJetOutputMinPt(0.150),
113   fMaxTrackPtInJet(100.),
114   fJetTriggerPtCut(0),
115   fVtxZCut(8),
116   fVtxR2Cut(1),
117   fCentCutUp(0),
118   fCentCutLo(0),
119   fNonStdBranch(""),
120   fBackgroundBranch(""),
121   fNonStdFile(""),
122   fRparam(1.0), 
123   fAlgorithm(fastjet::kt_algorithm),
124   fStrategy(fastjet::Best),
125   fRecombScheme(fastjet::BIpt_scheme),
126   fAreaType(fastjet::active_area), 
127   fGhostArea(0.01),
128   fActiveAreaRepeats(1),
129   fGhostEtamax(1.5),
130   fTCAJetsOut(0x0),
131   fTCAJetsOutRan(0x0),
132   fTCARandomConesOut(0x0),
133   fTCARandomConesOutRan(0x0),
134   fAODJetBackgroundOut(0x0),
135   fRandom(0),
136   fh1Xsec(0x0),
137   fh1Trials(0x0),
138   fh1PtHard(0x0),
139   fh1PtHardNoW(0x0),  
140   fh1PtHardTrials(0x0),
141   fh1NJetsRec(0x0),
142   fh1NConstRec(0x0),
143   fh1NConstLeadingRec(0x0),
144   fh1PtJetsRecIn(0x0),
145   fh1PtJetsLeadingRecIn(0x0),
146   fh1PtJetConstRec(0x0),
147   fh1PtJetConstLeadingRec(0x0),
148   fh1PtTracksRecIn(0x0),
149   fh1PtTracksLeadingRecIn(0x0),
150   fh1NJetsRecRan(0x0),
151   fh1NConstRecRan(0x0),
152   fh1PtJetsLeadingRecInRan(0x0),
153   fh1NConstLeadingRecRan(0x0),
154   fh1PtJetsRecInRan(0x0),
155   fh1PtTracksGenIn(0x0),
156   fh1Nch(0x0),
157   fh1CentralityPhySel(0x0), 
158   fh1Centrality(0x0), 
159   fh1CentralitySelect(0x0),
160   fh1ZPhySel(0x0), 
161   fh1Z(0x0), 
162   fh1ZSelect(0x0),
163   fh2NRecJetsPt(0x0),
164   fh2NRecTracksPt(0x0),
165   fh2NConstPt(0x0),
166   fh2NConstLeadingPt(0x0),
167   fh2JetPhiEta(0x0),
168   fh2LeadingJetPhiEta(0x0),
169   fh2JetEtaPt(0x0),
170   fh2LeadingJetEtaPt(0x0),
171   fh2TrackEtaPt(0x0),
172   fh2LeadingTrackEtaPt(0x0),
173   fh2JetsLeadingPhiEta(0x0),
174   fh2JetsLeadingPhiPt(0x0),
175   fh2TracksLeadingPhiEta(0x0),
176   fh2TracksLeadingPhiPt(0x0),
177   fh2TracksLeadingJetPhiPt(0x0),
178   fh2JetsLeadingPhiPtW(0x0),
179   fh2TracksLeadingPhiPtW(0x0),
180   fh2TracksLeadingJetPhiPtW(0x0),
181   fh2NRecJetsPtRan(0x0),
182   fh2NConstPtRan(0x0),
183   fh2NConstLeadingPtRan(0x0),
184   fh2PtNch(0x0),
185   fh2PtNchRan(0x0),
186   fh2PtNchN(0x0),
187   fh2PtNchNRan(0x0),
188   fh2TracksLeadingJetPhiPtRan(0x0),
189   fh2TracksLeadingJetPhiPtWRan(0x0),
190   fHistList(0x0)  
191 {
192   //
193   // Constructor
194   //
195
196   for(int i = 0;i<3;i++){
197     fh1BiARandomCones[i] = 0;
198     fh1BiARandomConesRan[i] = 0;
199   }
200   for(int i = 0;i<kMaxCent;i++){
201     fh2JetsLeadingPhiPtC[i] = 0;     
202     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
203     fh2TracksLeadingJetPhiPtC[i] = 0;
204     fh2TracksLeadingJetPhiPtWC[i] = 0;
205   }
206 }
207
208 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
209   AliAnalysisTaskSE(name),
210   fAOD(0x0),
211   fAODExtension(0x0),
212   fRef(new TRefArray),
213   fUseAODTrackInput(kFALSE),
214   fUseAODMCInput(kFALSE),
215   fUseBackgroundCalc(kFALSE),
216   fEventSelection(kFALSE),                                                      
217   fFilterMask(0),
218   fFilterType(0),
219   fJetTypes(kJet),
220   fTrackTypeRec(kTrackUndef),
221   fTrackTypeGen(kTrackUndef),
222   fNSkipLeadingRan(0),
223   fNSkipLeadingCone(0),
224   fNRandomCones(0),
225   fAvgTrials(1),
226   fExternalWeight(1),    
227   fTrackEtaWindow(0.9),    
228   fRecEtaWindow(0.5),
229   fTrackPtCut(0.),                                                      
230   fJetOutputMinPt(0.150),
231   fMaxTrackPtInJet(100.),
232   fJetTriggerPtCut(0),
233   fVtxZCut(8),
234   fVtxR2Cut(1),
235   fCentCutUp(0),
236   fCentCutLo(0),
237   fNonStdBranch(""),
238   fBackgroundBranch(""),
239   fNonStdFile(""),
240   fRparam(1.0), 
241   fAlgorithm(fastjet::kt_algorithm),
242   fStrategy(fastjet::Best),
243   fRecombScheme(fastjet::BIpt_scheme),
244   fAreaType(fastjet::active_area), 
245   fGhostArea(0.01),
246   fActiveAreaRepeats(1),
247   fGhostEtamax(1.5),
248   fTCAJetsOut(0x0),
249   fTCAJetsOutRan(0x0),
250   fTCARandomConesOut(0x0),
251   fTCARandomConesOutRan(0x0),
252   fAODJetBackgroundOut(0x0),
253   fRandom(0),
254   fh1Xsec(0x0),
255   fh1Trials(0x0),
256   fh1PtHard(0x0),
257   fh1PtHardNoW(0x0),  
258   fh1PtHardTrials(0x0),
259   fh1NJetsRec(0x0),
260   fh1NConstRec(0x0),
261   fh1NConstLeadingRec(0x0),
262   fh1PtJetsRecIn(0x0),
263   fh1PtJetsLeadingRecIn(0x0),
264   fh1PtJetConstRec(0x0),
265   fh1PtJetConstLeadingRec(0x0),
266   fh1PtTracksRecIn(0x0),
267   fh1PtTracksLeadingRecIn(0x0),
268   fh1NJetsRecRan(0x0),
269   fh1NConstRecRan(0x0),
270   fh1PtJetsLeadingRecInRan(0x0),
271   fh1NConstLeadingRecRan(0x0),
272   fh1PtJetsRecInRan(0x0),
273   fh1PtTracksGenIn(0x0),
274   fh1Nch(0x0),
275   fh1CentralityPhySel(0x0), 
276   fh1Centrality(0x0), 
277   fh1CentralitySelect(0x0),
278   fh1ZPhySel(0x0), 
279   fh1Z(0x0), 
280   fh1ZSelect(0x0),
281   fh2NRecJetsPt(0x0),
282   fh2NRecTracksPt(0x0),
283   fh2NConstPt(0x0),
284   fh2NConstLeadingPt(0x0),
285   fh2JetPhiEta(0x0),
286   fh2LeadingJetPhiEta(0x0),
287   fh2JetEtaPt(0x0),
288   fh2LeadingJetEtaPt(0x0),
289   fh2TrackEtaPt(0x0),
290   fh2LeadingTrackEtaPt(0x0),
291   fh2JetsLeadingPhiEta(0x0),
292   fh2JetsLeadingPhiPt(0x0),
293   fh2TracksLeadingPhiEta(0x0),
294   fh2TracksLeadingPhiPt(0x0),
295   fh2TracksLeadingJetPhiPt(0x0),
296   fh2JetsLeadingPhiPtW(0x0),
297   fh2TracksLeadingPhiPtW(0x0),
298   fh2TracksLeadingJetPhiPtW(0x0),
299   fh2NRecJetsPtRan(0x0),
300   fh2NConstPtRan(0x0),
301   fh2NConstLeadingPtRan(0x0),
302   fh2PtNch(0x0),
303   fh2PtNchRan(0x0),
304   fh2PtNchN(0x0),
305   fh2PtNchNRan(0x0),
306   fh2TracksLeadingJetPhiPtRan(0x0),
307   fh2TracksLeadingJetPhiPtWRan(0x0),
308   fHistList(0x0)
309 {
310   //
311   // named ctor
312   //
313   for(int i = 0;i<3;i++){
314     fh1BiARandomCones[i] = 0;
315     fh1BiARandomConesRan[i] = 0;
316   }
317   for(int i = 0;i<kMaxCent;i++){
318     fh2JetsLeadingPhiPtC[i] = 0;     
319     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
320     fh2TracksLeadingJetPhiPtC[i] = 0;
321     fh2TracksLeadingJetPhiPtWC[i] = 0;
322   }
323  DefineOutput(1, TList::Class());  
324 }
325
326
327
328 Bool_t AliAnalysisTaskJetCluster::Notify()
329 {
330   //
331   // Implemented Notify() to read the cross sections
332   // and number of trials from pyxsec.root
333   // 
334   return kTRUE;
335 }
336
337 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
338 {
339
340   //
341   // Create the output container
342   //
343
344   fRandom = new TRandom3(0);
345
346
347   // Connect the AOD
348
349
350   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
351
352   
353
354   if(fNonStdBranch.Length()!=0)
355     {
356       // only create the output branch if we have a name
357       // Create a new branch for jets...
358       //  -> cleared in the UserExec....
359       // here we can also have the case that the brnaches are written to a separate file
360       
361       if(fJetTypes&kJet){
362         fTCAJetsOut = new TClonesArray("AliAODJet", 0);
363         fTCAJetsOut->SetName(fNonStdBranch.Data());
364         AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
365       }
366    
367       if(fJetTypes&kJetRan){
368         fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
369         fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
370         AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
371       }
372
373       if(fUseBackgroundCalc){
374         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
375           fAODJetBackgroundOut = new AliAODJetEventBackground();
376           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
377           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());  
378         }
379       }
380       // create the branch for the random cones with the same R 
381       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
382   
383       if(fNRandomCones>0){
384         if(fJetTypes&kRC){
385           if(!AODEvent()->FindListObject(cName.Data())){
386             fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
387             fTCARandomConesOut->SetName(cName.Data());
388             AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
389           }
390         }
391         // create the branch with the random for the random cones on the random event
392         if(fJetTypes&kRCRan){
393           cName = Form("%sRandomCone_random",fNonStdBranch.Data());
394           if(!AODEvent()->FindListObject(cName.Data())){
395             fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
396             fTCARandomConesOutRan->SetName(cName.Data());
397             AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
398           }
399         }
400       }
401     
402       if(fNonStdFile.Length()!=0){
403         // 
404         // case that we have an AOD extension we need to fetch the jets from the extended output
405         // we identify the extension aod event by looking for the branchname
406         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
407         // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
408         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
409       }
410     }
411
412
413   if(!fHistList)fHistList = new TList();
414   fHistList->SetOwner();
415   PostData(1, fHistList); // post data in any case once
416
417   Bool_t oldStatus = TH1::AddDirectoryStatus();
418   TH1::AddDirectory(kFALSE);
419
420   //
421   //  Histogram
422     
423   const Int_t nBinPt = 100;
424   Double_t binLimitsPt[nBinPt+1];
425   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
426     if(iPt == 0){
427       binLimitsPt[iPt] = 0.0;
428     }
429     else {// 1.0
430       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
431     }
432   }
433   
434   const Int_t nBinPhi = 90;
435   Double_t binLimitsPhi[nBinPhi+1];
436   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
437     if(iPhi==0){
438       binLimitsPhi[iPhi] = -1.*TMath::Pi();
439     }
440     else{
441       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
442     }
443   }
444
445
446
447   const Int_t nBinEta = 40;
448   Double_t binLimitsEta[nBinEta+1];
449   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
450     if(iEta==0){
451       binLimitsEta[iEta] = -2.0;
452     }
453     else{
454       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
455     }
456   }
457
458   const int nChMax = 5000;
459
460   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
461   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
462
463   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
464   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
465
466
467   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
468   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
469
470   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
471   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
472   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
473   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
474
475
476   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
477   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
478   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
479
480   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
481   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
482   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
483   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
484   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
485   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
486   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
487   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
488   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
489   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
490
491   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
492   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
493   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
494
495   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
496   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
497   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
498
499   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
500   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
501   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
502   // 
503
504
505   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
506   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
507   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
508   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
509
510   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
511   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
512   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
513   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
514
515
516
517   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
518                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
519   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
520                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
521
522   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
523                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
524   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
525                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
526
527   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
528                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
529   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
530                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
531
532
533
534   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
535                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
536   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
537                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
538   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
539                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
540   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
541                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
542   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
543                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
544   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
545                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
546
547   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
548                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
549   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
550                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
551
552   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
553                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
554   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
555                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
556
557
558   if(fNRandomCones>0&&fUseBackgroundCalc){
559     for(int i = 0;i<3;i++){
560       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
561       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
562     }
563   }
564
565   for(int i = 0;i < kMaxCent;i++){
566     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
567     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
568     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
569     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
570   }
571
572   const Int_t saveLevel = 3; // large save level more histos
573   if(saveLevel>0){
574     fHistList->Add(fh1Xsec);
575     fHistList->Add(fh1Trials);
576
577     fHistList->Add(fh1NJetsRec);
578     fHistList->Add(fh1NConstRec);
579     fHistList->Add(fh1NConstLeadingRec);
580     fHistList->Add(fh1PtJetsRecIn);
581     fHistList->Add(fh1PtJetsLeadingRecIn);
582     fHistList->Add(fh1PtTracksRecIn);
583     fHistList->Add(fh1PtTracksLeadingRecIn);
584     fHistList->Add(fh1PtJetConstRec);
585     fHistList->Add(fh1PtJetConstLeadingRec);
586     fHistList->Add(fh1NJetsRecRan);
587     fHistList->Add(fh1NConstRecRan);
588     fHistList->Add(fh1PtJetsLeadingRecInRan);
589     fHistList->Add(fh1NConstLeadingRecRan);
590     fHistList->Add(fh1PtJetsRecInRan);
591     fHistList->Add(fh1Nch);
592     fHistList->Add(fh1Centrality);
593     fHistList->Add(fh1CentralitySelect);
594     fHistList->Add(fh1CentralityPhySel);
595     fHistList->Add(fh1Z);
596     fHistList->Add(fh1ZSelect);
597     fHistList->Add(fh1ZPhySel);
598     if(fNRandomCones>0&&fUseBackgroundCalc){
599       for(int i = 0;i<3;i++){
600         fHistList->Add(fh1BiARandomCones[i]);
601         fHistList->Add(fh1BiARandomConesRan[i]);
602       }
603     }
604   for(int i = 0;i < kMaxCent;i++){
605     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
606     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
607     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
608     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
609   }
610
611     fHistList->Add(fh2NRecJetsPt);
612     fHistList->Add(fh2NRecTracksPt);
613     fHistList->Add(fh2NConstPt);
614     fHistList->Add(fh2NConstLeadingPt);
615     fHistList->Add(fh2PtNch);
616     fHistList->Add(fh2PtNchRan);
617     fHistList->Add(fh2PtNchN);
618     fHistList->Add(fh2PtNchNRan);
619     fHistList->Add(fh2JetPhiEta);
620     fHistList->Add(fh2LeadingJetPhiEta);
621     fHistList->Add(fh2JetEtaPt);
622     fHistList->Add(fh2LeadingJetEtaPt);
623     fHistList->Add(fh2TrackEtaPt);
624     fHistList->Add(fh2LeadingTrackEtaPt);
625     fHistList->Add(fh2JetsLeadingPhiEta );
626     fHistList->Add(fh2JetsLeadingPhiPt);
627     fHistList->Add(fh2TracksLeadingPhiEta);
628     fHistList->Add(fh2TracksLeadingPhiPt);
629     fHistList->Add(fh2TracksLeadingJetPhiPt);
630     fHistList->Add(fh2JetsLeadingPhiPtW);
631     fHistList->Add(fh2TracksLeadingPhiPtW);
632     fHistList->Add(fh2TracksLeadingJetPhiPtW);
633     fHistList->Add(fh2NRecJetsPtRan);
634     fHistList->Add(fh2NConstPtRan);
635     fHistList->Add(fh2NConstLeadingPtRan);
636     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
637     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
638   }
639
640   // =========== Switch on Sumw2 for all histos ===========
641   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
642     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
643     if (h1){
644       h1->Sumw2();
645       continue;
646     }
647     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
648     if(hn)hn->Sumw2();
649   }
650   TH1::AddDirectory(oldStatus);
651 }
652
653 void AliAnalysisTaskJetCluster::Init()
654 {
655   //
656   // Initialization
657   //
658
659   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
660
661 }
662
663 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
664 {
665
666   // handle and reset the output jet branch 
667
668   if(fTCAJetsOut)fTCAJetsOut->Delete();
669   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
670   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
671   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
672   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
673
674   AliAODJetEventBackground* externalBackground = 0;
675   if(!externalBackground&&fBackgroundBranch.Length()){
676     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
677     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
678     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
679   }
680   //
681   // Execute analysis for current event
682   //
683   AliESDEvent *fESD = 0;
684   if(fUseAODTrackInput){    
685     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
686     if(!fAOD){
687       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
688       return;
689     }
690     // fethc the header
691   }
692   else{
693     //  assume that the AOD is in the general output...
694     fAOD  = AODEvent();
695     if(!fAOD){
696       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
697       return;
698     }
699     if(fDebug>0){
700       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
701     }
702   }
703   
704   Bool_t selectEvent =  false;
705   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
706
707   Float_t cent = 0;
708   Float_t zVtx  = 0;
709   Int_t cenClass = -1;
710   if(fAOD){
711     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
712     TString vtxTitle(vtxAOD->GetTitle());
713     zVtx = vtxAOD->GetZ();
714
715     cent = fAOD->GetHeader()->GetCentrality();
716     if(cent<10)cenClass = 0;
717     else if(cent<30)cenClass = 1;
718     else if(cent<50)cenClass = 2;
719     else if(cent<80)cenClass = 3;
720     if(physicsSelection){
721       fh1CentralityPhySel->Fill(cent);
722       fh1ZPhySel->Fill(zVtx);
723     }
724
725     if(fEventSelection){
726       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
727         Float_t yvtx = vtxAOD->GetY();
728         Float_t xvtx = vtxAOD->GetX();
729         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
730         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
731           if(physicsSelection){
732             selectEvent = true;
733           }
734         }
735       }
736       if(fCentCutUp>0){
737         if(cent<fCentCutLo||cent>fCentCutUp){
738           selectEvent = false;
739         }
740       }
741     }else{
742       selectEvent = true;
743     }
744   }
745   
746
747   if(!selectEvent){
748     PostData(1, fHistList);
749     return;
750   }
751   fh1Centrality->Fill(cent);  
752   fh1Z->Fill(zVtx);
753   fh1Trials->Fill("#sum{ntrials}",1);
754   
755
756   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
757
758   // ==== General variables needed
759
760
761
762   // we simply fetch the tracks/mc particles as a list of AliVParticles
763
764   TList recParticles;
765   TList genParticles;
766
767   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
768   Float_t nCh = recParticles.GetEntries(); 
769   fh1Nch->Fill(nCh);
770   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
771   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
772   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
773
774   // find the jets....
775
776   vector<fastjet::PseudoJet> inputParticlesRec;
777   vector<fastjet::PseudoJet> inputParticlesRecRan;
778   
779   // Generate the random cones
780   
781   AliAODJet vTmpRan(1,0,0,1);
782   for(int i = 0; i < recParticles.GetEntries(); i++){
783     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
784     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
785     // we take total momentum here
786     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
787     jInp.set_user_index(i);
788     inputParticlesRec.push_back(jInp);
789
790     // the randomized input changes eta and phi, but keeps the p_T
791     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
792       Double_t pT = vp->Pt();
793       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
794       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
795       
796       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
797       Double_t pZ = pT/TMath::Tan(theta);
798
799       Double_t pX = pT * TMath::Cos(phi);
800       Double_t pY = pT * TMath::Sin(phi);
801       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
802       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
803
804       jInpRan.set_user_index(i);
805       inputParticlesRecRan.push_back(jInpRan);
806       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
807     }
808
809     // fill the tref array, only needed when we write out jets
810     if(fTCAJetsOut){
811       if(i == 0){
812         fRef->Delete(); // make sure to delete before placement new...
813         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
814       }
815       fRef->Add(vp);
816     }
817   }// recparticles
818
819   if(inputParticlesRec.size()==0){
820     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
821     PostData(1, fHistList);
822     return;
823   }
824   
825   // run fast jet
826   // employ setters for these...
827
828  
829   // now create the object that holds info about ghosts                        
830   /*
831   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
832     // reduce CPU time...
833     fGhostArea = 0.5; 
834     fActiveAreaRepeats = 0; 
835   }
836   */
837
838   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
839   fastjet::AreaType areaType =   fastjet::active_area;
840   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
841   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
842   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
843   
844   //range where to compute background
845   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
846   phiMin = 0;
847   phiMax = 2*TMath::Pi();
848   rapMax = fGhostEtamax - fRparam;
849   rapMin = - fGhostEtamax + fRparam;
850   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
851  
852
853   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
854   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
855
856  
857   fh1NJetsRec->Fill(sortedJets.size());
858
859  // loop over all jets an fill information, first on is the leading jet
860
861   Int_t nRecOver = inclusiveJets.size();
862   Int_t nRec     = inclusiveJets.size();
863   if(inclusiveJets.size()>0){
864     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
865     Double_t area = clustSeq.area(sortedJets[0]);
866     leadingJet.SetEffArea(area,0);
867     Float_t pt = leadingJet.Pt();
868     Int_t nAodOutJets = 0;
869     Int_t nAodOutTracks = 0;
870     AliAODJet *aodOutJet = 0;
871
872     Int_t iCount = 0;
873     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
874       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
875       while(pt<ptCut&&iCount<nRec){
876         nRecOver--;
877         iCount++;
878         if(iCount<nRec){
879           pt = sortedJets[iCount].perp();
880         }
881       }
882       if(nRecOver<=0)break;
883       fh2NRecJetsPt->Fill(ptCut,nRecOver);
884     }
885     Float_t phi = leadingJet.Phi();
886     if(phi<0)phi+=TMath::Pi()*2.;    
887     Float_t eta = leadingJet.Eta();
888     Float_t pTback = 0;
889     if(externalBackground){
890       // carefull has to be filled in a task before
891       // todo, ReArrange to the botom
892       pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
893     }
894     pt = leadingJet.Pt() - pTback;
895     // correlation of leading jet with tracks
896     TIterator *recIter = recParticles.MakeIterator();
897     recIter->Reset();
898     AliVParticle *tmpRecTrack = 0;
899     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
900       Float_t tmpPt = tmpRecTrack->Pt();
901       // correlation
902       Float_t tmpPhi =  tmpRecTrack->Phi();     
903       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
904       Float_t dPhi = phi - tmpPhi;
905       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
906       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
907       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
908       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
909       if(cenClass>=0){
910         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
911         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
912       }
913
914     }  
915     
916    
917     TLorentzVector vecareab;
918     for(int j = 0; j < nRec;j++){
919       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
920       aodOutJet = 0;
921       nAodOutTracks = 0;
922       Float_t tmpPt = tmpRec.Pt();  
923       
924       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
925         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
926         aodOutJet->GetRefTracks()->Clear();
927         Double_t area1 = clustSeq.area(sortedJets[j]);
928         aodOutJet->SetEffArea(area1,0);
929         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);  
930         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());     
931         aodOutJet->SetVectorAreaCharged(&vecareab);
932       }
933
934
935       Float_t tmpPtBack = 0;
936       if(externalBackground){
937         // carefull has to be filled in a task before
938        // todo, ReArrange to the botom
939         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
940       }
941       tmpPt = tmpPt - tmpPtBack;
942       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
943       
944       fh1PtJetsRecIn->Fill(tmpPt);
945       // Fill Spectra with constituentsemacs
946       const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
947
948       fh1NConstRec->Fill(constituents.size());
949       fh2PtNch->Fill(nCh,tmpPt);
950       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
951       fh2NConstPt->Fill(tmpPt,constituents.size());
952       // loop over constiutents and fill spectrum
953    
954       for(unsigned int ic = 0; ic < constituents.size();ic++){
955         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
956         fh1PtJetConstRec->Fill(part->Pt());
957         if(aodOutJet){
958           aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
959           if(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
960         }
961         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
962       }
963       
964      // correlation
965      Float_t tmpPhi =  tmpRec.Phi();
966      Float_t tmpEta =  tmpRec.Eta();
967      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;        
968      if(j==0){
969        fh1PtJetsLeadingRecIn->Fill(tmpPt);
970        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
971        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
972        fh1NConstLeadingRec->Fill(constituents.size());
973        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
974        continue;
975      }
976      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
977      fh2JetEtaPt->Fill(tmpEta,tmpPt);
978      Float_t dPhi = phi - tmpPhi;
979      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
980      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
981      Float_t dEta = eta - tmpRec.Eta();
982      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
983      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
984      if(cenClass>=0){
985        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
986        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
987      }
988      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
989     }// loop over reconstructed jets
990    delete recIter;
991
992
993
994    // Add the random cones...
995    if(fNRandomCones>0&&fTCARandomConesOut){       
996      // create a random jet within the acceptance
997      Double_t etaMax = fTrackEtaWindow - fRparam;
998      Int_t nCone = 0;
999      Int_t nConeRan = 0;
1000      Double_t pTC = 1; // small number
1001      for(int ir = 0;ir < fNRandomCones;ir++){
1002        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1003        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1004        // massless jet
1005        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1006        Double_t pZC = pTC/TMath::Tan(thetaC);
1007        Double_t pXC = pTC * TMath::Cos(phiC);
1008        Double_t pYC = pTC * TMath::Sin(phiC);
1009        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1010        AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
1011        bool skip = false;
1012        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1013          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1014          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1015            skip = true;
1016            break;
1017          }
1018        }
1019        // test for overlap with previous cones to avoid double counting
1020        for(int iic = 0;iic<ir;iic++){
1021          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1022          if(iicone){
1023            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1024              skip = true;
1025              break;
1026            }
1027          }
1028        }
1029        if(skip)continue;
1030        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1031        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1032        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1033      }// loop over random cones creation
1034
1035   
1036      // loop over the reconstructed particles and add up the pT in the random cones
1037      // maybe better to loop over randomized particles not in the real jets...
1038      // but this by definition brings dow average energy in the whole  event
1039      AliAODJet vTmpRanR(1,0,0,1);
1040      for(int i = 0; i < recParticles.GetEntries(); i++){
1041        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1042        if(fTCARandomConesOut){
1043          for(int ir = 0;ir < fNRandomCones;ir++){
1044            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);  
1045            if(jC&&jC->DeltaR(vp)<fRparam){
1046              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1047              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1048            }
1049          }  
1050        }// add up energy in cone
1051       
1052        // the randomized input changes eta and phi, but keeps the p_T
1053        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1054          Double_t pTR = vp->Pt();
1055          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1056          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1057          
1058          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
1059          Double_t pZR = pTR/TMath::Tan(thetaR);
1060          
1061          Double_t pXR = pTR * TMath::Cos(phiR);
1062          Double_t pYR = pTR * TMath::Sin(phiR);
1063          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
1064          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1065          if(fTCARandomConesOutRan){
1066            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1067              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);  
1068              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1069                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1070                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1071              }
1072            }  
1073          }
1074        }
1075      }// loop over recparticles
1076     
1077      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1078      if(fTCARandomConesOut){
1079        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1080          // rescale the momntum vectors for the random cones
1081          
1082          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1083          if(rC){
1084            Double_t etaC = rC->Eta();
1085            Double_t phiC = rC->Phi();
1086            // massless jet, unit vector
1087            pTC = rC->ChargedBgEnergy();
1088            if(pTC<=0)pTC = 0.001; // for almost empty events
1089            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1090            Double_t pZC = pTC/TMath::Tan(thetaC);
1091            Double_t pXC = pTC * TMath::Cos(phiC);
1092            Double_t pYC = pTC * TMath::Sin(phiC);
1093            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1094            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1095            rC->SetBgEnergy(0,0);
1096            rC->SetEffArea(jetArea,0);
1097          }
1098        }
1099      }
1100      if(fTCARandomConesOutRan){
1101        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1102          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1103          // same wit random
1104          if(rC){
1105            Double_t etaC = rC->Eta();
1106            Double_t phiC = rC->Phi();
1107            // massless jet, unit vector
1108            pTC = rC->ChargedBgEnergy();
1109            if(pTC<=0)pTC = 0.001;// for almost empty events
1110            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1111            Double_t pZC = pTC/TMath::Tan(thetaC);
1112            Double_t pXC = pTC * TMath::Cos(phiC);
1113            Double_t pYC = pTC * TMath::Sin(phiC);
1114            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1115            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1116            rC->SetBgEnergy(0,0);
1117            rC->SetEffArea(jetArea,0);
1118          }
1119        }
1120      }
1121    }// if(fNRandomCones
1122   
1123    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1124  
1125
1126
1127
1128
1129    if(fAODJetBackgroundOut){
1130      vector<fastjet::PseudoJet> jets2=sortedJets;
1131      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
1132      Double_t bkg1=0;
1133      Double_t sigma1=0.;
1134      Double_t meanarea1=0.;
1135      Double_t bkg2=0;
1136      Double_t sigma2=0.;
1137      Double_t meanarea2=0.;
1138
1139      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1140      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1141
1142      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1143      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
1144      
1145      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1146      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1147      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1148      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1149
1150    }
1151   }
1152    
1153
1154   
1155   
1156  
1157   // fill track information
1158   Int_t nTrackOver = recParticles.GetSize();
1159   // do the same for tracks and jets
1160
1161   if(nTrackOver>0){
1162    TIterator *recIter = recParticles.MakeIterator();
1163    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1164    Float_t pt = tmpRec->Pt();
1165
1166    //    Printf("Leading track p_t %3.3E",pt);
1167    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1168      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1169      while(pt<ptCut&&tmpRec){
1170        nTrackOver--;
1171        tmpRec = (AliVParticle*)(recIter->Next()); 
1172        if(tmpRec){
1173          pt = tmpRec->Pt();
1174        }
1175      }
1176      if(nTrackOver<=0)break;
1177      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1178    }
1179    
1180    recIter->Reset();
1181    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1182    Float_t phi = leading->Phi();
1183    if(phi<0)phi+=TMath::Pi()*2.;    
1184    Float_t eta = leading->Eta();
1185    pt  = leading->Pt();
1186    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1187      Float_t tmpPt = tmpRec->Pt();
1188      Float_t tmpEta = tmpRec->Eta();
1189      fh1PtTracksRecIn->Fill(tmpPt);
1190      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1191      if(tmpRec==leading){
1192        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1193        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1194        continue;
1195      }
1196       // correlation
1197      Float_t tmpPhi =  tmpRec->Phi();
1198      
1199      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1200      Float_t dPhi = phi - tmpPhi;
1201      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1202      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1203      Float_t dEta = eta - tmpRec->Eta();
1204      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1205      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1206      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1207    }  
1208    delete recIter;
1209  }
1210
1211  // find the random jets
1212
1213  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1214   
1215  // fill the jet information from random track
1216  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1217  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1218
1219   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1220
1221  // loop over all jets an fill information, first on is the leading jet
1222
1223  Int_t nRecOverRan = inclusiveJetsRan.size();
1224  Int_t nRecRan     = inclusiveJetsRan.size();
1225
1226  if(inclusiveJetsRan.size()>0){
1227    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1228    Float_t pt = leadingJet.Pt();
1229    
1230    Int_t iCount = 0;
1231    TLorentzVector vecarearanb;
1232
1233    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1234      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1235       while(pt<ptCut&&iCount<nRecRan){
1236         nRecOverRan--;
1237         iCount++;
1238         if(iCount<nRecRan){
1239           pt = sortedJetsRan[iCount].perp();
1240         }
1241       }
1242       if(nRecOverRan<=0)break;
1243       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1244     }
1245     Float_t phi = leadingJet.Phi();
1246     if(phi<0)phi+=TMath::Pi()*2.;    
1247     pt = leadingJet.Pt();
1248
1249     // correlation of leading jet with random tracks
1250
1251     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1252       { 
1253         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1254         // correlation
1255         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1256         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1257         Float_t dPhi = phi - tmpPhi;
1258         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1259         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1260         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1261         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1262       }  
1263
1264     Int_t nAodOutJetsRan = 0;
1265      AliAODJet *aodOutJetRan = 0;
1266     for(int j = 0; j < nRecRan;j++){
1267       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1268       Float_t tmpPt = tmpRec.Pt();
1269       fh1PtJetsRecInRan->Fill(tmpPt);
1270       // Fill Spectra with constituents
1271       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1272       fh1NConstRecRan->Fill(constituents.size());
1273       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1274       fh2PtNchRan->Fill(nCh,tmpPt);
1275       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1276
1277
1278      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1279        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1280        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1281        aodOutJetRan->GetRefTracks()->Clear();
1282        aodOutJetRan->SetEffArea(arearan,0);
1283        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1284        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1285        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1286
1287      }
1288
1289       // correlation
1290       Float_t tmpPhi =  tmpRec.Phi();
1291       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1292
1293       if(j==0){
1294         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1295         fh1NConstLeadingRecRan->Fill(constituents.size());
1296         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1297         continue;
1298       }
1299     }  
1300
1301      
1302     if(fAODJetBackgroundOut){
1303      Double_t bkg3=0.;
1304      Double_t sigma3=0.;
1305      Double_t meanarea3=0.;
1306      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1307      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1308      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1309      /*
1310      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1311      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1312      */
1313     }
1314
1315
1316
1317  }
1318
1319
1320  // do the event selection if activated
1321  if(fJetTriggerPtCut>0){
1322    bool select = false;
1323    Float_t minPt = fJetTriggerPtCut;
1324    /*
1325    // hard coded for now ...
1326    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1327    if(cent<10)minPt = 50;
1328    else if(cent<30)minPt = 42;
1329    else if(cent<50)minPt = 28;
1330    else if(cent<80)minPt = 18;
1331    */
1332    float rho = 0;
1333    if(externalBackground)rho = externalBackground->GetBackground(2);
1334    if(fTCAJetsOut){
1335      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1336        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1337        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1338        if(ptSub>=minPt){
1339          select = true;
1340          break;
1341        }
1342      }
1343    }   
1344  
1345    if(select){
1346      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1347      fh1CentralitySelect->Fill(cent);
1348      fh1ZSelect->Fill(zVtx);
1349      aodH->SetFillAOD(kTRUE);
1350    }
1351  }
1352  if (fDebug > 2){
1353    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1354    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1355    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1356    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1357  }
1358  PostData(1, fHistList);
1359 }
1360
1361 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1362 {
1363   //
1364   // Terminate analysis
1365   //
1366   if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1367 }
1368
1369
1370 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1371
1372   //
1373   // get list of tracks/particles for different types
1374   //
1375
1376   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1377
1378   Int_t iCount = 0;
1379   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1380     if(type!=kTrackAODextraonly) {
1381       AliAODEvent *aod = 0;
1382       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1383       else aod = AODEvent();
1384       if(!aod){
1385         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1386         return iCount;
1387       }
1388       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1389         AliAODTrack *tr = aod->GetTrack(it);
1390         Bool_t bGood = false;
1391         if(fFilterType == 0)bGood = true;
1392         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1393         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1394         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1395           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
1396           continue;
1397         }
1398         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1399           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
1400           continue;
1401         }
1402         if(tr->Pt()<fTrackPtCut){
1403           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
1404           continue;
1405         }
1406         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
1407         list->Add(tr);
1408         iCount++;
1409       }
1410     }
1411     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1412       AliAODEvent *aod = 0;
1413       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1414       else aod = AODEvent();
1415       
1416       if(!aod){
1417         return iCount;
1418       }
1419       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1420       if(!aodExtraTracks)return iCount;
1421       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1422         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1423         if (!track) continue;
1424
1425         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1426         if(!trackAOD)continue;
1427         Bool_t bGood = false;
1428         if(fFilterType == 0)bGood = true;
1429         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1430         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1431         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1432         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1433         if(trackAOD->Pt()<fTrackPtCut) continue;
1434         list->Add(trackAOD);
1435         iCount++;
1436       }
1437     }
1438   }
1439   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1440     AliMCEvent* mcEvent = MCEvent();
1441     if(!mcEvent)return iCount;
1442     // we want to have alivpartilces so use get track
1443     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1444       if(!mcEvent->IsPhysicalPrimary(it))continue;
1445       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1446       if(type == kTrackKineAll){
1447         if(part->Pt()<fTrackPtCut)continue;
1448         list->Add(part);
1449         iCount++;
1450       }
1451       else if(type == kTrackKineCharged){
1452         if(part->Particle()->GetPDG()->Charge()==0)continue;
1453         if(part->Pt()<fTrackPtCut)continue;
1454         list->Add(part);
1455         iCount++;
1456       }
1457     }
1458   }
1459   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1460     AliAODEvent *aod = 0;
1461     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1462     else aod = AODEvent();
1463     if(!aod)return iCount;
1464     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1465     if(!tca)return iCount;
1466     for(int it = 0;it < tca->GetEntriesFast();++it){
1467       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1468       if(!part->IsPhysicalPrimary())continue;
1469       if(type == kTrackAODMCAll){
1470         if(part->Pt()<fTrackPtCut)continue;
1471         list->Add(part);
1472         iCount++;
1473       }
1474       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1475         if(part->Charge()==0)continue;
1476         if(part->Pt()<fTrackPtCut)continue;
1477         if(kTrackAODMCCharged){
1478           list->Add(part);
1479         }
1480         else {
1481           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1482           list->Add(part);
1483         }
1484         iCount++;
1485       }
1486     }
1487   }// AODMCparticle
1488   list->Sort();
1489   return iCount;
1490 }
1491
1492 /*
1493 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1494   for(int i = 0; i < particles.GetEntries(); i++){
1495     AliVParticle *vp = (AliVParticle*)particles.At(i);
1496     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1497     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1498     jInp.set_user_index(i);
1499     inputParticles.push_back(jInp);
1500   }
1501
1502   return 0;
1503
1504 }
1505 */