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   ,fJetRadius(0.4)
  82   ,fJetList(0)
  83   ,fTrackList(0)
  84   ,fTrackListUE(0)
  85   ,fTrackListJet(0)
  86   ,fCommonHistList(0)
  87   ,fh1EvtSelection(0)
  88   ,fh1VertexNContributors(0)
  89   ,fh1VertexZ(0)
  90   ,fh1Xsec(0)
  91   ,fh1Trials(0)
  92   ,fh1PtHard(0)
  93   ,fh1PtHardTrials(0)
  94   ,fh2EtaJet(0)
  95   ,fh2PhiJet(0)
  96   ,fh2PtJet(0)
  97   ,fh1PtJet(0)
  98   ,fh2NtracksJet(0)
  99   ,fProNtracksJet(0)
 100   ,fh2EtaTrack(0)
 101   ,fh2PhiTrack(0)
 102   ,fh2PtTrack(0)
 103   ,fh2FF(0)
 104   ,fh2DelEta(0)
 105   ,fh2DelPhi(0)
 106   ,fh2DelR(0)
 107
 108   ,fh1PtLeadingJet(0)
 109   ,fh2NtracksLeadingJet(0)
 110   ,fProNtracksLeadingJet(0)
 111   ,fh2DelR80pcNch(0)
 112   ,fProDelR80pcNch(0)
 113   ,fh2DelR80pcPt(0)
 114   ,fProDelR80pcPt(0)
 115   ,fh2AreaCh(0)
 116   ,fProAreaCh(0)
 117   ,fh3PtDelRNchSum(0)
 118   ,fh3PtDelRPtSum(0)
 119   ,fProDiffJetShape(0)
 120   ,fProIntJetShape(0)
 121
 122   ,fh1PtSumInJetConeUE(0)
 123   ,fh2NtracksLeadingJetUE(0)
 124   ,fProNtracksLeadingJetUE(0)
 125   ,fh2DelR80pcNchUE(0)
 126   ,fProDelR80pcNchUE(0)
 127   ,fh2DelR80pcPtUE(0)
 128   ,fProDelR80pcPtUE(0)
 129   ,fh2AreaChUE(0)
 130   ,fProAreaChUE(0)
 131   ,fh3PtDelRNchSumUE(0)
 132   ,fh3PtDelRPtSumUE(0)
 133   ,fProDiffJetShapeUE(0)
 134   ,fProIntJetShapeUE(0)
 135 {
 136   for(Int_t ii=0; ii<13; ii++){
 137     fProDelRNchSum[ii]    = NULL;
 138     fProDelRPtSum[ii]     = NULL;
 139     fProDiffJetShapeA[ii] = NULL;
 140     fProIntJetShapeA[ii]  = NULL;
 141
 142     fProDelRNchSumUE[ii]    = NULL;
 143     fProDelRPtSumUE[ii]     = NULL;
 144     fProDiffJetShapeAUE[ii] = NULL;
 145     fProIntJetShapeAUE[ii]  = NULL;
 146   }//ii loop
 147   // default constructor
 148 }
 149 //_________________________________________________________________________________//
 150
 151 AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties(const char *name)
 152   : AliAnalysisTaskSE(name)
 153   ,fESD(0)
 154   ,fAOD(0)
 155   ,fAODJets(0)
 156   ,fAODExtension(0)
 157   ,fNonStdFile("")
 158   ,fBranchJets("jets")
 159   ,fTrackType(kTrackAOD)
 160   ,fJetRejectType(0)
 161   ,fUseAODInputJets(kTRUE)
 162   ,fFilterMask(0)
 163   ,fUsePhysicsSelection(kTRUE)
 164   ,fMaxVertexZ(10)
 165   ,fNContributors(2)
 166   ,fTrackPtCut(0)
 167   ,fTrackEtaMin(0)
 168   ,fTrackEtaMax(0)
 169   ,fJetPtCut(0)
 170   ,fJetEtaMin(0)
 171   ,fJetEtaMax(0)
 172   ,fAvgTrials(0)
 173   ,fJetRadius(0.4)
 174   ,fJetList(0)
 175   ,fTrackList(0)
 176   ,fTrackListUE(0)
 177   ,fTrackListJet(0)
 178   ,fCommonHistList(0)
 179   ,fh1EvtSelection(0)
 180   ,fh1VertexNContributors(0)
 181   ,fh1VertexZ(0)
 182   ,fh1Xsec(0)
 183   ,fh1Trials(0)
 184   ,fh1PtHard(0)
 185   ,fh1PtHardTrials(0)
 186   ,fh2EtaJet(0)
 187   ,fh2PhiJet(0)
 188   ,fh2PtJet(0)
 189   ,fh1PtJet(0)
 190   ,fh2NtracksJet(0)
 191   ,fProNtracksJet(0)
 192   ,fh2EtaTrack(0)
 193   ,fh2PhiTrack(0)
 194   ,fh2PtTrack(0)
 195   ,fh2FF(0)
 196   ,fh2DelEta(0)
 197   ,fh2DelPhi(0)
 198   ,fh2DelR(0)
 199
 200   ,fh1PtLeadingJet(0)
 201   ,fh2NtracksLeadingJet(0)
 202   ,fProNtracksLeadingJet(0)
 203   ,fh2DelR80pcNch(0)
 204   ,fProDelR80pcNch(0)
 205   ,fh2DelR80pcPt(0)
 206   ,fProDelR80pcPt(0)
 207   ,fh2AreaCh(0)
 208   ,fProAreaCh(0)
 209   ,fh3PtDelRNchSum(0)
 210   ,fh3PtDelRPtSum(0)
 211   ,fProDiffJetShape(0)
 212   ,fProIntJetShape(0)
 213
 214   ,fh1PtSumInJetConeUE(0)
 215   ,fh2NtracksLeadingJetUE(0)
 216   ,fProNtracksLeadingJetUE(0)
 217   ,fh2DelR80pcNchUE(0)
 218   ,fProDelR80pcNchUE(0)
 219   ,fh2DelR80pcPtUE(0)
 220   ,fProDelR80pcPtUE(0)
 221   ,fh2AreaChUE(0)
 222   ,fProAreaChUE(0)
 223   ,fh3PtDelRNchSumUE(0)
 224   ,fh3PtDelRPtSumUE(0)
 225   ,fProDiffJetShapeUE(0)
 226   ,fProIntJetShapeUE(0)
 227 {
 228   for(Int_t ii=0; ii<13; ii++){
 229     fProDelRNchSum[ii]    = NULL;
 230     fProDelRPtSum[ii]     = NULL;
 231     fProDiffJetShapeA[ii] = NULL;
 232     fProIntJetShapeA[ii]  = NULL;
 233
 234     fProDelRNchSumUE[ii]    = NULL;
 235     fProDelRPtSumUE[ii]     = NULL;
 236     fProDiffJetShapeAUE[ii] = NULL;
 237     fProIntJetShapeAUE[ii]  = NULL;
 238   }//ii loop
 239   // constructor
 240   DefineOutput(1,TList::Class());
 241 }
 242 //_________________________________________________________________________________//
 243
 244 AliAnalysisTaskJetProperties::~AliAnalysisTaskJetProperties()
 245 {
 246   // destructor
 247   if(fJetList)   delete fJetList;
 248   if(fTrackList) delete fTrackList;
 249   if(fTrackListUE) delete fTrackListUE;
 250   if(fTrackListJet) delete fTrackListJet;
 251 }
 252 //_________________________________________________________________________________//
 253
 254 Bool_t AliAnalysisTaskJetProperties::Notify()
 255 {
 256   //
 257   // Implemented Notify() to read the cross sections
 258   // and number of trials from pyxsec.root
 259   // (taken from AliAnalysisTaskJetSpectrum2)
 260   //
 261   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Notify()");
 262
 263   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 264   Float_t xsection = 0;
 265   Float_t ftrials  = 1;
 266
 267   fAvgTrials = 1;
 268   if(tree){
 269     TFile *curfile = tree->GetCurrentFile();
 270     if (!curfile) {
 271       Error("Notify","No current file");
 272       return kFALSE;
 273     }
 274     if(!fh1Xsec||!fh1Trials){
 275       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 276       return kFALSE;
 277     }
 278     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 279     fh1Xsec->Fill("<#sigma>",xsection);
 280     // construct a poor man average trials
 281     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 282     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 283   }
 284   return kTRUE;
 285 }
 286 //_________________________________________________________________________________//
 287
 288 void AliAnalysisTaskJetProperties::UserCreateOutputObjects()
 289 {
 290   //(Here, event selection part is taken from AliAnalysisTaskFramentationFunction)
 291   // create output objects
 292   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserCreateOutputObjects()");
 293
 294   // create list of tracks and jets
 295   fJetList = new TList();
 296   fJetList->SetOwner(kFALSE);
 297   fTrackList = new TList();
 298   fTrackList->SetOwner(kFALSE);
 299   fTrackListUE = new TList();
 300   fTrackListUE->SetOwner(kFALSE);
 301   fTrackListJet = new TList();
 302   fTrackListJet->SetOwner(kFALSE);
 303
 304   // Create histograms / output container
 305   OpenFile(1);
 306   fCommonHistList = new TList();
 307   fCommonHistList->SetOwner();
 308
 309   Bool_t oldStatus = TH1::AddDirectoryStatus();
 310   TH1::AddDirectory(kFALSE);
 311   // Histograms/TProfiles
 312   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
 313   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 314   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
 315   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
 316   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
 317   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
 318   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
 319
 320
 321   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
 322   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
 323   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 324   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 325   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
 326   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 327   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 328   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 329
 330
 331   Int_t kNbinsPtSlice=20; Float_t xMinPtSlice=0.0; Float_t xMaxPtSlice=100.0;
 332   Int_t kNbinsEta=40;     Float_t xMinEta=-2.0;    Float_t xMaxEta=2.0;
 333   Int_t kNbinsPhi=90;     Float_t xMinPhi=0.0;     Float_t xMaxPhi=TMath::TwoPi();
 334   Int_t kNbinsPt=125;     Float_t xMinPt=0.0;      Float_t xMaxPt=250.0;
 335   Int_t kNbinsNtracks=50; Float_t xMinNtracks=0.0; Float_t xMaxNtracks=50.0;
 336   Int_t kNbinsFF=80;      Float_t xMinFF=0.0;      Float_t xMaxFF=2.0;
 337   Int_t kNbinsDelR1D=50;  Float_t xMinDelR1D=0.0;  Float_t xMaxDelR1D=1.0;
 338
 339
 340   fh2EtaJet      = new TH2F("EtaJet","EtaJet",
 341                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 342                             kNbinsEta,     xMinEta,     xMaxEta);
 343   fh2PhiJet      = new TH2F("PhiJet","PhiJet",
 344                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 345                             kNbinsPhi,     xMinPhi,     xMaxPhi);
 346   fh2PtJet       = new TH2F("PtJet","PtJet",
 347                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 348                             kNbinsPt,      xMinPt,      xMaxPt);
 349   fh1PtJet       = new TH1F("PtJet1D","PtJet1D",
 350                             kNbinsPt,      xMinPt,      xMaxPt);
 351   fh2NtracksJet  = new TH2F("NtracksJet","NtracksJet",
 352                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 353                             kNbinsNtracks, xMinNtracks, xMaxNtracks);
 354   fProNtracksJet = new TProfile("AvgNoOfTracksJet","AvgNoOfTracksJet",
 355                                 kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 356                                 xMinNtracks, xMaxNtracks);
 357   fh2EtaTrack    = new TH2F("EtaTrack","EtaTrack",
 358                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 359                             kNbinsEta,     xMinEta,     xMaxEta);
 360   fh2PhiTrack    = new TH2F("PhiTrack","PhiTrack",
 361                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 362                             kNbinsPhi,     xMinPhi,     xMaxPhi);
 363   fh2PtTrack     = new TH2F("PtTrack","PtTrack",
 364                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 365                             kNbinsPt,      xMinPt,      xMaxPt);
 366   fh2FF          = new TH2F("FF","FF",
 367                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 368                             kNbinsFF,      xMinFF,      xMaxFF);
 369   fh2DelEta      = new TH2F("DelEta","DelEta",
 370                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 371                             kNbinsEta,     xMinEta,     xMaxEta);
 372   fh2DelPhi      = new TH2F("DelPhi","DelPhi",
 373                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 374                             kNbinsPhi,     xMinPhi,     xMaxPhi);
 375   fh2DelR        = new TH2F("DelR","DelR",
 376                             kNbinsPtSlice, xMinPtSlice, xMaxPtSlice,
 377                             kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D);
 378
 379
 380
 381   Int_t kNbinsDelR=100;      Float_t xMinDelR=0.0;      Float_t xMaxDelR=1.0;
 382   Int_t kNbinsPtSliceJS=100; Float_t xMinPtSliceJS=0.0; Float_t xMaxPtSliceJS=250.0;
 383
 384   fh1PtLeadingJet       = new TH1F("PtLeadingJet","PtLeadingJet",
 385                                    kNbinsPt,      xMinPt,      xMaxPt);
 386   fh2NtracksLeadingJet  = new TH2F("NtracksLeadingJet","NtracksLeadingJet",
 387                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 388                                    kNbinsNtracks,   xMinNtracks,   xMaxNtracks);
 389   fProNtracksLeadingJet = new TProfile("AvgNoOfTracksLeadingJet","AvgNoOfTracksLeadingJet",
 390                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 391                                        xMinNtracks, xMaxNtracks);
 392   fh2DelR80pcNch        = new TH2F("delR80pcNch","delR80pcNch",
 393                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 394                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
 395   fProDelR80pcNch       = new TProfile("AvgdelR80pcNch","AvgdelR80pcNch",
 396                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 397                                        xMinDelR,  xMaxDelR);
 398   fh2DelR80pcPt         = new TH2F("delR80pcPt","delR80pcPt",
 399                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 400                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
 401   fProDelR80pcPt        = new TProfile("AvgdelR80pcPt","AvgdelR80pcPt",
 402                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 403                                        xMinDelR,  xMaxDelR);
 404   fh2AreaCh             = new TH2F("Area","Area",
 405                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 406                                    72, 0.0, TMath::Pi());
 407   fProAreaCh            = new TProfile("AvgArea","AvgArea",
 408                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 409                                    0.0, TMath::Pi());
 410   fh3PtDelRNchSum       = new TH3F("PtDelRNchSum","PtDelRNchSum",
 411                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 412                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
 413                                    kNbinsNtracks, xMinNtracks, xMaxNtracks);
 414   fh3PtDelRPtSum        = new TH3F("PtDelRPtSum","PtDelRPtSum",
 415                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 416                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
 417                                    kNbinsPt,        xMinPt,        xMaxPt);
 418   fProDiffJetShape      = new TProfile("DiffJetShape","DiffJetShape",
 419                                   10,0.0,1.0,0.0,250.0);
 420   fProIntJetShape       = new TProfile("IntJetShape","IntJetShape",
 421                                   10,0.0,1.0,0.0,250.0);
 422
 423
 424   fh1PtSumInJetConeUE       = new TH1F("PtSumInJetConeUE","PtSumInJetConeUE",
 425                                        500,      0.,      100.);
 426   fh2NtracksLeadingJetUE  = new TH2F("NtracksLeadingJetUE","NtracksLeadingJetUE",
 427                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 428                                    kNbinsNtracks,   xMinNtracks,   xMaxNtracks);
 429   fProNtracksLeadingJetUE = new TProfile("AvgNoOfTracksLeadingJetUE","AvgNoOfTracksLeadingJetUE",
 430                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 431                                        xMinNtracks, xMaxNtracks);
 432   fh2DelR80pcNchUE        = new TH2F("delR80pcNchUE","delR80pcNchUE",
 433                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 434                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
 435   fProDelR80pcNchUE       = new TProfile("AvgdelR80pcNchUE","AvgdelR80pcNchUE",
 436                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 437                                        xMinDelR,  xMaxDelR);
 438   fh2DelR80pcPtUE         = new TH2F("delR80pcPtUE","delR80pcPtUE",
 439                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 440                                    kNbinsDelR,      xMinDelR,      xMaxDelR);
 441   fProDelR80pcPtUE        = new TProfile("AvgdelR80pcPtUE","AvgdelR80pcPtUE",
 442                                        kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 443                                        xMinDelR,  xMaxDelR);
 444   fh2AreaChUE             = new TH2F("AreaUE","AreaUE",
 445                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 446                                    72, 0.0, TMath::Pi());
 447   fProAreaChUE            = new TProfile("AvgAreaUE","AvgAreaUE",
 448                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 449                                    0.0, TMath::Pi());
 450   fh3PtDelRNchSumUE       = new TH3F("PtDelRNchSumUE","PtDelRNchSumUE",
 451                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 452                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
 453                                    kNbinsNtracks, xMinNtracks, xMaxNtracks);
 454   fh3PtDelRPtSumUE        = new TH3F("PtDelRPtSumUE","PtDelRPtSumUE",
 455                                    kNbinsPtSliceJS, xMinPtSliceJS, xMaxPtSliceJS,
 456                                    kNbinsDelR1D,    xMinDelR1D,    xMaxDelR1D,
 457                                    kNbinsPt,        xMinPt,        xMaxPt);
 458   fProDiffJetShapeUE      = new TProfile("DiffJetShapeUE","DiffJetShapeUE",
 459                                   10,0.0,1.0,0.0,250.0);
 460   fProIntJetShapeUE       = new TProfile("IntJetShapeUE","IntJetShapeUE",
 461                                   10,0.0,1.0,0.0,250.0);
 462
 463   TString title;
 464   for(Int_t ii=0; ii<13; ii++){
 465     if(ii==0)title = "_JetPt20to25";
 466     if(ii==1)title = "_JetPt25to30";
 467     if(ii==2)title = "_JetPt20to30";
 468
 469     if(ii==3)title = "_JetPt30to35";
 470     if(ii==4)title = "_JetPt35to40";
 471     if(ii==5)title = "_JetPt30to40";
 472
 473     if(ii==6)title = "_JetPt40to50";
 474     if(ii==7)title = "_JetPt50to60";
 475     if(ii==8)title = "_JetPt40to60";
 476
 477     if(ii==9) title = "_JetPt60to70";
 478     if(ii==10)title = "_JetPt70to80";
 479     if(ii==11)title = "_JetPt60to80";
 480
 481     if(ii==12)title = "_JetPt80to100";
 482     fProDelRNchSum[ii] = new TProfile(Form("AvgNchSumDelR%s",title.Data()),Form("AvgNchSumDelR%s",title.Data()),
 483                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
 484                                       xMinNtracks, xMaxNtracks);
 485
 486     fProDelRPtSum[ii]  = new TProfile(Form("AvgPtSumDelR%s",title.Data()),Form("AvgPtSumDelR%s",title.Data()),
 487                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
 488                                       xMinPt,      xMaxPt);
 489     fProDelRNchSum[ii] ->GetXaxis()->SetTitle("R");
 490     fProDelRNchSum[ii] ->GetYaxis()->SetTitle("<NchSum>");
 491     fProDelRPtSum[ii]  ->GetXaxis()->SetTitle("R");
 492     fProDelRPtSum[ii]  ->GetYaxis()->SetTitle("<PtSum>");
 493
 494     fProDiffJetShapeA[ii] = new TProfile(Form("DiffJetShape%s",title.Data()),Form("DiffJetShape%s",title.Data()),
 495                                          10,0.0,1.0,0.0,250.0);
 496     fProIntJetShapeA[ii]  = new TProfile(Form("IntJetShape%s",title.Data()),Form("IntJetShape%s",title.Data()),
 497                                          10,0.0,1.0,0.0,250.0);
 498
 499     fProDiffJetShapeA[ii]->GetXaxis()->SetTitle("R");
 500     fProDiffJetShapeA[ii]->GetYaxis()->SetTitle("Diff jet shape");
 501     fProIntJetShapeA[ii]->GetXaxis()->SetTitle("R");
 502     fProIntJetShapeA[ii]->GetYaxis()->SetTitle("Integrated jet shape");
 503
 504     fCommonHistList->Add(fProDelRNchSum[ii]);
 505     fCommonHistList->Add(fProDelRPtSum[ii]);
 506     fCommonHistList->Add(fProDiffJetShapeA[ii]);
 507     fCommonHistList->Add(fProIntJetShapeA[ii]);
 508
 509     fProDelRNchSumUE[ii] = new TProfile(Form("AvgNchSumDelR%sUE",title.Data()),Form("AvgNchSumDelR%sUE",title.Data()),
 510                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
 511                                       xMinNtracks, xMaxNtracks);
 512
 513     fProDelRPtSumUE[ii]  = new TProfile(Form("AvgPtSumDelR%sUE",title.Data()),Form("AvgPtSumDelR%sUE",title.Data()),
 514                                       kNbinsDelR1D,  xMinDelR1D,  xMaxDelR1D,
 515                                       xMinPt,      xMaxPt);
 516     fProDelRNchSumUE[ii] ->GetXaxis()->SetTitle("R");
 517     fProDelRNchSumUE[ii] ->GetYaxis()->SetTitle("<N_{ch}^{Sum,UE}>");
 518     fProDelRPtSumUE[ii]  ->GetXaxis()->SetTitle("R");
 519     fProDelRPtSumUE[ii]  ->GetYaxis()->SetTitle("<p_{T}^{sum, UE}>");
 520
 521     fProDiffJetShapeAUE[ii] = new TProfile(Form("DiffJetShape%sUE",title.Data()),Form("DiffJetShape%sUE",title.Data()),
 522                                          10,0.0,1.0,0.0,250.0);
 523     fProIntJetShapeAUE[ii]  = new TProfile(Form("IntJetShape%sUE",title.Data()),Form("IntJetShape%sUE",title.Data()),
 524                                          10,0.0,1.0,0.0,250.0);
 525
 526     fProDiffJetShapeAUE[ii]->GetXaxis()->SetTitle("R");
 527     fProDiffJetShapeAUE[ii]->GetYaxis()->SetTitle("Diff jet shape UE");
 528     fProIntJetShapeAUE[ii]->GetXaxis()->SetTitle("R");
 529     fProIntJetShapeAUE[ii]->GetYaxis()->SetTitle("Integrated jet shape UE");
 530
 531     fCommonHistList->Add(fProDelRNchSumUE[ii]);
 532     fCommonHistList->Add(fProDelRPtSumUE[ii]);
 533     fCommonHistList->Add(fProDiffJetShapeAUE[ii]);
 534     fCommonHistList->Add(fProIntJetShapeAUE[ii]);
 535   }//ii loop
 536
 537   fh2EtaJet     ->GetXaxis()->SetTitle("JetPt");  fh2EtaJet     ->GetYaxis()->SetTitle("JetEta");
 538   fh2PhiJet     ->GetXaxis()->SetTitle("JetPt");  fh2PhiJet     ->GetYaxis()->SetTitle("JetPhi");
 539   fh2PtJet      ->GetXaxis()->SetTitle("JetPt");  fh2PtJet      ->GetYaxis()->SetTitle("JetPt");
 540   fh1PtJet      ->GetXaxis()->SetTitle("JetPt");  fh1PtJet      ->GetYaxis()->SetTitle("#jets");
 541   fh2NtracksJet ->GetXaxis()->SetTitle("JetPt");  fh2NtracksJet ->GetYaxis()->SetTitle("#tracks");
 542   fProNtracksJet->GetXaxis()->SetTitle("JetPt");  fProNtracksJet->GetYaxis()->SetTitle("AgvNtracks");
 543   fh2EtaTrack   ->GetXaxis()->SetTitle("JetPt");  fh2EtaTrack   ->GetYaxis()->SetTitle("TrackEta");
 544   fh2PhiTrack   ->GetXaxis()->SetTitle("JetPt");  fh2PhiTrack   ->GetYaxis()->SetTitle("TrackPhi");
 545   fh2PtTrack    ->GetXaxis()->SetTitle("JetPt");  fh2PtTrack    ->GetYaxis()->SetTitle("TrackPt");
 546   fh2FF         ->GetXaxis()->SetTitle("JetPt");  fh2FF         ->GetYaxis()->SetTitle("FF");
 547   fh2DelEta     ->GetXaxis()->SetTitle("JetPt");  fh2DelEta     ->GetYaxis()->SetTitle("DelEta");
 548   fh2DelPhi     ->GetXaxis()->SetTitle("JetPt");  fh2DelPhi     ->GetYaxis()->SetTitle("DelPhi");
 549   fh2DelR       ->GetXaxis()->SetTitle("JetPt");  fh2DelR       ->GetYaxis()->SetTitle("DelR");
 550
 551   fh1PtLeadingJet       ->GetXaxis()->SetTitle("JetPt");  fh1PtLeadingJet       ->GetYaxis()->SetTitle("#leading jets");
 552   fh2NtracksLeadingJet  ->GetXaxis()->SetTitle("JetPt");  fh2NtracksLeadingJet  ->GetYaxis()->SetTitle("#tracks leading jet");
 553   fProNtracksLeadingJet ->GetXaxis()->SetTitle("JetPt");  fProNtracksLeadingJet ->GetYaxis()->SetTitle("AvgNtracks leading jet");
 554   fh2DelR80pcNch        ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcNch        ->GetYaxis()->SetTitle("R containing 80% of tracks");
 555   fProDelR80pcNch       ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcNch       ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
 556   fh2DelR80pcPt         ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcPt         ->GetYaxis()->SetTitle("R containing 80% of pT");
 557   fProDelR80pcPt        ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcPt        ->GetYaxis()->SetTitle("<R> containing 80% of pT");
 558   fh2AreaCh             ->GetXaxis()->SetTitle("JetPt");  fh2AreaCh             ->GetYaxis()->SetTitle("Jet area");
 559   fProAreaCh            ->GetXaxis()->SetTitle("JetPt");  fProAreaCh            ->GetYaxis()->SetTitle("<jet area>");
 560   fh3PtDelRNchSum       ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRNchSum       ->GetYaxis()->SetTitle("R");  fh3PtDelRNchSum->GetZaxis()->SetTitle("NchSum");
 561   fh3PtDelRPtSum        ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRPtSum        ->GetYaxis()->SetTitle("R");  fh3PtDelRPtSum ->GetZaxis()->SetTitle("PtSum");
 562   fProDiffJetShape->GetXaxis()->SetTitle("R");
 563   fProDiffJetShape->GetYaxis()->SetTitle("Diff jet shape");
 564   fProIntJetShape->GetXaxis()->SetTitle("R");
 565   fProIntJetShape->GetYaxis()->SetTitle("Integrated jet shape");
 566
 567   fh1PtSumInJetConeUE     ->GetXaxis()->SetTitle("p_{T}^{sum, UE}(in cone R)");  fh1PtSumInJetConeUE     ->GetYaxis()->SetTitle("#leading jets");
 568   fh2NtracksLeadingJetUE  ->GetXaxis()->SetTitle("JetPt");  fh2NtracksLeadingJetUE  ->GetYaxis()->SetTitle("#tracks UE");
 569   fProNtracksLeadingJetUE ->GetXaxis()->SetTitle("JetPt");  fProNtracksLeadingJetUE ->GetYaxis()->SetTitle("AvgNtracks UE");
 570   fh2DelR80pcNchUE        ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcNchUE        ->GetYaxis()->SetTitle("R containing 80% of tracks");
 571   fProDelR80pcNchUE       ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcNchUE       ->GetYaxis()->SetTitle("<R> containing 80% of tracks");
 572   fh2DelR80pcPtUE         ->GetXaxis()->SetTitle("JetPt");  fh2DelR80pcPtUE         ->GetYaxis()->SetTitle("R containing 80% of pT");
 573   fProDelR80pcPtUE        ->GetXaxis()->SetTitle("JetPt");  fProDelR80pcPtUE        ->GetYaxis()->SetTitle("<R> containing 80% of pT");
 574   fh2AreaChUE             ->GetXaxis()->SetTitle("JetPt");  fh2AreaChUE             ->GetYaxis()->SetTitle("UE area");
 575   fProAreaChUE            ->GetXaxis()->SetTitle("JetPt");  fProAreaChUE            ->GetYaxis()->SetTitle("<UE area>");
 576   fh3PtDelRNchSumUE       ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRNchSumUE       ->GetYaxis()->SetTitle("R");  fh3PtDelRNchSumUE->GetZaxis()->SetTitle("NchSumUE");
 577   fh3PtDelRPtSumUE        ->GetXaxis()->SetTitle("JetPt");  fh3PtDelRPtSumUE        ->GetYaxis()->SetTitle("R");  fh3PtDelRPtSumUE ->GetZaxis()->SetTitle("PtSumUE");
 578   fProDiffJetShapeUE->GetXaxis()->SetTitle("R");
 579   fProDiffJetShapeUE->GetYaxis()->SetTitle("Diff jet shape UE");
 580   fProIntJetShapeUE->GetXaxis()->SetTitle("R");
 581   fProIntJetShapeUE->GetYaxis()->SetTitle("Integrated jet shape UE");
 582
 583   fCommonHistList->Add(fh1EvtSelection);
 584   fCommonHistList->Add(fh1VertexNContributors);
 585   fCommonHistList->Add(fh1VertexZ);
 586   fCommonHistList->Add(fh1Xsec);
 587   fCommonHistList->Add(fh1Trials);
 588   fCommonHistList->Add(fh1PtHard);
 589   fCommonHistList->Add(fh1PtHardTrials);
 590   fCommonHistList->Add(fh2EtaJet);
 591   fCommonHistList->Add(fh2PhiJet);
 592   fCommonHistList->Add(fh2PtJet);
 593   fCommonHistList->Add(fh1PtJet);
 594   fCommonHistList->Add(fh2NtracksJet);
 595   fCommonHistList->Add(fProNtracksJet);
 596   fCommonHistList->Add(fh2EtaTrack);
 597   fCommonHistList->Add(fh2PhiTrack);
 598   fCommonHistList->Add(fh2PtTrack);
 599   fCommonHistList->Add(fh2FF);
 600   fCommonHistList->Add(fh2DelEta);
 601   fCommonHistList->Add(fh2DelPhi);
 602   fCommonHistList->Add(fh2DelR);
 603
 604   fCommonHistList->Add(fh1PtLeadingJet);
 605   fCommonHistList->Add(fh2NtracksLeadingJet);
 606   fCommonHistList->Add(fProNtracksLeadingJet);
 607   fCommonHistList->Add(fh2DelR80pcNch);
 608   fCommonHistList->Add(fProDelR80pcNch);
 609   fCommonHistList->Add(fh2DelR80pcPt);
 610   fCommonHistList->Add(fProDelR80pcPt);
 611   fCommonHistList->Add(fh2AreaCh);
 612   fCommonHistList->Add(fProAreaCh);
 613   fCommonHistList->Add(fh3PtDelRNchSum);
 614   fCommonHistList->Add(fh3PtDelRPtSum);
 615   fCommonHistList->Add(fProDiffJetShape);
 616   fCommonHistList->Add(fProIntJetShape);
 617
 618
 619   fCommonHistList->Add(fh1PtSumInJetConeUE);
 620   fCommonHistList->Add(fh2NtracksLeadingJetUE);
 621   fCommonHistList->Add(fProNtracksLeadingJetUE);
 622   fCommonHistList->Add(fh2DelR80pcNchUE);
 623   fCommonHistList->Add(fProDelR80pcNchUE);
 624   fCommonHistList->Add(fh2DelR80pcPtUE);
 625   fCommonHistList->Add(fProDelR80pcPtUE);
 626   fCommonHistList->Add(fh2AreaChUE);
 627   fCommonHistList->Add(fProAreaChUE);
 628   fCommonHistList->Add(fh3PtDelRNchSumUE);
 629   fCommonHistList->Add(fh3PtDelRPtSumUE);
 630   fCommonHistList->Add(fProDiffJetShapeUE);
 631   fCommonHistList->Add(fProIntJetShapeUE);
 632
 633   // =========== Switch on Sumw2 for all histos ===========
 634   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
 635     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
 636     if (h1) h1->Sumw2();
 637     else{
 638       TProfile *hPro = dynamic_cast<TProfile*>(fCommonHistList->At(i));
 639       if(hPro) hPro->Sumw2();
 640     }
 641   }
 642
 643   TH1::AddDirectory(oldStatus);
 644   PostData(1, fCommonHistList);
 645 }
 646 //_________________________________________________________________________________//
 647
 648 void AliAnalysisTaskJetProperties::Init()
 649 {
 650   // Initialization
 651   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::Init()");
 652
 653 }
 654 //_________________________________________________________________________________//
 655
 656 void AliAnalysisTaskJetProperties::UserExec(Option_t *)
 657 {
 658   // Main loop
 659   // Called for each event
 660   if(fDebug > 1) Printf("AliAnalysisTaskJetProperties::UserExec()");
 661
 662   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
 663   // Trigger selection
 664
 665   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 666     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 667   if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
 668     if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
 669       fh1EvtSelection->Fill(1.);
 670       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
 671       PostData(1, fCommonHistList);
 672       return;
 673     }
 674   }
 675
 676   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 677   if(!fESD){
 678     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
 679   }
 680
 681   fMCEvent = MCEvent();
 682   if(!fMCEvent){
 683     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
 684   }
 685
 686   // get AOD event from input/ouput
 687   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 688   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
 689     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
 690     if(fUseAODInputJets) fAODJets = fAOD;
 691     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
 692   }
 693   else {
 694     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 695     if( handler && handler->InheritsFrom("AliAODHandler") ) {
 696       fAOD = ((AliAODHandler*)handler)->GetAOD();
 697       fAODJets = fAOD;
 698       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
 699     }
 700   }
 701
 702   if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
 703     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 704     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
 705       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
 706       if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
 707     }
 708   }
 709
 710   if(fNonStdFile.Length()!=0){
 711     // case we have an AOD extension - fetch the jets from the extended output
 712
 713     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 714     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 715     if(!fAODExtension){
 716       if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
 717     }
 718   }
 719
 720   if(!fAOD){
 721     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
 722     return;
 723   }
 724   if(!fAODJets){
 725     Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
 726     return;
 727   }
 728
 729
 730   // *** vertex cut ***
 731   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 732   Int_t nTracksPrim = primVtx->GetNContributors();
 733   fh1VertexNContributors->Fill(nTracksPrim);
 734
 735   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
 736   //if(!nTracksPrim){
 737   if(nTracksPrim<fNContributors){
 738     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
 739     fh1EvtSelection->Fill(3.);
 740     PostData(1, fCommonHistList);
 741     return;
 742   }
 743
 744   fh1VertexZ->Fill(primVtx->GetZ());
 745
 746   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
 747     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 748     fh1EvtSelection->Fill(4.);
 749     PostData(1, fCommonHistList);
 750     return;
 751   }
 752
 753   TString primVtxName(primVtx->GetName());
 754
 755   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
 756     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
 757     fh1EvtSelection->Fill(5.);
 758     PostData(1, fCommonHistList);
 759     return;
 760   }
 761
 762   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
 763   fh1EvtSelection->Fill(0.);
 764   //___ get MC information __________________________________________________________________
 765   Double_t ptHard = 0.;
 766   Double_t nTrials = 1; // trials for MC trigger weight for real data
 767
 768   if(fMCEvent){
 769     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
 770     if(genHeader){
 771       AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 772       AliGenHijingEventHeader*  hijingGenHeader = 0x0;
 773
 774       if(pythiaGenHeader){
 775         if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
 776         nTrials = pythiaGenHeader->Trials();
 777         ptHard  = pythiaGenHeader->GetPtHard();
 778
 779         fh1PtHard->Fill(ptHard);
 780         fh1PtHardTrials->Fill(ptHard,nTrials);
 781
 782
 783       } else { // no pythia, hijing?
 784
 785         if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
 786
 787         hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 788         if(!hijingGenHeader){
 789           Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
 790         } else {
 791           if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
 792         }
 793       }
 794
 795       fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 796     }
 797   }
 798
 799   //___ fetch jets __________________________________________________________________________
 800   Int_t nJ = GetListOfJets(fJetList);
 801   Int_t nJets = 0;
 802   if(nJ>=0) nJets = fJetList->GetEntries();
 803   if(fDebug>2){
 804     Printf("%s:%d Selected jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
 805     if(nJ != nJets) Printf("%s:%d Mismatch Selected Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nJets);
 806   }
 807   FillJetProperties(fJetList);
 808   FillJetShape(fJetList);
 809   FillJetShapeUE(fJetList);
 810   fJetList->Clear();
 811   //Post output data.
 812   PostData(1, fCommonHistList);
 813 }
 814 //_________________________________________________________________________________//
 815
 816 void AliAnalysisTaskJetProperties::Terminate(Option_t *)
 817 {
 818   // terminated
 819
 820   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::Terminate() \n");
 821 }
 822 //_________________________________________________________________________________//
 823
 824 Int_t AliAnalysisTaskJetProperties::GetListOfJets(TList *list)
 825 {
 826   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
 827   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJets() \n");
 828   // fill list of jets selected according to type
 829   if(!list){
 830     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
 831     return -1;
 832   }
 833
 834   if(fBranchJets.Length()==0){
 835     Printf("%s:%d no jet branch specified", (char*)__FILE__,__LINE__);
 836     if(fDebug>1)fAOD->Print();
 837     return 0;
 838   }
 839   TClonesArray *aodJets = 0;
 840   if(fBranchJets.Length())      aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchJets.Data()));
 841   if(!aodJets)                  aodJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchJets.Data()));
 842   if(fAODExtension&&!aodJets)   aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchJets.Data()));
 843
 844   if(!aodJets){
 845     if(fBranchJets.Length())Printf("%s:%d no jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchJets.Data());
 846     if(fDebug>1)fAOD->Print();
 847     return 0;
 848   }
 849
 850   Int_t nJets = 0;
 851   for(Int_t ijet=0; ijet<aodJets->GetEntries(); ijet++){
 852     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodJets->At(ijet));
 853     if(!tmp) continue;
 854     if( tmp->Pt() < fJetPtCut ) continue;
 855     if( tmp->Eta() < fJetEtaMin || tmp->Eta() > fJetEtaMax)continue;
 856     if(fJetRejectType==kReject1Track && tmp->GetRefTracks()->GetEntriesFast()==1)continue;//rejecting 1track jet if...
 857     list->Add(tmp);
 858     nJets++;
 859   }//ij loop
 860   list->Sort();
 861   return nJets;
 862 }
 863 //_________________________________________________________________________________//
 864
 865 Int_t AliAnalysisTaskJetProperties::GetListOfJetTracks(TList* list, const AliAODJet* jet)
 866 {
 867   //this functionality is motivated by AliAnalysisTaskFragmentationFunction
 868   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::GetListOfJetTracks() \n");
 869
 870   // list of jet tracks from trackrefs
 871   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
 872   Int_t NoOfTracks=0;
 873   for (Int_t itrack=0; itrack<nTracks; itrack++) {
 874     if(fTrackType==kTrackUndef){
 875       if(fDebug>2)Printf("%s:%d unknown track type %d in AOD", (char*)__FILE__,__LINE__,kTrackUndef);
 876       return 0;
 877     }//if
 878     else if(fTrackType == kTrackAOD){
 879       AliAODTrack* trackAod = dynamic_cast<AliAODTrack*>(jet->GetRefTracks()->At(itrack));
 880       if(!trackAod){
 881         AliError("expected ref track not found ");
 882         continue;
 883       }
 884       list->Add(trackAod);
 885       NoOfTracks++;
 886     }//if
 887
 888     else if(fTrackType==kTrackAODMC){
 889       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(jet->GetRefTracks()->At(itrack));
 890       if(!trackmc){
 891         AliError("expected ref trackmc not found ");
 892         continue;
 893       }
 894       list->Add(trackmc);
 895       NoOfTracks++;
 896     }//if
 897
 898     else if (fTrackType==kTrackKine){
 899       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
 900       if(!trackkine){
 901         AliError("expected ref trackkine not found ");
 902         continue;
 903       }
 904       list->Add(trackkine);
 905       NoOfTracks++;
 906     }//if
 907   }//itrack loop
 908   list->Sort();
 909   return NoOfTracks;
 910 }//initial loop
 911 //_________________________________________________________________________________//
 912
 913 void AliAnalysisTaskJetProperties::FillJetProperties(TList *jetlist){
 914   //filling up the histograms jets and tracks inside jet
 915   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetProperties() \n");
 916
 917   for(Int_t iJet=0; iJet < jetlist->GetEntries(); iJet++){
 918     Float_t JetPt;Float_t JetPhi;Float_t JetEta;Float_t JetE;
 919     AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(iJet));
 920     if(!jet)continue;
 921     JetEta = jet->Eta();
 922     JetPhi = jet->Phi();
 923     JetPt  = jet->Pt();
 924     JetE   = jet->E();
 925     fh2EtaJet ->Fill(JetPt,JetEta);
 926     fh2PhiJet ->Fill(JetPt,JetPhi);
 927     fh2PtJet  ->Fill(JetPt,JetPt);
 928     fh1PtJet  ->Fill(JetPt);
 929
 930     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
 931     Int_t nJetTracks = 0;
 932     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
 933     if(fDebug>2){
 934       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
 935       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
 936     }
 937
 938     fh2NtracksJet  ->Fill(JetPt,fTrackListJet->GetEntries());
 939     fProNtracksJet ->Fill(JetPt,fTrackListJet->GetEntries());
 940
 941     for (Int_t j =0; j< fTrackListJet->GetEntries(); j++){
 942       if(fTrackType==kTrackUndef)continue;
 943         Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
 944         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;
 945         if(fTrackType==kTrackAOD){
 946           AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
 947           if(!trackaod)continue;
 948           TrackEta = trackaod->Eta();
 949           TrackPhi = trackaod->Phi();
 950           TrackPt  = trackaod->Pt();
 951         }//if kTrackAOD
 952         else if(fTrackType==kTrackAODMC){
 953           AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
 954           if(!trackmc)continue;
 955           TrackEta = trackmc->Eta();
 956           TrackPhi = trackmc->Phi();
 957           TrackPt  = trackmc->Pt();
 958         }//if kTrackAODMC
 959         else if(fTrackType==kTrackKine){
 960           AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
 961           if(!trackkine)continue;
 962           TrackEta = trackkine->Eta();
 963           TrackPhi = trackkine->Phi();
 964           TrackPt  = trackkine->Pt();
 965         }//kTrackKine
 966         if(JetPt)FF = TrackPt/JetPt;
 967         DelEta      = TMath::Abs(JetEta - TrackEta);
 968         DelPhi      = TMath::Abs(JetPhi - TrackPhi);
 969         if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
 970         DelR        = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
 971         AreaJ       = TMath::Pi()*DelR*DelR;
 972         fh2EtaTrack ->Fill(JetPt,TrackEta);
 973         fh2PhiTrack ->Fill(JetPt,TrackPhi);
 974         fh2PtTrack  ->Fill(JetPt,TrackPt);
 975         fh2FF       ->Fill(JetPt,FF);
 976         fh2DelEta   ->Fill(JetPt,DelEta);
 977         fh2DelPhi   ->Fill(JetPt,DelPhi);
 978         fh2DelR     ->Fill(JetPt,DelR);
 979     }//track loop
 980     fTrackListJet->Clear();
 981   }//iJet loop
 982 }//FillJetProperties
 983 //_________________________________________________________________________________//
 984
 985 void AliAnalysisTaskJetProperties::FillJetShape(TList *jetlist){
 986   //filling up the histograms
 987   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
 988   Float_t JetEta; Float_t JetPhi; Float_t JetPt;
 989   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
 990   if(jet){
 991     JetEta = jet->Eta();
 992     JetPhi = jet->Phi();
 993     JetPt  = jet->Pt();
 994     fh1PtLeadingJet->Fill(JetPt);
 995     Float_t NchSumA[50]         = {0.};
 996     Float_t PtSumA[50]          = {0.};
 997     Float_t delRPtSum80pc       = 0;
 998     Float_t delRNtrkSum80pc     = 0;
 999     Float_t PtSumDiffShape[10]  = {0.0};
1000     Float_t PtSumIntShape[10]   = {0.0};
1001     Int_t kNbinsR               = 10;
1002
1003     Int_t nJT = GetListOfJetTracks(fTrackListJet,jet);
1004     Int_t nJetTracks = 0;
1005     if(nJT>=0) nJetTracks = fTrackListJet->GetEntries();
1006     if(fDebug>2){
1007       Printf("%s:%d Jet tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1008       if(nJT != nJetTracks) Printf("%s:%d Mismatch Jet Tracks: %d %d",(char*)__FILE__,__LINE__,nJT,nJetTracks);
1009     }
1010     fh2NtracksLeadingJet->Fill(JetPt,nJetTracks);
1011     fProNtracksLeadingJet->Fill(JetPt,nJetTracks);
1012     Int_t   *index     = new Int_t   [nJetTracks];//dynamic array containing index
1013     Float_t *delRA     = new Float_t [nJetTracks];//dynamic array containing delR
1014     Float_t *delEtaA   = new Float_t [nJetTracks];//dynamic array containing delEta
1015     Float_t *delPhiA   = new Float_t [nJetTracks];//dynamic array containing delPhi
1016     Float_t *trackPtA  = new Float_t [nJetTracks];//dynamic array containing pt-track
1017     Float_t *trackEtaA = new Float_t [nJetTracks];//dynamic array containing eta-track
1018     Float_t *trackPhiA = new Float_t [nJetTracks];//dynamic array containing phi-track
1019     for(Int_t ii=0; ii<nJetTracks; ii++){
1020       index[ii]     = 0;
1021       delRA[ii]     = 0.;
1022       delEtaA[ii]   = 0.;
1023       delPhiA[ii]   = 0.;
1024       trackPtA[ii]  = 0.;
1025       trackEtaA[ii] = 0.;
1026       trackPhiA[ii] = 0.;
1027     }//ii loop
1028
1029     for (Int_t j =0; j< nJetTracks; j++){
1030       if(fTrackType==kTrackUndef)continue;
1031
1032       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1033       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1034
1035       if(fTrackType==kTrackAOD){
1036         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListJet->At(j));
1037         if(!trackaod)continue;
1038         TrackEta = trackaod->Eta();
1039         TrackPhi = trackaod->Phi();
1040         TrackPt  = trackaod->Pt();
1041       }//if kTrackAOD
1042       else if(fTrackType==kTrackAODMC){
1043         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListJet->At(j));
1044         if(!trackmc)continue;
1045         TrackEta = trackmc->Eta();
1046         TrackPhi = trackmc->Phi();
1047         TrackPt  = trackmc->Pt();
1048       }//if kTrackAODMC
1049       else if(fTrackType==kTrackKine){
1050         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListJet->At(j));
1051         if(!trackkine)continue;
1052         TrackEta = trackkine->Eta();
1053         TrackPhi = trackkine->Phi();
1054         TrackPt  = trackkine->Pt();
1055       }//if kTrackKine
1056
1057       DelEta = TMath::Abs(JetEta - TrackEta);
1058       DelPhi = TMath::Abs(JetPhi - TrackPhi);
1059       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1060       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1061       AreaJ  = TMath::Pi()*DelR*DelR;
1062
1063       fh2AreaCh ->Fill(JetPt,AreaJ);
1064       fProAreaCh->Fill(JetPt,AreaJ);
1065       delRA[j]     = DelR;
1066       delEtaA[j]   = DelEta;
1067       delPhiA[j]   = DelPhi;
1068       trackPtA[j]  = TrackPt;
1069       trackEtaA[j] = TrackEta;
1070       trackPhiA[j] = TrackPhi;
1071
1072       //calculating diff and integrated jet shapes
1073       Float_t kDeltaR = 0.1;
1074       Float_t RMin    = kDeltaR/2.0;
1075       Float_t RMax    = kDeltaR/2.0;
1076       Float_t tmpR    = 0.05;
1077       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1078         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1079         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1080         tmpR += 0.1;
1081       }//ii1 loop
1082
1083       for(Int_t ibin=1; ibin<=50; ibin++){
1084         Float_t xlow = 0.02*(ibin-1);
1085         Float_t xup  = 0.02*ibin;
1086         if( xlow <= DelR && DelR < xup){
1087           NchSumA[ibin-1]++;
1088           PtSumA[ibin-1]+= TrackPt;
1089         }//if loop
1090       }//for ibin loop
1091     }//track loop
1092     fTrackListJet->Clear();
1093
1094     //---------------------
1095     Float_t tmp1R = 0.05;
1096     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1097       if(JetPt>20 && JetPt<=100){
1098         fProDiffJetShape->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1099         fProIntJetShape ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1100       }
1101       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1102       for(Int_t k=0; k<13; k++){
1103         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1104         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1105         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1106         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1107         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1108         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1109         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1110         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1111         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1112         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1113         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1114         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1115         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1116         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1117           fProDiffJetShapeA[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1118           fProIntJetShapeA[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1119         }//if
1120       }//k loop
1121       tmp1R +=0.1;
1122     }//jj1 loop
1123     //----------------------//
1124     Float_t PtSum = 0;
1125     Int_t NtrkSum = 0;
1126     Bool_t iflagPtSum   = kFALSE;
1127     Bool_t iflagNtrkSum = kFALSE;
1128     TMath::Sort(nJetTracks,delRA,index,0);
1129     for(Int_t ii=0; ii<nJetTracks; ii++){
1130       NtrkSum ++;
1131       PtSum += trackPtA[index[ii]];
1132       /*
1133         cout << index[ii] << "\t" <<
1134         delR[ii] << "\t" <<
1135         delEta[ii]<< "\t" <<
1136         delPhi[ii]<< "\t" <<
1137         trackPt[ii]<< "\t" <<
1138         trackEta[ii]<< "\t" <<
1139         trackPhi[ii]<< "\t DelR " <<
1140         delR[index[ii]] << endl;
1141       */
1142       if(!iflagNtrkSum){
1143         if((Float_t)NtrkSum/(Float_t)nJetTracks > 0.79){
1144           delRNtrkSum80pc = delRA[index[ii]];
1145           iflagNtrkSum = kTRUE;
1146         }//if NtrkSum
1147       }//if iflag
1148       if(!iflagPtSum){
1149         if(PtSum/JetPt >= 0.8000){
1150           delRPtSum80pc = delRA[index[ii]];
1151           iflagPtSum = kTRUE;
1152         }//if PtSum
1153       }//if iflag
1154     }//track loop 2nd
1155     delete [] index;
1156     delete [] delRA;
1157     delete [] delEtaA;
1158     delete [] delPhiA;
1159     delete [] trackPtA;
1160     delete [] trackEtaA;
1161     delete [] trackPhiA;
1162     fh2DelR80pcNch ->Fill(JetPt,delRNtrkSum80pc);
1163     fProDelR80pcNch->Fill(JetPt,delRNtrkSum80pc);
1164     fh2DelR80pcPt  ->Fill(JetPt,delRPtSum80pc);
1165     fProDelR80pcPt ->Fill(JetPt,delRPtSum80pc);
1166     for(Int_t ibin=0; ibin<50; ibin++){
1167       Float_t iR = 0.02*ibin + 0.01;
1168       fh3PtDelRNchSum->Fill(JetPt,iR,NchSumA[ibin]);
1169       fh3PtDelRPtSum ->Fill(JetPt,iR,PtSumA[ibin]);
1170       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1171       for(Int_t k=0; k<13; k++){
1172         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1173         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1174         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1175         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1176         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1177         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1178         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1179         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1180         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1181         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1182         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1183         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1184         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1185         if(JetPt>jetPtMin && JetPt<jetPtMax){
1186           fProDelRNchSum[k]->Fill(iR,NchSumA[ibin]);
1187           fProDelRPtSum[k]->Fill(iR,PtSumA[ibin]);
1188         }//if
1189       }//k loop
1190     }//ibin loop
1191   }//if(jet)
1192 }//FillJetShape()
1193 //_________________________________________________________________________________//
1194
1195 void AliAnalysisTaskJetProperties::FillJetShapeUE(TList *jetlist){
1196   //filling up the histograms
1197   if(fDebug > 1) printf("AliAnalysisTaskJetProperties::FillJetShape() \n");
1198   AliAODJet *jet = dynamic_cast<AliAODJet*>(jetlist->At(0));//Leading jet only
1199   if(jet){
1200     Double_t jetMom[3];
1201     jet->PxPyPz(jetMom);
1202
1203     TVector3 jet3mom(jetMom);
1204     // Rotate phi and keep eta unchanged
1205
1206     const Double_t alpha = TMath::Pi()/2.;
1207     Double_t etaTilted = jet3mom.Eta();
1208     Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
1209     if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1210
1211     Double_t JetPt  = jet->Pt();
1212
1213     Float_t NchSumA[50]         = {0.};
1214     Float_t PtSumA[50]          = {0.};
1215     Float_t delRPtSum80pc       = 0;
1216     Float_t delRNtrkSum80pc     = 0;
1217     Float_t PtSumDiffShape[10]  = {0.0};
1218     Float_t PtSumIntShape[10]   = {0.0};
1219     Int_t kNbinsR               = 10;
1220
1221     //Int_t nTracks = GetListOfTracks(fTrackList,fTrackType);
1222     GetListOfTracks(fTrackList,fTrackType);
1223     Double_t sumPtPerp = 0;
1224     if(fDebug > 100) printf("Cone radius for bckg. = %f Track type = %d\n",fJetRadius, fTrackType);
1225     GetTracksTiltedwrpJetAxis(alpha, fTrackList, fTrackListUE,jet,fJetRadius,sumPtPerp);
1226     fTrackList->Clear();
1227     fh1PtSumInJetConeUE->Fill(sumPtPerp);
1228
1229     Int_t nTracksUE = fTrackListUE->GetEntries();
1230     fh2NtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1231     fProNtracksLeadingJetUE->Fill(JetPt,nTracksUE);
1232
1233     Int_t   *index     = new Int_t   [nTracksUE];//dynamic array containing index
1234     Float_t *delRA     = new Float_t [nTracksUE];//dynamic array containing delR
1235     Float_t *delEtaA   = new Float_t [nTracksUE];//dynamic array containing delEta
1236     Float_t *delPhiA   = new Float_t [nTracksUE];//dynamic array containing delPhi
1237     Float_t *trackPtA  = new Float_t [nTracksUE];//dynamic array containing pt-track
1238     Float_t *trackEtaA = new Float_t [nTracksUE];//dynamic array containing eta-track
1239     Float_t *trackPhiA = new Float_t [nTracksUE];//dynamic array containing phi-track
1240     for(Int_t ii=0; ii<nTracksUE; ii++){
1241       index[ii]     = 0;
1242       delRA[ii]     = 0.;
1243       delEtaA[ii]   = 0.;
1244       delPhiA[ii]   = 0.;
1245       trackPtA[ii]  = 0.;
1246       trackEtaA[ii] = 0.;
1247       trackPhiA[ii] = 0.;
1248     }//ii loop
1249
1250     for (Int_t j =0; j< nTracksUE; j++){
1251       if(fTrackType==kTrackUndef)continue;
1252
1253       Float_t TrackEta=-99.0; Float_t TrackPt=-99.0; Float_t TrackPhi=-99.0;
1254       Float_t DelEta=-99.0;  Float_t DelPhi=-99.0; Float_t DelR=-99.0; Float_t AreaJ=-99.0;
1255
1256       if(fTrackType==kTrackAOD){
1257         AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(fTrackListUE->At(j));
1258         if(!trackaod)continue;
1259         TrackEta = trackaod->Eta();
1260         TrackPhi = trackaod->Phi();
1261         TrackPt  = trackaod->Pt();
1262         if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1263       }//if kTrackAOD
1264       else if(fTrackType==kTrackAODMC){
1265         AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(fTrackListUE->At(j));
1266         if(!trackmc)continue;
1267         TrackEta = trackmc->Eta();
1268         TrackPhi = trackmc->Phi();
1269         TrackPt  = trackmc->Pt();
1270         if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1271       }//if kTrackAODMC
1272       else if(fTrackType==kTrackKine){
1273         AliVParticle* trackkine = dynamic_cast<AliVParticle*>(fTrackListUE->At(j));
1274         if(!trackkine)continue;
1275         TrackEta = trackkine->Eta();
1276         TrackPhi = trackkine->Phi();
1277         TrackPt  = trackkine->Pt();
1278         if(fDebug > 100) printf("FillJetShapeUE itrack, trackPt %d,  %f\n",j,TrackPt);
1279       }//if kTrackKine
1280
1281       DelEta = TMath::Abs(etaTilted - TrackEta);
1282       DelPhi = TMath::Abs(phiTilted - TrackPhi);
1283       if(DelPhi>TMath::Pi())DelPhi = TMath::Abs(DelPhi-TMath::TwoPi());
1284       DelR   = TMath::Sqrt(DelEta*DelEta + DelPhi*DelPhi);
1285       AreaJ  = TMath::Pi()*DelR*DelR;
1286
1287       fh2AreaChUE ->Fill(JetPt,AreaJ);
1288       fProAreaChUE->Fill(JetPt,AreaJ);
1289
1290       delRA[j]     = DelR;
1291       delEtaA[j]   = DelEta;
1292       delPhiA[j]   = DelPhi;
1293       trackPtA[j]  = TrackPt;
1294       trackEtaA[j] = TrackEta;
1295       trackPhiA[j] = TrackPhi;
1296
1297       //calculating diff and integrated jet shapes
1298       Float_t kDeltaR = 0.1;
1299       Float_t RMin    = kDeltaR/2.0;
1300       Float_t RMax    = kDeltaR/2.0;
1301       Float_t tmpR    = 0.05;
1302       for(Int_t ii1=0; ii1<kNbinsR;ii1++){
1303         if((DelR > (tmpR-RMin)) && (DelR <=(tmpR+RMax)))PtSumDiffShape[ii1]+= TrackPt;
1304         if(DelR>0.0 && DelR <=(tmpR+RMax))PtSumIntShape[ii1]+= TrackPt;
1305         tmpR += 0.1;
1306       }//ii1 loop
1307
1308       for(Int_t ibin=1; ibin<=50; ibin++){
1309         Float_t xlow = 0.02*(ibin-1);
1310         Float_t xup  = 0.02*ibin;
1311         if( xlow <= DelR && DelR < xup){
1312           NchSumA[ibin-1]++;
1313           PtSumA[ibin-1]+= TrackPt;
1314         }//if loop
1315       }//for ibin loop
1316     }//track loop
1317     fTrackListUE->Clear();
1318
1319     //---------------------
1320     Float_t tmp1R = 0.05;
1321     for(Int_t jj1=0; jj1<kNbinsR;jj1++){
1322       if(JetPt>20 && JetPt<=100){
1323         fProDiffJetShapeUE->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1324         fProIntJetShapeUE ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1325       }
1326       Float_t jetPtMin0=20.0; Float_t jetPtMax0=30.0;
1327       for(Int_t k=0; k<13; k++){
1328         if(k==0) {jetPtMin0=20.0;jetPtMax0=25.0;}
1329         if(k==1) {jetPtMin0=25.0;jetPtMax0=30.0;}
1330         if(k==2) {jetPtMin0=20.0;jetPtMax0=30.0;}
1331         if(k==3) {jetPtMin0=30.0;jetPtMax0=35.0;}
1332         if(k==4) {jetPtMin0=35.0;jetPtMax0=40.0;}
1333         if(k==5) {jetPtMin0=30.0;jetPtMax0=40.0;}
1334         if(k==6) {jetPtMin0=40.0;jetPtMax0=50.0;}
1335         if(k==7) {jetPtMin0=50.0;jetPtMax0=60.0;}
1336         if(k==8) {jetPtMin0=40.0;jetPtMax0=60.0;}
1337         if(k==9) {jetPtMin0=60.0;jetPtMax0=70.0;}
1338         if(k==10){jetPtMin0=70.0;jetPtMax0=80.0;}
1339         if(k==11){jetPtMin0=60.0;jetPtMax0=80.0;}
1340         if(k==12){jetPtMin0=80.0;jetPtMax0=100.0;}
1341         if(JetPt>jetPtMin0 && JetPt<=jetPtMax0){
1342           fProDiffJetShapeAUE[k]->Fill(tmp1R,PtSumDiffShape[jj1]/JetPt);
1343           fProIntJetShapeAUE[k] ->Fill(tmp1R,PtSumIntShape[jj1]/JetPt);
1344         }//if
1345       }//k loop
1346       tmp1R +=0.1;
1347     }//jj1 loop
1348     //----------------------//
1349     Float_t PtSum = 0;
1350     Int_t NtrkSum = 0;
1351     Bool_t iflagPtSum   = kFALSE;
1352     Bool_t iflagNtrkSum = kFALSE;
1353     TMath::Sort(nTracksUE,delRA,index,0);
1354     for(Int_t ii=0; ii<nTracksUE; ii++){
1355       NtrkSum ++;
1356       PtSum += trackPtA[index[ii]];
1357       /*
1358         cout << index[ii] << "\t" <<
1359         delR[ii] << "\t" <<
1360         delEta[ii]<< "\t" <<
1361         delPhi[ii]<< "\t" <<
1362         trackPt[ii]<< "\t" <<
1363         trackEta[ii]<< "\t" <<
1364         trackPhi[ii]<< "\t DelR " <<
1365         delR[index[ii]] << endl;
1366       */
1367       if(!iflagNtrkSum){
1368         if((Float_t)NtrkSum/(Float_t)nTracksUE > 0.79){
1369           delRNtrkSum80pc = delRA[index[ii]];
1370           iflagNtrkSum = kTRUE;
1371         }//if NtrkSum
1372       }//if iflag
1373       if(!iflagPtSum){
1374         if(PtSum/JetPt >= 0.8000){
1375           delRPtSum80pc = delRA[index[ii]];
1376           iflagPtSum = kTRUE;
1377         }//if PtSum
1378       }//if iflag
1379     }//track loop 2nd
1380     delete [] index;
1381     delete [] delRA;
1382     delete [] delEtaA;
1383     delete [] delPhiA;
1384     delete [] trackPtA;
1385     delete [] trackEtaA;
1386     delete [] trackPhiA;
1387     fh2DelR80pcNchUE ->Fill(JetPt,delRNtrkSum80pc);
1388     fProDelR80pcNchUE->Fill(JetPt,delRNtrkSum80pc);
1389     fh2DelR80pcPtUE  ->Fill(JetPt,delRPtSum80pc);
1390     fProDelR80pcPtUE ->Fill(JetPt,delRPtSum80pc);
1391     for(Int_t ibin=0; ibin<50; ibin++){
1392       Float_t iR = 0.02*ibin + 0.01;
1393       fh3PtDelRNchSumUE->Fill(JetPt,iR,NchSumA[ibin]);
1394       fh3PtDelRPtSumUE ->Fill(JetPt,iR,PtSumA[ibin]);
1395       Float_t jetPtMin=20.0; Float_t jetPtMax=30.0;
1396       for(Int_t k=0; k<13; k++){
1397         if(k==0) {jetPtMin=20.0;jetPtMax=25.0;}
1398         if(k==1) {jetPtMin=25.0;jetPtMax=30.0;}
1399         if(k==2) {jetPtMin=20.0;jetPtMax=30.0;}
1400         if(k==3) {jetPtMin=30.0;jetPtMax=35.0;}
1401         if(k==4) {jetPtMin=35.0;jetPtMax=40.0;}
1402         if(k==5) {jetPtMin=30.0;jetPtMax=40.0;}
1403         if(k==6) {jetPtMin=40.0;jetPtMax=50.0;}
1404         if(k==7) {jetPtMin=50.0;jetPtMax=60.0;}
1405         if(k==8) {jetPtMin=40.0;jetPtMax=60.0;}
1406         if(k==9) {jetPtMin=60.0;jetPtMax=70.0;}
1407         if(k==10){jetPtMin=70.0;jetPtMax=80.0;}
1408         if(k==11){jetPtMin=60.0;jetPtMax=80.0;}
1409         if(k==12){jetPtMin=80.0;jetPtMax=100.0;}
1410         if(JetPt>jetPtMin && JetPt<jetPtMax){
1411           fProDelRNchSumUE[k]->Fill(iR,NchSumA[ibin]);
1412           fProDelRPtSumUE[k]->Fill(iR,PtSumA[ibin]);
1413         }//if
1414       }//k loop
1415     }//ibin loop
1416   }//if(jet)
1417 }//FillJetShapeUE()
1418 //_________________________________________________________________________________//
1419
1420 void AliAnalysisTaskJetProperties::GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt){
1421   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx
1422   // List of tracks in cone perpendicular to the jet azimuthal direction
1423   Double_t jetMom[3];
1424   jet->PxPyPz(jetMom);
1425
1426   TVector3 jet3mom(jetMom);
1427   // Rotate phi and keep eta unchanged
1428   Double_t etaTilted = jet3mom.Eta();//no change in eta
1429   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;//rotate phi by alpha
1430   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
1431
1432   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
1433     if(fTrackType==kTrackUndef)continue;
1434     Double_t trackMom[3];
1435
1436     if(fTrackType==kTrackAOD){
1437       AliAODTrack *trackaod = dynamic_cast<AliAODTrack*>(inputlist->At(itrack));
1438       if(!trackaod)continue;
1439       trackaod->PxPyPz(trackMom);
1440       TVector3 track3mom(trackMom);
1441
1442       Double_t deta = track3mom.Eta() - etaTilted;
1443       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1444       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1445       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1446
1447       if(dR<=radius){
1448         outputlist->Add(trackaod);
1449         sumPt += trackaod->Pt();
1450         if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1451       }//if dR<=radius
1452     }//if kTrackAOD
1453     else if(fTrackType==kTrackAODMC){
1454       AliAODMCParticle* trackmc = dynamic_cast<AliAODMCParticle*>(inputlist->At(itrack));
1455       if(!trackmc)continue;
1456       trackmc->PxPyPz(trackMom);
1457       TVector3 track3mom(trackMom);
1458
1459       Double_t deta = track3mom.Eta() - etaTilted;
1460       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1461       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1462       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1463
1464       if(dR<=radius){
1465         outputlist->Add(trackmc);
1466         sumPt += trackmc->Pt();
1467         if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1468       }//if dR<=radius
1469     }//if kTrackAODMC
1470     else if(fTrackType==kTrackKine){
1471       AliVParticle* trackkine = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
1472       if(!trackkine)continue;
1473       trackkine->PxPyPz(trackMom);
1474       TVector3 track3mom(trackMom);
1475
1476       Double_t deta = track3mom.Eta() - etaTilted;
1477       Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
1478       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
1479       Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
1480
1481       if(dR<=radius){
1482         outputlist->Add(trackkine);
1483         sumPt += trackkine->Pt();
1484         if(fDebug > 100) printf("GetTracksTiltewrpJetAxis itrack, trackPt %d,  %f\n",itrack,track3mom.Pt());
1485       }//if dR<=radius
1486     }//kTrackKine
1487   }//itrack
1488 }//initial loop
1489 //-----------------------------------------------//
1490 Int_t AliAnalysisTaskJetProperties::GetListOfTracks(TList *list, Int_t type)
1491 {
1492   //This part  is inherited from AliAnalysisTaskFragmentationFunction.cxx (and modified)
1493   // fill list of tracks selected according to type
1494
1495   if(fDebug > 2) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
1496
1497   if(!list){
1498     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1499     return -1;
1500   }
1501
1502   if(!fAOD) return -1;
1503
1504   if(!fAOD->GetTracks()) return 0;
1505
1506   if(type==kTrackUndef) return 0;
1507
1508   Int_t iCount = 0;
1509   if(type==kTrackAOD){
1510     // all rec. tracks, esd filter mask, within acceptance
1511     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
1512       AliAODTrack *tr = fAOD->GetTrack(it);
1513       if(!tr)continue;
1514       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;//selecting filtermask
1515       if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
1516       if(tr->Pt()  < fTrackPtCut) continue;
1517       list->Add(tr);
1518       iCount++;
1519     }//track loop
1520   }//if type
1521   else if (type==kTrackKine){
1522     // kine particles, primaries, charged only within acceptance
1523     if(!fMCEvent) return iCount;
1524     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
1525       if(!fMCEvent->IsPhysicalPrimary(it))continue;//selecting only primaries
1526       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
1527       if(!part)continue;
1528       if(part->Particle()->GetPDG()->Charge()==0)continue;//removing neutrals
1529       if(part->Eta() < fTrackEtaMin || part->Eta() > fTrackEtaMax) continue;//eta cut
1530       if(part->Pt()  < fTrackPtCut)continue;//pt cut
1531       list->Add(part);
1532       iCount++;
1533     }//track loop
1534   }//if type
1535   else if (type==kTrackAODMC) {
1536     // MC particles (from AOD), physical primaries, charged only within acceptance
1537     if(!fAOD) return -1;
1538
1539     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1540     if(!tca)return iCount;
1541
1542     for(int it=0; it<tca->GetEntriesFast(); ++it){
1543       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1544       if(!part)continue;
1545       if(!part->IsPhysicalPrimary())continue;
1546       if(part->Charge()==0) continue;
1547       if(part->Eta() > fTrackEtaMax || part->Eta() < fTrackEtaMin)continue;
1548       if(part->Pt()  < fTrackPtCut) continue;
1549       list->Add(part);
1550       iCount++;
1551     }//track loop
1552   }//if type
1553
1554   list->Sort();
1555   return iCount;
1556 }
1557 //_______________________________________________________________________________