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