Adding task to analyze properties of jet like clusters, using directly the fastjet...
[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 <TRandom.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 "AliJetReaderHeader.h"
46 #include "AliUA1JetHeaderV1.h"
47 #include "AliESDEvent.h"
48 #include "AliAODEvent.h"
49 #include "AliAODHandler.h"
50 #include "AliAODTrack.h"
51 #include "AliAODJet.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliStack.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliJetKineReaderHeader.h"
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliInputEventHandler.h"
60
61
62 #include "fastjet/PseudoJet.hh"
63 #include "fastjet/ClusterSequenceArea.hh"
64 #include "fastjet/AreaDefinition.hh"
65 #include "fastjet/JetDefinition.hh"
66 // get info on how fastjet was configured
67 #include "fastjet/config.h"
68
69
70 ClassImp(AliAnalysisTaskJetCluster)
71
72 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
73   fAOD(0x0),
74   fUseAODTrackInput(kFALSE),
75   fUseAODMCInput(kFALSE),
76   fUseGlobalSelection(kFALSE),
77   fFilterMask(0),
78   fTrackTypeRec(kTrackUndef),
79   fTrackTypeGen(kTrackUndef),
80   fAvgTrials(1),
81   fExternalWeight(1),    
82   fRecEtaWindow(0.5),
83   fRparam(1.0), 
84   fAlgorithm(fastjet::kt_algorithm),
85   fStrategy(fastjet::Best),
86   fRecombScheme(fastjet::BIpt_scheme),
87   fAreaType(fastjet::active_area), 
88   fh1Xsec(0x0),
89   fh1Trials(0x0),
90   fh1PtHard(0x0),
91   fh1PtHardNoW(0x0),  
92   fh1PtHardTrials(0x0),
93   fh1NJetsRec(0x0),
94   fh1NConstRec(0x0),
95   fh1NConstLeadingRec(0x0),
96   fh1PtJetsRecIn(0x0),
97   fh1PtJetsLeadingRecIn(0x0),
98   fh1PtJetConstRec(0x0),
99   fh1PtJetConstLeadingRec(0x0),
100   fh1PtTracksRecIn(0x0),
101   fh1PtTracksLeadingRecIn(0x0),
102   fh1NJetsRecRan(0x0),
103   fh1NConstRecRan(0x0),
104   fh1PtJetsLeadingRecInRan(0x0),
105   fh1NConstLeadingRecRan(0x0),
106   fh1PtJetsRecInRan(0x0),
107   fh1PtTracksGenIn(0x0),
108   fh2NRecJetsPt(0x0),
109   fh2NRecTracksPt(0x0),
110   fh2NConstPt(0x0),
111   fh2NConstLeadingPt(0x0),
112   fh2JetPhiEta(0x0),
113   fh2LeadingJetPhiEta(0x0),
114   fh2JetEtaPt(0x0),
115   fh2LeadingJetEtaPt(0x0),
116   fh2TrackEtaPt(0x0),
117   fh2LeadingTrackEtaPt(0x0),
118   fh2JetsLeadingPhiEta(0x0),
119   fh2JetsLeadingPhiPt(0x0),
120   fh2TracksLeadingPhiEta(0x0),
121   fh2TracksLeadingPhiPt(0x0),
122   fh2TracksLeadingJetPhiPt(0x0),
123   fh2JetsLeadingPhiPtW(0x0),
124   fh2TracksLeadingPhiPtW(0x0),
125   fh2TracksLeadingJetPhiPtW(0x0),
126   fh2NRecJetsPtRan(0x0),
127   fh2NConstPtRan(0x0),
128   fh2NConstLeadingPtRan(0x0),
129   fh2TracksLeadingJetPhiPtRan(0x0),
130   fh2TracksLeadingJetPhiPtWRan(0x0),
131   fHistList(0x0)  
132 {
133
134 }
135
136 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
137   AliAnalysisTaskSE(name),
138   fAOD(0x0),
139   fUseAODTrackInput(kFALSE),
140   fUseAODMCInput(kFALSE),
141   fUseGlobalSelection(kFALSE),
142   fFilterMask(0),
143   fTrackTypeRec(kTrackUndef),
144   fTrackTypeGen(kTrackUndef),
145   fAvgTrials(1),
146   fExternalWeight(1),    
147   fRecEtaWindow(0.5),
148   fRparam(1.0), 
149   fAlgorithm(fastjet::kt_algorithm),
150   fStrategy(fastjet::Best),
151   fRecombScheme(fastjet::BIpt_scheme),
152   fAreaType(fastjet::active_area), 
153   fh1Xsec(0x0),
154   fh1Trials(0x0),
155   fh1PtHard(0x0),
156   fh1PtHardNoW(0x0),  
157   fh1PtHardTrials(0x0),
158   fh1NJetsRec(0x0),
159   fh1NConstRec(0x0),
160   fh1NConstLeadingRec(0x0),
161   fh1PtJetsRecIn(0x0),
162   fh1PtJetsLeadingRecIn(0x0),
163   fh1PtJetConstRec(0x0),
164   fh1PtJetConstLeadingRec(0x0),
165   fh1PtTracksRecIn(0x0),
166   fh1PtTracksLeadingRecIn(0x0),
167   fh1NJetsRecRan(0x0),
168   fh1NConstRecRan(0x0),
169   fh1PtJetsLeadingRecInRan(0x0),
170   fh1NConstLeadingRecRan(0x0),
171   fh1PtJetsRecInRan(0x0),
172   fh1PtTracksGenIn(0x0),
173   fh2NRecJetsPt(0x0),
174   fh2NRecTracksPt(0x0),
175   fh2NConstPt(0x0),
176   fh2NConstLeadingPt(0x0),
177   fh2JetPhiEta(0x0),
178   fh2LeadingJetPhiEta(0x0),
179   fh2JetEtaPt(0x0),
180   fh2LeadingJetEtaPt(0x0),
181   fh2TrackEtaPt(0x0),
182   fh2LeadingTrackEtaPt(0x0),
183   fh2JetsLeadingPhiEta(0x0),
184   fh2JetsLeadingPhiPt(0x0),
185   fh2TracksLeadingPhiEta(0x0),
186   fh2TracksLeadingPhiPt(0x0),
187   fh2TracksLeadingJetPhiPt(0x0),
188   fh2JetsLeadingPhiPtW(0x0),
189   fh2TracksLeadingPhiPtW(0x0),
190   fh2TracksLeadingJetPhiPtW(0x0),
191   fh2NRecJetsPtRan(0x0),
192   fh2NConstPtRan(0x0),
193   fh2NConstLeadingPtRan(0x0),
194   fh2TracksLeadingJetPhiPtRan(0x0),
195   fh2TracksLeadingJetPhiPtWRan(0x0),
196   fHistList(0x0)
197 {
198   DefineOutput(1, TList::Class());  
199 }
200
201
202
203 Bool_t AliAnalysisTaskJetCluster::Notify()
204 {
205   //
206   // Implemented Notify() to read the cross sections
207   // and number of trials from pyxsec.root
208   // 
209   return kTRUE;
210 }
211
212 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
213 {
214
215   //
216   // Create the output container
217   //
218
219
220   // Connect the AOD
221
222
223   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
224
225   OpenFile(1);
226   if(!fHistList)fHistList = new TList();
227
228   Bool_t oldStatus = TH1::AddDirectoryStatus();
229   TH1::AddDirectory(kFALSE);
230
231   //
232   //  Histogram
233     
234   const Int_t nBinPt = 200;
235   Double_t binLimitsPt[nBinPt+1];
236   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
237     if(iPt == 0){
238       binLimitsPt[iPt] = 0.0;
239     }
240     else {// 1.0
241       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 0.5;
242     }
243   }
244   
245   const Int_t nBinPhi = 90;
246   Double_t binLimitsPhi[nBinPhi+1];
247   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
248     if(iPhi==0){
249       binLimitsPhi[iPhi] = -1.*TMath::Pi();
250     }
251     else{
252       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
253     }
254   }
255
256
257
258   const Int_t nBinEta = 40;
259   Double_t binLimitsEta[nBinEta+1];
260   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
261     if(iEta==0){
262       binLimitsEta[iEta] = -2.0;
263     }
264     else{
265       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
266     }
267   }
268
269
270   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
271   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
272
273   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
274   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
275
276
277   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
278   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
279
280   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
281   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
282   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
283   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
284
285
286   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
287   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
288   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
289
290   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
291   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
292   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
293   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
294   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
295   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
296   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
297   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
298   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
299
300   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
301   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
302   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
303   // 
304
305
306   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
307   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
308   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
309   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
310
311   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
312                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
313   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
314                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
315
316   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
317                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
318   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
319                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
320
321   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
322                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
323   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
324                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
325
326
327
328   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
329                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
330   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
331                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
332   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
333                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
334   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
335                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
336   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
337                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
338   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
339                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
340
341   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
342                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
343   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
344                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
345
346   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
347                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
348   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
349                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
350
351
352
353   const Int_t saveLevel = 3; // large save level more histos
354   if(saveLevel>0){
355     fHistList->Add(fh1Xsec);
356     fHistList->Add(fh1Trials);
357
358     fHistList->Add(fh1NJetsRec);
359     fHistList->Add(fh1NConstRec);
360     fHistList->Add(fh1NConstLeadingRec);
361     fHistList->Add(fh1PtJetsRecIn);
362     fHistList->Add(fh1PtJetsLeadingRecIn);
363     fHistList->Add(fh1PtTracksRecIn);
364     fHistList->Add(fh1PtTracksLeadingRecIn);
365     fHistList->Add(fh1PtJetConstRec);
366     fHistList->Add(fh1PtJetConstLeadingRec);
367     fHistList->Add(fh1NJetsRecRan);
368     fHistList->Add(fh1NConstRecRan);
369     fHistList->Add(fh1PtJetsLeadingRecInRan);
370     fHistList->Add(fh1NConstLeadingRecRan);
371     fHistList->Add(fh1PtJetsRecInRan);
372     fHistList->Add(fh2NRecJetsPt);
373     fHistList->Add(fh2NRecTracksPt);
374     fHistList->Add(fh2NConstPt);
375     fHistList->Add(fh2NConstLeadingPt);
376     fHistList->Add(fh2JetPhiEta);
377     fHistList->Add(fh2LeadingJetPhiEta);
378     fHistList->Add(fh2JetEtaPt);
379     fHistList->Add(fh2LeadingJetEtaPt);
380     fHistList->Add(fh2TrackEtaPt);
381     fHistList->Add(fh2LeadingTrackEtaPt);
382     fHistList->Add(fh2JetsLeadingPhiEta );
383     fHistList->Add(fh2JetsLeadingPhiPt);
384     fHistList->Add(fh2TracksLeadingPhiEta);
385     fHistList->Add(fh2TracksLeadingPhiPt);
386     fHistList->Add(fh2TracksLeadingJetPhiPt);
387     fHistList->Add(fh2JetsLeadingPhiPtW);
388     fHistList->Add(fh2TracksLeadingPhiPtW);
389     fHistList->Add(fh2TracksLeadingJetPhiPtW);
390     fHistList->Add(fh2NRecJetsPtRan);
391     fHistList->Add(fh2NConstPtRan);
392     fHistList->Add(fh2NConstLeadingPtRan);
393     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
394     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
395     }
396
397   // =========== Switch on Sumw2 for all histos ===========
398   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
399     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
400     if (h1){
401       h1->Sumw2();
402       continue;
403     }
404     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
405     if(hn)hn->Sumw2();
406   }
407   TH1::AddDirectory(oldStatus);
408 }
409
410 void AliAnalysisTaskJetCluster::Init()
411 {
412   //
413   // Initialization
414   //
415
416   Printf(">>> AnalysisTaskJetCluster::Init() debug level %d\n",fDebug);
417   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
418
419 }
420
421 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
422 {
423
424   if(fUseGlobalSelection){
425     // no selection by the service task, we continue
426     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
427     PostData(1, fHistList);
428     return;
429   }
430
431   //
432   // Execute analysis for current event
433   //
434   AliESDEvent *fESD = 0;
435   if(fUseAODTrackInput){    
436     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
437     if(!fAOD){
438       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
439       return;
440     }
441     // fethc the header
442   }
443   else{
444     //  assume that the AOD is in the general output...
445     fAOD  = AODEvent();
446     if(!fAOD){
447       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
448       return;
449     }
450     if(fDebug>0){
451       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
452     }
453   }
454   
455   Bool_t selectEvent =  false;
456   Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
457
458   if(fAOD){
459     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
460     if(vtxAOD->GetNContributors()>0){
461         Float_t zvtx = vtxAOD->GetZ();
462         Float_t yvtx = vtxAOD->GetY();
463         Float_t xvtx = vtxAOD->GetX();
464         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
465         if(TMath::Abs(zvtx)<8.&&r2<1.){
466           if(physicsSelection)selectEvent = true;
467         }
468     }
469   }
470   if(!selectEvent){
471     PostData(1, fHistList);
472     return;
473   }
474   
475   fh1Trials->Fill("#sum{ntrials}",1);
476   
477
478   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
479
480   // ==== General variables needed
481
482
483
484   // we simply fetch the tracks/mc particles as a list of AliVParticles
485
486   TList recParticles;
487   TList genParticles;
488
489   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
490   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
491   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
492   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
493
494   // find the jets....
495
496   vector<fastjet::PseudoJet> inputParticlesRec;
497   vector<fastjet::PseudoJet> inputParticlesRecRan;
498   for(int i = 0; i < recParticles.GetEntries(); i++){
499     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
500     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
501     // we talk total momentum here
502     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
503     jInp.set_user_index(i);
504     inputParticlesRec.push_back(jInp);
505
506     // the randomized input changes eta and phi, but keeps the p_T
507     Double_t pT = vp->Pt();
508     Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
509     Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
510
511
512     Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));  
513     Double_t pZ = pT/TMath::Tan(theta);
514
515     Double_t pX = pT * TMath::Cos(phi);
516     Double_t pY = pT * TMath::Sin(phi);
517     Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
518
519     fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
520     jInpRan.set_user_index(i);
521     inputParticlesRecRan.push_back(jInpRan);
522
523   }
524
525   if(inputParticlesRec.size()==0){
526     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
527     PostData(1, fHistList);
528     return;
529   }
530   
531   // run fast jet
532   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
533   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
534   fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
535   
536   inclusiveJets = clustSeq.inclusive_jets();
537   sortedJets    = sorted_by_pt(inclusiveJets);
538  
539   if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
540
541   fh1NJetsRec->Fill(sortedJets.size());
542
543  // loop over all jets an fill information, first on is the leading jet
544
545  Int_t nRecOver = inclusiveJets.size();
546  Int_t nRec     = inclusiveJets.size();
547  if(inclusiveJets.size()>0){
548    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
549     Float_t pt = leadingJet.Pt();
550
551     Int_t iCount = 0;
552     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
553       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
554       while(pt<ptCut&&iCount<nRec){
555         nRecOver--;
556         iCount++;
557         if(iCount<nRec){
558           pt = sortedJets[iCount].perp();
559         }
560       }
561       if(nRecOver<=0)break;
562       fh2NRecJetsPt->Fill(ptCut,nRecOver);
563     }
564     Float_t phi = leadingJet.Phi();
565     if(phi<0)phi+=TMath::Pi()*2.;    
566     Float_t eta = leadingJet.Eta();
567     pt = leadingJet.Pt();
568
569     // correlation of leading jet with tracks
570     TIterator *recIter = recParticles.MakeIterator();
571     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());  
572
573     recIter->Reset();
574     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
575       Float_t tmpPt = tmpRecTrack->Pt();
576       // correlation
577       Float_t tmpPhi =  tmpRecTrack->Phi();     
578       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
579       Float_t dPhi = phi - tmpPhi;
580       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
581       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
582       //      Float_t dEta = eta - tmpRecTrack->Eta();
583       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
584       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
585    }  
586
587
588
589     for(int j = 0; j < nRec;j++){
590       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
591       Float_t tmpPt = tmpRec.Pt();
592       fh1PtJetsRecIn->Fill(tmpPt);
593       // Fill Spectra with constituents
594       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
595       fh1NConstRec->Fill(constituents.size());
596       fh2NConstPt->Fill(tmpPt,constituents.size());
597       // loop over constiutents and fill spectrum
598       for(unsigned int ic = 0; ic < constituents.size();ic++){
599         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
600         fh1PtJetConstRec->Fill(part->Pt());
601         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
602       }
603
604       // correlation
605       Float_t tmpPhi =  tmpRec.Phi();
606       Float_t tmpEta =  tmpRec.Eta();
607       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
608
609       if(j==0){
610         fh1PtJetsLeadingRecIn->Fill(tmpPt);
611         fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
612         fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
613         fh1NConstLeadingRec->Fill(constituents.size());
614         fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
615         continue;
616       }
617       fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
618       fh2JetEtaPt->Fill(tmpEta,tmpPt);
619       Float_t dPhi = phi - tmpPhi;
620       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
621       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
622       Float_t dEta = eta - tmpRec.Eta();
623       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
624       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
625       fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
626     }  
627  }
628  
629
630  // fill track information
631  Int_t nTrackOver = recParticles.GetSize();
632   // do the same for tracks and jets
633  if(nTrackOver>0){
634    TIterator *recIter = recParticles.MakeIterator();
635    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
636    Float_t pt = tmpRec->Pt();
637    //    Printf("Leading track p_t %3.3E",pt);
638    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
639      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
640      while(pt<ptCut&&tmpRec){
641        nTrackOver--;
642        tmpRec = (AliVParticle*)(recIter->Next()); 
643        if(tmpRec){
644          pt = tmpRec->Pt();
645        }
646      }
647      if(nTrackOver<=0)break;
648      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
649    }
650    
651    recIter->Reset();
652    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
653    Float_t phi = leading->Phi();
654    if(phi<0)phi+=TMath::Pi()*2.;    
655    Float_t eta = leading->Eta();
656    pt  = leading->Pt();
657    while((tmpRec = (AliVParticle*)(recIter->Next()))){
658      Float_t tmpPt = tmpRec->Pt();
659      Float_t tmpEta = tmpRec->Eta();
660      fh1PtTracksRecIn->Fill(tmpPt);
661      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
662      if(tmpRec==leading){
663        fh1PtTracksLeadingRecIn->Fill(tmpPt);
664        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
665        continue;
666      }
667       // correlation
668      Float_t tmpPhi =  tmpRec->Phi();
669      
670      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
671      Float_t dPhi = phi - tmpPhi;
672      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
673      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
674      Float_t dEta = eta - tmpRec->Eta();
675      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
676      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
677      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
678    }  
679    delete recIter;
680  }
681
682  // find the random jets
683  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
684  fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
685   
686  inclusiveJetsRan = clustSeqRan.inclusive_jets();
687  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
688
689  // fill the jet information from random track
690
691   fh1NJetsRecRan->Fill(sortedJetsRan.size());
692
693  // loop over all jets an fill information, first on is the leading jet
694
695  Int_t nRecOverRan = inclusiveJetsRan.size();
696  Int_t nRecRan     = inclusiveJetsRan.size();
697  if(inclusiveJetsRan.size()>0){
698    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
699    Float_t pt = leadingJet.Pt();
700    
701    Int_t iCount = 0;
702    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
703      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
704       while(pt<ptCut&&iCount<nRecRan){
705         nRecOverRan--;
706         iCount++;
707         if(iCount<nRecRan){
708           pt = sortedJetsRan[iCount].perp();
709         }
710       }
711       if(nRecOverRan<=0)break;
712       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
713     }
714     Float_t phi = leadingJet.Phi();
715     if(phi<0)phi+=TMath::Pi()*2.;    
716     pt = leadingJet.Pt();
717
718     // correlation of leading jet with random tracks
719
720     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
721       { 
722         Float_t tmpPt = inputParticlesRecRan[ip].perp();
723         // correlation
724         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
725         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
726         Float_t dPhi = phi - tmpPhi;
727         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
728         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
729         //      Float_t dEta = eta - tmpRecTrack->Eta();
730         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
731         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
732       }  
733
734
735
736     for(int j = 0; j < nRecRan;j++){
737       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
738       Float_t tmpPt = tmpRec.Pt();
739       fh1PtJetsRecInRan->Fill(tmpPt);
740       // Fill Spectra with constituents
741       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
742       fh1NConstRecRan->Fill(constituents.size());
743       fh2NConstPtRan->Fill(tmpPt,constituents.size());
744
745       // correlation
746       Float_t tmpPhi =  tmpRec.Phi();
747       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
748
749       if(j==0){
750         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
751         fh1NConstLeadingRecRan->Fill(constituents.size());
752         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
753         continue;
754       }
755     }  
756  }
757  
758
759  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
760  PostData(1, fHistList);
761 }
762
763 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
764 {
765 // Terminate analysis
766 //
767     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
768 }
769
770
771 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
772
773   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
774
775   Int_t iCount = 0;
776   if(type==kTrackAOD){
777     AliAODEvent *aod = 0;
778     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
779     else aod = AODEvent();
780     if(!aod){
781       return iCount;
782     }
783     for(int it = 0;it < aod->GetNumberOfTracks();++it){
784       AliAODTrack *tr = aod->GetTrack(it);
785       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
786       if(TMath::Abs(tr->Eta())>0.9)continue;
787       //      if(tr->Pt()<0.3)continue;
788       list->Add(tr);
789       iCount++;
790     }
791   }
792   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
793     AliMCEvent* mcEvent = MCEvent();
794     if(!mcEvent)return iCount;
795     // we want to have alivpartilces so use get track
796     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
797       if(!mcEvent->IsPhysicalPrimary(it))continue;
798       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
799       if(type == kTrackKineAll){
800         list->Add(part);
801         iCount++;
802       }
803       else if(type == kTrackKineCharged){
804         if(part->Particle()->GetPDG()->Charge()==0)continue;
805         list->Add(part);
806         iCount++;
807       }
808     }
809   }
810   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
811     AliAODEvent *aod = 0;
812     if(fUseAODMCInput)dynamic_cast<AliAODEvent*>(InputEvent());
813     else aod = AODEvent();
814     if(!aod)return iCount;
815     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
816     if(!tca)return iCount;
817     for(int it = 0;it < tca->GetEntriesFast();++it){
818       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
819       if(!part->IsPhysicalPrimary())continue;
820       if(type == kTrackAODMCAll){
821         list->Add(part);
822         iCount++;
823       }
824       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
825         if(part->Charge()==0)continue;
826         if(kTrackAODMCCharged){
827           list->Add(part);
828         }
829         else {
830           if(TMath::Abs(part->Eta())>0.9)continue;
831           list->Add(part);
832         }
833         iCount++;
834       }
835     }
836   }// AODMCparticle
837   list->Sort();
838   return iCount;
839
840 }
841
842 /*
843 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
844   for(int i = 0; i < particles.GetEntries(); i++){
845     AliVParticle *vp = (AliVParticle*)particles.At(i);
846     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
847     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
848     jInp.set_user_index(i);
849     inputParticles.push_back(jInp);
850   }
851
852   return 0;
853
854 }
855 */