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