Fix compilation error on C++11 (H.Beck)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliCFVertexingHFCascade.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //-----------------------------------------------------------------------
17 // Class for HF corrections as a function of many variables and steps
18 // For D* and other cascades
19 //
20 // Author : A.Grelli a.grelli@uu.nl  UTECHT
21 //-----------------------------------------------------------------------
22
23 #include "AliAODRecoDecayHF2Prong.h"
24 #include "AliAODMCParticle.h"
25 #include "AliAODEvent.h"
26 #include "TClonesArray.h"
27 #include "AliCFVertexingHF.h"
28 #include "AliESDtrack.h"
29 #include "TDatabasePDG.h"
30 #include "AliAODRecoCascadeHF.h"
31 #include "AliCFVertexingHFCascade.h"
32 #include "AliCFContainer.h"
33 #include "AliCFTaskVertexingHF.h"
34 #include "AliPIDResponse.h"
35 #include "AliPID.h"
36
37 ClassImp(AliCFVertexingHFCascade)
38
39
40 //_________________________________________
41 AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
42 AliCFVertexingHF(mcArray, originDselection),
43   fPDGcascade(0),
44   fPDGbachelor(0),
45   fPDGneutrDaugh(0),
46   fPDGneutrDaughPositive(0),
47   fPDGneutrDaughNegative(0),
48   fPrimVtx(0x0)
49 {
50   // standard constructor
51
52   SetNProngs(3);
53   fPtAccCut = new Float_t[fProngs];
54   fEtaAccCut = new Float_t[fProngs];
55   // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
56   fPtAccCut[0] = 0.;
57   fEtaAccCut[0] = 0.;
58   for(Int_t iP=1; iP<fProngs; iP++){
59     fPtAccCut[iP] = 0.1;
60     fEtaAccCut[iP] = 0.9;
61   }
62
63 }
64
65
66 //_____________________________________
67 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
68 {
69   // operator =
70  
71   if  (this != &c) {
72
73     AliCFVertexingHF::operator=(c);
74
75   }
76   return *this;
77 }
78
79 //__________________________________________
80 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
81 {
82   // set the AliAODRecoDecay candidate
83
84   Bool_t bSignAssoc = kFALSE;
85
86   fRecoCandidate = recoCand;
87   AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
88
89   if (!fRecoCandidate) {
90     AliError("fRecoCandidate not found, problem in assignement\n");
91     return bSignAssoc;
92   }
93
94   if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
95
96   //Int_t pdgCand = 413;
97
98   Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
99   Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
100
101   Int_t nentries = fmcArray->GetEntriesFast();
102
103   AliDebug(3,Form("nentries = %d\n", nentries));
104  
105   Bool_t isV0 = kFALSE;
106   if (fPDGcascade == 4122) {
107     isV0 = kTRUE;
108     pdgDgCascade[0] = fPDGbachelor;
109     pdgDgCascade[1] = fPDGneutrDaugh;
110   }
111   AliDebug(3, Form("calling MatchToMC with: fPDGcascade = %d, fPDGneutrDaugh = %d, pdgDgCascade[0] = %d, pdgDgCascade[1] = %d, pdgDgNeutrDaugh[0] = %d, pdgDgNeutrDaugh[1] = %d, fmcArray = %p", fPDGcascade, fPDGneutrDaugh, pdgDgCascade[0], pdgDgCascade[1], pdgDgNeutrDaugh[0], pdgDgNeutrDaugh[1], fmcArray));
112  Int_t mcLabel = cascade->MatchToMC(fPDGcascade, fPDGneutrDaugh, pdgDgCascade, pdgDgNeutrDaugh, fmcArray, isV0); 
113   
114   if (mcLabel == -1) return bSignAssoc;
115
116   if (fRecoCandidate->NumberOfFakeDaughters()>0){
117     fFake = 0;    // fake candidate
118     if (fFakeSelection == 1) return bSignAssoc;
119   }
120   if (fRecoCandidate->NumberOfFakeDaughters()==0){
121     fFake = 2;    // non-fake candidate
122     if (fFakeSelection == 2) return bSignAssoc;
123   }
124   
125   SetMCLabel(mcLabel);
126   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel)); 
127
128   if (!fmcPartCandidate){
129     AliDebug(3,"No part candidate");
130     return bSignAssoc;
131   }
132  
133   bSignAssoc = kTRUE;
134   return bSignAssoc;
135 }
136
137 //______________________________________________
138 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) 
139 {
140   // 
141   // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
142   //
143
144   Bool_t bGenValues = kFALSE;
145
146   Int_t daughter0cascade = fmcPartCandidate->GetDaughter(0);
147   Int_t daughter1cascade = fmcPartCandidate->GetDaughter(1);
148
149   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0cascade));
150   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1cascade));
151   AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
152   AliAODMCParticle* mcPartDaughterBachelor = NULL;
153
154   // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
155   // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the 
156   // charge of the daughters to decide which is which
157   if (mcPartDaughter0->Charge()/3 == 0){
158     mcPartDaughterNeutrDaugh = mcPartDaughter0;
159     mcPartDaughterBachelor =  mcPartDaughter1;
160   }
161   else {
162     mcPartDaughterNeutrDaugh = mcPartDaughter1;
163     mcPartDaughterBachelor =  mcPartDaughter0;
164   }
165
166   if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
167
168   Double_t vtx1[3] = {0,0,0};   // primary vertex
169   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
170   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
171   fmcPartCandidate->XvYvZv(vtx1);  // cm
172
173   //Daughters of the neutral particle of the cascade
174   Int_t daughter0 = mcPartDaughterNeutrDaugh->GetDaughter(0); // this is the positive
175   Int_t daughter1 = mcPartDaughterNeutrDaugh->GetDaughter(1); // this is the negative
176
177   AliAODMCParticle* mcPartNeutrDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
178   AliAODMCParticle* mcPartNeutrDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
179
180   if (!mcPartNeutrDaughter0 || !mcPartNeutrDaughter1) return kFALSE;
181
182   // getting vertex from daughters
183   mcPartNeutrDaughter0->XvYvZv(vtx2daughter0);  // cm
184   mcPartNeutrDaughter1->XvYvZv(vtx2daughter1);  //cm
185   if (TMath::Abs(vtx2daughter0[0] - vtx2daughter1[0])>1E-5 || TMath::Abs(vtx2daughter0[1]- vtx2daughter1[1])>1E-5 || TMath::Abs(vtx2daughter0[2] - vtx2daughter1[2])>1E-5) {
186     AliError("Daughters have different secondary vertex, skipping the track");
187     return bGenValues;
188   }
189
190   Int_t nprongs = 2;
191   Short_t charge = 0;
192   // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
193   AliAODMCParticle* positiveDaugh = mcPartNeutrDaughter0;
194   AliAODMCParticle* negativeDaugh = mcPartNeutrDaughter1;
195   if (mcPartNeutrDaughter0->GetPdgCode() < 0 && mcPartNeutrDaughter1->GetPdgCode() > 0){
196     // inverting in case the positive daughter is the second one
197     positiveDaugh = mcPartNeutrDaughter1;
198     negativeDaugh = mcPartNeutrDaughter0;
199   }
200
201   // getting the momentum from the daughters
202   Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};          
203   Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};          
204   Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
205
206   Double_t d0[2] = {0.,0.};
207
208   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1, vtx2daughter0, nprongs, charge, px, py, pz, d0);
209         
210   Double_t cosThetaStar = 0.;
211   Double_t cosThetaStarNeutrDaugh = 0.;
212   Double_t cosThetaStarNeutrDaughBar = 0.;
213   cosThetaStarNeutrDaugh = decay->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
214   cosThetaStarNeutrDaughBar = decay->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughNegative, fPDGneutrDaughPositive);
215   if (mcPartDaughterNeutrDaugh->GetPdgCode() == fPDGneutrDaughForMC){  // neutral particle
216     AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
217     cosThetaStar = cosThetaStarNeutrDaugh;
218   }
219   else if (mcPartDaughterNeutrDaugh->GetPdgCode() == -fPDGneutrDaughForMC){  // neutral particle bar
220     AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
221     cosThetaStar = cosThetaStarNeutrDaughBar;
222   }
223   else{
224     AliWarning(Form("There are problems!! particle was expected to be either with pdg = %d or its antiparticle with pdg = %d, while we have a %d, check...", fPDGneutrDaughForMC, -fPDGneutrDaughForMC, mcPartDaughterNeutrDaugh->GetPdgCode()));
225     delete decay;
226     return vectorMC;
227   }
228   if (cosThetaStar < -1 || cosThetaStar > 1) {
229     AliWarning(Form("Invalid value for cosine Theta star %f, returning", cosThetaStar));
230     delete decay;
231     return bGenValues;
232   }
233         
234   Double_t vectorNeutrDaugh[2] = {0.,0.};
235
236   // evaluate the correct cascade
237   if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
238     AliDebug(2, "Error! the Neutral Daughter MC doesn't have correct daughters!!");
239     delete decay;
240     return bGenValues;  
241   }     
242
243   //ct
244   Double_t cT = decay->Ct(fPDGneutrDaugh);
245   // get the pT of the daughters
246   Double_t pTNeutrDaugh= 0.;
247   Double_t pTBachelor = 0.;
248         
249   if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
250     pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
251     pTBachelor = mcPartDaughterBachelor->Pt();
252   }
253
254   AliDebug(3, Form("The candidate has pt = %f, y = %f", fmcPartCandidate->Pt(), fmcPartCandidate->Y()));
255
256   switch (fConfiguration){
257   case AliCFTaskVertexingHF::kSnail:
258     vectorMC[0] = fmcPartCandidate->Pt();
259     vectorMC[1] = fmcPartCandidate->Y() ;
260     vectorMC[2] = cosThetaStar ;
261     vectorMC[3] = vectorNeutrDaugh[0];
262     vectorMC[4] = vectorNeutrDaugh[1];
263     vectorMC[5] = cT*1.E4 ;  // in micron
264     vectorMC[6] = 0.;   // dummy value, meaningless in MC
265     vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
266     vectorMC[8] = 1.01;    // dummy value, meaningless in MC
267     vectorMC[9] = fmcPartCandidate->Phi(); 
268     vectorMC[10] = fzMCVertex;    // z of reconstructed of primary vertex
269     vectorMC[11] = fCentValue; // reconstructed centrality
270     vectorMC[12] = 1.;           // always filling with 1 at MC level 
271     vectorMC[13] = 1.01; // dummy value for cosPointingXY  multiplicity
272     vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
273     vectorMC[15] = fMultiplicity; // reconstructed multiplicity
274     break;
275   case AliCFTaskVertexingHF::kCheetah:
276     vectorMC[0] = fmcPartCandidate->Pt();
277     vectorMC[1] = fmcPartCandidate->Y() ;
278     vectorMC[2] = cT*1.E4; // in micron
279     vectorMC[3] = fmcPartCandidate->Phi();
280     vectorMC[4] = fzMCVertex;
281     vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
282     vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level 
283     vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
284     break;
285   }
286
287   delete decay;
288   bGenValues = kTRUE;
289   return bGenValues;
290 }
291 //____________________________________________
292 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
293
294   // read the variables for the container
295
296   Bool_t bFillRecoValues = kFALSE;
297   
298   //Get the cascade and the neutral particle from it
299   AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
300   AliAODRecoDecay* neutrDaugh = NULL; 
301   if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
302   else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
303   else {
304     return kFALSE;
305   }
306
307   //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
308   //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
309
310   Double_t pt =  cascade->Pt();
311   Double_t rapidity =  cascade->Y(fPDGcascade);
312   Double_t invMass = 0.;
313   Double_t cosThetaStar = 9999.;
314   Double_t pTneutrDaughPos = 0.;
315   Double_t pTneutrDaughNeg = 0.;
316   Double_t dca = neutrDaugh->GetDCA();
317   Double_t d0neutrDaughPos = 0.;
318   Double_t d0neutrDaughNeg = 0.;
319   Double_t d0xd0 = neutrDaugh->Prodd0d0();
320   Double_t cosPointingAngle = neutrDaugh->CosPointingAngle(fPrimVtx);
321   Double_t phi = cascade->Phi();
322   Double_t cosPointingAngleXY = neutrDaugh->CosPointingAngleXY(fPrimVtx);
323   Double_t normDecayLengthXY = neutrDaugh->NormalizedDecayLengthXY(fPrimVtx);
324
325   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
326  
327   UInt_t pdgDaughCascade[2] = { static_cast<UInt_t>(fPDGbachelor),  static_cast<UInt_t>(fPDGneutrDaugh) };    // bachelor is first daughter of cascade
328   UInt_t pdgDaughBarCascade[2] = { static_cast<UInt_t>(fPDGneutrDaugh),  static_cast<UInt_t>(fPDGbachelor) }; // bachelor is second daughter in case of a cascade-bar 
329
330   if (pdgCode > 0){
331     cosThetaStar = neutrDaugh->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
332     pTneutrDaughPos = neutrDaugh->PtProng(0);
333     pTneutrDaughNeg = neutrDaugh->PtProng(1);
334     d0neutrDaughPos = neutrDaugh->Getd0Prong(0);
335     d0neutrDaughNeg = neutrDaugh->Getd0Prong(1);
336     invMass = neutrDaugh->InvMass(2, pdgDaughCascade);
337   }
338   else {
339     cosThetaStar = neutrDaugh->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
340     pTneutrDaughPos = neutrDaugh->PtProng(1);
341     pTneutrDaughNeg = neutrDaugh->PtProng(0);
342     d0neutrDaughPos = neutrDaugh->Getd0Prong(1);
343     d0neutrDaughNeg = neutrDaugh->Getd0Prong(0);
344     invMass = neutrDaugh->InvMass(2, pdgDaughBarCascade);
345   }
346   
347   Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
348   
349   switch (fConfiguration){
350   case AliCFTaskVertexingHF::kSnail:
351     vectorReco[0] = pt;
352     vectorReco[1] = rapidity;
353     vectorReco[2] = cosThetaStar;
354     vectorReco[3] = pTneutrDaughPos;
355     vectorReco[4] = pTneutrDaughNeg;
356     vectorReco[5] = cT*1.E4;  // in micron
357     vectorReco[6] = dca*1.E4;  // in micron
358     vectorReco[7] = d0xd0*1.E8;  // in micron^2
359     vectorReco[8] = cosPointingAngle;  // in micron
360     vectorReco[9] = phi;
361     vectorReco[10] = fzPrimVertex;    // z of reconstructed of primary vertex
362     vectorReco[11] = fCentValue;
363     vectorReco[12] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
364     vectorReco[13] = cosPointingAngleXY;
365     vectorReco[14] = normDecayLengthXY; // in cm
366     vectorReco[15] = fMultiplicity; // reconstructed multiplicity
367     break;
368   case AliCFTaskVertexingHF::kCheetah:
369     vectorReco[0] = pt;
370     vectorReco[1] = rapidity;
371     vectorReco[2] = cT*1.E4; // in micron
372     vectorReco[3] = phi;
373     vectorReco[4] = fzPrimVertex;
374     vectorReco[5] = fCentValue;
375     vectorReco[6] = fFake;
376     vectorReco[7] = fMultiplicity; 
377     break;
378   }
379
380   bFillRecoValues = kTRUE;
381
382   return bFillRecoValues;
383 }
384
385
386 //_____________________________________________________________
387 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
388
389   // check the required decay channel
390
391   Bool_t checkCD = kFALSE;
392   
393   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
394   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
395   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
396   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
397
398   if (!mcPartDaughter0 || !mcPartDaughter1) {
399     AliDebug (2,"Problems in the MC Daughters\n");
400     return checkCD;
401   }
402
403   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGneutrDaughForMC &&
404         TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGbachelor) && 
405       !(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGbachelor &&
406         TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGneutrDaughForMC)) {
407     AliDebug(2, Form("The cascade MC doesn't come from a the decay under study, skipping!! (Pdg codes of daughters = %d, %d)", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
408     return checkCD;  
409   }
410   
411   // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
412   AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
413
414   // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the 
415   // charge of teh daughters to decide which is which
416   AliDebug(3, Form("Charge0 = %d, Charge1 = %d", mcPartDaughter0->Charge()/3, mcPartDaughter1->Charge()/3));
417   if (mcPartDaughter0->Charge()/3 != 0){
418     mcPartDaughterNeutrDaugh = mcPartDaughter1;
419   }
420   else {
421     mcPartDaughterNeutrDaugh = mcPartDaughter0;
422   }
423
424   Double_t vectorNeutrDaugh[2] ={0., 0.};
425
426   // We are looking at a cascade ...evaluate the correct cascade
427   if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
428     AliDebug(2, "Error! the NeutrDaugh MC doesn't have correct daughters!!");
429     return checkCD;  
430   }
431    
432   checkCD = kTRUE;
433   return checkCD;
434   
435 }
436
437 //__________________________________________
438 Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
439 {  
440   //
441   // check wether D0 is decaing into kpi
442   //
443   
444   Bool_t isHadronic = kFALSE;
445   AliDebug(2, Form("neutralDaugh = %p, pdg = %d", neutralDaugh, neutralDaugh->GetPdgCode()));
446
447   if (fPDGcascade == 4122) {
448     Int_t labelresonanceDaugh = neutralDaugh->GetDaughter(0);
449     AliAODMCParticle* resonanceDaugh = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelresonanceDaugh));
450     if (!resonanceDaugh){
451       return kFALSE;
452     }
453     else {
454       if (TMath::Abs(resonanceDaugh->GetPdgCode()) != fPDGneutrDaugh){
455         return kFALSE;
456       }
457       else {
458         neutralDaugh = resonanceDaugh;
459       }
460     }
461   }
462
463   Int_t daughterNeutrDaugh0 = neutralDaugh->GetDaughter(0);
464   Int_t daughterNeutrDaugh1 = neutralDaugh->GetDaughter(1);
465   
466   AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
467   if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
468     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
469     return isHadronic;  
470   }
471   
472   Int_t numberOfExpectedDaughters = 2;
473   if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
474     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
475     return isHadronic;  
476   }
477   
478   AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
479   AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
480   if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
481     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping"); 
482     return isHadronic;  
483   }
484   
485   if (!(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughPositive &&
486         TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughNegative) && 
487       !(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughNegative &&
488         TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughPositive)) {
489     AliDebug(2, "The neutral particle (MC) doesn't come from the required decay, skipping!!");
490     return isHadronic;  
491   }
492   
493   Double_t sumPxDau = mcPartDaughterNeutrDaugh0->Px()+mcPartDaughterNeutrDaugh1->Px();
494   Double_t sumPyDau = mcPartDaughterNeutrDaugh0->Py()+mcPartDaughterNeutrDaugh1->Py();
495   Double_t sumPzDau = mcPartDaughterNeutrDaugh0->Pz()+mcPartDaughterNeutrDaugh1->Pz();
496   Double_t pxMother = neutralDaugh->Px();
497   Double_t pyMother = neutralDaugh->Py();
498   Double_t pzMother = neutralDaugh->Pz();
499   if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>0.00001 ||
500      TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>0.00001 ||
501      TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>0.00001){
502     AliDebug(2, "Momentum conservation violated, skipping!!");
503     return isHadronic;  
504   }
505
506   Double_t pTNeutrDaughPositive = 0;
507   Double_t pTNeutrDaughNegative = 0;
508   
509   if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
510     pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
511     pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
512   }
513   else {
514     pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
515     pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
516   }
517   
518   isHadronic = kTRUE;
519   
520   vectorNeutrDaugh[0] = pTNeutrDaughPositive;
521   vectorNeutrDaugh[1] = pTNeutrDaughNegative;
522  
523   return isHadronic;
524
525 }
526
527 //___________________________________________________________
528
529 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
530 {
531   //
532   // setting the pt cut to be used in the Acceptance steps (MC+Reco)
533   //
534
535   AliDebug(3, "The 3rd element of the pt cut array will correspond to the cut applied to the soft pion - please check that it is correct");
536   if (fProngs>0){
537     for (Int_t iP=0; iP<fProngs; iP++){
538       fPtAccCut[iP]=ptAccCut[iP];
539     }
540   }
541   return;
542 }               
543
544
545
546 //___________________________________________________________
547
548 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
549 {
550   //
551   // setting the eta cut to be used in the Acceptance steps (MC+Reco)
552   //
553
554   AliDebug(3, "The 3rd element of the eta cut array will correspond to the cut applied to the soft pion - please check that it is correct");
555   if (fProngs>0){
556     for (Int_t iP=0; iP<fProngs; iP++){
557       fEtaAccCut[iP] = etaAccCut[iP];
558     }
559   }
560   return;
561 }
562 //___________________________________________________________
563
564 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
565 {
566   //
567   // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
568   //
569
570   AliDebug(3, "The 3rd element of the pt and cut array will correspond to the cut applied to the soft pion - please check that they are correct");
571   if (fProngs>0){
572     for (Int_t iP=0; iP<fProngs; iP++){
573       fPtAccCut[iP]=ptAccCut[iP];
574       fEtaAccCut[iP]=etaAccCut[iP];
575     }
576   }
577   return;
578 }
579
580 //___________________________________________________________
581
582 void AliCFVertexingHFCascade::SetAccCut()
583 {
584   //
585   // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
586   //
587
588   Int_t bachelorPosition = 2;
589   if (fPDGcascade == 4122) bachelorPosition = 0;
590   AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[bachelorPosition]));  // should be the soft pion...  
591   if(!mcPartDaughter) return;
592   Int_t mother =  mcPartDaughter->GetMother();
593   AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); 
594   if(!mcMother) return;
595
596   if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
597     AliFatal("Apparently the soft pion is not in the third position, causing a crash!!");
598   }              
599   if (fProngs>0){
600     for (Int_t iP=0; iP<fProngs-1; iP++){
601       fPtAccCut[iP]=0.1;
602       fEtaAccCut[iP]=0.9;
603     }
604     fPtAccCut[2]=0.06;  // soft pion
605     fEtaAccCut[2]=0.9;  // soft pion
606   }
607   return;
608 }
609
610 //_____________________________________________________________
611 Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const 
612 {
613   //
614   // getting eta of the prong - overload the mother class method
615   //
616
617   if (fRecoCandidate){
618    
619     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
620
621     Double_t etaProng =-9999;
622     AliAODRecoDecay* neutrDaugh; 
623     if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
624     else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
625     if (iProng==0) etaProng = neutrDaugh->EtaProng(0);
626     if (iProng==1) etaProng = neutrDaugh->EtaProng(1);
627     if (iProng==2) etaProng = cascade->EtaProng(1);
628     
629     return etaProng;
630     
631   }
632   return 999999;    
633 }
634 //_____________________________________________________________
635 Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const 
636 {
637   //
638   // getting pt of the prong
639   //
640   
641   if (fRecoCandidate){
642
643     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
644     Double_t ptProng= -9999;
645     AliAODRecoDecay* neutrDaugh; 
646     if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
647     else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
648     if (iProng == 0) ptProng = neutrDaugh->PtProng(0);
649     if (iProng == 1) ptProng = neutrDaugh->PtProng(1);
650     if (iProng == 2) ptProng = cascade->PtProng(1);
651     
652     //  Double_t ptProng = fRecoCandidate->PtProng(iProng);  
653     return ptProng;
654     
655   }
656   return 999999;  
657
658 }
659 //_____________________________________________________________
660 Bool_t AliCFVertexingHFCascade::CheckAdditionalCuts(AliPIDResponse* pidResponse) const {
661
662   // function to check whether the candidate passes the additional cuts defined in the task to get the
663   // invariant mass spectra; these cuts are NOT pt-dependent
664
665   if (fPDGcascade == 4122){
666     // the method is implemented only in this case so far
667     AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
668     AliAODv0 * v0part = cascade->Getv0();
669     AliAODTrack *bachelor = (AliAODTrack*)cascade->GetBachelor();
670
671     Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
672     Double_t nSigmaTPCpr=-999.;
673     Double_t nSigmaTOFpr=-999.;
674     nSigmaTPCpr = pidResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
675     nSigmaTOFpr = pidResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
676     Double_t ptArm = v0part->PtArmV0();
677     Double_t invmassK0s = v0part->MassK0Short();
678     Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
679
680     Bool_t cutsForTMVA = ((nSigmaTOFpr > -800) && (nSigmaTOFpr < 3)) &&  
681       ((ptArm > 0.07) && (ptArm < 0.105)) &&
682       ((TMath::Abs(invmassK0s - mK0SPDG)) < 0.01);
683     
684     if (!fUseCutsForTMVA) cutsForTMVA = kTRUE;
685
686     Bool_t cutsForInvMassTask = !(onFlyV0) && 
687       (cascade->CosV0PointingAngle()>0.99) && 
688       (TMath::Abs(nSigmaTPCpr) <= 3) && 
689       (v0part->Getd0Prong(0) < 20) && 
690       (v0part->Getd0Prong(1) < 20);
691
692     if (cutsForTMVA && cutsForInvMassTask) {    
693       // K0Smass cut
694       // eta cut
695       // TOF PID cut
696       // Arm cut
697       return kTRUE;
698     }
699   }
700
701   return kFALSE;
702
703 }