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