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