Fixing coverty reports
[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)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
673   }
674   //
675   // Execute analysis for current event
676   //
677   AliESDEvent *fESD = 0;
678   if(fUseAODTrackInput){    
679     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
680     if(!fAOD){
681       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
682       return;
683     }
684     // fethc the header
685   }
686   else{
687     //  assume that the AOD is in the general output...
688     fAOD  = AODEvent();
689     if(!fAOD){
690       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
691       return;
692     }
693     if(fDebug>0){
694       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
695     }
696   }
697   
698   Bool_t selectEvent =  false;
699   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
700
701   Float_t cent = 0;
702   Float_t zVtx  = 0;
703   Int_t cenClass = -1;
704   if(fAOD){
705     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
706     TString vtxTitle(vtxAOD->GetTitle());
707     zVtx = vtxAOD->GetZ();
708
709     cent = fAOD->GetHeader()->GetCentrality();
710     if(cent<10)cenClass = 0;
711     else if(cent<30)cenClass = 1;
712     else if(cent<50)cenClass = 2;
713     else if(cent<80)cenClass = 3;
714     if(physicsSelection){
715       fh1CentralityPhySel->Fill(cent);
716       fh1ZPhySel->Fill(zVtx);
717     }
718
719
720     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
721         Float_t yvtx = vtxAOD->GetY();
722         Float_t xvtx = vtxAOD->GetX();
723         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
724         if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
725           if(physicsSelection){
726             selectEvent = true;
727           }
728         }
729     }
730     if(fCentCutUp>0){
731       if(cent<fCentCutLo||cent>fCentCutUp){
732         selectEvent = false;
733       }
734     }
735
736   }
737   if(!selectEvent){
738     PostData(1, fHistList);
739     return;
740   }
741   fh1Centrality->Fill(cent);  
742   fh1Z->Fill(zVtx);
743   fh1Trials->Fill("#sum{ntrials}",1);
744   
745
746   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
747
748   // ==== General variables needed
749
750
751
752   // we simply fetch the tracks/mc particles as a list of AliVParticles
753
754   TList recParticles;
755   TList genParticles;
756
757   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
758   Float_t nCh = recParticles.GetEntries(); 
759   fh1Nch->Fill(nCh);
760   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
761   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
762   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
763
764   // find the jets....
765
766   vector<fastjet::PseudoJet> inputParticlesRec;
767   vector<fastjet::PseudoJet> inputParticlesRecRan;
768   
769   // Generate the random cones
770   
771   AliAODJet vTmpRan(1,0,0,1);
772   for(int i = 0; i < recParticles.GetEntries(); i++){
773     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
774     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
775     // we take total momentum here
776     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
777     jInp.set_user_index(i);
778     inputParticlesRec.push_back(jInp);
779
780     // the randomized input changes eta and phi, but keeps the p_T
781     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
782       Double_t pT = vp->Pt();
783       Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
784       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
785       
786       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
787       Double_t pZ = pT/TMath::Tan(theta);
788
789       Double_t pX = pT * TMath::Cos(phi);
790       Double_t pY = pT * TMath::Sin(phi);
791       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
792       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
793
794       jInpRan.set_user_index(i);
795       inputParticlesRecRan.push_back(jInpRan);
796       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
797     }
798
799     // fill the tref array, only needed when we write out jets
800     if(jarray){
801       if(i == 0){
802         fRef->Delete(); // make sure to delete before placement new...
803         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
804       }
805       fRef->Add(vp);
806     }
807   }// recparticles
808
809   if(inputParticlesRec.size()==0){
810     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
811     PostData(1, fHistList);
812     return;
813   }
814   
815   // run fast jet
816   // employ setters for these...
817
818  
819   // now create the object that holds info about ghosts                        
820   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
821     // reduce CPU time...
822     fGhostArea = 0.5; 
823     fActiveAreaRepeats = 0; 
824   }
825    
826  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
827   fastjet::AreaType areaType =   fastjet::active_area;
828   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
829   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
830   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
831   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
832   
833   //range where to compute background
834   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
835   phiMin = 0;
836   phiMax = 2*TMath::Pi();
837   rapMax = fGhostEtamax - fRparam;
838   rapMin = - fGhostEtamax + fRparam;
839   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
840  
841
842
843   inclusiveJets = clustSeq.inclusive_jets();
844   sortedJets    = sorted_by_pt(inclusiveJets);
845  
846   fh1NJetsRec->Fill(sortedJets.size());
847
848  // loop over all jets an fill information, first on is the leading jet
849
850   Int_t nRecOver = inclusiveJets.size();
851   Int_t nRec     = inclusiveJets.size();
852   if(inclusiveJets.size()>0){
853     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
854     Double_t area = clustSeq.area(sortedJets[0]);
855     leadingJet.SetEffArea(area,0);
856     Float_t pt = leadingJet.Pt();
857     Int_t nAodOutJets = 0;
858     Int_t nAodOutTracks = 0;
859     AliAODJet *aodOutJet = 0;
860
861     Int_t iCount = 0;
862     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
863       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
864       while(pt<ptCut&&iCount<nRec){
865         nRecOver--;
866         iCount++;
867         if(iCount<nRec){
868           pt = sortedJets[iCount].perp();
869         }
870       }
871       if(nRecOver<=0)break;
872       fh2NRecJetsPt->Fill(ptCut,nRecOver);
873     }
874     Float_t phi = leadingJet.Phi();
875     if(phi<0)phi+=TMath::Pi()*2.;    
876     Float_t eta = leadingJet.Eta();
877     Float_t pTback = 0;
878     if(externalBackground){
879       // carefull has to be filled in a task before
880       // todo, ReArrange to the botom
881       pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
882     }
883     pt = leadingJet.Pt() - pTback;
884     // correlation of leading jet with tracks
885     TIterator *recIter = recParticles.MakeIterator();
886     recIter->Reset();
887     AliVParticle *tmpRecTrack = 0;
888     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
889       Float_t tmpPt = tmpRecTrack->Pt();
890       // correlation
891       Float_t tmpPhi =  tmpRecTrack->Phi();     
892       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
893       Float_t dPhi = phi - tmpPhi;
894       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
895       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
896       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
897       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
898       if(cenClass>=0){
899         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
900         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
901       }
902
903     }  
904     
905    
906     TLorentzVector vecareab;
907     for(int j = 0; j < nRec;j++){
908       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
909       aodOutJet = 0;
910       nAodOutTracks = 0;
911       Float_t tmpPt = tmpRec.Pt();  
912       
913       if(tmpPt>fJetOutputMinPt&&jarray){// cut on the non-background subtracted...
914         aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
915         Double_t area1 = clustSeq.area(sortedJets[j]);
916         aodOutJet->SetEffArea(area1,0);
917         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);  
918         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());     
919         aodOutJet->SetVectorAreaCharged(&vecareab);
920
921
922       }
923
924
925       Float_t tmpPtBack = 0;
926       if(externalBackground){
927         // carefull has to be filled in a task before
928        // todo, ReArrange to the botom
929         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
930       }
931       tmpPt = tmpPt - tmpPtBack;
932       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
933       
934       fh1PtJetsRecIn->Fill(tmpPt);
935       // Fill Spectra with constituents
936       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
937
938       fh1NConstRec->Fill(constituents.size());
939       fh2PtNch->Fill(nCh,tmpPt);
940       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
941       fh2NConstPt->Fill(tmpPt,constituents.size());
942       // loop over constiutents and fill spectrum
943       
944       // Add the jet information and the track references to the output AOD
945             
946
947       if(fNRandomCones>0){       
948         // create a random jet within the acceptance
949         Double_t etaMax = 0.8 - fRparam;
950         Int_t nCone = 0;
951         Int_t nConeRan = 0;
952         Double_t pTC = 1; // small number
953         for(int ir = 0;ir < fNRandomCones;ir++){
954           Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
955           Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
956           // massless jet
957           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
958           Double_t pZC = pTC/TMath::Tan(thetaC);
959           Double_t pXC = pTC * TMath::Cos(phiC);
960           Double_t pYC = pTC * TMath::Sin(phiC);
961           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
962           AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
963           bool skip = false;
964           for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
965             AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
966             if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
967               skip = true;
968               break;
969             }
970           }
971           // test for overlap with previous cones to avoid double counting
972           for(int iic = 0;iic<ir;iic++){
973             AliAODJet *iicone = (AliAODJet*)rConeArray->At(iic);
974             if(iicone){
975               if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
976                 skip = true;
977                 break;
978               }
979             }
980           }
981           if(skip)continue;
982           tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
983           if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
984           if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
985         }// random cones
986
987
988         // loop over the reconstructed particles and add up the pT in the random cones
989         // maybe better to loop over randomized particles not in the real jets...
990         // but this by definition brings dow average energy in the whole  event
991         AliAODJet vTmpRanR(1,0,0,1);
992         for(int i = 0; i < recParticles.GetEntries(); i++){
993           AliVParticle *vp = (AliVParticle*)recParticles.At(i);
994           if(rConeArray){
995             for(int ir = 0;ir < fNRandomCones;ir++){
996               AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);  
997               if(jC&&jC->DeltaR(vp)<fRparam){
998                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
999               }
1000             }  
1001           }// add up energy in cone
1002       
1003           // the randomized input changes eta and phi, but keeps the p_T
1004           if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1005             Double_t pTR = vp->Pt();
1006             Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
1007             Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1008             
1009             Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
1010             Double_t pZR = pTR/TMath::Tan(thetaR);
1011         
1012             Double_t pXR = pTR * TMath::Cos(phiR);
1013             Double_t pYR = pTR * TMath::Sin(phiR);
1014             Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
1015             vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1016             if(rConeArrayRan){
1017               for(int ir = 0;ir < fNRandomCones;ir++){
1018                 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);  
1019                 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1020                   jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1021                 }
1022               }  
1023             }
1024           }
1025         }// loop over recparticles
1026     
1027         Float_t jetArea = fRparam*fRparam*TMath::Pi();
1028         for(int ir = 0;ir < fNRandomCones;ir++){
1029           // rescale the momntum vectors for the random cones
1030           if(!rConeArray)continue;
1031           AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1032           if(rC){
1033             Double_t etaC = rC->Eta();
1034             Double_t phiC = rC->Phi();
1035             // massless jet, unit vector
1036             pTC = rC->ChargedBgEnergy();
1037             if(pTC<=0)pTC = 0.1; // for almost empty events
1038             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1039             Double_t pZC = pTC/TMath::Tan(thetaC);
1040             Double_t pXC = pTC * TMath::Cos(phiC);
1041             Double_t pYC = pTC * TMath::Sin(phiC);
1042             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1043             rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1044             rC->SetBgEnergy(0,0);
1045             rC->SetEffArea(jetArea,0);
1046           }
1047           rC = (AliAODJet*)rConeArrayRan->At(ir);
1048           // same wit random
1049           if(rC){
1050             Double_t etaC = rC->Eta();
1051             Double_t phiC = rC->Phi();
1052             // massless jet, unit vector
1053             pTC = rC->ChargedBgEnergy();
1054             if(pTC<=0)pTC = 0.1;// for almost empty events
1055             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1056             Double_t pZC = pTC/TMath::Tan(thetaC);
1057             Double_t pXC = pTC * TMath::Cos(phiC);
1058             Double_t pYC = pTC * TMath::Sin(phiC);
1059             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1060             rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1061             rC->SetBgEnergy(0,0);
1062             rC->SetEffArea(jetArea,0);
1063           }
1064         }
1065       }// if(fUseBackgroundCalc
1066      
1067      for(unsigned int ic = 0; ic < constituents.size();ic++){
1068        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1069        fh1PtJetConstRec->Fill(part->Pt());
1070        if(aodOutJet){
1071          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1072        }
1073        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1074      }
1075      
1076      // correlation
1077      Float_t tmpPhi =  tmpRec.Phi();
1078      Float_t tmpEta =  tmpRec.Eta();
1079      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;        
1080      if(j==0){
1081        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1082        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1083        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1084        fh1NConstLeadingRec->Fill(constituents.size());
1085        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1086        continue;
1087      }
1088      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1089      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1090      Float_t dPhi = phi - tmpPhi;
1091      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1092      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1093      Float_t dEta = eta - tmpRec.Eta();
1094      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1095      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1096      if(cenClass>=0){
1097        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1098        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1099      }
1100      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1101     }
1102    delete recIter;
1103   
1104    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1105  
1106
1107    if(evBkg){
1108      vector<fastjet::PseudoJet> jets2=sortedJets;
1109      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
1110      Double_t bkg1=0;
1111      Double_t sigma1=0.;
1112      Double_t meanarea1=0.;
1113      Double_t bkg2=0;
1114      Double_t sigma2=0.;
1115      Double_t meanarea2=0.;
1116
1117      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1118      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1119
1120      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1121      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
1122      
1123      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1124      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1125      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1126      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1127
1128    }
1129   }
1130    
1131
1132   
1133   
1134  
1135   // fill track information
1136   Int_t nTrackOver = recParticles.GetSize();
1137   // do the same for tracks and jets
1138
1139   if(nTrackOver>0){
1140    TIterator *recIter = recParticles.MakeIterator();
1141    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1142    Float_t pt = tmpRec->Pt();
1143
1144    //    Printf("Leading track p_t %3.3E",pt);
1145    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1146      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1147      while(pt<ptCut&&tmpRec){
1148        nTrackOver--;
1149        tmpRec = (AliVParticle*)(recIter->Next()); 
1150        if(tmpRec){
1151          pt = tmpRec->Pt();
1152        }
1153      }
1154      if(nTrackOver<=0)break;
1155      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1156    }
1157    
1158    recIter->Reset();
1159    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1160    Float_t phi = leading->Phi();
1161    if(phi<0)phi+=TMath::Pi()*2.;    
1162    Float_t eta = leading->Eta();
1163    pt  = leading->Pt();
1164    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1165      Float_t tmpPt = tmpRec->Pt();
1166      Float_t tmpEta = tmpRec->Eta();
1167      fh1PtTracksRecIn->Fill(tmpPt);
1168      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1169      if(tmpRec==leading){
1170        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1171        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1172        continue;
1173      }
1174       // correlation
1175      Float_t tmpPhi =  tmpRec->Phi();
1176      
1177      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1178      Float_t dPhi = phi - tmpPhi;
1179      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1180      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1181      Float_t dEta = eta - tmpRec->Eta();
1182      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1183      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1184      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1185    }  
1186    delete recIter;
1187  }
1188
1189  // find the random jets
1190  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1191  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1192   
1193  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1194  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1195
1196  // fill the jet information from random track
1197
1198   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1199
1200  // loop over all jets an fill information, first on is the leading jet
1201
1202  Int_t nRecOverRan = inclusiveJetsRan.size();
1203  Int_t nRecRan     = inclusiveJetsRan.size();
1204
1205  if(inclusiveJetsRan.size()>0){
1206    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1207    Float_t pt = leadingJet.Pt();
1208    
1209    Int_t iCount = 0;
1210    TLorentzVector vecarearanb;
1211
1212    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1213      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1214       while(pt<ptCut&&iCount<nRecRan){
1215         nRecOverRan--;
1216         iCount++;
1217         if(iCount<nRecRan){
1218           pt = sortedJetsRan[iCount].perp();
1219         }
1220       }
1221       if(nRecOverRan<=0)break;
1222       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1223     }
1224     Float_t phi = leadingJet.Phi();
1225     if(phi<0)phi+=TMath::Pi()*2.;    
1226     pt = leadingJet.Pt();
1227
1228     // correlation of leading jet with random tracks
1229
1230     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1231       { 
1232         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1233         // correlation
1234         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1235         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1236         Float_t dPhi = phi - tmpPhi;
1237         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1238         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1239         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1240         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1241       }  
1242
1243     Int_t nAodOutJetsRan = 0;
1244      AliAODJet *aodOutJetRan = 0;
1245     for(int j = 0; j < nRecRan;j++){
1246       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1247       Float_t tmpPt = tmpRec.Pt();
1248       fh1PtJetsRecInRan->Fill(tmpPt);
1249       // Fill Spectra with constituents
1250       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1251       fh1NConstRecRan->Fill(constituents.size());
1252       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1253       fh2PtNchRan->Fill(nCh,tmpPt);
1254       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1255
1256
1257      if(tmpPt>fJetOutputMinPt&&jarrayran){
1258        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1259        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1260        
1261        aodOutJetRan->SetEffArea(arearan,0);
1262        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1263        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1264        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1265
1266      }
1267
1268
1269
1270      
1271      for(unsigned int ic = 0; ic < constituents.size();ic++){
1272        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1273        // fh1PtJetConstRec->Fill(part->Pt());
1274        if(aodOutJetRan){
1275          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1276        }
1277      }
1278       
1279
1280
1281
1282       // correlation
1283       Float_t tmpPhi =  tmpRec.Phi();
1284       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1285
1286       if(j==0){
1287         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1288         fh1NConstLeadingRecRan->Fill(constituents.size());
1289         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1290         continue;
1291       }
1292     }  
1293
1294      
1295     if(evBkg){
1296      Double_t bkg3=0.;
1297      Double_t sigma3=0.;
1298      Double_t meanarea3=0.;
1299      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1300      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1301      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1302      /*
1303      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1304      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1305      */
1306     }
1307
1308
1309
1310  }
1311
1312
1313  // do the event selection if activated
1314  if(fJetTriggerPtCut>0){
1315    bool select = false;
1316    Float_t minPt = fJetTriggerPtCut;
1317    /*
1318    // hard coded for now ...
1319    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1320    if(cent<10)minPt = 50;
1321    else if(cent<30)minPt = 42;
1322    else if(cent<50)minPt = 28;
1323    else if(cent<80)minPt = 18;
1324    */
1325    float rho = 0;
1326    if(externalBackground)rho = externalBackground->GetBackground(2);
1327    if(jarray){
1328      for(int i = 0;i < jarray->GetEntriesFast();i++){
1329        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1330        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1331        if(ptSub>=minPt){
1332          select = true;
1333          break;
1334        }
1335      }
1336    }   
1337
1338    if(select){
1339      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1340      fh1CentralitySelect->Fill(cent);
1341      fh1ZSelect->Fill(zVtx);
1342      aodH->SetFillAOD(kTRUE);
1343    }
1344  }
1345
1346  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1347  PostData(1, fHistList);
1348 }
1349
1350 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1351 {
1352 // Terminate analysis
1353 //
1354     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1355 }
1356
1357
1358 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1359
1360   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1361
1362   Int_t iCount = 0;
1363   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1364     if(type!=kTrackAODextraonly) {
1365       AliAODEvent *aod = 0;
1366       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1367       else aod = AODEvent();
1368       if(!aod){
1369         return iCount;
1370       }
1371       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1372         AliAODTrack *tr = aod->GetTrack(it);
1373         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1374         if(TMath::Abs(tr->Eta())>0.9)continue;
1375         if(tr->Pt()<fTrackPtCut)continue;
1376         list->Add(tr);
1377         iCount++;
1378       }
1379     }
1380     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1381       AliAODEvent *aod = 0;
1382       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1383       else aod = AODEvent();
1384       
1385       if(!aod){
1386         return iCount;
1387       }
1388       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1389       if(!aodExtraTracks)return iCount;
1390       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1391         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1392         if (!track) continue;
1393
1394         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1395         if(!trackAOD)continue;
1396         if(trackAOD->Pt()<fTrackPtCut) continue;
1397         list->Add(trackAOD);
1398         iCount++;
1399       }
1400     }
1401   }
1402   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1403     AliMCEvent* mcEvent = MCEvent();
1404     if(!mcEvent)return iCount;
1405     // we want to have alivpartilces so use get track
1406     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1407       if(!mcEvent->IsPhysicalPrimary(it))continue;
1408       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1409       if(type == kTrackKineAll){
1410         if(part->Pt()<fTrackPtCut)continue;
1411         list->Add(part);
1412         iCount++;
1413       }
1414       else if(type == kTrackKineCharged){
1415         if(part->Particle()->GetPDG()->Charge()==0)continue;
1416         if(part->Pt()<fTrackPtCut)continue;
1417         list->Add(part);
1418         iCount++;
1419       }
1420     }
1421   }
1422   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1423     AliAODEvent *aod = 0;
1424     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1425     else aod = AODEvent();
1426     if(!aod)return iCount;
1427     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1428     if(!tca)return iCount;
1429     for(int it = 0;it < tca->GetEntriesFast();++it){
1430       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1431       if(!part->IsPhysicalPrimary())continue;
1432       if(type == kTrackAODMCAll){
1433         if(part->Pt()<fTrackPtCut)continue;
1434         list->Add(part);
1435         iCount++;
1436       }
1437       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1438         if(part->Charge()==0)continue;
1439         if(part->Pt()<fTrackPtCut)continue;
1440         if(kTrackAODMCCharged){
1441           list->Add(part);
1442         }
1443         else {
1444           if(TMath::Abs(part->Eta())>0.9)continue;
1445           list->Add(part);
1446         }
1447         iCount++;
1448       }
1449     }
1450   }// AODMCparticle
1451   list->Sort();
1452   return iCount;
1453 }
1454
1455 /*
1456 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1457   for(int i = 0; i < particles.GetEntries(); i++){
1458     AliVParticle *vp = (AliVParticle*)particles.At(i);
1459     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1460     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1461     jInp.set_user_index(i);
1462     inputParticles.push_back(jInp);
1463   }
1464
1465   return 0;
1466
1467 }
1468 */