Fixes for not filled histograms and calculation of Dijet binning
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskFragmentationFunction.cxx
index 3aa6101398cfd82ccf829150befbbcafba0a746f..ce780956c8062941eeeeada39b97bba5d0d6b280 100644 (file)
@@ -1,3 +1,10 @@
+/*************************************************************************
+ *                                                                       *
+ * Task for Fragmentation Function Analysis in PWG4 Jet Task Force Train *
+ *                                                                       *
+ *************************************************************************/
+
+
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-#include <Riostream.h>
+/* $Id: */
 
-#include <TFile.h>
-#include <TClonesArray.h>
-#include <TLorentzVector.h>
-#include <TProfile.h>
-#include <TList.h>
-#include <TH1.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TH3F.h>
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TString.h"
+#include "THnSparse.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TKey.h"
 
-#include "AliAnalysisTaskFragmentationFunction.h"
-#include "AliAnalysisManager.h"
-#include "AliInputEventHandler.h"
-#include "AliAODHandler.h"
-#include "AliAODTrack.h"
-#include "AliJetHeader.h"
-#include "AliAODEvent.h"
+#include "AliAODInputHandler.h" 
+#include "AliAODHandler.h" 
+#include "AliESDEvent.h"
+#include "AliAODMCParticle.h"
 #include "AliAODJet.h"
-#include "AliAODDiJet.h"
 #include "AliGenPythiaEventHeader.h"
+#include "AliInputEventHandler.h"
+
 #include "AliAnalysisHelperJetTasks.h"
-#include "AliMCEvent.h"
-#include "AliMCParticle.h"
-#include "AliAODMCParticle.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisTaskSE.h"
 
-ClassImp(AliAnalysisTaskFragmentationFunction)
+#include "AliAnalysisTaskFragmentationFunction.h"
 
-//#######################################################################
-AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction():
-  AliAnalysisTaskSE(),
-  fJetHeaderRec(0x0),
-  fJetHeaderGen(0x0),
-  fAOD(0x0),
-  fBranchRec(""),
-  fBranchGen(""),
-  fUseAODInput(kFALSE),
-  fUseAODJetInput(kFALSE),
-  fUseAODTrackInput(kFALSE),
-  fUseAODMCInput(kFALSE),
-  fUseGlobalSelection(kFALSE),
-  fUseExternalWeightOnly(kFALSE),
-  fLimitGenJetEta(kFALSE),
-  fFilterMask(0),
-  fAnalysisType(0),
-  fTrackTypeRec(kTrackUndef),  
-  fTrackTypeGen(kTrackUndef),
-  fAvgTrials(1.),
-  fExternalWeight(1.),    
-  fRecEtaWindow(0.5),
-  fR(0.),
-  fdRdNdxi(0.7),
-  fPartPtCut(0.),
-  fEfactor(1.),
-  fNff(5),
-  fNim(5),
-  fList(0x0),
-  fGlobVar(1),
-  // Number of energy bins
-  fnEBin(6),         
-  fEmin(10.),
-  fEmax(70.),
-  fnEInterval(6),
-  // Number of radius bins
-  fnRBin(10),        
-  fRmin(0.1),
-  fRmax(1.),   
-  fnRInterval(9),
-  // eta histograms
-  fnEtaHBin(50),
-  fEtaHBinMin(-0.9),
-  fEtaHBinMax(0.9),
-  // phi histograms
-  fnPhiHBin(60),
-  fPhiHBinMin(0.),
-  fPhiHBinMax(2*TMath::Pi()),
-  // pt histograms
-  fnPtHBin(300),
-  fPtHBinMin(0.),
-  fPtHBinMax(300.),
-  // E histograms
-  fnEHBin(300),
-  fEHBinMin(0.),
-  fEHBinMax(300.),
-  // Xi histograms
-  fnXiHBin(27),
-  fXiHBinMin(0.),
-  fXiHBinMax(9.),
-  // Pthad histograms
-  fnPthadHBin(60),
-  fPthadHBinMin(0.),
-  fPthadHBinMax(30.),
-  // z histograms
-  fnZHBin(30), 
-  fZHBinMin(0.),
-  fZHBinMax(1.),
-  // theta histograms
-  fnThetaHBin(200),
-  fThetaHBinMin(-0.5),
-  fThetaHBinMax(1.5),
-  fnCosThetaHBin(100),
-  fcosThetaHBinMin(0.),
-  fcosThetaHBinMax(1.),
-  // kT histograms
-  fnkTHBin(25), 
-  fkTHBinMin(0.),
-  fkTHBinMax(5.),
-  // kT histograms
-  fnRHBin(10),
-  fRHBinMin(0),
-  fRHBinMax(1),
-  // pt trig
-  fnPtTrigBin(10),
-  //Histograms
-  fEtaMonoJet1H(0x0),
-  fPhiMonoJet1H(0x0),
-  fPtMonoJet1H(0x0),
-  fEMonoJet1H(0x0),
-  fdNdXiMonoJet1H(0x0),
-  fdNdPtMonoJet1H(0x0),
-  fdNdZMonoJet1H(0x0),
-  fdNdThetaMonoJet1H(0x0),
-  fdNdcosThetaMonoJet1H(0x0),
-  fdNdkTMonoJet1H(0x0),
-  fdNdpTvsZMonoJet1H(0x0),
-  fShapeMonoJet1H(0x0),
-  fNMonoJet1sH(0x0),
-  fThetaPtPartMonoJet1H(0x0),
-  fcosThetaPtPartMonoJet1H(0x0),
-  fkTPtPartMonoJet1H(0x0),
-  fThetaPtJetMonoJet1H(0x0),
-  fcosThetaPtJetMonoJet1H(0x0),
-  fkTPtJetMonoJet1H(0x0),
-  fpTPtJetMonoJet1H(0x0),
-  farrayEmin(0x0),
-  farrayEmax(0x0),
-  farrayRadii(0x0),
-  farrayPtTrigmin(0x0),
-  farrayPtTrigmax(0x0),
-  // Track control plots
-  fptAllTracks(0x0),
-  fetaAllTracks(0x0),
-  fphiAllTracks(0x0),
-  fetaphiptAllTracks(0x0),
-  fetaphiAllTracks(0x0),
-  fptAllTracksCut(0x0),
-  fetaAllTracksCut(0x0),
-  fphiAllTracksCut(0x0),
-  fetaphiptAllTracksCut(0x0),
-  fetaphiAllTracksCut(0x0),
-  fptTracks(0x0),
-  fetaTracks(0x0),
-  fphiTracks(0x0),
-  fdetaTracks(0x0),
-  fdphiTracks(0x0),
-  fetaphiptTracks(0x0),
-  fetaphiTracks(0x0),
-  fdetadphiTracks(0x0),
-  fptTracksCut(0x0),
-  fetaTracksCut(0x0),
-  fphiTracksCut(0x0),
-  fdetaTracksCut(0x0),
-  fdphiTracksCut(0x0),
-  fetaphiptTracksCut(0x0),
-  fetaphiTracksCut(0x0),
-  fdetadphiTracksCut(0x0),
-  fNPtTrig(0x0),
-  fNPtTrigCut(0x0),
-  fvertexXY(0x0),
-  fvertexZ(0x0),
-  fEvtMult(0x0),
-  fEvtMultvsJetPt(0x0),
-  fPtvsEtaJet(0x0),
-  fNpvsEtaJet(0x0),
-  fNpevtvsEtaJet(0x0),
-  fPtvsPtJet(0x0),
-  fNpvsPtJet(0x0),
-  fNpevtvsPtJet(0x0),
-  fPtvsPtJet1D(0x0),
-  fNpvsPtJet1D(0x0),
-  fNpevtvsPtJet1D(0x0),
-  fptLeadingJet(0x0),
-  fetaLeadingJet(0x0),
-  fphiLeadingJet(0x0),
-  fptJet(0x0),
-  fetaJet(0x0),
-  fphiJet(0x0),
-  fHistList(0x0),
-  fNBadRuns(0),
-  fNBadRunsH(0x0)
- {
-  //
-  // Default constructor
-  //
-  /*
-  for(int i = 0;i < kMaxStep*2;++i){
-    fhnJetContainer[i] = 0;
-  }
-  */
-//   for(int ij  = 0;ij<kMaxJets;++ij){
-//     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
-//     fh1Eta[ij] = fh1Phi[ij] = 0;
-//   }
-}
 
+ClassImp(AliAnalysisTaskFragmentationFunction)
 
-//#######################################################################
-AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const char* name):
-  AliAnalysisTaskSE(name),
-  fJetHeaderRec(0x0),
-  fJetHeaderGen(0x0),
-  fAOD(0x0),
-  fBranchRec(""),
-  fBranchGen(""),
-  fUseAODInput(kFALSE),
-  fUseAODJetInput(kFALSE),
-  fUseAODTrackInput(kFALSE),
-  fUseAODMCInput(kFALSE),
-  fUseGlobalSelection(kFALSE),
-  fUseExternalWeightOnly(kFALSE),
-  fLimitGenJetEta(kFALSE),
-  fFilterMask(0),
-  fAnalysisType(0),
-  fTrackTypeRec(kTrackUndef), 
-  fTrackTypeGen(kTrackUndef),
-  fAvgTrials(1.),
-  fExternalWeight(1.),    
-  fRecEtaWindow(0.5),
-  fR(0.),
-  fdRdNdxi(0.7),
-  fPartPtCut(0.),
-  fEfactor(1.),
-  fNff(5),
-  fNim(5),
-  fList(0x0),
-  fGlobVar(1),
-  fCDFCut(1),
-  // Number of energy bins
-  fnEBin(6),         
-  fEmin(10.),
-  fEmax(70.),
-  fnEInterval(6),
-  // Number of radius bins
-  fnRBin(10),        
-  fRmin(0.1),
-  fRmax(1.),   
-  fnRInterval(9),
-  // eta histograms
-  fnEtaHBin(50),
-  fEtaHBinMin(-0.9),
-  fEtaHBinMax(0.9),
-  // phi histograms
-  fnPhiHBin(60),
-  fPhiHBinMin(0.),
-  fPhiHBinMax(2*TMath::Pi()),
-  // pt histograms
-  fnPtHBin(300),
-  fPtHBinMin(0.),
-  fPtHBinMax(300.),
-  // E histograms
-  fnEHBin(300),
-  fEHBinMin(0.),
-  fEHBinMax(300.),
-  // Xi histograms
-  fnXiHBin(27),
-  fXiHBinMin(0.),
-  fXiHBinMax(9.),
-  // Pthad histograms
-  fnPthadHBin(60),
-  fPthadHBinMin(0.),
-  fPthadHBinMax(30.),
-  // z histograms
-  fnZHBin(30),
-  fZHBinMin(0.),
-  fZHBinMax(1.),
-  fnPtTrigBin(10),
-  // theta histograms
-  fnThetaHBin(200),
-  fThetaHBinMin(-0.5),
-  fThetaHBinMax(1.5),
-  fnCosThetaHBin(100),
-  fcosThetaHBinMin(0.),
-  fcosThetaHBinMax(1.),
-  // kT histograms
-  fnkTHBin(25),
-  fkTHBinMin(0.),
-  fkTHBinMax(5.),
-  // kT histograms
-  fnRHBin(10),
-  fRHBinMin(0.),
-  fRHBinMax(1.),
-  //Histograms
-  fEtaMonoJet1H(0x0),
-  fPhiMonoJet1H(0x0),
-  fPtMonoJet1H(0x0),
-  fEMonoJet1H(0x0),
-  fdNdXiMonoJet1H(0x0),
-  fdNdPtMonoJet1H(0x0),
-  fdNdZMonoJet1H(0x0),
-  fdNdThetaMonoJet1H(0x0),
-  fdNdcosThetaMonoJet1H(0x0),
-  fdNdkTMonoJet1H(0x0),
-  fdNdpTvsZMonoJet1H(0x0),
-  fShapeMonoJet1H(0x0),
-  fNMonoJet1sH(0x0),
-  fThetaPtPartMonoJet1H(0x0),
-  fcosThetaPtPartMonoJet1H(0x0),
-  fkTPtPartMonoJet1H(0x0),
-  fThetaPtJetMonoJet1H(0x0),
-  fcosThetaPtJetMonoJet1H(0x0),
-  fkTPtJetMonoJet1H(0x0),
-  fpTPtJetMonoJet1H(0x0),
-  farrayEmin(0x0),
-  farrayEmax(0x0),
-  farrayRadii(0x0),
-  farrayPtTrigmin(0x0),
-  farrayPtTrigmax(0x0),
-  // Track control plots
-  fptAllTracks(0x0),
-  fetaAllTracks(0x0),
-  fphiAllTracks(0x0),
-  fetaphiptAllTracks(0x0),
-  fetaphiAllTracks(0x0),
-  fptAllTracksCut(0x0),
-  fetaAllTracksCut(0x0),
-  fphiAllTracksCut(0x0),
-  fetaphiptAllTracksCut(0x0),
-  fetaphiAllTracksCut(0x0),
-  fptTracks(0x0),
-  fetaTracks(0x0),
-  fphiTracks(0x0),
-  fdetaTracks(0x0),
-  fdphiTracks(0x0),
-  fetaphiptTracks(0x0),
-  fetaphiTracks(0x0),
-  fdetadphiTracks(0x0),
-  fptTracksCut(0x0),
-  fetaTracksCut(0x0),
-  fphiTracksCut(0x0),
-  fdetaTracksCut(0x0),
-  fdphiTracksCut(0x0),
-  fetaphiptTracksCut(0x0),
-  fetaphiTracksCut(0x0),
-  fdetadphiTracksCut(0x0),
-  fNPtTrig(0x0),
-  fNPtTrigCut(0x0),
-  fvertexXY(0x0),
-  fvertexZ(0x0),
-  fEvtMult(0x0),
-  fEvtMultvsJetPt(0x0),
-  fPtvsEtaJet(0x0),
-  fNpvsEtaJet(0x0),
-  fNpevtvsEtaJet(0x0),
-  fPtvsPtJet(0x0),
-  fNpvsPtJet(0x0),
-  fNpevtvsPtJet(0x0),
-  fPtvsPtJet1D(0x0),
-  fNpvsPtJet1D(0x0),
-  fNpevtvsPtJet1D(0x0),
-  fptLeadingJet(0x0),
-  fetaLeadingJet(0x0),
-  fphiLeadingJet(0x0),
-  fptJet(0x0),
-  fetaJet(0x0),
-  fphiJet(0x0),
-  fHistList(0x0),
-  fNBadRuns(0),
-  fNBadRunsH(0x0)
+//____________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
+   : AliAnalysisTaskSE()
+   ,fESD(0)
+   ,fAOD(0)
+   ,fMCEvent(0)
+   ,fBranchRecJets("jets")
+   ,fBranchGenJets("")
+   ,fTrackTypeGen(0)
+   ,fJetTypeGen(0)
+   ,fJetTypeRecEff(0)
+   ,fFilterMask(0)
+   ,fUsePhysicsSelection(kTRUE)
+   ,fTrackPtCut(0)
+   ,fTrackEtaMin(0)
+   ,fTrackEtaMax(0)
+   ,fTrackPhiMin(0)
+   ,fTrackPhiMax(0)
+   ,fJetPtCut(0)
+   ,fJetEtaMin(0)
+   ,fJetEtaMax(0)
+   ,fJetPhiMin(0)
+   ,fJetPhiMax(0)
+   ,fDiJetCut(0)
+   ,fDiJetDeltaPhiCut(0)
+   ,fDiJetPtFractionCut(0)
+   ,fDiJetCDFCut(0)
+   ,fDiJetKindBins(0)
+   ,fFFRadius(0)
+   ,fTracksRec(0)
+   ,fTracksRecCuts(0)
+   ,fTracksGen(0)
+   ,fTracksAODMCCharged(0)
+   ,fTracksRecQualityCuts(0)
+   ,fJetsRec(0)
+   ,fJetsRecCuts(0)
+   ,fJetsGen(0)
+   ,fJetsRecEff(0)
+   ,fQATrackHistosRec(0)
+   ,fQATrackHistosRecCuts(0)
+   ,fQATrackHistosGen(0)
+   ,fQAJetHistosRec(0)
+   ,fQAJetHistosRecCuts(0)
+   ,fQAJetHistosRecCutsLeading(0)
+   ,fQAJetHistosGen(0)
+   ,fQAJetHistosGenLeading(0)
+   ,fQAJetHistosRecEffLeading(0)
+   ,fFFHistosRecCuts(0)
+   ,fFFHistosRecLeading(0)
+   ,fFFHistosRecLeadingTrack(0)
+   ,fFFHistosGen(0)
+   ,fFFHistosGenLeading(0)
+   ,fFFHistosGenLeadingTrack(0)
+   ,fIJHistosRecCuts(0)
+   ,fIJHistosRecLeading(0)
+   ,fIJHistosRecLeadingTrack(0)
+   ,fIJHistosGen(0)
+   ,fIJHistosGenLeading(0)
+   ,fIJHistosGenLeadingTrack(0)
+   ,fFFDiJetHistosRecCuts(0)
+   ,fFFDiJetHistosRecLeading(0)
+   ,fFFDiJetHistosRecLeadingTrack(0)
+   ,fFFDiJetHistosGen(0)
+   ,fFFDiJetHistosGenLeading(0)
+   ,fFFDiJetHistosGenLeadingTrack(0)
+   ,fQADiJetHistosRecCuts(0)
+   ,fQADiJetHistosGen(0)
+   ,fQATrackHighPtThreshold(0)
+   ,fFFNBinsJetPt(0)    
+   ,fFFJetPtMin(0) 
+   ,fFFJetPtMax(0)
+   ,fFFNBinsPt(0)      
+   ,fFFPtMin(0)        
+   ,fFFPtMax(0)        
+   ,fFFNBinsXi(0)      
+   ,fFFXiMin(0)        
+   ,fFFXiMax(0)        
+   ,fFFNBinsZ(0)       
+   ,fFFZMin(0)         
+   ,fFFZMax(0)         
+   ,fQAJetNBinsPt(0)   
+   ,fQAJetPtMin(0)     
+   ,fQAJetPtMax(0)     
+   ,fQAJetNBinsEta(0)  
+   ,fQAJetEtaMin(0)    
+   ,fQAJetEtaMax(0)    
+   ,fQAJetNBinsPhi(0)  
+   ,fQAJetPhiMin(0)    
+   ,fQAJetPhiMax(0)    
+   ,fQATrackNBinsPt(0) 
+   ,fQATrackPtMin(0)   
+   ,fQATrackPtMax(0)   
+   ,fQATrackNBinsEta(0)
+   ,fQATrackEtaMin(0)  
+   ,fQATrackEtaMax(0)  
+   ,fQATrackNBinsPhi(0)
+   ,fQATrackPhiMin(0)  
+   ,fQATrackPhiMax(0)
+   ,fIJNBinsJetPt(0)
+   ,fIJJetPtMin(0)
+   ,fIJJetPtMax(0)
+   ,fIJNBinsPt(0)
+   ,fIJPtMin(0)
+   ,fIJPtMax(0)
+   ,fIJNBinsZ(0)
+   ,fIJZMin(0)
+   ,fIJZMax(0)
+   ,fIJNBinsCosTheta(0)
+   ,fIJCosThetaMin(0)
+   ,fIJCosThetaMax(0)
+   ,fIJNBinsTheta(0)
+   ,fIJThetaMin(0)
+   ,fIJThetaMax(0)
+   ,fIJNBinsJt(0)
+   ,fIJJtMin(0)
+   ,fIJJtMax(0)
+   ,fDiJetNBinsJetInvMass(0)
+   ,fDiJetJetInvMassMin(0)
+   ,fDiJetJetInvMassMax(0)
+   ,fDiJetNBinsJetPt(0)
+   ,fDiJetJetPtMin(0)
+   ,fDiJetJetPtMax(0)
+   ,fDiJetNBinsPt(0)
+   ,fDiJetPtMin(0)
+   ,fDiJetPtMax(0)
+   ,fDiJetNBinsXi(0)
+   ,fDiJetXiMin(0)
+   ,fDiJetXiMax(0)
+   ,fDiJetNBinsZ(0)
+   ,fDiJetZMin(0)
+   ,fDiJetZMax(0)
+   ,fQADiJetNBinsInvMass(0)
+   ,fQADiJetInvMassMin(0)
+   ,fQADiJetInvMassMax(0)
+   ,fQADiJetNBinsJetPt(0)
+   ,fQADiJetJetPtMin(0)
+   ,fQADiJetJetPtMax(0)
+   ,fQADiJetNBinsDeltaPhi(0)
+   ,fQADiJetDeltaPhiMin(0)
+   ,fQADiJetDeltaPhiMax(0)
+   ,fQADiJetNBinsDeltaEta(0)
+   ,fQADiJetDeltaEtaMin(0)
+   ,fQADiJetDeltaEtaMax(0)
+   ,fQADiJetNBinsDeltaPt(0)
+   ,fQADiJetDeltaPtMin(0)
+   ,fQADiJetDeltaPtMax(0)
+   ,fCommonHistList(0)
+   ,fh1EvtSelection(0)
+   ,fh1VertexNContributors(0)
+   ,fh1VertexZ(0)
+   ,fh1EvtMult(0)
+   ,fh1Xsec(0)
+   ,fh1Trials(0)
+   ,fh1PtHard(0)
+   ,fh1PtHardTrials(0)
+   ,fh1nRecJetsCuts(0)
+   ,fh1nGenJets(0)
+   ,fh1nRecEffJets(0)
+   ,fhnSingleTrackRecEffHisto(0)
+   ,fhnJetTrackRecEffHisto(0)
 {
-  //
-  // Default constructor
-  //
-  /*
-  for(int i = 0;i < kMaxStep*2;++i){
-    fhnJetContainer[i] = 0;
-  } 
-  */
-//   for(int ij  = 0;ij<kMaxJets;++ij){
-//     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
-//     fh1Eta[ij] = fh1Phi[ij] = 0;
-//   }
-  
-  DefineOutput(1, TList::Class());  
+   // default constructor
 }
 
-//////////////////////////////////////////////////////////////////////////////
-
-Bool_t AliAnalysisTaskFragmentationFunction::Notify()
+//__________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const char *name) 
+  : AliAnalysisTaskSE(name)
+  ,fESD(0)
+  ,fAOD(0)
+  ,fMCEvent(0)
+  ,fBranchRecJets("jets")
+  ,fBranchGenJets("")
+  ,fTrackTypeGen(0)
+  ,fJetTypeGen(0)
+  ,fJetTypeRecEff(0)
+  ,fFilterMask(0)
+  ,fUsePhysicsSelection(kTRUE)
+  ,fTrackPtCut(0)
+  ,fTrackEtaMin(0)
+  ,fTrackEtaMax(0)
+  ,fTrackPhiMin(0)
+  ,fTrackPhiMax(0)
+  ,fJetPtCut(0)
+  ,fJetEtaMin(0)
+  ,fJetEtaMax(0)
+  ,fJetPhiMin(0)
+  ,fJetPhiMax(0)
+  ,fDiJetCut(0)
+  ,fDiJetDeltaPhiCut(0)
+  ,fDiJetPtFractionCut(0)
+  ,fDiJetCDFCut(0)
+  ,fDiJetKindBins(0)
+  ,fFFRadius(0)
+  ,fTracksRec(0)
+  ,fTracksRecCuts(0)
+  ,fTracksGen(0)
+  ,fTracksAODMCCharged(0)
+  ,fTracksRecQualityCuts(0)
+  ,fJetsRec(0)
+  ,fJetsRecCuts(0)
+  ,fJetsGen(0)
+  ,fJetsRecEff(0)
+  ,fQATrackHistosRec(0)
+  ,fQATrackHistosRecCuts(0)
+  ,fQATrackHistosGen(0)
+  ,fQAJetHistosRec(0)
+  ,fQAJetHistosRecCuts(0)
+  ,fQAJetHistosRecCutsLeading(0)
+  ,fQAJetHistosGen(0)
+  ,fQAJetHistosGenLeading(0)
+  ,fQAJetHistosRecEffLeading(0)
+  ,fFFHistosRecCuts(0)
+  ,fFFHistosRecLeading(0)
+  ,fFFHistosRecLeadingTrack(0)
+  ,fFFHistosGen(0)
+  ,fFFHistosGenLeading(0)
+  ,fFFHistosGenLeadingTrack(0)
+  ,fIJHistosRecCuts(0)
+  ,fIJHistosRecLeading(0)
+  ,fIJHistosRecLeadingTrack(0)
+  ,fIJHistosGen(0)
+  ,fIJHistosGenLeading(0)
+  ,fIJHistosGenLeadingTrack(0)
+  ,fFFDiJetHistosRecCuts(0)
+  ,fFFDiJetHistosRecLeading(0)
+  ,fFFDiJetHistosRecLeadingTrack(0)
+  ,fFFDiJetHistosGen(0)
+  ,fFFDiJetHistosGenLeading(0)
+  ,fFFDiJetHistosGenLeadingTrack(0)
+  ,fQADiJetHistosRecCuts(0)
+  ,fQADiJetHistosGen(0)
+  ,fQATrackHighPtThreshold(0) 
+  ,fFFNBinsJetPt(0)    
+  ,fFFJetPtMin(0) 
+  ,fFFJetPtMax(0)
+  ,fFFNBinsPt(0)      
+  ,fFFPtMin(0)        
+  ,fFFPtMax(0)        
+  ,fFFNBinsXi(0)      
+  ,fFFXiMin(0)        
+  ,fFFXiMax(0)        
+  ,fFFNBinsZ(0)       
+  ,fFFZMin(0)         
+  ,fFFZMax(0)         
+  ,fQAJetNBinsPt(0)   
+  ,fQAJetPtMin(0)     
+  ,fQAJetPtMax(0)     
+  ,fQAJetNBinsEta(0)  
+  ,fQAJetEtaMin(0)    
+  ,fQAJetEtaMax(0)    
+  ,fQAJetNBinsPhi(0)  
+  ,fQAJetPhiMin(0)    
+  ,fQAJetPhiMax(0)    
+  ,fQATrackNBinsPt(0) 
+  ,fQATrackPtMin(0)   
+  ,fQATrackPtMax(0)   
+  ,fQATrackNBinsEta(0)
+  ,fQATrackEtaMin(0)  
+  ,fQATrackEtaMax(0)  
+  ,fQATrackNBinsPhi(0)
+  ,fQATrackPhiMin(0)  
+  ,fQATrackPhiMax(0)  
+  ,fIJNBinsJetPt(0)
+  ,fIJJetPtMin(0)
+  ,fIJJetPtMax(0)
+  ,fIJNBinsPt(0)
+  ,fIJPtMin(0)
+  ,fIJPtMax(0)
+  ,fIJNBinsZ(0)
+  ,fIJZMin(0)
+  ,fIJZMax(0)
+  ,fIJNBinsCosTheta(0)
+  ,fIJCosThetaMin(0)
+  ,fIJCosThetaMax(0)
+  ,fIJNBinsTheta(0)
+  ,fIJThetaMin(0)
+  ,fIJThetaMax(0)
+  ,fIJNBinsJt(0)
+  ,fIJJtMin(0)
+  ,fIJJtMax(0)
+  ,fDiJetNBinsJetInvMass(0)
+  ,fDiJetJetInvMassMin(0)
+  ,fDiJetJetInvMassMax(0)
+  ,fDiJetNBinsJetPt(0)
+  ,fDiJetJetPtMin(0)
+  ,fDiJetJetPtMax(0)
+  ,fDiJetNBinsPt(0)
+  ,fDiJetPtMin(0)
+  ,fDiJetPtMax(0)
+  ,fDiJetNBinsXi(0)
+  ,fDiJetXiMin(0)
+  ,fDiJetXiMax(0)
+  ,fDiJetNBinsZ(0)
+  ,fDiJetZMin(0)
+  ,fDiJetZMax(0)
+  ,fQADiJetNBinsInvMass(0)
+  ,fQADiJetInvMassMin(0)
+  ,fQADiJetInvMassMax(0)
+  ,fQADiJetNBinsJetPt(0)
+  ,fQADiJetJetPtMin(0)
+  ,fQADiJetJetPtMax(0)
+  ,fQADiJetNBinsDeltaPhi(0)
+  ,fQADiJetDeltaPhiMin(0)
+  ,fQADiJetDeltaPhiMax(0)
+  ,fQADiJetNBinsDeltaEta(0)
+  ,fQADiJetDeltaEtaMin(0)
+  ,fQADiJetDeltaEtaMax(0)
+  ,fQADiJetNBinsDeltaPt(0)
+  ,fQADiJetDeltaPtMin(0)
+  ,fQADiJetDeltaPtMax(0)
+  ,fCommonHistList(0)
+  ,fh1EvtSelection(0)
+  ,fh1VertexNContributors(0)
+  ,fh1VertexZ(0)
+  ,fh1EvtMult(0)
+  ,fh1Xsec(0)
+  ,fh1Trials(0)
+  ,fh1PtHard(0)
+  ,fh1PtHardTrials(0)
+  ,fh1nRecJetsCuts(0)
+  ,fh1nGenJets(0)
+  ,fh1nRecEffJets(0)
+  ,fhnSingleTrackRecEffHisto(0)
+  ,fhnJetTrackRecEffHisto(0)
 {
-  //
-  // Implemented Notify() to read the cross sections
-  // and number of trials from pyxsec.root
-  // 
+  // constructor
+  
+  DefineOutput(1,TList::Class());
+  
 
-//   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
-//   UInt_t ntrials  = 0;
-//   Float_t ftrials  = 0;
-//   if(tree){
-//     TFile *curfile = tree->GetCurrentFile();
-//     if (!curfile) {
-//       Error("Notify","No current file");
-//       return kFALSE;
-//     }
-
-//     if(!fh1Xsec||!fh1Trials){
-//       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
-//       return kFALSE;
-//     }
-
-//     TString fileName(curfile->GetName());
-//     if(fileName.Contains("AliESDs.root")){
-//         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
-//     }
-//     else if(fileName.Contains("AliAOD.root")){
-//         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
-//     }
-//     else if(fileName.Contains("AliAODs.root")){
-//         fileName.ReplaceAll("AliAODs.root", "");
-//     }
-//     else if(fileName.Contains("galice.root")){
-//         // for running with galice and kinematics alone...                      
-//         fileName.ReplaceAll("galice.root", "pyxsec.root");
-//     }
-//     TFile *fxsec = TFile::Open(fileName.Data());
-//     if(!fxsec){
-//       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
-//       // no a severe condition
-//       return kTRUE;
-//     }
-//     TTree *xtree = (TTree*)fxsec->Get("Xsection");
-//     if(!xtree){
-//       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
-//     }
-//     xtree->SetBranchAddress("xsection",&fXsection);
-//     xtree->SetBranchAddress("ntrials",&ntrials);
-//     ftrials = ntrials;
-//     xtree->GetEntry(0);
-
-//     fh1Xsec->Fill("<#sigma>",fXsection);
-//     fh1Trials->Fill("#sum{ntrials}",ftrials);
-
-//   }
+}
 
-  return kTRUE;
+//__________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const  AliAnalysisTaskFragmentationFunction &copy)
+  : AliAnalysisTaskSE()
+  ,fESD(copy.fESD)
+  ,fAOD(copy.fAOD)
+  ,fMCEvent(copy.fMCEvent)
+  ,fBranchRecJets(copy.fBranchRecJets)
+  ,fBranchGenJets(copy.fBranchGenJets)
+  ,fTrackTypeGen(copy.fTrackTypeGen)
+  ,fJetTypeGen(copy.fJetTypeGen)
+  ,fJetTypeRecEff(copy.fJetTypeRecEff)
+  ,fFilterMask(copy.fFilterMask)
+  ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
+  ,fTrackPtCut(copy.fTrackPtCut)
+  ,fTrackEtaMin(copy.fTrackEtaMin)
+  ,fTrackEtaMax(copy.fTrackEtaMax)
+  ,fTrackPhiMin(copy.fTrackPhiMin)
+  ,fTrackPhiMax(copy.fTrackPhiMax)
+  ,fJetPtCut(copy.fJetPtCut)
+  ,fJetEtaMin(copy.fJetEtaMin)
+  ,fJetEtaMax(copy.fJetEtaMax)
+  ,fJetPhiMin(copy.fJetPhiMin)
+  ,fJetPhiMax(copy.fJetPhiMax)
+  ,fDiJetCut(copy.fDiJetCut)
+  ,fDiJetDeltaPhiCut(copy.fDiJetDeltaPhiCut)
+  ,fDiJetPtFractionCut(copy.fDiJetPtFractionCut)
+  ,fDiJetCDFCut(copy.fDiJetCDFCut)
+  ,fDiJetKindBins(copy.fDiJetKindBins)
+  ,fFFRadius(copy.fFFRadius)
+  ,fTracksRec(copy.fTracksRec)
+  ,fTracksRecCuts(copy.fTracksRecCuts)
+  ,fTracksGen(copy.fTracksGen)
+  ,fTracksAODMCCharged(copy.fTracksAODMCCharged)
+  ,fTracksRecQualityCuts(copy.fTracksRecQualityCuts)
+  ,fJetsRec(copy.fJetsRec)
+  ,fJetsRecCuts(copy.fJetsRecCuts)
+  ,fJetsGen(copy.fJetsGen)
+  ,fJetsRecEff(copy.fJetsRecEff)
+  ,fQATrackHistosRec(copy.fQATrackHistosRec)
+  ,fQATrackHistosRecCuts(copy.fQATrackHistosRecCuts)
+  ,fQATrackHistosGen(copy.fQATrackHistosGen)
+  ,fQAJetHistosRec(copy.fQAJetHistosRec)
+  ,fQAJetHistosRecCuts(copy.fQAJetHistosRecCuts)
+  ,fQAJetHistosRecCutsLeading(copy.fQAJetHistosRecCutsLeading)
+  ,fQAJetHistosGen(copy.fQAJetHistosGen)
+  ,fQAJetHistosGenLeading(copy.fQAJetHistosGenLeading)
+  ,fQAJetHistosRecEffLeading(copy.fQAJetHistosRecEffLeading)
+  ,fFFHistosRecCuts(copy.fFFHistosRecCuts)
+  ,fFFHistosRecLeading(copy.fFFHistosRecLeading)
+  ,fFFHistosRecLeadingTrack(copy.fFFHistosRecLeadingTrack)
+  ,fFFHistosGen(copy.fFFHistosGen)
+  ,fFFHistosGenLeading(copy.fFFHistosGenLeading)
+  ,fFFHistosGenLeadingTrack(copy.fFFHistosGenLeadingTrack)
+  ,fIJHistosRecCuts(copy.fIJHistosRecCuts)
+  ,fIJHistosRecLeading(copy.fIJHistosRecLeading)
+  ,fIJHistosRecLeadingTrack(copy.fIJHistosRecLeadingTrack)
+  ,fIJHistosGen(copy.fIJHistosGen)
+  ,fIJHistosGenLeading(copy.fIJHistosGenLeading)
+  ,fIJHistosGenLeadingTrack(copy.fIJHistosGenLeadingTrack)
+  ,fFFDiJetHistosRecCuts(copy.fFFDiJetHistosRecCuts)
+  ,fFFDiJetHistosRecLeading(copy.fFFDiJetHistosRecLeading)
+  ,fFFDiJetHistosRecLeadingTrack(copy.fFFDiJetHistosRecLeadingTrack)
+  ,fFFDiJetHistosGen(copy.fFFDiJetHistosGen)
+  ,fFFDiJetHistosGenLeading(copy.fFFDiJetHistosGenLeading)
+  ,fFFDiJetHistosGenLeadingTrack(copy.fFFDiJetHistosGenLeadingTrack)
+  ,fQADiJetHistosRecCuts(copy.fQADiJetHistosRecCuts)
+  ,fQADiJetHistosGen(copy.fQADiJetHistosGen)
+  ,fQATrackHighPtThreshold(copy.fQATrackHighPtThreshold) 
+  ,fFFNBinsJetPt(copy.fFFNBinsJetPt)    
+  ,fFFJetPtMin(copy.fFFJetPtMin) 
+  ,fFFJetPtMax(copy.fFFJetPtMax)
+  ,fFFNBinsPt(copy.fFFNBinsPt)      
+  ,fFFPtMin(copy.fFFPtMin)        
+  ,fFFPtMax(copy.fFFPtMax)        
+  ,fFFNBinsXi(copy.fFFNBinsXi)      
+  ,fFFXiMin(copy.fFFXiMin)        
+  ,fFFXiMax(copy.fFFXiMax)        
+  ,fFFNBinsZ(copy.fFFNBinsZ)       
+  ,fFFZMin(copy.fFFZMin)         
+  ,fFFZMax(copy.fFFZMax)         
+  ,fQAJetNBinsPt(copy.fQAJetNBinsPt)   
+  ,fQAJetPtMin(copy.fQAJetPtMin)     
+  ,fQAJetPtMax(copy.fQAJetPtMax)     
+  ,fQAJetNBinsEta(copy.fQAJetNBinsEta)  
+  ,fQAJetEtaMin(copy.fQAJetEtaMin)    
+  ,fQAJetEtaMax(copy.fQAJetEtaMax)    
+  ,fQAJetNBinsPhi(copy.fQAJetNBinsPhi)  
+  ,fQAJetPhiMin(copy.fQAJetPhiMin)    
+  ,fQAJetPhiMax(copy.fQAJetPhiMax)    
+  ,fQATrackNBinsPt(copy.fQATrackNBinsPt) 
+  ,fQATrackPtMin(copy.fQATrackPtMin)   
+  ,fQATrackPtMax(copy.fQATrackPtMax)   
+  ,fQATrackNBinsEta(copy.fQATrackNBinsEta)
+  ,fQATrackEtaMin(copy.fQATrackEtaMin)  
+  ,fQATrackEtaMax(copy.fQATrackEtaMax)  
+  ,fQATrackNBinsPhi(copy.fQATrackNBinsPhi)
+  ,fQATrackPhiMin(copy.fQATrackPhiMin)  
+  ,fQATrackPhiMax(copy.fQATrackPhiMax)
+  ,fIJNBinsJetPt(copy.fIJNBinsJetPt)
+  ,fIJJetPtMin(copy.fIJJetPtMin)
+  ,fIJJetPtMax(copy.fIJJetPtMax)
+  ,fIJNBinsPt(copy.fIJNBinsPt)
+  ,fIJPtMin(copy.fIJPtMin)
+  ,fIJPtMax(copy.fIJPtMax)
+  ,fIJNBinsZ(copy.fIJNBinsZ)
+  ,fIJZMin(copy.fIJZMin)
+  ,fIJZMax(copy.fIJZMax)
+  ,fIJNBinsCosTheta(copy.fIJNBinsCosTheta)
+  ,fIJCosThetaMin(copy.fIJCosThetaMin)
+  ,fIJCosThetaMax(copy.fIJCosThetaMax)
+  ,fIJNBinsTheta(copy.fIJNBinsTheta)
+  ,fIJThetaMin(copy.fIJThetaMin)
+  ,fIJThetaMax(copy.fIJThetaMax)
+  ,fIJNBinsJt(copy.fIJNBinsJt)
+  ,fIJJtMin(copy.fIJJtMin)
+  ,fIJJtMax(copy.fIJJtMax)
+  ,fDiJetNBinsJetInvMass(copy.fDiJetNBinsJetInvMass)
+  ,fDiJetJetInvMassMin(copy.fDiJetJetInvMassMin)
+  ,fDiJetJetInvMassMax(copy.fDiJetJetInvMassMax)
+  ,fDiJetNBinsJetPt(copy.fDiJetNBinsJetPt)
+  ,fDiJetJetPtMin(copy.fDiJetJetPtMin)
+  ,fDiJetJetPtMax(copy.fDiJetJetPtMax)
+  ,fDiJetNBinsPt(copy.fDiJetNBinsPt)
+  ,fDiJetPtMin(copy.fDiJetPtMin)
+  ,fDiJetPtMax(copy.fDiJetPtMax)
+  ,fDiJetNBinsXi(copy.fDiJetNBinsXi)
+  ,fDiJetXiMin(copy.fDiJetXiMin)
+  ,fDiJetXiMax(copy.fDiJetXiMax)
+  ,fDiJetNBinsZ(copy.fDiJetNBinsZ)
+  ,fDiJetZMin(copy.fDiJetZMin)
+  ,fDiJetZMax(copy.fDiJetZMax)
+  ,fQADiJetNBinsInvMass(copy.fQADiJetNBinsInvMass)
+  ,fQADiJetInvMassMin(copy.fQADiJetInvMassMin)
+  ,fQADiJetInvMassMax(copy.fQADiJetInvMassMax)
+  ,fQADiJetNBinsJetPt(copy.fQADiJetNBinsJetPt)
+  ,fQADiJetJetPtMin(copy.fQADiJetJetPtMin)
+  ,fQADiJetJetPtMax(copy.fQADiJetJetPtMax)
+  ,fQADiJetNBinsDeltaPhi(copy.fQADiJetNBinsDeltaPhi)
+  ,fQADiJetDeltaPhiMin(copy.fQADiJetDeltaPhiMin)
+  ,fQADiJetDeltaPhiMax(copy.fQADiJetDeltaPhiMax)
+  ,fQADiJetNBinsDeltaEta(copy.fQADiJetNBinsDeltaEta)
+  ,fQADiJetDeltaEtaMin(copy.fQADiJetDeltaEtaMin)
+  ,fQADiJetDeltaEtaMax(copy.fQADiJetDeltaEtaMax)
+  ,fQADiJetNBinsDeltaPt(copy.fQADiJetNBinsDeltaPt)
+  ,fQADiJetDeltaPtMin(copy.fQADiJetDeltaPtMin)
+  ,fQADiJetDeltaPtMax(copy.fQADiJetDeltaPtMax)
+  ,fCommonHistList(copy.fCommonHistList)
+  ,fh1EvtSelection(copy.fh1EvtSelection)
+  ,fh1VertexNContributors(copy.fh1VertexNContributors)
+  ,fh1VertexZ(copy.fh1VertexZ)
+  ,fh1EvtMult(copy.fh1EvtMult)
+  ,fh1Xsec(copy.fh1Xsec)
+  ,fh1Trials(copy.fh1Trials)
+  ,fh1PtHard(copy.fh1PtHard)  
+  ,fh1PtHardTrials(copy.fh1PtHardTrials)  
+  ,fh1nRecJetsCuts(copy.fh1nRecJetsCuts)
+  ,fh1nGenJets(copy.fh1nGenJets)
+  ,fh1nRecEffJets(copy.fh1nRecEffJets)
+  ,fhnSingleTrackRecEffHisto(copy.fhnSingleTrackRecEffHisto)
+  ,fhnJetTrackRecEffHisto(copy.fhnJetTrackRecEffHisto)
+{
+  // copy constructor
 
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+// _________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::operator=(const AliAnalysisTaskFragmentationFunction& o)
+{
+  // assignment
+  
+  if(this!=&o){
+
+    AliAnalysisTaskSE::operator=(o);
+    fESD                          = o.fESD;
+    fAOD                          = o.fAOD;
+    fMCEvent                      = o.fMCEvent;
+    fBranchRecJets                = o.fBranchRecJets;
+    fBranchGenJets                = o.fBranchGenJets;
+    fTrackTypeGen                 = o.fTrackTypeGen;
+    fJetTypeGen                   = o.fJetTypeGen;
+    fJetTypeRecEff                = o.fJetTypeRecEff;
+    fFilterMask                   = o.fFilterMask;
+    fUsePhysicsSelection          = o.fUsePhysicsSelection;
+    fTrackPtCut                   = o.fTrackPtCut;
+    fTrackEtaMin                  = o.fTrackEtaMin;
+    fTrackEtaMax                  = o.fTrackEtaMax;
+    fTrackPhiMin                  = o.fTrackPhiMin;
+    fTrackPhiMax                  = o.fTrackPhiMax;
+    fJetPtCut                     = o.fJetPtCut;
+    fJetEtaMin                    = o.fJetEtaMin;
+    fJetEtaMax                    = o.fJetEtaMax;
+    fJetPhiMin                    = o.fJetPhiMin;
+    fJetPhiMax                    = o.fJetPhiMin;
+    fDiJetCut                     = o.fDiJetCut;
+    fDiJetDeltaPhiCut             = o.fDiJetDeltaPhiCut;
+    fDiJetPtFractionCut           = o.fDiJetPtFractionCut;
+    fDiJetCDFCut                  = o.fDiJetCDFCut;
+    fDiJetKindBins                = o.fDiJetKindBins;
+    fFFRadius                     = o.fFFRadius;
+    fTracksRec                    = o.fTracksRec;
+    fTracksRecCuts                = o.fTracksRecCuts;
+    fTracksGen                    = o.fTracksGen;
+    fTracksAODMCCharged           = o.fTracksAODMCCharged;
+    fTracksRecQualityCuts         = o.fTracksRecQualityCuts;
+    fJetsRec                      = o.fJetsRec;
+    fJetsRecCuts                  = o.fJetsRecCuts;
+    fJetsGen                      = o.fJetsGen;
+    fJetsRecEff                   = o.fJetsRecEff;
+    fQATrackHistosRec             = o.fQATrackHistosRec;
+    fQATrackHistosRecCuts         = o.fQATrackHistosRecCuts;
+    fQATrackHistosGen             = o.fQATrackHistosGen;
+    fQAJetHistosRec               = o.fQAJetHistosRec;
+    fQAJetHistosRecCuts           = o.fQAJetHistosRecCuts;
+    fQAJetHistosRecCutsLeading    = o.fQAJetHistosRecCutsLeading;
+    fQAJetHistosGen               = o.fQAJetHistosGen;
+    fQAJetHistosGenLeading        = o.fQAJetHistosGenLeading;
+    fQAJetHistosRecEffLeading     = o.fQAJetHistosRecEffLeading;
+    fFFHistosRecCuts              = o.fFFHistosRecCuts;
+    fFFHistosRecLeading           = o.fFFHistosRecLeading;
+    fFFHistosRecLeadingTrack      = o.fFFHistosRecLeadingTrack;
+    fFFHistosGen                  = o.fFFHistosGen;
+    fFFHistosGenLeading           = o.fFFHistosGenLeading;
+    fFFHistosGenLeadingTrack      = o.fFFHistosGenLeadingTrack;
+    fIJHistosRecCuts              = o.fIJHistosRecCuts;
+    fIJHistosRecLeading           = o.fIJHistosRecLeading;
+    fIJHistosRecLeadingTrack      = o.fIJHistosRecLeadingTrack;
+    fIJHistosGen                  = o.fIJHistosGen;
+    fIJHistosGenLeading           = o.fIJHistosGenLeading;
+    fIJHistosGenLeadingTrack      = o.fIJHistosGenLeadingTrack;
+    fFFDiJetHistosRecCuts         = o.fFFDiJetHistosRecCuts;
+    fFFDiJetHistosRecLeading      = o.fFFDiJetHistosRecLeading;
+    fFFDiJetHistosRecLeadingTrack = o.fFFDiJetHistosRecLeadingTrack;
+    fFFDiJetHistosGen             = o.fFFDiJetHistosGen;
+    fFFDiJetHistosGenLeading      = o.fFFDiJetHistosGenLeading;
+    fFFDiJetHistosGenLeadingTrack = o.fFFDiJetHistosGenLeadingTrack;
+    fQADiJetHistosRecCuts         = o.fQADiJetHistosRecCuts;
+    fQADiJetHistosGen             = o.fQADiJetHistosGen;
+    fQATrackHighPtThreshold       = o.fQATrackHighPtThreshold; 
+    fFFNBinsJetPt                 = o.fFFNBinsJetPt;    
+    fFFJetPtMin                   = o.fFFJetPtMin; 
+    fFFJetPtMax                   = o.fFFJetPtMax;
+    fFFNBinsPt                    = o.fFFNBinsPt;      
+    fFFPtMin                      = o.fFFPtMin;        
+    fFFPtMax                      = o.fFFPtMax;        
+    fFFNBinsXi                    = o.fFFNBinsXi;      
+    fFFXiMin                      = o.fFFXiMin;        
+    fFFXiMax                      = o.fFFXiMax;        
+    fFFNBinsZ                     = o.fFFNBinsZ;       
+    fFFZMin                       = o.fFFZMin;         
+    fFFZMax                       = o.fFFZMax;         
+    fQAJetNBinsPt                 = o.fQAJetNBinsPt;   
+    fQAJetPtMin                   = o.fQAJetPtMin;     
+    fQAJetPtMax                   = o.fQAJetPtMax;     
+    fQAJetNBinsEta                = o.fQAJetNBinsEta;  
+    fQAJetEtaMin                  = o.fQAJetEtaMin;    
+    fQAJetEtaMax                  = o.fQAJetEtaMax;    
+    fQAJetNBinsPhi                = o.fQAJetNBinsPhi;  
+    fQAJetPhiMin                  = o.fQAJetPhiMin;    
+    fQAJetPhiMax                  = o.fQAJetPhiMax;    
+    fQATrackNBinsPt               = o.fQATrackNBinsPt; 
+    fQATrackPtMin                 = o.fQATrackPtMin;   
+    fQATrackPtMax                 = o.fQATrackPtMax;   
+    fQATrackNBinsEta              = o.fQATrackNBinsEta;
+    fQATrackEtaMin                = o.fQATrackEtaMin;  
+    fQATrackEtaMax                = o.fQATrackEtaMax;  
+    fQATrackNBinsPhi              = o.fQATrackNBinsPhi;
+    fQATrackPhiMin                = o.fQATrackPhiMin;  
+    fQATrackPhiMax                = o.fQATrackPhiMax;  
+    fIJNBinsJetPt                 = o.fIJNBinsJetPt;
+    fIJJetPtMin                   = o.fIJJetPtMin;
+    fIJJetPtMax                   = o.fIJJetPtMax;
+    fIJNBinsPt                    = o.fIJNBinsPt;
+    fIJPtMin                      = o.fIJPtMin;
+    fIJPtMax                      = o.fIJPtMax;
+    fIJNBinsZ                     = o.fIJNBinsZ;
+    fIJZMin                       = o.fIJZMin;
+    fIJZMax                       = o.fIJZMax;
+    fIJNBinsCosTheta              = o.fIJNBinsCosTheta;
+    fIJCosThetaMin                = o.fIJCosThetaMin;
+    fIJCosThetaMax                = o.fIJCosThetaMax;
+    fIJNBinsTheta                 = o.fIJNBinsTheta;
+    fIJThetaMin                   = o.fIJThetaMin;
+    fIJThetaMax                   = o.fIJThetaMax;
+    fIJNBinsJt                    = o.fIJNBinsJt;
+    fIJJtMin                      = o.fIJJtMin;
+    fIJJtMax                      = o.fIJJtMax;
+    fDiJetNBinsJetInvMass         = o.fDiJetNBinsJetInvMass;
+    fDiJetJetInvMassMin           = o.fDiJetJetInvMassMin;
+    fDiJetJetInvMassMax           = o.fDiJetJetInvMassMax;
+    fDiJetNBinsJetPt              = o.fDiJetNBinsJetPt;
+    fDiJetJetPtMin                = o.fDiJetJetPtMin;
+    fDiJetJetPtMax                = o.fDiJetJetPtMax;
+    fDiJetNBinsPt                 = o.fDiJetNBinsPt;
+    fDiJetPtMin                   = o.fDiJetPtMin;
+    fDiJetPtMax                   = o.fDiJetPtMax;
+    fDiJetNBinsXi                 = o.fDiJetNBinsXi;
+    fDiJetXiMin                   = o.fDiJetXiMin;
+    fDiJetXiMax                   = o.fDiJetXiMax;
+    fDiJetNBinsZ                  = o.fDiJetNBinsZ;
+    fDiJetZMin                    = o.fDiJetZMin;
+    fDiJetZMax                    = o.fDiJetZMax;
+    fQADiJetNBinsInvMass          = o.fQADiJetNBinsInvMass;
+    fQADiJetInvMassMin            = o.fQADiJetInvMassMin;
+    fQADiJetInvMassMax            = o.fQADiJetInvMassMax;
+    fQADiJetNBinsJetPt            = o.fQADiJetNBinsJetPt;
+    fQADiJetJetPtMin              = o.fQADiJetJetPtMin;
+    fQADiJetJetPtMax              = o.fQADiJetJetPtMax;
+    fQADiJetNBinsDeltaPhi         = o.fQADiJetNBinsDeltaPhi;
+    fQADiJetDeltaPhiMin           = o.fQADiJetDeltaPhiMin;
+    fQADiJetDeltaPhiMax           = o.fQADiJetDeltaPhiMax;
+    fQADiJetNBinsDeltaEta         = o.fQADiJetNBinsDeltaEta;
+    fQADiJetDeltaEtaMin           = o.fQADiJetDeltaEtaMin;
+    fQADiJetDeltaEtaMax           = o.fQADiJetDeltaEtaMax;
+    fQADiJetNBinsDeltaPt          = o.fQADiJetNBinsDeltaPt;
+    fQADiJetDeltaPtMin            = o.fQADiJetDeltaPtMin;
+    fQADiJetDeltaPtMax            = o.fQADiJetDeltaPtMax;
+    fCommonHistList               = o.fCommonHistList;
+    fh1EvtSelection               = o.fh1EvtSelection;
+    fh1VertexNContributors        = o.fh1VertexNContributors;
+    fh1VertexZ                    = o.fh1VertexZ;
+    fh1EvtMult                    = o.fh1EvtMult;
+    fh1Xsec                       = o.fh1Xsec;
+    fh1Trials                     = o.fh1Trials;
+    fh1PtHard                     = o.fh1PtHard;
+    fh1PtHardTrials               = o.fh1PtHardTrials;
+    fh1nRecJetsCuts               = o.fh1nRecJetsCuts;
+    fh1nGenJets                   = o.fh1nGenJets; 
+    fh1nRecEffJets                = o.fh1nRecEffJets;
+    fhnSingleTrackRecEffHisto     = o.fhnSingleTrackRecEffHisto;
+    fhnJetTrackRecEffHisto        = o.fhnJetTrackRecEffHisto;
+  }
+    
+  return *this;
+}
 
-void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
+//___________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::~AliAnalysisTaskFragmentationFunction()
 {
-  //
-  // Create the output container
-  //
+  // destructor
+  
+  if(fTracksRec)            delete fTracksRec;
+  if(fTracksRecCuts)        delete fTracksRecCuts;
+  if(fTracksGen)            delete fTracksGen;
+  if(fTracksAODMCCharged)   delete fTracksAODMCCharged;  
+  if(fTracksRecQualityCuts) delete fTracksRecQualityCuts; 
+  if(fJetsRec)              delete fJetsRec;
+  if(fJetsRecCuts)          delete fJetsRecCuts;
+  if(fJetsGen)              delete fJetsGen;
+  if(fJetsRecEff)           delete fJetsRecEff;
+
+  //  if(fDiJetBins)     delete fDiJetBins;
 
-  //**** Connect the AOD
-  if(fUseAODInput) // Use AODs as input not ESDs
-  { 
-    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
-    if(!fAOD)
-    {
-      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
-      return;
-    }
+}
 
-    // fetch the header
-    fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
-    if(!fJetHeaderRec)
-    {
-      Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
-    }
-  }
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const char* name, 
+                                                        Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
+                                                        Int_t nPt, Float_t ptMin, Float_t ptMax,
+                                                        Int_t nXi, Float_t xiMin, Float_t xiMax,
+                                                        Int_t nZ , Float_t zMin , Float_t zMax )
+  : TObject()
+  ,fNBinsJetPt(nJetPt)
+  ,fJetPtMin(jetPtMin)
+  ,fJetPtMax(jetPtMax)
+  ,fNBinsPt(nPt) 
+  ,fPtMin(ptMin)   
+  ,fPtMax(ptMax)   
+  ,fNBinsXi(nXi) 
+  ,fXiMin(xiMin)   
+  ,fXiMax(xiMax)   
+  ,fNBinsZ(nZ)  
+  ,fZMin(zMin)    
+  ,fZMax(zMax)    
+  ,fh2TrackPt(0)
+  ,fh2Xi(0)
+  ,fh2Z(0)
+  ,fh1JetPt(0)
+  ,fName(name)
+{
+  // default constructor
 
+}
 
-  else // Use the AOD on the flight 
-  {
-    //  assume that the AOD is in the general output...
-    fAOD  = AODEvent();
-    if(!fAOD)
-    {
-      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
-      return;
-    }
+//___________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const AliFragFuncHistos& copy)
+  : TObject()
+  ,fNBinsJetPt(copy.fNBinsJetPt)
+  ,fJetPtMin(copy.fJetPtMin)
+  ,fJetPtMax(copy.fJetPtMax)
+  ,fNBinsPt(copy.fNBinsPt) 
+  ,fPtMin(copy.fPtMin)   
+  ,fPtMax(copy.fPtMax)   
+  ,fNBinsXi(copy.fNBinsXi) 
+  ,fXiMin(copy.fXiMin)   
+  ,fXiMax(copy.fXiMax)   
+  ,fNBinsZ(copy.fNBinsZ)  
+  ,fZMin(copy.fZMin)    
+  ,fZMax(copy.fZMax)    
+  ,fh2TrackPt(copy.fh2TrackPt)
+  ,fh2Xi(copy.fh2Xi)
+  ,fh2Z(copy.fh2Z)
+  ,fh1JetPt(copy.fh1JetPt)
+  ,fName(copy.fName)
+{
+  // copy constructor
+}
 
-    //    ((TList*)OutputTree()->GetUserInfo())->Dump();
-    fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));    
-    if(!fJetHeaderRec)
-    {
-      Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
-    }
-    else
-    {
-      if(fDebug>10)fJetHeaderRec->Dump();
-    }
+//_______________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsJetPt = o.fNBinsJetPt;
+    fJetPtMin   = o.fJetPtMin;
+    fJetPtMax   = o.fJetPtMax;
+    fNBinsPt    = o.fNBinsPt; 
+    fPtMin      = o.fPtMin;   
+    fPtMax      = o.fPtMax;   
+    fNBinsXi    = o.fNBinsXi; 
+    fXiMin      = o.fXiMin;   
+    fXiMax      = o.fXiMax;   
+    fNBinsZ     = o.fNBinsZ;  
+    fZMin       = o.fZMin;    
+    fZMax       = o.fZMax;    
+    fh2TrackPt  = o.fh2TrackPt;
+    fh2Xi       = o.fh2Xi;
+    fh2Z        = o.fh2Z;
+    fh1JetPt    = o.fh1JetPt;
+    fName       = o.fName;
   }
+    
+  return *this;
+}
 
-////////////////////////////
-
-  if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
+//_________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::~AliFragFuncHistos()
+{
+  // destructor 
 
-  OpenFile(1);
-  if(!fHistList)fHistList = new TList();
-  fHistList->SetOwner(kTRUE);
-  Bool_t oldStatus = TH1::AddDirectoryStatus();
-  TH1::AddDirectory(kFALSE);
+  if(fh1JetPt)   delete fh1JetPt;
+  if(fh2TrackPt) delete fh2TrackPt;
+  if(fh2Xi)      delete fh2Xi;
+  if(fh2Z)       delete fh2Z;
+}
 
-//////////////////////////////////////////////////
-//////////// HISTOGRAM DECLARATION ///////////////
-//////////////////////////////////////////////////
+//_________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::DefineHistos()
+{
+  // book FF histos
 
-  DefineJetH();
+  fh1JetPt   = new TH1F(Form("fh1FFJetPt%s", fName.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
+  fh2TrackPt = new TH2F(Form("fh2FFTrackPt%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsPt, fPtMin, fPtMax);
+  fh2Xi      = new TH2F(Form("fh2FFXi%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsXi, fXiMin, fXiMax);
+  fh2Z       = new TH2F(Form("fh2FFZ%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, fZMin, fZMax);
 
-//////////////////////////////////////////////////
-////////////// HISTOGRAM SAVING //////////////////
-//////////////////////////////////////////////////
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries"); 
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPt,"jet p_{T} [GeV/c]","p_{T} [GeV/c]","entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi,"jet p_{T} [GeV/c]","#xi", "entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z,"jet p_{T} [GeV/c]","z","entries");
+}
 
-  for (Int_t i3 = 0; i3 < fnEBin; i3++)
-  {
-    fHistList->Add(fEtaMonoJet1H[i3]);
-    fHistList->Add(fPhiMonoJet1H[i3]);
-    fHistList->Add(fPtMonoJet1H[i3]);
-    fHistList->Add(fEMonoJet1H[i3]);
-   
-    for(Int_t i4 = 0; i4 < fnRBin; i4++)
-    {
-      fHistList->Add(fdNdXiMonoJet1H[i3][i4]);
-      fHistList->Add(fdNdPtMonoJet1H[i3][i4]);
-      fHistList->Add(fdNdZMonoJet1H[i3][i4]);
-      fHistList->Add(fNMonoJet1sH[i3][i4]);
-    }
-  }
+//_______________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::FillFF(Float_t trackPt, Float_t jetPt, Bool_t incrementJetPt)
+{
+  // fill FF
  
-  // Theta, kT particles/jet
-  for (Int_t i3 = 0; i3 < fnEBin; i3++)
-  {
-    for(Int_t i4 = 0; i4 < fnRBin; i4++)
-    {
-      fHistList->Add(fdNdThetaMonoJet1H[i3][i4]);
-      fHistList->Add(fdNdcosThetaMonoJet1H[i3][i4]);
-      fHistList->Add(fdNdkTMonoJet1H[i3][i4]);
-      fHistList->Add(fdNdpTvsZMonoJet1H[i3][i4]);
-      fHistList->Add(fShapeMonoJet1H[i3][i4]);
-         
-      fHistList->Add(fThetaPtPartMonoJet1H[i3][i4]);
-      fHistList->Add(fcosThetaPtPartMonoJet1H[i3][i4]);
-      fHistList->Add(fkTPtPartMonoJet1H[i3][i4]);
-      fHistList->Add(fThetaPtJetMonoJet1H[i3][i4]);
-      fHistList->Add(fcosThetaPtJetMonoJet1H[i3][i4]);
-      fHistList->Add(fkTPtJetMonoJet1H[i3][i4]);
-      fHistList->Add(fpTPtJetMonoJet1H[i3][i4]);
-    }
-  }
+  if(incrementJetPt) fh1JetPt->Fill(jetPt);    
+  fh2TrackPt->Fill(jetPt,trackPt);
   
-  // Track QA - Correlations
-  for (Int_t iPtBin=0; iPtBin<fnPtTrigBin; iPtBin++)
-    {
-      fHistList->Add(fptTracks[iPtBin]);
-      fHistList->Add(fetaTracks[iPtBin]);
-      fHistList->Add(fphiTracks[iPtBin]);
-      fHistList->Add(fdetaTracks[iPtBin]);
-      fHistList->Add(fdphiTracks[iPtBin]);
-      fHistList->Add(fetaphiptTracks[iPtBin]);
-      fHistList->Add(fetaphiTracks[iPtBin]);
-      fHistList->Add(fdetadphiTracks[iPtBin]);
-      fHistList->Add(fptTracksCut[iPtBin]);
-      fHistList->Add(fetaTracksCut[iPtBin]);
-      fHistList->Add(fphiTracksCut[iPtBin]);
-      fHistList->Add(fdetaTracksCut[iPtBin]);
-      fHistList->Add(fdphiTracksCut[iPtBin]);
-      fHistList->Add(fetaphiptTracksCut[iPtBin]);
-      fHistList->Add(fetaphiTracksCut[iPtBin]);
-      fHistList->Add(fdetadphiTracksCut[iPtBin]);
-      fHistList->Add(fNPtTrig[iPtBin]);
-      fHistList->Add(fNPtTrigCut[iPtBin]);
-    }
+  Double_t z = 0.;
+  if(jetPt>0) z = trackPt / jetPt;
+  Double_t xi = 0;
+  if(z>0) xi = TMath::Log(1/z);
+  
+  fh2Xi->Fill(jetPt,xi);
+  fh2Z->Fill(jetPt,z);
+}
 
-  // Track QA 
-  fHistList->Add(fptAllTracks);
-  fHistList->Add(fetaAllTracks);
-  fHistList->Add(fphiAllTracks);
-  fHistList->Add(fetaphiptAllTracks);
-  fHistList->Add(fetaphiAllTracks);
-  fHistList->Add(fptAllTracksCut);
-  fHistList->Add(fetaAllTracksCut);
-  fHistList->Add(fphiAllTracksCut);
-  fHistList->Add(fetaphiptAllTracksCut);
-  fHistList->Add(fetaphiAllTracksCut);
-
-  // Event caracterisation QA
-  fHistList->Add(fvertexXY);
-  fHistList->Add(fvertexZ);
-  fHistList->Add(fEvtMult);
-  fHistList->Add(fEvtMultvsJetPt);
-  fHistList->Add(fPtvsEtaJet);
-  fHistList->Add(fNpvsEtaJet);
-  fHistList->Add(fNpevtvsEtaJet);
-  fHistList->Add(fPtvsPtJet);
-  fHistList->Add(fNpvsPtJet);
-  fHistList->Add(fNpevtvsPtJet);
-  fHistList->Add(fPtvsPtJet1D);
-  fHistList->Add(fNpvsPtJet1D);
-  fHistList->Add(fNpevtvsPtJet1D);
-  fHistList->Add(fptLeadingJet);
-  fHistList->Add(fetaLeadingJet);
-  fHistList->Add(fphiLeadingJet);
-  fHistList->Add(fptJet);
-  fHistList->Add(fetaJet);
-  fHistList->Add(fphiJet);
-  fHistList->Add(fNBadRunsH);
-
-//////////////////////////////////////////////////
-///////// END OF HISTOGRAM DECLARATION ///////////
-//////////////////////////////////////////////////
+//_________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AddToOutput(TList* list) const
+{
+  // add histos to list
 
-  // =========== Switch on Sumw2 for all histos ===========
-  for (Int_t i=0; i<fHistList->GetEntries(); ++i) 
-  {
-    TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
-    if (h1)
-    {
-      // Printf("%s ",h1->GetName());
-      h1->Sumw2();
-      continue;
-    }
-  }
+  list->Add(fh1JetPt);
   
-  TH1::AddDirectory(oldStatus);
+  list->Add(fh2TrackPt);
+  list->Add(fh2Xi);
+  list->Add(fh2Z);
+}
+
+//_________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const char* name,
+                                                              Int_t nPt,   Float_t ptMin,   Float_t ptMax,
+                                                              Int_t nEta,  Float_t etaMin,  Float_t etaMax,
+                                                              Int_t nPhi,  Float_t phiMin,  Float_t phiMax)
+  : TObject()
+  ,fNBinsPt(nPt)
+  ,fPtMin(ptMin)
+  ,fPtMax(ptMax)
+  ,fNBinsEta(nEta)
+  ,fEtaMin(etaMin)
+  ,fEtaMax(etaMax)
+  ,fNBinsPhi(nPhi)
+  ,fPhiMin(phiMin)
+  ,fPhiMax(phiMax)
+  ,fh2EtaPhi(0)
+  ,fh1Pt(0)
+  ,fName(name)
+{
+  // default constructor
+}
+
+//____________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const AliFragFuncQAJetHistos& copy)
+  : TObject()
+  ,fNBinsPt(copy.fNBinsPt)
+  ,fPtMin(copy.fPtMin)
+  ,fPtMax(copy.fPtMax)
+  ,fNBinsEta(copy.fNBinsEta)
+  ,fEtaMin(copy.fEtaMin)
+  ,fEtaMax(copy.fEtaMax)
+  ,fNBinsPhi(copy.fNBinsPhi)
+  ,fPhiMin(copy.fPhiMin)
+  ,fPhiMax(copy.fPhiMax)
+  ,fh2EtaPhi(copy.fh2EtaPhi)
+  ,fh1Pt(copy.fh1Pt)
+  ,fName(copy.fName)
+{
+  // copy constructor
+}
+
+//________________________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& o)
+{
+  // assignment
   
-  /*
-    if (fDebug > 1) printf("AnalysisTaskDiJets::CreateOutPutData() \n");
-    fDiJets = new TClonesArray("AliAODDiJet", 0);
-    fDiJets->SetName("Dijets");
-    AddAODBranch("TClonesArray", &fDiJets);
-  */
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsPt  = o.fNBinsPt;
+    fPtMin    = o.fPtMin;
+    fPtMax    = o.fPtMax;
+    fNBinsEta = o.fNBinsEta;
+    fEtaMin   = o.fEtaMin;
+    fEtaMax   = o.fEtaMax;
+    fNBinsPhi = o.fNBinsPhi;
+    fPhiMin   = o.fPhiMin;
+    fPhiMax   = o.fPhiMax;
+    fh2EtaPhi = o.fh2EtaPhi;
+    fh1Pt     = o.fh1Pt;
+    fName     = o.fName;
+  }
   
+  return *this;
+}
+
+//______________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::~AliFragFuncQAJetHistos()
+{
+  // destructor 
   
+  if(fh2EtaPhi) delete fh2EtaPhi;
+  if(fh1Pt)     delete fh1Pt;
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+//____________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::DefineHistos()
+{
+  // book jet QA histos
 
-void AliAnalysisTaskFragmentationFunction::Init()
+  fh2EtaPhi  = new TH2F(Form("fh2JetQAEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
+  fh1Pt      = new TH1F(Form("fh1JetQAPt%s", fName.Data()), Form("%s: p_{T} distribution", fName.Data()), fNBinsPt, fPtMin, fPtMax);
+       
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi"); 
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
+}
+
+//____________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::FillJetQA(Float_t eta, Float_t phi, Float_t pt)
 {
-  //
-  // Initialization
-  //
+  // fill jet QA histos 
 
-  Printf(">>> AnalysisTaskFragmentationFunction::Init() debug level %d\n",fDebug);
-  if (fDebug > 1) printf("AnalysisTaskDiJets::Init() \n");
+  fh2EtaPhi->Fill( eta, phi);
+  fh1Pt->Fill( pt );
 }
 
-/////////////////////////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////////
+//____________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AddToOutput(TList* list) const 
+{
+  // add histos to list
+
+  list->Add(fh2EtaPhi);
+  list->Add(fh1Pt);
+}
 
-void AliAnalysisTaskFragmentationFunction::UserExec(Option_t */*option*/) 
+//___________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const char* name,
+                                                                  Int_t nPt, Float_t ptMin, Float_t ptMax,
+                                                                  Int_t nEta, Float_t etaMin, Float_t etaMax,
+                                                                  Int_t nPhi, Float_t phiMin, Float_t phiMax,
+                                                                  Float_t ptThresh) 
+  : TObject()
+  ,fNBinsPt(nPt)
+  ,fPtMin(ptMin)
+  ,fPtMax(ptMax)
+  ,fNBinsEta(nEta)
+  ,fEtaMin(etaMin)
+  ,fEtaMax(etaMax)
+  ,fNBinsPhi(nPhi)
+  ,fPhiMin(phiMin)
+  ,fPhiMax(phiMax)
+  ,fHighPtThreshold(ptThresh)
+  ,fh2EtaPhi(0)
+  ,fh1Pt(0)
+  ,fh2HighPtEtaPhi(0)
+  ,fName(name)
 {
-  //
-  // Execute analysis for current event
-  //
+  // default constructor
+}
 
-  //****
-  //**** Check of input data
-  //****
+//__________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const AliFragFuncQATrackHistos& copy)
+  : TObject()
+  ,fNBinsPt(copy.fNBinsPt)
+  ,fPtMin(copy.fPtMin)
+  ,fPtMax(copy.fPtMax)
+  ,fNBinsEta(copy.fNBinsEta)
+  ,fEtaMin(copy.fEtaMin)
+  ,fEtaMax(copy.fEtaMax)
+  ,fNBinsPhi(copy.fNBinsPhi)
+  ,fPhiMin(copy.fPhiMin)
+  ,fPhiMax(copy.fPhiMax)
+  ,fHighPtThreshold(copy.fHighPtThreshold)
+  ,fh2EtaPhi(copy.fh2EtaPhi)
+  ,fh1Pt(copy.fh1Pt)
+  ,fh2HighPtEtaPhi(copy.fh2HighPtEtaPhi)
+  ,fName(copy.fName)
+{
+  // copy constructor
+}
 
-  printf("Analysing event # %5d\n", (Int_t) fEntry);
-  if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
+// _____________________________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& o)
+{
+  // assignment
   
-  AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
-  if(!aodH)
-  {
-    Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
-    return;
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsPt         = o.fNBinsPt;
+    fPtMin           = o.fPtMin;
+    fPtMax           = o.fPtMax;
+    fNBinsEta        = o.fNBinsEta;
+    fEtaMin          = o.fEtaMin;
+    fEtaMax          = o.fEtaMax;
+    fNBinsPhi        = o.fNBinsPhi;
+    fPhiMin          = o.fPhiMin;
+    fPhiMax          = o.fPhiMax;
+    fHighPtThreshold = o.fHighPtThreshold;
+    fh2EtaPhi        = o.fh2EtaPhi;
+    fh1Pt            = o.fh1Pt;
+    fh2HighPtEtaPhi  = o.fh2HighPtEtaPhi;
+    fName            = o.fName;
   }
   
-  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
-  if(!aodRecJets)
-  {
-    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
-    return;
-  }
-  if (fDebug > 10) Printf("%s:%d",(char*)__FILE__,__LINE__);
-  
-  //****
-  //**** Check of primary vertex
-  //****
-  AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
-  if( (!pvtx) ||
-      (pvtx->GetZ()<-10. || pvtx->GetZ()>10.) ||
-      (pvtx->GetNContributors()<0) )
-    { 
-       fNBadRuns++;
-       fNBadRunsH->Fill(0.5);
-       return;
-    }
+  return *this;
+}
 
-  //****
-  //**** Check number of tracks 
-  //****
-  TClonesArray* tracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
-
-  //****
-  //**** Declaration of arrays and variables 
-  //****
-  // We use static array, not to fragment the memory
-  AliAODJet recJets[kMaxJets];
-  Int_t nRecJets       = 0;
-  Int_t nTracks = 0;
-
-//////////////////////////////////////////////////
-///////// Get the reconstructed jets /////////////
-//////////////////////////////////////////////////
-
-  nRecJets = aodRecJets->GetEntries();
-  nRecJets = TMath::Min(nRecJets,kMaxJets);
-  nTracks  = fAOD->GetNumberOfTracks();
-    
-  for(int ir = 0;ir < nRecJets;++ir)
-  {
-    AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
-    if(!tmp)continue;
-    recJets[ir] = *tmp;
-    cout << "recJets[" << ir << "].Eta(): " << recJets[ir].Eta() << ", recJets[" << ir <<"].Phi(): " << recJets[ir].Phi() << ", recJets[" << ir << "].E(): " << recJets[ir].E() << endl;
-  }
-  if(nRecJets>1) {
-    Float_t detaJ = recJets[0].Eta() - recJets[1].Eta();
-    Float_t dphiJ = recJets[0].Phi() - recJets[1].Phi();
-    Float_t detJ = recJets[0].Pt() - recJets[1].Pt();
-    cout << "detaJ: " << detaJ << ", dphiJ: " << dphiJ << ", detJ: " << detJ << endl;
-  }
+//___________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::~AliFragFuncQATrackHistos()
+{
+  // destructor 
+  
+  if(fh2EtaPhi)       delete fh2EtaPhi;
+  if(fh2HighPtEtaPhi) delete fh2HighPtEtaPhi;
+  if(fh1Pt)           delete fh1Pt;
+}
 
-  // Get vertex information
-  fvertexXY->Fill(pvtx->GetX(),pvtx->GetY());
-  fvertexZ->Fill(pvtx->GetZ());
+//______________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::DefineHistos()
+{
+  // book track QA histos
 
-  //////////////////////////////////////////////////
-  ////////////// TRACK QUALITY ASSURANCE ///////////           TO BE OPTIMISED!!
-  //////////////////////////////////////////////////
-  Int_t evtMult = 0;
-  for(Int_t it=0; it<nTracks; it++)
-  {
-    AliAODTrack* aodTrack = (AliAODTrack*)tracks->At(it);
-    Float_t etaT = aodTrack->Eta();
-    Float_t phiT = aodTrack->Phi();
-    Float_t ptT = aodTrack->Pt(); 
-    phiT = ((phiT < 0) ? phiT + 2 * TMath::Pi() : phiT);
-    //    cout << "etaT: " << etaT << ", phiT: " << phiT << endl;
-    fptAllTracks->Fill(ptT);
-    fetaAllTracks->Fill(etaT);
-    fphiAllTracks->Fill(phiT);
-    fetaphiptAllTracks->Fill(etaT,phiT,ptT);
-    fetaphiAllTracks->Fill(etaT,phiT,1);
-    UInt_t status = aodTrack->GetStatus();
-
-    for(Int_t i=0; i<fnPtTrigBin; i++)
-    {
-      if(ptT>=farrayPtTrigmin[i] && ptT<farrayPtTrigmax[i])
-      {
-       fptTracks[i]->Fill(ptT);
-       fetaTracks[i]->Fill(etaT);
-       fphiTracks[i]->Fill(phiT);
-       fNPtTrig[i]->Fill(0.5);
-       // Compute deta/dphi
-       Float_t etaT2 = 0.; 
-       Float_t phiT2 = 0.;
-       Float_t ptT2 = 0.;
-       Float_t detaT = 0.;
-       Float_t dphiT = 0.;
-       for(Int_t it2 = 0; it2< nTracks; it2++)
-       {
-          AliAODTrack* aodTrack2 = (AliAODTrack*)tracks->At(it2);
-          etaT2 = aodTrack2->Eta(); phiT2 = aodTrack2->Phi(); ptT2 = aodTrack2->Pt();
-          phiT2 = ((phiT2 < 0) ? phiT2 + 2 * TMath::Pi() : phiT2);
-          //         cout << "etaT2: " << etaT2 << ", phiT2: " << phiT2 << endl;
-          if(ptT2 > 2.*i+4.) continue;
-          if(it2==it) continue;
-          detaT = etaT - etaT2; 
-          dphiT = phiT - phiT2; 
-          if (dphiT  >   TMath::Pi()) dphiT = (-TMath::Pi() +TMath::Abs(dphiT - TMath::Pi()));
-          if (dphiT  < -1.0*TMath::Pi())  dphiT = (TMath::Pi() -  TMath::Abs(dphiT + TMath::Pi()));
-             
-          fdetaTracks[i]->Fill(detaT);
-          fdphiTracks[i]->Fill(dphiT);
-          fdetadphiTracks[i]->Fill(detaT,dphiT,1);
-       }
-        fetaphiptTracks[i]->Fill(etaT,phiT,ptT);
-       fetaphiTracks[i]->Fill(etaT,phiT,1);
-      }
-    } // End loop over trigger ranges 
+  fh2EtaPhi       = new TH2F(Form("fh2TrackQAEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
+  fh2HighPtEtaPhi = new TH2F(Form("fh2TrackQAHighPtEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution for high-p_{T}", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
+  fh1Pt           = new TH1F(Form("fh1TrackQAPt%s", fName.Data()), Form("%s: p_{T} distribution", fName.Data()), fNBinsPt, fPtMin, fPtMax);
+  
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi"); 
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2HighPtEtaPhi, "#eta", "#phi");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
+}
 
-    if (status == 0) continue;
-    if((fFilterMask>0)&&!(aodTrack->TestFilterBit(fFilterMask))) continue;
-    fptAllTracksCut->Fill(ptT);
-    fetaAllTracksCut->Fill(etaT);
-    fphiAllTracksCut->Fill(phiT);
-    fetaphiptAllTracksCut->Fill(etaT,phiT,ptT);
-    fetaphiAllTracksCut->Fill(etaT,phiT,1);
-    if(ptT > 0.150 && TMath::Abs(etaT) < 0.9) evtMult++;
-  } // end loop over tracks
-  fEvtMult->Fill(evtMult);   
+//________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::FillTrackQA(Float_t eta, Float_t phi, Float_t pt)
+{
+  // fill track QA histos
+    
+  fh2EtaPhi->Fill( eta, phi);
+  if(pt > fHighPtThreshold) fh2HighPtEtaPhi->Fill( eta, phi);
+  fh1Pt->Fill( pt );   
+}
 
-//////////////////////////////////////////////////
-///////////////// MONOJET PART ///////////////////
-//////////////////////////////////////////////////
+//______________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AddToOutput(TList* list) const
+{
+  // add histos to list
 
-  if (nRecJets == 0) return;
+  list->Add(fh2EtaPhi);
+  list->Add(fh2HighPtEtaPhi);
+  list->Add(fh1Pt);
+}
 
-  Double_t jetEnergy = recJets[0].E();
-  Int_t    goodBin   = 999;
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AliFragFuncIntraJetHistos(const char* name, 
+                                                        Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
+                                                        Int_t nPt, Float_t ptMin, Float_t ptMax,
+                                                        Int_t nZ , Float_t zMin , Float_t zMax,
+                                                        Int_t nCosTheta , Float_t costhetaMin , Float_t costhetaMax,
+                                                        Int_t nTheta , Float_t thetaMin , Float_t thetaMax,
+                                                        Int_t nJt , Float_t jtMin , Float_t jtMax)
+  : TObject()
+  ,fNBinsJetPt(nJetPt)
+  ,fJetPtMin(jetPtMin)
+  ,fJetPtMax(jetPtMax)
+  ,fNBinsPt(nPt) 
+  ,fPtMin(ptMin)   
+  ,fPtMax(ptMax)   
+  ,fNBinsZ(nZ) 
+  ,fZMin(zMin)   
+  ,fZMax(zMax)   
+  ,fNBinsJt(nJt)
+  ,fJtMin(jtMin)
+  ,fJtMax(jtMax)
+  ,fNBinsTheta(nTheta)
+  ,fThetaMin(thetaMin)
+  ,fThetaMax(thetaMax)
+  ,fNBinsCosTheta(nCosTheta)
+  ,fCosThetaMin(costhetaMin)
+  ,fCosThetaMax(costhetaMax)
+  ,fh2Theta(0)
+  ,fh2CosTheta(0)
+  ,fh2Jt(0)
+  ,fh2PtvsZ(0)
+  ,fhnIntraJet(0)
+  ,fnDim(6)
+  ,fName(name)
+{
+  // default constructor
 
-  for (Int_t i1 = 0; i1 < fnEBin; i1++)
-  {
-    if (jetEnergy < farrayEmax[i1] && jetEnergy >= farrayEmin[i1]) 
-    {
-      goodBin = i1;
-      continue;
-    }
-  }
+}
 
-  fptLeadingJet->Fill(recJets[0].Pt());
-  fetaLeadingJet->Fill(recJets[0].Eta());
-  fphiLeadingJet->Fill(recJets[0].Phi());
+//___________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AliFragFuncIntraJetHistos(const AliFragFuncIntraJetHistos& copy)
+  : TObject()
+  ,fNBinsJetPt(copy.fNBinsJetPt)
+  ,fJetPtMin(copy.fJetPtMin)
+  ,fJetPtMax(copy.fJetPtMax)
+  ,fNBinsPt(copy.fNBinsPt) 
+  ,fPtMin(copy.fPtMin)   
+  ,fPtMax(copy.fPtMax)   
+  ,fNBinsZ(copy.fNBinsZ) 
+  ,fZMin(copy.fZMin)   
+  ,fZMax(copy.fZMax)   
+  ,fNBinsJt(copy.fNBinsJt)
+  ,fJtMin(copy.fJtMin)
+  ,fJtMax(copy.fJtMax)
+  ,fNBinsTheta(copy.fNBinsTheta)
+  ,fThetaMin(copy.fThetaMin)
+  ,fThetaMax(copy.fThetaMax)
+  ,fNBinsCosTheta(copy.fNBinsCosTheta)
+  ,fCosThetaMin(copy.fCosThetaMin)
+  ,fCosThetaMax(copy.fCosThetaMax)
+  ,fh2Theta(copy.fh2Theta)
+  ,fh2CosTheta(copy.fh2CosTheta)
+  ,fh2Jt(copy.fh2Jt)
+  ,fh2PtvsZ(copy.fh2PtvsZ)
+  ,fhnIntraJet(copy.fhnIntraJet)
+  ,fnDim(copy.fnDim)
+  ,fName(copy.fName)
+{
+  // copy constructor
+}
 
-  for(Int_t ij=0; ij<nRecJets; ij++)
-  {  
-    fptJet->Fill(recJets[ij].Pt());
-    fetaJet->Fill(recJets[ij].Eta());
-    fphiJet->Fill(recJets[ij].Phi());
+//_______________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsJetPt       = o.fNBinsJetPt;
+    fJetPtMin         = o.fJetPtMin;
+    fJetPtMax         = o.fJetPtMax;
+    fNBinsPt          = o.fNBinsPt; 
+    fPtMin            = o.fPtMin;   
+    fPtMax            = o.fPtMax;   
+    fNBinsZ           = o.fNBinsZ; 
+    fZMin             = o.fZMin;   
+    fZMax             = o.fZMax;   
+    fNBinsJt          = o.fNBinsJt;
+    fJtMin            = o.fJtMin;
+    fJtMax            = o.fJtMax;
+    fNBinsTheta       = o.fNBinsTheta;
+    fThetaMin         = o.fThetaMin;
+    fThetaMax         = o.fThetaMax;
+    fNBinsCosTheta    = o.fNBinsCosTheta;
+    fCosThetaMin      = o.fCosThetaMin;
+    fCosThetaMax      = o.fCosThetaMax;
+    fh2Theta          = o.fh2Theta;
+    fh2CosTheta       = o.fh2CosTheta;
+    fh2Jt             = o.fh2Jt;
+    fh2PtvsZ          = o.fh2PtvsZ;
+    fhnIntraJet       = o.fhnIntraJet;
+    fnDim             = o.fnDim;
+    fName             = o.fName;
   }
+    
+  return *this;
+}
 
-  // Get track ref 
-  TRefArray* ref = recJets[0].GetRefTracks();
-  for(Int_t it=0; it<ref->GetEntries(); it++)
-  {
-    Float_t ptTrack = ((AliVTrack*)ref->At(it))->Pt();
-    fPtvsEtaJet->Fill(recJets[0].Eta(),ptTrack);
-    fNpvsEtaJet->Fill(recJets[0].Eta(),ref->GetEntries());
-    fNpevtvsEtaJet->Fill(recJets[0].Eta(),evtMult);
-    fPtvsPtJet->Fill(recJets[0].Pt(),ptTrack);
-    fNpvsPtJet->Fill(recJets[0].Pt(),ref->GetEntries());
-    fNpevtvsPtJet->Fill(recJets[0].Pt(),evtMult);
-    fPtvsPtJet1D->Fill(recJets[0].Pt(),ptTrack);
-    fNpvsPtJet1D->Fill(recJets[0].Pt(),ref->GetEntries());
-    fNpevtvsPtJet1D->Fill(recJets[0].Pt(),evtMult);
-  }
+//_________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::~AliFragFuncIntraJetHistos()
+{
+  // destructor 
 
-  FillMonoJetH(goodBin, recJets, tracks);
 
-  //////////////////////////////////////////////////
-  ////////////////// DIJET PART ////////////////////       
-  //////////////////////////////////////////////////
+  if(fh2Theta)          delete fh2Theta;
+  if(fh2CosTheta)       delete fh2CosTheta;
+  if(fh2Jt)             delete fh2Jt;
+  if(fh2PtvsZ)          delete fh2PtvsZ;
+  if(fhnIntraJet)       delete fhnIntraJet;
 
-  // UNDER CONSTRUCTION
+}
 
-  PostData(1, fHistList);
+//_________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::DefineHistos()
+{
+  // book FF histos
+
+  fh2Theta    = new TH2F(Form("fh2IJTheta%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsTheta, fThetaMin, fThetaMax);
+  fh2CosTheta = new TH2F(Form("fh2IJcosTheta%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsCosTheta, fCosThetaMin, fCosThetaMax);
+  fh2Jt       = new TH2F(Form("fh2IJJt%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsJt, fJtMin, fJtMax);
+
+  // Create 3D histograms
+  Int_t    *iBin = new Int_t[fnDim];
+  Double_t *min  = new Double_t[fnDim];
+  Double_t *max  = new Double_t[fnDim];
+
+  iBin[0] = fNBinsJetPt; iBin[1] = fNBinsTheta; iBin[2] = fNBinsCosTheta; iBin[3] = fNBinsJt; iBin[4] = fNBinsZ; iBin[5] = fNBinsPt;
+  min[0]  = fJetPtMin; min[1] = fThetaMin; min[2] = fCosThetaMin; min[3] = fJtMin; min[4] = fZMin; min[5] = fPtMin; 
+  max[0]  = fJetPtMax; max[1] = fThetaMax; max[2] = fCosThetaMax; max[3] = fJtMax; max[4] = fZMax; max[5] = fPtMax;
+
+  const char* title = Form("fhnIntraJetPart%s",fName.Data());
+  const char* comment = "THnSparseF p_{T} jet [GeV/c] : #Theta : cos(#Theta) : j_{T} : Z : p_{T} part [GeV/c]";
+  fhnIntraJet = new THnSparseF(title,comment,fnDim,iBin,min,max);
+
+  const char** axisTitle = new const char*[fnDim];
+  axisTitle[0] = "p_{T}^{jet} [GeV/c]";
+  axisTitle[1] = "#Theta";
+  axisTitle[2] = "Cos(#Theta)";
+  axisTitle[3] = "j_{T} [GeV]";
+  axisTitle[4] = "z = p_{T}^{had}/p_{T}^{jet}";
+  axisTitle[5] = "p_{T}^{had} [GeV/c]";
+  
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Theta,"jet p_{T} [GeV/c]","#Theta","entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2CosTheta,"jet p_{T} [GeV/c]","cos(#Theta)", "entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Jt,"jet p_{T} [GeV/c]","j_{T}","entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fhnIntraJet,fnDim,axisTitle);
+  delete[] iBin;
+  delete[] min;
+  delete[] max;
+  delete[] axisTitle;
 }
 
-//#######################################################################
-void AliAnalysisTaskFragmentationFunction::Terminate(Option_t */*option*/)
+//_______________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::FillIntraJet(TLorentzVector* trackV, TLorentzVector* jetV)
 {
-// Terminate analysis
-//
-    if (fDebug > 1) printf("AnalysisDiJets: Terminate() \n");
-    printf("Number of events with vertex out of bound: %d", fNBadRuns);
+  // fill IntraJet histos
+  Float_t cosTheta = 0.; Float_t theta = 0.; 
+  Float_t jt = 0.; Float_t z = 0.; 
+  // For Theta distribution
+  Float_t pxT  = trackV->Px();
+  Float_t pyT  = trackV->Py();
+  Float_t pzT  = trackV->Pz();
+  Float_t ptT  = trackV->Pt();
+  Float_t pT   = trackV->P();
+  Float_t etaT = trackV->Eta();
+  Float_t phiT = trackV->Phi(); // Check the value returned
+  Float_t pxJ = jetV->Px();
+  Float_t pyJ = jetV->Py();
+  Float_t pzJ = jetV->Pz();
+  Float_t ptJ = jetV->Pt();
+  Float_t pJ  = jetV->P();
+
+  // Compute z
+  z = (Float_t)(ptT/ptJ);
+
+  // Compute theta
+  cosTheta = (pxT*pxJ+pyT*pyJ+pzT*pzJ)/(pT*pJ);
+  theta = TMath::ACos(cosTheta);
+
+  // Compute jt
+  TVector3 trackP; TVector3 jetP;
+  jetP[0] = pxJ;
+  jetP[1] = pyJ;
+  jetP[2] = pzJ;
+  trackP.SetPtEtaPhi(ptT,etaT,phiT);
+  jt = TMath::Sin(trackP.Angle(jetP))*trackP.Mag();
+
+  // Fill histos and THnSparse
+  fh2CosTheta->Fill(ptJ,cosTheta);
+  fh2Theta->Fill(ptJ,theta);
+  fh2Jt->Fill(ptJ,jt);
+
+  // Fill THnSparse
+  Double_t *content = new Double_t[fnDim];
+  content[0]= ptJ; content[1] = theta; content[2] = cosTheta; content[3] = jt; content[4] = z; content[5] = ptT; 
+
+  fhnIntraJet->Fill(content);
+
+  delete content;
+
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AliFragFuncDiJetHistos(const char* name, Int_t kindSlices,
+                                                        Int_t nJetInvMass, Float_t jetInvMassMin, Float_t jetInvMassMax,  
+                                                        Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
+                                                        Int_t nPt, Float_t ptMin, Float_t ptMax,
+                                                        Int_t nXi, Float_t xiMin, Float_t xiMax,
+                                                        Int_t nZ , Float_t zMin , Float_t zMax)
+  : TObject()
+  ,fKindSlices(kindSlices)
+  ,fNBinsJetInvMass(nJetInvMass)
+  ,fJetInvMassMin(jetInvMassMin)
+  ,fJetInvMassMax(jetInvMassMax)
+  ,fNBinsJetPt(nJetPt)
+  ,fJetPtMin(jetPtMin)
+  ,fJetPtMax(jetPtMax)
+  ,fNBinsPt(nPt) 
+  ,fPtMin(ptMin)   
+  ,fPtMax(ptMax)   
+  ,fNBinsXi(nXi) 
+  ,fXiMin(xiMin)   
+  ,fXiMax(xiMax)   
+  ,fNBinsZ(nZ)  
+  ,fZMin(zMin)    
+  ,fZMax(zMax)
+  ,fh2TrackPtJet1(0)
+  ,fh2TrackPtJet2(0)
+  ,fh2TrackPtJet(0)
+  ,fh1Jet1Pt(0)
+  ,fh1Jet2Pt(0)
+  ,fh1JetPt(0)
+  ,fh2Xi1(0)
+  ,fh2Xi2(0)
+  ,fh2Xi(0)
+  ,fh2Z1(0)
+  ,fh2Z2(0)
+  ,fh2Z(0)
+  ,fh2Pt1(0)
+  ,fh2Pt2(0)
+  ,fh2Pt(0)
+  ,fName(name)
+{
+  // default constructor
+
+}
 
-void AliAnalysisTaskFragmentationFunction::DefineJetH()
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AliFragFuncDiJetHistos(const AliFragFuncDiJetHistos& copy)
+  : TObject()
+  ,fKindSlices(copy.fKindSlices)
+  ,fNBinsJetInvMass(copy.fNBinsJetInvMass)
+  ,fJetInvMassMin(copy.fJetInvMassMin)
+  ,fJetInvMassMax(copy.fJetInvMassMax)
+  ,fNBinsJetPt(copy.fNBinsJetPt)
+  ,fJetPtMin(copy.fJetPtMin)
+  ,fJetPtMax(copy.fJetPtMax)
+  ,fNBinsPt(copy.fNBinsPt) 
+  ,fPtMin(copy.fPtMin)   
+  ,fPtMax(copy.fPtMax)   
+  ,fNBinsXi(copy.fNBinsXi) 
+  ,fXiMin(copy.fXiMin)   
+  ,fXiMax(copy.fXiMax)   
+  ,fNBinsZ(copy.fNBinsZ)  
+  ,fZMin(copy.fZMin)    
+  ,fZMax(copy.fZMax)
+  ,fh2TrackPtJet1(copy.fh2TrackPtJet1)
+  ,fh2TrackPtJet2(copy.fh2TrackPtJet2)
+  ,fh2TrackPtJet(copy.fh2TrackPtJet)
+  ,fh1Jet1Pt(copy.fh1Jet1Pt)
+  ,fh1Jet2Pt(copy.fh1Jet2Pt)
+  ,fh1JetPt(copy.fh1JetPt)
+  ,fh2Xi1(copy.fh2Xi1)
+  ,fh2Xi2(copy.fh2Xi2)
+  ,fh2Xi(copy.fh2Xi2)
+  ,fh2Z1(copy.fh2Z1)
+  ,fh2Z2(copy.fh2Z2)
+  ,fh2Z(copy.fh2Z)
+  ,fh2Pt1(copy.fh2Pt1)
+  ,fh2Pt2(copy.fh2Pt2)
+  ,fh2Pt(copy.fh2Pt)
+  ,fName(copy.fName)
 {
+  // default constructor
 
-  /////////////////////////////////////// HISTOGRAMS FIRST JET
-  fEtaMonoJet1H   = new TH1F*[fnEBin+1];
-  fPhiMonoJet1H   = new TH1F*[fnEBin+1];
-  fPtMonoJet1H    = new TH1F*[fnEBin+1];
-  fEMonoJet1H     = new TH1F*[fnEBin+1];
-
-  fdNdXiMonoJet1H = new TH1F**[fnEBin+1];
-  fdNdPtMonoJet1H = new TH1F**[fnEBin+1];
-  fdNdZMonoJet1H  = new TH1F**[fnEBin+1];
-  fdNdThetaMonoJet1H    = new TH1F**[fnEBin+1];
-  fdNdcosThetaMonoJet1H    = new TH1F**[fnEBin+1];
-  fdNdkTMonoJet1H       = new TH1F**[fnEBin+1];
-  fdNdpTvsZMonoJet1H    = new TH1F**[fnEBin+1];
-  fShapeMonoJet1H       = new TH1F**[fnEBin+1];
-  fNMonoJet1sH          = new TH1F**[fnEBin+1];
-
-  fThetaPtPartMonoJet1H = new TH2F**[fnEBin+1];
-  fcosThetaPtPartMonoJet1H = new TH2F**[fnEBin+1];
-  fkTPtPartMonoJet1H    = new TH2F**[fnEBin+1];
-  fThetaPtJetMonoJet1H = new TH2F**[fnEBin+1];
-  fcosThetaPtJetMonoJet1H = new TH2F**[fnEBin+1];
-  fkTPtJetMonoJet1H    = new TH2F**[fnEBin+1];
-  fpTPtJetMonoJet1H    = new TH2F**[fnEBin+1];
-
-  for (Int_t iEbin=0;iEbin<fnEBin+1;iEbin++)
-  {
-    fdNdXiMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
-    fdNdPtMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
-    fdNdZMonoJet1H[iEbin]        = new TH1F*[fnRBin+1];
-    fdNdThetaMonoJet1H[iEbin]    = new TH1F*[fnRBin+1];
-    fdNdcosThetaMonoJet1H[iEbin] = new TH1F*[fnRBin+1];
-    fdNdkTMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
-    fdNdpTvsZMonoJet1H[iEbin]    = new TH1F*[fnRBin+1];
-    fShapeMonoJet1H[iEbin]       = new TH1F*[fnRBin+1];
-    fNMonoJet1sH[iEbin]          = new TH1F*[fnRBin+1];
-
-    fThetaPtPartMonoJet1H[iEbin]    = new TH2F*[fnRBin+1];
-    fcosThetaPtPartMonoJet1H[iEbin] = new TH2F*[fnRBin+1];
-    fkTPtPartMonoJet1H[iEbin]       = new TH2F*[fnRBin+1];
-    fThetaPtJetMonoJet1H[iEbin]     = new TH2F*[fnRBin+1];
-    fcosThetaPtJetMonoJet1H[iEbin]  = new TH2F*[fnRBin+1];
-    fkTPtJetMonoJet1H[iEbin]        = new TH2F*[fnRBin+1];
-    fpTPtJetMonoJet1H[iEbin]        = new TH2F*[fnRBin+1];
+}
+
+//_______________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fKindSlices      = o.fKindSlices;
+    fNBinsJetInvMass = o.fNBinsJetInvMass;
+    fJetInvMassMin   = o.fJetInvMassMin;
+    fJetInvMassMax   = o.fJetInvMassMax;
+    fNBinsJetPt      = o.fNBinsJetPt;
+    fJetPtMin        = o.fJetPtMin;
+    fJetPtMax        = o.fJetPtMax;
+    fNBinsPt         = o.fNBinsPt; 
+    fPtMin           = o.fPtMin;   
+    fPtMax           = o.fPtMax;   
+    fNBinsXi         = o.fNBinsXi; 
+    fXiMin           = o.fXiMin;   
+    fXiMax           = o.fXiMax;   
+    fNBinsZ          = o.fNBinsZ;  
+    fZMin            = o.fZMin;    
+    fZMax            = o.fZMax;   
+    fh2TrackPtJet1   = o.fh2TrackPtJet1;
+    fh2TrackPtJet2   = o.fh2TrackPtJet2;
+    fh2TrackPtJet    = o.fh2TrackPtJet;
+    fh1Jet1Pt        = o.fh1Jet1Pt;
+    fh1Jet2Pt        = o.fh1Jet2Pt;
+    fh1JetPt         = o.fh1JetPt;
+    fh2Xi1           = o.fh2Xi1;
+    fh2Xi2           = o.fh2Xi2;
+    fh2Xi            = o.fh2Xi;
+    fh2Z1            = o.fh2Z1;
+    fh2Z2            = o.fh2Z2;
+    fh2Z             = o.fh2Z;
+    fh2Pt1           = o.fh2Pt1;
+    fh2Pt2           = o.fh2Pt2;
+    fh2Pt            = o.fh2Pt;
+    fName            = o.fName;
   }
+    
+  return *this;
+}
 
-  for (Int_t iEbin=0;iEbin<fnEBin+1;iEbin++)
-  {
-    fEtaMonoJet1H[iEbin]    = 0;
-    fPhiMonoJet1H[iEbin]    = 0;
-    fPtMonoJet1H[iEbin]     = 0;
-    fEMonoJet1H[iEbin]      = 0;
+//_________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::~AliFragFuncDiJetHistos()
+{
+  // destructor 
+
+  if(fh2TrackPtJet1) delete fh2TrackPtJet1;
+  if(fh2TrackPtJet2) delete fh2TrackPtJet2;
+  if(fh2TrackPtJet ) delete fh2TrackPtJet;
+  if(fh1Jet1Pt)      delete fh1Jet1Pt;
+  if(fh1Jet2Pt)      delete fh1Jet2Pt;
+  if(fh1JetPt)       delete fh1JetPt;
+  if(fh2Xi1)         delete fh2Xi1;
+  if(fh2Xi2)         delete fh2Xi2;
+  if(fh2Xi)          delete fh2Xi;
+  if(fh2Z1)          delete fh2Z1;
+  if(fh2Z2)          delete fh2Z2;
+  if(fh2Z)           delete fh2Z;
+  if(fh2Pt1)         delete fh2Pt1;
+  if(fh2Pt2)         delete fh2Pt2;
+  if(fh2Pt)          delete fh2Pt;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::DefineDiJetHistos()
+{
 
-    for (Int_t iRbin=0;iRbin<fnRBin+1;iRbin++)
+  Int_t nBins = 0;
+  Double_t min = 0.;
+  Double_t max = 0.;
+  const char *xaxis = "";
+  if(fKindSlices == 1)
     {
-      fdNdXiMonoJet1H[iEbin][iRbin]       = 0;
-      fdNdPtMonoJet1H[iEbin][iRbin]       = 0;
-      fdNdZMonoJet1H[iEbin][iRbin]        = 0;
-      fdNdThetaMonoJet1H[iEbin][iRbin]    = 0;
-      fdNdcosThetaMonoJet1H[iEbin][iRbin] = 0;
-      fdNdkTMonoJet1H[iEbin][iRbin]       = 0;
-      fdNdpTvsZMonoJet1H[iEbin][iRbin]    = 0;
-      fShapeMonoJet1H[iEbin][iRbin]       = 0;
-      fNMonoJet1sH[iEbin][iRbin]          = 0;
-
-      fThetaPtPartMonoJet1H[iEbin][iRbin]    = 0;
-      fcosThetaPtPartMonoJet1H[iEbin][iRbin] = 0;
-      fkTPtPartMonoJet1H[iEbin][iRbin]       = 0;
-      fThetaPtJetMonoJet1H[iEbin][iRbin]     = 0;
-      fcosThetaPtJetMonoJet1H[iEbin][iRbin]  = 0;
-      fkTPtJetMonoJet1H[iEbin][iRbin]        = 0;
-      fpTPtJetMonoJet1H[iEbin][iRbin]        = 0;
+      nBins = fNBinsJetInvMass;
+      min   = fJetInvMassMin;
+      max   = fJetInvMassMax;
+      xaxis = "M_{JJ} [GeV]";
     }
-  }
+  if(fKindSlices == 2 || fKindSlices == 3)
+    {
+      nBins = fNBinsJetPt;
+      min   = fJetPtMin;
+      max   = fJetPtMax;
+      if(fKindSlices == 2) xaxis = "E_{Tmean} [GeV]";
+      if(fKindSlices == 3) xaxis ="leading jet p_{T} [GeV/c]";
+    }
+  
+  fh1Jet1Pt      = new TH1F(Form("fh1DJJet1Pt%s", fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
+  fh1Jet2Pt      = new TH1F(Form("fh1DJJet2Pt%s", fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
+  fh1JetPt       = new TH1F(Form("fh1DJJetPt%s",  fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
+  
+  fh2TrackPtJet1 = new TH2F(Form("fh2DJTrackPtJet1%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+  fh2TrackPtJet2 = new TH2F(Form("fh2DJTrackPtJet2%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+  fh2TrackPtJet  = new TH2F(Form("fh2DJTrackPtJet%s", fName.Data()),  "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+  
+  fh2Xi1         = new TH2F(Form("fh2DJXi1%s", fName.Data()), "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
+  fh2Xi2         = new TH2F(Form("fh2DJXi2%s", fName.Data()), "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
+  fh2Xi          = new TH2F(Form("fh2DJXi%s", fName.Data()),  "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
+  
+  fh2Z1          = new TH2F(Form("fh2DJZ1%s", fName.Data()), "",nBins, min, max, fNBinsZ, fZMin, fZMax);
+  fh2Z2          = new TH2F(Form("fh2DJZ2%s", fName.Data()), "",nBins, min, max, fNBinsZ, fZMin, fZMax);
+  fh2Z           = new TH2F(Form("fh2DJZ%s", fName.Data()),  "",nBins, min, max, fNBinsZ, fZMin, fZMax);
+  
+  fh2Pt1         = new TH2F(Form("fh2DJPt1%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+  fh2Pt2         = new TH2F(Form("fh2DJPt2%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+  fh2Pt          = new TH2F(Form("fh2DJPtZ%s", fName.Data()),  "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
+      
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1Jet1Pt, "p_{T} [GeV/c]", "entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1Jet2Pt, "p_{T} [GeV/c]", "entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries");
+
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet1, xaxis, "p_{T} [GeV/c]", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet2, xaxis, "p_{T} [GeV/c]", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet, xaxis, "p_{T} [GeV/c]", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi1, xaxis, "#xi", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi2, xaxis, "#xi", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi, xaxis, "#xi", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z1, xaxis, "z", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z2, xaxis, "z", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z, xaxis, "z", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt1, xaxis, "p_{T} [GeV/c]", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt2, xaxis, "p_{T} [GeV/c]", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt, xaxis, "p_{T} [GeV/c]", "Entries");
+  
+}
 
-  fptTracks          = new TH1F*[fnPtTrigBin+1];
-  fetaTracks         = new TH1F*[fnPtTrigBin+1];
-  fphiTracks         = new TH1F*[fnPtTrigBin+1];
-  fdetaTracks        = new TH1F*[fnPtTrigBin+1];
-  fdphiTracks        = new TH1F*[fnPtTrigBin+1];
-  fetaphiptTracks    = new TH2F*[fnPtTrigBin+1];
-  fetaphiTracks      = new TH2F*[fnPtTrigBin+1];
-  fdetadphiTracks    = new TH2F*[fnPtTrigBin+1];
-  fptTracksCut       = new TH1F*[fnPtTrigBin+1];
-  fetaTracksCut      = new TH1F*[fnPtTrigBin+1];
-  fphiTracksCut      = new TH1F*[fnPtTrigBin+1];
-  fdetaTracksCut     = new TH1F*[fnPtTrigBin+1];
-  fdphiTracksCut     = new TH1F*[fnPtTrigBin+1];
-  fetaphiptTracksCut = new TH2F*[fnPtTrigBin+1];
-  fetaphiTracksCut   = new TH2F*[fnPtTrigBin+1];
-  fdetadphiTracksCut = new TH2F*[fnPtTrigBin+1];
-  fNPtTrig           = new TH1F*[fnPtTrigBin+1];
-  fNPtTrigCut        = new TH1F*[fnPtTrigBin+1];
-
-  for(Int_t iPtTrigBin=0; iPtTrigBin<fnPtTrigBin; iPtTrigBin++)
-  {
-    fptTracks[iPtTrigBin]          = 0;
-    fetaTracks[iPtTrigBin]         = 0;
-    fphiTracks[iPtTrigBin]         = 0;
-    fdetaTracks[iPtTrigBin]        = 0;
-    fdphiTracks[iPtTrigBin]        = 0;
-    fetaphiptTracks[iPtTrigBin]    = 0;
-    fetaphiTracks[iPtTrigBin]      = 0;
-    fdetadphiTracks[iPtTrigBin]    = 0;
-    fptTracksCut[iPtTrigBin]       = 0;
-    fetaTracksCut[iPtTrigBin]      = 0;
-    fphiTracksCut[iPtTrigBin]      = 0;
-    fdetaTracksCut[iPtTrigBin]     = 0;
-    fdphiTracksCut[iPtTrigBin]     = 0;
-    fetaphiptTracksCut[iPtTrigBin] = 0;
-    fetaphiTracksCut[iPtTrigBin]   = 0;
-    fdetadphiTracksCut[iPtTrigBin] = 0;
-    fNPtTrig[iPtTrigBin]           = 0;
-    fNPtTrigCut[iPtTrigBin]        = 0;
-  }
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::FillDiJetFF(Int_t jetType, Float_t trackPt, Float_t jetPt, Double_t jetBin, Bool_t incrementJetPt)
+{
+  if(jetType == 0)
+    {
+      if(incrementJetPt) fh1JetPt->Fill(jetPt);  
+      
+      fh2TrackPtJet->Fill(jetBin, trackPt);
+      
+      Double_t z = trackPt / jetPt;
+      Double_t xi = 0;
+      if(z!=0) xi = TMath::Log(1/z);
+      
+      fh2Xi->Fill(jetBin, xi);
+      fh2Z->Fill(jetBin, z);
+    }
+  if(jetType == 1)
+    {
+      if(incrementJetPt) fh1Jet1Pt->Fill(jetPt);
+      
+      fh2TrackPtJet1->Fill(jetBin, trackPt);
+      
+      Double_t z = trackPt / jetPt;
+      Double_t xi = 0;
+      if(z!=0) xi = TMath::Log(1/z);
+      
+      fh2Xi1->Fill(jetBin, xi);
+      fh2Z1->Fill(jetBin, z);
+    }
+  if(jetType == 2)
+    {
+      if(incrementJetPt) fh1Jet2Pt->Fill(jetPt);
+      
+      fh2TrackPtJet2->Fill(jetBin, trackPt);
+      
+      Double_t z = trackPt / jetPt;
+      Double_t xi = 0;
+      if(z!=0) xi = TMath::Log(1/z);
+      
+      fh2Xi2->Fill(jetBin, xi);
+      fh2Z2->Fill(jetBin, z);
+    }
 
-  farrayEmin = new Double_t[fnEBin];
-  farrayEmax = new Double_t[fnEBin];
 
-  farrayPtTrigmin = new Double_t[fnPtTrigBin];
-  farrayPtTrigmax = new Double_t[fnPtTrigBin];
+}
 
-  farrayRadii = new Double_t[fnRBin];
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AddToOutput(TList* list)const
+{
+  list->Add(fh1Jet1Pt);
+  list->Add(fh1Jet2Pt);
+  list->Add(fh1JetPt);
+  list->Add(fh2TrackPtJet1);
+  list->Add(fh2TrackPtJet2);
+  list->Add(fh2TrackPtJet);
+  list->Add(fh2Xi1);
+  list->Add(fh2Xi2);
+  list->Add(fh2Xi);
+  list->Add(fh2Z1);
+  list->Add(fh2Z2);
+  list->Add(fh2Z);
+}
 
-  Double_t Emin    = 0;
-  Double_t Emax    = 0;
-  Double_t pasE    = (Double_t)((fEmax-fEmin)/fnEInterval);
-  TString energy;
-  TString energy2;
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AliFragFuncQADiJetHistos(const char* name, Int_t kindSlices,
+                                           Int_t nInvMass, Float_t invMassMin, Float_t invMassMax, 
+                                            Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
+                                           Int_t nDeltaPhi, Float_t deltaPhiMin, Float_t deltaPhiMax, 
+                                          Int_t nDeltaEta, Float_t deltaEtaMin, Float_t deltaEtaMax, 
+                                          Int_t nDeltaPt, Float_t deltaPtMin, Float_t deltaPtMax)
+  : TObject()
+  ,fKindSlices(kindSlices)
+  ,fNBinsJetInvMass(nInvMass)
+  ,fJetInvMassMin(invMassMin)
+  ,fJetInvMassMax(invMassMax)
+  ,fNBinsJetPt(nJetPt)
+  ,fJetPtMin(jetPtMin)
+  ,fJetPtMax(jetPtMax)
+  ,fNBinsDeltaPhi(nDeltaPhi)
+  ,fDeltaPhiMin(deltaPhiMin)
+  ,fDeltaPhiMax(deltaPhiMax)
+  ,fNBinsDeltaEta(nDeltaEta)
+  ,fDeltaEtaMin(deltaEtaMin)
+  ,fDeltaEtaMax(deltaEtaMax)
+  ,fNBinsDeltaPt(nDeltaPt)
+  ,fDeltaPtMin(deltaPtMin)
+  ,fDeltaPtMax(deltaPtMax)
+  ,fh2InvMass(0)
+  ,fh2DeltaPhi(0)
+  ,fh2DeltaEta(0)
+  ,fh2DeltaPt(0)
+  ,fName(name)
+{
+  // default constructor
 
-  Double_t R    = 0.;
-  Double_t pasR = (Double_t)((fRmax-fRmin)/fnRInterval);
-  TString radius;
+}
 
-  for (Int_t i = 0; i < fnEBin; i++)
-  {
-    Emin       = 0;
-    Emax       = 0;
-
-    if (i==0) Emin = fEmin;
-    if (i!=0) Emin = fEmin + pasE*i;
-    Emax = Emin+pasE;
-    energy2 = "E_{jet1} : ";
-    energy = "E_{jet} : ";
-    energy += Emin;
-    energy2 += Emin;
-    energy += "-";
-    energy2 += "-";
-    energy +=Emax;
-    energy2 += Emax;
-    energy += "GeV";
-    energy2 += "GeV";
-
-    farrayEmin[i]      = Emin;
-    farrayEmax[i]      = Emax;
-
-    for (Int_t j = 0; j < fnRBin; j++)
-    {
-      if (j==0) R = fRmin;
-      if (j!=0) R = fRmin + pasR*j;
-      radius = ", R = ";
-      radius += R;    
-
-      farrayRadii[j] = R;
-
-      fEtaMonoJet1H[i]   = new TH1F("fEtaMonoJet1H,"+energy, "#eta_{jet1},"+energy, fnEtaHBin, fEtaHBinMin, fEtaHBinMax);
-      fPhiMonoJet1H[i]   = new TH1F("fPhiMonoJet1H,"+energy, "#phi_{jet1},"+energy, fnPhiHBin, fPhiHBinMin, fPhiHBinMax);
-      fPtMonoJet1H[i]    = new TH1F("fPtMonoJet1H,"+energy, "pT_{jet1},"+energy, fnPtHBin, fPtHBinMin, fPtHBinMax);
-      fEMonoJet1H[i]     = new TH1F("fEMonoJet1H,"+energy, "E_{jet1},"+energy, fnEHBin, fEHBinMin, fEHBinMax);
-
-      fdNdXiMonoJet1H[i][j] = new TH1F("fdNdXiMonoJet1H,"+energy+radius, "dN_{ch}/d#xi,"+energy+radius, fnXiHBin, fXiHBinMin, fXiHBinMax);
-      fdNdPtMonoJet1H[i][j] = new TH1F("fdNdPtMonoJet1H,"+energy+radius, "dN_{ch}/dPt_{had},"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax);
-      fdNdZMonoJet1H[i][j]  = new TH1F("fdNdZMonoJet1H,"+energy+radius, "dN_{ch}/dz,"+energy+radius, fnZHBin, fZHBinMin, fZHBinMax);
-
-      fdNdThetaMonoJet1H[i][j]    = new TH1F("fdNdThetaMonoJet1H,"+energy+radius, "dN_{ch}/d#Theta,"+energy+radius, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
-      fdNdcosThetaMonoJet1H[i][j] = new TH1F("fdNdcosThetaMonoJet1H,"+energy+radius, "dN_{ch}/dcos(#Theta),"+energy+radius, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
-      fdNdkTMonoJet1H[i][j]   = new TH1F("fdNdkTMonoJet1H,"+energy+radius, "dN_{ch}/dk_{T},"+energy+radius, fnkTHBin, fkTHBinMin, fkTHBinMax);
-      fdNdpTvsZMonoJet1H[i][j]= new TH1F("fdNdpTvsZMonoJet1H,"+energy+radius, "dN_{ch}/dk_{T} vs Z,"+energy+radius, fnZHBin, fZHBinMin, fZHBinMax);
-      fShapeMonoJet1H[i][j]   = new TH1F("fShapeMonoJet1H,"+energy+radius, "E(R=x)/E(R=1)"+energy+radius, fnRHBin, fRHBinMin, fRHBinMax);
-
-      fThetaPtPartMonoJet1H[i][j] = new TH2F("fThetaPtPartMonoJet1H,"+energy+radius, "#Theta vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
-      fcosThetaPtPartMonoJet1H[i][j] = new TH2F("fcosThetaPtPartMonoJet1H,"+energy+radius, "cos(#Theta) vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
-      fkTPtPartMonoJet1H[i][j] = new TH2F("fkTPtPartMonoJet1H,"+energy+radius, "kT vs Pt particle,"+energy+radius, fnPthadHBin, fPthadHBinMin, fPthadHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
-      fThetaPtJetMonoJet1H[i][j] = new TH2F("fThetaPtJetMonoJet1H,"+energy+radius, "#Theta vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnThetaHBin, fThetaHBinMin, fThetaHBinMax);
-      fcosThetaPtJetMonoJet1H[i][j] = new TH2F("fcosThetaPtJetMonoJet1H,"+energy+radius, "cos(#Theta) vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnCosThetaHBin, fcosThetaHBinMin, fcosThetaHBinMax);
-      fkTPtJetMonoJet1H[i][j] = new TH2F("fkTPtJetMonoJet1H,"+energy+radius, "kT vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
-      fpTPtJetMonoJet1H[i][j] = new TH2F("fpTPtJetMonoJet1H,"+energy+radius, "pT vs Pt jet,"+energy+radius, fnPtHBin, fPtHBinMin, fPtHBinMax, fnkTHBin, fkTHBinMin, fkTHBinMax);
-
-      fNMonoJet1sH[i][j]    = new TH1F("fNMonoJet1sH,"+energy+radius, "N_{jets1},"+energy+radius, 1, 0., 1.);
-      fNBadRunsH = new TH1F("fNBadRunsH","Number of events with Z vertex out of range", 1, 0., 1.);
-
-      SetProperties(fEtaMonoJet1H[i], "#eta_{jet1}", "Entries");
-      SetProperties(fPhiMonoJet1H[i], "#phi_{jet1}", "Entries");
-      SetProperties(fPtMonoJet1H[i], "p_{Tjet1} (GeV/c)", "Entries");
-      SetProperties(fEMonoJet1H[i], "E_{jet1} (GeV)", "Entries");
-
-      SetProperties(fdNdXiMonoJet1H[i][j], "#xi = ln(E_{jet1}/p_{Thad})", "dN_{had}/d#xi");
-      SetProperties(fdNdPtMonoJet1H[i][j], "p_{Thad} (GeV/c)", "dN_{had}/dp_{Thad}");
-      SetProperties(fdNdZMonoJet1H[i][j], "z = (p_{Thad}/E_{jet1})", "dN_{had}/dz");
-      SetProperties(fdNdThetaMonoJet1H[i][j], "#Theta", "dN_{had}/d#Theta");
-      SetProperties(fdNdcosThetaMonoJet1H[i][j], "cos(#Theta)", "dN_{had}/dcos(#Theta)");
-      SetProperties(fdNdkTMonoJet1H[i][j], "k_{Thad}", "dN_{had}/dk_{Thad}");
-      SetProperties(fdNdpTvsZMonoJet1H[i][j], "z = (p_{Thad}/E_{jet1})", "dN_{had}/dp_{T}");
-      SetProperties(fShapeMonoJet1H[i][j], "R", "#Psi(R)");
-      SetProperties(fNMonoJet1sH[i][j], "Bin", "N_{jets1}");
-
-      SetProperties(fThetaPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","#Theta");
-      SetProperties(fcosThetaPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","cos(#Theta)");
-      SetProperties(fkTPtPartMonoJet1H[i][j],"p_{Thad} [GeV/c]","k_{Thad}");
-      SetProperties(fThetaPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "#Theta");
-      SetProperties(fcosThetaPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "cos(#Theta)");
-      SetProperties(fkTPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "k_{Thad} [GeV/c]");
-      SetProperties(fpTPtJetMonoJet1H[i][j], "p_{Tjet} [GeV/c]", "p_{Thad} [GeV/c]");
-    }
-  }
+//______________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AliFragFuncQADiJetHistos(const AliFragFuncQADiJetHistos& copy)
+  : TObject()
+  ,fKindSlices(copy.fKindSlices)
+  ,fNBinsJetInvMass(copy.fNBinsJetInvMass)
+  ,fJetInvMassMin(copy.fJetInvMassMin)
+  ,fJetInvMassMax(copy.fJetInvMassMax)
+  ,fNBinsJetPt(copy.fNBinsJetPt)
+  ,fJetPtMin(copy.fJetPtMin)
+  ,fJetPtMax(copy.fJetPtMax)
+  ,fNBinsDeltaPhi(copy.fNBinsDeltaPhi)
+  ,fDeltaPhiMin(copy.fDeltaPhiMin)
+  ,fDeltaPhiMax(copy.fDeltaPhiMax)
+  ,fNBinsDeltaEta(copy.fNBinsDeltaEta)
+  ,fDeltaEtaMin(copy.fDeltaEtaMin)
+  ,fDeltaEtaMax(copy.fDeltaEtaMax)
+  ,fNBinsDeltaPt(copy.fNBinsDeltaPt)
+  ,fDeltaPtMin(copy.fDeltaPtMin)
+  ,fDeltaPtMax(copy.fDeltaPtMax)
+  ,fh2InvMass(copy.fh2InvMass)
+  ,fh2DeltaPhi(copy.fh2DeltaPhi)
+  ,fh2DeltaEta(copy.fh2DeltaEta)
+  ,fh2DeltaPt(copy.fh2DeltaPt)
+  ,fName(copy.fName)
+{
+  // default constructor
 
-  for(Int_t i=0; i<fnPtTrigBin; i++)
-  {
-    if(i==0) farrayPtTrigmin[i] = 1.; 
-    else farrayPtTrigmin[i] = i*5.;
-    farrayPtTrigmax[i] = i*5+5.;
-    
-    TString ptTrigRange;
-    ptTrigRange = "; p_{T} trig range: ";
-    ptTrigRange += farrayPtTrigmin[i];
-    ptTrigRange += "-";
-    ptTrigRange += farrayPtTrigmax[i];
-    ptTrigRange += " [GeV]";
-
-    fptTracks[i] = new TH1F("fptTracks"+ptTrigRange, "Track transverse momentum [GeV]"+ptTrigRange,300,0.,150.);
-    fetaTracks[i] = new TH1F("fetaTracks"+ptTrigRange, "#eta tracks"+ptTrigRange, 36, -0.9, 0.9);
-    fphiTracks[i] = new TH1F("fphiTracks"+ptTrigRange, "#phi tracks"+ptTrigRange, 60, 0., 2*TMath::Pi());
-    fdetaTracks[i] = new TH1F("fdetaTracks"+ptTrigRange, "#Delta #eta tracks"+ptTrigRange,80, -2., 2.);
-    fdphiTracks[i] = new TH1F("fdphiTracks"+ptTrigRange, "#Delta #phi tracks"+ptTrigRange, 120, -TMath::Pi(), TMath::Pi());
-    fetaphiptTracks[i] = new TH2F("fetaphiptTracks"+ptTrigRange,"#eta/#phi track p_{T} mapping"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-    fetaphiTracks[i] = new TH2F("fetaphiTracks"+ptTrigRange,"#eta/#phi track mapping"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-    fdetadphiTracks[i] = new TH2F("fdetadphiTracks"+ptTrigRange,"#Delta #eta/#Delta #phi track mapping"+ptTrigRange,80, -2., 2., 120, -TMath::Pi(), TMath::Pi());
-    fptTracksCut[i] = new TH1F("fptTracksCut"+ptTrigRange, "Track transverse momentum after cuts [GeV]"+ptTrigRange,300,0.,150.);
-    fetaTracksCut[i] = new TH1F("fetaTracksCut"+ptTrigRange, "#eta tracks after cuts"+ptTrigRange, 36, -0.9, 0.9);
-    fphiTracksCut[i] = new TH1F("fphiTracksCuts"+ptTrigRange, "#phi tracks after cuts"+ptTrigRange, 60, 0., 2*TMath::Pi());
-    fdetaTracksCut[i] = new TH1F("fdetaTracksCuts"+ptTrigRange, "#Delta #eta tracks after cuts"+ptTrigRange,80, -2., 2.);
-    fdphiTracksCut[i] = new TH1F("fdphiTracksCuts"+ptTrigRange, "#Delta #phi tracks after cuts"+ptTrigRange, 120, -TMath::Pi(), TMath::Pi());
-    fetaphiptTracksCut[i] = new TH2F("fetaphiptTracksCuts"+ptTrigRange,"#eta/#phi track p_{T} mapping after cuts"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-    fetaphiTracksCut[i] = new TH2F("fetaphiTracksCuts"+ptTrigRange,"#eta/#phi track mapping after cuts"+ptTrigRange,36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-    fdetadphiTracksCut[i] = new TH2F("fdetadphiTracksCuts"+ptTrigRange,"#Delta #eta/#Delta #phi track mapping after cuts"+ptTrigRange,80, -2., 2., 120, -TMath::Pi(), TMath::Pi());
-    fNPtTrig[i] = new TH1F("fNPtTrig"+ptTrigRange,"Number of triggers"+ptTrigRange,1,0.,1.);
-    fNPtTrigCut[i] = new TH1F("fNPtTrigCut"+ptTrigRange,"Number of triggers after cut"+ptTrigRange,1,0.,1.);
-
-    SetProperties(fptTracks[i], "Track p_{T} [GeV]", "dN/dp_{T}"); 
-    SetProperties(fetaTracks[i], "Track #eta", "dN/d#eta"); 
-    SetProperties(fphiTracks[i], "Track #phi", "dN/d#phi");
-    SetProperties(fdetaTracks[i], "#eta_{track} - #eta_{trig}", "dN/d#Delta#eta");
-    SetProperties(fdphiTracks[i], "#phi_{track} - #phi_{trig}", "dN/d#Delta#phi");
-    SetProperties(fetaphiptTracks[i], "#eta_{track}", "#phi_{track}"); 
-    SetProperties(fetaphiTracks[i], "#eta_{track}", "#phi_{track}");
-    SetProperties(fdetadphiTracks[i], "#Delta #eta_{track}", "#Delta #phi_{track}");
-    SetProperties(fptTracksCut[i], "p_{T}track [GeV]", "dN/dp_{T}");
-    SetProperties(fetaTracksCut[i], "#eta_{track}", "dN/d#eta"); 
-    SetProperties(fphiTracksCut[i], "#phi_{track}", "dN/d#phi"); 
-    SetProperties(fdetaTracksCut[i], "#eta_{track} - #eta_{trig}", "dN/d#Delta#eta");
-    SetProperties(fdphiTracksCut[i], "#phi_{track} - #phi_{trig}", "dN/d#Delta#phi");
-    SetProperties(fetaphiptTracksCut[i], "#eta_{track}", "#phi_{track}");
-    SetProperties(fetaphiTracksCut[i], "#eta_{track}", "#phi_{track}");
-    SetProperties(fdetadphiTracksCut[i], "#Delta #eta_{track}", "#Delta #phi_{track}");
-    SetProperties(fNPtTrig[i], "", "Number of triggers");
-    SetProperties(fNPtTrigCut[i], "", "Number of triggers");
+}
 
+//_______________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fKindSlices       = o.fKindSlices;
+    fNBinsJetInvMass  = o.fNBinsJetInvMass;
+    fJetInvMassMin    = o.fJetInvMassMin;
+    fJetInvMassMax    = o.fJetInvMassMax;
+    fNBinsJetPt       = o.fNBinsJetPt;
+    fJetPtMin         = o.fJetPtMin;
+    fJetPtMax         = o.fJetPtMax;
+    fNBinsDeltaPhi    = o.fNBinsDeltaPhi;
+    fDeltaPhiMin      = o.fDeltaPhiMin;
+    fDeltaPhiMax      = o.fDeltaPhiMax;
+    fNBinsDeltaEta    = o.fNBinsDeltaEta;
+    fDeltaEtaMin      = o.fDeltaEtaMin;
+    fDeltaEtaMax      = o.fDeltaEtaMax;
+    fNBinsDeltaPt     = o.fNBinsDeltaPt;
+    fDeltaPtMin       = o.fDeltaPtMin;
+    fDeltaPtMax       = o.fDeltaPtMax;
+    fh2InvMass        = o.fh2InvMass;
+    fh2DeltaPhi       = o.fh2DeltaPhi;
+    fh2DeltaEta       = o.fh2DeltaEta;
+    fh2DeltaPt        = o.fh2DeltaPt;
+    fName             = o.fName;
   }
-
-  fptAllTracks = new TH1F("fptAllTracks", "Track transverse momentum [GeV]",300,0.,150.);
-  fetaAllTracks = new TH1F("fetaAllTracks", "#eta tracks", 36, -0.9, 0.9);
-  fphiAllTracks = new TH1F("fphiAllTracks", "#phi tracks", 60, 0., 2*TMath::Pi());
-  fetaphiptAllTracks = new TH2F("fetaphiptAllTracks","#eta/#phi track p_{T} mapping",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-  fetaphiAllTracks = new TH2F("fetaphiAllTracks","#eta/#phi track mapping",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-  fptAllTracksCut = new TH1F("fptAllTracksCut", "Track transverse momentum after cuts [GeV]",300,0.,150.);
-  fetaAllTracksCut = new TH1F("fetaAllTracksCut", "#eta tracks after cuts", 36, -0.9, 0.9);
-  fphiAllTracksCut = new TH1F("fphiAllTracksCuts", "#phi tracks after cuts", 60, 0., 2*TMath::Pi());
-  fetaphiptAllTracksCut = new TH2F("fetaphiptAllTracksCuts","#eta/#phi track p_{T} mapping after cuts",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-  fetaphiAllTracksCut = new TH2F("fetaphiAllTracksCuts","#eta/#phi track mapping after cuts",36, -0.9, 0.9,60, 0., 2*TMath::Pi());
-
-  SetProperties(fptAllTracks, "Track p_{T} [GeV]", "dN/dp_{T}"); 
-  SetProperties(fetaAllTracks, "Track #eta", "dN/d#eta");
-  SetProperties(fphiAllTracks, "Track #phi", "dN/d#phi");
-  SetProperties(fetaphiptAllTracks, "#eta_{track}", "#phi_{track}");
-  SetProperties(fetaphiAllTracks, "#eta_{track}", "#phi_{track}");
-  SetProperties(fptAllTracksCut, "p_{T}track [GeV]", "dN/dp_{T}"); 
-  SetProperties(fetaAllTracksCut, "#eta_{track}", "dN/d#eta"); 
-  SetProperties(fphiAllTracksCut, "#phi_{track}", "dN/d#phi"); 
-  SetProperties(fetaphiptAllTracksCut, "#eta_{track}", "#phi_{track}");
-  SetProperties(fetaphiAllTracksCut, "#eta_{track}", "#phi_{track}");
-
-  fvertexXY = new TH2F("fvertexXY","X-Y vertex position",30,0.,10.,30,0.,10.);
-  fvertexZ = new TH1F("fvertexZ","Z vertex position",60,-30.,30.);
-  fEvtMult = new TH1F("fEvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",100,0.,100.);
-  fEvtMultvsJetPt = new TH2F("fEvtMultvsJetPt","Event multiplicity vs pT_{jet}",60,0.,60.,100,0.,100.);
-  fPtvsEtaJet = new TH2F("fPtvsEtaJet","Pt vs #eta_{jet}",20,-1.,1.,60,0.,60.);
-  fNpvsEtaJet = new TH2F("fNpvsEtaJet","N_{part} inside jet vs #eta_{jet}",20,-1.,1.,20,0,20); 
-  fNpevtvsEtaJet = new TH2F("fNpevtvsEtaJet","N_{part} in evt vs #eta_{jet}",20,-1.,1.,90,0,90); 
-  fPtvsPtJet = new TH2F("fPtvsPtJet","Pt vs #p_{Tjet}",60,0.,60.,60,0.,60.);
-  fNpvsPtJet = new TH2F("fNpvsPtJet","N_{part} inside jet vs #pt_{Tjet}",60,0.,60.,20,0,20); 
-  fNpevtvsPtJet = new TH2F("fNpevtvsPtJet","N_{part} in evt vs #pt_{Tjet}",60,0.,60.,90,0,90); 
-  fPtvsPtJet1D = new TH1F("fPtvsPtJet1D","Pt vs #p_{Tjet}",60,0.,60.);
-  fNpvsPtJet1D = new TH1F("fNpvsPtJet1D","N_{part} inside jet vs #pt_{Tjet}",60,0.,60.); 
-  fNpevtvsPtJet1D = new TH1F("fNpevtvsPtJet1D","N_{part} in evt vs #pt_{Tjet}",60,0.,60.); 
-  fptLeadingJet = new TH1F("fptLeadingJet","Pt leading Jet [GeV/c]",60,0.,60.); 
-  fetaLeadingJet = new TH1F("fetaLeadingJet","#eta leading jet",20,-1.,1.);
-  fphiLeadingJet = new TH1F("fphiLeadingJet","#phi leading jet",12,0.,2*TMath::Pi());
-  fptJet = new TH1F("fptJet","Pt Jets [GeV/c]",60,0.,60.); 
-  fetaJet = new TH1F("fetaJet","#eta jet",20,-1.,1.);
-  fphiJet = new TH1F("fphiJet","#phi jet",12,0.,2*TMath::Pi());
-  fNBadRunsH = new TH1F("fNBadRunsH","Number of events with Z vertex out of range", 1, 0., 1.);
-
-  SetProperties(fvertexXY, "vtx X", "vtx Y");
-  SetProperties(fvertexZ, "vtx Z", "Count");
-  SetProperties(fEvtMult, "N_{part} / event", "Count");
-  SetProperties(fEvtMultvsJetPt, "p_{T jet}", "Event multiplicity");
-  SetProperties(fPtvsEtaJet, "#eta_{leading jet}", "p_{T} part [GeV/c]");
-  SetProperties(fNpvsEtaJet, "#eta_{leading jet}", "N_{part} in leading jet");
-  SetProperties(fNpevtvsEtaJet, "#eta_{leading jet}", "N_{part} in event");
-  SetProperties(fPtvsPtJet, "#p_{T leading jet}", "p_{T} part [GeV/c]");
-  SetProperties(fNpvsPtJet, "#p_{T leading jet}", "N_{part} in leading jet");
-  SetProperties(fNpevtvsPtJet, "#p_{T leading jet}", "N_{part} in event");
-  SetProperties(fPtvsPtJet1D, "#p_{T leading jet}", "<p_{T}> part [GeV/c]");
-  SetProperties(fNpvsPtJet1D, "#p_{T leading jet}", "<N_{part}> in leading jet");
-  SetProperties(fNpevtvsPtJet1D, "#p_{T leading jet}", "<N_{part}> in event");
-  SetProperties(fptLeadingJet, "p_{T} leading jet", "dN/dp_{T} leading jet");
-  SetProperties(fetaLeadingJet, "#eta leading jet", "dN/d#eta leading jet");
-  SetProperties(fphiLeadingJet, "#phi leading jet", "dN/d#phi leading jet");
-  SetProperties(fptJet, "p_{T} jet [GeV/c]", "dN/dp_{T}");
-  SetProperties(fetaJet, "#eta jet", "dN/d#eta");
-  SetProperties(fphiJet, "#phi jet", "dN/d#phi");
-
+    
+  return *this;
 }
 
+//_________________________________________________________
+AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::~AliFragFuncQADiJetHistos()
+{
+  // destructor 
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+  if(fh2InvMass)  delete fh2InvMass;
+  if(fh2DeltaPhi) delete fh2DeltaPhi;
+  if(fh2DeltaEta) delete fh2DeltaEta;
+  if(fh2DeltaPt)  delete fh2DeltaPt;
+}
 
-void AliAnalysisTaskFragmentationFunction::FillMonoJetH(Int_t goodBin, AliAODJet* recJets, TClonesArray* tracks)
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::DefineQADiJetHistos()
 {
-  if (goodBin == 999) return;
-
-  Int_t nTracks = tracks->GetEntries();
-  Float_t xi,t1,ene,dr2,deta2,dphi2, z, cosTheta, theta, kt;
-  Bool_t jetOk1 = 0;
-  Bool_t jetOk2 = 0;
-  Bool_t jetOk3 = 0;
-  Bool_t jetOk4 = 0;
-  Bool_t jetOk5 = 0;
-  Bool_t jetOk6 = 0;
-  Bool_t jetOk7 = 0;
-  Bool_t jetOk8 = 0;
-  Bool_t jetOk9 = 0;
-  Bool_t jetOk10 = 0;
-
-  fEtaMonoJet1H[goodBin]->Fill(recJets[0].Eta());
-  fPhiMonoJet1H[goodBin]->Fill(recJets[0].Phi());
-  fPtMonoJet1H[goodBin]->Fill(recJets[0].Pt());
-  fEMonoJet1H[goodBin]->Fill(recJets[0].E());
-
-  Int_t mult = 0;
-  for(Int_t it=0; it<nTracks; it++)
-  {
-    AliAODTrack* aodTrack = (AliAODTrack*)tracks->At(it);
-    // Apply track cuts
-    UInt_t status = aodTrack->GetStatus();
-    if (status == 0) continue;
-    if((fFilterMask>0)&&!(aodTrack->TestFilterBit(fFilterMask)))continue;
-    mult++;
-    Float_t etaT = aodTrack->Eta();
-    Float_t phiT = aodTrack->Phi();
-    Float_t ptT = aodTrack->Pt();
-    // For Theta distribution
-    Float_t pxT = aodTrack->Px();
-    Float_t pyT = aodTrack->Py();
-    Float_t pzT = aodTrack->Pz();
-    Float_t pT = aodTrack->P();
-    // Compute theta
-    cosTheta = (pxT*recJets[0].Px()+pyT*recJets[0].Py()+pzT*recJets[0].Pz())/(pT*recJets[0].P());
-    theta = TMath::ACos(cosTheta);
-    // Compute xi
-    deta2 = etaT - recJets[0].Eta();
-    dphi2 = phiT - recJets[0].Phi();
-    if (dphi2 < -TMath::Pi()) dphi2= -dphi2 - 2.0 * TMath::Pi();
-    if (dphi2 > TMath::Pi()) dphi2 = 2.0 * TMath::Pi() - dphi2;
-    dr2 = TMath::Sqrt(deta2 * deta2 + dphi2 * dphi2);
-    t1  = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etaT)));
-    ene = ptT*TMath::Sqrt(1.+1./(t1*t1));
-    xi  = (Float_t) TMath::Log(recJets[0].E()/ptT);
-    // Compute z
-    z   = (Double_t)(ptT/recJets[0].E());
-    // Compute kT/jT
-    TVector3 partP; TVector3 jetP;
-    jetP[0] = recJets[0].Px();
-    jetP[1] = recJets[0].Py();
-    jetP[2] = recJets[0].Pz();
-    partP.SetPtEtaPhi(ptT,etaT,phiT);
-    kt = TMath::Sin(partP.Angle(jetP))*partP.Mag();
-    // Compute Jet shape
-
-    for(Int_t i2 = 0; i2 < fnRBin; i2++)
+
+  Int_t nBins = 0;
+  Double_t min = 0.;
+  Double_t max = 0.;
+  const char *xaxis = "";
+  if(fKindSlices == 1)
     {
-      if ((dr2<farrayRadii[i2]) && (ptT > fPartPtCut)) 
-      {
-        if (i2 == 0) jetOk1 = 1;
-        if (i2 == 1) jetOk2 = 1;
-        if (i2 == 2) jetOk3 = 1;
-        if (i2 == 3) jetOk4 = 1;
-        if (i2 == 4) jetOk5 = 1;
-        if (i2 == 5) jetOk6 = 1;
-        if (i2 == 6) jetOk7 = 1;
-        if (i2 == 7) jetOk8 = 1;
-        if (i2 == 8) jetOk9 = 1;
-        if (i2 == 9) jetOk10 = 1;
-
-        fdNdXiMonoJet1H[goodBin][i2]->Fill(xi);
-        fdNdPtMonoJet1H[goodBin][i2]->Fill(ptT);
-        fdNdZMonoJet1H[goodBin][i2]->Fill(z);
-       fdNdThetaMonoJet1H[goodBin][i2]->Fill(theta);
-       fdNdcosThetaMonoJet1H[goodBin][i2]->Fill(cosTheta);
-       fdNdkTMonoJet1H[goodBin][i2]->Fill(kt);
-       fdNdpTvsZMonoJet1H[goodBin][i2]->Fill(z,1/((fPthadHBinMax-fPthadHBinMin)/fnPthadHBin));
-
-       fThetaPtPartMonoJet1H[goodBin][i2]->Fill(ptT,theta);
-       fcosThetaPtPartMonoJet1H[goodBin][i2]->Fill(ptT,cosTheta);
-       fkTPtPartMonoJet1H[goodBin][i2]->Fill(ptT,kt);
-       fThetaPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),theta);
-       fcosThetaPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),cosTheta);
-       fkTPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),kt);
-       fpTPtJetMonoJet1H[goodBin][i2]->Fill(recJets[0].Pt(),ptT);
-      }
+      nBins = fNBinsJetInvMass;
+      min   = fJetInvMassMin;
+      max   = fJetInvMassMax;
+      xaxis = "M_{JJ} [GeV]";
     }
-  }
-  fEvtMultvsJetPt->Fill(recJets[0].Pt(),mult);
-
-  if (jetOk1)  fNMonoJet1sH[goodBin][0]->Fill(0.5);
-  if (jetOk2)  fNMonoJet1sH[goodBin][1]->Fill(0.5);
-  if (jetOk3)  fNMonoJet1sH[goodBin][2]->Fill(0.5);
-  if (jetOk4)  fNMonoJet1sH[goodBin][3]->Fill(0.5);
-  if (jetOk5)  fNMonoJet1sH[goodBin][4]->Fill(0.5);
-  if (jetOk6)  fNMonoJet1sH[goodBin][5]->Fill(0.5);
-  if (jetOk7)  fNMonoJet1sH[goodBin][6]->Fill(0.5);
-  if (jetOk8)  fNMonoJet1sH[goodBin][7]->Fill(0.5);
-  if (jetOk9)  fNMonoJet1sH[goodBin][8]->Fill(0.5);
-  if (jetOk10) fNMonoJet1sH[goodBin][9]->Fill(0.5);
+  if(fKindSlices == 2 || fKindSlices == 3)
+    {
+      nBins = fNBinsJetPt;
+      min   = fJetPtMin;
+      max   = fJetPtMax;
+      if(fKindSlices == 2) xaxis = "E_{Tmean} [GeV]";
+      if(fKindSlices == 3) xaxis ="leading jet p_{T} [GeV/c]";
+    }
+  
+  
+  fh2InvMass  = new TH2F(Form("fh2DJInvMassPositionCut%s",  fName.Data()), "",nBins, min, max, fNBinsJetInvMass, fJetInvMassMin, fJetInvMassMax);
+  fh2DeltaPhi = new TH2F(Form("fh2DJDeltaPhiPositionCut%s", fName.Data()), "",nBins, min, max, fNBinsDeltaPhi, fDeltaPhiMin, fDeltaPhiMax);
+  fh2DeltaEta = new TH2F(Form("fh2DJDeltaEtaPositionCut%s", fName.Data()), "",nBins, min, max, fNBinsDeltaEta, fDeltaEtaMin, fDeltaEtaMax);
+  fh2DeltaPt  = new TH2F(Form("fh2DJDeltaPtPositionCut%s",  fName.Data()), "",nBins, min, max, fNBinsDeltaPt, fDeltaPtMin, fDeltaPtMax);
+  
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2InvMass, xaxis, "Invariant Mass", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaPhi, xaxis, "#Delta #phi", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaEta, xaxis, "#Delta #eta", "Entries");
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaPt, xaxis, "#Delta p_{T}", "Entries");
 
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
-void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::FillDiJetQA(Double_t invMass, Double_t deltaPhi, Double_t deltaEta,Double_t deltaPt, Double_t jetBin)
 {
-  //Set properties of histos (x and y title and error propagation)
-  h->SetXTitle(x);
-  h->SetYTitle(y);
-  h->GetXaxis()->SetTitleColor(1);
-  h->GetYaxis()->SetTitleColor(1);
-  h->Sumw2();
+  fh2InvMass->Fill(jetBin, invMass);
+  fh2DeltaPhi->Fill(jetBin, deltaPhi);
+  fh2DeltaEta->Fill(jetBin, deltaEta);
+  fh2DeltaPt->Fill(jetBin, deltaPt);
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+//________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AddToOutput(TList* list)const
+{
+  list->Add(fh2InvMass);
+  list->Add(fh2DeltaPhi);
+  list->Add(fh2DeltaEta);
+  list->Add(fh2DeltaPt);
+}
 
-void AliAnalysisTaskFragmentationFunction::MakeJetContainer()
+//_________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AddToOutput(TList* list) const
 {
-  //
-  // Create the particle container for the correction framework manager and 
-  // link it
-  //
-  const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
-  const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
-  const Double_t kEtamin = -3.0, kEtamax = 3.0;
-  const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
-
-  // can we neglect migration in eta and phi?
-  // phi should be no problem since we cover full phi and are phi symmetric
-  // eta migration is more difficult  due to needed acceptance correction
-  // in limited eta range
-
-  //arrays for the number of bins in each dimension
-  Int_t iBin[kNvar];
-  iBin[0] = 100; //bins in pt
-  iBin[1] =  1; //bins in eta
-  iBin[2] = 1; // bins in phi
-
-  //arrays for lower bounds :
-  Double_t* binEdges[kNvar];
-  for(Int_t ivar = 0; ivar < kNvar; ivar++)
-    binEdges[ivar] = new Double_t[iBin[ivar] + 1];
-
-  //values for bin lower bounds
-  //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);  
-  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
-  for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
-  for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
-
-  /*
-  for(int i = 0;i < kMaxStep*2;++i){
-    fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
-    for (int k=0; k<kNvar; k++) {
-      fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
-    }
-  }
-  //create correlation matrix for unfolding
-  Int_t thnDim[2*kNvar];
-  for (int k=0; k<kNvar; k++) {
-    //first half  : reconstructed 
-    //second half : MC
-    thnDim[k]      = iBin[k];
-    thnDim[k+kNvar] = iBin[k];
-  }
+  // add histos to list
+
+  list->Add(fh2CosTheta);
+  list->Add(fh2Theta);
+  list->Add(fh2Jt);
+
+  list->Add(fhnIntraJet);
 
-  fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
-  for (int k=0; k<kNvar; k++) {
-    fhnCorrelation->SetBinEdges(k,binEdges[k]);
-    fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
-  }
-  fhnCorrelation->Sumw2();
-
-  // Add a histogram for Fake jets
-  //  thnDim[3] = AliPID::kSPECIES;
-  //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
-  // for(Int_t idim = 0; idim < kNvar; idim++)
-  //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
-  */
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+
+Bool_t AliAnalysisTaskFragmentationFunction::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // (taken from AliAnalysisTaskJetSpectrum)
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Double_t xsection = 0;
+  UInt_t   ntrials  = 0;
+  Float_t   ftrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "");
+    }
+    else if(fileName.Contains("AliAODs.root")){
+        fileName.ReplaceAll("AliAODs.root", "");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "");
+    }
+    TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
+    if(!fxsec){
+      if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
+      // next trial fetch the histgram file
+      fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+      if(!fxsec){
+       // not a severe condition
+       if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));        
+       return kTRUE;
+      }
+      else{
+       // find the tlist we want to be independtent of the name so use the Tkey
+       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
+       if(!key){
+         if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);       
+         return kTRUE;
+       }
+       TList *list = dynamic_cast<TList*>(key->ReadObj());
+       if(!list){
+         if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);      
+         return kTRUE;
+       }
+       xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+       ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+      }
+    }
+    else{
+      TTree *xtree = (TTree*)fxsec->Get("Xsection");
+      if(!xtree){
+       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+      }
+      xtree->SetBranchAddress("xsection",&xsection);
+      xtree->SetBranchAddress("ntrials",&ntrials);
+      ftrials = ntrials;
+      xtree->GetEntry(0);
+    }
+    fh1Xsec->Fill("<#sigma>",xsection);
+    fh1Trials->Fill("#sum{ntrials}",ftrials);
+  }
+  
+  return kTRUE;
+}
+
+
+
+//__________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
+{
+  // create output objects
+
+  if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
+  // create list of tracks and jets 
+  
+  fTracksRec = new TList();
+  fTracksRec->SetOwner(kFALSE);  
+
+  fTracksRecCuts = new TList();
+  fTracksRecCuts->SetOwner(kFALSE);  
+
+  fTracksGen = new TList();
+  fTracksGen->SetOwner(kFALSE);
+
+  fTracksAODMCCharged = new TList();
+  fTracksAODMCCharged->SetOwner(kFALSE);
+    
+  fTracksRecQualityCuts = new TList(); 
+  fTracksRecQualityCuts->SetOwner(kFALSE);
+
+  fJetsRec = new TList();
+  fJetsRec->SetOwner(kFALSE);
+
+  fJetsRecCuts = new TList();
+  fJetsRecCuts->SetOwner(kFALSE);
+
+  fJetsGen = new TList();
+  fJetsGen->SetOwner(kFALSE);
+
+  fJetsRecEff = new TList();
+  fJetsRecEff->SetOwner(kFALSE);
+
+  // fJetsKine = new TList();
+  // fJetsKine->SetOwner(kTRUE); // delete AOD jets using mom from Kine Tree via TList::Clear()
+
+
+  //
+  // Create histograms / output container
+  //
+
+  OpenFile(1);
+  fCommonHistList = new TList();
+  
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+  
+  
+  // Histograms        
+  fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
+  fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
+  fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
+  fh1EvtMult                = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,120.);
+
+  fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
+  fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
+
+  fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
+  fh1nGenJets                = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5);
+  fh1nRecEffJets             = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5);
+
+  // 5D single track eff histo: phi:eta:gen pt:rec pt:isReconstructed -  use binning as for track QA
+  Int_t   nBinsSingleTrackEffHisto[5]      = { fQATrackNBinsPhi, fQATrackNBinsEta, fQATrackNBinsPt, fQATrackNBinsPt, 2 };
+  Double_t binMinSingleTrackEffHisto[5]    = { fQATrackPhiMin, fQATrackEtaMin, fQATrackPtMin, fQATrackPtMin, 0 };
+  Double_t binMaxSingleTrackEffHisto[5]    = { fQATrackPhiMax, fQATrackEtaMax, fQATrackPtMax, fQATrackPtMax, 2 };
+  const char* labelsSingleTrackEffHisto[5] = {"#phi","#eta","gen p_{T} [GeV/c]", "rec p_{T} [GeV/c]", "isRec"};
+
+  fhnSingleTrackRecEffHisto = new THnSparseF("fhnSingleTrackRecEffHisto","generated tracks phi:eta:pt:isReconstructed",5,
+                                            nBinsSingleTrackEffHisto,binMinSingleTrackEffHisto,binMaxSingleTrackEffHisto);
+  
+  AliAnalysisTaskFragmentationFunction::SetProperties(fhnSingleTrackRecEffHisto,5,labelsSingleTrackEffHisto); 
+  
+  
+  // 7D jets track eff histo: jet phi:eta:pt:track pt:z:xi:isReconstructed - use binning as for track/jet QA
+  Int_t    nBinsJetTrackEffHisto[7]     = { fQAJetNBinsPhi, fQAJetNBinsEta, fFFNBinsJetPt, fFFNBinsPt, fFFNBinsZ, fFFNBinsXi, 2};
+  Double_t binMinJetTrackEffHisto[7]    = { fQAJetPhiMin, fQAJetEtaMin, fFFJetPtMin , fFFPtMin, fFFZMin ,  fFFXiMin, 0 };
+  Double_t binMaxJetTrackEffHisto[7]    = { fQAJetPhiMax, fQAJetEtaMax, fFFJetPtMax , fFFPtMax, fFFZMax ,  fFFXiMax, 2 };
+  const char* labelsJetTrackEffHisto[7] = {"jet #phi","jet #eta","jet p_{T} [GeV/c]","track p_{T} [GeV/c]","z","#xi","isRec"};
+
+  fhnJetTrackRecEffHisto       = new THnSparseF("fhnJetTrackRecEffHisto","generated tracks - jet phi:jet eta:jet pt:track pt:z:xi:isReconstructed",7,
+                                               nBinsJetTrackEffHisto,binMinJetTrackEffHisto,binMaxJetTrackEffHisto);
+
+  AliAnalysisTaskFragmentationFunction::SetProperties(fhnJetTrackRecEffHisto,7,labelsJetTrackEffHisto);
+
+
+  fQATrackHistosRec          = new AliFragFuncQATrackHistos("Rec", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
+                                                           fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
+                                                           fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
+                                                           fQATrackHighPtThreshold);
+  fQATrackHistosRecCuts      = new AliFragFuncQATrackHistos("RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
+                                                           fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
+                                                           fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
+                                                           fQATrackHighPtThreshold);
+  fQATrackHistosGen          = new AliFragFuncQATrackHistos("Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
+                                                           fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
+                                                           fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
+                                                           fQATrackHighPtThreshold);
+  
+
+  fQAJetHistosRec            = new AliFragFuncQAJetHistos("Rec", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
+                                                         fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+  fQAJetHistosRecCuts        = new AliFragFuncQAJetHistos("RecCuts", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
+                                                         fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+  fQAJetHistosRecCutsLeading = new AliFragFuncQAJetHistos("RecCutsLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
+                                                         fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+  fQAJetHistosGen            = new AliFragFuncQAJetHistos("Gen", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
+                                                         fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+  fQAJetHistosGenLeading     = new AliFragFuncQAJetHistos("GenLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
+                                                         fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);  
+  fQAJetHistosRecEffLeading  = new AliFragFuncQAJetHistos("RecEffLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
+                                                         fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+
+
+  fFFHistosRecCuts          = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  fFFHistosRecLeading        = new AliFragFuncHistos("RecLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  fFFHistosRecLeadingTrack   = new AliFragFuncHistos("RecLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  fFFHistosGen              = new AliFragFuncHistos("Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  fFFHistosGenLeading        = new AliFragFuncHistos("GenLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  fFFHistosGenLeadingTrack   = new AliFragFuncHistos("GenLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+
+
+  fIJHistosRecCuts          = new AliFragFuncIntraJetHistos("RecCuts", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax,
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+  fIJHistosRecLeading        = new AliFragFuncIntraJetHistos("RecLeading", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+  fIJHistosRecLeadingTrack   = new AliFragFuncIntraJetHistos("RecLeadingTrack", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+  fIJHistosGen              = new AliFragFuncIntraJetHistos("Gen", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+  fIJHistosGenLeading        = new AliFragFuncIntraJetHistos("GenLeading", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+  fIJHistosGenLeadingTrack   = new AliFragFuncIntraJetHistos("GenLeadingTrack", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
+                                                            fIJNBinsPt, fIJPtMin, fIJPtMax, 
+                                                            fIJNBinsZ, fIJZMin, fIJZMax,  
+                                                            fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
+                                                            fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
+                                                            fIJNBinsJt , fIJJtMin , fIJJtMax);
+
+
+  fFFDiJetHistosRecCuts         = new AliFragFuncDiJetHistos("RecCuts", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
+  fFFDiJetHistosRecLeading      = new AliFragFuncDiJetHistos("RecLeading", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax); 
+  fFFDiJetHistosRecLeadingTrack = new AliFragFuncDiJetHistos("RecLeadingTrack", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
+
+  fFFDiJetHistosGen             = new AliFragFuncDiJetHistos("Gen", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax,
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
+  fFFDiJetHistosGenLeading      = new AliFragFuncDiJetHistos("GenLeading", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax,
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax,
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
+  fFFDiJetHistosGenLeadingTrack = new AliFragFuncDiJetHistos("GenLeadingTrack", fDiJetKindBins, 
+                                                            fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
+                                                            fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
+                                                            fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
+                                                            fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
+                                                            fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
+
+  fQADiJetHistosRecCuts = new AliFragFuncQADiJetHistos("RecCuts", fDiJetKindBins, 
+                                                      fQADiJetNBinsInvMass, fQADiJetInvMassMin, fQADiJetInvMassMax,
+                                                      fQADiJetNBinsJetPt, fQADiJetJetPtMin, fQADiJetJetPtMax,
+                                                      fQADiJetNBinsDeltaPhi, fQADiJetDeltaPhiMin, fQADiJetDeltaPhiMax , 
+                                                      fQADiJetNBinsDeltaEta, fQADiJetDeltaEtaMin, fQADiJetDeltaEtaMax , 
+                                                      fQADiJetNBinsDeltaPt, fQADiJetDeltaPtMin, fQADiJetDeltaPtMax);
+  fQADiJetHistosGen     = new AliFragFuncQADiJetHistos("Gen", fDiJetKindBins, 
+                                                      fQADiJetNBinsInvMass, fQADiJetInvMassMin,  fQADiJetInvMassMax, 
+                                                      fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
+                                                      fQADiJetNBinsDeltaPhi, fQADiJetDeltaPhiMin, fQADiJetDeltaPhiMax,
+                                                      fQADiJetNBinsDeltaEta, fQADiJetDeltaEtaMin, fQADiJetDeltaEtaMax,
+                                                      fQADiJetNBinsDeltaPt, fQADiJetDeltaPtMin, fQADiJetDeltaPtMax);
+
+
+  fQATrackHistosRec->DefineHistos();
+  fQATrackHistosRecCuts->DefineHistos();
+  fQATrackHistosGen->DefineHistos();
+
+  fQAJetHistosRec->DefineHistos();
+  fQAJetHistosRecCuts->DefineHistos();
+  fQAJetHistosRecCutsLeading->DefineHistos();
+  fQAJetHistosGen->DefineHistos();
+  fQAJetHistosGenLeading->DefineHistos();
+  fQAJetHistosRecEffLeading->DefineHistos();
+
+  fFFHistosRecCuts->DefineHistos();
+  fFFHistosRecLeading->DefineHistos();
+  fFFHistosRecLeadingTrack->DefineHistos();
+  fFFHistosGen->DefineHistos();
+  fFFHistosGenLeading->DefineHistos();
+  fFFHistosGenLeadingTrack->DefineHistos();
+
+  fIJHistosRecCuts->DefineHistos();
+  fIJHistosRecLeading->DefineHistos();
+  fIJHistosRecLeadingTrack->DefineHistos();
+  fIJHistosGen->DefineHistos();
+  fIJHistosGenLeading->DefineHistos();
+  fIJHistosGenLeadingTrack->DefineHistos();
+  
+  fFFDiJetHistosRecCuts->DefineDiJetHistos();
+  fFFDiJetHistosRecLeading->DefineDiJetHistos();
+  fFFDiJetHistosRecLeadingTrack->DefineDiJetHistos();
+  fFFDiJetHistosGen->DefineDiJetHistos();
+  fFFDiJetHistosGenLeading->DefineDiJetHistos();
+  fFFDiJetHistosGenLeadingTrack->DefineDiJetHistos();
+  fQADiJetHistosRecCuts->DefineQADiJetHistos();
+  fQADiJetHistosGen->DefineQADiJetHistos();
+
+  Bool_t genJets    = (fJetTypeGen != kJetsUndef) ? kTRUE : kFALSE;
+  Bool_t genTracks  = (fTrackTypeGen != kTrackUndef) ? kTRUE : kFALSE;
+  Bool_t recJetsEff = (fJetTypeRecEff != kJetsUndef) ? kTRUE : kFALSE;
+
+
+  const Int_t saveLevel = 5;
+  if(saveLevel>0){
+    fCommonHistList->Add(fh1EvtSelection);
+    fFFHistosRecCuts->AddToOutput(fCommonHistList);
+    fFFHistosRecLeading->AddToOutput(fCommonHistList);
+    fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
+
+    if(genJets && genTracks){
+       fFFHistosGen->AddToOutput(fCommonHistList);
+       fFFHistosGenLeading->AddToOutput(fCommonHistList);
+       fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+       fCommonHistList->Add(fh1Xsec);
+       fCommonHistList->Add(fh1Trials);
+       fCommonHistList->Add(fh1PtHard);
+       fCommonHistList->Add(fh1PtHardTrials);
+    }
+  }
+  if(saveLevel>1){
+    fQATrackHistosRec->AddToOutput(fCommonHistList);
+    fQATrackHistosRecCuts->AddToOutput(fCommonHistList);
+    if(genTracks) fQATrackHistosGen->AddToOutput(fCommonHistList);
+    
+    fQAJetHistosRec->AddToOutput(fCommonHistList);
+    fQAJetHistosRecCuts->AddToOutput(fCommonHistList);
+    fQAJetHistosRecCutsLeading->AddToOutput(fCommonHistList);
+    if(recJetsEff) fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList); 
+
+    if(genJets){
+       fQAJetHistosGen->AddToOutput(fCommonHistList);
+       fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
+    }
+
+    fCommonHistList->Add(fh1EvtMult);
+    fCommonHistList->Add(fh1nRecJetsCuts);
+    if(genJets) fCommonHistList->Add(fh1nGenJets);
+  }
+  if(saveLevel>2){
+    fCommonHistList->Add(fh1VertexNContributors);
+    fCommonHistList->Add(fh1VertexZ);    
+  }
+  if(saveLevel>3){
+    fIJHistosRecCuts->AddToOutput(fCommonHistList);
+    fIJHistosRecLeading->AddToOutput(fCommonHistList);
+    fIJHistosRecLeadingTrack->AddToOutput(fCommonHistList);
+
+    if(genJets && genTracks){
+       fIJHistosGen->AddToOutput(fCommonHistList);
+       fIJHistosGenLeading->AddToOutput(fCommonHistList);
+       fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+    }
+  }
+  if(saveLevel>4){
+    fFFDiJetHistosRecCuts->AddToOutput(fCommonHistList);
+    fFFDiJetHistosRecLeading->AddToOutput(fCommonHistList);
+    fFFDiJetHistosRecLeadingTrack->AddToOutput(fCommonHistList);
+    fQADiJetHistosRecCuts->AddToOutput(fCommonHistList);
+
+    if(genJets && genTracks){
+        fFFDiJetHistosGen->AddToOutput(fCommonHistList);
+        fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
+        fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+        fQADiJetHistosGen->AddToOutput(fCommonHistList);
+    }
+
+    if(recJetsEff && genTracks){
+       fCommonHistList->Add(fhnSingleTrackRecEffHisto);
+       fCommonHistList->Add(fhnJetTrackRecEffHisto);
+       fCommonHistList->Add(fh1nRecEffJets);
+    }
+  }
+
+  // =========== Switch on Sumw2 for all histos ===========
+  for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
+    TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
+    if (h1) h1->Sumw2();
+    else{
+      THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
+      if(hnSparse) hnSparse->Sumw2();
+    }
+  }
+  
+  TH1::AddDirectory(oldStatus);
+}
+
+//_______________________________________________
+void AliAnalysisTaskFragmentationFunction::Init()
+{
+  // Initialization
+  if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::Init()");
+
+}
+
+//_____________________________________________________________
+void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserExec()");
+       
+
+  if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
+  // Trigger selection
+  
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+  if(inputHandler->IsEventSelected() & AliVEvent::kMB){
+    if(fDebug > 1)  Printf(" Trigger Selection: event ACCEPTED ... ");
+    fh1EvtSelection->Fill(1.);
+  } else {
+    fh1EvtSelection->Fill(0.);
+    if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
+      if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
+      PostData(1, fCommonHistList);
+      return;
+    }
+  }
+  
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if(!fESD){
+    if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
+  }
+  
+  fMCEvent = MCEvent();
+  if(!fMCEvent){
+    if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
+  }
+  
+  // get AOD event from input/ouput
+  TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+  if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
+    fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
+    if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
+  }
+  else {
+    handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
+    if( handler && handler->InheritsFrom("AliAODHandler") ) {
+      fAOD  = ((AliAODHandler*)handler)->GetAOD();
+      if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
+    }
+  }
+  
+  if(!fAOD){
+    Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
+    return;
+  }
+  
+
+  // event selection (vertex) *****************************************
+  
+  AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
+  Int_t nTracksPrim = primVtx->GetNContributors();
+  fh1VertexNContributors->Fill(nTracksPrim);
+  
+  if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
+  if(!nTracksPrim){
+    if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
+    fh1EvtSelection->Fill(2.);
+    PostData(1, fCommonHistList);
+    return;
+  }
+
+  fh1VertexZ->Fill(primVtx->GetZ());
+  
+  if(TMath::Abs(primVtx->GetZ())>10){
+    if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
+    fh1EvtSelection->Fill(3.);
+    PostData(1, fCommonHistList);
+    return; 
+  }
+
+  TString primVtxName(primVtx->GetName());
+
+  if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
+    if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
+    fh1EvtSelection->Fill(4.);
+    PostData(1, fCommonHistList);
+    return;
+  }
+  if (fDebug > 1) Printf("%s:%d primary vertex selection: event ACCEPTED ...",(char*)__FILE__,__LINE__); 
+  fh1EvtSelection->Fill(5.);
+
+
+  //___ get MC information __________________________________________________________________
+
+  Double_t ptHard = 0.;
+  Double_t nTrials = 1; // trials for MC trigger weight for real data
+  AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
+  if(!pythiaGenHeader){
+     if(fJetTypeGen != kJetsUndef && fTrackTypeGen != kTrackUndef){
+        Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
+        return;
+     }
+  } else {
+     nTrials = pythiaGenHeader->Trials();
+     ptHard  = pythiaGenHeader->GetPtHard();
+
+     fh1PtHard->Fill(ptHard);
+     fh1PtHardTrials->Fill(ptHard,nTrials);
+  }
+  
+  
+  //___ fetch jets __________________________________________________________________________
+  
+  Int_t nJ = GetListOfJets(fJetsRec, kJetsRec);
+  Int_t nRecJets = 0;
+  if(nJ>=0) nRecJets = fJetsRec->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Rec jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
+  if(nJ != nRecJets) Printf("%s:%d Mismatch Selected Rec Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
+
+  Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
+  Int_t nRecJetsCuts = 0;
+  if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
+  if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
+  fh1nRecJetsCuts->Fill(nRecJetsCuts);
+
+  
+  if(fJetTypeGen==kJetsKine || fJetTypeGen == kJetsKineAcceptance) fJetsGen->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear() 
+  Int_t nJGen  = GetListOfJets(fJetsGen, fJetTypeGen);
+  Int_t nGenJets = 0;
+  if(nJGen>=0) nGenJets = fJetsGen->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
+  if(nJGen != nGenJets) Printf("%s:%d Mismatch selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
+  fh1nGenJets->Fill(nGenJets);
+
+
+  if(fJetTypeRecEff==kJetsKine || fJetTypeRecEff == kJetsKineAcceptance) fJetsRecEff->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear() 
+  Int_t nJRecEff  = GetListOfJets(fJetsRecEff, fJetTypeRecEff);
+  Int_t nRecEffJets = 0;
+  if(nJRecEff>=0) nRecEffJets = fJetsRecEff->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
+  if(nJRecEff != nRecEffJets) Printf("%s:%d Mismatch selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
+  fh1nRecEffJets->Fill(nRecEffJets);
+
+
+  //____ fetch particles __________________________________________________________
+  
+  Int_t nT = GetListOfTracks(fTracksRec, kTrackAOD);
+  Int_t nRecPart = 0;
+  if(nT>=0) nRecPart = fTracksRec->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart);
+  if(nRecPart != nT) Printf("%s:%d Mismatch selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart);
+  
+
+  Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
+  Int_t nRecPartCuts = 0;
+  if(nTCuts>=0) nRecPartCuts = fTracksRecCuts->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
+  if(nRecPartCuts != nTCuts) Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
+  fh1EvtMult->Fill(nRecPartCuts);
+
+
+  Int_t nTGen = GetListOfTracks(fTracksGen,fTrackTypeGen);
+  Int_t nGenPart = 0;
+  if(nTGen>=0) nGenPart = fTracksGen->GetEntries();
+  if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
+  if(nGenPart != nTGen) Printf("%s:%d Mismatch selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
+  
+  
+  //____ analysis, fill histos ___________________________________________________
+  
+  // loop over tracks
+
+  for(Int_t it=0; it<nRecPart; ++it){
+    AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksRec->At(it));
+    fQATrackHistosRec->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
+  }
+  for(Int_t it=0; it<nRecPartCuts; ++it){
+    AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksRecCuts->At(it));
+    fQATrackHistosRecCuts->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
+  }
+  for(Int_t it=0; it<nGenPart; ++it){
+    AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksGen->At(it));
+    fQATrackHistosGen->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
+  }
+  
+  // loop over jets
+
+  for(Int_t ij=0; ij<nRecJets; ++ij){
+
+    AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRec->At(ij));
+    fQAJetHistosRec->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
+  }
+  
+
+  for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
+
+    AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
+    fQAJetHistosRecCuts->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
+
+    if(ij==0){ // leading jet
+      
+      fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
+      
+      TList* jettracklist = new TList();
+      Double_t sumPt = 0.;
+      Float_t leadTrackPx = 0.;
+      Float_t leadTrackPy = 0.;
+      Float_t leadTrackPz = 0.;
+      Float_t leadTrackP  = 0.;
+      Float_t leadTrackPt = 0.;
+      TLorentzVector* leadTrackV = new TLorentzVector();
+      
+      if(GetFFRadius()<=0){
+       GetJetTracksTrackrefs(jettracklist, jet);
+       } else {
+       GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
+      }
+      
+      for(Int_t it=0; it<jettracklist->GetSize(); ++it){
+       Float_t trackPx = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Px();
+       Float_t trackPy = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Py();
+       Float_t trackPz = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Pz();
+       Float_t trackP = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->P();
+       Float_t trackPt = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Pt();
+       Float_t jetPx = jet->Px();
+       Float_t jetPy = jet->Py();
+       Float_t jetPz = jet->Pz();
+       Float_t jetP  = jet->P();
+       Float_t jetPt = jet->Pt();
+       TLorentzVector* trackV = new TLorentzVector();
+       TLorentzVector *jetV = new TLorentzVector();
+       trackV->SetPxPyPzE(trackPx,trackPy,trackPz,trackP);
+       jetV->SetPxPyPzE(jetPx,jetPy,jetPz,jetP);
+
+       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+       
+       fFFHistosRecCuts->FillFF( trackPt, jetPt, incrementJetPt);
+       fIJHistosRecCuts->FillIntraJet( trackV, jetV );
+       
+       if(it==0){ 
+         leadTrackPx = trackPx;
+         leadTrackPy = trackPy;
+         leadTrackPz = trackPz;
+         leadTrackP  = trackP;
+         leadTrackPt = trackPt;
+         fFFHistosRecLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE);
+
+         leadTrackV->SetPxPyPzE(leadTrackPx,leadTrackPy,leadTrackPz,leadTrackP);
+         fIJHistosRecLeadingTrack->FillIntraJet( leadTrackV, jetV );
+       }
+       fFFHistosRecLeading->FillFF( trackPt, leadTrackPt , incrementJetPt);
+       fIJHistosRecLeading->FillIntraJet( trackV, leadTrackV );
+
+       delete trackV;
+       delete jetV;
+      }
+      
+      delete leadTrackV;
+      delete jettracklist;
+    }
+  }
+       
+
+  for(Int_t ij=0; ij<nGenJets; ++ij){
+
+    AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsGen->At(ij));
+    fQAJetHistosGen->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
+    
+    if(ij==0){ // leading jet
+    
+      fQAJetHistosGenLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
+      
+      TList* jettracklist = new TList();
+      Double_t sumPt = 0.;
+      Float_t leadTrackPx = 0.;
+      Float_t leadTrackPy = 0.;
+      Float_t leadTrackPz = 0.;
+      Float_t leadTrackP  = 0.;
+      Float_t leadTrackPt = 0.;
+      TLorentzVector* leadTrackV = new TLorentzVector();
+
+      if(GetFFRadius()<=0){
+       GetJetTracksTrackrefs(jettracklist, jet);
+      } else {
+       GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt);
+      }
+      
+      for(Int_t it=0; it<jettracklist->GetSize(); ++it){
+       Float_t trackPx = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Px();
+       Float_t trackPy = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Py();
+       Float_t trackPz = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Pz();
+       Float_t trackP  = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->P();
+       Float_t trackPt = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Pt();
+       Float_t jetPx = jet->Px();
+       Float_t jetPy = jet->Py();
+       Float_t jetPz = jet->Pz();
+       Float_t jetP  = jet->P();
+       Float_t jetPt = jet->Pt();
+       TLorentzVector* trackV = new TLorentzVector();
+       TLorentzVector  *jetV = new TLorentzVector();
+       trackV->SetPxPyPzE(trackPx,trackPy,trackPz,trackP);
+       jetV->SetPxPyPzE(jetPx,jetPy,jetPz,jetP);
+
+       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+
+       fFFHistosGen->FillFF( trackPt, jetPt, incrementJetPt);
+       fIJHistosGen->FillIntraJet( trackV, jetV );
+       
+       if(it==0){ 
+         leadTrackPx = trackPx;
+         leadTrackPy = trackPy;
+         leadTrackPz = trackPz;
+         leadTrackP  = trackP;
+         leadTrackPt = trackPt;
+         fFFHistosGenLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE);
+
+         leadTrackV->SetPxPyPzE(leadTrackPx,leadTrackPy,leadTrackPz,leadTrackP);
+         fIJHistosGenLeadingTrack->FillIntraJet( leadTrackV, jetV );
+       }
+       fFFHistosGenLeading->FillFF( trackPt, leadTrackPt, incrementJetPt);
+       fIJHistosGenLeading->FillIntraJet( trackV, leadTrackV );
+
+       delete trackV;
+       delete jetV;
+      }
+      
+      delete leadTrackV;
+      delete jettracklist;
+    }
+  }
+
+  //_______ DiJet part _____________________________________________________
+
+  if (nRecJetsCuts > 1)  // OB: debugged this
+  {
+
+    AliAODJet* jet1 = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(0));
+    AliAODJet* jet2 = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(1));
+    
+    // DiJet deltaphi calculation
+    Double_t phi1 = TVector2::Phi_0_2pi(jet1->Phi());
+    Double_t phi2 = TVector2::Phi_0_2pi(jet2->Phi());
+    Double_t deltaPhi = TMath::Abs(phi1-phi2); 
+    if (deltaPhi > TMath::Pi() && deltaPhi < 2*TMath::Pi()) deltaPhi = 2*TMath::Pi() - deltaPhi;
+    
+    // DiJet CDF cut calculation
+    Double_t Et1     = TMath::Abs(jet1->E()*TMath::Sin(jet1->Theta()));
+    Double_t Et2     = TMath::Abs(jet2->E()*TMath::Sin(jet2->Theta()));
+    Double_t sumEt   = Et1 + Et2;
+    Double_t normEt1PlusEt2   = TMath::Sqrt(Et1*Et1+Et2*Et2+2*Et1*Et2*TMath::Cos(deltaPhi));
+    Double_t ratio = (Double_t)(normEt1PlusEt2/sumEt);
+    
+    // DiJet events selection
+    Bool_t positionCut       = 0;
+    Bool_t positionEnergyCut = 0;
+    Bool_t cdfCut            = 0; 
+
+    // Position cut :
+    if (deltaPhi > fDiJetDeltaPhiCut) positionCut = 1;
+    // Position-Energy cut :
+    if ((deltaPhi > fDiJetDeltaPhiCut) && ((jet2->Pt()) >= fDiJetPtFractionCut*(jet1->Pt()))) positionEnergyCut = 1;
+    // CDF cut :
+    if (ratio < fDiJetCDFCut) cdfCut = 1;
+    
+    Int_t go = 0;
+    
+    if (fDiJetCut == 1 && positionCut == 1) go = 1;
+    if (fDiJetCut == 2 && positionEnergyCut == 1) go = 1;
+    if (fDiJetCut == 3 && cdfCut == 1) go = 1;
+
+    if (go)
+      {
+       Double_t deltaEta      = TMath::Abs(jet1->Eta()-jet2->Eta());
+       Double_t deltaPt       = TMath::Abs(jet1->Pt()-jet2->Pt());
+       Double_t meanEt        = (Double_t)((Et1+Et2)/2.);
+       Double_t invariantMass = (Double_t)InvMass(jet1,jet2);
+       
+       Double_t  jetBin = GetDiJetBin(invariantMass, jet1->Pt(), meanEt, fDiJetKindBins);
+
+       if (jetBin > 0)
+         {
+           fQADiJetHistosRecCuts->FillDiJetQA(invariantMass, deltaPhi, deltaEta, deltaPt, jetBin);
+           
+           TList* jettracklist1 = new TList();
+           Double_t sumPt1      = 0.;
+           Float_t leadTrackPt1 = 0;
+           
+           TList* jettracklist2 = new TList();
+           Double_t sumPt2      = 0.;
+           Float_t leadTrackPt2 = 0;
+           
+           if(GetFFRadius()<=0)
+             {
+               GetJetTracksTrackrefs(jettracklist1, jet1);
+               GetJetTracksTrackrefs(jettracklist2, jet2);
+             }
+           else
+             {
+               GetJetTracksPointing(fTracksRecCuts, jettracklist1, jet1, GetFFRadius(), sumPt1);
+               GetJetTracksPointing(fTracksRecCuts, jettracklist2, jet2, GetFFRadius(), sumPt2);
+             }
+           
+           Int_t nTracks = jettracklist1->GetSize(); 
+           if (jettracklist1->GetSize() < jettracklist2->GetSize()) nTracks = jettracklist2->GetSize();
+           
+           for(Int_t it=0; it<nTracks; ++it)
+             {
+               if (it < jettracklist1->GetSize())
+                 { 
+                   Float_t trackPt1 = (dynamic_cast<AliVParticle*> (jettracklist1->At(it)))->Pt();
+                   Float_t jetPt1   = jet1->Pt();
+                   
+                   Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+                   
+                   fFFDiJetHistosRecCuts->FillDiJetFF(1, trackPt1, jetPt1, jetBin, incrementJetPt);
+                   fFFDiJetHistosRecCuts->FillDiJetFF(0, trackPt1, jetPt1, jetBin, incrementJetPt);
+                   
+                   if (it == 0)
+                     {
+                       leadTrackPt1 = trackPt1;
+                       
+                       fFFDiJetHistosRecLeadingTrack->FillDiJetFF(1, leadTrackPt1, jetPt1, jetBin, kTRUE); 
+                       fFFDiJetHistosRecLeadingTrack->FillDiJetFF(0, leadTrackPt1, jetPt1, jetBin, kTRUE); 
+                     }
+                   
+                   fFFDiJetHistosRecLeading->FillDiJetFF(1, trackPt1, leadTrackPt1, jetBin, incrementJetPt); 
+                   fFFDiJetHistosRecLeading->FillDiJetFF(0, trackPt1, leadTrackPt1, jetBin, incrementJetPt); 
+                 }
+               
+               if (it < jettracklist2->GetSize())
+                 { 
+                   Float_t trackPt2   = (dynamic_cast<AliVParticle*>(jettracklist2->At(it)))->Pt();
+                   Float_t jetPt2     = jet2->Pt();
+                   
+                   Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+                   
+                   fFFDiJetHistosRecCuts->FillDiJetFF(2, trackPt2, jetPt2, jetBin, incrementJetPt);
+                   fFFDiJetHistosRecCuts->FillDiJetFF(0, trackPt2, jetPt2, jetBin, incrementJetPt);
+                   
+                   if (it == 0)
+                     {
+                       leadTrackPt2 = trackPt2;
+                       
+                       fFFDiJetHistosRecLeadingTrack->FillDiJetFF(2, leadTrackPt2, jetPt2, jetBin, kTRUE);
+                       fFFDiJetHistosRecLeadingTrack->FillDiJetFF(0, leadTrackPt2, jetPt2, jetBin, kTRUE);
+                     }
+                   
+                   fFFDiJetHistosRecLeading->FillDiJetFF(2, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
+                   fFFDiJetHistosRecLeading->FillDiJetFF(0, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
+                 }
+             } // End loop on tracks
+
+           delete jettracklist1;
+           delete jettracklist2;
+
+         } // End if(jetBin > 0)
+       else { Printf("Jet bins for di-jet studies not set !");}
+      } // End if(go)
+  } // End if(nRecJets > 1)
+
+  if (nGenJets > 1)
+  {
+    AliAODJet* jet1 = dynamic_cast<AliAODJet*>(fJetsGen->At(0));
+    AliAODJet* jet2 = dynamic_cast<AliAODJet*>(fJetsGen->At(1));
+
+    Double_t deltaPhi = 0;
+    Double_t phi1 = TVector2::Phi_0_2pi(jet1->Phi());
+    Double_t phi2 = TVector2::Phi_0_2pi(jet2->Phi());
+    deltaPhi      = TMath::Abs(phi1-phi2); 
+    if (deltaPhi > TMath::Pi() && deltaPhi < 2*TMath::Pi()) deltaPhi = 2*TMath::Pi() - deltaPhi;
+
+    Double_t Et1            = TMath::Abs(jet1->E()*TMath::Sin(jet1->Theta()));
+    Double_t Et2            = TMath::Abs(jet2->E()*TMath::Sin(jet2->Theta()));
+    Double_t sumEt          = Et1 + Et2;
+    Double_t normEt1PlusEt2 = TMath::Sqrt(Et1*Et1+Et2*Et2+2*Et1*Et2*TMath::Cos(deltaPhi));
+    Double_t ratio          = (Double_t)(normEt1PlusEt2/sumEt);
+
+    // DiJet events selection
+    Bool_t positionCut       = 0;
+    Bool_t positionEnergyCut = 0;
+    Bool_t cdfCut            = 0; 
+
+    // Position cut :
+    if (deltaPhi > fDiJetDeltaPhiCut) positionCut = 1;
+    // Position-Energy cut :
+    if ((deltaPhi > fDiJetDeltaPhiCut) && ((jet2->Pt()) >= fDiJetPtFractionCut*(jet1->Pt()))) positionEnergyCut = 1;
+    // CDF cut :
+    if (ratio < fDiJetCDFCut) cdfCut = 1;    
+
+    Int_t go = 0;
+
+    if (fDiJetCut == 1 && positionCut == 1) go = 1;
+    if (fDiJetCut == 2 && positionEnergyCut == 1) go = 1;
+    if (fDiJetCut == 3 && cdfCut == 1) go = 1;
+
+    if (go)
+    {
+      Double_t deltaEta      = TMath::Abs(jet1->Eta()-jet2->Eta());
+      Double_t deltaPt       = TMath::Abs(jet1->Pt()-jet2->Pt());
+      Double_t meanEt        = (Double_t)((Et1+Et2)/2.);
+      Double_t invariantMass = (Double_t)InvMass(jet1,jet2);
+
+      Double_t jetBin = GetDiJetBin(invariantMass, jet1->Pt(), meanEt, fDiJetKindBins);
+
+      if(jetBin > 0)
+      {
+        fQADiJetHistosGen->FillDiJetQA(invariantMass, deltaPhi, deltaEta, deltaPt, jetBin);
+
+        TList* jettracklist1 = new TList();
+        Double_t sumPt1 = 0.;
+        Float_t leadTrackPt1 = 0.;
+
+        TList* jettracklist2 = new TList();
+        Double_t sumPt2 = 0.;
+        Float_t leadTrackPt2 = 0.;
+      
+        if(GetFFRadius()<=0)
+        {
+         GetJetTracksTrackrefs(jettracklist1, jet1);
+         GetJetTracksTrackrefs(jettracklist2, jet2);
+        }
+        else
+        {
+         GetJetTracksPointing(fTracksGen, jettracklist1, jet1, GetFFRadius(), sumPt1);
+         GetJetTracksPointing(fTracksGen, jettracklist2, jet2, GetFFRadius(), sumPt2);
+        }
+      
+        Int_t nTracks = jettracklist1->GetSize(); 
+        if (jettracklist1->GetSize() < jettracklist2->GetSize()) nTracks = jettracklist2->GetSize();
+
+        for(Int_t it=0; it<nTracks; ++it)
+        {
+          if (it < jettracklist1->GetSize())
+         { 
+           Float_t trackPt1 = (dynamic_cast<AliAODTrack*>(jettracklist1->At(it)))->Pt();
+           Float_t jetPt1 = jet1->Pt();
+
+            Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+
+           fFFDiJetHistosGen->FillDiJetFF( 1, trackPt1, jetPt1, jetBin, incrementJetPt);
+           fFFDiJetHistosGen->FillDiJetFF( 0, trackPt1, jetPt1, jetBin, incrementJetPt);
+
+           if(it==0)
+            { 
+             leadTrackPt1 = trackPt1;
+
+             fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 1, leadTrackPt1, jetPt1, jetBin, kTRUE);
+             fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 0, leadTrackPt1, jetPt1, jetBin, kTRUE);
+           }
+
+           fFFDiJetHistosGenLeading->FillDiJetFF( 1, trackPt1, leadTrackPt1, jetBin, incrementJetPt);
+           fFFDiJetHistosGenLeading->FillDiJetFF( 0, trackPt1, leadTrackPt1, jetBin, incrementJetPt);
+         }
+      
+          if (it < jettracklist2->GetSize())
+         { 
+           Float_t trackPt2 = (dynamic_cast<AliAODTrack*>(jettracklist2->At(it)))->Pt();
+           Float_t jetPt2 = jet2->Pt();
+
+            Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+
+           fFFDiJetHistosGen->FillDiJetFF( 2, trackPt2, jetPt2, jetBin, incrementJetPt);
+            fFFDiJetHistosGen->FillDiJetFF( 0, trackPt2, jetPt2, jetBin, incrementJetPt);
+       
+            if (it==0)
+            { 
+             leadTrackPt2 = trackPt2;
+
+             fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 2, leadTrackPt2, jetPt2, jetBin, kTRUE);
+             fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 0, leadTrackPt2, jetPt2, jetBin, kTRUE);
+           }
+
+           fFFDiJetHistosGenLeading->FillDiJetFF( 2, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
+            fFFDiJetHistosGenLeading->FillDiJetFF( 0, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
+         }
+       } // End loop on tracks
+
+       delete jettracklist1;
+       delete jettracklist2;
+
+      } // End if(jetBin > 0)
+      else { Printf("Jet bins for di-jet studies not set !");}
+    } // End if (go)
+  } // End if(nGenJets > 1)
+
+  
+  // ____ efficiency _______________________________
+
+  // arrays for generated particles: reconstructed AOD track index, isPrimary flag
+  TArrayI indexAODTr; 
+  TArrayS isGenPrim; 
+  
+  // array for reconcstructed AOD tracks: generated particle index
+  TArrayI indexMCTr; 
+  
+  Int_t  nTracksAODMCCharged = GetListOfTracks(fTracksAODMCCharged, kTrackAODMCCharged);
+  if(fDebug>2)Printf("%s:%d selected AODMC tracks: %d ",(char*)__FILE__,__LINE__,nTracksAODMCCharged);
+  
+  Int_t  nTracksRecQualityCuts = GetListOfTracks(fTracksRecQualityCuts, kTrackAODQualityCuts);
+  if(fDebug>2)Printf("%s:%d selected rec tracks quality after cuts, full acceptance/pt : %d ",(char*)__FILE__,__LINE__,nTracksRecQualityCuts);
+  
+  // associate gen and rec tracks, store indices in TArrays 
+  AssociateGenRec(fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,indexMCTr,isGenPrim);
+  
+  // single track eff
+  FillSingleTrackRecEffHisto(fhnSingleTrackRecEffHisto,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim);
+
+  // jet track eff
+  for(Int_t ij=0; ij<nRecEffJets; ++ij){ 
+    
+    AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecEff->At(ij));
+    
+    if(ij==0){ // leading jet
+      
+      TList* jettracklist = new TList();
+      Double_t sumPt = 0.;
+      
+      GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt); // for efficiency: gen tracks from pointing with gen/rec jet
+      
+      Double_t jetEta = jet->Eta();
+      Double_t jetPhi = TVector2::Phi_0_2pi(jet->Phi());
+      Double_t jetPt  = sumPt;
+      
+      fQAJetHistosRecEffLeading->FillJetQA( jetEta, jetPhi, jetPt );
+      FillJetTrackRecEffHisto(fhnJetTrackRecEffHisto,jetPhi,jetEta,jetPt,jettracklist,fTracksAODMCCharged,indexAODTr,isGenPrim);
+      
+      delete jettracklist;
+    }
+  }
+   
+  //___________________
+  
+  fTracksRec->Clear();
+  fTracksRecCuts->Clear();
+  fTracksGen->Clear();
+  fTracksAODMCCharged->Clear();
+  fTracksRecQualityCuts->Clear();
+
+  fJetsRec->Clear();
+  fJetsRecCuts->Clear();
+  fJetsGen->Clear();
+  fJetsRecEff->Clear();
+
+  //Post output data.
+  PostData(1, fCommonHistList);
+  
+}
+
+//________________________________________________________________________________________
+Double_t AliAnalysisTaskFragmentationFunction::InvMass(AliAODJet* jet1, AliAODJet* jet2)
+{
+
+  Double_t invMass = 0.;
+  invMass = TMath::Sqrt(pow(jet1->E()+jet2->E(),2) - pow(jet1->Px()+jet2->Px(),2) - 
+                       pow(jet1->Py()+jet2->Py(),2) - pow(jet1->Pz()+jet2->Pz(),2));
+
+  return invMass;
+
+}
+
+//________________________________________________________________________________________
+Double_t AliAnalysisTaskFragmentationFunction::GetDiJetBin(Double_t invMass, Double_t leadingJetPt, Double_t EtMean, Int_t kindBins)
+{
+  Double_t jetBinOk = 0.;
+  Double_t jetBin = 0.;
+
+  Float_t stepInvMass = (fDiJetJetInvMassMax - fDiJetJetInvMassMin)/fDiJetNBinsJetInvMass;
+  Float_t stepPt = (fDiJetJetPtMax - fDiJetJetPtMin)/fDiJetNBinsJetPt;
+
+  if (kindBins == 1)
+    {
+      for(Int_t i=0; i<fDiJetNBinsJetInvMass; ++i)
+       {
+         jetBin = fDiJetJetInvMassMin + i*stepInvMass/2.;
+         if(((fDiJetJetInvMassMin+i*stepInvMass) <= invMass) &&
+            (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
+       }
+    }
+  else if (kindBins == 3)
+    {
+      for(Int_t i=0; i<fDiJetNBinsJetPt; ++i)
+       {
+         jetBin = fDiJetJetPtMin + i*stepPt/2.;
+         if(((fDiJetJetPtMin+i*stepPt) <= EtMean) &&
+            (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
+       }
+    }
+  else if (kindBins == 2)
+    {
+      for(Int_t i=0; i<fDiJetNBinsJetPt; ++i)
+       {
+         jetBin = fDiJetJetPtMin + i*stepPt/2.;
+         if(((fDiJetJetPtMin+i*stepPt) <= leadingJetPt) &&
+            (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
+       }
+    }
+  else {Printf("WARNING: kindBins wrongly set ! Please make sure to call SetKindSlices() and set the kind parameter to 1, 2 or 3.\n");}
+
+  return jetBinOk;
+
+}
+
+
+//______________________________________________________________
+void AliAnalysisTaskFragmentationFunction::Terminate(Option_t *) 
+{
+  // terminated
+
+  if(fDebug > 1) printf("AliAnalysisTaskFragmentationFunction::Terminate() \n");
+}  
+
+//_________________________________________________________________________________
+Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t type)
+{
+  // fill list of tracks selected according to type
+
+  if(fDebug > 2) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
+  
+  if(!list){
+    if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
+    return -1;
+  }
+
+  if(type==kTrackUndef) return 0;
+  
+  Int_t iCount = 0;
+  if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD){
+
+    // all rec. tracks, esd filter mask, eta range
+    if(!fAOD) return -1;
+    
+    for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
+      AliAODTrack *tr = fAOD->GetTrack(it);
+      
+      if(type == kTrackAODCuts || type==kTrackAODQualityCuts ){
+       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
+       if(type == kTrackAODCuts){
+         if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
+         if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
+         if(tr->Pt()  < fTrackPtCut) continue;
+       }
+      }
+      list->Add(tr);
+      iCount++;
+    }
+  }
+  else if (type==kTrackKineAll || type==kTrackKineCharged || type==kTrackKineChargedAcceptance){
+    // kine particles, all or rather charged
+    if(!fMCEvent) return iCount;
+    
+    for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
+      AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
+      
+      if(type == kTrackKineCharged || type == kTrackKineChargedAcceptance){
+       if(part->Charge()==0) continue;
+       
+       if(type == kTrackKineChargedAcceptance && 
+          (       part->Eta() < fTrackEtaMin
+               || part->Eta() > fTrackEtaMax
+               || part->Phi() < fTrackPhiMin
+               || part->Phi() > fTrackPhiMax 
+               || part->Pt()  < fTrackPtCut)) continue;
+      }
+      
+      list->Add(part);
+      iCount++;
+    }
+  }
+  else if (type==kTrackAODMCCharged || type==kTrackAODMCAll || type==kTrackAODMCChargedAcceptance) {
+    // MC particles (from AOD), physical primaries, all or rather charged or rather charged within acceptance
+    if(!fAOD) return -1;
+    
+    TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
+    if(!tca)return iCount;
+    
+    for(int it=0; it<tca->GetEntriesFast(); ++it){
+      AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
+      if(!part->IsPhysicalPrimary())continue;
+      
+      if (type==kTrackAODMCCharged || type==kTrackAODMCChargedAcceptance){
+       if(part->Charge()==0) continue;
+       if(type==kTrackAODMCChargedAcceptance && 
+          (     part->Eta() > fTrackEtaMax
+             || part->Eta() < fTrackEtaMin
+             || part->Phi() > fTrackPhiMax
+             || part->Phi() < fTrackPhiMin
+             || part->Pt()  < fTrackPtCut)) continue;
+      }
+      
+      list->Add(part);
+      iCount++;
+    }
+  }
+  
+  list->Sort();
+  return iCount;
+  
+}
+// _______________________________________________________________________________
+Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t type)
+{
+  // fill list of jets selected according to type
+  
+  if(!list){
+    if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
+    return -1;
+  }
+
+  if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
+
+    if(fBranchRecJets.Length()==0){
+      Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
+      if(fDebug>1)fAOD->Print();
+      return 0;
+    }
+
+    TClonesArray *aodRecJets = new TClonesArray();
+    if(fBranchRecJets.Length()) aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRecJets.Data()));\r
+    if(!aodRecJets)             aodRecJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fBranchRecJets.Data()));\r
+
+    if(!aodRecJets){
+      if(fBranchRecJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecJets.Data());
+
+      if(fDebug>1)fAOD->Print();
+      return 0;
+    }
+
+    Int_t nRecJets = 0;
+    
+    for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
+
+      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
+      if(!tmp) continue;
+       
+      if( tmp->Pt() < fJetPtCut ) continue;
+      if( type == kJetsRecAcceptance &&
+         (    tmp->Eta() < fJetEtaMin
+           || tmp->Eta() > fJetEtaMax
+           || tmp->Phi() < fJetPhiMin
+           || tmp->Phi() > fJetPhiMax )) continue;
+      
+      list->Add(tmp);
+         
+      nRecJets++;
+    }
+
+    list->Sort();
+    return nRecJets;
+    delete aodRecJets;
+  }
+  else if(type == kJetsKine || type == kJetsKineAcceptance){
+    
+    // generated jets
+    Int_t nGenJets = 0;
+    
+    if(!fMCEvent){
+      if(fDebug>1) Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
+      return 0;
+    }
+    
+    AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
+    if(!pythiaGenHeader){
+      Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
+      return 0;
+    }
+    
+    // fetch the pythia generated jets
+    for(int ip=0; ip<pythiaGenHeader->NTriggerJets(); ++ip){
+      
+      Float_t p[4];
+      AliAODJet *jet = new AliAODJet();
+      pythiaGenHeader->TriggerJet(ip, p);
+      jet->SetPxPyPzE(p[0], p[1], p[2], p[3]);
+
+      if( type == kJetsKineAcceptance &&
+          (    jet->Eta() < fJetEtaMin
+            || jet->Eta() > fJetEtaMax
+            || jet->Phi() < fJetPhiMin
+            || jet->Phi() > fJetPhiMax )) continue;
+      
+      list->Add(jet);
+      nGenJets++;
+    }
+    list->Sort();
+    return nGenJets;
+  }
+  else if(type == kJetsGen || type == kJetsGenAcceptance ){
+
+    if(fBranchGenJets.Length()==0){
+      if(fDebug>1) Printf("%s:%d no gen jet branch specified", (char*)__FILE__,__LINE__);
+      return 0;
+    }
+    
+    TClonesArray *aodGenJets = new TClonesArray();
+    if(fBranchGenJets.Length()) aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGenJets.Data()));\r
+    if(!aodGenJets)             aodGenJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fBranchGenJets.Data()));\r
+
+    if(!aodGenJets){
+      if(fDebug>0){
+       if(fBranchGenJets.Length())         Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data());
+      }
+      if(fDebug>1)fAOD->Print();
+      return 0;
+    }
+
+    Int_t nGenJets = 0;
+    
+    for(Int_t ig=0; ig<aodGenJets->GetEntries(); ++ig){
+         
+      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
+      if(!tmp) continue;
+         
+      if( tmp->Pt() < fJetPtCut ) continue;
+      if( type == kJetsGenAcceptance &&
+         (    tmp->Eta() < fJetEtaMin
+           || tmp->Eta() > fJetEtaMax
+           || tmp->Phi() < fJetPhiMin
+           || tmp->Phi() > fJetPhiMax )) continue;
+      
+      list->Add(tmp);
+      
+      nGenJets++;
+    }
+    list->Sort();
+    return nGenJets;
+    delete aodGenJets;
+  } 
+  else{
+    if(fDebug>0)Printf("%s:%d no such type %d",(char*)__FILE__,__LINE__,type);
+    return 0;
+  }
+}
+
+// _________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::SetProperties(THnSparse* h,const Int_t dim, const char** labels)
+{
+  //Set properties of THnSparse 
+
+  for(Int_t i=0; i<dim; i++){
+
+    h->GetAxis(i)->SetTitle(labels[i]);
+    h->GetAxis(i)->SetTitleColor(1);
+  }
+}
+
+// __________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
+{
+  //Set properties of histos (x and y title)
+
+  h->SetXTitle(x);
+  h->SetYTitle(y);
+  h->GetXaxis()->SetTitleColor(1);
+  h->GetYaxis()->SetTitleColor(1);
+}
+
+// _________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::SetProperties(TH2* h,const char* x, const char* y, const char* z)
+{
+  //Set properties of histos (x,y and z title)
+
+  h->SetXTitle(x);
+  h->SetYTitle(y);
+  h->SetZTitle(z);
+  h->GetXaxis()->SetTitleColor(1);
+  h->GetYaxis()->SetTitleColor(1);
+  h->GetZaxis()->SetTitleColor(1);
+}
+
+// ________________________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist, TList* outputlist, AliAODJet* jet, const Double_t radius,Double_t& sumPt)
+{
+  // fill list of tracks in cone around jet axis  
+
+  sumPt = 0;
+
+  Double_t jetMom[3];
+  jet->PxPyPz(jetMom);
+  TVector3 jet3mom(jetMom);
+
+  for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
+
+    AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
+
+    Double_t trackMom[3];
+    track->PxPyPz(trackMom);
+    TVector3 track3mom(trackMom);
+
+    Double_t dR = jet3mom.DeltaR(track3mom);
+
+    if(dR<radius){
+
+      outputlist->Add(track);
+      
+      sumPt += track->Pt();
+    }
+  }
+  
+  outputlist->Sort();
+}
+
+// ___________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, AliAODJet* jet)
+{
+  // list of jet tracks from trackrefs
+  
+  Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
+
+  for (Int_t itrack=0; itrack<nTracks; itrack++) {
+    
+    AliVParticle* track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
+    if(!track){
+      AliError("expected ref track not found ");
+      continue;
+    }
+        
+    list->Add(track);
+  }
+  
+  list->Sort();
+}
+
+// _ ________________________________________________________________________________________________________________________________
+void  AliAnalysisTaskFragmentationFunction::AssociateGenRec(TList* tracksAODMCCharged,TList* tracksRec, TArrayI& indexAODTr,TArrayI& indexMCTr,TArrayS& isGenPrim)
+{
+  // associate generated and reconstructed tracks, fill TArrays of list indices
+
+
+  Int_t nTracksRec  = tracksRec->GetSize();
+  Int_t nTracksGen  = tracksAODMCCharged->GetSize();
+  TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
+
+  if(!nTracksGen) return;
+  if(!tca)        return;
+  
+  // set size
+  indexAODTr.Set(nTracksGen);
+  indexMCTr.Set(nTracksRec);
+  isGenPrim.Set(nTracksGen);
+
+  indexAODTr.Reset(-1);
+  indexMCTr.Reset(-1);
+  isGenPrim.Reset(0);
+
+  // loop over reconstructed tracks, get generated track 
+
+  for(Int_t iRec=0; iRec<nTracksRec; iRec++){ 
+      
+    AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec)); 
+
+    Int_t label = rectrack->GetLabel();
+
+    // find MC track in our list
+    AliAODMCParticle* gentrack = 0x0;
+    if(label>=0) gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(label));
+
+    Int_t listIndex = -1;
+    if(gentrack) listIndex = tracksAODMCCharged->IndexOf(gentrack);
+
+    if(listIndex>=0){
+
+      indexAODTr[listIndex] = iRec;
+      indexMCTr[iRec]       = listIndex;
+    }
+  } 
+
+
+  // define primary sample for reconstruction efficiency
+
+  for(Int_t iGen=0; iGen<nTracksGen; iGen++){
+
+    AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksAODMCCharged->At(iGen));
+
+    Int_t pdg = gentrack->GetPdgCode();    
+
+    // 211 - pi, 2212 - proton, 321 - Kaon, 11 - electron, 13 - muon
+    if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || 
+       TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
+      
+      isGenPrim[iGen] = kTRUE;
+    }
+  }
+}
+
+// _____________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::FillSingleTrackRecEffHisto(THnSparse* histo, TList* tracksGen, TList* tracksRec, TArrayI& indexAODTr, TArrayS& isGenPrim){
+
+  // fill THnSparse for single track reconstruction efficiency
+
+  Int_t nTracksGen  = tracksGen->GetSize();
+
+  if(!nTracksGen) return;
+
+  for(Int_t iGen=0; iGen<nTracksGen; iGen++){
+
+    if(isGenPrim[iGen] != 1) continue; // select primaries
+
+    AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksGen->At(iGen));
+    
+    Double_t ptGen  = gentrack->Pt();
+    Double_t etaGen = gentrack->Eta();
+    Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
+
+    // apply same acc & pt cuts as for FF 
+    // could in principle also be done setting THNsparse axis limits before projecting, 
+    // but then the binning needs to be fine grained enough 
+
+    if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
+    if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
+    if(ptGen  < fTrackPtCut) continue;
+
+    Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track 
+    Double_t isRec =  0;
+    Double_t ptRec = -1;
+
+    if(iRec>=0){
+
+      AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec)); 
+      ptRec = rectrack->Pt();
+      isRec = 1;
+    }
+
+    Double_t entries[5] = {phiGen,etaGen,ptGen,ptRec,isRec};
+    histo->Fill(entries);
+  }
+}
+
+// ______________________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::FillJetTrackRecEffHisto(THnSparse* histo,Double_t jetPhi, Double_t jetEta, Double_t jetPt, TList* jetTrackList, 
+                                                                  TList* tracksGen, TArrayI& indexAODTr, TArrayS& isGenPrim)
+{
+  // fill THnSparse for jet track reconstruction efficiency
+
+  Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
+
+  if(!nTracksJet) return; 
+
+  for(Int_t iTr=0; iTr<nTracksJet; iTr++){
+
+    AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (jetTrackList->At(iTr));
+
+    // find jet track in gen tracks list
+    Int_t iGen = tracksGen->IndexOf(gentrack); 
+
+    if(iGen<0){
+      if(fDebug>0) Printf("%s:%d gen jet track not found ",(char*)__FILE__,__LINE__);
+      continue;
+    }
+
+    if(isGenPrim[iGen] != 1) continue; // select primaries
+    
+    Double_t ptGen  = gentrack->Pt();
+    Double_t etaGen = gentrack->Eta();
+    Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
+
+    // apply same acc & pt cuts as for FF 
+    // could in principle also be done setting THNsparse axis limits before projecting, 
+    // but then the binning needs to be fine grained enough 
+
+    if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
+    if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
+    if(ptGen  < fTrackPtCut) continue;
+
+    Double_t z = ptGen / jetPt;
+    Double_t xi = 0;
+    if(z>0) xi = TMath::Log(1/z);
+
+    Double_t isRec =  0;
+    Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track 
+    if(iRec>=0) isRec = 1;
+
+    Double_t entries[7] = {jetPhi,jetEta,jetPt,ptGen,z,xi,isRec};
+    histo->Fill(entries);
+  }
+}