AliPhysics  608b256 (608b256)
AliCFVertexingHFCascade.cxx
Go to the documentation of this file.
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 
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"
32 #include "AliCFContainer.h"
33 #include "AliCFTaskVertexingHF.h"
34 #include "AliPIDResponse.h"
35 #include "AliPID.h"
36 
38 ClassImp(AliCFVertexingHFCascade);
40 
41 //_________________________________________
44  fPDGcascade(0),
45  fPDGbachelor(0),
46  fPDGneutrDaugh(0),
47  fPDGneutrDaughForMC(0),
48  fPDGneutrDaughPositive(0),
49  fPDGneutrDaughNegative(0),
50  fPrimVtx(0x0),
51  fUseCutsForTMVA(kFALSE),
52  fCutOnMomConservation(0.00001)
53 {
54  // default constructor
55 
56  SetNProngs(3);
57  fPtAccCut = new Float_t[fProngs];
58  fEtaAccCut = new Float_t[fProngs];
59  // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
60  fPtAccCut[0] = 0.;
61  fEtaAccCut[0] = 0.;
62  for(Int_t iP=1; iP<fProngs; iP++){
63  fPtAccCut[iP] = 0.1;
64  fEtaAccCut[iP] = 0.9;
65  }
66 
67 }
68 
69 //_________________________________________
70 AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
71 AliCFVertexingHF(mcArray, originDselection),
72  fPDGcascade(0),
73  fPDGbachelor(0),
74  fPDGneutrDaugh(0),
78  fPrimVtx(0x0),
79  fUseCutsForTMVA(kFALSE),
80  fCutOnMomConservation(0.00001)
81 {
82  // standard constructor
83 
84  SetNProngs(3);
85  fPtAccCut = new Float_t[fProngs];
86  fEtaAccCut = new Float_t[fProngs];
87  // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
88  fPtAccCut[0] = 0.;
89  fEtaAccCut[0] = 0.;
90  for(Int_t iP=1; iP<fProngs; iP++){
91  fPtAccCut[iP] = 0.1;
92  fEtaAccCut[iP] = 0.9;
93  }
94 
95 }
96 
97 
98 //_____________________________________
100 {
101  // operator =
102 
103  if (this != &c) {
104 
106 
107  }
108  return *this;
109 }
110 
111 //__________________________________________
113 {
114  // set the AliAODRecoDecay candidate
115 
116  Bool_t bSignAssoc = kFALSE;
117 
118  fRecoCandidate = recoCand;
119  AliAODRecoCascadeHF* cascade = (AliAODRecoCascadeHF*)recoCand;
120 
121  if (!fRecoCandidate) {
122  AliError("fRecoCandidate not found, problem in assignement\n");
123  return bSignAssoc;
124  }
125 
126  if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
127 
128  //Int_t pdgCand = 413;
129 
130  Int_t pdgDgCascade[2] = {fPDGneutrDaugh, fPDGbachelor};
131  Int_t pdgDgNeutrDaugh[2] = {fPDGneutrDaughPositive, fPDGneutrDaughNegative};
132 
133  Int_t nentries = fmcArray->GetEntriesFast();
134 
135  AliDebug(3,Form("nentries = %d\n", nentries));
136 
137  Bool_t isV0 = kFALSE;
138  if (fPDGcascade == 4122) {
139  isV0 = kTRUE;
140  pdgDgCascade[0] = fPDGbachelor;
141  pdgDgCascade[1] = fPDGneutrDaugh;
142  }
143  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));
144  Int_t mcLabel = cascade->MatchToMC(fPDGcascade, fPDGneutrDaugh, pdgDgCascade, pdgDgNeutrDaugh, fmcArray, isV0);
145 
146  if (mcLabel == -1) return bSignAssoc;
147 
149  fFake = 0; // fake candidate
150  if (fFakeSelection == 1) return bSignAssoc;
151  }
153  fFake = 2; // non-fake candidate
154  if (fFakeSelection == 2) return bSignAssoc;
155  }
156 
157  SetMCLabel(mcLabel);
158  fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
159 
160  if (!fmcPartCandidate){
161  AliDebug(3,"No part candidate");
162  return bSignAssoc;
163  }
164 
165  bSignAssoc = kTRUE;
166  return bSignAssoc;
167 }
168 
169 //______________________________________________
171 {
172  //
173  // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
174  //
175 
176  Bool_t bGenValues = kFALSE;
177 
178  Int_t daughter0cascade = fmcPartCandidate->GetDaughterLabel(0);
179  Int_t daughter1cascade = fmcPartCandidate->GetDaughterLabel(1);
180 
181  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0cascade));
182  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1cascade));
183  AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
184  AliAODMCParticle* mcPartDaughterBachelor = NULL;
185 
186  // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
187  // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
188  // charge of the daughters to decide which is which
189  if (mcPartDaughter0->Charge()/3 == 0){
190  mcPartDaughterNeutrDaugh = mcPartDaughter0;
191  mcPartDaughterBachelor = mcPartDaughter1;
192  }
193  else {
194  mcPartDaughterNeutrDaugh = mcPartDaughter1;
195  mcPartDaughterBachelor = mcPartDaughter0;
196  }
197 
198  if (!mcPartDaughterNeutrDaugh || !mcPartDaughterBachelor) return kFALSE;
199 
200  Double_t vtx1[3] = {0,0,0}; // primary vertex
201  Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
202  Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
203  fmcPartCandidate->XvYvZv(vtx1); // cm
204 
205  //Daughters of the neutral particle of the cascade
206  Int_t daughter0 = mcPartDaughterNeutrDaugh->GetDaughterLabel(0); // this is the positive
207  Int_t daughter1 = mcPartDaughterNeutrDaugh->GetDaughterLabel(1); // this is the negative
208 
209  AliAODMCParticle* mcPartNeutrDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
210  AliAODMCParticle* mcPartNeutrDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
211 
212  if (!mcPartNeutrDaughter0 || !mcPartNeutrDaughter1) return kFALSE;
213 
214  // getting vertex from daughters
215  mcPartNeutrDaughter0->XvYvZv(vtx2daughter0); // cm
216  mcPartNeutrDaughter1->XvYvZv(vtx2daughter1); //cm
217  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) {
218  AliError("Daughters have different secondary vertex, skipping the track");
219  return bGenValues;
220  }
221 
222  Int_t nprongs = 2;
223  Short_t charge = 0;
224  // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
225  AliAODMCParticle* positiveDaugh = mcPartNeutrDaughter0;
226  AliAODMCParticle* negativeDaugh = mcPartNeutrDaughter1;
227  if (mcPartNeutrDaughter0->GetPdgCode() < 0 && mcPartNeutrDaughter1->GetPdgCode() > 0){
228  // inverting in case the positive daughter is the second one
229  positiveDaugh = mcPartNeutrDaughter1;
230  negativeDaugh = mcPartNeutrDaughter0;
231  }
232 
233  // getting the momentum from the daughters
234  Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
235  Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
236  Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
237 
238  Double_t d0[2] = {0.,0.};
239 
240  AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1, vtx2daughter0, nprongs, charge, px, py, pz, d0);
241 
242  Double_t cosThetaStar = 0.;
243  Double_t cosThetaStarNeutrDaugh = 0.;
244  Double_t cosThetaStarNeutrDaughBar = 0.;
245  cosThetaStarNeutrDaugh = decay->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
246  cosThetaStarNeutrDaughBar = decay->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughNegative, fPDGneutrDaughPositive);
247  if (mcPartDaughterNeutrDaugh->GetPdgCode() == fPDGneutrDaughForMC){ // neutral particle
248  AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d", mcPartDaughter0->GetPdgCode(), mcPartDaughter1->GetPdgCode()));
249  cosThetaStar = cosThetaStarNeutrDaugh;
250  }
251  else if (mcPartDaughterNeutrDaugh->GetPdgCode() == -fPDGneutrDaughForMC){ // neutral particle bar
252  AliDebug(3, Form("Neutral Daughter, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
253  cosThetaStar = cosThetaStarNeutrDaughBar;
254  }
255  else{
256  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()));
257  delete decay;
258  return vectorMC;
259  }
260  if (cosThetaStar < -1 || cosThetaStar > 1) {
261  AliWarning(Form("Invalid value for cosine Theta star %f, returning", cosThetaStar));
262  delete decay;
263  return bGenValues;
264  }
265 
266  Double_t vectorNeutrDaugh[2] = {0.,0.};
267 
268  // evaluate the correct cascade
269  if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
270  AliDebug(2, "Error! the Neutral Daughter MC doesn't have correct daughters!!");
271  delete decay;
272  return bGenValues;
273  }
274 
275  //ct
276  Double_t cT = decay->Ct(fPDGneutrDaugh);
277  // // get the pT of the daughters
278  // Double_t pTNeutrDaugh= 0.;
279  // Double_t pTBachelor = 0.;
280 
281  // if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == fPDGcascade) {
282  // pTNeutrDaugh = mcPartDaughterNeutrDaugh->Pt();
283  // pTBachelor = mcPartDaughterBachelor->Pt();
284  // }
285 
286  AliDebug(3, Form("The candidate has pt = %f, y = %f", fmcPartCandidate->Pt(), fmcPartCandidate->Y()));
287 
288  Int_t localmult = -1;
290  localmult = ComputeLocalMultiplicity(decay->Eta(), decay->Phi(), 0.4);
291  }
292 
293  switch (fConfiguration){
295  vectorMC[0] = fmcPartCandidate->Pt();
296  vectorMC[1] = fmcPartCandidate->Y() ;
297  vectorMC[2] = cosThetaStar ;
298  vectorMC[3] = vectorNeutrDaugh[0];
299  vectorMC[4] = vectorNeutrDaugh[1];
300  vectorMC[5] = cT*1.E4 ; // in micron
301  vectorMC[6] = 0.; // dummy value, meaningless in MC
302  vectorMC[7] = -100000.; // dummy value, meaningless in MC, in micron^2
303  vectorMC[8] = 1.01; // dummy value, meaningless in MC
304  vectorMC[9] = fmcPartCandidate->Phi();
305  vectorMC[10] = fzMCVertex; // z of reconstructed of primary vertex
306  vectorMC[11] = fCentValue; // reconstructed centrality
307  vectorMC[12] = 1.; // always filling with 1 at MC level
308  vectorMC[13] = 1.01; // dummy value for cosPointingXY multiplicity
309  vectorMC[14] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
310  vectorMC[15] = fMultiplicity; // reconstructed multiplicity
311  break;
313  vectorMC[0] = fmcPartCandidate->Pt();
314  vectorMC[1] = fmcPartCandidate->Y() ;
315  vectorMC[2] = cT*1.E4; // in micron
316  vectorMC[3] = fmcPartCandidate->Phi();
317  vectorMC[4] = fzMCVertex;
318  vectorMC[5] = fCentValue; // dummy value for dca, meaningless in MC
319  vectorMC[6] = 1. ; // fake: always filling with 1 at MC level
320  vectorMC[7] = fMultiplicity; // dummy value for d0pi, meaningless in MC, in micron
321  break;
323  vectorMC[0] = fmcPartCandidate->Pt();
324  vectorMC[1] = fmcPartCandidate->Y() ;
325  vectorMC[2] = fCentValue; // dummy value for dca, meaningless in MC
326  vectorMC[3] = fMultiplicity; // dummy value for d0pi, meaningless in MC, in micron
327  break;
329  vectorMC[0] = fmcPartCandidate->Pt();
330  vectorMC[1] = fmcPartCandidate->Y() ;
331  vectorMC[2] = fCentValue; // centrality
332  vectorMC[3] = fMultiplicity; // multiplicity (diff estimators can be used)
333  vectorMC[4] = localmult; // local multiplicity (Ntracks in R<0.4)
334  vectorMC[5] = fq2; // magnitude of reduced flow vector (computed using TPC tracks)
335  break;
336  }
337 
338  delete decay;
339  bGenValues = kTRUE;
340  return bGenValues;
341 }
342 //____________________________________________
344 {
345  // read the variables for the container
346 
347  Bool_t bFillRecoValues = kFALSE;
348 
349  //Get the cascade and the neutral particle from it
351  AliAODRecoDecay* neutrDaugh = NULL;
352  if (fPDGcascade == 413) neutrDaugh = cascade->Get2Prong();
353  else if (fPDGcascade == 4122) neutrDaugh = cascade->Getv0();
354  else {
355  return kFALSE;
356  }
357 
358  //if (cascade->GetPrimaryVtx())printf("cascade has primary vtx\n");
359  //if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
360 
361  Double_t pt = cascade->Pt();
362  Double_t rapidity = cascade->Y(fPDGcascade);
363  // Double_t invMass = 0.;
364  Double_t cosThetaStar = 9999.;
365  Double_t pTneutrDaughPos = 0.;
366  Double_t pTneutrDaughNeg = 0.;
367  Double_t dca = neutrDaugh->GetDCA();
368  // Double_t d0neutrDaughPos = 0.;
369  // Double_t d0neutrDaughNeg = 0.;
370  Double_t d0xd0 = neutrDaugh->Prodd0d0();
371  Double_t cosPointingAngle = neutrDaugh->CosPointingAngle(fPrimVtx);
372  Double_t phi = cascade->Phi();
373  Double_t cosPointingAngleXY = neutrDaugh->CosPointingAngleXY(fPrimVtx);
374  Double_t normDecayLengthXY = neutrDaugh->NormalizedDecayLengthXY(fPrimVtx);
375 
376  Int_t pdgCode = fmcPartCandidate->GetPdgCode();
377 
378  // UInt_t pdgDaughCascade[2] = { static_cast<UInt_t>(fPDGbachelor), static_cast<UInt_t>(fPDGneutrDaugh) }; // bachelor is first daughter of cascade
379  // UInt_t pdgDaughBarCascade[2] = { static_cast<UInt_t>(fPDGneutrDaugh), static_cast<UInt_t>(fPDGbachelor) }; // bachelor is second daughter in case of a cascade-bar
380 
381  if (pdgCode > 0){
382  cosThetaStar = neutrDaugh->CosThetaStar(1, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
383  pTneutrDaughPos = neutrDaugh->PtProng(0);
384  pTneutrDaughNeg = neutrDaugh->PtProng(1);
385  // d0neutrDaughPos = neutrDaugh->Getd0Prong(0);
386  // d0neutrDaughNeg = neutrDaugh->Getd0Prong(1);
387  // invMass = neutrDaugh->InvMass(2, pdgDaughCascade);
388  }
389  else {
390  cosThetaStar = neutrDaugh->CosThetaStar(0, fPDGneutrDaugh, fPDGneutrDaughPositive, fPDGneutrDaughNegative);
391  pTneutrDaughPos = neutrDaugh->PtProng(1);
392  pTneutrDaughNeg = neutrDaugh->PtProng(0);
393  // d0neutrDaughPos = neutrDaugh->Getd0Prong(1);
394  // d0neutrDaughNeg = neutrDaugh->Getd0Prong(0);
395  // invMass = neutrDaugh->InvMass(2, pdgDaughBarCascade);
396  }
397 
398  Double_t cT = neutrDaugh->Ct(fPDGneutrDaugh, fPrimVtx);
399 
400  Int_t localmult = -1;
402  localmult = ComputeLocalMultiplicity(cascade->Eta(), cascade->Phi(), 0.4);
403  }
404 
405  switch (fConfiguration){
407  vectorReco[0] = pt;
408  vectorReco[1] = rapidity;
409  vectorReco[2] = cosThetaStar;
410  vectorReco[3] = pTneutrDaughPos;
411  vectorReco[4] = pTneutrDaughNeg;
412  vectorReco[5] = cT*1.E4; // in micron
413  vectorReco[6] = dca*1.E4; // in micron
414  vectorReco[7] = d0xd0*1.E8; // in micron^2
415  vectorReco[8] = cosPointingAngle; // in micron
416  vectorReco[9] = phi;
417  vectorReco[10] = fzPrimVertex; // z of reconstructed of primary vertex
418  vectorReco[11] = fCentValue;
419  vectorReco[12] = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
420  vectorReco[13] = cosPointingAngleXY;
421  vectorReco[14] = normDecayLengthXY; // in cm
422  vectorReco[15] = fMultiplicity; // reconstructed multiplicity
423  break;
425  vectorReco[0] = pt;
426  vectorReco[1] = rapidity;
427  vectorReco[2] = cT*1.E4; // in micron
428  vectorReco[3] = phi;
429  vectorReco[4] = fzPrimVertex;
430  vectorReco[5] = fCentValue;
431  vectorReco[6] = fFake;
432  vectorReco[7] = fMultiplicity;
433  break;
435  vectorReco[0] = pt;
436  vectorReco[1] = rapidity;
437  vectorReco[2] = fCentValue;
438  vectorReco[3] = fMultiplicity;
439  break;
441  vectorReco[0] = pt;
442  vectorReco[1] = rapidity;
443  vectorReco[2] = fCentValue; // centrality
444  vectorReco[3] = fMultiplicity; // multiplicity (diff estimators can be used)
445  vectorReco[4] = localmult; // local multiplicity (Ntracks in DeltaEta<0.1 and DeltaPhi<0.1)
446  vectorReco[5] = fq2; // magnitude of reduced flow vector (computed using TPC tracks)
447  break;
448  }
449 
450  bFillRecoValues = kTRUE;
451 
452  return bFillRecoValues;
453 }
454 
455 
456 //_____________________________________________________________
458 {
459  // check the required decay channel
460 
461  Bool_t checkCD = kFALSE;
462 
463  Int_t daughter0 = fmcPartCandidate->GetDaughterLabel(0);
464  Int_t daughter1 = fmcPartCandidate->GetDaughterLabel(1);
465  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
466  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
467 
468  if (!mcPartDaughter0 || !mcPartDaughter1) {
469  AliDebug (2,"Problems in the MC Daughters\n");
470  return checkCD;
471  }
472 
473  if(TMath::Abs(fmcPartCandidate->GetPdgCode()) == 4122 && (daughter1 - daughter0 != 1)) {
474  AliDebug(2, Form("The MC particle doesn't have the correct daughters!!"));
475  return checkCD;
476  }
477 
478  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGneutrDaughForMC &&
479  TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGbachelor) &&
480  !(TMath::Abs(mcPartDaughter0->GetPdgCode()) == fPDGbachelor &&
481  TMath::Abs(mcPartDaughter1->GetPdgCode()) == fPDGneutrDaughForMC)) {
482  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()));
483  return checkCD;
484  }
485 
486  // the Neutral Particle (e.g. D0 for D*, K0S for Lc...)
487  AliAODMCParticle* mcPartDaughterNeutrDaugh = NULL;
488 
489  // for D* the D0 (the neutral) is the first daughter, while for Lc the V0 is the second, so we check the
490  // charge of teh daughters to decide which is which
491  AliDebug(3, Form("Charge0 = %d, Charge1 = %d", mcPartDaughter0->Charge()/3, mcPartDaughter1->Charge()/3));
492  if (mcPartDaughter0->Charge()/3 != 0){
493  mcPartDaughterNeutrDaugh = mcPartDaughter1;
494  }
495  else {
496  mcPartDaughterNeutrDaugh = mcPartDaughter0;
497  }
498 
499  Double_t vectorNeutrDaugh[2] ={0., 0.};
500 
501  // We are looking at a cascade ...evaluate the correct cascade
502  if (!EvaluateIfCorrectNeutrDaugh(mcPartDaughterNeutrDaugh, vectorNeutrDaugh)) {
503  AliDebug(2, "Error! the NeutrDaugh MC doesn't have correct daughters!!");
504  return checkCD;
505  }
506 
507  checkCD = kTRUE;
508  return checkCD;
509 
510 }
511 
512 //__________________________________________
513 Bool_t AliCFVertexingHFCascade::EvaluateIfCorrectNeutrDaugh(AliAODMCParticle* neutralDaugh, Double_t* vectorNeutrDaugh)const
514 {
515  //
516  // check wether D0 is decaing into kpi
517  //
518 
519  Bool_t isHadronic = kFALSE;
520  AliDebug(2, Form("neutralDaugh = %p, pdg = %d", neutralDaugh, neutralDaugh->GetPdgCode()));
521 
522  if (fPDGcascade == 4122) {
523  Int_t labelresonanceDaugh = neutralDaugh->GetDaughterLabel(0);
524  AliAODMCParticle* resonanceDaugh = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelresonanceDaugh));
525  if (!resonanceDaugh){
526  return kFALSE;
527  }
528  else {
529  AliDebug(3, Form("The daughter of the resonant particle is a %d (we are looking for a %d)", resonanceDaugh->GetPdgCode(), fPDGneutrDaugh));
530  if (TMath::Abs(resonanceDaugh->GetPdgCode()) != fPDGneutrDaugh){
531  return kFALSE;
532  }
533  else {
534  neutralDaugh = resonanceDaugh;
535  }
536  }
537  }
538 
539  Int_t daughterNeutrDaugh0 = neutralDaugh->GetDaughterLabel(0);
540  Int_t daughterNeutrDaugh1 = neutralDaugh->GetDaughterLabel(1);
541 
542  AliDebug(2, Form("daughter0 = %d and daughter1 = %d", daughterNeutrDaugh0, daughterNeutrDaugh1));
543  if (daughterNeutrDaugh0 == 0 || daughterNeutrDaugh1 == 0) {
544  AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
545  return isHadronic;
546  }
547 
548  Int_t numberOfExpectedDaughters = 2;
549  if (TMath::Abs(daughterNeutrDaugh1 - daughterNeutrDaugh0) != numberOfExpectedDaughters-1) { // should be everytime true - see PDGdatabooklet
550  AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
551  return isHadronic;
552  }
553 
554  AliAODMCParticle* mcPartDaughterNeutrDaugh0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh0));
555  AliAODMCParticle* mcPartDaughterNeutrDaugh1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterNeutrDaugh1));
556  if (!mcPartDaughterNeutrDaugh0 || !mcPartDaughterNeutrDaugh1) {
557  AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
558  return isHadronic;
559  }
560 
561  AliDebug(3, Form("Daughter 0 has pdg = %d, daughter 1 has pdg = %d", mcPartDaughterNeutrDaugh0->GetPdgCode(), mcPartDaughterNeutrDaugh1->GetPdgCode()));
562  if (!(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughPositive &&
563  TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughNegative) &&
564  !(TMath::Abs(mcPartDaughterNeutrDaugh0->GetPdgCode()) == fPDGneutrDaughNegative &&
565  TMath::Abs(mcPartDaughterNeutrDaugh1->GetPdgCode()) == fPDGneutrDaughPositive)) {
566  AliDebug(2, "The neutral particle (MC) doesn't come from the required decay, skipping!!");
567  return isHadronic;
568  }
569 
570  Double_t sumPxDau = mcPartDaughterNeutrDaugh0->Px()+mcPartDaughterNeutrDaugh1->Px();
571  Double_t sumPyDau = mcPartDaughterNeutrDaugh0->Py()+mcPartDaughterNeutrDaugh1->Py();
572  Double_t sumPzDau = mcPartDaughterNeutrDaugh0->Pz()+mcPartDaughterNeutrDaugh1->Pz();
573  Double_t pxMother = neutralDaugh->Px();
574  Double_t pyMother = neutralDaugh->Py();
575  Double_t pzMother = neutralDaugh->Pz();
576  AliDebug(3, Form("pxMother = %f, pyMother = %f, pzMother = %f", pxMother, pyMother, pzMother));
577  AliDebug(3, Form("sumPxDau = %f, sumPyDau = %f, sumPzDau = %f", sumPxDau, sumPyDau, sumPzDau));
578  if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>fCutOnMomConservation ||
579  TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>fCutOnMomConservation ||
580  TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>fCutOnMomConservation){
581  AliDebug(2, "Momentum conservation violated, skipping!!");
582  return isHadronic;
583  }
584 
585  Double_t pTNeutrDaughPositive = 0;
586  Double_t pTNeutrDaughNegative = 0;
587 
588  if (mcPartDaughterNeutrDaugh0->GetPdgCode() > 0 ) {
589  pTNeutrDaughPositive = mcPartDaughterNeutrDaugh0->Pt();
590  pTNeutrDaughNegative = mcPartDaughterNeutrDaugh1->Pt();
591  }
592  else {
593  pTNeutrDaughPositive = mcPartDaughterNeutrDaugh1->Pt();
594  pTNeutrDaughNegative = mcPartDaughterNeutrDaugh0->Pt();
595  }
596 
597  isHadronic = kTRUE;
598 
599  vectorNeutrDaugh[0] = pTNeutrDaughPositive;
600  vectorNeutrDaugh[1] = pTNeutrDaughNegative;
601 
602  return isHadronic;
603 
604 }
605 
606 //___________________________________________________________
607 
609 {
610  //
611  // setting the pt cut to be used in the Acceptance steps (MC+Reco)
612  //
613 
614  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");
615  if (fProngs>0){
616  for (Int_t iP=0; iP<fProngs; iP++){
617  fPtAccCut[iP]=ptAccCut[iP];
618  }
619  }
620  return;
621 }
622 
623 
624 
625 //___________________________________________________________
626 
628 {
629  //
630  // setting the eta cut to be used in the Acceptance steps (MC+Reco)
631  //
632 
633  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");
634  if (fProngs>0){
635  for (Int_t iP=0; iP<fProngs; iP++){
636  fEtaAccCut[iP] = etaAccCut[iP];
637  }
638  }
639  return;
640 }
641 //___________________________________________________________
642 
644 {
645  //
646  // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
647  //
648 
649  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");
650  if (fProngs>0){
651  for (Int_t iP=0; iP<fProngs; iP++){
652  fPtAccCut[iP]=ptAccCut[iP];
653  fEtaAccCut[iP]=etaAccCut[iP];
654  }
655  }
656  return;
657 }
658 
659 //___________________________________________________________
660 
662 {
663  //
664  // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
665  //
666 
667  Int_t bachelorPosition = 2;
668  if (fPDGcascade == 4122) bachelorPosition = 0;
669  AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[bachelorPosition])); // should be the soft pion...
670  if(!mcPartDaughter) return;
671  Int_t mother = mcPartDaughter->GetMother();
672  AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
673  if(!mcMother) return;
674 
675  if (TMath::Abs(mcPartDaughter->GetPdgCode()) != fPDGbachelor || TMath::Abs(mcMother->GetPdgCode()) != fPDGcascade){
676  AliError(Form("Apparently the expected bachelor is not in the third position, causing an error (pdg expected = %d, actual = %d)!!", fPDGbachelor, mcPartDaughter->GetPdgCode()));
677  AliError("This should be fixed when checking the MC part family in the CF task...");
678  return;
679  }
680  if (fProngs>0){
681  for (Int_t iP=0; iP<fProngs; iP++){
682  fPtAccCut[iP]=0.1;
683  fEtaAccCut[iP]=0.9;
684  }
685 
686  if (fPDGcascade != 4122){
687  fPtAccCut[2]=0.06; // soft pion
688  fEtaAccCut[2]=0.9; // soft pion
689  }
690  }
691  return;
692 }
693 
694 //_____________________________________________________________
696 {
697  //
698  // getting eta of the prong - overload the mother class method
699  //
700 
701  if (fRecoCandidate){
702 
704 
705  Double_t etaProng =-9999;
706  AliAODRecoDecay* neutrDaugh=0;
707  Int_t ibachelor = 0;
708  if (fPDGcascade == 413) {
709  neutrDaugh = cascade->Get2Prong();
710  }
711  else if (fPDGcascade == 4122) {
712  neutrDaugh = cascade->Getv0();
713  }
714  if (iProng==0) etaProng = neutrDaugh->EtaProng(0);
715  if (iProng==1) etaProng = neutrDaugh->EtaProng(1);
716  if (iProng==2) etaProng = cascade->EtaProng(ibachelor);
717 
718  return etaProng;
719 
720  }
721  return 999999;
722 }
723 //_____________________________________________________________
725 {
726  //
727  // getting pt of the prong
728  //
729 
730  if (fRecoCandidate){
731 
733  Double_t ptProng= -9999;
734  AliAODRecoDecay* neutrDaugh=0;
735  Int_t ibachelor = 0;
736  if (fPDGcascade == 413) {
737  neutrDaugh = cascade->Get2Prong();
738  }
739  else if (fPDGcascade == 4122) {
740  neutrDaugh = cascade->Getv0();
741  }
742  if (iProng == 0) ptProng = neutrDaugh->PtProng(0);
743  if (iProng == 1) ptProng = neutrDaugh->PtProng(1);
744  if (iProng == 2) ptProng = cascade->PtProng(ibachelor);
745 
746  // Double_t ptProng = fRecoCandidate->PtProng(iProng);
747  return ptProng;
748 
749  }
750  return 999999;
751 
752 }
753 //_____________________________________________________________
754 Bool_t AliCFVertexingHFCascade::CheckAdditionalCuts(AliPIDResponse* pidResponse) const {
755 
756  // function to check whether the candidate passes the additional cuts defined in the task to get the
757  // invariant mass spectra; these cuts are NOT pt-dependent
758 
759  if (fPDGcascade == 4122){
760  // the method is implemented only in this case so far
762  AliAODv0 * v0part = cascade->Getv0();
763  AliAODTrack *bachelor = (AliAODTrack*)cascade->GetBachelor();
764  Double_t bachelorEta = bachelor->Eta();
765  AliAODTrack *v0pos = (AliAODTrack*)v0part->GetDaughter(0);
766  AliAODTrack *v0neg = (AliAODTrack*)v0part->GetDaughter(1);
767  Double_t v0posEta = v0pos->Eta();
768  Double_t v0negEta = v0neg->Eta();
769 
770  Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
771  Double_t nSigmaTPCpr=-999.;
772  Double_t nSigmaTOFpr=-999.;
773  nSigmaTPCpr = pidResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
774  nSigmaTOFpr = pidResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
775  Double_t ptArm = v0part->PtArmV0();
776  Double_t invmassK0s = v0part->MassK0Short();
777  Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
778 
779  Bool_t cutsForTMVA = (TMath::Abs(bachelorEta) < 0.8 && TMath::Abs(v0posEta) < 0.8 && TMath::Abs(v0negEta) < 0.8) &&
780  ((nSigmaTOFpr < -800) || (TMath::Abs(nSigmaTOFpr) < 3)) &&
781  ((ptArm > 0.01 && ptArm < 0.07) || (ptArm > 0.105)) &&
782  ((TMath::Abs(invmassK0s - mK0SPDG)) < 0.01);
783 
784 
785  if (!fUseCutsForTMVA) cutsForTMVA = kTRUE;
786 
787  Bool_t cutsForInvMassTask = !(onFlyV0) &&
788  (cascade->CosV0PointingAngle()>0.99) &&
789  (TMath::Abs(nSigmaTPCpr) <= 3) &&
790  (v0part->Getd0Prong(0) < 20) &&
791  (v0part->Getd0Prong(1) < 20);
792 
793  if (cutsForTMVA && cutsForInvMassTask) {
794  // K0Smass cut
795  // eta cut
796  // TOF PID cut
797  // Arm cut
798  return kTRUE;
799  }
800  }
801 
802  return kFALSE;
803 
804 }
Int_t charge
fast configuration, only a subset of variables
AliCFVertexingHF & operator=(const AliCFVertexingHF &c)
Bool_t GetRecoValuesFromCandidate(Double_t *) const
double Double_t
Definition: External.C:58
Int_t fPDGneutrDaughPositive
pdg code of the V0
Double_t GetPtProng(Int_t iProng) const
TClonesArray * fmcArray
Class for HF corrections as a function of many variables and steps For D* and other cascades...
Double_t Ct(UInt_t pdg) const
Int_t ComputeLocalMultiplicity(Double_t etaD, Double_t phiD, Double_t R) const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliAODv0 * Getv0() const
AliAODVertex * fPrimVtx
pdg code of the negative daughter of the V0
AliAODMCParticle * fmcPartCandidate
Reconstructed HF candidate.
Int_t fPDGbachelor
pdg code of the cascade
AliCFVertexingHFCascade & operator=(const AliCFVertexingHFCascade &other)
TCanvas * c
Definition: TestFitELoss.C:172
Double_t CosPointingAngleXY() const
void SetEtaAccCut(Float_t *etaAccCut)
void SetNProngs(Int_t nProngs)
void SetPtAccCut(Float_t *ptAccCut)
Float_t * fPtAccCut
centrality value
Double_t GetEtaProng(Int_t iProng) const
Float_t fFake
fakes selection: 0 –> all, 1 –> non-fake, 2 –> fake
Int_t fConfiguration
array of tracks
int Int_t
Definition: External.C:63
Int_t fPDGneutrDaughForMC
pdg code of the V0
float Float_t
Definition: External.C:68
AliAODTrack * GetBachelor() const
slow configuration, all variables
Int_t fPDGneutrDaughNegative
pdg code of the positive daughter of the V0
rapidity
Definition: HFPtSpectrum.C:49
Bool_t GetGeneratedValuesFromMCParticle(Double_t *)
Double_t CosV0PointingAngle() const
short Short_t
Definition: External.C:23
Bool_t fUseCutsForTMVA
primaryVertex
Int_t * fLabelArray
n. of prongs
Int_t fPDGneutrDaugh
pdg code of the bachelor
AliAODRecoDecayHF * fRecoCandidate
mcArray candidate
Bool_t CheckAdditionalCuts(AliPIDResponse *pidResponse) const
decay
Definition: HFPtSpectrum.C:42
Int_t fmcLabel
flag to keep only the charm particles that comes from beauty decays
Double_t fzPrimVertex
get Number of variables for the container from the channel decay
Int_t NumberOfFakeDaughters() const
Int_t fProngs
results of the MatchToMC()
unsigned short UShort_t
Definition: External.C:28
AliAODVertex * GetPrimaryVtx() const
Double_t fq2
multiplicity of the event
bool Bool_t
Definition: External.C:53
void SetMCLabel(Int_t mcLabel)
AliAODRecoDecayHF2Prong * Get2Prong() const
Bool_t EvaluateIfCorrectNeutrDaugh(AliAODMCParticle *neutralDaugh, Double_t *VectorD0) const
Double_t fzMCVertex
Reco z primary vertex.
Double_t fMultiplicity
flag to remove events not geenrated with PYTHIA
Bool_t SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
Class for HF corrections as a function of many variables and step.
super fast configuration, only (pt,y,centrality)