PWGJE/AliAnalysisTaskFragmentationFunction.cxx

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