1 /**************************************************************************
2 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------------
17 // Class for HF corrections as a function of many variables and steps
18 // For D* and other cascades
20 // Author : A.Grelli a.grelli@uu.nl UTECHT
21 //-----------------------------------------------------------------------
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"
37 ClassImp(AliCFVertexingHFCascade)
40 //_________________________________________
41 AliCFVertexingHFCascade::AliCFVertexingHFCascade():
46 fPDGneutrDaughForMC(0),
47 fPDGneutrDaughPositive(0),
48 fPDGneutrDaughNegative(0),
50 fUseCutsForTMVA(kFALSE),
51 fCutOnMomConservation(0.00001)
53 // default constructor
56 fPtAccCut = new Float_t[fProngs];
57 fEtaAccCut = new Float_t[fProngs];
58 // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
61 for(Int_t iP=1; iP<fProngs; iP++){
68 //_________________________________________
69 AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
70 AliCFVertexingHF(mcArray, originDselection),
74 fPDGneutrDaughForMC(0),
75 fPDGneutrDaughPositive(0),
76 fPDGneutrDaughNegative(0),
78 fUseCutsForTMVA(kFALSE),
79 fCutOnMomConservation(0.00001)
81 // standard constructor
84 fPtAccCut = new Float_t[fProngs];
85 fEtaAccCut = new Float_t[fProngs];
86 // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
89 for(Int_t iP=1; iP<fProngs; iP++){
97 //_____________________________________
98 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
104 AliCFVertexingHF::operator=(c);
110 //__________________________________________
111 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
113 // set the AliAODRecoDecay candidate
115 Bool_t bSignAssoc = kFALSE;
117 fRecoCandidate = recoCand;
118 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
120 if (!fRecoCandidate) {
121 AliError("fRecoCandidate not found, problem in assignement\n");
125 if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
127 //Int_t pdgCand = 413;
129 Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
130 Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
132 Int_t nentries = fmcArray->GetEntriesFast();
134 AliDebug(3,Form("nentries = %d\n", nentries));
136 Bool_t isV0 = kFALSE;
137 if (fPDGcascade == 4122) {
139 pdgDgCascade[0] = fPDGbachelor;
140 pdgDgCascade[1] = fPDGneutrDaugh;
142 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));
143 Int_t mcLabel = cascade->MatchToMC(fPDGcascade, fPDGneutrDaugh, pdgDgCascade, pdgDgNeutrDaugh, fmcArray, isV0);
145 if (mcLabel == -1) return bSignAssoc;
147 if (fRecoCandidate->NumberOfFakeDaughters()>0){
148 fFake = 0; // fake candidate
149 if (fFakeSelection == 1) return bSignAssoc;
151 if (fRecoCandidate->NumberOfFakeDaughters()==0){
152 fFake = 2; // non-fake candidate
153 if (fFakeSelection == 2) return bSignAssoc;
157 fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
159 if (!fmcPartCandidate){
160 AliDebug(3,"No part candidate");
168 //______________________________________________
169 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
172 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
175 Bool_t bGenValues = kFALSE;
177 Int_t daughter0cascade = fmcPartCandidate->GetDaughter(0);
178 Int_t daughter1cascade = fmcPartCandidate->GetDaughter(1);
180 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0cascade));
181 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1cascade));
182 AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
183 AliAODMCParticle* mcPartDaughterBachelor = NULL;
185 // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
186 // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
187 // charge of the daughters to decide which is which
188 if (mcPartDaughter0->Charge()/3 == 0){
189 mcPartDaughterNeutrDaugh = mcPartDaughter0;
190 mcPartDaughterBachelor = mcPartDaughter1;
193 mcPartDaughterNeutrDaugh = mcPartDaughter1;
194 mcPartDaughterBachelor = mcPartDaughter0;
197 if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
199 Double_t vtx1[3] = {0,0,0}; // primary vertex
200 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
201 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
202 fmcPartCandidate->XvYvZv(vtx1); // cm
204 //Daughters of the neutral particle of the cascade
205 Int_t daughter0 = mcPartDaughterNeutrDaugh->GetDaughter(0); // this is the positive
206 Int_t daughter1 = mcPartDaughterNeutrDaugh->GetDaughter(1); // this is the negative
208 AliAODMCParticle* mcPartNeutrDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
209 AliAODMCParticle* mcPartNeutrDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
211 if (!mcPartNeutrDaughter0 || !mcPartNeutrDaughter1) return kFALSE;
213 // getting vertex from daughters
214 mcPartNeutrDaughter0->XvYvZv(vtx2daughter0); // cm
215 mcPartNeutrDaughter1->XvYvZv(vtx2daughter1); //cm
216 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) {
217 AliError("Daughters have different secondary vertex, skipping the track");
223 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
224 AliAODMCParticle* positiveDaugh = mcPartNeutrDaughter0;
225 AliAODMCParticle* negativeDaugh = mcPartNeutrDaughter1;
226 if (mcPartNeutrDaughter0->GetPdgCode() < 0 && mcPartNeutrDaughter1->GetPdgCode() > 0){
227 // inverting in case the positive daughter is the second one
228 positiveDaugh = mcPartNeutrDaughter1;
229 negativeDaugh = mcPartNeutrDaughter0;
232 // getting the momentum from the daughters
233 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
234 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
235 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
237 Double_t d0[2] = {0.,0.};
239 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1, vtx2daughter0, nprongs, charge, px, py, pz, d0);
241 Double_t cosThetaStar = 0.;
242 Double_t cosThetaStarNeutrDaugh = 0.;
243 Double_t cosThetaStarNeutrDaughBar = 0.;
244 cosThetaStarNeutrDaugh = decay->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
245 cosThetaStarNeutrDaughBar = decay->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughNegative, fPDGneutrDaughPositive);
246 if (mcPartDaughterNeutrDaugh->GetPdgCode() == fPDGneutrDaughForMC){ // neutral particle
247 AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
248 cosThetaStar = cosThetaStarNeutrDaugh;
250 else if (mcPartDaughterNeutrDaugh->GetPdgCode() == -fPDGneutrDaughForMC){ // neutral particle bar
251 AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
252 cosThetaStar = cosThetaStarNeutrDaughBar;
255 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()));
259 if (cosThetaStar < -1 || cosThetaStar > 1) {
260 AliWarning(Form("Invalid value for cosine Theta star %f, returning", cosThetaStar));
265 Double_t vectorNeutrDaugh[2] = {0.,0.};
267 // evaluate the correct cascade
268 if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
269 AliDebug(2, "Error! the Neutral Daughter MC doesn't have correct daughters!!");
275 Double_t cT = decay->Ct(fPDGneutrDaugh);
276 // // get the pT of the daughters
277 // Double_t pTNeutrDaugh= 0.;
278 // Double_t pTBachelor = 0.;
280 // if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
281 // pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
282 // pTBachelor = mcPartDaughterBachelor->Pt();
285 AliDebug(3, Form("The candidate has pt = %f, y = %f", fmcPartCandidate->Pt(), fmcPartCandidate->Y()));
287 switch (fConfiguration){
288 case AliCFTaskVertexingHF::kSnail:
289 vectorMC[0] = fmcPartCandidate->Pt();
290 vectorMC[1] = fmcPartCandidate->Y() ;
291 vectorMC[2] = cosThetaStar ;
292 vectorMC[3] = vectorNeutrDaugh[0];
293 vectorMC[4] = vectorNeutrDaugh[1];
294 vectorMC[5] = cT*1.E4 ; // in micron
295 vectorMC[6] = 0.; // dummy value, meaningless in MC
296 vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
297 vectorMC[8] = 1.01; // dummy value, meaningless in MC
298 vectorMC[9] = fmcPartCandidate->Phi();
299 vectorMC[10] = fzMCVertex; // z of reconstructed of primary vertex
300 vectorMC[11] = fCentValue; // reconstructed centrality
301 vectorMC[12] = 1.; // always filling with 1 at MC level
302 vectorMC[13] = 1.01; // dummy value for cosPointingXY multiplicity
303 vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
304 vectorMC[15] = fMultiplicity; // reconstructed multiplicity
306 case AliCFTaskVertexingHF::kCheetah:
307 vectorMC[0] = fmcPartCandidate->Pt();
308 vectorMC[1] = fmcPartCandidate->Y() ;
309 vectorMC[2] = cT*1.E4; // in micron
310 vectorMC[3] = fmcPartCandidate->Phi();
311 vectorMC[4] = fzMCVertex;
312 vectorMC[5] = fCentValue; // dummy value for dca, meaningless in MC
313 vectorMC[6] = 1. ; // fake: always filling with 1 at MC level
314 vectorMC[7] = fMultiplicity; // dummy value for d0pi, meaningless in MC, in micron
322 //____________________________________________
323 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
325 // read the variables for the container
327 Bool_t bFillRecoValues = kFALSE;
329 //Get the cascade and the neutral particle from it
330 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
331 AliAODRecoDecay* neutrDaugh = NULL;
332 if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
333 else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
338 //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
339 //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
341 Double_t pt = cascade->Pt();
342 Double_t rapidity = cascade->Y(fPDGcascade);
343 // Double_t invMass = 0.;
344 Double_t cosThetaStar = 9999.;
345 Double_t pTneutrDaughPos = 0.;
346 Double_t pTneutrDaughNeg = 0.;
347 Double_t dca = neutrDaugh->GetDCA();
348 // Double_t d0neutrDaughPos = 0.;
349 // Double_t d0neutrDaughNeg = 0.;
350 Double_t d0xd0 = neutrDaugh->Prodd0d0();
351 Double_t cosPointingAngle = neutrDaugh->CosPointingAngle(fPrimVtx);
352 Double_t phi = cascade->Phi();
353 Double_t cosPointingAngleXY = neutrDaugh->CosPointingAngleXY(fPrimVtx);
354 Double_t normDecayLengthXY = neutrDaugh->NormalizedDecayLengthXY(fPrimVtx);
356 Int_t pdgCode = fmcPartCandidate->GetPdgCode();
358 // UInt_t pdgDaughCascade[2] = { static_cast<UInt_t>(fPDGbachelor), static_cast<UInt_t>(fPDGneutrDaugh) }; // bachelor is first daughter of cascade
359 // UInt_t pdgDaughBarCascade[2] = { static_cast<UInt_t>(fPDGneutrDaugh), static_cast<UInt_t>(fPDGbachelor) }; // bachelor is second daughter in case of a cascade-bar
362 cosThetaStar = neutrDaugh->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
363 pTneutrDaughPos = neutrDaugh->PtProng(0);
364 pTneutrDaughNeg = neutrDaugh->PtProng(1);
365 // d0neutrDaughPos = neutrDaugh->Getd0Prong(0);
366 // d0neutrDaughNeg = neutrDaugh->Getd0Prong(1);
367 // invMass = neutrDaugh->InvMass(2, pdgDaughCascade);
370 cosThetaStar = neutrDaugh->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
371 pTneutrDaughPos = neutrDaugh->PtProng(1);
372 pTneutrDaughNeg = neutrDaugh->PtProng(0);
373 // d0neutrDaughPos = neutrDaugh->Getd0Prong(1);
374 // d0neutrDaughNeg = neutrDaugh->Getd0Prong(0);
375 // invMass = neutrDaugh->InvMass(2, pdgDaughBarCascade);
378 Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
380 switch (fConfiguration){
381 case AliCFTaskVertexingHF::kSnail:
383 vectorReco[1] = rapidity;
384 vectorReco[2] = cosThetaStar;
385 vectorReco[3] = pTneutrDaughPos;
386 vectorReco[4] = pTneutrDaughNeg;
387 vectorReco[5] = cT*1.E4; // in micron
388 vectorReco[6] = dca*1.E4; // in micron
389 vectorReco[7] = d0xd0*1.E8; // in micron^2
390 vectorReco[8] = cosPointingAngle; // in micron
392 vectorReco[10] = fzPrimVertex; // z of reconstructed of primary vertex
393 vectorReco[11] = fCentValue;
394 vectorReco[12] = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
395 vectorReco[13] = cosPointingAngleXY;
396 vectorReco[14] = normDecayLengthXY; // in cm
397 vectorReco[15] = fMultiplicity; // reconstructed multiplicity
399 case AliCFTaskVertexingHF::kCheetah:
401 vectorReco[1] = rapidity;
402 vectorReco[2] = cT*1.E4; // in micron
404 vectorReco[4] = fzPrimVertex;
405 vectorReco[5] = fCentValue;
406 vectorReco[6] = fFake;
407 vectorReco[7] = fMultiplicity;
411 bFillRecoValues = kTRUE;
413 return bFillRecoValues;
417 //_____________________________________________________________
418 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
420 // check the required decay channel
422 Bool_t checkCD = kFALSE;
424 Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
425 Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
426 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
427 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
429 if (!mcPartDaughter0 || !mcPartDaughter1) {
430 AliDebug (2,"Problems in the MC Daughters\n");
434 if(TMath::Abs(fmcPartCandidate->GetPdgCode()) == 4122 && (daughter1 - daughter0 != 1)) {
435 AliDebug(2, Form("The MC particle doesn't have the correct daughters!!"));
439 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGneutrDaughForMC &&
440 TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGbachelor) &&
441 !(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGbachelor &&
442 TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGneutrDaughForMC)) {
443 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()));
447 // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
448 AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
450 // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
451 // charge of teh daughters to decide which is which
452 AliDebug(3, Form("Charge0 = %d, Charge1 = %d", mcPartDaughter0->Charge()/3, mcPartDaughter1->Charge()/3));
453 if (mcPartDaughter0->Charge()/3 != 0){
454 mcPartDaughterNeutrDaugh = mcPartDaughter1;
457 mcPartDaughterNeutrDaugh = mcPartDaughter0;
460 Double_t vectorNeutrDaugh[2] ={0., 0.};
462 // We are looking at a cascade ...evaluate the correct cascade
463 if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
464 AliDebug(2, "Error! the NeutrDaugh MC doesn't have correct daughters!!");
473 //__________________________________________
474 Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
477 // check wether D0 is decaing into kpi
480 Bool_t isHadronic = kFALSE;
481 AliDebug(2, Form("neutralDaugh = %p, pdg = %d", neutralDaugh, neutralDaugh->GetPdgCode()));
483 if (fPDGcascade == 4122) {
484 Int_t labelresonanceDaugh = neutralDaugh->GetDaughter(0);
485 AliAODMCParticle* resonanceDaugh = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelresonanceDaugh));
486 if (!resonanceDaugh){
490 AliDebug(3, Form("The daughter of the resonant particle is a %d (we are looking for a %d)", resonanceDaugh->GetPdgCode(), fPDGneutrDaugh));
491 if (TMath::Abs(resonanceDaugh->GetPdgCode()) != fPDGneutrDaugh){
495 neutralDaugh = resonanceDaugh;
500 Int_t daughterNeutrDaugh0 = neutralDaugh->GetDaughter(0);
501 Int_t daughterNeutrDaugh1 = neutralDaugh->GetDaughter(1);
503 AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
504 if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
505 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
509 Int_t numberOfExpectedDaughters = 2;
510 if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
511 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
515 AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
516 AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
517 if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
518 AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
522 AliDebug(3, Form("Daughter 0 has pdg = %d, daughter 1 has pdg = %d", mcPartDaughterNeutrDaugh0->GetPdgCode(), mcPartDaughterNeutrDaugh1->GetPdgCode()));
523 if (!(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughPositive &&
524 TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughNegative) &&
525 !(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughNegative &&
526 TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughPositive)) {
527 AliDebug(2, "The neutral particle (MC) doesn't come from the required decay, skipping!!");
531 Double_t sumPxDau = mcPartDaughterNeutrDaugh0->Px()+mcPartDaughterNeutrDaugh1->Px();
532 Double_t sumPyDau = mcPartDaughterNeutrDaugh0->Py()+mcPartDaughterNeutrDaugh1->Py();
533 Double_t sumPzDau = mcPartDaughterNeutrDaugh0->Pz()+mcPartDaughterNeutrDaugh1->Pz();
534 Double_t pxMother = neutralDaugh->Px();
535 Double_t pyMother = neutralDaugh->Py();
536 Double_t pzMother = neutralDaugh->Pz();
537 AliDebug(3, Form("pxMother = %f, pyMother = %f, pzMother = %f", pxMother, pyMother, pzMother));
538 AliDebug(3, Form("sumPxDau = %f, sumPyDau = %f, sumPzDau = %f", sumPxDau, sumPyDau, sumPzDau));
539 if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>fCutOnMomConservation ||
540 TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>fCutOnMomConservation ||
541 TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>fCutOnMomConservation){
542 AliDebug(2, "Momentum conservation violated, skipping!!");
546 Double_t pTNeutrDaughPositive = 0;
547 Double_t pTNeutrDaughNegative = 0;
549 if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
550 pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
551 pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
554 pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
555 pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
560 vectorNeutrDaugh[0] = pTNeutrDaughPositive;
561 vectorNeutrDaugh[1] = pTNeutrDaughNegative;
567 //___________________________________________________________
569 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
572 // setting the pt cut to be used in the Acceptance steps (MC+Reco)
575 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");
577 for (Int_t iP=0; iP<fProngs; iP++){
578 fPtAccCut[iP]=ptAccCut[iP];
586 //___________________________________________________________
588 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
591 // setting the eta cut to be used in the Acceptance steps (MC+Reco)
594 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");
596 for (Int_t iP=0; iP<fProngs; iP++){
597 fEtaAccCut[iP] = etaAccCut[iP];
602 //___________________________________________________________
604 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
607 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
610 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");
612 for (Int_t iP=0; iP<fProngs; iP++){
613 fPtAccCut[iP]=ptAccCut[iP];
614 fEtaAccCut[iP]=etaAccCut[iP];
620 //___________________________________________________________
622 void AliCFVertexingHFCascade::SetAccCut()
625 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
628 Int_t bachelorPosition = 2;
629 if (fPDGcascade == 4122) bachelorPosition = 0;
630 AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[bachelorPosition])); // should be the soft pion...
631 if(!mcPartDaughter) return;
632 Int_t mother = mcPartDaughter->GetMother();
633 AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
634 if(!mcMother) return;
636 if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
637 AliError(Form("Apparently the expected bachelor is not in the third position, causing an error (pdg expected = %d, actual = %d)!!", fPDGbachelor, mcPartDaughter->GetPdgCode()));
638 AliError("This should be fixed when checking the MC part family in the CF task...");
642 for (Int_t iP=0; iP<fProngs; iP++){
647 if (fPDGcascade != 4122){
648 fPtAccCut[2]=0.06; // soft pion
649 fEtaAccCut[2]=0.9; // soft pion
655 //_____________________________________________________________
656 Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const
659 // getting eta of the prong - overload the mother class method
664 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
666 Double_t etaProng =-9999;
667 AliAODRecoDecay* neutrDaugh=0;
669 if (fPDGcascade == 413) {
670 neutrDaugh = cascade->Get2Prong();
672 else if (fPDGcascade == 4122) {
673 neutrDaugh = cascade->Getv0();
676 if (iProng==0) etaProng = neutrDaugh->EtaProng(0);
677 if (iProng==1) etaProng = neutrDaugh->EtaProng(1);
678 if (iProng==2) etaProng = cascade->EtaProng(ibachelor);
685 //_____________________________________________________________
686 Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const
689 // getting pt of the prong
694 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
695 Double_t ptProng= -9999;
696 AliAODRecoDecay* neutrDaugh=0;
697 if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
698 else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
699 if (iProng == 0) ptProng = neutrDaugh->PtProng(0);
700 if (iProng == 1) ptProng = neutrDaugh->PtProng(1);
701 if (iProng == 2) ptProng = cascade->PtProng(1);
703 // Double_t ptProng = fRecoCandidate->PtProng(iProng);
710 //_____________________________________________________________
711 Bool_t AliCFVertexingHFCascade::CheckAdditionalCuts(AliPIDResponse* pidResponse) const {
713 // function to check whether the candidate passes the additional cuts defined in the task to get the
714 // invariant mass spectra; these cuts are NOT pt-dependent
716 if (fPDGcascade == 4122){
717 // the method is implemented only in this case so far
718 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
719 AliAODv0 * v0part = cascade->Getv0();
720 AliAODTrack *bachelor = (AliAODTrack*)cascade->GetBachelor();
721 Double_t bachelorEta = bachelor->Eta();
722 AliAODTrack *v0pos = (AliAODTrack*)v0part->GetDaughter(0);
723 AliAODTrack *v0neg = (AliAODTrack*)v0part->GetDaughter(1);
724 Double_t v0posEta = v0pos->Eta();
725 Double_t v0negEta = v0neg->Eta();
727 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
728 Double_t nSigmaTPCpr=-999.;
729 Double_t nSigmaTOFpr=-999.;
730 nSigmaTPCpr = pidResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
731 nSigmaTOFpr = pidResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
732 Double_t ptArm = v0part->PtArmV0();
733 Double_t invmassK0s = v0part->MassK0Short();
734 Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
736 Bool_t cutsForTMVA = (TMath::Abs(bachelorEta) < 0.8 && TMath::Abs(v0posEta) < 0.8 && TMath::Abs(v0negEta) < 0.8) &&
737 ((nSigmaTOFpr < -800) || (TMath::Abs(nSigmaTOFpr) < 3)) &&
738 ((ptArm < 0.07) || (ptArm > 0.105)) &&
739 ((TMath::Abs(invmassK0s - mK0SPDG)) < 0.01);
742 if (!fUseCutsForTMVA) cutsForTMVA = kTRUE;
744 Bool_t cutsForInvMassTask = !(onFlyV0) &&
745 (cascade->CosV0PointingAngle()>0.99) &&
746 (TMath::Abs(nSigmaTPCpr) <= 3) &&
747 (v0part->Getd0Prong(0) < 20) &&
748 (v0part->Getd0Prong(1) < 20);
750 if (cutsForTMVA && cutsForInvMassTask) {