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