IdentifiedHighPt task for train usage
authorddobrigk <david.dobrigkeit.chinellato@cern.ch>
Sat, 4 Oct 2014 23:37:19 +0000 (20:37 -0300)
committerddobrigk <david.dobrigkeit.chinellato@cern.ch>
Sat, 4 Oct 2014 23:37:42 +0000 (20:37 -0300)
PWGLF/CMakelibPWGLFspectra.pkg
PWGLF/PWGLFspectraLinkDef.h
PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.cxx [new file with mode: 0755]
PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.h [new file with mode: 0755]
PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.cxx [new file with mode: 0755]
PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.h [new file with mode: 0755]
PWGLF/SPECTRA/IdentifiedHighPt/train/macros/AddTaskHighPtDeDx.C [new file with mode: 0755]

index aa4168a..a46b32b 100644 (file)
@@ -76,6 +76,8 @@ set ( SRCS SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAcceptanceCuts.cxx
        SPECTRA/XtAnalysis/AliJBaseTrack.cxx
        SPECTRA/XtAnalysis/AliJBaseCard.cxx
         SPECTRA/XtAnalysis/AliJCard.cxx
+        SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.cxx
+        SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.cxx
 
  )
 
@@ -85,4 +87,4 @@ set ( DHDR PWGLFspectraLinkDef.h)
 
 set ( EXPORT )
 
-set ( EINCLUDE  PWGLF/SPECTRA/Nuclei/B2 PWGLF/SPECTRA/PiKaPr/TPCTOF PWGLF/SPECTRA/PiKaPr/TPCTOFpA PWGLF/SPECTRA/PiKaPr/TOF/pp7 PWGLF/SPECTRA/ChargedHadrons/dNdPt PWGLF/SPECTRA/XtAnalysis TPC TOF PWGUD/base STEER/STEER STEER/ESD STEER/AOD STEER/CDB STEER/STEERBase ANALYSIS PWG/Tools)
+set ( EINCLUDE PWGLF/SPECTRA/IdentifiedHighPt/train PWGLF/SPECTRA/Nuclei/B2 PWGLF/SPECTRA/PiKaPr/TPCTOF PWGLF/SPECTRA/PiKaPr/TPCTOFpA PWGLF/SPECTRA/PiKaPr/TOF/pp7 PWGLF/SPECTRA/ChargedHadrons/dNdPt PWGLF/SPECTRA/XtAnalysis TPC TOF PWGUD/base STEER/STEER STEER/ESD STEER/AOD STEER/CDB STEER/STEERBase ANALYSIS PWG/Tools)
index 28bdad8..43bda68 100644 (file)
 #pragma link C++ class AliJBaseCard+;
 #pragma link C++ class AliJCard+;
 #pragma link C++ class AliJXtHistos+;
+#pragma link C++ class AliAnalysisTaskHighPtDeDx+; 
 
+#pragma link C++ class DeDxTrack+; 
+#pragma link C++ class VZEROCell+; 
+#pragma link C++ class DeDxV0+; 
+#pragma link C++ class DeDxTrackMC+; 
+#pragma link C++ class DeDxEvent+;
 
 #endif
diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.cxx b/PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.cxx
new file mode 100755 (executable)
index 0000000..622a432
--- /dev/null
@@ -0,0 +1,4510 @@
+/*
+  30/03/2012
+  Update: the bug on V0's was fixed, also an option to save VZERO cells was added
+  For the AOD analysis the filter 32 was replaced by the proper golden: 1024
+  // 1024 1<<10                        
+  trackFilter->AddCuts(esdTrackCutsH2); // add r_aa cuts
+
+  30/05/2012
+  Include neutral particles in MC truth
+  Add number of tpc shared clusters
+
+  10/01/13. 
+  Add a flag to see if the particle comes from material or WD
+  
+  20/02/13
+  Correct the bug on the multiplicity estimator. Before it was filled two times, one for the v0ana and other for the track ana. SPD vertex
+
+  10/04/13
+  Store all events with different vertex status.
+  Now the vertex is from TPC if not from SPD. Similar as in TOF
+
+  16/04/13
+  Add quality check to SPD vertex. Status -2
+
+  23/04/13
+  Store MC truth, |eta_lab|<2.4, also store y
+
+
+*/
+
+#include "AliAnalysisTaskHighPtDeDx.h"
+#include <bitset>
+// ROOT includes
+#include <TList.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TH1.h>
+#include <TParticle.h>
+#include <TFile.h>
+
+// AliRoot includes
+#include <AliAnalysisManager.h>
+#include <AliAnalysisFilter.h>
+#include <AliESDInputHandler.h>
+#include <AliESDEvent.h>
+#include <AliESDVertex.h>
+#include <AliLog.h>
+#include <AliExternalTrackParam.h>
+#include <AliESDtrackCuts.h>
+#include <AliESDVZERO.h>
+#include <AliAODVZERO.h>
+
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include <AliStack.h>
+
+#include <TTreeStream.h>
+
+#include <AliHeader.h>
+#include <AliGenPythiaEventHeader.h>
+#include <AliGenDPMjetEventHeader.h>
+
+#include "AliCentrality.h" 
+#include <AliESDv0.h>
+#include <AliKFVertex.h>
+#include <AliAODVertex.h>
+
+#include <AliAODTrack.h> 
+#include <AliAODPid.h> 
+#include <AliAODMCHeader.h> 
+#include <iostream>
+
+
+// STL includes
+#include <iostream>
+using namespace std;
+
+
+//
+// Responsible:
+// Alexandru Dobrin (Wayne State) 
+// Peter Christiansen (Lund)
+//
+
+/* 
+To do:
+
+Separate the code into two
+
+*/
+
+
+
+
+ClassImp(AliAnalysisTaskHighPtDeDx)
+
+const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
+
+//_____________________________________________________________________________
+AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
+  AliAnalysisTaskSE(),
+  fESD(0x0),
+  fAOD(0x0),
+  fMC(0x0),
+  fMCStack(0x0),
+  fMCArray(0x0),
+  fTrackFilter(0x0),
+  fTrackFilterGolden(0x0),
+  fTrackFilterTPC(0x0),
+  fV0ArrayGlobalPar(0x0), //
+  fV0ArrayTPCPar(0x0),//
+  fAnalysisType("ESD"),
+  fAnalysisMC(kFALSE),
+  fAnalysisPbPb(kFALSE),
+  fVZEROBranch(kFALSE),
+  fTPCBranch(kFALSE),
+  ftrigBit1(0x0),
+  ftrigBit2(0x0),
+  fRandom(0x0),
+  fEvent(0x0),
+  fTrackArrayGlobalPar(0x0),
+  fTrackArrayTPCPar(0x0),
+  fTrackArrayMC(0x0),
+  fVZEROArray(0x0),
+  fVtxCut(10.0),  
+  fEtaCut(0.9),  
+  fMinPt(0.1),
+  fMinPtV0(0.1),
+  fMinCent(0.0),
+  fMaxCent(100.0),
+  fMassCut(0.1),//
+  fLowPtFraction(0.01),
+  fTreeOption(0),
+  fRequireRecV0(kTRUE), //
+  fStoreMcIn(kFALSE),//
+  fMcProcessType(-999),
+  fTriggeredEventMB(-999),
+  fVtxStatus(-999),
+  fZvtx(-999),
+  fZvtxMC(-999),
+  fRun(-999),
+  fEventId(-999),
+  fListOfObjects(0), 
+  fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
+  fTree(0x0),
+  fn1(0),
+  fn2(0),
+  fcent(0)
+
+
+{
+  // Default constructor (should not be used)
+
+  //  fRandom = new TRandom(0); // 0 means random seed
+}
+
+//______________________________________________________________________________
+AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
+  AliAnalysisTaskSE(name),
+  fESD(0x0),
+  fAOD(0x0),
+  fMC(0x0),
+  fMCStack(0x0),
+  fMCArray(0x0),
+  fTrackFilter(0x0),
+  fTrackFilterGolden(0x0),
+  fTrackFilterTPC(0x0),
+  fV0ArrayGlobalPar(0x0), //
+  fV0ArrayTPCPar(0x0),//
+  fAnalysisType("ESD"),
+  fAnalysisMC(kFALSE),
+  fAnalysisPbPb(kFALSE),
+  fVZEROBranch(kFALSE),
+  fTPCBranch(kFALSE),
+  ftrigBit1(0x0),
+  ftrigBit2(0x0),
+  fRandom(0x0),
+  fEvent(0x0),
+  fTrackArrayGlobalPar(0x0),
+  fTrackArrayTPCPar(0x0),
+  fTrackArrayMC(0x0),
+  fVZEROArray(0x0),
+  fVtxCut(10.0),  
+  fEtaCut(0.9),  
+  fMinPt(0.1),
+  fMinPtV0(0.1),
+  fMinCent(0.0),
+  fMaxCent(100.0),
+  fMassCut(0.1),//
+  fLowPtFraction(0.01),
+  fTreeOption(0),
+  fRequireRecV0(kTRUE), //
+  fStoreMcIn(kFALSE),//
+  fMcProcessType(-999),
+  fTriggeredEventMB(-999),
+  fVtxStatus(-999),
+  fZvtx(-999),
+  fZvtxMC(-999),
+  fRun(-999),
+  fEventId(-999),
+  fListOfObjects(0), 
+  fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
+  fTree(0x0),
+  fn1(0),
+  fn2(0),
+  fcent(0)
+
+{
+
+  if(fTree)fTree=0;
+  // Output slot #1 writes into a TList
+  DefineOutput(1, TList::Class());
+
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
+{
+  // Destructor
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+  if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+    delete fListOfObjects;
+    fListOfObjects = 0;
+  }
+  if (fRandom) delete fRandom;
+  fRandom=0;
+  
+  // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse 
+  // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
+  
+  
+  
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
+{ 
+  // This method is called once per worker node
+  // Here we define the output: histograms and debug tree if requested 
+  // We also create the random generator here so it might get different seeds...
+  fRandom = new TRandom(0); // 0 means random seed
+
+  OpenFile(1);
+  fListOfObjects = new TList();
+  fListOfObjects->SetOwner();
+  
+  //
+  // Histograms
+  //  
+  fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
+  fListOfObjects->Add(fEvents);
+
+  fn1=new TH1F("fn1","fn1",5001,-1,5000);
+  fListOfObjects->Add(fn1);
+
+  fn2=new TH1F("fn2","fn2",5001,-1,5000);
+  fListOfObjects->Add(fn2);
+
+  fcent=new TH1F("fcent","fcent",104,-2,102);
+  fListOfObjects->Add(fcent);
+
+  fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
+  fListOfObjects->Add(fVtx);
+
+  fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
+  fListOfObjects->Add(fVtxBeforeCuts);
+  
+  fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
+  fListOfObjects->Add(fVtxAfterCuts);
+
+  if (fAnalysisMC) {    
+    fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
+    fListOfObjects->Add(fVtxMC);
+  }
+
+  if (fTreeOption) {
+
+    fTree = new TTree("tree","Event data");
+    fEvent = new DeDxEvent();
+    fTree->Branch("event", &fEvent);
+    //array with tracks
+    fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
+    fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
+    //array with v0's
+    fV0ArrayGlobalPar = new TClonesArray("DeDxV0", 1000);
+    fTree->Bronch("v0GlobalPar", "TClonesArray", &fV0ArrayGlobalPar);
+
+    if(fTPCBranch){
+      fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
+      fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
+      fV0ArrayTPCPar = new TClonesArray("DeDxV0", 1000);
+      fTree->Bronch("v0TPCPar", "TClonesArray", &fV0ArrayTPCPar);
+    }
+    if(fVZEROBranch){
+      fVZEROArray = new TClonesArray("VZEROCell", 1000);
+      fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
+    }
+
+
+    if (fAnalysisMC && fStoreMcIn) {    
+      fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
+      fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
+    }
+
+
+
+    fTree->SetDirectory(0);
+
+    fListOfObjects->Add(fTree);
+
+  }
+
+  // Post output data.
+  PostData(1, fListOfObjects);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *) 
+{
+  // Main loop
+
+  //
+  // First we make sure that we have valid input(s)!
+  //
+
+
+
+  AliVEvent *event = InputEvent();
+  if (!event) {
+     Error("UserExec", "Could not retrieve event");
+     return;
+  }
+
+
+  if (fAnalysisType == "ESD"){
+    fESD = dynamic_cast<AliESDEvent*>(event);
+    if(!fESD){
+      Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+      this->Dump();
+      return;
+    }    
+  } else {
+    fAOD = dynamic_cast<AliAODEvent*>(event);
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+      this->Dump();
+      return;
+    }    
+
+  }
+
+  
+  if (fAnalysisMC) {
+
+    if (fAnalysisType == "ESD"){
+      fMC = dynamic_cast<AliMCEvent*>(MCEvent());
+      if(!fMC){
+       Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+      }    
+
+      fMCStack = fMC->Stack();
+      
+      if(!fMCStack){
+       Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+      }    
+    } else { // AOD
+
+      fMC = dynamic_cast<AliMCEvent*>(MCEvent());
+      if(fMC)
+       fMC->Dump();
+
+      fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
+      if(!fMCArray){
+       Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
+       this->Dump();
+       return;
+      }    
+    }
+  }
+
+  
+  // Get trigger decision
+  fTriggeredEventMB = 0; //init
+
+  // cout << " trigger " << ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())) ->IsEventSelected() << "cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc"<<endl;
+    if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
+     ->IsEventSelected() & ftrigBit1 ){
+    fn1->Fill(1);
+    fTriggeredEventMB = 1;  //event triggered as minimum bias
+  }
+  if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
+     ->IsEventSelected() & ftrigBit2 ){
+    // From AliVEvent:
+    //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
+    fTriggeredEventMB += 2;  
+    fn2->Fill(1);
+  }
+  //if(fTriggeredEventMB == 0) fTriggeredEventMB = 1; //ignore trigger //hljunggr
+  //std::cout << "trigger " << fTriggeredEventMB << " ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc" << std::endl;
+  // Get process type for MC
+  fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+
+  // real data that are not triggered we skip
+  if(!fAnalysisMC && !fTriggeredEventMB)
+    return; 
+
+
+  //   std::cout << "trigger Martin: " << (bitset<20>) ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() << std::endl; print trigger
+
+
+  if (fAnalysisMC) {
+    
+
+
+    if (fAnalysisType == "ESD"){
+
+  
+
+      AliHeader* headerMC = fMC->Header();
+      if (headerMC) {
+       
+       AliGenEventHeader* genHeader = headerMC->GenEventHeader();
+       TArrayF vtxMC(3); // primary vertex  MC 
+       vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
+       if (genHeader) {
+         genHeader->PrimaryVertex(vtxMC);
+       }
+       fZvtxMC = vtxMC[2];
+       
+       // PYTHIA:
+       AliGenPythiaEventHeader* pythiaGenHeader =
+         dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
+       if (pythiaGenHeader) {  //works only for pythia
+         fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
+       }
+       // PHOJET:
+       AliGenDPMjetEventHeader* dpmJetGenHeader =
+         dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
+       if (dpmJetGenHeader) {
+         fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
+       }
+      }
+    } else { // AOD
+      
+  
+
+      AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); 
+
+
+      if(mcHeader) {
+       fZvtxMC = mcHeader->GetVtxZ();
+       
+
+
+       if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
+         fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
+       } else {
+         fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
+       }
+      }
+    }
+  }
+  
+  // There are 3 cases
+  // Vertex: NO  - status -1
+  // SDP Vertex wit bad quality: status -2
+  // Vertex: YES : outside cut - status 0
+  //             : inside cut  - status 1
+  // We have to be careful how we normalize because we probably want to
+  // normalize to:
+  // Nevents=(No vertex + outside + inside)/(out + in)*in
+  
+
+  if (fAnalysisType == "ESD"){
+    
+    const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
+    if(vtxESD->GetNContributors()<1) {
+      // SPD vertex
+      vtxESD = fESD->GetPrimaryVertexSPD();
+      /* quality checks on SPD-vertex */
+      TString vertexType = vtxESD->GetTitle();
+      if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))  
+       fZvtx  = -1599; //vertex = 0x0; //
+      else if (vtxESD->GetNContributors()<1) 
+       fZvtx  = -999; //vertex = 0x0; //
+      else
+       fZvtx = vtxESD->GetZ();
+    }  
+    else
+      fZvtx = vtxESD->GetZ();
+    
+  }
+  else // AOD
+    fZvtx = GetVertex(fAOD);
+    
+  fVtxStatus = -999;
+  
+
+
+  if(fZvtx<-1500) {
+    
+    fVtxStatus = -2;
+  } 
+  else if(fZvtx<-990&&fZvtx>-1500) {
+    
+    fVtxStatus = -1;
+    if(fTriggeredEventMB)
+      fVtx->Fill(0);
+    if(fAnalysisMC)
+      fVtxMC->Fill(0);
+  } else {
+    
+    if(fTriggeredEventMB)
+      fVtx->Fill(1);
+    if(fAnalysisMC)
+      fVtxMC->Fill(1);
+    fVtxBeforeCuts->Fill(fZvtx);
+    fVtxStatus = 0;
+    if (TMath::Abs(fZvtx) < fVtxCut) { 
+      fVtxAfterCuts->Fill(fZvtx);
+      fVtxStatus = 1;
+    }
+  }
+  
+  // store MC event data no matter what
+  //if(fAnalysisMC && fStoreMcIn) {
+  if(fAnalysisMC) {
+
+    if (fAnalysisType == "ESD") {
+      ProcessMCTruthESD();
+    } else { // AOD
+
+      ProcessMCTruthAOD();
+  
+    }
+  }      
+  
+
+
+  Float_t centralityV0M = -10;
+  Float_t centralityV0A = -10;
+  Float_t centralityZNA = -10;
+  Float_t centralityCL1 = -10;
+
+  // only analyze triggered events
+  //   cout << " fTriggeredEventMB " << fTriggeredEventMB << "fAnalysisType " << fAnalysisType << " fAnalysisPbPb " <<fAnalysisPbPb << endl;
+  if(fTriggeredEventMB) {
+    //cout << "inside triggered !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
+    if (fAnalysisType == "ESD"){
+      if(fAnalysisPbPb){
+       AliCentrality *centObject = fESD->GetCentrality();
+       centralityV0M = centObject->GetCentralityPercentile("V0M");
+       centralityV0A = centObject->GetCentralityPercentile("V0A");
+       centralityZNA = centObject->GetCentralityPercentile("ZNA");
+       centralityCL1 = centObject->GetCentralityPercentile("CL1");
+     
+       //if(  !(fMinCent>=20 && fMaxCent <=40) ) return; //take out //hljunggr
+       if((centralityV0M>fMaxCent)||(centralityV0M<fMinCent))return; //hljunggr comment out, put back in!!
+      }
+      fcent->Fill(centralityV0A);
+
+      AnalyzeESD(fESD);
+    } else { // AOD
+      if(fAnalysisPbPb){
+       AliCentrality *centObject = fAOD->GetCentrality();
+       if(centObject){
+         centralityV0M = centObject->GetCentralityPercentile("V0M"); 
+         centralityV0A = centObject->GetCentralityPercentile("V0A");
+         centralityZNA = centObject->GetCentralityPercentile("ZNA");
+         centralityCL1 = centObject->GetCentralityPercentile("CL1");
+       }
+
+       //if(  !(fMinCent>=20 && fMaxCent <=40) ) return; //take out //hljunggr
+       if((centralityV0M>fMaxCent)||(centralityV0M<fMinCent))return; //hljunggr comment out, put back in!!
+      }
+      fcent->Fill(centralityV0M);
+      AnalyzeAOD(fAOD);
+    }
+  }
+
+  if( fTreeOption) {
+    
+    fEvent->process = fMcProcessType;
+    fEvent->trig    = fTriggeredEventMB;
+    fEvent->zvtxMC  = fZvtxMC;
+    fEvent->cent      = centralityV0M;
+    fEvent->centV0A      = centralityV0A;
+    fEvent->centZNA      = centralityZNA;
+    fEvent->centCL1      = centralityCL1;
+
+
+    fTree->Fill();
+    if(fTrackArrayGlobalPar)fTrackArrayGlobalPar->Clear();
+    if(fV0ArrayGlobalPar)fV0ArrayGlobalPar->Clear();
+
+    if(fTPCBranch){
+      if(fTrackArrayTPCPar)fTrackArrayTPCPar->Clear();
+      if(fV0ArrayTPCPar)fV0ArrayTPCPar->Clear();
+    }
+    if(fVZEROArray)
+      fVZEROArray->Clear();
+
+    if (fAnalysisMC)    
+      fTrackArrayMC->Clear();
+      
+
+  }
+  
+  // Post output data.
+  PostData(1, fListOfObjects);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
+{
+  fRun  = esdEvent->GetRunNumber();
+  fEventId = 0;
+  if(esdEvent->GetHeader())
+    fEventId = GetEventIdAsLong(esdEvent->GetHeader());
+  
+  Short_t isPileup = esdEvent->IsPileupFromSPD();
+  /*
+  cout"  nTPCtracks="<<esdEvent->nTPCtracks<<"  nITStracks="<<esdEvent->nITStracks<<endl;
+  if(esdEvent->nTPCtracks==0 && esdEvent->nITStracks==0){
+    cout<<"Evento defectuoso!!"<<endl;
+    }*/
+
+
+  //  Int_t     event     = esdEvent->GetEventNumberInFile();
+  UInt_t    time      = esdEvent->GetTimeStamp();
+  //  ULong64_t trigger   = esdEvent->GetTriggerMask();
+  Float_t   magf      = esdEvent->GetMagneticField();
+
+
+
+
+
+  if(fTriggeredEventMB) {// Only MC case can we have not triggered events
+    
+    // accepted event
+    fEvents->Fill(0);
+    
+    
+    //Change, 10/04/13. Now accept all events to do a correct normalization
+    //if(fVtxStatus!=1) return; // accepted vertex
+    Int_t nESDTracks = esdEvent->GetNumberOfTracks();
+    
+    ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
+    ProduceArrayV0ESD( esdEvent, kGlobalTrk );//v0's
+
+    if(fTPCBranch){
+      ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
+      ProduceArrayV0ESD( esdEvent, kTPCTrk );
+    }
+    
+    fEvents->Fill(1);
+    
+    if(fVZEROBranch){
+      AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
+      for (Int_t iCh=0; iCh<64; ++iCh) { 
+       Float_t multv=esdV0->GetMultiplicity(iCh); 
+       Int_t intexv=iCh;
+       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
+       cellv0->cellindex=intexv;
+       cellv0->cellmult= multv;
+      }   
+    }
+
+
+
+  } // end if triggered
+  
+  if(fTreeOption) {
+
+    fEvent->run       = fRun;
+    fEvent->eventid   = fEventId;
+    fEvent->time      = time;
+    //fEvent->cent      = centrality;
+    fEvent->mag       = magf;
+    fEvent->zvtx      = fZvtx;
+    fEvent->vtxstatus = fVtxStatus;
+    fEvent->pileup    = isPileup;
+
+  }
+
+
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
+{
+  fRun  = aodEvent->GetRunNumber();
+  fEventId = 0;
+  if(aodEvent->GetHeader())
+    fEventId = GetEventIdAsLong(aodEvent->GetHeader());
+   
+  UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
+  Float_t   magf      = aodEvent->GetMagneticField();
+
+  //Int_t     trackmult = 0; // no pt cuts
+  //Int_t     nadded    = 0;
+  Short_t isPileup = aodEvent->IsPileupFromSPD();
+
+  
+
+  if(fTriggeredEventMB) {// Only MC case can we have not triggered events
+    
+    // accepted event
+    fEvents->Fill(0);
+    
+    //if(fVtxStatus!=1) return; // accepted vertex
+    //Int_t nAODTracks = aodEvent->GetNumberOfTracks();  
+
+    ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
+
+    ProduceArrayV0AOD( aodEvent, kGlobalTrk );
+    if(fTPCBranch){
+      ProduceArrayTrksAOD( aodEvent, kTPCTrk );
+      ProduceArrayV0AOD( aodEvent, kTPCTrk );
+    }
+    fEvents->Fill(1);
+
+    if(fVZEROBranch){
+      AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
+      for (Int_t iCh=0; iCh<64; ++iCh) { 
+       Float_t multv=esdV0->GetMultiplicity(iCh); 
+       Int_t intexv=iCh;
+       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
+       cellv0->cellindex=intexv;
+       cellv0->cellmult= multv;
+      }   
+    }
+
+
+
+  } // end if triggered
+  
+  if(fTreeOption) {
+
+    //Sort(fTrackArrayGlobalPar, kFALSE);
+
+    fEvent->run       = fRun;
+    fEvent->eventid   = fEventId;
+    fEvent->time      = time;
+    //fEvent->cent      = centrality;
+    fEvent->mag       = magf;
+    fEvent->zvtx      = fZvtx;
+    fEvent->vtxstatus = fVtxStatus;
+    //fEvent->trackmult = trackmult;
+    //fEvent->n         = nadded;
+    fEvent->pileup    = isPileup;
+  }
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
+{
+  Float_t zvtx = -999;
+  
+  const AliVVertex* primaryVertex = event->GetPrimaryVertex(); 
+  
+  if(primaryVertex->GetNContributors()>0)
+    zvtx = primaryVertex->GetZ();
+
+  return zvtx;
+}
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const 
+{
+  // return our internal code for pions, kaons, and protons
+  
+  Short_t pidCode = 6;
+  
+  switch (TMath::Abs(pdgCode)) {
+  case 211:
+    pidCode = 1; // pion
+    break;
+  case 321:
+    pidCode = 2; // kaon
+    break;
+  case 2212:
+    pidCode = 3; // proton
+    break;
+  case 11:
+    pidCode = 4; // electron
+    break;
+  case 13:
+    pidCode = 5; // muon
+    break;
+  default:
+    pidCode = 6;  // something else?
+  };
+  
+  return pidCode;
+}
+
+
+//_____________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD() 
+{
+  // Fill the special MC histogram with the MC truth info
+  
+  //cout << "mc truth " << " fAnalysisMC " << fAnalysisMC << " fstoremcin " << fStoreMcIn << endl;
+  Short_t trackmult = 0;
+  Short_t nadded    = 0;
+  const Int_t nTracksMC = fMCStack->GetNtrack();
+
+  Float_t  sphericityMC=GetSphericityTrue(fMCStack, 0.8, 0.5);
+  Float_t  spherocityMC=GetSpherocityTrue(fMCStack, 0.8, 0.5);
+
+
+  for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
+    
+   
+
+    //Cuts
+    if(!(fMCStack->IsPhysicalPrimary(iTracks)))
+      continue;
+    
+    TParticle* trackMC = fMCStack->Particle(iTracks);
+    
+    TParticlePDG* pdgPart = trackMC ->GetPDG();
+    Double_t chargeMC = pdgPart->Charge();
+
+    
+    if (TMath::Abs(trackMC->Eta()) > 2.4 )
+      continue;
+    
+    trackmult++;
+    
+    //   "charge:pt:p:eta:phi:pidCode"
+    Float_t ptMC      = trackMC->Pt();
+    Float_t pMC       = trackMC->P();
+    Float_t etaMC     = trackMC->Eta();
+    Float_t phiMC     = trackMC->Phi();
+    Float_t yMC       = trackMC->Y();
+
+    Int_t pdgCode = trackMC->GetPdgCode();
+    Short_t pidCodeMC = 0;
+    pidCodeMC = GetPidCode(pdgCode);
+    
+    // Here we want to add some of the MC histograms!
+    
+    // And therefore we first cut here!
+    if (trackMC->Pt() < fMinPt) {
+      
+      // Keep small fraction of low pT tracks
+      if(fRandom->Rndm() > fLowPtFraction)
+       continue;
+    } // else {
+    // Here we want to add the high pt part of the MC histograms!
+    //    }
+    
+    if(fTreeOption) {
+      
+      DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
+      nadded++;
+      
+      track->pMC   = pMC;
+      track->ptMC  = ptMC;
+      track->etaMC = etaMC;
+      track->phiMC = phiMC;
+      track->yMC   = yMC;
+      track->qMC   = Short_t(chargeMC);
+      track->pidMC = pidCodeMC;
+      track->pdgMC = pdgCode;
+    }
+    
+  }//MC track loop
+  
+  if(fTreeOption) {
+    
+    Sort(fTrackArrayMC, kTRUE);
+
+    fEvent->trackmultMC = trackmult;
+    fEvent->nMC         = nadded;
+    fEvent->sphericityMC         = sphericityMC;
+    fEvent->spherocityMC         = spherocityMC;
+
+  }
+  
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD() 
+{
+  // Fill the special MC histogram with the MC truth info
+  //cout << " debug 1 " << endl;
+  Short_t trackmult = 0;
+  Short_t nadded    = 0;
+  const Int_t nTracksMC = fMCArray->GetEntriesFast();
+  //cout << " debug 2 " << endl;
+
+  //cout << " debug 3 " << endl;
+  for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
+    //cout << " debug 31 " << endl;
+    AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
+    //cout << " debug 32 " << endl;
+    //Cuts
+    if(!(trackMC->IsPhysicalPrimary()))
+      continue;
+    //cout << " debug 33 " << endl;
+    Double_t chargeMC = trackMC->Charge();
+
+    
+    if (TMath::Abs(trackMC->Eta()) > 2.4 )
+      continue;
+    //cout << " debug 34 " << endl;
+    trackmult++;
+    
+    //   "charge:pt:p:eta:phi:pidCode"
+    Float_t ptMC      = trackMC->Pt();
+    Float_t pMC       = trackMC->P();
+    Float_t etaMC     = trackMC->Eta();
+    Float_t phiMC     = trackMC->Phi();
+    Float_t yMC       = trackMC->Y();
+    //cout << " debug 35 " << endl;
+    Int_t pdgCode = trackMC->PdgCode();
+    Short_t pidCodeMC = 0;
+    pidCodeMC = GetPidCode(pdgCode);
+   
+
+    // Here we want to add some of the MC histograms!
+    
+    // And therefore we first cut here!
+    if (trackMC->Pt() < fMinPt) {
+      
+      // Keep small fraction of low pT tracks
+      if(fRandom->Rndm() > fLowPtFraction)
+       continue;
+    } // else {
+    // Here we want to add the high pt part of the MC histograms!
+    //    }
+    //cout << " debug 36 " << endl;
+    //cout << "fTreeOption " << fTreeOption << endl;
+    if(fTreeOption) {
+      
+
+      //cout << " debug 361, nadded =  " << nadded << " fTrackArrayMC" << fTrackArrayMC <<  endl;
+      DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
+      //cout << " debug 3611, nadded =  " << nadded <<  endl;
+      nadded++;
+      //cout << " debug 3612, nadded =  " << nadded <<  endl;
+      //cout << " debug 362 " << endl;
+      track->pMC   = pMC;
+      //cout << " debug 363 " << endl;
+      track->ptMC  = ptMC;
+      track->etaMC = etaMC;
+      track->phiMC = phiMC;
+      track->yMC   = yMC;
+      track->qMC   = Short_t(chargeMC);
+      track->pidMC = pidCodeMC;
+      track->pdgMC = pdgCode;
+      //cout << " debug 364 " << endl;
+
+    }
+    //cout << " debug 37 " << endl;
+  }//MC track loop
+  //cout << " debug 4 " << endl;
+  if(fTreeOption) {
+    
+    Sort(fTrackArrayMC, kTRUE);
+
+    fEvent->trackmultMC = trackmult;
+    fEvent->nMC         = nadded;
+    fEvent->sphericityMC         = 0;
+    fEvent->spherocityMC         = 0;
+
+  }
+  //cout << " debug 5 " << endl;
+}
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
+  //
+  // Get the process type of the event.  PYTHIA
+  //
+  // source PWG0   dNdpt 
+
+  Short_t globalType = -1; //init
+      
+  if(pythiaType==92||pythiaType==93){
+    globalType = 2; //single diffractive
+  }
+  else if(pythiaType==94){
+    globalType = 3; //double diffractive
+  }
+  //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
+  else {
+    globalType = 1;  //non diffractive
+  }
+  return globalType;
+}
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
+  //
+  // get the process type of the event.  PHOJET
+  //
+  //source PWG0   dNdpt 
+  // can only read pythia headers, either directly or from cocktalil header
+  Short_t globalType = -1;
+  
+  if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
+    globalType = 1;
+  }
+  else if (dpmJetType==5 || dpmJetType==6) {
+    globalType = 2;
+  }
+  else if (dpmJetType==7) {
+    globalType = 3;
+  }
+  return globalType;
+}
+
+//_____________________________________________________________________________
+ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
+{
+  // To have a unique id for each event in a run!
+  // Modified from AliRawReader.h
+  return ((ULong64_t)header->GetBunchCrossNumber()+
+         (ULong64_t)header->GetOrbitNumber()*3564+
+         (ULong64_t)header->GetPeriodNumber()*16777215*3564);
+}
+
+
+//____________________________________________________________________
+TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // Taken from AliPWG0Helper class
+  //
+
+  Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
+  if (motherLabel < 0)
+    return 0;
+
+  return stack->Particle(motherLabel);
+}
+
+//____________________________________________________________________
+Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // returns its label
+  //
+  // Taken from AliPWG0Helper class
+  //
+  const Int_t nPrim  = stack->GetNprimary();
+  
+  while (label >= nPrim) {
+
+    //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle) {
+      
+      AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
+      return -1;
+    }
+    
+    // find mother
+    if (particle->GetMother(0) < 0) {
+
+      AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+      return -1;
+    }
+    
+    label = particle->GetMother(0);
+  }
+  
+  return label;
+}
+
+//____________________________________________________________________
+AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // Taken from AliPWG0Helper class
+  //
+
+  AliAODMCParticle* mcPart = startParticle;
+
+  while (mcPart)
+    {
+      
+      if(mcPart->IsPrimary())
+       return mcPart;
+      
+      Int_t mother = mcPart->GetMother();
+
+      mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
+    }
+
+  return 0;
+}
+
+
+//V0______________________________________
+//____________________________________________________________________
+TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // Taken from AliPWG0Helper class
+  //
+
+  Int_t nSteps = 0;
+
+  Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
+  if (motherLabel < 0)
+    return 0;
+
+  return stack->Particle(motherLabel);
+}
+
+//____________________________________________________________________
+Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // returns its label
+  //
+  // Taken from AliPWG0Helper class
+  //
+  nSteps = 0;
+  const Int_t nPrim  = stack->GetNprimary();
+  
+  while (label >= nPrim) {
+
+    //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+    nSteps++; // 1 level down
+    
+    TParticle* particle = stack->Particle(label);
+    if (!particle) {
+      
+      AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
+      return -1;
+    }
+    
+    // find mother
+    if (particle->GetMother(0) < 0) {
+
+      AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+      return -1;
+    }
+    
+    label = particle->GetMother(0);
+  }
+  
+  return label;
+}
+
+//____________________________________________________________________
+AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+  // Taken from AliPWG0Helper class
+  //
+
+  nSteps = 0;
+
+  AliAODMCParticle* mcPart = startParticle;
+
+  while (mcPart)
+    {
+      
+      if(mcPart->IsPrimary())
+       return mcPart;
+      
+      Int_t mother = mcPart->GetMother();
+
+      mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
+      nSteps++; // 1 level down
+    }
+
+  return 0;
+}
+
+
+
+
+//_____________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC) 
+{
+  const Int_t n = array->GetEntriesFast(); 
+  if(n==0) {
+
+    if(isMC) 
+      fEvent->ptmaxMC = 0;
+    else
+      fEvent->ptmax   = 0;
+      
+    return;
+  }
+
+  Float_t ptArray[n];
+  Int_t   indexArray[n];
+  
+  for(Int_t i = 0; i < n; i++) {
+
+    Float_t pt = 0;
+    if(isMC) {
+      DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
+      pt = track->ptMC;
+    } else {
+      DeDxTrack* track = (DeDxTrack*)array->At(i);
+      pt = track->pt;
+    }
+    ptArray[i]    = pt;
+    indexArray[i] = i;
+  }
+
+  TMath::Sort(n, ptArray, indexArray);
+  
+  // set max pt
+  if(isMC) {
+    DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
+    fEvent->ptmaxMC = track->ptMC;
+  } else {
+    DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
+    fEvent->ptmax   = track->pt;
+  }
+  
+  // set order of each track
+  for(Int_t i = 0; i < n; i++) {
+    
+    if(isMC) {
+      DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
+      track->orderMC = i;
+    } else {
+      DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
+      track->order = i;
+    }
+  }
+}
+//__________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
+  
+  const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
+  Int_t trackmult=0;
+
+
+  Int_t nadded=0;
+  if( analysisMode == kGlobalTrk ){
+    if(fTrackArrayGlobalPar)
+      fTrackArrayGlobalPar->Clear();
+  } else if( analysisMode == kTPCTrk ){
+    if(fTrackArrayTPCPar)
+      fTrackArrayTPCPar->Clear();
+  }
+
+  //Fix, for pPb LHC13b pass2 it does not work
+  
+  AliESDtrackCuts* esdTrackCutsGolden = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
+  esdTrackCutsGolden->SetMinNCrossedRowsTPC(120);
+  esdTrackCutsGolden->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+  esdTrackCutsGolden->SetMaxChi2PerClusterITS(36);
+  esdTrackCutsGolden->SetMaxFractionSharedTPCClusters(0.4);
+  esdTrackCutsGolden->SetMaxChi2TPCConstrainedGlobal(36);
+  
+  trackmult=esdTrackCutsGolden->GetReferenceMultiplicity(ESDevent, AliESDtrackCuts::MultEstTrackType(0), 0.5);
+  
+  //trackmult=AliESDtrackCuts::GetReferenceMultiplicity(ESDevent, AliESDtrackCuts::kTrackletsITSTPC, 0.5);
+
+
+  Float_t spherocity=-1;
+  Float_t sphericity=-1;
+  Float_t spherocityTPC=-1;
+  Float_t sphericityTPC=-1;
+
+
+  spherocity=GetSpherocity(ESDevent,fTrackFilterGolden,0.8, 0.5, kFALSE);
+  sphericity=GetSphericity(ESDevent,fTrackFilterGolden,0.8, 0.5, kFALSE);
+  //Lot of memory consumption solved by reducing the number of files per job
+  //spherocityTPC=GetSpherocity(ESDevent,fTrackFilterGolden,0.8, 0.5, kTRUE);
+  //sphericityTPC=GetSphericity(ESDevent,fTrackFilterGolden,0.8, 0.5, kTRUE);
+  
+  spherocityTPC=spherocity;
+  sphericityTPC=sphericity;
+
+
+  if( analysisMode==kGlobalTrk ){  
+    
+    for(Int_t iT = 0; iT < nESDTracks; iT++) {
+      
+      AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+      
+      
+      if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
+       continue;
+      
+      UShort_t filterFlag = 0;
+      
+      UInt_t selectDebug = 0;
+      if (fTrackFilterGolden) {
+       selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
+       if (selectDebug) {
+         filterFlag +=1;
+       }
+      }
+      
+      if (fTrackFilterTPC) {
+       
+       selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
+       if (selectDebug){//only tracks which pass the TPC-only track cuts
+         filterFlag +=2;
+         
+       }
+       
+      }
+      
+      if (fTrackFilter) {
+       selectDebug = fTrackFilter->IsSelected(esdTrack);
+       if (selectDebug) {
+         filterFlag +=4;
+       }
+      }
+      
+     
+      if(filterFlag==0)
+       continue;
+      
+      Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
+   
+
+      //
+      // Track was accepted
+      //      
+      
+      // Here we want to add histograms!
+      
+      if (esdTrack->Pt() < fMinPt) {
+       
+       // Keep small fraction of low pT tracks
+       if(fRandom->Rndm() > fLowPtFraction)
+         continue;
+      } // else {
+      // Here we want to add the high pt part of the histograms!
+      //    }
+    
+      Short_t charge  = esdTrack->Charge();
+      Float_t pt      = esdTrack->Pt();
+      Float_t p       = esdTrack->P(); 
+      Float_t eta     = esdTrack->Eta();
+      Float_t phi     = esdTrack->Phi();
+      Short_t ncl     = esdTrack->GetTPCsignalN();
+      Short_t neff    = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      //         Short_t nclf    = esdTrack->GetTPCNclsF();
+      Float_t dedx    = esdTrack->GetTPCsignal();
+      Float_t tpcchi  = 0;
+      if(esdTrack->GetTPCNcls() > 0)
+       tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
+      Float_t b[2];
+      Float_t bCov[3];
+      esdTrack->GetImpactParameters(b,bCov);
+      Float_t dcaxy   = b[0];
+      Float_t dcaz    = b[1];
+      Double_t p_con[3] = {0, 0, 0};
+      esdTrack->GetConstrainedPxPyPz(p_con);
+      
+      
+      //       Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
+      // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
+      // Float_t pttpc   = tpcParam->Pt();
+      // Float_t ptpc    = tpcParam->P();
+      
+      Float_t ptMC        = 0;
+      Short_t pidCode     = 0; // 0 = real data / no mc track!
+      Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   pdgMother   = 0;
+      
+      
+      //with Globaltrack parameters????
+      if(fAnalysisMC) {
+       
+       const Int_t label = TMath::Abs(esdTrack->GetLabel());
+       TParticle* mcTrack = fMCStack->Particle(label);     
+       if (mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(label))
+           primaryFlag = 1;
+         
+         //10/01/13. Add a flag to see if it is from material or WD
+
+         if(fMCStack->IsSecondaryFromWeakDecay(label))
+           primaryFlag = 2;
+
+         if(fMCStack->IsSecondaryFromMaterial(label))
+           primaryFlag = 3;
+
+
+
+
+         Int_t pdgCode = mcTrack->GetPdgCode();
+         pidCode = GetPidCode(pdgCode);
+         
+         ptMC      = mcTrack->Pt();
+         
+         TParticle* mother = FindPrimaryMother(fMCStack, label);
+         pdgMother = mother->GetPdgCode();
+       }
+      }
+      
+    
+      //TOF
+      /*
+       Float_t beta = -99;
+       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
+       if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
+         beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
+         }
+      */
+      Bool_t IsTOFout=kFALSE;
+      if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+       IsTOFout=kTRUE;
+      
+      Bool_t HasTOFTime=kTRUE;
+      if ((esdTrack->GetStatus()&AliESDtrack::kTIME)==0)
+       HasTOFTime=kFALSE;
+      
+      Bool_t IsMatched=kFALSE;
+      if (!(esdTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+       IsMatched=kFALSE;
+      else
+       IsMatched=kTRUE;
+
+      Float_t lengthtrack=esdTrack->GetIntegratedLength();
+
+      Float_t timeTOF=esdTrack->GetTOFsignal();
+
+      Double_t inttime[5]={0,0,0,0,0}; 
+      esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
+      Float_t fexptimeel=inttime[0];
+      Float_t fexptimemu=inttime[1];
+      Float_t fexptimepi=inttime[2];
+      Float_t fexptimeka=inttime[3];
+      Float_t fexptimepr=inttime[4];  
+
+      
+      if(fTreeOption) {
+       
+       DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
+       nadded++;
+       
+       track->p          = p;
+       track->pt         = pt;
+       //        track->ptcon   = pt_con;
+       //        track->tpcchi  = tpcchi;
+       track->eta        = eta;
+       track->phi        = phi;
+       track->q          = charge;
+       track->filter     = filterFlag;
+       track->ncl        = ncl;
+       track->neff       = neff;
+       track->dedx       = dedx;
+
+       track->isTOFout   = IsTOFout;
+       track->hasTOFtime = HasTOFTime;
+       track->isTOFmatched = IsMatched;
+       track->flength    = lengthtrack;
+       track->ftimetof   = timeTOF;
+       track->exptoftimeel = fexptimeel;
+       track->exptoftimemu = fexptimemu;
+       track->exptoftimepi = fexptimepi;
+       track->exptoftimeka = fexptimeka;
+       track->exptoftimepr = fexptimepr;
+
+       track->dcaxy      = dcaxy;
+       track->dcaz       = dcaz;
+       track->pid        = pidCode;
+       track->primary    = primaryFlag;
+       track->pttrue     = ptMC;
+       track->mother     = pdgMother;
+       track->tpcnclS    = sharedtpcclusters;
+       
+      }
+    }//end of track loop
+
+  }//end global
+  
+  else if( analysisMode==kTPCTrk ){  
+    const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
+    if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
+
+    for(Int_t iT = 0; iT < nESDTracks; iT++) {
+      
+      AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+
+      AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
+      
+      if(!trackTPC) continue;
+    
+      UInt_t selectDebug = 0;
+      selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
+      if(selectDebug==0) continue;
+    
+      if(trackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParam;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate = false;
+       relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
+                                          kVeryBig,&exParam);
+       if(!relate){
+         delete trackTPC;
+         continue;
+       }
+       trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
+                     exParam.GetCovariance());
+      }
+      else continue;     
+      
+      Int_t sharedtpcclusters = trackTPC->GetTPCnclsS();
+
+      if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
+       continue;
+      
+      //
+      // Track was accepted
+      //      
+      
+      // Here we want to add histograms!
+      
+      if (trackTPC->Pt() < fMinPt) {
+       
+       // Keep small fraction of low pT tracks
+       if(fRandom->Rndm() > fLowPtFraction)
+         continue;
+      } // else {
+      // Here we want to add the high pt part of the histograms!
+      //    }
+      
+      Short_t charge  = trackTPC->Charge();
+      Float_t pt      = trackTPC->Pt();
+      Float_t p       = trackTPC->P(); 
+      Float_t eta     = trackTPC->Eta();
+      Float_t phi     = trackTPC->Phi();
+      Short_t ncl     = trackTPC->GetTPCsignalN();
+      Short_t neff    = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      //         Short_t nclf    = esdTrack->GetTPCNclsF();
+      Float_t dedx    = trackTPC->GetTPCsignal();
+      Float_t tpcchi  = 0;
+      if(trackTPC->GetTPCNcls() > 0)
+       tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
+      //Float_t b[2];
+      //Float_t bCov[3];
+      //trackTPC->GetImpactParameters(b,bCov);
+      //Float_t dcaxy   = b[0];
+      //Float_t dcaz    = b[1];
+
+      Float_t dcaxy = 0.;
+      Float_t dcaz = 0.;
+      trackTPC->GetImpactParameters(dcaxy,dcaz);
+
+
+      Double_t p_con[3] = {0, 0, 0};
+      trackTPC->GetConstrainedPxPyPz(p_con);
+      
+    
+      // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
+      // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
+      // Float_t pttpc   = tpcParam->Pt();
+      // Float_t ptpc    = tpcParam->P();
+      
+      Float_t ptMC        = 0;
+      Short_t pidCode     = 0; // 0 = real data / no mc track!
+      Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   pdgMother   = 0;
+      
+      
+      //with Globaltrack parameters????
+      if(fAnalysisMC) {
+       
+       const Int_t label = TMath::Abs(esdTrack->GetLabel());
+       TParticle* mcTrack = fMCStack->Particle(label);     
+       if (mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(label))
+           primaryFlag = 1;
+         
+         
+         //10/01/13. Add a flag to see if it is from material or WD
+         
+         if(fMCStack->IsSecondaryFromWeakDecay(label))
+           primaryFlag = 2;
+         
+         if(fMCStack->IsSecondaryFromMaterial(label))
+           primaryFlag = 3;
+
+
+         Int_t pdgCode = mcTrack->GetPdgCode();
+         pidCode = GetPidCode(pdgCode);
+         
+         ptMC      = mcTrack->Pt();
+         
+         TParticle* mother = FindPrimaryMother(fMCStack, label);
+         pdgMother = mother->GetPdgCode();
+       }
+      }
+    
+    
+      //TOF
+      /*
+      Float_t beta = -99;
+      if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
+       if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
+         beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
+      }
+      */
+      if(fTreeOption) {
+       
+       DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
+       nadded++;
+       
+       tracktpc->p          = p;
+       tracktpc->pt         = pt;
+       //        track->ptcon   = pt_con;
+       //        track->tpcchi  = tpcchi;
+       tracktpc->eta        = eta;
+       tracktpc->phi        = phi;
+       tracktpc->q          = charge;
+       tracktpc->filter     = 1;
+       tracktpc->ncl        = ncl;
+       tracktpc->neff       = neff;
+       tracktpc->dedx       = dedx;
+
+       tracktpc->isTOFout   = 0;
+       tracktpc->hasTOFtime = 0;
+       tracktpc->isTOFmatched = 0;
+       tracktpc->flength    = 0;
+       tracktpc->ftimetof   = 0;
+       tracktpc->exptoftimeel = 0;
+       tracktpc->exptoftimemu = 0;
+       tracktpc->exptoftimepi = 0;
+       tracktpc->exptoftimeka = 0;
+       tracktpc->exptoftimepr = 0;
+
+
+       tracktpc->dcaxy      = dcaxy;
+       tracktpc->dcaz       = dcaz;
+       tracktpc->pid        = pidCode;
+       tracktpc->primary    = primaryFlag;
+       tracktpc->pttrue     = ptMC;
+       tracktpc->mother     = pdgMother;
+       tracktpc->tpcnclS    = sharedtpcclusters;
+      }
+
+
+      if(trackTPC)
+       delete trackTPC;
+
+
+    }//end of track loop
+  }//end of: if isglobal==kFALSE
+
+
+  if(fTreeOption) {
+
+    if( analysisMode==kGlobalTrk ){
+      Sort(fTrackArrayGlobalPar, kFALSE);
+      
+      fEvent->trackmult = trackmult;
+      fEvent->n         = nadded;
+
+      fEvent->sphericity         = sphericity;
+      fEvent->spherocity         = spherocity;
+      fEvent->sphericityTPC      = sphericityTPC;
+      fEvent->spherocityTPC      = spherocityTPC;
+
+    }
+    else if( analysisMode==kTPCTrk ){
+      Sort(fTrackArrayTPCPar, kFALSE);
+    }
+
+  }
+
+
+}
+//__________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
+
+
+  Int_t     trackmult = 0; // no pt cuts
+  Int_t     nadded    = 0;
+  Int_t nAODTracks = AODevent->GetNumberOfTracks();
+  if( analysisMode == kGlobalTrk ){
+    
+    if(fTrackArrayGlobalPar)
+      fTrackArrayGlobalPar->Clear();
+    
+  } 
+  if( analysisMode == kTPCTrk ){
+    if(fTrackArrayTPCPar)
+      fTrackArrayTPCPar->Clear();
+  }
+
+
+
+  //const AliAODVertex*        vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
+  //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
+
+
+  
+  if( analysisMode==kGlobalTrk ){  
+    
+     const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
+    for(Int_t iT = 0; iT < nAODTracks; iT++) {
+      
+      AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+      
+      //FOR AOD068, filterCut_Set2=KTRUE WHEN THE TRACK SATISFIES THE CUTS FROM JET ANALYSIS
+      
+      UShort_t filterFlag = 0;
+
+      
+      if (fTrackFilterGolden) {
+       
+       // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
+       if(aodTrack->TestFilterBit(1024)) {
+         filterFlag +=1;
+       }
+      }
+      
+      
+      if (fTrackFilterTPC) {
+       // TPC only cuts is bit 1 according to above macro
+       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+       if(aodTrack->TestFilterBit(1)){
+         filterFlag +=2;
+       }
+      }
+      
+    
+      
+      if (fTrackFilter) {
+       /*
+       // tighter cuts on primary particles for high pT tracks
+       // take the standard cuts, which include already 
+       // ITSrefit and use only primaries...
+       bit 32
+       */
+       if(aodTrack->TestFilterBit(32)) {
+         filterFlag +=4;
+       }
+
+      }
+      if(aodTrack->TestFilterBit(768)) {
+       filterFlag +=8; // 272 for 2012 hljunggr, included to get hybrid tracks
+       if(!aodTrack->IsHybridGlobalConstrainedGlobal()) filterFlag+=128; //just to test the selection, remove this line (/Martin)
+      }
+      
+      
+      // if(aodTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
+      //       // Minimum number of clusters
+      //       Float_t nCrossedRowsTPC = aodTrack->GetTPCClusterInfo(2,1);
+      //   if (nCrossedRowsTPC >= 70) {
+      //         Int_t findable=aodTrack->GetTPCNclsF();
+      //         if (findable > 0){ 
+      //           if (nCrossedRowsTPC/findable >= 0.8) filterFlag += 16;
+      //         }
+      //       }
+       
+
+
+      // }
+      
+
+      
+      if(filterFlag==0)
+       continue;
+      
+      
+      Double_t b[2], cov[3];
+      if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
+       filterFlag = 32; // propagation failed!!!!!;
+      Float_t dcaxy   = b[0];
+      Float_t dcaz    = b[1];
+      
+      TBits sharedTPC=aodTrack->GetTPCSharedMap();
+      Int_t sharedtpcclusters=sharedTPC.CountBits(0)-sharedTPC.CountBits(159);
+
+
+      // As I understand this routine recalculates the momentum so it should be called first!
+      //Double_t b[2], cov[3];
+  
+      
+      //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
+      //       filterFlag = 32; // propagation failed!!!!!
+      
+      if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+       continue;
+      if (aodTrack->Pt() < fMinPt) {
+       
+       // Keep small fraction of low pT tracks
+       if(fRandom->Rndm() > fLowPtFraction)
+         continue;
+      } // else {
+      // Here we want to add the high pt part of the histograms!
+      //    }
+     
+     
+
+      
+      Short_t charge  = aodTrack->Charge();
+      Float_t pt      = aodTrack->Pt();
+      Float_t p       = aodTrack->P(); 
+      Float_t eta     = aodTrack->Eta();
+      Float_t phi     = aodTrack->Phi();
+      AliAODPid* aodPid = aodTrack->GetDetPid();
+      Short_t ncl     = -10;
+      Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      //         Short_t nclf    = aodTrack->GetTPCNclsF();
+      Float_t dedx    = -10;
+      Bool_t IsTOFout=kFALSE;
+      Bool_t HasTOFTime=kTRUE;
+      Bool_t IsMatched=kFALSE;
+      Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
+       
+      Float_t timeTOF=-999;
+      Float_t fexptimeel=-999;
+      Float_t fexptimemu=-999;
+      Float_t fexptimepi=-999;
+      Float_t fexptimeka=-999;
+      Float_t fexptimepr=-999;  
+
+
+      //Float_t beta = -99;
+      if(aodPid) {
+       ncl     = aodPid->GetTPCsignalN();
+       dedx    = aodPid->GetTPCsignal();
+       //TOF
+       /*
+       if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
+         Double_t tof[5];
+         aodPid->GetIntegratedTimes(tof);
+         beta = tof[0]/aodPid->GetTOFsignal();
+         }*/
+
+
+
+       if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+         IsTOFout=kTRUE;
+       
+
+       if ((aodTrack->GetStatus()&AliESDtrack::kTIME)==0)
+         HasTOFTime=kFALSE;
+       
+
+       if (!(aodTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+         IsMatched=kFALSE;
+       else
+         IsMatched=kTRUE;
+       
+       lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
+       
+       timeTOF=aodPid->GetTOFsignal();
+       
+       Double_t inttime[5]={0,0,0,0,0}; 
+       aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
+       fexptimeel=inttime[0];
+       fexptimemu=inttime[1];
+       fexptimepi=inttime[2];
+       fexptimeka=inttime[3];
+       fexptimepr=inttime[4];  
+       
+
+
+      }
+      //       Float_t tpcchi = aodTrack->Chi2perNDF();
+      
+      // Double_t p_con[3] = {0, 0, 0};
+      // aodTrack->GetConstrainedPxPyPz(p_con);
+      //       Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
+      // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
+      // Float_t pttpc   = tpcParam->Pt();
+      // Float_t ptpc    = tpcParam->P();
+      
+      Float_t ptMC        = 0;
+      Short_t pidCode     = 0; // 0 = real data / no mc track!
+      Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   pdgMother   = 0;
+      
+      if(fAnalysisMC) {
+       
+       const Int_t label = TMath::Abs(aodTrack->GetLabel());
+       
+       AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
+       if (mcTrack){
+         
+         if(mcTrack->IsPhysicalPrimary())
+           primaryFlag = 1;
+         
+
+         Int_t pdgCode = mcTrack->GetPdgCode();
+         pidCode = GetPidCode(pdgCode);
+         
+         ptMC      = mcTrack->Pt();
+         
+         AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
+         pdgMother = mother->GetPdgCode();         
+       }
+      }
+    
+      
+      if(fTreeOption) {
+       
+       DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
+       nadded++;
+       
+       track->p          = p;
+       track->pt         = pt;
+       //        track->ptcon   = pt_con;
+       //        track->tpcchi  = tpcchi;
+       track->eta        = eta;
+       track->phi        = phi;
+       track->q          = charge;
+       track->filter     = filterFlag;
+       track->ncl        = ncl;
+       track->neff       = neff;
+       track->dedx       = dedx;
+       //track->beta       = beta;
+       track->isTOFout   = IsTOFout;
+       track->hasTOFtime = HasTOFTime;
+       track->isTOFmatched = IsMatched;
+       track->flength    = lengthtrack;
+       track->ftimetof   = timeTOF;
+       track->exptoftimeel = fexptimeel;
+       track->exptoftimemu = fexptimemu;
+       track->exptoftimepi = fexptimepi;
+       track->exptoftimeka = fexptimeka;
+       track->exptoftimepr = fexptimepr;
+
+       track->dcaxy      = dcaxy;
+       track->dcaz       = dcaz;
+       track->pid        = pidCode;
+       track->primary    = primaryFlag;
+       track->pttrue     = ptMC;
+       track->mother     = pdgMother;
+       track->tpcnclS    = sharedtpcclusters;
+      }
+    }//end of track loop
+
+
+  }//end of global track analysis
+  
+
+
+  else if( analysisMode==kTPCTrk ){  
+    const AliAODVertex*        vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
+    if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
+
+
+    //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
+    for(Int_t iT = 0; iT < nAODTracks; iT++) {
+      
+      AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+      
+      
+      UShort_t filterFlag = 0;
+  
+          
+      // TPC only cuts is bit 1 according to above macro
+      // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+      if(aodTrack->TestFilterBit(128)){
+       filterFlag +=2;
+       trackmult++;
+      }
+      
+      
+      if(filterFlag==0)
+       continue;
+      
+      TBits sharedTPC=aodTrack->GetTPCSharedMap();
+      Int_t sharedtpcclusters=sharedTPC.CountBits(0)-sharedTPC.CountBits(159);
+
+      
+      Double_t b[2], cov[3];
+      //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
+      //Double_t b[2] = {-99., -99.};
+      //Double_t bCov[3] = {-99., -99., -99.};
+      //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
+      if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
+       filterFlag = 32; // propagation failed!!!!!;
+      Float_t dcaxy   = b[0];
+      Float_t dcaz    = b[1];
+
+
+
+
+      // As I understand this routine recalculates the momentum so it should be called first!
+      //Double_t b[2], cov[3];
+  
+      
+      //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
+      //       filterFlag = 32; // propagation failed!!!!!
+      
+      if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+       continue;
+      if (aodTrack->Pt() < fMinPt) {
+       
+       // Keep small fraction of low pT tracks
+       if(fRandom->Rndm() > fLowPtFraction)
+         continue;
+      } // else {
+      // Here we want to add the high pt part of the histograms!
+      //    }
+     
+     
+
+      
+      Short_t charge  = aodTrack->Charge();
+      Float_t pt      = aodTrack->Pt();
+      Float_t p       = aodTrack->P(); 
+      Float_t eta     = aodTrack->Eta();
+      Float_t phi     = aodTrack->Phi();
+      AliAODPid* aodPid = aodTrack->GetDetPid();
+      Short_t ncl     = -10;
+      Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      //         Short_t nclf    = aodTrack->GetTPCNclsF();
+      Float_t dedx    = -10;
+      //Float_t beta = -99;
+      if(aodPid) {
+       ncl     = aodPid->GetTPCsignalN();
+       dedx    = aodPid->GetTPCsignal();
+       //TOF
+       /*
+       if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
+         Double_t tof[5];
+         aodPid->GetIntegratedTimes(tof);
+         beta = tof[0]/aodPid->GetTOFsignal();
+         }*/
+      }
+      //       Float_t tpcchi = aodTrack->Chi2perNDF();
+      
+      // Double_t p_con[3] = {0, 0, 0};
+      // aodTrack->GetConstrainedPxPyPz(p_con);
+      //       Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
+      // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
+      // Float_t pttpc   = tpcParam->Pt();
+      // Float_t ptpc    = tpcParam->P();
+      
+      Float_t ptMC        = 0;
+      Short_t pidCode     = 0; // 0 = real data / no mc track!
+      Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   pdgMother   = 0;
+      
+      if(fAnalysisMC) {
+       
+       const Int_t label = TMath::Abs(aodTrack->GetLabel());
+       
+       AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
+       if (mcTrack){
+         
+         if(mcTrack->IsPhysicalPrimary())
+           primaryFlag = 1;
+
+         
+         Int_t pdgCode = mcTrack->GetPdgCode();
+         pidCode = GetPidCode(pdgCode);
+         
+         ptMC      = mcTrack->Pt();
+         
+         AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
+         pdgMother = mother->GetPdgCode();         
+       }
+      }
+    
+    
+      if(fTreeOption) {
+       
+       DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
+       nadded++;
+       
+       tracktpc->p          = p;
+       tracktpc->pt         = pt;
+       //        track->ptcon   = pt_con;
+       //        track->tpcchi  = tpcchi;
+       tracktpc->eta        = eta;
+       tracktpc->phi        = phi;
+       tracktpc->q          = charge;
+       tracktpc->filter     = filterFlag;
+       tracktpc->ncl        = ncl;
+       tracktpc->neff       = neff;
+       tracktpc->dedx       = dedx;
+       //tracktpc->beta       = beta;
+
+       tracktpc->isTOFout   = 0;
+       tracktpc->hasTOFtime = 0;
+       tracktpc->isTOFmatched = 0;
+       tracktpc->flength    = 0;
+       tracktpc->ftimetof   = 0;
+       tracktpc->exptoftimeel = 0;
+       tracktpc->exptoftimemu = 0;
+       tracktpc->exptoftimepi = 0;
+       tracktpc->exptoftimeka = 0;
+       tracktpc->exptoftimepr = 0;
+
+
+       tracktpc->dcaxy      = dcaxy;
+       tracktpc->dcaz       = dcaz;
+       tracktpc->pid        = pidCode;
+       tracktpc->primary    = primaryFlag;
+       tracktpc->pttrue     = ptMC;
+       tracktpc->mother     = pdgMother;
+       tracktpc->tpcnclS       = sharedtpcclusters;
+      }
+    }//end of track loop
+
+
+  }//end of global track analysis
+
+
+  if(fTreeOption) {
+    
+    if( analysisMode==kGlobalTrk ){
+      Sort(fTrackArrayGlobalPar, kFALSE);
+      
+      fEvent->trackmult = trackmult;
+      fEvent->n         = nadded;
+      fEvent->sphericity         = 0;
+      fEvent->spherocity         = 0;
+
+    }
+    if( analysisMode==kTPCTrk ){
+      Sort(fTrackArrayTPCPar, kFALSE);
+    }
+    
+  }
+  
+}
+//__________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
+  Int_t nv0s = ESDevent->GetNumberOfV0s();
+  if(nv0s<1)return;
+  //Int_t     trackmult = 0; // no pt cuts
+  Int_t     nadded    = 0;
+  const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
+  if (!myBestPrimaryVertex) return;
+  if (!(myBestPrimaryVertex->GetStatus())) return;
+  Double_t  lPrimaryVtxPosition[3];
+  myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
+  
+  Double_t  lPrimaryVtxCov[6];
+  myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
+  Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
+  
+  AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
+
+
+  if( analysisMode == kGlobalTrk ){
+    if(fV0ArrayGlobalPar)
+      fV0ArrayGlobalPar->Clear();
+  } else if( analysisMode == kTPCTrk ){
+    if(fV0ArrayTPCPar)
+      fV0ArrayTPCPar->Clear();
+  }
+  
+  if( analysisMode==kGlobalTrk ){  
+
+
+
+    
+    // ################################
+    // #### BEGINNING OF V0 CODE ######
+    // ################################
+
+    
+    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+      
+      // This is the begining of the V0 loop  
+      AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
+      if (!esdV0) continue;
+      
+      // AliESDTrack (V0 Daughters)
+      UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
+      UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
+      
+      AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
+      AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
+      if (!pTrack || !nTrack) {
+       Printf("ERROR: Could not retreive one of the daughter track");
+       continue;
+      }
+      
+      // Remove like-sign
+      if (pTrack->GetSign() == nTrack->GetSign()) {
+       //cout<< "like sign, continue"<< endl;
+       continue;
+      } 
+      
+      // Eta cut on decay products
+      if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
+       continue;
+      
+      // Pt cut on decay products
+      if (esdV0->Pt() < fMinPtV0)
+       //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
+       continue;
+
+      //filter for positive track
+      UShort_t filterFlag_p = 0;
+      
+      UInt_t selectDebug_p = 0;
+      if (fTrackFilterGolden) {
+       selectDebug_p = fTrackFilterGolden->IsSelected(pTrack);
+       if (selectDebug_p) {
+         filterFlag_p +=1;
+       }
+      }
+      
+      if (fTrackFilterTPC) {
+       
+       selectDebug_p = fTrackFilterTPC->IsSelected(pTrack);
+       if (selectDebug_p){//only tracks which pass the TPC-only track cuts
+         filterFlag_p +=2;
+       }
+       
+      }
+      
+      if (fTrackFilter) {
+       selectDebug_p = fTrackFilter->IsSelected(pTrack);
+       if (selectDebug_p) {
+         filterFlag_p +=4;
+       }
+      }
+      
+
+    
+      
+
+      //filter for negative track
+      UShort_t filterFlag_n = 0;
+      
+      UInt_t selectDebug_n = 0;
+      if (fTrackFilterGolden) {
+       selectDebug_n = fTrackFilterGolden->IsSelected(nTrack);
+       if (selectDebug_n) {
+         filterFlag_n +=1;
+       }
+      }
+      
+      if (fTrackFilterTPC) {
+       
+       selectDebug_n = fTrackFilterTPC->IsSelected(nTrack);
+       if (selectDebug_n){//only tracks which pass the TPC-only track cuts
+         filterFlag_n +=2;
+       }
+       
+      }
+      
+      if (fTrackFilter) {
+       selectDebug_n = fTrackFilter->IsSelected(nTrack);
+       if (selectDebug_n) {
+         filterFlag_n +=4;
+       }
+      }
+      
+
+
+    
+
+
+      // Check if switch does anything!
+      Bool_t isSwitched = kFALSE;
+      if (pTrack->GetSign() < 0) { // switch
+       
+       isSwitched = kTRUE;
+       AliESDtrack* helpTrack = nTrack;
+       nTrack = pTrack;
+       pTrack = helpTrack;
+      }        
+      
+      Float_t alpha = esdV0->AlphaV0();
+      Float_t ptarm = esdV0->PtArmV0();
+      // Double_t pVtxPos= v0->PrimaryVtxPosition();      
+      
+      Double_t  lV0Position[3];
+      esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
+      
+      Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
+      Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
+                                           TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
+                                           TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
+      AliKFVertex primaryVtxKF( *myPrimaryVertex );
+      AliKFParticle::SetField(ESDevent->GetMagneticField());
+      
+      // Also implement switch here!!!!!!
+      AliKFParticle* negEKF  = 0; // e-
+      AliKFParticle* posEKF  = 0; // e+
+      AliKFParticle* negPiKF = 0; // pi -
+      AliKFParticle* posPiKF = 0; // pi +
+      AliKFParticle* posPKF  = 0; // p
+      AliKFParticle* negAPKF = 0; // p-bar
+      
+      if(!isSwitched) {
+       negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
+       posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
+       negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
+       posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
+       posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
+       negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
+      } else { // switch + and - 
+       negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
+       posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
+       negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
+       posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
+       posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
+       negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
+      }
+      
+      AliKFParticle v0GKF;  // Gamma e.g. from pi0
+      v0GKF+=(*negEKF);
+      v0GKF+=(*posEKF);
+      v0GKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0K0sKF; // K0 short
+      v0K0sKF+=(*negPiKF);
+      v0K0sKF+=(*posPiKF);
+      v0K0sKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0LambdaKF; // Lambda
+      v0LambdaKF+=(*negPiKF);
+      v0LambdaKF+=(*posPKF);   
+      v0LambdaKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0AntiLambdaKF; // Lambda-bar
+      v0AntiLambdaKF+=(*posPiKF);
+      v0AntiLambdaKF+=(*negAPKF);
+      v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
+      
+      Double_t deltaInvMassG     = v0GKF.GetMass();
+      Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
+      Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
+      Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
+
+
+      if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
+        TMath::Abs(deltaInvMassL) > fMassCut &&
+        TMath::Abs(deltaInvMassAntiL) > fMassCut &&
+       TMath::Abs(deltaInvMassG) > fMassCut )
+       continue;
+
+      
+      // Extract track information
+      
+      Short_t pcharge  = pTrack->Charge();
+      Float_t ppt      = pTrack->Pt();
+      Float_t pp       = pTrack->P(); 
+      Float_t peta     = pTrack->Eta();
+      Float_t pphi     = pTrack->Phi();
+      Short_t pncl     = pTrack->GetTPCsignalN();
+      Short_t pneff    = Short_t(pTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t pdedx    = pTrack->GetTPCsignal();
+      Int_t   psharedtpcclusters = pTrack->GetTPCnclsS();
+      
+      Float_t ptpcchi  = 0;
+      if(pTrack->GetTPCNcls() > 0)
+       ptpcchi = pTrack->GetTPCchi2()/Float_t(pTrack->GetTPCNcls());
+      
+
+      Bool_t pIsTOFout=kFALSE;
+      if ((pTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+       pIsTOFout=kTRUE;
+      
+      Bool_t pHasTOFTime=kTRUE;
+      if ((pTrack->GetStatus()&AliESDtrack::kTIME)==0)
+       pHasTOFTime=kFALSE;
+      
+      Bool_t pIsMatched=kFALSE;
+      if (!(pTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+       pIsMatched=kFALSE;
+      else
+       pIsMatched=kTRUE;
+
+      Float_t plengthtrack=pTrack->GetIntegratedLength();
+
+      Float_t ptimeTOF=pTrack->GetTOFsignal();
+
+      Double_t pinttime[5]={0,0,0,0,0}; 
+      pTrack->GetIntegratedTimes(pinttime);// Returns the array with integrated times for each particle hypothesis
+      Float_t pfexptimeel=pinttime[0];
+      Float_t pfexptimemu=pinttime[1];
+      Float_t pfexptimepi=pinttime[2];
+      Float_t pfexptimeka=pinttime[3];
+      Float_t pfexptimepr=pinttime[4];  
+
+
+
+
+      Short_t ncharge  = nTrack->Charge();
+      Float_t npt      = nTrack->Pt();
+      Float_t np       = nTrack->P(); 
+      Float_t neta     = nTrack->Eta();
+      Float_t nphi     = nTrack->Phi();
+      Short_t nncl     = nTrack->GetTPCsignalN();
+      Short_t nneff    = Short_t(nTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t ndedx    = nTrack->GetTPCsignal();
+      Int_t   nsharedtpcclusters = nTrack->GetTPCnclsS();
+
+      Float_t ntpcchi  = 0;
+      if(nTrack->GetTPCNcls() > 0)
+       ntpcchi = nTrack->GetTPCchi2()/Float_t(nTrack->GetTPCNcls());
+      
+
+
+      Bool_t nIsTOFout=kFALSE;
+      if ((nTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+       nIsTOFout=kTRUE;
+      
+      Bool_t nHasTOFTime=kTRUE;
+      if ((nTrack->GetStatus()&AliESDtrack::kTIME)==0)
+       nHasTOFTime=kFALSE;
+      
+      Bool_t nIsMatched=kFALSE;
+      if (!(nTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+       nIsMatched=kFALSE;
+      else
+       nIsMatched=kTRUE;
+
+      Float_t nlengthtrack=nTrack->GetIntegratedLength();
+
+      Float_t ntimeTOF=nTrack->GetTOFsignal();
+
+      Double_t ninttime[5]={0,0,0,0,0}; 
+      nTrack->GetIntegratedTimes(ninttime);// Returns the array with integrated times for each particle hypothesis
+      Float_t nfexptimeel=ninttime[0];
+      Float_t nfexptimemu=ninttime[1];
+      Float_t nfexptimepi=ninttime[2];
+      Float_t nfexptimeka=ninttime[3];
+      Float_t nfexptimepr=ninttime[4];  
+
+
+
+      Float_t b[2];
+      Float_t bCov[3];
+      pTrack->GetImpactParameters(b,bCov);
+      Float_t pdcaxy   = b[0];
+      Float_t pdcaz    = b[1];
+      nTrack->GetImpactParameters(b,bCov);
+      Float_t ndcaxy   = b[0];
+      Float_t ndcaz    = b[1];
+      
+      Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
+      Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
+      Float_t p_ptMC        = 0;
+      Short_t p_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   p_pdgMother   = 0;
+      Float_t n_ptMC        = 0;
+      Short_t n_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   n_pdgMother   = 0;
+      if(fAnalysisMC) {
+       
+       Int_t p_mother_label = 0;
+       Int_t p_mother_steps = 0;
+       Int_t n_mother_label = 0;
+       Int_t n_mother_steps = 0;
+       
+       // positive track
+       const Int_t p_label = TMath::Abs(pTrack->GetLabel());
+       TParticle* p_mcTrack = fMCStack->Particle(p_label);         
+       if (p_mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(p_label))
+           p_primaryFlag = 1;
+         
+         
+         //10/01/13. Add a flag to see if it is from material or WD
+
+         if(fMCStack->IsSecondaryFromWeakDecay(p_label))
+           p_primaryFlag = 2;
+
+         if(fMCStack->IsSecondaryFromMaterial(p_label))
+           p_primaryFlag = 3;
+
+         Int_t p_pdgCode = p_mcTrack->GetPdgCode();
+         p_pidCode = GetPidCode(p_pdgCode);
+         
+         p_ptMC      = p_mcTrack->Pt();
+         
+         p_mother_label = FindPrimaryMotherLabelV0(fMCStack, p_label, 
+                                                 p_mother_steps);
+         if(p_mother_label>0) {
+           TParticle* p_mother = fMCStack->Particle(p_mother_label);
+           p_pdgMother = p_mother->GetPdgCode();
+         }
+       }
+       
+       // negative track
+       const Int_t n_label = TMath::Abs(nTrack->GetLabel());
+       TParticle* n_mcTrack = fMCStack->Particle(n_label);         
+       if (n_mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(n_label))
+           n_primaryFlag = 1;
+         
+         //10/01/13. Add a flag to see if it is from material or WD
+
+         if(fMCStack->IsSecondaryFromWeakDecay(n_label))
+           n_primaryFlag = 2;
+
+         if(fMCStack->IsSecondaryFromMaterial(n_label))
+           n_primaryFlag = 3;
+
+         Int_t n_pdgCode = n_mcTrack->GetPdgCode();
+         n_pidCode = GetPidCode(n_pdgCode);
+         
+         n_ptMC      = n_mcTrack->Pt();
+         
+         n_mother_label = FindPrimaryMotherLabelV0(fMCStack, n_label, 
+                                                 n_mother_steps);
+         if(n_mother_label>0) {
+           TParticle* n_mother = fMCStack->Particle(n_mother_label);
+           n_pdgMother = n_mother->GetPdgCode();
+         }
+       }
+       
+       // Check if V0 is primary = first and the same mother of both partciles
+       if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
+         pdgV0 = p_pdgMother;
+         if(p_mother_steps == 1 && n_mother_steps == 1) {
+           primaryV0 = 1;
+         }
+       }
+      }
+
+    
+      if(fTreeOption) {
+       
+
+
+       DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
+       nadded++;
+       
+       // v0 data
+       v0data->p       = esdV0->P();
+       v0data->pt      = esdV0->Pt();
+       v0data->eta     = esdV0->Eta();
+       v0data->phi     = esdV0->Phi();
+       v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
+       v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
+       v0data->dmassG  = deltaInvMassG;
+       v0data->dmassK0 = deltaInvMassK0s;
+       v0data->dmassL  = deltaInvMassL;
+       v0data->dmassAL = deltaInvMassAntiL;
+       v0data->alpha   = alpha;
+       v0data->ptarm   = ptarm;
+       v0data->decayr  = lV0Radius;
+       v0data->decayl  = lV0DecayLength;
+       
+       // New parameters
+       v0data->status  = esdV0->GetOnFlyStatus();
+       v0data->chi2    = esdV0->GetChi2V0();
+       v0data->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
+       // cospt: as I understand this means that the pointing to the vertex
+       // is fine so I remove the dcaxy and dcaz for the V= class
+       v0data->dcadaughters = esdV0->GetDcaV0Daughters();
+       v0data->primary = primaryV0;
+       v0data->pdg     = pdgV0;
+       
+       // positive track
+       v0data->ptrack.p       = pp;
+       v0data->ptrack.pt      = ppt;
+       //        v0data->ptrack.ptcon   = ppt_con;
+       //        v0data->ptrack.tpcchi  = ptpcchi;
+       v0data->ptrack.eta     = peta;
+       v0data->ptrack.phi     = pphi;
+       v0data->ptrack.q       = pcharge;
+       v0data->ptrack.ncl     = pncl;
+       v0data->ptrack.neff    = pneff;
+       v0data->ptrack.dedx    = pdedx;
+
+
+       v0data->ptrack.isTOFout   = pIsTOFout;
+       v0data->ptrack.hasTOFtime = pHasTOFTime;
+       v0data->ptrack.isTOFmatched = pIsMatched;
+       v0data->ptrack.flength    =   plengthtrack;
+       v0data->ptrack.ftimetof   = ptimeTOF;
+       v0data->ptrack.exptoftimeel = pfexptimeel;
+       v0data->ptrack.exptoftimemu = pfexptimemu;
+       v0data->ptrack.exptoftimepi = pfexptimepi;
+       v0data->ptrack.exptoftimeka = pfexptimeka;
+       v0data->ptrack.exptoftimepr = pfexptimepr;
+
+       v0data->ptrack.dcaxy   = pdcaxy;
+       v0data->ptrack.dcaz    = pdcaz;
+       v0data->ptrack.pid     = p_pidCode;
+       v0data->ptrack.primary = p_primaryFlag;
+       v0data->ptrack.pttrue  = p_ptMC;
+       v0data->ptrack.mother  = p_pdgMother;
+       v0data->ptrack.filter  = filterFlag_p;
+       v0data->ptrack.tpcnclS    = psharedtpcclusters;
+
+       // negative track
+       v0data->ntrack.p       = np;
+       v0data->ntrack.pt      = npt;
+       //        v0data->ntrack.ptcon   = npt_con;
+       //        v0data->ntrack.tpcchi  = ntpcchi;
+       v0data->ntrack.eta     = neta;
+       v0data->ntrack.phi     = nphi;
+       v0data->ntrack.q       = ncharge;
+       v0data->ntrack.ncl     = nncl;
+       v0data->ntrack.neff    = nneff;
+       v0data->ntrack.dedx    = ndedx;
+
+       v0data->ntrack.isTOFout   = nIsTOFout;
+       v0data->ntrack.hasTOFtime = nHasTOFTime;
+       v0data->ntrack.isTOFmatched = nIsMatched;
+       v0data->ntrack.flength    =   nlengthtrack;
+       v0data->ntrack.ftimetof   = ntimeTOF;
+       v0data->ntrack.exptoftimeel = nfexptimeel;
+       v0data->ntrack.exptoftimemu = nfexptimemu;
+       v0data->ntrack.exptoftimepi = nfexptimepi;
+       v0data->ntrack.exptoftimeka = nfexptimeka;
+       v0data->ntrack.exptoftimepr = nfexptimepr;
+
+
+       v0data->ntrack.dcaxy   = ndcaxy;
+       v0data->ntrack.dcaz    = ndcaz;
+       v0data->ntrack.pid     = n_pidCode;
+       v0data->ntrack.primary = n_primaryFlag;
+       v0data->ntrack.pttrue  = n_ptMC;
+       v0data->ntrack.mother  = n_pdgMother;
+       v0data->ntrack.filter  = filterFlag_n;
+       v0data->ntrack.tpcnclS    = nsharedtpcclusters;
+
+      }
+      
+      // clean up loop over v0
+
+      delete negPiKF;
+      delete posPiKF;
+      delete posPKF;
+      delete negAPKF;
+  
+
+
+    }
+  
+    // clean up event
+    //delete myPrimaryVertex;
+  }
+  else if( analysisMode==kTPCTrk ){  
+    const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
+    if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
+    
+    
+    // ################################
+    // #### BEGINNING OF V0 CODE ######
+    // ################################
+
+    
+    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+      
+      // This is the begining of the V0 loop  
+      AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
+      if (!esdV0) continue;
+      
+      // AliESDTrack (V0 Daughters)
+      UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
+      UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
+      
+      AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
+      AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
+      if (!pTrack || !nTrack) {
+       Printf("ERROR: Could not retreive one of the daughter track");
+       continue;
+      }
+
+      // Remove like-sign
+      if (pTrack->GetSign() == nTrack->GetSign()) {
+       //cout<< "like sign, continue"<< endl;
+       continue;
+      } 
+
+      AliESDtrack *pTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),pTrack->GetID());
+      AliESDtrack *nTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),nTrack->GetID());
+
+      if (!pTrackTPC || !nTrackTPC) {
+       Printf("ERROR: Could not retreive one of the daughter TPC track");
+       continue;
+      }
+
+      //filter for positive track
+      UInt_t selectDebug_p = 0;
+      UShort_t filterFlag_p = 0;
+      selectDebug_p = fTrackFilterTPC->IsSelected(pTrackTPC);
+      //if(selectDebug_p==0) continue;
+
+      //filter for negative track
+      UInt_t selectDebug_n = 0;
+      UShort_t filterFlag_n = 0;
+      selectDebug_n = fTrackFilterTPC->IsSelected(nTrackTPC);
+      //if(selectDebug_n==0 || selectDebug_p==0) continue;
+
+      if(selectDebug_p)
+       filterFlag_p += 1;
+
+      if(pTrackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParamp;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate_p = false;
+       relate_p = pTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
+                                               kVeryBig,&exParamp);
+       if(!relate_p){
+         delete pTrackTPC;
+         continue;
+       }
+       pTrackTPC->Set(exParamp.GetX(),exParamp.GetAlpha(),exParamp.GetParameter(),
+                      exParamp.GetCovariance());
+      }
+      else continue;     
+      
+
+      //filter for negative track
+      if(selectDebug_n)
+       filterFlag_n += 1;
+
+  
+      if(nTrackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParamn;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate_n = false;
+       relate_n = nTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
+                                               kVeryBig,&exParamn);
+       if(!relate_n){
+         delete nTrackTPC;
+         continue;
+       }
+       nTrackTPC->Set(exParamn.GetX(),exParamn.GetAlpha(),exParamn.GetParameter(),
+                      exParamn.GetCovariance());
+      }
+      else continue;  
+      
+      
+      // Eta cut on decay products
+      if(TMath::Abs(pTrackTPC->Eta()) > fEtaCut || TMath::Abs(nTrackTPC->Eta()) > fEtaCut)
+       continue;
+      
+      // Pt cut on decay products
+      if (esdV0->Pt() < fMinPtV0)
+       //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
+       continue;
+      
+      // Check if switch does anything!
+      Bool_t isSwitched = kFALSE;
+      if (pTrackTPC->GetSign() < 0) { // switch
+       
+       isSwitched = kTRUE;
+       AliESDtrack* helpTrack = nTrack;
+       nTrackTPC = pTrackTPC;
+       pTrackTPC = helpTrack;
+      }        
+      
+
+      // Extract track information
+      
+      Short_t pcharge  = pTrackTPC->Charge();
+      Float_t ppt      = pTrackTPC->Pt();
+      Float_t pp       = pTrackTPC->P(); 
+      Float_t peta     = pTrackTPC->Eta();
+      Float_t pphi     = pTrackTPC->Phi();
+      Short_t pncl     = pTrackTPC->GetTPCsignalN();
+      Short_t pneff    = Short_t(pTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t pdedx    = pTrackTPC->GetTPCsignal();
+      Int_t   psharedtpcclusters = pTrackTPC->GetTPCnclsS();
+
+      Float_t ptpcchi  = 0;
+      if(pTrackTPC->GetTPCNcls() > 0)
+       ptpcchi = pTrackTPC->GetTPCchi2()/Float_t(pTrackTPC->GetTPCNcls());
+      
+      Short_t ncharge  = nTrackTPC->Charge();
+      Float_t npt      = nTrackTPC->Pt();
+      Float_t np       = nTrackTPC->P(); 
+      Float_t neta     = nTrackTPC->Eta();
+      Float_t nphi     = nTrackTPC->Phi();
+      Short_t nncl     = nTrackTPC->GetTPCsignalN();
+      Short_t nneff    = Short_t(nTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t ndedx    = nTrackTPC->GetTPCsignal();
+      Int_t   nsharedtpcclusters = nTrackTPC->GetTPCnclsS();
+
+      Float_t ntpcchi  = 0;
+      if(nTrackTPC->GetTPCNcls() > 0)
+       ntpcchi = nTrackTPC->GetTPCchi2()/Float_t(nTrackTPC->GetTPCNcls());
+      
+      Float_t bp[2]={0,0};
+      Float_t bCovp[3]={0,0,0};
+      pTrackTPC->GetImpactParameters(bp,bCovp);
+      Float_t pdcaxy   = bp[0];
+      Float_t pdcaz    = bp[1];
+
+      Float_t bn[2]={0,0};
+      Float_t bCovn[3]={0,0,0};
+      nTrackTPC->GetImpactParameters(bn,bCovn);
+      Float_t ndcaxy   = bn[0];
+      Float_t ndcaz    = bn[1];
+
+
+      Float_t alpha = esdV0->AlphaV0();
+      Float_t ptarm = esdV0->PtArmV0();
+      // Double_t pVtxPos= v0->PrimaryVtxPosition();      
+      
+      Double_t  lV0Position[3];
+      esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
+      
+      Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
+      Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
+                                           TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
+                                           TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
+      AliKFVertex primaryVtxKF( *myPrimaryVertex );
+      AliKFParticle::SetField(ESDevent->GetMagneticField());
+      
+      // Also implement switch here!!!!!!
+      AliKFParticle* negEKF  = 0; // e-
+      AliKFParticle* posEKF  = 0; // e+
+      AliKFParticle* negPiKF = 0; // pi -
+      AliKFParticle* posPiKF = 0; // pi +
+      AliKFParticle* posPKF  = 0; // p
+      AliKFParticle* negAPKF = 0; // p-bar
+      
+      
+      if(!isSwitched) {
+       /*      
+               negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
+               posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
+               negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
+               posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
+               posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
+               negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
+       */
+
+       negEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) , 11);
+       posEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) ,-11);
+       negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-211);
+       posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 211);
+       posPKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 2212);
+       negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-2212);
+
+
+
+      } else { // switch + and - 
+       /*
+       negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
+       posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
+       negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
+       posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
+       posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
+       negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
+       */
+
+       negEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)), 11);
+       posEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)),-11);
+       negPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-211);
+       posPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 211);
+       posPKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 2212);
+       negAPKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-2212);
+
+
+      }
+   
+
+
+
+      
+      AliKFParticle v0GKF;  // Gamma e.g. from pi0
+      v0GKF+=(*negEKF);
+      v0GKF+=(*posEKF);
+      v0GKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0K0sKF; // K0 short
+      v0K0sKF+=(*negPiKF);
+      v0K0sKF+=(*posPiKF);
+      v0K0sKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0LambdaKF; // Lambda
+      v0LambdaKF+=(*negPiKF);
+      v0LambdaKF+=(*posPKF);   
+      v0LambdaKF.SetProductionVertex(primaryVtxKF);
+      
+      AliKFParticle v0AntiLambdaKF; // Lambda-bar
+      v0AntiLambdaKF+=(*posPiKF);
+      v0AntiLambdaKF+=(*negAPKF);
+      v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
+      
+      Double_t deltaInvMassG     = v0GKF.GetMass();
+      Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
+      Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
+      Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
+      
+      
+      if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
+        TMath::Abs(deltaInvMassL) > fMassCut &&
+        TMath::Abs(deltaInvMassAntiL) > fMassCut &&
+        TMath::Abs(deltaInvMassG) > fMassCut )
+       continue;
+      
+
+
+      // TODO: Whe should these be different? Different mass hypothesis = energy loss
+      // This is not important for us as we focus on the decay products!
+      // Double_t ptK0s        = v0K0sKF.GetPt(); 
+      // Double_t ptL          = v0LambdaKF.GetPt();
+      // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
+      
+      Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
+      Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
+      Float_t p_ptMC        = 0;
+      Short_t p_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   p_pdgMother   = 0;
+      Float_t n_ptMC        = 0;
+      Short_t n_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   n_pdgMother   = 0;
+      if(fAnalysisMC) {
+       
+       Int_t p_mother_label = 0;
+       Int_t p_mother_steps = 0;
+       Int_t n_mother_label = 0;
+       Int_t n_mother_steps = 0;
+       
+       // positive track
+       const Int_t p_label = TMath::Abs(pTrackTPC->GetLabel());
+       TParticle* p_mcTrack = fMCStack->Particle(p_label);         
+       if (p_mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(p_label))
+           p_primaryFlag = 1;
+
+         //10/01/13. Add a flag to see if it is from material or WD
+         
+         if(fMCStack->IsSecondaryFromWeakDecay(p_label))
+           p_primaryFlag = 2;
+         
+         if(fMCStack->IsSecondaryFromMaterial(p_label))
+           p_primaryFlag = 3;
+
+
+         
+         Int_t p_pdgCode = p_mcTrack->GetPdgCode();
+         p_pidCode = GetPidCode(p_pdgCode);
+         
+         p_ptMC      = p_mcTrack->Pt();
+         
+         p_mother_label = FindPrimaryMotherLabelV0(fMCStack, p_label, 
+                                                 p_mother_steps);
+         if(p_mother_label>0) {
+           TParticle* p_mother = fMCStack->Particle(p_mother_label);
+           p_pdgMother = p_mother->GetPdgCode();
+         }
+       }
+       
+       // negative track
+       const Int_t n_label = TMath::Abs(nTrackTPC->GetLabel());
+       TParticle* n_mcTrack = fMCStack->Particle(n_label);         
+       if (n_mcTrack){
+         
+         if(fMCStack->IsPhysicalPrimary(n_label))
+           n_primaryFlag = 1;
+
+         //10/01/13. Add a flag to see if it is from material or WD
+         
+         if(fMCStack->IsSecondaryFromWeakDecay(n_label))
+           n_primaryFlag = 2;
+         
+         if(fMCStack->IsSecondaryFromMaterial(n_label))
+           n_primaryFlag = 3;
+
+         
+         Int_t n_pdgCode = n_mcTrack->GetPdgCode();
+         n_pidCode = GetPidCode(n_pdgCode);
+         
+         n_ptMC      = n_mcTrack->Pt();
+         
+         n_mother_label = FindPrimaryMotherLabelV0(fMCStack, n_label, 
+                                                 n_mother_steps);
+         if(n_mother_label>0) {
+           TParticle* n_mother = fMCStack->Particle(n_mother_label);
+           n_pdgMother = n_mother->GetPdgCode();
+         }
+       }
+       
+       // Check if V0 is primary = first and the same mother of both partciles
+       if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
+         pdgV0 = p_pdgMother;
+         if(p_mother_steps == 1 && n_mother_steps == 1) {
+           primaryV0 = 1;
+         }
+       }
+      }
+      
+
+      if(fTreeOption) {
+       
+       DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
+       nadded++;
+       
+       // v0 data
+       v0datatpc->p       = esdV0->P();
+       v0datatpc->pt      = esdV0->Pt();
+       v0datatpc->eta     = esdV0->Eta();
+       v0datatpc->phi     = esdV0->Phi();
+       v0datatpc->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
+       v0datatpc->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
+       v0datatpc->dmassG  = deltaInvMassG;
+       v0datatpc->dmassK0 = deltaInvMassK0s;
+       v0datatpc->dmassL  = deltaInvMassL;
+       v0datatpc->dmassAL = deltaInvMassAntiL;
+       v0datatpc->alpha   = alpha;
+       v0datatpc->ptarm   = ptarm;
+       v0datatpc->decayr  = lV0Radius;
+       v0datatpc->decayl  = lV0DecayLength;
+       
+       // New parameters
+       v0datatpc->status  = esdV0->GetOnFlyStatus();
+       v0datatpc->chi2    = esdV0->GetChi2V0();
+       v0datatpc->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
+       // cospt: as I understand this means that the pointing to the vertex
+       // is fine so I remove the dcaxy and dcaz for the V= class
+       v0datatpc->dcadaughters = esdV0->GetDcaV0Daughters();
+       v0datatpc->primary = primaryV0;
+       v0datatpc->pdg     = pdgV0;
+       
+       // positive track
+       v0datatpc->ptrack.p       = pp;
+       v0datatpc->ptrack.pt      = ppt;
+       //        v0data->ptrack.ptcon   = ppt_con;
+       //        v0data->ptrack.tpcchi  = ptpcchi;
+       v0datatpc->ptrack.eta     = peta;
+       v0datatpc->ptrack.phi     = pphi;
+       v0datatpc->ptrack.q       = pcharge;
+       v0datatpc->ptrack.ncl     = pncl;
+       v0datatpc->ptrack.neff    = pneff;
+       v0datatpc->ptrack.dedx    = pdedx;
+
+       v0datatpc->ptrack.isTOFout   = 0;
+       v0datatpc->ptrack.hasTOFtime = 0;
+       v0datatpc->ptrack.isTOFmatched = 0;
+       v0datatpc->ptrack.flength    =   0;
+       v0datatpc->ptrack.ftimetof   = 0;
+       v0datatpc->ptrack.exptoftimeel = 0;
+       v0datatpc->ptrack.exptoftimemu = 0;
+       v0datatpc->ptrack.exptoftimepi = 0;
+       v0datatpc->ptrack.exptoftimeka = 0;
+       v0datatpc->ptrack.exptoftimepr = 0;
+
+
+       v0datatpc->ptrack.dcaxy   = pdcaxy;
+       v0datatpc->ptrack.dcaz    = pdcaz;
+       v0datatpc->ptrack.pid     = p_pidCode;
+       v0datatpc->ptrack.primary = p_primaryFlag;
+       v0datatpc->ptrack.pttrue  = p_ptMC;
+       v0datatpc->ptrack.mother  = p_pdgMother;
+       v0datatpc->ptrack.filter  = filterFlag_p;
+       v0datatpc->ptrack.tpcnclS    = psharedtpcclusters;
+
+
+       // negative track
+       v0datatpc->ntrack.p       = np;
+       v0datatpc->ntrack.pt      = npt;
+       //        v0data->ntrack.ptcon   = npt_con;
+       //        v0data->ntrack.tpcchi  = ntpcchi;
+       v0datatpc->ntrack.eta     = neta;
+       v0datatpc->ntrack.phi     = nphi;
+       v0datatpc->ntrack.q       = ncharge;
+       v0datatpc->ntrack.ncl     = nncl;
+       v0datatpc->ntrack.neff    = nneff;
+       v0datatpc->ntrack.dedx    = ndedx;
+
+       v0datatpc->ntrack.isTOFout   = 0;
+       v0datatpc->ntrack.hasTOFtime = 0;
+       v0datatpc->ntrack.isTOFmatched = 0;
+       v0datatpc->ntrack.flength    =   0;
+       v0datatpc->ntrack.ftimetof   = 0;
+       v0datatpc->ntrack.exptoftimeel = 0;
+       v0datatpc->ntrack.exptoftimemu = 0;
+       v0datatpc->ntrack.exptoftimepi = 0;
+       v0datatpc->ntrack.exptoftimeka = 0;
+       v0datatpc->ntrack.exptoftimepr = 0;
+
+       v0datatpc->ntrack.dcaxy   = ndcaxy;
+       v0datatpc->ntrack.dcaz    = ndcaz;
+       v0datatpc->ntrack.pid     = n_pidCode;
+       v0datatpc->ntrack.primary = n_primaryFlag;
+       v0datatpc->ntrack.pttrue  = n_ptMC;
+       v0datatpc->ntrack.mother  = n_pdgMother;
+       v0datatpc->ntrack.filter  = filterFlag_n;
+       v0datatpc->ntrack.tpcnclS    = nsharedtpcclusters;
+      }
+      
+      // clean up loop over v0
+      
+      delete negPiKF;
+      delete posPiKF;
+      delete posPKF;
+      delete negAPKF;
+
+
+    }
+
+
+
+  }
+  delete myPrimaryVertex;
+
+  
+}
+//__________________________________________________________________________
+void AliAnalysisTaskHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
+  Int_t nv0s = AODevent->GetNumberOfV0s();
+  if(nv0s<1)return;
+  //Int_t     trackmult = 0; // no pt cuts
+  Int_t     nadded    = 0;
+  
+  AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
+  if (!myBestPrimaryVertex) return;
+  
+  
+  if( analysisMode == kGlobalTrk ){
+    if(fV0ArrayGlobalPar)
+      fV0ArrayGlobalPar->Clear();
+  } else if( analysisMode == kTPCTrk ){
+    if(fV0ArrayTPCPar)
+      fV0ArrayTPCPar->Clear();
+  }
+  
+  if( analysisMode==kGlobalTrk ){  
+    
+    
+    // ################################
+    // #### BEGINNING OF V0 CODE ######
+    // ################################
+    // This is the begining of the V0 loop  
+    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+      AliAODv0 *aodV0 = AODevent->GetV0(iV0);
+      if (!aodV0) continue;
+      
+      // common part
+      
+      // AliAODTrack (V0 Daughters)
+      AliAODVertex* vertex = aodV0->GetSecondaryVtx();
+      if (!vertex) {
+       Printf("ERROR: Could not retrieve vertex");
+       continue;
+      }
+      
+      AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
+      AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
+      if (!pTrack || !nTrack) {
+       Printf("ERROR: Could not retrieve one of the daughter track");
+       continue;
+      }
+      
+      // Remove like-sign
+      if (pTrack->Charge() == nTrack->Charge()) {
+       //cout<< "like sign, continue"<< endl;
+       continue;
+      } 
+      
+      // Make sure charge ordering is ok
+      if (pTrack->Charge() < 0) {
+       AliAODTrack* helpTrack = pTrack;
+       pTrack = nTrack;
+       nTrack = helpTrack;
+      } 
+      
+      // Eta cut on decay products
+      if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
+       continue;
+      
+      // Pt cut on decay products
+      if (aodV0->Pt() < fMinPtV0)
+       //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
+       continue;
+      
+
+
+      //check positive tracks
+      UShort_t filterFlag_p = 0;
+      
+      if (fTrackFilterGolden) {  
+       // For AOD068 
+       
+       if(pTrack->TestFilterBit(1024)) {
+         filterFlag_p +=1;
+       }
+      }
+      
+      
+      if (fTrackFilterTPC) {
+       // TPC only cuts is bit 1 according to above macro
+       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+       if(pTrack->TestFilterBit(1)){
+         filterFlag_p +=2;       
+       }
+      }
+
+
+      // if(pTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
+      //       // Minimum number of clusters
+      //   Float_t nCrossedRowsTPC = pTrack->GetTPCClusterInfo(2,1);
+      //   if (nCrossedRowsTPC >= 70) {
+      //         Int_t findable=pTrack->GetTPCNclsF();
+      //         if (findable > 0){ 
+      //           if (nCrossedRowsTPC/findable >= 0.8) 
+      //             filterFlag_p += 16;
+      //         }
+      //       }
+       
+      // }
+      
+
+      
+      //check negative tracks
+      UShort_t filterFlag_n = 0;
+
+      
+      if (fTrackFilterGolden) {
+       
+       // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
+       if(nTrack->TestFilterBit(1024)) {
+         filterFlag_n +=1;
+       }
+      }
+      
+      
+      if (fTrackFilterTPC) {
+       // TPC only cuts is bit 1 according to above macro
+       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+       if(nTrack->TestFilterBit(1)){
+         filterFlag_n +=2;
+       }
+      }
+      
+
+      // if(nTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
+      //       // Minimum number of clusters
+      //   Float_t nCrossedRowsTPC = pTrack->GetTPCClusterInfo(2,1);
+      //   if (nCrossedRowsTPC >= 70) {
+      //         Int_t findable=nTrack->GetTPCNclsF();
+      //         if (findable > 0){ 
+      //           if (nCrossedRowsTPC/findable >= 0.8) filterFlag_n += 16;
+      //         }
+      //       }
+      // }
+      
+      
+      Float_t alpha = aodV0->AlphaV0();
+      Float_t ptarm = aodV0->PtArmV0();
+      // Double_t pVtxPos= v0->PrimaryVtxPosition();      
+      
+      Double_t lV0Radius      = aodV0->RadiusV0();
+      Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
+      
+      Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
+      Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
+      Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
+      Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
+      
+      if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
+        TMath::Abs(deltaInvMassL) > fMassCut &&
+        TMath::Abs(deltaInvMassAntiL) > fMassCut &&
+        TMath::Abs(deltaInvMassG) > fMassCut )
+       continue;
+
+
+      // TODO: Why should these be different? Different mass hypothesis = energy loss
+      // This is not important for us as we focus on the decay products!
+      // Double_t ptK0s        = v0K0sKF.GetPt(); 
+      // Double_t ptL          = v0LambdaKF.GetPt();
+      // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
+      
+      // Extract track information
+      
+      Double_t b[2], cov[3];
+      if(!pTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
+       filterFlag_p += 32; // propagation failed!!!!!
+      
+      Float_t pdcaxy   = b[0];
+      Float_t pdcaz    = b[1];
+      if(!nTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
+       filterFlag_n += 32; // propagation failed!!!!!
+      Float_t ndcaxy   = b[0];
+      Float_t ndcaz    = b[1];
+      
+      Short_t pcharge  = pTrack->Charge();
+      Float_t ppt      = pTrack->Pt();
+      Float_t pp       = pTrack->P(); 
+      Float_t peta     = pTrack->Eta();
+      Float_t pphi     = pTrack->Phi();
+      //       Float_t ptpcchi  = pTrack->Chi2perNDF();
+      
+      AliAODPid* pPid = pTrack->GetDetPid();
+      Short_t pncl     = -10;
+      Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t pdedx    = -10;
+   
+      Bool_t pIsTOFout=kFALSE;
+      Bool_t pHasTOFTime=kTRUE;
+      Bool_t pIsMatched=kFALSE;
+      Float_t plengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
+      Float_t ptimeTOF=-999;
+      Float_t pfexptimeel=-999;
+      Float_t pfexptimemu=-999;
+      Float_t pfexptimepi=-999;
+      Float_t pfexptimeka=-999;
+      Float_t pfexptimepr=-999;  
+
+
+      if(pPid) {
+       pncl     = pPid->GetTPCsignalN();
+       pdedx    = pPid->GetTPCsignal();
+       //TOF
+
+       /*
+       if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
+         Double_t tof[5];
+         pPid->GetIntegratedTimes(tof);
+         pbeta = tof[0]/pPid->GetTOFsignal();
+       }
+       */
+
+
+       if ((pTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+         pIsTOFout=kTRUE;
+       
+
+       if ((pTrack->GetStatus()&AliESDtrack::kTIME)==0)
+         pHasTOFTime=kFALSE;
+       
+
+       if (!(pTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+         pIsMatched=kFALSE;
+       else
+         pIsMatched=kTRUE;
+       
+       ptimeTOF=pPid->GetTOFsignal();
+       
+       Double_t pinttime[5]={0,0,0,0,0}; 
+       pPid->GetIntegratedTimes(pinttime);// Returns the array with integrated times for each particle hypothesis
+       pfexptimeel=pinttime[0];
+       pfexptimemu=pinttime[1];
+       pfexptimepi=pinttime[2];
+       pfexptimeka=pinttime[3];
+       pfexptimepr=pinttime[4];  
+
+
+
+
+
+      }
+      TBits psharedTPC=pTrack->GetTPCSharedMap();
+      Int_t psharedtpcclusters=psharedTPC.CountBits(0)-psharedTPC.CountBits(159);
+
+      TBits nsharedTPC=nTrack->GetTPCSharedMap();
+      Int_t nsharedtpcclusters=nsharedTPC.CountBits(0)-nsharedTPC.CountBits(159);
+
+
+      
+      Short_t ncharge  = nTrack->Charge();
+      Float_t npt      = nTrack->Pt();
+      Float_t np       = nTrack->P(); 
+      Float_t neta     = nTrack->Eta();
+      Float_t nphi     = nTrack->Phi();
+      //       Float_t ntpcchi  = nTrack->Chi2perNDF();
+      
+      AliAODPid* nPid = nTrack->GetDetPid();
+      Short_t nncl     = -10;
+      Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+      Float_t ndedx    = -10;
+      //Float_t nbeta = -99;
+
+      Bool_t nIsTOFout=kFALSE;
+      Bool_t nHasTOFTime=kTRUE;
+      Bool_t nIsMatched=kFALSE;
+      Float_t nlengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
+      Float_t ntimeTOF=-999;
+      Float_t nfexptimeel=-999;
+      Float_t nfexptimemu=-999;
+      Float_t nfexptimepi=-999;
+      Float_t nfexptimeka=-999;
+      Float_t nfexptimepr=-999;  
+
+
+      if(pPid) {
+       nncl     = nPid->GetTPCsignalN();
+       ndedx    = nPid->GetTPCsignal();
+       //TOF
+       /*
+       if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
+         Double_t tof[5];
+         nPid->GetIntegratedTimes(tof);
+         nbeta = tof[0]/nPid->GetTOFsignal();
+         }*/
+
+
+       if ((nTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+         nIsTOFout=kTRUE;
+       
+
+       if ((nTrack->GetStatus()&AliESDtrack::kTIME)==0)
+         nHasTOFTime=kFALSE;
+       
+       if (!(nTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
+         nIsMatched=kFALSE;
+       else
+         nIsMatched=kTRUE;
+       
+       ntimeTOF=nPid->GetTOFsignal();
+       
+       Double_t ninttime[5]={0,0,0,0,0}; 
+       nPid->GetIntegratedTimes(ninttime);// Returns the array with integrated times for each particle hypothesis
+       nfexptimeel=ninttime[0];
+       nfexptimemu=ninttime[1];
+       nfexptimepi=ninttime[2];
+       nfexptimeka=ninttime[3];
+       nfexptimepr=ninttime[4];  
+
+      }
+      
+      Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
+      Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
+      Float_t p_ptMC        = 0;
+      Short_t p_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   p_pdgMother   = 0;
+      Float_t n_ptMC        = 0;
+      Short_t n_pidCode     = 0; // 0 = real data / no mc track!
+      Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
+      Int_t   n_pdgMother   = 0;
+      if(fAnalysisMC) {
+       
+       AliAODMCParticle* p_mother = 0;
+       Int_t p_mother_steps = 0;
+       AliAODMCParticle* n_mother = 0;
+       Int_t n_mother_steps = 0;
+       
+       // positive track
+       const Int_t p_label = TMath::Abs(pTrack->GetLabel());
+       
+       AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
+       if (p_mcTrack){
+         
+         if(p_mcTrack->IsPhysicalPrimary())
+           p_primaryFlag = 1;
+         
+
+
+         Int_t p_pdgCode = p_mcTrack->GetPdgCode();
+         p_pidCode = GetPidCode(p_pdgCode);
+         
+         p_ptMC      = p_mcTrack->Pt();
+         
+         p_mother = FindPrimaryMotherAODV0(p_mcTrack, p_mother_steps);
+         if(p_mother)
+           p_pdgMother = p_mother->GetPdgCode();
+       }
+       
+       // negative track
+       const Int_t n_label = TMath::Abs(pTrack->GetLabel());
+       
+       AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
+       if (n_mcTrack){
+         
+         if(n_mcTrack->IsPhysicalPrimary())
+           n_primaryFlag = 1;
+
+         
+         Int_t n_pdgCode = n_mcTrack->GetPdgCode();
+         n_pidCode = GetPidCode(n_pdgCode);
+         
+         n_ptMC      = n_mcTrack->Pt();
+         
+         n_mother = FindPrimaryMotherAODV0(n_mcTrack, n_mother_steps);
+         if(n_mother)
+           n_pdgMother = n_mother->GetPdgCode();
+       }
+       
+       // Check if V0 is primary = first and the same mother of both partciles
+       if(p_mother && n_mother && p_mother == n_mother) {
+         pdgV0 = p_pdgMother;
+         if(p_mother_steps == 1 && n_mother_steps == 1) {
+           primaryV0 = 1;
+         }
+       }
+      }
+      if(fTreeOption) {
+       
+       DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
+       nadded++;
+       
+       // v0 data
+       v0data->p       = aodV0->P();
+       v0data->pt      = aodV0->Pt();
+       v0data->eta     = aodV0->Eta();
+       v0data->phi     = aodV0->Phi();
+       v0data->pdca    = aodV0->DcaPosToPrimVertex();
+       v0data->ndca    = aodV0->DcaNegToPrimVertex();
+       v0data->dmassG  = deltaInvMassG;
+       v0data->dmassK0 = deltaInvMassK0s;
+       v0data->dmassL  = deltaInvMassL;
+       v0data->dmassAL = deltaInvMassAntiL;
+       v0data->alpha   = alpha;
+       v0data->ptarm   = ptarm;
+       v0data->decayr  = lV0Radius;
+       v0data->decayl  = lV0DecayLength;
+       
+       // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
+       // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
+       
+       // New parameters
+       v0data->status  = aodV0->GetOnFlyStatus();
+       v0data->chi2    = aodV0->Chi2V0();
+       v0data->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
+       // cospt: as I understand this means that the pointing to the vertex
+       // is fine so I remove the dcaxy and dcaz for the V= class
+       v0data->dcav0   = aodV0->DcaV0ToPrimVertex();
+       v0data->dcadaughters = aodV0->DcaV0Daughters();
+       v0data->primary = primaryV0;
+       v0data->pdg     = pdgV0;
+       
+       //hljunggr 
+       v0data->y = aodV0->Y(pdgV0);
+       
+       // positive track
+       v0data->ptrack.p       = pp;
+       v0data->ptrack.pt      = ppt;
+       //        v0data->ptrack.ptcon   = ppt_con;
+       //        v0data->ptrack.tpcchi  = ptpcchi;
+       v0data->ptrack.eta     = peta;
+       v0data->ptrack.phi     = pphi;
+       v0data->ptrack.q       = pcharge;
+       v0data->ptrack.ncl     = pncl;
+       v0data->ptrack.neff    = pneff;
+       v0data->ptrack.dedx    = pdedx;
+
+       v0data->ptrack.isTOFout   = pIsTOFout;
+       v0data->ptrack.hasTOFtime = pHasTOFTime;
+       v0data->ptrack.isTOFmatched = pIsMatched;
+       v0data->ptrack.flength    =   plengthtrack;
+       v0data->ptrack.ftimetof   = ptimeTOF;
+       v0data->ptrack.exptoftimeel = pfexptimeel;
+       v0data->ptrack.exptoftimemu = pfexptimemu;
+       v0data->ptrack.exptoftimepi = pfexptimepi;
+       v0data->ptrack.exptoftimeka = pfexptimeka;
+       v0data->ptrack.exptoftimepr = pfexptimepr;
+
+       v0data->ptrack.dcaxy   = pdcaxy;
+       v0data->ptrack.dcaz    = pdcaz;
+       v0data->ptrack.pid     = p_pidCode;
+       v0data->ptrack.primary = p_primaryFlag;
+       v0data->ptrack.pttrue  = p_ptMC;
+       v0data->ptrack.mother  = p_pdgMother;
+       v0data->ptrack.filter  = filterFlag_p;
+       v0data->ptrack.tpcnclS    = psharedtpcclusters;
+
+
+
+       
+       // negative track
+       v0data->ntrack.p       = np;
+       v0data->ntrack.pt      = npt;
+       //        v0data->ntrack.ptcon   = npt_con;
+       //        v0data->ntrack.tpcchi  = ntpcchi;
+       v0data->ntrack.eta     = neta;
+       v0data->ntrack.phi     = nphi;
+       v0data->ntrack.q       = ncharge;
+       v0data->ntrack.ncl     = nncl;
+       v0data->ntrack.neff    = nneff;
+       v0data->ntrack.dedx    = ndedx;
+
+       v0data->ntrack.isTOFout   = nIsTOFout;
+       v0data->ntrack.hasTOFtime = nHasTOFTime;
+       v0data->ntrack.isTOFmatched = nIsMatched;
+       v0data->ntrack.flength    =   nlengthtrack;
+       v0data->ntrack.ftimetof   = ntimeTOF;
+       v0data->ntrack.exptoftimeel = nfexptimeel;
+       v0data->ntrack.exptoftimemu = nfexptimemu;
+       v0data->ntrack.exptoftimepi = nfexptimepi;
+       v0data->ntrack.exptoftimeka = nfexptimeka;
+       v0data->ntrack.exptoftimepr = nfexptimepr;
+
+       v0data->ntrack.dcaxy   = ndcaxy;
+       v0data->ntrack.dcaz    = ndcaz;
+       v0data->ntrack.pid     = n_pidCode;
+       v0data->ntrack.primary = n_primaryFlag;
+       v0data->ntrack.pttrue  = n_ptMC;
+       v0data->ntrack.mother  = n_pdgMother;
+       v0data->ntrack.filter  = filterFlag_n;
+       v0data->ntrack.tpcnclS    = nsharedtpcclusters;
+       
+      }
+    }//end loop over v0's
+    
+    
+  }else if( analysisMode==kTPCTrk ){  
+    
+    const AliAODVertex*        vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
+    if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
+    
+    // ################################
+    // #### BEGINNING OF V0 CODE ######
+    // ################################
+    // This is the begining of the V0 loop  
+
+
+
+     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+       AliAODv0 *aodV0 = AODevent->GetV0(iV0);
+       if (!aodV0) continue;
+       
+       // common part
+       
+       // AliAODTrack (V0 Daughters)
+       AliAODVertex* vertex = aodV0->GetSecondaryVtx();
+       if (!vertex) {
+        Printf("ERROR: Could not retrieve vertex");
+        continue;
+       }
+       
+       AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
+       AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
+       if (!pTrack || !nTrack) {
+        Printf("ERROR: Could not retrieve one of the daughter track");
+        continue;
+       }
+       
+       
+       // Remove like-sign
+       if (pTrack->Charge() == nTrack->Charge()) {
+        //cout<< "like sign, continue"<< endl;
+        continue;
+       } 
+       
+       // Make sure charge ordering is ok
+       if (pTrack->Charge() < 0) {
+        AliAODTrack* helpTrack = pTrack;
+        pTrack = nTrack;
+        nTrack = helpTrack;
+       } 
+       
+       // Eta cut on decay products
+       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
+        continue;
+
+       cout<<"Eta positive track:"<<pTrack->Eta()<<endl;
+
+       
+       TBits psharedTPC=pTrack->GetTPCSharedMap();
+       Int_t psharedtpcclusters=psharedTPC.CountBits(0)-psharedTPC.CountBits(159);
+       
+       TBits nsharedTPC=nTrack->GetTPCSharedMap();
+       Int_t nsharedtpcclusters=nsharedTPC.CountBits(0)-nsharedTPC.CountBits(159);
+       
+
+       // Pt cut on decay products
+       if (aodV0->Pt() < fMinPt)
+        //     if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
+        continue;
+       
+       //check positive tracks
+       UShort_t filterFlag_p = 0;
+       
+       // TPC only cuts is bit 1 according to above macro
+       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+       if(pTrack->TestFilterBit(128)) {
+        cout<<"este track paso el corte bit 128"<<endl;
+        filterFlag_p +=1;
+       }
+       cout<<"filterFlag_p="<<filterFlag_p<<endl;      
+
+
+
+
+       //check negative tracks
+       UShort_t filterFlag_n = 0;
+       
+       // TPC only cuts is bit 1 according to above macro
+       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
+       if(nTrack->TestFilterBit(128)) {
+        filterFlag_n +=1;
+       }
+       
+       
+       
+       Float_t alpha = aodV0->AlphaV0();
+       Float_t ptarm = aodV0->PtArmV0();
+       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
+       
+       Double_t lV0Radius      = aodV0->RadiusV0();
+       Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
+       
+       Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
+       Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
+       Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
+       Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
+       
+
+       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
+         TMath::Abs(deltaInvMassL) > fMassCut &&
+         TMath::Abs(deltaInvMassAntiL) > fMassCut &&
+         TMath::Abs(deltaInvMassG) > fMassCut )
+        continue;
+       
+       
+       // TODO: Why should these be different? Different mass hypothesis = energy loss
+       // This is not important for us as we focus on the decay products!
+       // Double_t ptK0s        = v0K0sKF.GetPt(); 
+       // Double_t ptL          = v0LambdaKF.GetPt();
+       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
+       
+       // Extract track information
+       
+       Double_t b[2], cov[3];
+       if(!pTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
+        filterFlag_p += 32; // propagation failed!!!!!
+       
+       Float_t pdcaxy   = b[0];
+       Float_t pdcaz    = b[1];
+       if(!nTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
+        filterFlag_n += 32; // propagation failed!!!!!
+       Float_t ndcaxy   = b[0];
+       Float_t ndcaz    = b[1];
+       
+       Short_t pcharge  = pTrack->Charge();
+       Float_t ppt      = pTrack->Pt();
+       Float_t pp       = pTrack->P(); 
+       Float_t peta     = pTrack->Eta();
+       Float_t pphi     = pTrack->Phi();
+       //      Float_t ptpcchi  = pTrack->Chi2perNDF();
+       
+       AliAODPid* pPid = pTrack->GetDetPid();
+       Short_t pncl     = -10;
+       Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+       Float_t pdedx    = -10;
+       Float_t pbeta = -99;
+       if(pPid) {
+        pncl     = pPid->GetTPCsignalN();
+        pdedx    = pPid->GetTPCsignal();
+        //TOF
+        if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
+          Double_t tof[5];
+          pPid->GetIntegratedTimes(tof);
+          pbeta = tof[0]/pPid->GetTOFsignal();
+        }
+       }
+       
+       Short_t ncharge  = nTrack->Charge();
+       Float_t npt      = nTrack->Pt();
+       Float_t np       = nTrack->P(); 
+       Float_t neta     = nTrack->Eta();
+       Float_t nphi     = nTrack->Phi();
+       //      Float_t ntpcchi  = nTrack->Chi2perNDF();
+       
+       AliAODPid* nPid = nTrack->GetDetPid();
+       Short_t nncl     = -10;
+       Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
+       Float_t ndedx    = -10;
+       Float_t nbeta = -99;
+       if(pPid) {
+        nncl     = nPid->GetTPCsignalN();
+        ndedx    = nPid->GetTPCsignal();
+        //TOF
+        if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
+          Double_t tof[5];
+          nPid->GetIntegratedTimes(tof);
+          nbeta = tof[0]/nPid->GetTOFsignal();
+        }
+       }
+       
+       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
+       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
+       Float_t p_ptMC        = 0;
+       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
+       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
+       Int_t   p_pdgMother   = 0;
+       Float_t n_ptMC        = 0;
+       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
+       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
+       Int_t   n_pdgMother   = 0;
+       if(fAnalysisMC) {
+        
+        AliAODMCParticle* p_mother = 0;
+        Int_t p_mother_steps = 0;
+        AliAODMCParticle* n_mother = 0;
+        Int_t n_mother_steps = 0;
+        
+        // positive track
+        const Int_t p_label = TMath::Abs(pTrack->GetLabel());
+        
+        AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
+        if (p_mcTrack){
+          
+          if(p_mcTrack->IsPhysicalPrimary())
+            p_primaryFlag = 1;
+
+
+          
+          Int_t p_pdgCode = p_mcTrack->GetPdgCode();
+          p_pidCode = GetPidCode(p_pdgCode);
+          
+          p_ptMC      = p_mcTrack->Pt();
+          
+          p_mother = FindPrimaryMotherAODV0(p_mcTrack, p_mother_steps);
+          if(p_mother)
+            p_pdgMother = p_mother->GetPdgCode();
+        }
+        
+        // negative track
+        const Int_t n_label = TMath::Abs(pTrack->GetLabel());
+        
+        AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
+        if (n_mcTrack){
+          
+          if(n_mcTrack->IsPhysicalPrimary())
+            n_primaryFlag = 1;
+          
+          Int_t n_pdgCode = n_mcTrack->GetPdgCode();
+          n_pidCode = GetPidCode(n_pdgCode);
+          
+          n_ptMC      = n_mcTrack->Pt();
+          
+          n_mother = FindPrimaryMotherAODV0(n_mcTrack, n_mother_steps);
+          if(n_mother)
+            n_pdgMother = n_mother->GetPdgCode();
+        }
+        
+        // Check if V0 is primary = first and the same mother of both partciles
+        if(p_mother && n_mother && p_mother == n_mother) {
+          pdgV0 = p_pdgMother;
+          if(p_mother_steps == 1 && n_mother_steps == 1) {
+            primaryV0 = 1;
+          }
+        }
+       }
+       
+       if(fTreeOption) {
+        
+        DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
+        nadded++;
+        
+        // v0 data
+        v0datatpc->p       = aodV0->P();
+        v0datatpc->pt      = aodV0->Pt();
+        v0datatpc->eta     = aodV0->Eta();
+        v0datatpc->phi     = aodV0->Phi();
+        v0datatpc->pdca    = aodV0->DcaPosToPrimVertex();
+        v0datatpc->ndca    = aodV0->DcaNegToPrimVertex();
+        v0datatpc->dmassG  = deltaInvMassG;
+        v0datatpc->dmassK0 = deltaInvMassK0s;
+        v0datatpc->dmassL  = deltaInvMassL;
+        v0datatpc->dmassAL = deltaInvMassAntiL;
+        v0datatpc->alpha   = alpha;
+        v0datatpc->ptarm   = ptarm;
+        v0datatpc->decayr  = lV0Radius;
+        v0datatpc->decayl  = lV0DecayLength;
+        // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
+        // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
+        
+        // New parameters
+        v0datatpc->status  = aodV0->GetOnFlyStatus();
+        v0datatpc->chi2    = aodV0->Chi2V0();
+        v0datatpc->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
+        // cospt: as I understand this means that the pointing to the vertex
+        // is fine so I remove the dcaxy and dcaz for the V= class
+        v0datatpc->dcav0   = aodV0->DcaV0ToPrimVertex();
+        v0datatpc->dcadaughters = aodV0->DcaV0Daughters();
+        v0datatpc->primary = primaryV0;
+        v0datatpc->pdg     = pdgV0;
+        
+        // positive track
+        v0datatpc->ptrack.p       = pp;
+        v0datatpc->ptrack.pt      = ppt;
+        //       v0data->ptrack.ptcon   = ppt_con;
+        //       v0data->ptrack.tpcchi  = ptpcchi;
+        v0datatpc->ptrack.eta     = peta;
+        v0datatpc->ptrack.phi     = pphi;
+        v0datatpc->ptrack.q       = pcharge;
+        v0datatpc->ptrack.ncl     = pncl;
+        v0datatpc->ptrack.neff    = pneff;
+        v0datatpc->ptrack.dedx    = pdedx;
+
+       v0datatpc->ptrack.isTOFout   = 0;
+       v0datatpc->ptrack.hasTOFtime = 0;
+       v0datatpc->ptrack.isTOFmatched = 0;
+       v0datatpc->ptrack.flength    =   0;
+       v0datatpc->ptrack.ftimetof   = 0;
+       v0datatpc->ptrack.exptoftimeel = 0;
+       v0datatpc->ptrack.exptoftimemu = 0;
+       v0datatpc->ptrack.exptoftimepi = 0;
+       v0datatpc->ptrack.exptoftimeka = 0;
+       v0datatpc->ptrack.exptoftimepr = 0;
+
+        v0datatpc->ptrack.dcaxy   = pdcaxy;
+        v0datatpc->ptrack.dcaz    = pdcaz;
+        v0datatpc->ptrack.pid     = p_pidCode;
+        v0datatpc->ptrack.primary = p_primaryFlag;
+        v0datatpc->ptrack.pttrue  = p_ptMC;
+        v0datatpc->ptrack.mother  = p_pdgMother;
+        v0datatpc->ptrack.filter  = filterFlag_p;
+        v0datatpc->ptrack.tpcnclS    = psharedtpcclusters;
+        
+        // negative track
+        v0datatpc->ntrack.p       = np;
+        v0datatpc->ntrack.pt      = npt;
+        //       v0data->ntrack.ptcon   = npt_con;
+        //       v0data->ntrack.tpcchi  = ntpcchi;
+        v0datatpc->ntrack.eta     = neta;
+        v0datatpc->ntrack.phi     = nphi;
+        v0datatpc->ntrack.q       = ncharge;
+        v0datatpc->ntrack.ncl     = nncl;
+        v0datatpc->ntrack.neff    = nneff;
+        v0datatpc->ntrack.dedx    = ndedx;
+
+       v0datatpc->ntrack.isTOFout   = 0;
+       v0datatpc->ntrack.hasTOFtime = 0;
+       v0datatpc->ntrack.isTOFmatched = 0;
+       v0datatpc->ntrack.flength    =   0;
+       v0datatpc->ntrack.ftimetof   = 0;
+       v0datatpc->ntrack.exptoftimeel = 0;
+       v0datatpc->ntrack.exptoftimemu = 0;
+       v0datatpc->ntrack.exptoftimepi = 0;
+       v0datatpc->ntrack.exptoftimeka = 0;
+       v0datatpc->ntrack.exptoftimepr = 0;
+
+        v0datatpc->ntrack.dcaxy   = ndcaxy;
+        v0datatpc->ntrack.dcaz    = ndcaz;
+        v0datatpc->ntrack.pid     = n_pidCode;
+        v0datatpc->ntrack.primary = n_primaryFlag;
+        v0datatpc->ntrack.pttrue  = n_ptMC;
+        v0datatpc->ntrack.mother  = n_pdgMother;
+        v0datatpc->ntrack.filter  = filterFlag_n;
+        v0datatpc->ntrack.tpcnclS    = nsharedtpcclusters;
+        
+       }
+     }//end v0's loop
+    
+  }
+  
+
+  
+}
+//_________________________________________
+Float_t  AliAnalysisTaskHighPtDeDx::GetSpherocity(AliESDEvent *esdevent, AliAnalysisFilter* fTrackFilterCuts, Float_t Absetacut, Float_t ptCut, Bool_t useTPCtrack){
+
+  const Int_t nESDTracks = esdevent->GetNumberOfTracks();
+  Int_t ntracksapproved=0;
+  const AliESDVertex *vtxSPD = esdevent->GetPrimaryVertexSPD();
+  if(useTPCtrack)
+    if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) 
+      return -3;
+  
+  
+  for(Int_t iT = 0; iT < nESDTracks; iT++) {
+    
+
+    AliESDtrack* esdTrack = esdevent->GetTrack(iT);
+    Float_t eta=10;
+    Float_t pt=0;    
+
+    AliESDtrack *trackTPC = 0;
+
+
+    if(useTPCtrack){
+
+      //      AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      
+      if(!trackTPC) continue;
+      
+      if(fTrackFilterTPC->IsSelected(trackTPC)==kFALSE){
+       delete trackTPC;
+       continue;
+      }
+    
+      if(trackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParam;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate = false;
+       relate = trackTPC->RelateToVertexTPC(vtxSPD,esdevent->GetMagneticField(),
+                                          kVeryBig,&exParam);
+       if(!relate){
+         delete trackTPC;
+         continue;
+       }
+       trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
+                     exParam.GetCovariance());
+      }
+      else{
+       delete trackTPC;
+       continue;     
+      }
+  
+
+      eta=trackTPC->Eta();
+      pt=trackTPC->Pt();
+
+
+    }else{
+      if(fTrackFilterGolden->IsSelected(esdTrack)==kFALSE)
+       continue;
+      
+      eta=esdTrack->Eta();
+      pt=esdTrack->Pt();
+    }
+
+
+
+
+    if(TMath::Abs(eta)>Absetacut)
+      continue;
+
+    if(pt<ptCut)
+      continue;
+
+    ntracksapproved++;
+
+    if(useTPCtrack)
+      if(trackTPC)
+       delete trackTPC;
+
+  }
+
+  if(ntracksapproved<3)//events with more than 2 primary charged particles
+    return -2;
+  
+  
+  Float_t *pxA=new Float_t[ntracksapproved];
+  Float_t *pyA=new Float_t[ntracksapproved];
+  Float_t sumapt=0;
+  Int_t counter=0;
+  for(Int_t iT = 0; iT < nESDTracks; iT++) {
+    
+    AliESDtrack* esdTrack = esdevent->GetTrack(iT);
+    Float_t eta=10;
+    Float_t pt=0;    
+    Float_t px=0;    
+    Float_t py=0;    
+       
+    AliESDtrack *trackTPC = 0;
+
+    if(useTPCtrack){
+      
+      //      AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      
+      if(!trackTPC) continue;
+      
+      if(fTrackFilterTPC->IsSelected(trackTPC)==kFALSE){
+       delete trackTPC;
+       continue;
+      }
+      
+      
+      if(trackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParam;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate = false;
+       relate = trackTPC->RelateToVertexTPC(vtxSPD,esdevent->GetMagneticField(),
+                                            kVeryBig,&exParam);
+       if(!relate){
+         delete trackTPC;
+         continue;
+       }
+       trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
+                     exParam.GetCovariance());
+      }
+      else{
+       delete trackTPC;
+       continue;     
+      }
+    
+
+      eta=trackTPC->Eta();
+      pt=trackTPC->Pt();
+      px=trackTPC->Px();
+      py=trackTPC->Py();
+
+
+    }else{
+      if(fTrackFilterGolden->IsSelected(esdTrack)==kFALSE)
+       continue;
+      
+      eta=esdTrack->Eta();
+      pt=esdTrack->Pt();
+      px=esdTrack->Px();
+      py=esdTrack->Py();
+
+
+    }
+
+
+
+
+    if(TMath::Abs(eta)>Absetacut)
+      continue;
+
+    if(pt<ptCut)
+      continue;
+
+    pxA[counter]=px;
+    pyA[counter]=py;
+    sumapt+=pt;
+    counter++;
+
+    if(useTPCtrack)
+      if(trackTPC)
+       delete trackTPC;
+    
+    
+  }
+
+  Float_t pFull = 0;
+  Float_t Spherocity = 2;
+  //Getting thrust
+  for(Int_t i = 0; i < 360; ++i){
+    Float_t numerador = 0;
+    Float_t phiparam  = 0;
+    Float_t nx = 0;
+    Float_t ny = 0;
+    phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
+    nx = TMath::Cos(phiparam);            // x component of an unitary vector n
+    ny = TMath::Sin(phiparam);            // y component of an unitary vector n
+    for(Int_t i1 = 0; i1 < ntracksapproved; ++i1){
+      numerador += TMath::Abs(ny * pxA[i1] - nx * pyA[i1]);//product between momentum proyection in XY plane and the unitari vector.
+    }
+    pFull=TMath::Power( (numerador / sumapt),2 );
+    if(pFull < Spherocity)//maximization of pFull
+      {
+       Spherocity = pFull;
+      }
+  }
+
+
+  Float_t spherocity=((Spherocity)*TMath::Pi()*TMath::Pi())/4.0;
+
+  if(pxA){// clean up array memory used for TMath::Sort
+    delete[] pxA; 
+    pxA=0;
+  }
+  if(pyA){// clean up array memory used for TMath::Sort
+    delete[] pyA; 
+    pyA=0;
+  }
+
+  return spherocity;
+
+
+}
+//_________________________________________________________________________________________________________________________________________________
+Float_t  AliAnalysisTaskHighPtDeDx::GetSphericity(AliESDEvent *esdevent, AliAnalysisFilter* fTrackFilterCuts, Float_t Absetacut, Float_t ptCut, Bool_t useTPCtrack){
+
+  const Int_t nESDTracks = esdevent->GetNumberOfTracks();
+  Int_t ntracksapproved=0;
+
+  const AliESDVertex *vtxSPD = esdevent->GetPrimaryVertexSPD();
+  if(useTPCtrack)
+    if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) 
+      return -3;
+
+
+  Double_t s00gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s01gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s10gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s11gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t sumaptg=0;
+
+  for(Int_t iT = 0; iT < nESDTracks; iT++) {
+    
+    AliESDtrack* esdTrack = esdevent->GetTrack(iT);
+    Float_t eta=10;
+    Float_t pt=0;    
+    Float_t px=0;    
+    Float_t py=0;    
+    AliESDtrack *trackTPC = 0;
+
+    if(useTPCtrack){
+      
+      //AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdevent),esdTrack->GetID());
+      
+      if(!trackTPC) continue;
+      
+      if(fTrackFilterTPC->IsSelected(trackTPC)==kFALSE){
+       delete trackTPC;
+       continue;
+      }
+      
+      if(trackTPC->Pt()>0.){
+       // only constrain tracks above threshold
+       AliExternalTrackParam exParam;
+       // take the B-field from the ESD, no 3D fieldMap available at this point
+       Bool_t relate = false;
+       relate = trackTPC->RelateToVertexTPC(vtxSPD,esdevent->GetMagneticField(),
+                                            kVeryBig,&exParam);
+       if(!relate){
+         delete trackTPC;
+         continue;
+       }
+       trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
+                     exParam.GetCovariance());
+      }
+      else{
+       delete trackTPC;
+       continue;     
+      }
+      
+      eta=trackTPC->Eta();
+      pt=trackTPC->Pt();
+      px=trackTPC->Px();
+      py=trackTPC->Py();
+
+      
+    }else{
+      if(fTrackFilterGolden->IsSelected(esdTrack)==kFALSE)
+       continue;
+      
+      eta=esdTrack->Eta();
+      pt=esdTrack->Pt();
+      px=esdTrack->Px();
+      py=esdTrack->Py();
+      
+
+    }
+
+
+    if(TMath::Abs(eta)>Absetacut)
+      continue;
+
+    if(pt<ptCut)
+      continue;
+
+    ntracksapproved++;
+
+    sumaptg+=pt;
+
+    s00gl += (px * px)/ pt;//&&&&&&&&&&&&&&&&&&&&&
+    s01gl += (px * py)/ pt;//&&&&&&&&&&&&&&&&&&&&&
+    s10gl += (py * px)/ pt;//&&&&&&&&&&&&&&&&&&&&&
+    s11gl += (py * py)/ pt;//&&&&&&&&&&&&&&&&&&&&&
+
+    if(useTPCtrack)
+      if(trackTPC)
+       delete trackTPC;
+    
+    
+  }
+
+  if(ntracksapproved<3)
+    return -2;
+  
+
+  //SPHERICITY
+  Double_t S00gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S01gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S10gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S11gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  
+
+  S00gl=s00gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S01gl=s01gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S10gl=s10gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S11gl=s11gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  
+  Float_t sphericity=-2;//&&&&&&&&&&&&&&&&&&&&&
+  
+  Float_t lambda1gl=((S00gl+S11gl)+TMath::Sqrt((S00gl+S11gl)*(S00gl+S11gl)-4*(S00gl*S11gl-S01gl*S01gl)))/2;//&&&&&&&&&&&&&&&&&&&&&
+  Float_t lambda2gl=((S00gl+S11gl)-TMath::Sqrt((S00gl+S11gl)*(S00gl+S11gl)-4*(S00gl*S11gl-S01gl*S01gl)))/2;//&&&&&&&&&&&&&&&&&&&&&
+  if((lambda2gl==0)&&(lambda1gl==0))//&&&&&&&&&&&&&&&&&&&&&
+    sphericity=0;//&&&&&&&&&&&&&&&&&&&&&
+  if(lambda1gl+lambda2gl!=0) //&&&&&&&&&&&&&&&&&&&&&//&&&&&&&&&&&&&&&&&&&&&
+    sphericity=2*TMath::Min( lambda1gl,lambda2gl )/( lambda1gl+lambda2gl );//&&&&&&&&&&&&&&&&&&&&&
+  
+  return sphericity;
+
+
+}
+
+
+//_________________________________________
+Float_t  AliAnalysisTaskHighPtDeDx::GetSpherocityTrue(AliStack *stack, Float_t Absetacut, Float_t ptCut){
+  Int_t nPrim  = stack->GetNprimary();
+  Int_t ntracksapproved=0;
+
+  for(Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {
+    
+    TParticle* trackmc = stack->Particle(iMCTracks);
+    if (!trackmc) continue;
+    // Check if particle is charged, and primary
+    
+    Float_t etamc =trackmc ->Eta();
+    Double_t ptmc=trackmc->Pt();
+    Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
+    
+    Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
+    if(isprimary==0) continue;
+    TParticlePDG* pdgPart =trackmc ->GetPDG();
+    
+    if (pdgPart->Charge() == 0)
+      continue;
+  
+    if (TMath::Abs(etamc) > Absetacut) continue;
+    if(ptmc < ptCut) continue;
+    
+    ntracksapproved++;
+  }
+
+  if(ntracksapproved<3)
+    return -2;
+  
+  
+  Float_t *pxA=new Float_t[ntracksapproved];
+  Float_t *pyA=new Float_t[ntracksapproved];
+  Float_t sumapt=0;
+  Int_t counter=0;
+
+  for(Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {
+    
+    TParticle* trackmc = stack->Particle(iMCTracks);
+    if (!trackmc) continue;
+    // Check if particle is charged, and primary
+    
+    Float_t etamc =trackmc ->Eta();
+    Double_t ptmc=trackmc->Pt();
+    Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
+    
+    Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
+    if(isprimary==0) continue;
+    TParticlePDG* pdgPart =trackmc ->GetPDG();
+    
+    if (pdgPart->Charge() == 0)
+      continue;
+  
+    if (TMath::Abs(etamc) > Absetacut) continue;
+    if(ptmc < ptCut) continue;
+    
+    pxA[counter]=trackmc->Px();
+    pyA[counter]=trackmc->Py();
+    sumapt+=trackmc->Pt();
+    counter++;
+
+  } 
+
+  Float_t pFull = 0;
+  Float_t Spherocity = 2;
+  //Getting thrust
+  for(Int_t i = 0; i < 360; ++i){
+    Float_t numerador = 0;
+    Float_t phiparam  = 0;
+    Float_t nx = 0;
+    Float_t ny = 0;
+    phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
+    nx = TMath::Cos(phiparam);            // x component of an unitary vector n
+    ny = TMath::Sin(phiparam);            // y component of an unitary vector n
+    for(Int_t i1 = 0; i1 < ntracksapproved; ++i1){
+      numerador += TMath::Abs(ny * pxA[i1] - nx * pyA[i1]);//product between momentum proyection in XY plane and the unitari vector.
+    }
+    pFull=TMath::Power( (numerador / sumapt),2 );
+    if(pFull < Spherocity)//maximization of pFull
+      {
+       Spherocity = pFull;
+      }
+  }
+
+
+  Float_t spherocity=((Spherocity)*TMath::Pi()*TMath::Pi())/4.0;
+
+  if(pxA){// clean up array memory used for TMath::Sort
+    delete[] pxA; 
+    pxA=0;
+  }
+  if(pyA){// clean up array memory used for TMath::Sort
+    delete[] pyA; 
+    pyA=0;
+  }
+
+  return spherocity;
+
+
+}
+//_________________________________________________________________________________________________________________________________________________
+Float_t  AliAnalysisTaskHighPtDeDx::GetSphericityTrue(AliStack *stack, Float_t Absetacut, Float_t ptCut){
+
+  Int_t nPrim  = stack->GetNprimary();
+  Int_t ntracksapproved=0;
+
+  Double_t s00gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s01gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s10gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t s11gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t sumaptg=0;
+
+  for(Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {
+    
+    TParticle* trackmc = stack->Particle(iMCTracks);
+    if (!trackmc) continue;
+    // Check if particle is charged, and primary
+    
+    Float_t etamc =trackmc ->Eta();
+    Double_t ptmc=trackmc->Pt();
+    Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
+    
+    Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
+    if(isprimary==0) continue;
+    TParticlePDG* pdgPart =trackmc ->GetPDG();
+    
+    if (pdgPart->Charge() == 0)
+      continue;
+  
+    if (TMath::Abs(etamc) > Absetacut) continue;
+    if(ptmc < ptCut) continue;
+
+    ntracksapproved++;
+
+    sumaptg+=trackmc->Pt();
+    
+    s00gl += (trackmc->Px() * trackmc->Px())/trackmc->Pt();//&&&&&&&&&&&&&&&&&&&&&
+    s01gl += (trackmc->Px() * trackmc->Py())/trackmc->Pt();//&&&&&&&&&&&&&&&&&&&&&
+    s10gl += (trackmc->Py() * trackmc->Px())/trackmc->Pt();//&&&&&&&&&&&&&&&&&&&&&
+    s11gl += (trackmc->Py() * trackmc->Py())/trackmc->Pt();//&&&&&&&&&&&&&&&&&&&&&
+    
+  }
+
+  if(ntracksapproved<3)
+    return -2;
+  
+
+  //SPHERICITY
+  Double_t S00gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S01gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S10gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  Double_t S11gl=0;//&&&&&&&&&&&&&&&&&&&&&
+  
+
+  S00gl=s00gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S01gl=s01gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S10gl=s10gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  S11gl=s11gl/sumaptg;//&&&&&&&&&&&&&&&&&&&&&
+  
+  Float_t sphericity=-2;//&&&&&&&&&&&&&&&&&&&&&
+  
+  Float_t lambda1gl=((S00gl+S11gl)+TMath::Sqrt((S00gl+S11gl)*(S00gl+S11gl)-4*(S00gl*S11gl-S01gl*S01gl)))/2;//&&&&&&&&&&&&&&&&&&&&&
+  Float_t lambda2gl=((S00gl+S11gl)-TMath::Sqrt((S00gl+S11gl)*(S00gl+S11gl)-4*(S00gl*S11gl-S01gl*S01gl)))/2;//&&&&&&&&&&&&&&&&&&&&&
+  if((lambda2gl==0)&&(lambda1gl==0))//&&&&&&&&&&&&&&&&&&&&&
+    sphericity=0;//&&&&&&&&&&&&&&&&&&&&&
+  if(lambda1gl+lambda2gl!=0) //&&&&&&&&&&&&&&&&&&&&&//&&&&&&&&&&&&&&&&&&&&&
+    sphericity=2*TMath::Min( lambda1gl,lambda2gl )/( lambda1gl+lambda2gl );//&&&&&&&&&&&&&&&&&&&&&
+  
+  return sphericity;
+
+
+}
diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.h b/PWGLF/SPECTRA/IdentifiedHighPt/train/AliAnalysisTaskHighPtDeDx.h
new file mode 100755 (executable)
index 0000000..b723150
--- /dev/null
@@ -0,0 +1,170 @@
+#ifndef ALIANALYSISTASKHIGHPTDEDX_H
+#define ALIANALYSISTASKHIGHPTDEDX_H
+
+// ROOT includes
+#include <TList.h>
+#include <TH1.h>
+#include <TTreeStream.h>
+#include <TRandom.h>
+#include <TObject.h>
+
+// AliRoot includes
+#include <AliAnalysisTaskSE.h>
+#include <AliESDEvent.h>
+#include <AliAODEvent.h>
+#include <AliMCEvent.h>
+#include <AliAnalysisFilter.h>
+#include <AliStack.h>
+#include <AliGenEventHeader.h>
+#include <AliVHeader.h>
+#include <AliAODMCParticle.h> 
+#include <AliESDtrackCuts.h>
+#include "DebugClassesMultESA2013.h"
+
+
+
+class AliAnalysisTaskHighPtDeDx : public AliAnalysisTaskSE {
+ public:
+  enum AnalysisMode { kInvalid = -1, kGlobalTrk = 0x1, kTPCTrk = 0x2 }; 
+  AliAnalysisTaskHighPtDeDx();
+  AliAnalysisTaskHighPtDeDx(const char *name);
+
+  virtual ~AliAnalysisTaskHighPtDeDx();
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+
+  Bool_t   GetAnalysisMC() { return fAnalysisMC; }   
+  Double_t GetVtxCut() { return fVtxCut; }   
+  Double_t GetEtaCut() { return fEtaCut; }     
+  Double_t GetMinPt() { return fMinPt; }   
+  Double_t GetMinPtV0() { return fMinPtV0; }
+  Int_t    GetTreeOption() { return fTreeOption; }
+
+  virtual void  SetTrigger1(UInt_t ktriggerInt1) {ftrigBit1 = ktriggerInt1;}
+  virtual void  SetTrigger2(UInt_t ktriggerInt2) {ftrigBit2 = ktriggerInt2;}
+  virtual void  SetTrackFilter(AliAnalysisFilter* trackF) {fTrackFilter = trackF;}
+  virtual void  SetTrackFilterGolden(AliAnalysisFilter* trackF) {fTrackFilterGolden = trackF;}
+  virtual void  SetTrackFilterTPC(AliAnalysisFilter* trackF) {fTrackFilterTPC = trackF;}
+  virtual void  SetProduceVZEROBranch(Bool_t prodvzerob) {fVZEROBranch = prodvzerob;}
+  virtual void  SetProduceTPCBranch(Bool_t prodtpcb) {fTPCBranch = prodtpcb;}
+  virtual void  SetAnalysisType(const char* analysisType) {fAnalysisType = analysisType;}
+  virtual void  SetAnalysisMC(Bool_t isMC) {fAnalysisMC = isMC;}
+  virtual void  SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;}
+  virtual void  SetEtaCut(Double_t etaCut){fEtaCut = etaCut;}
+  virtual void  SetMinPt(Double_t value) {fMinPt = value;}   
+  virtual void  SetMinPtV0(Double_t value) {fMinPtV0 = value;}
+  virtual void  SetMinCent(Float_t minvalc) {fMinCent = minvalc;}
+  virtual void  SetMaxCent(Float_t maxvalc) {fMaxCent = maxvalc;}
+  virtual void  SetLowPtFraction(Double_t value) {fLowPtFraction = value;}   
+  virtual void  SetMassCut(Double_t massCut){fMassCut = massCut;}
+  virtual void  SetTreeOption(Int_t value) {fTreeOption = value;}    
+  virtual void  SetRequireRecV0(Bool_t value) {fRequireRecV0 = value;}
+  virtual void  SetStoreMcIn(Bool_t value) {fStoreMcIn = value;}
+  virtual void  SetAnalysisPbPb(Bool_t isanaPbPb) {fAnalysisPbPb = isanaPbPb;}
+
+ private:
+  virtual Float_t GetVertex(const AliVEvent* event) const;
+  virtual void AnalyzeESD(AliESDEvent* esd); 
+  virtual void AnalyzeAOD(AliAODEvent* aod); 
+  virtual void ProduceArrayTrksESD(AliESDEvent* event, AnalysisMode anamode );
+  virtual void ProduceArrayV0ESD(AliESDEvent* event, AnalysisMode anamode );
+  virtual void ProduceArrayTrksAOD(AliAODEvent* event, AnalysisMode anamode );
+  virtual void ProduceArrayV0AOD(AliAODEvent* event, AnalysisMode anamode );
+  Short_t   GetPidCode(Int_t pdgCode) const;
+  Float_t   GetSpherocity(AliESDEvent* event, AliAnalysisFilter* cuts, Float_t etacut, Float_t ptcut, Bool_t useTPCtrack);
+  Float_t   GetSphericity(AliESDEvent* event, AliAnalysisFilter* cuts, Float_t etacut, Float_t ptcut, Bool_t useTPCtrack);
+  Float_t   GetSpherocityTrue(AliStack *Stack, Float_t etacut, Float_t ptcut);
+  Float_t   GetSphericityTrue(AliStack *Stack, Float_t etacut, Float_t ptcut);
+
+  void      ProcessMCTruthESD();
+  void      ProcessMCTruthAOD(); 
+  void      Sort(TClonesArray* array, Bool_t isMC);
+  Short_t   GetPythiaEventProcessType(Int_t pythiaType);
+  Short_t   GetDPMjetEventProcessType(Int_t dpmJetType);
+  ULong64_t GetEventIdAsLong(AliVHeader* header) const;
+
+  TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
+  Int_t      FindPrimaryMotherLabel(AliStack* stack, Int_t label);
+
+  AliAODMCParticle* FindPrimaryMotherAOD(AliAODMCParticle* startParticle);
+
+  TParticle* FindPrimaryMotherV0(AliStack* stack, Int_t label);
+  Int_t      FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps);
+
+  AliAODMCParticle* FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps);
+
+
+
+  static const Double_t fgkClight;   // Speed of light (cm/ps)
+
+  AliESDEvent* fESD;                  //! ESD object
+  AliAODEvent* fAOD;                  //! AOD object
+  AliMCEvent*  fMC;                   //! MC object
+  AliStack*    fMCStack;              //! MC ESD stack
+  TClonesArray* fMCArray;             //! MC array for AOD
+  AliAnalysisFilter* fTrackFilter;    //  Track Filter, old cuts 2010
+  AliAnalysisFilter* fTrackFilterGolden;    //  Track Filter, set 2010 with golden cuts
+  AliAnalysisFilter* fTrackFilterTPC; // track filter for TPC only tracks
+  TString       fAnalysisType;        //  "ESD" or "AOD"
+  Bool_t        fAnalysisMC;          //  Real(kFALSE) or MC(kTRUE) flag
+  Bool_t        fAnalysisPbPb;        //  true you want to analyze PbPb data, false for pp
+  Bool_t        fVZEROBranch;         //true if you want to store VZERO cells information
+  Bool_t        fTPCBranch;           //tru if you want to produce the TPC branch
+  TRandom*      fRandom;              //! random number generator
+  DeDxEvent*    fEvent;               //! event pointer
+  TClonesArray* fTrackArrayGlobalPar;          //! track array pointer, global tracks
+  TClonesArray* fTrackArrayTPCPar;          //! track array pointer, tpc track parameters
+  TClonesArray* fV0ArrayGlobalPar;             //! V0 array pointer, global tracks
+  TClonesArray* fV0ArrayTPCPar;             //! V0 array pointer, tpc tracks
+  TClonesArray* fTrackArrayMC;        //! MC track array pointer
+  TClonesArray* fVZEROArray;          //! array of the v0 cells.
+
+
+  //
+  // Cuts and options
+  //
+  UInt_t       ftrigBit1;
+  UInt_t       ftrigBit2;
+  Double_t     fVtxCut;             // Vtx cut on z position in cm
+  Double_t     fEtaCut;             // Eta cut used to select particles
+  Double_t     fMinPt;              // Min pt - for histogram limits
+  Double_t     fMinPtV0;              // Min pt - for histogram limits
+  Double_t     fLowPtFraction;      // Fraction of tracks below min pt to keep
+  Double_t     fMassCut;            // Reject all v0 with all dmass > masscut!
+  Int_t        fTreeOption;         // 0: no tree, >0: enable debug tree
+  Float_t      fMinCent; //minimum centrality
+  Float_t      fMaxCent; //maximum centrality
+  Bool_t       fRequireRecV0;       // Require a v0 before updating tree
+                                    // For a spectra analysis we will need to
+                                    // keep track also of the empty events
+  Bool_t       fStoreMcIn;          // Store MC input tracks
+  //
+  // Help variables
+  //
+  Short_t      fMcProcessType;      // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+  Short_t      fTriggeredEventMB;   // 1 = triggered, 0 = not trigged (MC only)
+  Short_t      fVtxStatus;          // -1 = no vtx, 0 = outside cut, 1 = inside cut
+  Float_t      fZvtx;               // z vertex
+  Float_t      fZvtxMC;             // z vertex MC (truth)
+  Int_t        fRun;                // run no
+  ULong64_t    fEventId;            // unique event id
+              
+  //
+  // Output objects
+  //
+  TList*        fListOfObjects;     //! Output list of objects
+  TH1I*         fEvents;            //! No of accepted events
+  TH1I*         fVtx;               //! Event vertex info
+  TH1I*         fVtxMC;             //! Event vertex info for ALL MC events
+  TH1F*         fVtxBeforeCuts;     //! Vertex z dist before cuts
+  TH1F*         fVtxAfterCuts;      //! Vertex z dist after cuts
+  TH1F* fn1;
+  TH1F* fn2;
+  TH1F* fcent;
+  TTree*        fTree;              //! Debug tree 
+
+  ClassDef(AliAnalysisTaskHighPtDeDx, 1);    //Analysis task for high pt analysis 
+};
+
+#endif
diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.cxx b/PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.cxx
new file mode 100755 (executable)
index 0000000..bd2e09e
--- /dev/null
@@ -0,0 +1,291 @@
+// ROOT includes
+#include <TObject.h>
+
+#include "DebugClassesMultESA2013.h"
+
+//_____________________________________________________________________________
+ClassImp(DeDxTrack)
+
+DeDxTrack::DeDxTrack():
+TObject(),
+  p(-1),
+  pt(-1),
+//  ptcon(-1),
+  pttrue(-1),
+//  tpcchi(-1),
+  eta(-999),
+  phi(-999),
+  dedx(-999),
+
+  isTOFout(-1),
+  hasTOFtime(-1),
+  isTOFmatched(-1),
+  flength(-9999),
+  ftimetof(-9999),
+  exptoftimeel(-9999),
+  exptoftimemu(-9999),
+  exptoftimepi(-9999),
+  exptoftimeka(-9999),
+  exptoftimepr(-9999),
+  
+  dcaxy(-999),
+  dcaz(-999),
+  mother(0),
+  q(-999),
+  filter(-999),
+  ncl(-999),
+  neff(-999),
+  pid(-999),
+  primary(-999),
+  order(-1),
+//filterset1(0),
+// filterset2(0),
+//filterset3(0),
+  tpcnclS(0)
+
+{
+  // default constructor
+}
+
+void DeDxTrack::Copy(TObject& object) const
+{
+  TObject::Copy(object);
+
+  DeDxTrack* track = (DeDxTrack*)(&object);
+  if(!track)
+    return;
+  
+  track->p         = p;               
+  track->pt        = pt;              
+  track->pttrue     = pttrue;    
+  track->eta        = eta;       
+  track->phi        = phi;       
+  track->dedx       = dedx;    
+  track->isTOFout     = isTOFout;
+  track->hasTOFtime   = hasTOFtime;
+  track->isTOFmatched = isTOFmatched;
+
+  track->flength      = flength;
+  track->ftimetof     = ftimetof;
+  track->exptoftimeel = exptoftimeel;
+  track->exptoftimemu = exptoftimemu;
+  track->exptoftimepi = exptoftimepi;
+  track->exptoftimeka = exptoftimeka;
+  track->exptoftimepr = exptoftimepr;
+  track->dcaxy      = dcaxy;     
+  track->dcaz       = dcaz;      
+  track->mother     = mother;    
+  track->q         = q;               
+  track->filter     = filter;    
+  track->ncl        = ncl;       
+  track->neff       = neff;      
+  track->pid        = pid;       
+  track->primary    = primary;   
+  track->order      = order;     
+  track->tpcnclS    = tpcnclS;
+}
+
+//_____________________________________________________________________________
+ClassImp(VZEROCell)
+
+VZEROCell::VZEROCell():
+TObject(),
+  cellmult(-1),
+  cellindex(-999)
+
+{
+  // default constructor
+}
+
+void VZEROCell::Copy(TObject& object) const
+{
+  TObject::Copy(object);
+
+  VZEROCell* cellv0 = (VZEROCell*)(&object);
+  if(!cellv0)
+    return;
+  
+  cellv0->cellmult  = cellmult;               
+  cellv0->cellindex = cellindex;              
+}
+
+//_____________________________________________________________________________
+ClassImp(DeDxV0)
+
+DeDxV0::DeDxV0():
+TObject(),
+  p(-1),
+  pt(-1),
+  eta(-999),
+  phi(-999),
+  pdca(-1),
+  ndca(-1),
+  dmassG(-1),
+  dmassK0(-1),
+  dmassL(-1),
+  dmassAL(-1),
+  alpha(-999),
+  ptarm(-999),
+  decayr(-999),
+  decayl(-999),
+  chi2(-1),
+  cospt(-999),
+  dcav0(999),
+  dcadaughters(999),
+  pdg(0),
+  primary(-1),  
+  status(),  
+  ptrack(),
+  ntrack(),
+  y(-999)
+{
+  // default constructor
+}
+
+void DeDxV0::Copy(TObject& object) const
+{
+  TObject::Copy(object);
+
+  DeDxV0* v0 = (DeDxV0*)(&object);
+  if(!v0)
+    return;
+
+  v0->p                    = p;                 
+  v0->pt           = pt;                
+  v0->eta          = eta;               
+  v0->phi          = phi;               
+  v0->pdca         = pdca;              
+  v0->ndca         = ndca;              
+  v0->dmassG       = dmassG;    
+  v0->dmassK0      = dmassK0;   
+  v0->dmassL       = dmassL;    
+  v0->dmassAL      = dmassAL;   
+  v0->alpha        = alpha;     
+  v0->ptarm        = ptarm;     
+  v0->decayr       = decayr;    
+  v0->decayl       = decayl;    
+  v0->chi2         = chi2;              
+  v0->cospt        = cospt;     
+  v0->dcav0        = dcav0;     
+  v0->dcadaughters = dcadaughters;      
+  v0->pdg          = pdg;               
+  v0->primary       = primary;          
+  v0->status       = status;    
+
+  ptrack.Copy(v0->ptrack);      
+  ntrack.Copy(v0->ntrack);        
+  v0->y = y;
+}
+
+//_____________________________________________________________________________
+ClassImp(DeDxTrackMC)
+
+DeDxTrackMC::DeDxTrackMC():
+TObject(),
+  pMC(-1),
+  ptMC(-1),
+  etaMC(-999),
+  phiMC(-999),
+  yMC(-999),
+  qMC(-999),
+  pidMC(-999),
+  orderMC(-1),
+  pdgMC(0)
+{
+  // default constructor
+}
+
+void DeDxTrackMC::Copy(TObject& object) const
+{
+  TObject::Copy(object);
+
+  DeDxTrackMC* trackmc = (DeDxTrackMC*)(&object);
+  if(!trackmc)
+    return;
+  
+  trackmc->pMC     = pMC;             
+  trackmc->ptMC            = ptMC;            
+  trackmc->etaMC    = etaMC;
+  trackmc->phiMC    = phiMC;
+  trackmc->yMC      = yMC;
+  trackmc->qMC      = qMC;
+  trackmc->pidMC    = pidMC;
+  trackmc->orderMC  = orderMC;
+  trackmc->pdgMC    = pdgMC;
+}
+
+//_____________________________________________________________________________
+ClassImp(DeDxEvent)
+
+DeDxEvent::DeDxEvent():
+TObject(),
+  eventid(0),      // unique event id
+  run(-1),         // run number
+  time(-1),        // time of event
+  cent(1000),      // centrality
+  centV0A(1000),      // centrality
+  centZNA(1000),      // centrality
+  centCL1(1000),      // centrality
+  mag(+999),       // magnetic field
+  zvtx(+999),      // rec vertex
+  zvtxMC(+999),    // MC true vertes
+  ptmax(-1),       // Max pt of tracks for this event
+  ptmaxMC(-1),     // Max pt of MC tracks
+  vtxstatus(-2),   // Vtx status (-1=no vtx, 0 = outside, 1 = inside cuts)
+  trackmult(-1),   // Track mult (no cuts)
+  n(-1),           // Number of added tracks 
+  trackmultMC(-1), // MC track mult (primary tracks)
+  nMC(-1),         // MC number of added tracks 
+  process(-2),     // MC process: -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+  trig(-1),        // Was the event triggered
+  pileup(-1),       // Is the event marked as pileup?
+  sphericity(-2),  //|eta|<0.8, pt>0.5, Nch>3, rec
+  spherocity(-2),  //|eta|<0.8, pt>0.5, Nch>3, rec  
+  sphericityTPC(-2),  //|eta|<0.8, pt>0.5, Nch>3, rec
+  spherocityTPC(-2),  //|eta|<0.8, pt>0.5, Nch>3, rec  
+  sphericityMC(-2),  //|eta|<0.8, pt>0.5, Nch>3, true  
+  spherocityMC(-2)  //|eta|<0.8, pt>0.5, Nch>3, true  
+
+{
+  // default constructor
+}
+
+void DeDxEvent::Copy(TObject& object) const
+{
+  TObject::Copy(object);
+
+  DeDxEvent* eventIn = (DeDxEvent*)(&object);
+  if(!eventIn)
+    return;
+
+  eventIn->eventid       = eventid     ; 
+  eventIn->run           = run         ; 
+  eventIn->time          = time        ; 
+  eventIn->cent          = cent        ;
+  eventIn->centV0A       = centV0A     ;
+  eventIn->centZNA       = centZNA     ;
+  eventIn->centCL1       = centCL1     ;
+  eventIn->mag           = mag         ; 
+  eventIn->zvtx          = zvtx        ; 
+  eventIn->zvtxMC        = zvtxMC      ; 
+  eventIn->ptmax         = ptmax       ; 
+  eventIn->ptmaxMC       = ptmaxMC     ; 
+  eventIn->vtxstatus     = vtxstatus   ; 
+  eventIn->trackmult     = trackmult   ; 
+  eventIn->n             = n           ; 
+  eventIn->trackmultMC   = trackmultMC ; 
+  eventIn->nMC           = nMC         ; 
+  eventIn->process       = process     ; 
+  eventIn->trig          = trig        ; 
+  eventIn->pileup        = pileup      ;
+  eventIn->sphericity    = sphericity   ;
+  eventIn->spherocity    = spherocity   ;  
+  eventIn->sphericityTPC = sphericityTPC;
+  eventIn->spherocityTPC = spherocityTPC;  
+  eventIn->sphericityMC  = sphericityMC ; 
+  eventIn->spherocityMC  = spherocityMC ; 
+
+      
+}
diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.h b/PWGLF/SPECTRA/IdentifiedHighPt/train/DebugClassesMultESA2013.h
new file mode 100755 (executable)
index 0000000..da9ffc9
--- /dev/null
@@ -0,0 +1,159 @@
+#ifndef DEBUGCLASSESMULTESA2013_H
+#define DEBUGCLASSESMULTESA2013_H
+
+class DeDxTrack : public TObject
+{
+ public:
+  Float_t   p;
+  Float_t   pt;
+  //  Float_t   ptcon;
+  Float_t   pttrue;
+  //  Float_t   tpcchi;
+  Float_t   eta;
+  Float_t   phi;
+  Float_t   dedx;
+
+  Bool_t    isTOFout;
+  Bool_t    hasTOFtime;
+  Bool_t    isTOFmatched;
+
+  Float_t   flength;
+  Float_t   ftimetof;
+  Float_t   exptoftimeel;
+  Float_t   exptoftimemu;
+  Float_t   exptoftimepi;
+  Float_t   exptoftimeka;
+  Float_t   exptoftimepr;
+
+
+  Float_t   dcaxy;
+  Float_t   dcaz;
+  Int_t     mother; // pdg of mother (can be same particle)
+  Short_t   q;
+  Short_t   filter;
+  Short_t   ncl;
+  Short_t   neff;
+  Short_t   pid;
+  Short_t   primary;  
+  Short_t   order;
+  //Bool_t    filterset1;//TPC  
+  //Bool_t    filterset2;//2010 old
+  //Bool_t    filterset3;//2010 golden
+  Int_t     tpcnclS; //number of shared TPC clusters
+
+  DeDxTrack();
+  void Copy(TObject& object) const;
+
+  ClassDef(DeDxTrack, 2);    // Help class
+};
+//_________________________________________________________
+class VZEROCell : public TObject
+{
+ public:
+
+  Int_t   cellmult;
+  Float_t cellindex;
+  VZEROCell();
+  void Copy(TObject& object) const;
+  
+  ClassDef(VZEROCell, 2);    // Help class
+};
+
+
+//_____________________________________________________________________________
+class DeDxV0 : public TObject
+{
+ public:
+  Float_t   p;
+  Float_t   pt;
+  Float_t   eta;
+  Float_t   phi;
+  Float_t   pdca;     // Distance of Closest Approach for positive track
+  Float_t   ndca;     // Distance of Closest Approach for positive track
+  Float_t   dmassG;
+  Float_t   dmassK0;
+  Float_t   dmassL;
+  Float_t   dmassAL;
+  Float_t   alpha;
+  Float_t   ptarm;
+  Float_t   decayr;
+  Float_t   decayl;
+  // new
+  Float_t   chi2;
+  Float_t   cospt;
+  Float_t   dcav0;
+  Float_t   dcadaughters;
+  Int_t     pdg;
+  Short_t   primary;  
+  Short_t   status;  
+  // old
+  DeDxTrack ptrack;
+  DeDxTrack ntrack;
+  // incl. by hljunggr
+  Float_t   y;
+  DeDxV0();
+  void Copy(TObject& object) const;
+
+  
+  ClassDef(DeDxV0, 3);    // Help class
+};
+
+
+//_____________________________________________________________________________
+class DeDxTrackMC : public TObject
+{
+ public:
+  Float_t pMC;
+  Float_t ptMC;
+  Float_t etaMC;
+  Float_t phiMC;
+  Float_t yMC;
+  Short_t qMC;
+  Short_t pidMC;
+  Short_t orderMC;
+  Int_t   pdgMC;
+
+  DeDxTrackMC();
+  void Copy(TObject& object) const;
+
+  ClassDef(DeDxTrackMC, 2);    // Help class for MC track debug info
+};
+
+//_____________________________________________________________________________
+class DeDxEvent : public TObject
+{
+ public:
+  ULong64_t eventid;     // unique event id
+  Int_t     run;         // run number
+  UInt_t    time;        // time of event
+  Float_t   cent;        // centrality V0A+V0C, default
+  Float_t   centV0A;        // centrality V0A
+  Float_t   centZNA;        // centrality ZNA
+  Float_t   centCL1;        // centrality from number of clusters in layer 1, SPD
+  Float_t   mag;         // magnetic field
+  Float_t   zvtx;        // rec vertex
+  Float_t   zvtxMC;      // MC true vertes
+  Float_t   ptmax;       // Max pt of tracks for this event
+  Float_t   ptmaxMC;     // Max pt of MC tracks
+  Short_t   vtxstatus;   // Vtx status (-1=no vtx, 0 = outside, 1 = inside cuts)
+  Short_t   trackmult;   // Track mult (no cuts)
+  Short_t   n;           // Number of added tracks 
+  Short_t   trackmultMC; // MC track mult (primary tracks)
+  Short_t   nMC;         // MC number of added tracks 
+  Short_t   process;     // MC process: -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+  Short_t   trig;        // 0=untriggered, &1 = MB, &2=V0 AND
+  Short_t   pileup;      // Is the event marked as pileup?
+  Float_t   sphericity;  //|eta|<0.8, pt>0.5, Nch>3  
+  Float_t   spherocity;  //|eta|<0.8, pt>0.5, Nch>3
+  Float_t   sphericityTPC;  //|eta|<0.8, pt>0.5, Nch>3  
+  Float_t   spherocityTPC;  //|eta|<0.8, pt>0.5, Nch>3
+  Float_t   sphericityMC;  //|eta|<0.8, pt>0.5, Nch>3, true  
+  Float_t   spherocityMC;  //|eta|<0.8, pt>0.5, Nch>3, true  
+
+  DeDxEvent();
+  void Copy(TObject& object) const;
+
+  ClassDef(DeDxEvent, 4);    // Help class
+};
+
+#endif
diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/train/macros/AddTaskHighPtDeDx.C b/PWGLF/SPECTRA/IdentifiedHighPt/train/macros/AddTaskHighPtDeDx.C
new file mode 100755 (executable)
index 0000000..8142fcb
--- /dev/null
@@ -0,0 +1,277 @@
+/*
+  Last update: 26/03/2012, vzero branch, in PbPb macro execute the vzero code, a bug was fixed 
+*/
+
+AliAnalysisTask* AddTask(Bool_t AnalysisMC, const Char_t* taskname, Int_t typerun, UInt_t kTriggerInt[], Float_t minc[], Float_t maxc[] )
+{
+  // Creates a pid task and adds it to the analysis manager
+  
+  // Get the pointer to the existing analysis manager via the static
+  //access method
+  //=========================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskHighPtDeDx", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the
+  // analysis manager The availability of MC handler can also be
+  // checked here.
+  // =========================================================================
+  if (!mgr->GetInputEventHandler()) {
+    Error("AddTaskHighPtDeDx", "This task requires an input event handler");
+    return NULL;
+  }  
+  
+
+
+  //
+  // Add track filters, with Golden Cuts
+  //
+  /*
+  AliAnalysisFilter* trackFilterGolden = new AliAnalysisFilter("trackFilter");
+  AliESDtrackCuts* esdTrackCutsGolden = 
+    AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
+  trackFilterGolden->AddCuts(esdTrackCutsGolden);
+  */
+
+  AliAnalysisFilter* trackFilterGolden = new AliAnalysisFilter("trackFilter");
+  AliESDtrackCuts* esdTrackCutsGolden = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
+  esdTrackCutsGolden->SetMinNCrossedRowsTPC(120);
+  esdTrackCutsGolden->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+  esdTrackCutsGolden->SetMaxChi2PerClusterITS(36);
+  esdTrackCutsGolden->SetMaxFractionSharedTPCClusters(0.4);
+  esdTrackCutsGolden->SetMaxChi2TPCConstrainedGlobal(36);
+  esdTrackCutsGolden->SetMaxDCAToVertexXY(3.0);
+
+
+  trackFilterGolden->AddCuts(esdTrackCutsGolden);
+
+
+  
+  //old cuts without golden cut
+  AliAnalysisFilter* trackFilter0 = new AliAnalysisFilter("trackFilter");
+  Bool_t clusterCut=0;
+  Bool_t selPrimaries=kTRUE;
+  AliESDtrackCuts* esdTrackCutsL0 = new AliESDtrackCuts;
+ // TPC  
+  if(clusterCut == 0)  esdTrackCutsL0->SetMinNClustersTPC(70);
+  else if (clusterCut == 1) {
+    esdTrackCutsL0->SetMinNCrossedRowsTPC(70);
+    esdTrackCutsL0->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+  }
+  else {
+    AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut));
+    esdTrackCutsL0->SetMinNClustersTPC(70);
+  }
+  esdTrackCutsL0->SetMaxChi2PerClusterTPC(4);
+  esdTrackCutsL0->SetAcceptKinkDaughters(kFALSE);
+  esdTrackCutsL0->SetRequireTPCRefit(kTRUE);
+  // ITS
+  esdTrackCutsL0->SetRequireITSRefit(kTRUE);
+  esdTrackCutsL0->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                        AliESDtrackCuts::kAny);
+  if(selPrimaries) {
+    // 7*(0.0026+0.0050/pt^1.01)
+    esdTrackCutsL0->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
+  }
+  esdTrackCutsL0->SetMaxDCAToVertexZ(2);
+  esdTrackCutsL0->SetDCAToVertex2D(kFALSE);
+  esdTrackCutsL0->SetRequireSigmaToVertex(kFALSE);
+  esdTrackCutsL0->SetMaxChi2PerClusterITS(1e10);
+  esdTrackCutsL0->SetMaxChi2TPCConstrainedGlobal(1e10);
+
+
+  trackFilter0->AddCuts(esdTrackCutsL0);
+
+  
+  
+  AliAnalysisFilter* trackFilterTPC = new AliAnalysisFilter("trackFilterTPC");
+  AliESDtrackCuts* esdTrackCutsTPC = 
+    AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+  trackFilterTPC->AddCuts(esdTrackCutsTPC);
+
+
+
+  // Create the task and configure it 
+  //========================================================================
+  if(typerun==2){//heavy ion analysis
+    
+    
+    AliAnalysisTaskHighPtDeDx* taskHighPtDeDx[6];
+    for( Int_t i=0; i<6; ++i ){
+      taskHighPtDeDx[i]=0;
+      Char_t TaskName[256]={0};
+      sprintf(TaskName,"%s_%1.0f_%1.0f",taskname,minc[i],maxc[i]);
+      
+      taskHighPtDeDx[i] = new AliAnalysisTaskHighPtDeDx(TaskName);
+      TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+      taskHighPtDeDx[i]->SetAnalysisType(type);
+      taskHighPtDeDx[i]->SetAnalysisMC(AnalysisMC);
+      taskHighPtDeDx[i]->SetAnalysisPbPb(kTRUE);
+      taskHighPtDeDx[i]->SetProduceVZEROBranch(kTRUE);
+      taskHighPtDeDx[i]->SetProduceTPCBranch(kFALSE);
+      taskHighPtDeDx[i]->SetDebugLevel(0);
+      taskHighPtDeDx[i]->SetEtaCut(0.8);
+      taskHighPtDeDx[i]->SetVtxCut(10.0);
+      taskHighPtDeDx[i]->SetMinPt(0.0); // default 2.0
+      taskHighPtDeDx[i]->SetLowPtFraction(0.01); // keep 1% of tracks below min pt
+      taskHighPtDeDx[i]->SetTreeOption(1);
+      taskHighPtDeDx[i]->SetTrigger1(kTriggerInt[0]);
+      taskHighPtDeDx[i]->SetTrigger2(kTriggerInt[1]);
+      taskHighPtDeDx[i]->SetMinCent(minc[i]);
+      taskHighPtDeDx[i]->SetMaxCent(maxc[i]);
+      //Set Filtesr
+      taskHighPtDeDx[i]->SetTrackFilterGolden(trackFilterGolden);
+      taskHighPtDeDx[i]->SetTrackFilter(trackFilter0);
+      taskHighPtDeDx[i]->SetTrackFilterTPC(trackFilterTPC);
+      
+      //V0's cuts
+      taskHighPtDeDx[i]->SetMassCut(0.1);         // def: 0.03
+      taskHighPtDeDx[i]->SetRequireRecV0(kTRUE); // def: kTRUE
+      //taskHighPtDeDx[i]->SetStoreMcIn(kFALSE);     // def: kFALSE
+      taskHighPtDeDx[i]->SetStoreMcIn(kTRUE);     // def: kFALSE
+
+      mgr->AddTask(taskHighPtDeDx[i]);
+      
+    }
+    
+    // Create ONLY the output containers for the data produced by the
+    // task.  Get and connect other common input/output containers via
+    // the manager as below
+    //=======================================================================
+    AliAnalysisDataContainer *cout_hist[6];
+    for( Int_t i=0; i<6; ++i ){
+      
+      cout_hist[i]=0;
+      Char_t outFileName[256]={0};
+      sprintf(outFileName,"%s_Tree_%1.0f_%1.0f.root",taskname,minc[i],maxc[i]);
+      //AliAnalysisDataContainer *cout_hist    = 0;
+      
+      cout_hist[i] = mgr->CreateContainer(Form("output_%1.0f_%1.0f",minc[i],maxc[i]), TList::Class(), AliAnalysisManager::kOutputContainer, outFileName);
+      mgr->ConnectInput (taskHighPtDeDx[i], 0, mgr->GetCommonInputContainer());
+      mgr->ConnectOutput(taskHighPtDeDx[i], 1, cout_hist[i]);
+      
+    }  
+    
+    // Return task pointer at the end
+    return taskHighPtDeDx[0];
+  }
+  if(typerun==3){//pp analysis
+  
+    
+    AliAnalysisTaskHighPtDeDx* taskHighPtDeDx = new AliAnalysisTaskHighPtDeDx("taskHighPtDeDxpp");
+    TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+    taskHighPtDeDx->SetAnalysisType(type);
+    taskHighPtDeDx->SetAnalysisMC(AnalysisMC);
+    taskHighPtDeDx->SetAnalysisPbPb(kFALSE);
+    taskHighPtDeDx->SetProduceVZEROBranch(kTRUE);
+    taskHighPtDeDx->SetProduceTPCBranch(kFALSE);
+    taskHighPtDeDx->SetDebugLevel(0);
+    taskHighPtDeDx->SetEtaCut(0.8);
+    taskHighPtDeDx->SetVtxCut(10.0);
+    taskHighPtDeDx->SetMinPt(0.0); // default 2.0
+    taskHighPtDeDx->SetLowPtFraction(0.01); // keep 1% of tracks below min pt
+    taskHighPtDeDx->SetTreeOption(1);
+    taskHighPtDeDx->SetTrigger1(kTriggerInt[0]);
+    taskHighPtDeDx->SetTrigger2(kTriggerInt[1]);
+    //Set Filtesr
+    taskHighPtDeDx->SetTrackFilterGolden(trackFilterGolden);
+    taskHighPtDeDx->SetTrackFilter(trackFilter0);
+    taskHighPtDeDx->SetTrackFilterTPC(trackFilterTPC);
+    //V0's cuts
+    taskHighPtDeDx->SetMassCut(0.1);         // def: 0.03
+    taskHighPtDeDx->SetRequireRecV0(kTRUE); // def: kTRUE
+    taskHighPtDeDx->SetStoreMcIn(kFALSE);     // def: kFALSE
+
+    mgr->AddTask(taskHighPtDeDx);
+      
+    // Create ONLY the output containers for the data produced by the
+    // task.  Get and connect other common input/output containers via
+    // the manager as below
+    //=======================================================================
+  
+    AliAnalysisDataContainer *cout_histdedx;
+
+    cout_histdedx=0;
+    Char_t outFileName[256]={0};
+    sprintf(outFileName,"%s_Tree.root",taskname);
+    //AliAnalysisDataContainer *cout_hist    = 0;
+    cout_histdedx = mgr->CreateContainer("outputdedx", TList::Class(), AliAnalysisManager::kOutputContainer, outFileName);
+    mgr->ConnectInput (taskHighPtDeDx, 0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput(taskHighPtDeDx, 1, cout_histdedx);
+    // Return task pointer at the end
+    return taskHighPtDeDx;
+    
+
+    
+    
+    
+  }
+
+
+
+ if(typerun==4){//pPb analysis
+   cout<<"<<<<<<<<<<<<<<<<<<<<<<<<<<                 Runing pPb    <<<<<<<<<<<<<<"<<endl;
+    
+    AliAnalysisTaskHighPtDeDx* taskHighPtDeDx = new AliAnalysisTaskHighPtDeDx("taskHighPtDeDxpp");
+    TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+    taskHighPtDeDx->SetAnalysisType(type);
+    taskHighPtDeDx->SetAnalysisMC(AnalysisMC);
+    taskHighPtDeDx->SetAnalysisPbPb(kTRUE);
+    taskHighPtDeDx->SetMinCent(-200);
+    taskHighPtDeDx->SetMaxCent(200);
+    taskHighPtDeDx->SetProduceVZEROBranch(kTRUE);
+    taskHighPtDeDx->SetProduceTPCBranch(kFALSE);
+    taskHighPtDeDx->SetDebugLevel(0);
+    taskHighPtDeDx->SetEtaCut(0.8);
+    taskHighPtDeDx->SetVtxCut(10.0);
+    taskHighPtDeDx->SetMinPt(0.0); // default 2.0
+    taskHighPtDeDx->SetLowPtFraction(0.01); // keep 1% of tracks below min pt
+    taskHighPtDeDx->SetTreeOption(1);
+    taskHighPtDeDx->SetTrigger1(kTriggerInt[0]);
+    taskHighPtDeDx->SetTrigger2(kTriggerInt[1]);
+    //Set Filtesr
+    taskHighPtDeDx->SetTrackFilterGolden(trackFilterGolden);
+    taskHighPtDeDx->SetTrackFilter(trackFilter0);
+    taskHighPtDeDx->SetTrackFilterTPC(trackFilterTPC);
+    //V0's cuts
+    taskHighPtDeDx->SetMassCut(0.1);         // def: 0.03
+    taskHighPtDeDx->SetRequireRecV0(kTRUE); // def: kTRUE
+    taskHighPtDeDx->SetStoreMcIn(kTRUE);     // def: kFALSE
+
+    mgr->AddTask(taskHighPtDeDx);
+      
+    // Create ONLY the output containers for the data produced by the
+    // task.  Get and connect other common input/output containers via
+    // the manager as below
+    //=======================================================================
+  
+    AliAnalysisDataContainer *cout_histdedx;
+
+    cout_histdedx=0;
+    Char_t outFileName[256]={0};
+    sprintf(outFileName,"%s_Tree.root",taskname);
+    //AliAnalysisDataContainer *cout_hist    = 0;
+    cout_histdedx = mgr->CreateContainer("outputdedx", TList::Class(), AliAnalysisManager::kOutputContainer, outFileName);
+    mgr->ConnectInput (taskHighPtDeDx, 0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput(taskHighPtDeDx, 1, cout_histdedx);
+    // Return task pointer at the end
+    return taskHighPtDeDx;
+    
+
+    
+    
+    
+  }
+
+
+
+  
+  
+}