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