PWG3/hfe/AliHFEV0cuts.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //  * 20/04/2010 *
  16 // Class for optimising and applying V0 cuts to obtain clean V0 samples
  17 // Compatible with ESDs only
  18 //
  19 // Authors:
  20 //    Matus Kalisky <matus.kalisky@cern.ch>
  21 //
  22
  23 #include "TDatabasePDG.h"
  24
  25 #include "AliESDtrack.h"
  26 #include "AliMCEvent.h"
  27 #include "AliESDv0.h"
  28 #include "AliKFParticle.h"
  29 #include "AliKFVertex.h"
  30 #include "AliLog.h"
  31 #include "AliExternalTrackParam.h"
  32
  33 #include "AliHFEcollection.h"
  34
  35 #include "AliHFEV0cuts.h"
  36
  37 ClassImp(AliHFEV0cuts)
  38
  39 //________________________________________________________________
  40 AliHFEV0cuts::AliHFEV0cuts():
  41   fQA(NULL)
  42   , fQAmc(NULL)
  43   , fMCEvent(NULL)
  44   , fInputEvent(NULL)
  45   , fPrimaryVertex(NULL)
  46   , fCurrentV0id(0)
  47   , fPdaughterPDG(0)
  48   , fNdaughterPDG(0)
  49 {
  50
  51   //
  52   // Default constructor
  53   //
  54
  55
  56 }
  57 //________________________________________________________________
  58 AliHFEV0cuts::~AliHFEV0cuts()
  59 {
  60   //
  61   // destructor
  62   //
  63   if (fQA) delete fQA;
  64   if (fQAmc) delete fQAmc;
  65 }
  66
  67 //________________________________________________________________
  68 AliHFEV0cuts::AliHFEV0cuts(const AliHFEV0cuts &ref):
  69   TObject(ref)
  70   , fQA(NULL)
  71   , fQAmc(NULL)
  72   , fMCEvent(NULL)
  73   , fInputEvent(NULL)
  74   , fPrimaryVertex(NULL)
  75   , fCurrentV0id(0)
  76   , fPdaughterPDG(0)
  77   , fNdaughterPDG(0)
  78 {
  79   //
  80   // Copy constructor
  81   //
  82   ref.Copy(*this);
  83 }
  84 //________________________________________________________________
  85 AliHFEV0cuts &AliHFEV0cuts::operator=(const AliHFEV0cuts &ref){
  86   //
  87   // Assignment operator
  88   //
  89   if(this != &ref)
  90     ref.Copy(*this);
  91   return *this;
  92 }
  93 //________________________________________________________________
  94 void AliHFEV0cuts::Copy(TObject &ref) const{
  95   //
  96   // Copy function
  97   //
  98   AliHFEV0cuts &target = dynamic_cast<AliHFEV0cuts &>(ref);
  99
 100   if(fQA) target.fQA = dynamic_cast<AliHFEcollection *>(fQA->Clone());
 101
 102   if(fQAmc) target.fQAmc = dynamic_cast<AliHFEcollection *>(fQAmc->Clone());
 103
 104   if(target.fMCEvent) delete target.fMCEvent;
 105   target.fMCEvent = new AliMCEvent;
 106
 107   if(target.fPrimaryVertex) delete target.fPrimaryVertex;
 108   target.fPrimaryVertex = new AliKFVertex;
 109
 110   TObject::Copy(ref);
 111
 112 }
 113 //___________________________________________________________________
 114 void AliHFEV0cuts::Init(const char* name){
 115   //
 116   // initialize the output objects and create histograms
 117   //
 118
 119   //
 120   // all the "h_cut_XXX" histograms hare cut value distributions:
 121   // [0] for all candidates
 122   // [1] jus before the cut on given variable was applied, but after all the previous cuts
 123   //
 124
 125   fQA = new AliHFEcollection("fQA", name);
 126
 127   fQAmc = new AliHFEcollection("fQAmc", name);
 128
 129   // common for all V0s
 130   fQA->CreateTH2Fvector1(2, "h_all_AP", "armenteros plot for all V0 candidates", 200, -1, 1, 200, 0, 0.25);
 131
 132   // gammas
 133   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_CosPoint", "Gamma Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
 134   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_DCA", "DCA between the gamma daughters; dca (cm); counts", 100, 0, 2);
 135   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_VtxR_old", "*old* Radius of the gamma conversion vertex; r (cm); counts", 1000, 0, 100);
 136   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_VtxR", "Radius of the gamma conversion vertex; r (cm); counts", 1000, 0, 100);
 137   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_PP", "gamma psi pair angle; psi pairangle (rad); counts", 100, 0, 2);
 138   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_Chi2", "gamma Chi2/NDF; Chi2/NDF; counts", 100, 0, 50);
 139   fQA->CreateTH1Fvector1(2, "h_cut_Gamma_Sep", "gamma separation dist at TPC inned wall", 100, 0, 10);
 140   fQA->CreateTH1Fvector1(7, "h_Gamma_Mass", "Invariant mass of gammas; mass (GeV/c^{2}); counts", 100, 0, 0.2);
 141
 142
 143   // kaons
 144   fQA->CreateTH1Fvector1(2, "h_cut_K0_CosPoint", "K0 Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
 145   fQA->CreateTH1Fvector1(2, "h_cut_K0_DCA", "DCA between the K0 daughters; dca (cm); counts", 100, 0, 2);
 146   fQA->CreateTH1Fvector1(2, "h_cut_K0_VtxR", "Radius of the K0 decay vertex; r (cm); counts", 1000, 0, 100);
 147   fQA->CreateTH1Fvector1(2, "h_cut_K0_Chi2", "K0 Chi2/NDF; Chi2/NDF; counts", 100, 0, 50);
 148   fQA->CreateTH1Fvector1(5, "h_K0_Mass", "Invariant mass of K0; mass (GeV/c^{2}); counts", 125, 0.45, 0.55);
 149
 150   // lambda
 151   fQA->CreateTH1Fvector1(2, "h_cut_L_CosPoint", "L Cosine pointing angle; cos point. angle; counts", 100, 0, 0.1);
 152   fQA->CreateTH1Fvector1(2, "h_cut_L_DCA", "DCA between the L daughters; dca (cm); counts", 100, 0, 2);
 153   fQA->CreateTH1Fvector1(2, "h_cut_L_VtxR", "Radius of the L decay vertex; r (cm); counts", 1000, 0, 100);
 154   fQA->CreateTH1Fvector1(2, "h_cut_L_Chi2", "L Chi2/NDF; Chi2/NDF; counts", 100, 0, 50);
 155   fQA->CreateTH1Fvector1(5, "h_L_Mass", "Invariant mass of L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
 156   fQA->CreateTH1Fvector1(5, "h_AL_Mass", "Invariant mass of anti L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
 157
 158   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 );
 159
 160   // electrons
 161   fQA->CreateTH1Fvector1(7, "h_Electron_P", "Momenta of conversion electrons -cuts-; P (GeV/c); counts", 50, 0.1, 20, 0);
 162
 163   // K0 pions
 164   fQA->CreateTH1Fvector1(8, "h_PionK0_P", "Momenta of K0 pions -cuts-; P (GeV/c) counts;", 50, 0.1, 20, 0);
 165
 166   // L pions
 167   fQA->CreateTH1Fvector1(9, "h_PionL_P", "Momenta of L pions -cuts-; P (GeV/c) counts;", 50, 0.1, 20, 0);
 168
 169   // L protons
 170   fQA->CreateTH1Fvector1(9, "h_ProtonL_P", "Momenta of L protons -cuts-; P (GeV/c) counts;", 50, 0.1, 20, 0);
 171
 172   // single track cuts
 173   fQA->CreateTH1F("h_ST_NclsTPC", "Number of TPC clusters", 161, -1, 160);
 174   fQA->CreateTH1F("h_ST_TPCrefit", "TPC refit", 2, -0.5, 1.5);
 175   fQA->CreateTH1F("h_ST_chi2TPCcls", "chi2 per TPC cluster", 100, 0, 10);
 176   fQA->CreateTH1F("h_ST_TPCclsR", "TPC cluster ratio", 120, -0.1, 1.1);
 177   fQA->CreateTH1F("h_ST_kinks", "kinks", 2, -0.5, 1.5);
 178   fQA->CreateTH1F("h_ST_pt", "track pt", 100, 0.1, 20, 0);
 179   fQA->CreateTH1F("h_ST_eta", "track eta", 100, -1.5, 1.5);
 180
 181   //
 182   // possibly new cuts
 183   //
 184  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);
 185
 186   // THnSparse histograms
 187
 188   // THnSparse for the K0 mass
 189   // to be looked at after merging run by run
 190   // axes: mass, pt, theta, phi
 191   {
 192     Int_t nBin[4] = {100, 10, 10, 18};
 193     Double_t nMin[4] = {0.45, 0.1, 0., 0.};
 194     Double_t nMax[4] = {0.55, 10., TMath::Pi(), 2*TMath::Pi()};
 195     TString htitle = "K0 sparse; mass (GeV/c^{2}); p_{T} (GeV/c); theta (rad); phi(rad)";
 196     fQA->CreateTHnSparse("hK0", htitle, 4, nBin, nMin, nMax);
 197     fQA->BinLogAxis("hK0", 1);
 198   }
 199
 200
 201   //
 202   // MC plots for checking and tuning the V0 cuts
 203   //
 204
 205   const char *v0[4] = {"G", "K", "L"}; // to keep the names short
 206   // number of V0s left after each cut step - for signal and background - within given mass window
 207   for(Int_t i=0; i<3; ++i){
 208     fQAmc->CreateTH1F(Form("h_%s_cuts_S", v0[i]), Form("h_%s_cuts_S", v0[i]), 10, -0.5, 9.5);
 209     fQAmc->CreateTH1F(Form("h_%s_cuts_B", v0[i]), Form("h_%s_cuts_B", v0[i]), 10, -0.5, 9.5);
 210   }
 211
 212   //
 213   // cut distributions for signal and background
 214   //
 215
 216   const Float_t pMin = 0.1;
 217   const Float_t pMax = 10.;
 218   const Int_t pN = 12;
 219
 220
 221   // gamma signal
 222   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_CosPoint_S", "S - Gamma Cosine pointing angle; mom (GeV/c); cos point. angle",  pN, pMin, pMax, 50, 0, 0.1, 0);
 223   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_DCA_S", "S - DCA between the gamma daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 224   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_VtxR_S", "S - Radius of the gamma conversion vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 100, 0, 100, 0);
 225   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_PP_S", "S - gamma psi pair angle; mom (GeV/c); psi pairangle (rad)", pN, pMin, pMax, 50, 0, 0.5, 0);
 226   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_Chi2_S", "S - gamma Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 25, 0);
 227   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_Sep_S", "S - gamma separation TPC-inner; mom (GeV/c); tracks separatin (cm)", pN, pMin, pMax, 100, 0, 10, 0);
 228   // as a function of radius, not momentum
 229   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_SepR_S", "S - gamma separation TPC-inner; radius cm; tracks separatin (cm)", 20, 0, 100, 100, 0, 20);
 230
 231   fQAmc->CreateTH1Fvector1(9, "h_Gamma_Mass_S", "S - Invariant mass of gammas; mass (GeV/c^{2}); counts", 100, 0, 0.2);
 232   // gamma background
 233   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_CosPoint_B", "B - Gamma Cosine pointing angle; mom (GeV/c); cos point. angle", pN, pMin, pMax, 50, 0, 0.1, 0);
 234   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_DCA_B", "B - DCA between the gamma daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 235   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_VtxR_B", "B - Radius of the gamma conversion vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 100, 0, 100, 0);
 236   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_PP_B", "B - gamma psi pair angle; mom (GeV/c); psi pairangle (rad)", pN, pMin, pMax, 50, 0, 0.5, 0);
 237   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_Chi2_B", "B - gamma Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 25, 0);
 238   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_Sep_B", "B - gamma separation TPC-inner; mom (GeV/c); tracks separatin (cm)", pN, pMin, pMax, 100, 0, 50, 0);
 239   //
 240   fQAmc->CreateTH2Fvector1(2, "h_cut_Gamma_SepR_B", "S - gamma separation TPC-inner; radius cm; tracks separatin (cm)", 20, 0, 100, 100, 0, 20);
 241   fQAmc->CreateTH1Fvector1(9, "h_Gamma_Mass_B", "B - Invariant mass of gammas; mass (GeV/c^{2}); counts", 100, 0, 0.2);
 242
 243   // kaons signal
 244   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_CosPoint_S", "S - K0 Cosine pointing angle; mom (GeV/c); cos point. angle", pN, pMin, pMax, 50, 0, 0.1, 0);
 245   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_DCA_S", "S - DCA between the K0 daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 246   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_VtxR_S", "S - Radius of the K0 decay vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 50, 0, 100, 0);
 247   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_Chi2_S", "S - K0 Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 25, 0);
 248
 249   fQAmc->CreateTH1Fvector1(5, "h_K0_Mass_S", "S - Invariant mass of K0; mass (GeV/c^{2}); counts", 125, 0.45, 0.55);
 250   // kaons background
 251   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_CosPoint_B", "B - K0 Cosine pointing angle; mom (GeV/c); cos point. angle", pN, pMin, pMax, 50, 0, 0.1, 0);
 252   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_DCA_B", "B - DCA between the K0 daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 253   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_VtxR_B", "B - Radius of the K0 decay vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 50, 0, 100, 0);
 254   fQAmc->CreateTH2Fvector1(2, "h_cut_K0_Chi2_B", "B - K0 Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 50, 0);
 255
 256   fQAmc->CreateTH1Fvector1(5, "h_K0_Mass_B", "B - Invariant mass of K0; mass (GeV/c^{2}); counts", 125, 0.45, 0.55);
 257
 258   // lambda signal
 259   fQAmc->CreateTH2Fvector1(2, "h_cut_L_CosPoint_S", "S - L Cosine pointing angle; mom (GeV/c); cos point. angle", pN, pMin, pMax, 50, 0, 0.1, 0);
 260   fQAmc->CreateTH2Fvector1(2, "h_cut_L_DCA_S", "S - DCA between the L daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 261   fQAmc->CreateTH2Fvector1(2, "h_cut_L_VtxR_S", "S - Radius of the L decay vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 50, 0, 100, 0);
 262   fQAmc->CreateTH2Fvector1(2, "h_cut_L_Chi2_S", "S - L Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 50, 0);
 263
 264   fQAmc->CreateTH1Fvector1(5, "h_L_Mass_S", "S - Invariant mass of L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
 265   fQAmc->CreateTH1Fvector1(5, "h_AL_Mass_S", "S - Invariant mass of anti L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
 266   // lambda background
 267   fQAmc->CreateTH2Fvector1(2, "h_cut_L_CosPoint_B", "B - L Cosine pointing angle; mom (GeV/c); cos point. angle", pN, pMin, pMax, 50, 0, 0.1, 0);
 268   fQAmc->CreateTH2Fvector1(2, "h_cut_L_DCA_B", "B - DCA between the L daughters; mom (GeV/c); dca (cm)", pN, pMin, pMax, 50, 0, 2, 0);
 269   fQAmc->CreateTH2Fvector1(2, "h_cut_L_VtxR_B", "B - Radius of the L decay vertex; mom (GeV/c); r (cm)", pN, pMin, pMax, 50, 0, 100, 0);
 270   fQAmc->CreateTH2Fvector1(2, "h_cut_L_Chi2_B", "B - L Chi2/NDF; mom (GeV/c); Chi2/NDF", pN, pMin, pMax, 50, 0, 50, 0);
 271   fQAmc->CreateTH2Fvector1(2, "h_cut_L_rdp_v_mp_S", "S - relative L daughter mom -v- mother mom; L mom (GeV/c); relative daughter mom p2/p1", 100, 0.1, 10, 100, 0, 1);
 272   fQAmc->CreateTH2Fvector1(2, "h_cut_L_rdp_v_mp_B", "B - relative L daughter mom -v- mother mom; L mom (GeV/c); relative daughter mom p2/p1", 100, 0.1, 10, 100, 0, 1);
 273   fQAmc->CreateTH1Fvector1(5, "h_LAL_Mass_B", "B - Invariant mass of anti L; mass (GeV/c^{2}); counts", 60, 1.1, 1.13);
 274
 275
 276   // MC tagged daughter track momentum distribution after each cut step
 277 //   fQAmc->CreateTH1Fvector1(10, "h_electron_p_S", "h_electron_p_S", 20, 0.1, 20, 0);
 278 //   fQAmc->CreateTH1Fvector1(10, "h_K0pion_p_S", "h_K0pion_p_S", 20, 0.1, 20, 0);
 279 //   fQAmc->CreateTH1Fvector1(10, "h_Lpion_p_S", "h_Lpion_p_S", 20, 0.1, 20, 0);
 280 //   fQAmc->CreateTH1Fvector1(10, "h_proton_p_S", "h_proton_p_S", 20, 0.1, 20, 0);
 281
 282   // V0 momnetum distribution of MC tagged signal and backglound after all cuts
 283   fQAmc->CreateTH1F("h_gamma_p_S", "true gammas after all cuts", 20, 0.1, 10, 0);
 284   fQAmc->CreateTH1F("h_gamma_p_B", "true gamma BG after all cuts", 20, 0.1, 10, 0);
 285   fQAmc->CreateTH1F("h_K0_p_S", "true K0s after all cuts", 20, 0.1, 10, 0);
 286   fQAmc->CreateTH1F("h_K0_p_B", "true K0 BG after all cuts", 20, 0.1, 10, 0);
 287   fQAmc->CreateTH1F("h_lambda_p_S", "MC true lambdas after all cuts", 20, 0.1, 10, 0);
 288   fQAmc->CreateTH1F("h_lambda_p_B", "MC true lambda BG after all cuts", 20, 0.1, 10, 0);
 289   fQAmc->CreateTH1F("h_alambda_p_S", "MC true anti-lambdas after all cuts", 20, 0.1, 10, 0);
 290   fQAmc->CreateTH1F("h_alambda_p_B", "MC true anti-lambda BG after all cuts", 20, 0.1, 10, 0);
 291
 292   // invariant mass ditributions for the V0 for different hypoteses (gamma, K0, L, AL)
 293   fQAmc->CreateTH1F("h_Mass_gamma_as_K0","h_Mass_gamma_as_K0", 200, 0, 2);
 294   fQAmc->CreateTH1F("h_Mass_gamma_as_L","h_Mass_gamma_as_L", 200, 0, 2);
 295   fQAmc->CreateTH1F("h_Mass_K0_as_G", "h_Mass_K0_as_gamma", 200, 0, 2);
 296   fQAmc->CreateTH1F("h_Mass_K0_as_L", "h_Mass_K0_as_Lambda", 200, 0, 2);
 297   fQAmc->CreateTH1F("h_Mass_L_as_G", "h_Mass_L_as_gamma", 200, 0, 2);
 298   fQAmc->CreateTH1F("h_Mass_L_as_K0", "h_Mass_L_as_K0", 200, 0, 2);
 299
 300   // Invariant mass distribution of MC tagged signal for diffrent momenta
 301   fQAmc->CreateTH2F("h_gamma_MvP_S", "mc tagged gammas - signal; p (GeV/c); m (GeV/c^{2})", 12, 0.1, 20, 100, 0., 0.1, 0);
 302   fQAmc->CreateTH2F("h_K0_MvP_S", "mc tagged K0s - signal; p (GeV/c); m (GeV/c^{2})", 12, 0.1, 20, 100, 0.45, 0.55, 0);
 303   fQAmc->CreateTH2F("h_lambda_MvP_S", "mc tagged Lambdas - signal; p (GeV/c); m (GeV/c^{2})", 12, 0.1, 20, 100, 1.08, 1.14, 0);
 304
 305   // electrons
 306   fQAmc->CreateTH1Fvector1(8, "h_Electron_P_S", "MC-S momenta of conversion electrons -cuts-; P (GeV/c); counts", 20, 0.1, 20, 0);
 307   fQAmc->CreateTH1Fvector1(8, "h_Electron_P_B", "MC-B momenta of conversion electrons -cuts-; P (GeV/c); counts", 20, 0.1, 20, 0);
 308
 309   // K0 pions
 310   fQAmc->CreateTH1Fvector1(7, "h_PionK0_P_S", "MC-S momenta of K0 pions -cuts-; P (GeV/c) counts;", 20, 0.1, 20, 0);
 311   fQAmc->CreateTH1Fvector1(7, "h_PionK0_P_B", "MC-B momenta of K0 pions -cuts-; P (GeV/c) counts;", 20, 0.1, 20, 0);
 312
 313   // L pions
 314   fQAmc->CreateTH1Fvector1(8, "h_PionL_P_S", "MC-S momenta of L pions -cuts-; P (GeV/c) counts;", 20, 0.1, 50, 0);
 315   fQAmc->CreateTH1Fvector1(8, "h_PionL_P_B", "MC-B momenta of L pions -cuts-; P (GeV/c) counts;", 20, 0.1, 50, 0);
 316
 317   // L protons
 318   fQAmc->CreateTH1Fvector1(8, "h_ProtonL_P_S", "MC-S momenta of L protons -cuts-; P (GeV/c) counts;", 20, 0.1, 20, 0);
 319   fQAmc->CreateTH1Fvector1(8, "h_ProtonL_P_B", "MC-B momenta of L protons -cuts-; P (GeV/c) counts;", 20, 0.1, 20, 0);
 320
 321
 322
 323   // cut efficiencies
 324 }
 325 //________________________________________________________________
 326 Bool_t AliHFEV0cuts::TrackCutsCommon(AliESDtrack* track){
 327   //
 328   // singe track cuts commom for all particle candidates
 329   //
 330
 331   if(!track) return kFALSE;
 332
 333
 334   // status word
 335   ULong_t status = track->GetStatus();
 336
 337
 338   // No. of TPC clusters
 339   fQA->Fill("h_ST_NclsTPC", track->GetTPCNcls());
 340   if(track->GetTPCNcls() < 1) return kFALSE;   //
 341
 342   // TPC refit
 343   if((status & AliESDtrack::kTPCrefit)){
 344     fQA->Fill("h_ST_TPCrefit", 1);
 345   }
 346   if(!(status & AliESDtrack::kTPCrefit)){
 347     fQA->Fill("h_ST_TPCrefit", 0);
 348     return kFALSE;
 349   }
 350
 351   // Chi2 per TPC cluster
 352   Int_t nTPCclusters = track->GetTPCclusters(0);
 353   Float_t chi2perTPCcluster = track->GetTPCchi2()/Float_t(nTPCclusters);
 354   fQA->Fill("h_ST_chi2TPCcls", chi2perTPCcluster);
 355   if(chi2perTPCcluster > 4.0) return kFALSE;   // 4.0
 356
 357   // TPC cluster ratio
 358   Float_t cRatioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t> (track->GetTPCNclsF()) : 1.;
 359   fQA->Fill("h_ST_TPCclsR", cRatioTPC);
 360   if(cRatioTPC < 0.6) return kFALSE;
 361
 362   // kinks
 363   fQA->Fill("h_ST_kinks", track->GetKinkIndex(0));
 364   if(track->GetKinkIndex(0) != 0) return kFALSE;
 365
 366   // pt
 367   fQA->Fill("h_ST_pt",track->Pt());
 368   //if(track->Pt() < 0.1 || track->Pt() > 100) return kFALSE; //
 369
 370   // eta
 371   fQA->Fill("h_ST_eta", track->Eta());
 372   //if(TMath::Abs(track->Eta()) > 0.9) return kFALSE;
 373
 374   return kTRUE;
 375 }
 376 //________________________________________________________________
 377 Bool_t AliHFEV0cuts::V0CutsCommon(AliESDv0 *v0){
 378   //
 379   // V0 cuts common to all V0s
 380   //
 381
 382   AliESDtrack* dN, *dP;
 383
 384   dP = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(v0->GetPindex()));
 385   dN = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(v0->GetNindex()));
 386
 387   if(!dN || !dP) return kFALSE;
 388
 389   Int_t qP = dP->Charge();
 390   Int_t qN = dN->Charge();
 391
 392   if((qP*qN) != -1) return kFALSE;
 393
 394   return kTRUE;
 395 }
 396 //________________________________________________________________
 397 Bool_t AliHFEV0cuts::GammaCuts(AliESDv0 *v0){
 398   //
 399   // gamma cuts
 400   //
 401
 402   if(!v0) return kFALSE;
 403
 404   if(fMCEvent){
 405     if(1 == fCurrentV0id){
 406       fQAmc->Fill("h_Mass_gamma_as_K0",  v0->GetEffMass(2, 2));
 407       fQAmc->Fill("h_Mass_gamma_as_L",  v0->GetEffMass(2, 4));
 408       fQAmc->Fill("h_Mass_gamma_as_L",  v0->GetEffMass(4, 2));
 409     }
 410   }
 411
 412   // loose cuts first
 413   //if(LooseRejectK0(v0) || LooseRejectLambda(v0)) return kFALSE;
 414
 415   AliVTrack* daughter[2];
 416   Int_t pIndex = 0, nIndex = 0;
 417   if(CheckSigns(v0)){
 418     pIndex = v0->GetPindex();
 419     nIndex = v0->GetNindex();
 420   }
 421   else{
 422     pIndex = v0->GetNindex();
 423     nIndex = v0->GetPindex();
 424   }
 425   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 426   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 427   if(!daughter[0] || !daughter[1]) return kFALSE;
 428
 429   AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
 430   if(!kfMother) return kFALSE;
 431
 432   // production vertex is set in the 'CreateMotherParticle' function
 433   //kfMother->SetMassConstraint(0, 0.001);
 434
 435   AliESDtrack* d[2];
 436   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
 437   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
 438
 439   Float_t iMass = v0->GetEffMass(0, 0);
 440   Float_t iP = v0->P();
 441   Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
 442
 443   // Cut values
 444   const Double_t cutChi2NDF = 1.5;              // ORG [7.]
 445   const Double_t cutCosPoint[2] = {0., 0.007};  // ORG [0., 0.02]
 446   const Double_t cutDCA[2] = {0., 0.25};       // ORG [0., 0.25]
 447   const Double_t cutProdVtxR[2] = {6., 90.};   // ORG [3., 90]
 448   const Double_t cutPsiPair[2] = {0., 0.05};   // ORG [0. 0.05]
 449   // mass cut
 450   const Double_t cutMass = 0.05;               // ORG [0.05]
 451
 452   //
 453   // possible new cuts
 454   //
 455   // separation cut at the entrance to the TPC
 456   const Double_t cutSeparation = 0.;          // ORG 3.0 cm
 457
 458
 459
 460   // Values
 461
 462   // cos pointing angle
 463   Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
 464   cosPoint = TMath::ACos(cosPoint);
 465
 466   // DCA between daughters
 467   Double_t dca = v0->GetDcaV0Daughters();
 468
 469   // Production vertex
 470   Double_t x, y, z;
 471   v0->GetXYZ(x,y,z);
 472   Double_t r = TMath::Sqrt(x*x + y*y);
 473
 474   Double_t xy[2];
 475   Double_t r2 = -1.;
 476   if ( GetConvPosXY(d[0], d[1], xy) ){
 477     r2 = TMath::Sqrt(xy[0]*xy[0] + xy[1]*xy[1]);
 478   }
 479
 480   // psi pair
 481   Double_t psiPair = PsiPair(v0);
 482
 483   // V0 chi2/ndf
 484   Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
 485   if(kfMother) delete kfMother;
 486
 487   // Separation
 488   AliExternalTrackParam const *param[2];
 489   param[0] = d[0]->GetInnerParam();
 490   param[1] = d[1]->GetInnerParam();
 491   Double_t sep = 999.;
 492   if(param[0] && param[1]){
 493     TVector3 xyz[3];
 494     xyz[0].SetXYZ(param[0]->GetX(), param[0]->GetY(), param[0]->GetZ());
 495     xyz[1].SetXYZ(param[1]->GetX(), param[1]->GetY(), param[1]->GetZ());
 496     xyz[2] = xyz[0] - xyz[1];
 497     sep = xyz[2].Mag();
 498   }
 499
 500
 501
 502   //
 503   // Apply the cuts, produce QA plots (with mass cut)
 504   //
 505   fQA->Fill("h_Gamma_Mass", 0, iMass);
 506
 507   // MC
 508   if(fMCEvent){
 509     if(1 == fCurrentV0id){
 510       fQAmc->Fill("h_Gamma_Mass_S", 0, iMass);
 511       fQAmc->Fill("h_gamma_MvP_S", iP, iMass);
 512     }
 513     else if(-2 != fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_B", 0, iMass);
 514   }
 515   // cut distributions
 516   if(iMass < cutMass){
 517     fQA->Fill("h_Electron_P", 0, p[0]);
 518     fQA->Fill("h_Electron_P", 0, p[1]);
 519     fQA->Fill("h_cut_Gamma_CosPoint", 0, cosPoint);
 520     fQA->Fill("h_cut_Gamma_DCA", 0, dca);
 521     fQA->Fill("h_cut_Gamma_VtxR_old", 0, r);
 522     fQA->Fill("h_cut_Gamma_VtxR", 0, r2);
 523     fQA->Fill("h_cut_Gamma_PP", 0, psiPair);
 524     fQA->Fill("h_cut_Gamma_Chi2", 0, chi2ndf);
 525     fQA->Fill("h_cut_Gamma_Chi2", 1, chi2ndf, iP);
 526     fQA->Fill("h_cut_Gamma_Sep", 0, iP, sep);
 527
 528
 529     if(fMCEvent){
 530       // MC signal
 531       if(1 == fCurrentV0id){
 532         fQAmc->Fill("h_cut_Gamma_CosPoint_S", 0, iP, cosPoint);
 533         fQAmc->Fill("h_cut_Gamma_DCA_S", 0, iP, dca);
 534         fQAmc->Fill("h_cut_Gamma_VtxR_S", 0, iP, r2);
 535         fQAmc->Fill("h_cut_Gamma_PP_S", 0, iP, psiPair);
 536         fQAmc->Fill("h_cut_Gamma_Chi2_S", 0, iP, chi2ndf);
 537         fQAmc->Fill("h_cut_Gamma_Chi2_S", 1, iP, chi2ndf);
 538         fQAmc->Fill("h_cut_Gamma_Sep_S", 0, iP, sep);
 539         fQAmc->Fill("h_cut_Gamma_SepR_S", 0,r2, sep);
 540         fQAmc->Fill("h_Electron_P_S", 0, p[0]);
 541         fQAmc->Fill("h_Electron_P_S", 0, p[1]);
 542       }
 543       // MC background
 544       else if(-2 != fCurrentV0id){
 545         fQAmc->Fill("h_cut_Gamma_CosPoint_B", 0, iP, cosPoint);
 546         fQAmc->Fill("h_cut_Gamma_DCA_B", 0, iP, dca);
 547         fQAmc->Fill("h_cut_Gamma_VtxR_B", 0, iP, r2);
 548         fQAmc->Fill("h_cut_Gamma_PP_B", 0, iP, psiPair);
 549         fQAmc->Fill("h_cut_Gamma_Chi2_B", 0, iP, chi2ndf);
 550         fQAmc->Fill("h_cut_Gamma_Chi2_B", 1, iP, chi2ndf);
 551         fQAmc->Fill("h_cut_Gamma_Sep_B", 0, iP, sep);
 552         fQAmc->Fill("h_cut_Gamma_SepR_B", 0,r2, sep);
 553         fQAmc->Fill("h_Electron_P_B", 0, p[0]);
 554         fQAmc->Fill("h_Electron_P_B", 0, p[1]);
 555       }
 556     }
 557   }
 558
 559
 560   //
 561   // Chi2/NDF cut
 562   //
 563   if(chi2ndf > cutChi2NDF) return kFALSE;
 564   fQA->Fill("h_Gamma_Mass", 1, iMass);
 565   if(iMass < cutMass){
 566     fQA->Fill("h_cut_Gamma_CosPoint", 1, cosPoint);
 567     fQA->Fill("h_Electron_P", 1, p[0]);
 568     fQA->Fill("h_Electron_P", 1, p[1]);
 569   }
 570   if(fMCEvent){
 571     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 1, iMass);
 572     else if(-2 != fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_B", 1, iMass);
 573     if(iMass < cutMass){
 574       if(1 == fCurrentV0id){
 575         fQAmc->Fill("h_cut_Gamma_CosPoint_S", 1, iP, cosPoint);
 576         fQAmc->Fill("h_Electron_P_S", 1, p[0]);
 577         fQAmc->Fill("h_Electron_P_S", 1, p[1]);
 578       }
 579       else if(-2 != fCurrentV0id){
 580         fQAmc->Fill("h_cut_Gamma_CosPoint_B", 1, iP, cosPoint);
 581         fQAmc->Fill("h_Electron_P_B", 1, p[0]);
 582         fQAmc->Fill("h_Electron_P_B", 1, p[1]);
 583       }
 584     }
 585   }
 586
 587   //
 588   // Cos point cut
 589   //
 590   if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
 591   fQA->Fill("h_Gamma_Mass", 2, iMass);
 592   if(iMass < cutMass){
 593     fQA->Fill("h_Electron_P", 2, p[0]);
 594     fQA->Fill("h_Electron_P", 2, p[1]);
 595     fQA->Fill("h_cut_Gamma_DCA", 1, dca);
 596   }
 597   if(fMCEvent){
 598     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 2, iMass);
 599     else if(-2 != fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_B", 2, iMass);
 600     if(iMass < cutMass){
 601       if(1 == fCurrentV0id){
 602         fQAmc->Fill("h_cut_Gamma_DCA_S", 1, iP, dca);
 603         fQAmc->Fill("h_Electron_P_S", 2, p[0]);
 604         fQAmc->Fill("h_Electron_P_S", 2, p[1]);
 605
 606       }
 607       else if(-2 != fCurrentV0id){
 608         fQAmc->Fill("h_cut_Gamma_DCA_B", 1, iP, dca);
 609         fQAmc->Fill("h_Electron_P_B", 2, p[0]);
 610         fQAmc->Fill("h_Electron_P_B", 2, p[1]);
 611
 612       }
 613     }
 614   }
 615
 616   //
 617   // DCA cut
 618   //
 619   if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
 620   fQA->Fill("h_Gamma_Mass", 3, iMass);
 621   if(iMass < cutMass){
 622     fQA->Fill("h_Electron_P", 3, p[0]);
 623     fQA->Fill("h_Electron_P", 3, p[1]);
 624     fQA->Fill("h_cut_Gamma_VtxR_old", 1, r);
 625     fQA->Fill("h_cut_Gamma_VtxR", 1, r2);
 626   }
 627   if(fMCEvent){
 628     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 3, iMass);
 629     else if(-2 != fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_B", 3, iMass);
 630     if(iMass < cutMass){
 631       if(1 == fCurrentV0id){
 632         fQAmc->Fill("h_cut_Gamma_VtxR_S", 1, iP, r2);
 633         fQAmc->Fill("h_Electron_P_S", 3, p[0]);
 634         fQAmc->Fill("h_Electron_P_S", 3, p[1]);
 635
 636       }
 637       else if(-2 != fCurrentV0id){
 638         fQAmc->Fill("h_cut_Gamma_VtxR_B", 1, iP, r2);
 639         fQAmc->Fill("h_Electron_P_B", 3, p[0]);
 640         fQAmc->Fill("h_Electron_P_B", 3, p[1]);
 641       }
 642     }
 643   }
 644
 645   //
 646   // Vertex radius cut
 647   //
 648   if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
 649   fQA->Fill("h_Gamma_Mass", 4, iMass);
 650   if(iMass < cutMass){
 651     fQA->Fill("h_cut_Gamma_PP", 1, psiPair);
 652     fQA->Fill("h_Electron_P", 4, p[0]);
 653     fQA->Fill("h_Electron_P", 4, p[1]);
 654   }
 655   if(fMCEvent){
 656     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 4, iMass);
 657     else if(-2 != fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_B", 4, iMass);
 658     if(iMass < cutMass){
 659       if(1 == fCurrentV0id){
 660         fQAmc->Fill("h_cut_Gamma_PP_S", 1, iP, psiPair);
 661         fQAmc->Fill("h_Electron_P_S", 4, p[0]);
 662         fQAmc->Fill("h_Electron_P_S", 4, p[1]);
 663       }
 664       else if(-2 != fCurrentV0id){
 665         fQAmc->Fill("h_cut_Gamma_PP_B", 1, iP, psiPair);
 666         fQAmc->Fill("h_Electron_P_B", 4, p[0]);
 667         fQAmc->Fill("h_Electron_P_B", 4, p[1]);
 668       }
 669     }
 670   }
 671
 672
 673   //
 674   // PsiPair cut
 675   //
 676   if(psiPair < cutPsiPair[0] || psiPair > cutPsiPair[1]) return kFALSE;
 677   fQA->Fill("h_Gamma_Mass", 5, iMass);
 678   if(iMass < cutMass){
 679     fQA->Fill("h_cut_Gamma_Sep", 1, iP, sep);
 680     fQA->Fill("h_Electron_P", 5, p[0]);
 681     fQA->Fill("h_Electron_P", 5, p[1]);
 682   }
 683   if(fMCEvent){
 684     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 5, iMass);
 685     else if(-2 != fCurrentV0id)fQAmc->Fill("h_Gamma_Mass_B", 5, iMass);
 686
 687     if(iMass < cutMass){
 688       if(1 == fCurrentV0id){
 689         fQAmc->Fill("h_cut_Gamma_Sep_S", 1, iP, sep);
 690         fQAmc->Fill("h_cut_Gamma_SepR_S", 1,r2, sep);
 691         fQAmc->Fill("h_Electron_P_S", 5, p[0]);
 692         fQAmc->Fill("h_Electron_P_S", 5, p[1]);
 693       }
 694       else if(-2 != fCurrentV0id){
 695         fQAmc->Fill("h_cut_Gamma_Sep_B", 1, iP, sep);
 696         fQAmc->Fill("h_cut_Gamma_SepR_B", 1,r2, sep);
 697         fQAmc->Fill("h_Electron_P_B", 5, p[0]);
 698         fQAmc->Fill("h_Electron_P_B", 5, p[1]);
 699       }
 700     }
 701   }
 702
 703
 704   // TESTING NEW CUT
 705   //
 706   // distance of the tracks at the entrance of the TPC
 707   //
 708   if(sep < cutSeparation) return kFALSE;
 709   fQA->Fill("h_Gamma_Mass", 6, iMass);
 710   if(iMass < cutMass){
 711     fQA->Fill("h_Electron_P", 6, p[0]);
 712     fQA->Fill("h_Electron_P", 6, p[1]);
 713   }
 714   if(fMCEvent){
 715     if(1 == fCurrentV0id) fQAmc->Fill("h_Gamma_Mass_S", 6, iMass);
 716     else if(-2 != fCurrentV0id)fQAmc->Fill("h_Gamma_Mass_B", 6, iMass);
 717
 718     if(iMass < cutMass){
 719       if(1 == fCurrentV0id){
 720         fQAmc->Fill("h_Electron_P_S", 6, p[0]);
 721         fQAmc->Fill("h_Electron_P_S", 6, p[1]);
 722       }
 723       else if(-2 != fCurrentV0id){
 724         fQAmc->Fill("h_Electron_P_B", 6, p[0]);
 725         fQAmc->Fill("h_Electron_P_B", 6, p[1]);
 726       }
 727     }
 728   }
 729
 730   // .. test
 731
 732
 733   if(iMass > cutMass) return kFALSE;
 734
 735   // all cuts passed
 736
 737
 738   // some MC stuff
 739   //printf("**D: gamma V0id: %i, P: %i, N: %i \n", fCurrentV0id, fPdaughterPDG, fNdaughterPDG);
 740   if(1 == fCurrentV0id){
 741     fQAmc->Fill("h_gamma_p_S", iP);
 742     fQAmc->Fill("h_Electron_P_S", 7, p[0]);
 743     fQAmc->Fill("h_Electron_P_S", 7, p[1]);
 744   }
 745   else if (-2 != fCurrentV0id){
 746     fQAmc->Fill("h_gamma_p_B", iP);
 747     fQAmc->Fill("h_Electron_P_B", 7, p[0]);
 748     fQAmc->Fill("h_Electron_P_B", 7, p[1]);
 749   }
 750
 751
 752   return kTRUE;
 753 }
 754 //________________________________________________________________
 755 Bool_t AliHFEV0cuts::K0Cuts(AliESDv0 *v0){
 756   //
 757   // K0 cuts
 758   //
 759
 760   if(!v0) return kFALSE;
 761
 762   if(fMCEvent){
 763     if(2 == fCurrentV0id){
 764       fQAmc->Fill("h_Mass_K0_as_G",  v0->GetEffMass(0, 0));
 765       fQAmc->Fill("h_Mass_K0_as_L",  v0->GetEffMass(2, 4));
 766       fQAmc->Fill("h_Mass_K0_as_L",  v0->GetEffMass(4, 2));
 767     }
 768   }
 769
 770   //const Double_t cK0mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();  // PDG K0s mass
 771   AliVTrack* daughter[2];
 772   Int_t pIndex = 0, nIndex = 0;
 773   if(CheckSigns(v0)){
 774     pIndex = v0->GetPindex();
 775     nIndex = v0->GetNindex();
 776   }
 777   else{
 778     pIndex = v0->GetNindex();
 779     nIndex = v0->GetPindex();
 780   }
 781
 782   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 783   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 784   if(!daughter[0] || !daughter[1]) return kFALSE;
 785
 786   AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
 787   if(!kfMother) return kFALSE;
 788   // production vertex is set in the 'CreateMotherParticle' function
 789   //kfMother->SetMassConstraint(cK0mass, 0.);
 790
 791   AliESDtrack* d[2];
 792   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
 793   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
 794
 795   Float_t iMass = v0->GetEffMass(2, 2);
 796   Float_t iP = v0->P();
 797   Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
 798   Double_t theta = v0->Theta();
 799   Double_t phi = v0->Phi();
 800   Double_t pt = v0->Pt();
 801   Double_t data[4] = {0., 0., 0., 0.};
 802
 803   // Cut values
 804   const Double_t cutChi2NDF = 2.;              // ORG [7.]
 805   const Double_t cutCosPoint[2] = {0., 0.02};  // ORG [0., 0.03]
 806   const Double_t cutDCA[2] = {0., 0.2};        // ORG [0., 0.1]
 807   const Double_t cutProdVtxR[2] = {2.0, 30.};   // ORG [0., 8.1]
 808   const Double_t cutMass[2] = {0.486, 0.508};   // ORG [0.485, 0.51]
 809   // Values
 810
 811   // cos pointing angle
 812   Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
 813   cosPoint = TMath::ACos(cosPoint);
 814
 815   // DCA between daughters
 816   Double_t dca = v0->GetDcaV0Daughters();
 817
 818   // Production vertex
 819   Double_t x, y, z;
 820   v0->GetXYZ(x,y,z);
 821
 822   Double_t r = TMath::Sqrt(x*x + y*y);
 823
 824   // V0 chi2/ndf
 825   Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
 826
 827   if(kfMother) delete kfMother;
 828
 829   //
 830   // Apply the cuts, produce QA plots (with mass cut)
 831   //
 832
 833   fQA->Fill("h_K0_Mass", 0, iMass);
 834   // MC
 835   if(fMCEvent){
 836     if(2 == fCurrentV0id){
 837       fQAmc->Fill("h_K0_Mass_S", 0, iMass);
 838       fQAmc->Fill("h_K0_MvP_S", iP, iMass);
 839     }
 840     else if(-2 != fCurrentV0id) fQAmc->Fill("h_K0_Mass_B", 0, iMass);
 841   }
 842
 843   if(iMass > cutMass[0] && iMass < cutMass[1]){
 844     fQA->Fill("h_PionK0_P", 0, p[0]);
 845     fQA->Fill("h_PionK0_P", 0, p[1]);
 846     fQA->Fill("h_cut_K0_CosPoint", 0, cosPoint);
 847     fQA->Fill("h_cut_K0_DCA", 0, dca);
 848     fQA->Fill("h_cut_K0_VtxR", 0, r);
 849     fQA->Fill("h_cut_K0_Chi2", 0, chi2ndf);
 850     fQA->Fill("h_cut_K0_Chi2", 1, chi2ndf);
 851   }
 852
 853   // MC
 854   if(fMCEvent){
 855     if(iMass > cutMass[0] && iMass < cutMass[1]){
 856       if(2 == fCurrentV0id){
 857         fQAmc->Fill("h_cut_K0_CosPoint_S", 0, iP, cosPoint);
 858         fQAmc->Fill("h_cut_K0_DCA_S", 0, iP, dca);
 859         fQAmc->Fill("h_cut_K0_VtxR_S", 0, iP, r);
 860         fQAmc->Fill("h_cut_K0_Chi2_S", 0, iP, chi2ndf);
 861         fQAmc->Fill("h_cut_K0_Chi2_S", 1, iP, chi2ndf);
 862         fQAmc->Fill("h_PionK0_P_S", 0, p[0]);
 863         fQAmc->Fill("h_PionK0_P_S", 0, p[1]);
 864        }
 865       else if(-2 != fCurrentV0id){
 866         fQAmc->Fill("h_cut_K0_CosPoint_B", 0, iP, cosPoint);
 867         fQAmc->Fill("h_cut_K0_DCA_B", 0, iP, dca);
 868         fQAmc->Fill("h_cut_K0_VtxR_B", 0, iP, r);
 869         fQAmc->Fill("h_cut_K0_Chi2_B", 0, iP, chi2ndf);
 870         fQAmc->Fill("h_cut_K0_Chi2_B", 1, iP, chi2ndf);
 871         fQAmc->Fill("h_PionK0_P_B", 0, p[0]);
 872         fQAmc->Fill("h_PionK0_P_B", 0, p[1]);
 873       }
 874     }
 875   }
 876
 877   //
 878   // Chi2/NDF cut
 879   //
 880   if(chi2ndf > cutChi2NDF) return kFALSE;
 881   fQA->Fill("h_K0_Mass", 1, iMass);
 882   if(iMass > cutMass[0] && iMass < cutMass[1]){
 883     fQA->Fill("h_cut_K0_CosPoint", 1, cosPoint);
 884     fQA->Fill("h_PionK0_P", 1, p[0]);
 885     fQA->Fill("h_PionK0_P", 1, p[1]);
 886   }
 887   if(fMCEvent){
 888     if(2 == fCurrentV0id) fQAmc->Fill("h_K0_Mass_S", 1, iMass);
 889     else if(-2 != fCurrentV0id) fQAmc->Fill("h_K0_Mass_B", 1, iMass);
 890      if(iMass > cutMass[0] && iMass < cutMass[1]){
 891        if(2 == fCurrentV0id){
 892          fQAmc->Fill("h_cut_K0_CosPoint_S", 1, iP, cosPoint);
 893          fQAmc->Fill("h_PionK0_P_S", 1, p[0]);
 894          fQAmc->Fill("h_PionK0_P_S", 1, p[1]);
 895        }
 896        else if(-2 != fCurrentV0id){
 897          fQAmc->Fill("h_cut_K0_CosPoint_B", 1, iP, cosPoint);
 898          fQAmc->Fill("h_PionK0_P_B", 1, p[0]);
 899          fQAmc->Fill("h_PionK0_P_B", 1, p[1]);
 900        }
 901      }
 902   }
 903
 904   //
 905   // Cos point cut
 906   //
 907   if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
 908   fQA->Fill("h_K0_Mass", 2, iMass);
 909   if(iMass > cutMass[0] && iMass < cutMass[1]){
 910     fQA->Fill("h_PionK0_P", 2, p[0]);
 911     fQA->Fill("h_PionK0_P", 2, p[1]);
 912     fQA->Fill("h_cut_K0_DCA", 1, dca);
 913   }
 914   if(fMCEvent){
 915     if(2 == fCurrentV0id) fQAmc->Fill("h_K0_Mass_S", 2, iMass);
 916     else if(-2 != fCurrentV0id) fQAmc->Fill("h_K0_Mass_B", 2, iMass);
 917      if(iMass > cutMass[0] && iMass < cutMass[1]){
 918        if(2 == fCurrentV0id){
 919          fQAmc->Fill("h_cut_K0_DCA_S", 1, iP, dca);
 920          fQAmc->Fill("h_PionK0_P_S", 2, p[0]);
 921          fQAmc->Fill("h_PionK0_P_S", 2, p[1]);
 922        }
 923        else if(-2 != fCurrentV0id){
 924          fQAmc->Fill("h_cut_K0_DCA_B", 1, iP, dca);
 925          fQAmc->Fill("h_PionK0_P_B", 2, p[0]);
 926          fQAmc->Fill("h_PionK0_P_B", 2, p[1]);
 927        }
 928      }
 929   }
 930
 931
 932   //
 933   // DCA cut
 934   //
 935   if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
 936   fQA->Fill("h_K0_Mass", 3, iMass);
 937   if(iMass > cutMass[0] && iMass < cutMass[1]){
 938     fQA->Fill("h_PionK0_P", 3, p[0]);
 939     fQA->Fill("h_PionK0_P", 3, p[1]);
 940     fQA->Fill("h_cut_K0_VtxR", 1, r);
 941   }
 942   if(fMCEvent){
 943     if(2 == fCurrentV0id) fQAmc->Fill("h_K0_Mass_S", 3, iMass);
 944     else if(-2 != fCurrentV0id) fQAmc->Fill("h_K0_Mass_B", 3, iMass);
 945      if(iMass > cutMass[0] && iMass < cutMass[1]){
 946        if(2 == fCurrentV0id){
 947          fQAmc->Fill("h_cut_K0_VtxR_S", 1, iP, r);
 948          fQAmc->Fill("h_PionK0_P_S", 3, p[0]);
 949          fQAmc->Fill("h_PionK0_P_S", 3, p[1]);
 950        }
 951        else if(-2 != fCurrentV0id){
 952          fQAmc->Fill("h_cut_K0_VtxR_B", 1, iP, r);
 953          fQAmc->Fill("h_PionK0_P_B", 3, p[0]);
 954          fQAmc->Fill("h_PionK0_P_B", 3, p[1]);
 955        }
 956      }
 957   }
 958
 959
 960   //
 961   // Vertex R cut
 962   //
 963   if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
 964   fQA->Fill("h_K0_Mass", 4, iMass);
 965   if(iMass > cutMass[0] && iMass < cutMass[1]){
 966     fQA->Fill("h_PionK0_P", 4, p[0]);
 967     fQA->Fill("h_PionK0_P", 4, p[1]);
 968   }
 969   if(fMCEvent){
 970     if(2 == fCurrentV0id) fQAmc->Fill("h_K0_Mass_S", 4, iMass);
 971     else if(-2 != fCurrentV0id) fQAmc->Fill("h_K0_Mass_B", 4, iMass);
 972      if(iMass > cutMass[0] && iMass < cutMass[1]){
 973        if(2 == fCurrentV0id){
 974          fQAmc->Fill("h_PionK0_P_S", 4, p[0]);
 975          fQAmc->Fill("h_PionK0_P_S", 4, p[1]);
 976        }
 977        else if(-2 != fCurrentV0id){
 978          fQAmc->Fill("h_PionK0_P_B", 4, p[0]);
 979          fQAmc->Fill("h_PionK0_P_B", 4, p[1]);
 980        }
 981      }
 982   }
 983
 984   data[0] = iMass;
 985   data[1] = pt;
 986   data[2] = theta;
 987   data[3] = phi;
 988   //printf("-D: m: %f, pT: %f, theta: %f, phi: %f\n", invMass, mPt, theta, phi);
 989   fQA->Fill("hK0", data);
 990
 991
 992   if(iMass < cutMass[0] || iMass > cutMass[1]) return kFALSE;
 993
 994   // all cuts passed
 995
 996   // some MC stuff
 997   if(2 == fCurrentV0id){
 998     fQAmc->Fill("h_K0_p_S", iP);
 999     fQAmc->Fill("h_PionK0_P_S", 5, p[0]);
1000     fQAmc->Fill("h_PionK0_P_S", 5, p[1]);
1001   }
1002   else if (-2 != fCurrentV0id){
1003     fQAmc->Fill("h_K0_p_B", iP);
1004     fQAmc->Fill("h_PionK0_P_B", 5, p[0]);
1005     fQAmc->Fill("h_PionK0_P_B", 5, p[1]);
1006   }
1007
1008   return kTRUE;
1009 }
1010 //________________________________________________________________
1011 Bool_t AliHFEV0cuts::LambdaCuts(AliESDv0 *v0, Bool_t &isLambda ){
1012   //
1013   // Lambda cuts - decision on Lambda - AntiLambda is taken too
1014   //
1015   // discrimination between lambda and antilambda - correlation of the following variables necessary:
1016   // - momentum of the proton AND momentum of the pion (proton momentum is allways larger)
1017   // - mass of the mother particle
1018
1019   if(!v0) return kFALSE;
1020
1021   if(fMCEvent){
1022     if(4 == fCurrentV0id){
1023       fQAmc->Fill("h_Mass_L_as_G",  v0->GetEffMass(0, 0));
1024       fQAmc->Fill("h_Mass_L_as_K0",  v0->GetEffMass(2, 0));
1025     }
1026   }
1027   // loose cuts first
1028   //if(LooseRejectK0(v0) || LooseRejectGamma(v0)) return kFALSE;
1029
1030   const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();  // PDG lambda mass
1031
1032   AliVTrack* daughter[2];
1033   Int_t pIndex = 0, nIndex = 0;
1034   Float_t mMass[2] = {-1., -1.};
1035   if(CheckSigns(v0)){
1036     pIndex = v0->GetPindex();
1037     nIndex = v0->GetNindex();
1038     mMass[0] = v0->GetEffMass(4, 2);
1039     mMass[1] = v0->GetEffMass(2, 4);
1040   }
1041   else{
1042     pIndex = v0->GetNindex();
1043     nIndex = v0->GetPindex();
1044     mMass[0] = v0->GetEffMass(2, 4);
1045     mMass[1] = v0->GetEffMass(4, 2);
1046   }
1047
1048   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
1049   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
1050   if(!daughter[0] || !daughter[1]) return kFALSE;
1051
1052   AliKFParticle *kfMother[2] = {0x0, 0x0};
1053   // Lambda
1054   kfMother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
1055   if(!kfMother[0]) return kFALSE;
1056
1057   // production vertex is set in the 'CreateMotherParticle' function
1058   //kfMother[0]->SetMassConstraint(cL0mass, 0.);
1059
1060   // Anti Lambda
1061   kfMother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
1062   if(!kfMother[1]) return kFALSE;
1063   // production vertex is set in the 'CreateMotherParticle' function
1064   //kfMother[1]->SetMassConstraint(cL0mass, 0.);
1065
1066   Float_t dMass[2] = {TMath::Abs(mMass[0] - cL0mass), TMath::Abs(mMass[1] - cL0mass)};
1067
1068   AliESDtrack* d[2];
1069   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
1070   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
1071   if(!d[0] || !d[1])    return kFALSE;
1072
1073   Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
1074
1075   // check the 3 lambda - antilambda variables
1076   Int_t check[2] = {-1, -1};   // 0 : lambda, 1 : antilambda
1077   // 1) momentum of the daughter particles - proton is expected to have higher momentum than pion
1078   check[0] = (p[0] > p[1]) ? 0 : 1;
1079   // 2) mass of the mother particle
1080   check[1] = (dMass[0] < dMass[1]) ? 0 : 1;
1081   fQA->Fill("h_L_checks", check[0]*1.0, check[1]*1.0);
1082
1083   // if the two check do not agree
1084   if(check[0] != check[1]){
1085     if(kfMother[0]) delete kfMother[0];
1086     if(kfMother[1]) delete kfMother[1];
1087     return kFALSE;
1088   }
1089
1090   // now that the check[0] == check[1]
1091   const Int_t type = check[0];
1092
1093   Float_t iMass =0.;
1094   if(CheckSigns(v0)){
1095     iMass = (type == 0) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
1096   }
1097   else{
1098     iMass = (type == 0) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
1099   }
1100   Float_t iP = v0->P();
1101
1102    // Cuts
1103   const Double_t cutChi2NDF = 2.;              // ORG [5.]
1104   const Double_t cutCosPoint[2] = {0., 0.01};  // ORG [0., 0.02]
1105   const Double_t cutDCA[2] = {0., 0.2};        // ORG [0., 0.2]
1106   const Double_t cutProdVtxR[2] = {2., 40.};   // ORG [0., 24.]
1107   const Double_t cutMass[2] = {1.11, 1.12};   // ORG [1.11, 1.12]
1108   // cundidate cuts
1109   // relative daughter momentum versusu mother momentum
1110
1111   // compute the cut values
1112
1113   // cos pointing angle
1114   Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
1115   cosPoint = TMath::ACos(cosPoint);
1116
1117   // DCA between daughters
1118   Double_t dca = v0->GetDcaV0Daughters();
1119
1120   // Production vertex
1121   Double_t x, y, z;
1122   v0->GetXYZ(x,y,z);
1123   Double_t r = TMath::Sqrt(x*x + y*y);
1124
1125   // proton - pion indices
1126   Int_t ix[2] = {0, 1};
1127   if(1 == type){
1128     ix[0] = 1;
1129     ix[1] = 0;
1130   }
1131
1132   // proton - pion indices - based on MC truth
1133   // for background use the reconstructed indices
1134   Int_t ixMC[2] = {-1, -1}; // {proton, pion}
1135   if(fMCEvent){
1136     if(4 == fCurrentV0id){
1137       ixMC[0] = 0;
1138       ixMC[1] = 1;
1139     }
1140     else if(-4 == fCurrentV0id){
1141       ixMC[0] = 1;
1142       ixMC[1] = 0;
1143     }
1144     else{
1145       ixMC[0] = ix[0];
1146       ixMC[1] = ix[1];
1147     }
1148   }
1149
1150   // V0 chi2/ndf
1151   Double_t chi2ndf = kfMother[type]->GetChi2()/kfMother[type]->GetNDF();
1152
1153   if(kfMother[0]) delete kfMother[0];
1154   if(kfMother[1]) delete kfMother[1];
1155
1156   // Relative daughter momentum
1157   Double_t rP = (0 == check[0]) ? p[1]/p[0] : p[0]/p[1];
1158
1159
1160   //
1161   // Apply the cuts, produce QA plots (with mass cut)
1162   //
1163
1164   (type == 0) ?   fQA->Fill("h_L_Mass", 0, iMass) :  fQA->Fill("h_AL_Mass", 0, iMass);
1165
1166
1167
1168   // MC
1169   if(fMCEvent){
1170     if(4 == fCurrentV0id){
1171       fQAmc->Fill("h_L_Mass_S", 0, iMass);
1172       fQAmc->Fill("h_lambda_MvP_S", iP, iMass);
1173     }
1174     else if(-4 == fCurrentV0id){
1175       fQAmc->Fill("h_AL_Mass_S", 0, iMass);
1176       fQAmc->Fill("h_lambda_MvP_S", iP, iMass);
1177     }
1178     else if(-2 != fCurrentV0id) fQAmc->Fill("h_LAL_Mass_B", 0, iMass);
1179   }
1180
1181
1182   if(iMass > cutMass[0] && iMass < cutMass[1]){
1183     fQA->Fill("h_ProtonL_P", 0, p[ix[0]]);
1184     fQA->Fill("h_PionL_P", 0, p[ix[1]]);
1185     fQA->Fill("h_cut_L_Chi2", 0, chi2ndf);
1186     fQA->Fill("h_cut_L_Chi2", 1, chi2ndf);
1187     fQA->Fill("h_cut_L_CosPoint", 0, cosPoint);
1188     fQA->Fill("h_cut_L_DCA", 0, dca);
1189     fQA->Fill("h_cut_L_VtxR", 0, r);
1190     fQA->Fill("h_cut_L_rdp_v_mp", 0, iP, rP);
1191   }
1192   if(fMCEvent){
1193     if(iMass > cutMass[0] && iMass < cutMass[1]){
1194       if(4 == TMath::Abs(fCurrentV0id)){
1195         fQAmc->Fill("h_cut_L_Chi2_S", 0, iP, chi2ndf);
1196         fQAmc->Fill("h_cut_L_Chi2_S", 1, iP, chi2ndf);
1197         fQAmc->Fill("h_cut_L_CosPoint_S", 0, iP, cosPoint);
1198         fQAmc->Fill("h_cut_L_DCA_S", 0, iP, dca);
1199         fQAmc->Fill("h_cut_L_VtxR_S", 0, iP, r);
1200         fQAmc->Fill("h_cut_L_rdp_v_mp_S", 0, iP, rP);
1201         fQAmc->Fill("h_ProtonL_P_S", 0, p[ixMC[0]]);
1202         fQAmc->Fill("h_PionL_P_S", 0, p[ixMC[1]]);
1203       }
1204       else if(-2 != fCurrentV0id){
1205         fQAmc->Fill("h_cut_L_Chi2_B", 0, iP, chi2ndf);
1206         fQAmc->Fill("h_cut_L_Chi2_B", 1, iP, chi2ndf);
1207         fQAmc->Fill("h_cut_L_CosPoint_B", 0, iP, cosPoint);
1208         fQAmc->Fill("h_cut_L_DCA_B", 0, iP, dca);
1209         fQAmc->Fill("h_cut_L_VtxR_B", 0, iP, r);
1210         fQAmc->Fill("h_cut_L_rdp_v_mp_B", 0, iP, rP);
1211         fQAmc->Fill("h_ProtonL_P_B", 0, p[ixMC[0]]);
1212         fQAmc->Fill("h_PionL_P_B", 0, p[ixMC[1]]);
1213       }
1214     }
1215   }
1216   //
1217   // Chi2/NDF cut
1218   //
1219   if(chi2ndf > cutChi2NDF) return kFALSE;
1220   (type == 0) ?   fQA->Fill("h_L_Mass", 1, iMass) :  fQA->Fill("h_AL_Mass", 1, iMass);
1221   if(iMass > cutMass[0] && iMass < cutMass[1]){
1222     fQA->Fill("h_cut_L_CosPoint", 1, cosPoint);
1223     fQA->Fill("h_ProtonL_P", 1, p[ix[0]]);
1224     fQA->Fill("h_PionL_P", 1, p[ix[1]]);
1225   }
1226   if(fMCEvent){
1227     if(4 == fCurrentV0id) fQAmc->Fill("h_L_Mass_S", 1, iMass);
1228     else if(-4 == fCurrentV0id)  fQAmc->Fill("h_AL_Mass_S", 1, iMass);
1229     else if(-2 != fCurrentV0id)  fQAmc->Fill("h_LAL_Mass_B", 1, iMass);
1230     if(iMass > cutMass[0] && iMass < cutMass[1]){
1231       if(4 == TMath::Abs(fCurrentV0id)){
1232         fQAmc->Fill("h_cut_L_CosPoint_S", 1, iP, cosPoint);
1233         fQAmc->Fill("h_ProtonL_P_S", 1, p[ixMC[0]]);
1234         fQAmc->Fill("h_PionL_P_S", 1, p[ixMC[1]]);
1235       }
1236       else if(-2 != fCurrentV0id){
1237         fQAmc->Fill("h_cut_L_CosPoint_B", 1, iP, cosPoint);
1238         fQAmc->Fill("h_ProtonL_P_B", 1, p[ixMC[0]]);
1239         fQAmc->Fill("h_PionL_P_B", 1, p[ixMC[1]]);
1240       }
1241     }
1242   }
1243
1244   //
1245   // Cos point cut
1246   //
1247   if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
1248   (type == 0) ?   fQA->Fill("h_L_Mass", 2, iMass) :  fQA->Fill("h_AL_Mass", 2, iMass);
1249   if(iMass > cutMass[0] && iMass < cutMass[1]){
1250     fQA->Fill("h_ProtonL_P", 2, p[ix[0]]);
1251     fQA->Fill("h_PionL_P", 2, p[ix[1]]);
1252     fQA->Fill("h_cut_L_DCA", 1, dca);
1253   }
1254   if(fMCEvent){
1255     if(4 == fCurrentV0id) fQAmc->Fill("h_L_Mass_S", 2, iMass);
1256     else if(-4 == fCurrentV0id)  fQAmc->Fill("h_AL_Mass_S", 2, iMass);
1257     else if(-2 != fCurrentV0id)  fQAmc->Fill("h_LAL_Mass_B", 2, iMass);
1258     if(iMass > cutMass[0] && iMass < cutMass[1]){
1259       if(4 == TMath::Abs(fCurrentV0id)){
1260         fQAmc->Fill("h_cut_L_DCA_S", 1, iP, dca);
1261         fQAmc->Fill("h_ProtonL_P_S", 2, p[ixMC[0]]);
1262         fQAmc->Fill("h_PionL_P_S", 2, p[ixMC[1]]);
1263       }
1264       else if(-2 != fCurrentV0id){
1265         fQAmc->Fill("h_cut_L_DCA_B", 1, iP, dca);
1266         fQAmc->Fill("h_ProtonL_P_B", 2, p[ixMC[0]]);
1267         fQAmc->Fill("h_PionL_P_B", 2, p[ixMC[1]]);
1268       }
1269     }
1270   }
1271
1272   //
1273   // DCA cut
1274   //
1275   if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
1276   (type == 0) ?   fQA->Fill("h_L_Mass", 3, iMass) :  fQA->Fill("h_AL_Mass", 3, iMass);
1277   if(iMass > cutMass[0] && iMass < cutMass[1]){
1278     fQA->Fill("h_ProtonL_P", 3, p[ix[0]]);
1279     fQA->Fill("h_PionL_P", 3, p[ix[1]]);
1280     fQA->Fill("h_cut_L_VtxR", 1, r);
1281   }
1282   if(fMCEvent){
1283     if(4 == fCurrentV0id) fQAmc->Fill("h_L_Mass_S", 3, iMass);
1284     else if(-4 == fCurrentV0id)  fQAmc->Fill("h_AL_Mass_S", 3, iMass);
1285     else if(-2 != fCurrentV0id)  fQAmc->Fill("h_LAL_Mass_B", 3, iMass);
1286     if(iMass > cutMass[0] && iMass < cutMass[1]){
1287       if(4 == TMath::Abs(fCurrentV0id)){
1288         fQAmc->Fill("h_cut_L_VtxR_S", 1, iP, r);
1289         fQAmc->Fill("h_ProtonL_P_S", 3, p[ixMC[0]]);
1290         fQAmc->Fill("h_PionL_P_S", 3, p[ixMC[1]]);
1291       }
1292       else if(-2 != fCurrentV0id){
1293         fQAmc->Fill("h_cut_L_VtxR_B", 1, iP, r);
1294         fQAmc->Fill("h_ProtonL_P_B", 3, p[ixMC[0]]);
1295         fQAmc->Fill("h_PionL_P_B", 3, p[ixMC[1]]);
1296       }
1297     }
1298   }
1299
1300   //
1301   // Vertex radius cut
1302   //
1303   if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
1304   (type == 0) ?   fQA->Fill("h_L_Mass", 4, iMass) :  fQA->Fill("h_AL_Mass", 4, iMass);
1305   if(iMass > cutMass[0] && iMass < cutMass[1]){
1306     fQA->Fill("h_ProtonL_P", 4, p[ix[0]]);
1307     fQA->Fill("h_PionL_P", 4, p[ix[1]]);
1308   }
1309   if(fMCEvent){
1310     if(4 == fCurrentV0id) fQAmc->Fill("h_L_Mass_S", 4, iMass);
1311     else if(-4 == fCurrentV0id)  fQAmc->Fill("h_AL_Mass_S", 4, iMass);
1312     else if(-2 != fCurrentV0id)  fQAmc->Fill("h_LAL_Mass_B", 4, iMass);
1313     if(iMass > cutMass[0] && iMass < cutMass[1]){
1314       if(4 == TMath::Abs(fCurrentV0id)){
1315         fQAmc->Fill("h_ProtonL_P_S", 4, p[ixMC[0]]);
1316         fQAmc->Fill("h_PionL_P_S", 4, p[ixMC[1]]);
1317       }
1318       else if(-2 != fCurrentV0id){
1319         fQAmc->Fill("h_ProtonL_P_B", 4, p[ixMC[0]]);
1320         fQAmc->Fill("h_PionL_P_B", 4, p[ixMC[1]]);
1321       }
1322     }
1323   }
1324
1325   if(iMass < cutMass[0] || iMass > cutMass[1]) {
1326     return kFALSE;
1327   }
1328
1329   // all cuts passed
1330
1331   // assign the lambda type value: Lambda: kTRUE, Anti-Lambda: kFALSE
1332   isLambda = (0 == type) ? kTRUE : kFALSE;
1333
1334   // some MC stuff
1335   if(4 == fCurrentV0id){
1336     fQAmc->Fill("h_lambda_p_S", iP);
1337   }
1338   else if(-4 == fCurrentV0id){
1339     fQAmc->Fill("h_alambda_p_S", iP);
1340   }
1341   else if (-2 != fCurrentV0id && 0 == type){
1342     fQAmc->Fill("h_lambda_p_B", iP);
1343   }
1344   else if(-2 != fCurrentV0id && 0 != type ){
1345     fQAmc->Fill("h_alambda_p_B", iP);
1346   }
1347   //
1348   if(4 == TMath::Abs(fCurrentV0id)){
1349     fQAmc->Fill("h_ProtonL_P_S", 5, p[ixMC[0]]);
1350     fQAmc->Fill("h_PionL_P_S", 5, p[ixMC[1]]);
1351   }
1352   else if(-2 != fCurrentV0id){
1353     fQAmc->Fill("h_ProtonL_P_B", 5, p[ixMC[0]]);
1354     fQAmc->Fill("h_PionL_P_B", 5, p[ixMC[1]]);
1355   }
1356
1357   return kTRUE;
1358 }
1359 //________________________________________________________________
1360 Double_t AliHFEV0cuts::OpenAngle(AliESDv0 *v0) const {
1361   //
1362   // Opening angle between two daughter tracks
1363   //
1364   Double_t mn[3] = {0,0,0};
1365   Double_t mp[3] = {0,0,0};
1366
1367
1368   v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
1369   v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
1370
1371
1372   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])));
1373
1374   return TMath::Abs(openAngle);
1375 }
1376 //________________________________________________________________
1377 Double_t AliHFEV0cuts::PsiPair(AliESDv0 *v0) {
1378   //
1379   // Angle between daughter momentum plane and plane
1380   //
1381
1382   if(!fInputEvent) return -1.;
1383
1384   Float_t magField = fInputEvent->GetMagneticField();
1385
1386   Int_t pIndex = -1;
1387   Int_t nIndex = -1;
1388   if(CheckSigns(v0)){
1389     pIndex = v0->GetPindex();
1390     nIndex = v0->GetNindex();
1391   }
1392   else{
1393     pIndex = v0->GetNindex();
1394     nIndex = v0->GetPindex();
1395   }
1396
1397
1398   AliESDtrack* daughter[2];
1399
1400   daughter[0] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(pIndex));
1401   daughter[1] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(nIndex));
1402
1403   Double_t x, y, z;
1404   v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
1405
1406   Double_t mn[3] = {0,0,0};
1407   Double_t mp[3] = {0,0,0};
1408
1409
1410   v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
1411   v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
1412
1413
1414   Double_t deltat = 1.;
1415   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
1416
1417   Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
1418
1419   Double_t momPosProp[3];
1420   Double_t momNegProp[3];
1421
1422   AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
1423
1424   Double_t psiPair = 4.;
1425
1426   if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
1427     psiPair =  -5.;
1428   if(pt.PropagateTo(radiussum,magField) == 0)
1429     psiPair = -5.;
1430   pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
1431   nt.GetPxPyPz(momNegProp);
1432
1433   Double_t pEle =
1434     TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
1435   Double_t pPos =
1436     TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
1437
1438   Double_t scalarproduct =
1439     momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
1440
1441   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
1442
1443   psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
1444
1445   return psiPair;
1446 }
1447 //________________________________________________________________
1448 AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec){
1449   //
1450   // Creates a mother particle
1451   //
1452   AliKFParticle pkfdaughter(*pdaughter, pspec);
1453   AliKFParticle nkfdaughter(*ndaughter, nspec);
1454
1455   // - check if the daughter particles are coming from the primary vertex
1456   // - check the number of tracks that take part in the creaton of primary vertex.
1457   //   important: after removeal of candidate tracks there must be at least 2 tracks left
1458   //   otherwise the primary vertex will be corrupted
1459
1460
1461   // Create the mother particle
1462   AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);
1463   // important !!!
1464   m->SetField(fInputEvent->GetMagneticField());
1465   if(TMath::Abs(kElectron) == pspec && TMath::Abs(kElectron) == nspec) m->SetMassConstraint(0, 0.001);
1466   else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass(), 0.);
1467   else if(TMath::Abs(kProton) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
1468   else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kProton) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
1469   else{
1470     AliError("Wrong daughter ID - mass constraint can not be set");
1471   }
1472
1473   AliKFVertex improvedVertex = *fPrimaryVertex;
1474   improvedVertex += *m;
1475   m->SetProductionVertex(improvedVertex);
1476
1477   // update 15/06/2010
1478   // mother particle will not be added to primary vertex but only to its copy
1479   // as this confilcts with calling
1480   // m->SetPrimaryVertex() function and
1481   // subsequently removing the mother particle afterwards
1482   // Sourse: Sergey Gorbunov
1483
1484   return m;
1485 }
1486 //_________________________________________________
1487 Bool_t AliHFEV0cuts::LooseRejectK0(AliESDv0 * const v0) const {
1488   //
1489   // Reject K0 based on loose cuts
1490   //
1491   Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
1492   if(mass > 0.494 && mass < 0.501) return kTRUE;
1493   return kFALSE;
1494 }
1495
1496 //_________________________________________________
1497 Bool_t AliHFEV0cuts::LooseRejectLambda(AliESDv0 * const v0) const {
1498   //
1499   // Reject Lambda based on loose cuts
1500   //
1501   Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
1502   Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
1503
1504   if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
1505   if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
1506   return kFALSE;
1507 }
1508
1509 //_________________________________________________
1510 Bool_t AliHFEV0cuts::LooseRejectGamma(AliESDv0 * const v0) const {
1511   //
1512   // Reject Lambda based on loose cuts
1513   //
1514
1515   Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
1516
1517   if(mass < 0.02) return kTRUE;
1518   return kFALSE;
1519 }
1520 //___________________________________________________________________
1521 void  AliHFEV0cuts::Armenteros(AliESDv0 *v0, Float_t val[2]){
1522   //
1523   // computes the Armenteros variables for given V0
1524   // fills the histogram
1525   // returns the values via "val"
1526   //
1527
1528   Double_t mn[3] = {0,0,0};
1529   Double_t mp[3] = {0,0,0};
1530   Double_t mm[3] = {0,0,0};
1531
1532   if(CheckSigns(v0)){
1533     v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1534     v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
1535   }
1536   else{
1537     v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1538     v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
1539   }
1540   v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
1541
1542   TVector3 vecN(mn[0],mn[1],mn[2]);
1543   TVector3 vecP(mp[0],mp[1],mp[2]);
1544   TVector3 vecM(mm[0],mm[1],mm[2]);
1545
1546   Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
1547   Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
1548
1549   Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
1550     ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
1551   Double_t qt = vecP.Mag()*sin(thetaP);
1552
1553   val[0] = alfa;
1554   val[1] = qt;
1555
1556 }
1557 //___________________________________________________________________
1558 Bool_t AliHFEV0cuts::CheckSigns(AliESDv0* const v0){
1559   //
1560   // check wheter the sign was correctly applied to
1561   // V0 daughter tracks
1562   //
1563
1564   Bool_t correct = kFALSE;
1565
1566   Int_t pIndex = 0, nIndex = 0;
1567   pIndex = v0->GetPindex();
1568   nIndex = v0->GetNindex();
1569
1570   AliESDtrack* d[2];
1571   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
1572   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
1573
1574   Int_t sign[2];
1575   sign[0] = (int)d[0]->GetSign();
1576   sign[1] = (int)d[1]->GetSign();
1577
1578   if(-1 == sign[0] && 1 == sign[1]){
1579     correct = kFALSE;
1580     //v0->SetIndex(0, pIndex);  // set the index of the negative v0 track
1581     //v0->SetIndex(1, nIndex);  // set the index of the positive v0 track
1582   }
1583   else{
1584     correct = kTRUE;
1585   }
1586
1587   //pIndex = v0->GetPindex();
1588   //nIndex = v0->GetNindex();
1589   //printf("-D2: P: %i, N: %i\n", pIndex, nIndex);
1590
1591   return correct;
1592 }
1593 //___________________________________________________________________
1594 Bool_t AliHFEV0cuts::GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]){
1595   //
1596   // recalculate the gamma conversion XY postition
1597   //
1598
1599   const Double_t b = fInputEvent->GetMagneticField();
1600
1601   Double_t helixcenterpos[2];
1602   GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
1603
1604   Double_t helixcenterneg[2];
1605   GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
1606
1607   Double_t  poshelix[6];
1608   ptrack->GetHelixParameters(poshelix,b);
1609   Double_t posradius = TMath::Abs(1./poshelix[4]);
1610
1611   Double_t  neghelix[6];
1612   ntrack->GetHelixParameters(neghelix,b);
1613   Double_t negradius = TMath::Abs(1./neghelix[4]);
1614
1615   Double_t xpos = helixcenterpos[0];
1616   Double_t ypos = helixcenterpos[1];
1617   Double_t xneg = helixcenterneg[0];
1618   Double_t yneg = helixcenterneg[1];
1619
1620   convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1621   convpos[1] = (ypos*negradius+  yneg*posradius)/(negradius+posradius);
1622
1623   return 1;
1624 }
1625 //___________________________________________________________________
1626 Bool_t AliHFEV0cuts::GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]){
1627   // see header file for documentation
1628
1629   Double_t pi = TMath::Pi();
1630
1631   Double_t  helix[6];
1632   track->GetHelixParameters(helix,b);
1633
1634   Double_t xpos =  helix[5];
1635   Double_t ypos =  helix[0];
1636   Double_t radius = TMath::Abs(1./helix[4]);
1637   Double_t phi = helix[2];
1638
1639   if(phi < 0){
1640     phi = phi + 2*pi;
1641   }
1642
1643   phi -= pi/2.;
1644   Double_t xpoint =  radius * TMath::Cos(phi);
1645   Double_t ypoint =  radius * TMath::Sin(phi);
1646
1647   if(b<0){
1648     if(charge > 0){
1649       xpoint = - xpoint;
1650       ypoint = - ypoint;
1651     }
1652
1653     if(charge < 0){
1654       xpoint =  xpoint;
1655       ypoint =  ypoint;
1656     }
1657   }
1658   if(b>0){
1659     if(charge > 0){
1660       xpoint =  xpoint;
1661       ypoint =  ypoint;
1662     }
1663
1664     if(charge < 0){
1665       xpoint = - xpoint;
1666       ypoint = - ypoint;
1667     }
1668   }
1669   center[0] =  xpos + xpoint;
1670   center[1] =  ypos + ypoint;
1671
1672   return 1;
1673 }