fixed reading from OADB (M. Verweij)
[u/mrichter/AliRoot.git] / JETAN / DEV / 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 <TF1.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include  "TDatabasePDG.h"
40
41 #include "AliAnalysisTaskJetCluster.h"
42 #include "AliAnalysisManager.h"
43 #include "AliJetFinder.h"
44 #include "AliJetHeader.h"
45 #include "AliJetReader.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODExtension.h"
50 #include "AliAODTrack.h"
51 #include "AliAODJet.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliStack.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliJetKineReaderHeader.h"
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliInputEventHandler.h"
60 #include "AliAODJetEventBackground.h"
61
62 #include "fastjet/PseudoJet.hh"
63 #include "fastjet/ClusterSequenceArea.hh"
64 #include "fastjet/AreaDefinition.hh"
65 #include "fastjet/JetDefinition.hh"
66 // get info on how fastjet was configured
67 #include "fastjet/config.h"
68
69 using std::vector;
70
71 ClassImp(AliAnalysisTaskJetCluster)
72
73 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
74   //
75   // Destructor
76   //
77
78   delete fRef;
79   delete fRandom;
80
81   if(fTCAJetsOut)fTCAJetsOut->Delete();
82   delete fTCAJetsOut;
83   
84   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
85   delete fTCAJetsOutRan;
86   
87   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
88   delete fTCARandomConesOut;
89   
90   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
91   delete fTCARandomConesOutRan;
92   
93   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
94   delete fAODJetBackgroundOut;
95 }
96
97 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): 
98   AliAnalysisTaskSE(),
99   fAOD(0x0),
100   fAODExtension(0x0),
101   fRef(new TRefArray),
102   fUseAODTrackInput(kFALSE),
103   fUseAODMCInput(kFALSE),
104   fUseBackgroundCalc(kFALSE),
105   fEventSelection(kFALSE),     
106   fFilterMask(0),
107   fFilterMaskBestPt(0),
108   fFilterType(0),
109   fJetTypes(kJet),
110   fTrackTypeRec(kTrackUndef),
111   fTrackTypeGen(kTrackUndef),  
112   fNSkipLeadingRan(0),
113   fNSkipLeadingCone(0),
114   fNRandomCones(0),
115   fAvgTrials(1),
116   fExternalWeight(1),
117   fTrackEtaWindow(0.9),    
118   fRecEtaWindow(0.5),
119   fTrackPtCut(0.),                                                      
120   fJetOutputMinPt(0.150),
121   fMaxTrackPtInJet(100.),
122   fJetTriggerPtCut(0),
123   fVtxZCut(8),
124   fVtxR2Cut(1),
125   fCentCutUp(0),
126   fCentCutLo(0),
127   fNonStdBranch(""),
128   fBackgroundBranch(""),
129   fNonStdFile(""),
130   fMomResH1(0x0),
131   fMomResH2(0x0),
132   fMomResH3(0x0),
133   fMomResH1Fit(0x0),
134   fMomResH2Fit(0x0),
135   fMomResH3Fit(0x0),
136   fhEffH1(0x0),
137   fhEffH2(0x0),
138   fhEffH3(0x0),
139   fUseTrPtResolutionSmearing(kFALSE),
140   fUseDiceEfficiency(kFALSE),
141   fDiceEfficiencyMinPt(0.),
142   fUseTrPtResolutionFromOADB(kFALSE),
143   fUseTrEfficiencyFromOADB(kFALSE),
144   fPathTrPtResolution(""),
145   fPathTrEfficiency(""),
146   fChangeEfficiencyFraction(0.),
147   fRparam(1.0), 
148   fAlgorithm(fastjet::kt_algorithm),
149   fStrategy(fastjet::Best),
150   fRecombScheme(fastjet::BIpt_scheme),
151   fAreaType(fastjet::active_area), 
152   fGhostArea(0.01),
153   fActiveAreaRepeats(1),
154   fGhostEtamax(1.5),
155   fTCAJetsOut(0x0),
156   fTCAJetsOutRan(0x0),
157   fTCARandomConesOut(0x0),
158   fTCARandomConesOutRan(0x0),
159   fAODJetBackgroundOut(0x0),
160   fRandom(0),
161   fh1Xsec(0x0),
162   fh1Trials(0x0),
163   fh1PtHard(0x0),
164   fh1PtHardNoW(0x0),  
165   fh1PtHardTrials(0x0),
166   fh1NJetsRec(0x0),
167   fh1NConstRec(0x0),
168   fh1NConstLeadingRec(0x0),
169   fh1PtJetsRecIn(0x0),
170   fh1PtJetsLeadingRecIn(0x0),
171   fh1PtJetConstRec(0x0),
172   fh1PtJetConstLeadingRec(0x0),
173   fh1PtTracksRecIn(0x0),
174   fh1PtTracksLeadingRecIn(0x0),
175   fh1NJetsRecRan(0x0),
176   fh1NConstRecRan(0x0),
177   fh1PtJetsLeadingRecInRan(0x0),
178   fh1NConstLeadingRecRan(0x0),
179   fh1PtJetsRecInRan(0x0),
180   fh1PtTracksGenIn(0x0),
181   fh1Nch(0x0),
182   fh1CentralityPhySel(0x0), 
183   fh1Centrality(0x0), 
184   fh1CentralitySelect(0x0),
185   fh1ZPhySel(0x0), 
186   fh1Z(0x0), 
187   fh1ZSelect(0x0),
188   fh2NRecJetsPt(0x0),
189   fh2NRecTracksPt(0x0),
190   fh2NConstPt(0x0),
191   fh2NConstLeadingPt(0x0),
192   fh2JetPhiEta(0x0),
193   fh2LeadingJetPhiEta(0x0),
194   fh2JetEtaPt(0x0),
195   fh2LeadingJetEtaPt(0x0),
196   fh2TrackEtaPt(0x0),
197   fh2LeadingTrackEtaPt(0x0),
198   fh2JetsLeadingPhiEta(0x0),
199   fh2JetsLeadingPhiPt(0x0),
200   fh2TracksLeadingPhiEta(0x0),
201   fh2TracksLeadingPhiPt(0x0),
202   fh2TracksLeadingJetPhiPt(0x0),
203   fh2JetsLeadingPhiPtW(0x0),
204   fh2TracksLeadingPhiPtW(0x0),
205   fh2TracksLeadingJetPhiPtW(0x0),
206   fh2NRecJetsPtRan(0x0),
207   fh2NConstPtRan(0x0),
208   fh2NConstLeadingPtRan(0x0),
209   fh2PtNch(0x0),
210   fh2PtNchRan(0x0),
211   fh2PtNchN(0x0),
212   fh2PtNchNRan(0x0),
213   fh2TracksLeadingJetPhiPtRan(0x0),
214   fh2TracksLeadingJetPhiPtWRan(0x0),
215   fh2PtGenPtSmeared(0x0),
216   fp1Efficiency(0x0),
217   fp1PtResolution(0x0),
218   fHistList(0x0)  
219 {
220   //
221   // Constructor
222   //
223
224   for(int i = 0;i<3;i++){
225     fh1BiARandomCones[i] = 0;
226     fh1BiARandomConesRan[i] = 0;
227   }
228   for(int i = 0;i<kMaxCent;i++){
229     fh2JetsLeadingPhiPtC[i] = 0;     
230     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
231     fh2TracksLeadingJetPhiPtC[i] = 0;
232     fh2TracksLeadingJetPhiPtWC[i] = 0;
233   }
234 }
235
236 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
237   AliAnalysisTaskSE(name),
238   fAOD(0x0),
239   fAODExtension(0x0),
240   fRef(new TRefArray),
241   fUseAODTrackInput(kFALSE),
242   fUseAODMCInput(kFALSE),
243   fUseBackgroundCalc(kFALSE),
244   fEventSelection(kFALSE),                                                        fFilterMask(0),
245   fFilterMaskBestPt(0),
246   fFilterType(0),
247   fJetTypes(kJet),
248   fTrackTypeRec(kTrackUndef),
249   fTrackTypeGen(kTrackUndef),
250   fNSkipLeadingRan(0),
251   fNSkipLeadingCone(0),
252   fNRandomCones(0),
253   fAvgTrials(1),
254   fExternalWeight(1),    
255   fTrackEtaWindow(0.9),    
256   fRecEtaWindow(0.5),
257   fTrackPtCut(0.),                                                      
258   fJetOutputMinPt(0.150),
259   fMaxTrackPtInJet(100.),
260   fJetTriggerPtCut(0),
261   fVtxZCut(8),
262   fVtxR2Cut(1),
263   fCentCutUp(0),
264   fCentCutLo(0),
265   fNonStdBranch(""),
266   fBackgroundBranch(""),
267   fNonStdFile(""),
268   fMomResH1(0x0),
269   fMomResH2(0x0),
270   fMomResH3(0x0),
271   fMomResH1Fit(0x0),
272   fMomResH2Fit(0x0),
273   fMomResH3Fit(0x0),
274   fhEffH1(0x0),
275   fhEffH2(0x0),
276   fhEffH3(0x0),
277   fUseTrPtResolutionSmearing(kFALSE),
278   fUseDiceEfficiency(kFALSE),
279   fDiceEfficiencyMinPt(0.),
280   fUseTrPtResolutionFromOADB(kFALSE),
281   fUseTrEfficiencyFromOADB(kFALSE),
282   fPathTrPtResolution(""),
283   fPathTrEfficiency(""),
284   fChangeEfficiencyFraction(0.),
285   fRparam(1.0), 
286   fAlgorithm(fastjet::kt_algorithm),
287   fStrategy(fastjet::Best),
288   fRecombScheme(fastjet::BIpt_scheme),
289   fAreaType(fastjet::active_area), 
290   fGhostArea(0.01),
291   fActiveAreaRepeats(1),
292   fGhostEtamax(1.5),
293   fTCAJetsOut(0x0),
294   fTCAJetsOutRan(0x0),
295   fTCARandomConesOut(0x0),
296   fTCARandomConesOutRan(0x0),
297   fAODJetBackgroundOut(0x0),
298   fRandom(0),
299   fh1Xsec(0x0),
300   fh1Trials(0x0),
301   fh1PtHard(0x0),
302   fh1PtHardNoW(0x0),  
303   fh1PtHardTrials(0x0),
304   fh1NJetsRec(0x0),
305   fh1NConstRec(0x0),
306   fh1NConstLeadingRec(0x0),
307   fh1PtJetsRecIn(0x0),
308   fh1PtJetsLeadingRecIn(0x0),
309   fh1PtJetConstRec(0x0),
310   fh1PtJetConstLeadingRec(0x0),
311   fh1PtTracksRecIn(0x0),
312   fh1PtTracksLeadingRecIn(0x0),
313   fh1NJetsRecRan(0x0),
314   fh1NConstRecRan(0x0),
315   fh1PtJetsLeadingRecInRan(0x0),
316   fh1NConstLeadingRecRan(0x0),
317   fh1PtJetsRecInRan(0x0),
318   fh1PtTracksGenIn(0x0),
319   fh1Nch(0x0),
320   fh1CentralityPhySel(0x0), 
321   fh1Centrality(0x0), 
322   fh1CentralitySelect(0x0),
323   fh1ZPhySel(0x0), 
324   fh1Z(0x0), 
325   fh1ZSelect(0x0),
326   fh2NRecJetsPt(0x0),
327   fh2NRecTracksPt(0x0),
328   fh2NConstPt(0x0),
329   fh2NConstLeadingPt(0x0),
330   fh2JetPhiEta(0x0),
331   fh2LeadingJetPhiEta(0x0),
332   fh2JetEtaPt(0x0),
333   fh2LeadingJetEtaPt(0x0),
334   fh2TrackEtaPt(0x0),
335   fh2LeadingTrackEtaPt(0x0),
336   fh2JetsLeadingPhiEta(0x0),
337   fh2JetsLeadingPhiPt(0x0),
338   fh2TracksLeadingPhiEta(0x0),
339   fh2TracksLeadingPhiPt(0x0),
340   fh2TracksLeadingJetPhiPt(0x0),
341   fh2JetsLeadingPhiPtW(0x0),
342   fh2TracksLeadingPhiPtW(0x0),
343   fh2TracksLeadingJetPhiPtW(0x0),
344   fh2NRecJetsPtRan(0x0),
345   fh2NConstPtRan(0x0),
346   fh2NConstLeadingPtRan(0x0),
347   fh2PtNch(0x0),
348   fh2PtNchRan(0x0),
349   fh2PtNchN(0x0),
350   fh2PtNchNRan(0x0),
351   fh2TracksLeadingJetPhiPtRan(0x0),
352   fh2TracksLeadingJetPhiPtWRan(0x0),
353   fh2PtGenPtSmeared(0x0),
354   fp1Efficiency(0x0),
355   fp1PtResolution(0x0),
356   fHistList(0x0)
357 {
358   //
359   // named ctor
360   //
361
362   for(int i = 0;i<3;i++){
363     fh1BiARandomCones[i] = 0;
364     fh1BiARandomConesRan[i] = 0;
365   }
366   for(int i = 0;i<kMaxCent;i++){
367     fh2JetsLeadingPhiPtC[i] = 0;     
368     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
369     fh2TracksLeadingJetPhiPtC[i] = 0;
370     fh2TracksLeadingJetPhiPtWC[i] = 0;
371   }
372  DefineOutput(1, TList::Class());  
373 }
374
375
376
377 Bool_t AliAnalysisTaskJetCluster::Notify()
378 {
379   //
380   // Implemented Notify() to read the cross sections
381   // and number of trials from pyxsec.root
382   // 
383   return kTRUE;
384 }
385
386 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
387 {
388
389   //
390   // Create the output container
391   //
392
393   fRandom = new TRandom3(0);
394
395
396   // Connect the AOD
397
398
399   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
400
401   
402
403   if(fNonStdBranch.Length()!=0)
404     {
405       // only create the output branch if we have a name
406       // Create a new branch for jets...
407       //  -> cleared in the UserExec....
408       // here we can also have the case that the brnaches are written to a separate file
409       
410       if(fJetTypes&kJet){
411         fTCAJetsOut = new TClonesArray("AliAODJet", 0);
412         fTCAJetsOut->SetName(fNonStdBranch.Data());
413         AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
414       }
415
416       if(fJetTypes&kJetRan){
417         fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
418         fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
419         if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
420           fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
421         AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
422       }
423
424       if(fUseBackgroundCalc){
425         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
426           fAODJetBackgroundOut = new AliAODJetEventBackground();
427           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
428           if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
429             fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
430
431           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());  
432         }
433       }
434       // create the branch for the random cones with the same R 
435       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
436       if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
437         cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
438
439       if(fNRandomCones>0){
440         if(fJetTypes&kRC){
441           if(!AODEvent()->FindListObject(cName.Data())){
442             fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
443             fTCARandomConesOut->SetName(cName.Data());
444             AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
445           }
446         }
447         // create the branch with the random for the random cones on the random event
448         if(fJetTypes&kRCRan){
449           cName = Form("%sRandomCone_random",fNonStdBranch.Data());
450           if(!AODEvent()->FindListObject(cName.Data())){
451             fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
452             fTCARandomConesOutRan->SetName(cName.Data());
453             AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
454           }
455         }
456       }
457     
458       if(fNonStdFile.Length()!=0){
459         // 
460         // case that we have an AOD extension we need to fetch the jets from the extended output
461         // we identify the extension aod event by looking for the branchname
462         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
463         // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
464         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
465       }
466     }
467
468
469   if(!fHistList)fHistList = new TList();
470   fHistList->SetOwner();
471   PostData(1, fHistList); // post data in any case once
472
473   Bool_t oldStatus = TH1::AddDirectoryStatus();
474   TH1::AddDirectory(kFALSE);
475
476   //
477   //  Histogram
478     
479   const Int_t nBinPt = 100;
480   Double_t binLimitsPt[nBinPt+1];
481   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
482     if(iPt == 0){
483       binLimitsPt[iPt] = 0.0;
484     }
485     else {// 1.0
486       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
487     }
488   }
489   
490   const Int_t nBinPhi = 90;
491   Double_t binLimitsPhi[nBinPhi+1];
492   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
493     if(iPhi==0){
494       binLimitsPhi[iPhi] = -1.*TMath::Pi();
495     }
496     else{
497       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
498     }
499   }
500
501
502
503   const Int_t nBinEta = 40;
504   Double_t binLimitsEta[nBinEta+1];
505   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
506     if(iEta==0){
507       binLimitsEta[iEta] = -2.0;
508     }
509     else{
510       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
511     }
512   }
513
514   const int nChMax = 5000;
515
516   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
517   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
518
519   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
520   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
521
522
523   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
524   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
525
526   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
527   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
528   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
529   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
530
531
532   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
533   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
534   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
535
536   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
537   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
538   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
539   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
540   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
541   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
542   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
543   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
544   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
545   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
546
547   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
548   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
549   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
550
551   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
552   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
553   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
554
555   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
556   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
557   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
558   // 
559
560
561   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
562   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
563   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
564   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
565
566   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
567   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);
568   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);
569   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);
570
571
572
573   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
574                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
575   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
576                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
577
578   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
579                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
580   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
581                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
582
583   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
584                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
585   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
586                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
587
588
589
590   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
591                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
592   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
593                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
594   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
595                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
596   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
597                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
598   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
599                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
600   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
601                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
602
603   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
604                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
605   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
606                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
607
608   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
609                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
610   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
611                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
612
613   //Detector level effects histos
614   fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
615
616   fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
617   fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
618
619   fHistList->Add(fh2PtGenPtSmeared);
620   fHistList->Add(fp1Efficiency);
621   fHistList->Add(fp1PtResolution);
622
623   if(fNRandomCones>0&&fUseBackgroundCalc){
624     for(int i = 0;i<3;i++){
625       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
626       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
627     }
628   }
629
630   for(int i = 0;i < kMaxCent;i++){
631     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
632     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
633     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
634     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
635   }
636
637   const Int_t saveLevel = 3; // large save level more histos
638   if(saveLevel>0){
639     fHistList->Add(fh1Xsec);
640     fHistList->Add(fh1Trials);
641
642     fHistList->Add(fh1NJetsRec);
643     fHistList->Add(fh1NConstRec);
644     fHistList->Add(fh1NConstLeadingRec);
645     fHistList->Add(fh1PtJetsRecIn);
646     fHistList->Add(fh1PtJetsLeadingRecIn);
647     fHistList->Add(fh1PtTracksRecIn);
648     fHistList->Add(fh1PtTracksLeadingRecIn);
649     fHistList->Add(fh1PtJetConstRec);
650     fHistList->Add(fh1PtJetConstLeadingRec);
651     fHistList->Add(fh1NJetsRecRan);
652     fHistList->Add(fh1NConstRecRan);
653     fHistList->Add(fh1PtJetsLeadingRecInRan);
654     fHistList->Add(fh1NConstLeadingRecRan);
655     fHistList->Add(fh1PtJetsRecInRan);
656     fHistList->Add(fh1Nch);
657     fHistList->Add(fh1Centrality);
658     fHistList->Add(fh1CentralitySelect);
659     fHistList->Add(fh1CentralityPhySel);
660     fHistList->Add(fh1Z);
661     fHistList->Add(fh1ZSelect);
662     fHistList->Add(fh1ZPhySel);
663     if(fNRandomCones>0&&fUseBackgroundCalc){
664       for(int i = 0;i<3;i++){
665         fHistList->Add(fh1BiARandomCones[i]);
666         fHistList->Add(fh1BiARandomConesRan[i]);
667       }
668     }
669   for(int i = 0;i < kMaxCent;i++){
670     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
671     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
672     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
673     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
674   }
675
676     fHistList->Add(fh2NRecJetsPt);
677     fHistList->Add(fh2NRecTracksPt);
678     fHistList->Add(fh2NConstPt);
679     fHistList->Add(fh2NConstLeadingPt);
680     fHistList->Add(fh2PtNch);
681     fHistList->Add(fh2PtNchRan);
682     fHistList->Add(fh2PtNchN);
683     fHistList->Add(fh2PtNchNRan);
684     fHistList->Add(fh2JetPhiEta);
685     fHistList->Add(fh2LeadingJetPhiEta);
686     fHistList->Add(fh2JetEtaPt);
687     fHistList->Add(fh2LeadingJetEtaPt);
688     fHistList->Add(fh2TrackEtaPt);
689     fHistList->Add(fh2LeadingTrackEtaPt);
690     fHistList->Add(fh2JetsLeadingPhiEta );
691     fHistList->Add(fh2JetsLeadingPhiPt);
692     fHistList->Add(fh2TracksLeadingPhiEta);
693     fHistList->Add(fh2TracksLeadingPhiPt);
694     fHistList->Add(fh2TracksLeadingJetPhiPt);
695     fHistList->Add(fh2JetsLeadingPhiPtW);
696     fHistList->Add(fh2TracksLeadingPhiPtW);
697     fHistList->Add(fh2TracksLeadingJetPhiPtW);
698     fHistList->Add(fh2NRecJetsPtRan);
699     fHistList->Add(fh2NConstPtRan);
700     fHistList->Add(fh2NConstLeadingPtRan);
701     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
702     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
703   }
704
705   // =========== Switch on Sumw2 for all histos ===========
706   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
707     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
708     if (h1){
709       h1->Sumw2();
710       continue;
711     }
712     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
713     if(hn)hn->Sumw2();
714   }
715   TH1::AddDirectory(oldStatus);
716 }
717
718 void AliAnalysisTaskJetCluster::LocalInit()
719 {
720   //
721   // Initialization
722   //
723
724   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
725
726   if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
727   if(fUseTrEfficiencyFromOADB)   LoadTrEfficiencyRootFileFromOADB();
728
729
730   FitMomentumResolution();
731
732 }
733
734 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
735 {
736
737   // handle and reset the output jet branch 
738
739   if(fTCAJetsOut)fTCAJetsOut->Delete();
740   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
741   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
742   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
743   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
744
745   AliAODJetEventBackground* externalBackground = 0;
746   if(!externalBackground&&fBackgroundBranch.Length()){
747     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
748     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
749     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
750   }
751   //
752   // Execute analysis for current event
753   //
754   AliESDEvent *fESD = 0;
755   if(fUseAODTrackInput){    
756     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
757     if(!fAOD){
758       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
759       return;
760     }
761     // fetch the header
762   }
763   else{
764     //  assume that the AOD is in the general output...
765     fAOD  = AODEvent();
766     if(!fAOD){
767       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
768       return;
769     }
770     if(fDebug>0){
771       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
772     }
773   }
774
775   //Check if information is provided detector level effects
776   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
777   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
778   
779   Bool_t selectEvent =  false;
780   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
781
782   Float_t cent = 0;
783   Float_t zVtx  = 0;
784   Int_t cenClass = -1;
785   if(fAOD){
786     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
787     TString vtxTitle(vtxAOD->GetTitle());
788     zVtx = vtxAOD->GetZ();
789
790     cent = fAOD->GetHeader()->GetCentrality();
791     if(cent<10)cenClass = 0;
792     else if(cent<30)cenClass = 1;
793     else if(cent<50)cenClass = 2;
794     else if(cent<80)cenClass = 3;
795     if(physicsSelection){
796       fh1CentralityPhySel->Fill(cent);
797       fh1ZPhySel->Fill(zVtx);
798     }
799
800     if(fEventSelection){
801       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
802         Float_t yvtx = vtxAOD->GetY();
803         Float_t xvtx = vtxAOD->GetX();
804         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
805         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
806           if(physicsSelection){
807             selectEvent = true;
808           }
809         }
810       }
811       if(fCentCutUp>0){
812         if(cent<fCentCutLo||cent>fCentCutUp){
813           selectEvent = false;
814         }
815       }
816     }else{
817       selectEvent = true;
818     }
819   }
820   
821
822   if(!selectEvent){
823     PostData(1, fHistList);
824     return;
825   }
826   fh1Centrality->Fill(cent);  
827   fh1Z->Fill(zVtx);
828   fh1Trials->Fill("#sum{ntrials}",1);
829   
830
831   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
832
833   // ==== General variables needed
834
835
836
837   // we simply fetch the tracks/mc particles as a list of AliVParticles
838
839   TList recParticles;
840   TList genParticles;
841
842   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
843   Float_t nCh = recParticles.GetEntries(); 
844   fh1Nch->Fill(nCh);
845   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
846   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
847   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
848
849   // find the jets....
850
851   vector<fastjet::PseudoJet> inputParticlesRec;
852   vector<fastjet::PseudoJet> inputParticlesRecRan;
853   
854   // Generate the random cones
855   
856   AliAODJet vTmpRan(1,0,0,1);
857   for(int i = 0; i < recParticles.GetEntries(); i++){
858     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
859
860     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
861     // we take total momentum here
862
863     if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
864       //Add particles to fastjet in case we are not running toy model
865       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
866       jInp.set_user_index(i);
867       inputParticlesRec.push_back(jInp);
868     }
869     else if(fUseDiceEfficiency) {
870
871       // Dice to decide if particle is kept or not - toy  model for efficiency
872       //
873       Double_t rnd = fRandom->Uniform(1.);
874       Double_t pT = vp->Pt();
875       Double_t eff[3] = {0.};
876       Double_t pTtmp = pT;
877       if(pT>10.) pTtmp = 10.;
878       Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
879       Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
880       Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
881       Int_t cat[3] = {0};
882       //Sort efficiencies from large to small
883       if(eff1>eff2 && eff1>eff3) { 
884         eff[0] = eff1; 
885         cat[0] = 1;
886         if(eff2>eff3) {
887           eff[1] = eff2;
888           eff[2] = eff3; 
889           cat[1] = 2;
890           cat[2] = 3;
891         } else {
892           eff[1] = eff3;
893           eff[2] = eff2; 
894           cat[1] = 3;
895           cat[2] = 2;
896         }
897       }
898       else if(eff2>eff1 && eff2>eff3) {
899         eff[0] = eff2;
900         cat[0] = 2;
901         if(eff1>eff3) {
902           eff[1] = eff1;
903           eff[2] = eff3; 
904           cat[1] = 1;
905           cat[2] = 3;
906         } else {
907           eff[1] = eff3;
908           eff[2] = eff1; 
909           cat[1] = 3;
910           cat[2] = 1;
911         }
912       }
913       else if(eff3>eff1 && eff3>eff2) {
914         eff[0] = eff3;
915         cat[0] = 3;
916         if(eff1>eff2) {
917           eff[1] = eff1;
918           eff[2] = eff2; 
919           cat[1] = 1;
920           cat[2] = 2;
921         } else {
922           eff[1] = eff2;
923           eff[2] = eff1; 
924           cat[1] = 2;
925           cat[2] = 1;
926         }
927       }
928       
929       Double_t sumEff = eff[0]+eff[1]+eff[2];
930       fp1Efficiency->Fill(vp->Pt(),sumEff);
931       if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
932
933       if(fUseTrPtResolutionSmearing) {
934         //Smear momentum of generated particle
935         Double_t smear = 1.;
936         //Select hybrid track category
937         if(rnd<=eff[2]) 
938           smear = GetMomentumSmearing(cat[2],pT);
939         else if(rnd<=(eff[2]+eff[1])) 
940           smear = GetMomentumSmearing(cat[1],pT);
941         else 
942           smear = GetMomentumSmearing(cat[0],pT);
943
944         fp1PtResolution->Fill(vp->Pt(),smear);
945
946         Double_t sigma = vp->Pt()*smear;
947         Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
948         
949         Double_t phi   = vp->Phi();
950         Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
951         Double_t pX    = pTrec * TMath::Cos(phi);
952         Double_t pY    = pTrec * TMath::Sin(phi);
953         Double_t pZ    = pTrec/TMath::Tan(theta);
954         Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
955
956         fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
957
958         fastjet::PseudoJet jInp(pX,pY,pZ,p);
959         jInp.set_user_index(i);
960         inputParticlesRec.push_back(jInp);
961
962       }
963       else {
964         fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
965         jInp.set_user_index(i);
966         inputParticlesRec.push_back(jInp);
967
968       }
969
970     }
971
972     // the randomized input changes eta and phi, but keeps the p_T
973     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
974       Double_t pT = vp->Pt();
975       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
976       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
977       
978       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
979       Double_t pZ = pT/TMath::Tan(theta);
980
981       Double_t pX = pT * TMath::Cos(phi);
982       Double_t pY = pT * TMath::Sin(phi);
983       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
984       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
985
986       jInpRan.set_user_index(i);
987       inputParticlesRecRan.push_back(jInpRan);
988       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
989     }
990
991     // fill the tref array, only needed when we write out jets
992     if(fTCAJetsOut){
993       if(i == 0){
994         fRef->Delete(); // make sure to delete before placement new...
995         if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
996           new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
997         } 
998       }
999       if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp);  //TRefArray does not work with toy model ...
1000     }
1001   }// recparticles
1002
1003   if(inputParticlesRec.size()==0){
1004     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1005     PostData(1, fHistList);
1006     return;
1007   }
1008   
1009   // run fast jet
1010   // employ setters for these...
1011
1012  
1013   // now create the object that holds info about ghosts                        
1014   /*
1015   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1016     // reduce CPU time...
1017     fGhostArea = 0.5; 
1018     fActiveAreaRepeats = 0; 
1019   }
1020   */
1021
1022   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1023   fastjet::AreaType areaType =   fastjet::active_area;
1024   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1025   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1026   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1027   
1028   //range where to compute background
1029   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1030   phiMin = 0;
1031   phiMax = 2*TMath::Pi();
1032   rapMax = fGhostEtamax - fRparam;
1033   rapMin = - fGhostEtamax + fRparam;
1034   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1035  
1036
1037   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1038   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1039
1040  
1041   fh1NJetsRec->Fill(sortedJets.size());
1042
1043  // loop over all jets an fill information, first on is the leading jet
1044
1045   Int_t nRecOver = inclusiveJets.size();
1046   Int_t nRec     = inclusiveJets.size();
1047   if(inclusiveJets.size()>0){
1048     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1049     Double_t area = clustSeq.area(sortedJets[0]);
1050     leadingJet.SetEffArea(area,0);
1051     Float_t pt = leadingJet.Pt();
1052     Int_t nAodOutJets = 0;
1053     Int_t nAodOutTracks = 0;
1054     AliAODJet *aodOutJet = 0;
1055
1056     Int_t iCount = 0;
1057     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1058       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1059       while(pt<ptCut&&iCount<nRec){
1060         nRecOver--;
1061         iCount++;
1062         if(iCount<nRec){
1063           pt = sortedJets[iCount].perp();
1064         }
1065       }
1066       if(nRecOver<=0)break;
1067       fh2NRecJetsPt->Fill(ptCut,nRecOver);
1068     }
1069     Float_t phi = leadingJet.Phi();
1070     if(phi<0)phi+=TMath::Pi()*2.;    
1071     Float_t eta = leadingJet.Eta();
1072     Float_t pTback = 0;
1073     if(externalBackground){
1074       // carefull has to be filled in a task before
1075       // todo, ReArrange to the botom 
1076      pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1077     }
1078     pt = leadingJet.Pt() - pTback;
1079     // correlation of leading jet with tracks
1080     TIterator *recIter = recParticles.MakeIterator();
1081     recIter->Reset();
1082     AliVParticle *tmpRecTrack = 0;
1083     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1084       Float_t tmpPt = tmpRecTrack->Pt();
1085       // correlation
1086       Float_t tmpPhi =  tmpRecTrack->Phi();     
1087       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1088       Float_t dPhi = phi - tmpPhi;
1089       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1090       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1091       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1092       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1093       if(cenClass>=0){
1094         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1095         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1096       }
1097
1098     }  
1099     
1100    
1101     TLorentzVector vecareab;
1102     for(int j = 0; j < nRec;j++){
1103       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1104       aodOutJet = 0;
1105       nAodOutTracks = 0;
1106       Float_t tmpPt = tmpRec.Pt();  
1107       
1108       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1109         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1110         aodOutJet->GetRefTracks()->Clear();
1111         Double_t area1 = clustSeq.area(sortedJets[j]);
1112         aodOutJet->SetEffArea(area1,0);
1113         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);  
1114         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());     
1115         aodOutJet->SetVectorAreaCharged(&vecareab);
1116       }
1117
1118
1119       Float_t tmpPtBack = 0;
1120       if(externalBackground){
1121         // carefull has to be filled in a task before
1122        // todo, ReArrange to the botom
1123         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1124       }
1125       tmpPt = tmpPt - tmpPtBack;
1126       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1127       
1128       fh1PtJetsRecIn->Fill(tmpPt);
1129       // Fill Spectra with constituentsemacs
1130       const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1131
1132       fh1NConstRec->Fill(constituents.size());
1133       fh2PtNch->Fill(nCh,tmpPt);
1134       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1135       fh2NConstPt->Fill(tmpPt,constituents.size());
1136       // loop over constiutents and fill spectrum
1137    
1138       AliVParticle *partLead = 0;
1139       Float_t ptLead = -1;
1140
1141       for(unsigned int ic = 0; ic < constituents.size();ic++){
1142         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1143         if(!part) continue;
1144         fh1PtJetConstRec->Fill(part->Pt());
1145         if(aodOutJet){
1146           if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1147           if(part->Pt()>fMaxTrackPtInJet){
1148             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1149           }
1150         }
1151         if(part->Pt()>ptLead){
1152           ptLead = part->Pt();
1153           partLead = part;
1154         }
1155         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1156       }
1157
1158       AliAODTrack *aodT = 0;
1159       if(partLead&&aodOutJet){
1160         //set pT of leading constituent of jet
1161         aodOutJet->SetPtLeading(partLead->Pt());
1162         aodT = dynamic_cast<AliAODTrack*>(partLead);
1163         if(aodT){
1164           if(aodT->TestFilterBit(fFilterMaskBestPt)){
1165             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1166           }
1167         }
1168       }
1169       
1170      // correlation
1171      Float_t tmpPhi =  tmpRec.Phi();
1172      Float_t tmpEta =  tmpRec.Eta();
1173      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;        
1174      if(j==0){
1175        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1176        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1177        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1178        fh1NConstLeadingRec->Fill(constituents.size());
1179        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1180        continue;
1181      }
1182      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1183      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1184      Float_t dPhi = phi - tmpPhi;
1185      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1186      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1187      Float_t dEta = eta - tmpRec.Eta();
1188      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1189      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1190      if(cenClass>=0){
1191        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1192        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1193      }
1194      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1195     }// loop over reconstructed jets
1196    delete recIter;
1197
1198
1199
1200    // Add the random cones...
1201    if(fNRandomCones>0&&fTCARandomConesOut){       
1202      // create a random jet within the acceptance
1203      Double_t etaMax = fTrackEtaWindow - fRparam;
1204      Int_t nCone = 0;
1205      Int_t nConeRan = 0;
1206      Double_t pTC = 1; // small number
1207      for(int ir = 0;ir < fNRandomCones;ir++){
1208        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1209        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1210        // massless jet
1211        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1212        Double_t pZC = pTC/TMath::Tan(thetaC);
1213        Double_t pXC = pTC * TMath::Cos(phiC);
1214        Double_t pYC = pTC * TMath::Sin(phiC);
1215        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1216        AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
1217        bool skip = false;
1218        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1219          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1220          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1221            skip = true;
1222            break;
1223          }
1224        }
1225        // test for overlap with previous cones to avoid double counting
1226        for(int iic = 0;iic<ir;iic++){
1227          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1228          if(iicone){
1229            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1230              skip = true;
1231              break;
1232            }
1233          }
1234        }
1235        if(skip)continue;
1236        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1237        tmpRecC.SetPtLeading(-1.);
1238        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1239        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1240      }// loop over random cones creation
1241
1242   
1243      // loop over the reconstructed particles and add up the pT in the random cones
1244      // maybe better to loop over randomized particles not in the real jets...
1245      // but this by definition brings dow average energy in the whole  event
1246      AliAODJet vTmpRanR(1,0,0,1);
1247      for(int i = 0; i < recParticles.GetEntries(); i++){
1248        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1249        if(fTCARandomConesOut){
1250          for(int ir = 0;ir < fNRandomCones;ir++){
1251            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);  
1252            if(jC&&jC->DeltaR(vp)<fRparam){
1253              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1254              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1255              if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1256            }
1257          }  
1258        }// add up energy in cone
1259
1260        // the randomized input changes eta and phi, but keeps the p_T
1261        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1262          Double_t pTR = vp->Pt();
1263          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1264          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1265          
1266          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
1267          Double_t pZR = pTR/TMath::Tan(thetaR);
1268          
1269          Double_t pXR = pTR * TMath::Cos(phiR);
1270          Double_t pYR = pTR * TMath::Sin(phiR);
1271          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
1272          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1273          if(fTCARandomConesOutRan){
1274            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1275              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);  
1276              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1277                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1278                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1279                if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1280              }
1281            }  
1282          }
1283        }
1284      }// loop over recparticles
1285     
1286      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1287      if(fTCARandomConesOut){
1288        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1289          // rescale the momentum vectors for the random cones
1290          
1291          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1292          if(rC){
1293            Double_t etaC = rC->Eta();
1294            Double_t phiC = rC->Phi();
1295            // massless jet, unit vector
1296            pTC = rC->ChargedBgEnergy();
1297            if(pTC<=0)pTC = 0.001; // for almost empty events
1298            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1299            Double_t pZC = pTC/TMath::Tan(thetaC);
1300            Double_t pXC = pTC * TMath::Cos(phiC);
1301            Double_t pYC = pTC * TMath::Sin(phiC);
1302            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1303            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1304            rC->SetBgEnergy(0,0);
1305            rC->SetEffArea(jetArea,0);
1306          }
1307        }
1308      }
1309      if(fTCARandomConesOutRan){
1310        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1311          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1312          // same wit random
1313          if(rC){
1314            Double_t etaC = rC->Eta();
1315            Double_t phiC = rC->Phi();
1316            // massless jet, unit vector
1317            pTC = rC->ChargedBgEnergy();
1318            if(pTC<=0)pTC = 0.001;// for almost empty events
1319            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1320            Double_t pZC = pTC/TMath::Tan(thetaC);
1321            Double_t pXC = pTC * TMath::Cos(phiC);
1322            Double_t pYC = pTC * TMath::Sin(phiC);
1323            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1324            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1325            rC->SetBgEnergy(0,0);
1326            rC->SetEffArea(jetArea,0);
1327          }
1328        }
1329      }
1330    }// if(fNRandomCones
1331   
1332    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1333  
1334
1335
1336
1337
1338    if(fAODJetBackgroundOut){
1339      vector<fastjet::PseudoJet> jets2=sortedJets;
1340      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
1341      Double_t bkg1=0;
1342      Double_t sigma1=0.;
1343      Double_t meanarea1=0.;
1344      Double_t bkg2=0;
1345      Double_t sigma2=0.;
1346      Double_t meanarea2=0.;
1347
1348      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1349      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1350
1351      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1352      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
1353      
1354      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1355      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1356      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1357      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1358
1359    }
1360   }
1361    
1362
1363   
1364   
1365  
1366   // fill track information
1367   Int_t nTrackOver = recParticles.GetSize();
1368   // do the same for tracks and jets
1369
1370   if(nTrackOver>0){
1371    TIterator *recIter = recParticles.MakeIterator();
1372    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1373    Float_t pt = tmpRec->Pt();
1374
1375    //    Printf("Leading track p_t %3.3E",pt);
1376    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1377      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1378      while(pt<ptCut&&tmpRec){
1379        nTrackOver--;
1380        tmpRec = (AliVParticle*)(recIter->Next()); 
1381        if(tmpRec){
1382          pt = tmpRec->Pt();
1383        }
1384      }
1385      if(nTrackOver<=0)break;
1386      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1387    }
1388    
1389    recIter->Reset();
1390    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1391    Float_t phi = leading->Phi();
1392    if(phi<0)phi+=TMath::Pi()*2.;    
1393    Float_t eta = leading->Eta();
1394    pt  = leading->Pt();
1395    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1396      Float_t tmpPt = tmpRec->Pt();
1397      Float_t tmpEta = tmpRec->Eta();
1398      fh1PtTracksRecIn->Fill(tmpPt);
1399      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1400      if(tmpRec==leading){
1401        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1402        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1403        continue;
1404      }
1405       // correlation
1406      Float_t tmpPhi =  tmpRec->Phi();
1407      
1408      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1409      Float_t dPhi = phi - tmpPhi;
1410      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1411      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1412      Float_t dEta = eta - tmpRec->Eta();
1413      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1414      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1415      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1416    }  
1417    delete recIter;
1418  }
1419
1420  // find the random jets
1421
1422  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1423   
1424  // fill the jet information from random track
1425  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1426  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1427
1428   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1429
1430  // loop over all jets an fill information, first on is the leading jet
1431
1432  Int_t nRecOverRan = inclusiveJetsRan.size();
1433  Int_t nRecRan     = inclusiveJetsRan.size();
1434
1435  if(inclusiveJetsRan.size()>0){
1436    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1437    Float_t pt = leadingJet.Pt();
1438    
1439    Int_t iCount = 0;
1440    TLorentzVector vecarearanb;
1441
1442    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1443      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1444       while(pt<ptCut&&iCount<nRecRan){
1445         nRecOverRan--;
1446         iCount++;
1447         if(iCount<nRecRan){
1448           pt = sortedJetsRan[iCount].perp();
1449         }
1450       }
1451       if(nRecOverRan<=0)break;
1452       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1453     }
1454     Float_t phi = leadingJet.Phi();
1455     if(phi<0)phi+=TMath::Pi()*2.;    
1456     pt = leadingJet.Pt();
1457
1458     // correlation of leading jet with random tracks
1459
1460     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1461       { 
1462         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1463         // correlation
1464         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1465         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1466         Float_t dPhi = phi - tmpPhi;
1467         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1468         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1469         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1470         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1471       }  
1472
1473     Int_t nAodOutJetsRan = 0;
1474      AliAODJet *aodOutJetRan = 0;
1475     for(int j = 0; j < nRecRan;j++){
1476       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1477       Float_t tmpPt = tmpRec.Pt();
1478       fh1PtJetsRecInRan->Fill(tmpPt);
1479       // Fill Spectra with constituents
1480       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1481       fh1NConstRecRan->Fill(constituents.size());
1482       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1483       fh2PtNchRan->Fill(nCh,tmpPt);
1484       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1485
1486
1487      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1488        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1489        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1490        aodOutJetRan->GetRefTracks()->Clear();
1491        aodOutJetRan->SetEffArea(arearan,0);
1492        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1493        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1494        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1495
1496      }
1497
1498       // correlation
1499       Float_t tmpPhi =  tmpRec.Phi();
1500       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1501
1502       if(j==0){
1503         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1504         fh1NConstLeadingRecRan->Fill(constituents.size());
1505         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1506         continue;
1507       }
1508     }  
1509
1510      
1511     if(fAODJetBackgroundOut){
1512      Double_t bkg3=0.;
1513      Double_t sigma3=0.;
1514      Double_t meanarea3=0.;
1515      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1516      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1517      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1518      /*
1519      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1520      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1521      */
1522     }
1523
1524
1525
1526  }
1527
1528
1529  // do the event selection if activated
1530  if(fJetTriggerPtCut>0){
1531    bool select = false;
1532    Float_t minPt = fJetTriggerPtCut;
1533    /*
1534    // hard coded for now ...
1535    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1536    if(cent<10)minPt = 50;
1537    else if(cent<30)minPt = 42;
1538    else if(cent<50)minPt = 28;
1539    else if(cent<80)minPt = 18;
1540    */
1541    float rho = 0;
1542    if(externalBackground)rho = externalBackground->GetBackground(2);
1543    if(fTCAJetsOut){
1544      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1545        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1546        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1547        if(ptSub>=minPt){
1548          select = true;
1549          break;
1550        }
1551      }
1552    }   
1553  
1554    if(select){
1555      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1556      fh1CentralitySelect->Fill(cent);
1557      fh1ZSelect->Fill(zVtx);
1558      aodH->SetFillAOD(kTRUE);
1559    }
1560  }
1561  if (fDebug > 2){
1562    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1563    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1564    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1565    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1566  }
1567  PostData(1, fHistList);
1568 }
1569
1570 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1571 {
1572   //
1573   // Terminate analysis
1574   //
1575     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1576
1577     if(fMomResH1Fit) delete fMomResH1Fit;
1578     if(fMomResH2Fit) delete fMomResH2Fit;
1579     if(fMomResH3Fit) delete fMomResH3Fit;
1580     
1581 }
1582
1583
1584 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1585
1586   //
1587   // get list of tracks/particles for different types
1588   //
1589
1590   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1591
1592   Int_t iCount = 0;
1593   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1594     if(type!=kTrackAODextraonly) {
1595       AliAODEvent *aod = 0;
1596       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1597       else aod = AODEvent();
1598       if(!aod){
1599         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1600         return iCount;
1601       }
1602
1603       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1604         AliAODTrack *tr = aod->GetTrack(it);
1605         Bool_t bGood = false;
1606         if(fFilterType == 0)bGood = true;
1607         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1608         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1609         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1610           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
1611           continue;
1612         }
1613         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1614           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
1615           continue;
1616         }
1617         if(tr->Pt()<fTrackPtCut){
1618           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
1619           continue;
1620         }
1621         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
1622         list->Add(tr);
1623         iCount++;
1624       }
1625     }
1626     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1627       AliAODEvent *aod = 0;
1628       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1629       else aod = AODEvent();
1630       
1631       if(!aod){
1632         return iCount;
1633       }
1634       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1635       if(!aodExtraTracks)return iCount;
1636       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1637         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1638         if (!track) continue;
1639
1640         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1641         if(!trackAOD)continue;
1642         Bool_t bGood = false;
1643         if(fFilterType == 0)bGood = true;
1644         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1645         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1646         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1647         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1648         if(trackAOD->Pt()<fTrackPtCut) continue;
1649         list->Add(trackAOD);
1650         iCount++;
1651       }
1652     }
1653   }
1654   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1655     AliMCEvent* mcEvent = MCEvent();
1656     if(!mcEvent)return iCount;
1657     // we want to have alivpartilces so use get track
1658     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1659       if(!mcEvent->IsPhysicalPrimary(it))continue;
1660       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1661       if(type == kTrackKineAll){
1662         if(part->Pt()<fTrackPtCut)continue;
1663         list->Add(part);
1664         iCount++;
1665       }
1666       else if(type == kTrackKineCharged){
1667         if(part->Particle()->GetPDG()->Charge()==0)continue;
1668         if(part->Pt()<fTrackPtCut)continue;
1669         list->Add(part);
1670         iCount++;
1671       }
1672     }
1673   }
1674   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1675     AliAODEvent *aod = 0;
1676     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1677     else aod = AODEvent();
1678     if(!aod)return iCount;
1679     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1680     if(!tca)return iCount;
1681     for(int it = 0;it < tca->GetEntriesFast();++it){
1682       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1683       if(!part->IsPhysicalPrimary())continue;
1684       if(type == kTrackAODMCAll){
1685         if(part->Pt()<fTrackPtCut)continue;
1686         list->Add(part);
1687         iCount++;
1688       }
1689       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1690         if(part->Charge()==0)continue;
1691         if(part->Pt()<fTrackPtCut)continue;
1692         if(kTrackAODMCCharged){
1693           list->Add(part);
1694         }
1695         else {
1696           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1697           list->Add(part);
1698         }
1699         iCount++;
1700       }
1701     }
1702   }// AODMCparticle
1703   list->Sort();
1704   return iCount;
1705 }
1706
1707 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1708
1709   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1710
1711   if(!f)return;
1712
1713   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1714   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1715   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1716
1717   SetSmearResolution(kTRUE);
1718   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1719   
1720
1721 }
1722
1723 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1724
1725   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1726   if(!f)return;
1727
1728   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1729   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1730   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1731
1732   SetDiceEfficiency(kTRUE);
1733
1734   if(fChangeEfficiencyFraction>0.) {
1735
1736     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1737     
1738     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1739       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1740       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1741     }
1742
1743     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1744
1745   }
1746   else
1747     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1748
1749 }
1750
1751 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1752
1753   //
1754   // set mom res profiles
1755   //
1756
1757   if(fMomResH1) delete fMomResH1;
1758   if(fMomResH2) delete fMomResH2;
1759   if(fMomResH3) delete fMomResH3;
1760
1761   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1762   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1763   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1764
1765 }
1766
1767 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1768   //
1769   // set tracking efficiency histos
1770   //
1771
1772   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1773   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1774   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1775 }
1776
1777 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1778
1779   //
1780   // Get smearing on generated momentum
1781   //
1782
1783   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1784
1785   TProfile *fMomRes = 0x0;
1786   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1787   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1788   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1789
1790   if(!fMomRes) {
1791     return 0.;
1792   }
1793
1794
1795   Double_t smear = 0.;
1796
1797   if(pt>20.) {
1798     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1799     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1800     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1801   }
1802   else {
1803
1804     Int_t bin = fMomRes->FindBin(pt);
1805
1806     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1807
1808   }
1809
1810   if(fMomRes) delete fMomRes;
1811   
1812   return smear;
1813 }
1814
1815 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1816   //
1817   // Fit linear function on momentum resolution at high pT
1818   //
1819
1820   if(!fMomResH1Fit && fMomResH1) {
1821     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1822     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1823     fMomResH1Fit ->SetRange(5.,100.);
1824   }
1825
1826   if(!fMomResH2Fit && fMomResH2) {
1827     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1828     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1829     fMomResH2Fit ->SetRange(5.,100.);
1830   }
1831
1832   if(!fMomResH3Fit && fMomResH3) {
1833     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1834     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1835     fMomResH3Fit ->SetRange(5.,100.);
1836   }
1837
1838 }
1839
1840 /*
1841 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1842   for(int i = 0; i < particles.GetEntries(); i++){
1843     AliVParticle *vp = (AliVParticle*)particles.At(i);
1844     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1845     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1846     jInp.set_user_index(i);
1847     inputParticles.push_back(jInp);
1848   }
1849
1850   return 0;
1851
1852 }
1853 */