First preliminary tasks for PWGLF production QA (cont'd)
authordelia <delia@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 20:23:37 +0000 (20:23 +0000)
committerdelia <delia@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 20:23:37 +0000 (20:23 +0000)
PWGLF/CMakelibPWGLFQATasks.pkg [new file with mode: 0644]
PWGLF/PWGLFQATasksLinkDef.h [new file with mode: 0644]
PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx [new file with mode: 0644]
PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.h [new file with mode: 0644]
PWGLF/QATasks/AliQAProdMultistrange.cxx [new file with mode: 0644]
PWGLF/QATasks/AliQAProdMultistrange.h [new file with mode: 0644]

diff --git a/PWGLF/CMakelibPWGLFQATasks.pkg b/PWGLF/CMakelibPWGLFQATasks.pkg
new file mode 100644 (file)
index 0000000..5854660
--- /dev/null
@@ -0,0 +1,40 @@
+# -*- mode: CMake -*-
+#-----------------------------------------------------------------------#
+# Package File for PWG2spectra                                          #
+# Author : Johny Jose (johny.jose@cern.ch)                              #
+# Variables Defined :                                                   #
+#                                                                       #
+# SRCS - C++ source files                                               #
+# HDRS - C++ header files                                               #
+# DHDR - ROOT Dictionary Linkdef header file                            #
+# CSRCS - C source files                                                #
+# CHDRS - C header files                                                #
+# EINCLUDE - Include directories                                        #
+# EDEFINE - Compiler definitions                                        #
+# ELIBS - Extra libraries to link                                       #
+# ELIBSDIR - Extra library directories                                  #
+# PACKFFLAGS - Fortran compiler flags for package                       #
+# PACKCXXFLAGS - C++ compiler flags for package                         #
+# PACKCFLAGS - C compiler flags for package                             #
+# PACKSOFLAGS - Shared library linking flags                            #
+# PACKLDFLAGS - Module linker flags                                     #
+# PACKBLIBS - Libraries to link (Executables only)                      #
+# EXPORT - Header files to be exported                                  #
+# CINTHDRS - Dictionary header files                                    #
+# CINTAUTOLINK - Set automatic dictionary generation                    #
+# ARLIBS - Archive Libraries and objects for linking (Executables only) #
+# SHLIBS - Shared Libraries and objects for linking (Executables only)  #
+#-----------------------------------------------------------------------#
+
+set ( SRCS  
+  QATasks/AliAnalysisTaskQAHighPtDeDx.cxx
+  QATasks/AliQAProdMultistrange.cxx
+)
+
+string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
+
+set ( DHDR PWGLFQATasksLinkDef.h)
+
+set ( EXPORT )
+
+set ( EINCLUDE ANALYSIS PWGLF/QATasks STEER/ESD STEER/STEERBase)
diff --git a/PWGLF/PWGLFQATasksLinkDef.h b/PWGLF/PWGLFQATasksLinkDef.h
new file mode 100644 (file)
index 0000000..fe46219
--- /dev/null
@@ -0,0 +1,10 @@
+#ifdef __CINT__
+
+#pragma link off all glols;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnalysisTaskQAHighPtDeDx+;
+#pragma link C++ class AliQAProdMultistrange+;
+
+#endif
diff --git a/PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx b/PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx
new file mode 100644 (file)
index 0000000..d10501b
--- /dev/null
@@ -0,0 +1,1941 @@
+ /*************************************************************************
+* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  * **************************************************************************/
+
+
+#include "AliAnalysisTaskQAHighPtDeDx.h"
+
+// ROOT includes
+#include <TList.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TProfile.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> 
+
+
+// STL includes
+#include <iostream>
+using namespace std;
+
+
+//
+// Responsible:
+// Antonio Ortiz (Lund)
+// Peter Christiansen (Lund)
+//
+
+
+
+
+
+
+
+
+const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
+Float_t magf = -1;
+TF1* cutLow  = new TF1("StandardPhiCutLow",  "0.1/x/x+pi/18.0-0.025", 0, 50);
+TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
+Double_t DeDxMIPMin  = 30;
+Double_t DeDxMIPMax  = 65;
+const Int_t nHists = 9;
+
+Int_t etaLow[nHists]  = {-8, -8, -6, -4, -2, 0, 2, 4, 6};
+Int_t etaHigh[nHists] = { 8, -6, -4, -2,  0, 2, 4, 6, 8};
+
+Int_t nDeltaPiBins   = 60;
+Double_t deltaPiLow  = 40;
+Double_t deltaPiHigh = 100;
+const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
+ClassImp(AliAnalysisTaskQAHighPtDeDx)
+//_____________________________________________________________________________
+AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
+  AliAnalysisTaskSE(name),
+  fESD(0x0),
+  fAOD(0x0),
+  fMC(0x0),
+  fMCStack(0x0),
+  fMCArray(0x0),
+  fTrackFilterGolden(0x0),
+  fTrackFilterTPC(0x0),
+  fCentEst("V0M"),
+  fAnalysisType("ESD"),
+  fAnalysisMC(kFALSE),
+  fAnalysisPbPb(kFALSE),
+  ftrigBit(0x0),
+  fRandom(0x0),
+  fPileUpRej(kFALSE),
+  fVtxCut(10.0),  
+  fEtaCut(0.9),  
+  fMinCent(0.0),
+  fMaxCent(100.0),
+  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),
+  fn1(0x0),
+  fcent(0x0),
+  hMIPVsEta(0x0),
+  pMIPVsEta(0x0),
+  hMIPVsEtaV0s(0x0),
+  pMIPVsEtaV0s(0x0),
+  hPlateauVsEta(0x0),
+  pPlateauVsEta(0x0),
+  hPhi(0x0)
+
+
+{
+  // Default constructor (should not be used)
+  for(Int_t i=0;i<9;++i){
+    hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
+    pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
+    hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
+    pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
+    hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+    pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+    histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
+    histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
+    histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
+    histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
+    histEV0[i]=0;
+
+    for(Int_t pid=0;pid<7;++pid){
+      hMcIn[pid][i]=0;
+      hMcOut[pid][i]=0;
+    }
+
+  }
+  DefineOutput(1, TList::Class());//esto es nuevo
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::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",11,-1,10);
+  fListOfObjects->Add(fn1);
+
+  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);
+
+
+  const Int_t nPtBinsV0s = 25;
+  Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
+                                      1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
+                                      9.0 , 12.0, 15.0, 20.0 };
+  
+
+
+
+  const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
+  
+  const Char_t* LatexEta[nHists] = {
+    "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0", 
+    "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" 
+  };
+  hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4); 
+  //dE/dx vs phi, pions at the MIP
+  fListOfObjects->Add(hPhi);
+
+
+  
+
+  Int_t nPhiBins = 36;
+  
+  for(Int_t i = 0; i < nHists; i++) {
+    
+   
+    hMIPVsPhi[i]      = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(), 
+                                DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+    hMIPVsPhi[i]->Sumw2();
+    
+    pMIPVsPhi[i]     = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
+                                   DeDxMIPMin, DeDxMIPMax);
+    pMIPVsPhi[i]->SetMarkerStyle(20);
+    pMIPVsPhi[i]->SetMarkerColor(1);
+    pMIPVsPhi[i]->SetLineColor(1);
+    pMIPVsPhi[i]->Sumw2();
+
+    fListOfObjects->Add(hMIPVsPhi[i]);
+    fListOfObjects->Add(pMIPVsPhi[i]);
+
+    hPlateauVsPhi[i]  = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
+                                 95-DeDxMIPMax, DeDxMIPMax, 95);
+    hPlateauVsPhi[i]->Sumw2();
+    
+    pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
+                                    DeDxMIPMax, 95);
+    pPlateauVsPhi[i]->SetMarkerStyle(20);
+    pPlateauVsPhi[i]->SetMarkerColor(1);
+    pPlateauVsPhi[i]->SetLineColor(1);
+    pPlateauVsPhi[i]->Sumw2();
+
+    fListOfObjects->Add(hPlateauVsPhi[i]);
+    fListOfObjects->Add(pPlateauVsPhi[i]);
+
+
+    hMIPVsNch[i]      = new TH2D(Form("hMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001, 
+                                 DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+    hMIPVsNch[i]->Sumw2();
+    
+    pMIPVsNch[i] = new TProfile(Form("pMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001, DeDxMIPMin, DeDxMIPMax);
+    pMIPVsNch[i]->SetMarkerStyle(20);
+    pMIPVsNch[i]->SetMarkerColor(1);
+    pMIPVsNch[i]->SetLineColor(1);
+    pMIPVsNch[i]->Sumw2();
+
+    fListOfObjects->Add(hMIPVsNch[i]);
+    fListOfObjects->Add(pMIPVsNch[i]);
+
+
+    histPiV0[i]  = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+    histPiV0[i]->Sumw2();
+    histPV0[i]   = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+
+    fListOfObjects->Add(histPiV0[i]);
+    fListOfObjects->Add(histPV0[i]);
+
+    histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+    histPiTof[i]->Sumw2();
+
+    histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+    histAllCh[i]->Sumw2();
+
+    fListOfObjects->Add(histPiTof[i]);
+    fListOfObjects->Add(histAllCh[i]);
+
+    histEV0[i]   = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+    histEV0[i]->Sumw2();
+    fListOfObjects->Add(histEV0[i]);
+
+
+
+  }
+
+
+  hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+  pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
+  hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+  pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
+  
+  hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
+  pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
+
+  fListOfObjects->Add(hMIPVsEta);
+  fListOfObjects->Add(pMIPVsEta);
+  fListOfObjects->Add(hMIPVsEtaV0s);
+  fListOfObjects->Add(pMIPVsEtaV0s);
+  fListOfObjects->Add(hPlateauVsEta);
+  fListOfObjects->Add(pPlateauVsEta);
+
+
+
+
+
+
+  if (fAnalysisMC) {    
+    for(Int_t i = 0; i < nHists; i++) {
+      for(Int_t pid = 0; pid < 7; pid++) {
+       
+       hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]), 
+                                nPtBinsV0s, ptBinsV0s);
+       fListOfObjects->Add(hMcIn[pid][i]);
+       
+       hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]), 
+                                 nPtBinsV0s, ptBinsV0s);
+       fListOfObjects->Add(hMcOut[pid][i]);
+       
+       
+      }
+    }
+
+    fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
+    fListOfObjects->Add(fVtxMC);
+    
+  }
+  // Post output data.
+  PostData(1, fListOfObjects);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::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
+  
+  fn1->Fill(0);
+  
+  if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
+     ->IsEventSelected() & ftrigBit ){
+    fTriggeredEventMB = 1;  //event triggered as minimum bias
+  }
+
+  // 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; 
+  
+  fn1->Fill(1);
+
+
+  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());
+       }
+      }
+    }
+
+
+  }
+  
+  
+
+  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);
+    
+  fVtxBeforeCuts->Fill(fZvtx);
+  
+  //cut on the z position of vertex
+  if (TMath::Abs(fZvtx) > fVtxCut) {   
+    return;
+  }
+  fn1->Fill(2);
+  
+
+
+
+
+
+  Float_t centrality = -10;
+
+  // only analyze triggered events
+  if(fTriggeredEventMB) {
+    
+    if (fAnalysisType == "ESD"){
+      if(fAnalysisPbPb){
+       AliCentrality *centObject = fESD->GetCentrality();
+       centrality = centObject->GetCentralityPercentile(fCentEst);
+       
+       if((centrality>fMaxCent)||(centrality<fMinCent))return;
+      }
+      fcent->Fill(centrality);
+      fn1->Fill(3);
+      if(fAnalysisMC){
+       if(TMath::Abs(fZvtxMC)<fVtxCut){
+         ProcessMCTruthESD();
+         fVtxMC->Fill(fZvtxMC);
+       }
+      }
+      AnalyzeESD(fESD);
+    } else { // AOD
+      if(fAnalysisPbPb){
+       AliCentrality *centObject = fAOD->GetCentrality();
+       if(centObject){
+         centrality = centObject->GetCentralityPercentile(fCentEst); 
+       }
+       //cout<<"centrality="<<centrality<<endl;
+       if((centrality>fMaxCent)||(centrality<fMinCent))return;
+      }
+      fcent->Fill(centrality);
+      fn1->Fill(3);
+      if(fAnalysisMC){
+       if(TMath::Abs(fZvtxMC)<fVtxCut){
+
+         ProcessMCTruthAOD();
+         fVtxMC->Fill(fZvtxMC);
+       }
+      }
+      AnalyzeAOD(fAOD);
+    }
+  }
+
+  fVtxAfterCuts->Fill(fZvtx);
+
+  
+  
+
+  // Post output data.
+  PostData(1, fListOfObjects);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
+{
+  fRun  = esdEvent->GetRunNumber();
+  fEventId = 0;
+  if(esdEvent->GetHeader())
+    fEventId = GetEventIdAsLong(esdEvent->GetHeader());
+  
+  Bool_t isPileup = esdEvent->IsPileupFromSPD();
+  if(fPileUpRej)
+    if(isPileup)
+      return;
+  fn1->Fill(4);
+
+
+  //  Int_t     event     = esdEvent->GetEventNumberInFile();
+  //UInt_t    time      = esdEvent->GetTimeStamp();
+  //  ULong64_t trigger   = esdEvent->GetTriggerMask();
+  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 );//produce array with global track parameters
+    ProduceArrayV0ESD( esdEvent );//v0's
+
+    
+    fEvents->Fill(1);
+
+
+
+
+  } // end if triggered
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
+{
+  fRun  = aodEvent->GetRunNumber();
+  fEventId = 0;
+  if(aodEvent->GetHeader())
+    fEventId = GetEventIdAsLong(aodEvent->GetHeader());
+   
+  //UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
+  magf      = aodEvent->GetMagneticField();
+
+  //Int_t     trackmult = 0; // no pt cuts
+  //Int_t     nadded    = 0;
+
+  Bool_t isPileup = aodEvent->IsPileupFromSPD();
+  if(fPileUpRej)
+    if(isPileup)
+      return;
+  fn1->Fill(4);
+
+
+
+  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 );
+    ProduceArrayV0AOD( aodEvent );
+
+    fEvents->Fill(1);
+
+
+
+  } // end if triggered
+
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD() 
+{
+  // Fill the special MC histogram with the MC truth info
+  
+  const Int_t nTracksMC = fMCStack->GetNtrack();
+
+  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(chargeMC==0)
+      continue;
+    
+    if (TMath::Abs(trackMC->Eta()) > fEtaCut )
+      continue;
+
+    Double_t etaMC = trackMC->Eta();
+    Int_t pdgCode = trackMC->GetPdgCode();
+    Short_t pidCodeMC = 0;
+    pidCodeMC = GetPidCode(pdgCode);
+    
+    
+    for(Int_t nh = 0; nh < 9; nh++) {
+      
+      if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
+       continue;
+      
+      hMcIn[0][nh]->Fill(trackMC->Pt());
+      hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
+      
+      
+    }
+
+  }//MC track loop
+  
+  
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD() 
+{
+  // Fill the special MC histogram with the MC truth info
+  
+
+  const Int_t nTracksMC = fMCArray->GetEntriesFast();
+  
+  for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
+    
+    AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
+    
+    //Cuts
+    if(!(trackMC->IsPhysicalPrimary()))
+      continue;
+    
+   
+    Double_t chargeMC = trackMC->Charge();
+    if(chargeMC==0)
+      continue;
+
+    
+    if (TMath::Abs(trackMC->Eta()) > fEtaCut )
+      continue;
+    
+    Double_t etaMC = trackMC->Eta();
+    Int_t pdgCode = trackMC->GetPdgCode();
+    Short_t pidCodeMC = 0;
+    pidCodeMC = GetPidCode(pdgCode);
+    
+    //cout<<"pidcode="<<pidCodeMC<<endl;
+    for(Int_t nh = 0; nh < 9; nh++) {
+      
+      if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
+       continue;
+      
+      hMcIn[0][nh]->Fill(trackMC->Pt());
+      hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
+      
+      
+    }
+    
+  }//MC track loop
+  
+  
+}
+
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::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* AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::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* AliAnalysisTaskQAHighPtDeDx::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* AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::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* AliAnalysisTaskQAHighPtDeDx::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 AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
+  
+  const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
+  //Int_t trackmult=0;
+
+
+  Int_t multTPC = 0;
+  
+  //get multiplicity tpc only track cuts
+  for(Int_t iT = 0; iT < nESDTracks; iT++) {
+    
+    AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+    
+    
+    if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
+      continue;
+    
+    //only golden track cuts
+    UInt_t selectDebug = 0;
+    if (fTrackFilterTPC) {
+      selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
+      if (!selectDebug) {
+       //cout<<"this is not a golden track"<<endl;
+       continue;
+      }
+    }
+   
+    multTPC++;
+    
+  }
+
+
+
+  for(Int_t iT = 0; iT < nESDTracks; iT++) {
+    
+    AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+    
+      
+    if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
+      continue;
+    
+    //only golden track cuts
+    UInt_t selectDebug = 0;
+    if (fTrackFilterGolden) {
+      selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
+      if (!selectDebug) {
+       //cout<<"this is not a golden track"<<endl;
+       continue;
+      }
+    }
+    
+    
+    //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
+    Short_t ncl     = esdTrack->GetTPCsignalN();
+
+
+    if(ncl<70)
+      continue;
+    Double_t eta  = esdTrack->Eta();
+    Double_t phi  = esdTrack->Phi();
+    Double_t momentum = esdTrack->P();
+    Float_t  dedx    = esdTrack->GetTPCsignal();
+
+    //cout<<"magf="<<magf<<endl;
+
+    if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
+      continue;
+
+       
+
+
+    //TOF
+    
+    Bool_t IsTOFout=kFALSE;
+    if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+      IsTOFout=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 beta = -99;
+    if ( !IsTOFout ){
+      if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
+       beta = inttime[0] / timeTOF;
+    }
+
+    Short_t pidCode     = 0; 
+    
+    if(fAnalysisMC) {
+      
+      const Int_t label = TMath::Abs(esdTrack->GetLabel());
+      TParticle* mcTrack = fMCStack->Particle(label);      
+      if (mcTrack){
+       
+       Int_t pdgCode = mcTrack->GetPdgCode();
+       pidCode = GetPidCode(pdgCode);
+       
+      }
+      
+    }
+
+
+    if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
+      
+      if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){      
+       if(momentum<0.6&&momentum>0.4){
+         hMIPVsEta->Fill(eta,dedx);
+         pMIPVsEta->Fill(eta,dedx);
+       }
+      }
+      if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+       if(TMath::Abs(beta-1)<0.1){
+         hPlateauVsEta->Fill(eta,dedx);
+         pPlateauVsEta->Fill(eta,dedx); 
+       }
+      }
+    }
+    
+    
+    for(Int_t nh = 0; nh < 9; nh++) {
+      
+      if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
+       continue;
+      
+  
+      if(fAnalysisMC){
+       hMcOut[0][nh]->Fill(esdTrack->Pt());
+       hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
+      }
+
+      histAllCh[nh]->Fill(momentum, dedx);
+      if(beta>1)
+       histPiTof[nh]->Fill(momentum, dedx);
+      
+      if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
+       //Fill  pion MIP, before calibration
+       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){     
+         hMIPVsPhi[nh]->Fill(phi,dedx);
+         pMIPVsPhi[nh]->Fill(phi,dedx);
+         
+         
+         hMIPVsNch[nh]->Fill(multTPC,dedx);
+         pMIPVsNch[nh]->Fill(multTPC,dedx);
+         
+       }
+       
+       //Fill electrons, before calibration
+       if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+         if(TMath::Abs(beta-1)<0.1){
+           hPlateauVsPhi[nh]->Fill(phi,dedx);
+           pPlateauVsPhi[nh]->Fill(phi,dedx);
+         }
+       }
+
+      }
+
+
+    }//end loop over eta intervals
+
+    
+
+    
+    
+  }//end of track loop
+  
+  
+
+  
+}
+//__________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
+
+
+  Int_t nAODTracks = AODevent->GetNumberOfTracks();
+  Int_t multTPC = 0;
+  
+  //get multiplicity tpc only track cuts
+  for(Int_t iT = 0; iT < nAODTracks; iT++) {
+    
+    AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+    
+    if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+      continue;
+    
+
+    if (fTrackFilterTPC) {
+      // TPC only cuts is bit 1
+      if(!aodTrack->TestFilterBit(1))
+       continue;
+    }
+   
+    multTPC++;
+    
+  }
+
+
+  for(Int_t iT = 0; iT < nAODTracks; iT++) {
+    
+    AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+    
+    if (fTrackFilterGolden) {     
+      // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
+      if(!aodTrack->TestFilterBit(1024))
+       continue;
+    }
+    
+    
+    if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+      continue;
+    
+    
+    Double_t eta  = aodTrack->Eta();
+    Double_t phi  = aodTrack->Phi();
+    Double_t momentum = aodTrack->P();
+    
+    if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
+      continue;
+    
+
+    
+    AliAODPid* aodPid = aodTrack->GetDetPid();
+    Short_t ncl     = -10;
+    Float_t dedx    = -10;
+
+    //TOF    
+    Float_t beta = -99;
+
+
+    if(aodPid) {
+      ncl     = aodPid->GetTPCsignalN();
+      dedx    = aodPid->GetTPCsignal();
+      //TOF
+      Bool_t IsTOFout=kFALSE;
+      Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
+      Float_t timeTOF=-999;
+
+      if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+       IsTOFout=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
+
+      
+      if ( !IsTOFout ){
+       if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
+         beta = inttime[0] / timeTOF;
+      }
+      
+    }
+
+
+    if(ncl<70)
+      continue;
+
+
+    Short_t pidCode     = 0; 
+    
+    if(fAnalysisMC) {
+      
+      const Int_t label = TMath::Abs(aodTrack->GetLabel());
+      AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
+      
+      if (mcTrack){
+       
+       Int_t pdgCode = mcTrack->GetPdgCode();
+       pidCode = GetPidCode(pdgCode);
+       
+      }
+      
+    }
+
+
+
+
+    //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
+       //continue;
+    
+    //etaLow
+    //etaHigh
+    if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
+      if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){      
+       if(momentum<0.6&&momentum>0.4){
+         hMIPVsEta->Fill(eta,dedx);
+         pMIPVsEta->Fill(eta,dedx);
+       }
+      }
+      if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+       if(TMath::Abs(beta-1)<0.1){
+         hPlateauVsEta->Fill(eta,dedx);
+         pPlateauVsEta->Fill(eta,dedx); 
+       }
+      }
+    }
+    
+    
+    for(Int_t nh = 0; nh < 9; nh++) {
+      
+      if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
+       continue;
+      
+      
+      if(fAnalysisMC){
+       hMcOut[0][nh]->Fill(aodTrack->Pt());
+       hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
+      }
+      
+      histAllCh[nh]->Fill(momentum, dedx);
+      if(beta>1)
+       histPiTof[nh]->Fill(momentum, dedx);
+      
+      if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
+       //Fill  pion MIP, before calibration
+       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){     
+         hMIPVsPhi[nh]->Fill(phi,dedx);
+         pMIPVsPhi[nh]->Fill(phi,dedx);
+         
+         
+         hMIPVsNch[nh]->Fill(multTPC,dedx);
+         pMIPVsNch[nh]->Fill(multTPC,dedx);
+         
+       }
+       
+       //Fill electrons, before calibration
+       if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+         if(TMath::Abs(beta-1)<0.1){
+           hPlateauVsPhi[nh]->Fill(phi,dedx);
+           pPlateauVsPhi[nh]->Fill(phi,dedx);
+         }
+       }
+
+      }
+
+    }//end loop over eta intervals
+    
+    
+
+    
+    
+  }//end of track loop
+  
+  
+  
+  
+  
+}
+//__________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
+  Int_t nv0s = ESDevent->GetNumberOfV0s();
+  /*
+  if(nv0s<1){
+    return;
+    }*/
+  
+  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);
+
+
+    
+  //
+  // LOOP OVER V0s, K0s, L, AL
+  //
+  
+  
+  for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+    
+    // This is the begining of the V0 loop  
+    AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
+    if (!esdV0) continue;
+    
+    //check onfly status
+    if(esdV0->GetOnFlyStatus()!=0)
+      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;
+    
+    
+    // Check if switch does anything!
+    Bool_t isSwitched = kFALSE;
+    if (pTrack->GetSign() < 0) { // switch
+      
+      isSwitched = kTRUE;
+      AliESDtrack* helpTrack = nTrack;
+      nTrack = pTrack;
+      pTrack = helpTrack;
+    }  
+    
+    //Double_t  lV0Position[3];
+    //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[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 dmassG     = v0GKF.GetMass();
+    Double_t dmassK     = v0K0sKF.GetMass()-0.498;
+    Double_t dmassL     = v0LambdaKF.GetMass()-1.116;
+    Double_t dmassAL    = v0AntiLambdaKF.GetMass()-1.116;
+
+
+    for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
+
+
+      switch(case_v0){
+      case 0:{    
+       
+       Bool_t fillPos = kFALSE;
+       Bool_t fillNeg = kFALSE;
+       
+       if(dmassG < 0.1)
+         continue;
+       
+       if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
+         continue;
+       }
+       
+       if(dmassL<0.01){
+         fillPos = kTRUE;
+       }
+       if(dmassAL<0.01) {
+         if(fillPos)
+           continue;
+         fillNeg = kTRUE;
+       }
+       
+       if(dmassK<0.01) {
+         if(fillPos||fillNeg)
+           continue;
+         fillPos = kTRUE;
+         fillNeg = kTRUE;
+       }
+       
+       
+       for(Int_t j = 0; j < 2; j++) {
+         
+         AliESDtrack* track = 0;
+         
+         if(j==0) {
+           
+           if(fillNeg)
+             track = nTrack;
+           else
+             continue;
+         } else {
+           
+           if(fillPos)
+             track = pTrack;
+           else
+             continue;
+         }
+         
+         if(track->GetTPCsignalN()<=70)continue;
+         Double_t phi     = track->Phi();
+         
+         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+           continue;
+         
+         
+         //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
+         //    continue;
+         
+         Double_t eta  = track->Eta();
+         Double_t momentum = track->Pt();
+         Double_t dedx = track->GetTPCsignal();
+         
+         if(fillPos&&fillNeg){
+           
+           
+           if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
+             if(momentum<0.6&&momentum>0.4){
+               hMIPVsEtaV0s->Fill(eta,dedx);
+               pMIPVsEtaV0s->Fill(eta,dedx);
+             }
+           }
+           
+           
+         }
+         
+         for(Int_t nh = 0; nh < nHists; nh++) {
+           
+           
+           
+           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+             continue;
+           
+           if(fillPos&&fillNeg){
+             
+             histPiV0[nh]->Fill(momentum, dedx);           
+             
+           }
+           else{
+             
+             histPV0[nh]->Fill(momentum, dedx);
+             
+           }
+           
+         }
+                 
+       }//end loop over two tracks
+
+      };
+       break;
+
+      case 1:{//gammas    
+       
+       Bool_t fillPos = kFALSE;
+       Bool_t fillNeg = kFALSE;
+       
+
+
+       if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
+         if(dmassG<0.01 && dmassG>0.0001) {
+           
+
+           if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
+             fillPos = kTRUE;
+           if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
+             fillNeg = kTRUE;
+           
+         } else {
+           continue;
+         }
+       }
+
+       //cout<<"fillPos="<<fillPos<<endl;
+
+       if(fillPos == kTRUE && fillNeg == kTRUE)      
+         continue;
+       
+       
+       AliESDtrack* track = 0;
+       if(fillNeg)
+         track = nTrack;
+       else if(fillPos)
+         track = pTrack;
+       else
+         continue;
+
+       Double_t dedx  = track->GetTPCsignal();
+       Double_t eta  = track->Eta();
+       Double_t phi  = track->Phi();
+       Double_t momentum = track->P();
+
+       if(track->GetTPCsignalN()<=70)continue;
+
+       if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+         continue;
+       
+       for(Int_t nh = 0; nh < nHists; nh++) {
+      
+         if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+           continue;
+
+         histEV0[nh]->Fill(momentum, dedx);
+
+       }
+       
+      };
+       break;
+       
+       
+      }//end switch
+      
+    }//end loop over case V0
+
+    
+    // clean up loop over v0
+    
+    delete negPiKF;
+    delete posPiKF;
+    delete posPKF;
+    delete negAPKF;
+    
+    
+    
+  }
+  
+  
+  delete myPrimaryVertex;
+  
+  
+}
+//__________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
+  Int_t nv0s = AODevent->GetNumberOfV0s();
+  /*
+  if(nv0s<1){
+    return;
+    }*/
+  
+  AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
+  if (!myBestPrimaryVertex) 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;
+    
+    
+    //check onfly status
+    if(aodV0->GetOnFlyStatus()!=0)
+      continue;
+
+    // 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;
+    
+    
+    Double_t dmassG  = aodV0->InvMass2Prongs(0,1,11,11);
+    Double_t dmassK  = aodV0->MassK0Short()-0.498;
+    Double_t dmassL  = aodV0->MassLambda()-1.116;
+    Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
+    
+    for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
+      
+      
+      switch(case_v0){
+      case 0:{   
+       Bool_t fillPos = kFALSE;
+       Bool_t fillNeg = kFALSE;
+       
+       
+       if(dmassG < 0.1)
+         continue;
+       
+       if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
+         continue;
+       }
+       
+       if(dmassL<0.01){
+         fillPos = kTRUE;
+       }
+       if(dmassAL<0.01) {
+         if(fillPos)
+           continue;
+         fillNeg = kTRUE;
+       }
+       
+       if(dmassK<0.01) {
+         if(fillPos||fillNeg)
+           continue;
+         fillPos = kTRUE;
+         fillNeg = kTRUE;
+       }
+       
+       
+       
+       for(Int_t j = 0; j < 2; j++) {
+         
+         AliAODTrack* track = 0;
+         
+         if(j==0) {
+           
+           if(fillNeg)
+             track = nTrack;
+           else
+             continue;
+         } else {
+           
+           if(fillPos)
+             track = pTrack;
+           else
+             continue;
+         }
+         
+         if(track->GetTPCsignalN()<=70)continue;
+         
+         Double_t phi     = track->Phi();
+         
+         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+           continue;
+         
+         
+         //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
+         //    continue;
+         
+         Double_t eta  = track->Eta();
+         Double_t momentum = track->Pt();
+         Double_t dedx = track->GetTPCsignal();
+         
+         if(fillPos&&fillNeg){
+           
+           
+           if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
+             if(momentum<0.6&&momentum>0.4){
+               hMIPVsEtaV0s->Fill(eta,dedx);
+               pMIPVsEtaV0s->Fill(eta,dedx);
+             }
+           }
+           
+           
+         }
+         
+         for(Int_t nh = 0; nh < nHists; nh++) {
+           
+           
+           
+           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+             continue;
+           
+           if(fillPos&&fillNeg){
+             
+             histPiV0[nh]->Fill(momentum, dedx);           
+             
+           }
+           else{
+             
+             histPV0[nh]->Fill(momentum, dedx);
+             
+           }
+           
+         }
+         
+
+       }//end loop over two tracks
+      };
+       break;
+       
+      case 1:{//gammas    
+       
+       Bool_t fillPos = kFALSE;
+       Bool_t fillNeg = kFALSE;
+       
+       if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
+         if(dmassG<0.01 && dmassG>0.0001) {
+           
+           if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
+             fillPos = kTRUE;
+           if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
+             fillNeg = kTRUE;
+           
+         } else {
+           continue;
+         }
+       }
+
+       
+       if(fillPos == kTRUE && fillNeg == kTRUE)      
+         continue;
+       
+       
+       AliAODTrack* track = 0;
+       if(fillNeg)
+         track = nTrack;
+       else if(fillPos)
+         track = pTrack;
+       else
+         continue;
+
+       Double_t dedx  = track->GetTPCsignal();
+       Double_t eta  = track->Eta();
+       Double_t phi  = track->Phi();
+       Double_t momentum = track->P();
+
+       if(track->GetTPCsignalN()<=70)continue;
+
+       if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+         continue;
+       
+       for(Int_t nh = 0; nh < nHists; nh++) {
+      
+         if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+           continue;
+
+         histEV0[nh]->Fill(momentum, dedx);
+
+       }
+       
+      };
+       break;
+
+
+      }//end switch
+    }//end loop over V0s cases
+    
+  }//end loop over v0's
+  
+  
+  
+  
+}
+Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t   mag, TF1* phiCutLow, TF1* phiCutHigh)
+{
+  if(pt < 2.0)
+    return kTRUE;
+  
+  //Double_t phi = track->Phi();
+  if(mag < 0)    // for negatve polarity field
+    phi = TMath::TwoPi() - phi;
+  if(q < 0) // for negatve charge
+    phi = TMath::TwoPi()-phi;
+  
+  phi += TMath::Pi()/18.0; // to center gap in the middle
+  phi = fmod(phi, TMath::Pi()/9.0);
+    
+  if(phi<phiCutHigh->Eval(pt) 
+     && phi>phiCutLow->Eval(pt))
+  return kFALSE; // reject track
+
+  hPhi->Fill(pt, phi);
+
+  return kTRUE;
+}
diff --git a/PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.h b/PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.h
new file mode 100644 (file)
index 0000000..29a178c
--- /dev/null
@@ -0,0 +1,175 @@
+#ifndef ALIANALYSISTASKQAHIGHPTDEDX_H
+#define ALIANALYSISTASKQAHIGHPTDEDX_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id$ */
+
+
+// ROOT includes
+#include <TList.h>
+#include <TH1.h>
+#include <TProfile.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>
+
+
+
+class AliAnalysisTaskQAHighPtDeDx : public AliAnalysisTaskSE {
+ public:
+  //AliAnalysisTaskQAHighPtDeDx();
+  //AliAnalysisTaskQAHighPtDeDx(const char *name);
+  //virtual ~AliAnalysisTaskQAHighPtDeDx();
+  AliAnalysisTaskQAHighPtDeDx(const char *name="<default name>");
+  virtual ~AliAnalysisTaskQAHighPtDeDx() { /*if (fOutputList) delete fOutputList;*/}//;
+
+  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; }   
+  //Int_t    GetTreeOption() { return fTreeOption; }  
+
+  virtual void  SetTrigger(UInt_t ktriggerInt) {ftrigBit = ktriggerInt;}
+  virtual void  SetTrackFilterGolden(AliAnalysisFilter* trackF) {fTrackFilterGolden = trackF;}
+  virtual void  SetTrackFilterTPC(AliAnalysisFilter* trackF) {fTrackFilterTPC = trackF;}
+  virtual void  SetCentralityEstimator(const char * centEst) {fCentEst = centEst;}
+  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  SetPileUpRej(Bool_t isrej) {fPileUpRej = isrej;}   
+  virtual void  SetMinCent(Float_t minvalc) {fMinCent = minvalc;}
+  virtual void  SetMaxCent(Float_t maxvalc) {fMaxCent = maxvalc;}
+  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);
+  virtual void ProduceArrayV0ESD(AliESDEvent* event);
+  virtual void ProduceArrayTrksAOD(AliAODEvent* event);
+  virtual void ProduceArrayV0AOD(AliAODEvent* event);
+  Short_t   GetPidCode(Int_t pdgCode) const;
+  void      ProcessMCTruthESD();
+  void      ProcessMCTruthAOD(); 
+
+  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);
+  Bool_t PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t   mag, TF1* phiCutLow, TF1* phiCutHigh);
+
+
+  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* fTrackFilterGolden;    //  Track Filter, set 2010 with golden cuts
+  AliAnalysisFilter* fTrackFilterTPC; // track filter for TPC only tracks
+  TString       fCentEst;             // V0A , V0M, 
+  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
+  UInt_t       ftrigBit;
+  TRandom*      fRandom;              //! random number generator
+  Bool_t        fPileUpRej;           // kTRUE is pile-up is rejected
+
+
+  //
+  // Cuts and options
+  //
+
+  Double_t     fVtxCut;             // Vtx cut on z position in cm
+  Double_t     fEtaCut;             // Eta cut used to select particles
+  Float_t      fMinCent; //minimum centrality
+  Float_t      fMaxCent; //maximum centrality
+  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
+  TH1F*         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* fcent;
+
+  TH2D *hMIPVsEta;
+  TProfile *pMIPVsEta;
+  TH2D *hMIPVsEtaV0s;
+  TProfile *pMIPVsEtaV0s;
+  TH2D *hPlateauVsEta;
+  TProfile *pPlateauVsEta;
+  TH2D *hPhi;
+
+  TH2D     *hMIPVsNch[9];
+  TProfile *pMIPVsNch[9];
+
+  TH2D     *hMIPVsPhi[9];
+  TProfile *pMIPVsPhi[9];
+  TH2D     *hPlateauVsPhi[9];
+  TProfile *pPlateauVsPhi[9];
+
+  TH2D* histPiV0[9];
+  TH2D* histPV0[9];
+
+  TH2D* histAllCh[9];
+  TH2D* histPiTof[9];
+  TH2D* histEV0[9];
+  TH1D* hMcIn[7][9];
+  TH1D* hMcOut[7][9];
+
+
+
+  //TTree*        fTree;              //! Debug tree 
+
+  ClassDef(AliAnalysisTaskQAHighPtDeDx, 1);    //Analysis task for high pt analysis 
+};
+
+#endif
diff --git a/PWGLF/QATasks/AliQAProdMultistrange.cxx b/PWGLF/QATasks/AliQAProdMultistrange.cxx
new file mode 100644 (file)
index 0000000..e3e9f2f
--- /dev/null
@@ -0,0 +1,901 @@
+/**************************************************************************
+ *  Authors : Antonin Maire, Boris Hippolyte                              *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//            AliQAProdMultistrange class
+//
+//            Origin AliAnalysisTaskCheckCascade which has four roles :
+//              1. QAing the Cascades from ESD and AOD
+//                 Origin:  AliAnalysisTaskESDCheckV0 by Boris Hippolyte Nov2007, hippolyt@in2p3.fr
+//              2. Prepare the plots which stand as raw material for yield extraction (wi/wo PID)
+//              3. Supply an AliCFContainer meant to define the optimised topological selections
+//              4. Rough azimuthal correlation study (Eta, Phi)
+//              Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
+//              Modified :           A.Maire Mar2010 
+//
+//              Adapted to PbPb analysis: M. Nicassio, maria.nicassio@ba.infn.it
+//               Feb-August2011
+//                - Physics selection moved to the run.C macro
+//                - Centrality selection added (+ setters) 
+//                - setters added (vertex range)
+//                - histo added and histo/container binning changed 
+//                - protection in the destructor for CAF usage          
+//                - AliWarning disabled
+//                - automatic settings for PID
+//               September2011
+//                - proper time histos/container added (V0 and Cascades)
+//               November2011
+//                - re-run V0's and cascade's vertexers (SetCuts instead of SetDefaultCuts!!)
+//               Genuary2012 
+//                - AOD analysis part completed 
+//               March2012
+//                - min number of TPC clusters for track selection as a parameter       
+//               July2012
+//                - cut on min pt for daughter tracks as a parameter (+control histos)
+//               August2012 
+//                - cut on pseudorapidity for daughter tracks as a parameter (+control histos for Xi-)
+//-----------------------------------------------------------------
+
+class TTree;
+class TParticle;
+class TVector3;
+
+class AliESDVertex;
+class AliAODVertex;
+class AliESDv0;
+class AliAODv0;
+
+#include <Riostream.h>
+#include "THnSparse.h"
+#include "TVector3.h"
+#include "TMath.h"
+
+#include "AliLog.h"
+#include "AliCentrality.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliESDtrackCuts.h"
+#include "AliPIDResponse.h"
+
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliESDInputHandler.h" 
+#include "AliAODInputHandler.h"
+#include "AliCFContainer.h"
+#include "AliMultiplicity.h"
+
+#include "AliESDcascade.h"
+#include "AliAODcascade.h"
+#include "AliAODTrack.h"
+
+#include "AliQAProdMultistrange.h"
+
+ClassImp(AliQAProdMultistrange)
+
+
+
+//________________________________________________________________________
+AliQAProdMultistrange::AliQAProdMultistrange() 
+  : AliAnalysisTaskSE(), 
+    fAnalysisType               ("ESD"), 
+    fESDtrackCuts               (0),
+    fCollidingSystem            ("PbPb"),
+    fPIDResponse                (0),
+    fkSDDSelectionOn            (kTRUE),
+    fkQualityCutZprimVtxPos     (kTRUE),
+    fkQualityCutNoTPConlyPrimVtx(kTRUE),
+    fkQualityCutTPCrefit        (kTRUE),
+    fkQualityCutnTPCcls         (kTRUE),
+    fkQualityCutPileup          (kTRUE),
+    fwithSDD                    (kTRUE),
+    fMinnTPCcls                 (0),  
+    fCentrLowLim                (0),
+    fCentrUpLim                 (0),
+    fCentrEstimator             (0),
+    fkUseCleaning               (0),
+    fVtxRange                   (0),
+    fMinPtCutOnDaughterTracks   (0),
+    fEtaCutOnDaughterTracks     (0),
+
+    
+    fCFContCascadeCuts(0)
+    
+
+{
+  // Dummy Constructor
+}
+
+
+//________________________________________________________________________
+AliQAProdMultistrange::AliQAProdMultistrange(const char *name) 
+  : AliAnalysisTaskSE(name),
+    fAnalysisType               ("ESD"), 
+    fESDtrackCuts               (0),
+    fCollidingSystem            ("PbPb"),
+    fPIDResponse                (0),
+    fkSDDSelectionOn            (kTRUE),
+    fkQualityCutZprimVtxPos     (kTRUE),
+    fkQualityCutNoTPConlyPrimVtx(kTRUE),
+    fkQualityCutTPCrefit        (kTRUE),
+    fkQualityCutnTPCcls         (kTRUE),
+    fkQualityCutPileup          (kTRUE),
+    fwithSDD                    (kTRUE),
+    fMinnTPCcls                 (0),
+    fCentrLowLim                (0),
+    fCentrUpLim                 (0),
+    fCentrEstimator             (0),
+    fkUseCleaning               (0),
+    fVtxRange                   (0),
+    fMinPtCutOnDaughterTracks   (0),
+    fEtaCutOnDaughterTracks     (0),
+
+
+    fCFContCascadeCuts(0)
+    
+
+{
+  // Constructor
+  // Output slot #0 writes into a TList container (Cascade)
+  DefineOutput(1, AliCFContainer::Class());
+
+  AliLog::SetClassDebugLevel("AliQAProdMultistrange",1); // MN this should (?) enable only AliFatal
+}
+
+
+AliQAProdMultistrange::~AliQAProdMultistrange() {
+  //
+  // Destructor
+  //
+  // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
+  // They will be deleted when fListCascade is deleted by the TSelector dtor
+  // Because of TList::SetOwner() ...
+  if (fCFContCascadeCuts && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fCFContCascadeCuts;     fCFContCascadeCuts = 0x0;  }
+  if (fESDtrackCuts)         { delete fESDtrackCuts;        fESDtrackCuts = 0x0; }
+}
+
+
+
+//________________________________________________________________________
+void AliQAProdMultistrange::UserCreateOutputObjects() {
+  // Create histograms
+  // Called once
+
+ //-----------------------------------------------
+ // Particle Identification Setup (new PID object)
+ //-----------------------------------------------
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+ fPIDResponse = inputHandler->GetPIDResponse();
+
+
+ // Only used to get the number of primary reconstructed tracks
+ if (fAnalysisType == "ESD" && (! fESDtrackCuts )){
+   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
+   //Printf("CheckCascade - ESDtrackCuts set up to 2010 std ITS-TPC cuts...");
+ }
+
+
+ //---------------------------------------------------
+ // Define the container for the topological variables
+ //---------------------------------------------------
+  if(! fCFContCascadeCuts) {
+      // Container meant to store all the relevant distributions corresponding to the cut variables.
+      // NB: overflow/underflow of variables on which we want to cut later should be 0!!! 
+      const Int_t  lNbSteps      =  4 ;
+      const Int_t  lNbVariables  =  21 ;
+      //Array for the number of bins in each dimension :
+      Int_t lNbBinsPerVar[lNbVariables] = {0};
+      lNbBinsPerVar[0]  = 25;     //DcaCascDaughters             :  [0.0,2.4,3.0]       -> Rec.Cut = 2.0;
+      lNbBinsPerVar[1]  = 25;     //DcaBachToPrimVertex          :  [0.0,0.24,100.0]    -> Rec.Cut = 0.01; 
+      lNbBinsPerVar[2]  = 30;     //CascCosineOfPointingAngle    :  [0.97,1.0]          -> Rec.Cut = 0.98;
+      lNbBinsPerVar[3]  = 40;     //CascRadius                   :  [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
+      lNbBinsPerVar[4]  = 30;     //InvMassLambdaAsCascDghter    :  [1.1,1.3]           -> Rec.Cut = 0.008;
+      lNbBinsPerVar[5]  = 20;     //DcaV0Daughters               :  [0.0,2.0]           -> Rec.Cut = 1.5;
+      lNbBinsPerVar[6]  = 201;    //V0CosineOfPointingAngleToPV  :  [0.89,1.0]          -> Rec.Cut = 0.9;
+      lNbBinsPerVar[7]  = 40;     //V0Radius                     :  [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
+      lNbBinsPerVar[8]  = 40;     //DcaV0ToPrimVertex            :  [0.0,0.39,110.0]    -> Rec.Cut = 0.01;  
+      lNbBinsPerVar[9]  = 25;     //DcaPosToPrimVertex           :  [0.0,0.24,100.0]    -> Rec.Cut = 0.05;
+      lNbBinsPerVar[10] = 25;     //DcaNegToPrimVertex           :  [0.0,0.24,100.0]    -> Rec.Cut = 0.05
+      lNbBinsPerVar[11] = 150;    //InvMassXi                    :   2-MeV/c2 bins
+      lNbBinsPerVar[12] = 120;    //InvMassOmega                 :   2-MeV/c2 bins
+      lNbBinsPerVar[13] = 100;    //XiTransvMom                  :  [0.0,10.0]
+      lNbBinsPerVar[14] = 110;    //Y(Xi)                        :   0.02 in rapidity units
+      lNbBinsPerVar[15] = 110;    //Y(Omega)                     :   0.02 in rapidity units
+      lNbBinsPerVar[16] = 112;    //Proper lenght of cascade       
+      lNbBinsPerVar[17] = 112;    //Proper lenght of V0
+      lNbBinsPerVar[18] = 201;    //V0CosineOfPointingAngleToXiV
+      lNbBinsPerVar[19] = 11;     //Centrality
+      lNbBinsPerVar[20] = 100;    //ESD track multiplicity
+      //define the container
+      fCFContCascadeCuts = new AliCFContainer("fCFContCascadeCuts","Container for Cascade cuts", lNbSteps, lNbVariables, lNbBinsPerVar );
+      //Setting the bin limits 
+      //0 -  DcaXiDaughters
+      Double_t *lBinLim0  = new Double_t[ lNbBinsPerVar[0] + 1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[0]; i++) lBinLim0[i] = (Double_t)0.0 + (2.4 - 0.0)/(lNbBinsPerVar[0] - 1) * (Double_t)i;
+         lBinLim0[ lNbBinsPerVar[0] ] = 3.0;
+      fCFContCascadeCuts -> SetBinLimits(0, lBinLim0);
+      delete [] lBinLim0;
+      //1 - DcaToPrimVertexXi
+      Double_t *lBinLim1  = new Double_t[ lNbBinsPerVar[1] + 1 ];
+         for(Int_t i=0; i<lNbBinsPerVar[1]; i++) lBinLim1[i] = (Double_t)0.0 + (0.24  - 0.0)/(lNbBinsPerVar[1] - 1) * (Double_t)i;
+         lBinLim1[ lNbBinsPerVar[1] ] = 100.0;
+      fCFContCascadeCuts -> SetBinLimits(1, lBinLim1);
+      delete [] lBinLim1;
+      //2 - CascCosineOfPointingAngle 
+      fCFContCascadeCuts->SetBinLimits(2, 0.97, 1.);
+      //3 - CascRadius
+      Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[3]; i++)   lBinLim3[i]  = (Double_t)0.0   + (3.9  - 0.0 )/(lNbBinsPerVar[3] - 1)  * (Double_t)i ;
+         lBinLim3[ lNbBinsPerVar[3] ] = 1000.0;
+      fCFContCascadeCuts -> SetBinLimits(3,  lBinLim3 );
+      delete [] lBinLim3;
+      //4 - InvMassLambdaAsCascDghter
+      fCFContCascadeCuts->SetBinLimits(4, 1.1, 1.13);
+      //5 - DcaV0Daughters
+      fCFContCascadeCuts -> SetBinLimits(5, 0., 2.);
+      //6 - V0CosineOfPointingAngleToPV
+      fCFContCascadeCuts -> SetBinLimits(6, 0.8, 1.001);
+      //7 - V0Radius
+      Double_t *lBinLim7 = new Double_t[ lNbBinsPerVar[7] + 1];
+         for(Int_t i=0; i< lNbBinsPerVar[7];i++) lBinLim7[i] = (Double_t)0.0 + (3.9 - 0.0)/(lNbBinsPerVar[7] - 1) * (Double_t)i;
+         lBinLim7[ lNbBinsPerVar[7] ] = 1000.0;
+      fCFContCascadeCuts -> SetBinLimits(7, lBinLim7);
+      delete [] lBinLim7;
+      //8 - DcaV0ToPrimVertex
+      Double_t *lBinLim8  = new Double_t[ lNbBinsPerVar[8]+1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[8];i++)   lBinLim8[i]  = (Double_t)0.0   + (0.39  - 0.0 )/(lNbBinsPerVar[8]-1)  * (Double_t)i ;
+         lBinLim8[ lNbBinsPerVar[8]  ] = 100.0;
+      fCFContCascadeCuts -> SetBinLimits(8,  lBinLim8 );
+      delete [] lBinLim8;
+      //9 - DcaPosToPrimVertex
+      Double_t *lBinLim9  = new Double_t[ lNbBinsPerVar[9]+1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[9];i++)   lBinLim9[i]  = (Double_t)0.0   + (0.24  - 0.0 )/(lNbBinsPerVar[9]-1)  * (Double_t)i ;
+         lBinLim9[ lNbBinsPerVar[9]  ] = 100.0;
+      fCFContCascadeCuts -> SetBinLimits(9,  lBinLim9 );
+      delete [] lBinLim9;
+      //10 - DcaNegToPrimVertex
+      Double_t *lBinLim10  = new Double_t[ lNbBinsPerVar[10]+1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[10];i++)   lBinLim10[i]  = (Double_t)0.0   + (0.24  - 0.0 )/(lNbBinsPerVar[10]-1)  * (Double_t)i ;
+         lBinLim10[ lNbBinsPerVar[10]  ] = 100.0;
+      fCFContCascadeCuts -> SetBinLimits(10,  lBinLim10 );     
+      delete [] lBinLim10;
+      //11 - InvMassXi
+      fCFContCascadeCuts->SetBinLimits(11, 1.25, 1.40);
+      //12 - InvMassOmega
+      fCFContCascadeCuts->SetBinLimits(12, 1.62, 1.74);
+      //13 - XiTransvMom
+      fCFContCascadeCuts->SetBinLimits(13, 0.0, 10.0);
+      //14 - Y(Xi)
+      fCFContCascadeCuts->SetBinLimits(14, -1.1, 1.1);
+      //15 - Y(Omega)
+      fCFContCascadeCuts->SetBinLimits(15, -1.1, 1.1);
+      //16 - Proper time of cascade
+      Double_t *lBinLim16  = new Double_t[ lNbBinsPerVar[16]+1 ];
+         for(Int_t i=0; i< lNbBinsPerVar[16];i++) lBinLim16[i] = (Double_t) -1. + (110. + 1.0 ) / (lNbBinsPerVar[16] - 1) * (Double_t) i;
+         lBinLim16[ lNbBinsPerVar[16] ] = 2000.0;
+      fCFContCascadeCuts->SetBinLimits(16, lBinLim16);
+      //17 - Proper time of V0
+      fCFContCascadeCuts->SetBinLimits(17, lBinLim16);
+      //18 - V0CosineOfPointingAngleToXiV
+      fCFContCascadeCuts -> SetBinLimits(18, 0.8, 1.001);
+      //19
+      Double_t *lBinLim19  = new Double_t[ lNbBinsPerVar[19]+1 ];
+         for(Int_t i=3; i< lNbBinsPerVar[19]+1;i++)   lBinLim19[i]  = (Double_t)(i-1)*10.;
+         lBinLim19[0] = 0.0; 
+         lBinLim19[1] = 5.0; 
+         lBinLim19[2] = 10.0;
+      fCFContCascadeCuts->SetBinLimits(19,  lBinLim19 );     
+      delete [] lBinLim19;
+      //20
+      fCFContCascadeCuts->SetBinLimits(20, 0.0, 6000.0);
+      // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
+      fCFContCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates");
+      fCFContCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates");
+      fCFContCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates");
+      fCFContCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates");
+      // Setting the variable title, per axis
+      fCFContCascadeCuts->SetVarTitle(0,  "Dca(cascade daughters) (cm)");
+      fCFContCascadeCuts->SetVarTitle(1,  "ImpactParamToPV(bachelor) (cm)");
+      fCFContCascadeCuts->SetVarTitle(2,  "cos(cascade PA)");
+      fCFContCascadeCuts->SetVarTitle(3,  "R_{2d}(cascade decay) (cm)");
+      fCFContCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(as casc dghter) (GeV/c^{2})");
+      fCFContCascadeCuts->SetVarTitle(5,  "Dca(V0 daughters) in Xi (cm)");
+      fCFContCascadeCuts->SetVarTitle(6,  "cos(V0 PA) in cascade to PV");
+      fCFContCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
+      fCFContCascadeCuts->SetVarTitle(8,  "ImpactParamToPV(V0) (cm)");
+      fCFContCascadeCuts->SetVarTitle(9,  "ImpactParamToPV(Pos) (cm)");
+      fCFContCascadeCuts->SetVarTitle(10, "ImpactParamToPV(Neg) (cm)");
+      fCFContCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
+      fCFContCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
+      fCFContCascadeCuts->SetVarTitle(13, "pt(cascade) (GeV/c)");
+      fCFContCascadeCuts->SetVarTitle(14, "Y(Xi)");
+      fCFContCascadeCuts->SetVarTitle(15, "Y(Omega)");
+      fCFContCascadeCuts->SetVarTitle(16, "mL/p (cascade) (cm)");
+      fCFContCascadeCuts->SetVarTitle(17, "mL/p (V0) (cm)");
+      fCFContCascadeCuts->SetVarTitle(18,  "cos(V0 PA) in cascade to XiV");
+      fCFContCascadeCuts->SetVarTitle(19, "Centrality");
+      fCFContCascadeCuts->SetVarTitle(20, "ESD track multiplicity");
+  }
+
+PostData(1, fCFContCascadeCuts);
+
+}// end UserCreateOutputObjects
+
+
+//________________________________________________________________________
+void AliQAProdMultistrange::UserExec(Option_t *) {
+
+  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  // Main loop (called for each event)
+  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  //-----------------------
+  //Define ESD/AOD handlers 
+  AliESDEvent *lESDevent = 0x0;
+  AliAODEvent *lAODevent = 0x0;
+
+  //---------------------
+  //Check the PIDresponse
+  if(!fPIDResponse) {
+       AliError("Cannot get pid response");
+       return;
+  }
+
+  //__________________________________________________
+  // After these lines we should have an ESD/AOD event
+
+  //---------------------------------------------------------
+  //Load the InputEvent and check, before any event selection  
+  //---------------------------------------------------------
+  Float_t  lPrimaryTrackMultiplicity = -1.;
+  AliCentrality* centrality;
+  if (fAnalysisType == "ESD") {
+      lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
+      if (!lESDevent) {
+          AliWarning("ERROR: lESDevent not available \n");
+          return;
+      }
+      if (fCollidingSystem == "PbPb") lPrimaryTrackMultiplicity = fESDtrackCuts->CountAcceptedTracks(lESDevent);
+      if (fCollidingSystem == "PbPb") centrality = lESDevent->GetCentrality();
+      
+  } else if (fAnalysisType == "AOD") {
+      lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
+      if (!lAODevent) {
+          AliWarning("ERROR: lAODevent not available \n");
+          return;
+      }
+      if (fCollidingSystem == "PbPb") {
+          lPrimaryTrackMultiplicity = 0;
+          Int_t    nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
+          for (Int_t itrack = 0; itrack < nTrackMultiplicity; itrack++) {
+               AliAODTrack* track = lAODevent->GetTrack(itrack);
+               if (track->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) lPrimaryTrackMultiplicity++; 
+          }
+      }
+      if (fCollidingSystem == "PbPb") centrality = lAODevent->GetCentrality();
+  } else {
+    Printf("Analysis type (ESD or AOD) not specified \n");
+    return;
+  }
+
+  //-----------------------------------------
+  // Centrality selection for PbPb collisions
+  //-----------------------------------------
+  Float_t lcentrality = 0.;
+  if (fCollidingSystem == "PbPb") { 
+       if (fkUseCleaning) lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
+       else {
+           lcentrality = centrality->GetCentralityPercentileUnchecked(fCentrEstimator.Data());
+           if (centrality->GetQuality()>1) {
+               PostData(1, fCFContCascadeCuts);
+               return;
+           }
+       }
+       if (lcentrality<fCentrLowLim||lcentrality>=fCentrUpLim) { 
+           PostData(1, fCFContCascadeCuts);
+           return;
+       }
+  } else if (fCollidingSystem == "pp") lcentrality = 0.;
+
+
+  //----------------------------------------
+  // SDD selection for pp@2.76TeV collisions
+  //----------------------------------------
+  if (fCollidingSystem == "pp") {
+      if (fAnalysisType == "ESD") {
+          if (fkSDDSelectionOn) {
+              TString trcl = lESDevent->GetFiredTriggerClasses();
+              if      (fwithSDD) { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+              else if (!fwithSDD){ if((trcl.Contains("ALLNOTRD")))  { PostData(1, fCFContCascadeCuts); return; } }
+          }
+      } else if (fAnalysisType == "AOD") {
+          if (fkSDDSelectionOn) {
+              TString trcl = lAODevent->GetFiredTriggerClasses();
+              if      (fwithSDD)  { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+              else if (!fwithSDD) { if((trcl.Contains("ALLNOTRD")))  { PostData(1, fCFContCascadeCuts); return; } }
+          }
+      }
+  }
+
+  //--------------------------------------------
+  // Physics selection for pp@2.76TeV collisions
+  //--------------------------------------------
+  // - moved to the runGrid for the PbPb collisions
+  if (fCollidingSystem == "pp") {
+      if (fAnalysisType == "ESD") {
+          UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+          Bool_t isSelected = 0;
+          isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
+          if(! isSelected){ PostData(1, fCFContCascadeCuts); return; }
+      } else if (fAnalysisType == "AOD") {
+          UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+          Bool_t isSelected = 0;
+          isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
+          if(! isSelected){ PostData(1, fCFContCascadeCuts); return; }
+      }
+  }
+
+  //------------------------------
+  // Well-established PV selection
+  //------------------------------
+  Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0}; 
+  Double_t lMagneticField        = -10.;
+  if (fAnalysisType == "ESD") {
+       const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();   
+       const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex(); 
+       if (fkQualityCutNoTPConlyPrimVtx) {
+           const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
+           if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtx->GetStatus() ){
+               AliWarning(" No SPD prim. vertex nor prim. Tracking vertex ... return !");
+               PostData(1, fCFContCascadeCuts);
+               return;
+           }
+       }
+       lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
+       lMagneticField = lESDevent->GetMagneticField( );
+  } else if (fAnalysisType == "AOD") {
+       const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex();
+       if (!lPrimaryBestAODVtx){
+           AliWarning("No prim. vertex in AOD... return!");
+           PostData(1, fCFContCascadeCuts);
+           return;
+       }
+       lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
+       lMagneticField = lAODevent->GetMagneticField();  
+  }
+
+  //------------------------------------------
+  // Pilup selection for pp@2.76TeV collisions
+  //------------------------------------------
+  if (fCollidingSystem == "pp") { 
+      if (fAnalysisType == "ESD") {
+          if (fkQualityCutPileup) { if(lESDevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
+      } else if (fAnalysisType == "AOD") {
+          if (fkQualityCutPileup) { if(lAODevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
+      }
+  }
+
+  //----------------------------
+  // Vertex Z position selection
+  //----------------------------
+  if (fkQualityCutZprimVtxPos) {
+      if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
+          PostData(1, fCFContCascadeCuts);
+          return;
+      }
+  }
+
+
+
+  //////////////////////////////
+  // CASCADE RECONSTRUCTION PART
+  //////////////////////////////
+
+  //%%%%%%%%%%%%%
+  // Cascade loop
+  Int_t ncascades = 0;
+  if      (fAnalysisType == "ESD") ncascades = lESDevent->GetNumberOfCascades();
+  else if (fAnalysisType == "AOD") ncascades = lAODevent->GetNumberOfCascades();
+
+  for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop (ESD or AOD)
+          
+    // -------------------------------------
+    // - Initialisation of the local variables that will be needed for ESD/AOD
+    // -- Container variables (1st round)
+    Double_t lDcaXiDaughters              = -1. ;                   //[Container]
+    Double_t lXiCosineOfPointingAngle     = -1. ;                   //[Container]
+    Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };             //Useful to define other variables: radius fid. vol., ctau, etc. for cascade
+    Double_t lXiRadius                    = -1000. ;                //[Container]
+    UShort_t lPosTPCClusters              = -1;                     //To check the quality of the tracks. For ESD only ...
+    UShort_t lNegTPCClusters              = -1;                     //To check the quality of the tracks. For ESD only ...
+    UShort_t lBachTPCClusters             = -1;                     //To check the quality of the tracks. For ESD only ...
+    Double_t lInvMassLambdaAsCascDghter   = 0.;                     //[Container]
+    Double_t lDcaV0DaughtersXi            = -1.;                    //[Container]
+    Double_t lDcaBachToPrimVertexXi       = -1.;                    //[Container]
+    Double_t lDcaV0ToPrimVertexXi         = -1.;                    //[Container]
+    Double_t lDcaPosToPrimVertexXi        = -1.;                    //[Container]
+    Double_t lDcaNegToPrimVertexXi        = -1.;                    //[Container]
+    Double_t lV0CosineOfPointingAngle     = -1.;                    //[Container]
+    Double_t lV0toXiCosineOfPointingAngle = -1.;                    //[Container] 
+    Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. };             //Useful to define other variables: radius fid. vol., ctau, etc. for VO 
+    Double_t lV0RadiusXi                  = -1000.0;                //[Container]
+    Double_t lV0quality                   = 0.;                     //  ??
+    Double_t lInvMassXiMinus              = 0.;                     //[Container]
+    Double_t lInvMassXiPlus               = 0.;                     //[Container]
+    Double_t lInvMassOmegaMinus           = 0.;                     //[Container]
+    Double_t lInvMassOmegaPlus            = 0.;                     //[Container]
+    // -- PID treatment
+    Bool_t   lIsBachelorKaonForTPC = kFALSE; 
+    Bool_t   lIsBachelorPionForTPC = kFALSE; 
+    Bool_t   lIsNegPionForTPC      = kFALSE; 
+    Bool_t   lIsPosPionForTPC      = kFALSE; 
+    Bool_t   lIsNegProtonForTPC    = kFALSE; 
+    Bool_t   lIsPosProtonForTPC    = kFALSE; 
+    // -- More container variables and quality checks
+    Double_t lXiMomX           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
+    Double_t lXiMomY           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
+    Double_t lXiMomZ           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
+    Double_t lXiTransvMom      = 0.;                               //[Container]
+    Double_t lXiTotMom         = 0.;                               //Useful to define other variables: cTau
+    Double_t lV0PMomX          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+    Double_t lV0PMomY          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+    Double_t lV0PMomZ          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+    Double_t lV0NMomX          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+    Double_t lV0NMomY          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+    Double_t lV0NMomZ          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+    Double_t lV0TotMom         = 0.;                               //Useful to define other variables: lctauV0
+    Double_t lBachMomX         = 0.;                               //Useful to define other variables: lBachTransvMom
+    Double_t lBachMomY         = 0.;                               //Useful to define other variables: lBachTransvMom
+    Double_t lBachMomZ         = 0.;                               //Useful to define other variables: lBachTransvMom
+    Double_t lBachTransvMom    = 0.;                               //Selection on the min bachelor pT
+    Double_t lpTrackTransvMom  = 0.;                               //Selection on the min bachelor pT
+    Double_t lnTrackTransvMom  = 0.;                               //Selection on the min bachelor pT
+    Short_t  lChargeXi         = -2;                               //Useful to select the particles based on the charge
+    Double_t lRapXi            = -20.0;                            //[Container]
+    Double_t lRapOmega         = -20.0;                            //[Container]
+    Float_t  etaBach           = 0.;                               //Selection on the eta range
+    Float_t  etaPos            = 0.;                               //Selection on the eta range
+    Float_t  etaNeg            = 0.;                               //Selection on the eta range
+    // --  variables for the AliCFContainer dedicated to cascade cut optmisiation: ESD and AOD 
+    if (fAnalysisType == "ESD") { 
+  
+          // -------------------------------------
+          // - Load the cascades from the handler 
+          AliESDcascade *xi = lESDevent->GetCascade(iXi);
+          if (!xi) continue;
+        
+          // ---------------------------------------------------------------------------
+          // - Assigning the necessary variables for specific AliESDcascade data members       
+          lV0quality = 0.;
+          xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
+          lDcaXiDaughters         = xi->GetDcaXiDaughters();
+          lXiCosineOfPointingAngle = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
+                                      // Take care : the best available vertex should be used (like in AliCascadeVertexer)
+          xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] ); 
+          lXiRadius            = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
+               
+          // -------------------------------------------------------------------------------------------------------------------------------
+          // - Around the tracks : Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
+          UInt_t lIdxPosXi     = (UInt_t) TMath::Abs( xi->GetPindex() );
+          UInt_t lIdxNegXi     = (UInt_t) TMath::Abs( xi->GetNindex() );
+          UInt_t lBachIdx      = (UInt_t) TMath::Abs( xi->GetBindex() );
+                                    // Care track label can be negative in MC production (linked with the track quality)
+                                    // However = normally, not the case for track index ...
+          // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
+          if (lBachIdx == lIdxNegXi) continue;    
+          if (lBachIdx == lIdxPosXi) continue; 
+          // - Get the track for the daughters
+          AliESDtrack *pTrackXi                = lESDevent->GetTrack( lIdxPosXi );
+          AliESDtrack *nTrackXi                = lESDevent->GetTrack( lIdxNegXi );
+          AliESDtrack *bachTrackXi     = lESDevent->GetTrack( lBachIdx );
+          if (!pTrackXi || !nTrackXi || !bachTrackXi )  continue;
+          // - Get the TPCnumber of cluster for the daughters
+          lPosTPCClusters   = pTrackXi->GetTPCNcls();
+          lNegTPCClusters   = nTrackXi->GetTPCNcls();
+          lBachTPCClusters  = bachTrackXi->GetTPCNcls();
+      
+          // ------------------------------------
+          // - Rejection of a poor quality tracks
+          if (fkQualityCutTPCrefit) {
+                // 1 - Poor quality related to TPCrefit
+                ULong_t pStatus    = pTrackXi->GetStatus();
+                ULong_t nStatus    = nTrackXi->GetStatus();
+                ULong_t bachStatus = bachTrackXi->GetStatus();
+                if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
+                if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
+                if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning(" Bach.   track has no TPCrefit ... continue!"); continue; }
+          }
+          if (fkQualityCutnTPCcls) {
+                // 2 - Poor quality related to TPC clusters
+                if (lPosTPCClusters  < fMinnTPCcls) { AliWarning(" V0 Pos. track has less than minn TPC clusters ... continue!"); continue; }
+                if (lNegTPCClusters  < fMinnTPCcls) { AliWarning(" V0 Neg. track has less than minn TPC clusters ... continue!"); continue; }
+                if (lBachTPCClusters < fMinnTPCcls) { AliWarning(" Bach.   track has less than minn TPC clusters ... continue!"); continue; }
+          }
+
+          // ------------------------------
+          etaPos  = pTrackXi->Eta();             
+          etaNeg  = nTrackXi->Eta();
+          etaBach = bachTrackXi->Eta();
+          lInvMassLambdaAsCascDghter = xi->GetEffMass();
+          lDcaV0DaughtersXi          = xi->GetDcaV0Daughters(); 
+          lV0CosineOfPointingAngle   = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
+          lDcaV0ToPrimVertexXi       = xi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );        
+          lDcaBachToPrimVertexXi     = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1],    lMagneticField) ); 
+          xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] ); 
+          lV0RadiusXi               = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
+          lDcaPosToPrimVertexXi      = TMath::Abs( pTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) ); 
+          lDcaNegToPrimVertexXi      = TMath::Abs( nTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) ); 
+       
+           //----------------------------------------------------------------------------------------------------       
+           // - Around effective masses. Change mass hypotheses to cover all the possibilities:  Xi-/+, Omega -/+
+           if (bachTrackXi->Charge() < 0) {
+                // Calculate the effective mass of the Xi- candidate. pdg code 3312 = Xi-
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , 3312);   
+                lInvMassXiMinus = xi->GetEffMassXi();
+                // Calculate the effective mass of the Xi- candidate. pdg code 3334 = Omega-
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , 3334);   
+                lInvMassOmegaMinus = xi->GetEffMassXi();
+               // Back to default hyp.                 
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , 3312);
+           }// end if negative bachelor
+           if ( bachTrackXi->Charge() >  0 ) {
+                // Calculate the effective mass of the Xi+ candidate. pdg code -3312 = Xi+
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , -3312);  
+                lInvMassXiPlus = xi->GetEffMassXi();
+                // Calculate the effective mass of the Xi+ candidate. pdg code -3334  = Omega+
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , -3334);  
+                lInvMassOmegaPlus = xi->GetEffMassXi();
+               // Back to "default" hyp.
+                lV0quality = 0.;
+                xi->ChangeMassHypothesis(lV0quality , -3312); 
+           }// end if positive bachelor
+
+           // ----------------------------------------------   
+          // - TPC PID : 3-sigma bands on Bethe-Bloch curve
+           // Bachelor
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
+           // Negative V0 daughter
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
+           // Positive V0 daughter
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
+        
+           // ------------------------------
+           // - Miscellaneous pieces of info that may help regarding data quality assessment.
+           xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
+           lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
+           lXiTotMom   = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
+           xi->GetNPxPyPz(lV0NMomX,lV0NMomY,lV0NMomZ);
+           xi->GetPPxPyPz(lV0PMomX,lV0PMomY,lV0PMomZ);
+           lV0TotMom = TMath::Sqrt(TMath::Power(lV0NMomX+lV0PMomX,2)+TMath::Power(lV0NMomY+lV0PMomY,2)+TMath::Power(lV0NMomZ+lV0PMomZ,2));     
+           xi->GetBPxPyPz( lBachMomX, lBachMomY, lBachMomZ );
+           lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
+           lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX + lV0NMomY*lV0NMomY );
+           lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX + lV0PMomY*lV0PMomY );
+           lChargeXi = xi->Charge();
+           lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
+           lRapXi    = xi->RapXi();
+           lRapOmega = xi->RapOmega();
+       
+    } else if (fAnalysisType == "AOD") {
+
+           // -------------------------------------
+           // - Load the cascades from the handler     
+           const AliAODcascade *xi = lAODevent->GetCascade(iXi);
+           if (!xi) continue;
+               
+          //----------------------------------------------------------------------------        
+           // - Assigning the necessary variables for specific AliESDcascade data members  
+           lDcaXiDaughters             = xi->DcaXiDaughters();
+           lXiCosineOfPointingAngle    = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0], 
+                                                         lBestPrimaryVtxPos[1], 
+                                                         lBestPrimaryVtxPos[2] );
+           lPosXi[0] = xi->DecayVertexXiX();
+           lPosXi[1] = xi->DecayVertexXiY();
+           lPosXi[2] = xi->DecayVertexXiZ();
+           lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );     
+
+           //-------------------------------------------------------------------------------------------------------------------------------
+           // - Around the tracks: Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
+           AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
+           AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
+           AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
+           if (!pTrackXi || !nTrackXi || !bachTrackXi ) continue;
+           UInt_t lIdxPosXi  = (UInt_t) TMath::Abs( pTrackXi->GetID() );  
+           UInt_t lIdxNegXi  = (UInt_t) TMath::Abs( nTrackXi->GetID() );
+           UInt_t lBachIdx   = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
+                                // Care track label can be negative in MC production (linked with the track quality)
+                                // However = normally, not the case for track index ...
+
+           // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
+           if (lBachIdx == lIdxNegXi) continue; 
+           if (lBachIdx == lIdxPosXi) continue;
+           // - Get the TPCnumber of cluster for the daughters
+           lPosTPCClusters   = pTrackXi->GetTPCNcls(); // FIXME: Is this ok? or something like in LambdaK0PbPb task AOD?
+           lNegTPCClusters   = nTrackXi->GetTPCNcls();
+           lBachTPCClusters  = bachTrackXi->GetTPCNcls();
+
+           // ------------------------------------
+           // - Rejection of a poor quality tracks
+           if (fkQualityCutTPCrefit) {
+                // - Poor quality related to TPCrefit
+                if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit)))    { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
+                if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit)))    { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
+                if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning(" Bach.   track has no TPCrefit ... continue!"); continue; }
+           }
+           if (fkQualityCutnTPCcls) {
+                // - Poor quality related to TPC clusters
+                if (lPosTPCClusters  < fMinnTPCcls) continue; 
+                if (lNegTPCClusters  < fMinnTPCcls) continue; 
+                if (lBachTPCClusters < fMinnTPCcls) continue; 
+           }
+  
+           // ------------------------------------------------------------------------------------------------------------------------------
+           // - Around the tracks: Bach + V0 (AOD). Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
+           etaPos  = pTrackXi->Eta();
+           etaNeg  = nTrackXi->Eta();
+           etaBach = bachTrackXi->Eta();
+           lChargeXi                   = xi->ChargeXi();
+           if ( lChargeXi < 0)         lInvMassLambdaAsCascDghter      = xi->MassLambda();
+           else                 lInvMassLambdaAsCascDghter     = xi->MassAntiLambda();
+           lDcaV0DaughtersXi           = xi->DcaV0Daughters(); 
+           lDcaV0ToPrimVertexXi                = xi->DcaV0ToPrimVertex();
+           lDcaBachToPrimVertexXi              = xi->DcaBachToPrimVertex(); 
+           lPosV0Xi[0] = xi->DecayVertexV0X();
+           lPosV0Xi[1] = xi->DecayVertexV0Y();
+           lPosV0Xi[2] = xi->DecayVertexV0Z(); 
+           lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
+           lV0CosineOfPointingAngle     = xi->CosPointingAngle( lBestPrimaryVtxPos ); 
+           lDcaPosToPrimVertexXi       = xi->DcaPosToPrimVertex(); 
+           lDcaNegToPrimVertexXi       = xi->DcaNegToPrimVertex(); 
+
+           // ---------------------------------------------------------------------------------------------------       
+           // - Around effective masses. Change mass hypotheses to cover all the possibilities:  Xi-/+, Omega -/+
+           if ( lChargeXi < 0 )                lInvMassXiMinus         = xi->MassXi();
+           if ( lChargeXi > 0 )                lInvMassXiPlus          = xi->MassXi();
+           if ( lChargeXi < 0 )                lInvMassOmegaMinus      = xi->MassOmega();
+           if ( lChargeXi > 0 )                lInvMassOmegaPlus       = xi->MassOmega();
+
+           // ----------------------------------------------
+           // - TPC PID : 3-sigma bands on Bethe-Bloch curve
+           // Bachelor
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
+           // Negative V0 daughter
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
+           // Positive V0 daughter
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
+
+           //---------------------------------
+           // - Extra info for QA (AOD)
+           // Miscellaneous pieces of info that may help regarding data quality assessment.
+           // Cascade transverse and total momentum     
+           lXiMomX = xi->MomXiX();
+           lXiMomY = xi->MomXiY();
+           lXiMomZ = xi->MomXiZ();
+           lXiTransvMom        = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
+           lXiTotMom   = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
+           Double_t lV0MomX = xi->MomV0X();
+           Double_t lV0MomY = xi->MomV0Y();
+           Double_t lV0MomZ = xi->MomV0Z();
+           lV0TotMom = TMath::Sqrt(TMath::Power(lV0MomX,2)+TMath::Power(lV0MomY,2)+TMath::Power(lV0MomZ,2));
+           lBachMomX = xi->MomBachX();
+           lBachMomY = xi->MomBachY();
+           lBachMomZ = xi->MomBachZ();         
+           lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
+           lV0NMomX = xi->MomNegX();
+           lV0NMomY = xi->MomNegY();
+           lV0PMomX = xi->MomPosX();
+           lV0PMomY = xi->MomPosY(); 
+           lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX   + lV0NMomY*lV0NMomY );
+           lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX   + lV0PMomY*lV0PMomY );
+           lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
+           lRapXi    = xi->RapXi();
+           lRapOmega = xi->RapOmega();
+
+    }// end of AOD treatment
+
+    //---------------------------------------
+    // Cut on pt of the three daughter tracks
+    if (lBachTransvMom<fMinPtCutOnDaughterTracks) continue;
+    if (lpTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
+    if (lnTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
+    //---------------------------------------------------
+    // Cut on pseudorapidity of the three daughter tracks
+    if (TMath::Abs(etaBach)>fEtaCutOnDaughterTracks) continue;
+    if (TMath::Abs(etaPos)>fEtaCutOnDaughterTracks) continue;
+    if (TMath::Abs(etaNeg)>fEtaCutOnDaughterTracks) continue;
+
+    //----------------------------------
+    // Calculate proper time for cascade
+    Double_t cascadeMass = 0.;
+    if ( ( (lChargeXi<0) && lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) ||
+         ( (lChargeXi>0) && lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC )  ) cascadeMass = 1.321;
+    if ( ( (lChargeXi<0) && lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC ) ||
+         ( (lChargeXi>0) && lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )  ) cascadeMass = 1.672; 
+    Double_t lctau =  TMath::Sqrt(TMath::Power((lPosXi[0]-lBestPrimaryVtxPos[0]),2)+TMath::Power((lPosXi[1]-lBestPrimaryVtxPos[1]),2)+TMath::Power(( lPosXi[2]-lBestPrimaryVtxPos[2]),2));
+    if (lXiTotMom!=0)         lctau = lctau*cascadeMass/lXiTotMom;
+    else lctau = -1.;
+    // Calculate proper time for Lambda (reconstructed)
+    Float_t lambdaMass = 1.115683; // PDG mass
+    Float_t distV0Xi =  TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2)+TMath::Power((lPosV0Xi[2]-lPosXi[2]),2));
+    Float_t lctauV0 = -1.;
+    if (lV0TotMom!=0) lctauV0 = distV0Xi*lambdaMass/lV0TotMom;
+
+               
+    // Fill the AliCFContainer (optimisation of topological selections)
+    Double_t lContainerCutVars[21] = {0.0};
+    lContainerCutVars[0]  = lDcaXiDaughters;
+    lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
+    lContainerCutVars[2]  = lXiCosineOfPointingAngle;
+    lContainerCutVars[3]  = lXiRadius;
+    lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
+    lContainerCutVars[5]  = lDcaV0DaughtersXi;
+    lContainerCutVars[6]  = lV0CosineOfPointingAngle;
+    lContainerCutVars[7]  = lV0RadiusXi;
+    lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
+    lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
+    lContainerCutVars[10] = lDcaNegToPrimVertexXi;
+    lContainerCutVars[13] = lXiTransvMom;
+    lContainerCutVars[16] = lctau;
+    lContainerCutVars[17] = lctauV0;
+    lContainerCutVars[18] = lV0toXiCosineOfPointingAngle;
+    lContainerCutVars[19] = lcentrality;
+    lContainerCutVars[20] = lPrimaryTrackMultiplicity;
+    if ( lChargeXi < 0 ) {
+         lContainerCutVars[11] = lInvMassXiMinus;
+         lContainerCutVars[12] = lInvMassOmegaMinus;
+         lContainerCutVars[14] = lRapXi;
+         lContainerCutVars[15] = -1.;
+         if (lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
+         lContainerCutVars[11] = lInvMassXiMinus;
+         lContainerCutVars[12] = lInvMassOmegaMinus;
+         lContainerCutVars[14] = -1.;
+         lContainerCutVars[15] = lRapOmega;
+         if (lIsBachelorKaonForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
+    } else {
+         lContainerCutVars[11] = lInvMassXiPlus;
+         lContainerCutVars[12] = lInvMassOmegaPlus;
+         lContainerCutVars[14] = lRapXi;
+         lContainerCutVars[15] = -1.;
+         if (lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
+         lContainerCutVars[11] = lInvMassXiPlus;
+         lContainerCutVars[12] = lInvMassOmegaPlus;
+         lContainerCutVars[14] = -1.;
+         lContainerCutVars[15] = lRapOmega;
+         if (lIsBachelorKaonForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,3); // for Omega+ 
+    }
+
+
+  }// end of the Cascade loop (ESD or AOD)
+    
+  
+  // Post output data.
+  PostData(1, fCFContCascadeCuts); 
+}
+
+//________________________________________________________________________
+void AliQAProdMultistrange::Terminate(Option_t *) 
+{
+
+}
diff --git a/PWGLF/QATasks/AliQAProdMultistrange.h b/PWGLF/QATasks/AliQAProdMultistrange.h
new file mode 100644 (file)
index 0000000..d8b999b
--- /dev/null
@@ -0,0 +1,102 @@
+#ifndef ALIANALYSISTASKCHECKCASCADEPBPB_H
+#define ALIANALYSISTASKCHECKCASCADEPBPB_H
+
+/*  See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+//            AliQAProdMultistrange class
+//              Origin AliAnalysisTaskCheckCascade
+//              This task has four roles :
+//                1. QAing the Cascades from ESD and AOD
+//                   Origin:  AliAnalysisTaskESDCheckV0 by Boris Hippolyte Nov2007, hippolyt@in2p3.fr
+//                2. Prepare the plots which stand as raw material for yield extraction (wi/wo PID)
+//                3. Supply an AliCFContainer meant to define the optimised topological selections
+//                4. Rough azimuthal correlation study (Eta, Phi)
+//                Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
+//                Modified :           A.Maire Mar2010, antonin.maire@ires.in2p3.fr
+//                Modified for PbPb analysis: M. Nicassio Feb 2011, maria.nicassio@ba.infn.it
+//-----------------------------------------------------------------
+
+class TList;
+class TH1F;
+class TH2F;
+class TH3F;
+class TVector3;
+class THnSparse;
+class AliESDEvent;
+class AliPhysicsSelection;
+class AliCFContainer;
+class AliPIDResponse;
+
+#include "TString.h"
+
+#include "AliAnalysisTaskSE.h"
+
+class AliQAProdMultistrange : public AliAnalysisTaskSE {
+ public:
+  AliQAProdMultistrange();
+  AliQAProdMultistrange(const char *name);
+  virtual ~AliQAProdMultistrange();
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  
+  void SetAnalysisType               (const char* analysisType          = "ESD" ) { fAnalysisType                = analysisType;               }
+  void SetCollidingSystem            (const char* collidingSystem       = "PbPb") { fCollidingSystem             = collidingSystem;            }
+  void SetSDDselection               (Bool_t SDDSelection               = kFALSE) { fkSDDSelectionOn             = SDDSelection;               }
+  void SetQualityCutZprimVtxPos      (Bool_t qualityCutZprimVtxPos      = kTRUE ) { fkQualityCutZprimVtxPos      = qualityCutZprimVtxPos;      }
+  void SetQualityCutNoTPConlyPrimVtx (Bool_t qualityCutNoTPConlyPrimVtx = kTRUE ) { fkQualityCutNoTPConlyPrimVtx = qualityCutNoTPConlyPrimVtx; }
+  void SetQualityCutTPCrefit         (Bool_t qualityCutTPCrefit         = kTRUE ) { fkQualityCutTPCrefit         = qualityCutTPCrefit;         }
+  void SetQualityCutnTPCcls          (Bool_t qualityCutnTPCcls          = kTRUE ) { fkQualityCutnTPCcls          = qualityCutnTPCcls;          }
+  void SetQualityCutMinnTPCcls       (Int_t  minnTPCcls                 = 70    ) { fMinnTPCcls                  = minnTPCcls;                 }
+  void SetQualityCutPileup           (Bool_t qualitycutPileup           = kFALSE) { fkQualityCutPileup           = qualitycutPileup;           }
+  void SetwithSDD                    (Bool_t withSDD                    = kTRUE ) { fwithSDD                     = withSDD;                    } 
+  void SetCentralityLowLim           (Float_t centrlowlim               = 0.    ) { fCentrLowLim                 = centrlowlim;                }  
+  void SetCentralityUpLim            (Float_t centruplim                = 100.  ) { fCentrUpLim                  = centruplim;                 }
+  void SetCentralityEst              (TString   centrest                = "V0M" ) { fCentrEstimator              = centrest;                   }
+  void SetUseCleaning                (Bool_t   usecleaning              = kTRUE ) { fkUseCleaning                = usecleaning;                }
+  void SetVertexRange                (Float_t vtxrange                  = 0.    ) { fVtxRange                    = vtxrange;                   }
+  void SetMinptCutOnDaughterTracks   (Float_t minptdaughtrks            = 0.    ) { fMinPtCutOnDaughterTracks    = minptdaughtrks;             }
+  void SetEtaCutOnDaughterTracks     (Float_t etadaughtrks              = 0.    ) { fEtaCutOnDaughterTracks      = etadaughtrks;               }
+
+ private:
+        // Note : In ROOT, "//!" means "do not stream the data from Master node to Worker node" ...
+        // your data member object is created on the worker nodes and streaming is not needed.
+        // http://root.cern.ch/download/doc/11InputOutput.pdf, page 14
+
+
+        TString         fAnalysisType;                  // "ESD" or "AOD" analysis type        
+        AliESDtrackCuts *fESDtrackCuts;                 // ESD track cuts used for primary track definition
+        TString         fCollidingSystem;               // "PbPb" or "pp" colliding system
+        AliPIDResponse *fPIDResponse;                   //! PID response object
+        Bool_t          fkSDDSelectionOn;               // Boolean : kTRUE = apply the selection on SDD status
+        Bool_t          fkQualityCutZprimVtxPos;        // Boolean : kTRUE = cut on the prim.vtx  z-position
+        Bool_t          fkQualityCutNoTPConlyPrimVtx;   // Boolean : kTRUE = prim vtx should be SPD or Tracking vertex
+        Bool_t          fkQualityCutTPCrefit;           // Boolean : kTRUE = ask for TPCrefit for the 3 daughter tracks
+        Bool_t          fkQualityCutnTPCcls;            // Boolean : kTRUE = ask for at least n TPC clusters for each daughter track
+        Bool_t          fkQualityCutPileup;             // Boolean : kTRUE = ask for no pileup events
+        Bool_t          fwithSDD;                       // Boolean : kTRUE = Select the events that has and use the info from the SDD
+        Int_t           fMinnTPCcls;                    // minimum number of TPC cluster for daughter tracks
+        Float_t         fCentrLowLim;                   // Lower limit for centrality percentile selection
+        Float_t         fCentrUpLim;                    // Upper limit for centrality percentile selection
+        TString         fCentrEstimator;                // string for the centrality estimator
+        Bool_t          fkUseCleaning;                  // Boolean : kTRUE = uses all the cleaning criteria of centrality selections (vertex cut + outliers) otherwise only outliers
+        Float_t         fVtxRange;                      // to select events with |zvtx|<fVtxRange cm
+        Float_t         fMinPtCutOnDaughterTracks;      // minimum pt cut on daughter tracks
+        Float_t         fEtaCutOnDaughterTracks;        // pseudorapidity cut on daughter tracks
+       
+
+        
+       AliCFContainer  *fCFContCascadeCuts;            //! Container meant to store all the relevant distributions corresponding to the cut variables
+       
+       
+
+  AliQAProdMultistrange(const AliQAProdMultistrange&);            // not implemented
+  AliQAProdMultistrange& operator=(const AliQAProdMultistrange&); // not implemented
+  
+  ClassDef(AliQAProdMultistrange, 7);
+};
+
+#endif