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