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