PWGJE/AliAnalysisTaskFragmentationFunction.cxx

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