pileup rejection for Fragmentation Function and Jet Shape task
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetProperties.cxx
1 // **************************************************************************************
2 // *                                                                                    *
3 // * Task for Jet properties and jet shape analysis in PWG4 Jet Task Force Train for pp *
4 // *                                                                                    *
5 // **************************************************************************************
6
7
8 /**************************************************************************
9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10  *                                                                        *
11  * Author: The ALICE Off-line Project.                                    *
12  * Contributors are mentioned in the code where appropriate.              *
13  *                                                                        *
14  * Permission to use, copy, modify and distribute this software and its   *
15  * documentation strictly for non-commercial purposes is hereby granted   *
16  * without fee, provided that the above copyright notice appears in all   *
17  * copies and that both the copyright notice and this permission notice   *
18  * appear in the supporting documentation. The authors make no claims     *
19  * about the suitability of this software for any purpose. It is          *
20  * provided "as is" without express or implied warranty.                  *
21  **************************************************************************/
22
23 /* $Id: */
24
25 #include "TList.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TH3F.h"
29 #include "TString.h"
30 #include "THnSparse.h"
31 #include "TProfile.h"
32 #include "TFile.h"
33 #include "TKey.h"
34 #include "TRandom3.h"
35
36 #include "AliAODInputHandler.h" 
37 #include "AliAODHandler.h" 
38 #include "AliESDEvent.h"
39 #include "AliAODMCParticle.h"
40 #include "AliAODJet.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenHijingEventHeader.h"
43 #include "AliInputEventHandler.h"
44
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliAnalysisManager.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "AliVParticle.h"
49 #include "AliAODTrack.h"
50 #include "AliVEvent.h"
51
52
53 #include "AliAnalysisTaskJetProperties.h"
54
55 ClassImp(AliAnalysisTaskJetProperties)
56   
57 //_________________________________________________________________________________//
58   
59   AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties()
60     : AliAnalysisTaskSE()
61     ,fESD(0)
62     ,fAOD(0)
63     ,fAODJets(0)  
64     ,fAODExtension(0)
65     ,fNonStdFile("")
66     ,fBranchJets("jets")
67     ,fTrackType(kTrackAOD)
68     ,fJetRejectType(0)
69     ,fRejectPileup(1)
70     ,fUseAODInputJets(kTRUE)
71     ,fFilterMask(0)
72     ,fUsePhysicsSelection(kTRUE)
73     ,fMaxVertexZ(10)
74     ,fNContributors(2)
75     ,fTrackPtCut(0)
76     ,fTrackEtaMin(0)
77     ,fTrackEtaMax(0)
78     ,fJetPtCut(0)
79     ,fJetEtaMin(0)
80     ,fJetEtaMax(0)
81     ,fAvgTrials(0)
82     ,fJetRadius(0.4)
83     ,fJetList(0)
84     ,fTrackList(0)
85     ,fTrackListUE(0)
86     ,fTrackListJet(0)
87     ,fCommonHistList(0)
88     ,fh1EvtSelection(0)
89     ,fh1VertexNContributors(0)
90     ,fh1VertexZ(0)
91     ,fh1Xsec(0)
92     ,fh1Trials(0)
93     ,fh1PtHard(0)
94     ,fh1PtHardTrials(0)
95     ,fh2EtaJet(0)
96     ,fh2PhiJet(0)
97     ,fh2PtJet(0)
98     ,fh1PtJet(0)
99     ,fh2NtracksJet(0)
100     ,fProNtracksJet(0)
101     ,fh2EtaTrack(0)
102     ,fh2PhiTrack(0)
103     ,fh2PtTrack(0)
104     ,fh2FF(0)
105     ,fh2Ksi(0)
106     ,fh2DelEta(0)
107     ,fh2DelPhi(0)
108     ,fh2DelR(0)
109     
110     ,fh1PtLeadingJet(0)
111     ,fh2NtracksLeadingJet(0)
112     ,fProNtracksLeadingJet(0)
113     ,fh2DelR80pcNch(0)
114     ,fProDelR80pcNch(0)
115     ,fh2DelR80pcPt(0)
116     ,fProDelR80pcPt(0)
117     ,fh2AreaCh(0)
118     ,fProAreaCh(0)
119     ,fh3PtDelRNchSum(0)
120     ,fh3PtDelRPtSum(0)
121     ,fProDiffJetShape(0)
122     ,fProIntJetShape(0)
123     
124     ,fh1PtSumInJetConeUE(0)
125     ,fh2NtracksLeadingJetUE(0)
126     ,fProNtracksLeadingJetUE(0)
127     ,fh2DelR80pcNchUE(0)
128     ,fProDelR80pcNchUE(0)
129     ,fh2DelR80pcPtUE(0)
130     ,fProDelR80pcPtUE(0)
131     ,fh2AreaChUE(0)
132     ,fProAreaChUE(0)
133     ,fh3PtDelRNchSumUE(0)
134     ,fh3PtDelRPtSumUE(0)
135     ,fProDiffJetShapeUE(0)
136     ,fProIntJetShapeUE(0)
137     
138     ,fh1CorrJetPt(0)
139     ,fh2CorrPtTrack1(0)
140     ,fh2CorrFF1(0)
141     ,fh2CorrKsi1(0)
142     ,fh2CorrjT1(0)
143     ,fh1JetPtvsTrkSum(0)
144     
145 {
146   for(Int_t ii=0; ii<13; ii++){
147     if(ii<6){
148       fh2CorrPt1Pt2[ii]   = NULL;
149       fh2CorrZ1Z2[ii]     = NULL;
150       fh2CorrKsi1Ksi2[ii] = NULL;
151       fh2CorrjT1jT2[ii]   = NULL;
152     }
153     
154     fProDelRNchSum[ii]    = NULL;  
155     fProDelRPtSum[ii]     = NULL;
156     fProDiffJetShapeA[ii] = NULL;
157     fProIntJetShapeA[ii]  = NULL;
158     
159     fProDelRNchSumUE[ii]    = NULL;  
160     fProDelRPtSumUE[ii]     = NULL;
161     fProDiffJetShapeAUE[ii] = NULL;
162     fProIntJetShapeAUE[ii]  = NULL;
163   }//ii loop
164   // default constructor
165 }
166 //_________________________________________________________________________________//
167
168 AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties(const char *name) 
169   : AliAnalysisTaskSE(name)
170   ,fESD(0)
171   ,fAOD(0)
172   ,fAODJets(0)  
173   ,fAODExtension(0)
174   ,fNonStdFile("")
175   ,fBranchJets("jets")
176   ,fTrackType(kTrackAOD)
177   ,fJetRejectType(0)
178   ,fRejectPileup(1)
179   ,fUseAODInputJets(kTRUE)
180   ,fFilterMask(0)
181   ,fUsePhysicsSelection(kTRUE)
182   ,fMaxVertexZ(10)
183   ,fNContributors(2)
184   ,fTrackPtCut(0)
185   ,fTrackEtaMin(0)
186   ,fTrackEtaMax(0)
187   ,fJetPtCut(0)
188   ,fJetEtaMin(0)
189   ,fJetEtaMax(0)
190   ,fAvgTrials(0)
191   ,fJetRadius(0.4)
192   ,fJetList(0)
193   ,fTrackList(0)
194   ,fTrackListUE(0)
195   ,fTrackListJet(0)
196   ,fCommonHistList(0)
197   ,fh1EvtSelection(0)
198   ,fh1VertexNContributors(0)
199   ,fh1VertexZ(0)
200   ,fh1Xsec(0)
201   ,fh1Trials(0)
202   ,fh1PtHard(0)
203   ,fh1PtHardTrials(0)
204   ,fh2EtaJet(0)
205   ,fh2PhiJet(0)
206   ,fh2PtJet(0)
207   ,fh1PtJet(0)
208   ,fh2NtracksJet(0)
209   ,fProNtracksJet(0)
210   ,fh2EtaTrack(0)
211   ,fh2PhiTrack(0)
212   ,fh2PtTrack(0)
213   ,fh2FF(0)
214   ,fh2Ksi(0)
215   ,fh2DelEta(0)
216   ,fh2DelPhi(0)
217   ,fh2DelR(0)
218   
219   ,fh1PtLeadingJet(0)
220   ,fh2NtracksLeadingJet(0)
221   ,fProNtracksLeadingJet(0)
222   ,fh2DelR80pcNch(0)
223   ,fProDelR80pcNch(0)
224   ,fh2DelR80pcPt(0)
225   ,fProDelR80pcPt(0)
226   ,fh2AreaCh(0)
227   ,fProAreaCh(0)
228   ,fh3PtDelRNchSum(0)
229   ,fh3PtDelRPtSum(0)
230   ,fProDiffJetShape(0)
231   ,fProIntJetShape(0)
232   
233   ,fh1PtSumInJetConeUE(0)
234   ,fh2NtracksLeadingJetUE(0)
235   ,fProNtracksLeadingJetUE(0)
236   ,fh2DelR80pcNchUE(0)
237   ,fProDelR80pcNchUE(0)
238   ,fh2DelR80pcPtUE(0)
239   ,fProDelR80pcPtUE(0)
240   ,fh2AreaChUE(0)
241   ,fProAreaChUE(0)
242   ,fh3PtDelRNchSumUE(0)
243   ,fh3PtDelRPtSumUE(0)
244   ,fProDiffJetShapeUE(0)
245   ,fProIntJetShapeUE(0)
246
247   ,fh1CorrJetPt(0)
248   ,fh2CorrPtTrack1(0)
249   ,fh2CorrFF1(0)
250   ,fh2CorrKsi1(0)
251   ,fh2CorrjT1(0)
252   ,fh1JetPtvsTrkSum(0)
253 {
254   for(Int_t ii=0; ii<13; ii++){
255     if(ii<6){
256       fh2CorrPt1Pt2[ii]   = NULL;
257       fh2CorrZ1Z2[ii]     = NULL;
258       fh2CorrKsi1Ksi2[ii] = NULL;
259       fh2CorrjT1jT2[ii]   = NULL;
260     }
261     
262     fProDelRNchSum[ii]    = NULL;  
263     fProDelRPtSum[ii]     = NULL;  
264     fProDiffJetShapeA[ii] = NULL;
265     fProIntJetShapeA[ii]  = NULL;
266
267     fProDelRNchSumUE[ii]    = NULL;  
268     fProDelRPtSumUE[ii]     = NULL;  
269     fProDiffJetShapeAUE[ii] = NULL;
270     fProIntJetShapeAUE[ii]  = NULL;
271   }//ii loop
272   // constructor
273   DefineOutput(1,TList::Class());
274 }
275 //_________________________________________________________________________________//
276
277 AliAnalysisTaskJetProperties::~AliAnalysisTaskJetProperties()
278 {
279   // destructor
280   if(fJetList)   delete fJetList;
281   if(fTrackList)  delete fTrackList;
282   if(fTrackListUE) delete fTrackListUE;
283   if(fTrackListJet) delete fTrackListJet;
284 }
285 //_________________________________________________________________________________//
286
287 Bool_t AliAnalysisTaskJetProperties::Notify()
288 {
289   //
290   // Implemented Notify() to read the cross sections
291   // and number of trials from pyxsec.root
292   // (taken from AliAnalysisTaskJetSpectrum2)
293   // 
294   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Notify()");
295   
296   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
297   Float_t xsection = 0;
298   Float_t ftrials  = 1;
299   
300   fAvgTrials = 1;
301   if(tree){
302     TFile *curfile = tree->GetCurrentFile();
303     if (!curfile) {
304       Error("Notify","No current file");
305       return kFALSE;
306     }
307     if(!fh1Xsec||!fh1Trials){
308       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
309       return kFALSE;
310     }
311     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
312     fh1Xsec->Fill("<#sigma>",xsection);
313     // construct a poor man average trials 
314     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
315     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
316   }
317   return kTRUE;
318 }
319 //_________________________________________________________________________________//
320
321 void AliAnalysisTaskJetProperties::UserCreateOutputObjects()
322 {
323   //(Here, event selection part is taken from AliAnalysisTaskFramentationFunction)
324   // create output objects
325   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserCreateOutputObjects()");
326   
327   // create list of tracks and jets   
328   fJetList = new TList();
329   fJetList->SetOwner(kFALSE);
330   fTrackList = new TList();
331   fTrackList->SetOwner(kFALSE);
332   fTrackListUE = new TList();
333   fTrackListUE->SetOwner(kFALSE);
334   fTrackListJet = new TList();
335   fTrackListJet->SetOwner(kFALSE);
336   
337   // Create histograms / output container
338   OpenFile(1);
339   fCommonHistList = new TList();
340   fCommonHistList->SetOwner();
341   
342   Bool_t oldStatus = TH1::AddDirectoryStatus();
343   TH1::AddDirectory(kFALSE);
344   // Histograms/TProfiles       
345   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
346   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
347   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
348   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
349   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
350   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
351   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"pile up: rejected");
352   
353   
354   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
355   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
356   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
357   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
358   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
359   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
360   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
361   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
362   
363   
364   Int_t kNbinsPtSlice=20; Float_t xMinPtSlice=0.0; Float_t xMaxPtSlice=100.0;
365   Int_t kNbinsEta=40;     Float_t xMinEta=-2.0;    Float_t xMaxEta=2.0;
366   Int_t kNbinsPhi=90;     Float_t xMinPhi=0.0;     Float_t xMaxPhi=TMath::TwoPi();
367   Int_t kNbinsPt=250;     Float_t xMinPt=0.0;      Float_t xMaxPt=250.0;
368   Int_t kNbinsNtracks=50; Float_t xMinNtracks=0.0; Float_t xMaxNtracks=50.0;
369   Int_t kNbinsFF=200;     Float_t xMinFF=-0.05;    Float_t xMaxFF=1.95;
370   Int_t kNbinsKsi=80;     Float_t xMinKsi=0.;      Float_t xMaxKsi = 8.0;
371   Int_t kNbinsDelR1D=50;  Float_t xMinDelR1D=0.0;  Float_t xMaxDelR1D=1.0;
372   Int_t kNbinsjT=100;     Float_t xMinjT=0.0;      Float_t xMaxjT=10.0;
373   
374   fh2EtaJet      = new TH2F("EtaJet","EtaJet", 
375                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
376                             kNbinsEta,     xMinEta,     xMaxEta);
377   fh2PhiJet      = new TH2F("PhiJet","PhiJet",
378                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
379                             kNbinsPhi,     xMinPhi,     xMaxPhi);
380   fh2PtJet       = new TH2F("PtJet","PtJet",
381                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
382                             kNbinsPt,      xMinPt,      xMaxPt);
383   fh1PtJet       = new TH1F("PtJet1D","PtJet1D",
384                             kNbinsPt,      xMinPt,      xMaxPt);
385   fh2NtracksJet  = new TH2F("NtracksJet","NtracksJet",
386                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
387                             kNbinsNtracks, xMinNtracks, xMaxNtracks);
388   fProNtracksJet = new TProfile("AvgNoOfTracksJet","AvgNoOfTracksJet",
389                                 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
390                                 xMinNtracks, xMaxNtracks);
391   fh2EtaTrack    = new TH2F("EtaTrack","EtaTrack", 
392                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
393                             kNbinsEta,     xMinEta,     xMaxEta);
394   fh2PhiTrack    = new TH2F("PhiTrack","PhiTrack",
395                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
396                             kNbinsPhi,     xMinPhi,     xMaxPhi);
397   fh2PtTrack     = new TH2F("PtTrack","PtTrack",
398                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
399                             2*kNbinsPt,      xMinPt,      xMaxPt);
400   fh2FF          = new TH2F("FF","FF",
401                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
402                             kNbinsFF,      xMinFF,      xMaxFF);
403   fh2Ksi         = new TH2F("Ksi","Ksi",
404                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
405                             kNbinsKsi,      xMinKsi,      xMaxKsi);
406   fh2DelEta      = new TH2F("DelEta","DelEta", 
407                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
408                             kNbinsEta,     xMinEta,     xMaxEta);
409   fh2DelPhi      = new TH2F("DelPhi","DelPhi",
410                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
411                             kNbinsPhi,     xMinPhi,     xMaxPhi);
412   fh2DelR        = new TH2F("DelR","DelR",
413                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
414                             kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D);
415   
416
417   
418   Int_t kNbinsDelR=100;      Float_t xMinDelR=0.0;      Float_t xMaxDelR=1.0;
419   Int_t kNbinsPtSliceJS=100; Float_t xMinPtSliceJS=0.0; Float_t xMaxPtSliceJS=250.0;
420   
421   fh1PtLeadingJet       = new TH1F("PtLeadingJet","PtLeadingJet",
422                                    kNbinsPt,      xMinPt,      xMaxPt);
423   fh2NtracksLeadingJet  = new TH2F("NtracksLeadingJet","NtracksLeadingJet",
424                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
425                                    kNbinsNtracks,   xMinNtracks,   xMaxNtracks);
426   fProNtracksLeadingJet = new TProfile("AvgNoOfTracksLeadingJet","AvgNoOfTracksLeadingJet",
427                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
428                                        xMinNtracks, xMaxNtracks);
429   fh2DelR80pcNch        = new TH2F("delR80pcNch","delR80pcNch",
430                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
431                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
432   fProDelR80pcNch       = new TProfile("AvgdelR80pcNch","AvgdelR80pcNch",
433                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
434                                        xMinDelR,  xMaxDelR);
435   fh2DelR80pcPt         = new TH2F("delR80pcPt","delR80pcPt",
436                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
437                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
438   fProDelR80pcPt        = new TProfile("AvgdelR80pcPt","AvgdelR80pcPt",
439                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
440                                        xMinDelR,  xMaxDelR);
441   fh2AreaCh             = new TH2F("Area","Area",
442                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
443                                    72, 0.0, TMath::Pi());
444   fProAreaCh            = new TProfile("AvgArea","AvgArea",
445                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
446                                    0.0, TMath::Pi());
447   fh3PtDelRNchSum       = new TH3F("PtDelRNchSum","PtDelRNchSum",
448                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
449                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
450                                    kNbinsNtracks, xMinNtracks, xMaxNtracks);
451   fh3PtDelRPtSum        = new TH3F("PtDelRPtSum","PtDelRPtSum",
452                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
453                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
454                                    kNbinsPt,        xMinPt,        xMaxPt);
455   fProDiffJetShape      = new TProfile("DiffJetShape","DiffJetShape",
456                                   10,0.0,1.0,0.0,250.0);
457   fProIntJetShape       = new TProfile("IntJetShape","IntJetShape",
458                                   10,0.0,1.0,0.0,250.0);
459
460
461   fh1PtSumInJetConeUE       = new TH1F("PtSumInJetConeUE","PtSumInJetConeUE",
462                                        500,      0.,      100.);
463   fh2NtracksLeadingJetUE  = new TH2F("NtracksLeadingJetUE","NtracksLeadingJetUE",
464                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
465                                    kNbinsNtracks,   xMinNtracks,   xMaxNtracks);
466   fProNtracksLeadingJetUE = new TProfile("AvgNoOfTracksLeadingJetUE","AvgNoOfTracksLeadingJetUE",
467                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
468                                        xMinNtracks, xMaxNtracks);
469   fh2DelR80pcNchUE        = new TH2F("delR80pcNchUE","delR80pcNchUE",
470                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
471                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
472   fProDelR80pcNchUE       = new TProfile("AvgdelR80pcNchUE","AvgdelR80pcNchUE",
473                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
474                                        xMinDelR,  xMaxDelR);
475   fh2DelR80pcPtUE         = new TH2F("delR80pcPtUE","delR80pcPtUE",
476                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
477                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
478   fProDelR80pcPtUE        = new TProfile("AvgdelR80pcPtUE","AvgdelR80pcPtUE",
479                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
480                                        xMinDelR,  xMaxDelR);
481   fh2AreaChUE             = new TH2F("AreaUE","AreaUE",
482                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
483                                    72, 0.0, TMath::Pi());
484   fProAreaChUE            = new TProfile("AvgAreaUE","AvgAreaUE",
485                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
486                                    0.0, TMath::Pi());
487   fh3PtDelRNchSumUE       = new TH3F("PtDelRNchSumUE","PtDelRNchSumUE",
488                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
489                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
490                                    kNbinsNtracks, xMinNtracks, xMaxNtracks);
491   fh3PtDelRPtSumUE        = new TH3F("PtDelRPtSumUE","PtDelRPtSumUE",
492                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
493                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
494                                    kNbinsPt,        xMinPt,        xMaxPt);
495   fProDiffJetShapeUE      = new TProfile("DiffJetShapeUE","DiffJetShapeUE",
496                                   10,0.0,1.0,0.0,250.0);
497   fProIntJetShapeUE       = new TProfile("IntJetShapeUE","IntJetShapeUE",
498                                   10,0.0,1.0,0.0,250.0);
499
500   
501   fh1CorrJetPt    = new TH1F("JetPtCorr","JetPtCorr", 
502                              kNbinsPt,      xMinPt,      xMaxPt);
503   fh2CorrPtTrack1 = new TH2F("TrackPtCorr", "TrackPtCorr",
504                              kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
505                              2*kNbinsPt,      xMinPt,      xMaxPt);
506   fh2CorrFF1      = new TH2F("FFCorr","FFCorr",
507                              kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
508                              kNbinsFF,      xMinFF,      xMaxFF);
509   fh2CorrKsi1     = new TH2F("KsiCorr","KsiCorr",
510                              kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
511                              kNbinsKsi,      xMinKsi,      xMaxKsi);
512   fh2CorrjT1      = new TH2F("jTCorr","jTCorr",
513                              kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
514                              kNbinsjT,      xMinjT,      xMaxjT);
515
516   fh1JetPtvsTrkSum = new TH1F("diffJetPt_TrkPtSum","diffJetPt_TrkPtSum",
517                               500, -250., 250.);
518   TString title;
519   for(Int_t ii=0; ii<13; ii++){
520     if(ii==0)title = "_JetPt20to25";
521     if(ii==1)title = "_JetPt25to30";
522     if(ii==2)title = "_JetPt20to30";
523
524     if(ii==3)title = "_JetPt30to35";
525     if(ii==4)title = "_JetPt35to40";
526     if(ii==5)title = "_JetPt30to40";
527
528     if(ii==6)title = "_JetPt40to50";
529     if(ii==7)title = "_JetPt50to60";
530     if(ii==8)title = "_JetPt40to60";
531
532     if(ii==9) title = "_JetPt60to70";
533     if(ii==10)title = "_JetPt70to80";
534     if(ii==11)title = "_JetPt60to80";
535
536     if(ii==12)title = "_JetPt80to100";
537     fProDelRNchSum[ii] = new TProfile(Form("AvgNchSumDelR%s",title.Data()),Form("AvgNchSumDelR%s",title.Data()),
538                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
539                                       xMinNtracks, xMaxNtracks);
540     
541     fProDelRPtSum[ii]  = new TProfile(Form("AvgPtSumDelR%s",title.Data()),Form("AvgPtSumDelR%s",title.Data()),
542                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
543                                       xMinPt,      xMaxPt);
544     fProDelRNchSum[ii] ->GetXaxis()->SetTitle("R");  
545     fProDelRNchSum[ii] ->GetYaxis()->SetTitle("<NchSum>");
546     fProDelRPtSum[ii]  ->GetXaxis()->SetTitle("R");  
547     fProDelRPtSum[ii]  ->GetYaxis()->SetTitle("<PtSum>");
548
549     fProDiffJetShapeA[ii] = new TProfile(Form("DiffJetShape%s",title.Data()),Form("DiffJetShape%s",title.Data()),
550                                          10,0.0,1.0,0.0,250.0);
551     fProIntJetShapeA[ii]  = new TProfile(Form("IntJetShape%s",title.Data()),Form("IntJetShape%s",title.Data()),
552                                          10,0.0,1.0,0.0,250.0);
553
554     fProDiffJetShapeA[ii]->GetXaxis()->SetTitle("R"); 
555     fProDiffJetShapeA[ii]->GetYaxis()->SetTitle("Diff jet shape"); 
556     fProIntJetShapeA[ii]->GetXaxis()->SetTitle("R");  
557     fProIntJetShapeA[ii]->GetYaxis()->SetTitle("Integrated jet shape");  
558     
559     fCommonHistList->Add(fProDelRNchSum[ii]);
560     fCommonHistList->Add(fProDelRPtSum[ii]);
561     fCommonHistList->Add(fProDiffJetShapeA[ii]);
562     fCommonHistList->Add(fProIntJetShapeA[ii]);
563
564     fProDelRNchSumUE[ii] = new TProfile(Form("AvgNchSumDelR%sUE",title.Data()),Form("AvgNchSumDelR%sUE",title.Data()),
565                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
566                                       xMinNtracks, xMaxNtracks);
567     
568     fProDelRPtSumUE[ii]  = new TProfile(Form("AvgPtSumDelR%sUE",title.Data()),Form("AvgPtSumDelR%sUE",title.Data()),
569                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
570                                       xMinPt,      xMaxPt);
571     fProDelRNchSumUE[ii] ->GetXaxis()->SetTitle("R");  
572     fProDelRNchSumUE[ii] ->GetYaxis()->SetTitle("<N_{ch}^{Sum,UE}>");
573     fProDelRPtSumUE[ii]  ->GetXaxis()->SetTitle("R");  
574     fProDelRPtSumUE[ii]  ->GetYaxis()->SetTitle("<p_{T}^{sum, UE}>");
575
576     fProDiffJetShapeAUE[ii] = new TProfile(Form("DiffJetShape%sUE",title.Data()),Form("DiffJetShape%sUE",title.Data()),
577                                          10,0.0,1.0,0.0,250.0);
578     fProIntJetShapeAUE[ii]  = new TProfile(Form("IntJetShape%sUE",title.Data()),Form("IntJetShape%sUE",title.Data()),
579                                          10,0.0,1.0,0.0,250.0);
580
581     fProDiffJetShapeAUE[ii]->GetXaxis()->SetTitle("R"); 
582     fProDiffJetShapeAUE[ii]->GetYaxis()->SetTitle("Diff jet shape UE"); 
583     fProIntJetShapeAUE[ii]->GetXaxis()->SetTitle("R");  
584     fProIntJetShapeAUE[ii]->GetYaxis()->SetTitle("Integrated jet shape UE");  
585     
586     fCommonHistList->Add(fProDelRNchSumUE[ii]);
587     fCommonHistList->Add(fProDelRPtSumUE[ii]);
588     fCommonHistList->Add(fProDiffJetShapeAUE[ii]);
589     fCommonHistList->Add(fProIntJetShapeAUE[ii]);
590   }//ii loop
591   
592   fh2EtaJet     ->GetXaxis()->SetTitle("JetPt");  fh2EtaJet     ->GetYaxis()->SetTitle("JetEta");
593   fh2PhiJet     ->GetXaxis()->SetTitle("JetPt");  fh2PhiJet     ->GetYaxis()->SetTitle("JetPhi");
594   fh2PtJet      ->GetXaxis()->SetTitle("JetPt");  fh2PtJet      ->GetYaxis()->SetTitle("JetPt");
595   fh1PtJet      ->GetXaxis()->SetTitle("JetPt");  fh1PtJet      ->GetYaxis()->SetTitle("#jets");
596   fh2NtracksJet ->GetXaxis()->SetTitle("JetPt");  fh2NtracksJet ->GetYaxis()->SetTitle("#tracks");
597   fProNtracksJet->GetXaxis()->SetTitle("JetPt");  fProNtracksJet->GetYaxis()->SetTitle("AgvNtracks");
598   fh2EtaTrack   ->GetXaxis()->SetTitle("JetPt");  fh2EtaTrack   ->GetYaxis()->SetTitle("TrackEta");
599   fh2PhiTrack   ->GetXaxis()->SetTitle("JetPt");  fh2PhiTrack   ->GetYaxis()->SetTitle("TrackPhi");
600   fh2PtTrack    ->GetXaxis()->SetTitle("JetPt");  fh2PtTrack    ->GetYaxis()->SetTitle("TrackPt");
601   fh2FF         ->GetXaxis()->SetTitle("JetPt");  fh2FF         ->GetYaxis()->SetTitle("FF");
602   fh2Ksi        ->GetXaxis()->SetTitle("JetPt");  fh2Ksi        ->GetYaxis()->SetTitle("Ksi");
603   fh2DelEta     ->GetXaxis()->SetTitle("JetPt");  fh2DelEta     ->GetYaxis()->SetTitle("DelEta");
604   fh2DelPhi     ->GetXaxis()->SetTitle("JetPt");  fh2DelPhi     ->GetYaxis()->SetTitle("DelPhi");
605   fh2DelR       ->GetXaxis()->SetTitle("JetPt");  fh2DelR       ->GetYaxis()->SetTitle("DelR");
606
607   fh1PtLeadingJet       ->GetXaxis()->SetTitle("JetPt");  fh1PtLeadingJet       ->GetYaxis()->SetTitle("#leading jets");
608   fh2NtracksLeadingJet  ->GetXaxis()->SetTitle("JetPt");  fh2NtracksLeadingJet  ->GetYaxis()->SetTitle("#tracks leading jet");
609   fProNtracksLeadingJet ->GetXaxis()->SetTitle("JetPt");  fProNtracksLeadingJet ->GetYaxis()->SetTitle("AvgNtracks leading jet");
610   fh2DelR80pcNch        ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcNch        ->GetYaxis()->SetTitle("R containing 80% of tracks");
611   fProDelR80pcNch       ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcNch       ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
612   fh2DelR80pcPt         ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcPt         ->GetYaxis()->SetTitle("R containing 80% of pT");
613   fProDelR80pcPt        ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcPt        ->GetYaxis()->SetTitle("<R> containing 80% of pT");
614   fh2AreaCh             ->GetXaxis()->SetTitle("JetPt");  fh2AreaCh             ->GetYaxis()->SetTitle("Jet area");
615   fProAreaCh            ->GetXaxis()->SetTitle("JetPt");  fProAreaCh            ->GetYaxis()->SetTitle("<jet area>");
616   fh3PtDelRNchSum       ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRNchSum       ->GetYaxis()->SetTitle("R");  fh3PtDelRNchSum->GetZaxis()->SetTitle("NchSum");
617   fh3PtDelRPtSum        ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRPtSum        ->GetYaxis()->SetTitle("R");  fh3PtDelRPtSum ->GetZaxis()->SetTitle("PtSum");
618   fProDiffJetShape->GetXaxis()->SetTitle("R"); 
619   fProDiffJetShape->GetYaxis()->SetTitle("Diff jet shape"); 
620   fProIntJetShape->GetXaxis()->SetTitle("R");  
621   fProIntJetShape->GetYaxis()->SetTitle("Integrated jet shape");  
622
623   fh1PtSumInJetConeUE     ->GetXaxis()->SetTitle("p_{T}^{sum, UE}(in cone R)");  fh1PtSumInJetConeUE     ->GetYaxis()->SetTitle("#leading jets");
624   fh2NtracksLeadingJetUE  ->GetXaxis()->SetTitle("JetPt");  fh2NtracksLeadingJetUE  ->GetYaxis()->SetTitle("#tracks UE");
625   fProNtracksLeadingJetUE ->GetXaxis()->SetTitle("JetPt");  fProNtracksLeadingJetUE ->GetYaxis()->SetTitle("AvgNtracks UE");
626   fh2DelR80pcNchUE        ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcNchUE        ->GetYaxis()->SetTitle("R containing 80% of tracks");
627   fProDelR80pcNchUE       ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcNchUE       ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
628   fh2DelR80pcPtUE         ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcPtUE         ->GetYaxis()->SetTitle("R containing 80% of pT");
629   fProDelR80pcPtUE        ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcPtUE        ->GetYaxis()->SetTitle("<R> containing 80% of pT");
630   fh2AreaChUE             ->GetXaxis()->SetTitle("JetPt");  fh2AreaChUE             ->GetYaxis()->SetTitle("UE area");
631   fProAreaChUE            ->GetXaxis()->SetTitle("JetPt");  fProAreaChUE            ->GetYaxis()->SetTitle("<UE area>");
632   fh3PtDelRNchSumUE       ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRNchSumUE       ->GetYaxis()->SetTitle("R");  fh3PtDelRNchSumUE->GetZaxis()->SetTitle("NchSumUE");
633   fh3PtDelRPtSumUE        ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRPtSumUE        ->GetYaxis()->SetTitle("R");  fh3PtDelRPtSumUE ->GetZaxis()->SetTitle("PtSumUE");
634   fProDiffJetShapeUE->GetXaxis()->SetTitle("R"); 
635   fProDiffJetShapeUE->GetYaxis()->SetTitle("Diff jet shape UE"); 
636   fProIntJetShapeUE->GetXaxis()->SetTitle("R");  
637   fProIntJetShapeUE->GetYaxis()->SetTitle("Integrated jet shape UE");  
638   
639   fh1CorrJetPt    ->GetXaxis()->SetTitle("JetPt");  fh1CorrJetPt    ->GetYaxis()->SetTitle("#jets");
640   fh2CorrPtTrack1 ->GetXaxis()->SetTitle("JetPt");  fh2CorrPtTrack1 ->GetYaxis()->SetTitle("pt_track");
641   fh2CorrFF1      ->GetXaxis()->SetTitle("JetPt");  fh2CorrFF1      ->GetYaxis()->SetTitle("z_track");
642   fh2CorrKsi1     ->GetXaxis()->SetTitle("JetPt");  fh2CorrKsi1     ->GetYaxis()->SetTitle("ksi_track");
643   fh2CorrjT1      ->GetXaxis()->SetTitle("JetPt");  fh2CorrjT1      ->GetYaxis()->SetTitle("jt_track");
644   fh1JetPtvsTrkSum->GetXaxis()->SetTitle("JetPt-TrkPtSum");
645
646   fCommonHistList->Add(fh1EvtSelection);
647   fCommonHistList->Add(fh1VertexNContributors);
648   fCommonHistList->Add(fh1VertexZ);
649   fCommonHistList->Add(fh1Xsec);
650   fCommonHistList->Add(fh1Trials);
651   fCommonHistList->Add(fh1PtHard);
652   fCommonHistList->Add(fh1PtHardTrials);
653   fCommonHistList->Add(fh2EtaJet);
654   fCommonHistList->Add(fh2PhiJet);
655   fCommonHistList->Add(fh2PtJet);
656   fCommonHistList->Add(fh1PtJet);
657   fCommonHistList->Add(fh2NtracksJet);
658   fCommonHistList->Add(fProNtracksJet);
659   fCommonHistList->Add(fh2EtaTrack);
660   fCommonHistList->Add(fh2PhiTrack); 
661   fCommonHistList->Add(fh2PtTrack); 
662   fCommonHistList->Add(fh2FF);  
663   fCommonHistList->Add(fh2Ksi);  
664   fCommonHistList->Add(fh2DelEta);       
665   fCommonHistList->Add(fh2DelPhi);   
666   fCommonHistList->Add(fh2DelR); 
667   
668   fCommonHistList->Add(fh1PtLeadingJet); 
669   fCommonHistList->Add(fh2NtracksLeadingJet); 
670   fCommonHistList->Add(fProNtracksLeadingJet); 
671   fCommonHistList->Add(fh2DelR80pcNch); 
672   fCommonHistList->Add(fProDelR80pcNch); 
673   fCommonHistList->Add(fh2DelR80pcPt); 
674   fCommonHistList->Add(fProDelR80pcPt); 
675   fCommonHistList->Add(fh2AreaCh); 
676   fCommonHistList->Add(fProAreaCh); 
677   fCommonHistList->Add(fh3PtDelRNchSum); 
678   fCommonHistList->Add(fh3PtDelRPtSum); 
679   fCommonHistList->Add(fProDiffJetShape);
680   fCommonHistList->Add(fProIntJetShape);
681
682   fCommonHistList->Add(fh1PtSumInJetConeUE);
683   fCommonHistList->Add(fh2NtracksLeadingJetUE); 
684   fCommonHistList->Add(fProNtracksLeadingJetUE); 
685   fCommonHistList->Add(fh2DelR80pcNchUE); 
686   fCommonHistList->Add(fProDelR80pcNchUE); 
687   fCommonHistList->Add(fh2DelR80pcPtUE); 
688   fCommonHistList->Add(fProDelR80pcPtUE); 
689   fCommonHistList->Add(fh2AreaChUE); 
690   fCommonHistList->Add(fProAreaChUE); 
691   fCommonHistList->Add(fh3PtDelRNchSumUE); 
692   fCommonHistList->Add(fh3PtDelRPtSumUE); 
693   fCommonHistList->Add(fProDiffJetShapeUE);
694   fCommonHistList->Add(fProIntJetShapeUE);
695
696   fCommonHistList->Add(fh1CorrJetPt);
697   fCommonHistList->Add(fh2CorrPtTrack1);
698   fCommonHistList->Add(fh2CorrFF1);
699   fCommonHistList->Add(fh2CorrKsi1);
700   fCommonHistList->Add(fh2CorrjT1);
701   fCommonHistList->Add(fh1JetPtvsTrkSum);
702
703       
704
705   TString titleCorr;
706   for(Int_t jj=0; jj<6; jj++){
707     if(jj == 0)titleCorr = "_JetPt10to20";
708     if(jj == 1)titleCorr = "_JetPt20to30";
709     if(jj == 2)titleCorr = "_JetPt30to40";
710     if(jj == 3)titleCorr = "_JetPt40to60";
711     if(jj == 4)titleCorr = "_JetPt60to80";
712     if(jj == 5)titleCorr = "_JetPt80to100";
713     
714     fh2CorrPt1Pt2[jj]   = new TH2F(Form("CorrPt1Pt2%s",titleCorr.Data()),Form("CorrPt1Pt2%s",titleCorr.Data()),
715                                    2*kNbinsPt,      xMinPt,      xMaxPt,
716                                    2*kNbinsPt,      xMinPt,      xMaxPt);
717     fh2CorrZ1Z2[jj]     = new TH2F(Form("CorrZ1Z2%s",titleCorr.Data()),Form("CorrZ1Z2%s",titleCorr.Data()),
718                                    kNbinsFF,      xMinFF,      xMaxFF,
719                                    kNbinsFF,      xMinFF,      xMaxFF);
720     fh2CorrKsi1Ksi2[jj] = new TH2F(Form("CorrKsi1Ksi2%s",titleCorr.Data()),Form("CorrKsi1Ksi2%s",titleCorr.Data()),
721                                    kNbinsKsi,      xMinKsi,      xMaxKsi,
722                                    kNbinsKsi,      xMinKsi,      xMaxKsi);
723     fh2CorrjT1jT2[jj]   = new TH2F(Form("CorrjT1jT2%s",titleCorr.Data()),Form("CorrjT1jT2%s",titleCorr.Data()),
724                                    kNbinsjT,      xMinjT,      xMaxjT,
725                                    kNbinsjT,      xMinjT,      xMaxjT);
726     
727     
728     fh2CorrPt1Pt2[jj]   ->GetXaxis()->SetTitle("pt_track1");
729     fh2CorrZ1Z2[jj]     ->GetXaxis()->SetTitle("z_track1");
730     fh2CorrKsi1Ksi2[jj] ->GetXaxis()->SetTitle("ksi_track1");
731     fh2CorrjT1jT2[jj]   ->GetXaxis()->SetTitle("jt_track1");
732     
733     fh2CorrPt1Pt2[jj]   ->GetYaxis()->SetTitle("pt_track2");
734     fh2CorrZ1Z2[jj]     ->GetYaxis()->SetTitle("z_track2");
735     fh2CorrKsi1Ksi2[jj] ->GetYaxis()->SetTitle("ksi_track2");
736     fh2CorrjT1jT2[jj]   ->GetYaxis()->SetTitle("jt_track2");
737     
738     fCommonHistList->Add(fh2CorrPt1Pt2[jj]);
739     fCommonHistList->Add(fh2CorrZ1Z2[jj]);
740     fCommonHistList->Add(fh2CorrKsi1Ksi2[jj]);
741     fCommonHistList->Add(fh2CorrjT1jT2[jj]);
742   }//jj loop
743
744
745
746   // =========== Switch on Sumw2 for all histos ===========
747   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
748     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
749     if (h1) h1->Sumw2();
750     else{
751       TProfile *hPro = dynamic_cast<TProfile*>(fCommonHistList->At(i));
752       if(hPro) hPro->Sumw2();
753     }
754   }
755   
756   TH1::AddDirectory(oldStatus);
757   PostData(1, fCommonHistList);
758 }
759 //_________________________________________________________________________________//
760
761 void AliAnalysisTaskJetProperties::Init()
762 {
763   // Initialization
764   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Init()");
765
766 }
767 //_________________________________________________________________________________//
768
769 void AliAnalysisTaskJetProperties::UserExec(Option_t *) 
770 {
771   // Main loop
772   // Called for each event
773   if(fDebug > 2) Printf("AliAnalysisTaskJetProperties::UserExec()");
774   
775   //if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
776   // Trigger selection
777   
778   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
779     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
780   if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
781     if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
782       fh1EvtSelection->Fill(1.);
783       if (fDebug > 2 ) Printf(" Trigger Selection: event REJECTED ... ");
784       PostData(1, fCommonHistList);
785       return;
786     }
787   }
788   
789   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
790   if(!fESD){
791     if(fDebug>2) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
792   }
793   
794   fMCEvent = MCEvent();
795   if(!fMCEvent){
796     if(fDebug>2) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
797   }
798   
799   // get AOD event from input/ouput
800     TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
801     if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
802       fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
803       if(fUseAODInputJets) fAODJets = fAOD;
804       if (fDebug > 2)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
805     }
806     else {
807       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
808       if( handler && handler->InheritsFrom("AliAODHandler") ) {
809         fAOD = ((AliAODHandler*)handler)->GetAOD();
810         fAODJets = fAOD;
811         if (fDebug > 2)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
812       }
813     }
814     
815     if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
816     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
817     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
818       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
819       if (fDebug > 2)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
820     }
821     }
822     
823     if(fNonStdFile.Length()!=0){
824       // case we have an AOD extension - fetch the jets from the extended output
825       
826       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
827       fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
828       if(!fAODExtension){
829         if(fDebug>2)Printf("AODExtension not found for %s",fNonStdFile.Data());
830       }
831     }
832   
833     if(!fAOD){
834       Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
835       return;
836     }
837     if(!fAODJets){
838       Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
839       return;
840     }
841     
842     
843     // *** vertex cut ***
844     AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
845     Int_t nTracksPrim = primVtx->GetNContributors();
846     fh1VertexNContributors->Fill(nTracksPrim);
847     
848   if (fDebug > 2) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
849   //if(!nTracksPrim){
850   if(nTracksPrim<fNContributors){
851     if (fDebug > 2) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
852     fh1EvtSelection->Fill(3.);
853     PostData(1, fCommonHistList);
854     return;
855   }
856   
857   if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
858     if (fDebug > 2) Printf("%s:%d SPD pileup : event REJECTED...",(char*)__FILE__,__LINE__);
859     fh1EvtSelection->Fill(6.);
860     PostData(1, fCommonHistList);
861     return; 
862   }//pile up rejection
863   
864
865   fh1VertexZ->Fill(primVtx->GetZ());
866   
867   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
868     if (fDebug > 2) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
869     fh1EvtSelection->Fill(4.);
870     PostData(1, fCommonHistList);
871     return; 
872   }
873   
874   TString primVtxName(primVtx->GetName());
875   
876   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
877     if (fDebug > 2) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
878     fh1EvtSelection->Fill(5.);
879     PostData(1, fCommonHistList);
880     return;
881   }
882   
883   if (fDebug > 2) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
884   fh1EvtSelection->Fill(0.);
885   //___ get MC information __________________________________________________________________
886   Double_t ptHard = 0.;
887   Double_t nTrials = 1; // trials for MC trigger weight for real data
888   
889   if(fMCEvent){
890     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
891     if(genHeader){
892       AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
893       AliGenHijingEventHeader*  hijingGenHeader = 0x0;
894       
895       if(pythiaGenHeader){
896         if(fDebug>2) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
897         nTrials = pythiaGenHeader->Trials();
898         ptHard  = pythiaGenHeader->GetPtHard();
899         
900         fh1PtHard->Fill(ptHard);
901         fh1PtHardTrials->Fill(ptHard,nTrials);
902         
903         
904       } else { // no pythia, hijing?
905         
906         if(fDebug>2) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
907         
908         hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
909         if(!hijingGenHeader){
910           Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
911         } else {
912           if(fDebug>2) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
913         }
914       }
915       
916       fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
917     }
918   }
919   
920   //___ fetch jets __________________________________________________________________________
921   Int_t nJ = GetListOfJets(fJetList);
922   Int_t nJets = 0;
923   if(nJ>=0) nJets = fJetList->GetEntries();
924   if(fDebug>2){
925     Printf("%s:%d Selected jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
926     if(nJ != nJets) Printf("%s:%d Mismatch Selected Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
927   }
928   FillJetProperties(fJetList);
929   FillJetShape(fJetList);
930   FillJetShapeUE(fJetList);
931   FillFFCorr(fJetList);
932   fJetList->Clear();
933   //Post output data.
934   PostData(1, fCommonHistList);
935 }//UserExec
936 //_________________________________________________________________________________//
937
938 void AliAnalysisTaskJetProperties::Terminate(Option_t *) 
939 {
940   // terminated
941   
942   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::Terminate() \n");
943 }
944 //_________________________________________________________________________________//
945
946 Int_t AliAnalysisTaskJetProperties::GetListOfJets(TList *list)
947 {
948   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
949   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::GetListOfJets() \n");
950   // fill list of jets selected according to type
951   if(!list){
952     if(fDebug>2) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
953     return -1;
954   }
955   
956   if(fBranchJets.Length()==0){
957     Printf("%s:%d no jet branch specified", (char*)__FILE__,__LINE__);
958     if(fDebug>2)fAOD->Print();
959     return 0;
960   }
961   TClonesArray *aodJets = 0; 
962   if(fBranchJets.Length())      aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchJets.Data()));
963   if(!aodJets)                  aodJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchJets.Data()));
964   if(fAODExtension&&!aodJets)   aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchJets.Data()));
965   
966   if(!aodJets){
967     if(fBranchJets.Length())Printf("%s:%d no jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchJets.Data());
968     if(fDebug>2)fAOD->Print();
969     return 0;
970   }
971   
972   Int_t nJets = 0;
973   for(Int_t ijet=0; ijet<aodJets->GetEntries(); ijet++){
974     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodJets->At(ijet));
975     if(!tmp) continue;
976     if( tmp->Pt() < fJetPtCut ) continue;
977     if( tmp->Eta() < fJetEtaMin || tmp->Eta() > fJetEtaMax)continue;
978     if(fJetRejectType==kReject1Track && tmp->GetRefTracks()->GetEntriesFast()==1)continue;//rejecting 1track jet if...
979     list->Add(tmp);
980     nJets++;
981   }//ij loop
982   list->Sort();
983   return nJets;
984 }
985 //_________________________________________________________________________________//
986
987 Int_t AliAnalysisTaskJetProperties::GetListOfJetTracks(TList* list, const AliAODJet* jet)
988 {
989   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
990   if(fDebug > 3) printf("AliAnalysisTaskJetProperties::GetListOfJetTracks() \n");
991   
992   // list of jet tracks from trackrefs
993   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
994   Int_t NoOfTracks=0;
995   for (Int_t itrack=0; itrack<nTracks; itrack++) {
996     if(fTrackType==kTrackUndef){
997       if(fDebug>3)Printf("%s:%d unknown track type %d in AOD", (char*)__FILE__,__LINE__,kTrackUndef);
998       return 0;
999     }//if
1000     else if(fTrackType == kTrackAOD){
1001       AliAODTrack* trackAod = dynamic_cast<AliAODTrack*>(jet->GetRefTracks()->At(itrack));
1002       if(!trackAod){
1003         AliError("expected ref track not found ");
1004         continue;
1005       }
1006       list->Add(trackAod);
1007       NoOfTracks++;
1008     }//if
1009     
1010     else if(fTrackType==kTrackAODMC){
1011       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(jet->GetRefTracks()->At(itrack));
1012       if(!trackmc){
1013         AliError("expected ref trackmc not found ");
1014         continue;
1015       }
1016       list->Add(trackmc);
1017       NoOfTracks++;
1018     }//if
1019     
1020     else if (fTrackType==kTrackKine){
1021       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
1022       if(!trackkine){
1023         AliError("expected ref trackkine not found ");
1024         continue;
1025       }
1026       list->Add(trackkine);
1027       NoOfTracks++;
1028     }//if
1029   }//itrack loop
1030   list->Sort();
1031   return NoOfTracks;
1032 }//initial loop
1033 //_________________________________________________________________________________//
1034
1035 void AliAnalysisTaskJetProperties::FillJetProperties(TList *jetlist){
1036   //filling up the histograms jets and tracks inside jet
1037   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillJetProperties() \n");
1038   
1039   for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
1040     Float_t JetPt;Float_t JetPhi;Float_t JetEta;Float_t JetE;
1041     AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
1042     if(!jet)continue;
1043     JetEta = jet->Eta();
1044     JetPhi = jet->Phi();
1045     JetPt  = jet->Pt();
1046     JetE   = jet->E();
1047     fh2EtaJet ->Fill(JetPt,JetEta);
1048     fh2PhiJet ->Fill(JetPt,JetPhi);
1049     fh2PtJet  ->Fill(JetPt,JetPt);
1050     fh1PtJet  ->Fill(JetPt);
1051     fTrackListJet->Clear();
1052     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1053     Int_t nJetTracks = 0;
1054     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1055     if(fDebug>2){
1056       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1057       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1058     }
1059     
1060     fh2NtracksJet  ->Fill(JetPt,fTrackListJet->GetEntries());
1061     fProNtracksJet ->Fill(JetPt,fTrackListJet->GetEntries());
1062     
1063     for (Int_t j =0; j< fTrackListJet->GetEntries(); j++){
1064       if(fTrackType==kTrackUndef)continue;
1065         Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1066         Float_t FF=-99.0;       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; 
1067         Float_t DelR=-99.0;     Float_t AreaJ=-99.0;   Float_t Ksi=-99.0;
1068         if(fTrackType==kTrackAOD){
1069           AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1070           if(!trackaod)continue;
1071           TrackEta = trackaod->Eta();
1072           TrackPhi = trackaod->Phi();
1073           TrackPt  = trackaod->Pt();
1074         }//if kTrackAOD
1075         else if(fTrackType==kTrackAODMC){
1076           AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1077           if(!trackmc)continue;
1078           TrackEta = trackmc->Eta();
1079           TrackPhi = trackmc->Phi();
1080           TrackPt  = trackmc->Pt();
1081         }//if kTrackAODMC
1082         else if(fTrackType==kTrackKine){
1083           AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1084           if(!trackkine)continue;
1085           TrackEta = trackkine->Eta();
1086           TrackPhi = trackkine->Phi();
1087           TrackPt  = trackkine->Pt();
1088         }//kTrackKine
1089         if(JetPt)FF = TrackPt/JetPt;
1090         if(FF)Ksi = TMath::Log(1./FF);
1091         DelEta      = TMath::Abs(JetEta - TrackEta);
1092         DelPhi      = TMath::Abs(JetPhi - TrackPhi);
1093         if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1094         DelR        = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1095         AreaJ       = TMath::Pi()*DelR*DelR;
1096         fh2EtaTrack ->Fill(JetPt,TrackEta);
1097         fh2PhiTrack ->Fill(JetPt,TrackPhi);
1098         fh2PtTrack  ->Fill(JetPt,TrackPt);
1099         fh2FF       ->Fill(JetPt,FF);
1100         fh2Ksi      ->Fill(JetPt,Ksi);
1101         fh2DelEta   ->Fill(JetPt,DelEta);
1102         fh2DelPhi   ->Fill(JetPt,DelPhi);
1103         fh2DelR     ->Fill(JetPt,DelR);
1104     }//track loop
1105     fTrackListJet->Clear();
1106   }//iJet loop
1107 }//FillJetProperties
1108 //_________________________________________________________________________________//
1109 void AliAnalysisTaskJetProperties::FillFFCorr(TList *jetlist){
1110   //filling up the histograms jets and tracks inside jet
1111   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillFFCorr() \n");
1112   
1113   for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
1114     Float_t JetPt;Float_t JetTheta;
1115     AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
1116     if(!jet)continue;
1117     JetTheta = jet->Theta();
1118     JetPt  = jet->Pt();
1119     fh1CorrJetPt -> Fill(JetPt);
1120     
1121     fTrackListJet->Clear();
1122     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1123     Int_t nJetTracks = 0;
1124     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1125     if(fDebug>2){
1126       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1127       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1128     }//fDebug
1129     
1130     Float_t TotPt = 0.;
1131     for (Int_t j =0; j< fTrackListJet->GetEntries(); j++){
1132       if(fTrackType==kTrackUndef)continue;
1133       Float_t TrackPt1=-99.0;  Float_t TrackPt2=-99.0; 
1134       Float_t FF1=-99.0;       Float_t FF2=-99.0;
1135       Float_t Ksi1=-99.0;      Float_t Ksi2=-99.0;
1136       Float_t TrackTheta1=0.0; Float_t TrackTheta2=0.0; 
1137       Float_t DelTheta1=0.0;   Float_t DelTheta2=0.0; 
1138       Float_t jT1=-99.0;       Float_t jT2=-99.0;
1139       if(fTrackType==kTrackAOD){
1140         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1141         if(!trackaod)continue;
1142         TrackTheta1 = trackaod->Theta();
1143         TrackPt1    = trackaod->Pt();
1144       }//if kTrackAOD
1145       else if(fTrackType==kTrackAODMC){
1146         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1147         if(!trackmc)continue;
1148         TrackTheta1 = trackmc->Theta();
1149         TrackPt1    = trackmc->Pt();
1150       }//if kTrackAODMC
1151       else if(fTrackType==kTrackKine){
1152         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1153         if(!trackkine)continue;
1154         TrackTheta1 = trackkine->Theta();
1155         TrackPt1    = trackkine->Pt();
1156       }//kTrackKine
1157       TotPt += TrackPt1;
1158       if(JetPt)FF1 = TrackPt1/JetPt;
1159       if(FF1)Ksi1  = TMath::Log(1./FF1);
1160       DelTheta1    = TMath::Abs(JetTheta - TrackTheta1);
1161       jT1          = TrackPt1 * TMath::Sin(DelTheta1);
1162       fh2CorrPtTrack1  ->Fill(JetPt,TrackPt1);
1163       fh2CorrFF1       ->Fill(JetPt,FF1);
1164       fh2CorrKsi1      ->Fill(JetPt,Ksi1);
1165       fh2CorrjT1       ->Fill(JetPt,jT1);
1166       for (Int_t jj =j+1; jj< fTrackListJet->GetEntries(); jj++){
1167         if(fTrackType==kTrackUndef)continue;
1168         if(fTrackType==kTrackAOD){
1169           AliAODTrack *trackaod2 = dynamic_cast<AliAODTrack*>(fTrackListJet->At(jj));
1170           if(!trackaod2)continue;
1171           TrackTheta2 = trackaod2->Theta();
1172           TrackPt2    = trackaod2->Pt();
1173         }//if kTrackAOD
1174         else if(fTrackType==kTrackAODMC){
1175           AliAODMCParticle* trackmc2 = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(jj));
1176           if(!trackmc2)continue;
1177           TrackTheta2 = trackmc2->Theta();
1178           TrackPt2    = trackmc2->Pt();
1179         }//if kTrackAODMC
1180         else if(fTrackType==kTrackKine){
1181           AliVParticle* trackkine2 = dynamic_cast<AliVParticle*>(fTrackListJet->At(jj));
1182           if(!trackkine2)continue;
1183           TrackTheta2 = trackkine2->Theta();
1184           TrackPt2    = trackkine2->Pt();
1185         }//kTrackKine
1186         if(JetPt)FF2 = TrackPt2/JetPt;
1187         if(FF2)Ksi2  = TMath::Log(1./FF2);
1188         DelTheta2    = TMath::Abs(JetTheta - TrackTheta2);
1189         jT2          = TrackPt2 * TMath::Sin(DelTheta2);
1190         Float_t ptmin=10.; Float_t ptmax=20.0;
1191         for(Int_t iBin=0; iBin<6; iBin++){
1192           if(iBin == 0)ptmin=10.; ptmax=20.;
1193           if(iBin == 1)ptmin=20.; ptmax=30.;
1194           if(iBin == 2)ptmin=30.; ptmax=40.;
1195           if(iBin == 3)ptmin=40.; ptmax=60.;
1196           if(iBin == 4)ptmin=60.; ptmax=80.;
1197           if(iBin == 5)ptmin=80.; ptmax=100.;
1198           if(JetPt>ptmin && JetPt <= ptmax){
1199             fh2CorrPt1Pt2[iBin]  ->Fill(TrackPt1,TrackPt2);
1200             fh2CorrZ1Z2[iBin]    ->Fill(FF1,FF2);
1201             fh2CorrKsi1Ksi2[iBin]->Fill(Ksi1,Ksi2);
1202             fh2CorrjT1jT2[iBin]  ->Fill(jT1,jT2);
1203           }//if loop
1204         }//iBin loop
1205       }//inside track loop
1206     }//outer track loop
1207     Float_t diff_JetPt_TrkPtSum = JetPt - TotPt; 
1208     fh1JetPtvsTrkSum->Fill(diff_JetPt_TrkPtSum);
1209     TotPt = 0.;
1210     fTrackListJet->Clear();
1211   }//iJet loop
1212 }//FillFFCorr
1213 //_________________________________________________________________________________//
1214
1215 void AliAnalysisTaskJetProperties::FillJetShape(TList *jetlist){
1216   //filling up the histograms
1217   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
1218   Float_t JetEta; Float_t JetPhi; Float_t JetPt;
1219   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
1220   if(jet){
1221     JetEta = jet->Eta();
1222     JetPhi = jet->Phi();
1223     JetPt  = jet->Pt();
1224     fh1PtLeadingJet->Fill(JetPt);
1225     Float_t NchSumA[50]         = {0.};
1226     Float_t PtSumA[50]          = {0.};
1227     Float_t delRPtSum80pc       = 0;
1228     Float_t delRNtrkSum80pc     = 0;
1229     Float_t PtSumDiffShape[10]  = {0.0};
1230     Float_t PtSumIntShape[10]   = {0.0};
1231     Int_t kNbinsR               = 10;
1232     fTrackListJet->Clear();
1233     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1234     Int_t nJetTracks = 0;
1235     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1236     if(fDebug>3){
1237       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1238       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1239     }
1240     fh2NtracksLeadingJet->Fill(JetPt,nJetTracks);
1241     fProNtracksLeadingJet->Fill(JetPt,nJetTracks);
1242     Int_t   *index     = new Int_t   [nJetTracks];//dynamic array containing index
1243     Float_t *delRA     = new Float_t [nJetTracks];//dynamic array containing delR
1244     Float_t *delEtaA   = new Float_t [nJetTracks];//dynamic array containing delEta
1245     Float_t *delPhiA   = new Float_t [nJetTracks];//dynamic array containing delPhi
1246     Float_t *trackPtA  = new Float_t [nJetTracks];//dynamic array containing pt-track
1247     Float_t *trackEtaA = new Float_t [nJetTracks];//dynamic array containing eta-track
1248     Float_t *trackPhiA = new Float_t [nJetTracks];//dynamic array containing phi-track
1249     for(Int_t ii=0; ii<nJetTracks; ii++){
1250       index[ii]     = 0;
1251       delRA[ii]     = 0.;
1252       delEtaA[ii]   = 0.;
1253       delPhiA[ii]   = 0.;
1254       trackPtA[ii]  = 0.;
1255       trackEtaA[ii] = 0.;
1256       trackPhiA[ii] = 0.;
1257     }//ii loop
1258   
1259     for (Int_t j =0; j< nJetTracks; j++){
1260       if(fTrackType==kTrackUndef)continue;      
1261       
1262       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1263       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1264       
1265       if(fTrackType==kTrackAOD){
1266         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1267         if(!trackaod)continue;
1268         TrackEta = trackaod->Eta();
1269         TrackPhi = trackaod->Phi();
1270         TrackPt  = trackaod->Pt();
1271       }//if kTrackAOD
1272       else if(fTrackType==kTrackAODMC){
1273         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1274         if(!trackmc)continue;
1275         TrackEta = trackmc->Eta();
1276         TrackPhi = trackmc->Phi();
1277         TrackPt  = trackmc->Pt();
1278       }//if kTrackAODMC
1279       else if(fTrackType==kTrackKine){
1280         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1281         if(!trackkine)continue;
1282         TrackEta = trackkine->Eta();
1283         TrackPhi = trackkine->Phi();
1284         TrackPt  = trackkine->Pt();
1285       }//if kTrackKine
1286       
1287       DelEta = TMath::Abs(JetEta - TrackEta);
1288       DelPhi = TMath::Abs(JetPhi - TrackPhi);
1289       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1290       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1291       AreaJ  = TMath::Pi()*DelR*DelR;
1292       
1293       fh2AreaCh ->Fill(JetPt,AreaJ);
1294       fProAreaCh->Fill(JetPt,AreaJ);
1295       delRA[j]     = DelR;
1296       delEtaA[j]   = DelEta;
1297       delPhiA[j]   = DelPhi;
1298       trackPtA[j]  = TrackPt;
1299       trackEtaA[j] = TrackEta;
1300       trackPhiA[j] = TrackPhi;
1301       
1302       //calculating diff and integrated jet shapes
1303       Float_t kDeltaR = 0.1;
1304       Float_t RMin    = kDeltaR/2.0;
1305       Float_t RMax    = kDeltaR/2.0;
1306       Float_t tmpR    = 0.05;
1307       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1308         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1309         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1310         tmpR += 0.1;
1311       }//ii1 loop
1312       
1313       for(Int_t ibin=1; ibin<=50; ibin++){
1314         Float_t xlow = 0.02*(ibin-1);
1315         Float_t xup  = 0.02*ibin;
1316         if( xlow <= DelR && DelR < xup){
1317           NchSumA[ibin-1]++;
1318           PtSumA[ibin-1]+= TrackPt;
1319         }//if loop
1320       }//for ibin loop
1321     }//track loop
1322     fTrackListJet->Clear();
1323     
1324     //---------------------
1325     Float_t tmp1R = 0.05;
1326     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1327       if(JetPt>20 && JetPt<=100){
1328         fProDiffJetShape->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1329         fProIntJetShape ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1330       }
1331       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1332       for(Int_t k=0; k<13; k++){
1333         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1334         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1335         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1336         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1337         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1338         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1339         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1340         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1341         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1342         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1343         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1344         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1345         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1346         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1347           fProDiffJetShapeA[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1348           fProIntJetShapeA[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1349         }//if
1350       }//k loop
1351       tmp1R +=0.1;
1352     }//jj1 loop
1353     //----------------------//
1354     Float_t PtSum = 0;
1355     Int_t NtrkSum = 0;
1356     Bool_t iflagPtSum   = kFALSE;
1357     Bool_t iflagNtrkSum = kFALSE;
1358     TMath::Sort(nJetTracks,delRA,index,0);
1359     for(Int_t ii=0; ii<nJetTracks; ii++){
1360       NtrkSum ++;
1361       PtSum += trackPtA[index[ii]];
1362       /*
1363         cout << index[ii] << "\t" <<
1364         delR[ii] << "\t" <<
1365         delEta[ii]<< "\t" <<
1366         delPhi[ii]<< "\t" <<
1367         trackPt[ii]<< "\t" <<
1368         trackEta[ii]<< "\t" <<
1369         trackPhi[ii]<< "\t DelR " <<
1370         delR[index[ii]] << endl;
1371       */
1372       if(!iflagNtrkSum){
1373         if((Float_t)NtrkSum/(Float_t)nJetTracks > 0.79){
1374           delRNtrkSum80pc = delRA[index[ii]];
1375           iflagNtrkSum = kTRUE;
1376         }//if NtrkSum
1377       }//if iflag
1378       if(!iflagPtSum){
1379         if(PtSum/JetPt >= 0.8000){
1380           delRPtSum80pc = delRA[index[ii]];
1381           iflagPtSum = kTRUE;
1382         }//if PtSum
1383       }//if iflag
1384     }//track loop 2nd
1385     delete [] index;
1386     delete [] delRA;
1387     delete [] delEtaA;
1388     delete [] delPhiA;
1389     delete [] trackPtA;
1390     delete [] trackEtaA;
1391     delete [] trackPhiA;
1392     fh2DelR80pcNch ->Fill(JetPt,delRNtrkSum80pc);
1393     fProDelR80pcNch->Fill(JetPt,delRNtrkSum80pc);
1394     fh2DelR80pcPt  ->Fill(JetPt,delRPtSum80pc);
1395     fProDelR80pcPt ->Fill(JetPt,delRPtSum80pc);
1396     for(Int_t ibin=0; ibin<50; ibin++){
1397       Float_t iR = 0.02*ibin + 0.01;
1398       fh3PtDelRNchSum->Fill(JetPt,iR,NchSumA[ibin]);
1399       fh3PtDelRPtSum ->Fill(JetPt,iR,PtSumA[ibin]);
1400       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1401       for(Int_t k=0; k<13; k++){
1402         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1403         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1404         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1405         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1406         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1407         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1408         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1409         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1410         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1411         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1412         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1413         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1414         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1415         if(JetPt>jetPtMin && JetPt<jetPtMax){
1416           fProDelRNchSum[k]->Fill(iR,NchSumA[ibin]);
1417           fProDelRPtSum[k]->Fill(iR,PtSumA[ibin]);
1418         }//if
1419       }//k loop
1420     }//ibin loop
1421   }//if(jet)
1422 }//FillJetShape()
1423 //_________________________________________________________________________________//
1424
1425 void AliAnalysisTaskJetProperties::FillJetShapeUE(TList *jetlist){
1426   //filling up the histograms
1427   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
1428   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
1429   if(jet){
1430     Double_t jetMom[3];
1431     jet->PxPyPz(jetMom);
1432     
1433     TVector3 jet3mom(jetMom);
1434     // Rotate phi and keep eta unchanged
1435     
1436     const Double_t alpha = TMath::Pi()/2.;
1437     Double_t etaTilted = jet3mom.Eta();
1438     Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
1439     if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1440     
1441     Double_t JetPt  = jet->Pt();
1442     
1443     Float_t NchSumA[50]         = {0.};
1444     Float_t PtSumA[50]          = {0.};
1445     Float_t delRPtSum80pc       = 0;
1446     Float_t delRNtrkSum80pc     = 0;
1447     Float_t PtSumDiffShape[10]  = {0.0};
1448     Float_t PtSumIntShape[10]   = {0.0};
1449     Int_t kNbinsR               = 10;
1450     
1451     //Int_t nTracks = GetListOfTracks(fTrackList,fTrackType);
1452     fTrackList->Clear();
1453     GetListOfTracks(fTrackList,fTrackType);
1454     Double_t sumPtPerp = 0;
1455     //if(fDebug > 100) printf("Cone radius for bckg. = %f Track type = %d\n",fJetRadius, fTrackType);
1456     fTrackListUE->Clear();
1457     GetTracksTiltedwrpJetAxis(alpha, fTrackList, fTrackListUE,jet,fJetRadius,sumPtPerp);
1458     fTrackList->Clear();
1459     fh1PtSumInJetConeUE->Fill(sumPtPerp);
1460     Int_t nTracksUE = fTrackListUE->GetEntries();
1461     fh2NtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1462     fProNtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1463     
1464     Int_t   *index     = new Int_t   [nTracksUE];//dynamic array containing index
1465     Float_t *delRA     = new Float_t [nTracksUE];//dynamic array containing delR
1466     Float_t *delEtaA   = new Float_t [nTracksUE];//dynamic array containing delEta
1467     Float_t *delPhiA   = new Float_t [nTracksUE];//dynamic array containing delPhi
1468     Float_t *trackPtA  = new Float_t [nTracksUE];//dynamic array containing pt-track
1469     Float_t *trackEtaA = new Float_t [nTracksUE];//dynamic array containing eta-track
1470     Float_t *trackPhiA = new Float_t [nTracksUE];//dynamic array containing phi-track
1471     for(Int_t ii=0; ii<nTracksUE; ii++){
1472       index[ii]     = 0;
1473       delRA[ii]     = 0.;
1474       delEtaA[ii]   = 0.;
1475       delPhiA[ii]   = 0.;
1476       trackPtA[ii]  = 0.;
1477       trackEtaA[ii] = 0.;
1478       trackPhiA[ii] = 0.;
1479     }//ii loop
1480     
1481     for (Int_t j =0; j< nTracksUE; j++){
1482       if(fTrackType==kTrackUndef)continue;      
1483       
1484       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1485       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1486       
1487       if(fTrackType==kTrackAOD){
1488         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListUE->At(j));
1489         if(!trackaod)continue;
1490         TrackEta = trackaod->Eta();
1491         TrackPhi = trackaod->Phi();
1492         TrackPt  = trackaod->Pt();
1493         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1494       }//if kTrackAOD
1495       else if(fTrackType==kTrackAODMC){
1496         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListUE->At(j));
1497         if(!trackmc)continue;
1498         TrackEta = trackmc->Eta();
1499         TrackPhi = trackmc->Phi();
1500         TrackPt  = trackmc->Pt();
1501         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1502       }//if kTrackAODMC
1503       else if(fTrackType==kTrackKine){
1504         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListUE->At(j));
1505         if(!trackkine)continue;
1506         TrackEta = trackkine->Eta();
1507         TrackPhi = trackkine->Phi();
1508         TrackPt  = trackkine->Pt();
1509         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1510       }//if kTrackKine
1511       
1512       DelEta = TMath::Abs(etaTilted - TrackEta);
1513       DelPhi = TMath::Abs(phiTilted - TrackPhi);
1514       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1515       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1516       AreaJ  = TMath::Pi()*DelR*DelR;
1517       
1518       fh2AreaChUE ->Fill(JetPt,AreaJ);
1519       fProAreaChUE->Fill(JetPt,AreaJ);
1520       
1521       delRA[j]     = DelR;
1522       delEtaA[j]   = DelEta;
1523       delPhiA[j]   = DelPhi;
1524       trackPtA[j]  = TrackPt;
1525       trackEtaA[j] = TrackEta;
1526       trackPhiA[j] = TrackPhi;
1527       
1528       //calculating diff and integrated jet shapes
1529       Float_t kDeltaR = 0.1;
1530       Float_t RMin    = kDeltaR/2.0;
1531       Float_t RMax    = kDeltaR/2.0;
1532       Float_t tmpR    = 0.05;
1533       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1534         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1535         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1536         tmpR += 0.1;
1537       }//ii1 loop
1538       
1539       for(Int_t ibin=1; ibin<=50; ibin++){
1540         Float_t xlow = 0.02*(ibin-1);
1541         Float_t xup  = 0.02*ibin;
1542         if( xlow <= DelR && DelR < xup){
1543           NchSumA[ibin-1]++;
1544           PtSumA[ibin-1]+= TrackPt;
1545         }//if loop
1546       }//for ibin loop
1547     }//track loop
1548     fTrackListUE->Clear();
1549     
1550     //---------------------
1551     Float_t tmp1R = 0.05;
1552     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1553       if(JetPt>20 && JetPt<=100){
1554         fProDiffJetShapeUE->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1555         fProIntJetShapeUE ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1556       }
1557       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1558       for(Int_t k=0; k<13; k++){
1559         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1560         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1561         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1562         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1563         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1564         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1565         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1566         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1567         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1568         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1569         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1570         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1571         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1572         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1573           fProDiffJetShapeAUE[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1574           fProIntJetShapeAUE[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1575         }//if
1576       }//k loop
1577       tmp1R +=0.1;
1578     }//jj1 loop
1579     //----------------------//
1580     Float_t PtSum = 0;
1581     Int_t NtrkSum = 0;
1582     Bool_t iflagPtSum   = kFALSE;
1583     Bool_t iflagNtrkSum = kFALSE;
1584     TMath::Sort(nTracksUE,delRA,index,0);
1585     for(Int_t ii=0; ii<nTracksUE; ii++){
1586       NtrkSum ++;
1587       PtSum += trackPtA[index[ii]];
1588       /*
1589         cout << index[ii] << "\t" <<
1590         delR[ii] << "\t" <<
1591         delEta[ii]<< "\t" <<
1592         delPhi[ii]<< "\t" <<
1593         trackPt[ii]<< "\t" <<
1594         trackEta[ii]<< "\t" <<
1595         trackPhi[ii]<< "\t DelR " <<
1596         delR[index[ii]] << endl;
1597       */
1598       if(!iflagNtrkSum){
1599         if((Float_t)NtrkSum/(Float_t)nTracksUE > 0.79){
1600           delRNtrkSum80pc = delRA[index[ii]];
1601           iflagNtrkSum = kTRUE;
1602         }//if NtrkSum
1603       }//if iflag
1604       if(!iflagPtSum){
1605         if(PtSum/JetPt >= 0.8000){
1606           delRPtSum80pc = delRA[index[ii]];
1607           iflagPtSum = kTRUE;
1608         }//if PtSum
1609       }//if iflag
1610     }//track loop 2nd
1611     delete [] index;
1612     delete [] delRA;
1613     delete [] delEtaA;
1614     delete [] delPhiA;
1615     delete [] trackPtA;
1616     delete [] trackEtaA;
1617     delete [] trackPhiA;
1618     fh2DelR80pcNchUE ->Fill(JetPt,delRNtrkSum80pc);
1619     fProDelR80pcNchUE->Fill(JetPt,delRNtrkSum80pc);
1620     fh2DelR80pcPtUE  ->Fill(JetPt,delRPtSum80pc);
1621     fProDelR80pcPtUE ->Fill(JetPt,delRPtSum80pc);
1622     for(Int_t ibin=0; ibin<50; ibin++){
1623       Float_t iR = 0.02*ibin + 0.01;
1624       fh3PtDelRNchSumUE->Fill(JetPt,iR,NchSumA[ibin]);
1625       fh3PtDelRPtSumUE ->Fill(JetPt,iR,PtSumA[ibin]);
1626       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1627       for(Int_t k=0; k<13; k++){
1628         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1629         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1630         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1631         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1632         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1633         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1634         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1635         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1636         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1637         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1638         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1639         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1640         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1641         if(JetPt>jetPtMin && JetPt<jetPtMax){
1642           fProDelRNchSumUE[k]->Fill(iR,NchSumA[ibin]);
1643           fProDelRPtSumUE[k]->Fill(iR,PtSumA[ibin]);
1644         }//if
1645       }//k loop
1646     }//ibin loop
1647   }//if(jet)
1648 }//FillJetShapeUE()
1649 //_________________________________________________________________________________//
1650
1651 void AliAnalysisTaskJetProperties::GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt){
1652   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx
1653   // List of tracks in cone perpendicular to the jet azimuthal direction
1654   Double_t jetMom[3];
1655   jet->PxPyPz(jetMom);
1656
1657   TVector3 jet3mom(jetMom);
1658   // Rotate phi and keep eta unchanged
1659   Double_t etaTilted = jet3mom.Eta();//no change in eta
1660   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;//rotate phi by alpha
1661   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1662   
1663   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
1664     if(fTrackType==kTrackUndef)continue;    
1665     Double_t trackMom[3];
1666
1667     if(fTrackType==kTrackAOD){
1668       AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(inputlist->At(itrack));
1669       if(!trackaod)continue;
1670       trackaod->PxPyPz(trackMom);
1671       TVector3 track3mom(trackMom);
1672       
1673       Double_t deta = track3mom.Eta() - etaTilted;
1674       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1675       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1676       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1677       
1678       if(dR<=radius){
1679         outputlist->Add(trackaod);
1680         sumPt += trackaod->Pt();
1681         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1682       }//if dR<=radius
1683     }//if kTrackAOD
1684     else if(fTrackType==kTrackAODMC){
1685       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(inputlist->At(itrack));
1686       if(!trackmc)continue;
1687       trackmc->PxPyPz(trackMom);
1688       TVector3 track3mom(trackMom);
1689       
1690       Double_t deta = track3mom.Eta() - etaTilted;
1691       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1692       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1693       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1694       
1695       if(dR<=radius){
1696         outputlist->Add(trackmc);
1697         sumPt += trackmc->Pt();
1698         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1699       }//if dR<=radius
1700     }//if kTrackAODMC
1701     else if(fTrackType==kTrackKine){
1702       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
1703       if(!trackkine)continue;
1704       trackkine->PxPyPz(trackMom);
1705       TVector3 track3mom(trackMom);
1706       
1707       Double_t deta = track3mom.Eta() - etaTilted;
1708       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1709       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1710       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1711       
1712       if(dR<=radius){
1713         outputlist->Add(trackkine);
1714         sumPt += trackkine->Pt();
1715         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1716       }//if dR<=radius
1717     }//kTrackKine
1718   }//itrack
1719 }//initial loop
1720 //-----------------------------------------------//
1721 Int_t AliAnalysisTaskJetProperties::GetListOfTracks(TList *list, Int_t type)
1722 {
1723   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx (and modified)
1724   // fill list of tracks selected according to type
1725   
1726   if(fDebug > 3) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
1727   
1728   if(!list){
1729     if(fDebug>3) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1730     return -1;
1731   }
1732   
1733   if(!fAOD) return -1;
1734   
1735   if(!fAOD->GetTracks()) return 0;
1736   
1737   if(type==kTrackUndef) return 0;
1738   
1739   Int_t iCount = 0;
1740   if(type==kTrackAOD){
1741     // all rec. tracks, esd filter mask, within acceptance
1742     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
1743       AliAODTrack *tr = fAOD->GetTrack(it);
1744       if(!tr)continue;
1745       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;//selecting filtermask
1746       if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
1747       if(tr->Pt()  < fTrackPtCut) continue;
1748       list->Add(tr);
1749       iCount++;
1750     }//track loop
1751   }//if type
1752   else if (type==kTrackKine){
1753     // kine particles, primaries, charged only within acceptance
1754     if(!fMCEvent) return iCount;
1755     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
1756       if(!fMCEvent->IsPhysicalPrimary(it))continue;//selecting only primaries
1757       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
1758       if(!part)continue;
1759       if(part->Particle()->GetPDG()->Charge()==0)continue;//removing neutrals
1760       if(part->Eta() < fTrackEtaMin || part->Eta() > fTrackEtaMax) continue;//eta cut
1761       if(part->Pt()  < fTrackPtCut)continue;//pt cut
1762       list->Add(part);
1763       iCount++;
1764     }//track loop
1765   }//if type
1766   else if (type==kTrackAODMC) {
1767     // MC particles (from AOD), physical primaries, charged only within acceptance
1768     if(!fAOD) return -1;
1769     
1770     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1771     if(!tca)return iCount;
1772     
1773     for(int it=0; it<tca->GetEntriesFast(); ++it){
1774       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1775       if(!part)continue;
1776       if(!part->IsPhysicalPrimary())continue;
1777       if(part->Charge()==0) continue;
1778       if(part->Eta() > fTrackEtaMax || part->Eta() < fTrackEtaMin)continue;
1779       if(part->Pt()  < fTrackPtCut) continue;
1780       list->Add(part);
1781       iCount++;
1782     }//track loop
1783   }//if type
1784   
1785   list->Sort();
1786   return iCount;
1787 }
1788 //_______________________________________________________________________________