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 = 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) continue;
2441
2442       if(ptFractionEmbedded>=fCutFractionPtEmbedded){ // if no embedding: ptFraction = cutFraction = 0
2443
2444         for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2445
2446           AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
2447           if(!trackVP)continue;
2448
2449           AliAODTrack * aodtrack  = dynamic_cast<AliAODTrack*>(jettracklist->At(it));
2450           if(!aodtrack) continue;
2451
2452           TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
2453
2454           Float_t jetPt   = jet->Pt();
2455           if(fUseEmbeddedJetPt){
2456             if(embeddedJet) jetPt = embeddedJet->Pt();
2457             else jetPt = 0;
2458           }
2459           Float_t trackPt = trackV->Pt();
2460
2461           Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2462
2463           if(fFFMode && (ij==0)) fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
2464           if(fFFMode)            fFFHistosRecCutsInc->FillFF(trackPt, jetPt, incrementJetPt);
2465
2466           if(it==0){ // leading track
2467             if(fFFMode) fFFHistosRecLeadingTrack->FillFF( trackPt, jetPt, kTRUE);
2468           }
2469
2470           delete trackV;
2471         }
2472
2473         // background ff
2474         if(fBckgMode && (ij==0)){
2475           if(fBckgType[0]!=kBckgNone)
2476             FillBckgHistos(fBckgType[0], fTracksRecCuts, fJetsRecCuts, jet,
2477                            fFFBckgHisto0RecCuts,fQABckgHisto0RecCuts, fh1BckgMult0);
2478           if(fBckgType[1]!=kBckgNone)
2479             FillBckgHistos(fBckgType[1], fTracksRecCuts, fJetsRecCuts, jet,
2480                            fFFBckgHisto1RecCuts,fQABckgHisto1RecCuts, fh1BckgMult1);
2481           if(fBckgType[2]!=kBckgNone)
2482             FillBckgHistos(fBckgType[2], fTracksRecCuts, fJetsRecCuts, jet,
2483                            fFFBckgHisto2RecCuts,fQABckgHisto2RecCuts, fh1BckgMult2);
2484           if(fBckgType[3]!=kBckgNone)
2485             FillBckgHistos(fBckgType[3], fTracksRecCuts, fJetsRecCuts, jet,
2486                            fFFBckgHisto3RecCuts,fQABckgHisto3RecCuts, fh1BckgMult3);
2487           if(fBckgType[4]!=kBckgNone)
2488             FillBckgHistos(fBckgType[4], fTracksRecCuts, fJetsRecCuts, jet,
2489                            fFFBckgHisto4RecCuts,fQABckgHisto4RecCuts, fh1BckgMult4);
2490         } // end if(fBckgMode)
2491
2492
2493         if(fJSMode && (ij==0)) FillJetShape(jet, jettracklist, fProNtracksLeadingJet, fProDelRPtSum, fProDelR80pcPt);
2494
2495         delete jettracklist;
2496
2497       } // end: cut embedded ratio
2498     } // end: rec. jets after cuts
2499
2500     // generated jets
2501     for(Int_t ij=0; ij<nGenJets; ++ij){
2502
2503       AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsGen->At(ij));
2504       if(!jet)continue;
2505
2506       if(fQAMode&2) fQAJetHistosGen->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2507
2508       if(fQAMode&2 && (ij==0)) fQAJetHistosGenLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2509
2510       TList* jettracklist = new TList();
2511       Double_t sumPt      = 0.;
2512       Bool_t isBadJet     = kFALSE;
2513
2514       if(GetFFRadius()<=0){
2515         GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2516       } else {
2517         GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);
2518       }
2519
2520       if(GetFFMinNTracks()>0 && jettracklist->GetSize()<=GetFFMinNTracks()) isBadJet = kTRUE;;
2521
2522       if(isBadJet) continue;
2523
2524       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2525
2526         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
2527         if(!trackVP)continue;
2528         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
2529
2530         Float_t jetPt   = jet->Pt();
2531         Float_t trackPt = trackV->Pt();
2532
2533         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2534
2535         if(fFFMode && (ij==0)) fFFHistosGen->FillFF( trackPt, jetPt, incrementJetPt );
2536         if(fFFMode)            fFFHistosGenInc->FillFF( trackPt, jetPt, incrementJetPt );
2537
2538         if(it==0){ // leading track
2539           if(fFFMode) fFFHistosGenLeadingTrack->FillFF( trackPt, jetPt, kTRUE );
2540         }
2541
2542         delete trackV;
2543       }
2544
2545       if(fBckgMode && (ij==0)){
2546         if(fBckgType[0]!=kBckgNone)
2547           FillBckgHistos(fBckgType[0], fTracksGen, fJetsGen, jet,
2548                          fFFBckgHisto0Gen, fQABckgHisto0Gen);
2549         if(fBckgType[1]!=kBckgNone)
2550           FillBckgHistos(fBckgType[1], fTracksGen, fJetsGen, jet,
2551                          fFFBckgHisto1Gen, fQABckgHisto1Gen);
2552         if(fBckgType[2]!=kBckgNone)
2553           FillBckgHistos(fBckgType[2], fTracksGen, fJetsGen, jet,
2554                          fFFBckgHisto2Gen, fQABckgHisto2Gen);
2555         if(fBckgType[3]!=kBckgNone)
2556           FillBckgHistos(fBckgType[3], fTracksGen, fJetsGen, jet,
2557                          fFFBckgHisto3Gen, fQABckgHisto3Gen);
2558         if(fBckgType[4]!=kBckgNone)
2559           FillBckgHistos(fBckgType[4], fTracksGen, fJetsGen, jet,
2560                          fFFBckgHisto4Gen, fQABckgHisto4Gen);
2561       } // end if(fBckgMode)
2562
2563       if(fJSMode && (ij==0)) FillJetShape(jet, jettracklist, fProNtracksLeadingJetGen, fProDelRPtSumGen, fProDelR80pcPtGen);
2564
2565       delete jettracklist;
2566     }
2567   } // end: QA, FF and intra-jet
2568
2569
2570   // ____ efficiency _______________________________
2571
2572   if(fEffMode && (fJetTypeRecEff != kJetsUndef)){
2573
2574     // arrays holding for each generated particle the reconstructed AOD track index & isPrimary flag, are initialized in AssociateGenRec(...) function
2575     TArrayI indexAODTr;
2576     TArrayS isGenPrim;
2577
2578     // array holding for each reconstructed AOD track generated particle index, initialized in AssociateGenRec(...) function
2579     TArrayI indexMCTr;
2580
2581     // ... and another set for secondaries from strange/non strange mothers (secondary MC tracks are stored in different lists)
2582     TArrayI indexAODTrSecNS;
2583     TArrayS isGenSecNS;
2584     TArrayI indexMCTrSecNS;
2585
2586     TArrayI indexAODTrSecS;
2587     TArrayS isGenSecS;
2588     TArrayI indexMCTrSecS;
2589
2590     Int_t  nTracksAODMCCharged = GetListOfTracks(fTracksAODMCCharged, kTrackAODMCCharged);
2591     if(fDebug>2)Printf("%s:%d selected AODMC tracks: %d ",(char*)__FILE__,__LINE__,nTracksAODMCCharged);
2592
2593     Int_t  nTracksAODMCChargedSecNS = GetListOfTracks(fTracksAODMCChargedSecNS, kTrackAODMCChargedSecNS);
2594     if(fDebug>2)Printf("%s:%d selected AODMC secondary tracks NS: %d ",(char*)__FILE__,__LINE__,nTracksAODMCChargedSecNS);
2595
2596     Int_t  nTracksAODMCChargedSecS = GetListOfTracks(fTracksAODMCChargedSecS, kTrackAODMCChargedSecS);
2597     if(fDebug>2)Printf("%s:%d selected AODMC secondary tracks S: %d ",(char*)__FILE__,__LINE__,nTracksAODMCChargedSecS);
2598
2599     Int_t  nTracksRecQualityCuts = GetListOfTracks(fTracksRecQualityCuts, kTrackAODQualityCuts);
2600     if(fDebug>2)Printf("%s:%d selected rec tracks quality after cuts, full acceptance/pt : %d ",(char*)__FILE__,__LINE__,nTracksRecQualityCuts);
2601
2602     // associate gen and rec tracks, store indices in TArrays
2603     AssociateGenRec(fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,indexMCTr,isGenPrim,fh2PtRecVsGenPrim);
2604     AssociateGenRec(fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,indexMCTrSecNS,isGenSecNS,fh2PtRecVsGenSec);
2605     AssociateGenRec(fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,indexMCTrSecS,isGenSecS,fh2PtRecVsGenSec);
2606
2607     // single track eff
2608     if(fQAMode&1) FillSingleTrackHistosRecGen(fQATrackHistosRecEffGen,fQATrackHistosRecEffRec,fTracksAODMCCharged,indexAODTr,isGenPrim);
2609
2610     // secondaries
2611     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecNS,fTracksAODMCChargedSecNS,indexAODTrSecNS,isGenSecNS);
2612     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecS,fTracksAODMCChargedSecS,indexAODTrSecS,isGenSecS);
2613     if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRecSsc,fTracksAODMCChargedSecS,indexAODTrSecS,isGenSecS,kTRUE);
2614
2615
2616     // jet track eff
2617     Double_t sumPtGenLeadingJetRecEff = 0;
2618     Double_t sumPtGenLeadingJetSec    = 0;
2619     Double_t sumPtRecLeadingJetRecEff = 0;
2620
2621     for(Int_t ij=0; ij<nRecEffJets; ++ij){ // jet loop
2622
2623       AliAODJet* jet = (AliAODJet*)(fJetsRecEff->At(ij));
2624
2625       Bool_t isBadJetGenPrim = kFALSE;
2626       Bool_t isBadJetGenSec  = kFALSE;
2627       Bool_t isBadJetRec     = kFALSE;
2628
2629
2630       if(ij==0){ // leading jet
2631
2632         // for efficiency: gen tracks from pointing with gen/rec jet
2633         TList* jettracklistGenPrim = new TList();
2634
2635         // if radius<0 -> trackRefs: collect gen tracks in wide radius + fill FF recEff rec histos with tracks contained in track refs
2636         // note : FF recEff gen histos will be somewhat useless in this approach
2637
2638         if(GetFFRadius() >0)
2639           GetJetTracksPointing(fTracksAODMCCharged, jettracklistGenPrim, jet, GetFFRadius(), sumPtGenLeadingJetRecEff, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetGenPrim);
2640         else
2641           GetJetTracksPointing(fTracksAODMCCharged, jettracklistGenPrim, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetRecEff, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetGenPrim);
2642
2643         TList* jettracklistGenSecNS = new TList();
2644         if(GetFFRadius() >0)
2645           GetJetTracksPointing(fTracksAODMCChargedSecNS, jettracklistGenSecNS, jet, GetFFRadius(), sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2646         else
2647           GetJetTracksPointing(fTracksAODMCChargedSecNS, jettracklistGenSecNS, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2648
2649         TList* jettracklistGenSecS = new TList();
2650         if(GetFFRadius() >0)
2651           GetJetTracksPointing(fTracksAODMCChargedSecS, jettracklistGenSecS, jet, GetFFRadius(), sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2652         else
2653           GetJetTracksPointing(fTracksAODMCChargedSecS, jettracklistGenSecS, jet, TMath::Abs(GetFFRadius())+0.2, sumPtGenLeadingJetSec, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetGenSec);
2654
2655
2656         // bin efficiency in jet pt bins using rec tracks
2657         TList* jettracklistRec = new TList();
2658         if(GetFFRadius() >0) GetJetTracksPointing(fTracksRecCuts,jettracklistRec, jet, GetFFRadius(), sumPtRecLeadingJetRecEff, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetRec);
2659         else                 GetJetTracksTrackrefs(jettracklistRec, jet, GetFFMinLTrackPt() , GetFFMaxTrackPt(), isBadJetRec);
2660
2661
2662         Double_t jetEta   = jet->Eta();
2663         Double_t jetPhi   = TVector2::Phi_0_2pi(jet->Phi());
2664
2665         if(GetFFMinNTracks()>0 && jettracklistGenPrim->GetSize()<=GetFFMinNTracks())   isBadJetGenPrim = kTRUE;
2666         if(GetFFMinNTracks()>0 && jettracklistGenSecNS->GetSize()<=GetFFMinNTracks())  isBadJetGenSec  = kTRUE;
2667         if(GetFFMinNTracks()>0 && jettracklistRec->GetSize()<=GetFFMinNTracks())       isBadJetRec     = kTRUE;
2668
2669         if(isBadJetRec) continue;
2670
2671         if(fQAMode&2) fQAJetHistosRecEffLeading->FillJetQA( jetEta, jetPhi, sumPtGenLeadingJetRecEff );
2672
2673         if(fFFMode){
2674
2675           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosRecEffRec,jet,
2676                                                     jettracklistGenPrim,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim,
2677                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecPrim,fProDelRPtSumRecPrim,fProDelR80pcPtRecPrim);
2678
2679           else                FillJetTrackHistosRec(fFFHistosRecEffRec,jet,
2680                                                     jettracklistGenPrim,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim,
2681                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecPrim,fProDelRPtSumRecPrim,fProDelR80pcPtRecPrim);
2682
2683
2684           // secondaries: use jet pt from primaries
2685           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecNS,jet,
2686                                                     jettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts, indexAODTrSecNS,isGenSecNS,
2687                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecSecNS,fProDelRPtSumRecSecNS);
2688
2689           else                FillJetTrackHistosRec(fFFHistosSecRecNS,jet,
2690                                                     jettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,isGenSecNS,
2691                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecSecNS,fProDelRPtSumRecSecNS);
2692
2693           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecS,jet,
2694                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2695                                                     0,kFALSE,fJSMode,fProNtracksLeadingJetRecSecS,fProDelRPtSumRecSecS);
2696
2697           else                FillJetTrackHistosRec(fFFHistosSecRecS,jet,
2698                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2699                                                     jettracklistRec,kFALSE,fJSMode,fProNtracksLeadingJetRecSecS,fProDelRPtSumRecSecS);
2700
2701           if(GetFFRadius()>0) FillJetTrackHistosRec(fFFHistosSecRecSsc,jet,
2702                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2703                                                     0,kTRUE,fJSMode,fProNtracksLeadingJetRecSecSsc,fProDelRPtSumRecSecSsc);
2704
2705           else                FillJetTrackHistosRec(fFFHistosSecRecSsc,jet,
2706                                                     jettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,
2707                                                     jettracklistRec,kTRUE,fJSMode,fProNtracksLeadingJetRecSecSsc,fProDelRPtSumRecSecSsc);
2708         }
2709
2710         delete jettracklistGenPrim;
2711         delete jettracklistGenSecNS;
2712         delete jettracklistGenSecS;
2713         delete jettracklistRec;
2714
2715
2716         if(fBckgMode && fFFMode){
2717
2718           TList* perpjettracklistGen  = new TList();
2719           TList* perpjettracklistGen1 = new TList();
2720           TList* perpjettracklistGen2 = new TList();
2721
2722           Double_t sumPtGenPerp  = 0.;
2723           Double_t sumPtGenPerp1 = 0.;
2724           Double_t sumPtGenPerp2 = 0.;
2725           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCCharged, perpjettracklistGen1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerp1);
2726           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCCharged, perpjettracklistGen2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerp2);
2727
2728           perpjettracklistGen->AddAll(perpjettracklistGen1);
2729           perpjettracklistGen->AddAll(perpjettracklistGen2);
2730           sumPtGenPerp = 0.5*(sumPtGenPerp1+sumPtGenPerp2);
2731
2732           TList* perpjettracklistGenSecNS  = new TList();
2733           TList* perpjettracklistGenSecNS1 = new TList();
2734           TList* perpjettracklistGenSecNS2 = new TList();
2735
2736           Double_t sumPtGenPerpNS;
2737           Double_t sumPtGenPerpNS1;
2738           Double_t sumPtGenPerpNS2;
2739           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS1);
2740           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS2);
2741
2742           perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS1);
2743           perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS2);
2744           sumPtGenPerpNS = 0.5*(sumPtGenPerpNS1+sumPtGenPerpNS2);
2745
2746
2747           TList* perpjettracklistGenSecS  = new TList();
2748           TList* perpjettracklistGenSecS1 = new TList();
2749           TList* perpjettracklistGenSecS2 = new TList();
2750
2751           Double_t sumPtGenPerpS;
2752           Double_t sumPtGenPerpS1;
2753           Double_t sumPtGenPerpS2;
2754           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS1);
2755           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS2);
2756
2757           perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS1);
2758           perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS2);
2759           sumPtGenPerpS = 0.5*(sumPtGenPerpS1+sumPtGenPerpS2);
2760
2761
2762           if(perpjettracklistGen->GetSize() != perpjettracklistGen1->GetSize() + perpjettracklistGen2->GetSize()){
2763             cout<<" ERROR: perpjettracklistGen size "<<perpjettracklistGen->GetSize()<<" perp1 "<<perpjettracklistGen1->GetSize()
2764                 <<" perp2 "<<perpjettracklistGen2->GetSize()<<endl;
2765             exit(0);
2766           }
2767
2768           if(perpjettracklistGenSecNS->GetSize() != perpjettracklistGenSecNS1->GetSize() + perpjettracklistGenSecNS2->GetSize()){
2769             cout<<" ERROR: perpjettracklistGenSecNS size "<<perpjettracklistGenSecNS->GetSize()<<" perp1 "<<perpjettracklistGenSecNS1->GetSize()
2770                 <<" perp2 "<<perpjettracklistGenSecNS2->GetSize()<<endl;
2771             exit(0);
2772           }
2773
2774           if(perpjettracklistGenSecS->GetSize() != perpjettracklistGenSecS1->GetSize() + perpjettracklistGenSecS2->GetSize()){
2775             cout<<" ERROR: perpjettracklistGenSecS size "<<perpjettracklistGenSecS->GetSize()<<" perp1 "<<perpjettracklistGenSecS1->GetSize()
2776                 <<" perp2 "<<perpjettracklistGenSecS2->GetSize()<<endl;
2777             exit(0);
2778           }
2779
2780
2781           FillJetTrackHistosRec(fFFBckgHisto0RecEffRec,jet,
2782                                 perpjettracklistGen,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim);
2783
2784           FillJetTrackHistosRec(fFFBckgHisto0SecRecNS,jet,
2785                                 perpjettracklistGenSecNS,fTracksAODMCChargedSecNS,fTracksRecQualityCuts,indexAODTrSecNS,isGenSecNS);
2786
2787           FillJetTrackHistosRec(fFFBckgHisto0SecRecS,jet,
2788                                 perpjettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS);
2789
2790           FillJetTrackHistosRec(fFFBckgHisto0SecRecSsc,jet,
2791                                 perpjettracklistGenSecS,fTracksAODMCChargedSecS,fTracksRecQualityCuts,indexAODTrSecS,isGenSecS,0,kTRUE);
2792
2793
2794           delete perpjettracklistGen;
2795           delete perpjettracklistGen1;
2796           delete perpjettracklistGen2;
2797
2798           delete perpjettracklistGenSecNS;
2799           delete perpjettracklistGenSecNS1;
2800           delete perpjettracklistGenSecNS2;
2801
2802           delete perpjettracklistGenSecS;
2803           delete perpjettracklistGenSecS1;
2804           delete perpjettracklistGenSecS2;
2805
2806         }
2807       }
2808     }
2809   }
2810
2811   //___________________
2812
2813   fTracksRecCuts->Clear();
2814   fTracksGen->Clear();
2815   fTracksAODMCCharged->Clear();
2816   fTracksAODMCChargedSecNS->Clear();
2817   fTracksAODMCChargedSecS->Clear();
2818   fTracksRecQualityCuts->Clear();
2819
2820   fJetsRec->Clear();
2821   fJetsRecCuts->Clear();
2822   fJetsGen->Clear();
2823   fJetsRecEff->Clear();
2824   fJetsEmbedded->Clear();
2825
2826
2827   if(fBckgMode &&
2828      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
2829       fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading ||
2830       fBckgType[3]==kBckgClustersOutLeading || fBckgType[4]==kBckgClustersOutLeading)){
2831
2832     fBckgJetsRec->Clear();
2833     fBckgJetsRecCuts->Clear();
2834     fBckgJetsGen->Clear();
2835   }
2836
2837
2838   //Post output data.
2839   PostData(1, fCommonHistList);
2840 }
2841
2842 //______________________________________________________________
2843 void AliAnalysisTaskFragmentationFunction::Terminate(Option_t *)
2844 {
2845   // terminated
2846
2847   if(fDebug > 1) printf("AliAnalysisTaskFragmentationFunction::Terminate() \n");
2848 }
2849
2850 //_________________________________________________________________________________
2851 Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t type)
2852 {
2853   // fill list of tracks selected according to type
2854
2855   if(fDebug > 2) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
2856
2857   if(!list){
2858     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
2859     return -1;
2860   }
2861
2862   if(!fAOD) return -1;
2863
2864   if(!fAOD->GetTracks()) return 0;
2865
2866   if(type==kTrackUndef) return 0;
2867
2868   Int_t iCount = 0;
2869
2870   if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts || type==kTrackAODExtra || type==kTrackAODExtraonly){
2871
2872     TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(fAOD->FindListObject("aodExtraTracks"));
2873     if(!aodExtraTracks)return iCount;
2874     for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
2875       AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
2876       if (!track) continue;
2877
2878       AliAODTrack *tr = dynamic_cast<AliAODTrack*> (track);
2879       if(!tr)continue;
2880
2881       if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts){
2882
2883         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
2884
2885         if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
2886         if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
2887         if(tr->Pt()  < fTrackPtCut) continue;
2888       }
2889
2890       list->Add(tr);
2891       iCount++;
2892     }
2893   }
2894
2895   if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD || type==kTrackAODExtraCuts || type==kTrackAODExtra){
2896
2897     // all rec. tracks, esd filter mask, eta range
2898
2899     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
2900       AliAODTrack *tr = fAOD->GetTrack(it);
2901
2902       if(type == kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAODExtraCuts){
2903
2904         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
2905         if(type == kTrackAODCuts){
2906           if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
2907           if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
2908           if(tr->Pt()  < fTrackPtCut) continue;
2909         }
2910       }
2911       list->Add(tr);
2912       iCount++;
2913     }
2914   }
2915   else if (type==kTrackKineAll || type==kTrackKineCharged || type==kTrackKineChargedAcceptance){
2916     // kine particles, all or rather charged
2917     if(!fMCEvent) return iCount;
2918
2919     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
2920       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
2921
2922       if(type == kTrackKineCharged || type == kTrackKineChargedAcceptance){
2923         if(part->Charge()==0) continue;
2924
2925         if(type == kTrackKineChargedAcceptance &&
2926            (       part->Eta() < fTrackEtaMin
2927                 || part->Eta() > fTrackEtaMax
2928                 || part->Phi() < fTrackPhiMin
2929                 || part->Phi() > fTrackPhiMax
2930                 || part->Pt()  < fTrackPtCut)) continue;
2931       }
2932
2933       list->Add(part);
2934       iCount++;
2935     }
2936   }
2937   else if (type==kTrackAODMCCharged || type==kTrackAODMCAll || type==kTrackAODMCChargedAcceptance || type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS)  {
2938     // MC particles (from AOD), physical primaries, all or rather charged or rather charged within acceptance
2939     if(!fAOD) return -1;
2940
2941     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
2942     if(!tca)return iCount;
2943
2944     for(int it=0; it<tca->GetEntriesFast(); ++it){
2945       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
2946       if(!part)continue;
2947       if(type != kTrackAODMCChargedSecNS && type != kTrackAODMCChargedSecS  && !part->IsPhysicalPrimary())continue;
2948       if((type == kTrackAODMCChargedSecNS || type == kTrackAODMCChargedSecS) && part->IsPhysicalPrimary())continue;
2949
2950       if (type==kTrackAODMCCharged || type==kTrackAODMCChargedAcceptance || type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS){
2951         if(part->Charge()==0) continue;
2952
2953         if(type==kTrackAODMCChargedSecNS || type==kTrackAODMCChargedSecS){
2954           Bool_t isFromStrange = kFALSE;
2955           Int_t iMother = part->GetMother();
2956
2957           if(iMother < 0) continue; // throw out PYTHIA stack partons + incoming protons
2958
2959           AliAODMCParticle *partM = dynamic_cast<AliAODMCParticle*>(tca->At(iMother));
2960           if(!partM) continue;
2961
2962           Int_t codeM =  TMath::Abs(partM->GetPdgCode());
2963           Int_t mfl = Int_t (codeM/ TMath::Power(10, Int_t(TMath::Log10(codeM))));
2964           if  (mfl == 3 && codeM != 3) isFromStrange = kTRUE;
2965
2966           if(codeM == 130) isFromStrange = kTRUE; // K0 long
2967           if(part->IsSecondaryFromMaterial()) isFromStrange = kFALSE; // strange resonances from hadronic showers ?
2968
2969           // if(mfl ==3){
2970           //   cout<<" mfl "<<mfl<<" codeM "<<partM->GetPdgCode()<<" code this track "<<part->GetPdgCode()<<endl;
2971           //   cout<<" index this track "<<it<<" index daughter 0 "<<partM->GetDaughter(0)<<" 1 "<<partM->GetDaughter(1)<<endl;
2972           // }
2973
2974           if(type==kTrackAODMCChargedSecNS && isFromStrange) continue;
2975           if(type==kTrackAODMCChargedSecS  && !isFromStrange) continue;
2976         }
2977
2978
2979         if(type==kTrackAODMCChargedAcceptance &&
2980            (     part->Eta() > fTrackEtaMax
2981               || part->Eta() < fTrackEtaMin
2982               || part->Phi() > fTrackPhiMax
2983               || part->Phi() < fTrackPhiMin
2984               || part->Pt()  < fTrackPtCut)) continue;
2985       }
2986
2987       list->Add(part);
2988       iCount++;
2989     }
2990   }
2991
2992   list->Sort();
2993   return iCount;
2994
2995 }
2996 // _______________________________________________________________________________
2997 Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t type)
2998 {
2999   // fill list of jets selected according to type
3000
3001   if(!list){
3002     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
3003     return -1;
3004   }
3005
3006   if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
3007
3008     if(fBranchRecJets.Length()==0){
3009       Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
3010       if(fDebug>1)fAOD->Print();
3011       return 0;
3012     }
3013
3014     TClonesArray *aodRecJets = 0;
3015     if(fBranchRecJets.Length())      aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchRecJets.Data()));
3016     if(!aodRecJets)                  aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchRecJets.Data()));
3017     if(fAODExtension&&!aodRecJets)   aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRecJets.Data()));
3018
3019     if(!aodRecJets){
3020       if(fBranchRecJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecJets.Data());
3021       if(fDebug>1)fAOD->Print();
3022       return 0;
3023     }
3024
3025     // Reorder jet pt and fill new temporary AliAODJet objects
3026     Int_t nRecJets = 0;
3027
3028     for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
3029
3030       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
3031       if(!tmp) continue;
3032
3033       if( tmp->Pt() < fJetPtCut ) continue;
3034       if( tmp->EffectiveAreaCharged() < fJetMinArea ) continue;
3035       if( type == kJetsRecAcceptance &&
3036           (    tmp->Eta() < fJetEtaMin
3037             || tmp->Eta() > fJetEtaMax
3038             || tmp->Phi() < fJetPhiMin
3039             || tmp->Phi() > fJetPhiMax )) continue;
3040
3041
3042       list->Add(tmp);
3043       nRecJets++;
3044     }
3045
3046     list->Sort();
3047
3048     return nRecJets;
3049   }
3050   else if(type == kJetsKine || type == kJetsKineAcceptance){
3051
3052     // generated jets
3053     Int_t nGenJets = 0;
3054
3055     if(!fMCEvent){
3056       if(fDebug>1) Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3057       return 0;
3058     }
3059
3060     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
3061     AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
3062     AliGenHijingEventHeader*  hijingGenHeader = 0x0;
3063
3064     if(!pythiaGenHeader){
3065       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
3066
3067       if(!hijingGenHeader){
3068          Printf("%s:%d no pythiaGenHeader or hijingGenHeader found", (char*)__FILE__,__LINE__);
3069          return 0;
3070       }else{
3071          TLorentzVector mom[4];
3072          AliAODJet* jet[4];
3073          hijingGenHeader->GetJets(mom[0], mom[1], mom[2], mom[3]);
3074
3075          for(Int_t i=0; i<2; ++i){
3076             if(!mom[i].Pt()) continue;
3077             jet[i] = new AliAODJet(mom[i]);
3078
3079             if( type == kJetsKineAcceptance &&
3080                 (    jet[i]->Eta() < fJetEtaMin
3081                   || jet[i]->Eta() > fJetEtaMax
3082                   || jet[i]->Phi() < fJetPhiMin
3083                   || jet[i]->Phi() > fJetPhiMax )) continue;
3084
3085             list->Add(jet[i]);
3086             nGenJets++;
3087          }
3088          list->Sort();
3089          return nGenJets;
3090       }
3091     }
3092
3093     // fetch the pythia generated jets
3094     for(int ip=0; ip<pythiaGenHeader->NTriggerJets(); ++ip){
3095
3096       Float_t p[4];
3097       AliAODJet *jet = new AliAODJet();
3098       pythiaGenHeader->TriggerJet(ip, p);
3099       jet->SetPxPyPzE(p[0], p[1], p[2], p[3]);
3100
3101       if( type == kJetsKineAcceptance &&
3102           (    jet->Eta() < fJetEtaMin
3103             || jet->Eta() > fJetEtaMax
3104             || jet->Phi() < fJetPhiMin
3105             || jet->Phi() > fJetPhiMax )) continue;
3106
3107         list->Add(jet);
3108         nGenJets++;
3109     }
3110     list->Sort();
3111     return nGenJets;
3112   }
3113   else if(type == kJetsGen || type == kJetsGenAcceptance ){
3114
3115     if(fBranchGenJets.Length()==0){
3116       if(fDebug>1) Printf("%s:%d no gen jet branch specified", (char*)__FILE__,__LINE__);
3117       return 0;
3118     }
3119
3120     TClonesArray *aodGenJets = 0;
3121     if(fBranchGenJets.Length()) aodGenJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchGenJets.Data()));
3122     if(!aodGenJets)             aodGenJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchGenJets.Data()));
3123     if(fAODExtension&&!aodGenJets)   aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGenJets.Data()));
3124
3125     if(!aodGenJets){
3126       if(fDebug>0){
3127         if(fBranchGenJets.Length()) Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data());
3128       }
3129       if(fDebug>1)fAOD->Print();
3130       return 0;
3131     }
3132
3133     Int_t nGenJets = 0;
3134
3135     for(Int_t ig=0; ig<aodGenJets->GetEntries(); ++ig){
3136
3137       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
3138       if(!tmp) continue;
3139
3140       if( tmp->Pt() < fJetPtCut ) continue;
3141       if( tmp->EffectiveAreaCharged() < fJetMinArea ) continue;
3142       if( type == kJetsGenAcceptance &&
3143           (    tmp->Eta() < fJetEtaMin
3144             || tmp->Eta() > fJetEtaMax
3145             || tmp->Phi() < fJetPhiMin
3146             || tmp->Phi() > fJetPhiMax )) continue;
3147
3148         list->Add(tmp);
3149         nGenJets++;
3150     }
3151     list->Sort();
3152     return nGenJets;
3153   }
3154   else if(type == kJetsEmbedded){ // embedded jets
3155
3156     if(fBranchEmbeddedJets.Length()==0){
3157       Printf("%s:%d no embedded jet branch specified", (char*)__FILE__,__LINE__);
3158       if(fDebug>1)fAOD->Print();
3159       return 0;
3160     }
3161
3162     TClonesArray *aodEmbeddedJets = 0;
3163     if(fBranchEmbeddedJets.Length())      aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchEmbeddedJets.Data()));
3164     if(!aodEmbeddedJets)                  aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchEmbeddedJets.Data()));
3165     if(fAODExtension&&!aodEmbeddedJets)   aodEmbeddedJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchEmbeddedJets.Data()));
3166
3167     if(!aodEmbeddedJets){
3168       if(fBranchEmbeddedJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchEmbeddedJets.Data());
3169       if(fDebug>1)fAOD->Print();
3170       return 0;
3171     }
3172
3173     // Reorder jet pt and fill new temporary AliAODJet objects
3174     Int_t nEmbeddedJets = 0;
3175
3176     for(Int_t ij=0; ij<aodEmbeddedJets->GetEntries(); ++ij){
3177
3178       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodEmbeddedJets->At(ij));
3179       if(!tmp) continue;
3180
3181       if( tmp->Pt() < fJetPtCut ) continue;
3182       if( tmp->EffectiveAreaCharged() < fJetMinArea ) continue;
3183       if(    tmp->Eta() < fJetEtaMin
3184           || tmp->Eta() > fJetEtaMax
3185           || tmp->Phi() < fJetPhiMin
3186           || tmp->Phi() > fJetPhiMax ) continue;
3187
3188       list->Add(tmp);
3189       nEmbeddedJets++;
3190     }
3191
3192     list->Sort();
3193
3194     return nEmbeddedJets;
3195   }
3196   else{
3197     if(fDebug>0)Printf("%s:%d no such type %d",(char*)__FILE__,__LINE__,type);
3198     return 0;
3199   }
3200 }
3201
3202 // ___________________________________________________________________________________
3203 Int_t AliAnalysisTaskFragmentationFunction::GetListOfBckgJets(TList *list, Int_t type)
3204 {
3205   // fill list of bgr clusters selected according to type
3206
3207   if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
3208
3209     if(fBranchRecBckgClusters.Length()==0){
3210       Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
3211       if(fDebug>1)fAOD->Print();
3212       return 0;
3213     }
3214
3215     TClonesArray *aodRecJets = 0;
3216     if(fBranchRecBckgClusters.Length()) aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fBranchRecBckgClusters.Data()));
3217     if(!aodRecJets)                     aodRecJets = dynamic_cast<TClonesArray*>(fAODJets->GetList()->FindObject(fBranchRecBckgClusters.Data()));
3218     if(fAODExtension&&!aodRecJets)      aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRecBckgClusters.Data()));
3219
3220     if(!aodRecJets){
3221       if(fBranchRecBckgClusters.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecBckgClusters.Data());
3222       if(fDebug>1)fAOD->Print();
3223       return 0;
3224     }
3225
3226     // Reorder jet pt and fill new temporary AliAODJet objects
3227     Int_t nRecJets = 0;
3228
3229     for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
3230
3231       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
3232       if(!tmp) continue;
3233
3234       // if( tmp->Pt() < fJetPtCut ) continue; // no pt cut on bckg clusters !
3235       if( type == kJetsRecAcceptance &&
3236           (    tmp->Eta() < fJetEtaMin
3237                || tmp->Eta() > fJetEtaMax
3238                || tmp->Phi() < fJetPhiMin
3239                || tmp->Phi() > fJetPhiMax )) continue;
3240
3241       list->Add(tmp);
3242
3243       nRecJets++;
3244
3245     }
3246
3247     list->Sort();
3248
3249     return nRecJets;
3250   }
3251
3252   //  /*
3253   // MC clusters still Under construction
3254   //  */
3255
3256   return 0;
3257 }
3258
3259 // _________________________________________________________________________________________________________
3260 void AliAnalysisTaskFragmentationFunction::SetProperties(THnSparse* h, Int_t dim, const char** labels)
3261 {
3262   // Set properties of THnSparse
3263
3264   for(Int_t i=0; i<dim; i++){
3265     h->GetAxis(i)->SetTitle(labels[i]);
3266     h->GetAxis(i)->SetTitleColor(1);
3267   }
3268 }
3269
3270 // __________________________________________________________________________________________
3271 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
3272 {
3273   //Set properties of histos (x and y title)
3274
3275   h->SetXTitle(x);
3276   h->SetYTitle(y);
3277   h->GetXaxis()->SetTitleColor(1);
3278   h->GetYaxis()->SetTitleColor(1);
3279 }
3280
3281 // _________________________________________________________________________________________________________
3282 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y, const char* z)
3283 {
3284   //Set properties of histos (x,y and z title)
3285
3286   h->SetXTitle(x);
3287   h->SetYTitle(y);
3288   h->SetZTitle(z);
3289   h->GetXaxis()->SetTitleColor(1);
3290   h->GetYaxis()->SetTitleColor(1);
3291   h->GetZaxis()->SetTitleColor(1);
3292 }
3293
3294 // ________________________________________________________________________________________________________________________________________________________
3295 void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist, TList* outputlist, const AliAODJet* jet,
3296                                                                    Double_t radius, Double_t& sumPt, Double_t minPtL, Double_t maxPt, Bool_t& isBadPt)
3297 {
3298   // fill list of tracks in cone around jet axis
3299
3300   sumPt = 0;
3301   Bool_t isBadMaxPt = kFALSE;
3302   Bool_t isBadMinPt = kTRUE;
3303
3304   Double_t jetMom[3];
3305   jet->PxPyPz(jetMom);
3306   TVector3 jet3mom(jetMom);
3307
3308   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3309
3310     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3311     if(!track)continue;
3312     Double_t trackMom[3];
3313     track->PxPyPz(trackMom);
3314     TVector3 track3mom(trackMom);
3315
3316     Double_t dR = jet3mom.DeltaR(track3mom);
3317
3318     if(dR<radius){
3319
3320       outputlist->Add(track);
3321
3322       sumPt += track->Pt();
3323
3324       if(maxPt>0  && track->Pt()>maxPt)  isBadMaxPt = kTRUE;
3325       if(minPtL>0 && track->Pt()>minPtL) isBadMinPt = kFALSE;
3326     }
3327   }
3328
3329   isBadPt = kFALSE;
3330   if(minPtL>0 && isBadMinPt) isBadPt = kTRUE;
3331   if(maxPt>0  && isBadMaxPt) isBadPt = kTRUE;
3332
3333   outputlist->Sort();
3334 }
3335
3336 // _________________________________________________________________________________________________________________________________________________________________
3337 void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, const AliAODJet* jet, Double_t minPtL, Double_t maxPt, Bool_t& isBadPt)
3338 {
3339   // list of jet tracks from trackrefs
3340
3341   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
3342
3343   Bool_t isBadMaxPt = kFALSE;
3344   Bool_t isBadMinPt = kTRUE;
3345
3346   for(Int_t itrack=0; itrack<nTracks; itrack++) {
3347
3348     AliVParticle* track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
3349     if(!track){
3350       AliError("expected ref track not found ");
3351       continue;
3352     }
3353
3354     if(track->Pt()  < fTrackPtCut) continue; // track refs may contain low pt cut (bug in AliFastJetInput)
3355     if(maxPt>0 && track->Pt()>maxPt)   isBadMaxPt = kTRUE;
3356     if(minPtL>0 && track->Pt()>minPtL) isBadMinPt = kFALSE;
3357
3358     list->Add(track);
3359   }
3360
3361   isBadPt = kFALSE;
3362   if(minPtL>0 && isBadMinPt) isBadPt = kTRUE;
3363   if(maxPt>0 && isBadMaxPt)  isBadPt = kTRUE;
3364
3365   list->Sort();
3366 }
3367
3368 // _ ________________________________________________________________________________________________________________________________
3369 void  AliAnalysisTaskFragmentationFunction::AssociateGenRec(TList* tracksAODMCCharged,TList* tracksRec, TArrayI& indexAODTr,TArrayI& indexMCTr,
3370                                                             TArrayS& isRefGen,TH2F* fh2PtRecVsGen)
3371 {
3372   // associate generated and reconstructed tracks, fill TArrays of list indices
3373
3374   Int_t nTracksRec  = tracksRec->GetSize();
3375   Int_t nTracksGen  = tracksAODMCCharged->GetSize();
3376   TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
3377
3378
3379   if(!nTracksGen) return;
3380   if(!tca)        return;
3381
3382   // set size
3383   indexAODTr.Set(nTracksGen);
3384   indexMCTr.Set(nTracksRec);
3385   isRefGen.Set(nTracksGen);
3386
3387   indexAODTr.Reset(-1);
3388   indexMCTr.Reset(-1);
3389   isRefGen.Reset(0);
3390
3391   // loop over reconstructed tracks, get generated track
3392
3393   for(Int_t iRec=0; iRec<nTracksRec; iRec++){
3394
3395     AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
3396     if(!rectrack)continue;
3397     Int_t label = TMath::Abs(rectrack->GetLabel());
3398
3399     // find MC track in our list
3400     AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(label));
3401
3402     Int_t listIndex = -1;
3403     if(gentrack) listIndex = tracksAODMCCharged->IndexOf(gentrack);
3404
3405     if(listIndex>=0){
3406
3407       indexAODTr[listIndex] = iRec;
3408       indexMCTr[iRec]       = listIndex;
3409     }
3410   }
3411
3412
3413   // define reference sample of primaries/secondaries (for reconstruction efficiency / contamination)
3414
3415   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3416
3417     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksAODMCCharged->At(iGen));
3418     if(!gentrack)continue;
3419     Int_t pdg = gentrack->GetPdgCode();
3420
3421     // 211 - pi, 2212 - proton, 321 - Kaon, 11 - electron, 13 - muon
3422     if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 ||
3423        TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
3424
3425       isRefGen[iGen] = kTRUE;
3426
3427       Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3428
3429       if(iRec>=0){
3430         Float_t genPt = gentrack->Pt();
3431         AliAODTrack* vt = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
3432         if(vt){
3433           Float_t recPt = vt->Pt();
3434           fh2PtRecVsGen->Fill(genPt,recPt);
3435         }
3436       }
3437     }
3438   }
3439 }
3440
3441 // _____________________________________________________________________________________________________________________________________________
3442 void AliAnalysisTaskFragmentationFunction::FillSingleTrackHistosRecGen(AliFragFuncQATrackHistos* trackQAGen, AliFragFuncQATrackHistos* trackQARec, TList* tracksGen,
3443                                                                        const TArrayI& indexAODTr, const TArrayS& isRefGen, Bool_t scaleStrangeness){
3444
3445   // fill QA for single track reconstruction efficiency
3446
3447   Int_t nTracksGen  = tracksGen->GetSize();
3448
3449   if(!nTracksGen) return;
3450
3451   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3452
3453     if(isRefGen[iGen] != 1) continue; // select primaries
3454
3455     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksGen->At(iGen));
3456     if(!gentrack) continue;
3457     Double_t ptGen  = gentrack->Pt();
3458     Double_t etaGen = gentrack->Eta();
3459     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3460
3461     // apply same acc & pt cuts as for FF
3462
3463     if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
3464     if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
3465     if(ptGen  < fTrackPtCut) continue;
3466
3467     if(trackQAGen) trackQAGen->FillTrackQA(etaGen, phiGen, ptGen);
3468
3469     Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3470
3471     if(iRec>=0 && trackQARec){
3472       if(scaleStrangeness){
3473         //Double_t weight = GetMCStrangenessFactor(ptGen);
3474         Double_t weight = GetMCStrangenessFactorCMS(gentrack);
3475         trackQARec->FillTrackQA(etaGen, phiGen, ptGen, kFALSE, 0, kTRUE, weight);
3476       }
3477       else trackQARec->FillTrackQA(etaGen, phiGen, ptGen);
3478     }
3479   }
3480 }
3481
3482 // ______________________________________________________________________________________________________________________________________________________
3483
3484 void  AliAnalysisTaskFragmentationFunction::FillJetTrackHistosRec(AliFragFuncHistos* ffhistRec, AliAODJet* jet,
3485                                                                   TList* jetTrackList, const TList* tracksGen, const TList* tracksRec, const TArrayI& indexAODTr,
3486                                                                   const TArrayS& isRefGen, TList* jetTrackListTR, Bool_t scaleStrangeness,
3487                                                                   Bool_t fillJS, TProfile* hProNtracksLeadingJet, TProfile** hProDelRPtSum, TProfile* hProDelR80pcPt)
3488 {
3489   // fill objects for jet track reconstruction efficiency or secondaries contamination
3490   // arguments histGen/histRec can be of different type: AliFragFuncHistos*, AliFragFuncIntraJetHistos*, ...
3491   // jetTrackListTR pointer: track refs if not NULL
3492
3493
3494   // ensure proper normalization, even for secondaries
3495   Double_t jetPtRec = jet->Pt();
3496   ffhistRec->FillFF(-1, jetPtRec, kTRUE);
3497
3498   Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
3499   if(nTracksJet == 0) return;
3500
3501   TList* listRecTracks = new TList();
3502   listRecTracks->Clear();
3503
3504   for(Int_t iTr=0; iTr<nTracksJet; iTr++){ // jet tracks loop
3505
3506     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (jetTrackList->At(iTr));
3507     if(!gentrack)continue;
3508     // find jet track in gen tracks list
3509     Int_t iGen = tracksGen->IndexOf(gentrack);
3510
3511     if(iGen<0){
3512       if(fDebug>0) Printf("%s:%d gen jet track not found ",(char*)__FILE__,__LINE__);
3513       continue;
3514     }
3515
3516     if(isRefGen[iGen] != 1) continue; // select primaries
3517
3518     Double_t ptGen  = gentrack->Pt();
3519     Double_t etaGen = gentrack->Eta();
3520     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3521
3522     // gen level acc & pt cuts - skip in case of track refs
3523     if(!jetTrackListTR && (etaGen < fTrackEtaMin || etaGen > fTrackEtaMax)) continue;
3524     if(!jetTrackListTR && (phiGen < fTrackPhiMin || phiGen > fTrackPhiMax)) continue;
3525     if(!jetTrackListTR &&  ptGen  < fTrackPtCut) continue;
3526
3527
3528     Double_t ptRec = -1;
3529
3530     Int_t iRec   = indexAODTr[iGen]; // can be -1 if no good reconstructed track
3531     Bool_t isRec = (iRec>=0) ? kTRUE : kFALSE;
3532
3533     Bool_t isJetTrack = kFALSE;
3534     if(!jetTrackListTR) isJetTrack = kTRUE; // skip trackRefs check for tracks in ideal cone
3535
3536     if(isRec){
3537
3538       AliAODTrack* rectrack = dynamic_cast<AliAODTrack*> (tracksRec->At(iRec));
3539       if(!rectrack) continue;
3540
3541       ptRec = rectrack->Pt();
3542
3543       if(jetTrackListTR){
3544         Int_t iRecTR = jetTrackListTR->IndexOf(rectrack);
3545         if(iRecTR >=0 ) isJetTrack = kTRUE; // rec tracks assigned to jet
3546       }
3547
3548       if(isJetTrack){
3549
3550         Double_t trackPt = ptRec;
3551         Bool_t incrementJetPt = kFALSE;
3552
3553         if(scaleStrangeness){
3554           //Double_t weight = GetMCStrangenessFactor(ptGen);
3555           Double_t weight = GetMCStrangenessFactorCMS(gentrack);
3556
3557           ffhistRec->FillFF( trackPt, jetPtRec, incrementJetPt, 0, kTRUE, weight );
3558         }
3559         else{
3560           ffhistRec->FillFF( trackPt, jetPtRec, incrementJetPt );
3561         }
3562
3563         listRecTracks->Add(rectrack);
3564
3565       }
3566     }
3567   }
3568
3569
3570   if(fillJS) FillJetShape(jet,listRecTracks,hProNtracksLeadingJet, hProDelRPtSum, hProDelR80pcPt,0,0,scaleStrangeness);
3571
3572   delete listRecTracks;
3573
3574 }
3575
3576 // _____________________________________________________________________________________________________________________________________________________________________
3577 void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt)
3578 {
3579   // List of tracks in cone perpendicular to the jet azimuthal direction
3580
3581   Double_t jetMom[3];
3582   jet->PxPyPz(jetMom);
3583
3584   TVector3 jet3mom(jetMom);
3585   // Rotate phi and keep eta unchanged
3586   Double_t etaTilted = jet3mom.Eta();
3587   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
3588   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
3589
3590   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3591
3592     // embedded tracks
3593     if( fUseExtraTracksBgr != 1){
3594       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3595         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3596         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3597       }
3598     }
3599
3600     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3601     if(!track)continue;
3602     Double_t trackMom[3];
3603     track->PxPyPz(trackMom);
3604     TVector3 track3mom(trackMom);
3605
3606     Double_t deta = track3mom.Eta() - etaTilted;
3607     Double_t dphi = TMath::Abs(track3mom.Phi() - phiTilted);
3608     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
3609     Double_t dR = TMath::Sqrt(deta * deta + dphi * dphi);
3610
3611
3612     if(dR<=radius){
3613       outputlist->Add(track);
3614       sumPt += track->Pt();
3615     }
3616   }
3617
3618 }
3619
3620 // ________________________________________________________________________________________________________________________________________________________
3621 void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxisWindow(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius,Double_t& sumPt,Double_t &normFactor)
3622 {
3623   // List of tracks in cone perpendicular to the jet azimuthal direction
3624
3625   Double_t jetMom[3];
3626   jet->PxPyPz(jetMom);
3627
3628   TVector3 jet3mom(jetMom);
3629   // Rotate phi and keep eta unchanged
3630   Double_t etaTilted = jet3mom.Eta();
3631   Double_t phiTilted = TVector2::Phi_0_2pi(jet3mom.Phi()) + alpha;
3632   if(phiTilted > 2*TMath::Pi()) phiTilted = phiTilted - 2*TMath::Pi();
3633
3634   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++)
3635     {
3636
3637       // embedded tracks
3638       if( fUseExtraTracksBgr != 1){
3639         if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3640           if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3641           if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3642         }
3643       }
3644
3645       AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3646       if(!track)continue;
3647       Float_t trackEta = track->Eta();
3648       Float_t trackPhi = track->Phi();
3649
3650       if( ( phiTilted-radius >= 0 ) && ( phiTilted+radius <= 2*TMath::Pi()))
3651         {
3652           if((trackPhi<=phiTilted+radius) &&
3653              (trackPhi>=phiTilted-radius) &&
3654              (trackEta<=fTrackEtaMax)&&(trackEta>=fTrackEtaMin)) // 0.9 and - 0.9
3655             {
3656               outputlist->Add(track);
3657               sumPt += track->Pt();
3658             }
3659         }
3660       else if( phiTilted-radius < 0 )
3661         {
3662           if((( trackPhi < phiTilted+radius ) ||
3663               ( trackPhi > 2*TMath::Pi()-(radius-phiTilted) )) &&
3664              (( trackEta <= fTrackEtaMax ) && ( trackEta >= fTrackEtaMin )))
3665             {
3666               outputlist->Add(track);
3667               sumPt += track->Pt();
3668             }
3669         }
3670       else if( phiTilted+radius > 2*TMath::Pi() )
3671         {
3672           if((( trackPhi > phiTilted-radius ) ||
3673               ( trackPhi < phiTilted+radius-2*TMath::Pi() )) &&
3674              (( trackEta <= fTrackEtaMax ) && ( trackEta >= fTrackEtaMin )))
3675             {
3676               outputlist->Add(track);
3677               sumPt += track->Pt();
3678             }
3679         }
3680     }
3681
3682   // Jet area - Temporarily added should be modified with the proper jet area value
3683   Float_t areaJet = CalcJetArea(etaTilted,radius);
3684   Float_t areaTilted = 2*radius*(fTrackEtaMax-fTrackEtaMin);
3685
3686   normFactor = (Float_t) 1. / (areaJet / areaTilted);
3687
3688 }
3689
3690
3691 // ________________________________________________________________________________________________________________________________________________________
3692 void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJets(Int_t nCases, TList* inputlist, TList* outputlist, TList* jetlist, Double_t& sumPt)
3693 {
3694   // List of tracks outside cone around N jet axis
3695   // Particles taken randomly
3696
3697   sumPt = 0;
3698   //  Int_t   nj  = jetlist->GetSize();
3699   Float_t rc  = TMath::Abs(GetFFRadius());
3700   Float_t rcl = GetFFBckgRadius();
3701
3702   // Estimate jet and background areas
3703   Float_t* areaJet = new Float_t[nCases];
3704   memset(areaJet, 0, sizeof(Float_t) * nCases);
3705   Float_t* areaJetLarge = new Float_t[nCases];
3706   memset(areaJetLarge, 0, sizeof(Float_t) * nCases);
3707   Float_t areaFull = (fTrackEtaMax-fTrackEtaMin)*(fTrackPhiMax-fTrackPhiMin);
3708   Float_t areaOut = areaFull;
3709
3710   //estimate jets and background areas
3711   Int_t nOut = 0;
3712   Int_t ijet = 0;
3713   TList* templist = new TList();
3714   TClonesArray *vect3Jet = new TClonesArray("TVector3",nCases);
3715
3716   for(Int_t ij=0; ij<nCases; ++ij)
3717     {
3718       // Get jet information
3719       AliAODJet* jet = dynamic_cast<AliAODJet*>(jetlist->At(ij));
3720       if(!jet)continue;
3721       TVector3 jet3mom;
3722       jet3mom.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
3723       new((*vect3Jet)[ijet]) TVector3((TVector3)jet3mom);
3724       Float_t etaJet = (Float_t)((TVector3*) vect3Jet->At(ij))->Eta();
3725
3726       // Jet area
3727       areaJet[ij] = CalcJetArea(etaJet,rc);
3728
3729       // Area jet larger angle
3730       areaJetLarge[ij] = CalcJetArea(etaJet,rcl);
3731
3732       // Outside jet area
3733       areaOut = areaOut - areaJetLarge[ij];
3734       ijet++;
3735     }
3736
3737   // List of all tracks outside jet areas
3738   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3739
3740     // embedded tracks
3741     if( fUseExtraTracksBgr != 1){
3742       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3743         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3744         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3745       }
3746     }
3747
3748     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3749
3750     if(!track)continue;
3751     Double_t trackMom[3];
3752     track->PxPyPz(trackMom);
3753     TVector3 track3mom(trackMom);
3754
3755     Double_t *dR = new Double_t[nCases];
3756     for(Int_t ij=0; ij<nCases; ij++)
3757       dR[ij] = (Double_t)((TVector3*) vect3Jet->At(ij))->DeltaR(track3mom);
3758
3759     if((nCases==1 && (dR[0]>rcl)) ||
3760        (nCases==2 && (dR[0]>rcl && dR[1]>rcl)) ||
3761        (nCases==3 && (dR[0]>rcl && dR[1]>rcl && dR[2]>rcl)))
3762       {
3763         templist->Add(track);
3764         nOut++;
3765       }
3766     delete [] dR;
3767   }
3768
3769   // Take tracks randomly
3770   Int_t nScaled = (Int_t) (nOut * areaJet[0] / areaOut + 0.5);
3771   TArrayI* ar = new TArrayI(nOut);
3772
3773   for(Int_t init=0; init<nOut; init++)
3774     (*ar)[init] = init;
3775
3776   Int_t *randIndex = new Int_t[nScaled];
3777   for(Int_t init2=0; init2<nScaled; init2++)
3778     randIndex[init2] = -1;
3779
3780   // Select nScaled different random numbers in nOut
3781   for(Int_t i=0; i<nScaled; i++)
3782     {
3783       Int_t* tmpArr = new Int_t[nOut-i];
3784       Int_t temp = fRandom->Integer(nOut-i);
3785       for(Int_t ind = 0; ind< ar->GetSize()-1; ind++)
3786         {
3787           if(ind<temp) tmpArr[ind] = (*ar)[ind];
3788           else tmpArr[ind] = (*ar)[ind+1];
3789         }
3790       randIndex[i] = (*ar)[temp];
3791
3792       ar->Set(nOut-i-1,tmpArr);
3793
3794       delete [] tmpArr;
3795
3796     }
3797
3798   for(Int_t ipart=0; ipart<nScaled; ipart++)
3799     {
3800       AliVParticle* track = (AliVParticle*)(templist->At(randIndex[ipart]));
3801       outputlist->Add(track);
3802       sumPt += track->Pt();
3803     }
3804
3805   outputlist->Sort();
3806
3807   delete vect3Jet;
3808   delete templist;
3809   delete [] areaJetLarge;
3810   delete [] areaJet;
3811   delete ar;
3812   delete [] randIndex;
3813
3814 }
3815
3816 // ________________________________________________________________________________________________________________________________________________________
3817 void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJetsStat(Int_t nCases, TList* inputlist, TList* outputlist, TList* jetlist, Double_t& sumPt, Double_t &normFactor)
3818 {
3819   // List of tracks outside cone around N jet axis
3820   // All particles taken + final scaling factor
3821
3822   sumPt = 0;
3823   Float_t rc  = TMath::Abs(GetFFRadius());
3824   Float_t rcl = GetFFBckgRadius();
3825
3826   // Estimate jet and background areas
3827   Float_t* areaJet = new Float_t[nCases];
3828   memset(areaJet, 0, sizeof(Float_t) * nCases);
3829   Float_t* areaJetLarge = new Float_t[nCases];
3830   memset(areaJetLarge, 0, sizeof(Float_t) * nCases);
3831   Float_t areaFull = (fTrackEtaMax-fTrackEtaMin)*(fTrackPhiMax-fTrackPhiMin);
3832   Float_t areaOut = areaFull;
3833
3834   //estimate jets and background areas
3835   Int_t nOut = 0;
3836   Int_t ijet = 0;
3837   TClonesArray *vect3Jet = new TClonesArray("TVector3",nCases);
3838
3839   for(Int_t ij=0; ij<nCases; ++ij)
3840     {
3841       // Get jet information
3842       AliAODJet* jet = dynamic_cast<AliAODJet*>(jetlist->At(ij));
3843       if(!jet)continue;
3844       TVector3 jet3mom;
3845       jet3mom.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
3846       new((*vect3Jet)[ijet]) TVector3((TVector3)jet3mom);
3847       Float_t etaJet = (Float_t)((TVector3*) vect3Jet->At(ij))->Eta();
3848
3849       // Jet area
3850       areaJet[ij] = CalcJetArea(etaJet,rc);
3851
3852       // Area jet larger angle
3853       areaJetLarge[ij] = CalcJetArea(etaJet,rcl);
3854
3855       // Outside jets area
3856       areaOut = areaOut - areaJetLarge[ij];
3857       ijet++;
3858     }
3859
3860   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3861
3862     // embedded tracks
3863     if( fUseExtraTracksBgr != 1){
3864       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (inputlist->At(itrack))){
3865         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3866         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3867       }
3868     }
3869
3870     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3871     if(!track)continue;
3872     Double_t trackMom[3];
3873     track->PxPyPz(trackMom);
3874     TVector3 track3mom(trackMom);
3875
3876     Double_t *dR = new Double_t[nCases];
3877     for(Int_t ij=0; ij<nCases; ij++)
3878         dR[ij] = (Double_t)((TVector3*) vect3Jet->At(ij))->DeltaR(track3mom);
3879
3880     if((nCases==0) ||
3881        (nCases==1 && (dR[0]>rcl)) ||
3882        (nCases==2 && (dR[0]>rcl && dR[1]>rcl)) ||
3883        (nCases==3 && (dR[0]>rcl && dR[1]>rcl && dR[2]>rcl)))
3884       {
3885         outputlist->Add(track);
3886         sumPt += track->Pt();
3887         nOut++;
3888       }
3889     delete [] dR;
3890   }
3891
3892   if(nCases==0) areaJet[0] = TMath::Pi()*rc*rc;
3893   normFactor = (Float_t) 1./(areaJet[0] / areaOut);
3894
3895   outputlist->Sort();
3896
3897   delete vect3Jet;
3898   delete [] areaJetLarge;
3899   delete [] areaJet;
3900
3901 }
3902
3903 // ______________________________________________________________________________________________________________________________________________________
3904 Float_t AliAnalysisTaskFragmentationFunction::CalcJetArea(Float_t etaJet, Float_t rc) const
3905 {
3906   // calculate area of jet with eta etaJet and radius rc
3907
3908   Float_t detamax = etaJet + rc;
3909   Float_t detamin = etaJet - rc;
3910   Float_t accmax = 0.0; Float_t accmin = 0.0;
3911   if(detamax > fTrackEtaMax){ // sector outside etamax
3912     Float_t h = fTrackEtaMax - etaJet;
3913     accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
3914   }
3915   if(detamin < fTrackEtaMin){ // sector outside etamin
3916     Float_t h = fTrackEtaMax + etaJet;
3917     accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
3918   }
3919   Float_t areaJet = rc*rc*TMath::Pi() - accmax - accmin;
3920
3921   return areaJet;
3922
3923 }
3924
3925 // ___________________________________________________________________________________________________________________________
3926 void AliAnalysisTaskFragmentationFunction::GetClusterTracksOutOf1Jet(AliAODJet* jet, TList* outputlist, Double_t &normFactor)
3927 {
3928   // fill tracks from bckgCluster branch in list,
3929   // for all clusters outside jet cone
3930   // sum up total area of clusters
3931
3932   Double_t rc  = GetFFRadius();
3933   Double_t rcl = GetFFBckgRadius();
3934
3935   Double_t areaTotal   = 0;
3936   Double_t sumPtTotal  = 0;
3937
3938   for(Int_t ij=0; ij<fBckgJetsRec->GetEntries(); ++ij){
3939
3940     AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij)); // not 'recCuts': use all clusters in full eta range
3941
3942     Double_t dR = jet->DeltaR(bgrCluster);
3943
3944     if(dR<rcl) continue;
3945
3946     Double_t clusterPt = bgrCluster->Pt();
3947     Double_t area      = bgrCluster->EffectiveAreaCharged();
3948     areaTotal  += area;
3949     sumPtTotal += clusterPt;
3950
3951     Int_t nTracksJet = bgrCluster->GetRefTracks()->GetEntries();
3952
3953     for(Int_t it = 0; it<nTracksJet; it++){
3954
3955       // embedded tracks - note: using ref tracks here, fBranchRecBckgClusters has to be consistent
3956       if( fUseExtraTracksBgr != 1){
3957         if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (bgrCluster->GetTrack(it))){
3958           if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3959           if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
3960         }
3961       }
3962
3963       AliVParticle*   track = dynamic_cast<AliVParticle*>(bgrCluster->GetTrack(it));
3964       if(!track) continue;
3965
3966       Float_t trackPt  = track->Pt();
3967       Float_t trackEta = track->Eta();
3968       Float_t trackPhi = TVector2::Phi_0_2pi(track->Phi());
3969
3970       if(trackEta < fTrackEtaMin || trackEta > fTrackEtaMax) continue;
3971       if(trackPhi < fTrackPhiMin || trackPhi > fTrackPhiMax) continue;
3972       if(trackPt  < fTrackPtCut) continue;
3973
3974       outputlist->Add(track);
3975     }
3976   }
3977
3978   Double_t areaJet = TMath::Pi()*rc*rc;
3979   if(areaTotal) normFactor = (Float_t) 1./(areaJet / areaTotal);
3980
3981   outputlist->Sort();
3982 }
3983
3984 // _______________________________________________________________________________________________________________________
3985 void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputlist, Double_t &normFactor)
3986 {
3987   // fill tracks from bckgCluster branch,
3988   // using cluster with median density (odd number of clusters)
3989   // or picking randomly one of the two closest to median (even number)
3990
3991   normFactor = 0;
3992
3993   Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range
3994
3995   if(nBckgClusters<3) return; // need at least 3 clusters (skipping 2 highest)
3996
3997   Double_t* bgrDensity = new Double_t[nBckgClusters];
3998   Int_t*    indices    = new Int_t[nBckgClusters];
3999
4000   for(Int_t ij=0; ij<nBckgClusters; ++ij){
4001
4002     AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij));
4003     Double_t clusterPt    = bgrCluster->Pt();
4004     Double_t area         = bgrCluster->EffectiveAreaCharged();
4005
4006     Double_t density = 0;
4007     if(area>0) density = clusterPt/area;
4008
4009     bgrDensity[ij] = density;
4010     indices[ij]    = ij;
4011   }
4012
4013   TMath::Sort(nBckgClusters, bgrDensity, indices);
4014
4015   // get median cluster
4016
4017   AliAODJet* medianCluster = 0;
4018   Double_t   medianDensity = 0;
4019
4020   if(TMath::Odd(nBckgClusters)){
4021
4022     Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters-1))];
4023     medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex));
4024
4025     Double_t clusterPt = medianCluster->Pt();
4026     Double_t area      = medianCluster->EffectiveAreaCharged();
4027
4028     if(area>0) medianDensity = clusterPt/area;
4029   }
4030   else{
4031
4032     Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters-1)];
4033     Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters)];
4034
4035     AliAODJet* medianCluster1 = (AliAODJet*)(fBckgJetsRec->At(medianIndex1));
4036     AliAODJet* medianCluster2 = (AliAODJet*)(fBckgJetsRec->At(medianIndex2));
4037
4038     Double_t density1 = 0;
4039     Double_t clusterPt1 = medianCluster1->Pt();
4040     Double_t area1      = medianCluster1->EffectiveAreaCharged();
4041     if(area1>0) density1 = clusterPt1/area1;
4042
4043     Double_t density2 = 0;
4044     Double_t clusterPt2 = medianCluster2->Pt();
4045     Double_t area2      = medianCluster2->EffectiveAreaCharged();
4046     if(area2>0) density2 = clusterPt2/area2;
4047
4048     medianDensity = 0.5*(density1+density2);
4049
4050     medianCluster = ( (fRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 );  // select one randomly to avoid adding areas
4051   }
4052
4053   Int_t nTracksJet = medianCluster->GetRefTracks()->GetEntries();
4054
4055   for(Int_t it = 0; it<nTracksJet; it++){
4056
4057     // embedded tracks - note: using ref tracks here, fBranchRecBckgClusters has to be consistent
4058     if( fUseExtraTracksBgr != 1){
4059       if(AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*> (medianCluster->GetTrack(it))){
4060         if(fUseExtraTracksBgr == 0  &&  ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
4061         if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue;
4062       }
4063     }
4064
4065     AliVParticle* track = dynamic_cast<AliVParticle*>(medianCluster->GetTrack(it));
4066     if(!track) continue;
4067
4068     Float_t trackPt  = track->Pt();
4069     Float_t trackEta = track->Eta();
4070     Float_t trackPhi = TVector2::Phi_0_2pi(track->Phi());
4071
4072     if(trackEta < fTrackEtaMin || trackEta > fTrackEtaMax) continue;
4073     if(trackPhi < fTrackPhiMin || trackPhi > fTrackPhiMax) continue;
4074     if(trackPt  < fTrackPtCut) continue;
4075
4076     outputlist->Add(track);
4077   }
4078
4079   Double_t areaMedian = medianCluster->EffectiveAreaCharged();
4080   Double_t areaJet = TMath::Pi()*GetFFRadius()*GetFFRadius();
4081
4082   if(areaMedian) normFactor = (Float_t) 1./(areaJet / areaMedian);
4083
4084   outputlist->Sort();
4085
4086   delete[] bgrDensity;
4087   delete[] indices;
4088 }
4089
4090 // ______________________________________________________________________________________________________________________________________________________
4091 void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inputtracklist, TList* inputjetlist, AliAODJet* jet,
4092                                                              AliFragFuncHistos* ffbckghistocuts, AliFragFuncQATrackHistos* qabckghistocuts, TH1F* fh1Mult){
4093
4094   // List of tracks outside jets for background study
4095   TList* tracklistout2jets     = new TList();
4096   TList* tracklistout3jets     = new TList();
4097   TList* tracklistout2jetsStat = new TList();
4098   TList* tracklistout3jetsStat = new TList();
4099   Double_t sumPtOut2Jets       = 0.;
4100   Double_t sumPtOut3Jets       = 0.;
4101   Double_t sumPtOut2JetsStat   = 0.;
4102   Double_t sumPtOut3JetsStat   = 0.;
4103   Double_t normFactor2Jets     = 0.;
4104   Double_t normFactor3Jets     = 0.;
4105
4106   Int_t nRecJetsCuts = inputjetlist->GetEntries();
4107
4108   if(nRecJetsCuts>1) {
4109     GetTracksOutOfNJets(2,inputtracklist, tracklistout2jets, inputjetlist, sumPtOut2Jets);
4110     GetTracksOutOfNJetsStat(2,inputtracklist, tracklistout2jetsStat, inputjetlist,sumPtOut2JetsStat, normFactor2Jets);
4111
4112   }
4113   if(nRecJetsCuts>2) {
4114     GetTracksOutOfNJets(3,inputtracklist, tracklistout3jets, inputjetlist, sumPtOut3Jets);
4115     GetTracksOutOfNJetsStat(3,inputtracklist, tracklistout3jetsStat, inputjetlist, sumPtOut3JetsStat, normFactor3Jets);
4116   }
4117
4118   if(type==kBckgOutLJ || type==kBckgOutAJ)
4119     {
4120       TList* tracklistoutleading   = new TList();
4121       Double_t sumPtOutLeading     = 0.;
4122       GetTracksOutOfNJets(1,inputtracklist, tracklistoutleading, inputjetlist, sumPtOutLeading);
4123       if(type==kBckgOutLJ && fh1Mult) fh1Mult->Fill(tracklistoutleading->GetSize());
4124
4125       for(Int_t it=0; it<tracklistoutleading->GetSize(); ++it){
4126
4127         AliVParticle* trackVP   = (AliVParticle*)(tracklistoutleading->At(it));
4128         if(!trackVP) continue;
4129         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4130
4131         Float_t jetPt   = jet->Pt();
4132         Float_t trackPt = trackV->Pt();
4133
4134         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4135
4136         if(type==kBckgOutLJ)
4137           {
4138             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt);
4139
4140             // Fill track QA for background
4141             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4142           }
4143
4144         // All cases included
4145         if(nRecJetsCuts==1 && type==kBckgOutAJ)
4146           {
4147             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4148           }
4149         delete trackV;
4150       }
4151       // Increment jet pt with one entry in case #tracks outside jets = 0
4152       if(tracklistoutleading->GetSize()==0) {
4153          Float_t jetPt = jet->Pt();
4154          Bool_t incrementJetPt = kTRUE;
4155          if(type==kBckgOutLJ)
4156           {
4157             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4158           }
4159         // All cases included
4160         if(nRecJetsCuts==1 && type==kBckgOutAJ)
4161           {
4162             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4163           }
4164       }
4165       delete tracklistoutleading;
4166     }
4167   if(type==kBckgOutLJStat || type==kBckgOutAJStat)
4168     {
4169       TList* tracklistoutleadingStat   = new TList();
4170       Double_t sumPtOutLeadingStat = 0.;
4171       Double_t normFactorLeading   = 0.;
4172
4173       GetTracksOutOfNJetsStat(1,inputtracklist, tracklistoutleadingStat, inputjetlist, sumPtOutLeadingStat, normFactorLeading);
4174       if(type==kBckgOutLJStat && fh1Mult) fh1Mult->Fill(tracklistoutleadingStat->GetSize());
4175
4176       for(Int_t it=0; it<tracklistoutleadingStat->GetSize(); ++it){
4177
4178         AliVParticle* trackVP   = dynamic_cast<AliVParticle*>(tracklistoutleadingStat->At(it));
4179         if(!trackVP) continue;
4180         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4181
4182         Float_t jetPt   = jet->Pt();
4183         Float_t trackPt = trackV->Pt();
4184         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4185
4186         // Stat plots
4187         if(type==kBckgOutLJStat)
4188           {
4189             if(fFFMode)ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorLeading);
4190
4191             // Fill track QA for background
4192             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt); // OB added bgr QA
4193           }
4194
4195         // All cases included
4196         if(nRecJetsCuts==1 && type==kBckgOutAJStat)
4197           {
4198             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorLeading);
4199             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4200
4201           }
4202         delete trackV;
4203       }
4204       // Increment jet pt with one entry in case #tracks outside jets = 0
4205       if(tracklistoutleadingStat->GetSize()==0) {
4206          Float_t jetPt = jet->Pt();
4207          Bool_t incrementJetPt = kTRUE;
4208          if(type==kBckgOutLJStat)
4209           {
4210             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorLeading);
4211           }
4212         // All cases included
4213         if(nRecJetsCuts==1 && type==kBckgOutLJStat)
4214           {
4215             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorLeading);
4216           }
4217       }
4218
4219       delete tracklistoutleadingStat;
4220     }
4221
4222   if(type==kBckgPerp || type==kBckgPerp2 || type==kBckgPerp2Area)
4223     {
4224       Double_t sumPtPerp1 = 0.;
4225       Double_t sumPtPerp2 = 0.;
4226       TList* tracklistperp  = new TList();
4227       TList* tracklistperp1 = new TList();
4228       TList* tracklistperp2 = new TList();
4229
4230       Double_t norm = 1;
4231       if(type == kBckgPerp2)     norm = 2; // in FillFF() scaleFac = 1/norm = 0.5 - account for double area;
4232       if(type == kBckgPerp2Area) norm = 2*TMath::Pi()*TMath::Abs(GetFFRadius())*TMath::Abs(GetFFRadius()) / jet->EffectiveAreaCharged(); // in FillFF() scaleFac = 1/norm;
4233
4234       GetTracksTiltedwrpJetAxis(TMath::Pi()/2., inputtracklist,tracklistperp1,jet,TMath::Abs(GetFFRadius()),sumPtPerp1);
4235       if(type==kBckgPerp2 || type==kBckgPerp2Area) GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2., inputtracklist,tracklistperp2,jet,TMath::Abs(GetFFRadius()),sumPtPerp2);
4236
4237       tracklistperp->AddAll(tracklistperp1);
4238       tracklistperp->AddAll(tracklistperp2);
4239
4240       if(tracklistperp->GetSize() != tracklistperp1->GetSize() + tracklistperp2->GetSize()){
4241         cout<<" ERROR: tracklistperp size "<<tracklistperp->GetSize()<<" perp1 "<<tracklistperp1->GetSize()<<" perp2 "<<tracklistperp2->GetSize()<<endl;
4242         exit(0);
4243       }
4244
4245       if(fh1Mult) fh1Mult->Fill(tracklistperp->GetSize());
4246
4247       for(Int_t it=0; it<tracklistperp->GetSize(); ++it){
4248
4249         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistperp->At(it));
4250         if(!trackVP)continue;
4251         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4252
4253         Float_t jetPt   = jet->Pt();
4254         Float_t trackPt = trackV->Pt();
4255
4256         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4257
4258         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, norm );
4259
4260         // Fill track QA for background
4261         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4262
4263         delete trackV;
4264       }
4265       // Increment jet pt with one entry in case #tracks outside jets = 0
4266       if(tracklistperp->GetSize()==0) {
4267          Float_t jetPt = jet->Pt();
4268          Bool_t incrementJetPt = kTRUE;
4269          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4270       }
4271
4272
4273       if(fJSMode){
4274         // fill for tracklistperp1/2 separately, divide norm by 2
4275         if(type==kBckgPerp){
4276           FillJetShape(jet, tracklistperp, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, TMath::Pi()/2., 0., kFALSE);
4277         }
4278         if(type==kBckgPerp2){
4279           FillJetShape(jet, tracklistperp1, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0,    TMath::Pi()/2., 0., kFALSE);
4280           FillJetShape(jet, tracklistperp2, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, -1*TMath::Pi()/2., 0., kFALSE);
4281         }
4282         if(type==kBckgPerp2Area){ // divide norm by 2: listperp1/2 filled separately
4283           FillJetShape(jet, tracklistperp1, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0,    TMath::Pi()/2., 0.5*norm, kFALSE);
4284           FillJetShape(jet, tracklistperp2, fProNtracksLeadingJetBgrPerp2, fProDelRPtSumBgrPerp2, 0, -1*TMath::Pi()/2., 0.5*norm, kFALSE);
4285         }
4286       }
4287
4288       delete tracklistperp;
4289       delete tracklistperp1;
4290       delete tracklistperp2;
4291     }
4292
4293  if(type==kBckgASide)
4294     {
4295       Double_t sumPtASide = 0.;
4296       TList* tracklistaside = new TList();
4297       GetTracksTiltedwrpJetAxis(TMath::Pi(),inputtracklist,tracklistaside,jet,TMath::Abs(GetFFRadius()),sumPtASide);
4298       if(fh1Mult) fh1Mult->Fill(tracklistaside->GetSize());
4299
4300       for(Int_t it=0; it<tracklistaside->GetSize(); ++it){
4301
4302         AliVParticle*   trackVP = (AliVParticle*)(tracklistaside->At(it));
4303         if(!trackVP) continue;
4304         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4305
4306         Float_t jetPt   = jet->Pt();
4307         Float_t trackPt = trackV->Pt();
4308
4309         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4310
4311         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4312
4313         // Fill track QA for background
4314         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4315
4316         delete trackV;
4317       }
4318       if(tracklistaside->GetSize()==0) {
4319          Float_t jetPt = jet->Pt();
4320          Bool_t incrementJetPt = kTRUE;
4321          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4322       }
4323
4324       delete tracklistaside;
4325     }
4326
4327   if(type==kBckgASideWindow)
4328     {
4329       Double_t normFactorASide = 0.;
4330       Double_t sumPtASideW = 0.;
4331       TList* tracklistasidew = new TList();
4332       GetTracksTiltedwrpJetAxisWindow(TMath::Pi(),inputtracklist,tracklistasidew,jet,TMath::Abs(GetFFRadius()),sumPtASideW,normFactorASide);
4333       if(fh1Mult) fh1Mult->Fill(tracklistasidew->GetSize());
4334
4335       for(Int_t it=0; it<tracklistasidew->GetSize(); ++it){
4336
4337         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistasidew->At(it));
4338         if(!trackVP) continue;
4339         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4340
4341         Float_t jetPt   = jet->Pt();
4342         Float_t trackPt = trackV->Pt();
4343         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4344
4345         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorASide);
4346
4347         // Fill track QA for background
4348         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt, kFALSE, normFactorASide);
4349
4350         delete trackV;
4351       }
4352       if(tracklistasidew->GetSize()==0) {
4353          Float_t jetPt = jet->Pt();
4354          Bool_t incrementJetPt = kTRUE;
4355          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorASide);
4356       }
4357
4358       delete tracklistasidew;
4359     }
4360
4361   if(type==kBckgPerpWindow)
4362     {
4363       Double_t normFactorPerp = 0.;
4364       Double_t sumPtPerpW = 0.;
4365       TList* tracklistperpw = new TList();
4366       GetTracksTiltedwrpJetAxisWindow(TMath::Pi()/2.,inputtracklist,tracklistperpw,jet,TMath::Abs(GetFFRadius()),sumPtPerpW,normFactorPerp);
4367       if(fh1Mult) fh1Mult->Fill(tracklistperpw->GetSize());
4368
4369       for(Int_t it=0; it<tracklistperpw->GetSize(); ++it){
4370
4371         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistperpw->At(it));
4372         if(!trackVP) continue;
4373         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4374
4375         Float_t jetPt   = jet->Pt();
4376         Float_t trackPt = trackV->Pt();
4377         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4378
4379         if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorPerp);
4380
4381         // Fill track QA for background
4382         if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt, kFALSE, normFactorPerp);
4383
4384         delete trackV;
4385       }
4386       if(tracklistperpw->GetSize()==0) {
4387          Float_t jetPt = jet->Pt();
4388          Bool_t incrementJetPt = kTRUE;
4389          if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactorPerp);
4390       }
4391
4392       delete tracklistperpw;
4393     }
4394
4395
4396   if(type==kBckgOut2J || type==kBckgOutAJ)
4397     {
4398       if(type==kBckgOut2J && fh1Mult) fh1Mult->Fill(tracklistout2jets->GetSize());
4399       for(Int_t it=0; it<tracklistout2jets->GetSize(); ++it){
4400
4401         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout2jets->At(it));
4402         if(!trackVP) continue;
4403         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4404
4405         Float_t jetPt   = jet->Pt();
4406         Float_t trackPt = trackV->Pt();
4407
4408         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4409
4410         if(type==kBckgOut2J)
4411           {
4412             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4413             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4414           }
4415
4416         // All cases included
4417         if(nRecJetsCuts==2 && type==kBckgOutAJ)
4418           {
4419             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4420
4421           }
4422         delete trackV;
4423       }
4424       // Increment jet pt with one entry in case #tracks outside jets = 0
4425       if(tracklistout2jets->GetSize()==0) {
4426          Float_t jetPt = jet->Pt();
4427          Bool_t incrementJetPt = kTRUE;
4428          if(type==kBckgOut2J)
4429           {
4430             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4431           }
4432         // All cases included
4433         if(nRecJetsCuts==2 && type==kBckgOutAJ)
4434           {
4435             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4436           }
4437       }
4438     }
4439
4440   if(type==kBckgOut2JStat || type==kBckgOutAJStat)
4441     {
4442       for(Int_t it=0; it<tracklistout2jetsStat->GetSize(); ++it){
4443
4444         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout2jetsStat->At(it));
4445         if(!trackVP) continue;
4446         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4447
4448         Float_t jetPt   = jet->Pt();
4449         Float_t trackPt = trackV->Pt();
4450         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4451
4452         if(type==kBckgOut2JStat)
4453           {
4454             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor2Jets);
4455
4456             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4457           }
4458
4459         // All cases included
4460         if(nRecJetsCuts==2 && type==kBckgOutAJStat)
4461           {
4462             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor2Jets);
4463
4464             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt ); // OB added bgr QA
4465           }
4466         delete trackV;
4467       }
4468       // Increment jet pt with one entry in case #tracks outside jets = 0
4469       if(tracklistout2jetsStat->GetSize()==0) {
4470          Float_t jetPt = jet->Pt();
4471          Bool_t incrementJetPt = kTRUE;
4472          if(type==kBckgOut2JStat)
4473            {
4474             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor2Jets);
4475           }
4476         // All cases included
4477         if(nRecJetsCuts==2 && type==kBckgOutAJStat)
4478           {
4479             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor2Jets);
4480           }
4481       }
4482
4483     }
4484
4485   if(type==kBckgOut3J || type==kBckgOutAJ)
4486     {
4487       if(type==kBckgOut3J && fh1Mult) fh1Mult->Fill(tracklistout3jets->GetSize());
4488
4489       for(Int_t it=0; it<tracklistout3jets->GetSize(); ++it){
4490
4491         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout3jets->At(it));
4492         if(!trackVP) continue;
4493         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4494
4495         Float_t jetPt   = jet->Pt();
4496         Float_t trackPt = trackV->Pt();
4497
4498         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4499
4500         if(type==kBckgOut3J)
4501           {
4502             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4503
4504             qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt);
4505           }
4506
4507         // All cases included
4508         if(nRecJetsCuts==3 && type==kBckgOutAJ)
4509           {
4510             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt );
4511
4512           }
4513         delete trackV;
4514       }
4515       // Increment jet pt with one entry in case #tracks outside jets = 0
4516       if(tracklistout3jets->GetSize()==0) {
4517          Float_t jetPt = jet->Pt();
4518          Bool_t incrementJetPt = kTRUE;
4519          if(type==kBckgOut3J)
4520           {
4521             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4522           }
4523         // All cases included
4524         if(nRecJetsCuts==3 && type==kBckgOutAJ)
4525           {
4526             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt );
4527           }
4528       }
4529     }
4530
4531   if(type==kBckgOut3JStat || type==kBckgOutAJStat)
4532     {
4533       for(Int_t it=0; it<tracklistout3jetsStat->GetSize(); ++it){
4534
4535         AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistout3jetsStat->At(it));
4536         if(!trackVP) continue;
4537         TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4538
4539         Float_t jetPt   = jet->Pt();
4540         Float_t trackPt = trackV->Pt();
4541         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4542
4543         if(type==kBckgOut3JStat)
4544           {
4545             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor3Jets);
4546
4547             //if(fQAMode&1)     qabckghistocuts->FillTrackQA( trackEta, TVector2::Phi_0_2pi(trackPhi), trackPt);
4548           }
4549
4550         // All cases included
4551         if(nRecJetsCuts==3 && type==kBckgOutAJStat)
4552           {
4553             if(fFFMode) ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactor3Jets );
4554
4555             if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4556
4557           }
4558         delete trackV;
4559       }
4560       // Increment jet pt with one entry in case #tracks outside jets = 0
4561       if(tracklistout3jetsStat->GetSize()==0) {
4562          Float_t jetPt = jet->Pt();
4563          Bool_t incrementJetPt = kTRUE;
4564          if(type==kBckgOut3JStat)
4565           {
4566             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor3Jets);
4567           }
4568         // All cases included
4569         if(nRecJetsCuts==3 && type==kBckgOutAJStat)
4570           {
4571             if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor3Jets);
4572           }
4573       }
4574
4575     }
4576
4577   if(type==kBckgClustersOutLeading){ // clusters bgr: all tracks in clusters out of leading jet
4578
4579     TList* tracklistClustersOutLeading = new TList();
4580     Double_t normFactorClusters = 0;
4581     Float_t jetPt   = jet->Pt();
4582
4583     GetClusterTracksOutOf1Jet(jet, tracklistClustersOutLeading, normFactorClusters);
4584     if(fh1Mult) fh1Mult->Fill(tracklistClustersOutLeading->GetSize());
4585
4586     for(Int_t it=0; it<tracklistClustersOutLeading->GetSize(); ++it){
4587
4588       AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistClustersOutLeading->At(it));
4589       if(!trackVP) continue;
4590       TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4591
4592       Float_t trackPt = trackVP->Pt();
4593
4594       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4595
4596       if(fFFMode)   ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorClusters );
4597       if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4598
4599       delete trackV;
4600     }
4601
4602     delete tracklistClustersOutLeading;
4603
4604   }
4605
4606   if(type == kBckgClusters){ // clusters bgr: all tracks in 'median cluster'
4607
4608     TList* tracklistClustersMedian = new TList();
4609     Double_t normFactorClusters = 0;
4610     Float_t jetPt = jet->Pt();
4611
4612     GetClusterTracksMedian(tracklistClustersMedian, normFactorClusters);
4613     if(fh1Mult) fh1Mult->Fill(tracklistClustersMedian->GetSize());
4614
4615     for(Int_t it=0; it<tracklistClustersMedian->GetSize(); ++it){
4616
4617       AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(tracklistClustersMedian->At(it));
4618       if(!trackVP) continue;
4619       TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
4620
4621       Float_t trackPt = trackVP->Pt();
4622
4623       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4624
4625       if(fFFMode)   ffbckghistocuts->FillFF( trackPt, jetPt, incrementJetPt, normFactorClusters );
4626       if(fQAMode&1) qabckghistocuts->FillTrackQA( trackV->Eta(), TVector2::Phi_0_2pi(trackV->Phi()), trackPt );
4627
4628       delete trackV;
4629     }
4630
4631     delete tracklistClustersMedian;
4632   }
4633
4634   delete tracklistout2jets;
4635   delete tracklistout3jets;
4636   delete tracklistout2jetsStat;
4637   delete tracklistout3jetsStat;
4638 }
4639
4640 //_____________________________________________________________________________________
4641 Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactor(Double_t pt)
4642 {
4643   // factor strangeness data/MC as function of pt from UE analysis (Sara Vallero)
4644
4645   Double_t alpha = 1;
4646
4647   if(0.150<pt && pt<0.200) alpha = 3.639;
4648   if(0.200<pt && pt<0.250) alpha = 2.097;
4649   if(0.250<pt && pt<0.300) alpha = 1.930;
4650   if(0.300<pt && pt<0.350) alpha = 1.932;
4651   if(0.350<pt && pt<0.400) alpha = 1.943;
4652   if(0.400<pt && pt<0.450) alpha = 1.993;
4653   if(0.450<pt && pt<0.500) alpha = 1.989;
4654   if(0.500<pt && pt<0.600) alpha = 1.963;
4655   if(0.600<pt && pt<0.700) alpha = 1.917;
4656   if(0.700<pt && pt<0.800) alpha = 1.861;
4657   if(0.800<pt && pt<0.900) alpha = 1.820;
4658   if(0.900<pt && pt<1.000) alpha = 1.741;
4659   if(1.000<pt && pt<1.500) alpha = 0.878;
4660
4661   return alpha;
4662 }
4663
4664 //__________________________________________________________________________________________________
4665 Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactorCMS(AliAODMCParticle* daughter)
4666 {
4667   // strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
4668
4669   TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
4670   if(!tca) return 1;
4671
4672   AliAODMCParticle* currentMother   = daughter;
4673   AliAODMCParticle* currentDaughter = daughter;
4674
4675
4676   // find first primary mother K0s, Lambda or Xi
4677   while(1){
4678
4679     Int_t daughterPDG   = currentDaughter->GetPdgCode();
4680
4681     Int_t motherLabel   = currentDaughter->GetMother();
4682     if(motherLabel >= tca->GetEntriesFast()){ // protection
4683       currentMother = currentDaughter;
4684       break;
4685     }
4686
4687     currentMother     = (AliAODMCParticle*) tca->At(motherLabel);
4688
4689     if(!currentMother){
4690       currentMother = currentDaughter;
4691       break;
4692     }
4693
4694     Int_t motherPDG   = currentMother->GetPdgCode();
4695
4696     // phys. primary found ?
4697     if(currentMother->IsPhysicalPrimary()) break;
4698
4699     if(TMath::Abs(daughterPDG) == 321){ // K+/K- e.g. from phi (ref data not feeddown corrected)
4700       currentMother = currentDaughter; break;
4701     }
4702     if(TMath::Abs(motherPDG) == 310 ){ // K0s e.g. from phi (ref data not feeddown corrected)
4703       break;
4704     }
4705     if(TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122){ // mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
4706       currentMother = currentDaughter; break;
4707     }
4708
4709     currentDaughter = currentMother;
4710   }
4711
4712
4713   Int_t motherPDG   = currentMother->GetPdgCode();
4714   Double_t motherPt = currentMother->Pt();
4715
4716   Double_t fac = 1;
4717
4718   if(TMath::Abs(motherPDG) == 310 || TMath::Abs(motherPDG)==321){ // K0s / K+ / K-
4719
4720     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.768049;
4721     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.732933;
4722     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.650298;
4723     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.571332;
4724     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.518734;
4725     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.492543;
4726     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.482704;
4727     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.488056;
4728     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.488861;
4729     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.492862;
4730     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.504332;
4731     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.501858;
4732     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.512970;
4733     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.524131;
4734     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.539130;
4735     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.554101;
4736     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.560348;
4737     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.568869;
4738     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.583310;
4739     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.604818;
4740     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.632630;
4741     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.710070;
4742     else if(6.00 <= motherPt && motherPt < 8.00) fac = 0.736365;
4743     else if(8.00 <= motherPt && motherPt < 10.00) fac = 0.835865;
4744   }
4745
4746   if(TMath::Abs(motherPDG) == 3122){ // Lambda
4747
4748     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.645162;
4749     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.627431;
4750     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.457136;
4751     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.384369;
4752     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.330597;
4753     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.309571;
4754     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.293620;
4755     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.283709;
4756     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.282047;
4757     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.277261;
4758     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.275772;
4759     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.280726;
4760     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.288540;
4761     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.288315;
4762     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.296619;
4763     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.302993;
4764     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.338121;
4765     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.349800;
4766     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.356802;
4767     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.391202;
4768     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.422573;
4769     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.573815;
4770     else if(6.00 <= motherPt && motherPt < 8.00) fac = 0.786984;
4771     else if(8.00 <= motherPt && motherPt < 10.00) fac = 1.020021;
4772   }
4773
4774   if(TMath::Abs(motherPDG) == 3312 || TMath::Abs(motherPDG) == 3322){ // xi
4775
4776     if(0.00 <= motherPt && motherPt < 0.20) fac = 0.666620;
4777     else if(0.20 <= motherPt && motherPt < 0.40) fac = 0.575908;
4778     else if(0.40 <= motherPt && motherPt < 0.60) fac = 0.433198;
4779     else if(0.60 <= motherPt && motherPt < 0.80) fac = 0.340901;
4780     else if(0.80 <= motherPt && motherPt < 1.00) fac = 0.290896;
4781     else if(1.00 <= motherPt && motherPt < 1.20) fac = 0.236074;
4782     else if(1.20 <= motherPt && motherPt < 1.40) fac = 0.218681;
4783     else if(1.40 <= motherPt && motherPt < 1.60) fac = 0.207763;
4784     else if(1.60 <= motherPt && motherPt < 1.80) fac = 0.222848;
4785     else if(1.80 <= motherPt && motherPt < 2.00) fac = 0.208806;
4786     else if(2.00 <= motherPt && motherPt < 2.20) fac = 0.197275;
4787     else if(2.20 <= motherPt && motherPt < 2.40) fac = 0.183645;
4788     else if(2.40 <= motherPt && motherPt < 2.60) fac = 0.188788;
4789     else if(2.60 <= motherPt && motherPt < 2.80) fac = 0.188282;
4790     else if(2.80 <= motherPt && motherPt < 3.00) fac = 0.207442;
4791     else if(3.00 <= motherPt && motherPt < 3.20) fac = 0.240388;
4792     else if(3.20 <= motherPt && motherPt < 3.40) fac = 0.241916;
4793     else if(3.40 <= motherPt && motherPt < 3.60) fac = 0.208276;
4794     else if(3.60 <= motherPt && motherPt < 3.80) fac = 0.234550;
4795     else if(3.80 <= motherPt && motherPt < 4.00) fac = 0.251689;
4796     else if(4.00 <= motherPt && motherPt < 5.00) fac = 0.310204;
4797     else if(5.00 <= motherPt && motherPt < 6.00) fac = 0.343492;
4798   }
4799
4800   Double_t weight = 1;
4801   if(fac > 0) weight = 1/fac;
4802
4803   return weight;
4804 }
4805
4806 // _________________________________________________________________________________
4807 void  AliAnalysisTaskFragmentationFunction::FillJetShape(AliAODJet* jet, TList* list,
4808                                                          TProfile* hProNtracksLeadingJet, TProfile** hProDelRPtSum, TProfile* hProDelR80pcPt,
4809                                                          Double_t dPhiUE, Double_t normUE, Bool_t scaleStrangeness){
4810
4811   const Int_t   kNbinsR    = 50;
4812   const Float_t kBinWidthR = 0.02;
4813
4814   Int_t nJetTracks = list->GetEntries();
4815
4816   Float_t PtSumA[kNbinsR]     = {0.0};
4817
4818   Float_t *delRA     = new Float_t[nJetTracks];
4819   Float_t *trackPtA  = new Float_t[nJetTracks];
4820   Int_t   *index     = new Int_t[nJetTracks];
4821
4822   for(Int_t i=0; i<nJetTracks; i++){
4823     delRA[i]    = 0;
4824     trackPtA[i] = 0;
4825     index[i]    = 0;
4826   }
4827
4828   Double_t jetMom[3];
4829   jet->PxPyPz(jetMom);
4830   TVector3 jet3mom(jetMom);
4831
4832   if(TMath::Abs(dPhiUE)>0){
4833     Double_t phiTilted = jet3mom.Phi();
4834     phiTilted += dPhiUE;
4835     phiTilted = TVector2::Phi_0_2pi(phiTilted);
4836     jet3mom.SetPhi(phiTilted);
4837   }
4838
4839   Double_t jetPt = jet->Pt();
4840   Double_t sumWeights = 0;
4841
4842   for (Int_t j =0; j<nJetTracks; j++){
4843
4844     AliVParticle* track = dynamic_cast<AliVParticle*>(list->At(j));
4845     if(!track)continue;
4846
4847     Double_t trackMom[3];
4848     track->PxPyPz(trackMom);
4849     TVector3 track3mom(trackMom);
4850
4851     Double_t dR = jet3mom.DeltaR(track3mom);
4852
4853     delRA[j]    = dR;
4854     trackPtA[j] = track->Pt();
4855
4856     Double_t weight = GetMCStrangenessFactor(track->Pt()); // more correctly should be gen pt
4857     sumWeights += weight;
4858
4859     for(Int_t ibin=1; ibin<=kNbinsR; ibin++){
4860       Float_t xlow = kBinWidthR*(ibin-1);
4861       Float_t xup  = kBinWidthR*ibin;
4862       if(xlow <= dR && dR < xup){
4863
4864         if(scaleStrangeness) PtSumA[ibin-1]     += track->Pt()*weight;
4865         else                 PtSumA[ibin-1]     += track->Pt();
4866       }
4867     }
4868   } // track loop
4869
4870   Float_t jetPtMin=0;
4871   Float_t jetPtMax=0;
4872
4873   for(Int_t ibin=0; ibin<kNbinsR; ibin++){
4874     Float_t fR =  kBinWidthR*(ibin+0.5);
4875
4876     for(Int_t k=0; k<5; k++){
4877       if(k==0){jetPtMin=20.0;jetPtMax=30.0;}
4878       if(k==1){jetPtMin=30.0;jetPtMax=40.0;}
4879       if(k==2){jetPtMin=40.0;jetPtMax=60.0;}
4880       if(k==3){jetPtMin=60.0;jetPtMax=80.0;}
4881       if(k==4){jetPtMin=80.0;jetPtMax=100.0;}
4882       if(jetPt>jetPtMin && jetPt<jetPtMax){
4883
4884         hProDelRPtSum[k]->Fill(fR,PtSumA[ibin]);
4885
4886       }
4887     }
4888   }
4889
4890   if(scaleStrangeness) hProNtracksLeadingJet->Fill(jetPt,sumWeights);
4891   else                 hProNtracksLeadingJet->Fill(jetPt,nJetTracks);
4892
4893   if(normUE)           hProNtracksLeadingJet->Fill(jetPt,nJetTracks/normUE);
4894
4895   if(hProDelR80pcPt){
4896
4897     Float_t PtSum = 0;
4898     Float_t delRPtSum80pc = 0;
4899
4900     TMath::Sort(nJetTracks,delRA,index,0);
4901
4902     for(Int_t ii=0; ii<nJetTracks; ii++){
4903
4904       if(scaleStrangeness){
4905         Double_t weight = GetMCStrangenessFactor(trackPtA[index[ii]]); // more correctly should be gen pt
4906         PtSum += weight*trackPtA[index[ii]];
4907       }
4908       else PtSum += trackPtA[index[ii]];
4909
4910
4911       if(PtSum/jetPt >= 0.8000){
4912         delRPtSum80pc = delRA[index[ii]];
4913         break;
4914       }
4915     }
4916     hProDelR80pcPt->Fill(jetPt,delRPtSum80pc);
4917   }
4918
4919   delete[] delRA;
4920   delete[] trackPtA;
4921   delete[] index;
4922 }