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