PWGJE/AliAnalysisTaskFragmentationFunction.cxx

   1 // *************************************************************************
   2 // *                                                                       *
   3 // * Task for Fragmentation Function Analysis in PWG4 Jet Task Force Train *
   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 #include "TAxis.h"
  36
  37 #include "AliAODInputHandler.h"
  38 #include "AliAODHandler.h"
  39 #include "AliESDEvent.h"
  40 #include "AliAODMCParticle.h"
  41 #include "AliAODJet.h"
  42 #include "AliAODJetEventBackground.h"
  43 #include "AliGenPythiaEventHeader.h"
  44 #include "AliGenHijingEventHeader.h"
  45 #include "AliInputEventHandler.h"
  46
  47 #include "AliAnalysisHelperJetTasks.h"
  48 #include "AliAnalysisManager.h"
  49 #include "AliAnalysisTaskSE.h"
  50 #include "AliVParticle.h"
  51 #include "AliVEvent.h"
  52
  53 #include "AliAnalysisTaskFragmentationFunction.h"
  54 using std::cout;
  55 using std::endl;
  56 using std::cerr;
  57
  58 ClassImp(AliAnalysisTaskFragmentationFunction)
  59
  60 //____________________________________________________________________________
  61 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
  62    : AliAnalysisTaskSE()
  63    ,fESD(0)
  64    ,fAOD(0)
  65    ,fAODJets(0)
  66    ,fAODExtension(0)
  67    ,fNonStdFile("")
  68    ,fBranchRecJets("jets")
  69    ,fBranchRecBckgClusters("")
  70    ,fBranchGenJets("")
  71    ,fBranchEmbeddedJets("")
  72    ,fTrackTypeGen(0)
  73    ,fJetTypeGen(0)
  74    ,fJetTypeRecEff(0)
  75    ,fUseAODInputJets(kTRUE)
  76    ,fFilterMask(0)
  77    ,fUsePhysicsSelection(kTRUE)
  78    ,fEvtSelectionMask(0)
  79    ,fEventClass(0)
  80    ,fMaxVertexZ(10)
  81    ,fTrackPtCut(0)
  82    ,fTrackEtaMin(0)
  83    ,fTrackEtaMax(0)
  84    ,fTrackPhiMin(0)
  85    ,fTrackPhiMax(0)
  86    ,fUseExtraTracks(0)
  87    ,fUseExtraTracksBgr(0)
  88    ,fCutFractionPtEmbedded(0)
  89    ,fUseEmbeddedJetAxis(0)
  90    ,fUseEmbeddedJetPt(0)
  91    ,fJetPtCut(0)
  92    ,fJetEtaMin(0)
  93    ,fJetEtaMax(0)
  94    ,fJetPhiMin(0)
  95    ,fJetPhiMax(0)
  96    ,fFFRadius(0)
  97    ,fFFMinLTrackPt(-1)
  98    ,fFFMaxTrackPt(-1)
  99    ,fFFMinnTracks(0)
 100    ,fFFBckgRadius(0)
 101    ,fBckgMode(0)
 102    ,fQAMode(0)
 103    ,fFFMode(0)
 104    ,fEffMode(0)
 105    ,fJSMode(0)
 106    ,fAvgTrials(0)
 107    ,fTracksRecCuts(0)
 108    ,fTracksGen(0)
 109    ,fTracksAODMCCharged(0)
 110    ,fTracksAODMCChargedSecNS(0)
 111    ,fTracksAODMCChargedSecS(0)
 112    ,fTracksRecQualityCuts(0)
 113    ,fJetsRec(0)
 114    ,fJetsRecCuts(0)
 115    ,fJetsGen(0)
 116    ,fJetsRecEff(0)
 117    ,fJetsEmbedded(0)
 118    ,fBckgJetsRec(0)
 119    ,fBckgJetsRecCuts(0)
 120    ,fBckgJetsGen(0)
 121    ,fQATrackHistosRecCuts(0)
 122    ,fQATrackHistosGen(0)
 123    ,fQAJetHistosRec(0)
 124    ,fQAJetHistosRecCuts(0)
 125    ,fQAJetHistosRecCutsLeading(0)
 126    ,fQAJetHistosGen(0)
 127    ,fQAJetHistosGenLeading(0)
 128    ,fQAJetHistosRecEffLeading(0)
 129    ,fFFHistosRecCuts(0)
 130    ,fFFHistosRecCutsInc(0)
 131    ,fFFHistosRecLeadingTrack(0)
 132    ,fFFHistosGen(0)
 133    ,fFFHistosGenInc(0)
 134    ,fFFHistosGenLeadingTrack(0)
 135    ,fQATrackHighPtThreshold(0)
 136    ,fFFNBinsJetPt(0)
 137    ,fFFJetPtMin(0)
 138    ,fFFJetPtMax(0)
 139    ,fFFNBinsPt(0)
 140    ,fFFPtMin(0)
 141    ,fFFPtMax(0)
 142    ,fFFNBinsXi(0)
 143    ,fFFXiMin(0)
 144    ,fFFXiMax(0)
 145    ,fFFNBinsZ(0)
 146    ,fFFZMin(0)
 147    ,fFFZMax(0)
 148    ,fQAJetNBinsPt(0)
 149    ,fQAJetPtMin(0)
 150    ,fQAJetPtMax(0)
 151    ,fQAJetNBinsEta(0)
 152    ,fQAJetEtaMin(0)
 153    ,fQAJetEtaMax(0)
 154    ,fQAJetNBinsPhi(0)
 155    ,fQAJetPhiMin(0)
 156    ,fQAJetPhiMax(0)
 157    ,fQATrackNBinsPt(0)
 158    ,fQATrackPtMin(0)
 159    ,fQATrackPtMax(0)
 160    ,fQATrackNBinsEta(0)
 161    ,fQATrackEtaMin(0)
 162    ,fQATrackEtaMax(0)
 163    ,fQATrackNBinsPhi(0)
 164    ,fQATrackPhiMin(0)
 165    ,fQATrackPhiMax(0)
 166    ,fCommonHistList(0)
 167    ,fh1EvtSelection(0)
 168    ,fh1VertexNContributors(0)
 169    ,fh1VertexZ(0)
 170    ,fh1EvtMult(0)
 171    ,fh1EvtCent(0)
 172    ,fh1Xsec(0)
 173    ,fh1Trials(0)
 174    ,fh1PtHard(0)
 175    ,fh1PtHardTrials(0)
 176    ,fh1nRecJetsCuts(0)
 177    ,fh1nGenJets(0)
 178    ,fh1nRecEffJets(0)
 179    ,fh1nEmbeddedJets(0)
 180    ,fh1nRecBckgJetsCuts(0)
 181    ,fh1nGenBckgJets(0)
 182    ,fh2PtRecVsGenPrim(0)
 183    ,fh2PtRecVsGenSec(0)
 184    ,fQATrackHistosRecEffGen(0)
 185    ,fQATrackHistosRecEffRec(0)
 186    ,fQATrackHistosSecRecNS(0)
 187    ,fQATrackHistosSecRecS(0)
 188    ,fQATrackHistosSecRecSsc(0)
 189    ,fFFHistosRecEffRec(0)
 190    ,fFFHistosSecRecNS(0)
 191    ,fFFHistosSecRecS(0)
 192    ,fFFHistosSecRecSsc(0)
 193    // Background
 194    ,fh1BckgMult0(0)
 195    ,fh1BckgMult1(0)
 196    ,fh1BckgMult2(0)
 197    ,fh1BckgMult3(0)
 198    ,fh1BckgMult4(0)
 199    ,fh1FractionPtEmbedded(0)
 200    ,fh1IndexEmbedded(0)
 201    ,fh2DeltaPtVsJetPtEmbedded(0)
 202    ,fh2DeltaPtVsRecJetPtEmbedded(0)
 203    ,fh1DeltaREmbedded(0)
 204    ,fQABckgHisto0RecCuts(0)
 205    ,fQABckgHisto0Gen(0)
 206    ,fQABckgHisto1RecCuts(0)
 207    ,fQABckgHisto1Gen(0)
 208    ,fQABckgHisto2RecCuts(0)
 209    ,fQABckgHisto2Gen(0)
 210    ,fQABckgHisto3RecCuts(0)
 211    ,fQABckgHisto3Gen(0)
 212    ,fQABckgHisto4RecCuts(0)
 213    ,fQABckgHisto4Gen(0)
 214    ,fFFBckgHisto0RecCuts(0)
 215    ,fFFBckgHisto0Gen(0)
 216    ,fFFBckgHisto1RecCuts(0)
 217    ,fFFBckgHisto1Gen(0)
 218    ,fFFBckgHisto2RecCuts(0)
 219    ,fFFBckgHisto2Gen(0)
 220    ,fFFBckgHisto3RecCuts(0)
 221    ,fFFBckgHisto3Gen(0)
 222    ,fFFBckgHisto4RecCuts(0)
 223    ,fFFBckgHisto4Gen(0)
 224    ,fFFBckgHisto0RecEffRec(0)
 225    ,fFFBckgHisto0SecRecNS(0)
 226    ,fFFBckgHisto0SecRecS(0)
 227    ,fFFBckgHisto0SecRecSsc(0)
 228     // jet shape
 229    ,fProNtracksLeadingJet(0)
 230    ,fProDelR80pcPt(0)
 231    ,fProNtracksLeadingJetGen(0)
 232    ,fProDelR80pcPtGen(0)
 233    ,fProNtracksLeadingJetBgrPerp2(0)
 234    ,fProNtracksLeadingJetRecPrim(0)
 235    ,fProDelR80pcPtRecPrim(0)
 236    ,fProNtracksLeadingJetRecSecNS(0)
 237    ,fProNtracksLeadingJetRecSecS(0)
 238    ,fProNtracksLeadingJetRecSecSsc(0)
 239
 240    ,fRandom(0)
 241 {
 242    // default constructor
 243   fBckgType[0] = 0;
 244   fBckgType[1] = 0;
 245   fBckgType[2] = 0;
 246   fBckgType[3] = 0;
 247   fBckgType[4] = 0;
 248
 249   for(Int_t ii=0; ii<5; ii++){
 250     fProDelRPtSum[ii]          = 0;
 251     fProDelRPtSumGen[ii]       = 0;
 252     fProDelRPtSumBgrPerp2[ii]  = 0;
 253     fProDelRPtSumRecPrim[ii]   = 0;
 254     fProDelRPtSumRecSecNS[ii]  = 0;
 255     fProDelRPtSumRecSecS[ii]   = 0;
 256     fProDelRPtSumRecSecSsc[ii] = 0;
 257   }
 258 }
 259
 260 //_______________________________________________________________________________________________
 261 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const char *name)
 262   : AliAnalysisTaskSE(name)
 263   ,fESD(0)
 264   ,fAOD(0)
 265   ,fAODJets(0)
 266   ,fAODExtension(0)
 267   ,fNonStdFile("")
 268   ,fBranchRecJets("jets")
 269   ,fBranchRecBckgClusters("")
 270   ,fBranchGenJets("")
 271   ,fBranchEmbeddedJets("")
 272   ,fTrackTypeGen(0)
 273   ,fJetTypeGen(0)
 274   ,fJetTypeRecEff(0)
 275   ,fUseAODInputJets(kTRUE)
 276   ,fFilterMask(0)
 277   ,fUsePhysicsSelection(kTRUE)
 278   ,fEvtSelectionMask(0)
 279   ,fEventClass(0)
 280   ,fMaxVertexZ(10)
 281   ,fTrackPtCut(0)
 282   ,fTrackEtaMin(0)
 283   ,fTrackEtaMax(0)
 284   ,fTrackPhiMin(0)
 285   ,fTrackPhiMax(0)
 286   ,fUseExtraTracks(0)
 287   ,fUseExtraTracksBgr(0)
 288   ,fCutFractionPtEmbedded(0)
 289   ,fUseEmbeddedJetAxis(0)
 290   ,fUseEmbeddedJetPt(0)
 291   ,fJetPtCut(0)
 292   ,fJetEtaMin(0)
 293   ,fJetEtaMax(0)
 294   ,fJetPhiMin(0)
 295   ,fJetPhiMax(0)
 296   ,fFFRadius(0)
 297   ,fFFMinLTrackPt(-1)
 298   ,fFFMaxTrackPt(-1)
 299   ,fFFMinnTracks(0)
 300   ,fFFBckgRadius(0)
 301   ,fBckgMode(0)
 302   ,fQAMode(0)
 303   ,fFFMode(0)
 304   ,fEffMode(0)
 305   ,fJSMode(0)
 306   ,fAvgTrials(0)
 307   ,fTracksRecCuts(0)
 308   ,fTracksGen(0)
 309   ,fTracksAODMCCharged(0)
 310   ,fTracksAODMCChargedSecNS(0)
 311   ,fTracksAODMCChargedSecS(0)
 312   ,fTracksRecQualityCuts(0)
 313   ,fJetsRec(0)
 314   ,fJetsRecCuts(0)
 315   ,fJetsGen(0)
 316   ,fJetsRecEff(0)
 317   ,fJetsEmbedded(0)
 318   ,fBckgJetsRec(0)
 319   ,fBckgJetsRecCuts(0)
 320   ,fBckgJetsGen(0)
 321   ,fQATrackHistosRecCuts(0)
 322   ,fQATrackHistosGen(0)
 323   ,fQAJetHistosRec(0)
 324   ,fQAJetHistosRecCuts(0)
 325   ,fQAJetHistosRecCutsLeading(0)
 326   ,fQAJetHistosGen(0)
 327   ,fQAJetHistosGenLeading(0)
 328   ,fQAJetHistosRecEffLeading(0)
 329   ,fFFHistosRecCuts(0)
 330   ,fFFHistosRecCutsInc(0)
 331   ,fFFHistosRecLeadingTrack(0)
 332   ,fFFHistosGen(0)
 333   ,fFFHistosGenInc(0)
 334   ,fFFHistosGenLeadingTrack(0)
 335   ,fQATrackHighPtThreshold(0)
 336   ,fFFNBinsJetPt(0)
 337   ,fFFJetPtMin(0)
 338   ,fFFJetPtMax(0)
 339   ,fFFNBinsPt(0)
 340   ,fFFPtMin(0)
 341   ,fFFPtMax(0)
 342   ,fFFNBinsXi(0)
 343   ,fFFXiMin(0)
 344   ,fFFXiMax(0)
 345   ,fFFNBinsZ(0)
 346   ,fFFZMin(0)
 347   ,fFFZMax(0)
 348   ,fQAJetNBinsPt(0)
 349   ,fQAJetPtMin(0)
 350   ,fQAJetPtMax(0)
 351   ,fQAJetNBinsEta(0)
 352   ,fQAJetEtaMin(0)
 353   ,fQAJetEtaMax(0)
 354   ,fQAJetNBinsPhi(0)
 355   ,fQAJetPhiMin(0)
 356   ,fQAJetPhiMax(0)
 357   ,fQATrackNBinsPt(0)
 358   ,fQATrackPtMin(0)
 359   ,fQATrackPtMax(0)
 360   ,fQATrackNBinsEta(0)
 361   ,fQATrackEtaMin(0)
 362   ,fQATrackEtaMax(0)
 363   ,fQATrackNBinsPhi(0)
 364   ,fQATrackPhiMin(0)
 365   ,fQATrackPhiMax(0)
 366   ,fCommonHistList(0)
 367   ,fh1EvtSelection(0)
 368   ,fh1VertexNContributors(0)
 369   ,fh1VertexZ(0)
 370   ,fh1EvtMult(0)
 371   ,fh1EvtCent(0)
 372   ,fh1Xsec(0)
 373   ,fh1Trials(0)
 374   ,fh1PtHard(0)
 375   ,fh1PtHardTrials(0)
 376   ,fh1nRecJetsCuts(0)
 377   ,fh1nGenJets(0)
 378   ,fh1nRecEffJets(0)
 379   ,fh1nEmbeddedJets(0)
 380   ,fh1nRecBckgJetsCuts(0)
 381   ,fh1nGenBckgJets(0)
 382   ,fh2PtRecVsGenPrim(0)
 383   ,fh2PtRecVsGenSec(0)
 384   ,fQATrackHistosRecEffGen(0)
 385   ,fQATrackHistosRecEffRec(0)
 386   ,fQATrackHistosSecRecNS(0)
 387   ,fQATrackHistosSecRecS(0)
 388   ,fQATrackHistosSecRecSsc(0)
 389   ,fFFHistosRecEffRec(0)
 390   ,fFFHistosSecRecNS(0)
 391   ,fFFHistosSecRecS(0)
 392   ,fFFHistosSecRecSsc(0)
 393   // Background
 394   ,fh1BckgMult0(0)
 395   ,fh1BckgMult1(0)
 396   ,fh1BckgMult2(0)
 397   ,fh1BckgMult3(0)
 398   ,fh1BckgMult4(0)
 399   ,fh1FractionPtEmbedded(0)
 400   ,fh1IndexEmbedded(0)
 401   ,fh2DeltaPtVsJetPtEmbedded(0)
 402   ,fh2DeltaPtVsRecJetPtEmbedded(0)
 403   ,fh1DeltaREmbedded(0)
 404   ,fQABckgHisto0RecCuts(0)
 405   ,fQABckgHisto0Gen(0)
 406   ,fQABckgHisto1RecCuts(0)
 407   ,fQABckgHisto1Gen(0)
 408   ,fQABckgHisto2RecCuts(0)
 409   ,fQABckgHisto2Gen(0)
 410   ,fQABckgHisto3RecCuts(0)
 411   ,fQABckgHisto3Gen(0)
 412   ,fQABckgHisto4RecCuts(0)
 413   ,fQABckgHisto4Gen(0)
 414   ,fFFBckgHisto0RecCuts(0)
 415   ,fFFBckgHisto0Gen(0)
 416   ,fFFBckgHisto1RecCuts(0)
 417   ,fFFBckgHisto1Gen(0)
 418   ,fFFBckgHisto2RecCuts(0)
 419   ,fFFBckgHisto2Gen(0)
 420   ,fFFBckgHisto3RecCuts(0)
 421   ,fFFBckgHisto3Gen(0)
 422   ,fFFBckgHisto4RecCuts(0)
 423   ,fFFBckgHisto4Gen(0)
 424   ,fFFBckgHisto0RecEffRec(0)
 425   ,fFFBckgHisto0SecRecNS(0)
 426   ,fFFBckgHisto0SecRecS(0)
 427   ,fFFBckgHisto0SecRecSsc(0)
 428   // jet shape
 429   ,fProNtracksLeadingJet(0)
 430   ,fProDelR80pcPt(0)
 431   ,fProNtracksLeadingJetGen(0)
 432   ,fProDelR80pcPtGen(0)
 433   ,fProNtracksLeadingJetBgrPerp2(0)
 434   ,fProNtracksLeadingJetRecPrim(0)
 435   ,fProDelR80pcPtRecPrim(0)
 436   ,fProNtracksLeadingJetRecSecNS(0)
 437   ,fProNtracksLeadingJetRecSecS(0)
 438   ,fProNtracksLeadingJetRecSecSsc(0)
 439   ,fRandom(0)
 440 {
 441   // constructor
 442   fBckgType[0] = 0;
 443   fBckgType[1] = 0;
 444   fBckgType[2] = 0;
 445   fBckgType[3] = 0;
 446   fBckgType[4] = 0;
 447
 448   for(Int_t ii=0; ii<5; ii++){
 449     fProDelRPtSum[ii]          = 0;
 450     fProDelRPtSumGen[ii]       = 0;
 451     fProDelRPtSumBgrPerp2[ii]  = 0;
 452     fProDelRPtSumRecPrim[ii]   = 0;
 453     fProDelRPtSumRecSecNS[ii]  = 0;
 454     fProDelRPtSumRecSecS[ii]   = 0;
 455     fProDelRPtSumRecSecSsc[ii] = 0;
 456   }
 457
 458   DefineOutput(1,TList::Class());
 459 }
 460
 461 //__________________________________________________________________________________________________________________________
 462 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const  AliAnalysisTaskFragmentationFunction &copy)
 463   : AliAnalysisTaskSE()
 464   ,fESD(copy.fESD)
 465   ,fAOD(copy.fAOD)
 466   ,fAODJets(copy.fAODJets)
 467   ,fAODExtension(copy.fAODExtension)
 468   ,fNonStdFile(copy.fNonStdFile)
 469   ,fBranchRecJets(copy.fBranchRecJets)
 470   ,fBranchRecBckgClusters(copy.fBranchRecBckgClusters)
 471   ,fBranchGenJets(copy.fBranchGenJets)
 472   ,fBranchEmbeddedJets(copy.fBranchEmbeddedJets)
 473   ,fTrackTypeGen(copy.fTrackTypeGen)
 474   ,fJetTypeGen(copy.fJetTypeGen)
 475   ,fJetTypeRecEff(copy.fJetTypeRecEff)
 476   ,fUseAODInputJets(copy.fUseAODInputJets)
 477   ,fFilterMask(copy.fFilterMask)
 478   ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
 479   ,fEvtSelectionMask(copy.fEvtSelectionMask)
 480   ,fEventClass(copy.fEventClass)
 481   ,fMaxVertexZ(copy.fMaxVertexZ)
 482   ,fTrackPtCut(copy.fTrackPtCut)
 483   ,fTrackEtaMin(copy.fTrackEtaMin)
 484   ,fTrackEtaMax(copy.fTrackEtaMax)
 485   ,fTrackPhiMin(copy.fTrackPhiMin)
 486   ,fTrackPhiMax(copy.fTrackPhiMax)
 487   ,fUseExtraTracks(copy.fUseExtraTracks)
 488   ,fUseExtraTracksBgr(copy.fUseExtraTracksBgr)
 489   ,fCutFractionPtEmbedded(copy.fCutFractionPtEmbedded)
 490   ,fUseEmbeddedJetAxis(copy.fUseEmbeddedJetAxis)
 491   ,fUseEmbeddedJetPt(copy.fUseEmbeddedJetPt)
 492   ,fJetPtCut(copy.fJetPtCut)
 493   ,fJetEtaMin(copy.fJetEtaMin)
 494   ,fJetEtaMax(copy.fJetEtaMax)
 495   ,fJetPhiMin(copy.fJetPhiMin)
 496   ,fJetPhiMax(copy.fJetPhiMax)
 497   ,fFFRadius(copy.fFFRadius)
 498   ,fFFMinLTrackPt(copy.fFFMinLTrackPt)
 499   ,fFFMaxTrackPt(copy.fFFMaxTrackPt)
 500   ,fFFMinnTracks(copy.fFFMinnTracks)
 501   ,fFFBckgRadius(copy.fFFBckgRadius)
 502   ,fBckgMode(copy.fBckgMode)
 503   ,fQAMode(copy.fQAMode)
 504   ,fFFMode(copy.fFFMode)
 505   ,fEffMode(copy.fEffMode)
 506   ,fJSMode(copy.fJSMode)
 507   ,fAvgTrials(copy.fAvgTrials)
 508   ,fTracksRecCuts(copy.fTracksRecCuts)
 509   ,fTracksGen(copy.fTracksGen)
 510   ,fTracksAODMCCharged(copy.fTracksAODMCCharged)
 511   ,fTracksAODMCChargedSecNS(copy.fTracksAODMCChargedSecNS)
 512   ,fTracksAODMCChargedSecS(copy.fTracksAODMCChargedSecS)
 513   ,fTracksRecQualityCuts(copy.fTracksRecQualityCuts)
 514   ,fJetsRec(copy.fJetsRec)
 515   ,fJetsRecCuts(copy.fJetsRecCuts)
 516   ,fJetsGen(copy.fJetsGen)
 517   ,fJetsRecEff(copy.fJetsRecEff)
 518   ,fJetsEmbedded(copy.fJetsEmbedded)
 519   ,fBckgJetsRec(copy.fBckgJetsRec)
 520   ,fBckgJetsRecCuts(copy.fBckgJetsRecCuts)
 521   ,fBckgJetsGen(copy.fBckgJetsGen)
 522   ,fQATrackHistosRecCuts(copy.fQATrackHistosRecCuts)
 523   ,fQATrackHistosGen(copy.fQATrackHistosGen)
 524   ,fQAJetHistosRec(copy.fQAJetHistosRec)
 525   ,fQAJetHistosRecCuts(copy.fQAJetHistosRecCuts)
 526   ,fQAJetHistosRecCutsLeading(copy.fQAJetHistosRecCutsLeading)
 527   ,fQAJetHistosGen(copy.fQAJetHistosGen)
 528   ,fQAJetHistosGenLeading(copy.fQAJetHistosGenLeading)
 529   ,fQAJetHistosRecEffLeading(copy.fQAJetHistosRecEffLeading)
 530   ,fFFHistosRecCuts(copy.fFFHistosRecCuts)
 531   ,fFFHistosRecCutsInc(copy.fFFHistosRecCutsInc)
 532   ,fFFHistosRecLeadingTrack(copy.fFFHistosRecLeadingTrack)
 533   ,fFFHistosGen(copy.fFFHistosGen)
 534   ,fFFHistosGenInc(copy.fFFHistosGenInc)
 535   ,fFFHistosGenLeadingTrack(copy.fFFHistosGenLeadingTrack)
 536   ,fQATrackHighPtThreshold(copy.fQATrackHighPtThreshold)
 537   ,fFFNBinsJetPt(copy.fFFNBinsJetPt)
 538   ,fFFJetPtMin(copy.fFFJetPtMin)
 539   ,fFFJetPtMax(copy.fFFJetPtMax)
 540   ,fFFNBinsPt(copy.fFFNBinsPt)
 541   ,fFFPtMin(copy.fFFPtMin)
 542   ,fFFPtMax(copy.fFFPtMax)
 543   ,fFFNBinsXi(copy.fFFNBinsXi)
 544   ,fFFXiMin(copy.fFFXiMin)
 545   ,fFFXiMax(copy.fFFXiMax)
 546   ,fFFNBinsZ(copy.fFFNBinsZ)
 547   ,fFFZMin(copy.fFFZMin)
 548   ,fFFZMax(copy.fFFZMax)
 549   ,fQAJetNBinsPt(copy.fQAJetNBinsPt)
 550   ,fQAJetPtMin(copy.fQAJetPtMin)
 551   ,fQAJetPtMax(copy.fQAJetPtMax)
 552   ,fQAJetNBinsEta(copy.fQAJetNBinsEta)
 553   ,fQAJetEtaMin(copy.fQAJetEtaMin)
 554   ,fQAJetEtaMax(copy.fQAJetEtaMax)
 555   ,fQAJetNBinsPhi(copy.fQAJetNBinsPhi)
 556   ,fQAJetPhiMin(copy.fQAJetPhiMin)
 557   ,fQAJetPhiMax(copy.fQAJetPhiMax)
 558   ,fQATrackNBinsPt(copy.fQATrackNBinsPt)
 559   ,fQATrackPtMin(copy.fQATrackPtMin)
 560   ,fQATrackPtMax(copy.fQATrackPtMax)
 561   ,fQATrackNBinsEta(copy.fQATrackNBinsEta)
 562   ,fQATrackEtaMin(copy.fQATrackEtaMin)
 563   ,fQATrackEtaMax(copy.fQATrackEtaMax)
 564   ,fQATrackNBinsPhi(copy.fQATrackNBinsPhi)
 565   ,fQATrackPhiMin(copy.fQATrackPhiMin)
 566   ,fQATrackPhiMax(copy.fQATrackPhiMax)
 567   ,fCommonHistList(copy.fCommonHistList)
 568   ,fh1EvtSelection(copy.fh1EvtSelection)
 569   ,fh1VertexNContributors(copy.fh1VertexNContributors)
 570   ,fh1VertexZ(copy.fh1VertexZ)
 571   ,fh1EvtMult(copy.fh1EvtMult)
 572   ,fh1EvtCent(copy.fh1EvtCent)
 573   ,fh1Xsec(copy.fh1Xsec)
 574   ,fh1Trials(copy.fh1Trials)
 575   ,fh1PtHard(copy.fh1PtHard)
 576   ,fh1PtHardTrials(copy.fh1PtHardTrials)
 577   ,fh1nRecJetsCuts(copy.fh1nRecJetsCuts)
 578   ,fh1nGenJets(copy.fh1nGenJets)
 579   ,fh1nRecEffJets(copy.fh1nRecEffJets)
 580   ,fh1nEmbeddedJets(copy.fh1nEmbeddedJets)
 581   ,fh1nRecBckgJetsCuts(copy.fh1nRecBckgJetsCuts)
 582   ,fh1nGenBckgJets(copy.fh1nGenBckgJets)
 583   ,fh2PtRecVsGenPrim(copy.fh2PtRecVsGenPrim)
 584   ,fh2PtRecVsGenSec(copy.fh2PtRecVsGenSec)
 585   ,fQATrackHistosRecEffGen(copy.fQATrackHistosRecEffGen)
 586   ,fQATrackHistosRecEffRec(copy.fQATrackHistosRecEffRec)
 587   ,fQATrackHistosSecRecNS(copy.fQATrackHistosSecRecNS)
 588   ,fQATrackHistosSecRecS(copy.fQATrackHistosSecRecS)
 589   ,fQATrackHistosSecRecSsc(copy.fQATrackHistosSecRecSsc)
 590   ,fFFHistosRecEffRec(copy.fFFHistosRecEffRec)
 591   ,fFFHistosSecRecNS(copy.fFFHistosSecRecNS)
 592   ,fFFHistosSecRecS(copy.fFFHistosSecRecS)
 593   ,fFFHistosSecRecSsc(copy.fFFHistosSecRecSsc)
 594   // Background
 595   ,fh1BckgMult0(copy.fh1BckgMult0)
 596   ,fh1BckgMult1(copy.fh1BckgMult1)
 597   ,fh1BckgMult2(copy.fh1BckgMult2)
 598   ,fh1BckgMult3(copy.fh1BckgMult3)
 599   ,fh1BckgMult4(copy.fh1BckgMult4)
 600   ,fh1FractionPtEmbedded(copy.fh1FractionPtEmbedded)
 601   ,fh1IndexEmbedded(copy.fh1IndexEmbedded)
 602   ,fh2DeltaPtVsJetPtEmbedded(copy.fh2DeltaPtVsJetPtEmbedded)
 603   ,fh2DeltaPtVsRecJetPtEmbedded(copy.fh2DeltaPtVsRecJetPtEmbedded)
 604   ,fh1DeltaREmbedded(copy.fh1DeltaREmbedded)
 605   ,fQABckgHisto0RecCuts(copy.fQABckgHisto0RecCuts)
 606   ,fQABckgHisto0Gen(copy.fQABckgHisto0Gen)
 607   ,fQABckgHisto1RecCuts(copy.fQABckgHisto1RecCuts)
 608   ,fQABckgHisto1Gen(copy.fQABckgHisto1Gen)
 609   ,fQABckgHisto2RecCuts(copy.fQABckgHisto2RecCuts)
 610   ,fQABckgHisto2Gen(copy.fQABckgHisto2Gen)
 611   ,fQABckgHisto3RecCuts(copy.fQABckgHisto3RecCuts)
 612   ,fQABckgHisto3Gen(copy.fQABckgHisto3Gen)
 613   ,fQABckgHisto4RecCuts(copy.fQABckgHisto4RecCuts)
 614   ,fQABckgHisto4Gen(copy.fQABckgHisto4Gen)
 615   ,fFFBckgHisto0RecCuts(copy.fFFBckgHisto0RecCuts)
 616   ,fFFBckgHisto0Gen(copy.fFFBckgHisto0Gen)
 617   ,fFFBckgHisto1RecCuts(copy.fFFBckgHisto1RecCuts)
 618   ,fFFBckgHisto1Gen(copy.fFFBckgHisto1Gen)
 619   ,fFFBckgHisto2RecCuts(copy.fFFBckgHisto2RecCuts)
 620   ,fFFBckgHisto2Gen(copy.fFFBckgHisto2Gen)
 621   ,fFFBckgHisto3RecCuts(copy.fFFBckgHisto3RecCuts)
 622   ,fFFBckgHisto3Gen(copy.fFFBckgHisto3Gen)
 623   ,fFFBckgHisto4RecCuts(copy.fFFBckgHisto4RecCuts)
 624   ,fFFBckgHisto4Gen(copy.fFFBckgHisto4Gen)
 625   ,fFFBckgHisto0RecEffRec(copy.fFFBckgHisto0RecEffRec)
 626   ,fFFBckgHisto0SecRecNS(copy.fFFBckgHisto0SecRecNS)
 627   ,fFFBckgHisto0SecRecS(copy.fFFBckgHisto0SecRecS)
 628   ,fFFBckgHisto0SecRecSsc(copy.fFFBckgHisto0SecRecSsc)
 629   // jet shape
 630   ,fProNtracksLeadingJet(copy.fProNtracksLeadingJet)
 631   ,fProDelR80pcPt(copy.fProDelR80pcPt)
 632   ,fProNtracksLeadingJetGen(copy.fProNtracksLeadingJetGen)
 633   ,fProDelR80pcPtGen(copy.fProDelR80pcPtGen)
 634   ,fProNtracksLeadingJetBgrPerp2(copy.fProNtracksLeadingJetBgrPerp2)
 635   ,fProNtracksLeadingJetRecPrim(copy.fProNtracksLeadingJetRecPrim)
 636   ,fProDelR80pcPtRecPrim(copy.fProDelR80pcPtRecPrim)
 637   ,fProNtracksLeadingJetRecSecNS(copy.fProNtracksLeadingJetRecSecNS)
 638   ,fProNtracksLeadingJetRecSecS(copy.fProNtracksLeadingJetRecSecS)
 639   ,fProNtracksLeadingJetRecSecSsc(copy.fProNtracksLeadingJetRecSecSsc)
 640   ,fRandom(copy.fRandom)
 641 {
 642   // copy constructor
 643   fBckgType[0] = copy.fBckgType[0];
 644   fBckgType[1] = copy.fBckgType[1];
 645   fBckgType[2] = copy.fBckgType[2];
 646   fBckgType[3] = copy.fBckgType[3];
 647   fBckgType[4] = copy.fBckgType[4];
 648
 649
 650   for(Int_t ii=0; ii<5; ii++){
 651     fProDelRPtSum[ii]          = copy.fProDelRPtSum[ii];
 652     fProDelRPtSumGen[ii]       = copy.fProDelRPtSumGen[ii];
 653     fProDelRPtSumBgrPerp2[ii]  = copy.fProDelRPtSumBgrPerp2[ii];
 654     fProDelRPtSumRecPrim[ii]   = copy.fProDelRPtSumRecPrim[ii];
 655     fProDelRPtSumRecSecNS[ii]  = copy.fProDelRPtSumRecSecNS[ii];
 656     fProDelRPtSumRecSecS[ii]   = copy.fProDelRPtSumRecSecS[ii];
 657     fProDelRPtSumRecSecSsc[ii] = copy.fProDelRPtSumRecSecSsc[ii];
 658   }
 659 }
 660
 661 // _________________________________________________________________________________________________________________________________
 662 AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::operator=(const AliAnalysisTaskFragmentationFunction& o)
 663 {
 664   // assignment
 665
 666   if(this!=&o){
 667
 668     AliAnalysisTaskSE::operator=(o);
 669     fESD                           = o.fESD;
 670     fAOD                           = o.fAOD;
 671     fAODJets                       = o.fAODJets;
 672     fAODExtension                  = o.fAODExtension;
 673     fNonStdFile                    = o.fNonStdFile;
 674     fBranchRecJets                 = o.fBranchRecJets;
 675     fBranchRecBckgClusters         = o.fBranchRecBckgClusters;
 676     fBranchGenJets                 = o.fBranchGenJets;
 677     fBranchEmbeddedJets            = o.fBranchEmbeddedJets;
 678     fTrackTypeGen                  = o.fTrackTypeGen;
 679     fJetTypeGen                    = o.fJetTypeGen;
 680     fJetTypeRecEff                 = o.fJetTypeRecEff;
 681     fUseAODInputJets               = o.fUseAODInputJets;
 682     fFilterMask                    = o.fFilterMask;
 683     fUsePhysicsSelection           = o.fUsePhysicsSelection;
 684     fEvtSelectionMask              = o.fEvtSelectionMask;
 685     fEventClass                    = o.fEventClass;
 686     fMaxVertexZ                    = o.fMaxVertexZ;
 687     fTrackPtCut                    = o.fTrackPtCut;
 688     fTrackEtaMin                   = o.fTrackEtaMin;
 689     fTrackEtaMax                   = o.fTrackEtaMax;
 690     fTrackPhiMin                   = o.fTrackPhiMin;
 691     fTrackPhiMax                   = o.fTrackPhiMax;
 692     fUseExtraTracks                = o.fUseExtraTracks;
 693     fUseExtraTracksBgr             = o.fUseExtraTracksBgr;
 694     fCutFractionPtEmbedded         = o.fCutFractionPtEmbedded;
 695     fUseEmbeddedJetAxis            = o.fUseEmbeddedJetAxis;
 696     fUseEmbeddedJetPt              = o.fUseEmbeddedJetPt;
 697     fJetPtCut                      = o.fJetPtCut;
 698     fJetEtaMin                     = o.fJetEtaMin;
 699     fJetEtaMax                     = o.fJetEtaMax;
 700     fJetPhiMin                     = o.fJetPhiMin;
 701     fJetPhiMax                     = o.fJetPhiMin;
 702     fFFRadius                      = o.fFFRadius;
 703     fFFMinLTrackPt                 = o.fFFMinLTrackPt;
 704     fFFMaxTrackPt                  = o.fFFMaxTrackPt;
 705     fFFMinnTracks                  = o.fFFMinnTracks;
 706     fFFBckgRadius                  = o.fFFBckgRadius;
 707     fBckgMode                      = o.fBckgMode;
 708     fQAMode                        = o.fQAMode;
 709     fFFMode                        = o.fFFMode;
 710     fEffMode                       = o.fEffMode;
 711     fJSMode                        = o.fJSMode;
 712     fBckgType[0]                   = o.fBckgType[0];
 713     fBckgType[1]                   = o.fBckgType[1];
 714     fBckgType[2]                   = o.fBckgType[2];
 715     fBckgType[3]                   = o.fBckgType[3];
 716     fBckgType[4]                   = o.fBckgType[4];
 717     fAvgTrials                     = o.fAvgTrials;
 718     fTracksRecCuts                 = o.fTracksRecCuts;
 719     fTracksGen                     = o.fTracksGen;
 720     fTracksAODMCCharged            = o.fTracksAODMCCharged;
 721     fTracksAODMCChargedSecNS       = o.fTracksAODMCChargedSecNS;
 722     fTracksAODMCChargedSecS        = o.fTracksAODMCChargedSecS;
 723     fTracksRecQualityCuts          = o.fTracksRecQualityCuts;
 724     fJetsRec                       = o.fJetsRec;
 725     fJetsRecCuts                   = o.fJetsRecCuts;
 726     fJetsGen                       = o.fJetsGen;
 727     fJetsRecEff                    = o.fJetsRecEff;
 728     fJetsEmbedded                  = o.fJetsEmbedded;
 729     fBckgJetsRec                   = o.fBckgJetsRec;
 730     fBckgJetsRecCuts               = o.fBckgJetsRecCuts;
 731     fBckgJetsGen                   = o.fBckgJetsGen;
 732     fQATrackHistosRecCuts          = o.fQATrackHistosRecCuts;
 733     fQATrackHistosGen              = o.fQATrackHistosGen;
 734     fQAJetHistosRec                = o.fQAJetHistosRec;
 735     fQAJetHistosRecCuts            = o.fQAJetHistosRecCuts;
 736     fQAJetHistosRecCutsLeading     = o.fQAJetHistosRecCutsLeading;
 737     fQAJetHistosGen                = o.fQAJetHistosGen;
 738     fQAJetHistosGenLeading         = o.fQAJetHistosGenLeading;
 739     fQAJetHistosRecEffLeading      = o.fQAJetHistosRecEffLeading;
 740     fFFHistosRecCuts               = o.fFFHistosRecCuts;
 741     fFFHistosRecCutsInc            = o.fFFHistosRecCutsInc;
 742     fFFHistosRecLeadingTrack       = o.fFFHistosRecLeadingTrack;
 743     fFFHistosGen                   = o.fFFHistosGen;
 744     fFFHistosGenInc                = o.fFFHistosGenInc;
 745     fFFHistosGenLeadingTrack       = o.fFFHistosGenLeadingTrack;
 746     fQATrackHighPtThreshold        = o.fQATrackHighPtThreshold;
 747     fFFNBinsJetPt                  = o.fFFNBinsJetPt;
 748     fFFJetPtMin                    = o.fFFJetPtMin;
 749     fFFJetPtMax                    = o.fFFJetPtMax;
 750     fFFNBinsPt                     = o.fFFNBinsPt;
 751     fFFPtMin                       = o.fFFPtMin;
 752     fFFPtMax                       = o.fFFPtMax;
 753     fFFNBinsXi                     = o.fFFNBinsXi;
 754     fFFXiMin                       = o.fFFXiMin;
 755     fFFXiMax                       = o.fFFXiMax;
 756     fFFNBinsZ                      = o.fFFNBinsZ;
 757     fFFZMin                        = o.fFFZMin;
 758     fFFZMax                        = o.fFFZMax;
 759     fQAJetNBinsPt                  = o.fQAJetNBinsPt;
 760     fQAJetPtMin                    = o.fQAJetPtMin;
 761     fQAJetPtMax                    = o.fQAJetPtMax;
 762     fQAJetNBinsEta                 = o.fQAJetNBinsEta;
 763     fQAJetEtaMin                   = o.fQAJetEtaMin;
 764     fQAJetEtaMax                   = o.fQAJetEtaMax;
 765     fQAJetNBinsPhi                 = o.fQAJetNBinsPhi;
 766     fQAJetPhiMin                   = o.fQAJetPhiMin;
 767     fQAJetPhiMax                   = o.fQAJetPhiMax;
 768     fQATrackNBinsPt                = o.fQATrackNBinsPt;
 769     fQATrackPtMin                  = o.fQATrackPtMin;
 770     fQATrackPtMax                  = o.fQATrackPtMax;
 771     fQATrackNBinsEta               = o.fQATrackNBinsEta;
 772     fQATrackEtaMin                 = o.fQATrackEtaMin;
 773     fQATrackEtaMax                 = o.fQATrackEtaMax;
 774     fQATrackNBinsPhi               = o.fQATrackNBinsPhi;
 775     fQATrackPhiMin                 = o.fQATrackPhiMin;
 776     fQATrackPhiMax                 = o.fQATrackPhiMax;
 777     fCommonHistList                = o.fCommonHistList;
 778     fh1EvtSelection                = o.fh1EvtSelection;
 779     fh1VertexNContributors         = o.fh1VertexNContributors;
 780     fh1VertexZ                     = o.fh1VertexZ;
 781     fh1EvtMult                     = o.fh1EvtMult;
 782     fh1EvtCent                     = o.fh1EvtCent;
 783     fh1Xsec                        = o.fh1Xsec;
 784     fh1Trials                      = o.fh1Trials;
 785     fh1PtHard                      = o.fh1PtHard;
 786     fh1PtHardTrials                = o.fh1PtHardTrials;
 787     fh1nRecJetsCuts                = o.fh1nRecJetsCuts;
 788     fh1nGenJets                    = o.fh1nGenJets;
 789     fh1nRecEffJets                 = o.fh1nRecEffJets;
 790     fh1nEmbeddedJets               = o.fh1nEmbeddedJets;
 791     fh2PtRecVsGenPrim              = o.fh2PtRecVsGenPrim;
 792     fh2PtRecVsGenSec               = o.fh2PtRecVsGenSec;
 793     fQATrackHistosRecEffGen        = o.fQATrackHistosRecEffGen;
 794     fQATrackHistosRecEffRec        = o.fQATrackHistosRecEffRec;
 795     fQATrackHistosSecRecNS         = o.fQATrackHistosSecRecNS;
 796     fQATrackHistosSecRecS          = o.fQATrackHistosSecRecS;
 797     fQATrackHistosSecRecSsc        = o.fQATrackHistosSecRecSsc;
 798     fFFHistosRecEffRec             = o.fFFHistosRecEffRec;
 799     fFFHistosSecRecNS              = o.fFFHistosSecRecNS;
 800     fFFHistosSecRecS               = o.fFFHistosSecRecS;
 801     fFFHistosSecRecSsc             = o.fFFHistosSecRecSsc;
 802     // Background
 803     fh1BckgMult0                   = o.fh1BckgMult0;
 804     fh1BckgMult1                   = o.fh1BckgMult1;
 805     fh1BckgMult2                   = o.fh1BckgMult2;
 806     fh1BckgMult3                   = o.fh1BckgMult3;
 807     fh1BckgMult4                   = o.fh1BckgMult4;
 808     fh1FractionPtEmbedded          = o.fh1FractionPtEmbedded;
 809     fh1IndexEmbedded               = o.fh1IndexEmbedded;
 810     fh2DeltaPtVsJetPtEmbedded      = o.fh2DeltaPtVsJetPtEmbedded;
 811     fh2DeltaPtVsRecJetPtEmbedded   = o.fh2DeltaPtVsRecJetPtEmbedded;
 812     fh1DeltaREmbedded              = o.fh1DeltaREmbedded;
 813     fQABckgHisto0RecCuts           = o.fQABckgHisto0RecCuts;
 814     fQABckgHisto0Gen               = o.fQABckgHisto0Gen;
 815     fQABckgHisto1RecCuts           = o.fQABckgHisto1RecCuts;
 816     fQABckgHisto1Gen               = o.fQABckgHisto1Gen;
 817     fQABckgHisto2RecCuts           = o.fQABckgHisto2RecCuts;
 818     fQABckgHisto2Gen               = o.fQABckgHisto2Gen;
 819     fQABckgHisto3RecCuts           = o.fQABckgHisto3RecCuts;
 820     fQABckgHisto3Gen               = o.fQABckgHisto3Gen;
 821     fQABckgHisto4RecCuts           = o.fQABckgHisto4RecCuts;
 822     fQABckgHisto4Gen               = o.fQABckgHisto4Gen;
 823     fFFBckgHisto0RecCuts           = o.fFFBckgHisto0RecCuts;
 824     fFFBckgHisto0Gen               = o.fFFBckgHisto0Gen;
 825     fFFBckgHisto1RecCuts           = o.fFFBckgHisto1RecCuts;
 826     fFFBckgHisto1Gen               = o.fFFBckgHisto1Gen;
 827     fFFBckgHisto2RecCuts           = o.fFFBckgHisto2RecCuts;
 828     fFFBckgHisto2Gen               = o.fFFBckgHisto2Gen;
 829     fFFBckgHisto3RecCuts           = o.fFFBckgHisto3RecCuts;
 830     fFFBckgHisto3Gen               = o.fFFBckgHisto3Gen;
 831     fFFBckgHisto4RecCuts           = o.fFFBckgHisto4RecCuts;
 832     fFFBckgHisto4Gen               = o.fFFBckgHisto4Gen;
 833     fFFBckgHisto0RecEffRec         = o.fFFBckgHisto0RecEffRec;
 834     fFFBckgHisto0SecRecNS          = o.fFFBckgHisto0SecRecNS;
 835     fFFBckgHisto0SecRecS           = o.fFFBckgHisto0SecRecS;
 836     fFFBckgHisto0SecRecSsc         = o.fFFBckgHisto0SecRecSsc;
 837     fProNtracksLeadingJet          = o.fProNtracksLeadingJet;
 838     fProDelR80pcPt                 = o.fProDelR80pcPt;
 839     fProNtracksLeadingJetGen       = o.fProNtracksLeadingJetGen;
 840     fProDelR80pcPtGen              = o.fProDelR80pcPtGen;
 841     fProNtracksLeadingJetBgrPerp2  = o.fProNtracksLeadingJetBgrPerp2;
 842     fProNtracksLeadingJetRecPrim   = o.fProNtracksLeadingJetRecPrim;
 843     fProDelR80pcPtRecPrim          = o.fProDelR80pcPtRecPrim;
 844     fProNtracksLeadingJetRecSecNS  = o.fProNtracksLeadingJetRecSecNS;
 845     fProNtracksLeadingJetRecSecS   = o.fProNtracksLeadingJetRecSecS;
 846     fProNtracksLeadingJetRecSecSsc = o.fProNtracksLeadingJetRecSecSsc;
 847     fRandom                        = o.fRandom;
 848
 849     for(Int_t ii=0; ii<5; ii++){
 850       fProDelRPtSum[ii]           = o.fProDelRPtSum[ii];
 851       fProDelRPtSumGen[ii]        = o.fProDelRPtSumGen[ii];
 852       fProDelRPtSumBgrPerp2[ii]   = o.fProDelRPtSumBgrPerp2[ii];
 853       fProDelRPtSumRecPrim[ii]    = o.fProDelRPtSumRecPrim[ii];
 854       fProDelRPtSumRecSecNS[ii]   = o.fProDelRPtSumRecSecNS[ii];
 855       fProDelRPtSumRecSecS[ii]    = o.fProDelRPtSumRecSecS[ii];
 856       fProDelRPtSumRecSecSsc[ii]  = o.fProDelRPtSumRecSecSsc[ii];
 857     }
 858   }
 859
 860   return *this;
 861 }
 862
 863 //___________________________________________________________________________
 864 AliAnalysisTaskFragmentationFunction::~AliAnalysisTaskFragmentationFunction()
 865 {
 866   // destructor
 867
 868   if(fTracksRecCuts)           delete fTracksRecCuts;
 869   if(fTracksGen)               delete fTracksGen;
 870   if(fTracksAODMCCharged)      delete fTracksAODMCCharged;
 871   if(fTracksAODMCChargedSecNS) delete fTracksAODMCChargedSecNS;
 872   if(fTracksAODMCChargedSecS)  delete fTracksAODMCChargedSecS;
 873   if(fTracksRecQualityCuts)    delete fTracksRecQualityCuts;
 874   if(fJetsRec)                 delete fJetsRec;
 875   if(fJetsRecCuts)             delete fJetsRecCuts;
 876   if(fJetsGen)                 delete fJetsGen;
 877   if(fJetsRecEff)              delete fJetsRecEff;
 878   if(fJetsEmbedded)            delete fJetsEmbedded;
 879
 880   if(fBckgMode &&
 881      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
 882       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
 883       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)){
 884
 885     if(fBckgJetsRec)          delete fBckgJetsRec;
 886     if(fBckgJetsRecCuts)      delete fBckgJetsRecCuts;
 887     if(fBckgJetsGen)          delete fBckgJetsGen;
 888   }
 889   if(fRandom)               delete fRandom;
 890 }
 891
 892 //______________________________________________________________________________________________________
 893 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const char* name,
 894                                                          Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
 895                                                          Int_t nPt, Float_t ptMin, Float_t ptMax,
 896                                                          Int_t nXi, Float_t xiMin, Float_t xiMax,
 897                                                          Int_t nZ , Float_t zMin , Float_t zMax)
 898   : TObject()
 899   ,fNBinsJetPt(nJetPt)
 900   ,fJetPtMin(jetPtMin)
 901   ,fJetPtMax(jetPtMax)
 902   ,fNBinsPt(nPt)
 903   ,fPtMin(ptMin)
 904   ,fPtMax(ptMax)
 905   ,fNBinsXi(nXi)
 906   ,fXiMin(xiMin)
 907   ,fXiMax(xiMax)
 908   ,fNBinsZ(nZ)
 909   ,fZMin(zMin)
 910   ,fZMax(zMax)
 911   ,fh2TrackPt(0)
 912   ,fh2Xi(0)
 913   ,fh2Z(0)
 914   ,fh1JetPt(0)
 915   ,fNameFF(name)
 916 {
 917   // default constructor
 918
 919 }
 920
 921 //___________________________________________________________________________
 922 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const AliFragFuncHistos& copy)
 923   : TObject()
 924   ,fNBinsJetPt(copy.fNBinsJetPt)
 925   ,fJetPtMin(copy.fJetPtMin)
 926   ,fJetPtMax(copy.fJetPtMax)
 927   ,fNBinsPt(copy.fNBinsPt)
 928   ,fPtMin(copy.fPtMin)
 929   ,fPtMax(copy.fPtMax)
 930   ,fNBinsXi(copy.fNBinsXi)
 931   ,fXiMin(copy.fXiMin)
 932   ,fXiMax(copy.fXiMax)
 933   ,fNBinsZ(copy.fNBinsZ)
 934   ,fZMin(copy.fZMin)
 935   ,fZMax(copy.fZMax)
 936   ,fh2TrackPt(copy.fh2TrackPt)
 937   ,fh2Xi(copy.fh2Xi)
 938   ,fh2Z(copy.fh2Z)
 939   ,fh1JetPt(copy.fh1JetPt)
 940   ,fNameFF(copy.fNameFF)
 941 {
 942   // copy constructor
 943 }
 944
 945 //_______________________________________________________________________________________________________________________________________________________________
 946 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& o)
 947 {
 948   // assignment
 949
 950   if(this!=&o){
 951     TObject::operator=(o);
 952     fNBinsJetPt = o.fNBinsJetPt;
 953     fJetPtMin   = o.fJetPtMin;
 954     fJetPtMax   = o.fJetPtMax;
 955     fNBinsPt    = o.fNBinsPt;
 956     fPtMin      = o.fPtMin;
 957     fPtMax      = o.fPtMax;
 958     fNBinsXi    = o.fNBinsXi;
 959     fXiMin      = o.fXiMin;
 960     fXiMax      = o.fXiMax;
 961     fNBinsZ     = o.fNBinsZ;
 962     fZMin       = o.fZMin;
 963     fZMax       = o.fZMax;
 964     fh2TrackPt  = o.fh2TrackPt;
 965     fh2Xi       = o.fh2Xi;
 966     fh2Z        = o.fh2Z;
 967     fh1JetPt    = o.fh1JetPt;
 968     fNameFF     = o.fNameFF;
 969   }
 970
 971   return *this;
 972 }
 973
 974 //_________________________________________________________
 975 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::~AliFragFuncHistos()
 976 {
 977   // destructor
 978
 979   if(fh1JetPt)   delete fh1JetPt;
 980   if(fh2TrackPt) delete fh2TrackPt;
 981   if(fh2Xi)      delete fh2Xi;
 982   if(fh2Z)       delete fh2Z;
 983 }
 984
 985 //_________________________________________________________________
 986 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::DefineHistos()
 987 {
 988   // book FF histos
 989
 990   fh1JetPt   = new TH1F(Form("fh1FFJetPt%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
 991   fh2TrackPt = new TH2F(Form("fh2FFTrackPt%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsPt, fPtMin, fPtMax);
 992   fh2Z       = new TH2F(Form("fh2FFZ%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, fZMin, fZMax);
 993   fh2Xi      = new TH2F(Form("fh2FFXi%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsXi, fXiMin, fXiMax);
 994
 995   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries");
 996   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPt,"jet p_{T} [GeV/c]","p_{T} [GeV/c]","entries");
 997   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi,"jet p_{T} [GeV/c]","#xi", "entries");
 998   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z,"jet p_{T} [GeV/c]","z","entries");
 999 }
1000
1001 //_______________________________________________________________________________________________________________
1002 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::FillFF(Float_t trackPt, Float_t jetPt, Bool_t incrementJetPt, Float_t norm,
1003                                                                         Bool_t scaleStrangeness, Float_t scaleFacStrangeness)
1004 {
1005   // fill FF
1006
1007   if(incrementJetPt && norm) fh1JetPt->Fill(jetPt,1/norm);
1008   else if(incrementJetPt) fh1JetPt->Fill(jetPt);
1009
1010  // Added for proper normalization of FF background estimation
1011   // when zero track are found in the background region
1012   if((int)trackPt==-1) return;
1013
1014   if(norm)fh2TrackPt->Fill(jetPt,trackPt,1/norm);
1015   else if(scaleStrangeness) fh2TrackPt->Fill(jetPt,trackPt,scaleFacStrangeness);
1016   else fh2TrackPt->Fill(jetPt,trackPt);
1017
1018   Double_t z = 0.;
1019   if(jetPt>0) z = trackPt / jetPt;
1020   Double_t xi = 0;
1021   if(z>0) xi = TMath::Log(1/z);
1022
1023   if(trackPt>(1-1e-06)*jetPt && trackPt<(1+1e-06)*jetPt){ // case z=1 : move entry to last histo bin <1
1024     z  = 1-1e-06;
1025     xi = 1e-06;
1026   }
1027
1028
1029   if(norm){
1030     fh2Xi->Fill(jetPt,xi,1/norm);
1031     fh2Z->Fill(jetPt,z,1/norm);
1032   }
1033   else if(scaleStrangeness){
1034     fh2Xi->Fill(jetPt,xi,scaleFacStrangeness);
1035     fh2Z->Fill(jetPt,z,scaleFacStrangeness);
1036   }
1037   else {
1038     fh2Xi->Fill(jetPt,xi);
1039     fh2Z->Fill(jetPt,z);
1040   }
1041 }
1042
1043 //_________________________________________________________________________________
1044 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AddToOutput(TList* list) const
1045 {
1046   // add histos to list
1047
1048   list->Add(fh1JetPt);
1049
1050   list->Add(fh2TrackPt);
1051   list->Add(fh2Xi);
1052   list->Add(fh2Z);
1053 }
1054
1055 //_________________________________________________________________________________________________________
1056 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const char* name,
1057                                                                Int_t nPt,   Float_t ptMin,   Float_t ptMax,
1058                                                                Int_t nEta,  Float_t etaMin,  Float_t etaMax,
1059                                                                Int_t nPhi,  Float_t phiMin,  Float_t phiMax)
1060   : TObject()
1061   ,fNBinsPt(nPt)
1062   ,fPtMin(ptMin)
1063   ,fPtMax(ptMax)
1064   ,fNBinsEta(nEta)
1065   ,fEtaMin(etaMin)
1066   ,fEtaMax(etaMax)
1067   ,fNBinsPhi(nPhi)
1068   ,fPhiMin(phiMin)
1069   ,fPhiMax(phiMax)
1070   ,fh2EtaPhi(0)
1071   ,fh1Pt(0)
1072   ,fNameQAJ(name)
1073 {
1074   // default constructor
1075 }
1076
1077 //____________________________________________________________________________________
1078 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const AliFragFuncQAJetHistos& copy)
1079   : TObject()
1080   ,fNBinsPt(copy.fNBinsPt)
1081   ,fPtMin(copy.fPtMin)
1082   ,fPtMax(copy.fPtMax)
1083   ,fNBinsEta(copy.fNBinsEta)
1084   ,fEtaMin(copy.fEtaMin)
1085   ,fEtaMax(copy.fEtaMax)
1086   ,fNBinsPhi(copy.fNBinsPhi)
1087   ,fPhiMin(copy.fPhiMin)
1088   ,fPhiMax(copy.fPhiMax)
1089   ,fh2EtaPhi(copy.fh2EtaPhi)
1090   ,fh1Pt(copy.fh1Pt)
1091   ,fNameQAJ(copy.fNameQAJ)
1092 {
1093   // copy constructor
1094 }
1095
1096 //________________________________________________________________________________________________________________________________________________________________________
1097 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& o)
1098 {
1099   // assignment
1100
1101   if(this!=&o){
1102     TObject::operator=(o);
1103     fNBinsPt  = o.fNBinsPt;
1104     fPtMin    = o.fPtMin;
1105     fPtMax    = o.fPtMax;
1106     fNBinsEta = o.fNBinsEta;
1107     fEtaMin   = o.fEtaMin;
1108     fEtaMax   = o.fEtaMax;
1109     fNBinsPhi = o.fNBinsPhi;
1110     fPhiMin   = o.fPhiMin;
1111     fPhiMax   = o.fPhiMax;
1112     fh2EtaPhi = o.fh2EtaPhi;
1113     fh1Pt     = o.fh1Pt;
1114     fNameQAJ  = o.fNameQAJ;
1115   }
1116
1117   return *this;
1118 }
1119
1120 //______________________________________________________________
1121 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::~AliFragFuncQAJetHistos()
1122 {
1123   // destructor
1124
1125   if(fh2EtaPhi) delete fh2EtaPhi;
1126   if(fh1Pt)     delete fh1Pt;
1127 }
1128
1129 //____________________________________________________________________
1130 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::DefineHistos()
1131 {
1132   // book jet QA histos
1133
1134   fh2EtaPhi  = new TH2F(Form("fh2JetQAEtaPhi%s", fNameQAJ.Data()), Form("%s: #eta - #phi distribution", fNameQAJ.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
1135   fh1Pt      = new TH1F(Form("fh1JetQAPt%s", fNameQAJ.Data()), Form("%s: p_{T} distribution", fNameQAJ.Data()), fNBinsPt, fPtMin, fPtMax);
1136
1137   AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi");
1138   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
1139 }
1140
1141 //____________________________________________________________________________________________________
1142 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::FillJetQA(Float_t eta, Float_t phi, Float_t pt)
1143 {
1144   // fill jet QA histos
1145
1146   fh2EtaPhi->Fill( eta, phi);
1147   fh1Pt->Fill( pt );
1148 }
1149
1150 //____________________________________________________________________________________
1151 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AddToOutput(TList* list) const
1152 {
1153   // add histos to list
1154
1155   list->Add(fh2EtaPhi);
1156   list->Add(fh1Pt);
1157 }
1158
1159 //___________________________________________________________________________________________________________
1160 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const char* name,
1161                                                                    Int_t nPt, Float_t ptMin, Float_t ptMax,
1162                                                                    Int_t nEta, Float_t etaMin, Float_t etaMax,
1163                                                                    Int_t nPhi, Float_t phiMin, Float_t phiMax,
1164                                                                    Float_t ptThresh)
1165   : TObject()
1166   ,fNBinsPt(nPt)
1167   ,fPtMin(ptMin)
1168   ,fPtMax(ptMax)
1169   ,fNBinsEta(nEta)
1170   ,fEtaMin(etaMin)
1171   ,fEtaMax(etaMax)
1172   ,fNBinsPhi(nPhi)
1173   ,fPhiMin(phiMin)
1174   ,fPhiMax(phiMax)
1175   ,fHighPtThreshold(ptThresh)
1176   ,fh2EtaPhi(0)
1177   ,fh1Pt(0)
1178   ,fh2HighPtEtaPhi(0)
1179   ,fh2PhiPt(0)
1180   ,fNameQAT(name)
1181 {
1182   // default constructor
1183 }
1184
1185 //__________________________________________________________________________________________
1186 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const AliFragFuncQATrackHistos& copy)
1187   : TObject()
1188   ,fNBinsPt(copy.fNBinsPt)
1189   ,fPtMin(copy.fPtMin)
1190   ,fPtMax(copy.fPtMax)
1191   ,fNBinsEta(copy.fNBinsEta)
1192   ,fEtaMin(copy.fEtaMin)
1193   ,fEtaMax(copy.fEtaMax)
1194   ,fNBinsPhi(copy.fNBinsPhi)
1195   ,fPhiMin(copy.fPhiMin)
1196   ,fPhiMax(copy.fPhiMax)
1197   ,fHighPtThreshold(copy.fHighPtThreshold)
1198   ,fh2EtaPhi(copy.fh2EtaPhi)
1199   ,fh1Pt(copy.fh1Pt)
1200   ,fh2HighPtEtaPhi(copy.fh2HighPtEtaPhi)
1201   ,fh2PhiPt(copy.fh2PhiPt)
1202   ,fNameQAT(copy.fNameQAT)
1203 {
1204   // copy constructor
1205 }
1206
1207 // _____________________________________________________________________________________________________________________________________________________________________________
1208 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& o)
1209 {
1210   // assignment
1211
1212   if(this!=&o){
1213     TObject::operator=(o);
1214     fNBinsPt         = o.fNBinsPt;
1215     fPtMin           = o.fPtMin;
1216     fPtMax           = o.fPtMax;
1217     fNBinsEta        = o.fNBinsEta;
1218     fEtaMin          = o.fEtaMin;
1219     fEtaMax          = o.fEtaMax;
1220     fNBinsPhi        = o.fNBinsPhi;
1221     fPhiMin          = o.fPhiMin;
1222     fPhiMax          = o.fPhiMax;
1223     fHighPtThreshold = o.fHighPtThreshold;
1224     fh2EtaPhi        = o.fh2EtaPhi;
1225     fh1Pt            = o.fh1Pt;
1226     fh2HighPtEtaPhi  = o.fh2HighPtEtaPhi;
1227     fh2PhiPt         = o.fh2PhiPt;
1228     fNameQAT         = o.fNameQAT;
1229   }
1230
1231   return *this;
1232 }
1233
1234 //___________________________________________________________________
1235 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::~AliFragFuncQATrackHistos()
1236 {
1237   // destructor
1238
1239   if(fh2EtaPhi)       delete fh2EtaPhi;
1240   if(fh2HighPtEtaPhi) delete fh2HighPtEtaPhi;
1241   if(fh1Pt)           delete fh1Pt;
1242   if(fh2PhiPt)        delete fh2PhiPt;
1243 }
1244
1245 //______________________________________________________________________
1246 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::DefineHistos()
1247 {
1248   // book track QA histos
1249
1250   fh2EtaPhi       = new TH2F(Form("fh2TrackQAEtaPhi%s", fNameQAT.Data()), Form("%s: #eta - #phi distribution", fNameQAT.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
1251   fh2HighPtEtaPhi = new TH2F(Form("fh2TrackQAHighPtEtaPhi%s", fNameQAT.Data()), Form("%s: #eta - #phi distribution for high-p_{T}", fNameQAT.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
1252   fh1Pt           = new TH1F(Form("fh1TrackQAPt%s", fNameQAT.Data()), Form("%s: p_{T} distribution", fNameQAT.Data()), fNBinsPt, fPtMin, fPtMax);
1253     fh2PhiPt        = new TH2F(Form("fh2TrackQAPhiPt%s", fNameQAT.Data()), Form("%s: #phi - p_{T} distribution", fNameQAT.Data()), fNBinsPhi, fPhiMin, fPhiMax, fNBinsPt, fPtMin, fPtMax);
1254
1255   AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi");
1256   AliAnalysisTaskFragmentationFunction::SetProperties(fh2HighPtEtaPhi, "#eta", "#phi");
1257   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
1258   AliAnalysisTaskFragmentationFunction::SetProperties(fh2PhiPt, "#phi", "p_{T} [GeV/c]");
1259 }
1260
1261 //________________________________________________________________________________________________________
1262 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::FillTrackQA(Float_t eta, Float_t phi, Float_t pt, Bool_t weightPt, Float_t norm,
1263                                                                                     Bool_t scaleStrangeness, Float_t scaleFacStrangeness)
1264 {
1265   // fill track QA histos
1266   Float_t weight = 1.;
1267   if(weightPt) weight = pt;
1268   fh2EtaPhi->Fill( eta, phi, weight);
1269   if(scaleStrangeness) fh2EtaPhi->Fill( eta, phi, scaleFacStrangeness);
1270   if(pt > fHighPtThreshold) fh2HighPtEtaPhi->Fill( eta, phi, weight);
1271   if(pt > fHighPtThreshold && scaleStrangeness) fh2HighPtEtaPhi->Fill( eta, phi, weight);
1272   if(norm) fh1Pt->Fill( pt, 1/norm );
1273   else if(scaleStrangeness) fh1Pt->Fill(pt,scaleFacStrangeness);
1274   else  fh1Pt->Fill( pt );
1275
1276   if(scaleFacStrangeness) fh2PhiPt->Fill(phi, pt, scaleFacStrangeness);
1277   else fh2PhiPt->Fill(phi, pt);
1278 }
1279
1280 //______________________________________________________________________________________
1281 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AddToOutput(TList* list) const
1282 {
1283   // add histos to list
1284
1285   list->Add(fh2EtaPhi);
1286   list->Add(fh2HighPtEtaPhi);
1287   list->Add(fh1Pt);
1288   list->Add(fh2PhiPt);
1289 }
1290
1291 //_________________________________________________________________________________
1292 Bool_t AliAnalysisTaskFragmentationFunction::Notify()
1293 {
1294   //
1295   // Implemented Notify() to read the cross sections
1296   // and number of trials from pyxsec.root
1297   // (taken from AliAnalysisTaskJetSpectrum2)
1298   //
1299   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1300   Float_t xsection = 0;
1301   Float_t ftrials  = 1;
1302
1303   fAvgTrials = 1;
1304   if(tree){
1305     TFile *curfile = tree->GetCurrentFile();
1306     if (!curfile) {
1307       Error("Notify","No current file");
1308       return kFALSE;
1309     }
1310     if(!fh1Xsec||!fh1Trials){
1311       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1312       return kFALSE;
1313     }
1314     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1315     fh1Xsec->Fill("<#sigma>",xsection);
1316     // construct a poor man average trials
1317     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1318     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1319   }
1320
1321   // Set seed for backg study
1322   fRandom = new TRandom3();
1323   fRandom->SetSeed(0);
1324
1325   return kTRUE;
1326 }
1327
1328 //__________________________________________________________________
1329 void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
1330 {
1331   // create output objects
1332
1333   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
1334
1335   // create list of tracks and jets
1336
1337   fTracksRecCuts = new TList();
1338   fTracksRecCuts->SetOwner(kFALSE);
1339
1340   fTracksGen = new TList();
1341   fTracksGen->SetOwner(kFALSE);
1342
1343   fTracksAODMCCharged = new TList();
1344   fTracksAODMCCharged->SetOwner(kFALSE);
1345
1346   fTracksAODMCChargedSecNS = new TList();
1347   fTracksAODMCChargedSecNS->SetOwner(kFALSE);
1348
1349   fTracksAODMCChargedSecS = new TList();
1350   fTracksAODMCChargedSecS->SetOwner(kFALSE);
1351
1352   fTracksRecQualityCuts = new TList();
1353   fTracksRecQualityCuts->SetOwner(kFALSE);
1354
1355   fJetsRec = new TList();
1356   fJetsRec->SetOwner(kFALSE);
1357
1358   fJetsRecCuts = new TList();
1359   fJetsRecCuts->SetOwner(kFALSE);
1360
1361   fJetsGen = new TList();
1362   fJetsGen->SetOwner(kFALSE);
1363
1364   fJetsRecEff = new TList();
1365   fJetsRecEff->SetOwner(kFALSE);
1366
1367   fJetsEmbedded = new TList();
1368   fJetsEmbedded->SetOwner(kFALSE);
1369
1370
1371   if(fBckgMode &&
1372      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters ||  fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
1373       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
1374       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)){
1375
1376     fBckgJetsRec = new TList();
1377     fBckgJetsRec->SetOwner(kFALSE);
1378
1379     fBckgJetsRecCuts = new TList();
1380     fBckgJetsRecCuts->SetOwner(kFALSE);
1381
1382     fBckgJetsGen = new TList();
1383     fBckgJetsGen->SetOwner(kFALSE);
1384   }
1385
1386   //
1387   // Create histograms / output container
1388   //
1389
1390   OpenFile(1);
1391   fCommonHistList = new TList();
1392   fCommonHistList->SetOwner(kTRUE);
1393
1394   Bool_t oldStatus = TH1::AddDirectoryStatus();
1395   TH1::AddDirectory(kFALSE);
1396
1397
1398   // Histograms
1399   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
1400   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1401   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
1402   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
1403   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
1404   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
1405   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
1406
1407   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
1408   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
1409   fh1EvtMult                 = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,12000.);
1410   fh1EvtCent                 = new TH1F("fh1EvtCent","centrality",100,0.,100.);
1411
1412   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1413   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1414   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1415   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1416   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
1417   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
1418
1419   fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
1420   fh1nGenJets                = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5);
1421   fh1nRecEffJets             = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5);
1422   fh1nEmbeddedJets           = new TH1F("fh1nEmbeddedJets","embedded jets per event",10,-0.5,9.5);
1423
1424   fh2PtRecVsGenPrim          = new TH2F("fh2PtRecVsGenPrim","rec vs gen pt",fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax,fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax);
1425   fh2PtRecVsGenSec           = new TH2F("fh2PtRecVsGenSec","rec vs gen pt",fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax,fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax);
1426
1427   // embedding
1428   if(fBranchEmbeddedJets.Length()){
1429     fh1FractionPtEmbedded         = new TH1F("fh1FractionPtEmbedded","",200,0,2);
1430     fh1IndexEmbedded              = new TH1F("fh1IndexEmbedded","",11,-1,10);
1431     fh2DeltaPtVsJetPtEmbedded     = new TH2F("fh2DeltaPtVsJetPtEmbedded","",250,0,250,200,-100,100);
1432     fh2DeltaPtVsRecJetPtEmbedded  = new TH2F("fh2DeltaPtVsRecJetPtEmbedded","",250,0,250,200,-100,100);
1433     fh1DeltaREmbedded             = new TH1F("fh1DeltaREmbedded","",50,0,0.5);
1434     fh1nEmbeddedJets              = new TH1F("fh1nEmbeddedJets","embedded jets per event",10,-0.5,9.5);
1435   }
1436
1437
1438   if(fQAMode){
1439     if(fQAMode&1){ // track QA
1440        fQATrackHistosRecCuts      = new AliFragFuncQATrackHistos("RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1441                                                                 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1442                                                                 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1443                                                                 fQATrackHighPtThreshold);
1444       fQATrackHistosGen          = new AliFragFuncQATrackHistos("Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1445                                                                 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1446                                                                 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1447                                                                 fQATrackHighPtThreshold);
1448     }
1449
1450     if(fQAMode&2){ // jet QA
1451       fQAJetHistosRec            = new AliFragFuncQAJetHistos("Rec", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1452                                                               fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1453                                                               fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1454       fQAJetHistosRecCuts        = new AliFragFuncQAJetHistos("RecCuts", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1455                                                               fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1456                                                               fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1457       fQAJetHistosRecCutsLeading = new AliFragFuncQAJetHistos("RecCutsLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1458                                                               fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1459                                                               fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1460       fQAJetHistosGen            = new AliFragFuncQAJetHistos("Gen", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1461                                                               fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1462                                                               fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1463       fQAJetHistosGenLeading     = new AliFragFuncQAJetHistos("GenLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1464                                                               fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1465                                                               fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1466       if(fEffMode) fQAJetHistosRecEffLeading  = new AliFragFuncQAJetHistos("RecEffLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
1467                                                                            fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1468     }
1469   } // end: QA
1470
1471   if(fFFMode){
1472
1473     fFFHistosRecCuts         = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1474                                                      fFFNBinsPt, fFFPtMin, fFFPtMax,
1475                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,
1476                                                      fFFNBinsZ , fFFZMin , fFFZMax );
1477
1478
1479     fFFHistosRecCutsInc      = new AliFragFuncHistos("RecCutsInc", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1480                                                      fFFNBinsPt, fFFPtMin, fFFPtMax,
1481                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,
1482                                                      fFFNBinsZ , fFFZMin , fFFZMax );
1483
1484
1485     fFFHistosRecLeadingTrack  = new AliFragFuncHistos("RecLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1486                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1487                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1488                                                       fFFNBinsZ , fFFZMin , fFFZMax );
1489
1490     fFFHistosGen              = new AliFragFuncHistos("Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1491                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1492                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1493                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1494
1495     fFFHistosGenInc           = new AliFragFuncHistos("GenInc", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1496                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1497                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1498                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1499
1500     fFFHistosGenLeadingTrack  = new AliFragFuncHistos("GenLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1501                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1502                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1503                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1504
1505   } // end: FF
1506
1507   // efficiency
1508
1509   if(fEffMode){
1510     if(fQAMode&1){
1511       fQATrackHistosRecEffGen = new AliFragFuncQATrackHistos("RecEffGen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1512                                                              fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1513                                                              fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1514                                                              fQATrackHighPtThreshold);
1515
1516       fQATrackHistosRecEffRec = new AliFragFuncQATrackHistos("RecEffRec", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1517                                                              fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1518                                                              fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1519                                                              fQATrackHighPtThreshold);
1520
1521       fQATrackHistosSecRecNS   = new AliFragFuncQATrackHistos("SecRecNS", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1522                                                              fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1523                                                              fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1524                                                              fQATrackHighPtThreshold);
1525
1526       fQATrackHistosSecRecS    = new AliFragFuncQATrackHistos("SecRecS", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1527                                                              fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1528                                                              fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1529                                                              fQATrackHighPtThreshold);
1530
1531       fQATrackHistosSecRecSsc = new AliFragFuncQATrackHistos("SecRecSsc", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1532                                                                fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1533                                                                fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1534                                                                fQATrackHighPtThreshold);
1535
1536     }
1537     if(fFFMode){
1538       fFFHistosRecEffRec      = new AliFragFuncHistos("RecEffRec", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1539                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1540                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1541                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1542
1543       fFFHistosSecRecNS       = new AliFragFuncHistos("SecRecNS", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1544                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1545                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1546                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1547
1548       fFFHistosSecRecS        = new AliFragFuncHistos("SecRecS", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1549                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1550                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1551                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1552
1553       fFFHistosSecRecSsc      = new AliFragFuncHistos("SecRecSsc", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1554                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1555                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1556                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1557
1558     }
1559   } // end: efficiency
1560
1561   // Background
1562   if(fBckgMode){
1563     if(fBckgType[0]==kBckgNone){
1564       AliError("no bgr method selected !");
1565     }
1566
1567     TString title[5];
1568     for(Int_t i=0; i<5; i++){
1569       if(fBckgType[i]==kBckgPerp) title[i]="Perp";
1570       else if(fBckgType[i]==kBckgPerp2) title[i]="Perp2";
1571       else if(fBckgType[i]==kBckgPerp2Area) title[i]="Perp2Area";
1572       else if(fBckgType[i]==kBckgPerpWindow) title[i]="PerpW";
1573       else if(fBckgType[i]==kBckgASide) title[i]="ASide";
1574       else if(fBckgType[i]==kBckgASideWindow) title[i]="ASideW";
1575       else if(fBckgType[i]==kBckgOutLJ) title[i]="OutLeadingJet";
1576       else if(fBckgType[i]==kBckgOut2J) title[i]="Out2Jets";
1577       else if(fBckgType[i]==kBckgOut3J) title[i]="Out3Jets";
1578       else if(fBckgType[i]==kBckgOutAJ) title[i]="AllJets";
1579       else if(fBckgType[i]==kBckgOutLJStat) title[i]="OutLeadingJetStat";
1580       else if(fBckgType[i]==kBckgOut2JStat) title[i]="Out2JetsStat";
1581       else if(fBckgType[i]==kBckgOut3JStat) title[i]="Out3JetsStat";
1582       else if(fBckgType[i]==kBckgOutAJStat) title[i]="AllJetsStat";
1583       else if(fBckgType[i]==kBckgClustersOutLeading) title[i]="OutClusters";
1584       else if(fBckgType[i]==kBckgClusters) title[i]="MedianClusters";
1585       else if(fBckgType[i]==kBckgNone)  title[i]="";
1586       else printf("Please chose background method number %d!",i);
1587     }
1588
1589
1590     if(fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
1591        fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
1592        fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading){
1593
1594       fh1nRecBckgJetsCuts        = new TH1F("fh1nRecBckgJetsCuts","reconstructed background jets per event",10,-0.5,9.5);
1595       fh1nGenBckgJets            = new TH1F("fh1nGenBckgJets","generated background jets per event",10,-0.5,9.5);
1596     }
1597
1598
1599     fh1BckgMult0 = new TH1F("fh1BckgMult0","bckg mult "+title[0],500,0,500);
1600     if(fBckgType[1] != kBckgNone) fh1BckgMult1 = new TH1F("fh1BckgMult1","bckg mult "+title[1],500,0,500);
1601     if(fBckgType[2] != kBckgNone) fh1BckgMult2 = new TH1F("fh1BckgMult2","bckg mult "+title[2],500,0,500);
1602     if(fBckgType[3] != kBckgNone) fh1BckgMult3 = new TH1F("fh1BckgMult3","bckg mult "+title[3],500,0,500);
1603     if(fBckgType[4] != kBckgNone) fh1BckgMult4 = new TH1F("fh1BckgMult4","bckg mult "+title[4],500,0,500);
1604
1605
1606     if(fQAMode&1){
1607       fQABckgHisto0RecCuts      = new AliFragFuncQATrackHistos("Bckg"+title[0]+"RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1608                                                                fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1609                                                                fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1610                                                                fQATrackHighPtThreshold);
1611       fQABckgHisto0Gen          = new AliFragFuncQATrackHistos("Bckg"+title[0]+"Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1612                                                                fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1613                                                                fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1614                                                                fQATrackHighPtThreshold);
1615
1616       if(fBckgType[1] != kBckgNone){
1617         fQABckgHisto1RecCuts      = new AliFragFuncQATrackHistos("Bckg"+title[1]+"RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1618                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1619                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1620                                                                  fQATrackHighPtThreshold);
1621         fQABckgHisto1Gen          = new AliFragFuncQATrackHistos("Bckg"+title[1]+"Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1622                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1623                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1624                                                                  fQATrackHighPtThreshold);
1625       }
1626       if(fBckgType[2] != kBckgNone){
1627         fQABckgHisto2RecCuts      = new AliFragFuncQATrackHistos("Bckg"+title[2]+"RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1628                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1629                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1630                                                                  fQATrackHighPtThreshold);
1631         fQABckgHisto2Gen          = new AliFragFuncQATrackHistos("Bckg"+title[2]+"Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1632                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1633                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1634                                                                  fQATrackHighPtThreshold);
1635       }
1636       if(fBckgType[3] != kBckgNone){
1637         fQABckgHisto3RecCuts      = new AliFragFuncQATrackHistos("Bckg"+title[3]+"RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1638                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1639                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1640                                                                  fQATrackHighPtThreshold);
1641         fQABckgHisto3Gen          = new AliFragFuncQATrackHistos("Bckg"+title[3]+"Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1642                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1643                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1644                                                                  fQATrackHighPtThreshold);
1645       }
1646       if(fBckgType[4] != kBckgNone){
1647         fQABckgHisto4RecCuts      = new AliFragFuncQATrackHistos("Bckg"+title[4]+"RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1648                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1649                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1650                                                                  fQATrackHighPtThreshold);
1651         fQABckgHisto4Gen          = new AliFragFuncQATrackHistos("Bckg"+title[4]+"Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1652                                                                  fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1653                                                                  fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1654                                                                  fQATrackHighPtThreshold);
1655       }
1656     } // end: background QA
1657
1658     if(fFFMode){
1659       fFFBckgHisto0RecCuts    = new AliFragFuncHistos("Bckg"+title[0]+"RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1660                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1661                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1662                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1663
1664       fFFBckgHisto0Gen        = new AliFragFuncHistos("Bckg"+title[0]+"Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1665                                                       fFFNBinsPt, fFFPtMin, fFFPtMax,
1666                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,
1667                                                       fFFNBinsZ , fFFZMin , fFFZMax);
1668
1669       if(fBckgType[1] != kBckgNone){
1670         fFFBckgHisto1RecCuts    = new AliFragFuncHistos("Bckg"+title[1]+"RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1671                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1672                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1673                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1674         fFFBckgHisto1Gen        = new AliFragFuncHistos("Bckg"+title[1]+"Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1675                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1676                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1677                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1678       }
1679       if(fBckgType[2] != kBckgNone){
1680         fFFBckgHisto2RecCuts    = new AliFragFuncHistos("Bckg"+title[2]+"RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1681                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1682                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1683                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1684
1685         fFFBckgHisto2Gen        = new AliFragFuncHistos("Bckg"+title[2]+"Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1686                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1687                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1688                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1689       }
1690       if(fBckgType[3] != kBckgNone){
1691         fFFBckgHisto3RecCuts    = new AliFragFuncHistos("Bckg"+title[3]+"RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1692                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1693                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1694                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1695
1696         fFFBckgHisto3Gen        = new AliFragFuncHistos("Bckg"+title[3]+"Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1697                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1698                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1699                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1700       }
1701       if(fBckgType[4] != kBckgNone){
1702         fFFBckgHisto4RecCuts    = new AliFragFuncHistos("Bckg"+title[4]+"RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1703                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1704                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1705                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1706
1707         fFFBckgHisto4Gen        = new AliFragFuncHistos("Bckg"+title[4]+"Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1708                                                         fFFNBinsPt, fFFPtMin, fFFPtMax,
1709                                                         fFFNBinsXi, fFFXiMin, fFFXiMax,
1710                                                         fFFNBinsZ , fFFZMin , fFFZMax);
1711       }
1712       if(fEffMode){
1713         fFFBckgHisto0RecEffRec      = new AliFragFuncHistos("Bckg"+title[0]+"RecEffRec", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1714                                                             fFFNBinsPt, fFFPtMin, fFFPtMax,
1715                                                             fFFNBinsXi, fFFXiMin, fFFXiMax,
1716                                                             fFFNBinsZ , fFFZMin , fFFZMax);
1717
1718         fFFBckgHisto0SecRecNS       = new AliFragFuncHistos("Bckg"+title[0]+"SecRecNS", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1719                                                             fFFNBinsPt, fFFPtMin, fFFPtMax,
1720                                                             fFFNBinsXi, fFFXiMin, fFFXiMax,
1721                                                             fFFNBinsZ , fFFZMin , fFFZMax);
1722
1723         fFFBckgHisto0SecRecS        = new AliFragFuncHistos("Bckg"+title[0]+"SecRecS", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1724                                                             fFFNBinsPt, fFFPtMin, fFFPtMax,
1725                                                             fFFNBinsXi, fFFXiMin, fFFXiMax,
1726                                                             fFFNBinsZ , fFFZMin , fFFZMax);
1727
1728         fFFBckgHisto0SecRecSsc      = new AliFragFuncHistos("Bckg"+title[0]+"SecRecSsc", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1729                                                             fFFNBinsPt, fFFPtMin, fFFPtMax,
1730                                                             fFFNBinsXi, fFFXiMin, fFFXiMax,
1731                                                             fFFNBinsZ , fFFZMin , fFFZMax);
1732
1733       }
1734     } // end: background FF
1735
1736
1737   } // end: background
1738
1739
1740   // ____________ define histograms ____________________
1741
1742   if(fQAMode){
1743     if(fQAMode&1){ // track QA
1744       fQATrackHistosRecCuts->DefineHistos();
1745       fQATrackHistosGen->DefineHistos();
1746     }
1747
1748     if(fQAMode&2){ // jet QA
1749       fQAJetHistosRec->DefineHistos();
1750       fQAJetHistosRecCuts->DefineHistos();
1751       fQAJetHistosRecCutsLeading->DefineHistos();
1752       fQAJetHistosGen->DefineHistos();
1753       fQAJetHistosGenLeading->DefineHistos();
1754       if(fEffMode) fQAJetHistosRecEffLeading->DefineHistos();
1755     }
1756   }
1757
1758   if(fFFMode){
1759     fFFHistosRecCuts->DefineHistos();
1760     fFFHistosRecCutsInc->DefineHistos();
1761     fFFHistosRecLeadingTrack->DefineHistos();
1762     fFFHistosGen->DefineHistos();
1763     fFFHistosGenInc->DefineHistos();
1764     fFFHistosGenLeadingTrack->DefineHistos();
1765   }
1766
1767   if(fEffMode){
1768     if(fQAMode&1){
1769       fQATrackHistosRecEffGen->DefineHistos();
1770       fQATrackHistosRecEffRec->DefineHistos();
1771       fQATrackHistosSecRecNS->DefineHistos();
1772       fQATrackHistosSecRecS->DefineHistos();
1773       fQATrackHistosSecRecSsc->DefineHistos();
1774     }
1775     if(fFFMode){
1776       fFFHistosRecEffRec->DefineHistos();
1777       fFFHistosSecRecNS->DefineHistos();
1778       fFFHistosSecRecS->DefineHistos();
1779       fFFHistosSecRecSsc->DefineHistos();
1780     }
1781   } // end: efficiency
1782
1783   // Background
1784   if(fBckgMode){
1785     if(fFFMode){
1786       fFFBckgHisto0RecCuts->DefineHistos();
1787       fFFBckgHisto0Gen->DefineHistos();
1788       if(fBckgType[1] != kBckgNone) fFFBckgHisto1RecCuts->DefineHistos();
1789       if(fBckgType[1] != kBckgNone) fFFBckgHisto1Gen->DefineHistos();
1790       if(fBckgType[2] != kBckgNone) fFFBckgHisto2RecCuts->DefineHistos();
1791       if(fBckgType[2] != kBckgNone) fFFBckgHisto2Gen->DefineHistos();
1792       if(fBckgType[3] != kBckgNone) fFFBckgHisto3RecCuts->DefineHistos();
1793       if(fBckgType[3] != kBckgNone) fFFBckgHisto3Gen->DefineHistos();
1794       if(fBckgType[4] != kBckgNone) fFFBckgHisto4RecCuts->DefineHistos();
1795       if(fBckgType[4] != kBckgNone) fFFBckgHisto4Gen->DefineHistos();
1796
1797      if(fEffMode){
1798         fFFBckgHisto0RecEffRec->DefineHistos();
1799         fFFBckgHisto0SecRecNS->DefineHistos();
1800         fFFBckgHisto0SecRecS->DefineHistos();
1801         fFFBckgHisto0SecRecSsc->DefineHistos();
1802       }
1803     }
1804
1805     if(fQAMode&1){
1806       fQABckgHisto0RecCuts->DefineHistos();
1807       fQABckgHisto0Gen->DefineHistos();
1808       if(fBckgType[1] != kBckgNone) fQABckgHisto1RecCuts->DefineHistos();
1809       if(fBckgType[1] != kBckgNone) fQABckgHisto1Gen->DefineHistos();
1810       if(fBckgType[2] != kBckgNone) fQABckgHisto2RecCuts->DefineHistos();
1811       if(fBckgType[2] != kBckgNone) fQABckgHisto2Gen->DefineHistos();
1812       if(fBckgType[3] != kBckgNone) fQABckgHisto3RecCuts->DefineHistos();
1813       if(fBckgType[3] != kBckgNone) fQABckgHisto3Gen->DefineHistos();
1814       if(fBckgType[4] != kBckgNone) fQABckgHisto4RecCuts->DefineHistos();
1815       if(fBckgType[4] != kBckgNone) fQABckgHisto4Gen->DefineHistos();
1816     }
1817   } // end: background
1818
1819
1820   Bool_t genJets    = (fJetTypeGen != kJetsUndef) ? kTRUE : kFALSE;
1821   Bool_t genTracks  = (fTrackTypeGen != kTrackUndef) ? kTRUE : kFALSE;
1822   Bool_t recJetsEff = (fJetTypeRecEff != kJetsUndef) ? kTRUE : kFALSE;
1823
1824   fCommonHistList->Add(fh1EvtSelection);
1825   fCommonHistList->Add(fh1EvtMult);
1826   fCommonHistList->Add(fh1EvtCent);
1827   fCommonHistList->Add(fh1VertexNContributors);
1828   fCommonHistList->Add(fh1VertexZ);
1829   fCommonHistList->Add(fh1nRecJetsCuts);
1830   fCommonHistList->Add(fh1Xsec);
1831   fCommonHistList->Add(fh1Trials);
1832   fCommonHistList->Add(fh1PtHard);
1833   fCommonHistList->Add(fh1PtHardTrials);
1834
1835   if(genJets) fCommonHistList->Add(fh1nGenJets);
1836
1837   // FF histograms
1838   if(fFFMode){
1839     fFFHistosRecCuts->AddToOutput(fCommonHistList);
1840     fFFHistosRecCutsInc->AddToOutput(fCommonHistList);
1841     fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
1842
1843     if(genJets && genTracks){
1844       fFFHistosGen->AddToOutput(fCommonHistList);
1845       fFFHistosGenInc->AddToOutput(fCommonHistList);
1846       fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
1847     }
1848   }
1849
1850   // Background
1851   if(fBckgMode){
1852     if(fFFMode){
1853       fFFBckgHisto0RecCuts->AddToOutput(fCommonHistList);
1854       if(fBckgType[1] != kBckgNone) fFFBckgHisto1RecCuts->AddToOutput(fCommonHistList);
1855       if(fBckgType[2] != kBckgNone) fFFBckgHisto2RecCuts->AddToOutput(fCommonHistList);
1856       if(fBckgType[3] != kBckgNone) fFFBckgHisto3RecCuts->AddToOutput(fCommonHistList);
1857       if(fBckgType[4] != kBckgNone) fFFBckgHisto4RecCuts->AddToOutput(fCommonHistList);
1858
1859       if(genJets && genTracks){
1860         fFFBckgHisto0Gen->AddToOutput(fCommonHistList);
1861         if(fBckgType[1] != kBckgNone) fFFBckgHisto1Gen->AddToOutput(fCommonHistList);
1862         if(fBckgType[2] != kBckgNone) fFFBckgHisto2Gen->AddToOutput(fCommonHistList);
1863         if(fBckgType[3] != kBckgNone) fFFBckgHisto3Gen->AddToOutput(fCommonHistList);
1864         if(fBckgType[4] != kBckgNone) fFFBckgHisto4Gen->AddToOutput(fCommonHistList);
1865       }
1866
1867       if(fEffMode){
1868         fFFBckgHisto0RecEffRec->AddToOutput(fCommonHistList);
1869         fFFBckgHisto0SecRecNS->AddToOutput(fCommonHistList);
1870         fFFBckgHisto0SecRecS->AddToOutput(fCommonHistList);
1871         fFFBckgHisto0SecRecSsc->AddToOutput(fCommonHistList);
1872       }
1873     }
1874
1875     if(fQAMode&1){
1876       fQABckgHisto0RecCuts->AddToOutput(fCommonHistList);
1877       if(fBckgType[1] != kBckgNone) fQABckgHisto1RecCuts->AddToOutput(fCommonHistList);
1878       if(fBckgType[2] != kBckgNone) fQABckgHisto2RecCuts->AddToOutput(fCommonHistList);
1879       if(fBckgType[3] != kBckgNone) fQABckgHisto3RecCuts->AddToOutput(fCommonHistList);
1880       if(fBckgType[4] != kBckgNone) fQABckgHisto4RecCuts->AddToOutput(fCommonHistList);
1881       if(genJets && genTracks){
1882         fQABckgHisto0Gen->AddToOutput(fCommonHistList);
1883         if(fBckgType[1] != kBckgNone) fQABckgHisto1Gen->AddToOutput(fCommonHistList);
1884         if(fBckgType[2] != kBckgNone) fQABckgHisto2Gen->AddToOutput(fCommonHistList);
1885         if(fBckgType[3] != kBckgNone) fQABckgHisto3Gen->AddToOutput(fCommonHistList);
1886         if(fBckgType[4] != kBckgNone) fQABckgHisto4Gen->AddToOutput(fCommonHistList);
1887       }
1888     }
1889
1890     if(fh1BckgMult0) fCommonHistList->Add(fh1BckgMult0);
1891     if(fBckgType[1] != kBckgNone)  fCommonHistList->Add(fh1BckgMult1);
1892     if(fBckgType[2] != kBckgNone)  fCommonHistList->Add(fh1BckgMult2);
1893     if(fBckgType[3] != kBckgNone)  fCommonHistList->Add(fh1BckgMult3);
1894     if(fBckgType[4] != kBckgNone)  fCommonHistList->Add(fh1BckgMult4);
1895   }
1896
1897
1898   if(fBranchEmbeddedJets.Length()){
1899     fCommonHistList->Add(fh1FractionPtEmbedded);
1900     fCommonHistList->Add(fh1IndexEmbedded);
1901     fCommonHistList->Add(fh2DeltaPtVsJetPtEmbedded);
1902     fCommonHistList->Add(fh2DeltaPtVsRecJetPtEmbedded);
1903     fCommonHistList->Add(fh1DeltaREmbedded);
1904     fCommonHistList->Add(fh1nEmbeddedJets);
1905   }
1906
1907
1908   // QA
1909   if(fQAMode){
1910     if(fQAMode&1){ // track QA
1911       fQATrackHistosRecCuts->AddToOutput(fCommonHistList);
1912       if(genTracks) fQATrackHistosGen->AddToOutput(fCommonHistList);
1913     }
1914
1915     if(fQAMode&2){ // jet QA
1916       fQAJetHistosRec->AddToOutput(fCommonHistList);
1917       fQAJetHistosRecCuts->AddToOutput(fCommonHistList);
1918       fQAJetHistosRecCutsLeading->AddToOutput(fCommonHistList);
1919       if(recJetsEff && fEffMode) fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList);
1920       if(genJets){
1921         fQAJetHistosGen->AddToOutput(fCommonHistList);
1922         fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
1923       }
1924     }
1925   }
1926
1927   if(fBckgMode &&
1928      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
1929       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
1930       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)) {
1931     fCommonHistList->Add(fh1nRecBckgJetsCuts);
1932     if(genJets) fCommonHistList->Add(fh1nGenBckgJets);
1933   }
1934
1935
1936   if(fEffMode && recJetsEff && genTracks){
1937     if(fQAMode&1){
1938       fQATrackHistosRecEffGen->AddToOutput(fCommonHistList);
1939       fQATrackHistosRecEffRec->AddToOutput(fCommonHistList);
1940       fQATrackHistosSecRecNS->AddToOutput(fCommonHistList);
1941       fQATrackHistosSecRecS->AddToOutput(fCommonHistList);
1942       fQATrackHistosSecRecSsc->AddToOutput(fCommonHistList);
1943     }
1944     if(fFFMode){
1945       fFFHistosRecEffRec->AddToOutput(fCommonHistList);
1946       fFFHistosSecRecNS->AddToOutput(fCommonHistList);
1947       fFFHistosSecRecS->AddToOutput(fCommonHistList);
1948       fFFHistosSecRecSsc->AddToOutput(fCommonHistList);
1949     }
1950     fCommonHistList->Add(fh1nRecEffJets);
1951     fCommonHistList->Add(fh2PtRecVsGenPrim);
1952     fCommonHistList->Add(fh2PtRecVsGenSec);
1953   }
1954
1955   // jet shape
1956   if(fJSMode){
1957
1958     fProNtracksLeadingJet          = new TProfile("AvgNoOfTracksLeadingJet","AvgNoOfTracksLeadingJet",100,0,250,0,50);
1959     fProDelR80pcPt                 = new TProfile("AvgdelR80pcPt","AvgdelR80pcPt",100,0,250,0,1);
1960
1961     if(genJets && genTracks){
1962       fProNtracksLeadingJetGen       = new TProfile("AvgNoOfTracksLeadingJetGen","AvgNoOfTracksLeadingJetGen",100,0,250,0,50);
1963       fProDelR80pcPtGen              = new TProfile("AvgdelR80pcPtGen","AvgdelR80pcPtGen",100,0,250,0,1);
1964     }
1965
1966     if(fBckgMode)
1967       fProNtracksLeadingJetBgrPerp2  = new TProfile("AvgNoOfTracksLeadingJetBgrPerp2","AvgNoOfTracksLeadingJetBgrPerp2",100,0,250,0,50);
1968
1969     if(fEffMode){
1970       fProNtracksLeadingJetRecPrim   = new TProfile("AvgNoOfTracksLeadingJetRecPrim","AvgNoOfTracksLeadingJetRecPrim",100,0,250,0,50);
1971       fProDelR80pcPtRecPrim          = new TProfile("AvgdelR80pcPtRecPrim","AvgdelR80pcPtRecPrim",100,0,250,0,1);
1972       fProNtracksLeadingJetRecSecNS  = new TProfile("AvgNoOfTracksLeadingJetRecSecNS","AvgNoOfTracksLeadingJetRecSecNS",100,0,250,0,50);
1973       fProNtracksLeadingJetRecSecS   = new TProfile("AvgNoOfTracksLeadingJetRecSecS","AvgNoOfTracksLeadingJetRecSecS",100,0,250,0,50);
1974       fProNtracksLeadingJetRecSecSsc = new TProfile("AvgNoOfTracksLeadingJetRecSecSsc","AvgNoOfTracksLeadingJetRecSecSsc",100,0,250,0,50);
1975     }
1976
1977     TString strTitJS;
1978     for(Int_t ii=0; ii<5; ii++){
1979       if(ii==0)strTitJS = "_JetPt20to30";
1980       if(ii==1)strTitJS = "_JetPt30to40";
1981       if(ii==2)strTitJS = "_JetPt40to60";
1982       if(ii==3)strTitJS = "_JetPt60to80";
1983       if(ii==4)strTitJS = "_JetPt80to100";
1984
1985       fProDelRPtSum[ii]            = new TProfile(Form("AvgPtSumDelR%s",strTitJS.Data()),Form("AvgPtSumDelR%s",strTitJS.Data()),50,0,1,0,250);
1986       if(genJets && genTracks)
1987         fProDelRPtSumGen[ii]       = new TProfile(Form("AvgPtSumDelRGen%s",strTitJS.Data()),Form("AvgPtSumDelRGen%s",strTitJS.Data()),50,0,1,0,250);
1988       if(fBckgMode)
1989         fProDelRPtSumBgrPerp2[ii]  = new TProfile(Form("AvgPtSumDelRBgrPerp2%s",strTitJS.Data()),Form("AvgPtSumDelRBgrPerp2%s",strTitJS.Data()),50,0,1,0,250);
1990       if(fEffMode){
1991         fProDelRPtSumRecPrim[ii]   = new TProfile(Form("AvgPtSumDelRRecPrim%s",strTitJS.Data()),Form("AvgPtSumDelRRecPrim%s",strTitJS.Data()),50,0,1,0,250);
1992         fProDelRPtSumRecSecNS[ii]  = new TProfile(Form("AvgPtSumDelRRecSecNS%s",strTitJS.Data()),Form("AvgPtSumDelRRecSecNS%s",strTitJS.Data()),50,0,1,0,250);
1993         fProDelRPtSumRecSecS[ii]   = new TProfile(Form("AvgPtSumDelRRecSecS%s",strTitJS.Data()),Form("AvgPtSumDelRRecSecS%s",strTitJS.Data()),50,0,1,0,250);
1994         fProDelRPtSumRecSecSsc[ii] = new TProfile(Form("AvgPtSumDelRRecSecSsc%s",strTitJS.Data()),Form("AvgPtSumDelRRecSecSsc%s",strTitJS.Data()),50,0,1,0,250);
1995       }
1996     }
1997
1998     fCommonHistList->Add(fProNtracksLeadingJet);
1999     fCommonHistList->Add(fProDelR80pcPt);
2000     for(int ii=0; ii<5; ii++) fCommonHistList->Add(fProDelRPtSum[ii]);
2001
2002     if(genJets && genTracks){
2003       fCommonHistList->Add(fProNtracksLeadingJetGen);
2004       fCommonHistList->Add(fProDelR80pcPtGen);
2005       for(Int_t ii=0; ii<5; ii++)  fCommonHistList->Add(fProDelRPtSumGen[ii]);
2006     }
2007
2008     if(fBckgMode){
2009       fCommonHistList->Add(fProNtracksLeadingJetBgrPerp2);
2010       for(Int_t ii=0; ii<5; ii++) fCommonHistList->Add(fProDelRPtSumBgrPerp2[ii]);
2011     }
2012
2013     if(fEffMode){
2014       fCommonHistList->Add(fProNtracksLeadingJetRecPrim);
2015       fCommonHistList->Add(fProDelR80pcPtRecPrim);
2016       for(Int_t ii=0; ii<5; ii++)  fCommonHistList->Add(fProDelRPtSumRecPrim[ii]);
2017
2018       fCommonHistList->Add(fProNtracksLeadingJetRecSecNS);
2019       for(Int_t ii=0; ii<5; ii++)  fCommonHistList->Add(fProDelRPtSumRecSecNS[ii]);
2020
2021       fCommonHistList->Add(fProNtracksLeadingJetRecSecS);
2022       for(Int_t ii=0; ii<5; ii++)  fCommonHistList->Add(fProDelRPtSumRecSecS[ii]);
2023
2024       fCommonHistList->Add(fProNtracksLeadingJetRecSecSsc);
2025       for(Int_t ii=0; ii<5; ii++)  fCommonHistList->Add(fProDelRPtSumRecSecSsc[ii]);
2026     }
2027   }
2028
2029   // =========== Switch on Sumw2 for all histos ===========
2030   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
2031     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
2032     if (h1) h1->Sumw2();
2033     else{
2034       THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
2035       if(hnSparse) hnSparse->Sumw2();
2036     }
2037   }
2038
2039   TH1::AddDirectory(oldStatus);
2040
2041   PostData(1, fCommonHistList);
2042 }
2043
2044 //_______________________________________________
2045 void AliAnalysisTaskFragmentationFunction::Init()
2046 {
2047   // Initialization
2048   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::Init()");
2049
2050 }
2051
2052 //_____________________________________________________________
2053 void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
2054 {
2055   // Main loop
2056   // Called for each event
2057   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserExec()");
2058
2059
2060   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
2061
2062   // Trigger selection
2063   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2064     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2065
2066   if(!(inputHandler->IsEventSelected() & fEvtSelectionMask)){
2067     fh1EvtSelection->Fill(1.);
2068     if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
2069     PostData(1, fCommonHistList);
2070     return;
2071   }
2072
2073   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
2074   if(!fESD){
2075     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
2076   }
2077
2078   fMCEvent = MCEvent();
2079   if(!fMCEvent){
2080     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
2081   }
2082
2083   // get AOD event from input/ouput
2084   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
2085   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
2086     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
2087     if(fUseAODInputJets) fAODJets = fAOD;
2088     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
2089   }
2090   else {
2091     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
2092     if( handler && handler->InheritsFrom("AliAODHandler") ) {
2093       fAOD = ((AliAODHandler*)handler)->GetAOD();
2094       fAODJets = fAOD;
2095       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
2096     }
2097   }
2098
2099   if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
2100     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
2101     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
2102       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
2103       if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
2104     }
2105   }
2106
2107   if(fNonStdFile.Length()!=0){
2108     // case we have an AOD extension - fetch the jets from the extended output
2109
2110     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
2111     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
2112     if(!fAODExtension){
2113       if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
2114     }
2115   }
2116
2117   if(!fAOD){
2118     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
2119     return;
2120   }
2121   if(!fAODJets){
2122     Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
2123     return;
2124   }
2125
2126
2127   // event selection **************************************************
2128   // *** event class ***
2129   Double_t centPercent = -1;
2130   if(fEventClass>0){
2131     Int_t cl = 0;
2132     if(handler->InheritsFrom("AliAODInputHandler")){
2133       // since it is not supported by the helper task define own classes
2134       centPercent = fAOD->GetHeader()->GetCentrality();
2135       cl = 1;
2136       if(centPercent>10) cl = 2;
2137       if(centPercent>30) cl = 3;
2138       if(centPercent>50) cl = 4;
2139     }
2140     else {
2141       cl = AliAnalysisHelperJetTasks::EventClass();
2142       if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // retrieve value 'by hand'
2143     }
2144
2145     if(cl!=fEventClass){
2146       // event not in selected event class, reject event
2147       if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
2148       fh1EvtSelection->Fill(2.);
2149       PostData(1, fCommonHistList);
2150       return;
2151     }
2152   }
2153
2154   // *** vertex cut ***
2155   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
2156   Int_t nTracksPrim = primVtx->GetNContributors();
2157   fh1VertexNContributors->Fill(nTracksPrim);
2158
2159
2160   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
2161   if(!nTracksPrim){
2162     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
2163     fh1EvtSelection->Fill(3.);
2164     PostData(1, fCommonHistList);
2165     return;
2166   }
2167
2168   fh1VertexZ->Fill(primVtx->GetZ());
2169
2170   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
2171     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2172     fh1EvtSelection->Fill(4.);
2173     PostData(1, fCommonHistList);
2174     return;
2175   }
2176
2177   TString primVtxName(primVtx->GetName());
2178
2179   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
2180     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
2181     fh1EvtSelection->Fill(5.);
2182     PostData(1, fCommonHistList);
2183     return;
2184   }
2185
2186   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
2187   fh1EvtSelection->Fill(0.);
2188   fh1EvtCent->Fill(centPercent);
2189
2190
2191   //___ get MC information __________________________________________________________________
2192
2193   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
2194
2195   Double_t ptHard = 0.;
2196   Double_t nTrials = 1; // trials for MC trigger weight for real data
2197
2198   if(fMCEvent){
2199     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
2200
2201     if(genHeader){
2202
2203       AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
2204       AliGenHijingEventHeader*  hijingGenHeader = 0x0;
2205
2206       if(pythiaGenHeader){
2207         if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
2208         nTrials = pythiaGenHeader->Trials();
2209         ptHard  = pythiaGenHeader->GetPtHard();
2210
2211         fh1PtHard->Fill(ptHard);
2212         fh1PtHardTrials->Fill(ptHard,nTrials);
2213
2214       } else { // no pythia, hijing?
2215
2216         if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
2217
2218         hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
2219         if(!hijingGenHeader){
2220           Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
2221         } else {
2222           if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
2223         }
2224       }
2225
2226       //fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
2227     }
2228   }
2229
2230   //___ fetch jets __________________________________________________________________________
2231
2232   Int_t nJ = GetListOfJets(fJetsRec, kJetsRec);
2233   Int_t nRecJets = 0;
2234   if(nJ>=0) nRecJets = fJetsRec->GetEntries();
2235   if(fDebug>2)Printf("%s:%d Selected Rec jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
2236   if(nJ != nRecJets) Printf("%s:%d Mismatch Selected Rec Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
2237
2238   Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
2239   Int_t nRecJetsCuts = 0;
2240   if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
2241   if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2242   if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2243   fh1nRecJetsCuts->Fill(nRecJetsCuts);
2244
2245   if(fJetTypeGen==kJetsKine || fJetTypeGen == kJetsKineAcceptance) fJetsGen->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear()
2246
2247   Int_t nJGen  = GetListOfJets(fJetsGen, fJetTypeGen);
2248   Int_t nGenJets = 0;
2249   if(nJGen>=0) nGenJets = fJetsGen->GetEntries();
2250   if(fDebug>2)Printf("%s:%d Selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
2251
2252   if(nJGen != nGenJets) Printf("%s:%d Mismatch selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
2253   fh1nGenJets->Fill(nGenJets);
2254
2255
2256   if(fJetTypeRecEff==kJetsKine || fJetTypeRecEff == kJetsKineAcceptance) fJetsRecEff->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear()
2257   Int_t nJRecEff  = GetListOfJets(fJetsRecEff, fJetTypeRecEff);
2258   Int_t nRecEffJets = 0;
2259   if(nJRecEff>=0) nRecEffJets = fJetsRecEff->GetEntries();
2260   if(fDebug>2)Printf("%s:%d Selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
2261   if(nJRecEff != nRecEffJets) Printf("%s:%d Mismatch selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
2262   fh1nRecEffJets->Fill(nRecEffJets);
2263
2264
2265   Int_t nEmbeddedJets =  0;
2266   TArrayI iEmbeddedMatchIndex;
2267   TArrayF fEmbeddedPtFraction;
2268
2269
2270   if(fBranchEmbeddedJets.Length()){
2271     Int_t nJEmbedded = GetListOfJets(fJetsEmbedded, kJetsEmbedded);
2272     if(nJEmbedded>=0) nEmbeddedJets = fJetsEmbedded->GetEntries();
2273     if(fDebug>2)Printf("%s:%d Selected Embedded jets: %d %d",(char*)__FILE__,__LINE__,nJEmbedded,nEmbeddedJets);
2274     if(nJEmbedded != nEmbeddedJets) Printf("%s:%d Mismatch Selected Embedded Jets: %d %d",(char*)__FILE__,__LINE__,nJEmbedded,nEmbeddedJets);
2275     fh1nEmbeddedJets->Fill(nEmbeddedJets);
2276
2277     Float_t maxDist = 0.3;
2278
2279     iEmbeddedMatchIndex.Set(nEmbeddedJets);
2280     fEmbeddedPtFraction.Set(nEmbeddedJets);
2281
2282     iEmbeddedMatchIndex.Reset(-1);
2283     fEmbeddedPtFraction.Reset(0);
2284
2285     AliAnalysisHelperJetTasks::GetJetMatching(fJetsEmbedded, nEmbeddedJets,
2286                                               fJetsRecCuts, nRecJetsCuts,
2287                                               iEmbeddedMatchIndex, fEmbeddedPtFraction,
2288                                               fDebug, maxDist);
2289
2290   }
2291
2292   //____ fetch background clusters ___________________________________________________
2293   if(fBckgMode &&
2294      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
2295       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
2296       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)){
2297
2298     Int_t nBJ = GetListOfBckgJets(fBckgJetsRec, kJetsRec);
2299     Int_t nRecBckgJets = 0;
2300     if(nBJ>=0) nRecBckgJets = fBckgJetsRec->GetEntries();
2301     if(fDebug>2)Printf("%s:%d Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
2302     if(nBJ != nRecBckgJets) Printf("%s:%d Mismatch Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
2303
2304     Int_t nBJCuts = GetListOfBckgJets(fBckgJetsRecCuts, kJetsRecAcceptance);
2305     Int_t nRecBckgJetsCuts = 0;
2306     if(nBJCuts>=0) nRecBckgJetsCuts = fBckgJetsRecCuts->GetEntries();
2307     if(fDebug>2)Printf("%s:%d Selected Rec background jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2308     if(nRecBckgJetsCuts != nBJCuts) Printf("%s:%d Mismatch selected Rec background jets after cuts: %d %d",(char*)__FILE__,__LINE__,nBJCuts,nRecBckgJetsCuts);
2309     fh1nRecBckgJetsCuts->Fill(nRecBckgJetsCuts);
2310
2311     if(0){ // protection OB - not yet implemented
2312       if(fJetTypeGen==kJetsKine || fJetTypeGen == kJetsKineAcceptance) fBckgJetsGen->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear()
2313       Int_t nBJGen  = GetListOfBckgJets(fBckgJetsGen, fJetTypeGen);
2314       Int_t nGenBckgJets = 0;
2315       if(nBJGen>=0) nGenBckgJets = fBckgJetsGen->GetEntries();
2316       if(fDebug>2)Printf("%s:%d Selected Gen background jets: %d %d",(char*)__FILE__,__LINE__,nBJGen,nGenBckgJets);
2317       if(nBJGen != nGenBckgJets) Printf("%s:%d Mismatch selected Gen background jets: %d %d",(char*)__FILE__,__LINE__,nBJGen,nGenBckgJets);
2318       fh1nGenBckgJets->Fill(nGenBckgJets);
2319     }
2320   }
2321
2322
2323   //____ fetch particles __________________________________________________________
2324
2325   Int_t nTCuts;
2326   if(fUseExtraTracks ==  1)      nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODExtraCuts);
2327   else if(fUseExtraTracks == -1) nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODExtraonlyCuts);
2328   else                           nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
2329
2330   Int_t nRecPartCuts = 0;
2331   if(nTCuts>=0) nRecPartCuts = fTracksRecCuts->GetEntries();
2332   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
2333   if(nRecPartCuts != nTCuts) Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
2334   fh1EvtMult->Fill(nRecPartCuts);
2335
2336
2337   Int_t nTGen = GetListOfTracks(fTracksGen,fTrackTypeGen);
2338   Int_t nGenPart = 0;
2339   if(nTGen>=0) nGenPart = fTracksGen->GetEntries();
2340   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
2341   if(nGenPart != nTGen) Printf("%s:%d Mismatch selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
2342
2343
2344   //____ analysis, fill histos ___________________________________________________
2345
2346   if(fQAMode){
2347     // loop over tracks
2348     if(fQAMode&1){
2349       for(Int_t it=0; it<nRecPartCuts; ++it){
2350         AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksRecCuts->At(it));
2351         if(part)fQATrackHistosRecCuts->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt() );
2352       }
2353       for(Int_t it=0; it<nGenPart; ++it){
2354         AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksGen->At(it));
2355         if(part)fQATrackHistosGen->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
2356       }
2357     }
2358
2359     // loop over jets
2360     if(fQAMode&2){
2361       for(Int_t ij=0; ij<nRecJets; ++ij){
2362         AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRec->At(ij));
2363         if(jet)fQAJetHistosRec->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2364       }
2365     }
2366   }
2367
2368   if(fQAMode || fFFMode){
2369     for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
2370
2371       AliAODJet* jet = (AliAODJet*)(fJetsRecCuts->At(ij));
2372       if(fQAMode&2) fQAJetHistosRecCuts->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2373
2374       if(fQAMode&2 && ij==0) fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
2375
2376       Double_t ptFractionEmbedded = 0;
2377       AliAODJet* embeddedJet = 0;
2378
2379       if(fBranchEmbeddedJets.Length()){ // find embedded jet
2380
2381         Int_t indexEmbedded = -1;
2382         for(Int_t i=0; i<nEmbeddedJets; i++){
2383           if(iEmbeddedMatchIndex[i] == ij){
2384             indexEmbedded      = i;
2385             ptFractionEmbedded = fEmbeddedPtFraction[i];
2386           }
2387         }
2388
2389         fh1IndexEmbedded->Fill(indexEmbedded);
2390         fh1FractionPtEmbedded->Fill(ptFractionEmbedded);
2391
2392         if(indexEmbedded>-1){
2393
2394           embeddedJet = dynamic_cast<AliAODJet*>(fJetsEmbedded->At(indexEmbedded));
2395           if(!embeddedJet) continue;
2396
2397           Double_t deltaPt = jet->Pt() - embeddedJet->Pt();
2398           Double_t deltaR  = jet->DeltaR((AliVParticle*) (embeddedJet));
2399
2400           fh2DeltaPtVsJetPtEmbedded->Fill(embeddedJet->Pt(),deltaPt);
2401           fh2DeltaPtVsRecJetPtEmbedded->Fill(jet->Pt(),deltaPt);
2402           fh1DeltaREmbedded->Fill(deltaR);
2403         }
2404       }
2405
2406       // get tracks in jet
2407       TList* jettracklist = new TList();
2408       Double_t sumPt      = 0.;
2409       Bool_t isBadJet     = kFALSE;
2410
2411       if(GetFFRadius()<=0){
2412         GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2413       } else {
2414         if(fUseEmbeddedJetAxis){
2415           if(embeddedJet) GetJetTracksPointing(fTracksRecCuts, jettracklist, embeddedJet,
2416                                                GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2417           }
2418         else              GetJetTracksPointing(fTracksRecCuts, jettracklist, jet,
2419                                                GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2420       }
2421
2422       if(GetFFMinNTracks()>0 && jettracklist->GetSize()<=GetFFMinNTracks()) isBadJet = kTRUE;
2423
2424       if(isBadJet) continue;
2425
2426       if(ptFractionEmbedded>=fCutFractionPtEmbedded){ // if no embedding: ptFraction = cutFraction = 0
2427
2428         for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2429
2430           AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
2431           if(!trackVP)continue;
2432
2433           AliAODTrack * aodtrack  = dynamic_cast<AliAODTrack*>(jettracklist->At(it));
2434           if(!aodtrack) continue;
2435
2436           TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
2437
2438           Float_t jetPt   = jet->Pt();
2439           if(fUseEmbeddedJetPt){
2440             if(embeddedJet) jetPt = embeddedJet->Pt();
2441             else jetPt = 0;
2442           }
2443           Float_t trackPt = trackV->Pt();
2444
2445           Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2446
2447           if(fFFMode && (ij==0)) fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
2448           if(fFFMode)            fFFHistosRecCutsInc->FillFF(trackPt, jetPt, incrementJetPt);
2449
2450           if(it==0){ // leading track
2451             if(fFFMode) fFFHistosRecLeadingTrack->FillFF( trackPt, jetPt, kTRUE);
2452           }
2453
2454           delete trackV;
2455         }
2456
2457         // background ff
2458         if(fBckgMode && (ij==0)){
2459           if(fBckgType[0]!=kBckgNone)
2460             FillBckgHistos(fBckgType[0], fTracksRecCuts, fJetsRecCuts, jet,
2461                            fFFBckgHisto0RecCuts,fQABckgHisto0RecCuts, fh1BckgMult0);
2462           if(fBckgType[1]!=kBckgNone)
2463             FillBckgHistos(fBckgType[1], fTracksRecCuts, fJetsRecCuts, jet,
2464                            fFFBckgHisto1RecCuts,fQABckgHisto1RecCuts, fh1BckgMult1);
2465           if(fBckgType[2]!=kBckgNone)
2466             FillBckgHistos(fBckgType[2], fTracksRecCuts, fJetsRecCuts, jet,
2467                            fFFBckgHisto2RecCuts,fQABckgHisto2RecCuts, fh1BckgMult2);
2468           if(fBckgType[3]!=kBckgNone)
2469             FillBckgHistos(fBckgType[3], fTracksRecCuts, fJetsRecCuts, jet,
2470                            fFFBckgHisto3RecCuts,fQABckgHisto3RecCuts, fh1BckgMult3);
2471           if(fBckgType[4]!=kBckgNone)
2472             FillBckgHistos(fBckgType[4], fTracksRecCuts, fJetsRecCuts, jet,
2473                            fFFBckgHisto4RecCuts,fQABckgHisto4RecCuts, fh1BckgMult4);
2474         } // end if(fBckgMode)
2475
2476
2477         if(fJSMode && (ij==0)) FillJetShape(jet, jettracklist, fProNtracksLeadingJet, fProDelRPtSum, fProDelR80pcPt);
2478
2479         delete jettracklist;
2480
2481       } // end: cut embedded ratio
2482     } // end: rec. jets after cuts
2483
2484     // generated jets
2485     for(Int_t ij=0; ij<nGenJets; ++ij){
2486
2487       AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsGen->At(ij));
2488       if(!jet)continue;
2489
2490       if(fQAMode&2) fQAJetHistosGen->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2491
2492       if(fQAMode&2 && (ij==0)) fQAJetHistosGenLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2493
2494       TList* jettracklist = new TList();
2495       Double_t sumPt      = 0.;
2496       Bool_t isBadJet     = kFALSE;
2497
2498       if(GetFFRadius()<=0){
2499         GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2500       } else {
2501         GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2502       }
2503
2504       if(GetFFMinNTracks()>0 && jettracklist->GetSize()<=GetFFMinNTracks()) isBadJet = kTRUE;;
2505
2506       if(isBadJet) continue;
2507
2508       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2509
2510         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
2511         if(!trackVP)continue;
2512         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
2513
2514         Float_t jetPt   = jet->Pt();
2515         Float_t trackPt = trackV->Pt();
2516
2517         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2518
2519         if(fFFMode && (ij==0)) fFFHistosGen->FillFF( trackPt, jetPt, incrementJetPt );
2520         if(fFFMode)            fFFHistosGenInc->FillFF( trackPt, jetPt, incrementJetPt );
2521
2522         if(it==0){ // leading track
2523           if(fFFMode) fFFHistosGenLeadingTrack->FillFF( trackPt, jetPt, kTRUE );
2524         }
2525
2526         delete trackV;
2527       }
2528
2529       if(fBckgMode && (ij==0)){
2530         if(fBckgType[0]!=kBckgNone)
2531           FillBckgHistos(fBckgType[0], fTracksGen, fJetsGen, jet,
2532                          fFFBckgHisto0Gen, fQABckgHisto0Gen);
2533         if(fBckgType[1]!=kBckgNone)
2534           FillBckgHistos(fBckgType[1], fTracksGen, fJetsGen, jet,
2535                          fFFBckgHisto1Gen, fQABckgHisto1Gen);
2536         if(fBckgType[2]!=kBckgNone)
2537           FillBckgHistos(fBckgType[2], fTracksGen, fJetsGen, jet,
2538                          fFFBckgHisto2Gen, fQABckgHisto2Gen);
2539         if(fBckgType[3]!=kBckgNone)
2540           FillBckgHistos(fBckgType[3], fTracksGen, fJetsGen, jet,
2541                          fFFBckgHisto3Gen, fQABckgHisto3Gen);
2542         if(fBckgType[4]!=kBckgNone)
2543           FillBckgHistos(fBckgType[4], fTracksGen, fJetsGen, jet,
2544                          fFFBckgHisto4Gen, fQABckgHisto4Gen);
2545       } // end if(fBckgMode)
2546
2547       if(fJSMode && (ij==0)) FillJetShape(jet, jettracklist, fProNtracksLeadingJetGen, fProDelRPtSumGen, fProDelR80pcPtGen);
2548
2549       delete jettracklist;
2550     }
2551   } // end: QA, FF and intra-jet
2552
2553
2554   // ____ efficiency _______________________________
2555
2556   if(fEffMode && (fJetTypeRecEff != kJetsUndef)){
2557
2558     // arrays holding for each generated particle the reconstructed AOD track index & isPrimary flag, are initialized in AssociateGenRec(...) function
2559     TArrayI indexAODTr;
2560     TArrayS isGenPrim;
2561
2562     // array holding for each reconstructed AOD track generated particle index, initialized in AssociateGenRec(...) function
2563     TArrayI indexMCTr;
2564
2565     // ... and another set for secondaries from strange/non strange mothers (secondary MC tracks are stored in different lists)
2566     TArrayI indexAODTrSecNS;
2567     TArrayS isGenSecNS;
2568     TArrayI indexMCTrSecNS;
2569
2570     TArrayI indexAODTrSecS;
2571     TArrayS isGenSecS;
2572     TArrayI indexMCTrSecS;
2573
2574     Int_t  nTracksAODMCCharged = GetListOfTracks(fTracksAODMCCharged, kTrackAODMCCharged);
2575     if(fDebug>2)Printf("%s:%d selected AODMC tracks: %d ",(char*)__FILE__,__LINE__,nTracksAODMCCharged);
2576
2577     Int_t  nTracksAODMCChargedSecNS = GetListOfTracks(fTracksAODMCChargedSecNS, kTrackAODMCChargedSecNS);
2578     if(fDebug>2)Printf("%s:%d selected AODMC secondary tracks NS: %d ",(char*)__FILE__,__LINE__,nTracksAODMCChargedSecNS);
2579
2580     Int_t  nTracksAODMCChargedSecS = GetListOfTracks(fTracksAODMCChargedSecS, kTrackAODMCChargedSecS);
2581     if(fDebug>2)Printf("%s:%d selected AODMC secondary tracks S: %d ",(char*)__FILE__,__LINE__,nTracksAODMCChargedSecS);
2582
2583     Int_t  nTracksRecQualityCuts = GetListOfTracks(fTracksRecQualityCuts, kTrackAODQualityCuts);
2584     if(fDebug>2)Printf("%s:%d selected rec tracks quality after cuts, full acceptance/pt : %d ",(char*)__FILE__,__LINE__,nTracksRecQualityCuts);
2585
2586     // associate gen and rec tracks, store indices in TArrays
2587     AssociateGenRec(fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,indexMCTr,isGenPrim,fh2PtRecVsGenPrim);
2588     AssociateGenRec(fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,indexMCTrSecNS,isGenSecNS,fh2PtRecVsGenSec);
2589     AssociateGenRec(fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,indexMCTrSecS,isGenSecS,fh2PtRecVsGenSec);
2590
2591     // single track eff
2592     if(fQAMode&1) FillSingleTrackHistosRecGen(fQATrackHistosRecEffGen,fQATrackHistosRecEffRec,fTracksAODMCCharged,indexAODTr,isGenPrim);
2593
2594     // secondaries
2595     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecNS,fTracksAODMCChargedSecNS,indexAODTrSecNS,isGenSecNS);
2596     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecS,fTracksAODMCChargedSecS,indexAODTrSecS,isGenSecS);
2597     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecSsc,fTracksAODMCChargedSecS,indexAODTrSecS,isGenSecS,kTRUE);
2598
2599
2600     // jet track eff
2601     Double_t sumPtGenLeadingJetRecEff = 0;
2602     Double_t sumPtGenLeadingJetSec    = 0;
2603     Double_t sumPtRecLeadingJetRecEff = 0;
2604
2605     for(Int_t ij=0; ij<nRecEffJets; ++ij){ // jet loop
2606
2607       AliAODJet* jet = (AliAODJet*)(fJetsRecEff->At(ij));
2608
2609       Bool_t isBadJetGenPrim = kFALSE;
2610       Bool_t isBadJetGenSec  = kFALSE;
2611       Bool_t isBadJetRec     = kFALSE;
2612
2613
2614       if(ij==0){ // leading jet
2615
2616         // for efficiency: gen tracks from pointing with gen/rec jet
2617         TList* jettracklistGenPrim = new TList();
2618
2619         // if radius<0 -> trackRefs: collect gen tracks in wide radius + fill FF recEff rec histos with tracks contained in track refs
2620         // note : FF recEff gen histos will be somewhat useless in this approach
2621
2622         if(GetFFRadius() >0)
2623           GetJetTracksPointing(fTracksAODMCCharged, jettracklistGenPrim, jet, GetFFRadius(), sumPtGenLeadingJetRecEff, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetGenPrim);
2624         else
2625           GetJetTracksPointing(fTracksAODMCCharged, jettracklistGenPrim, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetRecEff, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetGenPrim);
2626
2627         TList* jettracklistGenSecNS = new TList();
2628         if(GetFFRadius() >0)
2629           GetJetTracksPointing(fTracksAODMCChargedSecNS, jettracklistGenSecNS, jet, GetFFRadius(), sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2630         else
2631           GetJetTracksPointing(fTracksAODMCChargedSecNS, jettracklistGenSecNS, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2632
2633         TList* jettracklistGenSecS = new TList();
2634         if(GetFFRadius() >0)
2635           GetJetTracksPointing(fTracksAODMCChargedSecS, jettracklistGenSecS, jet, GetFFRadius(), sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2636         else
2637           GetJetTracksPointing(fTracksAODMCChargedSecS, jettracklistGenSecS, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2638
2639
2640         // bin efficiency in jet pt bins using rec tracks
2641         TList* jettracklistRec = new TList();
2642         if(GetFFRadius() >0) GetJetTracksPointing(fTracksRecCuts,jettracklistRec, jet, GetFFRadius(), sumPtRecLeadingJetRecEff, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetRec);
2643         else                 GetJetTracksTrackrefs(jettracklistRec, jet, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetRec);
2644
2645
2646         Double_t jetEta   = jet->Eta();
2647         Double_t jetPhi   = TVector2::Phi_0_2pi(jet->Phi());
2648
2649         if(GetFFMinNTracks()>0 && jettracklistGenPrim->GetSize()<=GetFFMinNTracks())   isBadJetGenPrim = kTRUE;
2650         if(GetFFMinNTracks()>0 && jettracklistGenSecNS->GetSize()<=GetFFMinNTracks())  isBadJetGenSec  = kTRUE;
2651         if(GetFFMinNTracks()>0 && jettracklistRec->GetSize()<=GetFFMinNTracks())       isBadJetRec     = kTRUE;
2652
2653         if(isBadJetRec) continue;
2654
2655         if(fQAMode&2) fQAJetHistosRecEffLeading->FillJetQA( jetEta, jetPhi, sumPtGenLeadingJetRecEff );
2656
2657         if(fFFMode){
2658
2659           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosRecEffRec,jet,
2660                                                     jettracklistGenPrim,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim,
2661                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecPrim,fProDelRPtSumRecPrim,fProDelR80pcPtRecPrim);
2662
2663           else                FillJetTrackHistosRec(fFFHistosRecEffRec,jet,
2664                                                     jettracklistGenPrim,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim,
2665                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecPrim,fProDelRPtSumRecPrim,fProDelR80pcPtRecPrim);
2666
2667
2668           // secondaries: use jet pt from primaries
2669           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecNS,jet,
2670                                                     jettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts, indexAODTrSecNS,isGenSecNS,
2671                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecSecNS,fProDelRPtSumRecSecNS);
2672
2673           else                FillJetTrackHistosRec(fFFHistosSecRecNS,jet,
2674                                                     jettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,isGenSecNS,
2675                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecSecNS,fProDelRPtSumRecSecNS);
2676
2677           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecS,jet,
2678                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2679                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecSecS,fProDelRPtSumRecSecS);
2680
2681           else                FillJetTrackHistosRec(fFFHistosSecRecS,jet,
2682                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2683                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecSecS,fProDelRPtSumRecSecS);
2684
2685           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecSsc,jet,
2686                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2687                                                     0,kTRUE,fJSMode,fProNtracksLeadingJetRecSecSsc,fProDelRPtSumRecSecSsc);
2688
2689           else                FillJetTrackHistosRec(fFFHistosSecRecSsc,jet,
2690                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2691                                                     jettracklistRec,kTRUE,fJSMode,fProNtracksLeadingJetRecSecSsc,fProDelRPtSumRecSecSsc);
2692         }
2693
2694         delete jettracklistGenPrim;
2695         delete jettracklistGenSecNS;
2696         delete jettracklistGenSecS;
2697         delete jettracklistRec;
2698
2699
2700         if(fBckgMode && fFFMode){
2701
2702           TList* perpjettracklistGen  = new TList();
2703           TList* perpjettracklistGen1 = new TList();
2704           TList* perpjettracklistGen2 = new TList();
2705
2706           Double_t sumPtGenPerp1 = 0.;
2707           Double_t sumPtGenPerp2 = 0.;
2708           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCCharged, perpjettracklistGen1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerp1);
2709           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCCharged, perpjettracklistGen2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerp2);
2710
2711           perpjettracklistGen->AddAll(perpjettracklistGen1);
2712           perpjettracklistGen->AddAll(perpjettracklistGen2);
2713
2714           TList* perpjettracklistGenSecNS  = new TList();
2715           TList* perpjettracklistGenSecNS1 = new TList();
2716           TList* perpjettracklistGenSecNS2 = new TList();
2717
2718           Double_t sumPtGenPerpNS1 = 0;
2719           Double_t sumPtGenPerpNS2 = 0;
2720           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS1);
2721           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS2);
2722
2723           perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS1);
2724           perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS2);
2725
2726
2727           TList* perpjettracklistGenSecS  = new TList();
2728           TList* perpjettracklistGenSecS1 = new TList();
2729           TList* perpjettracklistGenSecS2 = new TList();
2730
2731           Double_t sumPtGenPerpS1 = 0;
2732           Double_t sumPtGenPerpS2 = 0;
2733           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS1);
2734           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS2);
2735
2736           perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS1);
2737           perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS2);
2738
2739
2740           if(perpjettracklistGen->GetSize() != perpjettracklistGen1->GetSize() + perpjettracklistGen2->GetSize()){
2741             cout<<" ERROR: perpjettracklistGen size "<<perpjettracklistGen->GetSize()<<" perp1 "<<perpjettracklistGen1->GetSize()
2742                 <<" perp2 "<<perpjettracklistGen2->GetSize()<<endl;
2743             exit(0);
2744           }
2745
2746           if(perpjettracklistGenSecNS->GetSize() != perpjettracklistGenSecNS1->GetSize() + perpjettracklistGenSecNS2->GetSize()){
2747             cout<<" ERROR: perpjettracklistGenSecNS size "<<perpjettracklistGenSecNS->GetSize()<<" perp1 "<<perpjettracklistGenSecNS1->GetSize()
2748                 <<" perp2 "<<perpjettracklistGenSecNS2->GetSize()<<endl;
2749             exit(0);
2750           }
2751
2752           if(perpjettracklistGenSecS->GetSize() != perpjettracklistGenSecS1->GetSize() + perpjettracklistGenSecS2->GetSize()){
2753             cout<<" ERROR: perpjettracklistGenSecS size "<<perpjettracklistGenSecS->GetSize()<<" perp1 "<<perpjettracklistGenSecS1->GetSize()
2754                 <<" perp2 "<<perpjettracklistGenSecS2->GetSize()<<endl;
2755             exit(0);
2756           }
2757
2758
2759           FillJetTrackHistosRec(fFFBckgHisto0RecEffRec,jet,
2760                                 perpjettracklistGen,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim);
2761
2762           FillJetTrackHistosRec(fFFBckgHisto0SecRecNS,jet,
2763                                 perpjettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,isGenSecNS);
2764
2765           FillJetTrackHistosRec(fFFBckgHisto0SecRecS,jet,
2766                                 perpjettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS);
2767
2768           FillJetTrackHistosRec(fFFBckgHisto0SecRecSsc,jet,
2769                                 perpjettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,0,kTRUE);
2770
2771
2772           delete perpjettracklistGen;
2773           delete perpjettracklistGen1;
2774           delete perpjettracklistGen2;
2775
2776           delete perpjettracklistGenSecNS;
2777           delete perpjettracklistGenSecNS1;
2778           delete perpjettracklistGenSecNS2;
2779
2780           delete perpjettracklistGenSecS;
2781           delete perpjettracklistGenSecS1;
2782           delete perpjettracklistGenSecS2;
2783
2784         }
2785       }
2786     }
2787   }
2788
2789   //___________________
2790
2791   fTracksRecCuts->Clear();
2792   fTracksGen->Clear();
2793   fTracksAODMCCharged->Clear();
2794   fTracksAODMCChargedSecNS->Clear();
2795   fTracksAODMCChargedSecS->Clear();
2796   fTracksRecQualityCuts->Clear();
2797
2798   fJetsRec->Clear();
2799   fJetsRecCuts->Clear();
2800   fJetsGen->Clear();
2801   fJetsRecEff->Clear();
2802   fJetsEmbedded->Clear();
2803
2804
2805   if(fBckgMode &&
2806      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
2807       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
2808       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)){
2809
2810     fBckgJetsRec->Clear();
2811     fBckgJetsRecCuts->Clear();
2812     fBckgJetsGen->Clear();
2813   }
2814
2815
2816   //Post output data.
2817   PostData(1, fCommonHistList);
2818 }
2819
2820 //______________________________________________________________
2821 void AliAnalysisTaskFragmentationFunction::Terminate(Option_t *)
2822 {
2823   // terminated
2824
2825   if(fDebug > 1) printf("AliAnalysisTaskFragmentationFunction::Terminate() \n");
2826 }
2827
2828 //_________________________________________________________________________________
2829 Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t type)
2830 {
2831   // fill list of tracks selected according to type
2832
2833   if(fDebug > 2) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
2834
2835   if(!list){
2836     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
2837     return -1;
2838   }
2839
2840   if(!fAOD) return -1;
2841
2842   if(!fAOD->GetTracks()) return 0;
2843
2844   if(type==kTrackUndef) return 0;
2845
2846   Int_t iCount = 0;
2847
2848   if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts || type==kTrackAODExtra || type==kTrackAODExtraonly){
2849
2850     TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(fAOD->FindListObject("aodExtraTracks"));
2851     if(!aodExtraTracks)return iCount;
2852     for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
2853       AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
2854       if (!track) continue;
2855
2856       AliAODTrack *tr = dynamic_cast<AliAODTrack*> (track);
2857       if(!tr)continue;
2858
2859       if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts){
2860
2861         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
2862
2863         if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
2864         if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
2865         if(tr->Pt()  < fTrackPtCut) continue;
2866       }
2867
2868       list->Add(tr);
2869       iCount++;
2870     }
2871   }
2872
2873   if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD || type==kTrackAODExtraCuts || type==kTrackAODExtra){
2874
2875     // all rec. tracks, esd filter mask, eta range
2876
2877     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
2878       AliAODTrack *tr = fAOD->GetTrack(it);
2879
2880       if(type == kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAODExtraCuts){
2881
2882         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
2883         if(type == kTrackAODCuts){
2884           if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
2885           if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
2886           if(tr->Pt()  < fTrackPtCut) continue;
2887         }
2888       }
2889       list->Add(tr);
2890       iCount++;
2891     }
2892   }
2893   else if (type==kTrackKineAll || type==kTrackKineCharged || type==kTrackKineChargedAcceptance){
2894     // kine particles, all or rather charged
2895     if(!fMCEvent) return iCount;
2896
2897     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
2898       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
2899
2900       if(type == kTrackKineCharged || type == kTrackKineChargedAcceptance){
2901         if(part->Charge()==0) continue;
2902
2903         if(type == kTrackKineChargedAcceptance &&
2904            (       part->Eta() < fTrackEtaMin
2905                 || part->Eta() > fTrackEtaMax
2906                 || part->Phi() < fTrackPhiMin
2907                 || part->Phi() > fTrackPhiMax
2908                 || part->Pt()  < fTrackPtCut)) continue;
2909       }
2910
2911       list->Add(part);
2912       iCount++;
2913     }
2914   }
2915   else if (type==kTrackAODMCCharged || type==kTrackAODMCAll || type==kTrackAODMCChargedAcceptance || type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS)  {
2916     // MC particles (from AOD), physical primaries, all or rather charged or rather charged within acceptance
2917     if(!fAOD) return -1;
2918
2919     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
2920     if(!tca)return iCount;
2921
2922     for(int it=0; it<tca->GetEntriesFast(); ++it){
2923       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
2924       if(!part)continue;
2925       if(type != kTrackAODMCChargedSecNS && type != kTrackAODMCChargedSecS  && !part->IsPhysicalPrimary())continue;
2926       if((type == kTrackAODMCChargedSecNS || type == kTrackAODMCChargedSecS) && part->IsPhysicalPrimary())continue;
2927
2928       if (type==kTrackAODMCCharged || type==kTrackAODMCChargedAcceptance || type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS){
2929         if(part->Charge()==0) continue;
2930
2931         if(type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS){
2932           Bool_t isFromStrange = kFALSE;
2933           Int_t iMother = part->GetMother();
2934
2935           if(iMother < 0) continue; // throw out PYTHIA stack partons + incoming protons
2936
2937           AliAODMCParticle *partM = dynamic_cast<AliAODMCParticle*>(tca->At(iMother));
2938           if(!partM) continue;
2939
2940           Int_t codeM =  TMath::Abs(partM->GetPdgCode());
2941           Int_t mfl = Int_t (codeM/ TMath::Power(10, Int_t(TMath::Log10(codeM))));
2942           if  (mfl == 3 && codeM != 3) isFromStrange = kTRUE;
2943
2944           if(codeM == 130) isFromStrange = kTRUE; // K0 long
2945           if(part->IsSecondaryFromMaterial()) isFromStrange = kFALSE; // strange resonances from hadronic showers ?
2946
2947           // if(mfl ==3){
2948           //   cout<<" mfl "<<mfl<<" codeM "<<partM->GetPdgCode()<<" code this track "<<part->GetPdgCode()<<endl;
2949           //   cout<<" index this track "<<it<<" index daughter 0 "<<partM->GetDaughter(0)<<" 1 "<<partM->GetDaughter(1)<<endl;
2950           // }
2951
2952           if(type==kTrackAODMCChargedSecNS && isFromStrange) continue;
2953           if(type==kTrackAODMCChargedSecS  && !isFromStrange) continue;
2954         }
2955
2956
2957         if(type==kTrackAODMCChargedAcceptance &&
2958            (     part->Eta() > fTrackEtaMax
2959               || part->Eta() < fTrackEtaMin
2960               || part->Phi() > fTrackPhiMax
2961               || part->Phi() < fTrackPhiMin
2962               || part->Pt()  < fTrackPtCut)) continue;
2963       }
2964
2965       list->Add(part);
2966       iCount++;
2967     }
2968   }
2969
2970   list->Sort();
2971   return iCount;
2972
2973 }
2974 // _______________________________________________________________________________
2975 Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t type)
2976 {
2977   // fill list of jets selected according to type
2978
2979   if(!list){
2980     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
2981     return -1;
2982   }
2983
2984   if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
2985
2986     if(fBranchRecJets.Length()==0){
2987       Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
2988       if(fDebug>1)fAOD->Print();
2989       return 0;
2990     }
2991
2992     TClonesArray *aodRecJets = 0;
2993     if(fBranchRecJets.Length())      aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchRecJets.Data()));
2994     if(!aodRecJets)                  aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchRecJets.Data()));
2995     if(fAODExtension&&!aodRecJets)   aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRecJets.Data()));
2996
2997     if(!aodRecJets){
2998       if(fBranchRecJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecJets.Data());
2999       if(fDebug>1)fAOD->Print();
3000       return 0;
3001     }
3002
3003     // Reorder jet pt and fill new temporary AliAODJet objects
3004     Int_t nRecJets = 0;
3005
3006     for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
3007
3008       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
3009       if(!tmp) continue;
3010
3011       if( tmp->Pt() < fJetPtCut ) continue;
3012       if( type == kJetsRecAcceptance &&
3013           (    tmp->Eta() < fJetEtaMin
3014             || tmp->Eta() > fJetEtaMax
3015             || tmp->Phi() < fJetPhiMin
3016             || tmp->Phi() > fJetPhiMax )) continue;
3017
3018
3019       list->Add(tmp);
3020       nRecJets++;
3021     }
3022
3023     list->Sort();
3024
3025     return nRecJets;
3026   }
3027   else if(type == kJetsKine || type == kJetsKineAcceptance){
3028
3029     // generated jets
3030     Int_t nGenJets = 0;
3031
3032     if(!fMCEvent){
3033       if(fDebug>1) Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3034       return 0;
3035     }
3036
3037     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
3038     AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
3039     AliGenHijingEventHeader*  hijingGenHeader = 0x0;
3040
3041     if(!pythiaGenHeader){
3042       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
3043
3044       if(!hijingGenHeader){
3045          Printf("%s:%d no pythiaGenHeader or hijingGenHeader found", (char*)__FILE__,__LINE__);
3046          return 0;
3047       }else{
3048          TLorentzVector mom[4];
3049          AliAODJet* jet[4];
3050          hijingGenHeader->GetJets(mom[0], mom[1], mom[2], mom[3]);
3051
3052          for(Int_t i=0; i<2; ++i){
3053             if(!mom[i].Pt()) continue;
3054             jet[i] = new AliAODJet(mom[i]);
3055
3056             if( type == kJetsKineAcceptance &&
3057                 (    jet[i]->Eta() < fJetEtaMin
3058                   || jet[i]->Eta() > fJetEtaMax
3059                   || jet[i]->Phi() < fJetPhiMin
3060                   || jet[i]->Phi() > fJetPhiMax )) continue;
3061
3062             list->Add(jet[i]);
3063             nGenJets++;
3064          }
3065          list->Sort();
3066          return nGenJets;
3067       }
3068     }
3069
3070     // fetch the pythia generated jets
3071     for(int ip=0; ip<pythiaGenHeader->NTriggerJets(); ++ip){
3072
3073       Float_t p[4];
3074       AliAODJet *jet = new AliAODJet();
3075       pythiaGenHeader->TriggerJet(ip, p);
3076       jet->SetPxPyPzE(p[0], p[1], p[2], p[3]);
3077
3078       if( type == kJetsKineAcceptance &&
3079           (    jet->Eta() < fJetEtaMin
3080             || jet->Eta() > fJetEtaMax
3081             || jet->Phi() < fJetPhiMin
3082             || jet->Phi() > fJetPhiMax )) continue;
3083
3084         list->Add(jet);
3085         nGenJets++;
3086     }
3087     list->Sort();
3088     return nGenJets;
3089   }
3090   else if(type == kJetsGen || type == kJetsGenAcceptance ){
3091
3092     if(fBranchGenJets.Length()==0){
3093       if(fDebug>1) Printf("%s:%d no gen jet branch specified", (char*)__FILE__,__LINE__);
3094       return 0;
3095     }
3096
3097     TClonesArray *aodGenJets = 0;
3098     if(fBranchGenJets.Length()) aodGenJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchGenJets.Data()));
3099     if(!aodGenJets)             aodGenJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchGenJets.Data()));
3100     if(fAODExtension&&!aodGenJets)   aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGenJets.Data()));
3101
3102     if(!aodGenJets){
3103       if(fDebug>0){
3104         if(fBranchGenJets.Length()) Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data());
3105       }
3106       if(fDebug>1)fAOD->Print();
3107       return 0;
3108     }
3109
3110     Int_t nGenJets = 0;
3111
3112     for(Int_t ig=0; ig<aodGenJets->GetEntries(); ++ig){
3113
3114       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
3115       if(!tmp) continue;
3116
3117       if( tmp->Pt() < fJetPtCut ) continue;
3118       if( type == kJetsGenAcceptance &&
3119           (    tmp->Eta() < fJetEtaMin
3120             || tmp->Eta() > fJetEtaMax
3121             || tmp->Phi() < fJetPhiMin
3122             || tmp->Phi() > fJetPhiMax )) continue;
3123
3124         list->Add(tmp);
3125         nGenJets++;
3126     }
3127     list->Sort();
3128     return nGenJets;
3129   }
3130   else if(type == kJetsEmbedded){ // embedded jets
3131
3132     if(fBranchEmbeddedJets.Length()==0){
3133       Printf("%s:%d no embedded jet branch specified", (char*)__FILE__,__LINE__);
3134       if(fDebug>1)fAOD->Print();
3135       return 0;
3136     }
3137
3138     TClonesArray *aodEmbeddedJets = 0;
3139     if(fBranchEmbeddedJets.Length())      aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchEmbeddedJets.Data()));
3140     if(!aodEmbeddedJets)                  aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchEmbeddedJets.Data()));
3141     if(fAODExtension&&!aodEmbeddedJets)   aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchEmbeddedJets.Data()));
3142
3143     if(!aodEmbeddedJets){
3144       if(fBranchEmbeddedJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchEmbeddedJets.Data());
3145       if(fDebug>1)fAOD->Print();
3146       return 0;
3147     }
3148
3149     // Reorder jet pt and fill new temporary AliAODJet objects
3150     Int_t nEmbeddedJets = 0;
3151
3152     for(Int_t ij=0; ij<aodEmbeddedJets->GetEntries(); ++ij){
3153
3154       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodEmbeddedJets->At(ij));
3155       if(!tmp) continue;
3156
3157       if( tmp->Pt() < fJetPtCut ) continue;
3158       if(    tmp->Eta() < fJetEtaMin
3159           || tmp->Eta() > fJetEtaMax
3160           || tmp->Phi() < fJetPhiMin
3161           || tmp->Phi() > fJetPhiMax ) continue;
3162
3163       list->Add(tmp);
3164       nEmbeddedJets++;
3165     }
3166
3167     list->Sort();
3168
3169     return nEmbeddedJets;
3170   }
3171   else{
3172     if(fDebug>0)Printf("%s:%d no such type %d",(char*)__FILE__,__LINE__,type);
3173     return 0;
3174   }
3175 }
3176
3177 // ___________________________________________________________________________________
3178 Int_t AliAnalysisTaskFragmentationFunction::GetListOfBckgJets(TList *list, Int_t type)
3179 {
3180   // fill list of bgr clusters selected according to type
3181
3182   if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
3183
3184     if(fBranchRecBckgClusters.Length()==0){
3185       Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
3186       if(fDebug>1)fAOD->Print();
3187       return 0;
3188     }
3189
3190     TClonesArray *aodRecJets = 0;
3191     if(fBranchRecBckgClusters.Length()) aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchRecBckgClusters.Data()));
3192     if(!aodRecJets)                     aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchRecBckgClusters.Data()));
3193     if(fAODExtension&&!aodRecJets)      aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRecBckgClusters.Data()));
3194
3195     if(!aodRecJets){
3196       if(fBranchRecBckgClusters.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecBckgClusters.Data());
3197       if(fDebug>1)fAOD->Print();
3198       return 0;
3199     }
3200
3201     // Reorder jet pt and fill new temporary AliAODJet objects
3202     Int_t nRecJets = 0;
3203
3204     for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
3205
3206       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
3207       if(!tmp) continue;
3208
3209       // if( tmp->Pt() < fJetPtCut ) continue; // no pt cut on bckg clusters !
3210       if( type == kJetsRecAcceptance &&
3211           (    tmp->Eta() < fJetEtaMin
3212                || tmp->Eta() > fJetEtaMax
3213                || tmp->Phi() < fJetPhiMin
3214                || tmp->Phi() > fJetPhiMax )) continue;
3215
3216       list->Add(tmp);
3217
3218       nRecJets++;
3219
3220     }
3221
3222     list->Sort();
3223
3224     return nRecJets;
3225   }
3226
3227   //  /*
3228   // MC clusters still Under construction
3229   //  */
3230
3231   return 0;
3232 }
3233
3234 // _________________________________________________________________________________________________________
3235 void AliAnalysisTaskFragmentationFunction::SetProperties(THnSparse* h,const Int_t dim, const char** labels)
3236 {
3237   // Set properties of THnSparse
3238
3239   for(Int_t i=0; i<dim; i++){
3240     h->GetAxis(i)->SetTitle(labels[i]);
3241     h->GetAxis(i)->SetTitleColor(1);
3242   }
3243 }
3244
3245 // __________________________________________________________________________________________
3246 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
3247 {
3248   //Set properties of histos (x and y title)
3249
3250   h->SetXTitle(x);
3251   h->SetYTitle(y);
3252   h->GetXaxis()->SetTitleColor(1);
3253   h->GetYaxis()->SetTitleColor(1);
3254 }
3255
3256 // _________________________________________________________________________________________________________
3257 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y, const char* z)
3258 {
3259   //Set properties of histos (x,y and z title)
3260
3261   h->SetXTitle(x);
3262   h->SetYTitle(y);
3263   h->SetZTitle(z);
3264   h->GetXaxis()->SetTitleColor(1);
3265   h->GetYaxis()->SetTitleColor(1);
3266   h->GetZaxis()->SetTitleColor(1);
3267 }
3268
3269 // ________________________________________________________________________________________________________________________________________________________
3270 void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist, TList* outputlist, const AliAODJet* jet,
3271                                                                 const Double_t& radius, Double_t& sumPt, const Double_t& minPtL, const Double_t& maxPt, Bool_t& isBadPt)
3272 {
3273   // fill list of tracks in cone around jet axis
3274
3275   sumPt = 0;
3276   Bool_t isBadMaxPt = kFALSE;
3277   Bool_t isBadMinPt = kTRUE;
3278
3279   Double_t jetMom[3];
3280   jet->PxPyPz(jetMom);
3281   TVector3 jet3mom(jetMom);
3282
3283   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3284
3285     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3286     if(!track)continue;
3287     Double_t trackMom[3];
3288     track->PxPyPz(trackMom);
3289     TVector3 track3mom(trackMom);
3290
3291     Double_t dR = jet3mom.DeltaR(track3mom);
3292
3293     if(dR<radius){
3294
3295       outputlist->Add(track);
3296
3297       sumPt += track->Pt();
3298
3299       if(maxPt>0  && track->Pt()>maxPt)  isBadMaxPt = kTRUE;
3300       if(minPtL>0 && track->Pt()>minPtL) isBadMinPt = kFALSE;
3301     }
3302   }
3303
3304   isBadPt = kFALSE;
3305   if(minPtL>0 && isBadMinPt) isBadPt = kTRUE;
3306   if(maxPt>0  && isBadMaxPt) isBadPt = kTRUE;
3307
3308   outputlist->Sort();
3309 }
3310
3311 // _________________________________________________________________________________________________________________________________________________________________
3312 void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, const AliAODJet* jet, const Double_t& minPtL,
3313                                                                  const Double_t& maxPt, Bool_t& isBadPt)
3314 {
3315   // list of jet tracks from trackrefs
3316
3317   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
3318
3319   Bool_t isBadMaxPt = kFALSE;
3320   Bool_t isBadMinPt = kTRUE;
3321
3322   for(Int_t itrack=0; itrack<nTracks; itrack++) {
3323
3324     AliVParticle* track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
3325     if(!track){
3326       AliError("expected ref track not found ");
3327       continue;
3328     }
3329
3330     if(track->Pt()  < fTrackPtCut) continue; // track refs may contain low pt cut (bug in AliFastJetInput)
3331     if(maxPt>0 && track->Pt()>maxPt)   isBadMaxPt = kTRUE;
3332     if(minPtL>0 && track->Pt()>minPtL) isBadMinPt = kFALSE;
3333
3334     list->Add(track);
3335   }
3336
3337   isBadPt = kFALSE;
3338   if(minPtL>0 && isBadMinPt) isBadPt = kTRUE;
3339   if(maxPt>0 && isBadMaxPt)  isBadPt = kTRUE;
3340
3341   list->Sort();
3342 }
3343
3344 // _ ________________________________________________________________________________________________________________________________
3345 void  AliAnalysisTaskFragmentationFunction::AssociateGenRec(TList* tracksAODMCCharged,TList* tracksRec, TArrayI& indexAODTr,TArrayI& indexMCTr,
3346                                                             TArrayS& isRefGen,TH2F* fh2PtRecVsGen)
3347 {
3348   // associate generated and reconstructed tracks, fill TArrays of list indices
3349
3350   Int_t nTracksRec  = tracksRec->GetSize();
3351   Int_t nTracksGen  = tracksAODMCCharged->GetSize();
3352   TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
3353
3354
3355   if(!nTracksGen) return;
3356   if(!tca)        return;
3357
3358   // set size
3359   indexAODTr.Set(nTracksGen);
3360   indexMCTr.Set(nTracksRec);
3361   isRefGen.Set(nTracksGen);
3362
3363   indexAODTr.Reset(-1);
3364   indexMCTr.Reset(-1);
3365   isRefGen.Reset(0);
3366
3367   // loop over reconstructed tracks, get generated track
3368
3369   for(Int_t iRec=0; iRec<nTracksRec; iRec++){
3370
3371     AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
3372     if(!rectrack)continue;
3373     Int_t label = TMath::Abs(rectrack->GetLabel());
3374
3375     // find MC track in our list
3376     AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(label));
3377
3378     Int_t listIndex = -1;
3379     if(gentrack) listIndex = tracksAODMCCharged->IndexOf(gentrack);
3380
3381     if(listIndex>=0){
3382
3383       indexAODTr[listIndex] = iRec;
3384       indexMCTr[iRec]       = listIndex;
3385     }
3386   }
3387
3388
3389   // define reference sample of primaries/secondaries (for reconstruction efficiency / contamination)
3390
3391   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3392
3393     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksAODMCCharged->At(iGen));
3394     if(!gentrack)continue;
3395     Int_t pdg = gentrack->GetPdgCode();
3396
3397     // 211 - pi, 2212 - proton, 321 - Kaon, 11 - electron, 13 - muon
3398     if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 ||
3399        TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
3400
3401       isRefGen[iGen] = kTRUE;
3402
3403       Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3404
3405       if(iRec>=0){
3406         Float_t genPt = gentrack->Pt();
3407         AliAODTrack* vt = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
3408         if(vt){
3409           Float_t recPt = vt->Pt();
3410           fh2PtRecVsGen->Fill(genPt,recPt);
3411         }
3412       }
3413     }
3414   }
3415 }
3416
3417 // _____________________________________________________________________________________________________________________________________________
3418 void AliAnalysisTaskFragmentationFunction::FillSingleTrackHistosRecGen(AliFragFuncQATrackHistos* trackQAGen, AliFragFuncQATrackHistos* trackQARec, TList* tracksGen,
3419                                                                        const TArrayI& indexAODTr, const TArrayS& isRefGen, Bool_t scaleStrangeness){
3420
3421   // fill QA for single track reconstruction efficiency
3422
3423   Int_t nTracksGen  = tracksGen->GetSize();
3424
3425   if(!nTracksGen) return;
3426
3427   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3428
3429     if(isRefGen[iGen] != 1) continue; // select primaries
3430
3431     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksGen->At(iGen));
3432     if(!gentrack) continue;
3433     Double_t ptGen  = gentrack->Pt();
3434     Double_t etaGen = gentrack->Eta();
3435     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3436
3437     // apply same acc & pt cuts as for FF
3438
3439     if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
3440     if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
3441     if(ptGen  < fTrackPtCut) continue;
3442
3443     if(trackQAGen) trackQAGen->FillTrackQA(etaGen, phiGen, ptGen);
3444
3445     Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3446
3447     if(iRec>=0 && trackQARec){
3448       if(scaleStrangeness){
3449         //Double_t weight = GetMCStrangenessFactor(ptGen);
3450         Double_t weight = GetMCStrangenessFactorCMS(gentrack);
3451         trackQARec->FillTrackQA(etaGen, phiGen, ptGen, kFALSE, 0, kTRUE, weight);
3452       }
3453       else trackQARec->FillTrackQA(etaGen, phiGen, ptGen);
3454     }
3455   }
3456 }
3457
3458 // ______________________________________________________________________________________________________________________________________________________
3459
3460 void  AliAnalysisTaskFragmentationFunction::FillJetTrackHistosRec(AliFragFuncHistos* ffhistRec, AliAODJet* jet,
3461                                                                   TList* jetTrackList, const TList* tracksGen, const TList* tracksRec, const TArrayI& indexAODTr,
3462                                                                   const TArrayS& isRefGen, TList* jetTrackListTR, Bool_t scaleStrangeness,
3463                                                                   Bool_t fillJS, TProfile* hProNtracksLeadingJet, TProfile** hProDelRPtSum, TProfile* hProDelR80pcPt)
3464 {
3465   // fill objects for jet track reconstruction efficiency or secondaries contamination
3466   // arguments histGen/histRec can be of different type: AliFragFuncHistos*, AliFragFuncIntraJetHistos*, ...
3467   // jetTrackListTR pointer: track refs if not NULL
3468
3469
3470   // ensure proper normalization, even for secondaries
3471   Double_t jetPtRec = jet->Pt();
3472   ffhistRec->FillFF(-1, jetPtRec, kTRUE);
3473
3474   Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
3475   if(nTracksJet == 0) return;
3476
3477   TList* listRecTracks = new TList();
3478   listRecTracks->Clear();
3479
3480   for(Int_t iTr=0; iTr<nTracksJet; iTr++){ // jet tracks loop
3481
3482     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (jetTrackList->At(iTr));
3483     if(!gentrack)continue;
3484     // find jet track in gen tracks list
3485     Int_t iGen = tracksGen->IndexOf(gentrack);
3486
3487     if(iGen<0){
3488       if(fDebug>0) Printf("%s:%d gen jet track not found ",(char*)__FILE__,__LINE__);
3489       continue;
3490     }
3491
3492     if(isRefGen[iGen] != 1) continue; // select primaries
3493
3494     Double_t ptGen  = gentrack->Pt();
3495     Double_t etaGen = gentrack->Eta();
3496     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3497
3498     // gen level acc & pt cuts - skip in case of track refs
3499     if(!jetTrackListTR && (etaGen < fTrackEtaMin || etaGen > fTrackEtaMax)) continue;
3500     if(!jetTrackListTR && (phiGen < fTrackPhiMin || phiGen > fTrackPhiMax)) continue;
3501     if(!jetTrackListTR &&  ptGen  < fTrackPtCut) continue;
3502
3503
3504     Double_t ptRec = -1;
3505
3506     Int_t iRec   = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3507     Bool_t isRec = (iRec>=0) ? kTRUE : kFALSE;
3508
3509     Bool_t isJetTrack = kFALSE;
3510     if(!jetTrackListTR) isJetTrack = kTRUE; // skip trackRefs check for tracks in ideal cone
3511
3512     if(isRec){
3513
3514       AliAODTrack* rectrack = dynamic_cast<AliAODTrack*> (tracksRec->At(iRec));
3515       if(!rectrack) continue;
3516
3517       ptRec = rectrack->Pt();
3518
3519       if(jetTrackListTR){
3520         Int_t iRecTR = jetTrackListTR->IndexOf(rectrack);
3521         if(iRecTR >=0 ) isJetTrack = kTRUE; // rec tracks assigned to jet
3522       }
3523
3524       if(isJetTrack){
3525
3526         Double_t trackPt = ptRec;
3527         Bool_t incrementJetPt = kFALSE;
3528
3529         if(scaleStrangeness){
3530           //Double_t weight = GetMCStrangenessFactor(ptGen);
3531           Double_t weight = GetMCStrangenessFactorCMS(gentrack);
3532
3533           ffhistRec->FillFF( trackPt, jetPtRec, incrementJetPt, 0, kTRUE, weight );
3534         }
3535         else{
3536           ffhistRec->FillFF( trackPt, jetPtRec, incrementJetPt );
3537         }
3538
3539         listRecTracks->Add(rectrack);
3540
3541       }
3542     }
3543   }
3544
3545
3546   if(fillJS) FillJetShape(jet,listRecTracks,hProNtracksLeadingJet, hProDelRPtSum, hProDelR80pcPt,0,0,scaleStrangeness);
3547
3548   delete listRecTracks;
3549
3550 }
3551
3552 // _____________________________________________________________________________________________________________________________________________________________________
3553 void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt)
3554 {
3555   // List of tracks in cone perpendicular to the jet azimuthal direction
3556
3557   Double_t jetMom[3];
3558   jet->PxPyPz(jetMom);
3559
3560   TVector3 jet3mom(jetMom);
3561   // Rotate phi and keep eta unchanged
3562   Double_t etaTilted = jet3mom.Eta();
3563   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
3564   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
3565
3566   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3567
3568     // embedded tracks
3569     if( fUseExtraTracksBgr != 1){
3570       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3571         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3572         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3573       }
3574     }
3575
3576     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3577     if(!track)continue;
3578     Double_t trackMom[3];
3579     track->PxPyPz(trackMom);
3580     TVector3 track3mom(trackMom);
3581
3582     Double_t deta = track3mom.Eta() - etaTilted;
3583     Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
3584     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
3585     Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
3586
3587
3588     if(dR<=radius){
3589       outputlist->Add(track);
3590       sumPt += track->Pt();
3591     }
3592   }
3593
3594 }
3595
3596 // ________________________________________________________________________________________________________________________________________________________
3597 void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxisWindow(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt,Double_t &normFactor)
3598 {
3599   // List of tracks in cone perpendicular to the jet azimuthal direction
3600
3601   Double_t jetMom[3];
3602   jet->PxPyPz(jetMom);
3603
3604   TVector3 jet3mom(jetMom);
3605   // Rotate phi and keep eta unchanged
3606   Double_t etaTilted = jet3mom.Eta();
3607   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
3608   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
3609
3610   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++)
3611     {
3612
3613       // embedded tracks
3614       if( fUseExtraTracksBgr != 1){
3615         if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3616           if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3617           if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3618         }
3619       }
3620
3621       AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3622       if(!track)continue;
3623       Float_t trackEta = track->Eta();
3624       Float_t trackPhi = track->Phi();
3625
3626       if( ( phiTilted-radius >= 0 ) && ( phiTilted+radius <= 2*TMath::Pi()))
3627         {
3628           if((trackPhi<=phiTilted+radius) &&
3629              (trackPhi>=phiTilted-radius) &&
3630              (trackEta<=fTrackEtaMax)&&(trackEta>=fTrackEtaMin)) // 0.9 and - 0.9
3631             {
3632               outputlist->Add(track);
3633               sumPt += track->Pt();
3634             }
3635         }
3636       else if( phiTilted-radius < 0 )
3637         {
3638           if((( trackPhi < phiTilted+radius ) ||
3639               ( trackPhi > 2*TMath::Pi()-(radius-phiTilted) )) &&
3640              (( trackEta <= fTrackEtaMax ) && ( trackEta >= fTrackEtaMin )))
3641             {
3642               outputlist->Add(track);
3643               sumPt += track->Pt();
3644             }
3645         }
3646       else if( phiTilted+radius > 2*TMath::Pi() )
3647         {
3648           if((( trackPhi > phiTilted-radius ) ||
3649               ( trackPhi < phiTilted+radius-2*TMath::Pi() )) &&
3650              (( trackEta <= fTrackEtaMax ) && ( trackEta >= fTrackEtaMin )))
3651             {
3652               outputlist->Add(track);
3653               sumPt += track->Pt();
3654             }
3655         }
3656     }
3657
3658   // Jet area - Temporarily added should be modified with the proper jet area value
3659   Float_t areaJet = CalcJetArea(etaTilted,radius);
3660   Float_t areaTilted = 2*radius*(fTrackEtaMax-fTrackEtaMin);
3661
3662   normFactor = (Float_t) 1. / (areaJet / areaTilted);
3663
3664 }
3665
3666
3667 // ________________________________________________________________________________________________________________________________________________________
3668 void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJets(Int_t nCases, TList* inputlist, TList* outputlist, TList* jetlist, Double_t& sumPt)
3669 {
3670   // List of tracks outside cone around N jet axis
3671   // Particles taken randomly
3672
3673   sumPt = 0;
3674   //  Int_t   nj  = jetlist->GetSize();
3675   Float_t rc  = TMath::Abs(GetFFRadius());
3676   Float_t rcl = GetFFBckgRadius();
3677
3678   // Estimate jet and background areas
3679   Float_t* areaJet = new Float_t[nCases];
3680   memset(areaJet, 0, sizeof(Float_t) * nCases);
3681   Float_t* areaJetLarge = new Float_t[nCases];
3682   memset(areaJetLarge, 0, sizeof(Float_t) * nCases);
3683   Float_t areaFull = (fTrackEtaMax-fTrackEtaMin)*(fTrackPhiMax-fTrackPhiMin);
3684   Float_t areaOut = areaFull;
3685
3686   //estimate jets and background areas
3687   Int_t nOut = 0;
3688   Int_t ijet = 0;
3689   TList* templist = new TList();
3690   TClonesArray *vect3Jet = new TClonesArray("TVector3",nCases);
3691
3692   for(Int_t ij=0; ij<nCases; ++ij)
3693     {
3694       // Get jet information
3695       AliAODJet* jet = dynamic_cast<AliAODJet*>(jetlist->At(ij));
3696       if(!jet)continue;
3697       TVector3 jet3mom;
3698       jet3mom.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
3699       new((*vect3Jet)[ijet]) TVector3((TVector3)jet3mom);
3700       Float_t etaJet = (Float_t)((TVector3*) vect3Jet->At(ij))->Eta();
3701
3702       // Jet area
3703       areaJet[ij] = CalcJetArea(etaJet,rc);
3704
3705       // Area jet larger angle
3706       areaJetLarge[ij] = CalcJetArea(etaJet,rcl);
3707
3708       // Outside jet area
3709       areaOut = areaOut - areaJetLarge[ij];
3710       ijet++;
3711     }
3712
3713   // List of all tracks outside jet areas
3714   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3715
3716     // embedded tracks
3717     if( fUseExtraTracksBgr != 1){
3718       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3719         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3720         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3721       }
3722     }
3723
3724     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3725
3726     if(!track)continue;
3727     Double_t trackMom[3];
3728     track->PxPyPz(trackMom);
3729     TVector3 track3mom(trackMom);
3730
3731     Double_t *dR = new Double_t[nCases];
3732     for(Int_t ij=0; ij<nCases; ij++)
3733       dR[ij] = (Double_t)((TVector3*) vect3Jet->At(ij))->DeltaR(track3mom);
3734
3735     if((nCases==1 && (dR[0]>rcl)) ||
3736        (nCases==2 && (dR[0]>rcl && dR[1]>rcl)) ||
3737        (nCases==3 && (dR[0]>rcl && dR[1]>rcl && dR[2]>rcl)))
3738       {
3739         templist->Add(track);
3740         nOut++;
3741       }
3742     delete [] dR;
3743   }
3744
3745   // Take tracks randomly
3746   Int_t nScaled = (Int_t) (nOut * areaJet[0] / areaOut + 0.5);
3747   TArrayI* ar = new TArrayI(nOut);
3748
3749   for(Int_t init=0; init<nOut; init++)
3750     (*ar)[init] = init;
3751
3752   Int_t *randIndex = new Int_t[nScaled];
3753   for(Int_t init2=0; init2<nScaled; init2++)
3754     randIndex[init2] = -1;
3755
3756   // Select nScaled different random numbers in nOut
3757   for(Int_t i=0; i<nScaled; i++)
3758     {
3759       Int_t* tmpArr = new Int_t[nOut-i];
3760       Int_t temp = fRandom->Integer(nOut-i);
3761       for(Int_t ind = 0; ind< ar->GetSize()-1; ind++)
3762         {
3763           if(ind<temp) tmpArr[ind] = (*ar)[ind];
3764           else tmpArr[ind] = (*ar)[ind+1];
3765         }
3766       randIndex[i] = (*ar)[temp];
3767
3768       ar->Set(nOut-i-1,tmpArr);
3769
3770       delete [] tmpArr;
3771
3772     }
3773
3774   for(Int_t ipart=0; ipart<nScaled; ipart++)
3775     {
3776       AliVParticle* track = (AliVParticle*)(templist->At(randIndex[ipart]));
3777       outputlist->Add(track);
3778       sumPt += track->Pt();
3779     }
3780
3781   outputlist->Sort();
3782
3783   delete vect3Jet;
3784   delete templist;
3785   delete [] areaJetLarge;
3786   delete [] areaJet;
3787   delete ar;
3788   delete [] randIndex;
3789
3790 }
3791
3792 // ________________________________________________________________________________________________________________________________________________________
3793 void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJetsStat(Int_t nCases, TList* inputlist, TList* outputlist, TList* jetlist, Double_t& sumPt, Double_t &normFactor)
3794 {
3795   // List of tracks outside cone around N jet axis
3796   // All particles taken + final scaling factor
3797
3798   sumPt = 0;
3799   Float_t rc  = TMath::Abs(GetFFRadius());
3800   Float_t rcl = GetFFBckgRadius();
3801
3802   // Estimate jet and background areas
3803   Float_t* areaJet = new Float_t[nCases];
3804   memset(areaJet, 0, sizeof(Float_t) * nCases);
3805   Float_t* areaJetLarge = new Float_t[nCases];
3806   memset(areaJetLarge, 0, sizeof(Float_t) * nCases);
3807   Float_t areaFull = (fTrackEtaMax-fTrackEtaMin)*(fTrackPhiMax-fTrackPhiMin);
3808   Float_t areaOut = areaFull;
3809
3810   //estimate jets and background areas
3811   Int_t nOut = 0;
3812   Int_t ijet = 0;
3813   TClonesArray *vect3Jet = new TClonesArray("TVector3",nCases);
3814
3815   for(Int_t ij=0; ij<nCases; ++ij)
3816     {
3817       // Get jet information
3818       AliAODJet* jet = dynamic_cast<AliAODJet*>(jetlist->At(ij));
3819       if(!jet)continue;
3820       TVector3 jet3mom;
3821       jet3mom.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
3822       new((*vect3Jet)[ijet]) TVector3((TVector3)jet3mom);
3823       Float_t etaJet = (Float_t)((TVector3*) vect3Jet->At(ij))->Eta();
3824
3825       // Jet area
3826       areaJet[ij] = CalcJetArea(etaJet,rc);
3827
3828       // Area jet larger angle
3829       areaJetLarge[ij] = CalcJetArea(etaJet,rcl);
3830
3831       // Outside jets area
3832       areaOut = areaOut - areaJetLarge[ij];
3833       ijet++;
3834     }
3835
3836   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3837
3838     // embedded tracks
3839     if( fUseExtraTracksBgr != 1){
3840       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3841         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3842         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3843       }
3844     }
3845
3846     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3847     if(!track)continue;
3848     Double_t trackMom[3];
3849     track->PxPyPz(trackMom);
3850     TVector3 track3mom(trackMom);
3851
3852     Double_t *dR = new Double_t[nCases];
3853     for(Int_t ij=0; ij<nCases; ij++)
3854         dR[ij] = (Double_t)((TVector3*) vect3Jet->At(ij))->DeltaR(track3mom);
3855
3856     if((nCases==0) ||
3857        (nCases==1 && (dR[0]>rcl)) ||
3858        (nCases==2 && (dR[0]>rcl && dR[1]>rcl)) ||
3859        (nCases==3 && (dR[0]>rcl && dR[1]>rcl && dR[2]>rcl)))
3860       {
3861         outputlist->Add(track);
3862         sumPt += track->Pt();
3863         nOut++;
3864       }
3865     delete [] dR;
3866   }
3867
3868   if(nCases==0) areaJet[0] = TMath::Pi()*rc*rc;
3869   normFactor = (Float_t) 1./(areaJet[0] / areaOut);
3870
3871   outputlist->Sort();
3872
3873   delete vect3Jet;
3874   delete [] areaJetLarge;
3875   delete [] areaJet;
3876
3877 }
3878
3879 // ______________________________________________________________________________________________________________________________________________________
3880 Float_t AliAnalysisTaskFragmentationFunction::CalcJetArea(const Float_t etaJet, const Float_t rc) const
3881 {
3882   // calculate area of jet with eta etaJet and radius rc
3883
3884   Float_t detamax = etaJet + rc;
3885   Float_t detamin = etaJet - rc;
3886   Float_t accmax = 0.0; Float_t accmin = 0.0;
3887   if(detamax > fTrackEtaMax){ // sector outside etamax
3888     Float_t h = fTrackEtaMax - etaJet;
3889     accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
3890   }
3891   if(detamin < fTrackEtaMin){ // sector outside etamin
3892     Float_t h = fTrackEtaMax + etaJet;
3893     accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
3894   }
3895   Float_t areaJet = rc*rc*TMath::Pi() - accmax - accmin;
3896
3897   return areaJet;
3898
3899 }
3900
3901 // ___________________________________________________________________________________________________________________________
3902 void AliAnalysisTaskFragmentationFunction::GetClusterTracksOutOf1Jet(AliAODJet* jet, TList* outputlist, Double_t &normFactor)
3903 {
3904   // fill tracks from bckgCluster branch in list,
3905   // for all clusters outside jet cone
3906   // sum up total area of clusters
3907
3908   Double_t rc  = GetFFRadius();
3909   Double_t rcl = GetFFBckgRadius();
3910
3911   Double_t areaTotal   = 0;
3912   Double_t sumPtTotal  = 0;
3913
3914   for(Int_t ij=0; ij<fBckgJetsRec->GetEntries(); ++ij){
3915
3916     AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij)); // not 'recCuts': use all clusters in full eta range
3917
3918     Double_t dR = jet->DeltaR(bgrCluster);
3919
3920     if(dR<rcl) continue;
3921
3922     Double_t clusterPt = bgrCluster->Pt();
3923     Double_t area      = bgrCluster->EffectiveAreaCharged();
3924     areaTotal  += area;
3925     sumPtTotal += clusterPt;
3926
3927     Int_t nTracksJet = bgrCluster->GetRefTracks()->GetEntries();
3928
3929     for(Int_t it = 0; it<nTracksJet; it++){
3930
3931       // embedded tracks - note: using ref tracks here, fBranchRecBckgClusters has to be consistent
3932       if( fUseExtraTracksBgr != 1){
3933         if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (bgrCluster->GetTrack(it))){
3934           if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3935           if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3936         }
3937       }
3938
3939       AliVParticle*   track = dynamic_cast<AliVParticle*>(bgrCluster->GetTrack(it));
3940       if(!track) continue;
3941
3942       Float_t trackPt  = track->Pt();
3943       Float_t trackEta = track->Eta();
3944       Float_t trackPhi = TVector2::Phi_0_2pi(track->Phi());
3945
3946       if(trackEta < fTrackEtaMin || trackEta > fTrackEtaMax) continue;
3947       if(trackPhi < fTrackPhiMin || trackPhi > fTrackPhiMax) continue;
3948       if(trackPt  < fTrackPtCut) continue;
3949
3950       outputlist->Add(track);
3951     }
3952   }
3953
3954   Double_t areaJet = TMath::Pi()*rc*rc;
3955   if(areaTotal) normFactor = (Float_t) 1./(areaJet / areaTotal);
3956
3957   outputlist->Sort();
3958 }
3959
3960 // _______________________________________________________________________________________________________________________
3961 void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputlist, Double_t &normFactor)
3962 {
3963   // fill tracks from bckgCluster branch,
3964   // using cluster with median density (odd number of clusters)
3965   // or picking randomly one of the two closest to median (even number)
3966
3967   normFactor = 0;
3968
3969   Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range
3970
3971   if(nBckgClusters<3) return; // need at least 3 clusters (skipping 2 highest)
3972
3973   Double_t* bgrDensity = new Double_t[nBckgClusters];
3974   Int_t*    indices    = new Int_t[nBckgClusters];
3975
3976   for(Int_t ij=0; ij<nBckgClusters; ++ij){
3977
3978     AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij));
3979     Double_t clusterPt    = bgrCluster->Pt();
3980     Double_t area         = bgrCluster->EffectiveAreaCharged();
3981
3982     Double_t density = 0;
3983     if(area>0) density = clusterPt/area;
3984
3985     bgrDensity[ij] = density;
3986     indices[ij]    = ij;
3987   }
3988
3989   TMath::Sort(nBckgClusters, bgrDensity, indices);
3990
3991   // get median cluster
3992
3993   AliAODJet* medianCluster = 0;
3994   //Double_t   medianDensity = 0;
3995
3996   if(TMath::Odd(nBckgClusters)){
3997
3998     Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters-1))];
3999     medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex));
4000
4001     Double_t clusterPt = medianCluster->Pt();
4002     Double_t area      = medianCluster->EffectiveAreaCharged();
4003
4004     //if(area>0) medianDensity = clusterPt/area;
4005   }
4006   else{
4007
4008     Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters-1)];
4009     Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters)];
4010
4011     AliAODJet* medianCluster1 = (AliAODJet*)(fBckgJetsRec->At(medianIndex1));
4012     AliAODJet* medianCluster2 = (AliAODJet*)(fBckgJetsRec->At(medianIndex2));
4013
4014     Double_t density1 = 0;
4015     Double_t clusterPt1 = medianCluster1->Pt();
4016     Double_t area1      = medianCluster1->EffectiveAreaCharged();
4017     if(area1>0) density1 = clusterPt1/area1;
4018
4019     Double_t density2 = 0;
4020     Double_t clusterPt2 = medianCluster2->Pt();
4021     Double_t area2      = medianCluster2->EffectiveAreaCharged();
4022     if(area2>0) density2 = clusterPt2/area2;
4023
4024     //medianDensity = 0.5*(density1+density2);
4025
4026     medianCluster = ( (fRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 );  // select one randomly to avoid adding areas
4027   }
4028
4029   Int_t nTracksJet = medianCluster->GetRefTracks()->GetEntries();
4030
4031   for(Int_t it = 0; it<nTracksJet; it++){
4032
4033     // embedded tracks - note: using ref tracks here, fBranchRecBckgClusters has to be consistent
4034     if( fUseExtraTracksBgr != 1){
4035       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (medianCluster->GetTrack(it))){
4036         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
4037         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
4038       }
4039     }
4040
4041     AliVParticle* track = dynamic_cast<AliVParticle*>(medianCluster->GetTrack(it));
4042     if(!track) continue;
4043
4044     Float_t trackPt  = track->Pt();
4045     Float_t trackEta = track->Eta();
4046     Float_t trackPhi = TVector2::Phi_0_2pi(track->Phi());
4047
4048     if(trackEta < fTrackEtaMin || trackEta > fTrackEtaMax) continue;
4049     if(trackPhi < fTrackPhiMin || trackPhi > fTrackPhiMax) continue;
4050     if(trackPt  < fTrackPtCut) continue;
4051
4052     outputlist->Add(track);
4053   }
4054
4055   Double_t areaMedian = medianCluster->EffectiveAreaCharged();
4056   Double_t areaJet = TMath::Pi()*GetFFRadius()*GetFFRadius();
4057
4058   if(areaMedian) normFactor = (Float_t) 1./(areaJet / areaMedian);
4059
4060   outputlist->Sort();
4061
4062   delete[] bgrDensity;
4063   delete[] indices;
4064 }
4065
4066 // ______________________________________________________________________________________________________________________________________________________
4067 void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inputtracklist, TList* inputjetlist, AliAODJet* jet,
4068                                                              AliFragFuncHistos* ffbckghistocuts, AliFragFuncQATrackHistos* qabckghistocuts, TH1F* fh1Mult){
4069
4070   // List of tracks outside jets for background study
4071   TList* tracklistout2jets     = new TList();
4072   TList* tracklistout3jets     = new TList();
4073   TList* tracklistout2jetsStat = new TList();
4074   TList* tracklistout3jetsStat = new TList();
4075   Double_t sumPtOut2Jets       = 0.;
4076   Double_t sumPtOut3Jets       = 0.;
4077   Double_t sumPtOut2JetsStat   = 0.;
4078   Double_t sumPtOut3JetsStat   = 0.;
4079   Double_t normFactor2Jets     = 0.;
4080   Double_t normFactor3Jets     = 0.;
4081
4082   Int_t nRecJetsCuts = inputjetlist->GetEntries();
4083
4084   if(nRecJetsCuts>1) {
4085     GetTracksOutOfNJets(2,inputtracklist, tracklistout2jets, inputjetlist, sumPtOut2Jets);
4086     GetTracksOutOfNJetsStat(2,inputtracklist, tracklistout2jetsStat, inputjetlist,sumPtOut2JetsStat, normFactor2Jets);
4087
4088   }
4089   if(nRecJetsCuts>2) {
4090     GetTracksOutOfNJets(3,inputtracklist, tracklistout3jets, inputjetlist, sumPtOut3Jets);
4091     GetTracksOutOfNJetsStat(3,inputtracklist, tracklistout3jetsStat, inputjetlist, sumPtOut3JetsStat, normFactor3Jets);
4092   }
4093
4094   if(type==kBckgOutLJ || type==kBckgOutAJ)
4095     {
4096       TList* tracklistoutleading   = new TList();
4097       Double_t sumPtOutLeading     = 0.;
4098       GetTracksOutOfNJets(1,inputtracklist, tracklistoutleading, inputjetlist, sumPtOutLeading);
4099       if(type==kBckgOutLJ && fh1Mult) fh1Mult->Fill(tracklistoutleading->GetSize());
4100
4101       for(Int_t it=0; it<tracklistoutleading->GetSize(); ++it){
4102
4103         AliVParticle* trackVP   = (AliVParticle*)(tracklistoutleading->At(it));
4104         if(!trackVP) continue;
4105         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4106
4107         Float_t jetPt   = jet->Pt();
4108         Float_t trackPt = trackV->Pt();
4109
4110         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4111
4112         if(type==kBckgOutLJ)
4113           {
4114             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt);
4115
4116             // Fill track QA for background
4117             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4118           }
4119
4120         // All cases included
4121         if(nRecJetsCuts==1 && type==kBckgOutAJ)
4122           {
4123             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4124           }
4125         delete trackV;
4126       }
4127       // Increment jet pt with one entry in case #tracks outside jets = 0
4128       if(tracklistoutleading->GetSize()==0) {
4129          Float_t jetPt = jet->Pt();
4130          Bool_t incrementJetPt = kTRUE;
4131          if(type==kBckgOutLJ)
4132           {
4133             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4134           }
4135         // All cases included
4136         if(nRecJetsCuts==1 && type==kBckgOutAJ)
4137           {
4138             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4139           }
4140       }
4141       delete tracklistoutleading;
4142     }
4143   if(type==kBckgOutLJStat || type==kBckgOutAJStat)
4144     {
4145       TList* tracklistoutleadingStat   = new TList();
4146       Double_t sumPtOutLeadingStat = 0.;
4147       Double_t normFactorLeading   = 0.;
4148
4149       GetTracksOutOfNJetsStat(1,inputtracklist, tracklistoutleadingStat, inputjetlist, sumPtOutLeadingStat, normFactorLeading);
4150       if(type==kBckgOutLJStat && fh1Mult) fh1Mult->Fill(tracklistoutleadingStat->GetSize());
4151
4152       for(Int_t it=0; it<tracklistoutleadingStat->GetSize(); ++it){
4153
4154         AliVParticle* trackVP   = dynamic_cast<AliVParticle*>(tracklistoutleadingStat->At(it));
4155         if(!trackVP) continue;
4156         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4157
4158         Float_t jetPt   = jet->Pt();
4159         Float_t trackPt = trackV->Pt();
4160         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4161
4162         // Stat plots
4163         if(type==kBckgOutLJStat)
4164           {
4165             if(fFFMode)ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorLeading);
4166
4167             // Fill track QA for background
4168             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt); // OB added bgr QA
4169           }
4170
4171         // All cases included
4172         if(nRecJetsCuts==1 && type==kBckgOutAJStat)
4173           {
4174             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorLeading);
4175             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4176
4177           }
4178         delete trackV;
4179       }
4180       // Increment jet pt with one entry in case #tracks outside jets = 0
4181       if(tracklistoutleadingStat->GetSize()==0) {
4182          Float_t jetPt = jet->Pt();
4183          Bool_t incrementJetPt = kTRUE;
4184          if(type==kBckgOutLJStat)
4185           {
4186             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorLeading);
4187           }
4188         // All cases included
4189         if(nRecJetsCuts==1 && type==kBckgOutLJStat)
4190           {
4191             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorLeading);
4192           }
4193       }
4194
4195       delete tracklistoutleadingStat;
4196     }
4197
4198   if(type==kBckgPerp || type==kBckgPerp2 || type==kBckgPerp2Area)
4199     {
4200       Double_t sumPtPerp1 = 0.;
4201       Double_t sumPtPerp2 = 0.;
4202       TList* tracklistperp  = new TList();
4203       TList* tracklistperp1 = new TList();
4204       TList* tracklistperp2 = new TList();
4205
4206       Double_t norm = 1;
4207       if(type == kBckgPerp2)     norm = 2; // in FillFF() scaleFac = 1/norm = 0.5 - account for double area;
4208       if(type == kBckgPerp2Area) norm = 2*TMath::Pi()*TMath::Abs(GetFFRadius())*TMath::Abs(GetFFRadius()) / jet->EffectiveAreaCharged(); // in FillFF() scaleFac = 1/norm;
4209
4210       GetTracksTiltedwrpJetAxis(TMath::Pi()/2., inputtracklist,tracklistperp1,jet,TMath::Abs(GetFFRadius()),sumPtPerp1);
4211       if(type==kBckgPerp2 || type==kBckgPerp2Area) GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2., inputtracklist,tracklistperp2,jet,TMath::Abs(GetFFRadius()),sumPtPerp2);
4212
4213       tracklistperp->AddAll(tracklistperp1);
4214       tracklistperp->AddAll(tracklistperp2);
4215
4216       if(tracklistperp->GetSize() != tracklistperp1->GetSize() + tracklistperp2->GetSize()){
4217         cout<<" ERROR: tracklistperp size "<<tracklistperp->GetSize()<<" perp1 "<<tracklistperp1->GetSize()<<" perp2 "<<tracklistperp2->GetSize()<<endl;
4218         exit(0);
4219       }
4220
4221       if(fh1Mult) fh1Mult->Fill(tracklistperp->GetSize());
4222
4223       for(Int_t it=0; it<tracklistperp->GetSize(); ++it){
4224
4225         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistperp->At(it));
4226         if(!trackVP)continue;
4227         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4228
4229         Float_t jetPt   = jet->Pt();
4230         Float_t trackPt = trackV->Pt();
4231
4232         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4233
4234         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, norm );
4235
4236         // Fill track QA for background
4237         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4238
4239         delete trackV;
4240       }
4241       // Increment jet pt with one entry in case #tracks outside jets = 0
4242       if(tracklistperp->GetSize()==0) {
4243          Float_t jetPt = jet->Pt();
4244          Bool_t incrementJetPt = kTRUE;
4245          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4246       }
4247
4248
4249       if(fJSMode){
4250         // fill for tracklistperp1/2 separately, divide norm by 2
4251         if(type==kBckgPerp){
4252           FillJetShape(jet, tracklistperp, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, TMath::Pi()/2., 0., kFALSE);
4253         }
4254         if(type==kBckgPerp2){
4255           FillJetShape(jet, tracklistperp1, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0,    TMath::Pi()/2., 0., kFALSE);
4256           FillJetShape(jet, tracklistperp2, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, -1*TMath::Pi()/2., 0., kFALSE);
4257         }
4258         if(type==kBckgPerp2Area){ // divide norm by 2: listperp1/2 filled separately
4259           FillJetShape(jet, tracklistperp1, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0,    TMath::Pi()/2., 0.5*norm, kFALSE);
4260           FillJetShape(jet, tracklistperp2, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, -1*TMath::Pi()/2., 0.5*norm, kFALSE);
4261         }
4262       }
4263
4264       delete tracklistperp;
4265       delete tracklistperp1;
4266       delete tracklistperp2;
4267     }
4268
4269  if(type==kBckgASide)
4270     {
4271       Double_t sumPtASide = 0.;
4272       TList* tracklistaside = new TList();
4273       GetTracksTiltedwrpJetAxis(TMath::Pi(),inputtracklist,tracklistaside,jet,TMath::Abs(GetFFRadius()),sumPtASide);
4274       if(fh1Mult) fh1Mult->Fill(tracklistaside->GetSize());
4275
4276       for(Int_t it=0; it<tracklistaside->GetSize(); ++it){
4277
4278         AliVParticle*   trackVP = (AliVParticle*)(tracklistaside->At(it));
4279         if(!trackVP) continue;
4280         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4281
4282         Float_t jetPt   = jet->Pt();
4283         Float_t trackPt = trackV->Pt();
4284
4285         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4286
4287         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4288
4289         // Fill track QA for background
4290         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4291
4292         delete trackV;
4293       }
4294       if(tracklistaside->GetSize()==0) {
4295          Float_t jetPt = jet->Pt();
4296          Bool_t incrementJetPt = kTRUE;
4297          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4298       }
4299
4300       delete tracklistaside;
4301     }
4302
4303   if(type==kBckgASideWindow)
4304     {
4305       Double_t normFactorASide = 0.;
4306       Double_t sumPtASideW = 0.;
4307       TList* tracklistasidew = new TList();
4308       GetTracksTiltedwrpJetAxisWindow(TMath::Pi(),inputtracklist,tracklistasidew,jet,TMath::Abs(GetFFRadius()),sumPtASideW,normFactorASide);
4309       if(fh1Mult) fh1Mult->Fill(tracklistasidew->GetSize());
4310
4311       for(Int_t it=0; it<tracklistasidew->GetSize(); ++it){
4312
4313         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistasidew->At(it));
4314         if(!trackVP) continue;
4315         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4316
4317         Float_t jetPt   = jet->Pt();
4318         Float_t trackPt = trackV->Pt();
4319         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4320
4321         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorASide);
4322
4323         // Fill track QA for background
4324         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt, kFALSE, normFactorASide);
4325
4326         delete trackV;
4327       }
4328       if(tracklistasidew->GetSize()==0) {
4329          Float_t jetPt = jet->Pt();
4330          Bool_t incrementJetPt = kTRUE;
4331          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorASide);
4332       }
4333
4334       delete tracklistasidew;
4335     }
4336
4337   if(type==kBckgPerpWindow)
4338     {
4339       Double_t normFactorPerp = 0.;
4340       Double_t sumPtPerpW = 0.;
4341       TList* tracklistperpw = new TList();
4342       GetTracksTiltedwrpJetAxisWindow(TMath::Pi()/2.,inputtracklist,tracklistperpw,jet,TMath::Abs(GetFFRadius()),sumPtPerpW,normFactorPerp);
4343       if(fh1Mult) fh1Mult->Fill(tracklistperpw->GetSize());
4344
4345       for(Int_t it=0; it<tracklistperpw->GetSize(); ++it){
4346
4347         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistperpw->At(it));
4348         if(!trackVP) continue;
4349         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4350
4351         Float_t jetPt   = jet->Pt();
4352         Float_t trackPt = trackV->Pt();
4353         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4354
4355         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorPerp);
4356
4357         // Fill track QA for background
4358         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt, kFALSE, normFactorPerp);
4359
4360         delete trackV;
4361       }
4362       if(tracklistperpw->GetSize()==0) {
4363          Float_t jetPt = jet->Pt();
4364          Bool_t incrementJetPt = kTRUE;
4365          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorPerp);
4366       }
4367
4368       delete tracklistperpw;
4369     }
4370
4371
4372   if(type==kBckgOut2J || type==kBckgOutAJ)
4373     {
4374       if(type==kBckgOut2J && fh1Mult) fh1Mult->Fill(tracklistout2jets->GetSize());
4375       for(Int_t it=0; it<tracklistout2jets->GetSize(); ++it){
4376
4377         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout2jets->At(it));
4378         if(!trackVP) continue;
4379         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4380
4381         Float_t jetPt   = jet->Pt();
4382         Float_t trackPt = trackV->Pt();
4383
4384         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4385
4386         if(type==kBckgOut2J)
4387           {
4388             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4389             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4390           }
4391
4392         // All cases included
4393         if(nRecJetsCuts==2 && type==kBckgOutAJ)
4394           {
4395             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4396
4397           }
4398         delete trackV;
4399       }
4400       // Increment jet pt with one entry in case #tracks outside jets = 0
4401       if(tracklistout2jets->GetSize()==0) {
4402          Float_t jetPt = jet->Pt();
4403          Bool_t incrementJetPt = kTRUE;
4404          if(type==kBckgOut2J)
4405           {
4406             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4407           }
4408         // All cases included
4409         if(nRecJetsCuts==2 && type==kBckgOutAJ)
4410           {
4411             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4412           }
4413       }
4414     }
4415
4416   if(type==kBckgOut2JStat || type==kBckgOutAJStat)
4417     {
4418       for(Int_t it=0; it<tracklistout2jetsStat->GetSize(); ++it){
4419
4420         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout2jetsStat->At(it));
4421         if(!trackVP) continue;
4422         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4423
4424         Float_t jetPt   = jet->Pt();
4425         Float_t trackPt = trackV->Pt();
4426         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4427
4428         if(type==kBckgOut2JStat)
4429           {
4430             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor2Jets);
4431
4432             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4433           }
4434
4435         // All cases included
4436         if(nRecJetsCuts==2 && type==kBckgOutAJStat)
4437           {
4438             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor2Jets);
4439
4440             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4441           }
4442         delete trackV;
4443       }
4444       // Increment jet pt with one entry in case #tracks outside jets = 0
4445       if(tracklistout2jetsStat->GetSize()==0) {
4446          Float_t jetPt = jet->Pt();
4447          Bool_t incrementJetPt = kTRUE;
4448          if(type==kBckgOut2JStat)
4449            {
4450             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor2Jets);
4451           }
4452         // All cases included
4453         if(nRecJetsCuts==2 && type==kBckgOutAJStat)
4454           {
4455             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor2Jets);
4456           }
4457       }
4458
4459     }
4460
4461   if(type==kBckgOut3J || type==kBckgOutAJ)
4462     {
4463       if(type==kBckgOut3J && fh1Mult) fh1Mult->Fill(tracklistout3jets->GetSize());
4464
4465       for(Int_t it=0; it<tracklistout3jets->GetSize(); ++it){
4466
4467         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout3jets->At(it));
4468         if(!trackVP) continue;
4469         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4470
4471         Float_t jetPt   = jet->Pt();
4472         Float_t trackPt = trackV->Pt();
4473
4474         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4475
4476         if(type==kBckgOut3J)
4477           {
4478             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4479
4480             qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4481           }
4482
4483         // All cases included
4484         if(nRecJetsCuts==3 && type==kBckgOutAJ)
4485           {
4486             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4487
4488           }
4489         delete trackV;
4490       }
4491       // Increment jet pt with one entry in case #tracks outside jets = 0
4492       if(tracklistout3jets->GetSize()==0) {
4493          Float_t jetPt = jet->Pt();
4494          Bool_t incrementJetPt = kTRUE;
4495          if(type==kBckgOut3J)
4496           {
4497             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4498           }
4499         // All cases included
4500         if(nRecJetsCuts==3 && type==kBckgOutAJ)
4501           {
4502             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4503           }
4504       }
4505     }
4506
4507   if(type==kBckgOut3JStat || type==kBckgOutAJStat)
4508     {
4509       for(Int_t it=0; it<tracklistout3jetsStat->GetSize(); ++it){
4510
4511         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout3jetsStat->At(it));
4512         if(!trackVP) continue;
4513         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4514
4515         Float_t jetPt   = jet->Pt();
4516         Float_t trackPt = trackV->Pt();
4517         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4518
4519         if(type==kBckgOut3JStat)
4520           {
4521             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor3Jets);
4522
4523             //if(fQAMode&1)     qabckghistocuts->FillTrackQA( trackEta, TVector2::Phi_0_2pi(trackPhi), trackPt);
4524           }
4525
4526         // All cases included
4527         if(nRecJetsCuts==3 && type==kBckgOutAJStat)
4528           {
4529             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor3Jets );
4530
4531             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4532
4533           }
4534         delete trackV;
4535       }
4536       // Increment jet pt with one entry in case #tracks outside jets = 0
4537       if(tracklistout3jetsStat->GetSize()==0) {
4538          Float_t jetPt = jet->Pt();
4539          Bool_t incrementJetPt = kTRUE;
4540          if(type==kBckgOut3JStat)
4541           {
4542             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor3Jets);
4543           }
4544         // All cases included
4545         if(nRecJetsCuts==3 && type==kBckgOutAJStat)
4546           {
4547             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor3Jets);
4548           }
4549       }
4550
4551     }
4552
4553   if(type==kBckgClustersOutLeading){ // clusters bgr: all tracks in clusters out of leading jet
4554
4555     TList* tracklistClustersOutLeading = new TList();
4556     Double_t normFactorClusters = 0;
4557     Float_t jetPt   = jet->Pt();
4558
4559     GetClusterTracksOutOf1Jet(jet, tracklistClustersOutLeading, normFactorClusters);
4560     if(fh1Mult) fh1Mult->Fill(tracklistClustersOutLeading->GetSize());
4561
4562     for(Int_t it=0; it<tracklistClustersOutLeading->GetSize(); ++it){
4563
4564       AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistClustersOutLeading->At(it));
4565       if(!trackVP) continue;
4566       TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4567
4568       Float_t trackPt = trackVP->Pt();
4569
4570       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4571
4572       if(fFFMode)   ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorClusters );
4573       if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4574
4575       delete trackV;
4576     }
4577
4578     delete tracklistClustersOutLeading;
4579
4580   }
4581
4582   if(type == kBckgClusters){ // clusters bgr: all tracks in 'median cluster'
4583
4584     TList* tracklistClustersMedian = new TList();
4585     Double_t normFactorClusters = 0;
4586     Float_t jetPt = jet->Pt();
4587
4588     GetClusterTracksMedian(tracklistClustersMedian, normFactorClusters);
4589     if(fh1Mult) fh1Mult->Fill(tracklistClustersMedian->GetSize());
4590
4591     for(Int_t it=0; it<tracklistClustersMedian->GetSize(); ++it){
4592
4593       AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistClustersMedian->At(it));
4594       if(!trackVP) continue;
4595       TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4596
4597       Float_t trackPt = trackVP->Pt();
4598
4599       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4600
4601       if(fFFMode)   ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorClusters );
4602       if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4603
4604       delete trackV;
4605     }
4606
4607     delete tracklistClustersMedian;
4608   }
4609
4610   delete tracklistout2jets;
4611   delete tracklistout3jets;
4612   delete tracklistout2jetsStat;
4613   delete tracklistout3jetsStat;
4614 }
4615
4616 //_____________________________________________________________________________________
4617 Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactor(const Double_t& pt)
4618 {
4619   // factor strangeness data/MC as function of pt from UE analysis (Sara Vallero)
4620
4621   Double_t alpha = 1;
4622
4623   if(0.150<pt && pt<0.200) alpha = 3.639;
4624   if(0.200<pt && pt<0.250) alpha = 2.097;
4625   if(0.250<pt && pt<0.300) alpha = 1.930;
4626   if(0.300<pt && pt<0.350) alpha = 1.932;
4627   if(0.350<pt && pt<0.400) alpha = 1.943;
4628   if(0.400<pt && pt<0.450) alpha = 1.993;
4629   if(0.450<pt && pt<0.500) alpha = 1.989;
4630   if(0.500<pt && pt<0.600) alpha = 1.963;
4631   if(0.600<pt && pt<0.700) alpha = 1.917;
4632   if(0.700<pt && pt<0.800) alpha = 1.861;
4633   if(0.800<pt && pt<0.900) alpha = 1.820;
4634   if(0.900<pt && pt<1.000) alpha = 1.741;
4635   if(1.000<pt && pt<1.500) alpha = 0.878;
4636
4637   return alpha;
4638 }
4639
4640 //__________________________________________________________________________________________________
4641 Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactorCMS(AliAODMCParticle* daughter)
4642 {
4643   // strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
4644
4645   TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
4646   if(!tca) return 1;
4647
4648   AliAODMCParticle* currentMother   = daughter;
4649   AliAODMCParticle* currentDaughter = daughter;
4650
4651
4652   // find first primary mother K0s, Lambda or Xi
4653   while(1){
4654
4655     Int_t daughterPDG   = currentDaughter->GetPdgCode();
4656
4657     Int_t motherLabel   = currentDaughter->GetMother();
4658     if(motherLabel >= tca->GetEntriesFast()){ // protection
4659       currentMother = currentDaughter;
4660       break;
4661     }
4662
4663     currentMother     = (AliAODMCParticle*) tca->At(motherLabel);
4664
4665     if(!currentMother){
4666       currentMother = currentDaughter;
4667       break;
4668     }
4669
4670     Int_t motherPDG   = currentMother->GetPdgCode();
4671
4672     // phys. primary found ?
4673     if(currentMother->IsPhysicalPrimary()) break;
4674
4675     if(TMath::Abs(daughterPDG) == 321){ // K+/K- e.g. from phi (ref data not feeddown corrected)
4676       currentMother = currentDaughter; break;
4677     }
4678     if(TMath::Abs(motherPDG) == 310 ){ // K0s e.g. from phi (ref data not feeddown corrected)
4679       break;
4680     }
4681     if(TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122){ // mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
4682       currentMother = currentDaughter; break;
4683     }
4684
4685     currentDaughter = currentMother;
4686   }
4687
4688
4689   Int_t motherPDG   = currentMother->GetPdgCode();
4690   Double_t motherPt = currentMother->Pt();
4691
4692   Double_t fac = 1;
4693
4694   if(TMath::Abs(motherPDG) == 310 || TMath::Abs(motherPDG)==321){ // K0s / K+ / K-
4695
4696     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.768049;
4697     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.732933;
4698     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.650298;
4699     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.571332;
4700     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.518734;
4701     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.492543;
4702     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.482704;
4703     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.488056;
4704     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.488861;
4705     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.492862;
4706     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.504332;
4707     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.501858;
4708     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.512970;
4709     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.524131;
4710     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.539130;
4711     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.554101;
4712     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.560348;
4713     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.568869;
4714     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.583310;
4715     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.604818;
4716     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.632630;
4717     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.710070;
4718     else if(6.00 <= motherPt && motherPt < 8.00) fac = 0.736365;
4719     else if(8.00 <= motherPt && motherPt < 10.00) fac = 0.835865;
4720   }
4721
4722   if(TMath::Abs(motherPDG) == 3122){ // Lambda
4723
4724     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.645162;
4725     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.627431;
4726     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.457136;
4727     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.384369;
4728     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.330597;
4729     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.309571;
4730     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.293620;
4731     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.283709;
4732     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.282047;
4733     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.277261;
4734     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.275772;
4735     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.280726;
4736     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.288540;
4737     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.288315;
4738     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.296619;
4739     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.302993;
4740     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.338121;
4741     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.349800;
4742     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.356802;
4743     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.391202;
4744     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.422573;
4745     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.573815;
4746     else if(6.00 <= motherPt && motherPt < 8.00) fac = 0.786984;
4747     else if(8.00 <= motherPt && motherPt < 10.00) fac = 1.020021;
4748   }
4749
4750   if(TMath::Abs(motherPDG) == 3312 || TMath::Abs(motherPDG) == 3322){ // xi
4751
4752     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.666620;
4753     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.575908;
4754     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.433198;
4755     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.340901;
4756     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.290896;
4757     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.236074;
4758     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.218681;
4759     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.207763;
4760     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.222848;
4761     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.208806;
4762     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.197275;
4763     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.183645;
4764     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.188788;
4765     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.188282;
4766     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.207442;
4767     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.240388;
4768     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.241916;
4769     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.208276;
4770     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.234550;
4771     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.251689;
4772     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.310204;
4773     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.343492;
4774   }
4775
4776   Double_t weight = 1;
4777   if(fac > 0) weight = 1/fac;
4778
4779   return weight;
4780 }
4781
4782 // _________________________________________________________________________________
4783 void  AliAnalysisTaskFragmentationFunction::FillJetShape(AliAODJet* jet, TList* list,
4784                                                          TProfile* hProNtracksLeadingJet, TProfile** hProDelRPtSum, TProfile* hProDelR80pcPt,
4785                                                          Double_t dPhiUE, Double_t normUE, Bool_t scaleStrangeness){
4786
4787   const Int_t   kNbinsR    = 50;
4788   const Float_t kBinWidthR = 0.02;
4789
4790   Int_t nJetTracks = list->GetEntries();
4791
4792   Float_t PtSumA[kNbinsR]     = {0.0};
4793
4794   Float_t *delRA     = new Float_t[nJetTracks];
4795   Float_t *trackPtA  = new Float_t[nJetTracks];
4796   Int_t   *index     = new Int_t[nJetTracks];
4797
4798   for(Int_t i=0; i<nJetTracks; i++){
4799     delRA[i]    = 0;
4800     trackPtA[i] = 0;
4801     index[i]    = 0;
4802   }
4803
4804   Double_t jetMom[3];
4805   jet->PxPyPz(jetMom);
4806   TVector3 jet3mom(jetMom);
4807
4808   if(TMath::Abs(dPhiUE)>0){
4809     Double_t phiTilted = jet3mom.Phi();
4810     phiTilted += dPhiUE;
4811     phiTilted = TVector2::Phi_0_2pi(phiTilted);
4812     jet3mom.SetPhi(phiTilted);
4813   }
4814
4815   Double_t jetPt = jet->Pt();
4816   Double_t sumWeights = 0;
4817
4818   for (Int_t j =0; j<nJetTracks; j++){
4819
4820     AliVParticle* track = dynamic_cast<AliVParticle*>(list->At(j));
4821     if(!track)continue;
4822
4823     Double_t trackMom[3];
4824     track->PxPyPz(trackMom);
4825     TVector3 track3mom(trackMom);
4826
4827     Double_t dR = jet3mom.DeltaR(track3mom);
4828
4829     delRA[j]    = dR;
4830     trackPtA[j] = track->Pt();
4831
4832     Double_t weight = GetMCStrangenessFactor(track->Pt()); // more correctly should be gen pt
4833     sumWeights += weight;
4834
4835     for(Int_t ibin=1; ibin<=kNbinsR; ibin++){
4836       Float_t xlow = kBinWidthR*(ibin-1);
4837       Float_t xup  = kBinWidthR*ibin;
4838       if(xlow <= dR && dR < xup){
4839
4840         if(scaleStrangeness) PtSumA[ibin-1]     += track->Pt()*weight;
4841         else                 PtSumA[ibin-1]     += track->Pt();
4842       }
4843     }
4844   } // track loop
4845
4846   Float_t jetPtMin=0;
4847   Float_t jetPtMax=0;
4848
4849   for(Int_t ibin=0; ibin<kNbinsR; ibin++){
4850     Float_t fR =  kBinWidthR*(ibin+0.5);
4851
4852     for(Int_t k=0; k<5; k++){
4853       if(k==0){jetPtMin=20.0;jetPtMax=30.0;}
4854       if(k==1){jetPtMin=30.0;jetPtMax=40.0;}
4855       if(k==2){jetPtMin=40.0;jetPtMax=60.0;}
4856       if(k==3){jetPtMin=60.0;jetPtMax=80.0;}
4857       if(k==4){jetPtMin=80.0;jetPtMax=100.0;}
4858       if(jetPt>jetPtMin && jetPt<jetPtMax){
4859
4860         hProDelRPtSum[k]->Fill(fR,PtSumA[ibin]);
4861
4862       }
4863     }
4864   }
4865
4866   if(scaleStrangeness) hProNtracksLeadingJet->Fill(jetPt,sumWeights);
4867   else                 hProNtracksLeadingJet->Fill(jetPt,nJetTracks);
4868
4869   if(normUE)           hProNtracksLeadingJet->Fill(jetPt,nJetTracks/normUE);
4870
4871   if(hProDelR80pcPt){
4872
4873     Float_t PtSum = 0;
4874     Float_t delRPtSum80pc = 0;
4875
4876     TMath::Sort(nJetTracks,delRA,index,0);
4877
4878     for(Int_t ii=0; ii<nJetTracks; ii++){
4879
4880       if(scaleStrangeness){
4881         Double_t weight = GetMCStrangenessFactor(trackPtA[index[ii]]); // more correctly should be gen pt
4882         PtSum += weight*trackPtA[index[ii]];
4883       }
4884       else PtSum += trackPtA[index[ii]];
4885
4886
4887       if(PtSum/jetPt >= 0.8000){
4888         delRPtSum80pc = delRA[index[ii]];
4889         break;
4890       }
4891     }
4892     hProDelR80pcPt->Fill(jetPt,delRPtSum80pc);
4893   }
4894
4895   delete[] delRA;
4896   delete[] trackPtA;
4897   delete[] index;
4898 }