Fix compilation error on C++11 (H.Beck)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliCFVertexingHFCascade.cxx
CommitLineData
a9f238cf 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
4539d15b 19//
a9f238cf 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"
ec5288c3 33#include "AliCFTaskVertexingHF.h"
4539d15b 34#include "AliPIDResponse.h"
35#include "AliPID.h"
a9f238cf 36
37ClassImp(AliCFVertexingHFCascade)
38
39
40//_________________________________________
4539d15b 41AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
42AliCFVertexingHF(mcArray, originDselection),
43 fPDGcascade(0),
44 fPDGbachelor(0),
45 fPDGneutrDaugh(0),
46 fPDGneutrDaughPositive(0),
47 fPDGneutrDaughNegative(0),
48 fPrimVtx(0x0)
49{
a9f238cf 50 // standard constructor
51
52 SetNProngs(3);
4539d15b 53 fPtAccCut = new Float_t[fProngs];
54 fEtaAccCut = new Float_t[fProngs];
2bf2e62b 55 // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
4539d15b 56 fPtAccCut[0] = 0.;
57 fEtaAccCut[0] = 0.;
2bf2e62b 58 for(Int_t iP=1; iP<fProngs; iP++){
4539d15b 59 fPtAccCut[iP] = 0.1;
60 fEtaAccCut[iP] = 0.9;
2bf2e62b 61 }
62
a9f238cf 63}
64
65
66//_____________________________________
67AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
68{
69 // operator =
70
71 if (this != &c) {
72
73 AliCFVertexingHF::operator=(c);
4539d15b 74
a9f238cf 75 }
76 return *this;
77}
78
79//__________________________________________
80Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
81{
82 // set the AliAODRecoDecay candidate
4539d15b 83
a9f238cf 84 Bool_t bSignAssoc = kFALSE;
4539d15b 85
a9f238cf 86 fRecoCandidate = recoCand;
4539d15b 87 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
a9f238cf 88
89 if (!fRecoCandidate) {
90 AliError("fRecoCandidate not found, problem in assignement\n");
91 return bSignAssoc;
92 }
4539d15b 93
a9f238cf 94 if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
4539d15b 95
a9f238cf 96 //Int_t pdgCand = 413;
97
4539d15b 98 Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
99 Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
a9f238cf 100
101 Int_t nentries = fmcArray->GetEntriesFast();
102
103 AliDebug(3,Form("nentries = %d\n", nentries));
104
4539d15b 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);
a9f238cf 113
fbec9fa9 114 if (mcLabel == -1) return bSignAssoc;
115
1f780958 116 if (fRecoCandidate->NumberOfFakeDaughters()>0){
4539d15b 117 fFake = 0; // fake candidate
118 if (fFakeSelection == 1) return bSignAssoc;
1f780958 119 }
120 if (fRecoCandidate->NumberOfFakeDaughters()==0){
4539d15b 121 fFake = 2; // non-fake candidate
122 if (fFakeSelection == 2) return bSignAssoc;
1f780958 123 }
fbec9fa9 124
a9f238cf 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//______________________________________________
138Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
139{
140 //
141 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
142 //
a9f238cf 143
4539d15b 144 Bool_t bGenValues = kFALSE;
a9f238cf 145
4539d15b 146 Int_t daughter0cascade = fmcPartCandidate->GetDaughter(0);
147 Int_t daughter1cascade = fmcPartCandidate->GetDaughter(1);
a828d135 148
4539d15b 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;
a9f238cf 153
4539d15b 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 }
a9f238cf 165
4539d15b 166 if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
a9f238cf 167
4539d15b 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
a828d135 172
4539d15b 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);
a9f238cf 209
4539d15b 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 }
a9f238cf 233
4539d15b 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.;
a9f238cf 248
4539d15b 249 if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
250 pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
251 pTBachelor = mcPartDaughterBachelor->Pt();
252 }
a9f238cf 253
4539d15b 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;
a9f238cf 290}
291//____________________________________________
292Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
293{
294 // read the variables for the container
295
4539d15b 296 Bool_t bFillRecoValues = kFALSE;
a9f238cf 297
4539d15b 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 }
a9f238cf 306
4539d15b 307 //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
c9b41591 308 //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
a9f238cf 309
4539d15b 310 Double_t pt = cascade->Pt();
311 Double_t rapidity = cascade->Y(fPDGcascade);
312 Double_t invMass = 0.;
a9f238cf 313 Double_t cosThetaStar = 9999.;
4539d15b 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);
a9f238cf 324
325 Int_t pdgCode = fmcPartCandidate->GetPdgCode();
326
bd666732 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
4539d15b 329
a9f238cf 330 if (pdgCode > 0){
4539d15b 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);
a9f238cf 337 }
338 else {
4539d15b 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);
a9f238cf 345 }
346
4539d15b 347 Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
a9f238cf 348
4539d15b 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
a9f238cf 382 return bFillRecoValues;
383}
384
385
386//_____________________________________________________________
387Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
388{
389 // check the required decay channel
390
391 Bool_t checkCD = kFALSE;
a9f238cf 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
4539d15b 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()));
a9f238cf 408 return checkCD;
409 }
410
4539d15b 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 }
a9f238cf 423
4539d15b 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!!");
a9f238cf 429 return checkCD;
430 }
431
432 checkCD = kTRUE;
433 return checkCD;
434
435}
436
437//__________________________________________
4539d15b 438Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
a9f238cf 439{
440 //
4539d15b 441 // check wether D0 is decaing into kpi
a9f238cf 442 //
443
444 Bool_t isHadronic = kFALSE;
4539d15b 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);
a9f238cf 465
4539d15b 466 AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
467 if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
a9f238cf 468 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
469 return isHadronic;
470 }
471
4539d15b 472 Int_t numberOfExpectedDaughters = 2;
473 if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
a9f238cf 474 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
475 return isHadronic;
476 }
477
4539d15b 478 AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
479 AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
480 if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
a9f238cf 481 AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
482 return isHadronic;
483 }
484
4539d15b 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!!");
a9f238cf 490 return isHadronic;
491 }
492
4539d15b 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();
88bc191b 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
4539d15b 506 Double_t pTNeutrDaughPositive = 0;
507 Double_t pTNeutrDaughNegative = 0;
a9f238cf 508
4539d15b 509 if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
510 pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
511 pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
a9f238cf 512 }
513 else {
4539d15b 514 pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
515 pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
a9f238cf 516 }
517
518 isHadronic = kTRUE;
519
4539d15b 520 vectorNeutrDaugh[0] = pTNeutrDaughPositive;
521 vectorNeutrDaugh[1] = pTNeutrDaughNegative;
a9f238cf 522
523 return isHadronic;
524
525}
526
2bf2e62b 527//___________________________________________________________
528
529void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
530{
4539d15b 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;
2bf2e62b 542}
543
544
545
546//___________________________________________________________
547
548void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
549{
4539d15b 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}
2bf2e62b 562//___________________________________________________________
563
564void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
565{
4539d15b 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}
2bf2e62b 579
580//___________________________________________________________
581
582void AliCFVertexingHFCascade::SetAccCut()
583{
4539d15b 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...
367e9aa3 591 if(!mcPartDaughter) return;
592 Int_t mother = mcPartDaughter->GetMother();
593 AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
594 if(!mcMother) return;
595
4539d15b 596 if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
367e9aa3 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;
4539d15b 608}
a9f238cf 609
367e9aa3 610//_____________________________________________________________
611Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const
612{
4539d15b 613 //
614 // getting eta of the prong - overload the mother class method
615 //
a9f238cf 616
4539d15b 617 if (fRecoCandidate){
367e9aa3 618
4539d15b 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);
367e9aa3 628
629 return etaProng;
630
631 }
632 return 999999;
633}
634//_____________________________________________________________
635Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const
636{
4539d15b 637 //
638 // getting pt of the prong
639 //
367e9aa3 640
641 if (fRecoCandidate){
642
4539d15b 643 AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)fRecoCandidate;
367e9aa3 644 Double_t ptProng= -9999;
4539d15b 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);
367e9aa3 651
652 // Double_t ptProng = fRecoCandidate->PtProng(iProng);
653 return ptProng;
654
655 }
656 return 999999;
4539d15b 657
658}
659//_____________________________________________________________
660Bool_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
367e9aa3 703}