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