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