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