Adapt add macro and particle/cluster containers for the analysis on jets
[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;   
1068         Float_t Ksi=-99.0;
1069         if(fTrackType==kTrackAOD){
1070           AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1071           if(!trackaod)continue;
1072           TrackEta = trackaod->Eta();
1073           TrackPhi = trackaod->Phi();
1074           TrackPt  = trackaod->Pt();
1075         }//if kTrackAOD
1076         else if(fTrackType==kTrackAODMC){
1077           AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1078           if(!trackmc)continue;
1079           TrackEta = trackmc->Eta();
1080           TrackPhi = trackmc->Phi();
1081           TrackPt  = trackmc->Pt();
1082         }//if kTrackAODMC
1083         else if(fTrackType==kTrackKine){
1084           AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1085           if(!trackkine)continue;
1086           TrackEta = trackkine->Eta();
1087           TrackPhi = trackkine->Phi();
1088           TrackPt  = trackkine->Pt();
1089         }//kTrackKine
1090         if(JetPt)FF = TrackPt/JetPt;
1091         if(FF)Ksi = TMath::Log(1./FF);
1092         DelEta      = TMath::Abs(JetEta - TrackEta);
1093         DelPhi      = TMath::Abs(JetPhi - TrackPhi);
1094         if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1095         DelR        = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1096         // AreaJ       = TMath::Pi()*DelR*DelR;
1097         fh2EtaTrack ->Fill(JetPt,TrackEta);
1098         fh2PhiTrack ->Fill(JetPt,TrackPhi);
1099         fh2PtTrack  ->Fill(JetPt,TrackPt);
1100         fh2FF       ->Fill(JetPt,FF);
1101         fh2Ksi      ->Fill(JetPt,Ksi);
1102         fh2DelEta   ->Fill(JetPt,DelEta);
1103         fh2DelPhi   ->Fill(JetPt,DelPhi);
1104         fh2DelR     ->Fill(JetPt,DelR);
1105     }//track loop
1106     fTrackListJet->Clear();
1107   }//iJet loop
1108 }//FillJetProperties
1109 //_________________________________________________________________________________//
1110 void AliAnalysisTaskJetProperties::FillFFCorr(TList *jetlist){
1111   //filling up the histograms jets and tracks inside jet
1112   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillFFCorr() \n");
1113   
1114   for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
1115     Float_t JetPt;Float_t JetTheta;
1116     AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
1117     if(!jet)continue;
1118     JetTheta = jet->Theta();
1119     JetPt  = jet->Pt();
1120     fh1CorrJetPt -> Fill(JetPt);
1121     
1122     fTrackListJet->Clear();
1123     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1124     Int_t nJetTracks = 0;
1125     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1126     if(fDebug>2){
1127       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1128       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1129     }//fDebug
1130     
1131     Float_t TotPt = 0.;
1132     for (Int_t j =0; j< fTrackListJet->GetEntries(); j++){
1133       if(fTrackType==kTrackUndef)continue;
1134       Float_t TrackPt1=-99.0;  Float_t TrackPt2=-99.0; 
1135       Float_t FF1=-99.0;       Float_t FF2=-99.0;
1136       Float_t Ksi1=-99.0;      Float_t Ksi2=-99.0;
1137       Float_t TrackTheta1=0.0; Float_t TrackTheta2=0.0; 
1138       Float_t DelTheta1=0.0;   Float_t DelTheta2=0.0; 
1139       Float_t jT1=-99.0;       Float_t jT2=-99.0;
1140       if(fTrackType==kTrackAOD){
1141         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1142         if(!trackaod)continue;
1143         TrackTheta1 = trackaod->Theta();
1144         TrackPt1    = trackaod->Pt();
1145       }//if kTrackAOD
1146       else if(fTrackType==kTrackAODMC){
1147         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1148         if(!trackmc)continue;
1149         TrackTheta1 = trackmc->Theta();
1150         TrackPt1    = trackmc->Pt();
1151       }//if kTrackAODMC
1152       else if(fTrackType==kTrackKine){
1153         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1154         if(!trackkine)continue;
1155         TrackTheta1 = trackkine->Theta();
1156         TrackPt1    = trackkine->Pt();
1157       }//kTrackKine
1158       TotPt += TrackPt1;
1159       if(JetPt)FF1 = TrackPt1/JetPt;
1160       if(FF1)Ksi1  = TMath::Log(1./FF1);
1161       DelTheta1    = TMath::Abs(JetTheta - TrackTheta1);
1162       jT1          = TrackPt1 * TMath::Sin(DelTheta1);
1163       fh2CorrPtTrack1  ->Fill(JetPt,TrackPt1);
1164       fh2CorrFF1       ->Fill(JetPt,FF1);
1165       fh2CorrKsi1      ->Fill(JetPt,Ksi1);
1166       fh2CorrjT1       ->Fill(JetPt,jT1);
1167       for (Int_t jj =j+1; jj< fTrackListJet->GetEntries(); jj++){
1168         if(fTrackType==kTrackUndef)continue;
1169         if(fTrackType==kTrackAOD){
1170           AliAODTrack *trackaod2 = dynamic_cast<AliAODTrack*>(fTrackListJet->At(jj));
1171           if(!trackaod2)continue;
1172           TrackTheta2 = trackaod2->Theta();
1173           TrackPt2    = trackaod2->Pt();
1174         }//if kTrackAOD
1175         else if(fTrackType==kTrackAODMC){
1176           AliAODMCParticle* trackmc2 = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(jj));
1177           if(!trackmc2)continue;
1178           TrackTheta2 = trackmc2->Theta();
1179           TrackPt2    = trackmc2->Pt();
1180         }//if kTrackAODMC
1181         else if(fTrackType==kTrackKine){
1182           AliVParticle* trackkine2 = dynamic_cast<AliVParticle*>(fTrackListJet->At(jj));
1183           if(!trackkine2)continue;
1184           TrackTheta2 = trackkine2->Theta();
1185           TrackPt2    = trackkine2->Pt();
1186         }//kTrackKine
1187         if(JetPt)FF2 = TrackPt2/JetPt;
1188         if(FF2)Ksi2  = TMath::Log(1./FF2);
1189         DelTheta2    = TMath::Abs(JetTheta - TrackTheta2);
1190         jT2          = TrackPt2 * TMath::Sin(DelTheta2);
1191         Float_t ptmin=10.; Float_t ptmax=20.0;
1192         for(Int_t iBin=0; iBin<6; iBin++){
1193           if(iBin == 0)ptmin=10.; ptmax=20.;
1194           if(iBin == 1)ptmin=20.; ptmax=30.;
1195           if(iBin == 2)ptmin=30.; ptmax=40.;
1196           if(iBin == 3)ptmin=40.; ptmax=60.;
1197           if(iBin == 4)ptmin=60.; ptmax=80.;
1198           if(iBin == 5)ptmin=80.; ptmax=100.;
1199           if(JetPt>ptmin && JetPt <= ptmax){
1200             fh2CorrPt1Pt2[iBin]  ->Fill(TrackPt1,TrackPt2);
1201             fh2CorrZ1Z2[iBin]    ->Fill(FF1,FF2);
1202             fh2CorrKsi1Ksi2[iBin]->Fill(Ksi1,Ksi2);
1203             fh2CorrjT1jT2[iBin]  ->Fill(jT1,jT2);
1204           }//if loop
1205         }//iBin loop
1206       }//inside track loop
1207     }//outer track loop
1208     Float_t diff_JetPt_TrkPtSum = JetPt - TotPt; 
1209     fh1JetPtvsTrkSum->Fill(diff_JetPt_TrkPtSum);
1210     TotPt = 0.;
1211     fTrackListJet->Clear();
1212   }//iJet loop
1213 }//FillFFCorr
1214 //_________________________________________________________________________________//
1215
1216 void AliAnalysisTaskJetProperties::FillJetShape(TList *jetlist){
1217   //filling up the histograms
1218   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
1219   Float_t JetEta; Float_t JetPhi; Float_t JetPt;
1220   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
1221   if(jet){
1222     JetEta = jet->Eta();
1223     JetPhi = jet->Phi();
1224     JetPt  = jet->Pt();
1225     fh1PtLeadingJet->Fill(JetPt);
1226     Float_t NchSumA[50]         = {0.};
1227     Float_t PtSumA[50]          = {0.};
1228     Float_t delRPtSum80pc       = 0;
1229     Float_t delRNtrkSum80pc     = 0;
1230     Float_t PtSumDiffShape[10]  = {0.0};
1231     Float_t PtSumIntShape[10]   = {0.0};
1232     Int_t kNbinsR               = 10;
1233     fTrackListJet->Clear();
1234     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1235     Int_t nJetTracks = 0;
1236     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1237     if(fDebug>3){
1238       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1239       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1240     }
1241     fh2NtracksLeadingJet->Fill(JetPt,nJetTracks);
1242     fProNtracksLeadingJet->Fill(JetPt,nJetTracks);
1243     Int_t   *index     = new Int_t   [nJetTracks];//dynamic array containing index
1244     Float_t *delRA     = new Float_t [nJetTracks];//dynamic array containing delR
1245     Float_t *delEtaA   = new Float_t [nJetTracks];//dynamic array containing delEta
1246     Float_t *delPhiA   = new Float_t [nJetTracks];//dynamic array containing delPhi
1247     Float_t *trackPtA  = new Float_t [nJetTracks];//dynamic array containing pt-track
1248     Float_t *trackEtaA = new Float_t [nJetTracks];//dynamic array containing eta-track
1249     Float_t *trackPhiA = new Float_t [nJetTracks];//dynamic array containing phi-track
1250     for(Int_t ii=0; ii<nJetTracks; ii++){
1251       index[ii]     = 0;
1252       delRA[ii]     = 0.;
1253       delEtaA[ii]   = 0.;
1254       delPhiA[ii]   = 0.;
1255       trackPtA[ii]  = 0.;
1256       trackEtaA[ii] = 0.;
1257       trackPhiA[ii] = 0.;
1258     }//ii loop
1259   
1260     for (Int_t j =0; j< nJetTracks; j++){
1261       if(fTrackType==kTrackUndef)continue;      
1262       
1263       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1264       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1265       
1266       if(fTrackType==kTrackAOD){
1267         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1268         if(!trackaod)continue;
1269         TrackEta = trackaod->Eta();
1270         TrackPhi = trackaod->Phi();
1271         TrackPt  = trackaod->Pt();
1272       }//if kTrackAOD
1273       else if(fTrackType==kTrackAODMC){
1274         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1275         if(!trackmc)continue;
1276         TrackEta = trackmc->Eta();
1277         TrackPhi = trackmc->Phi();
1278         TrackPt  = trackmc->Pt();
1279       }//if kTrackAODMC
1280       else if(fTrackType==kTrackKine){
1281         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1282         if(!trackkine)continue;
1283         TrackEta = trackkine->Eta();
1284         TrackPhi = trackkine->Phi();
1285         TrackPt  = trackkine->Pt();
1286       }//if kTrackKine
1287       
1288       DelEta = TMath::Abs(JetEta - TrackEta);
1289       DelPhi = TMath::Abs(JetPhi - TrackPhi);
1290       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1291       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1292       AreaJ  = TMath::Pi()*DelR*DelR;
1293       
1294       fh2AreaCh ->Fill(JetPt,AreaJ);
1295       fProAreaCh->Fill(JetPt,AreaJ);
1296       delRA[j]     = DelR;
1297       delEtaA[j]   = DelEta;
1298       delPhiA[j]   = DelPhi;
1299       trackPtA[j]  = TrackPt;
1300       trackEtaA[j] = TrackEta;
1301       trackPhiA[j] = TrackPhi;
1302       
1303       //calculating diff and integrated jet shapes
1304       Float_t kDeltaR = 0.1;
1305       Float_t RMin    = kDeltaR/2.0;
1306       Float_t RMax    = kDeltaR/2.0;
1307       Float_t tmpR    = 0.05;
1308       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1309         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1310         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1311         tmpR += 0.1;
1312       }//ii1 loop
1313       
1314       for(Int_t ibin=1; ibin<=50; ibin++){
1315         Float_t xlow = 0.02*(ibin-1);
1316         Float_t xup  = 0.02*ibin;
1317         if( xlow <= DelR && DelR < xup){
1318           NchSumA[ibin-1]++;
1319           PtSumA[ibin-1]+= TrackPt;
1320         }//if loop
1321       }//for ibin loop
1322     }//track loop
1323     fTrackListJet->Clear();
1324     
1325     //---------------------
1326     Float_t tmp1R = 0.05;
1327     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1328       if(JetPt>20 && JetPt<=100){
1329         fProDiffJetShape->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1330         fProIntJetShape ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1331       }
1332       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1333       for(Int_t k=0; k<13; k++){
1334         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1335         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1336         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1337         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1338         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1339         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1340         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1341         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1342         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1343         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1344         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1345         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1346         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1347         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1348           fProDiffJetShapeA[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1349           fProIntJetShapeA[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1350         }//if
1351       }//k loop
1352       tmp1R +=0.1;
1353     }//jj1 loop
1354     //----------------------//
1355     Float_t PtSum = 0;
1356     Int_t NtrkSum = 0;
1357     Bool_t iflagPtSum   = kFALSE;
1358     Bool_t iflagNtrkSum = kFALSE;
1359     TMath::Sort(nJetTracks,delRA,index,0);
1360     for(Int_t ii=0; ii<nJetTracks; ii++){
1361       NtrkSum ++;
1362       PtSum += trackPtA[index[ii]];
1363       /*
1364         cout << index[ii] << "\t" <<
1365         delR[ii] << "\t" <<
1366         delEta[ii]<< "\t" <<
1367         delPhi[ii]<< "\t" <<
1368         trackPt[ii]<< "\t" <<
1369         trackEta[ii]<< "\t" <<
1370         trackPhi[ii]<< "\t DelR " <<
1371         delR[index[ii]] << endl;
1372       */
1373       if(!iflagNtrkSum){
1374         if((Float_t)NtrkSum/(Float_t)nJetTracks > 0.79){
1375           delRNtrkSum80pc = delRA[index[ii]];
1376           iflagNtrkSum = kTRUE;
1377         }//if NtrkSum
1378       }//if iflag
1379       if(!iflagPtSum){
1380         if(PtSum/JetPt >= 0.8000){
1381           delRPtSum80pc = delRA[index[ii]];
1382           iflagPtSum = kTRUE;
1383         }//if PtSum
1384       }//if iflag
1385     }//track loop 2nd
1386     delete [] index;
1387     delete [] delRA;
1388     delete [] delEtaA;
1389     delete [] delPhiA;
1390     delete [] trackPtA;
1391     delete [] trackEtaA;
1392     delete [] trackPhiA;
1393     fh2DelR80pcNch ->Fill(JetPt,delRNtrkSum80pc);
1394     fProDelR80pcNch->Fill(JetPt,delRNtrkSum80pc);
1395     fh2DelR80pcPt  ->Fill(JetPt,delRPtSum80pc);
1396     fProDelR80pcPt ->Fill(JetPt,delRPtSum80pc);
1397     for(Int_t ibin=0; ibin<50; ibin++){
1398       Float_t iR = 0.02*ibin + 0.01;
1399       fh3PtDelRNchSum->Fill(JetPt,iR,NchSumA[ibin]);
1400       fh3PtDelRPtSum ->Fill(JetPt,iR,PtSumA[ibin]);
1401       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1402       for(Int_t k=0; k<13; k++){
1403         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1404         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1405         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1406         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1407         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1408         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1409         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1410         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1411         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1412         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1413         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1414         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1415         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1416         if(JetPt>jetPtMin && JetPt<jetPtMax){
1417           fProDelRNchSum[k]->Fill(iR,NchSumA[ibin]);
1418           fProDelRPtSum[k]->Fill(iR,PtSumA[ibin]);
1419         }//if
1420       }//k loop
1421     }//ibin loop
1422   }//if(jet)
1423 }//FillJetShape()
1424 //_________________________________________________________________________________//
1425
1426 void AliAnalysisTaskJetProperties::FillJetShapeUE(TList *jetlist){
1427   //filling up the histograms
1428   if(fDebug > 2) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
1429   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
1430   if(jet){
1431     Double_t jetMom[3];
1432     jet->PxPyPz(jetMom);
1433     
1434     TVector3 jet3mom(jetMom);
1435     // Rotate phi and keep eta unchanged
1436     
1437     const Double_t alpha = TMath::Pi()/2.;
1438     Double_t etaTilted = jet3mom.Eta();
1439     Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
1440     if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1441     
1442     Double_t JetPt  = jet->Pt();
1443     
1444     Float_t NchSumA[50]         = {0.};
1445     Float_t PtSumA[50]          = {0.};
1446     Float_t delRPtSum80pc       = 0;
1447     Float_t delRNtrkSum80pc     = 0;
1448     Float_t PtSumDiffShape[10]  = {0.0};
1449     Float_t PtSumIntShape[10]   = {0.0};
1450     Int_t kNbinsR               = 10;
1451     
1452     //Int_t nTracks = GetListOfTracks(fTrackList,fTrackType);
1453     fTrackList->Clear();
1454     GetListOfTracks(fTrackList,fTrackType);
1455     Double_t sumPtPerp = 0;
1456     //if(fDebug > 100) printf("Cone radius for bckg. = %f Track type = %d\n",fJetRadius, fTrackType);
1457     fTrackListUE->Clear();
1458     GetTracksTiltedwrpJetAxis(alpha, fTrackList, fTrackListUE,jet,fJetRadius,sumPtPerp);
1459     fTrackList->Clear();
1460     fh1PtSumInJetConeUE->Fill(sumPtPerp);
1461     Int_t nTracksUE = fTrackListUE->GetEntries();
1462     fh2NtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1463     fProNtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1464     
1465     Int_t   *index     = new Int_t   [nTracksUE];//dynamic array containing index
1466     Float_t *delRA     = new Float_t [nTracksUE];//dynamic array containing delR
1467     Float_t *delEtaA   = new Float_t [nTracksUE];//dynamic array containing delEta
1468     Float_t *delPhiA   = new Float_t [nTracksUE];//dynamic array containing delPhi
1469     Float_t *trackPtA  = new Float_t [nTracksUE];//dynamic array containing pt-track
1470     Float_t *trackEtaA = new Float_t [nTracksUE];//dynamic array containing eta-track
1471     Float_t *trackPhiA = new Float_t [nTracksUE];//dynamic array containing phi-track
1472     for(Int_t ii=0; ii<nTracksUE; ii++){
1473       index[ii]     = 0;
1474       delRA[ii]     = 0.;
1475       delEtaA[ii]   = 0.;
1476       delPhiA[ii]   = 0.;
1477       trackPtA[ii]  = 0.;
1478       trackEtaA[ii] = 0.;
1479       trackPhiA[ii] = 0.;
1480     }//ii loop
1481     
1482     for (Int_t j =0; j< nTracksUE; j++){
1483       if(fTrackType==kTrackUndef)continue;      
1484       
1485       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1486       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1487       
1488       if(fTrackType==kTrackAOD){
1489         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListUE->At(j));
1490         if(!trackaod)continue;
1491         TrackEta = trackaod->Eta();
1492         TrackPhi = trackaod->Phi();
1493         TrackPt  = trackaod->Pt();
1494         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1495       }//if kTrackAOD
1496       else if(fTrackType==kTrackAODMC){
1497         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListUE->At(j));
1498         if(!trackmc)continue;
1499         TrackEta = trackmc->Eta();
1500         TrackPhi = trackmc->Phi();
1501         TrackPt  = trackmc->Pt();
1502         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1503       }//if kTrackAODMC
1504       else if(fTrackType==kTrackKine){
1505         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListUE->At(j));
1506         if(!trackkine)continue;
1507         TrackEta = trackkine->Eta();
1508         TrackPhi = trackkine->Phi();
1509         TrackPt  = trackkine->Pt();
1510         //if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1511       }//if kTrackKine
1512       
1513       DelEta = TMath::Abs(etaTilted - TrackEta);
1514       DelPhi = TMath::Abs(phiTilted - TrackPhi);
1515       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1516       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1517       AreaJ  = TMath::Pi()*DelR*DelR;
1518       
1519       fh2AreaChUE ->Fill(JetPt,AreaJ);
1520       fProAreaChUE->Fill(JetPt,AreaJ);
1521       
1522       delRA[j]     = DelR;
1523       delEtaA[j]   = DelEta;
1524       delPhiA[j]   = DelPhi;
1525       trackPtA[j]  = TrackPt;
1526       trackEtaA[j] = TrackEta;
1527       trackPhiA[j] = TrackPhi;
1528       
1529       //calculating diff and integrated jet shapes
1530       Float_t kDeltaR = 0.1;
1531       Float_t RMin    = kDeltaR/2.0;
1532       Float_t RMax    = kDeltaR/2.0;
1533       Float_t tmpR    = 0.05;
1534       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1535         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1536         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1537         tmpR += 0.1;
1538       }//ii1 loop
1539       
1540       for(Int_t ibin=1; ibin<=50; ibin++){
1541         Float_t xlow = 0.02*(ibin-1);
1542         Float_t xup  = 0.02*ibin;
1543         if( xlow <= DelR && DelR < xup){
1544           NchSumA[ibin-1]++;
1545           PtSumA[ibin-1]+= TrackPt;
1546         }//if loop
1547       }//for ibin loop
1548     }//track loop
1549     fTrackListUE->Clear();
1550     
1551     //---------------------
1552     Float_t tmp1R = 0.05;
1553     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1554       if(JetPt>20 && JetPt<=100){
1555         fProDiffJetShapeUE->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1556         fProIntJetShapeUE ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1557       }
1558       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1559       for(Int_t k=0; k<13; k++){
1560         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1561         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1562         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1563         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1564         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1565         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1566         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1567         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1568         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1569         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1570         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1571         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1572         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1573         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1574           fProDiffJetShapeAUE[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1575           fProIntJetShapeAUE[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1576         }//if
1577       }//k loop
1578       tmp1R +=0.1;
1579     }//jj1 loop
1580     //----------------------//
1581     Float_t PtSum = 0;
1582     Int_t NtrkSum = 0;
1583     Bool_t iflagPtSum   = kFALSE;
1584     Bool_t iflagNtrkSum = kFALSE;
1585     TMath::Sort(nTracksUE,delRA,index,0);
1586     for(Int_t ii=0; ii<nTracksUE; ii++){
1587       NtrkSum ++;
1588       PtSum += trackPtA[index[ii]];
1589       /*
1590         cout << index[ii] << "\t" <<
1591         delR[ii] << "\t" <<
1592         delEta[ii]<< "\t" <<
1593         delPhi[ii]<< "\t" <<
1594         trackPt[ii]<< "\t" <<
1595         trackEta[ii]<< "\t" <<
1596         trackPhi[ii]<< "\t DelR " <<
1597         delR[index[ii]] << endl;
1598       */
1599       if(!iflagNtrkSum){
1600         if((Float_t)NtrkSum/(Float_t)nTracksUE > 0.79){
1601           delRNtrkSum80pc = delRA[index[ii]];
1602           iflagNtrkSum = kTRUE;
1603         }//if NtrkSum
1604       }//if iflag
1605       if(!iflagPtSum){
1606         if(PtSum/JetPt >= 0.8000){
1607           delRPtSum80pc = delRA[index[ii]];
1608           iflagPtSum = kTRUE;
1609         }//if PtSum
1610       }//if iflag
1611     }//track loop 2nd
1612     delete [] index;
1613     delete [] delRA;
1614     delete [] delEtaA;
1615     delete [] delPhiA;
1616     delete [] trackPtA;
1617     delete [] trackEtaA;
1618     delete [] trackPhiA;
1619     fh2DelR80pcNchUE ->Fill(JetPt,delRNtrkSum80pc);
1620     fProDelR80pcNchUE->Fill(JetPt,delRNtrkSum80pc);
1621     fh2DelR80pcPtUE  ->Fill(JetPt,delRPtSum80pc);
1622     fProDelR80pcPtUE ->Fill(JetPt,delRPtSum80pc);
1623     for(Int_t ibin=0; ibin<50; ibin++){
1624       Float_t iR = 0.02*ibin + 0.01;
1625       fh3PtDelRNchSumUE->Fill(JetPt,iR,NchSumA[ibin]);
1626       fh3PtDelRPtSumUE ->Fill(JetPt,iR,PtSumA[ibin]);
1627       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1628       for(Int_t k=0; k<13; k++){
1629         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1630         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1631         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1632         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1633         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1634         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1635         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1636         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1637         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1638         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1639         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1640         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1641         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1642         if(JetPt>jetPtMin && JetPt<jetPtMax){
1643           fProDelRNchSumUE[k]->Fill(iR,NchSumA[ibin]);
1644           fProDelRPtSumUE[k]->Fill(iR,PtSumA[ibin]);
1645         }//if
1646       }//k loop
1647     }//ibin loop
1648   }//if(jet)
1649 }//FillJetShapeUE()
1650 //_________________________________________________________________________________//
1651
1652 void AliAnalysisTaskJetProperties::GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt){
1653   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx
1654   // List of tracks in cone perpendicular to the jet azimuthal direction
1655   Double_t jetMom[3];
1656   jet->PxPyPz(jetMom);
1657
1658   TVector3 jet3mom(jetMom);
1659   // Rotate phi and keep eta unchanged
1660   Double_t etaTilted = jet3mom.Eta();//no change in eta
1661   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;//rotate phi by alpha
1662   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1663   
1664   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
1665     if(fTrackType==kTrackUndef)continue;    
1666     Double_t trackMom[3];
1667
1668     if(fTrackType==kTrackAOD){
1669       AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(inputlist->At(itrack));
1670       if(!trackaod)continue;
1671       trackaod->PxPyPz(trackMom);
1672       TVector3 track3mom(trackMom);
1673       
1674       Double_t deta = track3mom.Eta() - etaTilted;
1675       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1676       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1677       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1678       
1679       if(dR<=radius){
1680         outputlist->Add(trackaod);
1681         sumPt += trackaod->Pt();
1682         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1683       }//if dR<=radius
1684     }//if kTrackAOD
1685     else if(fTrackType==kTrackAODMC){
1686       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(inputlist->At(itrack));
1687       if(!trackmc)continue;
1688       trackmc->PxPyPz(trackMom);
1689       TVector3 track3mom(trackMom);
1690       
1691       Double_t deta = track3mom.Eta() - etaTilted;
1692       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1693       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1694       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1695       
1696       if(dR<=radius){
1697         outputlist->Add(trackmc);
1698         sumPt += trackmc->Pt();
1699         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1700       }//if dR<=radius
1701     }//if kTrackAODMC
1702     else if(fTrackType==kTrackKine){
1703       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
1704       if(!trackkine)continue;
1705       trackkine->PxPyPz(trackMom);
1706       TVector3 track3mom(trackMom);
1707       
1708       Double_t deta = track3mom.Eta() - etaTilted;
1709       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1710       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1711       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1712       
1713       if(dR<=radius){
1714         outputlist->Add(trackkine);
1715         sumPt += trackkine->Pt();
1716         //if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1717       }//if dR<=radius
1718     }//kTrackKine
1719   }//itrack
1720 }//initial loop
1721 //-----------------------------------------------//
1722 Int_t AliAnalysisTaskJetProperties::GetListOfTracks(TList *list, Int_t type)
1723 {
1724   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx (and modified)
1725   // fill list of tracks selected according to type
1726   
1727   if(fDebug > 3) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
1728   
1729   if(!list){
1730     if(fDebug>3) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1731     return -1;
1732   }
1733   
1734   if(!fAOD) return -1;
1735   
1736   if(!fAOD->GetTracks()) return 0;
1737   
1738   if(type==kTrackUndef) return 0;
1739   
1740   Int_t iCount = 0;
1741   if(type==kTrackAOD){
1742     // all rec. tracks, esd filter mask, within acceptance
1743     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
1744       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(it));
1745       if(!tr) AliFatal("Not a standard AOD");
1746       if(!tr)continue;
1747       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;//selecting filtermask
1748       if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
1749       if(tr->Pt()  < fTrackPtCut) continue;
1750       list->Add(tr);
1751       iCount++;
1752     }//track loop
1753   }//if type
1754   else if (type==kTrackKine){
1755     // kine particles, primaries, charged only within acceptance
1756     if(!fMCEvent) return iCount;
1757     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
1758       if(!fMCEvent->IsPhysicalPrimary(it))continue;//selecting only primaries
1759       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
1760       if(!part)continue;
1761       if(part->Particle()->GetPDG()->Charge()==0)continue;//removing neutrals
1762       if(part->Eta() < fTrackEtaMin || part->Eta() > fTrackEtaMax) continue;//eta cut
1763       if(part->Pt()  < fTrackPtCut)continue;//pt cut
1764       list->Add(part);
1765       iCount++;
1766     }//track loop
1767   }//if type
1768   else if (type==kTrackAODMC) {
1769     // MC particles (from AOD), physical primaries, charged only within acceptance
1770     if(!fAOD) return -1;
1771     
1772     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1773     if(!tca)return iCount;
1774     
1775     for(int it=0; it<tca->GetEntriesFast(); ++it){
1776       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1777       if(!part)continue;
1778       if(!part->IsPhysicalPrimary())continue;
1779       if(part->Charge()==0) continue;
1780       if(part->Eta() > fTrackEtaMax || part->Eta() < fTrackEtaMin)continue;
1781       if(part->Pt()  < fTrackPtCut) continue;
1782       list->Add(part);
1783       iCount++;
1784     }//track loop
1785   }//if type
1786   
1787   list->Sort();
1788   return iCount;
1789 }
1790 //_______________________________________________________________________________