]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJetChem.cxx
adding usage of vtx axis as a parameter in the addtask macro
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetChem.cxx
index cecb209ef4c13b5103562ababa85673a2caf400c..7161b2f4e480af964ad213de4d3f51a08f2be3c4 100644 (file)
@@ -1,3 +1,11 @@
+/*************************************************************************
+ *                                                                       *
+ *                                                                       *
+ *      Task for Jet Chemistry 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.                  *
  **************************************************************************/
 
-/* $Id:$ */
-
+/* $Id: */
 
-#include <TChain.h>
-#include <TFile.h>
-#include <TList.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TProfile.h>
-#include <TVector3.h>
-#include <TLorentzVector.h>
-#include <TMath.h>
-#include <TTree.h>
-#include <TDatabasePDG.h> 
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TProfile.h"
+#include "THnSparse.h"
 
-#include "AliAnalysisTaskJetChem.h"
+#include "AliAnalysisHelperJetTasks.h"
 #include "AliAnalysisManager.h"
-#include "AliMCEventHandler.h"
-#include "AliAODInputHandler.h"
-#include "AliAODHandler.h"
+#include "AliAODHandler.h" 
+#include "AliAODInputHandler.h" 
+#include "AliESDEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+
+#include "AliAODEvent.h"
 #include "AliAODJet.h"
+#include "AliAODv0.h"
 #include "AliAODTrack.h"
-#include "AliAODMCParticle.h"
 
-#include "AliGenPythiaEventHeader.h"
-#include "AliAnalysisHelperJetTasks.h"
-#include "AliESDtrack.h"
+
+#include "AliPID.h" 
 #include "AliExternalTrackParam.h"
 
-//#include <iostream> // TEST!!!
+#include "AliAnalysisTaskJetChem.h"
 
+ClassImp(AliAnalysisTaskJetChem)
 
-//
-// Analysis class for jet chemistry studies
-// based on AliAnalysisTaskUE by Arian Abrahantes Quintana and Enesto Lopez
-// contact: Oliver Busch, o.busch@gsi.de
+//____________________________________________________________________________
+AliAnalysisTaskJetChem::AliAnalysisTaskJetChem()
+   : AliAnalysisTaskFragmentationFunction()
+   ,fK0Type(0)
+   ,fFilterMaskK0(0)
+   ,fListK0s(0)
+   ,fV0QAK0(0)
+   ,fFFHistosRecCutsK0Evt(0)      
+   ,fFFHistosIMK0AllEvt(0)        
+   ,fFFHistosIMK0Jet(0)           
+   ,fFFHistosIMK0Cone(0)
+   ,fFFHistosPhiCorrIMK0(0)
+   ,fFFIMNBinsJetPt(0)    
+   ,fFFIMJetPtMin(0) 
+   ,fFFIMJetPtMax(0)
+   ,fFFIMNBinsInvM(0) 
+   ,fFFIMInvMMin(0)   
+   ,fFFIMInvMMax(0)   
+   ,fFFIMNBinsPt(0)      
+   ,fFFIMPtMin(0)        
+   ,fFFIMPtMax(0)        
+   ,fFFIMNBinsXi(0)      
+   ,fFFIMXiMin(0)        
+   ,fFFIMXiMax(0)        
+   ,fFFIMNBinsZ(0)       
+   ,fFFIMZMin(0)         
+   ,fFFIMZMax(0)
+   ,fPhiCorrIMNBinsPt(0)
+   ,fPhiCorrIMPtMin(0)
+   ,fPhiCorrIMPtMax(0)
+   ,fPhiCorrIMNBinsPhi(0)
+   ,fPhiCorrIMPhiMin(0)
+   ,fPhiCorrIMPhiMax(0)
+   ,fPhiCorrIMNBinsInvM(0)
+   ,fPhiCorrIMInvMMin(0)
+   ,fPhiCorrIMInvMMax(0)
+   ,fh1K0Mult(0)
+   ,fh1dPhiJetK0(0)
+{
+   // default constructor
+}
 
-// This class needs as input AOD with track and Jets
-// the output is a list of histograms
-//
-// AOD can be either connected to the InputEventHandler  
-// for a chain of AOD files 
-// or 
-// to the OutputEventHandler
-// for a chain of ESD files, so this case class should be 
-// in the train after the Jet finder
-//
+//__________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name) 
+  : AliAnalysisTaskFragmentationFunction(name)
+  ,fK0Type(0)  
+  ,fFilterMaskK0(0)
+  ,fListK0s(0)
+  ,fV0QAK0(0)
+  ,fFFHistosRecCutsK0Evt(0)      
+  ,fFFHistosIMK0AllEvt(0)        
+  ,fFFHistosIMK0Jet(0)           
+  ,fFFHistosIMK0Cone(0)
+  ,fFFHistosPhiCorrIMK0(0)
+  ,fFFIMNBinsJetPt(0)    
+  ,fFFIMJetPtMin(0) 
+  ,fFFIMJetPtMax(0)
+  ,fFFIMNBinsInvM(0) 
+  ,fFFIMInvMMin(0)   
+  ,fFFIMInvMMax(0)   
+  ,fFFIMNBinsPt(0)      
+  ,fFFIMPtMin(0)        
+  ,fFFIMPtMax(0)        
+  ,fFFIMNBinsXi(0)      
+  ,fFFIMXiMin(0)        
+  ,fFFIMXiMax(0)        
+  ,fFFIMNBinsZ(0)       
+  ,fFFIMZMin(0)         
+  ,fFFIMZMax(0)
+  ,fPhiCorrIMNBinsPt(0)
+  ,fPhiCorrIMPtMin(0)
+  ,fPhiCorrIMPtMax(0)
+  ,fPhiCorrIMNBinsPhi(0)
+  ,fPhiCorrIMPhiMin(0)
+  ,fPhiCorrIMPhiMax(0)
+  ,fPhiCorrIMNBinsInvM(0)
+  ,fPhiCorrIMInvMMin(0)
+  ,fPhiCorrIMInvMMax(0)
+  ,fh1K0Mult(0)
+  ,fh1dPhiJetK0(0)     
+{
+  // constructor
+  
+  DefineOutput(1,TList::Class());  
+}
 
+//__________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const  AliAnalysisTaskJetChem &copy)
+  : AliAnalysisTaskFragmentationFunction()
+  ,fK0Type(copy.fK0Type)
+  ,fFilterMaskK0(copy.fFilterMaskK0)
+  ,fListK0s(copy.fListK0s)
+  ,fV0QAK0(copy.fV0QAK0)
+  ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)      
+  ,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)        
+  ,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)           
+  ,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)          
+  ,fFFHistosPhiCorrIMK0(copy.fFFHistosPhiCorrIMK0)          
+  ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)   
+  ,fFFIMJetPtMin(copy.fFFIMJetPtMin)     
+  ,fFFIMJetPtMax(copy.fFFIMJetPtMax)     
+  ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)  
+  ,fFFIMInvMMin(copy.fFFIMInvMMin)    
+  ,fFFIMInvMMax(copy.fFFIMInvMMax)    
+  ,fFFIMNBinsPt(copy.fFFIMNBinsPt)      
+  ,fFFIMPtMin(copy.fFFIMPtMin)        
+  ,fFFIMPtMax(copy.fFFIMPtMax)        
+  ,fFFIMNBinsXi(copy.fFFIMNBinsXi)      
+  ,fFFIMXiMin(copy.fFFIMXiMin)        
+  ,fFFIMXiMax(copy.fFFIMXiMax)        
+  ,fFFIMNBinsZ(copy.fFFIMNBinsZ)       
+  ,fFFIMZMin(copy.fFFIMZMin)         
+  ,fFFIMZMax(copy.fFFIMZMax) 
+  ,fPhiCorrIMNBinsPt(copy.fPhiCorrIMNBinsPt)
+  ,fPhiCorrIMPtMin(copy.fPhiCorrIMPtMin)
+  ,fPhiCorrIMPtMax(copy.fPhiCorrIMPtMax)
+  ,fPhiCorrIMNBinsPhi(copy.fPhiCorrIMNBinsPhi)
+  ,fPhiCorrIMPhiMin(copy.fPhiCorrIMPhiMin)
+  ,fPhiCorrIMPhiMax(copy.fPhiCorrIMPhiMax)
+  ,fPhiCorrIMNBinsInvM(copy.fPhiCorrIMNBinsInvM)
+  ,fPhiCorrIMInvMMin(copy.fPhiCorrIMInvMMin)
+  ,fPhiCorrIMInvMMax(copy.fPhiCorrIMInvMMax)
+  ,fh1K0Mult(copy.fh1K0Mult)
+  ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
+{
+  // copy constructor
+  
+}
 
+// _________________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    AliAnalysisTaskFragmentationFunction::operator=(o);
+
+    fK0Type                         = o.fK0Type;
+    fFilterMaskK0                   = o.fFilterMaskK0;
+    fListK0s                        = o.fListK0s;
+    fV0QAK0                         = o.fV0QAK0;
+    fFFHistosRecCutsK0Evt           = o.fFFHistosRecCutsK0Evt;      
+    fFFHistosIMK0AllEvt             = o.fFFHistosIMK0AllEvt;        
+    fFFHistosIMK0Jet                = o.fFFHistosIMK0Jet;           
+    fFFHistosIMK0Cone               = o.fFFHistosIMK0Cone;          
+    fFFHistosPhiCorrIMK0            = o.fFFHistosPhiCorrIMK0;
+    fFFIMNBinsJetPt                 = o.fFFIMNBinsJetPt;    
+    fFFIMJetPtMin                   = o.fFFIMJetPtMin; 
+    fFFIMJetPtMax                   = o.fFFIMJetPtMax;
+    fFFIMNBinsPt                    = o.fFFIMNBinsPt;      
+    fFFIMPtMin                      = o.fFFIMPtMin;        
+    fFFIMPtMax                      = o.fFFIMPtMax;        
+    fFFIMNBinsXi                    = o.fFFIMNBinsXi;      
+    fFFIMXiMin                      = o.fFFIMXiMin;        
+    fFFIMXiMax                      = o.fFFIMXiMax;        
+    fFFIMNBinsZ                     = o.fFFIMNBinsZ;       
+    fFFIMZMin                       = o.fFFIMZMin;         
+    fFFIMZMax                       = o.fFFIMZMax;         
+    fPhiCorrIMNBinsPt               = o.fPhiCorrIMNBinsPt;
+    fPhiCorrIMPtMin                 = o.fPhiCorrIMPtMin;
+    fPhiCorrIMPtMax                 = o.fPhiCorrIMPtMax;
+    fPhiCorrIMNBinsPhi              = o.fPhiCorrIMNBinsPhi;
+    fPhiCorrIMPhiMin                = o.fPhiCorrIMPhiMin;
+    fPhiCorrIMPhiMax                = o.fPhiCorrIMPhiMax;
+    fPhiCorrIMNBinsInvM             = o.fPhiCorrIMNBinsInvM;
+    fPhiCorrIMInvMMin               = o.fPhiCorrIMInvMMin;
+    fPhiCorrIMInvMMax               = o.fPhiCorrIMInvMMax;
+    fh1K0Mult                       = o.fh1K0Mult;
+    fh1dPhiJetK0                    = o.fh1dPhiJetK0;
+  }
+    
+  return *this;
+}
 
+//_______________________________________________
+AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
+{
+  // destructor  
 
-ClassImp(AliAnalysisTaskJetChem)
+  
+  if(fTracksRecCuts) delete fTracksRecCuts;
+  if(fJetsRecCuts) delete fJetsRecCuts;
+  if(fListK0s) delete fListK0s;
+  if(fCommonHistList) delete fCommonHistList;
 
-////////////////////////////////////////////////////////////////////////
-
-
-//____________________________________________________________________
-AliAnalysisTaskJetChem:: AliAnalysisTaskJetChem(const char* name): AliAnalysisTaskSE(name),
-fDebug(kFALSE),
-fDeltaAOD(kFALSE),
-fDeltaAODBranch(""),
-fAODBranch("jets"),
-fDeltaAODBranchMC(""),
-fAODBranchMC("jetsMC"),
-fArrayJetsAOD(0x0), 
-fArrayJetsMC(0x0),
-fAOD(0x0),            
-fAODjets(0x0),
-fListOfHistos(0x0),
-fJetsOnFly(kFALSE),
-fUseLOConeJets(kFALSE),
-fUseLOConeMCJets(kFALSE),
-fUsePythiaJets(kFALSE),
-fConeRadius(0.),
-fTrackPtCutJF(0),
-fFilterBitJF(0x01),
-fRequireITSRefitJF(kFALSE),
-fRejectK0TracksJF(kFALSE),
-fJetPtCut(5.0),
-fJetEtaCut(0.5),
-fFilterBit(0x10),
-fTrackPtCut(0.),
-fTrackEtaCut(0.9),
-fUseOnFlyV0s(kFALSE),
-fCutnSigdEdx(0.), 
-fUseAODMCTracksForUE(kFALSE),
-fAreaReg(1.0), 
-fAvgTrials(1),
-fhPrimVertexNCont(0x0),
-fhPrimVertexRho(0x0),
-fhPrimVertexZ(0x0),
-fhNJets(0x0),
-fhNJetsMC(0x0),
-fhLeadingEta(0x0),
-fhLeadingNTracksVsEta(0x0),
-fhLeadingPtVsEta(0x0),
-fhLeadingPhi(0x0), 
-fhLeadingPt(0x0),
-fhLeadingPtDiffr(0x0),
-fhLeadingEtaMC(0x0), 
-fhLeadingPhiMC(0x0), 
-fhLeadingPtMC(0x0),
-fhLeadingPtMCDiffr(0x0),
-fhPhiEtaTracksNoCut(0x0),
-fhPtTracksNoCut(0x0),
-fhPhiEtaTracks(0x0),
-fhPtTracks(0x0),
-fhTrackMult(0x0),
-fhEtaMCTracks(0x0),            
-fhPhiMCTracks(0x0),            
-fhPtMCTracks(0x0),
-fhnTracksVsPtLeading(0x0),
-fhdNdEtaPhiDist(0x0),        
-fhRegionSumPtMaxVsEt(0x0),
-fhRegionMultMaxVsEt(0x0),     
-fhRegionSumPtMinVsEt(0x0),
-fhRegionMultMinVsEt(0x0),     
-fhNV0s(0x0),
-fhV0onFly(0x0),
-fhV0DCADaughters(0x0),
-fhV0Radius(0x0),
-fhV0DCAToVertex(0x0),
-fhV0DCAToVertexK0(0x0),
-fhV0InvMassK0(0x0),
-fhV0PtVsInvMassK0(0x0),
-fhV0InvMassK0JetEvt(0x0),
-fhV0InvMassLambda(0x0),
-fhV0InvMassAntiLambda(0x0),
-fhV0InvMassLambdaJetEvt(0x0),
-fhV0InvMassAntiLambdaJetEvt(0x0),
-fhdROpanK0VsPt(0x0),
-fhdPhiJetV0(0x0),
-fhdPhiJetK0(0x0),
-fhdRJetK0(0x0),
-fhdNdptV0(0x0),
-fhdNdptK0(0x0),
-fhPtVsEtaK0(0x0),
-fhV0InvMassK0DCA(0x0),
-fhV0InvMassK0DCAdEdx(0x0),
-fhV0PtVsInvMassK0DCAdEdx(0x0),
-fhV0InvMassK0DCAPID(0x0),
-fhV0InvMassLambdaDCAdEdx(0x0), 
-fhV0InvMassAntiLambdaDCAdEdx(0x0),
-fhdNdptK0DCA(0x0),
-fhdNdptK0DCAdEdx(0x0),
-fhV0InvMassK0Min(0x0),
-fhV0InvMassLambdaMin(0x0),
-fhV0InvMassAntiLambdaMin(0x0),
-fhV0InvMassK0Max(0x0),
-fhV0InvMassLambdaMax(0x0),
-fhV0InvMassAntiLambdaMax(0x0),
-fhV0InvMassK0Jet(0x0),
-fhV0InvMassLambdaJet(0x0),
-fhV0InvMassAntiLambdaJet(0x0),
-fhV0InvMassK0Lambda(0x0), 
-fhdNdptK0JetEvt(0x0),
-fhdNdzK0(0x0),  
-fhdNdzK05to10(0x0),  
-fhdNdzK010to20(0x0), 
-fhdNdzK020to30(0x0), 
-fhdNdzK030to40(0x0), 
-fhdNdzK040to60(0x0), 
-fhdNdxiK0(0x0),
-fhdNdzLambda(0x0),
-fhdNdzAntiLambda(0x0),
-fhdNdzK0Max(0x0),     
-fhdNdxiK0Max(0x0),    
-fhdNdzLambdaMax(0x0), 
-fhdNdxiLambdaMax(0x0),
-fhdNdptK0Max(0x0), 
-fhdNdptLambdaMax(0x0),
-fhdNdzK0Min(0x0),     
-fhdNdxiK0Min(0x0),    
-fhdNdzLambdaMin(0x0), 
-fhdNdxiLambdaMin(0x0),
-fhdNdptK0Min(0x0),
-fhdNdptLambdaMin(0x0),
-fhdNdzK0Jet(0x0),     
-fhdNdxiK0Jet(0x0),    
-fhdNdzLambdaJet(0x0), 
-fhdNdxiLambdaJet(0x0),
-fhdNdptK0Jet(0x0),
-fhdNdptLambdaJet(0x0),
-fhdEdxVsMomV0(0x0),
-fhdEdxVsMomV0pidEdx(0x0),
-fhdEdxVsMomV0piPID(0x0),
-fhdPhiJetK0MC(0x0),
-fhdRJetK0MC(0x0),
-fhdRV0MC(0x0),
-fhdNdptchPiMCMax(0x0),    
-fhdNdptK0MCMax(0x0),      
-fhdNdptchKMCMax(0x0),     
-fhdNdptpMCMax(0x0),       
-fhdNdptpBarMCMax(0x0),    
-fhdNdptLambdaMCMax(0x0),  
-fhdNdptLambdaBarMCMax(0x0),
-fhdNdptchPiMCMin(0x0),    
-fhdNdptK0MCMin(0x0),      
-fhdNdptchKMCMin(0x0),     
-fhdNdptpMCMin(0x0),       
-fhdNdptpBarMCMin(0x0),    
-fhdNdptLambdaMCMin(0x0),  
-fhdNdptLambdaBarMCMin(0x0),
-fhdNdptOmegaMCMin(0x0),   
-fhdNdptOmegaBarMCMin(0x0),
-fhdNdptchPiMCJet(0x0),    
-fhdNdptK0MCJet(0x0),      
-fhdNdptchKMCJet(0x0),     
-fhdNdptpMCJet(0x0),       
-fhdNdptpBarMCJet(0x0),    
-fhdNdptLambdaMCJet(0x0),  
-fhdNdptLambdaBarMCJet(0x0),
-fhPIDMC(0x0),
-fhPIDMC_quarkEv(0x0),
-fhPIDMC_gluonEv(0x0),
-fhPIDMCAll(0x0),
-fhPIDMCMin(0x0),
-fhPIDMCJet(0x0),
-fhPIDMCMotherK0(0x0),
-fhPIDMCGrandMotherK0(0x0),
-fhPIDMCMotherChK(0x0),
-fhPIDMCMotherK0Trans(0x0),
-fhPIDMCGrandMotherK0Trans(0x0),
-fhPIDMCMotherChKTrans(0x0),
-fhdNdptgammaMC(0x0),
-fhdNdptchPiMC(0x0),
-fhdNdptpi0MC(0x0),
-fhdNdptK0MC(0x0),
-fhdNdptchKMC(0x0),
-fhdNdptpMC(0x0),
-fhdNdptpBarMC(0x0),
-fhdNdptLambdaMC(0x0),
-fhdNdptLambdaBarMC(0x0),
-fhdNdptOmegaMC(0x0),
-fhdNdptOmegaBarMC(0x0),
-fhdNdxiMC(0x0),
-fhdNdxiK0MC(0x0),
-fhdNdxiK0MCJet(0x0),
-fhdNdzK0MC(0x0), 
-fhdNdzK0MCJet(0x0), 
-fhdNdptK0MCJetEvt(0x0),
-fhnJetsAODvsMC(0x0),
-fhLeadingPtAODvsMC(0x0),
-fhLeadingEtaAODvsMC(0x0),
-fhLeadingPhiAODvsMC(0x0),
-fhnTracksLeadingAODvsMC(0x0),
-fhLeadingdRAODMC(0x0),
-fhLeadingPtAODvsMCdRcut(0x0),
-fhdnTracksVsdPtLeadingAODMC(0x0),
-fhnTracksJetVsPtAOD(0x0),
-fhnTracksJetVsPtAODquarkEv(0x0),
-fhRadiusJetVsPtAOD(0x0), 
-fhnTracksJetVsPtMC(0x0),
-fhnTracksJetVsPtMCquarkEv(0x0),
-fhRadiusJetVsPtMC(0x0),
-fhnTracksJetVsPtMCK0(0x0),
-fhnTracksJetVsPtMCK0quarkEv(0x0),
-fhRadiusJetVsPtMCK0(0x0),
-fhnTracksJetVsPtAODK0(0x0),
-fhnTracksJetVsPtAODK0quarkEv(0x0),
-fhRadiusJetVsPtAODK0(0x0),
-fhnTracksJetVsPtAODpKch(0x0),
-fhRadiusJetVsPtAODpKch(0x0),
-fhPythiaProcess(0x0), 
-fhPythiaProcessK0(0x0),
-fhPythiaProcessKch(0x0),
-fhPythiaProcessp(0x0),
-fhPythiaProcesspbar(0x0),
-fhdNdzJets5to10(0x0),   
-fhdNdzJets10to20(0x0),  
-fhdNdzJets20to30(0x0),  
-fhdNdzJets30to40(0x0),  
-fhdNdzJets40to60(0x0),  
-fhdNdxiJets5to10(0x0),  
-fhdNdxiJets10to20(0x0),
-fhdNdxiJets20to30(0x0), 
-fhdNdxiJets30to40(0x0), 
-fhdNdxiJets40to60(0x0), 
-fhdNdptTracksJetPt5to10(0x0),
-fhdNdptTracksJetPt10to20(0x0),
-fhdNdptTracksJetPt20to30(0x0),
-fhdNdptTracksJetPt30to40(0x0),
-fhdNdptTracksJetPt40to60(0x0),
-fh1Xsec(0x0),
-fh1Trials(0x0),
-fpdgdb(0x0){
-  // Default constructor 
-
-  fAreaReg = 2*TMath::Pi()/6.0 * 2*fTrackEtaCut;
-  fpdgdb = TDatabasePDG::Instance(); 
-
-  // Output slot #1 writes into a TList container, 0 reserved for std AOD output
-  DefineOutput(1, TList::Class());
 }
 
-//______________________________________________________________
-Bool_t AliAnalysisTaskJetChem::UserNotify()
+//________________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name, 
+                                                                          Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
+                                                                          Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
+                                                                          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)
+  ,fNBinsInvMass(nInvMass)
+  ,fInvMassMin(invMassMin)  
+  ,fInvMassMax(invMassMax)
+  ,fNBinsPt(nPt) 
+  ,fPtMin(ptMin)   
+  ,fPtMax(ptMax)   
+  ,fNBinsXi(nXi) 
+  ,fXiMin(xiMin)   
+  ,fXiMax(xiMax)   
+  ,fNBinsZ(nZ)  
+  ,fZMin(zMin)    
+  ,fZMax(zMax)    
+  ,fh3TrackPt(0)
+  ,fh3Xi(0)
+  ,fh3Z(0)
+  ,fh1JetPt(0)
+  ,fNameFF(name)
 {
-  //
-  // read the cross sections
-  // and number of trials from pyxsec.root
-  // 
-
-  fAvgTrials = 1;
-  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
-  Float_t xsection = 0;
-  Float_t ftrials  = 1;
-  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;
-    }
-    AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
-    fh1Xsec->Fill("<#sigma>",xsection);
-    // construct average trials 
-    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
-    if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; 
-  }  
-  return kTRUE;
+  // default constructor
+
 }
 
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::UserCreateOutputObjects()
+//______________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
+  : TObject()
+  ,fNBinsJetPt(copy.fNBinsJetPt)
+  ,fJetPtMin(copy.fJetPtMin)
+  ,fJetPtMax(copy.fJetPtMax)
+  ,fNBinsInvMass(copy.fNBinsInvMass)
+  ,fInvMassMin(copy.fInvMassMin)  
+  ,fInvMassMax(copy.fInvMassMax)
+  ,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)    
+  ,fh3TrackPt(copy.fh3TrackPt)
+  ,fh3Xi(copy.fh3Xi)
+  ,fh3Z(copy.fh3Z)
+  ,fh1JetPt(copy.fh1JetPt)
+  ,fNameFF(copy.fNameFF)
 {
-  // Create the output container
-  //
-  AliInfo("UserCreateOutPutObjects()");
-  //
-  //  Histograms
-
-  OpenFile(1);
-  CreateHistos();
-  PostData(1, fListOfHistos); // PostData at least once
-
+  // copy constructor
 }
 
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::UserExec(Option_t */*option*/)
+//______________________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
 {
-  // Execute analysis for current event
-  //
-  if (fDebug > 3) AliInfo( " Processing event..." );
-
-  TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
-  
-  if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
-    fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
-    if (fDebug > 1) AliInfo("  ==== Tracks from AliAODInputHandler");
-    // Case when jets are reconstructed on the fly from AOD tracks
-    // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
-    // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
-    // different parameters to default ones stored in the AOD or to use a different algorithm
-    if( fJetsOnFly ) {
-      handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
-      if( handler && handler->InheritsFrom("AliAODHandler") ) {
-        fAODjets = ((AliAODHandler*)handler)->GetAOD();
-        if (fDebug > 1) AliInfo("  ==== Jets from AliAODHandler");
-      }
-    } else {
-      fAODjets = fAOD;
-      if (fDebug > 1) AliInfo("  ==== Jets from AliAODInputHandler");
-    }
-  } else {
-    handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
-    if( handler && handler->InheritsFrom("AliAODHandler") ) {
-      fAOD  = ((AliAODHandler*)handler)->GetAOD();
-      fAODjets = fAOD;
-      if (fDebug > 1) AliInfo("  ==== Tracks and Jets from AliAODHandler");
-    } else {
-      AliFatal("I can't get any AOD Event Handler");
-      return;
-    }
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsJetPt   = o.fNBinsJetPt;
+    fJetPtMin     = o.fJetPtMin;
+    fJetPtMax     = o.fJetPtMax;
+    fNBinsInvMass = o.fNBinsInvMass;
+    fInvMassMin   = o.fInvMassMin;  
+    fInvMassMax   = o.fInvMassMax;
+    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;    
+    fh3TrackPt    = o.fh3TrackPt;
+    fh3Xi         = o.fh3Xi;
+    fh3Z          = o.fh3Z;
+    fh1JetPt      = o.fh1JetPt;
+    fNameFF       = o.fNameFF;
   }
+    
+  return *this;
+}
 
-  // -------------
-
-  // fetch the pythia header info and get the trials
-  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  Float_t nTrials = 1;
-  if (mcHandler) {  
-    AliMCEvent* mcEvent = mcHandler->MCEvent();
-    if (mcEvent) {
-      AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
-      if(pythiaGenHeader){
-       nTrials = pythiaGenHeader->Trials();
-      }
-    }
-  }
+//___________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
+{ 
+  // destructor 
 
-  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
-    
-  AnalyseEvent();
-  
-  // Post the data
-  PostData(1, fListOfHistos);
+  if(fh1JetPt)   delete fh1JetPt;
+  if(fh3TrackPt) delete fh3TrackPt;
+  if(fh3Xi)      delete fh3Xi;
+  if(fh3Z)       delete fh3Z;
 }
 
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::AnalyseEvent()
+//_________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
 {
-  // check Trigger 
+  // book FF histos
 
-  // TString firedTriggerClasses = fAOD->GetFiredTriggerClasses(); 
-  // AliInfo(Form("firedTriggerClasses %s",firedTriggerClasses.Data()));
-  // if(firedTriggerClasses.Length() > 0 && !firedTriggerClasses.Contains("CINT1B")) return;  
+  fh1JetPt   = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
+  fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
+  fh3Xi      = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
+  fh3Z       = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
 
-  AliAODVertex *primVertex = fAOD->GetPrimaryVertex();  // is this from tracks or from SPD ? SPD should only be used if no track vertex
-  if(!primVertex){
-    AliInfo("no prim Vertex found - skip event");
-    fhPrimVertexNCont->Fill(-1);
-    return; 
-  }
-    
-  TString primVtxName(primVertex->GetName());
-  
-  if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
+  AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries"); 
+  AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
+  AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
+  AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
+}
 
-    AliInfo("found TPC prim Vertex  - skip event");
-    fhPrimVertexNCont->Fill(-1);
-    return; 
-  }
-    
-  Int_t nVertContributors = primVertex->GetNContributors();
-  fhPrimVertexNCont->Fill(nVertContributors);
-  
-  //cout<<" prim vertex name "<<primVertex->GetName()<<" nCont "<<nVertContributors<<endl;
-  
-  if(nVertContributors<1){ // eventually check if not SPD vertex ??? 
-    AliInfo("prim vertex no contributors - skip event");
-    return;
-  }
-    
-  Double_t vertX = primVertex->GetX();
-  Double_t vertY = primVertex->GetY();
-  Double_t vertZ = primVertex->GetZ();
-  
-  Double_t vertRho = TMath::Sqrt(vertX*vertX+vertY*vertY);
-  
-  fhPrimVertexRho->Fill(vertRho);
-  fhPrimVertexZ->Fill(vertZ);
-  
-  if(TMath::Abs(vertZ)>10){
-    AliInfo(Form("prim vertex z=%f - skip event",vertZ));
-    return; 
-  }
-  
-  
-  Int_t pythiaPID = GetPythiaProcessID();
-  
-  // ------------------------------------------------
-  // Find Leading Jets 1,2,3 
-  // (could be skipped if Jets are sort by Pt...)
-
-  //Int_t    index1 = -1;
-  //Int_t nTracksLeading = 0;
-
-  Int_t nJetsAOD = 0;
-  AliAODJet* leadingJetAOD = NULL; // non-zero if any jet in acc and leading jet survives pt cut 
-  Int_t indexLeadingAOD   = -1;
-  Double_t ptLeadingAOD   = 0.; 
-  Int_t indexMaxRegionAOD = 0; // initialize with '0', '+-1' transverse regions 
-
-  Int_t nJetsMC = 0;
-  AliAODJet* leadingJetMC = NULL; // non-zero if leading jet survives acc, pt cut 
-  Int_t indexLeadingMC   = -1;
-  Double_t ptLeadingMC   = 0.; 
-  Int_t indexMaxRegionMC = 0; // initialize with '0', '+-1' transverse regions 
-
-  Double_t   ptLeadingAODAllEta  = 0; 
-  AliAODJet* leadingJetAODAllEta = NULL; 
+//________________________________________________________________________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
+{
+  // fill FF
+  if(incrementJetPt) fh1JetPt->Fill(jetPt);    
+  fh3TrackPt->Fill(jetPt,invM,trackPt);
   
-  if(fUseLOConeJets)   fArrayJetsAOD  = FindChargedParticleJets(); 
-  else{ // Use jets on AOD/DeltaAOD
-
-    if(fDeltaAOD){
-      if (fDebug > 1) AliInfo(" ==== Jets From  Delta-AODs !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fDeltaAODBranch.Data()));
-      fArrayJetsAOD = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data());
-      if (!fArrayJetsAOD){
-       AliFatal(" No jet-array! ");
-       return;
-      }
-    }
-    else{
-      if (fDebug > 1) AliInfo(" ==== AOD jets: Read Standard-AODs  !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fAODBranch.Data()));
-      
-      fArrayJetsAOD = (TClonesArray*)fAODjets->FindListObject(fAODBranch.Data());
-    }
-  }
+  Double_t z = 0.;
+  if(jetPt>0) z = trackPt / jetPt;
+  Double_t xi = 0;
+  if(z>0) xi = TMath::Log(1/z);
   
+  fh3Xi->Fill(jetPt,invM,xi);
+  fh3Z->Fill(jetPt,invM,z);
+}
 
-  if(fUseLOConeMCJets)    fArrayJetsMC = FindChargedParticleJetsMC(); 
-  else if(fUsePythiaJets) fArrayJetsMC = GetPythiaJets();
-  else{
-
-    if(fDeltaAOD){
-      if (fDebug > 1) AliInfo(" ==== MC Jets From  Delta-AODs !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fDeltaAODBranchMC.Data()));
-      fArrayJetsMC = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranchMC.Data());
-      if (!fArrayJetsMC){
-       AliFatal(" No jet-array! ");
-       return;
-      }
-    }
-    else{
-      if (fDebug > 1) AliInfo(" ==== MC jets: Read Standard-AODs  !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fAODBranchMC.Data()));
-      
-      fArrayJetsMC = (TClonesArray*)fAODjets->FindListObject(fAODBranchMC.Data());
-    }
-  }
-
-  if(fArrayJetsAOD) nJetsAOD = fArrayJetsAOD->GetEntries();
-  if(fArrayJetsMC)  nJetsMC  = fArrayJetsMC->GetEntries();
-
-  fhNJets->Fill(nJetsAOD);
-  fhNJetsMC->Fill(nJetsMC);
-  fhnJetsAODvsMC->Fill(nJetsMC,nJetsAOD);
+//___________________________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
+{
+  // add histos to list
 
-  if(fDebug>1) AliInfo(Form("AOD %d jets",nJetsAOD)); 
+  list->Add(fh1JetPt);
+  
+  list->Add(fh3TrackPt);
+  list->Add(fh3Xi);
+  list->Add(fh3Z);
+}
 
+// ---
+
+
+//_______________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const char* name,
+                                                                                        Int_t nPt,   Float_t ptMin,   Float_t ptMax,
+                                                                                        Int_t nPhi,  Float_t phiMin,  Float_t phiMax,
+                                                                                        Int_t nInvMass,  Float_t invMassMin,  Float_t invMassMax)
+  : TObject()
+  ,fNBinsPt(nPt)
+  ,fPtMin(ptMin)
+  ,fPtMax(ptMax)
+  ,fNBinsPhi(nPhi)
+  ,fPhiMin(phiMin)
+  ,fPhiMax(phiMax)
+  ,fNBinsInvMass(nInvMass)
+  ,fInvMassMin(invMassMin)  
+  ,fInvMassMax(invMassMax)
+  ,fh3PhiCorr(0) 
+  ,fNamePhiCorr(name) 
+{
+  // default constructor
+}
 
-  // for xcheck: find leading jet in large eta range
-  for(Int_t i=0; i<nJetsAOD; i++){
+//____________________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const AliFragFuncHistosPhiCorrInvMass& copy)
+  : TObject()
+  ,fNBinsPt(copy.fNBinsPt)
+  ,fPtMin(copy.fPtMin)
+  ,fPtMax(copy.fPtMax)
+  ,fNBinsPhi(copy.fNBinsPhi)
+  ,fPhiMin(copy.fPhiMin)
+  ,fPhiMax(copy.fPhiMax)
+  ,fNBinsInvMass(copy.fNBinsInvMass)
+  ,fInvMassMin(copy.fInvMassMin)  
+  ,fInvMassMax(copy.fInvMassMax)
+  ,fh3PhiCorr(copy.fh3PhiCorr)
+  ,fNamePhiCorr(copy.fNamePhiCorr)
+{
+  // copy constructor
+}
 
-    AliAODJet* jet  = (AliAODJet*) fArrayJetsAOD->At(i);
-    Double_t jetPt  = jet->Pt();
+//________________________________________________________________________________________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& o)
+{
+  // assignment
+  
+  if(this!=&o){
+    TObject::operator=(o);
+    fNBinsPt      = o.fNBinsPt;
+    fPtMin        = o.fPtMin;
+    fPtMax        = o.fPtMax;
+    fNBinsPhi     = o.fNBinsPhi;
+    fPhiMin       = o.fPhiMin;
+    fPhiMax       = o.fPhiMax;
+    fNBinsInvMass = o.fNBinsInvMass;
+    fInvMassMin   = o.fInvMassMin;  
+    fInvMassMax   = o.fInvMassMax;
     
-    if(jetPt > ptLeadingAODAllEta){ 
-      ptLeadingAODAllEta  = jetPt;
-      leadingJetAODAllEta = jet;
-    }
+    fh3PhiCorr    = o.fh3PhiCorr;
+    fNamePhiCorr  = o.fNamePhiCorr;
   }
+  
+  return *this;
+}
 
-  if(leadingJetAODAllEta){
-    //cout<<" trackRefs entries "<<leadingJetAODAllEta->GetRefTracks()->GetEntriesFast()<<endl;
-    fhLeadingNTracksVsEta->Fill(leadingJetAODAllEta->Eta(),leadingJetAODAllEta->GetRefTracks()->GetEntriesFast());
-    fhLeadingPtVsEta->Fill(leadingJetAODAllEta->Eta(),leadingJetAODAllEta->Pt()); 
-    
-  }
+//_________________________________________________________________________________________
+AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::~AliFragFuncHistosPhiCorrInvMass()
+{
+  // destructor 
+  
+  if(fh3PhiCorr) delete fh3PhiCorr;
+}
+
+//__________________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::DefineHistos()
+{
+  // book jet QA histos
 
-  // find leading jet AOD
-  for(Int_t i=0; i<nJetsAOD; ++i){
+  fh3PhiCorr  = new TH3F(Form("fh3PhiCorrIM%s", fNamePhiCorr.Data()), 
+                        Form("%s: p_{t} - #phi - m_{inv} distribution",fNamePhiCorr.Data()), 
+                        fNBinsPt, fPtMin, fPtMax, 
+                        fNBinsPhi, fPhiMin, fPhiMax,
+                        fNBinsInvMass, fInvMassMin, fInvMassMax);
+  
+  AliAnalysisTaskJetChem::SetProperties(fh3PhiCorr, "p_{t} (GeV/c)", "#phi", "m_{inv} (GeV/c^2)"); 
+}
 
-    AliAODJet* jet  = (AliAODJet*) fArrayJetsAOD->At(i); 
-    Double_t jetPt  = jet->Pt();
-    Double_t jetEta = jet->Eta();
+//___________________________________________________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::FillPhiCorr(Float_t pt, Float_t phi, Float_t invM)
+{
+  // fill jet QA histos 
 
-    if((jetPt > ptLeadingAOD) && (TMath::Abs(jetEta)<fJetEtaCut)){ 
-      ptLeadingAOD    = jetPt; 
-      indexLeadingAOD = i;
-    }
-  }
+  fh3PhiCorr->Fill(pt, phi, invM);
+}
 
-  // find leading jet MC
-  for(Int_t i=0; i<nJetsMC; ++i){
+//______________________________________________________________________________________________
+void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AddToOutput(TList* list) const 
+{
+  // add histos to list
 
-    AliAODJet* jet  = (AliAODJet*) fArrayJetsMC->At(i); 
-    Double_t jetPt  = jet->Pt();
-    Double_t jetEta = jet->Eta();
+  list->Add(fh3PhiCorr);
+}
 
-    if((jetPt > ptLeadingMC) && (TMath::Abs(jetEta)<fJetEtaCut)){
-      ptLeadingMC    = jetPt; 
-      indexLeadingMC = i;
-    }
-  }
+//____________________________________________________
+void AliAnalysisTaskJetChem::UserCreateOutputObjects()
+{
+  // create output objects
 
+  if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
  
-  //cout<<" ptLeadingAOD "<<ptLeadingAOD<<" MC "<<ptLeadingMC<<endl;
-
-  if(indexLeadingAOD>=0){ // event with jet  
+  // create list of tracks and jets 
+  
+  fTracksRecCuts = new TList();
+  fTracksRecCuts->SetOwner(kFALSE);  
 
-    Double_t etaLeadingAOD  = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->Eta();
-    Double_t phiLeadingAOD  = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->Phi();
-    Int_t nTracksLeadingAOD = ((AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD))->GetRefTracks()->GetEntriesFast();
-    
-    if(fDebug>1) AliInfo(Form("\n Pt Leading AOD Jet = %6.1f eta=%5.3f, nTracks %d",ptLeadingAOD,etaLeadingAOD,nTracksLeadingAOD));
+  fJetsRecCuts = new TList();
+  fJetsRecCuts->SetOwner(kFALSE);
 
-    fhLeadingEta->Fill(etaLeadingAOD);
-    if(TMath::Abs(etaLeadingAOD)<fJetEtaCut){ // leading jet eta cut
+  fListK0s = new TList(); 
+  fListK0s->SetOwner(kFALSE);
 
-      fhnTracksVsPtLeading->Fill(ptLeadingAOD,nTracksLeadingAOD);
-      fhLeadingPt->Fill(ptLeadingAOD);
-      if(IsDiffractiveEvent(pythiaPID))  fhLeadingPtDiffr->Fill(ptLeadingAOD);
+  //
+  // Create histograms / output container
+  //
 
-      if(ptLeadingAOD>fJetPtCut){ // leading jet pt cut
-       
-       leadingJetAOD = (AliAODJet*) fArrayJetsAOD->At(indexLeadingAOD);
-       
-       fhLeadingPhi->Fill(phiLeadingAOD);
-       
-       // ----------------------------------------------
-       // Find max and min regions
-       Double_t sumPtRegionPosit  = 0.;
-       Double_t sumPtRegionNegat  = 0.;
-       Int_t    nTrackRegionPosit = 0;
-       Int_t    nTrackRegionNegat = 0;
-      
-       Int_t nTracks = fAOD->GetNTracks();
-         
-       for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-           
-         AliAODTrack* part = fAOD->GetTrack(ipart);
-         if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
-         if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point
-         // PID Selection: Reject everything but hadrons
-         Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || 
-           part->GetMostProbablePID()==AliAODTrack::kKaon || 
-           part->GetMostProbablePID()==AliAODTrack::kProton;
-         if (!isHadron ) continue;
-         if (!part->Charge() ) continue; //Only charged
-         if (part->Pt() < fTrackPtCut ) continue;
-         if(TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
-         
-         TVector3 partVect(part->Px(), part->Py(), part->Pz());
-           
-         Int_t region = IsTrackInsideRegion(leadingJetAOD,&partVect );  
-
-         if (region > 0) {
-           sumPtRegionPosit += part->Pt();
-           nTrackRegionPosit++;
-         }
-         if (region < 0) {
-           sumPtRegionNegat += part->Pt();
-           nTrackRegionNegat++;
-         }
-       } // tracks loop 
-
-       // fill sumPt and mult density for transverse regions
-       
-       if( sumPtRegionPosit > sumPtRegionNegat ) {
-         FillSumPtRegion( ptLeadingAOD, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
-       } 
-       else {
-         FillSumPtRegion( ptLeadingAOD, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
-       }
-       if (nTrackRegionPosit > nTrackRegionNegat ) {
-         FillMultRegion( ptLeadingAOD, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg);
-       } 
-       else {
-         FillMultRegion( ptLeadingAOD, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg);
-       }
-        
-       indexMaxRegionAOD = (sumPtRegionPosit > sumPtRegionNegat) ? 1 : -1;
-       
-      } // leading jet pt cut
-    } // leading jet eta cut  
-  } // jet event
+  OpenFile(1);
+  fCommonHistList = new TList();
   
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+       
 
-  if(indexLeadingMC>=0){ // event with MC jet  
+  // histograms inherited from AliAnalysisTaskFragmentationFunction
 
-    Double_t etaLeadingMC  = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->Eta();
-    Double_t phiLeadingMC  = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->Phi();
-    Int_t nTracksLeadingMC = ((AliAODJet*) fArrayJetsMC->At(indexLeadingMC))->GetRefTracks()->GetEntriesFast();
-    
-    if(fDebug>1) AliInfo(Form("\n Pt Leading MC Jet = %6.1f eta=%5.3f, nTracks %d",ptLeadingMC,etaLeadingMC,nTracksLeadingMC));
+  fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
+  fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
 
-    fhLeadingEtaMC->Fill(etaLeadingMC);
-  
-    if(TMath::Abs(etaLeadingMC)<fJetEtaCut){ // leading jet eta cut
+  fh1EvtCent                = new TH1F("fh1EvtCent","centrality",100,0.,100.);
+  fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
+  fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
+  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);
 
-      fhLeadingPtMC->Fill(ptLeadingMC);
-      if(IsDiffractiveEvent(pythiaPID)) fhLeadingPtMCDiffr->Fill(ptLeadingMC);
+  fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
 
-      if(ptLeadingMC>fJetPtCut){ // leading jet pt cut
-       
-       leadingJetMC = (AliAODJet*) fArrayJetsMC->At(indexLeadingMC);
-   
-       fhLeadingPhiMC->Fill(phiLeadingMC); // -pi to pi
 
-       // ----------------------------------------------
-       // Find max and min regions
-       Double_t sumPtRegionPosit  = 0;
-       Double_t sumPtRegionNegat  = 0;
-       Int_t    nTrackRegionPosit = 0;
-       Int_t    nTrackRegionNegat = 0;
+  // histograms jetChem proper 
 
-       TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-       
-       Int_t ntrks = farray->GetEntries();
-       if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
-         
-       for(Int_t i =0 ; i < ntrks; i++){   
-           
-         AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
-         //Cuts
-         if (!(mctrk->IsPhysicalPrimary())) continue;
-         //if (!(mctrk->IsPrimary())) continue;
-         
-         if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue;
-           
-         if (mctrk->Pt() < fTrackPtCut ) continue;
-         if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue;
-         
-         Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
-           TMath::Abs(mctrk->GetPdgCode())==2212 ||
-           TMath::Abs(mctrk->GetPdgCode())==321;
-         
-         if (!isHadron) continue;
-           
-         TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
-         
-         Int_t region = IsTrackInsideRegion(leadingJetMC,&partVect );  
-         
-         if (region > 0) {
-           sumPtRegionPosit += mctrk->Pt();
-           nTrackRegionPosit++;
-         }
-         if (region < 0) {
-           sumPtRegionNegat += mctrk->Pt();
-           nTrackRegionNegat++;
-         }
-       } // AliAODMCParticle loop
-       
-       indexMaxRegionMC = (sumPtRegionPosit > sumPtRegionNegat) ? 1 : -1;
-       
-      } // leading jet pt
-    } // leading jet eta 
-  } // found jet
+  fh1EvtMult                = new TH1F("fh1EvtMult","multiplicity",1200,0.,12000.);
+  fh1K0Mult                 = new TH1F("fh1K0Mult","K0 multiplicity",500,0.,500.);
+  fh1dPhiJetK0               = new TH1F("fh1dPhiJetK0","",640,-1,5.4);
 
-  Bool_t foundK0AOD   = kFALSE; // event with K0 ? 
-  Bool_t foundK0MC    = kFALSE; // event with K0 ? 
-
-  CheckV0s(leadingJetAOD,indexMaxRegionAOD,foundK0AOD); // here leadingJetAOD/MC nonzero if jet passes eta & pt cut
-  CheckMCParticles(leadingJetMC,indexMaxRegionMC,foundK0MC);
-  CompLeadingJets(leadingJetAOD,leadingJetMC,pythiaPID,foundK0AOD,foundK0MC);
-  FillReferencePlotsTracks();
-  FillReferenceFF(leadingJetAOD);
+
+
+  fFFHistosRecCuts          = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  
+  fV0QAK0                   = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
+                                                           fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
+                                                           fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
+                                                           fQATrackHighPtThreshold);
+  
+  fFFHistosRecCutsK0Evt      = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
+                                                    fFFNBinsPt, fFFPtMin, fFFPtMax, 
+                                                    fFFNBinsXi, fFFXiMin, fFFXiMax,  
+                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+  
+  
+  fFFHistosIMK0AllEvt        = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
+                                                           fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
+                                                           fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
+                                                           fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
+                                                           fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
+  
+  fFFHistosIMK0Jet           = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
+                                                           fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
+                                                           fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
+                                                           fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
+                                                           fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
   
+  
+  fFFHistosIMK0Cone          = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
+                                                           fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
+                                                           fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
+                                                           fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
+                                                           fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
+  
+  fFFHistosPhiCorrIMK0       = new AliFragFuncHistosPhiCorrInvMass("K0",fPhiCorrIMNBinsPt, fPhiCorrIMPtMin, fPhiCorrIMPtMax, 
+                                                                  fPhiCorrIMNBinsPhi, fPhiCorrIMPhiMin, fPhiCorrIMPhiMax,  
+                                                                  fPhiCorrIMNBinsInvM , fPhiCorrIMInvMMin , fPhiCorrIMInvMMax);
+  
+  
+  fV0QAK0->DefineHistos();
+  fFFHistosRecCuts->DefineHistos();
+  fFFHistosRecCutsK0Evt->DefineHistos();
+  fFFHistosIMK0AllEvt->DefineHistos();
+  fFFHistosIMK0Jet->DefineHistos();
+  fFFHistosIMK0Cone->DefineHistos();
+  fFFHistosPhiCorrIMK0->DefineHistos();
 
-  if(fUseLOConeJets && fArrayJetsAOD){
-    fArrayJetsAOD->Delete(); // no 'Clear': AliAODjet contains TMomentum and TRefArray
-    delete fArrayJetsAOD;
+  const Int_t saveLevel = 5;
+  if(saveLevel>0){
+    
+    fCommonHistList->Add(fh1EvtSelection);
+    fCommonHistList->Add(fh1EvtCent);
+    fCommonHistList->Add(fh1VertexNContributors);
+    fCommonHistList->Add(fh1VertexZ);
+    fCommonHistList->Add(fh1Xsec);
+    fCommonHistList->Add(fh1Trials);
+    fCommonHistList->Add(fh1PtHard);
+    fCommonHistList->Add(fh1PtHardTrials);
+    fCommonHistList->Add(fh1nRecJetsCuts);
+    fCommonHistList->Add(fh1EvtMult);
+    fCommonHistList->Add(fh1K0Mult);
+    fCommonHistList->Add(fh1dPhiJetK0);
+
+    fV0QAK0->AddToOutput(fCommonHistList);
+    fFFHistosRecCuts->AddToOutput(fCommonHistList);
+    fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
+
+    fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
+    fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
+    fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
+    fFFHistosPhiCorrIMK0->AddToOutput(fCommonHistList);
   }
-  if(fUseLOConeMCJets && fArrayJetsMC){
-    fArrayJetsMC->Delete(); // no 'Clear': AliAODjet contains TMomentum and TRefArray
-    delete fArrayJetsMC;
+  
+  // =========== 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 AliAnalysisTaskJetChem::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
 
+  if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
 
-Double_t AliAnalysisTaskJetChem::GetJetRadius(const AliAODJet* jet, const Double_t energyFrac){
+  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;
+    }
+  }
 
-  // calc jet radius containing fraction energyFrac of full jet pt  
+  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(fUseAODInputJets) fAODJets = fAOD;
+    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();
+      fAODJets = fAOD;
+      if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
+    }
+  }
+  
+  if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
+    TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
+    if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
+      fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
+      if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
+    }
+  }
+  
+  if(fNonStdFile.Length()!=0){
+    // case we have an AOD extension - fetch the jets from the extended output
+    
+    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+    fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
+    if(!fAODExtension){
+      if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
+    }
+  }
+  
+  if(!fAOD){
+    Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
+    return;
+  }
+  if(!fAODJets){
+    Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
+    return;
+  }
   
-  const Int_t kArraySize = 1000;
-  const Int_t kInitVal   = -999;
-
-  Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
 
-  if(nTracks>kArraySize){
-    AliError(Form("nTracks in jet %d exceeds max array size",nTracks));
-    return -1;
+  // event selection  *****************************************
+  
+  Double_t centPercent = -1;
+  if(fEventClass>0){
+    Int_t cl = 0;
+    if(handler && handler->InheritsFrom("AliAODInputHandler")){ 
+      // since it is not supported by the helper task define own classes
+      centPercent = fAOD->GetHeader()->GetCentrality();
+      cl = 1;
+      if(centPercent>10) cl = 2;
+      if(centPercent>30) cl = 3;
+      if(centPercent>50) cl = 4;
+    }
+    else {
+      cl = AliAnalysisHelperJetTasks::EventClass();
+      if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // OB added
+    }
+    
+    if(cl!=fEventClass){
+      // event not in selected event class, reject event
+      if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
+      fh1EvtSelection->Fill(2.);
+      PostData(1, fCommonHistList);
+      return;
+    }
+  }
+  
+  // *** vertex cut ***
+  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(3.);
+    PostData(1, fCommonHistList);
+    return;
+  }
+  
+  fh1VertexZ->Fill(primVtx->GetZ());
+  
+  if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
+    if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
+    fh1EvtSelection->Fill(4.);
+    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(5.);
+    PostData(1, fCommonHistList);
+    return;
   }
+  
+  if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
+  fh1EvtSelection->Fill(0.);
+  fh1EvtCent->Fill(centPercent);
+    
 
+  //___ get MC information __________________________________________________________________
 
-  Double_t deltaR[kArraySize];
-  Double_t pt[kArraySize];
-  Int_t index[kArraySize];
-  for(int i=0; i<kArraySize; i++) index[i] = kInitVal;
+  Double_t ptHard = 0.;
+  Double_t nTrials = 1; // trials for MC trigger weight for real data
   
+  if(fMCEvent){
+     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
+     AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
+     AliGenHijingEventHeader*  hijingGenHeader = 0x0;
 
-  Double_t ptTot = 0;
-
-  for(int i=0; i<nTracks; i++){
+     if(pythiaGenHeader){
+        if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
+        nTrials = pythiaGenHeader->Trials();
+        ptHard  = pythiaGenHeader->GetPtHard();
 
-    AliAODTrack* track = (AliAODTrack*) jet->GetRefTracks()->At(i);
+        fh1PtHard->Fill(ptHard);
+        fh1PtHardTrials->Fill(ptHard,nTrials);
 
-    TLorentzVector *mom4Jet   = jet->MomentumVector();
-    TVector3 mom3Jet = mom4Jet->Vect();
 
-    Double_t trackMom[3];
-    track->PxPyPz(trackMom);
-    TVector3 mom3Track(trackMom);
+     } else { // no pythia, hijing?
 
-    Double_t dR = mom3Jet.DeltaR(mom3Track);
+        if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
 
-    deltaR[i] = dR;
-    pt[i]     = track->Pt();
+         hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
+         if(!hijingGenHeader){
+            Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
+         } else {
+           if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
+        }
+     }
 
-    ptTot += pt[i];
+     fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
   }
 
-  //cout<<" ptTot "<<ptTot<<" jetPt "<<jet->Pt()<<endl; // Xcheck 
-
-  TMath::Sort(nTracks,deltaR,index,kFALSE); // sort in decreasing order 
-
-  Double_t ptSum = 0;
 
-  for(int i=0; i<nTracks; i++){
+  //____ fetch jets ______________________________________________________________
 
-    Int_t ind = index[i];
-    
-    Double_t ptTrack = pt[ind];
-    Double_t dR      = deltaR[ind];
-
-    ptSum += ptTrack;
-
-    if(ptSum >= ptTot*energyFrac) return dR;
-  }
+  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);
+  
 
-  return -1;
+  //____ fetch particles __________________________________________________________
+  Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
+  if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
+  if(fTracksRecCuts->GetEntries() != nTCuts) 
+    Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
+  fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
 
-}
+  Int_t nK0s = GetListOfK0s(fListK0s,fK0Type);
+  if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
+  if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
+  fh1K0Mult->Fill(fListK0s->GetEntries());
 
-//____________________________________________________________________
-void AliAnalysisTaskJetChem::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
-{
-  // Fill sumPt of control regions
+   // ___ V0 QA + K0 pt spectra all events _______________________________________________
   
-  fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
-  fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
-}
+  if(fListK0s->GetEntries()>0){
+    
+    for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s 
+      
+      AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
+      if(!v0) continue;
+    
+      Float_t trackPt       = v0->Pt();
+      Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+      Double_t invM         = v0->MassK0Short();
+      Double_t jetPt        = fFFIMJetPtMin; // assign pro forma jet energy
 
-//____________________________________________________________________
-void AliAnalysisTaskJetChem::FillMultRegion(Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin)
-{
-  // Fill Nch multiplicity of control regions
-  
-  fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
-  fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
+      fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt()); 
+      fFFHistosIMK0AllEvt->FillFF(trackPt, invM, jetPt, incrementJetPt);
+    }
+  }
 
-}
 
-//____________________________________________________________________
-Int_t AliAnalysisTaskJetChem::IsTrackInsideRegion(const AliAODJet* aodjetVect,const TVector3 *partVect) 
-{  
-  // return de region in delta phi
-  // -1 negative delta phi 
-  //  1 positive delta phi
-  //  0 outside region
+  //____ fill FF histos  __________________________________________________________
 
-  TLorentzVector* jetVectLorentz = aodjetVect->MomentumVector(); 
-  TVector3 jetVect = jetVectLorentz->Vect();
+  for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
 
-  static const Double_t k60rad  = 60.*TMath::Pi()/180.;
-  static const Double_t k120rad = 120.*TMath::Pi()/180.;
-  
-  Int_t region = 0;
-  if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
-  // transverse regions
-  if (jetVect.DeltaPhi(*partVect) < -k60rad && jetVect.DeltaPhi(*partVect) > -k120rad ) region = -1;
-  if (jetVect.DeltaPhi(*partVect) > k60rad && jetVect.DeltaPhi(*partVect) < k120rad ) region = 1;
+    AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
     
-  return region;
-}
+    if(ij==0){ // leading jet
+      
+      TList* jettracklist = new TList();
+      Double_t sumPt      = 0.;
+      
+      if(GetFFRadius()<=0){
+       GetJetTracksTrackrefs(jettracklist, jet);
+      } else {
+       GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
+      }
+      
+      for(Int_t it=0; it<jettracklist->GetSize(); ++it){
 
-//____________________________________________________________________
-TClonesArray* AliAnalysisTaskJetChem::FindChargedParticleJetsMC()
-{
-  // CDF jet finder: 
-  // loop over pt-ordered list of tracks, combine tracks within jet cone, recalc cone axis after each step
-  // based on implementation by Arian Abrahantes Quintana and Enesto Lopez
-  TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-  if(!farray){
-    AliInfo("no mcparticles branch"); 
-    return 0;
+       AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));      
+       if(!trackVP)continue;
+
+       Float_t jetPt   = jet->Pt();
+       Float_t trackPt = trackVP->Pt();
+       
+       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+       
+       fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
+       if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);
+      }
+      
+      delete jettracklist;
+    }
   }
 
-  Int_t nTracks = farray->GetEntries();
 
-  if( !nTracks ) return 0;
-  TObjArray tracks(nTracks);
+  for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // jet loop
 
-  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
+    AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
+    if(!jet) continue;
+    
+    if(ij==0){ // leading jet
+      
+      //fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
+      
+      Float_t jetPt   = jet->Pt();
 
-    AliAODMCParticle* part = (AliAODMCParticle*)farray->At(ipart);
+      for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s 
 
-    if(!part->IsPhysicalPrimary()) continue;
+       AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
+       if(!v0) continue;
+       Double_t v0Mom[3];
+       v0->PxPyPz(v0Mom);
+       TVector3 v0MomVect(v0Mom);
 
-    // exclude neutrinos
-    Int_t   pdg  = TMath::Abs(part->GetPdgCode());
-    if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
+       Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
+       
+       Float_t trackPt       = v0->Pt();
+       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
+       Double_t invM         = v0->MassK0Short();
+       
+       fFFHistosIMK0Jet->FillFF(trackPt, invM, jetPt, incrementJetPt);
+       fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invM);
 
-    if( !part->Charge() ) continue; // comment / uncomment here
-    fhEtaMCTracks->Fill(part->Eta());
+       if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
+       fh1dPhiJetK0->Fill(dPhiJetK0);
 
-    if(TMath::Abs(part->Eta()) > fTrackEtaCut) continue; 
-    fhPtMCTracks->Fill(part->Pt());
+      }
+      if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum 
 
-    if( part->Pt() < fTrackPtCutJF ) continue;
-    fhPhiMCTracks->Fill(part->Phi());
+       Bool_t incrementJetPt = kTRUE;
+       fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
+      }
+      
 
-    tracks.AddLast(part);
-  
-  }
-
-  QSortTracks( tracks, 0, tracks.GetEntriesFast() );
-  
-  nTracks = tracks.GetEntriesFast();
-
-  if( !nTracks ) return 0;
-  TObjArray *jets = new TObjArray(nTracks);
-  TIter itrack(&tracks);
-  while( nTracks ) {
-    // 2- Start with the highest pT particle ...
-    Float_t px,py,pz,pt; 
-
-    AliAODMCParticle* track = (AliAODMCParticle*)itrack.Next();
-
-    if( !track ) continue;
-    px = track->Px();
-    py = track->Py();
-    pz = track->Pz();
-    pt = track->Pt(); 
-    jets->AddLast( new AliAODJet(px,py,pz,pt) ); // Use the energy member to store Pt
-    AliAODJet* jet = (AliAODJet*)jets->Last();
-    jet->AddTrack(track);
-    tracks.Remove(track);
-   
-
-    TVector2 jetVect2(jet->Px(),jet->Py());
-
-    // 3- Go to the next highest pT particle not already included...
-    AliAODMCParticle* track1;
-    while ( (track1  = (AliAODMCParticle*)(itrack.Next())) ) {
-      TVector2 trackVect2(track1->Px(),track1->Py());
-      Double_t dphi = trackVect2.DeltaPhi(jetVect2);
-
-      Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
-                               dphi*dphi );
-
-      if( r < fConeRadius ) {
-        Double_t fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
-        // recalculating the centroid
-        Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
-       jetVect2.SetMagPhi(jet->E()/fPt,jetVect2.Phi());
-       trackVect2.SetMagPhi(track1->Pt()/fPt,trackVect2.Phi());
-       
-       TVector2 sumVect2 = jetVect2+trackVect2;
-       Double_t phi = sumVect2.Phi();
-
-        //jet->SetPtEtaPhiE( 1., eta, phi, fPt );
-       ((TLorentzVector*) jet->MomentumVector())->SetPtEtaPhiE(1,eta,phi,fPt);
-
-       jet->AddTrack(track1);
-        tracks.Remove(track1);
-      }
-    }
-    
-    tracks.Compress();
-
-    nTracks = tracks.GetEntries();
-    //   4- Continue until all particles are in a jet.
-    itrack.Reset();
-  } // end while nTracks
-  
-  // Convert to AODjets....
-  Int_t njets = jets->GetEntriesFast();
-
-  TClonesArray* aodjets = new TClonesArray("AliAODJet",njets);
-  aodjets->SetOwner(kTRUE);
-
-  Int_t count = 0;
-  for(Int_t ijet=0; ijet<njets; ++ijet) {
-    AliAODJet* jet = (AliAODJet*)jets->At(ijet);
-
-    if (jet->E() < fJetPtCut) continue;
-    Float_t px, py,pz,en; // convert to 4-vector
-    px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
-    py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
-    pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
-    en = TMath::Sqrt(px * px + py * py + pz * pz);
-    jet->SetPxPyPzE(px,py,pz,en);
-
-
-    TClonesArray &tmpaodjets = *aodjets;
-    new(tmpaodjets[count++]) AliAODJet(*jet);
-    //aodjets->AddLast( new AliAODJet(*jet));    
-    //aodjets->AddLast( new AliAODJet(px, py, pz, en) );
-  }
-  // Order jets according to their pT .
-  QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
-  
-  // debug
-  //if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
-  
-  jets->Delete(); // OB - should I cleanup or leave it to garbage collection ? 
-  delete jets;
-
-  return aodjets;
-}
-// ___________________________________________________________________
-
-Bool_t AliAnalysisTaskJetChem::IsTrackFromK0(const Int_t indexTrack){
-
-  // check wether track with index is from K0 (= V0 with proper inv mass)
-
-  TClonesArray* tracks = fAOD->GetTracks();
-
-  for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
-    
-    AliAODv0* v0 = fAOD->GetV0(i);
-    
-    Bool_t isOnFly = v0->GetOnFlyStatus();
-
-    if( (fUseOnFlyV0s && !isOnFly) || (!fUseOnFlyV0s && isOnFly) ) continue;
-
-    Double_t massK0 = v0->MassK0Short();
-    
-    // FIXME: here should go cuts (at least V0 quality)
-
-    if(IsK0InvMass(massK0)){
+      TList* jetConeK0list = new TList();
+      Double_t sumPt   = 0.;
       
-      AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
-      AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));   
-
-      Int_t indexPos = tracks->IndexOf(trackPos);
-      Int_t indexNeg = tracks->IndexOf(trackNeg);
-     
-
-      if(indexPos == indexTrack){ return kTRUE;}
-      if(indexNeg == indexTrack){ return kTRUE;}
-
-    }
-  }
-
-  return kFALSE;
-
-}
-
-//____________________________________________________________________
-TClonesArray*  AliAnalysisTaskJetChem::FindChargedParticleJets()
-{
-
-  // CDF jet finder: 
-  // loop over pt-ordered list of tracks, combine tracks within jet cone, recalc cone axis after each step 
-  // based on implementation by Arian Abrahantes Quintana and Enesto Lopez
-  // ref: PHYSICAL REVIEW D 65 092002, CDF Collaboration
+      GetJetTracksPointing(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPt);
  
-
-  //  1 - Order all charged particles according to their pT .
-  Int_t nTracks = fAOD->GetNTracks();
-
-  if( !nTracks ) return 0;
-  TObjArray tracks(nTracks);
-
-  // standard cuts + ITS refit = stdrd PWG4 cut (compose them for productions with old ESDFilter task) 
-
-  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-    AliAODTrack* part = fAOD->GetTrack( ipart );
-
-    UInt_t status = part->GetStatus();
-
-    if( !part->TestFilterBit(fFilterBitJF) ) continue; // track cut selection
-    if(fRequireITSRefitJF &&  ((status&AliESDtrack::kITSrefit)==0)) continue;  
-
-    if(!part->Charge() ) continue;
-
-    if(TMath::Abs(part->Eta()) > fTrackEtaCut) continue; 
-    
-    if(fRejectK0TracksJF && IsTrackFromK0(ipart)) continue; 
-
-    if( part->Pt() < fTrackPtCutJF ) continue;
-
-    tracks.AddLast(part);
-  }
-
-  QSortTracks( tracks, 0, tracks.GetEntriesFast() );
-  
-  nTracks = tracks.GetEntriesFast();
-  
-  if( !nTracks ) return 0;
-  TObjArray *jets = new TObjArray(nTracks);
-  TIter itrack(&tracks);
-  while( nTracks ) {
-    // 2- Start with the highest pT particle ...
-    Float_t px,py,pz,pt; 
-    AliAODTrack* track = (AliAODTrack*)itrack.Next();
-    if( !track ) continue;
-    px = track->Px();
-    py = track->Py();
-    pz = track->Pz();
-    pt = track->Pt(); 
-    //jets->AddLast( new TLorentzVector(px, py, pz, pt) ); // Use the energy member to store Pt
-    jets->AddLast( new AliAODJet(px,py,pz,pt) ); // Use the energy member to store Pt
-   //TLorentzVector* jet = (TLorentzVector*)jets->Last();
-    AliAODJet* jet = (AliAODJet*)jets->Last();
-    jet->AddTrack(track);
-    tracks.Remove( track );
-
-    TVector2 jetVect2(jet->Px(),jet->Py()); // OB
-
-    // 3- Go to the next highest pT particle not already included...
-    AliAODTrack* track1;
-    while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
-      //Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi()); // OB remove
-      TVector2 trackVect2(track1->Px(),track1->Py());
-      Double_t dphi = trackVect2.DeltaPhi(jetVect2);
-
-      Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
-                               dphi*dphi );
-
-      if( r < fConeRadius ) {
-        Double_t fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
-        // recalculating the centroid
-        Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
-
-        //Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()*track1->Pt()/fPt; // OB - remove
-       // OB - recalc phi via weighted 2vectors
-       jetVect2.SetMagPhi(jet->E()/fPt,jetVect2.Phi());
-       trackVect2.SetMagPhi(track1->Pt()/fPt,trackVect2.Phi());
+      if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
+      
+      for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop K0s in jet cone
        
-       TVector2 sumVect2 = jetVect2+trackVect2;
-       Double_t phi = sumVect2.Phi();
+       AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
+       if(!v0) continue;
 
-        //jet->SetPtEtaPhiE( 1., eta, phi, fPt );
-       ((TLorentzVector*) jet->MomentumVector())->SetPtEtaPhiE(1,eta,phi,fPt);
-
-       jet->AddTrack(track1);
-        tracks.Remove(track1);
+       Double_t invM           = v0->MassK0Short();
+       Float_t  trackPt        = v0->Pt();
+       Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
+       
+       //std::cout<<" trackPt "<<trackPt<<" invM "<<invM<<std::endl;
+       fFFHistosIMK0Cone->FillFF(trackPt, invM, jetPt, incrementJetPt);
       }
-    }
-    
-    tracks.Compress();
-
-    nTracks = tracks.GetEntries();
-    //   4- Continue until all particles are in a jet.
-    itrack.Reset();
-  } // end while nTracks
-  
-  // Convert to AODjets....
-  Int_t njets = jets->GetEntriesFast();
-
-  TClonesArray* aodjets = new TClonesArray("AliAODJet",njets);
-  aodjets->SetOwner(kTRUE);
-  
-  Int_t count = 0;
-
-  for(Int_t ijet=0; ijet<njets; ++ijet) {
-    //TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
-    AliAODJet* jet = (AliAODJet*)jets->At(ijet);
-
-    if (jet->E() < fJetPtCut) continue;
-    Float_t px, py,pz,en; // convert to 4-vector
-    px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
-    py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
-    pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
-    en = TMath::Sqrt(px * px + py * py + pz * pz);
-    jet->SetPxPyPzE(px,py,pz,en);
-
-    TClonesArray &tmpaodjets = *aodjets;
-    new(tmpaodjets[count++]) AliAODJet(*jet);
-    //aodjets->AddLast( new AliAODJet(*jet));    
-    //aodjets->AddLast( new AliAODJet(px, py, pz, en) );
-  }
-  // Order jets according to their pT .
-  //QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() ); // not for TClonesArray
-  
-  if (fDebug>3) AliInfo(Form(" %d Charged jets found - after cuts %d \n",njets,count));
-  
-  jets->Delete(); // OB: cleanup
-  delete jets;
-
-
-  return aodjets;
-}
-
-//____________________________________________________________________
-
-TClonesArray*  AliAnalysisTaskJetChem::GetPythiaJets(){
-
-  // return Pythia jets in jet acc
-
-  // note: present verion of AliMCEventHandler expects "AliAOD.root" (not "AliAODs.root"), otherwise galice.root will not be found in proper dir
-  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  if (!mcHandler) {
-    Printf("ERROR: Could not retrieve MC event handler");
-    return 0;
-  }
-  
-  AliMCEvent* mcEvent = mcHandler->MCEvent();
-  if (!mcEvent) {
-    Printf("ERROR: Could not retrieve MC event");
-    return 0;
-  }
-  
-  
-  AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
-  if(!pythiaGenHeader){
-    Printf("ERROR: Could not retrieve pythiaGenHeader");
-    return 0;
-  }
-  
-  //Get Jets from MC header
-  Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
-  
-  TClonesArray* aodjets = new TClonesArray("AliAODJet",nPythiaGenJets);
-  aodjets->SetOwner(kTRUE);
-
-  int count = 0;
-  for(int ip = 0; ip<nPythiaGenJets; ++ip){
-    Float_t p[4];
-    pythiaGenHeader->TriggerJet(ip,p);
-    TVector3 tempVect(p[0],p[1],p[2]);
-    if ( TMath::Abs(tempVect.Eta())>fJetEtaCut ) continue;
-
-    TClonesArray &tmpaodjets = *aodjets;
-    new(tmpaodjets[count++]) AliAODJet(p[0],p[1],p[2],p[3]);
-    //aodjets->AddLast(new AliAODJet(p[0],p[1],p[2],p[3]));
-   }
-  
-  return aodjets;
-}
+     if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum 
 
+       Bool_t incrementJetPt = kTRUE;
+       fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
+      }
 
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::QSortTracks(TObjArray &a, Int_t first, Int_t last)
-{
-  // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
-  
-  static TObject *tmp;
-  static int i;           // "static" to save stack space
-  int j;
-  
-  while (last - first > 1) {
-    i = first;
-    j = last;
-    for (;;) {
-      while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
-        ;
-      while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
-        ;
-      if (i >= j)
-        break;
-      
-      tmp  = a[i];
-      a[i] = a[j];
-      a[j] = tmp;
-    }
-    if (j == first) {
-      ++first;
-      continue;
+      delete jetConeK0list;
     }
-    tmp = a[first];
-    a[first] = a[j];
-    a[j] = tmp;
-    if (j - first < last - (j + 1)) {
-      QSortTracks(a, first, j);
-      first = j + 1;   // QSortTracks(j + 1, last);
-    } else {
-      QSortTracks(a, j + 1, last);
-      last = j;        // QSortTracks(first, j);
-    }
-  }
-}
-   
-//------------------------------------------------------------------
-
-TH1F* AliAnalysisTaskJetChem::CreatePIDhisto(const char* name){
-  
-  // create histogram
-  
-  TH1F* result = new TH1F(name,"",60,0,60);
-  result->SetOption("E");
-
-  // bin equal Geant ID
-
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(1),"photon");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(2),"e+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(3),"e-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(4),"e-neutrino"); 
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(5),"mu+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(6),"mu-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(7),"pi0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(8),"pi+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(9),"pi-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(10),"K long");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(11),"K+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(12),"K-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(13),"n");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(14),"p");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(15),"anti-proton");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(16),"K short");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(17),"eta");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(18),"Lambda");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(19),"Sigma+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(20),"Sigma0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(21),"Sigma-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(22),"Xi0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(23),"Xi-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(24),"Omega-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(25),"anti-neutron");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(26),"anti-Lambda");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(27),"anti-Sigma-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(28),"anti-Sigma0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(29),"anti-Sigma+"); 
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(30),"anti-Xi0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(31),"anti-Xi+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(32),"anti-Omega+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(33),"tau+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(34),"tau-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(35),"D+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(36),"D-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(37),"D0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(38),"anti-D0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(39),"Ds+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(40),"anti Ds-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(41),"Lamba_c+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(42),"W+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(43),"W-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(44),"Z0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(45),"d");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(46),"t");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(47),"alpha");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(48),"G_nu");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGpm311Bin),"K0/#bar{K0}");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDG333Bin),"phi");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGpm313Bin),"K*(892)0");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGp323Bin),"K*(892)+");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGm323Bin),"K*(892)-");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGNeutrinoBin),"nu");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGCharmedBaryonBin),"charmed baryon");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGQuarkBin),"q/#bar{q}");
-  result->GetXaxis()->SetBinLabel(result->GetXaxis()->FindBin(kPDGDiQuarkBin),"q #bar{q}");
-  result->GetXaxis()->LabelsOption("v"); // "u" ?
-  
-  return result;
-}
-
-//------------------------------------------------------------------
-
-TH1F*  AliAnalysisTaskJetChem::CreatePythiaIDhisto(const char* name){
-  
-  // create histogram
-  
-  TH1F* result = new TH1F(name,"",22,0,22);
-  result->SetOption("E");
-
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP11Bin,"qq #rightarrow qq");             // ISUB = 11 
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP12Bin,"q#bar{q} #rightarrow q#bar{q}"); // ISUB = 12
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP13Bin,"q#bar{q} #rightarrow gg");       // ISUB = 13
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP28Bin,"qg #rightarrow qg ");            // ISUB = 28 
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP53Bin,"gg #rightarrow q#bar{q}");       // ISUB = 53
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP68Bin,"gg #rightarrow gg");             // ISUB = 68
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP92Bin,"SD");                            // ISUB = 92
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP93Bin,"SD");                            // ISUB = 93
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP94Bin,"DD");                            // ISUB = 94
-  result->GetXaxis()->SetBinLabel(kPythiaPIDP95Bin,"low pt (MPI)");                  // ISUB = 95
-  result->GetXaxis()->SetBinLabel(kPythiaPIDPOtherBin,"other");                      // ISUB = XX
-
-  result->GetXaxis()->LabelsOption("v"); // "u" ?
-  
-  return result;
-}
-
-
-//------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::FillPythiaIDhisto(TH1F* h, const Int_t PID){
-
-  // fill Pyhtia PID histogram
-
-  Int_t bin = -1; 
-  if(PID == 11)      bin = kPythiaPIDP11Bin;
-  else if(PID == 12) bin = kPythiaPIDP12Bin;
-  else if(PID == 13) bin = kPythiaPIDP13Bin;
-  else if(PID == 28) bin = kPythiaPIDP28Bin;
-  else if(PID == 53) bin = kPythiaPIDP53Bin;
-  else if(PID == 68) bin = kPythiaPIDP68Bin;
-  else if(PID == 92) bin = kPythiaPIDP92Bin;
-  else if(PID == 93) bin = kPythiaPIDP93Bin;
-  else if(PID == 94) bin = kPythiaPIDP94Bin;
-  else if(PID == 95) bin = kPythiaPIDP95Bin;
-  else{ 
-    if(PID != -1) AliInfo(Form("unknown PID %d",PID));
-    return;
   }
   
-  h->Fill(h->GetBinCenter(bin),1);
-}
-
-
-//____________________________________________________________________
-
-
-void  AliAnalysisTaskJetChem::FillPIDhisto(TH1F* hist, Int_t pdg, Float_t weight){
-  
-  // convert pdg code to Geant ID and fill corresponding bin 
-
-  Int_t fGID = fpdgdb->ConvertPdgToGeant3(pdg);
-
-  //cout<<" pdg "<<pdg<<" fGID "<<fGID<<endl;
-
-  if(TMath::Abs(pdg) == 311) fGID = kPDGpm311Bin; 
-  if(pdg == 333)             fGID = kPDG333Bin; 
-  if(TMath::Abs(pdg) == 313) fGID = kPDGpm313Bin; 
-  if(pdg == 323)             fGID = kPDGp323Bin; 
-  if(pdg == -323)            fGID = kPDGm323Bin; 
-
-  if(TMath::Abs(pdg)==12 || TMath::Abs(pdg)==14 || TMath::Abs(pdg)==16) fGID = kPDGNeutrinoBin;
-
-  if(TMath::Abs(pdg)==4122) fGID = kPDGCharmedBaryonBin;
-
-  if(1<=TMath::Abs(pdg) && TMath::Abs(pdg)<=6) fGID = kPDGQuarkBin; 
-
-  if(TMath::Abs(pdg)==1103 || TMath::Abs(pdg)==2101 || TMath::Abs(pdg)==2103 || TMath::Abs(pdg)==2203 || 
-     TMath::Abs(pdg)==3101 || TMath::Abs(pdg)==3103 || TMath::Abs(pdg)==3201 || TMath::Abs(pdg)==3203 || 
-     TMath::Abs(pdg)==3303 || TMath::Abs(pdg)==4101 || TMath::Abs(pdg)==4103 || TMath::Abs(pdg)==4201 || 
-     TMath::Abs(pdg)==4203 || TMath::Abs(pdg)==4301 || TMath::Abs(pdg)==4303 || TMath::Abs(pdg)==4403 || 
-     TMath::Abs(pdg)==5101 || TMath::Abs(pdg)==5103 || TMath::Abs(pdg)==5201 || TMath::Abs(pdg)==5203 || 
-     TMath::Abs(pdg)==5301 || TMath::Abs(pdg)==5303 || TMath::Abs(pdg)==5401 || TMath::Abs(pdg)==5403 || 
-     TMath::Abs(pdg)==5503)  fGID = kPDGDiQuarkBin;
-    
-
-  hist->Fill(fGID,weight);
-
-  if(fGID>hist->GetBinLowEdge(hist->GetNbinsX()+1) || 
-     fGID<hist->GetBinLowEdge(1)){
-    
-    AliError(Form("fGID %d for pdg %d exceeding histo limits ",fGID,pdg));
-  }
-
-  if(fGID == 0){
-    AliInfo(Form("fGID 0 for pdg %d ",pdg));
-  }
-
-}
-
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::CreateHistos(){
-
-  // make histos
-  
-  fListOfHistos = new TList();
-  fListOfHistos->SetOwner(kTRUE);  
-
-  fhPrimVertexNCont = new TH1F("hPrimVertexNCont","",52,-2,50);
-  fhPrimVertexNCont->SetXTitle("");
-  fhPrimVertexNCont->Sumw2();
-  fListOfHistos->Add( fhPrimVertexNCont ); 
-
-  fhPrimVertexRho = new TH1F("hPrimVertexRho","",100,0,1);
-  fhPrimVertexRho->SetXTitle("");
-  fhPrimVertexRho->Sumw2();
-  fListOfHistos->Add( fhPrimVertexRho ); 
-
-  fhPrimVertexZ = new TH1F("hPrimVertexZ","",40,-20,20);
-  fhPrimVertexZ->SetXTitle("");
-  fhPrimVertexZ->Sumw2();
-  fListOfHistos->Add( fhPrimVertexZ ); 
-
-  fhNJets = new TH1F("hNJets","",10, 0, 10);
-  fhNJets->SetXTitle("# of jets");
-  fhNJets->Sumw2();
-  fListOfHistos->Add( fhNJets ); 
-
-  fhNJetsMC = new TH1F("hNJetsMC","",10,0,10);
-  fhNJetsMC->SetXTitle("# of jets");
-  fhNJetsMC->Sumw2();
-  fListOfHistos->Add( fhNJetsMC ); 
-
-  fhLeadingEta = new TH1F("hLeadingEta","Leading Jet eta",12,-0.6,0.6);
-  fhLeadingEta->SetXTitle("eta");
-  fhLeadingEta->SetYTitle("dN/deta");
-  fhLeadingEta->Sumw2();
-  fListOfHistos->Add(fhLeadingEta); 
-
-  fhLeadingNTracksVsEta = new TH2F("hLeadingNTracksVsEta","",20,-1.0,1.0,20,0,20);
-  fhLeadingNTracksVsEta->SetXTitle("eta");
-  fhLeadingNTracksVsEta->SetYTitle("# of tracks");
-  fhLeadingNTracksVsEta->Sumw2();
-  fListOfHistos->Add(fhLeadingNTracksVsEta); 
-
-  fhLeadingPtVsEta = new TH2F("hLeadingPtVsEta","",20,-1.0,1.0,50,0,50);
-  fhLeadingPtVsEta->SetXTitle("eta");
-  fhLeadingPtVsEta->SetYTitle("# of tracks");
-  fhLeadingPtVsEta->Sumw2();
-  fListOfHistos->Add(fhLeadingPtVsEta); 
-
-
-  fhLeadingPhi = new TH1F("hLeadingPhi","Leading Jet phi",63,0,6.3);
-  fhLeadingPhi->SetXTitle("phi");
-  fhLeadingPhi->SetYTitle("dN/dphi");
-  fhLeadingPhi->Sumw2();
-  fListOfHistos->Add(fhLeadingPhi);
-
-
-  fhLeadingPt  = new TH1F("hLeadingPt","leading Jet p_{T}",50,0,50);
-  fhLeadingPt->SetXTitle("p_{T} (GeV/c)");
-  fhLeadingPt->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhLeadingPt->Sumw2();
-  fListOfHistos->Add( fhLeadingPt ); 
-
-  fhLeadingPtDiffr  = new TH1F("hLeadingPtDiffr","leading Jet p_{T}",50,0,50);
-  fhLeadingPtDiffr->SetXTitle("P_{T} (GeV/c)");
-  fhLeadingPtDiffr->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhLeadingPtDiffr->Sumw2();
-  fListOfHistos->Add( fhLeadingPtDiffr ); 
-
-  fhLeadingEtaMC = new TH1F("hLeadingEtaMC","Leading Jet eta",12,-0.6,0.6);
-  fhLeadingEtaMC->SetXTitle("eta");
-  fhLeadingEtaMC->SetYTitle("dN/deta");
-  fhLeadingEtaMC->Sumw2();
-  fListOfHistos->Add(fhLeadingEtaMC); 
-
-  fhLeadingPhiMC = new TH1F("hLeadingPhiMC","Leading Jet phi",63,0,6.3);
-  fhLeadingPhiMC->SetXTitle("phi");
-  fhLeadingPhiMC->SetYTitle("dN/dphi");
-  fhLeadingPhiMC->Sumw2();
-  fListOfHistos->Add(fhLeadingPhiMC);
-
-  fhLeadingPtMC  = new TH1F("hLeadingPtMC","Leading Jet p_{T}",50,0,50);
-  fhLeadingPtMC->SetXTitle("p_{T} (GeV/c)");
-  fhLeadingPtMC->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhLeadingPtMC->Sumw2();
-  fListOfHistos->Add( fhLeadingPtMC ); 
-
-  fhLeadingPtMCDiffr  = new TH1F("hLeadingPtMCDiffr","Leading Jet p_{T}",50,0,50);
-  fhLeadingPtMCDiffr->SetXTitle("p_{T} (GeV/c)");
-  fhLeadingPtMCDiffr->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhLeadingPtMCDiffr->Sumw2();
-  fListOfHistos->Add( fhLeadingPtMCDiffr ); 
-
-  fhPhiEtaTracksNoCut = new TH2F("hPhiEtaTracksNoCut","phi vs eta tracks",20,-1.0,1.0,63,0,6.3);
-  fhPhiEtaTracksNoCut->SetXTitle("eta");
-  fhPhiEtaTracksNoCut->SetYTitle("phi");
-  fhPhiEtaTracksNoCut->Sumw2();
-  fListOfHistos->Add(fhPhiEtaTracksNoCut);
-
-  fhPtTracksNoCut = new TH1F("hPtTracksNoCut","p_{T} tracks",150,0,150);
-  fhPtTracksNoCut->SetXTitle("p_{T} (GeV)");
-  fhPtTracksNoCut->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhPtTracksNoCut->Sumw2();
-  fListOfHistos->Add(fhPtTracksNoCut);
-
-  fhPhiEtaTracks = new TH2F("hPhiEtaTracks","phi vs eta tracks",20,-1.0,1.0,63,0,6.3);
-  fhPhiEtaTracks->SetXTitle("eta");
-  fhPhiEtaTracks->SetYTitle("phi");
-  fhPhiEtaTracks->Sumw2();
-  fListOfHistos->Add(fhPhiEtaTracks);
-
-  fhPtTracks = new TH1F("hPtTracks","p_{T} tracks",150,0,150);
-  fhPtTracks->SetXTitle("P{T} (GeV)");
-  fhPtTracks->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhPtTracks->Sumw2();
-  fListOfHistos->Add(fhPtTracks);
-
-  fhTrackMult = new TH1F("hTrackMult","",150,0,150);
-  fhTrackMult->SetXTitle("n Tracks");
-  fhTrackMult->SetYTitle("counts");
-  fhTrackMult->Sumw2();
-  fListOfHistos->Add(fhTrackMult);
-
-  fhEtaMCTracks = new TH1F("hEtaMCTracks","eta tracks",30,-1.5,1.5);
-  fhEtaMCTracks->SetXTitle("eta");
-  fhEtaMCTracks->SetYTitle("dN/deta");
-  fhEtaMCTracks->Sumw2();
-  fListOfHistos->Add(fhEtaMCTracks);
-
-  fhPhiMCTracks = new TH1F("hPhiMCTracks","phi tracks",63,0,6.3);
-  fhPhiMCTracks->SetXTitle("phi");
-  fhPhiMCTracks->SetYTitle("dN/dphi");
-  fhPhiMCTracks->Sumw2();
-  fListOfHistos->Add(fhPhiMCTracks);
-
-  fhPtMCTracks = new TH1F("hPtMCTracks","p_{T} tracks",50,0,50);
-  fhPtMCTracks->SetXTitle("p_{T} (GeV)");
-  fhPtMCTracks->SetYTitle("dN/dp_{T} (1/GeV)");
-  fhPtMCTracks->Sumw2();
-  fListOfHistos->Add(fhPtMCTracks);
-
-  fhnTracksVsPtLeading  = new TH2F("hnTracksVsPtLeading","",50,0.,50.,20,-0.5,19.5);
-  fhnTracksVsPtLeading->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksVsPtLeading->SetYTitle("n tracks");
-  fhnTracksVsPtLeading->Sumw2();
-  fListOfHistos->Add( fhnTracksVsPtLeading );
-
-  fhdNdEtaPhiDist  = new TH1F("hdNdEtaPhiDist","Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
-  fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
-  fhdNdEtaPhiDist->SetYTitle("dN{ch}/d#etad#phi");
-  fhdNdEtaPhiDist->Sumw2();
-  fListOfHistos->Add( fhdNdEtaPhiDist );
-
-  fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt","P{T}^{90, max} vs Leading Jet P{T}",50,0.,50.);
-  fhRegionSumPtMaxVsEt->SetXTitle("P{T} (GeV/c)");
-  fhRegionSumPtMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionSumPtMaxVsEt ); 
-
-  fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt","N{ch}^{90, max} vs Leading Jet P{T}",50,0.,50.);
-  fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
-  fhRegionMultMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionMultMaxVsEt );
-
-  fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt","P{T}^{90, min} vs Leading Jet P{T}",50,0.,50.);
-  fhRegionSumPtMinVsEt->SetXTitle("P{T} (GeV/c)");
-  fhRegionSumPtMinVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionSumPtMinVsEt );
-    
-  fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt","N{ch}^{90, min} vs Leading Jet P{T}",50,0.,50.);
-  fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
-  fhRegionMultMinVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionMultMinVsEt );
-  
-  // V0s 
-
-  fhNV0s = new TH1F("hNV0s","n V0s",50,0,50);
-  fhNV0s->SetXTitle("n V0s");
-  fhNV0s->Sumw2();
-  fListOfHistos->Add(fhNV0s);
-
-  fhV0onFly = new TH1F("hV0onFly","on-the-fly V0",5,0,5);
-  fhV0onFly->SetXTitle("is on-the-fly V0");
-  fhV0onFly->Sumw2();
-  fListOfHistos->Add(fhV0onFly);
-
-  fhV0DCADaughters = new TH1F("hV0DCADaughters","V0 DCA daughters",200,0,2.0);
-  fhV0DCADaughters->SetXTitle("V0 DCA daughters");
-  fhV0DCADaughters->Sumw2();
-  fListOfHistos->Add(fhV0DCADaughters); 
-
-  fhV0Radius = new TH1F("hV0Radius","V0 radius",2500,0,250);
-  fhV0Radius->SetXTitle("V0 radius");
-  fhV0Radius->Sumw2();
-  fListOfHistos->Add(fhV0Radius); 
-
-  fhV0DCAToVertex = new TH1F("hV0DCAToVertex","",100,0,10);
-  fhV0DCAToVertex->SetXTitle("V0 DCA (cm)");
-  fhV0DCAToVertex->Sumw2();
-  fListOfHistos->Add(fhV0DCAToVertex);
-
-  fhV0DCAToVertexK0 = new TH1F("hV0DCAToVertexK0","",100,0,10); 
-  fhV0DCAToVertexK0->SetXTitle("V0 DCA (cm)");
-  fhV0DCAToVertexK0->Sumw2();
-  fListOfHistos->Add(fhV0DCAToVertexK0);
-
-  fhV0InvMassK0 = new TH1F("hV0InvMassK0","",2000,0,2);
-  fhV0InvMassK0->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0); 
-
-  fhV0PtVsInvMassK0 = new TH2F("hV0PtVsInvMassK0","",100,0.45,0.55,100,0.,50.);
-  fhV0PtVsInvMassK0->SetXTitle("inv mass (GeV)");
-  fhV0PtVsInvMassK0->SetYTitle("p_{T} (GeV)");
-  fhV0PtVsInvMassK0->Sumw2();
-  fListOfHistos->Add(fhV0PtVsInvMassK0); 
-
-  fhV0InvMassK0JetEvt = new TH1F("hV0InvMassK0JetEvt","",2000,0,2);
-  fhV0InvMassK0JetEvt->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0JetEvt->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0JetEvt); 
-
-  fhV0InvMassLambda = new TH1F("hV0InvMassLambda","",2000,0,2);
-  fhV0InvMassLambda->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambda->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambda);
-
-  fhV0InvMassAntiLambda = new TH1F("hV0InvMassAntiLambda","",2000,0,2);
-  fhV0InvMassAntiLambda->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambda->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambda);
-
-  fhV0InvMassLambdaJetEvt = new TH1F("hV0InvMassLambdaJetEvt","",2000,0,2);
-  fhV0InvMassLambdaJetEvt->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambdaJetEvt->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambdaJetEvt);
-
-  fhV0InvMassAntiLambdaJetEvt = new TH1F("hV0InvMassAntiLambdaJetEvt","",2000,0,2);
-  fhV0InvMassAntiLambdaJetEvt->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambdaJetEvt->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambdaJetEvt);
-
-  fhdROpanK0VsPt = new TH2F("hdROpanK0VsPt","V0 dR vs pt",100,0,10,100,0,10);
-  fhdROpanK0VsPt->SetXTitle("opening angle R (rad)");
-  fhdROpanK0VsPt->Sumw2();
-  fListOfHistos->Add(fhdROpanK0VsPt);
-
-  fhdPhiJetV0 = new TH1F("hdPhiJetV0","",640,-3.2,3.2);
-  fhdPhiJetV0->SetXTitle("#Delta #phi V0 - jet");
-  fhdPhiJetV0->Sumw2();
-  fListOfHistos->Add(fhdPhiJetV0); 
-
-  fhdPhiJetK0 = new TH1F("hdPhiJetK0","",640,-3.2,3.2);
-  fhdPhiJetK0->SetXTitle("#Delta #phi V0 - jet");
-  fhdPhiJetK0->Sumw2();
-  fListOfHistos->Add(fhdPhiJetK0); 
-
-  fhdRJetK0 = new TH1F("hdRJetK0","dN/dR K0-jet",500,0,5);
-  fhdRJetK0->SetXTitle("#Delta R K0 - jet");
-  fhdRJetK0->SetYTitle("1/N_{jet} dN/dR");
-  fhdRJetK0->Sumw2();
-  fListOfHistos->Add(fhdRJetK0); 
-
-  fhdNdptV0 = new TH1F("hdNdptV0","dN/dpt V0",100,0,10);
-  fhdNdptV0->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptV0->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptV0->Sumw2();
-  fListOfHistos->Add(fhdNdptV0);
-
-  fhdNdptK0 = new TH1F("hdNdptK0","",100,0,10);
-  fhdNdptK0->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0->Sumw2();
-  fListOfHistos->Add(fhdNdptK0); 
-
-  fhPtVsEtaK0 = new TH2F("hPtVsEtaK0","",20,-1,1,100,0,10);
-  fhPtVsEtaK0->SetXTitle("#eta");  
-  fhPtVsEtaK0->SetYTitle("p_{T} (GeV/c)");
-  fhPtVsEtaK0->Sumw2();
-  fListOfHistos->Add(fhPtVsEtaK0); 
-
-  fhV0InvMassK0DCA = new TH1F("hV0InvMassK0DCA","",2000,0,2);
-  fhV0InvMassK0DCA->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0DCA->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0DCA); 
-
-  fhV0InvMassK0DCAdEdx = new TH1F("hV0InvMassK0DCAdEdx","",2000,0,2);
-  fhV0InvMassK0DCAdEdx->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0DCAdEdx->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0DCAdEdx); 
-
-  fhV0PtVsInvMassK0DCAdEdx = new TH2F("hV0PtVsInvMassK0DCAdEdx","",100,0.45,0.55,100,0.,50.);
-  fhV0PtVsInvMassK0DCAdEdx->SetXTitle("inv mass (GeV)");
-  fhV0PtVsInvMassK0DCAdEdx->SetYTitle("p_{T} (GeV)");
-  fhV0PtVsInvMassK0DCAdEdx->Sumw2();
-  fListOfHistos->Add(fhV0PtVsInvMassK0DCAdEdx); 
-
-  fhV0InvMassK0DCAPID = new TH1F("hV0InvMassK0DCAPID","",2000,0,2);
-  fhV0InvMassK0DCAPID->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0DCAPID->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0DCAPID); 
-
-  fhV0InvMassLambdaDCAdEdx = new TH1F("hV0InvMassLambdaDCAdEdx","",2000,0,2);
-  fhV0InvMassLambdaDCAdEdx->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambdaDCAdEdx->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambdaDCAdEdx); 
-
-  fhV0InvMassAntiLambdaDCAdEdx = new TH1F("hV0InvMassAntiLambdaDCAdEdx","",2000,0,2);
-  fhV0InvMassAntiLambdaDCAdEdx->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambdaDCAdEdx->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambdaDCAdEdx); 
-
-  fhdNdptK0DCA = new TH1F("hdNdptK0DCA","",100,0,10);
-  fhdNdptK0DCA->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0DCA->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0DCA->Sumw2();
-  fListOfHistos->Add(fhdNdptK0DCA); 
-
-  fhdNdptK0DCAdEdx = new TH1F("hdNdptK0DCAdEdx","",100,0,10);
-  fhdNdptK0DCAdEdx->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0DCAdEdx->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0DCAdEdx->Sumw2();
-  fListOfHistos->Add(fhdNdptK0DCAdEdx); 
-
-  fhV0InvMassK0Min = new TH1F("hV0InvMassK0Min","",2000,0,2);
-  fhV0InvMassK0Min->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0Min->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0Min); 
-
-  fhV0InvMassLambdaMin = new TH1F("hV0InvMassLambdaMin","",2000,0,2);
-  fhV0InvMassLambdaMin->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambdaMin->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambdaMin); 
-
-  fhV0InvMassAntiLambdaMin = new TH1F("hV0InvMassAntiLambdaMin","",2000,0,2);
-  fhV0InvMassAntiLambdaMin->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambdaMin->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambdaMin);
-
-  fhV0InvMassK0Max = new TH1F("hV0InvMassK0Max","",2000,0,2);
-  fhV0InvMassK0Max->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0Max->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0Max); 
-
-  fhV0InvMassLambdaMax = new TH1F("hV0InvMassLambdaMax","",2000,0,2);
-  fhV0InvMassLambdaMax->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambdaMax->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambdaMax); 
-
-  fhV0InvMassAntiLambdaMax = new TH1F("hV0InvMassAntiLambdaMax","",2000,0,2);
-  fhV0InvMassAntiLambdaMax->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambdaMax->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambdaMax); 
-
-  fhV0InvMassK0Jet = new TH1F("hV0InvMassK0Jet","",2000,0,2);
-  fhV0InvMassK0Jet->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0Jet->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0Jet); 
-
-  fhV0InvMassLambdaJet = new TH1F("hV0InvMassLambdaJet","",2000,0,2);
-  fhV0InvMassLambdaJet->SetXTitle("inv mass (GeV)");
-  fhV0InvMassLambdaJet->Sumw2();
-  fListOfHistos->Add(fhV0InvMassLambdaJet); 
-
-  fhV0InvMassAntiLambdaJet = new TH1F("hV0InvMassAntiLambdaJet","",2000,0,2);
-  fhV0InvMassAntiLambdaJet->SetXTitle("inv mass (GeV)");
-  fhV0InvMassAntiLambdaJet->Sumw2();
-  fListOfHistos->Add(fhV0InvMassAntiLambdaJet); 
-
-  fhV0InvMassK0Lambda = new TH1F("hV0InvMassK0Lambda","",2000,0,2);
-  fhV0InvMassK0Lambda->SetXTitle("inv mass (GeV)");
-  fhV0InvMassK0Lambda->Sumw2();
-  fListOfHistos->Add(fhV0InvMassK0Lambda); 
-
-  fhdNdptK0JetEvt = new TH1F("hdNdptK0JetEvt","",100,0,10);
-  fhdNdptK0JetEvt->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0JetEvt->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0JetEvt->Sumw2();
-  fListOfHistos->Add(fhdNdptK0JetEvt); 
-  
-  fhdNdzK0 = new TH1F("hdNdzK0","",150,0,1.5);
-  fhdNdzK0->SetXTitle("z");
-  fhdNdzK0->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0->Sumw2();
-  fListOfHistos->Add(fhdNdzK0);
-
-  fhdNdzK05to10 = new TH1F("hdNdzK05to10","",150,0,1.5);
-  fhdNdzK05to10->SetXTitle("z");
-  fhdNdzK05to10->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK05to10->Sumw2();
-  fListOfHistos->Add(fhdNdzK05to10);
-
-  fhdNdzK010to20 = new TH1F("hdNdzK010to20","",150,0,1.5);
-  fhdNdzK010to20->SetXTitle("z");
-  fhdNdzK010to20->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK010to20->Sumw2();
-  fListOfHistos->Add(fhdNdzK010to20);
-
-  fhdNdzK020to30 = new TH1F("hdNdzK020to30","",150,0,1.5);
-  fhdNdzK020to30->SetXTitle("z");
-  fhdNdzK020to30->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK020to30->Sumw2();
-  fListOfHistos->Add(fhdNdzK020to30);
-
-  fhdNdzK030to40 = new TH1F("hdNdzK030to40","",150,0,1.5);
-  fhdNdzK030to40->SetXTitle("z");
-  fhdNdzK030to40->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK030to40->Sumw2();
-  fListOfHistos->Add(fhdNdzK030to40);
-
-  fhdNdzK040to60 = new TH1F("hdNdzK040to60","",150,0,1.5);
-  fhdNdzK040to60->SetXTitle("z");
-  fhdNdzK040to60->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK040to60->Sumw2();
-  fListOfHistos->Add(fhdNdzK040to60);
-
-  fhdNdxiK0 = new TH1F("hdNdxiK0","",100,0,10);
-  fhdNdxiK0->SetXTitle("xi");
-  fhdNdxiK0->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiK0->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0);
-  
-  fhdNdzLambda = new TH1F("hdNdzLambda","",150,0,1.5);
-  fhdNdzLambda->SetXTitle("z");
-  fhdNdzLambda->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzLambda->Sumw2();
-  fListOfHistos->Add(fhdNdzLambda);
-
-  fhdNdzAntiLambda = new TH1F("hdNdzAntiLambda","",150,0,1.5);
-  fhdNdzAntiLambda->SetXTitle("z");
-  fhdNdzAntiLambda->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzAntiLambda->Sumw2();
-  fListOfHistos->Add(fhdNdzAntiLambda);
-
-  fhdNdzK0Max = new TH1F("hdNdzK0Max","",150,0,1.5);
-  fhdNdzK0Max->SetXTitle("z");
-  fhdNdzK0Max->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0Max->Sumw2();
-  fListOfHistos->Add(fhdNdzK0Max);
-
-  fhdNdxiK0Max = new TH1F("hdNdxiK0Max","",100,0,10);
-  fhdNdxiK0Max->SetXTitle("xi");
-  fhdNdxiK0Max->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiK0Max->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0Max);
-
-  fhdNdzLambdaMax = new TH1F("hdNdzLambdaMax","",150,0,1.5);
-  fhdNdzLambdaMax->SetXTitle("z");
-  fhdNdzLambdaMax->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzLambdaMax->Sumw2();
-  fListOfHistos->Add(fhdNdzLambdaMax); 
-
-  fhdNdxiLambdaMax = new TH1F("hdNdxiLambdaMax","",700,0,7);
-  fhdNdxiLambdaMax->SetXTitle("xi");
-  fhdNdxiLambdaMax->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiLambdaMax->Sumw2();
-  fListOfHistos->Add(fhdNdxiLambdaMax); 
-
-  fhdNdptK0Max = new TH1F("hdNdptK0Max","",100,0,10);
-  fhdNdptK0Max->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0Max->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0Max->Sumw2();
-  fListOfHistos->Add(fhdNdptK0Max); 
-
-  fhdNdptLambdaMax = new TH1F("hdNdptLambdaMax","",100,0,10);
-  fhdNdptLambdaMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaMax->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMax); 
-
-  fhdNdzK0Min = new TH1F("hdNdzK0Min","",150,0,1.5);
-  fhdNdzK0Min->SetXTitle("z");
-  fhdNdzK0Min->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0Min->Sumw2();
-  fListOfHistos->Add(fhdNdzK0Min); 
-
-  fhdNdxiK0Min = new TH1F("hdNdxiK0Min","",100,0,10);
-  fhdNdxiK0Min->SetXTitle("xi");
-  fhdNdxiK0Min->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiK0Min->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0Min); 
-
-  fhdNdzLambdaMin = new TH1F("hdNdzLambdaMin","",150,0,1.5);
-  fhdNdzLambdaMin->SetXTitle("z");
-  fhdNdzLambdaMin->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzLambdaMin->Sumw2();
-  fListOfHistos->Add(fhdNdzLambdaMin); 
-
-  fhdNdxiLambdaMin = new TH1F("hdNdxiLambdaMin","",700,0,7);
-  fhdNdxiLambdaMin->SetXTitle("xi");
-  fhdNdxiLambdaMin->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiLambdaMin->Sumw2();
-  fListOfHistos->Add(fhdNdxiLambdaMin); 
-
-  fhdNdptK0Min = new TH1F("hdNdptK0Min","",100,0,10);
-  fhdNdptK0Min->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0Min->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0Min->Sumw2();
-  fListOfHistos->Add(fhdNdptK0Min); 
-
-  fhdNdptLambdaMin = new TH1F("hdNdptLambdaMin","",100,0,10);
-  fhdNdptLambdaMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaMin->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMin); 
-
-  fhdNdzK0Jet = new TH1F("hdNdzK0Jet","",150,0,1.5);
-  fhdNdzK0Jet->SetXTitle("z");
-  fhdNdzK0Jet->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0Jet->Sumw2();
-  fListOfHistos->Add(fhdNdzK0Jet);
-
-  fhdNdxiK0Jet = new TH1F("hdNdxiK0Jet","",100,0,10);
-  fhdNdxiK0Jet->SetXTitle("xi");
-  fhdNdxiK0Jet->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiK0Jet->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0Jet); 
-
-  fhdNdzLambdaJet = new TH1F("hdNdzLambdaJet","",150,0,1.5);
-  fhdNdzLambdaJet->SetXTitle("z");
-  fhdNdzLambdaJet->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzLambdaJet->Sumw2();
-  fListOfHistos->Add(fhdNdzLambdaJet); 
-
-  fhdNdxiLambdaJet = new TH1F("hdNdxiLambdaJet","",700,0,7);
-  fhdNdxiLambdaJet->SetXTitle("xi");
-  fhdNdxiLambdaJet->SetYTitle("1/N_{jet} dN/dxi");
-  fhdNdxiLambdaJet->Sumw2();
-  fListOfHistos->Add(fhdNdxiLambdaJet); 
-
-  fhdNdptK0Jet = new TH1F("hdNdptK0Jet","",100,0,10);
-  fhdNdptK0Jet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0Jet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0Jet->Sumw2();
-  fListOfHistos->Add(fhdNdptK0Jet); 
-
-  fhdNdptLambdaJet = new TH1F("hdNdptLambdaJet","",100,0,10);
-  fhdNdptLambdaJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaJet->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaJet); 
-
-  fhdEdxVsMomV0 = new TH2F("hdEdxVsMomV0","",200,-10,10,200,0,2000);
-  fhdEdxVsMomV0->SetXTitle("mom (GeV/c)");
-  fhdEdxVsMomV0->SetYTitle("dE/dx");
-  fhdEdxVsMomV0->Sumw2();
-  fListOfHistos->Add(fhdEdxVsMomV0); 
-
-  fhdEdxVsMomV0pidEdx = new TH2F("hdEdxVsMomV0pidEdx","",200,-10,10,200,0,2000);
-  fhdEdxVsMomV0pidEdx->SetXTitle("mom (GeV/c)");
-  fhdEdxVsMomV0pidEdx->SetYTitle("dE/dx");
-  fhdEdxVsMomV0pidEdx->Sumw2();
-  fListOfHistos->Add(fhdEdxVsMomV0pidEdx); 
-
-  fhdEdxVsMomV0piPID = new TH2F("hdEdxVsMomV0piPID","",200,-10,10,200,0,2000);
-  fhdEdxVsMomV0piPID->SetXTitle("mom (GeV/c)");
-  fhdEdxVsMomV0piPID->SetYTitle("dE/dx");
-  fhdEdxVsMomV0piPID->Sumw2();
-  fListOfHistos->Add(fhdEdxVsMomV0piPID); 
-
-  fhdPhiJetK0MC = new TH1F("hdPhiJetK0MC","",640,-3.2,3.2);
-  fhdPhiJetK0MC->SetXTitle("#Delta #phi K0 - jet");
-  fhdPhiJetK0MC->SetYTitle("1/N_{jet} dN/dphi");
-  fhdPhiJetK0MC->Sumw2();
-  fListOfHistos->Add(fhdPhiJetK0MC); 
-
-  fhdRJetK0MC   = new TH1F("hdRJetK0MC","dN/R K0-jet",500,0,5);
-  fhdRJetK0MC->SetXTitle("#Delta R K0 - jet");
-  fhdRJetK0MC->SetYTitle("1/N_{jet} dN/dR");
-  fhdRJetK0MC->Sumw2();
-  fListOfHistos->Add(fhdRJetK0MC); 
-
-  fhdRV0MC =  new TH1F("hdRV0MC","",500,0.,1.);
-  fhdRV0MC->SetXTitle("#Delta R");
-  fhdRV0MC->SetYTitle("");
-  fhdRV0MC->Sumw2();
-  fListOfHistos->Add(fhdRV0MC); 
-
-  fhdNdptchPiMCMax = new TH1F("hdNdptchPiMCMax","",100,0,10);
-  fhdNdptchPiMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchPiMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchPiMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptchPiMCMax); 
-
-  fhdNdptK0MCMax = new TH1F("hdNdptK0MCMax","",100,0,10);
-  fhdNdptK0MCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0MCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0MCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptK0MCMax); 
-
-  fhdNdptchKMCMax = new TH1F("hdNdptchKMCMax","",100,0,10);
-  fhdNdptchKMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchKMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchKMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptchKMCMax); 
-
-  fhdNdptpMCMax = new TH1F("hdNdptpMCMax","",100,0,10);
-  fhdNdptpMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptpMCMax); 
-
-  fhdNdptpBarMCMax = new TH1F("hdNdptpBarMCMax","",100,0,10);
-  fhdNdptpBarMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpBarMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpBarMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptpBarMCMax); 
-
-  fhdNdptLambdaMCMax = new TH1F("hdNdptLambdaMCMax","",100,0,10);
-  fhdNdptLambdaMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMCMax); 
-
-  fhdNdptLambdaBarMCMax = new TH1F("hdNdptLambdaBarMCMax","",100,0,10);
-  fhdNdptLambdaBarMCMax->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaBarMCMax->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaBarMCMax->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaBarMCMax); 
-
-
-  fhdNdptchPiMCMin = new TH1F("hdNdptchPiMCMin","",100,0,10);
-  fhdNdptchPiMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchPiMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchPiMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptchPiMCMin); 
-  fhdNdptK0MCMin = new TH1F("hdNdptK0MCMin","",100,0,10);
-  fhdNdptK0MCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0MCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0MCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptK0MCMin); 
-
-  fhdNdptchKMCMin = new TH1F("hdNdptchKMCMin","",100,0,10);
-  fhdNdptchKMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchKMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchKMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptchKMCMin); 
-
-  fhdNdptpMCMin = new TH1F("hdNdptpMCMin","",100,0,10);
-  fhdNdptpMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptpMCMin); 
-
-  fhdNdptpBarMCMin = new TH1F("hdNdptpBarMCMin","",100,0,10);
-  fhdNdptpBarMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpBarMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpBarMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptpBarMCMin); 
-
-  fhdNdptLambdaMCMin = new TH1F("hdNdptLambdaMCMin","",100,0,10);
-  fhdNdptLambdaMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMCMin); 
-
-  fhdNdptLambdaBarMCMin = new TH1F("hdNdptLambdaBarMCMin","",100,0,10);
-  fhdNdptLambdaBarMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaBarMCMin->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaBarMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaBarMCMin); 
-
-  fhdNdptOmegaMCMin = new TH1F("hdNdptOmegaMCMin","",100,0,10);;
-  fhdNdptOmegaMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptOmegaMCMin->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptOmegaMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptOmegaMCMin); 
-
-  fhdNdptOmegaBarMCMin = new TH1F("hdNdptOmegaBarMCMin","",100,0,10);;
-  fhdNdptOmegaBarMCMin->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptOmegaBarMCMin->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptOmegaBarMCMin->Sumw2();
-  fListOfHistos->Add(fhdNdptOmegaBarMCMin); 
-
-  fhdNdptchPiMCJet = new TH1F("hdNdptchPiMCJet","",100,0,10);
-  fhdNdptchPiMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchPiMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchPiMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptchPiMCJet); 
-
-  fhdNdptK0MCJet = new TH1F("hdNdptK0MCJet","",100,0,10);
-  fhdNdptK0MCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0MCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptK0MCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptK0MCJet); 
-
-  fhdNdptchKMCJet = new TH1F("hdNdptchKMCJet","",100,0,10);
-  fhdNdptchKMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchKMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptchKMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptchKMCJet);
-
-  fhdNdptpMCJet = new TH1F("hdNdptpMCJet","",100,0,10);
-  fhdNdptpMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptpMCJet); 
-
-  fhdNdptpBarMCJet = new TH1F("hdNdptpBarMCJet","",100,0,10);
-  fhdNdptpBarMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpBarMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptpBarMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptpBarMCJet);
-
-  fhdNdptLambdaMCJet = new TH1F("hdNdptLambdaMCJet","",100,0,10);
-  fhdNdptLambdaMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMCJet); 
-
-  fhdNdptLambdaBarMCJet = new TH1F("hdNdptLambdaBarMCJet","",100,0,10);
-  fhdNdptLambdaBarMCJet->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaBarMCJet->SetYTitle("1/N_{jet} dN/dpt (1/GeV)");
-  fhdNdptLambdaBarMCJet->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaBarMCJet); 
-
-  fhPIDMC = CreatePIDhisto("hPIDMC");
-  fhPIDMC->Sumw2();
-  fListOfHistos->Add(fhPIDMC); 
-
-  fhPIDMC_quarkEv = CreatePIDhisto("hPIDMC_quarkEv");
-  fhPIDMC_quarkEv->Sumw2();
-  fListOfHistos->Add(fhPIDMC_quarkEv); 
-
-  fhPIDMC_gluonEv = CreatePIDhisto("hPIDMC_gluonEv");
-  fhPIDMC_gluonEv->Sumw2();
-  fListOfHistos->Add(fhPIDMC_gluonEv); 
-
-  fhPIDMCAll = CreatePIDhisto("hPIDMCAll");
-  fhPIDMCAll->Sumw2();
-  fListOfHistos->Add(fhPIDMCAll); 
-
-  fhPIDMCMin = CreatePIDhisto("hPIDMCMin");
-  fhPIDMCMin->Sumw2();
-  fListOfHistos->Add(fhPIDMCMin); 
-
-  fhPIDMCJet = CreatePIDhisto("hPIDMCJet");
-  fhPIDMCJet->Sumw2();
-  fListOfHistos->Add(fhPIDMCJet); 
-
-  fhPIDMCMotherK0 = CreatePIDhisto("hPIDMCMotherK0"); 
-  fhPIDMCMotherK0->Sumw2();
-  fListOfHistos->Add(fhPIDMCMotherK0); 
-
-  fhPIDMCGrandMotherK0 = CreatePIDhisto("hPIDMCGrandMotherK0"); 
-  fhPIDMCGrandMotherK0->Sumw2();
-  fListOfHistos->Add(fhPIDMCGrandMotherK0); 
-
-  fhPIDMCMotherChK = CreatePIDhisto("hPIDMCMotherChK"); 
-  fhPIDMCMotherChK->Sumw2();
-  fListOfHistos->Add(fhPIDMCMotherChK); 
-
-  fhPIDMCMotherK0Trans = CreatePIDhisto("hPIDMCMotherK0Trans"); 
-  fhPIDMCMotherK0Trans->Sumw2();
-  fListOfHistos->Add(fhPIDMCMotherK0Trans); 
-
-  fhPIDMCGrandMotherK0Trans = CreatePIDhisto("hPIDMCGrandMotherK0Trans"); 
-  fhPIDMCGrandMotherK0Trans->Sumw2();
-  fListOfHistos->Add(fhPIDMCGrandMotherK0Trans); 
-
-  fhPIDMCMotherChKTrans = CreatePIDhisto("hPIDMCMotherChKTrans"); 
-  fhPIDMCMotherChKTrans->Sumw2();
-  fListOfHistos->Add(fhPIDMCMotherChKTrans); 
-
-  fhdNdptgammaMC = new TH1F("hdNdptgammaMC","",100,0,10);
-  fhdNdptgammaMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptgammaMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptgammaMC->Sumw2();
-  fListOfHistos->Add(fhdNdptgammaMC); 
-
-  fhdNdptchPiMC  = new TH1F("hdNdptchPiMC","",100,0,10);;
-  fhdNdptchPiMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchPiMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptchPiMC->Sumw2();
-  fListOfHistos->Add(fhdNdptchPiMC); 
-
-  fhdNdptpi0MC    = new TH1F("hdNdptpi0MC","",100,0,10);;
-  fhdNdptpi0MC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpi0MC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptpi0MC->Sumw2();
-  fListOfHistos->Add(fhdNdptpi0MC); 
-
-  fhdNdptK0MC = new TH1F("hdNdptK0MC","",100,0,10);;
-  fhdNdptK0MC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0MC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptK0MC->Sumw2();
-  fListOfHistos->Add(fhdNdptK0MC); 
-
-  fhdNdptchKMC    = new TH1F("hdNdptchKMC","",100,0,10);;
-  fhdNdptchKMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptchKMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptchKMC->Sumw2();
-  fListOfHistos->Add(fhdNdptchKMC); 
-
-  fhdNdptpMC    = new TH1F("hdNdptpMC","",100,0,10);;
-  fhdNdptpMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptpMC->Sumw2();
-  fListOfHistos->Add(fhdNdptpMC); 
-
-  fhdNdptpBarMC    = new TH1F("hdNdptpBarMC","",100,0,10);;
-  fhdNdptpBarMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptpBarMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptpBarMC->Sumw2();
-  fListOfHistos->Add(fhdNdptpBarMC); 
-
-  fhdNdptLambdaMC    = new TH1F("hdNdptLambdaMC","",100,0,10);;
-  fhdNdptLambdaMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptLambdaMC->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaMC); 
-
-  fhdNdptLambdaBarMC    = new TH1F("hdNdptLambdaBarMC","",100,0,10);;
-  fhdNdptLambdaBarMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptLambdaBarMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptLambdaBarMC->Sumw2();
-  fListOfHistos->Add(fhdNdptLambdaBarMC); 
-
-  fhdNdptOmegaMC = new TH1F("hdNdptOmegaMC","",100,0,10);;
-  fhdNdptOmegaMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptOmegaMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptOmegaMC->Sumw2();
-  fListOfHistos->Add(fhdNdptOmegaMC); 
-
-  fhdNdptOmegaBarMC = new TH1F("hdNdptOmegaBarMC","",100,0,10);;
-  fhdNdptOmegaBarMC->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptOmegaBarMC->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptOmegaBarMC->Sumw2();
-  fListOfHistos->Add(fhdNdptOmegaBarMC); 
-
-  fhdNdxiMC = new TH1F("hdNdxiMC","",100,0,10); 
-  fhdNdxiMC->SetXTitle("#xi");
-  fhdNdxiMC->Sumw2();
-  fListOfHistos->Add(fhdNdxiMC);            
-
-  fhdNdxiK0MC = new TH1F("hdNdxiK0MC","",100,0,10); 
-  fhdNdxiK0MC->SetXTitle("#xi");
-  fhdNdxiK0MC->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0MC);            
-
-  fhdNdxiK0MCJet = new TH1F("hdNdxiK0MCJet","",100,0,10); 
-  fhdNdxiK0MCJet->SetXTitle("#xi");
-  fhdNdxiK0MCJet->Sumw2();
-  fListOfHistos->Add(fhdNdxiK0MCJet);
-
-  fhdNdzK0MC = new TH1F("hdNdzK0MC","",150,0,1.5);
-  fhdNdzK0MC->SetXTitle("z");
-  fhdNdzK0MC->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0MC->Sumw2();
-  fListOfHistos->Add(fhdNdzK0MC);
-
-  fhdNdzK0MCJet = new TH1F("hdNdzK0MCJet","",150,0,1.5);
-  fhdNdzK0MCJet->SetXTitle("z");
-  fhdNdzK0MCJet->SetYTitle("1/N_{jet} dN/dz");
-  fhdNdzK0MCJet->Sumw2();
-  fListOfHistos->Add(fhdNdzK0MCJet);
-
-  fhdNdptK0MCJetEvt = new TH1F("hdNdptK0MCJetEvt","",100,0,10);;
-  fhdNdptK0MCJetEvt->SetXTitle("p_{T} (GeV/c)");
-  fhdNdptK0MCJetEvt->SetYTitle("1/N_{event} dN/dpt (1/GeV)");
-  fhdNdptK0MCJetEvt->Sumw2();
-  fListOfHistos->Add(fhdNdptK0MCJetEvt);
-
-  fhnJetsAODvsMC = new TH2F("hnJetsAODvsMC","",20,0,20,20,0.,20.); 
-  fhnJetsAODvsMC->SetXTitle("n jets MC");
-  fhnJetsAODvsMC->SetXTitle("n jets AOD");
-  fhnJetsAODvsMC->Sumw2();
-  fListOfHistos->Add(fhnJetsAODvsMC);
-
-  fhLeadingPtAODvsMC = new TH2F("hLeadingPtAODvsMC","",20,0,20,20,0.,20.); 
-  fhLeadingPtAODvsMC->SetXTitle("p_{T} MC (GeV/c)");
-  fhLeadingPtAODvsMC->SetYTitle("p_{T} AOD (GeV/c)");
-  fhLeadingPtAODvsMC->Sumw2();
-  fListOfHistos->Add(fhLeadingPtAODvsMC);
-
-  fhLeadingEtaAODvsMC = new TH2F("hLeadingEtaAODvsMC","",20,-1,1,20,-1.,1.); 
-  fhLeadingEtaAODvsMC->SetXTitle("#eta MC");
-  fhLeadingEtaAODvsMC->SetYTitle("#eta AOD");
-  fhLeadingEtaAODvsMC->Sumw2();
-  fListOfHistos->Add(fhLeadingEtaAODvsMC);
-
-  fhLeadingPhiAODvsMC = new TH2F("hLeadingPhiAODvsMC","",63,0,6.3,63,0.,6.3); 
-  fhLeadingPhiAODvsMC->SetXTitle("#phi MC");
-  fhLeadingPhiAODvsMC->SetYTitle("#phi AOD");
-  fhLeadingPhiAODvsMC->Sumw2();
-  fListOfHistos->Add(fhLeadingPhiAODvsMC);
-
-
-  fhnTracksLeadingAODvsMC = new TH2F("hnTracksLeadingAODvsMC","",20,0.,20,20,-0.5,19.5);
-  fhnTracksLeadingAODvsMC->SetXTitle("nTracks MC");
-  fhnTracksLeadingAODvsMC->SetYTitle("nTracks AOD");
-  fhnTracksLeadingAODvsMC->Sumw2();
-  fListOfHistos->Add(fhnTracksLeadingAODvsMC);
-
-  fhLeadingdRAODMC = new TH1F("hLeadingdRAODMC","",40,0.,4);
-  fhLeadingdRAODMC->SetXTitle("#Delta R");
-  fhLeadingdRAODMC->Sumw2();
-  fListOfHistos->Add(fhLeadingdRAODMC);
-
-
-  fhLeadingPtAODvsMCdRcut = new TH2F("hLeadingPtAODvsMCdRcut","",40,0,20,40,0.,20.); 
-  fhLeadingPtAODvsMCdRcut->SetXTitle("p_{T} MC (GeV/c)");
-  fhLeadingPtAODvsMCdRcut->SetYTitle("p_{T} AOD (GeV/c)");
-  fhLeadingPtAODvsMCdRcut->Sumw2();
-  fListOfHistos->Add(fhLeadingPtAODvsMCdRcut);
-
-
-  fhdnTracksVsdPtLeadingAODMC = new TH2F("hdnTracksVsdPtLeadingAODMC","",40,-10.,10,20,-10.5,9.5);
-  fhdnTracksVsdPtLeadingAODMC->SetXTitle("#Delta Pt AOD-MC (GeV/c)");
-  fhdnTracksVsdPtLeadingAODMC->SetYTitle("#Delta N AOD-MC");
-  fhdnTracksVsdPtLeadingAODMC->Sumw2();
-  fListOfHistos->Add(fhdnTracksVsdPtLeadingAODMC);
-
-  fhnTracksJetVsPtAOD = new TH2F("hnTracksJetVsPtAOD","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtAOD->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtAOD->SetYTitle("nTracks");
-  fhnTracksJetVsPtAOD->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtAOD);
-
-  fhnTracksJetVsPtAODquarkEv = new TH2F("hnTracksJetVsPtAODquarkEv","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtAODquarkEv->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtAODquarkEv->SetYTitle("nTracks");
-  fhnTracksJetVsPtAODquarkEv->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtAODquarkEv);
-
-  fhRadiusJetVsPtAOD = new TH2F("hRadiusJetVsPtAOD","",50,0,50,10,0.,1);
-  fhRadiusJetVsPtAOD->SetXTitle("p_{T} (GeV/c)");
-  fhRadiusJetVsPtAOD->SetYTitle("radius");
-  fhRadiusJetVsPtAOD->Sumw2();
-  fListOfHistos->Add(fhRadiusJetVsPtAOD);
-
-  fhnTracksJetVsPtMC = new TH2F("hnTracksJetVsPtMC","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtMC->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtMC->SetYTitle("nTracks");
-  fhnTracksJetVsPtMC->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtMC);
-
-  fhnTracksJetVsPtMCquarkEv = new TH2F("hnTracksJetVsPtMCquarkEv","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtMCquarkEv->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtMCquarkEv->SetYTitle("nTracks");
-  fhnTracksJetVsPtMCquarkEv->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtMCquarkEv);
-
-  fhRadiusJetVsPtMC = new TH2F("hRadiusJetVsPtMC","",50,0,50,10,0.,1);
-  fhRadiusJetVsPtMC->SetXTitle("p_{T} (GeV/c)");
-  fhRadiusJetVsPtMC->SetYTitle("radius");
-  fhRadiusJetVsPtMC->Sumw2();
-  fListOfHistos->Add(fhRadiusJetVsPtMC);
-
-  fhnTracksJetVsPtMCK0 = new TH2F("hnTracksJetVsPtMCK0","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtMCK0->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtMCK0->SetYTitle("nTracks");
-  fhnTracksJetVsPtMCK0->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtMCK0);
-
-
-  fhnTracksJetVsPtMCK0quarkEv = new TH2F("hnTracksJetVsPtMCK0quarkEv","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtMCK0quarkEv->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtMCK0quarkEv->SetYTitle("nTracks");
-  fhnTracksJetVsPtMCK0quarkEv->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtMCK0quarkEv);
-
-  fhRadiusJetVsPtMCK0 = new TH2F("hRadiusJetVsPtMCK0","",50,0,50,10,0.,1);
-  fhRadiusJetVsPtMCK0->SetXTitle("p_{T} (GeV/c)");
-  fhRadiusJetVsPtMCK0->SetYTitle("radius");
-  fhRadiusJetVsPtMCK0->Sumw2();
-  fListOfHistos->Add(fhRadiusJetVsPtMCK0);
-
-  fhnTracksJetVsPtAODK0 = new TH2F("hnTracksJetVsPtAODK0","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtAODK0->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtAODK0->SetYTitle("nTracks AODK0");
-  fhnTracksJetVsPtAODK0->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtAODK0);
-
-  fhnTracksJetVsPtAODK0quarkEv = new TH2F("hnTracksJetVsPtAODK0quarkEv","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtAODK0quarkEv->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtAODK0quarkEv->SetYTitle("nTracks AODK0quarkEv");
-  fhnTracksJetVsPtAODK0quarkEv->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtAODK0quarkEv);
-
-  fhRadiusJetVsPtAODK0 = new TH2F("hRadiusJetVsPtAODK0","",50,0,50,10,0.,1);
-  fhRadiusJetVsPtAODK0->SetXTitle("p_{T} (GeV/c)");
-  fhRadiusJetVsPtAODK0->SetYTitle("radius");
-  fhRadiusJetVsPtAODK0->Sumw2();
-  fListOfHistos->Add(fhRadiusJetVsPtAODK0);
-
-  fhnTracksJetVsPtAODpKch = new TH2F("hnTracksJetVsPtAODpKch","",50,0,50,20,-0.5,19.5);
-  fhnTracksJetVsPtAODpKch->SetXTitle("p_{T} (GeV/c)");
-  fhnTracksJetVsPtAODpKch->SetYTitle("nTracks AODpKch");
-  fhnTracksJetVsPtAODpKch->Sumw2();
-  fListOfHistos->Add(fhnTracksJetVsPtAODpKch);
-
-  fhRadiusJetVsPtAODpKch = new TH2F("hRadiusJetVsPtAODpKch","",50,0,50,20,-0.5,19.5);
-  fhRadiusJetVsPtAODpKch->SetXTitle("p_{T} (GeV/c)");
-  fhRadiusJetVsPtAODpKch->SetYTitle("Radius AODpKch");
-  fhRadiusJetVsPtAODpKch->Sumw2();
-  fListOfHistos->Add(fhRadiusJetVsPtAODpKch);
-
-  fhPythiaProcess      = CreatePythiaIDhisto("hPythiaProcess");
-  fListOfHistos->Add(fhPythiaProcess);  
-
-  fhPythiaProcessK0   = CreatePythiaIDhisto("hPythiaProcessK0");
-  fListOfHistos->Add(fhPythiaProcessK0);  
-
-  fhPythiaProcessKch   = CreatePythiaIDhisto("hPythiaProcessKch");
-  fListOfHistos->Add(fhPythiaProcessKch);  
-
-  fhPythiaProcessp    = CreatePythiaIDhisto("hPythiaProcessp");
-  fListOfHistos->Add(fhPythiaProcessp);  
-  
-  fhPythiaProcesspbar = CreatePythiaIDhisto("hPythiaProcesspbar");
-  fListOfHistos->Add(fhPythiaProcesspbar);  
-  fhdNdzJets5to10 = new TH1F("hdNdzJets5to10","",25,0,1.25);
-  fhdNdzJets5to10->SetXTitle("z");
-  fhdNdzJets5to10->SetYTitle("dN/dz");
-  fhdNdzJets5to10->Sumw2();
-  fListOfHistos->Add(fhdNdzJets5to10);
-
-  fhdNdzJets10to20 = new TH1F("hdNdzJets10to20","",25,0,1.25);
-  fhdNdzJets10to20->SetXTitle("z");
-  fhdNdzJets10to20->SetYTitle("dN/dz");
-  fhdNdzJets10to20->Sumw2();
-  fListOfHistos->Add(fhdNdzJets10to20);
-
-  fhdNdzJets20to30 = new TH1F("hdNdzJets20to30","",25,0,1.25);
-  fhdNdzJets20to30->SetXTitle("z");
-  fhdNdzJets20to30->SetYTitle("dN/dz");
-  fhdNdzJets20to30->Sumw2();
-  fListOfHistos->Add(fhdNdzJets20to30);
-
-  fhdNdzJets30to40 = new TH1F("hdNdzJets30to40","",25,0,1.25);
-  fhdNdzJets30to40->SetXTitle("z");
-  fhdNdzJets30to40->SetYTitle("dN/dz");
-  fhdNdzJets30to40->Sumw2();
-  fListOfHistos->Add(fhdNdzJets30to40);
-
-  fhdNdzJets40to60 = new TH1F("hdNdzJets40to60","",25,0,1.25);
-  fhdNdzJets40to60->SetXTitle("z");
-  fhdNdzJets40to60->SetYTitle("dN/dz");
-  fhdNdzJets40to60->Sumw2();
-  fListOfHistos->Add(fhdNdzJets40to60);
-  
-
-  fhdNdxiJets5to10 = new TH1F("hdNdxiJets5to10","",70,0,7);
-  fhdNdxiJets5to10->SetXTitle("z");
-  fhdNdxiJets5to10->SetYTitle("dN/dz");
-  fhdNdxiJets5to10->Sumw2();
-  fListOfHistos->Add(fhdNdxiJets5to10);
-
-  fhdNdxiJets10to20 = new TH1F("hdNdxiJets10to20","",70,0,7);
-  fhdNdxiJets10to20->SetXTitle("z");
-  fhdNdxiJets10to20->SetYTitle("dN/dz");
-  fhdNdxiJets10to20->Sumw2();
-  fListOfHistos->Add(fhdNdxiJets10to20);
-
-  fhdNdxiJets20to30 = new TH1F("hdNdxiJets20to30","",70,0,7);
-  fhdNdxiJets20to30->SetXTitle("z");
-  fhdNdxiJets20to30->SetYTitle("dN/dz");
-  fhdNdxiJets20to30->Sumw2();
-  fListOfHistos->Add(fhdNdxiJets20to30);
-
-  fhdNdxiJets30to40 = new TH1F("hdNdxiJets30to40","",70,0,7);
-  fhdNdxiJets30to40->SetXTitle("z");
-  fhdNdxiJets30to40->SetYTitle("dN/dz");
-  fhdNdxiJets30to40->Sumw2();
-  fListOfHistos->Add(fhdNdxiJets30to40);
-
-  fhdNdxiJets40to60 = new TH1F("hdNdxiJets40to60","",70,0,7);
-  fhdNdxiJets40to60->SetXTitle("z");
-  fhdNdxiJets40to60->SetYTitle("dN/dz");
-  fhdNdxiJets40to60->Sumw2();
-  fListOfHistos->Add(fhdNdxiJets40to60);
-
-  fhdNdptTracksJetPt5to10 = new TH1F("hdNdptTracksJetPt5to10","",250,0,25);
-  fhdNdptTracksJetPt5to10->SetXTitle("p_{T} (GeV)");
-  fhdNdptTracksJetPt5to10->SetYTitle("dN/dp_{T} 1/GeV");
-  fhdNdptTracksJetPt5to10->Sumw2();
-  fListOfHistos->Add(fhdNdptTracksJetPt5to10);
-
-  fhdNdptTracksJetPt10to20 = new TH1F("hdNdptTracksJetPt10to20","",25,0,25);
-  fhdNdptTracksJetPt10to20->SetXTitle("p_{T} (GeV)");
-  fhdNdptTracksJetPt10to20->SetYTitle("dN/dp_{T} 1/GeV");
-  fhdNdptTracksJetPt10to20->Sumw2();
-  fListOfHistos->Add(fhdNdptTracksJetPt10to20);
-
-  fhdNdptTracksJetPt20to30 = new TH1F("hdNdptTracksJetPt20to30","",25,0,25);
-  fhdNdptTracksJetPt20to30->SetXTitle("p_{T} (GeV)");
-  fhdNdptTracksJetPt20to30->SetYTitle("dN/dp_{T} 1/GeV");
-  fhdNdptTracksJetPt20to30->Sumw2();
-  fListOfHistos->Add(fhdNdptTracksJetPt20to30);
-
-  fhdNdptTracksJetPt30to40 = new TH1F("hdNdptTracksJetPt30to40","",25,0,25);
-  fhdNdptTracksJetPt30to40->SetXTitle("p_{T} (GeV)");
-  fhdNdptTracksJetPt30to40->SetYTitle("dN/dp_{T} 1/GeV");
-  fhdNdptTracksJetPt30to40->Sumw2();
-  fListOfHistos->Add(fhdNdptTracksJetPt30to40);
-
-  fhdNdptTracksJetPt40to60 = new TH1F("hdNdptTracksJetPt40to60","",25,0,25);
-  fhdNdptTracksJetPt40to60->SetXTitle("p_{T} (GeV)");
-  fhdNdptTracksJetPt40to60->SetYTitle("dN/dp_{T} 1/GeV");
-  fhdNdptTracksJetPt40to60->Sumw2();
-  fListOfHistos->Add(fhdNdptTracksJetPt40to60);
-
-
-  fh1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1); 
-  fh1Xsec->SetXTitle("<#sigma>");
-  fh1Xsec->Sumw2();
-  fListOfHistos->Add( fh1Xsec ); 
-  
-  fh1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
-  fh1Trials->SetXTitle("#sum{ntrials}");
-  fh1Trials->Sumw2();
-  fListOfHistos->Add( fh1Trials ); 
-//   fSettingsTree   = new TTree("JetChemAnalysisSettings","Analysis Settings");
-//   fSettingsTree->Branch("fUseLOConeJets",&fUseLOConeJets,"UseLOConeJets/O");
-//   fSettingsTree->Branch("fUseLOConeMCJets",&fUseLOConeMCJets,"UseLOConeMCJets/O");
-//   fSettingsTree->Branch("fUsePythiaJets",&fUsePythiaJets,"UsePythiaJets/O");
-//   fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
-//   fSettingsTree->Branch("fTrackPtCutJF", &fTrackPtCutJF,"TrackPt/D");
-//   fSettingsTree->Branch("fFilterBitJF",&fFilterBitJF,"FilterBitJF/i");
-//   fSettingsTree->Branch("fRequireITSRefitJF",&fRequireITSRefitJF,"RequireITSRefitJF/O");
-//   fSettingsTree->Branch("fRejectK0TracksJF",&fRejectK0TracksJF,"RejectK0TracksJF/O");
-//   fSettingsTree->Branch("fJetPtCut",&fJetPtCut,"JetPtCut/D");
-//   fSettingsTree->Branch("fJetEtaCut",&fJetEtaCut,"JetEtaCut/D");
-//   fSettingsTree->Branch("fFilterBit",&fFilterBit,"FilterBit/i");
-//   fSettingsTree->Branch("fTrackPtCut",&fTrackPtCut,"TrackPtCut/D");
-//   fSettingsTree->Branch("fTrackEtaCut",&fTrackEtaCut,"TrackEtaCut/D");
-//   fSettingsTree->Branch("fUseOnFlyV0s",&fUseOnFlyV0s,"UseOnFlyV0s/O");
-//   fSettingsTree->Branch("fCutnSigdEdx",&fCutnSigdEdx,"CutnSigdEdx/D");
-//   fSettingsTree->Branch("fUseAODMCTracksForUE",&fUseAODMCTracksForUE,"UseAODMCTracksForUE/O");
-//   fSettingsTree->Branch("fAreaReg",&fAreaReg,"AreaReg/D");
-//   fSettingsTree->Branch("fAvgTrials",&fAvgTrials,"AvgTrials/D");
-  // fListOfHistos->Add(fSettingsTree);
-  
-}
-
-//____________________________________________________________________
-void  AliAnalysisTaskJetChem::Terminate(Option_t */*option*/){
-
-  // Terminate analysis
-  //
-  
-  WriteSettings();
-  
-  // Update pointers reading them from the output slot
-  fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
-  if( !fListOfHistos ) {
-    AliError("Histogram List is not available");
-    return;
-  }
-
-  fhLeadingPt          = (TH1F*)fListOfHistos->FindObject("hLeadingPt");
-  fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->FindObject("hRegionSumPtMaxVsEt");
-  fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->FindObject("hRegionSumPtMinVsEt");
-  fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->FindObject("hRegionMultMaxVsEt");
-  fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->FindObject("hRegionMultMinVsEt");
-    
-   
-  fhRegionSumPtMaxVsEt->Divide(fhLeadingPt);
-  fhRegionSumPtMinVsEt->Divide(fhLeadingPt);
-  fhRegionMultMaxVsEt->Divide(fhLeadingPt);
-  fhRegionMultMinVsEt->Divide(fhLeadingPt);
-
-    
-  //Get Normalization
-  fh1Xsec                = (TProfile*) fListOfHistos->FindObject("fh1Xsec");
-  fh1Trials              = (TH1F*)fListOfHistos->FindObject("fh1Trials");
-  
-  //std::cout<<" fh1Xsec "<<fh1Xsec<<"  fh1Trials "<<fh1Trials<<std::endl;
-
-  if(fh1Xsec && fh1Trials){
-
-    Double_t xsec = fh1Xsec->GetBinContent(1);
-    Double_t ntrials = fh1Trials->GetBinContent(1);
-    Double_t normFactor = xsec/ntrials;
-    Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
-    fhLeadingPt->Scale(normFactor);
-  }
-
-  
-  if(fDebug > 1) AliInfo("End analysis");
-}
-
-// ----------------------------------------------------------------------------
-
-void  AliAnalysisTaskJetChem::WriteSettings(){ 
-
-  // fill settings tree
-  // not on GRID !
-
-  //fSettingsTree->Fill(); 
-
-}
-
-// ---------------------------------------------------------------------------
+  fTracksRecCuts->Clear();
+  fJetsRecCuts->Clear();
+  fListK0s->Clear();
 
-Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const {
-
-  // K0 mass ?  Use STAR params for bin counting & mass fit 
-
-  const Double_t massK0    = 0.497;        // from fits
-  const Double_t sigmaK0   = 0.0046;
-  const Double_t nSigmaSignal = 3.5; // STAR parameters for bin counting
-
-  if((massK0-nSigmaSignal*sigmaK0)<=mass &&  
-     mass<=(massK0 + nSigmaSignal*sigmaK0)) return kTRUE;
-
-  return kFALSE;
-}
-
-// ---------------------------------------------------------------------------
-
-Bool_t AliAnalysisTaskJetChem::IsLambdaInvMass(const Double_t mass) const{
-
-  // Lambda mass ?  
-
-
-  if(1.1<mass && mass<1.13) return kTRUE; // FIXME - adjust range from fit
-  return kFALSE;
-}
-
-// ---------------------------------------------------------------------------
-
-Bool_t AliAnalysisTaskJetChem::IsAcceptedDCAK0(/*const Double_t dca*/) const{
-
-  // DCA cut 
-
-  return kTRUE; // FIXME - adjust cut
+  //Post output data.
+  PostData(1, fCommonHistList);    
 }
 
-
-// ---------------------------------------------------------------------------
-
-Bool_t AliAnalysisTaskJetChem::IsAcceptedDCALambda(/*const Double_t dca*/) const {
-
-  // DCA cut
-
-  return kTRUE; // FIXME - adjust cut
+// ____________________________________________________________________________________________
+void AliAnalysisTaskJetChem::SetProperties(TH3F* 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);
 }
 
-// ---------------------------------------------------------------------------
-
+// ____________________________________________________________________________________________________________________________________________
 Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
 
   // apply TPC dE/dx cut similar as in AliTPCpidESD 
@@ -2729,50 +1015,59 @@ Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t
   return kTRUE;
 }
 
-// ----------------------------------------------------------------------------
+//___________________________________________________________________
+Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
+{
+  // K0 mass ? Use FF histo limits
+  
+  if(fFFIMInvMMin <= mass && mass <= fFFIMInvMMax) return kTRUE;
+
+  return kFALSE;
+}
 
-void AliAnalysisTaskJetChem::CheckV0s(AliAODJet* jetVect, Int_t maxPtRegionIndex,Bool_t& foundK0){
+//_____________________________________________________________________________________
+Int_t AliAnalysisTaskJetChem::GetListOfK0s(TList *list, const Int_t type)
+{
+  // fill list of V0s selected according to type
 
-  // loop over AOD V0s, fill masses etc
+  if(!list){
+    if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
+    return -1;
+  }
 
-  Int_t nV0 = fAOD->GetNumberOfV0s();
-  fhNV0s->Fill(nV0);  
+  if(type==kTrackUndef) return 0;
 
   for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
     
     AliAODv0* v0 = fAOD->GetV0(i);
 
     Bool_t isOnFly = v0->GetOnFlyStatus();
-    isOnFly ? fhV0onFly->Fill(1) : fhV0onFly->Fill(0);
     
-    if((fUseOnFlyV0s && isOnFly)  || (!fUseOnFlyV0s && !isOnFly) ){ // 'offline' V0s :  using vertex tracks, on-fly: take decay vertex into account during mom fit 
-      
-      Double_t massK0         = v0->MassK0Short();
-      Double_t massLambda     = v0->MassLambda();
-      Double_t massAntiLambda = v0->MassAntiLambda();
-      Double_t radiusV0       = v0->RadiusV0();
+    if(!isOnFly &&  (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue; 
+    if( isOnFly &&  (type == kOffl  || type == kOfflPID  || type == kOffldEdx  || type == kOfflPrim))  continue; 
 
-      Double_t etaV0          = v0->Eta();
-      Double_t fDCAV0ToVertex = v0->DcaV0ToPrimVertex();
-      Double_t fDCADaughters  = v0->DcaV0Daughters();
+    Double_t massK0 = v0->MassK0Short();
 
+    if(!(IsK0InvMass(massK0))) continue; // moved invMass cut for HI - otherwise too slow
 
-      Double_t v0Mom[3];
-      v0->PxPyPz(v0Mom);
-      TVector3 v0MomVect(v0Mom);
+    if(type == kOnFlyPID || type == kOfflPID){
       
-      AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
-      AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));   
+      AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
+      AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
+     
+      // AOD pid - cuts strongly into signal
+
+      AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
+      AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
+       
+      if(!( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ) ) continue;
+    }
+   
+    if(type == kOnFlydEdx || type == kOffldEdx){
+
+      AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
+      AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
 
-      Double_t posMom[3], negMom[3];
-      trackPos->PxPyPz(posMom);
-      trackNeg->PxPyPz(negMom);
-      
-      TVector3 posMomVec(posMom);
-      TVector3 negMomVec(negMom);
-      
-      Double_t dROpanV0 = posMomVec.DeltaR(negMomVec);
-      
       AliAODPid*  aodPidPos = trackPos->GetDetPid();
       AliAODPid*  aodPidNeg = trackNeg->GetDetPid();
 
@@ -2781,1094 +1076,29 @@ void AliAnalysisTaskJetChem::CheckV0s(AliAODJet* jetVect, Int_t maxPtRegionIndex
 
       Double_t momPos  = trackPos->P();
       Double_t momNeg  = trackNeg->P();
-           
-      fhdEdxVsMomV0->Fill(momPos,dEdxPos);
-      fhdEdxVsMomV0->Fill(momNeg,dEdxNeg);
-      Double_t pV0       = TMath::Sqrt(v0->Ptot2V0());
-      Double_t ptV0      = TMath::Sqrt(v0->Pt2V0());
-     
-      fhV0InvMassK0->Fill(massK0);
-      fhV0InvMassLambda->Fill(massLambda);
-      fhV0InvMassAntiLambda->Fill(massAntiLambda);
-    
-      fhV0DCADaughters->Fill(fDCADaughters);
-      fhV0Radius->Fill(radiusV0);
-      fhV0DCAToVertex->Fill(fDCAV0ToVertex);
-      fhdNdptV0->Fill(ptV0);
-
-      fhV0PtVsInvMassK0->Fill(massK0,ptV0); 
-
-
-      // K0 signal before cuts 
 
-      if(IsK0InvMass(massK0)){
-       
-       fhdNdptK0->Fill(ptV0);
-       
-       fhPtVsEtaK0->Fill(etaV0,ptV0);
-
-       Double_t dRV0MC = AssociateV0MC(&v0MomVect,310); // K0
-       if(dRV0MC < 0) dRV0MC = 0.99;
-       fhdRV0MC->Fill(dRV0MC);
-      }
-      if(IsAcceptedDCAK0(/*fDCAV0ToVertex*/)){
-       
-       fhV0InvMassK0DCA->Fill(massK0);
-       if(IsK0InvMass(massK0)) fhdNdptK0DCA->Fill(ptV0);
-
-       
-       if(IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)) fhdEdxVsMomV0pidEdx->Fill(momPos,dEdxPos);
-       if(IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) fhdEdxVsMomV0pidEdx->Fill(-1*momNeg,dEdxNeg);
-
-       if(IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) && 
-          IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){ 
-         
-         fhV0InvMassK0DCAdEdx->Fill(massK0);
-         fhV0PtVsInvMassK0DCAdEdx->Fill(massK0,ptV0); 
-
-         if(IsK0InvMass(massK0)){
-           fhdNdptK0DCAdEdx->Fill(ptV0);
-           fhV0DCAToVertexK0->Fill(fDCAV0ToVertex);
-           fhdROpanK0VsPt->Fill(ptV0,dROpanV0);
-           foundK0 = kTRUE;
-         }
-       }
-       
-
-       // AOD pid - cuts strongly into signal
-
-       AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
-       AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
-       
-       if(mpPIDNeg == AliAODTrack::kPion) fhdEdxVsMomV0piPID->Fill(momPos,dEdxPos);
-       if(mpPIDPos == AliAODTrack::kPion) fhdEdxVsMomV0piPID->Fill(-1*momNeg,dEdxNeg);
-
-
-       if( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ){
-
-         fhV0InvMassK0DCAPID->Fill(massK0);
-       }       
-      }
+      Int_t cutnSigdEdx = 2;
+      if(! (IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,cutnSigdEdx)) && (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,cutnSigdEdx)) ) continue;
       
-
-      if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
-       
-       fhV0InvMassK0Lambda->Fill(massK0);      
-      }
-
-
-      // Lambda 
-
-      if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-        IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-        IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
+    }   
     
-       fhV0InvMassLambdaDCAdEdx->Fill(massLambda);
-      }
+    if(type == kOnFlyPrim || type == kOfflPrim){
       
-      if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-        IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-        IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
-    
-       fhV0InvMassAntiLambdaDCAdEdx->Fill(massAntiLambda);
-      }
-
-
-      // jet events 
-
-      if(jetVect){
-       Int_t regionV0vect = IsTrackInsideRegion(jetVect,&v0MomVect);
-
-       // calc xi 
-       Double_t jetpt     = jetVect->Pt();
-       Double_t dPhiJetV0 = (jetVect->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
-       Double_t dRJetV0   = (jetVect->MomentumVector()->Vect()).DeltaR(v0MomVect);             
-
-       Double_t z         = pV0/jetpt;
-       Double_t xi        = TMath::Log(1/z);
-
-       fhdPhiJetV0->Fill(dPhiJetV0);
-
-       // K0
-
-       if(IsAcceptedDCAK0(/*fDCAV0ToVertex*/) && 
-          IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) && 
-          IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-         
-         fhV0InvMassK0JetEvt->Fill(massK0);
-
-         if(IsK0InvMass(massK0)){
-           fhdNdptK0JetEvt->Fill(ptV0);
-           fhdNdzK0->Fill(z); 
-           fhdNdxiK0->Fill(xi);
-
-           fhdPhiJetK0->Fill(dPhiJetV0);
-           fhdRJetK0->Fill(dRJetV0);
-
-           if(5<jetpt && jetpt<=10)       fhdNdzK05to10->Fill(z);
-           else if(10<jetpt && jetpt<=20) fhdNdzK010to20->Fill(z);
-           else if(20<jetpt && jetpt<=30) fhdNdzK020to30->Fill(z);
-           else if(30<jetpt && jetpt<=40) fhdNdzK030to40->Fill(z);
-           else if(40<jetpt && jetpt<=60) fhdNdzK040to60->Fill(z);
-         }
-       }
-       
-
-       // Lambda
-       
-       if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-          IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-          IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-         
-         fhV0InvMassLambdaJetEvt->Fill(massLambda);
-
-         if(IsLambdaInvMass(massLambda)){
-           fhdNdzLambda->Fill(z); 
-         }
-       }
-       
-       if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-          IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-          IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
-         
-         fhV0InvMassAntiLambdaJetEvt->Fill(massAntiLambda);
-
-         if(IsLambdaInvMass(massAntiLambda)){
-           fhdNdzAntiLambda->Fill(z); 
-         }
-       }
-       
-       // fill histos max region 
-       
-       if(regionV0vect != 0 && regionV0vect == maxPtRegionIndex){ // max region
-       
-         // K0
-
-         if(IsAcceptedDCAK0(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){ // K0 cuts
-           
-           fhV0InvMassK0Max->Fill(massK0);
-           
-           if(IsK0InvMass(massK0)){
-             fhdNdzK0Max->Fill(z);
-             fhdNdxiK0Max->Fill(xi);
-             fhdNdptK0Max->Fill(ptV0);
-           }
-         }
-
-         // Lambda
-         
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-           
-           fhV0InvMassLambdaMax->Fill(massLambda);
-         }
-
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
-           
-           fhV0InvMassAntiLambdaMax->Fill(massAntiLambda);
-         }
-         
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            ((IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
-             (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)))){
-              
-           if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
-             fhdNdzLambdaMax->Fill(z);
-             fhdNdxiLambdaMax->Fill(xi);
-             fhdNdptLambdaMax->Fill(ptV0);
-           }
-         }
-       }
-       
-       // fill histos min region
-
-       if(regionV0vect != 0 && regionV0vect != maxPtRegionIndex){ // min region 
-       
-         // K0
-  
-         if(IsAcceptedDCAK0(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-           
-           fhV0InvMassK0Min->Fill(massK0);
-           
-           if(IsK0InvMass(massK0)){
-             fhdNdzK0Min->Fill(z);
-             fhdNdxiK0Min->Fill(xi);
-             fhdNdptK0Min->Fill(ptV0);
-           }
-         }
-       
-
-         // Lambda
-         
-         
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-           
-           fhV0InvMassLambdaMin->Fill(massLambda);
-         }
-
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-            IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
-           
-           fhV0InvMassAntiLambdaMin->Fill(massAntiLambda);
-         }
-         
-         if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-            ((IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
-             (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)))){
-           
-           if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
-             fhdNdzLambdaMin->Fill(z);
-             fhdNdxiLambdaMin->Fill(xi);
-             fhdNdptLambdaMin->Fill(ptV0);
-           }
-         }
-       }
-       
-
-       // jet region 
-
-       if(regionV0vect == 0){ // jet region
-       
-         //Double_t dRJetV0 = jetVect->DeltaR(v0MomVect); 
-         
-         if(dRJetV0 <= fConeRadius){ 
-         
-           // fill histos jet region  
-
-           // K0
-
-           if(IsAcceptedDCAK0(/*fDCAV0ToVertex*/) && 
-              IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx) && 
-              IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){ // K0 cuts
-             
-             fhV0InvMassK0Jet->Fill(massK0);
-           
-             if(IsK0InvMass(massK0)){
-               fhdNdzK0Jet->Fill(z);
-               fhdNdxiK0Jet->Fill(xi);
-               fhdNdptK0Jet->Fill(ptV0);
-             }
-           }
-
-
-           // Lambda 
-         
-           if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-              IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)){
-           
-             fhV0InvMassLambdaJet->Fill(massLambda);
-           }     
-           if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-              IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-              IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)){
-             
-             fhV0InvMassAntiLambdaJet->Fill(massAntiLambda);
-           }
-           
-           if(IsAcceptedDCALambda(/*fDCAV0ToVertex*/) && 
-              ((IsAccepteddEdx(momPos,dEdxPos,AliPID::kProton,fCutnSigdEdx) && 
-                IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,fCutnSigdEdx)) ||
-               (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kProton,fCutnSigdEdx) && 
-                IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,fCutnSigdEdx)))){
-             
-             if(IsLambdaInvMass(massLambda) || (IsLambdaInvMass(massAntiLambda))){
-               fhdNdzLambdaJet->Fill(z);
-               fhdNdxiLambdaJet->Fill(xi);
-               fhdNdptLambdaJet->Fill(ptV0);
-             }
-           }
-         }
-       }
-      }
-    }
-  }
-}
-
-// ----------------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::CheckMCParticles(AliAODJet* jetVect,Int_t maxPtRegionIndex,Bool_t& isK0Event){
-
-  // histos for generated MC  
-
-  TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-       
-  if(!farray){
-    AliInfo("no mcparticles branch"); 
-    return;
-  }
-
-  Int_t ntrks = farray->GetEntries();
-  if (fDebug>1) AliInfo(Form("check MC particles, tracks %d \n",ntrks));
-
-  Int_t pythiaPID = GetPythiaProcessID();
-
-  isK0Event          = kFALSE;
-  Bool_t ispEvent    = kFALSE;
-  Bool_t ispBarEvent = kFALSE;
-  Bool_t isKchEvent  = kFALSE;
-
-  Bool_t isQuarkHardScatteringEvent = IsQuarkHardScatteringEvent(pythiaPID);
-  Bool_t isGluonHardScatteringEvent = IsGluonHardScatteringEvent(pythiaPID);
-
-  for(Int_t i =0 ; i<ntrks; i++){ // mc tracks loop
-
-    AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
-    
-    Double_t trackPt = mctrk->Pt(); 
-
-    //Cuts
-    if (!(mctrk->IsPhysicalPrimary())) continue;
-    if ((trackPt < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut)) continue; 
-
-    // OB PID histo
-    Int_t pdg = mctrk->GetPdgCode();
-    
-
-    FillPIDhisto(fhPIDMC,pdg);
-    if(isQuarkHardScatteringEvent) FillPIDhisto(fhPIDMC_quarkEv,pdg);
-    if(isGluonHardScatteringEvent) FillPIDhisto(fhPIDMC_gluonEv,pdg);
-
-    if(pdg == 111)              fhdNdptpi0MC->Fill(trackPt);
-    if(pdg == 22)               fhdNdptgammaMC->Fill(trackPt);
-    if(TMath::Abs(pdg) == 211)  fhdNdptchPiMC->Fill(trackPt);
-    if(pdg == 310)              fhdNdptK0MC->Fill(trackPt);
-    if(TMath::Abs(pdg) == 321)  fhdNdptchKMC->Fill(trackPt);
-    if(pdg == 2212)             fhdNdptpMC->Fill(trackPt);
-    if(pdg == -2212)            fhdNdptpBarMC->Fill(trackPt);
-    if(pdg == 3122)             fhdNdptLambdaMC->Fill(trackPt);
-    if(pdg == -3122)            fhdNdptLambdaBarMC->Fill(trackPt);
-    if(pdg == 3332)             fhdNdptOmegaMC->Fill(trackPt);
-    if(pdg == -3332)            fhdNdptOmegaBarMC->Fill(trackPt);
-    
-    if(pdg == 310)             isK0Event   = kTRUE;
-    if(TMath::Abs(pdg) == 321) isKchEvent  = kTRUE;
-    if(pdg == 2212)            ispEvent    = kTRUE;
-    if(pdg == -2212)           ispBarEvent = kTRUE;
-    
-    
-    Int_t pdgMotherChK = -1;
-    Int_t pdgMotherK0  = -1;
-    Int_t pdgGrandMotherK0 = -1;
-
-    if(TMath::Abs(pdg) == 321){ // chK
-      Int_t labelMother = mctrk->GetMother();
-      AliAODMCParticle* mctrkMother    = NULL;
-      
-      for(Int_t k=0 ; k<ntrks; k++){
-      
-       AliAODMCParticle* mctrk2 = (AliAODMCParticle*)farray->At(k);
-       if(mctrk2->GetLabel() == labelMother) mctrkMother = mctrk2;
-      }
-      
-      if(mctrkMother) pdgMotherChK = mctrkMother->GetPdgCode();
-      FillPIDhisto(fhPIDMCMotherChK,pdgMotherChK);
-
-      //printf("pdgMotherChK %d \n",pdgMotherChK);
-    }
-    
-    if(pdg == 310){ // K0
-      Int_t labelMother = mctrk->GetMother();
-      AliAODMCParticle* mctrkMother = NULL;
-      
-      for(Int_t k=0 ; k<ntrks; k++){
-       AliAODMCParticle* mctrk2 = (AliAODMCParticle*)farray->At(k);
-       if(mctrk2->GetLabel() == labelMother) mctrkMother = mctrk2;
-      }
-      
-      if(mctrkMother) pdgMotherK0 = mctrkMother->GetPdgCode();
-      FillPIDhisto(fhPIDMCMotherK0,pdgMotherK0);
-
-      Int_t labelGrandMother = -1; 
-      if(mctrkMother) mctrkMother->GetMother();
-      AliAODMCParticle* mctrkGrandMother = NULL;
-      
-      for(Int_t k=0 ; k<ntrks; k++){
-       AliAODMCParticle* mctrk2 = (AliAODMCParticle*)farray->At(k);
-       if(mctrk2->GetLabel() == labelGrandMother) mctrkGrandMother = mctrk2;
-      }
-      
-      if(mctrkGrandMother) pdgGrandMotherK0 = mctrkGrandMother->GetPdgCode();
-      FillPIDhisto(fhPIDMCGrandMotherK0,pdgGrandMotherK0);
-    }
-    
-    if(jetVect){ // jet event 
-
-      FillPIDhisto(fhPIDMCAll,pdg);
-    
-      TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
-      
-      Int_t region = IsTrackInsideRegion(jetVect, &partVect );  
-      
-      //Double_t deltaPhi = jetVect->DeltaPhi(partVect); //+k270rad;
-      Double_t deltaPhi = (jetVect->MomentumVector()->Vect()).DeltaPhi(partVect);
-      if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
-      //Double_t deltaR = jetVect->DeltaR(partVect); //+k270rad;
-      Double_t deltaR = (jetVect->MomentumVector()->Vect()).DeltaR(partVect);
-
-      // calc xi 
-      Double_t jetpt = jetVect->Pt(); 
-      //Double_t pV0   = partVect.Mag(); 
-      //Double_t ptV0  = partVect.Pt();
-   
-
-      Double_t z  = trackPt / jetpt;
-      Double_t xi = TMath::Log(1/z);
-
-      if(!(mctrk->Charge() == 0 || mctrk->Charge()==-99)) fhdNdxiMC->Fill(xi);
-
-      if(pdg == 310){ // K0
-       fhdPhiJetK0MC->Fill(deltaPhi);
-       fhdRJetK0MC->Fill(deltaR);
-       fhdNdxiK0MC->Fill(xi);
-       fhdNdzK0MC->Fill(z);
-       fhdNdptK0MCJetEvt->Fill(trackPt);
-      }
-
+      AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
+      AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
       
-      if(region != 0 && region == maxPtRegionIndex){ // max region
-       
-       if(TMath::Abs(pdg) == 211)  fhdNdptchPiMCMax->Fill(trackPt);
-       if(pdg == 310)              fhdNdptK0MCMax->Fill(trackPt);
-       if(TMath::Abs(pdg) == 321)  fhdNdptchKMCMax->Fill(trackPt);
-       if(pdg == 2212)             fhdNdptpMCMax->Fill(trackPt);
-       if(pdg == -2212)            fhdNdptpBarMCMax->Fill(trackPt);
-       if(pdg == 3122)             fhdNdptLambdaMCMax->Fill(trackPt);
-       if(pdg == -3122)            fhdNdptLambdaBarMCMax->Fill(trackPt);
-      }
+      //std::cout<<"  filer map trackPos "<<trackPos->GetFilterMap()<<" trackNeg "<<trackNeg->GetFilterMap()<<std::endl;
       
-      if(region != 0 && region != maxPtRegionIndex){ // min region
-       
-       FillPIDhisto(fhPIDMCMin,pdg);
-       
-       if(TMath::Abs(pdg) == 211)  fhdNdptchPiMCMin->Fill(trackPt);
-       if(pdg == 310)              fhdNdptK0MCMin->Fill(trackPt);
-       if(TMath::Abs(pdg) == 321)  fhdNdptchKMCMin->Fill(trackPt);
-       if(pdg == 2212)             fhdNdptpMCMin->Fill(trackPt);
-       if(pdg == -2212)            fhdNdptpBarMCMin->Fill(trackPt);
-       if(pdg == 3122)             fhdNdptLambdaMCMin->Fill(trackPt);
-       if(pdg == -3122)            fhdNdptLambdaBarMCMin->Fill(trackPt);
-
-       if(pdg == 3332)  fhdNdptOmegaMCMin->Fill(trackPt);
-       if(pdg == -3332) fhdNdptOmegaBarMCMin->Fill(trackPt);
-      }
-
-      // trans region
-      if(region != 0){
-       if(TMath::Abs(pdg) == 321) FillPIDhisto(fhPIDMCMotherChKTrans,pdgMotherChK);
-       if(pdg == 310){
-         FillPIDhisto(fhPIDMCMotherK0Trans,pdgMotherK0);
-         FillPIDhisto(fhPIDMCGrandMotherK0Trans,pdgGrandMotherK0);
-       }
-      }
-   
-      if(region == 0){ //  jet region ?
-       
-       FillPIDhisto(fhPIDMCJet,pdg);
-       
-       //Double_t dRJetV0 = jetVect->DeltaR(partVect); 
-       Double_t dRJetV0 =  (jetVect->MomentumVector()->Vect()).DeltaR(partVect);
-       
-       if(dRJetV0 <= fConeRadius){ 
-         
-         if(pdg == 310){ // K0
-           
-           fhdNdptK0MCJet->Fill(trackPt);
-           fhdNdxiK0MCJet->Fill(xi);
-           fhdNdzK0MCJet->Fill(z);
-         }
-       
-         if(TMath::Abs(pdg) == 211)  fhdNdptchPiMCJet->Fill(trackPt);
-         if(TMath::Abs(pdg) == 321)  fhdNdptchKMCJet->Fill(trackPt);
-         if(pdg == 2212)             fhdNdptpMCJet->Fill(trackPt);
-         if(pdg == -2212)            fhdNdptpBarMCJet->Fill(trackPt);
-         if(pdg == 3122)             fhdNdptLambdaMCJet->Fill(trackPt);
-         if(pdg == -3122)            fhdNdptLambdaBarMCJet->Fill(trackPt);  
-       }
-      }
+      if((fFilterMaskK0>0) && !(trackPos->TestFilterBit(fFilterMaskK0)))   continue;
+      if((fFilterMaskK0>0) && !(trackNeg->TestFilterBit(fFilterMaskK0)))   continue;
     }
-  }
-
-
-  FillPythiaIDhisto(fhPythiaProcess,pythiaPID);
-  if(isK0Event)   FillPythiaIDhisto(fhPythiaProcessK0,pythiaPID);
-  if(isKchEvent)  FillPythiaIDhisto(fhPythiaProcessKch,pythiaPID);
-  if(ispEvent)    FillPythiaIDhisto(fhPythiaProcessp,pythiaPID);
-  if(ispBarEvent) FillPythiaIDhisto(fhPythiaProcesspbar,pythiaPID);
-
-}
-       
-// ----------------------------------------------------------------------------
-
-Double_t AliAnalysisTaskJetChem::AssociateV0MC(const TVector3* V0Mom,const Int_t pdg){
-
-  // find closest MC gen. particle for V0 vector
-
-  TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-       
-  if(!farray){
-    AliInfo("no mcparticles branch"); 
-    return -1;
-  }
-
-  Double_t dRmin = -1;
-
-  Int_t ntrks = farray->GetEntries();
-
-  for(Int_t i =0 ; i<ntrks; i++){
-
-    AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
-    
-    //Cuts
-    if (!(mctrk->IsPhysicalPrimary())) continue;
-
-    Int_t pdgtrk = mctrk->GetPdgCode();
     
-    if(pdgtrk != pdg) continue;
-
-    TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
-
-    Double_t dR = V0Mom->DeltaR(partVect);
-
-    if(dRmin<0) dRmin = dR; // initialize
-    
-    if(dR < dRmin) dRmin = dR;
-      
+    list->Add(v0);
   }
   
-  return dRmin;
-}
-
-// ----------------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::CompLeadingJets(AliAODJet* jetLeadingAOD,AliAODJet* jetLeadingMC,const Int_t pythiaPID,
-                                            const Bool_t foundK0AOD,const Bool_t foundK0MC){
-
-  // leading jet properties
-
-  Double_t ptLeadingAOD      = -1;
-  Double_t etaLeadingAOD     = -1;
-  Double_t phiLeadingAOD     = -1;
-  Int_t    nTracksLeadingAOD = -1;
-
-  Double_t ptLeadingMC       = -1; 
-  Double_t etaLeadingMC      = -1;
-  Double_t phiLeadingMC      = -1;
-  Int_t    nTracksLeadingMC  = -1;
-
-  if(jetLeadingAOD){
-    if(jetLeadingAOD->Pt()>fJetPtCut && TMath::Abs(jetLeadingAOD->Eta())<fJetEtaCut){
-
-      ptLeadingAOD      = jetLeadingAOD->Pt();
-      etaLeadingAOD     = jetLeadingAOD->Eta();
-      phiLeadingAOD     = jetLeadingAOD->Phi();
-      nTracksLeadingAOD = jetLeadingAOD->GetRefTracks()->GetEntriesFast();
-
-      Double_t radiusAOD = GetJetRadius(jetLeadingAOD,0.8);
-      fhnTracksJetVsPtAOD->Fill(ptLeadingAOD,nTracksLeadingAOD);
-      fhRadiusJetVsPtAOD->Fill(ptLeadingAOD,radiusAOD);
-      if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtAODquarkEv->Fill(ptLeadingAOD,nTracksLeadingAOD);
-      
-      if(foundK0AOD){
-       
-       fhnTracksJetVsPtAODK0->Fill(ptLeadingAOD,nTracksLeadingAOD);
-       if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtAODK0quarkEv->Fill(ptLeadingAOD,nTracksLeadingAOD);
-       fhRadiusJetVsPtAODK0->Fill(ptLeadingAOD,radiusAOD);
-      }
-      
-
-      // check if p/Kch in jet 
-
-      Bool_t foundpKch = kFALSE; 
-      Int_t nTracksJet = jetLeadingAOD->GetRefTracks()->GetEntriesFast();
-
-      for(int i=0; i<nTracksJet; i++){
-
-       AliAODTrack* track = (AliAODTrack*) jetLeadingAOD->GetRefTracks()->At(i);
-
-       Double_t mom  = track->P();
-       
-       AliAODPid* aodPid = track->GetDetPid();
-       Double_t   dEdx   = aodPid->GetTPCsignal();
-
-       if(IsAccepteddEdx(mom,dEdx,AliPID::kKaon,fCutnSigdEdx) ||  
-          IsAccepteddEdx(mom,dEdx,AliPID::kProton,fCutnSigdEdx)){
-
-         foundpKch = kTRUE; 
-       }
-      } // track loop
-      
-      if(foundpKch){
-       fhnTracksJetVsPtAODpKch->Fill(ptLeadingAOD,nTracksLeadingAOD);
-       fhRadiusJetVsPtAODpKch->Fill(ptLeadingAOD,radiusAOD);
-      }
-    }
-  }
-
-
-  if(jetLeadingMC){
-    if(jetLeadingMC->Pt()>fJetPtCut && TMath::Abs(jetLeadingMC->Eta())<fJetEtaCut){
-
-      ptLeadingMC      = jetLeadingMC->Pt();
-      etaLeadingMC     = jetLeadingMC->Eta();
-      phiLeadingMC     = jetLeadingMC->Phi();
-      nTracksLeadingMC = jetLeadingMC->GetRefTracks()->GetEntriesFast();    
-      
-      Double_t radiusMC  = GetJetRadius(jetLeadingMC,0.8);
-      fhnTracksJetVsPtMC->Fill(ptLeadingMC,nTracksLeadingMC);
-      fhRadiusJetVsPtMC->Fill(ptLeadingMC,radiusMC);
-      if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtMCquarkEv->Fill(ptLeadingMC,nTracksLeadingMC);
-      
-      if(foundK0MC){
-       
-       fhnTracksJetVsPtMCK0->Fill(ptLeadingMC,nTracksLeadingMC);
-       if(IsQuarkHardScatteringEvent(pythiaPID)) fhnTracksJetVsPtMCK0quarkEv->Fill(ptLeadingMC,nTracksLeadingMC);
-       fhRadiusJetVsPtMCK0->Fill(ptLeadingMC,radiusMC);
-      }
-    }
-  }
-  
-  if(jetLeadingAOD && jetLeadingMC){
-    
-    //std::cout<<" comp: leading jetPt AOD "<<ptLeadingAOD<<" MC "<<ptLeadingMC<<std::endl;
-    //if(jetLeadingAOD && jetLeadingMC) 
-    //std::cout<<" leading jet eta AOD "<<jetLeadingAOD->Eta()<<" MC "<<jetLeadingMC->Eta()<<std::endl;
-
-    
-    if(jetLeadingMC->Pt()>fJetPtCut && jetLeadingAOD->Pt()>fJetPtCut && 
-       TMath::Abs(jetLeadingMC->Eta())<fJetEtaCut && TMath::Abs(jetLeadingAOD->Eta())<fJetEtaCut){
-      
-      fhLeadingPtAODvsMC->Fill(ptLeadingMC,ptLeadingAOD);
-      fhLeadingEtaAODvsMC->Fill(etaLeadingMC,etaLeadingAOD);
-      fhLeadingPhiAODvsMC->Fill(phiLeadingMC,phiLeadingAOD);
-      fhnTracksLeadingAODvsMC->Fill(nTracksLeadingMC,nTracksLeadingAOD);
-      
-      TLorentzVector *mom4MC  = jetLeadingMC->MomentumVector();
-      TLorentzVector *mom4AOD = jetLeadingAOD->MomentumVector();
-      
-      Double_t dR = mom4MC->DeltaR(*mom4AOD);
-      fhLeadingdRAODMC->Fill(dR);
-      
-      if(dR<0.4) fhLeadingPtAODvsMCdRcut->Fill(ptLeadingMC,ptLeadingAOD);
-      
-      Double_t dPt = ptLeadingAOD - ptLeadingMC;
-      Double_t dnTracks = nTracksLeadingAOD - nTracksLeadingMC;
-      
-      fhdnTracksVsdPtLeadingAODMC->Fill(dPt,dnTracks);  
-    }
-  }
+  Int_t nK0s = list->GetSize();
   
+  return nK0s;
 }
 
-// ---------------------------------------------------------------------------
-
-// void AliAnalysisTaskJetChem::CheckK0(){
-
-//   TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-  
-//   if(!farray){
-//     AliInfo("no mcparticles branch"); 
-//     return;
-//   }
-  
-  
-//   Int_t ntrks = farray->GetEntries();
-    
-//   for(Int_t i=0; i<ntrks; i++){ // trk loop
-    
-//     AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
-    
-//     Int_t pdg = mctrk->GetPdgCode();
-    
-//     if(pdg != 310) continue; // K0
-
-//     cout<<" MC K0: physPrim "<<mctrk->IsPhysicalPrimary()<<endl;
-
-//     Int_t labelMother = mctrk->GetMother();
-
-//     cout<<" found K00, label mother "<<labelMother<<endl;
-
-//     AliAODMCParticle* mctrkMother = NULL;
-//     Int_t pdgMother = -1;
-    
-//     for(Int_t k=0 ; k<ntrks; k++){
-      
-//       mctrkMother = (AliAODMCParticle*)farray->At(k);
-//       if(mctrkMother->GetLabel() == labelMother) break;
-//     }
-    
-//     pdgMother = mctrkMother->GetPdgCode();
-//     cout<<" K0 mother pdg "<<pdgMother<<" GID "<<fpdgdb->ConvertPdgToGeant3(pdgMother)<<" isPrimary "<<mctrkMother->IsPrimary()<<endl;
-//     //cout<<" mother name "<<mctrkMother->GetName()<<endl;
-    
-//   }
-// }
-
-// ---------------------------------------------------------------------------
-
-// void AliAnalysisTaskJetChem::CheckK0Stack(){
-  
-//   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-//   if (!mcHandler) {
-//     Printf("ERROR: Could not retrieve MC event handler");
-//     return;
-//   }
-  
-//   AliMCEvent* mcEvent = mcHandler->MCEvent();
-//   if (!mcEvent) {
-//     Printf("ERROR: Could not retrieve MC event");
-//     return;
-//   }
-  
-//    AliStack* mcStack = mcEvent->Stack();//Load Stack
-//    AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
-
-//    Int_t nTracksMC = mcStack->GetNtrack();
-
-//    Bool_t foundK0 = kFALSE;
-
-//    for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
-//      //Cuts
-//      if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
-     
-//      TParticle* mctrk = mcStack->Particle(iTracks);
-        
-//      Int_t pdg = mctrk->GetPdgCode();
-
-//     if ((mctrk->Pt() < 2*fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue; 
-
-//      if(pdg == 310) foundK0 = kTRUE;
-
-//      if(pdg != 310) continue; // K0 short
-
-//      cout<<" check K0 "<<endl;
-
-//      Int_t indexMother = -999;
-//      TParticle* mctrkMother = mctrk;
-//      TParticle* mctrkFirstGenMother = NULL; // last mother which is not the primary proton
-
-//      Int_t nAncestors = 0;
-
-//      while(indexMother != -1){
-
-//        indexMother = mctrkMother->GetFirstMother();
-     
-//        if(indexMother != -1){
-//      mctrkFirstGenMother = mctrkMother;
-//      mctrkMother  =  mcStack->Particle(indexMother);
-//      nAncestors++;
-//        }
-
-//        cout<<" nAncestors "<<nAncestors<<" pdg mother "<<mctrkMother->GetPdgCode()<<" name "<<mctrkMother->GetName()<<endl;
-//      }
-//      cout<<" pdg firstGenMother "<<mctrkFirstGenMother->GetPdgCode()<<" name "<<mctrkFirstGenMother->GetName()<<endl;
-
-
-
-//      cout<<" pythiaGenHeader "<<pythiaGenHeader<<endl;
-//      cout<<" Pythia Process type "<<pythiaGenHeader->ProcessType()<<" ptHard "<<pythiaGenHeader->GetPtHard()<<endl;
-
-
-//      fhPythiaProcess->Fill(pythiaGenHeader->ProcessType());
-//      if(foundK0) fhPythiaProcess_K0->Fill(pythiaGenHeader->ProcessType());
-
-//      //Int_t indexGrandMother      = mctrkMother->GetFirstMother();
-//      //cout<<" indexGrandMother "<<indexGrandMother<<endl;
-     
-//      //if(indexGrandMother>-1){
-//      //  TParticle* mctrkGrandMother = mcStack->Particle(indexGrandMother);
-//      //  cout<<" pdg grandMother "<<mctrkGrandMother->GetPdgCode()<<" name "<<mctrkGrandMother->GetName()<<endl;
-//      // }
-
-//    }
-// }
-
-// ---------------------------------------------------------------------------------
-
-Bool_t  AliAnalysisTaskJetChem::IsQuarkHardScatteringEvent(const Int_t PID){
-
-  // Pythia Manual sec. 8.2.1 : 
-  // if Pythia PID = 92,93,94   event is diffractive 
-  // if Pythia PID = 13, 28, 68 hard scattering products are gg, qg 
-  // if Pythia PID = 11, 12, 53 hard scattering products are qq,  q\bar{q}  
-
-  if(PID == 92 || PID == 93 || PID == 94 || PID==95) return kFALSE;
-  else if(PID == 13 || PID == 28 || PID == 68) return kFALSE;
-  else if(PID == 11 || PID == 12 || PID == 53) return kTRUE;
-  else{
-    AliInfo(Form("unknown Pythia PID %d",PID));
-  }
-
-  return kFALSE;
-}
-
-
-// ---------------------------------------------------------------------------------
-
-Bool_t  AliAnalysisTaskJetChem::IsGluonHardScatteringEvent(const Int_t PID){
-
-  // Pythia Manual sec. 8.2.1 : 
-  // if Pythia PID = 92,93,94   event is diffractive 
-  // if Pythia PID = 95: low pt event (MPI) 
-  // if Pythia PID = 13, 28, 68 hard scattering products are gg, qg 
-  // if Pythia PID = 11, 12, 53 hard scattering products are qq,  q\bar{q}  
-
-
-  if(PID == 92 || PID == 93 || PID == 94 || PID == 95) return kFALSE;
-  else if(PID == 13 || PID == 68) return kTRUE;
-  else if(PID == 28) return kFALSE; // mixed gq final state
-  else if(PID == 11 || PID == 12 || PID == 53) return kFALSE;
-  else{
-    AliInfo(Form("unknown Pythia PID %d",PID));
-  }
-
-  return kFALSE;
-}
-
-// ---------------------------------------------------------------------------------
-
-Bool_t  AliAnalysisTaskJetChem::IsDiffractiveEvent(const Int_t PID){
-
-  // Pythia Manual sec. 8.2.1 : 
-  // if Pythia PID = 92,93,94   event is diffractive 
-  // if Pythia PID = 13, 28, 68 hard scattering products are gg, qg 
-  // if Pythia PID = 11, 12, 53 hard scattering products are qq,  q\bar{q}  
-
-  if(PID == -1) return kFALSE;
-
-  if(PID == 13 || PID == 28 || PID == 68) return kFALSE;
-  else if(PID == 11 || PID == 12 || PID == 53) return kFALSE;
-  else if(PID == 92 || PID == 93 || PID == 94 || PID==95) return kTRUE;
-  else{
-    AliInfo(Form("unknown Pythia PID %d",PID));
-  }
-
-  return kFALSE;
-
-}
-
-
-// ----------------------------------------------------------------------------------
-
-Int_t AliAnalysisTaskJetChem::GetPythiaProcessID(){
-
-  // Pythia PID for this event 
-  
-  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  if (!mcHandler) {
-    //Printf("ERROR: Could not retrieve MC event handler");
-    return -1;
-  }
-  
-  AliMCEvent* mcEvent = mcHandler->MCEvent();
-  if (!mcEvent) {
-    AliInfo("could not retrieve MC event");
-    return -1;
-  }
-  
-  AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
-  
-  if (!pythiaGenHeader) {
-    AliInfo("Could not retrieve pythiaEventHeader");
-    return -1;
-  }
-
-  Int_t  pid = pythiaGenHeader->ProcessType(); 
-
-  return pid;
-
-}
-
-// ----------------------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::GetJetTracksResum(TList* list, AliAODJet* jet, const Double_t radius){
-  
-
-  if(!jet) return; // no jet in acc in event
-
-  // list of AOD tracks in jet cone, using cone axis and distance axis-track (and not trackrefs)
-
-  Int_t nTracks = fAOD->GetNTracks();
-
-  if(!nTracks) return;
-
-  Double_t jetMom[3];
-  jet->PxPyPz(jetMom);
-  TVector3 jet3mom(jetMom);
-    
-  
-  // standard cuts + ITS refit = stdrd PWG4 cut (compose them for productions with old ESDFilter task) 
-
-  for (Int_t itrack=0; itrack<nTracks; itrack++) {
-
-    AliAODTrack* track = fAOD->GetTrack(itrack);
-    
-    UInt_t status = track->GetStatus();
-
-    if(!track->TestFilterBit(fFilterBitJF)) continue; // track cut selection
-    if(fRequireITSRefitJF &&  ((status&AliESDtrack::kITSrefit)==0)) continue;  
-
-    if(TMath::Abs(track->Eta()) > fTrackEtaCut) continue; 
-    if( track->Pt() < fTrackPtCutJF ) continue;
-
-    Double_t trackMom[3];
-    track->PxPyPz(trackMom);
-    TVector3 track3mom(trackMom);
-
-    Double_t dR = jet3mom.DeltaR(track3mom);
-
-    if(dR<radius){
-
-      list->Add(track);
-      
-    }
-  }
-}
-
-// ----------------------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::GetJetTracksTrackrefs(TList* list, AliAODJet* jet){
-  
-  if(!jet) return; // no jet in acc in event
-
-  // list of AOD tracks in jet cone, using trackrefs
-  
-  Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
-
-  if(!nTracks) return;
-  
-  // standard cuts + ITS refit = stdrd PWG4 cut (compose them for productions with old ESDFilter task) 
-
-  for (Int_t itrack=0; itrack<nTracks; itrack++) {
-
-    AliAODTrack* track = (AliAODTrack*) jet->GetRefTracks()->At(itrack);
-
-    UInt_t status = track->GetStatus();
-
-    if(!track->TestFilterBit(fFilterBitJF)) continue; // track cut selection
-    if(fRequireITSRefitJF &&  ((status&AliESDtrack::kITSrefit)==0)) continue;  
-
-    if(TMath::Abs(track->Eta()) > fTrackEtaCut) continue; 
-    if( track->Pt() < fTrackPtCutJF ) continue;
-        
-    list->Add(track);
-  }
-
-  //cout<<" trackrefs Size "<<nTracks<<" acc track list size "<<list->GetEntries()<<endl;
-
-}
-
-
-// ----------------------------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::FillReferenceFF(AliAODJet* jet){
-  
-
-  if(!jet) return;
-
-  TList* jetTracks = new TList(); // FIXME - avoid new/delete
-  //GetJetTracksResum(jetTracks,jet,0.7); 
-  GetJetTracksTrackrefs(jetTracks,jet); 
-
-  Double_t jetpt = jet->Pt();
-
-  TIter next(jetTracks);
-  while(AliAODTrack* track = static_cast<AliAODTrack*>(next())){
-    
-    Double_t trackpt = track->Pt();
-    Double_t z       = trackpt/jetpt;
-    Double_t xi      = TMath::Log(1/z);
-    
-    //cout<<" trackpt "<<trackpt<<" jetpt "<<jetpt<<" z "<<z<<" xi "<<xi<<endl;
-
-    if(5<jetpt && jetpt<=10){
-      fhdNdptTracksJetPt5to10->Fill(trackpt);
-      fhdNdzJets5to10->Fill(z);
-      fhdNdxiJets5to10->Fill(xi);
-    }
-    else if(10<jetpt && jetpt<=20){
-      fhdNdptTracksJetPt10to20->Fill(trackpt);
-      fhdNdzJets10to20->Fill(z);
-      fhdNdxiJets10to20->Fill(xi);
-    }
-    else if(20<jetpt && jetpt<=30){
-      fhdNdptTracksJetPt20to30->Fill(trackpt);
-      fhdNdzJets20to30->Fill(z);
-      fhdNdxiJets20to30->Fill(xi);
-    }
-    else if(30<jetpt && jetpt<=40){
-      fhdNdptTracksJetPt30to40->Fill(trackpt);
-      fhdNdzJets30to40->Fill(z);
-      fhdNdxiJets30to40->Fill(xi);
-    }
-    else if(40<jetpt && jetpt<=60){
-      fhdNdptTracksJetPt40to60->Fill(trackpt);
-      fhdNdzJets40to60->Fill(z);
-      fhdNdxiJets40to60->Fill(xi);
-    }
-  }
-
-  
-  delete jetTracks;
-  
-}
-
-// -------------------------------------------------------------
-
-void AliAnalysisTaskJetChem::FillReferencePlotsTracks(){
-
-  // eta/phi & pt tracks before/after cuts
-  // track multiplicity / evt
-
-  Int_t nTracks     = fAOD->GetNTracks();
-  Int_t countTracks = 0;
-
-  // standard cuts + ITS refit = stdrd PWG4 cut (compose them for productions with old ESDFilter task) 
-
-  for(Int_t itrack=0; itrack<nTracks; itrack++){
-    
-    AliAODTrack* track = fAOD->GetTrack(itrack);
-    
-    Double_t trackPt  = track->Pt();
-    Double_t trackPhi = track->Phi();
-    Double_t trackEta = track->Eta();
-
-    fhPhiEtaTracksNoCut->Fill(trackEta,trackPhi);
-    fhPtTracksNoCut->Fill(trackPt);
-    
-    UInt_t status = track->GetStatus();
-    
-    if(!track->TestFilterBit(fFilterBitJF)) continue; // track cut selection
-    if(fRequireITSRefitJF &&  ((status&AliESDtrack::kITSrefit)==0)) continue;  
-
-    fhPhiEtaTracks->Fill(trackEta,trackPhi);
-    fhPtTracks->Fill(trackPt);
-
-    if(TMath::Abs(track->Eta()) > fTrackEtaCut) continue; 
-    if( track->Pt() < fTrackPtCutJF ) continue;
-
-    countTracks++;
-    
-  }
-
-  fhTrackMult->Fill(countTracks);
-}