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