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