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