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