Avoid comparing floats with != (ChiaraZ)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFVertexingHF2Prong.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 step 
18 // Author : C. Zampolli, CERN
19 // D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
20 // Base class for HF Unfolding - agrelli@uu.nl
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
31 #include "AliCFVertexingHF2Prong.h"
32 #include "AliCFContainer.h"
33
34 ClassImp(AliCFVertexingHF2Prong)
35
36
37 //_________________________________________
38   AliCFVertexingHF2Prong::AliCFVertexingHF2Prong(TClonesArray *mcArray, UShort_t originDselection):
39           AliCFVertexingHF(mcArray, originDselection)
40 {       
41         //
42         // constructor
43         //
44
45         SetNProngs(2);
46 }
47
48
49 //_____________________________________
50 AliCFVertexingHF2Prong& AliCFVertexingHF2Prong::operator=(const AliCFVertexingHF2Prong& c)
51 {
52         //
53         // copy constructor     
54         //
55
56         if  (this != &c) {              
57                 AliCFVertexingHF::operator=(c);         
58         }
59         return *this;
60 }
61
62 //__________________________________________
63 Bool_t AliCFVertexingHF2Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
64 {  
65         //
66         // setting the recontructed candidate
67         //
68         
69         Bool_t bSignAssoc = kFALSE;     
70         fRecoCandidate = recoCand;
71         if (!fRecoCandidate) {
72                 AliError("fRecoCandidate not found, problem in assignement\n");
73                 return bSignAssoc;
74         }
75         
76         if (fRecoCandidate->GetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n");
77         if (recoCand->GetPrimaryVtx()) AliDebug(3,"Reco Cand has a pointer to PrimVtx\n");
78         
79         Int_t pdgCand = 421;
80         Int_t pdgDgD0toKpi[2]={321,211};
81         Int_t nentries = fmcArray->GetEntriesFast();    
82
83         AliDebug(3,Form("nentries = %d\n", nentries));
84  
85         Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,2,pdgDgD0toKpi);
86         if (mcLabel == -1) return bSignAssoc;
87         SetMCLabel(mcLabel);
88         fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
89         if (!fmcPartCandidate){
90                 AliDebug(3,"No part candidate");
91                 return bSignAssoc;
92         }       
93         bSignAssoc = kTRUE;
94         return bSignAssoc;
95 }
96
97 //______________________________________________
98 Bool_t AliCFVertexingHF2Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) 
99 {
100         // 
101         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
102         //
103         
104         Bool_t bGenValues = kFALSE;
105         Double_t vtx1[3] = {0,0,0};   // primary vertex         
106         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
107         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
108         fmcPartCandidate->XvYvZv(vtx1);  // cm
109
110         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
111         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
112         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
113         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
114
115         // getting vertex from daughters
116         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
117         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
118         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) {
119                 AliError("Daughters have different secondary vertex, skipping the track");
120                 return bGenValues;
121         }
122         
123         Int_t nprongs = 2;
124         Short_t charge = 0;
125         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
126         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
127         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
128         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
129                 // inverting in case the positive daughter is the second one
130                 positiveDaugh = mcPartDaughter1;
131                 negativeDaugh = mcPartDaughter0;
132         }
133         // getting the momentum from the daughters
134         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};            
135         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};            
136         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
137         
138         Double_t d0[2] = {0.,0.};               
139         
140         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
141         
142         Double_t cosThetaStar = 0.;
143         Double_t cosThetaStarD0 = 0.;
144         Double_t cosThetaStarD0bar = 0.;
145         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
146         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
147         if (fmcPartCandidate->GetPdgCode() == 421){  // D0
148                 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
149                 cosThetaStar = cosThetaStarD0;
150         }
151         else if (fmcPartCandidate->GetPdgCode() == -421){  // D0bar{
152                 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
153                 cosThetaStar = cosThetaStarD0bar;
154         }
155         else{
156                 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
157                 return vectorMC;
158         }
159         if (cosThetaStar < -1 || cosThetaStar > 1) {
160                 AliWarning("Invalid value for cosine Theta star, returning");
161                 return bGenValues;
162         }
163                 
164         //ct
165         Double_t cT = decay->Ct(421);
166         // get the pT of the daughters
167         Double_t pTpi = 0.;
168         Double_t pTK = 0.;
169         
170         if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
171                 pTpi = mcPartDaughter0->Pt();
172                 pTK = mcPartDaughter1->Pt();
173         }
174         else {
175                 pTpi = mcPartDaughter1->Pt();
176                 pTK = mcPartDaughter0->Pt();
177         }
178         
179         vectorMC[0] = fmcPartCandidate->Pt();
180         vectorMC[1] = fmcPartCandidate->Y() ;
181         vectorMC[2] = cosThetaStar ;
182         vectorMC[3] = pTpi ;
183         vectorMC[4] = pTK ;
184         vectorMC[5] = cT*1.E4 ;  // in micron
185         vectorMC[6] = 0.;   // dummy value, meaningless in MC
186         vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
187         vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
188         vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
189         vectorMC[10] = 1.01;    // dummy value, meaningless in MC
190         vectorMC[11] = fmcPartCandidate->Phi(); 
191         vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
192         delete decay;
193         bGenValues = kTRUE;
194         return bGenValues;
195 }
196 //____________________________________________
197 Bool_t AliCFVertexingHF2Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
198 {
199         //
200         // Getting the reconstructed values from the candidate
201         // 
202         
203         Bool_t bFillRecoValues=kFALSE;
204         
205         AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)fRecoCandidate;
206         
207         if (d0toKpi->GetPrimaryVtx())AliDebug(3,"d0toKpi has primary vtx\n");
208         if (fRecoCandidate->GetPrimaryVtx())AliDebug(3,"fRecoCandidate has primary vtx\n");
209         
210         Double_t pt = d0toKpi->Pt();
211         Double_t rapidity = d0toKpi->YD0();
212         Double_t invMass=0.;
213         Double_t cosThetaStar = 9999.;
214         Double_t pTpi = 0.;
215         Double_t pTK = 0.;
216         Double_t dca = d0toKpi->GetDCA();
217         Double_t d0pi = 0.;
218         Double_t d0K = 0.;
219         Double_t d0xd0 = d0toKpi->Prodd0d0();
220         Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
221         Double_t phi = d0toKpi->Phi();
222         Int_t pdgCode = fmcPartCandidate->GetPdgCode();
223        
224         if (pdgCode > 0){
225                 cosThetaStar = d0toKpi->CosThetaStarD0();
226                 pTpi = d0toKpi->PtProng(0);
227                 pTK = d0toKpi->PtProng(1);
228                 d0pi = d0toKpi->Getd0Prong(0);
229                 d0K = d0toKpi->Getd0Prong(1);
230                 invMass=d0toKpi->InvMassD0();
231         }
232         else {
233                 cosThetaStar = d0toKpi->CosThetaStarD0bar();
234                 pTpi = d0toKpi->PtProng(1);
235                 pTK = d0toKpi->PtProng(0);
236                 d0pi = d0toKpi->Getd0Prong(1);
237                 d0K = d0toKpi->Getd0Prong(0);
238                 invMass= d0toKpi->InvMassD0bar();
239         }
240         
241         Double_t cT = d0toKpi->CtD0();  
242         
243         vectorReco[0] = pt;
244         vectorReco[1] = rapidity;
245         vectorReco[2] = cosThetaStar;
246         vectorReco[3] = pTpi;
247         vectorReco[4] = pTK;
248         vectorReco[5] = cT*1.E4;  // in micron
249         vectorReco[6] = dca*1.E4;  // in micron
250         vectorReco[7] = d0pi*1.E4;  // in micron
251         vectorReco[8] = d0K*1.E4;  // in micron
252         vectorReco[9] = d0xd0*1.E8;  // in micron^2
253         vectorReco[10] = cosPointingAngle;  // in micron
254         vectorReco[11] = phi;  
255         vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
256         bFillRecoValues = kTRUE;
257
258         return bFillRecoValues;
259 }
260
261 //_____________________________________________________________
262 Bool_t AliCFVertexingHF2Prong::CheckMCChannelDecay() const
263
264         //
265         // checking the MC decay channel
266         //
267         Bool_t checkCD = kFALSE;
268         Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
269         Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
270         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
271         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
272         
273         if (!mcPartDaughter0 || !mcPartDaughter1) {
274                 AliDebug (2,"Problems in the MC Daughters\n");
275                 return checkCD;
276         }
277         
278         if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
279               TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
280             !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
281               TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
282                 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
283                 return checkCD;  
284         }
285         
286         checkCD = kTRUE;
287         return checkCD;
288         
289 }
290