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