Add missing classes (first try)
authorsma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 Aug 2010 13:07:39 +0000 (13:07 +0000)
committersma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 Aug 2010 13:07:39 +0000 (13:07 +0000)
PWG3/hfe/AliHFEV0cuts.cxx [new file with mode: 0644]
PWG3/hfe/AliHFEV0cuts.h [new file with mode: 0644]
PWG3/hfe/AliHFEV0info.cxx [new file with mode: 0644]
PWG3/hfe/AliHFEV0info.h [new file with mode: 0644]
PWG3/hfe/AliHFEspectrum.cxx [new file with mode: 0644]
PWG3/hfe/AliHFEspectrum.h [new file with mode: 0644]
PWG3/hfe/AliHFEtrdPIDqa.cxx [new file with mode: 0644]
PWG3/hfe/AliHFEtrdPIDqa.h [new file with mode: 0644]

diff --git a/PWG3/hfe/AliHFEV0cuts.cxx b/PWG3/hfe/AliHFEV0cuts.cxx
new file mode 100644 (file)
index 0000000..3001426
--- /dev/null
@@ -0,0 +1,1104 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//  * 20/04/2010 *
+// Class for optimising and applying V0 cuts to obtain clean V0 samples
+// Compatible with ESDs only
+//
+// Authors:
+//    Matus Kalisky <matus.kalisky@cern.ch>
+//
+
+#include "TDatabasePDG.h"
+
+#include "AliESDtrack.h"
+#include "AliMCEvent.h"
+#include "AliESDv0.h"
+#include "AliKFParticle.h"
+#include "AliKFVertex.h"
+
+#include "AliHFEcollection.h"
+
+#include "AliHFEV0cuts.h"
+
+ClassImp(AliHFEV0cuts)
+
+//________________________________________________________________
+AliHFEV0cuts::AliHFEV0cuts():
+  fQA(NULL)
+  , fMCEvent(NULL)
+  , fInputEvent(NULL)
+  , fPrimaryVertex(NULL)
+  , fIsMC(kFALSE)
+{
+  //
+  // Default constructor
+  //
+  
+
+}
+//________________________________________________________________
+AliHFEV0cuts::~AliHFEV0cuts()
+{
+  //
+  // destructor
+  //
+  if (fQA) delete fQA;
+}
+
+//________________________________________________________________
+AliHFEV0cuts::AliHFEV0cuts(const AliHFEV0cuts &ref):
+  TObject(ref)
+  , fQA(NULL)
+  , fMCEvent(NULL)
+  , fInputEvent(NULL)
+  , fPrimaryVertex(NULL)
+  , fIsMC(kFALSE)
+{
+  //
+  // Copy constructor
+  //
+  ref.Copy(*this);  
+}
+//________________________________________________________________
+AliHFEV0cuts &AliHFEV0cuts::operator=(const AliHFEV0cuts &ref){
+  //
+  // Assignment operator
+  //
+  if(this != &ref)
+    ref.Copy(*this);
+  return *this;  
+}
+//________________________________________________________________
+void AliHFEV0cuts::Copy(TObject &ref) const{
+  //
+  // Copy function
+  // 
+  AliHFEV0cuts &target = dynamic_cast<AliHFEV0cuts &>(ref);
+
+  if(fQA) target.fQA = dynamic_cast<AliHFEcollection *>(fQA->Clone());  
+
+  if(target.fMCEvent) delete target.fMCEvent;
+  target.fMCEvent = new AliMCEvent;
+  
+  if(target.fPrimaryVertex) delete target.fPrimaryVertex;
+  target.fPrimaryVertex = new AliKFVertex;
+
+  TObject::Copy(ref);
+  
+}
+//___________________________________________________________________
+void AliHFEV0cuts::Init(const char* name){
+  //
+  // initialize the output objects and create histograms
+  //
+
+  //
+  // all the "h_cut_XXX" histograms hare cut value distributions:
+  // [0] for all candidates
+  // [1] jus before the cut on given variable was applied, but after all the previous cuts
+  //
+
+  fQA = new AliHFEcollection("fQA", name);
+
+
+  // common for all V0s
+  fQA->CreateTH2Fvector1(2, "h_all_AP", "armenteros plot for all V0 candidates", 200, -1, 1, 200, 0, 0.25);
+
+  // gammas
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_CosPoint", "Gamma Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_DCA", "DCA between the gamma daughters; dca (cm); counts", 100, 0, 2);
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_VtxR", "Radius of the gamma conversion vertex; r (cm); counts", 1000, 0, 100);
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_OA", "opening angle of the gamma products; opening angle (rad); counts", 100, 0, 1);
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_PP", "gamma psi pair angle; psi pairangle (rad); counts", 100, 0, 2);
+  fQA->CreateTH1Fvector1(2, "h_cut_Gamma_Chi2", "gamma Chi2/NDF; Chi2/NDF; counts", 100, 0, 25);
+  fQA->CreateTH1Fvector1(9, "h_Gamma_Mass", "Invariant mass of gammas; mass (GeV/c^{2}); counts", 100, 0, 0.2);
+
+  // kaons
+  fQA->CreateTH1Fvector1(2, "h_cut_K0_CosPoint", "K0 Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
+  fQA->CreateTH1Fvector1(2, "h_cut_K0_DCA", "DCA between the K0 daughters; dca (cm); counts", 100, 0, 2);
+  fQA->CreateTH1Fvector1(2, "h_cut_K0_VtxR", "Radius of the K0 decay vertex; r (cm); counts", 1000, 0, 100);
+  fQA->CreateTH1Fvector1(2, "h_cut_K0_Chi2", "K0 Chi2/NDF; Chi2/NDF; counts", 100, 0, 25);
+  fQA->CreateTH2Fvector1(2, "h_cut_K0_AP", "Armenteros plot for K0 candidates; #alpha; q_{T} (GeV/c)", 100, -1, 1, 100, 0, 0.3);
+  fQA->CreateTH1Fvector1(8, "h_K0_Mass", "Invariant mass of K0; mass (GeV/c^{2}); counts", 125, 0.45, 0.55);
+
+  // lambda
+  fQA->CreateTH1Fvector1(2, "h_cut_L_CosPoint", "L Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
+  fQA->CreateTH1Fvector1(2, "h_cut_L_DCA", "DCA between the L daughters; dca (cm); counts", 100, 0, 2);
+  fQA->CreateTH1Fvector1(2, "h_cut_L_VtxR", "Radius of the L decay vertex; r (cm); counts", 1000, 0, 100);
+  fQA->CreateTH1Fvector1(2, "h_cut_L_Chi2", "L Chi2/NDF; Chi2/NDF; counts", 100, 0, 25);
+  fQA->CreateTH2Fvector1(2, "h_cut_L_AP", "Armenteros plot for Lambda candidates; #alpha; q_{T} (GeV/c)", 100, -1, 1, 100, 0, 0.3);
+  fQA->CreateTH2Fvector1(2, "h_cut_AL_AP", "Armenteros plot for anti Lambda candidates; #alpha; q_{T} (GeV/c)", 100, -1, 1, 100, 0, 0.3);
+  fQA->CreateTH2Fvector1(2, "h_cut_Gamma_AP", "Armenteros plot for gamma candidates; #alpha; q_{T} (GeV/c)", 100, -1, 1, 100, 0, 0.3);
+  fQA->CreateTH1Fvector1(9, "h_L_Mass", "Invariant mass of L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
+  fQA->CreateTH1Fvector1(9, "h_AL_Mass", "Invariant mass of anti L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
+
+  fQA->CreateTH2F("h_L_checks", "Lambda candidate check[0] -v- check[1]; check[0]; check[1]", 5, -0.75, 1.75, 6, -0.75, 1.75 );
+  
+  // electrons
+  fQA->CreateTH1Fvector1(9, "h_Electron_P", "Momenta of conversion electrons -cuts-; P (GeV/c); counts", 50, 0.1, 20, 0);
+
+  // K0 pions
+  fQA->CreateTH1Fvector1(8, "h_PionK0_P", "Momenta of K0 pions -cuts-; P (GeV/c) counts;", 50, 0.1, 20, 0);
+  
+  // L pions
+  fQA->CreateTH1Fvector1(9, "h_PionL_P", "Momenta of L pions -cuts-; P (GeV/c) counts;", 50, 0.1, 50, 0);
+  
+  // L protons
+  fQA->CreateTH1Fvector1(9, "h_ProtonL_P", "Momenta of L protons -cuts-; P (GeV/c) counts;", 50, 0.1, 20, 0);    
+
+  // single track cuts
+  fQA->CreateTH1F("h_ST_NclsTPC", "Number of TPC clusters", 161, -1, 160);
+  fQA->CreateTH1F("h_ST_TPCrefit", "TPC refit", 2, -0.5, 1.5);
+  fQA->CreateTH1F("h_ST_chi2TPCcls", "chi2 per TPC cluster", 100, 0, 10);
+  fQA->CreateTH1F("h_ST_TPCclsR", "TPC cluster ratio", 120, -0.1, 1.1);
+  fQA->CreateTH1F("h_ST_kinks", "kinks", 2, -0.5, 1.5);
+  fQA->CreateTH1F("h_ST_pt", "track pt", 100, 0.1, 20, 0);
+  fQA->CreateTH1F("h_ST_eta", "track eta", 100, -1.5, 1.5);
+
+  //
+  // possibly new cuts
+  //
+  
+  // Gamma
+  fQA->CreateTH2Fvector1(2, "h_cut_Gamma_OAvP", "open. ang. of the Gamma daughters versus Gamma mom; Gamma p (GeV/c); opening angle (pions) (rad)", 100, 0.1, 10, 200, 0., 0.2);
+  // K0
+  fQA->CreateTH2Fvector1(2, "h_cut_K0_OAvP", "open. ang. of the K0 daughters versus K0 momentum; K0 p (GeV/c); opening angle (pions) (rad)", 100, 0.1, 10, 100, 0, 3.5);
+  // Lambda
+  fQA->CreateTH2Fvector1(2, "h_cut_L_OAvP", "open. ang. of the L daughters versus L momentum; Lambda p (GeV/c); openeing angle pion-proton (rad)", 100, 0.1, 10, 100, 0, 3.5);
+  fQA->CreateTH2Fvector1(2, "h_cut_L_rdp_v_mp", "relative L daughter mom -v- mother mom; L mom (GeV/c); relative daughter mom p2/p1", 100, 0.1, 10, 100, 0, 1);
+  fQA->CreateTH2Fvector1(2, "h_cut_L_qT_v_mp", "A-P q_{T} -v- Lambda momentum; mom (GeV/c); q_{T} GeV/c", 100, 0.1, 10, 50, 0., 0.12);
+
+
+  // THnSparse histograms
+  
+  // THnSparse for the K0 mass
+  // to be looked at after merging run by run
+  // axes: mass, pt, theta, phi
+  {
+    Int_t nBin[4] = {100, 10, 10, 18};
+    Double_t nMin[4] = {0.45, 0.1, 0., 0.};
+    Double_t nMax[4] = {0.55, 10., TMath::Pi(), 2*TMath::Pi()};
+    TString htitle = "K0 sparse; mass (GeV/c^{2}); p_{T} (GeV/c); theta (rad); phi(rad)";
+    fQA->CreateTHnSparse("hK0", htitle, 4, nBin, nMin, nMax);
+    fQA->BinLogAxis("hK0", 1);
+  }
+
+}
+//________________________________________________________________
+Bool_t AliHFEV0cuts::TrackCutsCommon(AliESDtrack* track){
+  //
+  // singe track cuts commom for all particle candidates
+  //
+
+  if(!track) return kFALSE;
+  
+  // status word
+  ULong_t status = track->GetStatus();
+
+
+  // No. of TPC clusters
+  fQA->Fill("h_ST_NclsTPC", track->GetTPCNcls());
+  if(track->GetTPCNcls() < 80) return kFALSE;
+
+  // TPC refit
+  if((status & AliESDtrack::kTPCrefit)){
+    fQA->Fill("h_ST_TPCrefit", 1);
+  }
+  if(!(status & AliESDtrack::kTPCrefit)){
+    fQA->Fill("h_ST_TPCrefit", 0);
+    return kFALSE;
+  }
+
+  // Chi2 per TPC cluster
+  Int_t nTPCclusters = track->GetTPCclusters(0);
+  Float_t chi2perTPCcluster = track->GetTPCchi2()/Float_t(nTPCclusters);
+  fQA->Fill("h_ST_chi2TPCcls", chi2perTPCcluster);
+  if(chi2perTPCcluster > 3.5) return kFALSE;
+
+  // TPC cluster ratio
+  Float_t cRatioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t> (track->GetTPCNclsF()) : 1.;
+  fQA->Fill("h_ST_TPCclsR", cRatioTPC);
+  if(cRatioTPC < 0.6) return kFALSE;
+
+  // kinks
+  fQA->Fill("h_ST_kinks", track->GetKinkIndex(0));
+  if(track->GetKinkIndex(0) != 0) return kFALSE;
+
+  // pt
+  fQA->Fill("h_ST_pt",track->Pt());
+  if(track->Pt() < 0.1 || track->Pt() > 100) return kFALSE;
+
+  // eta
+  fQA->Fill("h_ST_eta", track->Eta());
+  //if(TMath::Abs(track->Eta()) > 0.9) return kFALSE;  
+
+  return kTRUE;
+}
+//________________________________________________________________
+Bool_t AliHFEV0cuts::V0CutsCommon(AliESDv0 *v0){
+  //
+  // V0 cuts common to all V0s
+  //
+  
+  AliESDtrack* dN, *dP; 
+  dP = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(v0->GetPindex()));
+  dN = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(v0->GetNindex())); 
+  
+  if(!dN || !dP) return kFALSE;
+
+  Int_t qP = dP->Charge();
+  Int_t qN = dN->Charge();
+
+  if((qP*qN) != -1) return kFALSE;
+
+  return kTRUE;
+}
+//________________________________________________________________
+Bool_t AliHFEV0cuts::GammaCuts(AliESDv0 *v0){
+  //
+  // gamma cuts 
+  //
+  
+  if(!v0) return kFALSE;
+
+  // loose cuts first
+  if(LooseRejectK0(v0) || LooseRejectLambda(v0)) return kFALSE;
+
+  AliVTrack* daughter[2];
+  Int_t pIndex = 0, nIndex = 0;
+  if(CheckSigns(v0)){
+    pIndex = v0->GetPindex();
+    nIndex = v0->GetNindex();
+  }
+  else{
+    pIndex = v0->GetNindex();
+    nIndex = v0->GetPindex();    
+  }
+  daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
+  daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
+  if(!daughter[0] || !daughter[1]) return kFALSE;
+
+  AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
+  if(!kfMother) return kFALSE;
+  
+  // production vertex is set in the 'CreateMotherParticle' function
+  kfMother->SetMassConstraint(0, 0.001);
+
+  AliESDtrack* d[2];
+  d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
+  d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
+
+  Float_t iMass = v0->GetEffMass(0, 0);
+  Float_t iP = v0->P();
+  Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
+
+  // Armenteros
+  Float_t ap[2];
+  Armenteros(v0, ap);
+
+  // Cut values
+  const Double_t cutCosPoint[2] = {0., 0.03};  // ORG [0., 0.03]
+  const Double_t cutDCA[2] = {0., 0.25};       // ORG [0., 0.25]
+  const Double_t cutProdVtxR[2] = {6., 50.};   // ORG [6., 9999]
+  const Double_t cutOAngle[2] = {0, 0.1};      // ORG [0., 0.1]
+  const Double_t cutPsiPair[2] = {0., 0.05};   // ORG [0. 0.05]
+  const Double_t cutMass = 0.05;               // ORG [0.05]
+  const Double_t cutChi2NDF = 7.;              // ORG [7.]
+  // armenteros cuts
+  const Double_t cutAlpha[2] = {0.35, 0.45};   // [0.35, 0.45]
+  const Double_t cutQT = 0.015;
+  
+  // Values
+   
+  // cos pointing angle
+  Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
+  cosPoint = TMath::ACos(cosPoint);
+
+  // DCA between daughters
+  Double_t dca = v0->GetDcaV0Daughters();
+
+  // Production vertex
+  Double_t x, y, z; 
+  v0->GetXYZ(x,y,z);
+  Double_t r = TMath::Sqrt(x*x + y*y);
+
+  // Opening angle
+  Double_t oAngle = OpenAngle(v0);
+
+  // psi pair 
+  Double_t psiPair = PsiPair(v0);
+  
+  // V0 chi2/ndf
+  Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
+
+  if(kfMother) delete kfMother; 
+  //
+  // Apply the cuts, produce QA plots (with mass cut)
+  //
+  fQA->Fill("h_Gamma_Mass", 0, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_all_AP", 0, ap[0], ap[1]);
+    fQA->Fill("h_Electron_P", 0, p[0]);
+    fQA->Fill("h_Electron_P", 0, p[1]);
+    fQA->Fill("h_cut_Gamma_CosPoint", 0, cosPoint);
+    fQA->Fill("h_cut_Gamma_CosPoint", 1, cosPoint);
+    fQA->Fill("h_cut_Gamma_DCA", 0, dca);
+    fQA->Fill("h_cut_Gamma_VtxR", 0, r);
+    fQA->Fill("h_cut_Gamma_OA", 0, oAngle);
+    fQA->Fill("h_cut_Gamma_PP", 0, psiPair);
+    fQA->Fill("h_cut_Gamma_Chi2", 0, chi2ndf);
+    fQA->Fill("h_cut_Gamma_OAvP", 0, iP, oAngle);      
+    fQA->Fill("h_cut_Gamma_AP", 0, ap[0], ap[1]);
+ }
+
+  if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 1, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 1, p[0]);
+    fQA->Fill("h_Electron_P", 1, p[1]);
+    fQA->Fill("h_cut_Gamma_DCA", 1, dca);
+  }
+  
+  if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 2, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 2, p[0]);
+    fQA->Fill("h_Electron_P", 2, p[1]);
+    fQA->Fill("h_cut_Gamma_VtxR", 1, r);
+  }
+
+  if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 3, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 3, p[0]);
+    fQA->Fill("h_Electron_P", 3, p[1]);
+    fQA->Fill("h_cut_Gamma_OA", 1, oAngle);
+  }
+
+  if(oAngle < cutOAngle[0] || oAngle > cutOAngle[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 4, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 4, p[0]);
+    fQA->Fill("h_Electron_P", 4, p[1]);
+    fQA->Fill("h_cut_Gamma_PP", 1, psiPair);
+  }
+
+  if(psiPair < cutPsiPair[0] || psiPair > cutPsiPair[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 5, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 5, p[0]);
+    fQA->Fill("h_Electron_P", 5, p[1]);
+    fQA->Fill("h_cut_Gamma_Chi2", 1, chi2ndf);
+  }
+
+  if(chi2ndf > cutChi2NDF) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 6, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 6, p[0]);
+    fQA->Fill("h_Electron_P", 6, p[1]);
+    fQA->Fill("h_cut_Gamma_OAvP", 1, iP, oAngle);      
+    fQA->Fill("h_all_AP", 1, ap[0], ap[1]);
+    fQA->Fill("h_cut_Gamma_AP", 1, ap[0], ap[1]);
+  }
+
+  if(TMath::Abs(ap[0]) > cutAlpha[0] && TMath::Abs(ap[0]) < cutAlpha[1]) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 7, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 7, p[0]);
+    fQA->Fill("h_Electron_P", 7, p[1]);
+  }
+
+  if(ap[1] > cutQT) return kFALSE;
+  fQA->Fill("h_Gamma_Mass", 8, iMass);
+  if(iMass < cutMass){
+    fQA->Fill("h_Electron_P", 8, p[0]);
+    fQA->Fill("h_Electron_P", 8, p[1]);
+  }
+  
+
+  if(iMass > cutMass) return kFALSE;
+
+  // all cuts passed
+    
+  return kTRUE;
+}
+//________________________________________________________________
+Bool_t AliHFEV0cuts::K0Cuts(AliESDv0 *v0){
+  //
+  // K0 cuts
+  //
+
+  if(!v0) return kFALSE;
+
+  const Double_t cK0mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();  // PDG K0s mass
+  AliVTrack* daughter[2];
+  Int_t pIndex = 0, nIndex = 0;
+  if(CheckSigns(v0)){
+    pIndex = v0->GetPindex();
+    nIndex = v0->GetNindex();
+  }
+  else{
+    pIndex = v0->GetNindex();
+    nIndex = v0->GetPindex();    
+  }
+  daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
+  daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
+  if(!daughter[0] || !daughter[1]) return kFALSE;
+
+  AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
+  if(!kfMother) return kFALSE;
+  // production vertex is set in the 'CreateMotherParticle' function
+  kfMother->SetMassConstraint(cK0mass, 0.);
+  
+  AliESDtrack* d[2];
+  d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
+  d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
+
+  Float_t iMass = v0->GetEffMass(2, 2);
+  Float_t iP = v0->P();
+  Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
+  Double_t theta = v0->Theta();
+  Double_t phi = v0->Phi();
+  Double_t pt = v0->Pt();
+  Double_t data[4] = {0., 0., 0., 0.};
+  // Cut values
+  const Double_t cutCosPoint[2] = {0., 0.03};  // ORG [0., 0.03]
+  const Double_t cutDCA[2] = {0., 0.2};        // ORG [0., 0.1]
+  const Double_t cutProdVtxR[2] = {2.0, 30.};   // ORG [0., 8.1]
+  const Double_t cutMass[2] = {0.49, 0.51};   // ORG [0.485, 0.51]
+  const Double_t cutChi2NDF = 5.;              // ORG [7.]
+  const Double_t cutOAngleP = (1.0/(iP + 0.3) - 0.1); // momentum dependent min. OAngle ~ 1/x
+  // cundidate cuts
+  // armenteros plot
+  const Double_t cutQT = 0.1075;
+  // elipse cut - see bellow
+
+  // Values
+
+  // cos pointing angle
+  Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
+  cosPoint = TMath::ACos(cosPoint);
+
+  // DCA between daughters
+  Double_t dca = v0->GetDcaV0Daughters();
+
+  // Production vertex
+  Double_t x, y, z; 
+  v0->GetXYZ(x,y,z);
+
+  Double_t r = TMath::Sqrt(x*x + y*y);  
+
+  // V0 chi2/ndf
+  Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
+  
+  if(kfMother) delete kfMother; 
+
+  // Opening angle
+  Double_t oAngle = OpenAngle(v0);
+  
+  // Armenteros
+  Float_t ap[2];
+  Armenteros(v0, ap);
+  const Double_t cutAP = 0.22 * TMath::Sqrt( TMath::Abs( (1-ap[0]*ap[0]/(0.92*0.92)) ) );
+  
+
+  //
+  // Apply the cuts, produce QA plots (with mass cut)
+  //
+
+  fQA->Fill("h_K0_Mass", 0, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_all_AP", 0, ap[0], ap[1]);
+    fQA->Fill("h_PionK0_P", 0, p[0]);
+    fQA->Fill("h_PionK0_P", 0, p[1]);
+    fQA->Fill("h_cut_K0_CosPoint", 0, cosPoint);
+    fQA->Fill("h_cut_K0_CosPoint", 1, cosPoint);
+    fQA->Fill("h_cut_K0_DCA", 0, dca);
+    fQA->Fill("h_cut_K0_VtxR", 0, r);
+    fQA->Fill("h_cut_K0_Chi2", 0, chi2ndf);
+    fQA->Fill("h_cut_K0_OAvP", 0, iP, oAngle);
+    fQA->Fill("h_cut_K0_AP", 0, ap[0], ap[1]);
+  }
+
+  if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
+  fQA->Fill("h_K0_Mass", 1, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 1, p[0]);
+    fQA->Fill("h_PionK0_P", 1, p[1]);
+    fQA->Fill("h_cut_K0_DCA", 1, dca);
+  }
+  
+  if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
+  fQA->Fill("h_K0_Mass", 2, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 2, p[0]);
+    fQA->Fill("h_PionK0_P", 2, p[1]);
+    fQA->Fill("h_cut_K0_VtxR", 1, r);
+  }
+
+  
+  if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
+  fQA->Fill("h_K0_Mass", 3, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 3, p[0]);
+    fQA->Fill("h_PionK0_P", 3, p[1]);
+    fQA->Fill("h_cut_K0_Chi2", 1, chi2ndf);
+  }
+
+  if(chi2ndf > cutChi2NDF) return kFALSE;
+  fQA->Fill("h_K0_Mass", 4, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 4, p[0]);
+    fQA->Fill("h_PionK0_P", 4, p[1]);
+    fQA->Fill("h_cut_K0_OAvP", 1, iP, oAngle);
+  }  
+
+  if(oAngle < cutOAngleP) return kFALSE;
+  fQA->Fill("h_K0_Mass", 5, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 5, p[0]);
+    fQA->Fill("h_PionK0_P", 5, p[1]);
+    fQA->Fill("h_cut_K0_AP", 1, ap[0], ap[1]);
+    fQA->Fill("h_all_AP", 1, ap[0], ap[1]);
+  }
+
+  if(ap[1] < cutQT) return kFALSE;
+  fQA->Fill("h_K0_Mass", 6, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 6, p[0]);
+    fQA->Fill("h_PionK0_P", 6, p[1]);
+  }
+    
+  if(ap[1] > cutAP) return kFALSE;
+  fQA->Fill("h_K0_Mass", 7, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_PionK0_P", 7, p[0]);
+    fQA->Fill("h_PionK0_P", 7, p[1]);
+  }
+
+  data[0] = iMass;
+  data[1] = pt;
+  data[2] = theta;
+  data[3] = phi;
+  //printf("-D: m: %f, pT: %f, theta: %f, phi: %f\n", invMass, mPt, theta, phi);
+  fQA->Fill("hK0", data);
+  
+
+  if(iMass < cutMass[0] || iMass > cutMass[1]) return kFALSE;
+
+  // all cuts passed
+  
+  return kTRUE;
+}
+//________________________________________________________________
+Bool_t AliHFEV0cuts::LambdaCuts(AliESDv0 *v0, Bool_t &isLambda ){
+  //
+  // Lambda cuts - decision on Lambda - AntiLambda is taken too
+  //
+  // discrimination between lambda and antilambda - correlation of the following variables necessary:
+  // - momentum of the proton AND momentum of the pion (proton momentum is allways larger)
+  // - mass of the mother particle
+
+  if(!v0) return kFALSE;
+
+  // loose cuts first
+  if(LooseRejectK0(v0) || LooseRejectGamma(v0)) return kFALSE;
+  
+  const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();  // PDG lambda mass
+
+  AliVTrack* daughter[2];
+  Int_t pIndex = 0, nIndex = 0;
+  Float_t mMass[2] = {-1., -1.};
+  if(CheckSigns(v0)){
+    pIndex = v0->GetPindex();
+    nIndex = v0->GetNindex();
+    mMass[0] = v0->GetEffMass(4, 2);
+    mMass[1] = v0->GetEffMass(2, 4);
+  }
+  else{
+    pIndex = v0->GetNindex();
+    nIndex = v0->GetPindex();    
+    mMass[0] = v0->GetEffMass(2, 4);
+    mMass[1] = v0->GetEffMass(4, 2);
+  }
+  daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
+  daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
+  if(!daughter[0] || !daughter[1]) return kFALSE;
+
+  AliKFParticle *kfMother[2] = {0x0, 0x0};
+  // Lambda
+  kfMother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
+  if(!kfMother[0]) return kFALSE;
+  
+  // production vertex is set in the 'CreateMotherParticle' function
+  kfMother[0]->SetMassConstraint(cL0mass, 0.);
+
+  // Anti Lambda
+  kfMother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
+  if(!kfMother[1]) return kFALSE;
+  // production vertex is set in the 'CreateMotherParticle' function
+  kfMother[1]->SetMassConstraint(cL0mass, 0.);
+
+  Float_t dMass[2] = {TMath::Abs(mMass[0] - cL0mass), TMath::Abs(mMass[1] - cL0mass)};
+  
+  AliESDtrack* d[2];
+  d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
+  d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
+  if(!d[0] || !d[1])    return kFALSE;
+  
+  Float_t p[2] = {d[0]->GetP(), d[1]->GetP()}; 
+
+  // check the 3 lambda - antilambda variables
+  Int_t check[2] = {-1, -1};   // 0 : lambda, 1 : antilambda
+  // 1) momentum of the daughter particles - proton is expected to have higher momentum than pion
+  check[0] = (p[0] > p[1]) ? 0 : 1;
+  // 2) mass of the mother particle
+  check[1] = (dMass[0] < dMass[1]) ? 0 : 1;
+  fQA->Fill("h_L_checks", check[0]*1.0, check[1]*1.0);
+  // if the two check do not agree
+  if(check[0] != check[1]){
+    if(kfMother[0]) delete kfMother[0]; 
+    if(kfMother[1]) delete kfMother[1]; 
+    return kFALSE;
+  }
+
+  // now that the check[0] == check[1]
+  const Int_t type = check[0];
+
+  Float_t iMass =0.;
+  if(CheckSigns(v0)){
+    iMass = (type == 0) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
+  }
+  else{
+    iMass = (type == 0) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
+  }
+  Float_t iP = v0->P();
+
+   // Cuts
+  const Double_t cutCosPoint[2] = {0., 0.03};  // ORG [0., 0.03]
+  const Double_t cutDCA[2] = {0., 0.2};        // ORG [0., 0.2]
+  const Double_t cutProdVtxR[2] = {2., 30.};   // ORG [0., 24.]
+  const Double_t cutMass[2] = {1.11, 1.12};   // ORG [1.11, 1.12]
+  const Double_t cutChi2NDF = 5.;              // ORG [5.]
+  // cundidate cuts
+  // opening angle as a function of L momentum
+  const Double_t cutOAngleP = 0.3 - 0.2*iP; // momentum dependent min. OAngle linear cut
+  // relative daughter momentum versusu mother momentum
+  // armenteros plot cuts
+  const Double_t cutQT = 0.03;
+  const Double_t cutAlpha = 0.7;  // VERY strong - should supress the overlap with K0
+  // next cut see below
+
+  // compute the cut values
+  
+  // cos pointing angle
+  Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
+  cosPoint = TMath::ACos(cosPoint);
+
+  // DCA between daughters
+  Double_t dca = v0->GetDcaV0Daughters();
+  
+  // Production vertex
+  Double_t x, y, z; 
+  v0->GetXYZ(x,y,z);
+  Double_t r = TMath::Sqrt(x*x + y*y);
+
+  Int_t ix[2] = {0, 1};
+  if(1 == type){
+    ix[0] = 1;
+    ix[1] = 0;
+  }
+  
+  // V0 chi2/ndf
+  Double_t chi2ndf = kfMother[type]->GetChi2()/kfMother[type]->GetNDF();
+
+  if(kfMother[0]) delete kfMother[0]; 
+  if(kfMother[1]) delete kfMother[1]; 
+
+  // Opening angle
+  Double_t oAngle = OpenAngle(v0);
+
+  // Relative daughter momentum
+  Double_t rP = (0 == check[0]) ? p[1]/p[0] : p[0]/p[1];
+  
+  // Armenteros
+  Float_t ap[2];
+  Armenteros(v0, ap);
+  
+  Double_t cutAP[2]; // a bit of woodoo :-)
+  cutAP[0] = 1.0 - (ap[0]-0.7 * ap[0]-0.7)*1.1 - 0.87;
+  cutAP[1] = 1.0 - (ap[0]+0.7 * ap[0]+0.7)*1.1 - 0.87;
+
+
+  //
+  // Apply the cuts, produce QA plots (with mass cut)
+  //
+  
+  (type == 0) ?   fQA->Fill("h_L_Mass", 0, iMass) :  fQA->Fill("h_AL_Mass", 0, iMass);
+
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_all_AP", 0, ap[0], ap[1]);
+    fQA->Fill("h_ProtonL_P", 0, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 0, p[ix[1]]);
+    fQA->Fill("h_cut_L_Chi2", 0, chi2ndf);
+    fQA->Fill("h_cut_L_CosPoint", 0, cosPoint);
+    fQA->Fill("h_cut_L_CosPoint", 1, cosPoint);
+    fQA->Fill("h_cut_L_DCA", 0, dca);
+    fQA->Fill("h_cut_L_VtxR", 0, r);
+    fQA->Fill("h_cut_L_OAvP", 0, iP, oAngle);
+    fQA->Fill("h_cut_L_rdp_v_mp", 0, iP, rP);
+    if(0 ==type)  fQA->Fill("h_cut_L_AP", 0, ap[0], ap[1]);
+    else fQA->Fill("h_cut_AL_AP", 0, ap[0], ap[1]);
+    fQA->Fill("h_cut_L_qT_v_mp", 0, iP, ap[1]);
+  }
+
+  if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 1, iMass) :  fQA->Fill("h_AL_Mass", 1, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 1, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 1, p[ix[1]]);
+    fQA->Fill("h_cut_L_DCA", 1, dca);
+  }
+  
+  if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 2, iMass) :  fQA->Fill("h_AL_Mass", 2, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 2, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 2, p[ix[1]]);
+     fQA->Fill("h_cut_L_VtxR", 1, r);
+  }
+
+  if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 3, iMass) :  fQA->Fill("h_AL_Mass", 3, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 3, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 3, p[ix[1]]);
+    fQA->Fill("h_cut_L_Chi2", 1, chi2ndf);
+  }
+
+
+  if(chi2ndf > cutChi2NDF) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 4, iMass) :  fQA->Fill("h_AL_Mass", 4, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 4, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 4, p[ix[1]]);
+    fQA->Fill("h_cut_L_OAvP", 1, iP, oAngle);
+  }
+
+  if(oAngle < cutOAngleP) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 5, iMass) :  fQA->Fill("h_AL_Mass", 5, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 5, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 5, p[ix[1]]);
+    fQA->Fill("h_cut_L_rdp_v_mp", 1, iP, rP);
+    if(0 == type)  fQA->Fill("h_cut_L_AP", 1, ap[0], ap[1]);
+    else fQA->Fill("h_cut_AL_AP", 1, ap[0], ap[1]);
+    fQA->Fill("h_cut_L_qT_v_mp", 1, iP, ap[1]);
+    fQA->Fill("h_all_AP", 1, ap[0], ap[1]);
+  }
+  
+  if(ap[1] < cutQT) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 6, iMass) :  fQA->Fill("h_AL_Mass", 6, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 6, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 6, p[ix[1]]);
+  }
+
+  if(ap[1] > cutAP[type]) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 7, iMass) :  fQA->Fill("h_AL_Mass", 7, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 7, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 7, p[ix[1]]);
+  }
+  if(TMath::Abs(ap[0]) > cutAlpha) return kFALSE;
+  (type == 0) ?   fQA->Fill("h_L_Mass", 8, iMass) :  fQA->Fill("h_AL_Mass", 8, iMass);
+  if(iMass > cutMass[0] && iMass < cutMass[1]){
+    fQA->Fill("h_ProtonL_P", 8, p[ix[0]]);
+    fQA->Fill("h_PionL_P", 8, p[ix[1]]);
+  }
+
+  if(iMass < cutMass[0] || iMass > cutMass[1]) {
+    return kFALSE;
+  }
+
+  // all cuts passed
+
+  // assign the lambda type value: Lambda: kTRUE, Anti-Lambda: kFALSE
+  isLambda = (0 == type) ? kTRUE : kFALSE;
+
+
+  return kTRUE;
+}
+//________________________________________________________________
+Double_t AliHFEV0cuts::OpenAngle(AliESDv0 *v0) const {
+  //
+  // Opening angle between two daughter tracks
+  //
+  Double_t mn[3] = {0,0,0};
+  Double_t mp[3] = {0,0,0};
+    
+  
+  v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+  v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
+
+  
+  Double_t openAngle = TMath::ACos((mp[0]*mn[0] + mp[1]*mn[1] + mp[2]*mn[2])/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] + mp[2]*mp[2])*TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] + mn[2]*mn[2])));
+  
+  return TMath::Abs(openAngle);
+}
+//________________________________________________________________
+Double_t AliHFEV0cuts::PsiPair(AliESDv0 *v0) {
+  //
+  // Angle between daughter momentum plane and plane 
+  // 
+
+  if(!fInputEvent) return -1.;
+
+  Float_t magField = fInputEvent->GetMagneticField();
+
+  Int_t pIndex = -1;
+  Int_t nIndex = -1;
+  if(CheckSigns(v0)){
+    pIndex = v0->GetPindex();
+    nIndex = v0->GetNindex();
+  }
+  else{
+    pIndex = v0->GetNindex();
+    nIndex = v0->GetPindex();    
+  }
+
+  AliESDtrack* daughter[2];
+
+  daughter[0] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(pIndex));
+  daughter[1] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(nIndex));
+
+  Double_t x, y, z;
+  v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
+  
+  Double_t mn[3] = {0,0,0};
+  Double_t mp[3] = {0,0,0};
+  
+
+  v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+  v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter; 
+
+
+  Double_t deltat = 1.;
+  deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) -  TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
+
+  Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
+
+  Double_t momPosProp[3];
+  Double_t momNegProp[3];
+    
+  AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
+    
+  Double_t psiPair = 4.;
+
+  if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
+    psiPair =  -5.;
+  if(pt.PropagateTo(radiussum,magField) == 0)
+    psiPair = -5.;
+  pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
+  nt.GetPxPyPz(momNegProp);
+  
+  Double_t pEle =
+    TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
+  Double_t pPos =
+    TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
+    
+  Double_t scalarproduct =
+    momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
+    
+  Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
+
+  psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));  
+
+  return psiPair; 
+}
+//________________________________________________________________
+AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec){
+  //
+  // Creates a mother particle
+  //
+  AliKFParticle pkfdaughter(*pdaughter, pspec);
+  AliKFParticle nkfdaughter(*ndaughter, nspec);
+  
+  // - check if the daughter particles are coming from the primary vertex
+  // - check the number of tracks that take part in the creaton of primary vertex.
+  //   important: after removeal of candidate tracks there must be at least 2 tracks left
+  //   otherwise the primary vertex will be corrupted  
+  
+  // ESD Analyis
+  //const AliESDVertex *esdvertex = dynamic_cast<const AliESDVertex *>(fInputEvent->GetPrimaryVertex());
+  //if(!esdvertex) return NULL;
+  //UShort_t *contrib = esdvertex->GetIndices();
+  
+  //
+  // not using the removal of the daughter track now
+  //
+  //   Int_t nTracks = esdvertex->GetNIndices();
+  //   printf(" -D: N Vertex tracks: %i\n", nTracks);
+  //   printf(" -D: N Contributors: %i\n", fPrimaryVertex->GetNContributors()); 
+  //   Int_t nfound = 0;
+  //   for(Int_t id = 0; id < esdvertex->GetNIndices(); id++){
+  //     if(contrib[id] == pdaughter->GetID()){
+  //       if( (nTracks - nfound) <= 2 ) return NULL;
+  //       *fPrimaryVertex -= pkfdaughter;
+  //       removed[0] = kTRUE;
+  //       nfound++;
+  //     }
+  //     if(contrib[id] == ndaughter->GetID()){
+  //       if( (nTracks - nfound) <=2 ) return NULL;
+  //       *fPrimaryVertex -= nkfdaughter;
+  //       removed[1] = kTRUE;
+  //       nfound++;
+  //     }
+  //     if(nfound == 2) break;
+  //   }
+  
+  //  printf(" -D: n removed: %i\n", nfound);
+  
+  // Create the mother particle 
+  AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);
+
+  AliKFVertex improvedVertex = *fPrimaryVertex;
+  improvedVertex += *m;
+  m->SetProductionVertex(improvedVertex);
+  
+  // update 15/06/2010
+  // mother particle will not be added to primary vertex but only to its copy 
+  // as this confilcts with calling
+  // m->SetPrimaryVertex() function and
+  // subsequently removing the mother particle afterwards
+  // Sourse: Sergey Gorbunov
+
+  return m;
+}
+//_________________________________________________
+Bool_t AliHFEV0cuts::LooseRejectK0(AliESDv0 * const v0) const {
+  //
+  // Reject K0 based on loose cuts
+  //
+  Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
+  if(mass > 0.494 && mass < 0.501) return kTRUE;
+  return kFALSE;
+}
+
+//_________________________________________________
+Bool_t AliHFEV0cuts::LooseRejectLambda(AliESDv0 * const v0) const {
+  //
+  // Reject Lambda based on loose cuts
+  //
+  Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
+  Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
+  
+  if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
+  if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
+  return kFALSE;
+}
+
+//_________________________________________________
+Bool_t AliHFEV0cuts::LooseRejectGamma(AliESDv0 * const v0) const {
+  //
+  // Reject Lambda based on loose cuts
+  //
+  Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
+  
+  if(mass < 0.02) return kTRUE;
+  return kFALSE;
+}
+//___________________________________________________________________
+void  AliHFEV0cuts::Armenteros(AliESDv0 *v0, Float_t val[2]){
+  //
+  // computes the Armenteros variables for given V0
+  // fills the histogram
+  // returns the values via "val"
+  //
+  
+  Double_t mn[3] = {0,0,0};
+  Double_t mp[3] = {0,0,0};  
+  Double_t mm[3] = {0,0,0};  
+
+  if(CheckSigns(v0)){
+    v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
+    v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
+  }
+  else{
+    v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
+    v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
+  }
+  v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
+
+  TVector3 vecN(mn[0],mn[1],mn[2]);
+  TVector3 vecP(mp[0],mp[1],mp[2]);
+  TVector3 vecM(mm[0],mm[1],mm[2]);
+  
+  Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
+  Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
+  
+  Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
+    ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
+  Double_t qt = vecP.Mag()*sin(thetaP);
+
+  val[0] = alfa;
+  val[1] = qt;
+
+}
+//___________________________________________________________________
+Bool_t AliHFEV0cuts::CheckSigns(AliESDv0* const v0){
+  //
+  // check wheter the sign was correctly applied to 
+  // V0 daughter tracks
+  //
+  
+  Bool_t correct = kFALSE;
+
+  Int_t pIndex = 0, nIndex = 0;
+  pIndex = v0->GetPindex();
+  nIndex = v0->GetNindex();
+  
+  AliESDtrack* d[2];
+  d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
+  d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
+
+  Int_t sign[2];
+  sign[0] = (int)d[0]->GetSign();
+  sign[1] = (int)d[1]->GetSign();
+  
+  if(-1 == sign[0] && 1 == sign[1]){
+    correct = kFALSE;
+    //v0->SetIndex(0, pIndex);  // set the index of the negative v0 track
+    //v0->SetIndex(1, nIndex);  // set the index of the positive v0 track
+  }
+  else{
+    correct = kTRUE;
+  }
+  
+  //pIndex = v0->GetPindex();
+  //nIndex = v0->GetNindex();
+  //printf("-D2: P: %i, N: %i\n", pIndex, nIndex);
+
+  return correct;
+}
diff --git a/PWG3/hfe/AliHFEV0cuts.h b/PWG3/hfe/AliHFEV0cuts.h
new file mode 100644 (file)
index 0000000..65a8dc9
--- /dev/null
@@ -0,0 +1,87 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Class for the V0 cuts - tuned to obtain clean eletron, pion and proton samples.
+// NOT suitable for V0 analysis
+//
+#ifndef ALIHFEV0CUTS_H
+#define ALIHFEV0CUTS_H
+
+#include "AliHFEcollection.h"
+
+class TList;
+
+class AliMCParticle;
+class AliVEvent;
+class AliMCEvent;
+class AliESDtraack;
+class AliESDv0;
+class AliKFVertex;
+class AliKFParticle;
+class AliVTrack;
+
+class AliHFEV0cuts : public TObject {
+ public:
+  AliHFEV0cuts();
+  ~AliHFEV0cuts();
+  AliHFEV0cuts(const AliHFEV0cuts &ref);
+  AliHFEV0cuts &operator=(const AliHFEV0cuts &ref);
+
+  void Init(const char* name);
+  
+  void RunQA();
+  void SetMC(Bool_t b)                  { fIsMC = b; };
+  void SetMCEvent(AliMCEvent* const mce)      { fMCEvent = mce; };
+  void SetInputEvent(AliVEvent* const e)      { fInputEvent = e; };
+  void SetPrimaryVertex(AliKFVertex* const v) { fPrimaryVertex = v; };
+  
+  TList* GetList() { return fQA->GetList(); };
+  
+  Bool_t   TrackCutsCommon(AliESDtrack* track);
+  Bool_t   V0CutsCommon(AliESDv0 *v0);
+  Bool_t   GammaCuts(AliESDv0 *v0);
+  Bool_t   K0Cuts(AliESDv0 *v0);
+  Bool_t   LambdaCuts(AliESDv0 *v0, Bool_t &isLambda);
+  Bool_t LooseRejectK0(AliESDv0 * const v0) const;
+  Bool_t LooseRejectLambda(AliESDv0 * const v0) const;
+  Bool_t LooseRejectGamma(AliESDv0 * const v0) const;
+
+  void     Armenteros(AliESDv0 *v0, Float_t val[2]);
+
+  Double_t OpenAngle(AliESDv0 *v0) const;//opening angle between V0 daughters; close to zero for conversions
+  Double_t PsiPair(AliESDv0 *v0);
+  
+  Bool_t   CheckSigns(AliESDv0* const v0);
+  
+  AliKFParticle *CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec);
+
+ private:
+  void Copy(TObject &ref) const;
+      
+ private:
+  
+  AliHFEcollection     *fQA;            // store QA cut histograms
+  AliMCEvent           *fMCEvent;       // MC event
+  AliVEvent            *fInputEvent;    // Input Event
+  AliKFVertex          *fPrimaryVertex; // primary vertex
+
+  Bool_t                      fIsMC; // availability of MC information
+
+  ClassDef(AliHFEV0cuts, 1)
+};
+    
+
+#endif
diff --git a/PWG3/hfe/AliHFEV0info.cxx b/PWG3/hfe/AliHFEV0info.cxx
new file mode 100644 (file)
index 0000000..8ed8b0f
--- /dev/null
@@ -0,0 +1,44 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Stores aditional information about the V0 candidates
+// author: M.Fasel@gsi.de
+//
+//
+#include "AliHFEV0info.h"
+
+ClassImp(AliHFEV0info)
+
+AliHFEV0info::AliHFEV0info():
+  TObject(),
+  fTrack(NULL),
+  fIDpartnerTrack(-1),
+  fIDV0(-1)
+{
+  //
+  // Dummy constructor
+  //
+}
+
+AliHFEV0info::AliHFEV0info(AliVParticle *track, Int_t idPartnerTrack, Int_t v0id):
+  TObject(),
+  fTrack(track),
+  fIDpartnerTrack(idPartnerTrack),
+  fIDV0(v0id)
+{
+  //
+  // Default constructor
+  //
+}
diff --git a/PWG3/hfe/AliHFEV0info.h b/PWG3/hfe/AliHFEV0info.h
new file mode 100644 (file)
index 0000000..522ed65
--- /dev/null
@@ -0,0 +1,46 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Stores aditional information about the V0 candidates
+// author: M.Fasel@gsi.de
+//
+#ifndef ALIHFEV0INF0_H
+#define ALIHFEV0INFO_H
+
+#ifndef ROOT_TObject
+#include <TObject.h>
+#endif
+
+class AliVParticle;
+
+class AliHFEV0info : public TObject{
+  public:
+    AliHFEV0info();
+    AliHFEV0info(AliVParticle *track, Int_t idPartnerTrack, Int_t v0id);
+    ~AliHFEV0info() {};
+
+    AliVParticle *GetTrack() const { return fTrack; }
+    Int_t GetPartnerID() const { return fIDpartnerTrack; }
+    Int_t GetV0ID() const { return fIDV0; }
+  private:
+    AliHFEV0info(const AliHFEV0info &);
+    AliHFEV0info &operator=(const AliHFEV0info &);
+    AliVParticle *fTrack;       // The track itself
+    Int_t fIDpartnerTrack;      // The ID of the parter
+    Int_t fIDV0;                // the V0 ID
+
+    ClassDef(AliHFEV0info, 1)
+};
+#endif
diff --git a/PWG3/hfe/AliHFEspectrum.cxx b/PWG3/hfe/AliHFEspectrum.cxx
new file mode 100644 (file)
index 0000000..b747afe
--- /dev/null
@@ -0,0 +1,822 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: Markus Fasel <M.Fasel@gsi.de>
+//
+
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+
+
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFEspectrum.h"
+
+ClassImp(AliHFEspectrum)
+
+//____________________________________________________________
+AliHFEspectrum::AliHFEspectrum(const char *name):
+  TNamed(name, ""),
+  fCFContainers(NULL),
+  fTemporaryObjects(NULL),
+  fCorrelation(NULL),
+  fBackground(NULL),
+  fBackgroundSource(kMCbackground),
+  fDumpToFile(kFALSE),
+  fNEvents(0),
+  fStepMC(-1),
+  fStepTrue(-1),
+  fStepData(-1),
+  fStepGuessedUnfolding(-1),
+  fNumberOfIterations(5)
+{
+  //
+  // Default constructor
+  //
+}
+
+//____________________________________________________________
+AliHFEspectrum::~AliHFEspectrum(){
+  //
+  // Destructor
+  //
+  if(fCFContainers) delete fCFContainers;
+  if(fTemporaryObjects){
+    fTemporaryObjects->Clear();
+    delete fTemporaryObjects;
+  }
+}
+
+//____________________________________________________________
+void AliHFEspectrum::Correct(AliCFContainer *datacontainer,AliCFContainer *mccontainer,THnSparseF *mccorrelation){
+  //
+  // Correct the spectrum for efficiency and unfolding
+  // with both method and compare
+  //
+  
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(1111);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadLeftMargin(0.13);
+  gStyle->SetPadRightMargin(0.13);
+  
+  //////////////////
+  // Take only pt
+  /////////////////
+  
+  AliCFDataGrid *dataspectrum = (AliCFDataGrid *)GetSpectrum(datacontainer,20);
+  
+  //
+
+  SetCorrelation(mccorrelation);
+  SetContainer(datacontainer,AliHFEspectrum::kDataContainer);
+  SetContainer(mccontainer,AliHFEspectrum::kMCContainer);
+  
+  SetMCEffStep(8);
+  SetMCTruthStep(0);
+  SetNumberOfIteration(5);
+  SetStepGuessedUnfolding(16);
+  SetNumberOfIteration(5);
+  SetStepToCorrect(20);
+  
+  // Unfold
+
+  TList *listunfolded = Unfold(1);
+  if(!listunfolded){
+    printf("Unfolded failed\n");
+    return;
+  }
+  THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
+  THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
+  if(!correctedspectrum){
+    AliError("No corrected spectrum\n");
+    return;
+  }
+ if(!residualspectrum){
+    AliError("No residul spectrum\n");
+    return;
+  }   
+
+  // Simply correct
+ SetMCEffStep(20);  
+ AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(1);
+ //////////
+  // Plot
+  //////////
+
+  TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
+  cresidual->Divide(2,1);
+  cresidual->cd(1);
+  gPad->SetLogy();
+  TGraphErrors* residualspectrumD = Normalize(residualspectrum);
+  if(!residualspectrumD) {
+    AliError("Number of Events not set for the normalization");
+    return;
+  }
+  residualspectrumD->SetTitle("");
+  residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  residualspectrumD->SetMarkerStyle(26);
+  residualspectrumD->SetMarkerColor(kBlue);
+  residualspectrumD->SetLineColor(kBlue);
+  residualspectrumD->Draw("AP");
+  TGraphErrors* measuredspectrumD = Normalize(dataspectrum);
+  measuredspectrumD->SetTitle("");  
+  measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  measuredspectrumD->SetMarkerStyle(25);
+  measuredspectrumD->SetMarkerColor(kBlack);
+  measuredspectrumD->SetLineColor(kBlack);
+  measuredspectrumD->Draw("P");
+  TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
+  legres->AddEntry(residualspectrumD,"Residual","p");
+  legres->AddEntry(measuredspectrumD,"Measured","p");
+  legres->Draw("same");
+  cresidual->cd(2);
+  TH1D *residualTH1D = residualspectrum->Projection(0);
+  TH1D *measuredTH1D = dataspectrum->Project(0);
+  TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
+  ratioresidual->SetName("ratioresidual");
+  ratioresidual->SetTitle("");
+  ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
+  ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
+  ratioresidual->SetStats(0);
+  ratioresidual->Draw();
+
+  //
+
+  
+  TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
+  ccorrected->Divide(2,1);
+  ccorrected->cd(1);
+  gPad->SetLogy();
+  TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
+  correctedspectrumD->SetTitle("");
+  correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  correctedspectrumD->SetMarkerStyle(26);
+  correctedspectrumD->SetMarkerColor(kBlue);
+  correctedspectrumD->SetLineColor(kBlue);
+  correctedspectrumD->Draw("AP");
+  TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
+  alltogetherspectrumD->SetTitle("");
+  alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  alltogetherspectrumD->SetMarkerStyle(25);
+  alltogetherspectrumD->SetMarkerColor(kBlack);
+  alltogetherspectrumD->SetLineColor(kBlack);
+  alltogetherspectrumD->Draw("P");
+  TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
+  legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
+  legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
+  legcorrected->Draw("same");
+  ccorrected->cd(2);
+  TH1D *correctedTH1D = correctedspectrum->Projection(0);
+  TH1D *alltogetherTH1D = alltogetherCorrection->Project(0);
+  TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
+  ratiocorrected->SetName("ratiocorrected");
+  ratiocorrected->SetTitle("");
+  ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
+  ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
+  ratiocorrected->SetStats(0);
+  ratiocorrected->Draw();
+  
+  // Efficiency correction
+
+  AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(mccontainer,8,7);
+  AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(mccontainer,7,0);
+  AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(mccontainer,8,0);
+  
+  AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(mccontainer,20,0);
+  
+  TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
+  cefficiency->cd(1);
+  TH1D* efficiencymcPIDD = efficiencymcPID->Project(0);
+  efficiencymcPIDD->SetTitle("");
+  efficiencymcPIDD->SetStats(0);
+  efficiencymcPIDD->SetMarkerStyle(25);
+  efficiencymcPIDD->Draw();
+  TH1D* efficiencymctrackinggeoD = efficiencymctrackinggeo->Project(0);
+  efficiencymctrackinggeoD->SetTitle("");
+  efficiencymctrackinggeoD->SetStats(0);
+  efficiencymctrackinggeoD->SetMarkerStyle(26);
+  efficiencymctrackinggeoD->Draw("same");
+  TH1D* efficiencymcallD = efficiencymcall->Project(0);
+  efficiencymcallD->SetTitle("");
+  efficiencymcallD->SetStats(0);
+  efficiencymcallD->SetMarkerStyle(27);
+  efficiencymcallD->Draw("same");
+  TH1D* efficiencyesdallD = efficiencyesdall->Project(0);
+  efficiencyesdallD->SetTitle("");
+  efficiencyesdallD->SetStats(0);
+  efficiencyesdallD->SetMarkerStyle(24);
+  efficiencyesdallD->Draw("same");
+  TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
+  legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
+  legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
+  legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
+  legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
+  legeff->Draw("same");
+
+  // Dump to file if needed
+
+  if(fDumpToFile) {
+    
+    TFile *out = new TFile("finalSpectrum.root","recreate");
+    out->cd();
+    //
+    residualspectrumD->SetName("UnfoldingResidualSpectrum");
+    residualspectrumD->Write();
+    measuredspectrumD->SetName("MeasuredSpectrum");
+    measuredspectrumD->Write();
+    ratioresidual->SetName("RatioResidualSpectrum");
+    ratioresidual->Write();
+    //
+    correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
+    correctedspectrumD->Write();
+    alltogetherspectrumD->SetName("AlltogetherSpectrum");
+    alltogetherspectrumD->Write();
+    ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
+    ratiocorrected->Write();
+    //
+    correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
+    correctedspectrum->Write();
+    alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
+    alltogetherCorrection->Write();
+    //
+    out->Close(); delete out;
+    
+  }
+
+
+}
+
+//____________________________________________________________
+AliCFDataGrid* AliHFEspectrum::SubtractBackground(Int_t dimensions, Bool_t setBackground){
+  //
+  // Apply background subtraction
+  //
+
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  if(!dataContainer){
+    AliError("Data Container not available");
+    return NULL;
+  }
+  // Get the data grid in the requested format
+  Bool_t initContainer = kFALSE;
+  Int_t dims[3];
+  switch(dimensions){
+    case 1:   dims[0] = 0;
+              initContainer = kTRUE;
+              break;
+    case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
+              initContainer = kTRUE;
+              break;
+    default:
+              AliError("Container with this number of dimensions not foreseen (yet)");
+  };
+  if(!initContainer){
+    AliError("Creation of the sliced containers failed. Background estimation not possible");
+    return NULL;
+  }
+  AliCFContainer *spectrumSliced = GetSlicedContainer(dataContainer, dimensions, dims);
+
+  // Make Background Estimate
+  AliCFDataGrid *backgroundGrid = NULL;
+  if(fBackgroundSource == kDataBackground)
+    backgroundGrid = MakeBackgroundEstimateFromData(dimensions);
+  else
+    backgroundGrid = MakeBackgroundEstimateFromMC(dimensions);
+  if(!backgroundGrid){
+    AliError("Background Estimate not created");
+    ClearObject(spectrumSliced);
+    return NULL;
+  }
+
+
+  AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction");
+  spectrumSubtracted->ApplyBGCorrection(*backgroundGrid);
+  if(setBackground){
+    if(fBackground) delete fBackground;
+    fBackground = backgroundGrid;
+  } else delete backgroundGrid;
+
+  return spectrumSubtracted;
+}
+
+//____________________________________________________________
+TList *AliHFEspectrum::Unfold(Int_t dimensions,AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Unfold and eventually correct for efficiency the bgsubspectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainer);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  if(!fCorrelation){
+    AliError("No Correlation map available");
+    return NULL;
+  }
+
+  // Get the mc grid in the requested format
+  Bool_t initContainer = kFALSE;
+  Int_t dims[3];
+  switch(dimensions){
+    case 1:   dims[0] = 0;
+              initContainer = kTRUE;
+              break;
+    case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
+              initContainer = kTRUE;
+              break;
+    default:
+              AliError("Container with this number of dimensions not foreseen (yet)");
+  };
+  if(!initContainer){
+    AliError("Creation of the sliced containers failed. Background estimation not possible");
+    return NULL;
+  }
+  AliCFContainer *mcContainerD = GetSlicedContainer(mcContainer, dimensions, dims);
+  THnSparse *correlationD = GetSlicedCorrelation(dimensions, dims);
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    if(bgsubpectrum->GetNVar()!= dimensions) {
+      AliError("Not the expected number of dimensions for the AliCFDataGrid\n");
+      return NULL;
+    }
+    dataGrid = bgsubpectrum;
+    
+  }
+  else {
+
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    AliCFContainer *dataContainerD = GetSlicedContainer(dataContainer, dimensions, dims);
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainerD);
+    dataGrid->SetMeasured(fStepData);
+    
+  } 
+  
+  
+  // Guessed
+  AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainerD);
+  guessedGrid->SetMeasured(fStepGuessedUnfolding);
+  THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainerD);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+  // Unfold 
+  
+  AliCFUnfolding unfolding("unfolding","",dimensions,correlationD,efficiencyD->GetGrid(),((AliCFGridSparse*)dataGrid->GetData())->GetGrid(),guessedTHnSparse);
+  unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
+  unfolding.UseSmoothing();
+  unfolding.Unfold();
+
+  // Results
+  THnSparse* result = unfolding.GetUnfolded();
+  THnSparse* residual = unfolding.GetEstMeasured();
+  TList *listofresults = new TList;
+  listofresults->SetOwner();
+  listofresults->AddAt((THnSparse*)result->Clone(),0);
+  listofresults->AddAt((THnSparse*)residual->Clone(),1);
+
+  return listofresults;
+
+}
+
+//____________________________________________________________
+AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(Int_t dimensions,AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Apply unfolding and efficiency correction together to bgsubspectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainer);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  // Get the data grid in the requested format
+  Bool_t initContainer = kFALSE;
+  Int_t dims[3];
+  switch(dimensions){
+    case 1:   dims[0] = 0;
+              initContainer = kTRUE;
+              break;
+    case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
+              initContainer = kTRUE;
+              break;
+    default:
+              AliError("Container with this number of dimensions not foreseen (yet)");
+  };
+  if(!initContainer){
+    AliError("Creation of the sliced containers failed. Background estimation not possible");
+    return NULL;
+  }
+  AliCFContainer *mcContainerD = GetSlicedContainer(mcContainer, dimensions, dims);
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainerD);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    if(bgsubpectrum->GetNVar()!= dimensions) {
+      AliError("Not the expected number of dimensions for the AliCFDataGrid\n");
+      return NULL;
+    }
+    dataGrid = bgsubpectrum;
+    
+  }
+  else {
+
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    AliCFContainer *dataContainerD = GetSlicedContainer(dataContainer, dimensions, dims);
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainerD);
+    dataGrid->SetMeasured(fStepData);
+    
+  } 
+
+  // Correct
+  AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
+  result->ApplyEffCorrection(*efficiencyD);
+  
+
+  return result;
+
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+  if(fNEvents > 0) {
+    
+    TH1D* projection = spectrum->Projection(0);
+    CorrectFromTheWidth(projection);
+    TGraphErrors *graphError = NormalizeTH1(projection);
+    return graphError;
+  
+  }
+    
+  return 0x0;
+  
+
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+  if(fNEvents > 0) {
+
+    TH1D* projection = spectrum->Project(0);
+    CorrectFromTheWidth(projection);
+    TGraphErrors *graphError = NormalizeTH1(projection);
+    return graphError;
+    
+  }
+
+  return 0x0;
+  
+
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+  if(fNEvents > 0) {
+
+    TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
+    Double_t p = 0, dp = 0; Int_t point = 1;
+    Double_t n = 0, dN = 0;
+    Double_t nCorr = 0, dNcorr = 0;
+    Double_t errdN = 0, errdp = 0;
+    for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
+      point = ibin - input->GetXaxis()->GetFirst();
+      p = input->GetXaxis()->GetBinCenter(ibin);
+      dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
+      n = input->GetBinContent(ibin);
+      dN = input->GetBinError(ibin);
+
+      // New point
+      nCorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * 1/(2. * TMath::Pi() * p) * n;
+      errdN = 1./(2. * TMath::Pi() * p);
+      errdp = 1./(2. * TMath::Pi() * p*p) * n;
+      dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+      
+      spectrumNormalized->SetPoint(point, p, nCorr);
+      spectrumNormalized->SetPointError(point, dp, dNcorr);
+    }
+    spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
+    spectrumNormalized->SetMarkerStyle(22);
+    spectrumNormalized->SetMarkerColor(kBlue);
+    spectrumNormalized->SetLineColor(kBlue);
+    
+    return spectrumNormalized;
+    
+  }
+
+  return 0x0;
+  
+
+}
+//____________________________________________________________
+void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
+  //
+  // Set the container for a given step to the 
+  //
+  if(!fCFContainers) fCFContainers = new TList;
+  fCFContainers->AddAt(cont, type);
+}
+
+//____________________________________________________________
+AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
+  //
+  // Get Correction Framework Container for given type
+  //
+  if(!fCFContainers) return NULL;
+  return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
+}
+
+//____________________________________________________________
+AliCFDataGrid *AliHFEspectrum::MakeBackgroundEstimateFromData(Int_t nDim){
+  // 
+  // Make Background Estimate from Data
+  // Container having background source from 
+  // Correction has to be done for background selection Efficiency
+  //
+  AliCFContainer *backgroundContainer = GetContainer(kBackgroundMC);
+  AliCFContainer *dataContainer = GetContainer(kBackgroundData);
+  if(!backgroundContainer){
+    AliError("MC background container not found");
+    return NULL;
+  }
+  if(!dataContainer){
+    AliError("Data container not found");
+    return NULL;
+  }
+  AliInfo(Form("Making background Estimate from Data in %d Dimensions", nDim));
+
+  Bool_t initContainer = kFALSE;
+  Int_t dimensions[3];
+  switch(nDim){
+    case 1:   dimensions[0] = 1;
+              initContainer = kTRUE;
+              break;
+    case 3:   for(Int_t i = 0; i < 3; i++) dimensions[i] = i;
+              initContainer = kTRUE;
+              break;
+    default:
+              AliError("Container with this number of dimensions not foreseen (yet)");
+  };
+  if(!initContainer){
+    AliError("Creation of the sliced containers failed. Background estimation not possible");
+    return NULL;
+  }
+  AliCFContainer *slicedBackground = GetSlicedContainer(backgroundContainer, nDim, dimensions);
+  AliCFContainer *slicedData = GetSlicedContainer(dataContainer, nDim, dimensions);
+  
+  // Create Efficiency Grid and data grid
+  AliCFEffGrid backgroundRatio("backgroundRatio", "BackgroundRatio", *slicedBackground);
+  backgroundRatio.CalculateEfficiency(2, 1);
+  AliCFDataGrid *backgroundEstimate = new AliCFDataGrid("backgroundEstimate", "Grid for Background Estimate", *slicedData);
+  backgroundEstimate->SetMeasured(fStepData);
+  backgroundEstimate->ApplyEffCorrection(backgroundRatio);
+
+  return backgroundEstimate;
+}
+
+//____________________________________________________________
+AliCFDataGrid *AliHFEspectrum::MakeBackgroundEstimateFromMC(Int_t nDim){
+  //
+  // Make Background Estimate using MC
+  // Calculate ratio of hadronic background from MC and 
+  // apply this on data
+  //
+  AliCFContainer *backgroundContainer = GetContainer(kBackgroundMC);
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  if(!backgroundContainer){
+    AliError("MC background container not found");
+    return NULL;
+  }
+  if(!dataContainer){
+    AliError("Data container not found");
+    return NULL;
+  }
+  AliInfo(Form("Making background Estimate from MC in %d Dimensions", nDim));
+
+  Bool_t initContainer = kFALSE;
+  Int_t dimensions[3];
+  switch(nDim){
+    case 1:   dimensions[0] = 1;
+              initContainer = kTRUE;
+              break;
+    case 3:   for(Int_t i = 0; i < 3; i++) dimensions[i] = i;
+              initContainer = kTRUE;
+              break;
+    default:
+              AliError("Container with this number of dimensions not foreseen (yet)");
+  };
+  if(!initContainer){
+    AliError("Creation of the sliced containers failed. Background estimation not possible");
+    return NULL;
+  }
+  AliCFContainer *slicedBackground = GetSlicedContainer(backgroundContainer, nDim, dimensions);
+  AliCFContainer *slicedData = GetSlicedContainer(dataContainer, nDim, dimensions);
+  
+  // Create Efficiency Grid and data grid
+  AliCFEffGrid backgroundRatio("backgroundRatio", "BackgroundRatio", *slicedBackground);
+  backgroundRatio.CalculateEfficiency(1, 0);
+  AliCFDataGrid *backgroundEstimate = new AliCFDataGrid("backgroundEstimate", "Grid for Background Estimate", *slicedData);
+  backgroundEstimate->SetMeasured(fStepData);
+  backgroundEstimate->ApplyEffCorrection(backgroundRatio);
+
+  return backgroundEstimate;
+}
+
+//____________________________________________________________
+AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions) {
+  //
+  // Slice Pt bin
+  //
+  
+  Double_t *varMin = new Double_t[container->GetNVar()],
+           *varMax = new Double_t[container->GetNVar()];
+
+  Double_t *binLimits;
+  for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
+    
+    binLimits = new Double_t[container->GetNBins(ivar)+1];
+    container->GetBinLimits(ivar,binLimits);
+    varMin[ivar] = binLimits[0];
+    varMax[ivar] = binLimits[container->GetNBins(ivar)];
+    delete[] binLimits;
+  }
+  
+  AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
+  AddTemporaryObject(k);
+  delete[] varMin; delete[] varMax;
+
+  return k;
+
+}
+
+//_________________________________________________________________________
+THnSparse *AliHFEspectrum::GetSlicedCorrelation(Int_t nDim, Int_t *dimensions) const {
+  //
+  // Slice Pt correlation
+  //
+
+  Int_t ndimensions = fCorrelation->GetNdimensions();
+  printf("Number of dimension %d correlation map\n",ndimensions);
+  if(ndimensions < (2*nDim)) {
+    AliError("Problem in the dimensions");
+    return NULL;
+  }
+  Int_t ndimensionsContainer = (Int_t) ndimensions/2;
+  printf("Number of dimension %d container\n",ndimensionsContainer);
+
+  Int_t *dim = new Int_t[nDim*2];
+  for(Int_t iter=0; iter < nDim; iter++){
+    dim[iter] = dimensions[iter];
+    dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
+    printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
+  }
+    
+  THnSparse *k = fCorrelation->Projection(nDim*2,dim);
+
+  delete[] dim; 
+  return k;
+  
+}
+//___________________________________________________________________________
+void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
+  //
+  // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
+  //
+
+  TAxis *axis = h1->GetXaxis();
+  Int_t nbinX = h1->GetNbinsX();
+
+  for(Int_t i = 1; i <= nbinX; i++) {
+
+    Double_t width = axis->GetBinWidth(i);
+    Double_t content = h1->GetBinContent(i);
+    Double_t error = h1->GetBinError(i); 
+    h1->SetBinContent(i,content/width); 
+    h1->SetBinError(i,error/width);
+  }
+
+}
+//____________________________________________________________
+void AliHFEspectrum::AddTemporaryObject(TObject *o){
+  // 
+  // Emulate garbage collection on class level
+  // 
+  if(!fTemporaryObjects) {
+    fTemporaryObjects= new TList;
+    fTemporaryObjects->SetOwner();
+  }
+  fTemporaryObjects->Add(o);
+}
+
+//____________________________________________________________
+void AliHFEspectrum::ClearObject(TObject *o){
+  //
+  // Do a safe deletion
+  //
+  if(fTemporaryObjects){
+    if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
+    delete o;
+  } else{
+    // just delete
+    delete o;
+  }
+}
+//_________________________________________________________________________
+TObject* AliHFEspectrum::GetSpectrum(AliCFContainer * const c, Int_t step) {
+  AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
+  data->SetMeasured(step);
+  return data;
+}
+//_________________________________________________________________________
+TObject* AliHFEspectrum::GetEfficiency(AliCFContainer * const c, Int_t step, Int_t step0) {
+  // 
+  // Create efficiency grid and calculate efficiency
+  // of step to step0
+  //
+  TString name("eff");
+  name += step;
+  name+= step0;
+  AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
+  eff->CalculateEfficiency(step,step0);
+  return eff;
+}
diff --git a/PWG3/hfe/AliHFEspectrum.h b/PWG3/hfe/AliHFEspectrum.h
new file mode 100644 (file)
index 0000000..35dfb81
--- /dev/null
@@ -0,0 +1,113 @@
+        /**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFESPECTRUM_H
+#define ALIHFESPECTRUM_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+class TGraphErrors;
+class TObject;
+class TH1;
+class TList;
+class AliCFContainer;
+class AliCFDataGrid;
+class AliCFEffGrid;
+
+class AliHFEspectrum : public TNamed{
+  public:
+    enum CFContainer_t{
+      kDataContainer = 0,
+      kMCContainer = 1,
+      kBackgroundData = 2,
+      kBackgroundMC = 3
+    };
+    enum BackgroundSource_t{
+      kMCbackground = 0,
+      kDataBackground = 1
+    };
+    AliHFEspectrum(const char* name);
+    ~AliHFEspectrum();
+
+    void Correct(AliCFContainer *datacontainer,AliCFContainer *mccontainer,THnSparseF *mccorrelation);
+    AliCFDataGrid *SubtractBackground(Int_t dimensions, Bool_t setBackground = kFALSE);
+    
+    TList *Unfold(Int_t dimensions, AliCFDataGrid* const bgsubpectrum = 0x0);
+    AliCFDataGrid *CorrectForEfficiency(Int_t dimensions, AliCFDataGrid* const bgsubpectrum = 0x0);
+    TGraphErrors *Normalize(THnSparse * const spectrum) const;
+    TGraphErrors *Normalize(AliCFDataGrid * const spectrum) const;
+    void SetCorrelation(THnSparseF * const correlation) {fCorrelation = correlation; };
+    void SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type);
+    void SetNumberOfEvents(Int_t nEvents) { fNEvents = nEvents; }
+    void SetBackgroundSource(BackgroundSource_t source) { fBackgroundSource = source; };
+    void SetMCEffStep(Int_t step) { fStepMC = step; };
+    void SetMCTruthStep(Int_t step) { fStepTrue = step; };
+    void SetStepToCorrect(Int_t step) { fStepData = step; };
+    void SetStepGuessedUnfolding(Int_t stepGuessedUnfolding) { fStepGuessedUnfolding = stepGuessedUnfolding; };
+    void SetNumberOfIteration(Int_t numberOfIteration) { fNumberOfIterations = numberOfIteration; };
+    void SetDumpToFile(Bool_t dumpToFile) { fDumpToFile=dumpToFile; }; 
+  
+
+  protected:
+    AliCFDataGrid *MakeBackgroundEstimateFromData(Int_t nDimensions);
+    AliCFDataGrid *MakeBackgroundEstimateFromMC(Int_t nDimensions);
+    
+    AliCFContainer *GetContainer(AliHFEspectrum::CFContainer_t contt);
+    AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions);
+    THnSparse *GetSlicedCorrelation(Int_t nDim, Int_t *dimensions) const;
+    TObject* GetSpectrum(AliCFContainer * const c, Int_t step);
+    TObject* GetEfficiency(AliCFContainer * const c, Int_t step, Int_t step0);
+    void AddTemporaryObject(TObject *cont);
+    void ClearObject(TObject *o);
+    
+    void CorrectFromTheWidth(TH1D *h1) const;
+    TGraphErrors *NormalizeTH1(TH1 *input) const;
+
+
+  private:
+    AliHFEspectrum(const AliHFEspectrum &);
+    AliHFEspectrum &operator=(const AliHFEspectrum &);
+    TList *fCFContainers;         // List of Correction Framework Containers
+    TList *fTemporaryObjects;     // Emulate garbage collection
+    THnSparseF *fCorrelation;     // Correlation Matrices
+    AliCFDataGrid *fBackground;   // Background Grid
+
+    BackgroundSource_t fBackgroundSource;   // Source for the background estimate
+
+    Bool_t fDumpToFile;           // Write Result in a file
+
+    Int_t fNEvents;               // Number of Events
+    Int_t fStepMC;                // MC step (for unfolding)
+    Int_t fStepTrue;              // MC step of the final spectrum
+    Int_t fStepData;              // Data Step (various applications)
+    Int_t fStepGuessedUnfolding;  // Step for first guessed unfolding
+    Int_t fNumberOfIterations;    // Number of iterations
+
+
+    ClassDef(AliHFEspectrum, 1) 
+};
+#endif
+
diff --git a/PWG3/hfe/AliHFEtrdPIDqa.cxx b/PWG3/hfe/AliHFEtrdPIDqa.cxx
new file mode 100644 (file)
index 0000000..95a522b
--- /dev/null
@@ -0,0 +1,755 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// QA class for TRD PID
+// Plot Pion Efficiency at x electron efficiency
+// Calculate the threshold parametrisation and save
+// them in a root file
+//
+// Author:
+//   Markus Fasel <M.Fasel@gsi.de>
+//
+#include <TAxis.h>
+#include <TClass.h>
+#include <TCanvas.h>
+#include <TF1.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <THnSparse.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TIterator.h>
+#include <TLegend.h>
+#include <TList.h>
+#include <TMath.h>
+#include <TObjArray.h>
+#include <TString.h>
+
+#include "AliAODTrack.h"
+#include "AliAODPid.h"
+#include "AliESDtrack.h"
+#include "AliHFEtrdPIDqa.h"
+#include "AliHFEtools.h"
+#include "AliHFEpidTRD.h"
+#include "AliLog.h"
+
+ClassImp(AliHFEtrdPIDqa)
+
+const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
+//_______________________________________________________________
+// Definition of the common binning 
+const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
+  AliPID::kSPECIES + 1,         // species
+  40,                           // p-bins
+  AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
+};
+const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
+  -1,      // species
+  0.1,     // p-bins
+  0        // tracklets including 0
+};
+
+const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
+  AliPID::kSPECIES,               // species
+  10.,                            // p-bins
+  AliESDtrack::kTRDnPlanes + 1    // tracklets including 0
+};
+//_______________________________________________________________
+
+//__________________________________________________________________
+AliHFEtrdPIDqa::AliHFEtrdPIDqa():
+  TNamed("trdPIDqa", ""),
+  fTRDpid(NULL),
+  fLikeTRD(NULL),
+  fQAtrack(NULL),
+  fQAdEdx(NULL),
+  fTRDtruncMean(NULL),
+  fPionEfficiencies(NULL),
+  fProtonEfficiencies(NULL),
+  fKaonEfficiencies(NULL),
+  fThresholds(NULL)
+{
+  //
+  // Default Constructor
+  //
+}
+
+//__________________________________________________________________
+AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
+  TNamed(name, ""),
+  fTRDpid(NULL),
+  fLikeTRD(NULL),
+  fQAtrack(NULL),
+  fQAdEdx(NULL),
+  fTRDtruncMean(NULL),
+  fPionEfficiencies(NULL),
+  fProtonEfficiencies(NULL),
+  fKaonEfficiencies(NULL),
+  fThresholds(NULL)
+{
+  //
+  // Main Constructor
+  //
+}
+
+//__________________________________________________________________
+AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
+  TNamed(ref),
+  fTRDpid(NULL),
+  fLikeTRD(NULL),
+  fQAtrack(NULL),
+  fQAdEdx(NULL),
+  fTRDtruncMean(NULL),
+  fPionEfficiencies(NULL),
+  fProtonEfficiencies(NULL),
+  fKaonEfficiencies(NULL),
+  fThresholds(NULL)
+{
+  //
+  // Copy constructor
+  //
+  ref.Copy(*this);
+}
+
+//__________________________________________________________________
+AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
+  //
+  // Assignment operator
+  //
+  if(this != &ref)
+    ref.Copy(*this);
+  return *this;
+}
+
+//__________________________________________________________________
+AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
+  //
+  // Destructor
+  //
+  if(fTRDpid) delete fTRDpid;
+  if(fLikeTRD) delete fLikeTRD;
+  if(fQAtrack) delete fQAtrack;
+  if(fQAdEdx) delete fQAdEdx;
+  if(fTRDtruncMean) delete fTRDtruncMean;
+  
+  if(fPionEfficiencies) delete fPionEfficiencies;
+  if(fProtonEfficiencies) delete fProtonEfficiencies;
+  if(fKaonEfficiencies) delete fKaonEfficiencies;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::Copy(TObject &ref) const{
+  //
+  // Copies content of this object into object ref
+  //
+  TNamed::Copy(ref);
+
+  AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
+  target.fTRDpid = fTRDpid;
+  target.fLikeTRD = dynamic_cast<THnSparseF *>(fLikeTRD->Clone());
+  target.fQAtrack = dynamic_cast<THnSparseF *>(fQAtrack->Clone());
+  target.fQAdEdx = dynamic_cast<THnSparseF *>(fQAdEdx->Clone());
+  target.fTRDtruncMean = dynamic_cast<THnSparseF *>(fTRDtruncMean->Clone());
+}
+
+//__________________________________________________________________
+Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
+  //
+  // Merge objects
+  //
+  if(!coll) return 0;
+  if(coll->IsEmpty()) return 1;
+  
+  TIterator *it = coll->MakeIterator();
+  TObject *o = NULL;
+  Long64_t count = 0;
+  while((o = it->Next())){
+    AliHFEtrdPIDqa *refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
+    if(!refQA) continue;
+
+    if(fLikeTRD) fLikeTRD->Add(refQA->fLikeTRD);
+    if(fQAtrack) fQAtrack->Add(refQA->fQAtrack);
+    if(fQAdEdx) fQAdEdx->Add(refQA->fQAdEdx);
+    if(fTRDtruncMean) fTRDtruncMean->Add(refQA->fTRDtruncMean);
+    count++;
+  }
+  delete it;
+  return count+1;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::Init(){
+  //
+  // Initialize Object
+  //
+
+  CreateLikelihoodHistogram();
+  CreateQAHistogram();
+  CreatedEdxHistogram();
+  CreateHistoTruncatedMean();
+
+  fTRDpid = new AliHFEpidTRD("QAtrdPID");
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
+  //
+  // Create Histogram for TRD Likelihood Studies
+  //
+  Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
+  nbins[kElectronLike] = 100;
+  Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMin[kElectronLike] = 0.;      
+  binMax[kElectronLike] = 1.;
+
+  fLikeTRD = new THnSparseF("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
+  Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
+  fLikeTRD->GetAxis(kP)->Set(nbins[kP], binLog);
+  delete[] binLog;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::CreateQAHistogram(){
+  //
+  // Create Histogram for Basic TRD PID QA
+  //
+  AliDebug(1, "Called");
+  Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
+  nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
+  nbins[kNClusters] = 200;
+  Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMin[kNonZeroTrackletCharge] = 0.;      
+  binMin[kNClusters] = 0.;      
+  Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
+  binMax[kNClusters] = 200.;
+
+  fQAtrack = new THnSparseF("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
+  Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
+  fQAtrack->GetAxis(kP)->Set(nbins[kP], binLog);
+  delete[] binLog;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::CreatedEdxHistogram(){
+  //
+  // Create QA histogram for dEdx investigations
+  //
+  AliDebug(1, "Called");
+  Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
+  nbins[kdEdx] = 100;
+  Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMin[kdEdx] = 0.;      
+  Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMax[kdEdx] = 100000.;
+
+  fQAdEdx = new THnSparseF("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
+  Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
+  fQAdEdx->GetAxis(kP)->Set(nbins[kP], binLog);
+  delete[] binLog;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
+  //
+  // Create Histogram for Basic TRD PID QA
+  //
+  AliDebug(1, "Called");
+  Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
+  nbins[kTPCdEdx] = 600;
+  nbins[kTRDdEdxMethod1] = 1000;
+  nbins[kTRDdEdxMethod2] = 1000;
+  Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMin[kTPCdEdx] = 0.;      
+  binMin[kTRDdEdxMethod1] = 0.;      
+  binMin[kTRDdEdxMethod2] = 0.;      
+  Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
+  binMax[kTPCdEdx] = 600;
+  binMax[kTRDdEdxMethod1] = 20000.;
+  binMax[kTRDdEdxMethod2] = 20000.;
+
+  fTRDtruncMean = new THnSparseF("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
+  Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
+  fTRDtruncMean->GetAxis(kP)->Set(nbins[kP], binLog);
+  delete[] binLog;
+}
+
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::ProcessTracks(TObjArray * const tracks, Int_t species){
+  //
+  // Process Electron Tracks
+  //
+  if(species < -1 || species >= AliPID::kSPECIES) return;
+  TIterator *it = tracks->MakeIterator();
+  AliVTrack *track = NULL;
+  while((track = dynamic_cast<AliVTrack *>(it->Next()))){
+    if(track) ProcessTrack(track, species);
+  }
+  delete it;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::ProcessTrack(AliVTrack *track, Int_t species){
+  //
+  // Process Single Track
+  //
+  if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
+    ProcessTrackESD(dynamic_cast<AliESDtrack *>(track), species);
+  else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
+    ProcessTrackAOD(dynamic_cast<AliAODTrack *>(track), species);
+}
+
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::ProcessTrackESD(AliESDtrack *track, Int_t species){
+  //
+  // Process single ESD track
+  //
+  if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return;  // require TRD track
+  FillTRDLikelihoods(track, species);
+  FillTRDQAplots(track, species);
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::ProcessTrackAOD(AliAODTrack * const track, Int_t /*species*/){
+  //
+  // Process single AOD track
+  // AOD PID object required
+  //
+  AliAODPid *trackPID = track->GetDetPid();
+  if(!trackPID) return;
+
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::FillTRDLikelihoods(AliESDtrack *track, Int_t species){
+  //
+  // Fill TRD Likelihood Histogram
+  //
+  Double_t trdLike[AliPID::kSPECIES];
+  track->GetTRDpid(trdLike);
+  const AliExternalTrackParam *outerPars = track->GetOuterParam();
+
+  Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
+  // we store:
+  // species
+  // p
+  // ntracklets
+  // Electron Likelihood
+  quantities[kSpecies] = species;
+  quantities[kP] = outerPars ? outerPars->P() : track->P();
+  quantities[kNTracklets] = track->GetTRDntrackletsPID();
+  quantities[kElectronLike] = trdLike[AliPID::kElectron];
+
+  fLikeTRD->Fill(quantities);
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::FillTRDQAplots(AliESDtrack *track, Int_t species){
+  //
+  // Fill QA Plots containing further information
+  //
+  const AliExternalTrackParam *outerPars = track->GetOuterParam();
+
+  Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
+  // we store:
+  // 1. QA
+  // species
+  // p
+  // ntracklets
+  // Non-zero tracklet charges
+  // Number of clusters / full track
+  // 2. dEdx
+  // species
+  // p
+  // ntracklets
+  // dEdx
+  // 3. Truncated Mean
+  // ...
+  // TPC dEdx
+  // TRD dEdx Method 1
+  // TRD dEdx Method 2
+  quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies] 
+                          = quantitiesTruncMean[kSpecies] 
+                          = species;
+  quantitiesQA[kP]  = quantitiesdEdx[kP] 
+                    = quantitiesTruncMean[kP]
+                    = outerPars ? outerPars->P() : track->P();
+  quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets] 
+                            = quantitiesTruncMean[kNTracklets]
+                            = track->GetTRDntrackletsPID();
+  quantitiesQA[kNClusters] = track->GetTRDncls();
+
+  Double_t dEdxSum;
+  Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices();
+  for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
+    dEdxSum = 0.;
+    for(Int_t islice = 0; islice < nSlices; islice++)
+      dEdxSum += track->GetTRDslice(iplane, islice);
+    quantitiesdEdx[kdEdx] = dEdxSum;
+    if(dEdxSum) ntrackletsNonZero++;
+    // Fill dEdx histogram
+    fQAdEdx->Fill(quantitiesdEdx);
+  }
+  quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
+  fQAtrack->Fill(quantitiesQA);
+
+  quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
+  quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, -1);
+  quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, -1);
+  fTRDtruncMean->Fill(quantitiesTruncMean);
+}
+
+/////////////////////////////////////////////////////////
+//
+// Code for Post Processing
+//
+// //////////////////////////////////////////////////////
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::FinishAnalysis(){
+  //
+  // Finish Analysis:
+  // Calculate Electron Efficiency for ntracklets = 4...6
+  // Calculate thresholds for ntracklets = 4...6
+  //
+  
+  if(!fPionEfficiencies){
+    fPionEfficiencies = new TList;
+    fPionEfficiencies->SetName("pionEfficiencies");
+  }
+  if(!fProtonEfficiencies){
+    fProtonEfficiencies = new TList;
+    fProtonEfficiencies->SetName("protonEfficiencies");
+  }
+  if(!fThresholds){
+    fThresholds = new TList;
+    fThresholds->SetName("thresholds");
+  }
+
+  for(Int_t itr = 4; itr <= 6; itr++){
+    printf("========================================\n");
+    printf("Analysing %d trackltes\n", itr);
+    printf("========================================\n");
+    AnalyseNTracklets(itr);
+  }
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
+  //
+  // Store histos into a root file
+  //
+  TFile *outfile = new TFile(filename, "RECREATE");
+  outfile->cd();
+  fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
+  fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
+  fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
+  outfile->Close();
+  delete outfile;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
+  //
+  // Fit the threshold histograms with the given parametrisation
+  // and store the TF1 in the file
+  //
+
+  if(!fThresholds){
+    AliError("Threshold histograms have to be created first");
+    return;
+  }
+
+  printf("========================================\n");
+  printf("Calculating threshold parameters\n");
+  printf("========================================\n");
+
+  TList *outlist = new TList;
+  outlist->SetName("thresholdTRD");
+  
+  TGraph *threshhist = NULL;
+  TF1 *threshparam = NULL;
+  TList *lHistos = NULL, *lFormulas = NULL;
+  for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
+  
+    printf("-------------------------------\n");
+    printf("Processing %d tracklets\n", itracklet);
+    printf("-------------------------------\n");
+
+    lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
+    if(!lHistos){
+      AliError(Form("Threshold histograms for the case %d tracklets not found"));
+      continue;
+    }
+    lFormulas = new TList;
+    lFormulas->SetName(Form("%dTracklets", itracklet));
+    outlist->Add(lFormulas);
+    
+    for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
+      
+      printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
+      printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
+      printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
+
+      threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
+      threshparam = MakeThresholds(threshhist);
+      threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
+    }
+  }
+
+  // store the output
+  TFile *outfile = new TFile(name, "RECREATE");
+  outfile->cd();
+  outlist->Write(outlist->GetName(), kSingleKey);
+  outfile->Close();
+  delete outfile;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
+  //
+  // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
+  // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
+  // electron efficiencies
+  //
+  Int_t binTracklets = fLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
+  fLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
+  
+  Int_t binElectrons = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron); 
+  AliDebug(1, Form("BinElectrons %d", binElectrons));
+  Int_t binPions = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
+  AliDebug(1, Form("BinPions %d", binPions));
+  Int_t binProtons =  fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
+  AliDebug(1, Form("BinProtons %d", binProtons));
+  fLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
+  TH2 *likeElectron = fLikeTRD->Projection(kElectronLike, kP);
+  likeElectron->SetName("likeElectron");
+  fLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
+  TH2 *likePion = fLikeTRD->Projection(kElectronLike, kP);
+  likePion->SetName("likePion");
+  fLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
+  TH2 *likeProton = fLikeTRD->Projection(kElectronLike, kP);
+  likeProton->SetName("likeProton");
+  
+  // Undo ranges
+  fLikeTRD->GetAxis(kSpecies)->SetRange(0, fLikeTRD->GetAxis(kSpecies)->GetNbins());
+  fLikeTRD->GetAxis(kNTracklets)->SetRange(0, fLikeTRD->GetAxis(kNTracklets)->GetNbins());
+
+  // Prepare List for output
+  TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
+  TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
+  TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
+  fPionEfficiencies->Add(listPions);
+  fProtonEfficiencies->Add(listProtons);
+  fThresholds->Add(listThresholds);
+
+  TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
+  TGraphErrors *effPi = NULL, *effPr = NULL; TGraph *thresholds = NULL;
+  Double_t p = 0, dp = 0;
+  Int_t threshbin = 0;
+  Double_t noElEff[2]; // value and error
+  for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
+    printf("-----------------------------------------\n");
+    printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
+    printf("-----------------------------------------\n");
+    effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
+    effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
+    effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
+    effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
+    thresholds = new TGraph(likeElectron->GetXaxis()->GetNbins());
+    thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
+  
+    // Add to lists
+    listPions->Add(effPi);
+    listProtons->Add(effPr);
+    listThresholds->Add(thresholds);
+
+    for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
+      p = likeElectron->GetXaxis()->GetBinCenter(imom);
+      dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
+
+      probsEl = likeElectron->ProjectionY("el", imom);
+      if(!probsEl->GetEntries()) continue;
+      probsEl->Scale(1./probsEl->Integral());
+      probsPi = likePion->ProjectionY("pi", imom);
+      if(!probsPi->GetEntries()) continue;
+      probsPi->Scale(1./probsPi->Integral());
+      probsPr = likeProton->ProjectionY("pr", imom);
+      if(!probsPr->GetEntries()) continue;
+      probsPr->Scale(1./probsPr->Integral());
+      AliDebug(1, Form("Calculating Values for p = %f", p));
+
+      // Calculare threshold we need to achive the x% electron Efficiency
+      threshbin = GetThresholdBin(probsEl, fgkElectronEff[ieff]);
+      thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
+      AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
+
+      // Calculate non-electronEfficiency and error
+      CalculateEfficiency(probsPi, threshbin, noElEff);
+      AliDebug(1, Form("Pion Efficiency %f", noElEff[0]));
+      effPi->SetPoint(imom - 1, p, noElEff[0]);
+      effPi->SetPointError(imom - 1, dp, noElEff[1]);
+      CalculateEfficiency(probsPr, threshbin, noElEff);
+      effPr->SetPoint(imom - 1, p, noElEff[0]);
+      effPr->SetPointError(imom - 1, dp, noElEff[1]);
+      AliDebug(1, Form("Proton Efficiency %f", noElEff[0]));
+      // cleanup
+      delete probsEl;
+      delete probsPi;
+      delete probsPr;
+    }
+  }
+
+  // remove temporary histograms
+  delete likeElectron;
+  delete likePion;
+  delete likeProton;
+}
+
+//__________________________________________________________________
+Int_t AliHFEtrdPIDqa::GetThresholdBin(TH1 * const input, Double_t eff){
+  //
+  // Calculate the threshold bin  
+  //
+  Double_t integralEff = 0.;
+  Int_t currentBin = 0;
+  for(Int_t ibin = input->GetXaxis()->GetLast(); ibin >= input->GetXaxis()->GetFirst(); ibin--){
+    currentBin = ibin;
+    integralEff += input->GetBinContent(ibin);
+    if(integralEff >= eff){
+      // we found the matching bin, break the loop
+      break;
+    }
+  }
+  return currentBin;
+}
+
+//__________________________________________________________________
+Bool_t AliHFEtrdPIDqa::CalculateEfficiency(TH1 * const input, Int_t threshbin, Double_t *par){
+  // 
+  // Calculate non-electron efficiency
+  //
+  Double_t integralEff = 0; 
+  for(Int_t ibin = threshbin; ibin <= input->GetXaxis()->GetLast(); ibin++) 
+    integralEff += input->GetBinContent(ibin);
+  par[0] = integralEff;
+
+  // @TODO: Error calculation
+  par[1] = 0;
+
+  return kTRUE;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet){
+  //
+  // Draw efficiencies and threshold as function of p
+  //
+  if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
+    AliError("No graphs to draw available");  
+    return;
+  }
+
+  TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
+  TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
+  
+  TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet))); 
+
+  TGraphErrors *pi, *pr;
+  TGraph *tr;
+  TLegend *leg;
+  TCanvas *c1 = new TCanvas(Form("tracklet%d", itracklet), Form("Tracklet %d", itracklet));
+  c1->Divide(3,2);
+  for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
+    c1->cd(ieff + 1);
+    leg = new TLegend(0.6, 0.7, 0.89, 0.89);
+    leg->SetBorderSize(0);
+    leg->SetFillStyle(0);
+    pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
+    pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
+    tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
+
+    // Define nice plot, draw
+    // Axis Title
+    pi->GetXaxis()->SetTitle("p / GeV/c");
+    pi->GetYaxis()->SetTitle("Efficiency / %");
+    pr->GetXaxis()->SetTitle("p / GeV/c");
+    pr->GetYaxis()->SetTitle("Efficiency / %");
+    tr->GetXaxis()->SetTitle("p / GeV/c");
+    tr->GetYaxis()->SetTitle("Efficiency / %");
+    // Axis Range
+    pi->GetYaxis()->SetRangeUser(0., 1.);
+    pr->GetYaxis()->SetRangeUser(0., 1.);
+    tr->GetYaxis()->SetRangeUser(0., 1.);
+    // Marker
+    pi->SetMarkerColor(kRed);
+    pi->SetMarkerStyle(20);
+    pr->SetMarkerColor(kBlue);
+    pr->SetMarkerStyle(21);
+    tr->SetMarkerColor(kBlack);
+    tr->SetMarkerStyle(22);
+    // Title
+    pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
+    pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
+    tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
+    // Draw
+    pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("psame");
+
+    // Add entries to legend
+    leg->AddEntry(pi, "Pion Efficiency", "lp");
+    leg->AddEntry(pr, "Proton Efficiency", "lp");
+    leg->AddEntry(tr, "Thresholds", "lp");
+    leg->Draw();
+    c1->Update();
+  }
+}
+
+//__________________________________________________________________
+TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
+  //
+  // Create TF1 containing the threshold parametrisation
+  //
+
+  TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
+  threshist->Fit(threshparam, "NE", "", 0, 10);
+  return threshparam;
+}
+
+//__________________________________________________________________
+void AliHFEtrdPIDqa::ClearLists(){
+  //
+  // Clear lists for particle efficiencies and thresholds
+  //
+  if(fPionEfficiencies){
+    fPionEfficiencies->Delete();
+    delete fPionEfficiencies;
+    fPionEfficiencies = NULL;
+  }
+  if(fProtonEfficiencies){
+    fProtonEfficiencies->Delete();
+    delete fProtonEfficiencies;
+    fProtonEfficiencies = NULL;
+  }
+  if(fThresholds){
+    fThresholds->Delete();
+    delete fThresholds;
+    fThresholds = NULL;
+  }
+}
diff --git a/PWG3/hfe/AliHFEtrdPIDqa.h b/PWG3/hfe/AliHFEtrdPIDqa.h
new file mode 100644 (file)
index 0000000..45b3170
--- /dev/null
@@ -0,0 +1,143 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// QA class for TRD PID
+// Evaluate TRD PID using well identified reference tracks
+// For more information see implementation file
+//
+#ifndef ALIHFETRDPIDQA_H
+#define ALIHFETRDPIDQA_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+#ifndef ROOT_THnSparse
+#include <THnSparse.h>
+#endif
+
+class TCollection;
+class TF1;
+class TGraph;
+class TH1;
+class TList;
+class TObjArray;
+
+class AliAODTrack;
+class AliESDtrack;
+class AliVTrack;
+class AliHFEpidTRD;
+
+class AliHFEtrdPIDqa : public TNamed{
+  public:
+    AliHFEtrdPIDqa();
+    AliHFEtrdPIDqa(const Char_t *name);
+    AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref);
+    AliHFEtrdPIDqa &operator=(const AliHFEtrdPIDqa &ref);
+    virtual void Copy(TObject &o) const;
+    virtual Long64_t Merge(TCollection *coll);
+    virtual ~AliHFEtrdPIDqa();
+
+    void ProcessTracks(TObjArray * const  l, Int_t species);
+    void ProcessTrack(AliVTrack *track, Int_t species);
+
+    void Init();
+    void FinishAnalysis();
+    void StoreResults(const Char_t *filename = "HFEtrdPIDqa.root");
+    void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root");
+
+    void DrawTracklet(Int_t tracklet);
+    void ClearLists();
+
+    //---------------------------------------------------
+    // Getters for Histograms
+    THnSparseF *GetLikelihoodHistogram() const { return fLikeTRD; }
+    THnSparseF *GetQAHistogram() const { return fQAtrack; }
+    THnSparseF *GetdEdxHistogram() const { return fQAdEdx; }
+    THnSparseF *GetHistoTruncMean() const { return fTRDtruncMean; }
+    //---------------------------------------------------
+  protected:
+    // Description of the containers we use to store basic information
+    // we access in the Post Processing. For all containers we have a 
+    // common part containing species, momentum and number of tracklets,
+    // and a specific part. For both containers we define the number of 
+    // variables too
+    enum QuantitiesCommon_t{
+      kSpecies = 0,
+      kP = 1,
+      kNTracklets = 2,
+      kQuantitiesCommon = 3
+    };
+    enum QuantitiesLike_t{
+      kElectronLike = 3,
+      kQuantitiesLike = 4
+    };
+    enum QuantitiesQAtrack_t{
+      kNonZeroTrackletCharge = 3,
+      kNClusters = 4,
+      kQuantitiesQA = 5
+    };
+    enum QuantitiesdEdx_t{
+      kdEdx = 3,
+      kQuantitiesdEdx = 4
+    };
+    enum QuantitiesTruncMean_t{
+      kTPCdEdx = 3,
+      kTRDdEdxMethod1 = 4,
+      kTRDdEdxMethod2 = 5,
+      kQuantitiesTruncMean = 6
+    };
+
+    void ProcessTrackESD(AliESDtrack *track, Int_t species);
+    void ProcessTrackAOD(AliAODTrack * const track, Int_t species);
+
+    void FillTRDLikelihoods(AliESDtrack *track, Int_t species);
+    void FillTRDQAplots(AliESDtrack *track, Int_t species);
+
+    void AnalyseNTracklets(Int_t nTracklets);
+    Int_t GetThresholdBin(TH1 * const input, Double_t efficiency);
+    Bool_t CalculateEfficiency(TH1 * const input, Int_t threshbin, Double_t *params);
+    TF1 *MakeThresholds(TGraph *input);
+
+    void CreateLikelihoodHistogram();
+    void CreateQAHistogram();
+    void CreatedEdxHistogram();
+    void CreateHistoTruncatedMean();
+
+  private:
+    enum{
+      kNElectronEffs = 6
+    };
+    static const Double_t fgkElectronEff[kNElectronEffs];       // Electron efficiency bins
+    static const Int_t    fgkNBinsCommon[kQuantitiesCommon];    // Number of bins for common quantities
+    static const Double_t fgkMinBinCommon[kQuantitiesCommon];   // Bin Limits for common quantities (lower limit)
+    static const Double_t fgkMaxBinCommon[kQuantitiesCommon];   // Bin Limits for common quantities (upper limit)
+    AliHFEpidTRD *fTRDpid;        // HFE PID for TRD
+    THnSparseF *fLikeTRD;         // Histo for Likelihoods
+    THnSparseF *fQAtrack;         // QA histo for quantities based on track level
+    THnSparseF *fQAdEdx;          // QA for tracklet charge
+    THnSparseF *fTRDtruncMean;    // QA for truncated mean
+
+    // List for Histograms:
+    TList *fPionEfficiencies;     //! List for Pion efficiencies
+    TList *fProtonEfficiencies;   //! List for Proton efficiencies
+    TList *fKaonEfficiencies;     //! List for Kaon efficiencies
+
+    TList *fThresholds;           //! List for Threshold Graphs
+  
+  ClassDef(AliHFEtrdPIDqa, 2)     // QA class for TRD PID 
+};
+#endif
+