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