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