1326951d111c47304a3d1ee47a8558611c6b388c
[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){
1160         if(aodOutJet){
1161           //set pT of leading constituent of jet
1162           aodOutJet->SetPtLeading(partLead->Pt());
1163           aodT = dynamic_cast<AliAODTrack*>(partLead);
1164           if(aodT){
1165             if(aodT->TestFilterBit(fFilterMaskBestPt)){
1166               aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1167             }
1168           }
1169         }
1170       }
1171     
1172      // correlation
1173      Float_t tmpPhi =  tmpRec.Phi();
1174      Float_t tmpEta =  tmpRec.Eta();
1175      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;        
1176      if(j==0){
1177        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1178        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1179        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1180        fh1NConstLeadingRec->Fill(constituents.size());
1181        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1182        continue;
1183      }
1184      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1185      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1186      Float_t dPhi = phi - tmpPhi;
1187      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1188      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1189      Float_t dEta = eta - tmpRec.Eta();
1190      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1191      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1192      if(cenClass>=0){
1193        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1194        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1195      }
1196      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1197     }// loop over reconstructed jets
1198    delete recIter;
1199
1200
1201
1202    // Add the random cones...
1203    if(fNRandomCones>0&&fTCARandomConesOut){       
1204      // create a random jet within the acceptance
1205      Double_t etaMax = fTrackEtaWindow - fRparam;
1206      Int_t nCone = 0;
1207      Int_t nConeRan = 0;
1208      Double_t pTC = 1; // small number
1209      for(int ir = 0;ir < fNRandomCones;ir++){
1210        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1211        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1212        // massless jet
1213        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1214        Double_t pZC = pTC/TMath::Tan(thetaC);
1215        Double_t pXC = pTC * TMath::Cos(phiC);
1216        Double_t pYC = pTC * TMath::Sin(phiC);
1217        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1218        AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
1219        bool skip = false;
1220        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1221          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1222          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1223            skip = true;
1224            break;
1225          }
1226        }
1227        // test for overlap with previous cones to avoid double counting
1228        for(int iic = 0;iic<ir;iic++){
1229          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1230          if(iicone){
1231            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1232              skip = true;
1233              break;
1234            }
1235          }
1236        }
1237        if(skip)continue;
1238        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1239        tmpRecC.SetPtLeading(-1.);
1240        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1241        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1242      }// loop over random cones creation
1243
1244   
1245      // loop over the reconstructed particles and add up the pT in the random cones
1246      // maybe better to loop over randomized particles not in the real jets...
1247      // but this by definition brings dow average energy in the whole  event
1248      AliAODJet vTmpRanR(1,0,0,1);
1249      for(int i = 0; i < recParticles.GetEntries(); i++){
1250        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1251        if(fTCARandomConesOut){
1252          for(int ir = 0;ir < fNRandomCones;ir++){
1253            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);  
1254            if(jC&&jC->DeltaR(vp)<fRparam){
1255              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1256              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1257              if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1258            }
1259          }  
1260        }// add up energy in cone
1261
1262        // the randomized input changes eta and phi, but keeps the p_T
1263        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1264          Double_t pTR = vp->Pt();
1265          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1266          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1267          
1268          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
1269          Double_t pZR = pTR/TMath::Tan(thetaR);
1270          
1271          Double_t pXR = pTR * TMath::Cos(phiR);
1272          Double_t pYR = pTR * TMath::Sin(phiR);
1273          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
1274          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1275          if(fTCARandomConesOutRan){
1276            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1277              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);  
1278              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1279                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1280                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1281                if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1282              }
1283            }  
1284          }
1285        }
1286      }// loop over recparticles
1287     
1288      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1289      if(fTCARandomConesOut){
1290        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1291          // rescale the momentum vectors for the random cones
1292          
1293          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1294          if(rC){
1295            Double_t etaC = rC->Eta();
1296            Double_t phiC = rC->Phi();
1297            // massless jet, unit vector
1298            pTC = rC->ChargedBgEnergy();
1299            if(pTC<=0)pTC = 0.001; // for almost empty events
1300            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1301            Double_t pZC = pTC/TMath::Tan(thetaC);
1302            Double_t pXC = pTC * TMath::Cos(phiC);
1303            Double_t pYC = pTC * TMath::Sin(phiC);
1304            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1305            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1306            rC->SetBgEnergy(0,0);
1307            rC->SetEffArea(jetArea,0);
1308          }
1309        }
1310      }
1311      if(fTCARandomConesOutRan){
1312        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1313          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1314          // same wit random
1315          if(rC){
1316            Double_t etaC = rC->Eta();
1317            Double_t phiC = rC->Phi();
1318            // massless jet, unit vector
1319            pTC = rC->ChargedBgEnergy();
1320            if(pTC<=0)pTC = 0.001;// for almost empty events
1321            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1322            Double_t pZC = pTC/TMath::Tan(thetaC);
1323            Double_t pXC = pTC * TMath::Cos(phiC);
1324            Double_t pYC = pTC * TMath::Sin(phiC);
1325            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1326            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1327            rC->SetBgEnergy(0,0);
1328            rC->SetEffArea(jetArea,0);
1329          }
1330        }
1331      }
1332    }// if(fNRandomCones
1333   
1334    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1335  
1336
1337
1338
1339
1340    if(fAODJetBackgroundOut){
1341      vector<fastjet::PseudoJet> jets2=sortedJets;
1342      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
1343      Double_t bkg1=0;
1344      Double_t sigma1=0.;
1345      Double_t meanarea1=0.;
1346      Double_t bkg2=0;
1347      Double_t sigma2=0.;
1348      Double_t meanarea2=0.;
1349
1350      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1351      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1352
1353      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1354      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
1355      
1356      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1357      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1358      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1359      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1360
1361    }
1362   }
1363    
1364
1365   
1366   
1367  
1368   // fill track information
1369   Int_t nTrackOver = recParticles.GetSize();
1370   // do the same for tracks and jets
1371
1372   if(nTrackOver>0){
1373    TIterator *recIter = recParticles.MakeIterator();
1374    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1375    Float_t pt = tmpRec->Pt();
1376
1377    //    Printf("Leading track p_t %3.3E",pt);
1378    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1379      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1380      while(pt<ptCut&&tmpRec){
1381        nTrackOver--;
1382        tmpRec = (AliVParticle*)(recIter->Next()); 
1383        if(tmpRec){
1384          pt = tmpRec->Pt();
1385        }
1386      }
1387      if(nTrackOver<=0)break;
1388      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1389    }
1390    
1391    recIter->Reset();
1392    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1393    Float_t phi = leading->Phi();
1394    if(phi<0)phi+=TMath::Pi()*2.;    
1395    Float_t eta = leading->Eta();
1396    pt  = leading->Pt();
1397    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1398      Float_t tmpPt = tmpRec->Pt();
1399      Float_t tmpEta = tmpRec->Eta();
1400      fh1PtTracksRecIn->Fill(tmpPt);
1401      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1402      if(tmpRec==leading){
1403        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1404        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1405        continue;
1406      }
1407       // correlation
1408      Float_t tmpPhi =  tmpRec->Phi();
1409      
1410      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1411      Float_t dPhi = phi - tmpPhi;
1412      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1413      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1414      Float_t dEta = eta - tmpRec->Eta();
1415      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1416      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1417      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1418    }  
1419    delete recIter;
1420  }
1421
1422  // find the random jets
1423
1424  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1425   
1426  // fill the jet information from random track
1427  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1428  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1429
1430   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1431
1432  // loop over all jets an fill information, first on is the leading jet
1433
1434  Int_t nRecOverRan = inclusiveJetsRan.size();
1435  Int_t nRecRan     = inclusiveJetsRan.size();
1436
1437  if(inclusiveJetsRan.size()>0){
1438    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1439    Float_t pt = leadingJet.Pt();
1440    
1441    Int_t iCount = 0;
1442    TLorentzVector vecarearanb;
1443
1444    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1445      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1446       while(pt<ptCut&&iCount<nRecRan){
1447         nRecOverRan--;
1448         iCount++;
1449         if(iCount<nRecRan){
1450           pt = sortedJetsRan[iCount].perp();
1451         }
1452       }
1453       if(nRecOverRan<=0)break;
1454       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1455     }
1456     Float_t phi = leadingJet.Phi();
1457     if(phi<0)phi+=TMath::Pi()*2.;    
1458     pt = leadingJet.Pt();
1459
1460     // correlation of leading jet with random tracks
1461
1462     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1463       { 
1464         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1465         // correlation
1466         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1467         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1468         Float_t dPhi = phi - tmpPhi;
1469         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1470         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1471         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1472         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1473       }  
1474
1475     Int_t nAodOutJetsRan = 0;
1476      AliAODJet *aodOutJetRan = 0;
1477     for(int j = 0; j < nRecRan;j++){
1478       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1479       Float_t tmpPt = tmpRec.Pt();
1480       fh1PtJetsRecInRan->Fill(tmpPt);
1481       // Fill Spectra with constituents
1482       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1483       fh1NConstRecRan->Fill(constituents.size());
1484       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1485       fh2PtNchRan->Fill(nCh,tmpPt);
1486       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1487
1488
1489      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1490        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1491        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1492        aodOutJetRan->GetRefTracks()->Clear();
1493        aodOutJetRan->SetEffArea(arearan,0);
1494        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1495        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1496        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1497
1498      }
1499
1500       // correlation
1501       Float_t tmpPhi =  tmpRec.Phi();
1502       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1503
1504       if(j==0){
1505         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1506         fh1NConstLeadingRecRan->Fill(constituents.size());
1507         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1508         continue;
1509       }
1510     }  
1511
1512      
1513     if(fAODJetBackgroundOut){
1514      Double_t bkg3=0.;
1515      Double_t sigma3=0.;
1516      Double_t meanarea3=0.;
1517      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1518      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1519      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1520      /*
1521      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1522      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1523      */
1524     }
1525
1526
1527
1528  }
1529
1530
1531  // do the event selection if activated
1532  if(fJetTriggerPtCut>0){
1533    bool select = false;
1534    Float_t minPt = fJetTriggerPtCut;
1535    /*
1536    // hard coded for now ...
1537    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1538    if(cent<10)minPt = 50;
1539    else if(cent<30)minPt = 42;
1540    else if(cent<50)minPt = 28;
1541    else if(cent<80)minPt = 18;
1542    */
1543    float rho = 0;
1544    if(externalBackground)rho = externalBackground->GetBackground(2);
1545    if(fTCAJetsOut){
1546      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1547        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1548        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1549        if(ptSub>=minPt){
1550          select = true;
1551          break;
1552        }
1553      }
1554    }   
1555  
1556    if(select){
1557      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1558      fh1CentralitySelect->Fill(cent);
1559      fh1ZSelect->Fill(zVtx);
1560      aodH->SetFillAOD(kTRUE);
1561    }
1562  }
1563  if (fDebug > 2){
1564    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1565    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1566    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1567    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1568  }
1569  PostData(1, fHistList);
1570 }
1571
1572 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1573 {
1574   //
1575   // Terminate analysis
1576   //
1577     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1578
1579     if(fMomResH1Fit) delete fMomResH1Fit;
1580     if(fMomResH2Fit) delete fMomResH2Fit;
1581     if(fMomResH3Fit) delete fMomResH3Fit;
1582     
1583 }
1584
1585
1586 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1587
1588   //
1589   // get list of tracks/particles for different types
1590   //
1591
1592   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1593
1594   Int_t iCount = 0;
1595   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1596     if(type!=kTrackAODextraonly) {
1597       AliAODEvent *aod = 0;
1598       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1599       else aod = AODEvent();
1600       if(!aod){
1601         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1602         return iCount;
1603       }
1604
1605       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1606         AliAODTrack *tr = aod->GetTrack(it);
1607         Bool_t bGood = false;
1608         if(fFilterType == 0)bGood = true;
1609         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1610         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1611         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1612           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
1613           continue;
1614         }
1615         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1616           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
1617           continue;
1618         }
1619         if(tr->Pt()<fTrackPtCut){
1620           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
1621           continue;
1622         }
1623         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
1624         list->Add(tr);
1625         iCount++;
1626       }
1627     }
1628     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1629       AliAODEvent *aod = 0;
1630       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1631       else aod = AODEvent();
1632       
1633       if(!aod){
1634         return iCount;
1635       }
1636       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1637       if(!aodExtraTracks)return iCount;
1638       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1639         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1640         if (!track) continue;
1641
1642         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1643         if(!trackAOD)continue;
1644         Bool_t bGood = false;
1645         if(fFilterType == 0)bGood = true;
1646         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1647         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1648         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1649         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1650         if(trackAOD->Pt()<fTrackPtCut) continue;
1651         list->Add(trackAOD);
1652         iCount++;
1653       }
1654     }
1655   }
1656   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1657     AliMCEvent* mcEvent = MCEvent();
1658     if(!mcEvent)return iCount;
1659     // we want to have alivpartilces so use get track
1660     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1661       if(!mcEvent->IsPhysicalPrimary(it))continue;
1662       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1663       if(type == kTrackKineAll){
1664         if(part->Pt()<fTrackPtCut)continue;
1665         list->Add(part);
1666         iCount++;
1667       }
1668       else if(type == kTrackKineCharged){
1669         if(part->Particle()->GetPDG()->Charge()==0)continue;
1670         if(part->Pt()<fTrackPtCut)continue;
1671         list->Add(part);
1672         iCount++;
1673       }
1674     }
1675   }
1676   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1677     AliAODEvent *aod = 0;
1678     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1679     else aod = AODEvent();
1680     if(!aod)return iCount;
1681     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1682     if(!tca)return iCount;
1683     for(int it = 0;it < tca->GetEntriesFast();++it){
1684       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1685       if(!part->IsPhysicalPrimary())continue;
1686       if(type == kTrackAODMCAll){
1687         if(part->Pt()<fTrackPtCut)continue;
1688         list->Add(part);
1689         iCount++;
1690       }
1691       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1692         if(part->Charge()==0)continue;
1693         if(part->Pt()<fTrackPtCut)continue;
1694         if(kTrackAODMCCharged){
1695           list->Add(part);
1696         }
1697         else {
1698           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1699           list->Add(part);
1700         }
1701         iCount++;
1702       }
1703     }
1704   }// AODMCparticle
1705   list->Sort();
1706   return iCount;
1707 }
1708
1709 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1710
1711   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1712
1713   if(!f)return;
1714
1715   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1716   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1717   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1718
1719   SetSmearResolution(kTRUE);
1720   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1721   
1722
1723 }
1724
1725 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1726
1727   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1728   if(!f)return;
1729
1730   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1731   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1732   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1733
1734   SetDiceEfficiency(kTRUE);
1735
1736   if(fChangeEfficiencyFraction>0.) {
1737
1738     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1739     
1740     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1741       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1742       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1743     }
1744
1745     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1746
1747   }
1748   else
1749     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1750
1751 }
1752
1753 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1754
1755   //
1756   // set mom res profiles
1757   //
1758
1759   if(fMomResH1) delete fMomResH1;
1760   if(fMomResH2) delete fMomResH2;
1761   if(fMomResH3) delete fMomResH3;
1762
1763   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1764   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1765   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1766
1767 }
1768
1769 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1770   //
1771   // set tracking efficiency histos
1772   //
1773
1774   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1775   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1776   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1777 }
1778
1779 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1780
1781   //
1782   // Get smearing on generated momentum
1783   //
1784
1785   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1786
1787   TProfile *fMomRes = 0x0;
1788   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1789   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1790   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1791
1792   if(!fMomRes) {
1793     return 0.;
1794   }
1795
1796
1797   Double_t smear = 0.;
1798
1799   if(pt>20.) {
1800     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1801     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1802     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1803   }
1804   else {
1805
1806     Int_t bin = fMomRes->FindBin(pt);
1807
1808     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1809
1810   }
1811
1812   if(fMomRes) delete fMomRes;
1813   
1814   return smear;
1815 }
1816
1817 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1818   //
1819   // Fit linear function on momentum resolution at high pT
1820   //
1821
1822   if(!fMomResH1Fit && fMomResH1) {
1823     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1824     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1825     fMomResH1Fit ->SetRange(5.,100.);
1826   }
1827
1828   if(!fMomResH2Fit && fMomResH2) {
1829     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1830     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1831     fMomResH2Fit ->SetRange(5.,100.);
1832   }
1833
1834   if(!fMomResH3Fit && fMomResH3) {
1835     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1836     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1837     fMomResH3Fit ->SetRange(5.,100.);
1838   }
1839
1840 }
1841
1842 /*
1843 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1844   for(int i = 0; i < particles.GetEntries(); i++){
1845     AliVParticle *vp = (AliVParticle*)particles.At(i);
1846     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1847     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1848     jInp.set_user_index(i);
1849     inputParticles.push_back(jInp);
1850   }
1851
1852   return 0;
1853
1854 }
1855 */