removing print out in Init()
[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   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
417
418 }
419
420 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
421 {
422
423   if(fUseGlobalSelection){
424     // no selection by the service task, we continue
425     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
426     PostData(1, fHistList);
427     return;
428   }
429
430   //
431   // Execute analysis for current event
432   //
433   AliESDEvent *fESD = 0;
434   if(fUseAODTrackInput){    
435     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
436     if(!fAOD){
437       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
438       return;
439     }
440     // fethc the header
441   }
442   else{
443     //  assume that the AOD is in the general output...
444     fAOD  = AODEvent();
445     if(!fAOD){
446       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
447       return;
448     }
449     if(fDebug>0){
450       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
451     }
452   }
453   
454   Bool_t selectEvent =  false;
455   Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
456
457   if(fAOD){
458     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
459     if(vtxAOD->GetNContributors()>0){
460         Float_t zvtx = vtxAOD->GetZ();
461         Float_t yvtx = vtxAOD->GetY();
462         Float_t xvtx = vtxAOD->GetX();
463         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
464         if(TMath::Abs(zvtx)<8.&&r2<1.){
465           if(physicsSelection)selectEvent = true;
466         }
467     }
468   }
469   if(!selectEvent){
470     PostData(1, fHistList);
471     return;
472   }
473   
474   fh1Trials->Fill("#sum{ntrials}",1);
475   
476
477   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
478
479   // ==== General variables needed
480
481
482
483   // we simply fetch the tracks/mc particles as a list of AliVParticles
484
485   TList recParticles;
486   TList genParticles;
487
488   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
489   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
490   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
491   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
492
493   // find the jets....
494
495   vector<fastjet::PseudoJet> inputParticlesRec;
496   vector<fastjet::PseudoJet> inputParticlesRecRan;
497   for(int i = 0; i < recParticles.GetEntries(); i++){
498     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
499     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
500     // we talk total momentum here
501     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
502     jInp.set_user_index(i);
503     inputParticlesRec.push_back(jInp);
504
505     // the randomized input changes eta and phi, but keeps the p_T
506     Double_t pT = vp->Pt();
507     Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
508     Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
509
510
511     Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));  
512     Double_t pZ = pT/TMath::Tan(theta);
513
514     Double_t pX = pT * TMath::Cos(phi);
515     Double_t pY = pT * TMath::Sin(phi);
516     Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
517
518     fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
519     jInpRan.set_user_index(i);
520     inputParticlesRecRan.push_back(jInpRan);
521
522   }
523
524   if(inputParticlesRec.size()==0){
525     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
526     PostData(1, fHistList);
527     return;
528   }
529   
530   // run fast jet
531   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
532   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
533   fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
534   
535   inclusiveJets = clustSeq.inclusive_jets();
536   sortedJets    = sorted_by_pt(inclusiveJets);
537  
538   if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
539
540   fh1NJetsRec->Fill(sortedJets.size());
541
542  // loop over all jets an fill information, first on is the leading jet
543
544  Int_t nRecOver = inclusiveJets.size();
545  Int_t nRec     = inclusiveJets.size();
546  if(inclusiveJets.size()>0){
547    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
548     Float_t pt = leadingJet.Pt();
549
550     Int_t iCount = 0;
551     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
552       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
553       while(pt<ptCut&&iCount<nRec){
554         nRecOver--;
555         iCount++;
556         if(iCount<nRec){
557           pt = sortedJets[iCount].perp();
558         }
559       }
560       if(nRecOver<=0)break;
561       fh2NRecJetsPt->Fill(ptCut,nRecOver);
562     }
563     Float_t phi = leadingJet.Phi();
564     if(phi<0)phi+=TMath::Pi()*2.;    
565     Float_t eta = leadingJet.Eta();
566     pt = leadingJet.Pt();
567
568     // correlation of leading jet with tracks
569     TIterator *recIter = recParticles.MakeIterator();
570     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());  
571
572     recIter->Reset();
573     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
574       Float_t tmpPt = tmpRecTrack->Pt();
575       // correlation
576       Float_t tmpPhi =  tmpRecTrack->Phi();     
577       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
578       Float_t dPhi = phi - tmpPhi;
579       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
580       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
581       //      Float_t dEta = eta - tmpRecTrack->Eta();
582       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
583       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
584    }  
585
586
587
588     for(int j = 0; j < nRec;j++){
589       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
590       Float_t tmpPt = tmpRec.Pt();
591       fh1PtJetsRecIn->Fill(tmpPt);
592       // Fill Spectra with constituents
593       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
594       fh1NConstRec->Fill(constituents.size());
595       fh2NConstPt->Fill(tmpPt,constituents.size());
596       // loop over constiutents and fill spectrum
597       for(unsigned int ic = 0; ic < constituents.size();ic++){
598         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
599         fh1PtJetConstRec->Fill(part->Pt());
600         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
601       }
602
603       // correlation
604       Float_t tmpPhi =  tmpRec.Phi();
605       Float_t tmpEta =  tmpRec.Eta();
606       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
607
608       if(j==0){
609         fh1PtJetsLeadingRecIn->Fill(tmpPt);
610         fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
611         fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
612         fh1NConstLeadingRec->Fill(constituents.size());
613         fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
614         continue;
615       }
616       fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
617       fh2JetEtaPt->Fill(tmpEta,tmpPt);
618       Float_t dPhi = phi - tmpPhi;
619       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
620       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
621       Float_t dEta = eta - tmpRec.Eta();
622       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
623       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
624       fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
625     }
626     delete recIter;
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 */