setting of eta reange (Sidharth
[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    ,fUseAODInputJets(kTRUE)
70    ,fFilterMask(0)
71    ,fUsePhysicsSelection(kTRUE)
72    ,fMaxVertexZ(10)
73    ,fNContributors(2)
74    ,fTrackPtCut(0)
75    ,fTrackEtaMin(0)
76    ,fTrackEtaMax(0)
77    ,fJetPtCut(0)
78    ,fJetEtaMin(0)
79    ,fJetEtaMax(0)
80    ,fAvgTrials(0)
81    ,fJetList(0)
82    ,fTrackList(0)
83    ,fCommonHistList(0)
84    ,fh1EvtSelection(0)
85    ,fh1VertexNContributors(0)
86    ,fh1VertexZ(0)
87    ,fh1Xsec(0)
88    ,fh1Trials(0)
89    ,fh1PtHard(0)
90    ,fh1PtHardTrials(0)
91    ,fh2EtaJet(0)
92    ,fh2PhiJet(0)
93    ,fh2PtJet(0)
94    ,fh1PtJet(0)
95    ,fh2NtracksJet(0)
96    ,fProNtracksJet(0)
97    ,fh2EtaTrack(0)
98    ,fh2PhiTrack(0)
99    ,fh2PtTrack(0)
100    ,fh2FF(0)
101    ,fh2DelEta(0)
102    ,fh2DelPhi(0)
103    ,fh2DelR(0)
104    ,fh1PtLeadingJet(0)
105    ,fh2NtracksLeadingJet(0)
106    ,fProNtracksLeadingJet(0)
107    ,fh2DelR80pcNch(0)
108    ,fProDelR80pcNch(0)
109    ,fh2DelR80pcPt(0)
110    ,fProDelR80pcPt(0)
111    ,fh2AreaCh(0)
112    ,fProAreaCh(0)
113    ,fh3PtDelRNchSum(0)
114    ,fh3PtDelRPtSum(0)
115    ,fProDiffJetShape(0)
116    ,fProIntJetShape(0)
117 {
118   for(Int_t ii=0; ii<5; ii++){
119     fProDelRNchSum[ii]    = NULL;  
120     fProDelRPtSum[ii]     = NULL;
121     fProDiffJetShapeA[ii] = NULL;
122     fProIntJetShapeA[ii]  = NULL;
123   }//ii loop
124   // default constructor
125 }
126 //_________________________________________________________________________________//
127
128 AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties(const char *name) 
129   : AliAnalysisTaskSE(name)
130   ,fESD(0)
131   ,fAOD(0)
132   ,fAODJets(0)  
133   ,fAODExtension(0)
134   ,fNonStdFile("")
135   ,fBranchJets("jets")
136   ,fTrackType(kTrackAOD)
137   ,fJetRejectType(0)
138   ,fUseAODInputJets(kTRUE)
139   ,fFilterMask(0)
140   ,fUsePhysicsSelection(kTRUE)
141   ,fMaxVertexZ(10)
142   ,fNContributors(2)
143   ,fTrackPtCut(0)
144   ,fTrackEtaMin(0)
145   ,fTrackEtaMax(0)
146   ,fJetPtCut(0)
147   ,fJetEtaMin(0)
148   ,fJetEtaMax(0)
149   ,fAvgTrials(0)
150   ,fJetList(0)
151   ,fTrackList(0)
152   ,fCommonHistList(0)
153   ,fh1EvtSelection(0)
154   ,fh1VertexNContributors(0)
155   ,fh1VertexZ(0)
156   ,fh1Xsec(0)
157   ,fh1Trials(0)
158   ,fh1PtHard(0)
159   ,fh1PtHardTrials(0)
160   ,fh2EtaJet(0)
161   ,fh2PhiJet(0)
162   ,fh2PtJet(0)
163   ,fh1PtJet(0)
164   ,fh2NtracksJet(0)
165   ,fProNtracksJet(0)
166   ,fh2EtaTrack(0)
167   ,fh2PhiTrack(0)
168   ,fh2PtTrack(0)
169   ,fh2FF(0)
170   ,fh2DelEta(0)
171   ,fh2DelPhi(0)
172   ,fh2DelR(0)
173   ,fh1PtLeadingJet(0)
174   ,fh2NtracksLeadingJet(0)
175   ,fProNtracksLeadingJet(0)
176   ,fh2DelR80pcNch(0)
177   ,fProDelR80pcNch(0)
178   ,fh2DelR80pcPt(0)
179   ,fProDelR80pcPt(0)
180   ,fh2AreaCh(0)
181   ,fProAreaCh(0)
182   ,fh3PtDelRNchSum(0)
183   ,fh3PtDelRPtSum(0)
184   ,fProDiffJetShape(0)
185   ,fProIntJetShape(0)
186 {
187   for(Int_t ii=0; ii<5; ii++){
188     fProDelRNchSum[ii]    = NULL;  
189     fProDelRPtSum[ii]     = NULL;  
190     fProDiffJetShapeA[ii] = NULL;
191     fProIntJetShapeA[ii]  = NULL;
192   }//ii loop
193   // constructor
194   DefineOutput(1,TList::Class());
195 }
196 //_________________________________________________________________________________//
197
198 AliAnalysisTaskJetProperties::~AliAnalysisTaskJetProperties()
199 {
200   // destructor
201   if(fJetList)   delete fJetList;
202   if(fTrackList) delete fTrackList;
203 }
204 //_________________________________________________________________________________//
205
206 Bool_t AliAnalysisTaskJetProperties::Notify()
207 {
208   //
209   // Implemented Notify() to read the cross sections
210   // and number of trials from pyxsec.root
211   // (taken from AliAnalysisTaskJetSpectrum2)
212   // 
213   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Notify()");
214
215   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
216   Float_t xsection = 0;
217   Float_t ftrials  = 1;
218   
219   fAvgTrials = 1;
220   if(tree){
221     TFile *curfile = tree->GetCurrentFile();
222     if (!curfile) {
223       Error("Notify","No current file");
224       return kFALSE;
225     }
226     if(!fh1Xsec||!fh1Trials){
227       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
228       return kFALSE;
229     }
230     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
231     fh1Xsec->Fill("<#sigma>",xsection);
232     // construct a poor man average trials 
233     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
234     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
235   }
236   return kTRUE;
237 }
238 //_________________________________________________________________________________//
239
240 void AliAnalysisTaskJetProperties::UserCreateOutputObjects()
241 {
242   //(Here, event selection part is taken from AliAnalysisTaskFramentationFunction)
243   // create output objects
244   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserCreateOutputObjects()");
245
246   // create list of tracks and jets   
247   fJetList = new TList();
248   fJetList->SetOwner(kFALSE);
249   fTrackList = new TList();
250   fTrackList->SetOwner(kFALSE);
251   
252   // Create histograms / output container
253   OpenFile(1);
254   fCommonHistList = new TList();
255   fCommonHistList->SetOwner();
256   
257   Bool_t oldStatus = TH1::AddDirectoryStatus();
258   TH1::AddDirectory(kFALSE);
259   // Histograms/TProfiles       
260   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
261   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
262   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
263   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
264   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
265   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
266   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
267
268   
269   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
270   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
271   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
272   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
273   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
274   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
275   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
276   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
277   
278
279   Int_t kNbinsPtSlice=20; Float_t xMinPtSlice=0.0; Float_t xMaxPtSlice=100.0;
280   Int_t kNbinsEta=40;     Float_t xMinEta=-2.0;    Float_t xMaxEta=2.0;
281   Int_t kNbinsPhi=90;     Float_t xMinPhi=0.0;     Float_t xMaxPhi=TMath::TwoPi();
282   Int_t kNbinsPt=125;     Float_t xMinPt=0.0;      Float_t xMaxPt=250.0;
283   Int_t kNbinsNtracks=50; Float_t xMinNtracks=0.0; Float_t xMaxNtracks=50.0;
284   Int_t kNbinsFF=80;      Float_t xMinFF=0.0;      Float_t xMaxFF=2.0;
285   Int_t kNbinsDelR1D=50;  Float_t xMinDelR1D=0.0;  Float_t xMaxDelR1D=1.0;
286   
287   
288   fh2EtaJet      = new TH2F("EtaJet","EtaJet", 
289                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
290                             kNbinsEta,     xMinEta,     xMaxEta);
291   fh2PhiJet      = new TH2F("PhiJet","PhiJet",
292                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
293                             kNbinsPhi,     xMinPhi,     xMaxPhi);
294   fh2PtJet       = new TH2F("PtJet","PtJet",
295                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
296                             kNbinsPt,      xMinPt,      xMaxPt);
297   fh1PtJet       = new TH1F("PtJet1D","PtJet1D",
298                             kNbinsPt,      xMinPt,      xMaxPt);
299   fh2NtracksJet  = new TH2F("NtracksJet","NtracksJet",
300                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
301                             kNbinsNtracks, xMinNtracks, xMaxNtracks);
302   fProNtracksJet = new TProfile("AvgNoOfTracksJet","AvgNoOfTracksJet",
303                                 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
304                                 xMinNtracks, xMaxNtracks);
305   fh2EtaTrack    = new TH2F("EtaTrack","EtaTrack", 
306                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
307                             kNbinsEta,     xMinEta,     xMaxEta);
308   fh2PhiTrack    = new TH2F("PhiTrack","PhiTrack",
309                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
310                             kNbinsPhi,     xMinPhi,     xMaxPhi);
311   fh2PtTrack     = new TH2F("PtTrack","PtTrack",
312                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
313                             kNbinsPt,      xMinPt,      xMaxPt);
314   fh2FF          = new TH2F("FF","FF",
315                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
316                             kNbinsFF,      xMinFF,      xMaxFF);
317   fh2DelEta      = new TH2F("DelEta","DelEta", 
318                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
319                             kNbinsEta,     xMinEta,     xMaxEta);
320   fh2DelPhi      = new TH2F("DelPhi","DelPhi",
321                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
322                             kNbinsPhi,     xMinPhi,     xMaxPhi);
323   fh2DelR        = new TH2F("DelR","DelR",
324                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
325                             kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D);
326   
327
328   
329   Int_t kNbinsDelR=100;      Float_t xMinDelR=0.0;      Float_t xMaxDelR=1.0;
330   Int_t kNbinsPtSliceJS=100; Float_t xMinPtSliceJS=0.0; Float_t xMaxPtSliceJS=250.0;
331   
332   fh1PtLeadingJet       = new TH1F("PtLeadingJet","PtLeadingJet",
333                                    kNbinsPt,      xMinPt,      xMaxPt);
334   fh2NtracksLeadingJet  = new TH2F("NtracksLeadingJet","NtracksLeadingJet",
335                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
336                                    kNbinsNtracks,   xMinNtracks,   xMaxNtracks);
337   fProNtracksLeadingJet = new TProfile("AvgNoOfTracksLeadingJet","AvgNoOfTracksLeadingJet",
338                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
339                                        xMinNtracks, xMaxNtracks);
340   fh2DelR80pcNch        = new TH2F("delR80pcNch","delR80pcNch",
341                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
342                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
343   fProDelR80pcNch       = new TProfile("AvgdelR80pcNch","AvgdelR80pcNch",
344                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
345                                        xMinDelR,  xMaxDelR);
346   fh2DelR80pcPt         = new TH2F("delR80pcPt","delR80pcPt",
347                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
348                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
349   fProDelR80pcPt        = new TProfile("AvgdelR80pcPt","AvgdelR80pcPt",
350                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
351                                        xMinDelR,  xMaxDelR);
352   fh2AreaCh             = new TH2F("Area","Area",
353                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
354                                    72, 0.0, TMath::Pi());
355   fProAreaCh            = new TProfile("AvgArea","AvgArea",
356                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
357                                    0.0, TMath::Pi());
358   fh3PtDelRNchSum       = new TH3F("PtDelRNchSum","PtDelRNchSum",
359                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
360                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
361                                    kNbinsNtracks, xMinNtracks, xMaxNtracks);
362   fh3PtDelRPtSum        = new TH3F("PtDelRPtSum","PtDelRPtSum",
363                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
364                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
365                                    kNbinsPt,        xMinPt,        xMaxPt);
366   fProDiffJetShape      = new TProfile("DiffJetShape","DiffJetShape",
367                                   10,0.0,1.0,0.0,250.0);
368   fProIntJetShape       = new TProfile("IntJetShape","IntJetShape",
369                                   10,0.0,1.0,0.0,250.0);
370
371   TString title;
372   for(Int_t ii=0; ii<5; ii++){
373     if(ii==0)title = "_JetPt20to30";
374     if(ii==1)title = "_JetPt30to40";
375     if(ii==2)title = "_JetPt40to60";
376     if(ii==3)title = "_JetPt60to80";
377     if(ii==4)title = "_JetPt80to100";
378     fProDelRNchSum[ii] = new TProfile(Form("AvgNchSumDelR%s",title.Data()),Form("AvgNchSumDelR%s",title.Data()),
379                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
380                                       xMinNtracks, xMaxNtracks);
381     
382     fProDelRPtSum[ii]  = new TProfile(Form("AvgPtSumDelR%s",title.Data()),Form("AvgPtSumDelR%s",title.Data()),
383                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
384                                       xMinPt,      xMaxPt);
385     fProDelRNchSum[ii] ->GetXaxis()->SetTitle("R");  
386     fProDelRNchSum[ii] ->GetYaxis()->SetTitle("<NchSum>");
387     fProDelRPtSum[ii]  ->GetXaxis()->SetTitle("R");  
388     fProDelRPtSum[ii]  ->GetYaxis()->SetTitle("<PtSum>");
389
390     fProDiffJetShapeA[ii] = new TProfile(Form("DiffJetShape%s",title.Data()),Form("DiffJetShape%s",title.Data()),
391                                          10,0.0,1.0,0.0,250.0);
392     fProIntJetShapeA[ii]  = new TProfile(Form("IntJetShape%s",title.Data()),Form("IntJetShape%s",title.Data()),
393                                          10,0.0,1.0,0.0,250.0);
394
395     fProDiffJetShapeA[ii]->GetXaxis()->SetTitle("R"); 
396     fProDiffJetShapeA[ii]->GetYaxis()->SetTitle("Diff jet shape"); 
397     fProIntJetShapeA[ii]->GetXaxis()->SetTitle("R");  
398     fProIntJetShapeA[ii]->GetYaxis()->SetTitle("Integrated jet shape");  
399     
400     fCommonHistList->Add(fProDelRNchSum[ii]);
401     fCommonHistList->Add(fProDelRPtSum[ii]);
402     fCommonHistList->Add(fProDiffJetShapeA[ii]);
403     fCommonHistList->Add(fProIntJetShapeA[ii]);
404   
405 }//ii loop
406   
407   fh2EtaJet     ->GetXaxis()->SetTitle("JetPt");  fh2EtaJet     ->GetYaxis()->SetTitle("JetEta");
408   fh2PhiJet     ->GetXaxis()->SetTitle("JetPt");  fh2PhiJet     ->GetYaxis()->SetTitle("JetPhi");
409   fh2PtJet      ->GetXaxis()->SetTitle("JetPt");  fh2PtJet      ->GetYaxis()->SetTitle("JetPt");
410   fh1PtJet      ->GetXaxis()->SetTitle("JetPt");  fh1PtJet      ->GetYaxis()->SetTitle("#jets");
411   fh2NtracksJet ->GetXaxis()->SetTitle("JetPt");  fh2NtracksJet ->GetYaxis()->SetTitle("#tracks");
412   fProNtracksJet->GetXaxis()->SetTitle("JetPt");  fProNtracksJet->GetYaxis()->SetTitle("AgvNtracks");
413   fh2EtaTrack   ->GetXaxis()->SetTitle("JetPt");  fh2EtaTrack   ->GetYaxis()->SetTitle("TrackEta");
414   fh2PhiTrack   ->GetXaxis()->SetTitle("JetPt");  fh2PhiTrack   ->GetYaxis()->SetTitle("TrackPhi");
415   fh2PtTrack    ->GetXaxis()->SetTitle("JetPt");  fh2PtTrack    ->GetYaxis()->SetTitle("TrackPt");
416   fh2FF         ->GetXaxis()->SetTitle("JetPt");  fh2FF         ->GetYaxis()->SetTitle("FF");
417   fh2DelEta     ->GetXaxis()->SetTitle("JetPt");  fh2DelEta     ->GetYaxis()->SetTitle("DelEta");
418   fh2DelPhi     ->GetXaxis()->SetTitle("JetPt");  fh2DelPhi     ->GetYaxis()->SetTitle("DelPhi");
419   fh2DelR       ->GetXaxis()->SetTitle("JetPt");  fh2DelR       ->GetYaxis()->SetTitle("DelR");
420
421   fh1PtLeadingJet       ->GetXaxis()->SetTitle("JetPt");  fh1PtLeadingJet       ->GetYaxis()->SetTitle("#leading jets");
422   fh2NtracksLeadingJet  ->GetXaxis()->SetTitle("JetPt");  fh2NtracksLeadingJet  ->GetYaxis()->SetTitle("#tracks leading jet");
423   fProNtracksLeadingJet ->GetXaxis()->SetTitle("JetPt");  fProNtracksLeadingJet ->GetYaxis()->SetTitle("AvgNtracks leading jet");
424   fh2DelR80pcNch        ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcNch        ->GetYaxis()->SetTitle("R containing 80% of tracks");
425   fProDelR80pcNch       ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcNch       ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
426   fh2DelR80pcPt         ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcPt         ->GetYaxis()->SetTitle("R containing 80% of pT");
427   fProDelR80pcPt        ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcPt        ->GetYaxis()->SetTitle("<R> containing 80% of pT");
428   fh2AreaCh             ->GetXaxis()->SetTitle("JetPt");  fh2AreaCh             ->GetYaxis()->SetTitle("Jet area");
429   fProAreaCh            ->GetXaxis()->SetTitle("JetPt");  fProAreaCh            ->GetYaxis()->SetTitle("<jet area>");
430   fh3PtDelRNchSum       ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRNchSum       ->GetYaxis()->SetTitle("R");  fh3PtDelRNchSum->GetZaxis()->SetTitle("NchSum");
431   fh3PtDelRPtSum        ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRPtSum        ->GetYaxis()->SetTitle("R");  fh3PtDelRPtSum ->GetZaxis()->SetTitle("PtSum");
432   fProDiffJetShape->GetXaxis()->SetTitle("R"); 
433   fProDiffJetShape->GetYaxis()->SetTitle("Diff jet shape"); 
434   fProIntJetShape->GetXaxis()->SetTitle("R");  
435   fProIntJetShape->GetYaxis()->SetTitle("Integrated jet shape");  
436   
437   fCommonHistList->Add(fh1EvtSelection);
438   fCommonHistList->Add(fh1VertexNContributors);
439   fCommonHistList->Add(fh1VertexZ);
440   fCommonHistList->Add(fh1Xsec);
441   fCommonHistList->Add(fh1Trials);
442   fCommonHistList->Add(fh1PtHard);
443   fCommonHistList->Add(fh1PtHardTrials);
444   fCommonHistList->Add(fh2EtaJet);
445   fCommonHistList->Add(fh2PhiJet);
446   fCommonHistList->Add(fh2PtJet);
447   fCommonHistList->Add(fh1PtJet);
448   fCommonHistList->Add(fh2NtracksJet);
449   fCommonHistList->Add(fProNtracksJet);
450   fCommonHistList->Add(fh2EtaTrack);
451   fCommonHistList->Add(fh2PhiTrack); 
452   fCommonHistList->Add(fh2PtTrack); 
453   fCommonHistList->Add(fh2FF);  
454   fCommonHistList->Add(fh2DelEta);       
455   fCommonHistList->Add(fh2DelPhi);   
456   fCommonHistList->Add(fh2DelR); 
457   
458   fCommonHistList->Add(fh1PtLeadingJet); 
459   fCommonHistList->Add(fh2NtracksLeadingJet); 
460   fCommonHistList->Add(fProNtracksLeadingJet); 
461   fCommonHistList->Add(fh2DelR80pcNch); 
462   fCommonHistList->Add(fProDelR80pcNch); 
463   fCommonHistList->Add(fh2DelR80pcPt); 
464   fCommonHistList->Add(fProDelR80pcPt); 
465   fCommonHistList->Add(fh2AreaCh); 
466   fCommonHistList->Add(fProAreaCh); 
467   fCommonHistList->Add(fh3PtDelRNchSum); 
468   fCommonHistList->Add(fh3PtDelRPtSum); 
469   fCommonHistList->Add(fProDiffJetShape);
470   fCommonHistList->Add(fProIntJetShape);
471   // =========== Switch on Sumw2 for all histos ===========
472   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
473     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
474     if (h1) h1->Sumw2();
475     else{
476       TProfile *hPro = dynamic_cast<TProfile*>(fCommonHistList->At(i));
477       if(hPro) hPro->Sumw2();
478     }
479   }
480   
481   TH1::AddDirectory(oldStatus);
482   PostData(1, fCommonHistList);
483 }
484 //_________________________________________________________________________________//
485
486 void AliAnalysisTaskJetProperties::Init()
487 {
488   // Initialization
489   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Init()");
490
491 }
492 //_________________________________________________________________________________//
493
494 void AliAnalysisTaskJetProperties::UserExec(Option_t *) 
495 {
496   // Main loop
497   // Called for each event
498   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserExec()");
499   
500   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
501   // Trigger selection
502   
503   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
504     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
505   if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
506     if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
507       fh1EvtSelection->Fill(1.);
508       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
509       PostData(1, fCommonHistList);
510       return;
511     }
512   }
513   
514   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
515   if(!fESD){
516     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
517   }
518   
519   fMCEvent = MCEvent();
520   if(!fMCEvent){
521     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
522   }
523   
524   // get AOD event from input/ouput
525   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
526   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
527     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
528     if(fUseAODInputJets) fAODJets = fAOD;
529     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
530   }
531   else {
532     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
533     if( handler && handler->InheritsFrom("AliAODHandler") ) {
534       fAOD = ((AliAODHandler*)handler)->GetAOD();
535       fAODJets = fAOD;
536       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
537     }
538   }
539   
540   if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
541     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
542     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
543       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
544       if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
545     }
546   }
547   
548   if(fNonStdFile.Length()!=0){
549     // case we have an AOD extension - fetch the jets from the extended output
550     
551     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
552     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
553     if(!fAODExtension){
554       if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
555     }
556   }
557   
558   if(!fAOD){
559     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
560     return;
561   }
562   if(!fAODJets){
563     Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
564     return;
565   }
566   
567   
568   // *** vertex cut ***
569   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
570   Int_t nTracksPrim = primVtx->GetNContributors();
571   fh1VertexNContributors->Fill(nTracksPrim);
572   
573   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
574   //if(!nTracksPrim){
575   if(nTracksPrim<fNContributors){
576     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
577     fh1EvtSelection->Fill(3.);
578     PostData(1, fCommonHistList);
579     return;
580   }
581   
582   fh1VertexZ->Fill(primVtx->GetZ());
583   
584   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
585     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
586     fh1EvtSelection->Fill(4.);
587     PostData(1, fCommonHistList);
588     return; 
589   }
590   
591   TString primVtxName(primVtx->GetName());
592   
593   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
594     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
595     fh1EvtSelection->Fill(5.);
596     PostData(1, fCommonHistList);
597     return;
598   }
599   
600   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
601   fh1EvtSelection->Fill(0.);
602   //___ get MC information __________________________________________________________________
603   Double_t ptHard = 0.;
604   Double_t nTrials = 1; // trials for MC trigger weight for real data
605   
606   if(fMCEvent){
607     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
608     if(genHeader){
609       AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
610       AliGenHijingEventHeader*  hijingGenHeader = 0x0;
611       
612       if(pythiaGenHeader){
613         if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
614         nTrials = pythiaGenHeader->Trials();
615         ptHard  = pythiaGenHeader->GetPtHard();
616         
617         fh1PtHard->Fill(ptHard);
618         fh1PtHardTrials->Fill(ptHard,nTrials);
619         
620         
621       } else { // no pythia, hijing?
622         
623         if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
624         
625         hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
626         if(!hijingGenHeader){
627           Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
628         } else {
629           if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
630         }
631       }
632       
633       fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
634     }
635   }
636   
637   //___ fetch jets __________________________________________________________________________
638   Int_t nJ = GetListOfJets(fJetList);
639   Int_t nJets = 0;
640   if(nJ>=0) nJets = fJetList->GetEntries();
641   if(fDebug>2){
642     Printf("%s:%d Selected jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
643     if(nJ != nJets) Printf("%s:%d Mismatch Selected Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
644   }
645   FillJetProperties(fJetList);
646   FillJetShape(fJetList);
647   fJetList->Clear();
648   //Post output data.
649   PostData(1, fCommonHistList);
650 }
651 //_________________________________________________________________________________//
652
653 void AliAnalysisTaskJetProperties::Terminate(Option_t *) 
654 {
655   // terminated
656   
657   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::Terminate() \n");
658 }
659 //_________________________________________________________________________________//
660
661 Int_t AliAnalysisTaskJetProperties::GetListOfJets(TList *list)
662 {
663   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
664   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJets() \n");
665   // fill list of jets selected according to type
666   if(!list){
667     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
668     return -1;
669   }
670   
671   if(fBranchJets.Length()==0){
672     Printf("%s:%d no jet branch specified", (char*)__FILE__,__LINE__);
673     if(fDebug>1)fAOD->Print();
674     return 0;
675   }
676   TClonesArray *aodJets = 0; 
677   if(fBranchJets.Length())      aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchJets.Data()));
678   if(!aodJets)                  aodJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchJets.Data()));
679   if(fAODExtension&&!aodJets)   aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchJets.Data()));
680   
681   if(!aodJets){
682     if(fBranchJets.Length())Printf("%s:%d no jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchJets.Data());
683     if(fDebug>1)fAOD->Print();
684     return 0;
685   }
686   
687   Int_t nJets = 0;
688   for(Int_t ijet=0; ijet<aodJets->GetEntries(); ijet++){
689     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodJets->At(ijet));
690     if(!tmp) continue;
691     if( tmp->Pt() < fJetPtCut ) continue;
692     if( tmp->Eta() < fJetEtaMin || tmp->Eta() > fJetEtaMax)continue;
693     if(fJetRejectType==kReject1Track && tmp->GetRefTracks()->GetEntriesFast()==1)continue;//rejecting 1track jet if...
694     list->Add(tmp);
695     nJets++;
696   }//ij loop
697   list->Sort();
698   return nJets;
699 }
700 //_________________________________________________________________________________//
701
702 Int_t AliAnalysisTaskJetProperties::GetListOfJetTracks(TList* list, const AliAODJet* jet)
703 {
704   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
705   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJetTracks() \n");
706   
707   // list of jet tracks from trackrefs
708   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
709   Int_t NoOfTracks=0;
710   for (Int_t itrack=0; itrack<nTracks; itrack++) {
711     if(fTrackType==kTrackUndef){
712       if(fDebug>2)Printf("%s:%d unknown track type %d in AOD", (char*)__FILE__,__LINE__,kTrackUndef);
713       return 0;
714     }//if
715     else if(fTrackType == kTrackAOD){
716       AliAODTrack* trackAod = dynamic_cast<AliAODTrack*>(jet->GetRefTracks()->At(itrack));
717       if(!trackAod){
718         AliError("expected ref track not found ");
719         continue;
720       }
721       list->Add(trackAod);
722       NoOfTracks++;
723     }//if
724     
725     else if(fTrackType==kTrackAODMC){
726       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(jet->GetRefTracks()->At(itrack));
727       if(!trackmc){
728         AliError("expected ref trackmc not found ");
729         continue;
730       }
731       list->Add(trackmc);
732       NoOfTracks++;
733     }//if
734     
735     else if (fTrackType==kTrackKine){
736       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
737       if(!trackkine){
738         AliError("expected ref trackkine not found ");
739         continue;
740       }
741       list->Add(trackkine);
742       NoOfTracks++;
743     }//if
744   }//itrack loop
745   list->Sort();
746   return NoOfTracks;
747 }//initial loop
748 //_________________________________________________________________________________//
749
750 void AliAnalysisTaskJetProperties::FillJetProperties(TList *jetlist){
751   //filling up the histograms jets and tracks inside jet
752   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetProperties() \n");
753   
754   for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
755     Float_t JetPt;Float_t JetPhi;Float_t JetEta;Float_t JetE;
756     AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
757     if(!jet)continue;
758     JetEta = jet->Eta();
759     JetPhi = jet->Phi();
760     JetPt  = jet->Pt();
761     JetE   = jet->E();
762     fh2EtaJet ->Fill(JetPt,JetEta);
763     fh2PhiJet ->Fill(JetPt,JetPhi);
764     fh2PtJet  ->Fill(JetPt,JetPt);
765     fh1PtJet  ->Fill(JetPt);
766     
767     Int_t nJT = GetListOfJetTracks(fTrackList,jet);
768     Int_t nJetTracks = 0;
769     if(nJT>=0) nJetTracks = fTrackList->GetEntries();
770     if(fDebug>2){
771       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
772       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
773     }
774     
775     fh2NtracksJet  ->Fill(JetPt,fTrackList->GetEntries());
776     fProNtracksJet ->Fill(JetPt,fTrackList->GetEntries());
777     
778     for (Int_t j =0; j< fTrackList->GetEntries(); j++){
779       if(fTrackType==kTrackUndef)continue;
780         Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
781         Float_t FF=-99.0;       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
782         if(fTrackType==kTrackAOD){
783           AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackList->At(j));
784           if(!trackaod)continue;
785           TrackEta = trackaod->Eta();
786           TrackPhi = trackaod->Phi();
787           TrackPt  = trackaod->Pt();
788         }//if kTrackAOD
789         else if(fTrackType==kTrackAODMC){
790           AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackList->At(j));
791           if(!trackmc)continue;
792           TrackEta = trackmc->Eta();
793           TrackPhi = trackmc->Phi();
794           TrackPt  = trackmc->Pt();
795         }//if kTrackAODMC
796         else if(fTrackType==kTrackKine){
797           AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackList->At(j));
798           if(!trackkine)continue;
799           TrackEta = trackkine->Eta();
800           TrackPhi = trackkine->Phi();
801           TrackPt  = trackkine->Pt();
802         }//kTrackKine
803         if(JetPt)FF = TrackPt/JetPt;
804         DelEta      = TMath::Abs(JetEta - TrackEta);
805         DelPhi      = TMath::Abs(JetPhi - TrackPhi);
806         if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
807         DelR        = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
808         AreaJ       = TMath::Pi()*DelR*DelR;
809         fh2EtaTrack ->Fill(JetPt,TrackEta);
810         fh2PhiTrack ->Fill(JetPt,TrackPhi);
811         fh2PtTrack  ->Fill(JetPt,TrackPt);
812         fh2FF       ->Fill(JetPt,FF);
813         fh2DelEta   ->Fill(JetPt,DelEta);
814         fh2DelPhi   ->Fill(JetPt,DelPhi);
815         fh2DelR     ->Fill(JetPt,DelR);
816     }//track loop
817     fTrackList->Clear();
818   }//iJet loop
819 }//FillJetProperties
820 //_________________________________________________________________________________//
821
822 void AliAnalysisTaskJetProperties::FillJetShape(TList *jetlist){
823   //filling up the histograms
824   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
825   Float_t JetEta; Float_t JetPhi; Float_t JetPt;
826   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet
827   if(jet){
828     JetEta = jet->Eta();
829     JetPhi = jet->Phi();
830     JetPt  = jet->Pt();
831     fh1PtLeadingJet->Fill(JetPt);
832     Float_t NchSumA[50]         = {0.};
833     Float_t PtSumA[50]          = {0.};
834     Float_t delRPtSum80pc       = 0;
835     Float_t delRNtrkSum80pc     = 0;
836     Float_t PtSumDiffShape[10]  = {0.0};
837     Float_t PtSumIntShape[10]   = {0.0};
838     Int_t kNbinsR               = 10;
839
840     Int_t nJT = GetListOfJetTracks(fTrackList,jet);
841     Int_t nJetTracks = 0;
842     if(nJT>=0) nJetTracks = fTrackList->GetEntries();
843     if(fDebug>2){
844       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
845       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
846     }
847     fh2NtracksLeadingJet->Fill(JetPt,nJetTracks);
848     fProNtracksLeadingJet->Fill(JetPt,nJetTracks);
849     Int_t   *index     = new Int_t   [nJetTracks];//dynamic array containing index
850     Float_t *delRA     = new Float_t [nJetTracks];//dynamic array containing delR
851     Float_t *delEtaA   = new Float_t [nJetTracks];//dynamic array containing delEta
852     Float_t *delPhiA   = new Float_t [nJetTracks];//dynamic array containing delPhi
853     Float_t *trackPtA  = new Float_t [nJetTracks];//dynamic array containing pt-track
854     Float_t *trackEtaA = new Float_t [nJetTracks];//dynamic array containing eta-track
855     Float_t *trackPhiA = new Float_t [nJetTracks];//dynamic array containing phi-track
856     for(Int_t ii=0; ii<nJetTracks; ii++){
857       index[ii]     = 0;
858       delRA[ii]     = 0.;
859       delEtaA[ii]   = 0.;
860       delPhiA[ii]   = 0.;
861       trackPtA[ii]  = 0.;
862       trackEtaA[ii] = 0.;
863       trackPhiA[ii] = 0.;
864     }//ii loop
865   
866     for (Int_t j =0; j< nJetTracks; j++){
867       if(fTrackType==kTrackUndef)continue;      
868       
869       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
870       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
871       
872       if(fTrackType==kTrackAOD){
873         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackList->At(j));
874         if(!trackaod)continue;
875         TrackEta = trackaod->Eta();
876         TrackPhi = trackaod->Phi();
877         TrackPt  = trackaod->Pt();
878       }//if kTrackAOD
879       else if(fTrackType==kTrackAODMC){
880         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackList->At(j));
881         if(!trackmc)continue;
882         TrackEta = trackmc->Eta();
883         TrackPhi = trackmc->Phi();
884         TrackPt  = trackmc->Pt();
885       }//if kTrackAODMC
886       else if(fTrackType==kTrackKine){
887         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackList->At(j));
888         if(!trackkine)continue;
889         TrackEta = trackkine->Eta();
890         TrackPhi = trackkine->Phi();
891         TrackPt  = trackkine->Pt();
892       }//if kTrackKine
893       
894       DelEta = TMath::Abs(JetEta - TrackEta);
895       DelPhi = TMath::Abs(JetPhi - TrackPhi);
896       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
897       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
898       AreaJ  = TMath::Pi()*DelR*DelR;
899       
900       fh2AreaCh ->Fill(JetPt,AreaJ);
901       fProAreaCh->Fill(JetPt,AreaJ);
902       delRA[j]     = DelR;
903       delEtaA[j]   = DelEta;
904       delPhiA[j]   = DelPhi;
905       trackPtA[j]  = TrackPt;
906       trackEtaA[j] = TrackEta;
907       trackPhiA[j] = TrackPhi;
908       
909       //calculating diff and integrated jet shapes
910       Float_t kDeltaR = 0.1;
911       Float_t RMin    = kDeltaR/2.0;
912       Float_t RMax    = kDeltaR/2.0;
913       Float_t tmpR    = 0.05;
914       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
915         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
916         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
917         tmpR += 0.1;
918       }//ii1 loop
919       
920       for(Int_t ibin=1; ibin<=50; ibin++){
921         Float_t xlow = 0.02*(ibin-1);
922         Float_t xup  = 0.02*ibin;
923         if( xlow <= DelR && DelR < xup){
924           NchSumA[ibin-1]++;
925           PtSumA[ibin-1]+= TrackPt;
926         }//if loop
927       }//for ibin loop
928     }//track loop
929     fTrackList->Clear();
930     
931     //---------------------
932     Float_t tmp1R = 0.05;
933     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
934       if(JetPt>20 && JetPt<=100){
935         fProDiffJetShape->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
936         fProIntJetShape ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
937       }
938       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
939       for(Int_t k=0; k<5; k++){
940         if(k==0){jetPtMin0=20.0;jetPtMax0=30.0;}
941         if(k==1){jetPtMin0=30.0;jetPtMax0=40.0;}
942         if(k==2){jetPtMin0=40.0;jetPtMax0=60.0;}
943         if(k==3){jetPtMin0=60.0;jetPtMax0=80.0;}
944         if(k==4){jetPtMin0=80.0;jetPtMax0=100.0;}
945         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
946           fProDiffJetShapeA[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
947           fProIntJetShapeA[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
948         }//if
949       }//k loop
950       tmp1R +=0.1;
951     }//jj1 loop
952     //----------------------//
953     Float_t PtSum = 0;
954     Int_t NtrkSum = 0;
955     Bool_t iflagPtSum   = kFALSE;
956     Bool_t iflagNtrkSum = kFALSE;
957     TMath::Sort(nJetTracks,delRA,index,0);
958     for(Int_t ii=0; ii<nJetTracks; ii++){
959       NtrkSum ++;
960       PtSum += trackPtA[index[ii]];
961       /*
962         cout << index[ii] << "\t" <<
963         delR[ii] << "\t" <<
964         delEta[ii]<< "\t" <<
965         delPhi[ii]<< "\t" <<
966         trackPt[ii]<< "\t" <<
967         trackEta[ii]<< "\t" <<
968         trackPhi[ii]<< "\t DelR " <<
969         delR[index[ii]] << endl;
970       */
971       if(!iflagNtrkSum){
972         if((Float_t)NtrkSum/(Float_t)nJetTracks > 0.79){
973           delRNtrkSum80pc = delRA[index[ii]];
974           iflagNtrkSum = kTRUE;
975         }//if NtrkSum
976       }//if iflag
977       if(!iflagPtSum){
978         if(PtSum/JetPt >= 0.8000){
979           delRPtSum80pc = delRA[index[ii]];
980           iflagPtSum = kTRUE;
981         }//if PtSum
982       }//if iflag
983     }//track loop 2nd
984     delete [] index;
985     delete [] delRA;
986     delete [] delEtaA;
987     delete [] delPhiA;
988     delete [] trackPtA;
989     delete [] trackEtaA;
990     delete [] trackPhiA;
991     fh2DelR80pcNch ->Fill(JetPt,delRNtrkSum80pc);
992     fProDelR80pcNch->Fill(JetPt,delRNtrkSum80pc);
993     fh2DelR80pcPt  ->Fill(JetPt,delRPtSum80pc);
994     fProDelR80pcPt ->Fill(JetPt,delRPtSum80pc);
995     for(Int_t ibin=0; ibin<50; ibin++){
996       Float_t iR = 0.02*ibin + 0.01;
997       fh3PtDelRNchSum->Fill(JetPt,iR,NchSumA[ibin]);
998       fh3PtDelRPtSum ->Fill(JetPt,iR,PtSumA[ibin]);
999       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1000       for(Int_t k=0; k<5; k++){
1001         if(k==0){jetPtMin=20.0;jetPtMax=30.0;}
1002         if(k==1){jetPtMin=30.0;jetPtMax=40.0;}
1003         if(k==2){jetPtMin=40.0;jetPtMax=60.0;}
1004         if(k==3){jetPtMin=60.0;jetPtMax=80.0;}
1005         if(k==4){jetPtMin=80.0;jetPtMax=100.0;}
1006         if(JetPt>jetPtMin && JetPt<jetPtMax){
1007           fProDelRNchSum[k]->Fill(iR,NchSumA[ibin]);
1008           fProDelRPtSum[k]->Fill(iR,PtSumA[ibin]);
1009         }//if
1010       }//k loop
1011     }//ibin loop
1012   }//if(jet)
1013 }//FillJetShape()
1014 //_________________________________________________________________________________//
1015
1016