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