99c31415ff0b8987775417251e640a0a6f2ab263
[u/mrichter/AliRoot.git] / PWGHF / 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 #include "AliCFTaskVertexingHF.h"
38
39 ClassImp(AliCFVertexingHF3Prong)
40
41 //_________________________________________
42 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay, UInt_t resonantDecay):
43 AliCFVertexingHF(),
44   fDecay(decay),
45   fGenDsOption(kCountResonant),
46   fResonantDecay(resonantDecay)
47  {
48   // 
49   SetNProngs(3);
50
51   fPtAccCut=new Float_t[fProngs];
52   fEtaAccCut=new Float_t[fProngs];
53   for(Int_t iP=0; iP<fProngs; iP++){
54           fPtAccCut[iP]=0.1;
55           fEtaAccCut[iP]=0.9;
56   }
57
58 }
59 //_________________________________________
60 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay):
61 AliCFVertexingHF(),
62   fDecay(decay),
63   fGenDsOption(kCountResonant),
64   fResonantDecay(0)
65  {
66   //
67   SetNProngs(3);
68
69   fPtAccCut=new Float_t[fProngs];
70   fEtaAccCut=new Float_t[fProngs];
71   for(Int_t iP=0; iP<fProngs; iP++){
72           fPtAccCut[iP]=0.1;
73           fEtaAccCut[iP]=0.9;
74   }
75
76 }
77 //_________________________________________
78 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay, UInt_t resonantDecay):
79   AliCFVertexingHF(mcArray, originDselection),
80   fDecay(decay),
81   fGenDsOption(kCountResonant),
82   fResonantDecay(resonantDecay)
83 {
84   //
85   SetNProngs(3);
86   fPtAccCut=new Float_t[fProngs];
87   fEtaAccCut=new Float_t[fProngs];
88   for(Int_t iP=0; iP<fProngs; iP++){
89           fPtAccCut[iP]=0.1;
90           fEtaAccCut[iP]=0.9;
91   }
92 }
93
94 //_________________________________________
95 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
96   AliCFVertexingHF(mcArray, originDselection),
97   fDecay(decay),
98   fGenDsOption(kCountResonant),
99   fResonantDecay(0)
100 {
101   //
102   SetNProngs(3);
103   fPtAccCut=new Float_t[fProngs];
104   fEtaAccCut=new Float_t[fProngs];
105   for(Int_t iP=0; iP<fProngs; iP++){
106           fPtAccCut[iP]=0.1;
107           fEtaAccCut[iP]=0.9;
108   }
109 }
110 //_____________________________________
111 AliCFVertexingHF3Prong& AliCFVertexingHF3Prong::operator=(const AliCFVertexingHF3Prong& c){
112   //
113   if  (this != &c) {
114
115     AliCFVertexingHF::operator=(c);
116    
117   }
118     return *this;
119 }
120
121 //__________________________________________
122 Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
123   // Checks if candidate is signal and D meson is present in MC array
124   
125   Bool_t bSignAssoc = kFALSE;
126   fRecoCandidate = recoCand;
127
128   if (!fRecoCandidate) {
129     AliError("fRecoCandidate not found, problem in assignement\n");
130     return bSignAssoc;
131   }
132   
133   Int_t pdgCand = -1;
134   Int_t pdgDaughter[3]={-1,-1,-1};
135   if(fDecay==kDplustoKpipi){
136     pdgCand=411;
137     pdgDaughter[0]=321;
138     pdgDaughter[1]=211;
139     pdgDaughter[2]=211;
140   }else if(fDecay==kDstoKKpi){
141     pdgCand=431;
142     pdgDaughter[0]=321;
143     pdgDaughter[1]=321;
144     pdgDaughter[2]=211;
145   }else if(fDecay==kLctopKpi){
146           pdgCand=4122;
147           pdgDaughter[0]=2212;
148           pdgDaughter[1]=321;
149           pdgDaughter[2]=211;     
150   }else{
151     AliError("WRONG DECAY SETTING");
152     return bSignAssoc;    
153   }
154
155   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);  
156   if (mcLabel == -1) return bSignAssoc;
157
158   if (fRecoCandidate->NumberOfFakeDaughters()>0){
159           fFake = 0;    // fake candidate
160           if (fFakeSelection==1) return bSignAssoc;
161   }
162   if (fRecoCandidate->NumberOfFakeDaughters()==0){
163           fFake = 2;    // non-fake candidate
164           if (fFakeSelection==2) return bSignAssoc;
165   }
166   
167   SetMCLabel(mcLabel);
168   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
169     
170   if (!fmcPartCandidate){
171     AliDebug(3,"No part candidate");
172     return bSignAssoc;
173   }
174
175   if(fDecay==kDstoKKpi && fGenDsOption!=kCountAllDsKKpi){
176     if(!CheckMCChannelDecay()){
177       AliDebug(3,"Ds not from the selected resonant channel");
178       return bSignAssoc;
179     }
180   }
181
182   if (fDecay==kLctopKpi && fResonantDecay != AliCFTaskVertexingHF::kAll) {
183     if (!CheckLc3Prong()) return bSignAssoc;
184   }
185
186   bSignAssoc = kTRUE;
187   return bSignAssoc;
188 }
189
190 //______________________________________________
191 Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
192         // 
193         // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
194         // pt_D
195         // y_D
196         // phi_D
197         // ctau
198         // cos point
199         // pt_1
200         // pt_2
201         // pt_3
202         // d0_1
203         // d0_2
204         // d0_3
205         // zPrimVert
206         // centrality
207         
208         Bool_t bGenValues = kFALSE;
209         
210         Int_t pdgCand = -1;
211         if(fDecay==kDplustoKpipi){
212                 pdgCand=411;
213         }else if(fDecay==kDstoKKpi){
214                 pdgCand=431;
215         }else if(fDecay==kLctopKpi){
216                 pdgCand=4122;
217         }else{
218                 AliError("WRONG DECAY SETTING");
219                 return bGenValues;
220         }
221         
222         Double_t vertD[3] = {0,0,0};   // D origin
223         fmcPartCandidate->XvYvZv(vertD);  // cm
224         
225         Int_t nprongs = 3;
226         Int_t daughter[3];
227         Short_t charge = fmcPartCandidate->Charge();
228         
229         // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
230         // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
231         Int_t index=0;
232         Int_t nDauLS=0;
233         Int_t nDauOS=0;
234         
235         
236         Int_t nDau=fmcPartCandidate->GetNDaughters();
237         Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0); 
238         if(nDau==3){
239                 for(Int_t iDau=0; iDau<3; iDau++){
240                         Int_t ind = labelFirstDau+iDau;
241                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
242                         if(!part){
243                                 AliError("Daughter particle not found in MC array");
244                                 return bGenValues;
245                         }
246                         Short_t signDau=part->Charge();
247                         if(signDau==charge){
248                                 nDauLS++;
249                                 daughter[index] = ind;
250                                 index=2;
251                         }else{
252                                 daughter[1] = ind;
253                                 nDauOS++;
254                         }
255                 }
256         }else if(nDau==2){
257                 for(Int_t iDau=0; iDau<2; iDau++){
258                         Int_t ind = labelFirstDau+iDau;
259                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
260                         if(!part){
261                                 AliError("Daughter particle not found in MC array");
262                                 return bGenValues;
263                         }
264                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
265                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
266                                 Short_t signDau=part->Charge();
267                                 if(signDau==charge){
268                                         nDauLS++;
269                                         daughter[index] = ind;
270                                         index=2;
271                                 }else{
272                                         daughter[1] = ind;
273                                         nDauOS++;
274                                 }
275                         }else{
276                                 Int_t nDauRes=part->GetNDaughters();
277                                 if(nDauRes!=2){
278                                         AliError("Wrong resonant decay");
279                                         return bGenValues;
280                                 }
281                                 Int_t labelFirstDauRes = part->GetDaughter(0);  
282                                 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
283                                         Int_t indDR = labelFirstDauRes+iDauRes;
284                                         AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
285                                         if(!partDR){
286                                                 AliError("Daughter particle not found in MC array");
287                                                 return bGenValues;
288                                         }
289                                         Short_t signDau=partDR->Charge();
290                                         if(signDau==charge){
291                                                 nDauLS++;
292                                                 daughter[index] = ind;
293                                                 index=2;
294                                         }else{
295                                                 daughter[1] = ind;
296                                                 nDauOS++;
297                                         }
298                                 }
299                         }
300                 }
301         }else{
302                 AliError(Form("Wrong number of daughters %d",nDau));
303                 return bGenValues;
304         }
305         
306         if(nDauLS!=2 || nDauOS!=1){
307                 AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
308                 return bGenValues;
309         }
310         if(daughter[0]>daughter[2]){
311                 Int_t tmp=daughter[0];
312                 daughter[0]=daughter[2];
313                 daughter[2]=tmp;
314         }
315         
316         // getting the momentum from the daughters and decay vertex
317         Double_t px[3],py[3],pz[3],pt[3];
318         Double_t vertDec[3] = {0,0,0};   // decay vertex                
319         for(Int_t iDau=0; iDau<3; iDau++){
320                 AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
321                 if(!part){
322                         AliError("Daughter particle not found in MC array");
323                         return bGenValues;
324                 }
325                 px[iDau]=part->Px();
326                 py[iDau]=part->Py();
327                 pz[iDau]=part->Pz();
328                 pt[iDau]=part->Pt();
329                 if(iDau==0) part->XvYvZv(vertDec);
330         }
331         
332         Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
333         
334         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
335         Double_t cT = decay->Ct(pdgCand);
336         
337         switch (fConfiguration){
338         case AliCFTaskVertexingHF::kSnail:
339                 vectorMC[0] = fmcPartCandidate->Pt();
340                 vectorMC[1] = fmcPartCandidate->Y() ;
341                 vectorMC[2] = fmcPartCandidate->Phi();
342                 vectorMC[3] = cT*1.E4 ;  // in micron
343                 vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
344                 vectorMC[5] = pt[0];
345                 vectorMC[6] = pt[1];
346                 vectorMC[7] = pt[2];
347                 vectorMC[8] = fzMCVertex;    // z of reconstructed of primary vertex
348                 vectorMC[9] = fCentValue; // reconstructed centrality value 
349                 vectorMC[10] = 1.;           // fake: always filling with 1 at MC level 
350                 vectorMC[11] = 1.01; // dummy value for cosPointingXY  multiplicity
351                 vectorMC[12] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
352                 vectorMC[13] = fMultiplicity; // reconstructed multiplicity
353                 
354                 if (fDecay==kLctopKpi){
355                         vectorMC[11] = 0.; //dist12
356                         vectorMC[12] = 0.; //dist23
357                         vectorMC[13] = 0.; //sigmaVtx
358                         vectorMC[14] = 0.; //sumd02
359                         vectorMC[15] = 1.01; // dummy value for cosPointingXY  multiplicity
360                         vectorMC[16] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
361                         vectorMC[17] = fMultiplicity; // reconstructed multiplicity
362                 }
363                 break;
364                 
365         case AliCFTaskVertexingHF::kCheetah:
366                 vectorMC[0] = fmcPartCandidate->Pt();
367                 vectorMC[1] = fmcPartCandidate->Y() ;
368                 vectorMC[2] = cT*1.E4; // in micron
369                 vectorMC[3] = fmcPartCandidate->Phi();
370                 vectorMC[4] = fzMCVertex;
371                 vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
372                 vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level 
373                 vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
374                 break;
375         }
376         
377         bGenValues = kTRUE;
378         return bGenValues;
379 }
380
381
382 //____________________________________________
383 Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
384
385         // Fill vector (see above) with reconstructed quantities
386         Bool_t bFillRecoValues=kFALSE;
387         
388         Int_t pdgCand = -1;
389         if(fDecay==kDplustoKpipi){
390                 pdgCand=411;
391         }else if(fDecay==kDstoKKpi){
392                 pdgCand=431;
393         }else if(fDecay==kLctopKpi){
394                 pdgCand=4122;
395                 // AliError("LambdaC not yet implemented");
396                 // return bFillRecoValues;
397         }else{
398                 AliError("WRONG DECAY SETTING");
399                 return bFillRecoValues;
400         }
401         
402         AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
403         Short_t charge=decay3->Charge();
404         Double_t rapidity=decay3->Y(pdgCand);
405         Double_t cT=decay3->Ct(pdgCand); 
406         Double_t pt = decay3->Pt();
407         Double_t cosPointingAngle = decay3->CosPointingAngle();
408         Double_t phi = decay3->Phi();
409         Double_t dist12= decay3->GetDist12toPrim();
410         Double_t dist23 = decay3->GetDist23toPrim();
411         Double_t sigmVert = decay3->GetSigmaVert();
412         Double_t cosPointingAngleXY = decay3->CosPointingAngleXY();
413         Double_t normDecayLengthXY = decay3->NormalizedDecayLengthXY();
414                 
415         Int_t daughtSorted[3];
416         Int_t tmpIndex=0;
417         Int_t nDauLS=0;
418         Int_t nDauOS=0;
419         for(Int_t iDau=0; iDau<3; iDau++){
420                 AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
421                 Int_t label = TMath::Abs(trk->GetLabel());
422                 Short_t chargedau=trk->Charge();
423                 if(chargedau==charge){
424                         daughtSorted[tmpIndex]=label;
425                         tmpIndex=2;
426                         nDauLS++;
427                 }else{
428                         daughtSorted[1]=label;
429                         nDauOS++;
430                 }
431         }
432         
433         if(nDauLS!=2 || nDauOS!=1){
434                 AliError("Wrong decay channel: number of OS and LS tracks not OK");
435                 return bFillRecoValues;
436         }
437         
438         if(daughtSorted[0]>daughtSorted[2]){
439                 Int_t tmp=daughtSorted[0];
440                 daughtSorted[0]=daughtSorted[2];
441                 daughtSorted[2]=tmp;
442         }
443         
444         Double_t d0prong0 = decay3->Getd0Prong(daughtSorted[0]);
445         Double_t d0prong1 = decay3->Getd0Prong(daughtSorted[1]);
446         Double_t d0prong2 = decay3->Getd0Prong(daughtSorted[2]);
447         
448         switch (fConfiguration){
449         case AliCFTaskVertexingHF::kSnail:
450                 vectorReco[0] = pt;
451                 vectorReco[1] = rapidity;
452                 vectorReco[2] = phi;
453                 vectorReco[3] = cT*1.E4;  // in micron
454                 vectorReco[4] = cosPointingAngle;  // in micron
455                 vectorReco[5] = decay3->PtProng(daughtSorted[0]);
456                 vectorReco[6] = decay3->PtProng(daughtSorted[1]);
457                 vectorReco[7] = decay3->PtProng(daughtSorted[2]);
458                 vectorReco[8] = fzPrimVertex;    // z of reconstructed of primary vertex
459                 vectorReco[9] = fCentValue; //reconstructed centrality value
460                 vectorReco[10] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) 
461                 vectorReco[11] = cosPointingAngleXY; 
462                 vectorReco[12] = normDecayLengthXY; // in cm
463                 vectorReco[13] = fMultiplicity; // reconstructed multiplicity
464                 
465                 if(fDecay==kLctopKpi){  
466                         Double_t sumd02 =(d0prong0*d0prong0 + d0prong1*d0prong1 + d0prong2*d0prong2); 
467                         vectorReco[11] = dist12*1.E4;
468                         vectorReco[12] = dist23*1.E4;
469                         vectorReco[13] = sigmVert*1.E4;
470                         vectorReco[14] = sumd02*1.E8;
471                         vectorReco[15] = cosPointingAngleXY; 
472                         vectorReco[16] = normDecayLengthXY; // in cm
473                         vectorReco[17] = fMultiplicity; // reconstructed multiplicity
474                 }
475                 break;
476         case AliCFTaskVertexingHF::kCheetah:
477                 vectorReco[0] = pt;
478                 vectorReco[1] = rapidity ;
479                 vectorReco[2] = cT*1.E4; // in micron
480                 vectorReco[3] = phi; 
481                 vectorReco[4] = fzPrimVertex;
482                 vectorReco[5] = fCentValue;   
483                 vectorReco[6] = fFake ; 
484                 vectorReco[7] = fMultiplicity;  
485                 break;
486         }
487
488         bFillRecoValues = kTRUE;
489         return bFillRecoValues;
490 }
491
492
493 //_____________________________________________________________
494 Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
495
496   // Check the pdg codes of the daughters
497   Bool_t checkCD = kFALSE;
498
499   Int_t pdgCand = -1;
500   Int_t pdgDaughter[3]={-1,-1,-1};
501   if(fDecay==kDplustoKpipi){
502     pdgCand=411;
503     pdgDaughter[0]=321;
504     pdgDaughter[1]=211;
505     pdgDaughter[2]=211;
506   }else if(fDecay==kDstoKKpi){
507     pdgCand=431;
508     pdgDaughter[0]=321;
509     pdgDaughter[1]=321;
510     pdgDaughter[2]=211;
511   }else if(fDecay==kLctopKpi){
512           pdgCand=4122;
513           pdgDaughter[0]=2212;
514           pdgDaughter[1]=321;
515           pdgDaughter[2]=211;
516           
517   //  AliError("LambdaC not yet implemented");
518   //  return checkCD;
519   }else{
520     AliError("WRONG DECAY SETTING");
521     return checkCD;    
522   }
523
524
525   Int_t daughter[3];
526   Double_t sumPxDau=0.;
527   Double_t sumPyDau=0.;
528   Double_t sumPzDau=0.;
529
530   Int_t nDau=fmcPartCandidate->GetNDaughters();
531   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
532
533   if(fDecay==kLctopKpi && fResonantDecay!=AliCFTaskVertexingHF::kAll){if(!CheckLc3Prong()) return checkCD;}
534  
535   if(nDau==3){
536     if(fDecay==kDstoKKpi && !(fGenDsOption==kCountAllDsKKpi || fGenDsOption==kCountNonResonant)){
537       AliDebug(3,"Decay channel in direct KKpi, should be skipped");
538       return checkCD;
539     }
540     for(Int_t iDau=0; iDau<3; iDau++){
541       Int_t ind = labelFirstDau+iDau;
542       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
543       if(!part){
544         AliError("Daughter particle not found in MC array");
545         return checkCD;
546       }
547       daughter[iDau]=TMath::Abs(part->GetPdgCode());
548       sumPxDau+=part->Px();
549       sumPyDau+=part->Py();
550       sumPzDau+=part->Pz();
551     }
552   }else if(nDau==2){
553     if(fDecay==kDstoKKpi && fGenDsOption==kCountNonResonant) return checkCD;
554     Int_t nDauFound=0;
555     for(Int_t iDau=0; iDau<2; iDau++){
556       Int_t ind = labelFirstDau+iDau;
557       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
558       if(!part){
559         AliError("Daughter particle not found in MC array");
560         return checkCD;
561       }
562       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
563       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
564         if(nDauFound>=3) return checkCD;
565         daughter[nDauFound]=pdgCode;
566         sumPxDau+=part->Px();
567         sumPyDau+=part->Py();
568         sumPzDau+=part->Pz();
569         nDauFound++;
570       }else{
571         if(fDecay==kDstoKKpi && fGenDsOption!=3){
572           Int_t pdgCodeRes=TMath::Abs(part->GetPdgCode());
573           if(fGenDsOption==kCountPhipi && pdgCodeRes!=333) return checkCD;
574           else if(fGenDsOption==kCountK0stK && pdgCodeRes!=313) return checkCD;
575         }
576         Int_t nDauRes=part->GetNDaughters();
577         if(nDauRes!=2) return checkCD;
578         Int_t labelFirstDauRes = part->GetDaughter(0);  
579         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
580           Int_t indDR = labelFirstDauRes+iDauRes;
581           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
582           if(!partDR){
583             AliError("Daughter particle not found in MC array");
584             return checkCD;
585           }
586           Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
587           if(nDauFound>=3) return checkCD;
588           daughter[nDauFound]=pdgCodeDR;
589           sumPxDau+=partDR->Px();
590           sumPyDau+=partDR->Py();
591           sumPzDau+=partDR->Pz();
592           nDauFound++;
593         }
594       }
595     }
596   }else{
597     return checkCD;
598   }
599   for(Int_t iDau1=0; iDau1<3; iDau1++){
600     for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
601       if(daughter[iDau1]<daughter[iDau2]){
602         Int_t tmp=daughter[iDau1];
603         daughter[iDau1]=daughter[iDau2];
604         daughter[iDau2]=tmp;
605       }
606     }
607   }
608   for(Int_t iDau=0; iDau<3; iDau++){
609     if(daughter[iDau]!=pdgDaughter[iDau]){
610       AliDebug(2, "Wrong decay channel from MC, skipping!!");
611       return checkCD;  
612     }
613   }
614   
615   Double_t pxMother=fmcPartCandidate->Px();
616   Double_t pyMother=fmcPartCandidate->Py();
617   Double_t pzMother=fmcPartCandidate->Pz();
618   if(TMath::Abs(pxMother-sumPxDau)/(TMath::Abs(pxMother)+1.e-13)>0.00001 ||
619      TMath::Abs(pyMother-sumPyDau)/(TMath::Abs(pyMother)+1.e-13)>0.00001 ||
620      TMath::Abs(pzMother-sumPzDau)/(TMath::Abs(pzMother)+1.e-13)>0.00001){
621     AliDebug(2, "Momentum conservation violated, skipping!!");
622     return checkCD;  
623   }
624
625   checkCD = kTRUE;
626   return checkCD;
627   
628 }
629 //---------------------------------
630 Bool_t AliCFVertexingHF3Prong::CheckLc3Prong() const
631 {
632   Int_t numberOfLambdac=0;
633   if(TMath::Abs(fmcPartCandidate->GetPdgCode())!=4122) return kFALSE;
634   Int_t nDaugh = (Int_t)fmcPartCandidate->GetNDaughters();
635   if(nDaugh<2) return kFALSE;
636   if(nDaugh>3) return kFALSE;
637   AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)fmcArray->At(fmcPartCandidate->GetDaughter(0));
638   if(!pdaugh1) {return kFALSE;}
639   Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
640   AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)fmcArray->At(fmcPartCandidate->GetDaughter(1));
641   if(!pdaugh2) {return kFALSE;}
642   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
643   if(nDaugh==3){
644     if(fResonantDecay!=AliCFTaskVertexingHF::kNonResonant && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
645     Int_t thirdDaugh=fmcPartCandidate->GetDaughter(1)-1;
646     AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)fmcArray->At(thirdDaugh);
647     if(!pdaugh3) return kFALSE;
648     Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
649     if((number1==321 && number2==211 && number3==2212) || (number1==211 && number2==321 && number3==2212) || (number1==211 && number2==2212 && number3==321) || (number1==321 && number2==2212 && number3==211) || (number1==2212 && number2==321 && number3==211) || (number1==2212 && number2==211 && number3==321)) numberOfLambdac++;
650   }
651   if(nDaugh==2){
652     if(fResonantDecay==AliCFTaskVertexingHF::kNonResonant)return kFALSE; 
653     Int_t nfiglieK=0;
654     if((number1==2212 && number2==313)){
655       if(fResonantDecay!=AliCFTaskVertexingHF::kKstar && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
656       nfiglieK=pdaugh2->GetNDaughters();
657       if(nfiglieK!=2) return kFALSE;
658       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(0));
659       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(1));
660       if(!pdaughK1) return kFALSE;
661       if(!pdaughK2) return kFALSE;
662       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
663     }
664     if((number1==313 && number2==2212)){
665       if(fResonantDecay!=AliCFTaskVertexingHF::kKstar && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
666       nfiglieK=pdaugh1->GetNDaughters();
667       if(nfiglieK!=2) return kFALSE;
668       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(0));
669       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(1));
670       if(!pdaughK1) return kFALSE;
671       if(!pdaughK2) return kFALSE;
672       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
673     }
674     Int_t nfiglieDelta=0;
675     if(number1==321 && number2==2224){
676       if(fResonantDecay!=AliCFTaskVertexingHF::kDelta && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
677       nfiglieDelta=pdaugh2->GetNDaughters();
678       if(nfiglieDelta!=2) return kFALSE;
679       AliAODMCParticle *pdaughD1=(AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(0));
680       AliAODMCParticle *pdaughD2=(AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(1));
681       if(!pdaughD1) return kFALSE;
682       if(!pdaughD2) return kFALSE;
683       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
684     }
685     if(number1==2224 && number2==321){
686       if(fResonantDecay!=AliCFTaskVertexingHF::kDelta && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
687       nfiglieDelta=pdaugh1->GetNDaughters();
688       if(nfiglieDelta!=2) return kFALSE;
689       AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(0));
690       AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(1)); 
691       if(!pdaughD1) return kFALSE;
692       if(!pdaughD2) return kFALSE;
693       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
694     }
695     
696     Int_t nfiglieLa=0;
697     if(number1==3124 && number2==211){
698       if(fResonantDecay!=AliCFTaskVertexingHF::kL1520 && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
699       nfiglieLa=pdaugh1->GetNDaughters();
700       if(nfiglieLa!=2) return kFALSE;
701       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(0));
702       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)fmcArray->At(pdaugh1->GetDaughter(1));
703       if(!pdaughL1) return kFALSE;
704       if(!pdaughL2) return kFALSE;
705       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
706     }
707     if(number1==211 && number2==3124){
708       if(fResonantDecay!=AliCFTaskVertexingHF::kL1520 && fResonantDecay!=AliCFTaskVertexingHF::kAll)return kFALSE; 
709       nfiglieLa=pdaugh2->GetNDaughters();
710       if(nfiglieLa!=2) return kFALSE;
711       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(0));
712       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)fmcArray->At(pdaugh2->GetDaughter(1));
713       if(!pdaughL1) return kFALSE;
714       if(!pdaughL2) return kFALSE;
715       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
716       
717     }
718   }
719   if(numberOfLambdac>0) return kTRUE;
720   return kFALSE;
721 }