Update (Andrea)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFVertexingHF3Prong.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////
19 //                                                               //
20 // Class to compute variables for correction framework           //  
21 // for 3-body decays of D mesons (D+, Ds, Lc)                    //
22 // in bins of cut variables                                      //
23 // Origin:       Francesco Prino (prino@to.infn.it)              //
24 //               Renu Bala       (bala@to.infn.it)               //
25 //               Davide Caffarri (cafarri@pd.infn.it)            //
26 ///////////////////////////////////////////////////////////////////
27
28 #include "AliAODMCParticle.h"
29 #include "AliAODEvent.h"
30 #include "TClonesArray.h"
31 #include "AliCFVertexingHF.h"
32 #include "AliESDtrack.h"
33 #include "TDatabasePDG.h"
34
35 #include "AliCFVertexingHF3Prong.h"
36 #include "AliCFContainer.h"
37
38 ClassImp(AliCFVertexingHF3Prong)
39
40
41 //_________________________________________
42 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay):
43 AliCFVertexingHF(),
44   fDecay(decay)
45  {
46   // 
47   SetNProngs(3);
48
49   fPtAccCut=new Float_t[fProngs];
50   fEtaAccCut=new Float_t[fProngs];
51   for(Int_t iP=0; iP<fProngs; iP++){
52           fPtAccCut[iP]=0.1;
53           fEtaAccCut[iP]=0.9;
54   }
55
56 }
57 //_________________________________________
58 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
59   AliCFVertexingHF(mcArray, originDselection),
60   fDecay(decay)
61  
62 {
63   //
64   SetNProngs(3);
65   fPtAccCut=new Float_t[fProngs];
66   fEtaAccCut=new Float_t[fProngs];
67   for(Int_t iP=0; iP<fProngs; iP++){
68           fPtAccCut[iP]=0.1;
69           fEtaAccCut[iP]=0.9;
70   }
71 }
72
73
74 //_____________________________________
75 AliCFVertexingHF3Prong& AliCFVertexingHF3Prong::operator=(const AliCFVertexingHF3Prong& c){
76   //
77   if  (this != &c) {
78
79     AliCFVertexingHF::operator=(c);
80    
81   }
82     return *this;
83 }
84
85 //__________________________________________
86 Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
87   // Checks if candidate is signal and D meson is present in MC array
88   
89   Bool_t bSignAssoc = kFALSE;
90   fRecoCandidate = recoCand;
91
92   if (!fRecoCandidate) {
93     AliError("fRecoCandidate not found, problem in assignement\n");
94     return bSignAssoc;
95   }
96   
97   Int_t pdgCand = -1;
98   Int_t pdgDaughter[3]={-1,-1,-1};
99   if(fDecay==kDplustoKpipi){
100     pdgCand=411;
101     pdgDaughter[0]=321;
102     pdgDaughter[1]=211;
103     pdgDaughter[2]=211;
104   }else if(fDecay==kDstoKKpi){
105     pdgCand=431;
106     pdgDaughter[0]=321;
107     pdgDaughter[1]=321;
108     pdgDaughter[2]=211;
109   }else if(fDecay==kLctopKpi){
110           pdgCand=4122;
111           pdgDaughter[0]=2212;
112           pdgDaughter[1]=321;
113           pdgDaughter[2]=211;     
114   }else{
115     AliError("WRONG DECAY SETTING");
116     return bSignAssoc;    
117   }
118
119   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);  
120   if (mcLabel == -1) return bSignAssoc;
121
122   if (fRecoCandidate->NumberOfFakeDaughters()>0){
123           fFake = 0;    // fake candidate
124           if (fFakeSelection==1) return bSignAssoc;
125   }
126   if (fRecoCandidate->NumberOfFakeDaughters()==0){
127           fFake = 2;    // non-fake candidate
128           if (fFakeSelection==2) return bSignAssoc;
129   }
130   
131   SetMCLabel(mcLabel);
132   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
133     
134   if (!fmcPartCandidate){
135     AliDebug(3,"No part candidate");
136     return bSignAssoc;
137   }
138
139   bSignAssoc = kTRUE;
140   return bSignAssoc;
141 }
142
143 //______________________________________________
144 Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
145   // 
146   // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
147   // pt_D
148   // y_D
149   // phi_D
150   // ctau
151   // cos point
152   // pt_1
153   // pt_2
154   // pt_3
155   // d0_1
156   // d0_2
157   // d0_3
158   // zPrimVert
159   // centrality
160
161   Bool_t bGenValues = kFALSE;
162                 
163   Int_t pdgCand = -1;
164   if(fDecay==kDplustoKpipi){
165     pdgCand=411;
166   }else if(fDecay==kDstoKKpi){
167     pdgCand=431;
168   }else if(fDecay==kLctopKpi){
169         pdgCand=4122;
170   }else{
171     AliError("WRONG DECAY SETTING");
172     return bGenValues;
173   }
174
175   Double_t vertD[3] = {0,0,0};   // D origin
176   fmcPartCandidate->XvYvZv(vertD);  // cm
177
178   Int_t nprongs = 3;
179   Int_t daughter[3];
180   Short_t charge = fmcPartCandidate->Charge();
181
182   // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
183   // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
184   Int_t index=0;
185   Int_t nDauLS=0;
186   Int_t nDauOS=0;
187
188
189   Int_t nDau=fmcPartCandidate->GetNDaughters();
190   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0); 
191   if(nDau==3){
192     for(Int_t iDau=0; iDau<3; iDau++){
193       Int_t ind = labelFirstDau+iDau;
194       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
195       if(!part){
196         AliError("Daughter particle not found in MC array");
197         return bGenValues;
198       }
199       Short_t signDau=part->Charge();
200       if(signDau==charge){
201         nDauLS++;
202         daughter[index] = ind;
203         index=2;
204       }else{
205         daughter[1] = ind;
206         nDauOS++;
207       }
208     }
209   }else if(nDau==2){
210     for(Int_t iDau=0; iDau<2; iDau++){
211       Int_t ind = labelFirstDau+iDau;
212       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
213       if(!part){
214         AliError("Daughter particle not found in MC array");
215         return bGenValues;
216       }
217       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
218       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
219         Short_t signDau=part->Charge();
220         if(signDau==charge){
221           nDauLS++;
222           daughter[index] = ind;
223           index=2;
224         }else{
225           daughter[1] = ind;
226           nDauOS++;
227         }
228       }else{
229         Int_t nDauRes=part->GetNDaughters();
230         if(nDauRes!=2){
231           AliError("Wrong resonant decay");
232           return bGenValues;
233         }
234         Int_t labelFirstDauRes = part->GetDaughter(0);  
235         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
236           Int_t indDR = labelFirstDauRes+iDauRes;
237           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
238           if(!partDR){
239             AliError("Daughter particle not found in MC array");
240             return bGenValues;
241           }
242           Short_t signDau=partDR->Charge();
243           if(signDau==charge){
244             nDauLS++;
245             daughter[index] = ind;
246             index=2;
247           }else{
248             daughter[1] = ind;
249             nDauOS++;
250           }
251         }
252       }
253     }
254   }else{
255     AliError(Form("Wrong number of daughters %d",nDau));
256     return bGenValues;
257   }
258
259   if(nDauLS!=2 || nDauOS!=1){
260     AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
261     return bGenValues;
262   }
263   if(daughter[0]>daughter[2]){
264     Int_t tmp=daughter[0];
265     daughter[0]=daughter[2];
266     daughter[2]=tmp;
267   }
268   
269   // getting the momentum from the daughters and decay vertex
270   Double_t px[3],py[3],pz[3],pt[3];
271   Double_t vertDec[3] = {0,0,0};   // decay vertex              
272   for(Int_t iDau=0; iDau<3; iDau++){
273     AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
274     if(!part){
275       AliError("Daughter particle not found in MC array");
276       return bGenValues;
277     }
278     px[iDau]=part->Px();
279     py[iDau]=part->Py();
280     pz[iDau]=part->Pz();
281     pt[iDau]=part->Pt();
282     if(iDau==0) part->XvYvZv(vertDec);
283   }
284
285   Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
286   
287   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
288   Double_t cT = decay->Ct(pdgCand);
289                 
290   vectorMC[0] = fmcPartCandidate->Pt();
291   vectorMC[1] = fmcPartCandidate->Y() ;
292   vectorMC[2] = fmcPartCandidate->Phi();
293   vectorMC[3] = cT*1.E4 ;  // in micron
294   vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
295   vectorMC[5] = pt[0];
296   vectorMC[6] = pt[1];
297   vectorMC[7] = pt[2];
298   vectorMC[8] = 0.;   // imppar0, dummy value, meaningless in MC
299   vectorMC[9] = 0.;   // imppar1, dummy value, meaningless in MC, in micron
300   vectorMC[10] = 0.;   // imppar2, dummy value, meaningless in MC, in micron
301   vectorMC[11] = fzMCVertex;    // z of reconstructed of primary vertex
302   vectorMC[12] = fCentValue; // reconstructed centrality value 
303   vectorMC[13] = 1.;           // always filling with 1 at MC level 
304
305         if (fDecay==kLctopKpi){
306                 vectorMC[14] = 0.;
307                 vectorMC[15] = 0.;
308                 vectorMC[16] = 0.;
309                 vectorMC[17] = 0.;
310         }
311         
312
313   bGenValues = kTRUE;
314   return bGenValues;
315 }
316
317
318 //____________________________________________
319 Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
320
321   // Fill vector (see above) with reconstructed quantities
322   Bool_t bFillRecoValues=kFALSE;
323   
324   Int_t pdgCand = -1;
325   if(fDecay==kDplustoKpipi){
326     pdgCand=411;
327   }else if(fDecay==kDstoKKpi){
328     pdgCand=431;
329   }else if(fDecay==kLctopKpi){
330         pdgCand=4122;
331    // AliError("LambdaC not yet implemented");
332    // return bFillRecoValues;
333   }else{
334     AliError("WRONG DECAY SETTING");
335     return bFillRecoValues;
336   }
337
338   AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
339   Short_t charge=decay3->Charge();
340   Double_t rapidity=decay3->Y(pdgCand);
341   Double_t cT=decay3->Ct(pdgCand); 
342   Double_t pt = decay3->Pt();
343   Double_t cosPointingAngle = decay3->CosPointingAngle();
344   Double_t phi = decay3->Phi();
345   Double_t dist12= decay3->GetDist12toPrim();
346   Double_t dist23 = decay3->GetDist23toPrim();
347   Double_t sigmVert = decay3->GetSigmaVert();
348         
349
350   Int_t daughtSorted[3];
351   Int_t tmpIndex=0;
352   Int_t nDauLS=0;
353   Int_t nDauOS=0;
354   for(Int_t iDau=0; iDau<3; iDau++){
355     AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
356     Int_t label = TMath::Abs(trk->GetLabel());
357     Short_t chargedau=trk->Charge();
358     if(chargedau==charge){
359       daughtSorted[tmpIndex]=label;
360       tmpIndex=2;
361       nDauLS++;
362     }else{
363       daughtSorted[1]=label;
364       nDauOS++;
365     }
366   }
367
368   if(nDauLS!=2 || nDauOS!=1){
369     AliError("Wrong decay channel: number of OS and LS tracks not OK");
370     return bFillRecoValues;
371   }
372   
373   if(daughtSorted[0]>daughtSorted[2]){
374     Int_t tmp=daughtSorted[0];
375     daughtSorted[0]=daughtSorted[2];
376     daughtSorted[2]=tmp;
377   }
378   
379         Double_t d0prong0 = decay3->Getd0Prong(daughtSorted[0]);
380         Double_t d0prong1 = decay3->Getd0Prong(daughtSorted[1]);
381         Double_t d0prong2 = decay3->Getd0Prong(daughtSorted[2]);
382         
383   vectorReco[0] = pt;
384   vectorReco[1] = rapidity;
385   vectorReco[2] = phi;
386   vectorReco[3] = cT*1.E4;  // in micron
387   vectorReco[4] = cosPointingAngle;  // in micron
388   vectorReco[5] = decay3->PtProng(daughtSorted[0]);
389   vectorReco[6] = decay3->PtProng(daughtSorted[1]);
390   vectorReco[7] = decay3->PtProng(daughtSorted[2]);
391   vectorReco[8] = decay3->Getd0Prong(daughtSorted[0]);
392   vectorReco[9] = decay3->Getd0Prong(daughtSorted[1]);
393   vectorReco[10] = decay3->Getd0Prong(daughtSorted[2]);
394   vectorReco[11] = fzPrimVertex;    // z of reconstructed of primary vertex
395   vectorReco[12] = fCentValue; //reconstructed centrality value
396   vectorReco[13] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
397
398         if(fDecay==kLctopKpi){  
399                 Double_t sumd02 =(d0prong0*d0prong0 + d0prong1*d0prong1 + d0prong2*d0prong2); 
400                 vectorReco[14] = dist12*1.E4;
401                 vectorReco[15] = dist23*1.E4;
402                 vectorReco[16] = sigmVert*1.E4;
403                 vectorReco[17] = sumd02*1.E8;
404         }
405         
406
407   bFillRecoValues = kTRUE;
408   return bFillRecoValues;
409 }
410
411
412 //_____________________________________________________________
413 Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
414
415   // Check the pdg codes of the daughters
416   Bool_t checkCD = kFALSE;
417
418   Int_t pdgCand = -1;
419   Int_t pdgDaughter[3]={-1,-1,-1};
420   if(fDecay==kDplustoKpipi){
421     pdgCand=411;
422     pdgDaughter[0]=321;
423     pdgDaughter[1]=211;
424     pdgDaughter[2]=211;
425   }else if(fDecay==kDstoKKpi){
426     pdgCand=431;
427     pdgDaughter[0]=321;
428     pdgDaughter[1]=321;
429     pdgDaughter[2]=211;
430   }else if(fDecay==kLctopKpi){
431           pdgCand=4122;
432           pdgDaughter[0]=2212;
433           pdgDaughter[1]=321;
434           pdgDaughter[2]=211;
435           
436   //  AliError("LambdaC not yet implemented");
437   //  return checkCD;
438   }else{
439     AliError("WRONG DECAY SETTING");
440     return checkCD;    
441   }
442
443
444   Int_t daughter[3];
445
446   Int_t nDau=fmcPartCandidate->GetNDaughters();
447   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0); 
448   if(nDau==3){
449     for(Int_t iDau=0; iDau<3; iDau++){
450       Int_t ind = labelFirstDau+iDau;
451       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
452       if(!part){
453         AliError("Daughter particle not found in MC array");
454         return checkCD;
455       }
456       daughter[iDau]=TMath::Abs(part->GetPdgCode());
457     }
458   }else if(nDau==2){
459     Int_t nDauFound=0;
460     for(Int_t iDau=0; iDau<2; iDau++){
461       Int_t ind = labelFirstDau+iDau;
462       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
463       if(!part){
464         AliError("Daughter particle not found in MC array");
465         return checkCD;
466       }
467       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
468       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
469         if(nDauFound>=3) return checkCD;
470         daughter[nDauFound]=pdgCode;
471         nDauFound++;
472       }else{
473         Int_t nDauRes=part->GetNDaughters();
474         if(nDauRes!=2) return checkCD;
475         Int_t labelFirstDauRes = part->GetDaughter(0);  
476         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
477           Int_t indDR = labelFirstDauRes+iDauRes;
478           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
479           if(!partDR){
480             AliError("Daughter particle not found in MC array");
481             return checkCD;
482           }
483           Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
484           if(nDauFound>=3) return checkCD;
485           daughter[nDauFound]=pdgCodeDR;
486           nDauFound++;
487         }
488       }
489     }
490   }else{
491     return checkCD;
492   }
493   for(Int_t iDau1=0; iDau1<3; iDau1++){
494     for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
495       if(daughter[iDau1]<daughter[iDau2]){
496         Int_t tmp=daughter[iDau1];
497         daughter[iDau1]=daughter[iDau2];
498         daughter[iDau2]=tmp;
499       }
500     }
501   }
502   for(Int_t iDau=0; iDau<3; iDau++){
503     if(daughter[iDau]!=pdgDaughter[iDau]){
504       AliDebug(2, "Wrong decay channel from MC, skipping!!");
505       return checkCD;  
506     }
507   }
508   
509   checkCD = kTRUE;
510   return checkCD;
511   
512 }