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"
34 ClassImp(AliCFVertexingHFCascade)
37 //_________________________________________
38 AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
39 AliCFVertexingHF(mcArray, originDselection)
41 // standard constructor
44 fPtAccCut=new Float_t[fProngs];
45 fEtaAccCut=new Float_t[fProngs];
46 // element 0 in the cut arrays corresponds to the soft pion!!!!!!!! Careful when setting the values...
49 for(Int_t iP=1; iP<fProngs; iP++){
57 //_____________________________________
58 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
64 AliCFVertexingHF::operator=(c);
70 //__________________________________________
71 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
73 // set the AliAODRecoDecay candidate
75 Bool_t bSignAssoc = kFALSE;
77 fRecoCandidate = recoCand;
78 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)recoCand;
80 if (!fRecoCandidate) {
81 AliError("fRecoCandidate not found, problem in assignement\n");
85 if ( fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
86 //if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
88 //Int_t pdgCand = 413;
90 Int_t pdgDgDStartoD0pi[2]={421,211};
91 Int_t pdgDgD0toKpi[2]={321,211};
93 Int_t nentries = fmcArray->GetEntriesFast();
95 AliDebug(3,Form("nentries = %d\n", nentries));
97 Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
99 if (mcLabel < 0) return bSignAssoc;
101 fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
103 if (!fmcPartCandidate){
104 AliDebug(3,"No part candidate");
112 //______________________________________________
113 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC)
116 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
119 Bool_t bGenValues = kFALSE;
122 Int_t daughter0ds = fmcPartCandidate->GetDaughter(0);
123 Int_t daughter1ds = fmcPartCandidate->GetDaughter(1);
126 AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0ds));
127 AliAODMCParticle* mcPartDaughterPis = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1ds));
129 Double_t vtx1[3] = {0,0,0}; // primary vertex
130 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
131 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
132 fmcPartCandidate->XvYvZv(vtx1); // cm
135 Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
136 Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);
138 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
139 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis
141 // getting vertex from daughters
142 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
143 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
144 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) {
145 AliError("Daughters have different secondary vertex, skipping the track");
151 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
152 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
153 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
154 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
155 // inverting in case the positive daughter is the second one
156 positiveDaugh = mcPartDaughter1;
157 negativeDaugh = mcPartDaughter0;
159 // getting the momentum from the daughters
160 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
161 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
162 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
164 Double_t d0[2] = {0.,0.};
166 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
168 Double_t cosThetaStar = 0.;
169 Double_t cosThetaStarD0 = 0.;
170 Double_t cosThetaStarD0bar = 0.;
171 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
172 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
173 if (mcPartDaughterD0->GetPdgCode() == 421){ // D0
174 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
175 cosThetaStar = cosThetaStarD0;
177 else if (mcPartDaughterD0->GetPdgCode() == -421){ // D0bar{
178 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
179 cosThetaStar = cosThetaStarD0bar;
182 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
185 if (cosThetaStar < -1 || cosThetaStar > 1) {
186 AliWarning("Invalid value for cosine Theta star, returning");
190 Double_t vectorD0[2] ={0.,0.};
192 // evaluate the correct cascade
193 if (!EvaluateIfD0toKpi(mcPartDaughterD0,vectorD0)) {
194 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
199 Double_t cT = decay->Ct(421);
200 // get the pT of the daughters
204 if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
205 pTD0 = mcPartDaughterD0->Pt();
206 pTpi = mcPartDaughterPis->Pt();
210 vectorMC[0] = fmcPartCandidate->Pt();
211 vectorMC[1] = fmcPartCandidate->Y() ;
212 vectorMC[2] = cosThetaStar ;
213 vectorMC[3] = vectorD0[0];
214 vectorMC[4] = vectorD0[1];
215 vectorMC[5] = cT*1.E4 ; // in micron
216 vectorMC[6] = 0.; // dummy value, meaningless in MC
217 vectorMC[7] = 0.; // dummy value, meaningless in MC, in micron
218 vectorMC[8] = 0.; // dummy value, meaningless in MC, in micron
219 vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
220 vectorMC[10] = 1.01; // dummy value, meaningless in MC
221 vectorMC[11] = fmcPartCandidate->Phi();
222 vectorMC[12] = fzMCVertex; // z of reconstructed of primary vertex
223 vectorMC[13] = fCentValue; // reconstructed centrality
229 //____________________________________________
230 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
232 // read the variables for the container
234 Bool_t bFillRecoValues=kFALSE;
236 //Get the D* and the D0 from D*
237 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
238 AliAODRecoDecayHF2Prong* d0toKpi = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
241 if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi has primary vtx\n");
242 if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
244 Double_t pt = dstarD0pi->Pt();
245 Double_t rapidity = dstarD0pi->YDstar();
247 Double_t cosThetaStar = 9999.;
250 Double_t dca = d0toKpi->GetDCA();
253 Double_t d0xd0 = d0toKpi->Prodd0d0();
254 Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
255 Double_t phi = dstarD0pi->Phi();
257 Int_t pdgCode = fmcPartCandidate->GetPdgCode();
260 cosThetaStar = d0toKpi->CosThetaStarD0();
261 pTpi = d0toKpi->PtProng(0);
262 pTK = d0toKpi->PtProng(1);
263 d0pi = d0toKpi->Getd0Prong(0);
264 d0K = d0toKpi->Getd0Prong(1);
265 invMass=d0toKpi->InvMassD0();
268 cosThetaStar = d0toKpi->CosThetaStarD0bar();
269 pTpi = d0toKpi->PtProng(1);
270 pTK = d0toKpi->PtProng(0);
271 d0pi = d0toKpi->Getd0Prong(1);
272 d0K = d0toKpi->Getd0Prong(0);
273 invMass= d0toKpi->InvMassD0bar();
276 Double_t cT = d0toKpi->CtD0();
279 vectorReco[1] = rapidity;
280 vectorReco[2] = cosThetaStar;
281 vectorReco[3] = pTpi;
283 vectorReco[5] = cT*1.E4; // in micron
284 vectorReco[6] = dca*1.E4; // in micron
285 vectorReco[7] = d0pi*1.E4; // in micron
286 vectorReco[8] = d0K*1.E4; // in micron
287 vectorReco[9] = d0xd0*1.E8; // in micron^2
288 vectorReco[10] = cosPointingAngle; // in micron
289 vectorReco[11] = phi;
290 vectorReco[12] = fzPrimVertex; // z of reconstructed of primary vertex
291 vectorReco[13] = fCentValue;
293 bFillRecoValues = kTRUE;
295 return bFillRecoValues;
299 //_____________________________________________________________
300 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
302 // check the required decay channel
304 Bool_t checkCD = kFALSE;
307 Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
308 Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
309 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
310 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
312 if (!mcPartDaughter0 || !mcPartDaughter1) {
313 AliDebug (2,"Problems in the MC Daughters\n");
317 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
318 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
319 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
320 TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
321 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
325 Double_t vectorD0[2] ={0.,0.};
327 // D* is a cascade ...evaluate the correct cascade
328 if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
329 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
338 //__________________________________________
339 Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
342 // chack wether D0 is decaing into kpi
345 Bool_t isHadronic = kFALSE;
347 Int_t daughterD00 = neutralDaugh->GetDaughter(0);
348 Int_t daughterD01 = neutralDaugh->GetDaughter(1);
350 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
351 if (daughterD00 == 0 || daughterD01 == 0) {
352 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
356 if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
357 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
361 AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD00));
362 AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD01));
363 if (!mcPartDaughterD00 || !mcPartDaughterD01) {
364 AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
368 if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
369 TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
370 !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
371 TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
372 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
380 if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
381 pTD0pi = mcPartDaughterD00->Pt();
382 pTD0K = mcPartDaughterD01->Pt();
385 pTD0pi = mcPartDaughterD01->Pt();
386 pTD0K = mcPartDaughterD00->Pt();
391 vectorD0[0] = pTD0pi;
398 //___________________________________________________________
400 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
403 // setting the pt cut to be used in the Acceptance steps (MC+Reco)
406 AliInfo("The 3rd element of the pt cut array will correspond to the cut applied to the soft pion - please check that it is correct");
408 for (Int_t iP=0; iP<fProngs; iP++){
409 fPtAccCut[iP]=ptAccCut[iP];
417 //___________________________________________________________
419 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
422 // setting the eta cut to be used in the Acceptance steps (MC+Reco)
425 AliInfo("The 3rd element of the eta cut array will correspond to the cut applied to the soft pion - please check that it is correct");
427 for (Int_t iP=0; iP<fProngs; iP++){
428 fEtaAccCut[iP]=etaAccCut[iP];
433 //___________________________________________________________
435 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
438 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
441 AliInfo("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");
443 for (Int_t iP=0; iP<fProngs; iP++){
444 fPtAccCut[iP]=ptAccCut[iP];
445 fEtaAccCut[iP]=etaAccCut[iP];
451 //___________________________________________________________
453 void AliCFVertexingHFCascade::SetAccCut()
456 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
459 AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(2)); // should be the soft pion...
460 if (TMath::Abs(fLabelArray[0]-fLabelArray[1] != 1) || TMath::Abs(fLabelArray[1]-fLabelArray[2] == 1) || TMath::Abs(fLabelArray[0]-fLabelArray[2] == 1) || TMath::Abs(mcPartDaughter->GetPdgCode())){
461 AliFatal("Apparently the soft pion is not in the third position, causing a crash!!");
464 for (Int_t iP=0; iP<fProngs-1; iP++){
468 fPtAccCut[3]=0.; // soft pion
469 fEtaAccCut[3]=0.9; // soft pion