d46f142364dcb675d607cb5285c9c4f34029b69d
[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 <TRandom.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 }
73
74 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
75   fAOD(0x0),
76   fAODExtension(0x0),
77   fRef(new TRefArray),
78   fUseAODTrackInput(kFALSE),
79   fUseAODMCInput(kFALSE),
80   fUseGlobalSelection(kFALSE),
81   fFilterMask(0),
82   fTrackTypeRec(kTrackUndef),
83   fTrackTypeGen(kTrackUndef),  
84   fNSkipLeadingRan(0),
85   fAvgTrials(1),
86   fExternalWeight(1),    
87   fRecEtaWindow(0.5),
88   fTrackPtCut(0.),                                                      
89   fJetOutputMinPt(1),
90   fNonStdBranch(""),
91   fNonStdFile(""),
92   fRparam(1.0), 
93   fAlgorithm(fastjet::kt_algorithm),
94   fStrategy(fastjet::Best),
95   fRecombScheme(fastjet::BIpt_scheme),
96   fAreaType(fastjet::active_area), 
97   fGhostArea(0.01),
98   fActiveAreaRepeats(1),
99   fGhostEtamax(1.5),
100   fh1Xsec(0x0),
101   fh1Trials(0x0),
102   fh1PtHard(0x0),
103   fh1PtHardNoW(0x0),  
104   fh1PtHardTrials(0x0),
105   fh1NJetsRec(0x0),
106   fh1NConstRec(0x0),
107   fh1NConstLeadingRec(0x0),
108   fh1PtJetsRecIn(0x0),
109   fh1PtJetsLeadingRecIn(0x0),
110   fh1PtJetConstRec(0x0),
111   fh1PtJetConstLeadingRec(0x0),
112   fh1PtTracksRecIn(0x0),
113   fh1PtTracksLeadingRecIn(0x0),
114   fh1NJetsRecRan(0x0),
115   fh1NConstRecRan(0x0),
116   fh1PtJetsLeadingRecInRan(0x0),
117   fh1NConstLeadingRecRan(0x0),
118   fh1PtJetsRecInRan(0x0),
119   fh1PtTracksGenIn(0x0),
120   fh1Nch(0x0),
121   fh2NRecJetsPt(0x0),
122   fh2NRecTracksPt(0x0),
123   fh2NConstPt(0x0),
124   fh2NConstLeadingPt(0x0),
125   fh2JetPhiEta(0x0),
126   fh2LeadingJetPhiEta(0x0),
127   fh2JetEtaPt(0x0),
128   fh2LeadingJetEtaPt(0x0),
129   fh2TrackEtaPt(0x0),
130   fh2LeadingTrackEtaPt(0x0),
131   fh2JetsLeadingPhiEta(0x0),
132   fh2JetsLeadingPhiPt(0x0),
133   fh2TracksLeadingPhiEta(0x0),
134   fh2TracksLeadingPhiPt(0x0),
135   fh2TracksLeadingJetPhiPt(0x0),
136   fh2JetsLeadingPhiPtW(0x0),
137   fh2TracksLeadingPhiPtW(0x0),
138   fh2TracksLeadingJetPhiPtW(0x0),
139   fh2NRecJetsPtRan(0x0),
140   fh2NConstPtRan(0x0),
141   fh2NConstLeadingPtRan(0x0),
142   fh2PtNch(0x0),
143   fh2PtNchRan(0x0),
144   fh2PtNchN(0x0),
145   fh2PtNchNRan(0x0),
146   fh2TracksLeadingJetPhiPtRan(0x0),
147   fh2TracksLeadingJetPhiPtWRan(0x0),
148   fHistList(0x0)  
149 {
150   for(int i = 0;i<3;i++){
151     fh1BiARandomCones[i] = 0;
152     fh1BiARandomConesRan[i] = 0;
153   }
154 }
155
156 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
157   AliAnalysisTaskSE(name),
158   fAOD(0x0),
159   fAODExtension(0x0),
160   fRef(new TRefArray),
161   fUseAODTrackInput(kFALSE),
162   fUseAODMCInput(kFALSE),
163   fUseGlobalSelection(kFALSE),
164   fFilterMask(0),
165   fTrackTypeRec(kTrackUndef),
166   fTrackTypeGen(kTrackUndef),
167   fNSkipLeadingRan(0),
168   fAvgTrials(1),
169   fExternalWeight(1),    
170   fRecEtaWindow(0.5),
171   fTrackPtCut(0.),                                                      
172   fJetOutputMinPt(1),
173   fNonStdBranch(""),
174   fNonStdFile(""),
175   fRparam(1.0), 
176   fAlgorithm(fastjet::kt_algorithm),
177   fStrategy(fastjet::Best),
178   fRecombScheme(fastjet::BIpt_scheme),
179   fAreaType(fastjet::active_area), 
180   fGhostArea(0.01),
181   fActiveAreaRepeats(1),
182   fGhostEtamax(1.5),
183   fh1Xsec(0x0),
184   fh1Trials(0x0),
185   fh1PtHard(0x0),
186   fh1PtHardNoW(0x0),  
187   fh1PtHardTrials(0x0),
188   fh1NJetsRec(0x0),
189   fh1NConstRec(0x0),
190   fh1NConstLeadingRec(0x0),
191   fh1PtJetsRecIn(0x0),
192   fh1PtJetsLeadingRecIn(0x0),
193   fh1PtJetConstRec(0x0),
194   fh1PtJetConstLeadingRec(0x0),
195   fh1PtTracksRecIn(0x0),
196   fh1PtTracksLeadingRecIn(0x0),
197   fh1NJetsRecRan(0x0),
198   fh1NConstRecRan(0x0),
199   fh1PtJetsLeadingRecInRan(0x0),
200   fh1NConstLeadingRecRan(0x0),
201   fh1PtJetsRecInRan(0x0),
202   fh1PtTracksGenIn(0x0),
203   fh1Nch(0x0),
204   fh2NRecJetsPt(0x0),
205   fh2NRecTracksPt(0x0),
206   fh2NConstPt(0x0),
207   fh2NConstLeadingPt(0x0),
208   fh2JetPhiEta(0x0),
209   fh2LeadingJetPhiEta(0x0),
210   fh2JetEtaPt(0x0),
211   fh2LeadingJetEtaPt(0x0),
212   fh2TrackEtaPt(0x0),
213   fh2LeadingTrackEtaPt(0x0),
214   fh2JetsLeadingPhiEta(0x0),
215   fh2JetsLeadingPhiPt(0x0),
216   fh2TracksLeadingPhiEta(0x0),
217   fh2TracksLeadingPhiPt(0x0),
218   fh2TracksLeadingJetPhiPt(0x0),
219   fh2JetsLeadingPhiPtW(0x0),
220   fh2TracksLeadingPhiPtW(0x0),
221   fh2TracksLeadingJetPhiPtW(0x0),
222   fh2NRecJetsPtRan(0x0),
223   fh2NConstPtRan(0x0),
224   fh2NConstLeadingPtRan(0x0),
225   fh2PtNch(0x0),
226   fh2PtNchRan(0x0),
227   fh2PtNchN(0x0),
228   fh2PtNchNRan(0x0),
229   fh2TracksLeadingJetPhiPtRan(0x0),
230   fh2TracksLeadingJetPhiPtWRan(0x0),
231   fHistList(0x0)
232 {
233   for(int i = 0;i<3;i++){
234     fh1BiARandomCones[i] = 0;
235     fh1BiARandomConesRan[i] = 0;
236   }
237  DefineOutput(1, TList::Class());  
238 }
239
240
241
242 Bool_t AliAnalysisTaskJetCluster::Notify()
243 {
244   //
245   // Implemented Notify() to read the cross sections
246   // and number of trials from pyxsec.root
247   // 
248   return kTRUE;
249 }
250
251 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
252 {
253
254   //
255   // Create the output container
256   //
257
258
259
260
261   // Connect the AOD
262
263
264   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
265
266   
267
268   if(fNonStdBranch.Length()!=0)
269     {
270       // only create the output branch if we have a name
271       // Create a new branch for jets...
272       //  -> cleared in the UserExec....
273       // here we can also have the case that the brnaches are written to a separate file
274
275       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
276       tca->SetName(fNonStdBranch.Data());
277       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
278
279    
280       TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
281       tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
282       AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
283
284
285      if(!AODEvent() || !AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
286         AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
287         evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
288         AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
289        
290      }
291
292       if(fNonStdFile.Length()!=0){
293         // 
294         // case that we have an AOD extension we need to fetch the jets from the extended output
295         // we identifay the extension aod event by looking for the branchname
296         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
297         TObjArray* extArray = aodH->GetExtensions();
298         if (extArray) {
299           TIter next(extArray);
300           while ((fAODExtension=(AliAODExtension*)next())){
301             TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
302             if(fDebug>10){
303               Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
304               fAODExtension->GetAOD()->Dump();
305             }
306             if(obj){
307               if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
308               break;
309             }
310             fAODExtension = 0;
311           }
312         }
313       }
314     }
315
316
317   OpenFile(1);
318   if(!fHistList)fHistList = new TList();
319   fHistList->SetOwner();
320
321   Bool_t oldStatus = TH1::AddDirectoryStatus();
322   TH1::AddDirectory(kFALSE);
323
324   //
325   //  Histogram
326     
327   const Int_t nBinPt = 200;
328   Double_t binLimitsPt[nBinPt+1];
329   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
330     if(iPt == 0){
331       binLimitsPt[iPt] = 0.0;
332     }
333     else {// 1.0
334       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 0.5;
335     }
336   }
337   
338   const Int_t nBinPhi = 90;
339   Double_t binLimitsPhi[nBinPhi+1];
340   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
341     if(iPhi==0){
342       binLimitsPhi[iPhi] = -1.*TMath::Pi();
343     }
344     else{
345       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
346     }
347   }
348
349
350
351   const Int_t nBinEta = 40;
352   Double_t binLimitsEta[nBinEta+1];
353   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
354     if(iEta==0){
355       binLimitsEta[iEta] = -2.0;
356     }
357     else{
358       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
359     }
360   }
361
362   const int nChMax = 100;
363
364   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
365   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
366
367   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
368   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
369
370
371   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
372   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
373
374   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
375   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
376   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
377   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
378
379
380   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
381   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
382   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
383
384   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
385   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
386   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
387   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
388   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
389   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
390   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
391   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
392   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
393   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
394
395
396   for(int i = 0;i<3;i++){
397     fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
398     fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
399   }
400
401   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
402   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
403   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
404   // 
405
406
407   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
408   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
409   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
410   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
411
412   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
413   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);
414   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);
415   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);
416
417
418
419   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
420                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
421   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
422                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
423
424   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
425                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
426   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
427                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
428
429   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
430                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
431   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
432                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
433
434
435
436   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
437                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
438   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
439                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
440   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
441                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
442   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
443                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
444   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
445                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
446   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
447                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
448
449   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
450                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
451   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
452                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
453
454   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
455                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
456   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
457                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
458
459
460
461   const Int_t saveLevel = 3; // large save level more histos
462   if(saveLevel>0){
463     fHistList->Add(fh1Xsec);
464     fHistList->Add(fh1Trials);
465
466     fHistList->Add(fh1NJetsRec);
467     fHistList->Add(fh1NConstRec);
468     fHistList->Add(fh1NConstLeadingRec);
469     fHistList->Add(fh1PtJetsRecIn);
470     fHistList->Add(fh1PtJetsLeadingRecIn);
471     fHistList->Add(fh1PtTracksRecIn);
472     fHistList->Add(fh1PtTracksLeadingRecIn);
473     fHistList->Add(fh1PtJetConstRec);
474     fHistList->Add(fh1PtJetConstLeadingRec);
475     fHistList->Add(fh1NJetsRecRan);
476     fHistList->Add(fh1NConstRecRan);
477     fHistList->Add(fh1PtJetsLeadingRecInRan);
478     fHistList->Add(fh1NConstLeadingRecRan);
479     fHistList->Add(fh1PtJetsRecInRan);
480     fHistList->Add(fh1Nch);
481     for(int i = 0;i<3;i++){
482       fHistList->Add(fh1BiARandomCones[i]);
483       fHistList->Add(fh1BiARandomConesRan[i]);
484     }
485     fHistList->Add(fh2NRecJetsPt);
486     fHistList->Add(fh2NRecTracksPt);
487     fHistList->Add(fh2NConstPt);
488     fHistList->Add(fh2NConstLeadingPt);
489     fHistList->Add(fh2PtNch);
490     fHistList->Add(fh2PtNchRan);
491     fHistList->Add(fh2PtNchN);
492     fHistList->Add(fh2PtNchNRan);
493     fHistList->Add(fh2JetPhiEta);
494     fHistList->Add(fh2LeadingJetPhiEta);
495     fHistList->Add(fh2JetEtaPt);
496     fHistList->Add(fh2LeadingJetEtaPt);
497     fHistList->Add(fh2TrackEtaPt);
498     fHistList->Add(fh2LeadingTrackEtaPt);
499     fHistList->Add(fh2JetsLeadingPhiEta );
500     fHistList->Add(fh2JetsLeadingPhiPt);
501     fHistList->Add(fh2TracksLeadingPhiEta);
502     fHistList->Add(fh2TracksLeadingPhiPt);
503     fHistList->Add(fh2TracksLeadingJetPhiPt);
504     fHistList->Add(fh2JetsLeadingPhiPtW);
505     fHistList->Add(fh2TracksLeadingPhiPtW);
506     fHistList->Add(fh2TracksLeadingJetPhiPtW);
507     fHistList->Add(fh2NRecJetsPtRan);
508     fHistList->Add(fh2NConstPtRan);
509     fHistList->Add(fh2NConstLeadingPtRan);
510     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
511     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
512     }
513
514   // =========== Switch on Sumw2 for all histos ===========
515   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
516     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
517     if (h1){
518       h1->Sumw2();
519       continue;
520     }
521     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
522     if(hn)hn->Sumw2();
523   }
524   TH1::AddDirectory(oldStatus);
525 }
526
527 void AliAnalysisTaskJetCluster::Init()
528 {
529   //
530   // Initialization
531   //
532
533   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
534
535 }
536
537 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
538 {
539
540   if(fUseGlobalSelection){
541     // no selection by the service task, we continue
542     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
543     PostData(1, fHistList);
544     return;
545   }
546
547   // handle and reset the output jet branch 
548   // only need this once
549   TClonesArray* jarray = 0;  
550   TClonesArray* jarrayran = 0;  
551   AliAODJetEventBackground*  evBkg = 0;
552   if(fNonStdBranch.Length()!=0) {
553     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
554     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
555     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
556     if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
557     if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
558     if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
559
560     if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
561     if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
562     if(evBkg)evBkg->Reset(); 
563
564   }
565
566
567  
568   //
569   // Execute analysis for current event
570   //
571   AliESDEvent *fESD = 0;
572   if(fUseAODTrackInput){    
573     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
574     if(!fAOD){
575       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
576       return;
577     }
578     // fethc the header
579   }
580   else{
581     //  assume that the AOD is in the general output...
582     fAOD  = AODEvent();
583     if(!fAOD){
584       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
585       return;
586     }
587     if(fDebug>0){
588       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
589     }
590   }
591   
592   Bool_t selectEvent =  false;
593   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
594
595   if(fAOD){
596     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
597     TString vtxTitle(vtxAOD->GetTitle());
598
599     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
600         Float_t zvtx = vtxAOD->GetZ();
601         Float_t yvtx = vtxAOD->GetY();
602         Float_t xvtx = vtxAOD->GetX();
603         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
604         if(TMath::Abs(zvtx)<8.&&r2<1.){
605           if(physicsSelection)selectEvent = true;
606         }
607     }
608   }
609   if(!selectEvent){
610     PostData(1, fHistList);
611     return;
612   }
613   
614   fh1Trials->Fill("#sum{ntrials}",1);
615   
616
617   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
618
619   // ==== General variables needed
620
621
622
623   // we simply fetch the tracks/mc particles as a list of AliVParticles
624
625   TList recParticles;
626   TList genParticles;
627
628   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
629   Float_t nCh = recParticles.GetEntries(); 
630   fh1Nch->Fill(nCh);
631   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
632   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
633   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
634
635   // find the jets....
636
637   vector<fastjet::PseudoJet> inputParticlesRec;
638   vector<fastjet::PseudoJet> inputParticlesRecRan;
639
640
641   // create a random jet within the acceptance
642   const float rRandomCone2 = 0.4*0.4;
643   float etaRandomCone = gRandom->Rndm()-0.5; // +- 0.5
644   float phiRandomCone = gRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
645   float ptRandomCone = 0;
646   float ptRandomConeRan = 0;
647
648   for(int i = 0; i < recParticles.GetEntries(); i++){
649     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
650     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
651     // we talk total momentum here
652     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
653     jInp.set_user_index(i);
654     inputParticlesRec.push_back(jInp);
655
656     if(evBkg){
657       float deta = vp->Eta()-etaRandomCone;
658       float dphi = vp->Phi()-phiRandomCone;
659       if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
660       if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
661       float dr2 = dphi*dphi+deta*deta;
662       if(dr2<rRandomCone2){
663         // within random jet
664         ptRandomCone += vp->Pt();
665       }
666     }
667
668     // the randomized input changes eta and phi, but keeps the p_T
669     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
670       Double_t pT = vp->Pt();
671       Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
672       Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
673       
674       Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));  
675       Double_t pZ = pT/TMath::Tan(theta);
676
677       Double_t pX = pT * TMath::Cos(phi);
678       Double_t pY = pT * TMath::Sin(phi);
679       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
680       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
681       jInpRan.set_user_index(i);
682       inputParticlesRecRan.push_back(jInpRan);
683
684       if(evBkg){
685         float deta = eta-etaRandomCone;
686         float dphi = phi-phiRandomCone;
687         if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
688         if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
689         float dr2 = dphi*dphi+deta*deta;
690         if(dr2<rRandomCone2){
691           // within random jet
692           ptRandomConeRan += pT;
693         }
694       }
695     }
696
697     // fill the tref array, only needed when we write out jets
698     if(jarray){
699       if(i == 0){
700         fRef->Delete(); // make sure to delete before placement new...
701         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
702       }
703       fRef->Add(vp);
704     }
705   }
706
707   if(inputParticlesRec.size()==0){
708     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
709     PostData(1, fHistList);
710     return;
711   }
712   
713   // run fast jet
714   // employ setters for these...
715
716  
717   // now create the object that holds info about ghosts                         
718   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
719   fastjet::AreaType areaType =   fastjet::active_area;
720   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
721   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
722   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
723   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
724   
725   //range where to compute background
726   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
727   phiMin = 0;
728   phiMax = 2*TMath::Pi();
729   rapMax = fGhostEtamax - fRparam;
730   rapMin = - fGhostEtamax + fRparam;
731   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
732  
733
734
735   inclusiveJets = clustSeq.inclusive_jets();
736   sortedJets    = sorted_by_pt(inclusiveJets);
737  
738   fh1NJetsRec->Fill(sortedJets.size());
739
740  // loop over all jets an fill information, first on is the leading jet
741
742   Int_t nRecOver = inclusiveJets.size();
743   Int_t nRec     = inclusiveJets.size();
744   if(inclusiveJets.size()>0){
745     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
746     Float_t pt = leadingJet.Pt();
747     Int_t nAodOutJets = 0;
748     Int_t nAodOutTracks = 0;
749     AliAODJet *aodOutJet = 0;
750
751     Int_t iCount = 0;
752     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
753       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
754       while(pt<ptCut&&iCount<nRec){
755         nRecOver--;
756         iCount++;
757         if(iCount<nRec){
758           pt = sortedJets[iCount].perp();
759         }
760       }
761       if(nRecOver<=0)break;
762       fh2NRecJetsPt->Fill(ptCut,nRecOver);
763     }
764     Float_t phi = leadingJet.Phi();
765     if(phi<0)phi+=TMath::Pi()*2.;    
766     Float_t eta = leadingJet.Eta();
767     pt = leadingJet.Pt();
768     
769     // correlation of leading jet with tracks
770     TIterator *recIter = recParticles.MakeIterator();
771     recIter->Reset();
772     AliVParticle *tmpRecTrack = 0;
773     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
774       Float_t tmpPt = tmpRecTrack->Pt();
775       // correlation
776       Float_t tmpPhi =  tmpRecTrack->Phi();     
777       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
778       Float_t dPhi = phi - tmpPhi;
779       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
780       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
781       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
782       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
783     }  
784     
785    
786     
787    for(int j = 0; j < nRec;j++){
788      AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
789      aodOutJet = 0;
790      nAodOutTracks = 0;
791      Float_t tmpPt = tmpRec.Pt();  
792      fh1PtJetsRecIn->Fill(tmpPt);
793      // Fill Spectra with constituents
794      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
795      fh1NConstRec->Fill(constituents.size());
796      fh2PtNch->Fill(nCh,tmpPt);
797      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
798      fh2NConstPt->Fill(tmpPt,constituents.size());
799      // loop over constiutents and fill spectrum
800
801      // Add the jet information and the track references to the output AOD
802      
803      if(tmpPt>fJetOutputMinPt&&jarray){
804        aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
805        Double_t area=clustSeq.area(sortedJets[j]);
806        
807        aodOutJet->SetEffArea(area,0);
808      }
809
810      
811      for(unsigned int ic = 0; ic < constituents.size();ic++){
812        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
813        fh1PtJetConstRec->Fill(part->Pt());
814        if(aodOutJet){
815          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
816        }
817        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
818      }
819
820
821
822      
823      
824
825
826      
827      // correlation
828      Float_t tmpPhi =  tmpRec.Phi();
829      Float_t tmpEta =  tmpRec.Eta();
830      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
831      
832      if(j==0){
833        fh1PtJetsLeadingRecIn->Fill(tmpPt);
834        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
835        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
836        fh1NConstLeadingRec->Fill(constituents.size());
837        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
838        continue;
839      }
840      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
841      fh2JetEtaPt->Fill(tmpEta,tmpPt);
842      Float_t dPhi = phi - tmpPhi;
843      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
844      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
845      Float_t dEta = eta - tmpRec.Eta();
846      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
847      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
848      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
849    }
850    delete recIter;
851
852       //background estimates:all bckg jets(0) & wo the 2 hardest(1)
853  
854    if(evBkg){
855      vector<fastjet::PseudoJet> jets2=sortedJets;
856      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
857      Double_t bkg1=0;
858      Double_t sigma1=0.;
859      Double_t meanarea1=0.;
860      Double_t bkg2=0;
861      Double_t sigma2=0.;
862      Double_t meanarea2=0.;
863
864      float areaRandomCone = rRandomCone2 *TMath::Pi();         
865      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
866      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
867      fh1BiARandomCones[0]->Fill(ptRandomCone-(bkg1*areaRandomCone));    
868      fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
869      
870      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
871      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
872      fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
873      fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
874    }
875   }
876   
877
878   
879   
880  
881   // fill track information
882   Int_t nTrackOver = recParticles.GetSize();
883   // do the same for tracks and jets
884
885   if(nTrackOver>0){
886    TIterator *recIter = recParticles.MakeIterator();
887    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
888    Float_t pt = tmpRec->Pt();
889
890    //    Printf("Leading track p_t %3.3E",pt);
891    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
892      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
893      while(pt<ptCut&&tmpRec){
894        nTrackOver--;
895        tmpRec = (AliVParticle*)(recIter->Next()); 
896        if(tmpRec){
897          pt = tmpRec->Pt();
898        }
899      }
900      if(nTrackOver<=0)break;
901      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
902    }
903    
904    recIter->Reset();
905    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
906    Float_t phi = leading->Phi();
907    if(phi<0)phi+=TMath::Pi()*2.;    
908    Float_t eta = leading->Eta();
909    pt  = leading->Pt();
910    while((tmpRec = (AliVParticle*)(recIter->Next()))){
911      Float_t tmpPt = tmpRec->Pt();
912      Float_t tmpEta = tmpRec->Eta();
913      fh1PtTracksRecIn->Fill(tmpPt);
914      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
915      if(tmpRec==leading){
916        fh1PtTracksLeadingRecIn->Fill(tmpPt);
917        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
918        continue;
919      }
920       // correlation
921      Float_t tmpPhi =  tmpRec->Phi();
922      
923      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
924      Float_t dPhi = phi - tmpPhi;
925      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
926      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
927      Float_t dEta = eta - tmpRec->Eta();
928      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
929      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
930      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
931    }  
932    delete recIter;
933  }
934
935  // find the random jets
936  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
937  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
938   
939  inclusiveJetsRan = clustSeqRan.inclusive_jets();
940  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
941
942  // fill the jet information from random track
943
944   fh1NJetsRecRan->Fill(sortedJetsRan.size());
945
946  // loop over all jets an fill information, first on is the leading jet
947
948  Int_t nRecOverRan = inclusiveJetsRan.size();
949  Int_t nRecRan     = inclusiveJetsRan.size();
950  if(inclusiveJetsRan.size()>0){
951    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
952    Float_t pt = leadingJet.Pt();
953    
954    Int_t iCount = 0;
955    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
956      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
957       while(pt<ptCut&&iCount<nRecRan){
958         nRecOverRan--;
959         iCount++;
960         if(iCount<nRecRan){
961           pt = sortedJetsRan[iCount].perp();
962         }
963       }
964       if(nRecOverRan<=0)break;
965       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
966     }
967     Float_t phi = leadingJet.Phi();
968     if(phi<0)phi+=TMath::Pi()*2.;    
969     pt = leadingJet.Pt();
970
971     // correlation of leading jet with random tracks
972
973     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
974       { 
975         Float_t tmpPt = inputParticlesRecRan[ip].perp();
976         // correlation
977         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
978         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
979         Float_t dPhi = phi - tmpPhi;
980         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
981         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
982         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
983         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
984       }  
985
986     Int_t nAodOutJetsRan = 0;
987      AliAODJet *aodOutJetRan = 0;
988     for(int j = 0; j < nRecRan;j++){
989       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
990       Float_t tmpPt = tmpRec.Pt();
991       fh1PtJetsRecInRan->Fill(tmpPt);
992       // Fill Spectra with constituents
993       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
994       fh1NConstRecRan->Fill(constituents.size());
995       fh2NConstPtRan->Fill(tmpPt,constituents.size());
996       fh2PtNchRan->Fill(nCh,tmpPt);
997       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
998
999
1000      if(tmpPt>fJetOutputMinPt&&jarrayran){
1001        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1002        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1003        
1004        aodOutJetRan->SetEffArea(arearan,0);    }
1005
1006
1007
1008      
1009      for(unsigned int ic = 0; ic < constituents.size();ic++){
1010        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1011        // fh1PtJetConstRec->Fill(part->Pt());
1012        if(aodOutJetRan){
1013          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1014        }}
1015       
1016
1017
1018
1019       // correlation
1020       Float_t tmpPhi =  tmpRec.Phi();
1021       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1022
1023       if(j==0){
1024         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1025         fh1NConstLeadingRecRan->Fill(constituents.size());
1026         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1027         continue;
1028       }
1029     }  
1030
1031      
1032     if(evBkg){
1033      Double_t bkg3=0.;
1034      Double_t sigma3=0.;
1035      Double_t meanarea3=0.;
1036      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1037      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1038      float areaRandomCone = rRandomCone2 *TMath::Pi();         
1039      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1040      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1041     }
1042
1043
1044
1045  }
1046  
1047
1048  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1049  PostData(1, fHistList);
1050 }
1051
1052 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1053 {
1054 // Terminate analysis
1055 //
1056     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1057 }
1058
1059
1060 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1061
1062   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1063
1064   Int_t iCount = 0;
1065   if(type==kTrackAOD){
1066     AliAODEvent *aod = 0;
1067     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1068     else aod = AODEvent();
1069     if(!aod){
1070       return iCount;
1071     }
1072     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1073       AliAODTrack *tr = aod->GetTrack(it);
1074       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1075       if(TMath::Abs(tr->Eta())>0.9)continue;
1076       //      if(tr->Pt()<0.3)continue;
1077       list->Add(tr);
1078       iCount++;
1079     }
1080   }
1081   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1082     AliMCEvent* mcEvent = MCEvent();
1083     if(!mcEvent)return iCount;
1084     // we want to have alivpartilces so use get track
1085     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1086       if(!mcEvent->IsPhysicalPrimary(it))continue;
1087       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1088       if(type == kTrackKineAll){
1089         list->Add(part);
1090         iCount++;
1091       }
1092       else if(type == kTrackKineCharged){
1093         if(part->Particle()->GetPDG()->Charge()==0)continue;
1094         list->Add(part);
1095         iCount++;
1096       }
1097     }
1098   }
1099   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1100     AliAODEvent *aod = 0;
1101     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1102     else aod = AODEvent();
1103     if(!aod)return iCount;
1104     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1105     if(!tca)return iCount;
1106     for(int it = 0;it < tca->GetEntriesFast();++it){
1107       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1108       if(!part->IsPhysicalPrimary())continue;
1109       if(type == kTrackAODMCAll){
1110         list->Add(part);
1111         iCount++;
1112       }
1113       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1114         if(part->Charge()==0)continue;
1115         if(kTrackAODMCCharged){
1116           list->Add(part);
1117         }
1118         else {
1119           if(TMath::Abs(part->Eta())>0.9)continue;
1120           list->Add(part);
1121         }
1122         iCount++;
1123       }
1124     }
1125   }// AODMCparticle
1126   list->Sort();
1127   return iCount;
1128
1129 }
1130
1131 /*
1132 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1133   for(int i = 0; i < particles.GetEntries(); i++){
1134     AliVParticle *vp = (AliVParticle*)particles.At(i);
1135     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1136     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1137     jInp.set_user_index(i);
1138     inputParticles.push_back(jInp);
1139   }
1140
1141   return 0;
1142
1143 }
1144 */