Updates in PbPb cuts (Andrea Rossi)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFVertexingHFCascade.cxx
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
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
34 ClassImp(AliCFVertexingHFCascade)
35
36
37 //_________________________________________
38   AliCFVertexingHFCascade::AliCFVertexingHFCascade(TClonesArray *mcArray, UShort_t originDselection):
39   AliCFVertexingHF(mcArray, originDselection)
40 {       
41   // standard constructor
42
43   SetNProngs(3);
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...
47   fPtAccCut[0]=0.;
48   fEtaAccCut[0]=0.;
49   for(Int_t iP=1; iP<fProngs; iP++){
50           fPtAccCut[iP]=0.1;
51           fEtaAccCut[iP]=0.9;
52   }
53
54 }
55
56
57 //_____________________________________
58 AliCFVertexingHFCascade& AliCFVertexingHFCascade::operator=(const AliCFVertexingHFCascade& c)
59 {
60   // operator =
61  
62   if  (this != &c) {
63
64     AliCFVertexingHF::operator=(c);
65    
66   }
67   return *this;
68 }
69
70 //__________________________________________
71 Bool_t AliCFVertexingHFCascade::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand)
72 {
73   // set the AliAODRecoDecay candidate
74   
75   Bool_t bSignAssoc = kFALSE;
76  
77   fRecoCandidate = recoCand;
78   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)recoCand;
79
80   if (!fRecoCandidate) {
81     AliError("fRecoCandidate not found, problem in assignement\n");
82     return bSignAssoc;
83   }
84   
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");
87   
88   //Int_t pdgCand = 413;
89
90   Int_t pdgDgDStartoD0pi[2]={421,211};
91   Int_t pdgDgD0toKpi[2]={321,211};
92
93   Int_t nentries = fmcArray->GetEntriesFast();
94
95   AliDebug(3,Form("nentries = %d\n", nentries));
96  
97   Int_t mcLabel =  dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray); 
98   
99   if (mcLabel == -1) return bSignAssoc;
100
101   if (fRecoCandidate->NumberOfFakeDaughters()>0){
102           fFake = 0;    // fake candidate
103           if (fFakeSelection==1) return bSignAssoc;
104   }
105   if (fRecoCandidate->NumberOfFakeDaughters()==0){
106           fFake = 2;    // non-fake candidate
107           if (fFakeSelection==2) return bSignAssoc;
108   }
109   
110   SetMCLabel(mcLabel);
111   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel)); 
112
113   if (!fmcPartCandidate){
114     AliDebug(3,"No part candidate");
115     return bSignAssoc;
116   }
117  
118   bSignAssoc = kTRUE;
119   return bSignAssoc;
120 }
121
122 //______________________________________________
123 Bool_t AliCFVertexingHFCascade::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) 
124 {
125   // 
126   // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
127   //
128         
129         Bool_t bGenValues = kFALSE;
130         
131         
132         Int_t daughter0ds = fmcPartCandidate->GetDaughter(0);
133         Int_t daughter1ds = fmcPartCandidate->GetDaughter(1);
134
135         //the D0
136         AliAODMCParticle* mcPartDaughterD0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0ds));
137         AliAODMCParticle* mcPartDaughterPis = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1ds));
138
139         if(!mcPartDaughterD0 || !mcPartDaughterPis) return kFALSE;
140
141         Double_t vtx1[3] = {0,0,0};   // primary vertex         
142         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
143         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
144         fmcPartCandidate->XvYvZv(vtx1);  // cm
145
146         //D0 daughters
147         Int_t daughter0 = mcPartDaughterD0->GetDaughter(0);
148         Int_t daughter1 = mcPartDaughterD0->GetDaughter(1);
149
150         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0)); //D0
151         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1)); //pis
152
153         if(!mcPartDaughter0 || !mcPartDaughter1) return kFALSE;
154
155         // getting vertex from daughters
156         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
157         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
158         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) {
159           AliError("Daughters have different secondary vertex, skipping the track");
160           return bGenValues;
161         }
162         
163         Int_t nprongs = 2;
164         Short_t charge = 0;
165         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
166         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
167         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
168         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
169                 // inverting in case the positive daughter is the second one
170                 positiveDaugh = mcPartDaughter1;
171                 negativeDaugh = mcPartDaughter0;
172         }
173         // getting the momentum from the daughters
174         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};            
175         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};            
176         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
177         
178         Double_t d0[2] = {0.,0.};               
179         
180         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
181         
182         Double_t cosThetaStar = 0.;
183         Double_t cosThetaStarD0 = 0.;
184         Double_t cosThetaStarD0bar = 0.;
185         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
186         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
187         if (mcPartDaughterD0->GetPdgCode() == 421){  // D0
188           AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
189           cosThetaStar = cosThetaStarD0;
190         }
191         else if (mcPartDaughterD0->GetPdgCode() == -421){  // D0bar{
192           AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
193           cosThetaStar = cosThetaStarD0bar;
194         }
195         else{
196           AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
197           return vectorMC;
198         }
199         if (cosThetaStar < -1 || cosThetaStar > 1) {
200           AliWarning("Invalid value for cosine Theta star, returning");
201           return bGenValues;
202         }
203         
204         Double_t vectorD0[2] ={0.,0.};
205
206         // evaluate the correct cascade
207         if (!EvaluateIfD0toKpi(mcPartDaughterD0,vectorD0)) {
208           AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
209           return bGenValues;  
210         }       
211
212         //ct
213         Double_t cT = decay->Ct(421);
214         // get the pT of the daughters
215         Double_t pTD0 = 0.;
216         Double_t pTpi = 0.;
217         
218         if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 413) {
219                 pTD0 = mcPartDaughterD0->Pt();
220                 pTpi = mcPartDaughterPis->Pt();
221         }
222
223         
224         vectorMC[0] = fmcPartCandidate->Pt();
225         vectorMC[1] = fmcPartCandidate->Y() ;
226         vectorMC[2] = cosThetaStar ;
227         vectorMC[3] = vectorD0[0]; 
228         vectorMC[4] = vectorD0[1];
229         vectorMC[5] = cT*1.E4 ;  // in micron
230         vectorMC[6] = 0.;   // dummy value, meaningless in MC
231         vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
232         vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
233         vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
234         vectorMC[10] = 1.01;    // dummy value, meaningless in MC
235         vectorMC[11] = fmcPartCandidate->Phi(); 
236         vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
237         vectorMC[13] = fCentValue; // reconstructed centrality
238         vectorMC[14] = 1.;           // always filling with 1 at MC level 
239
240         delete decay;
241         bGenValues = kTRUE;
242         return bGenValues;
243 }
244 //____________________________________________
245 Bool_t AliCFVertexingHFCascade::GetRecoValuesFromCandidate(Double_t *vectorReco) const
246
247   // read the variables for the container
248
249   Bool_t bFillRecoValues=kFALSE;
250   
251   //Get the D* and the D0 from D*
252   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
253   AliAODRecoDecayHF2Prong* d0toKpi = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
254   
255
256   if (dstarD0pi->GetPrimaryVtx())printf("dstarD0pi has primary vtx\n");
257   if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidateDstar has primary vtx\n");
258
259   Double_t pt =  dstarD0pi->Pt();
260   Double_t rapidity =  dstarD0pi->YDstar();
261   Double_t invMass=0.;
262   Double_t cosThetaStar = 9999.;
263   Double_t pTpi = 0.;
264   Double_t pTK = 0.;
265   Double_t dca = d0toKpi->GetDCA();
266   Double_t d0pi = 0.;
267   Double_t d0K = 0.;
268   Double_t d0xd0 = d0toKpi->Prodd0d0();
269   Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
270   Double_t phi = dstarD0pi->Phi();
271
272   Int_t pdgCode = fmcPartCandidate->GetPdgCode();
273  
274   if (pdgCode > 0){
275     cosThetaStar = d0toKpi->CosThetaStarD0();
276     pTpi = d0toKpi->PtProng(0);
277     pTK = d0toKpi->PtProng(1);
278     d0pi = d0toKpi->Getd0Prong(0);
279     d0K = d0toKpi->Getd0Prong(1);
280     invMass=d0toKpi->InvMassD0();
281   }
282   else {
283     cosThetaStar = d0toKpi->CosThetaStarD0bar();
284     pTpi = d0toKpi->PtProng(1);
285     pTK = d0toKpi->PtProng(0);
286     d0pi = d0toKpi->Getd0Prong(1);
287     d0K = d0toKpi->Getd0Prong(0);
288     invMass= d0toKpi->InvMassD0bar();
289   }
290   
291   Double_t cT = d0toKpi->CtD0();
292   
293   vectorReco[0] = pt;
294   vectorReco[1] = rapidity;
295   vectorReco[2] = cosThetaStar;
296   vectorReco[3] = pTpi;
297   vectorReco[4] = pTK;
298   vectorReco[5] = cT*1.E4;  // in micron
299   vectorReco[6] = dca*1.E4;  // in micron
300   vectorReco[7] = d0pi*1.E4;  // in micron
301   vectorReco[8] = d0K*1.E4;  // in micron
302   vectorReco[9] = d0xd0*1.E8;  // in micron^2
303   vectorReco[10] = cosPointingAngle;  // in micron
304   vectorReco[11] = phi;  
305   vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
306   vectorReco[13] = fCentValue;
307   vectorReco[14] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
308
309   bFillRecoValues = kTRUE;
310
311   return bFillRecoValues;
312 }
313
314
315 //_____________________________________________________________
316 Bool_t AliCFVertexingHFCascade::CheckMCChannelDecay() const
317
318   // check the required decay channel
319
320   Bool_t checkCD = kFALSE;
321
322   
323   Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
324   Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
325   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
326   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
327
328   if (!mcPartDaughter0 || !mcPartDaughter1) {
329     AliDebug (2,"Problems in the MC Daughters\n");
330     return checkCD;
331   }
332
333   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
334         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
335       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
336         TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
337     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
338     return checkCD;  
339   }
340   
341   Double_t vectorD0[2] ={0.,0.};
342
343   // D* is a cascade ...evaluate the correct cascade
344   if (!EvaluateIfD0toKpi(mcPartDaughter0,vectorD0)) {
345     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
346     return checkCD;  
347   }
348    
349   checkCD = kTRUE;
350   return checkCD;
351   
352 }
353
354 //__________________________________________
355 Bool_t AliCFVertexingHFCascade::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, Double_t* vectorD0)const
356 {  
357   //
358   // chack wether D0 is decaing into kpi
359   //
360   
361   Bool_t isHadronic = kFALSE;
362   
363   Int_t daughterD00 = neutralDaugh->GetDaughter(0);
364   Int_t daughterD01 = neutralDaugh->GetDaughter(1);
365   
366   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
367   if (daughterD00 == 0 || daughterD01 == 0) {
368     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
369     return isHadronic;  
370   }
371   
372   if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
373     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
374     return isHadronic;  
375   }
376   
377   AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD00));
378   AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughterD01));
379   if (!mcPartDaughterD00 || !mcPartDaughterD01) {
380     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping"); 
381     return isHadronic;  
382   }
383   
384   if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
385         TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) && 
386       !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
387         TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
388     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
389     return isHadronic;  
390   }
391   
392   Double_t pTD0pi = 0;
393   Double_t pTD0K = 0;
394   
395   
396   if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
397     pTD0pi = mcPartDaughterD00->Pt();
398     pTD0K = mcPartDaughterD01->Pt();
399   }
400   else {
401     pTD0pi = mcPartDaughterD01->Pt();
402     pTD0K  = mcPartDaughterD00->Pt();
403   }
404   
405   isHadronic = kTRUE;
406   
407   vectorD0[0] = pTD0pi;
408   vectorD0[1] = pTD0K;
409  
410   return isHadronic;
411
412 }
413
414 //___________________________________________________________
415
416 void AliCFVertexingHFCascade::SetPtAccCut(Float_t* ptAccCut)
417 {
418         //
419         // setting the pt cut to be used in the Acceptance steps (MC+Reco)
420         //
421
422         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");
423         if (fProngs>0){
424                 for (Int_t iP=0; iP<fProngs; iP++){
425                         fPtAccCut[iP]=ptAccCut[iP];
426                 }
427         }
428         return;
429 }               
430
431
432
433 //___________________________________________________________
434
435 void AliCFVertexingHFCascade::SetEtaAccCut(Float_t* etaAccCut)
436 {
437         //
438         // setting the eta cut to be used in the Acceptance steps (MC+Reco)
439         //
440
441         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");
442         if (fProngs>0){
443                 for (Int_t iP=0; iP<fProngs; iP++){
444                         fEtaAccCut[iP]=etaAccCut[iP];
445                 }
446         }
447         return;
448 }       
449 //___________________________________________________________
450
451 void AliCFVertexingHFCascade::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
452 {
453         //
454         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
455         //
456
457         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");
458         if (fProngs>0){
459                 for (Int_t iP=0; iP<fProngs; iP++){
460                         fPtAccCut[iP]=ptAccCut[iP];
461                         fEtaAccCut[iP]=etaAccCut[iP];
462                 }
463         }
464         return;
465 }               
466
467 //___________________________________________________________
468
469 void AliCFVertexingHFCascade::SetAccCut()
470 {
471         //
472         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
473         //
474   
475   AliAODMCParticle* mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[2]));  // should be the soft pion...  
476   if(!mcPartDaughter) return;
477   Int_t mother =  mcPartDaughter->GetMother();
478   AliAODMCParticle* mcMother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); 
479   if(!mcMother) return;
480
481   if (TMath::Abs(mcPartDaughter->GetPdgCode())!= 211 || TMath::Abs(mcMother->GetPdgCode())!=413){
482     AliFatal("Apparently the soft pion is not in the third position, causing a crash!!");
483   }              
484   if (fProngs>0){
485     for (Int_t iP=0; iP<fProngs-1; iP++){
486       fPtAccCut[iP]=0.1;
487       fEtaAccCut[iP]=0.9;
488     }
489     fPtAccCut[2]=0.06;  // soft pion
490     fEtaAccCut[2]=0.9;  // soft pion
491   }
492   return;
493 }               
494
495 //_____________________________________________________________
496 Double_t AliCFVertexingHFCascade::GetEtaProng(Int_t iProng) const 
497 {
498         //
499         // getting eta of the prong - overload the mother class method
500         //
501
502  if (fRecoCandidate){
503    
504    AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
505
506    Double_t etaProng =-9999;
507     if(iProng==0) etaProng =dstarD0pi->Get2Prong()->EtaProng(0);
508     if(iProng==1) etaProng =dstarD0pi->Get2Prong()->EtaProng(1);
509     if(iProng==2) etaProng =dstarD0pi->EtaProng(1);
510     
511     return etaProng;
512     
513   }
514   return 999999;    
515 }
516 //_____________________________________________________________
517 Double_t AliCFVertexingHFCascade::GetPtProng(Int_t iProng) const 
518 {
519         //
520         // getting pt of the prong
521         //
522   
523   if (fRecoCandidate){
524
525     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)fRecoCandidate;
526     Double_t ptProng= -9999;
527     if(iProng==0) ptProng =dstarD0pi->Get2Prong()->PtProng(0);
528     if(iProng==1) ptProng =dstarD0pi->Get2Prong()->PtProng(1);
529     if(iProng==2) ptProng =dstarD0pi->PtProng(1);
530     
531     //  Double_t ptProng = fRecoCandidate->PtProng(iProng);  
532     return ptProng;
533     
534   }
535   return 999999;  
536   
537 }