PWGJE/AliAnalysisTaskFragmentationFunction.cxx

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