AliPhysics  master (3d17d9d)
AliAODRecoDecayHF2Prong.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2006, 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 
19 //
20 // Base class for AOD reconstructed heavy-flavour 2-prong decay
21 //
22 // Author: A.Dainese, andrea.dainese@lnl.infn.it
24 
25 #include <TDatabasePDG.h>
26 #include "AliAODRecoDecayHF.h"
29 #include "AliAODMCParticle.h"
30 
32 ClassImp(AliAODRecoDecayHF2Prong);
34 
35 //--------------------------------------------------------------------------
38 {
39  //
41  //
42 }
43 //--------------------------------------------------------------------------
45  Double_t *px,Double_t *py,Double_t *pz,
46  Double_t *d0,Double_t *d0err,Float_t dca) :
47  AliAODRecoDecayHF(vtx2,2,0,px,py,pz,d0,d0err)
48 {
49  //
51  //
52  SetDCA(dca);
53 }
54 //--------------------------------------------------------------------------
56  Double_t *d0,Double_t *d0err,Float_t dca) :
57  AliAODRecoDecayHF(vtx2,2,0,d0,d0err)
58 {
59  //
61  //
62  SetDCA(dca);
63 }
64 //--------------------------------------------------------------------------
66  AliAODRecoDecayHF(source)
67 {
68  //
70  //
71 }
72 //--------------------------------------------------------------------------
74 {
75  //
77  //
78  if(&source == this) return *this;
79 
81 
82  return *this;
83 }
84 //--------------------------------------------------------------------------
86  const {
103  Double_t mD0,mD0bar,ctsD0,ctsD0bar;
104  okD0=1; okD0bar=1;
105 
106  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
107 
108  if(PtProng(1) < cuts[3] || PtProng(0) < cuts[4]) okD0 = 0;
109  if(PtProng(0) < cuts[3] || PtProng(1) < cuts[4]) okD0bar = 0;
110  if(!okD0 && !okD0bar) return kFALSE;
111 
112  if(TMath::Abs(Getd0Prong(1)) > cuts[5] ||
113  TMath::Abs(Getd0Prong(0)) > cuts[6]) okD0 = 0;
114  if(TMath::Abs(Getd0Prong(0)) > cuts[6] ||
115  TMath::Abs(Getd0Prong(1)) > cuts[5]) okD0bar = 0;
116  if(!okD0 && !okD0bar) return kFALSE;
117 
118  if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
119 
120  InvMassD0(mD0,mD0bar);
121  if(TMath::Abs(mD0-mD0PDG) > cuts[0]) okD0 = 0;
122  if(TMath::Abs(mD0bar-mD0PDG) > cuts[0]) okD0bar = 0;
123  if(!okD0 && !okD0bar) return kFALSE;
124 
125  CosThetaStarD0(ctsD0,ctsD0bar);
126  if(TMath::Abs(ctsD0) > cuts[2]) okD0 = 0;
127  if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
128  if(!okD0 && !okD0bar) return kFALSE;
129 
130  if(Prodd0d0() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
131 
132  if(CosPointingAngle() < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
133 
134  return kTRUE;
135 }
136 //-----------------------------------------------------------------------------
138  const {
155  Double_t mJPsi,ctsJPsi;
156  okB=1;
157 
158  Double_t mJPSIPDG = TDatabasePDG::Instance()->GetParticle(443)->Mass();
159 
160  if(PtProng(1) < cuts[3] || PtProng(0) < cuts[4]) okB = 0;
161  if(!okB) return kFALSE;
162 
163  if(TMath::Abs(Getd0Prong(1)) > cuts[5] ||
164  TMath::Abs(Getd0Prong(0)) > cuts[6]) okB = 0;
165  if(!okB) return kFALSE;
166 
167  if(GetDCA() > cuts[1]) { okB = 0; return kFALSE; }
168 
169  mJPsi=InvMassJPSIee();
170  if(TMath::Abs(mJPsi-mJPSIPDG) > cuts[0]) okB = 0;
171  if(!okB) return kFALSE;
172 
173  ctsJPsi=CosThetaStarJPSI();
174  if(TMath::Abs(ctsJPsi) > cuts[2]) okB = 0;
175  if(!okB) return kFALSE;
176 
177  if(Prodd0d0() > cuts[7]) { okB = 0; return kFALSE; }
178 
179  if(CosPointingAngle() < cuts[8]) { okB = 0; return kFALSE; }
180 
181  return kTRUE;
182 }
183 //-----------------------------------------------------------------------------
184 Int_t AliAODRecoDecayHF2Prong::MatchToMCB2Prong(Int_t pdgabs,Int_t pdgabs2prong,Int_t *pdgDg,Int_t *pdgDg2prong,TClonesArray *mcArray) const
185 {
186  //
187  // Check if this candidate is matched to a MC signal
188  // If no, return -1
189  // If yes, return label (>=0) of the AliAODMCParticle
190  // If yes, but rejected by momentum conservation check, return -1*label (to allow checks at task level)
191  //
192  // NB: This function is only for non-resonant decays (for both 2 prongs)
193  // NB: Loosened cut on mom. conserv. (needed because of small issue in ITS Upgrade productions)
194 
195  Int_t signlabMother = 1;
196 
197  // Check number of daughters. Candidate is AliAODRecoDecayHF2Prong, so only continue when 2 daughters
198  Int_t ndg = GetNDaughters();
199  if (!ndg) { AliError("HF2Prong: No daughters available"); return -1;}
200  if (ndg != 2) { AliError(Form("HF2Prong: %d daughters instead of 2",ndg)); return -1;}
201 
202  // loop on daughters and write the labels
203  Int_t dgLabels[2] = {-1};
204  if(pdgabs2prong == 0){
205  for(Int_t i=0; i<ndg; i++) {
206  AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
207  dgLabels[i] = trk->GetLabel();
208  }
209  } else {
210  AliAODTrack *trk = (AliAODTrack*)GetDaughter(0);
211  dgLabels[0] = trk->GetLabel();
212  AliAODRecoDecayHF2Prong* daug2prong = (AliAODRecoDecayHF2Prong*)GetDaughter(1);
213  //Daughter prong (independent on their decay time) will also be far from PV, so call this function again.
214  Int_t pdgempty[2]={0,0};
215  dgLabels[1] = daug2prong->MatchToMCB2Prong(pdgabs2prong, 0, pdgDg2prong, pdgempty, mcArray);
216  if (dgLabels[1] < -1) signlabMother = -1;
217  }
218  if (dgLabels[0] == -1) return -1;
219  if (dgLabels[1] == -1) return -1;
220 
221  Int_t labMom[2] = {0, 0};
222  Int_t i, j, lab, labMother, pdgMother, pdgPart;
223  AliAODMCParticle *part = 0;
224  AliAODMCParticle *mother = 0;
225  Double_t pxSumDgs = 0., pySumDgs = 0., pzSumDgs = 0.;
226  Bool_t pdgUsed[2] = {kFALSE, kFALSE};
227 
228  // loop on daughter labels
229  for (i = 0; i < ndg; i++){
230  labMom[i] = -1;
231  lab = TMath::Abs(dgLabels[i]);
232  if (lab < 0){
233  printf("daughter with negative label %d\n", lab);
234  return -1;
235  }
236  part = (AliAODMCParticle*)mcArray->At(lab);
237  if (!part){
238  printf("no MC particle\n");
239  return -1;
240  }
241 
242  // check the PDG of the daughter, if requested
243  pdgPart = TMath::Abs(part->GetPdgCode());
244  for (j = 0; j < ndg; j++){
245  if (!pdgUsed[j] && pdgPart == pdgDg[j]){
246  pdgUsed[j] = kTRUE;
247  break;
248  }
249  }
250 
251  mother = part;
252  while (mother->GetMother() >= 0){
253  labMother = mother->GetMother();
254  mother = (AliAODMCParticle*)mcArray->At(labMother);
255  if (!mother){
256  printf("no MC mother particle\n");
257  break;
258  }
259  pdgMother = TMath::Abs(mother->GetPdgCode());
260  if (pdgMother == pdgabs){
261  labMom[i] = labMother;
262  // keep sum of daughters' momenta, to check for mom conservation
263  pxSumDgs += part->Px();
264  pySumDgs += part->Py();
265  pzSumDgs += part->Pz();
266  break;
267  } else break;
268  }
269  if (labMom[i] == -1) return -1; // mother PDG not ok for this daughter
270  } // end loop on daughters
271 
272  // check if the candidate is signal
273  labMother = labMom[0];
274  // all labels have to be the same and !=-1
275  for (i = 0; i < ndg; i++){
276  if (labMom[i] == -1) return -1;
277  if (labMom[i] != labMother) return -1;
278  }
279 
280  // check that all daughter PDGs are matched
281  for (i = 0; i < ndg; i++){
282  if (pdgUsed[i] == kFALSE) return -1;
283  }
284 
285  // check the number of daughters (we are not looking at resonant decay)
286  if(mother->GetNDaughters() != 2) return -1;
287 
288  // Check for mom conservation
289  mother = (AliAODMCParticle*)mcArray->At(labMother);
290  Double_t pxMother = mother->Px();
291  Double_t pyMother = mother->Py();
292  Double_t pzMother = mother->Pz();
293  // within 2% (fix for ITSUpgrade MC's: 0.00001 -> 0.02)
294  if ((TMath::Abs(pxMother - pxSumDgs) / (TMath::Abs(pxMother) + 1.e-13)) > 0.02 &&
295  (TMath::Abs(pyMother - pySumDgs) / (TMath::Abs(pyMother) + 1.e-13)) > 0.02 &&
296  (TMath::Abs(pzMother - pzSumDgs) / (TMath::Abs(pzMother) + 1.e-13)) > 0.02)
297  {
298  AliWarning(Form("Mom. cons. not within 2.0%% perc for decay pdgabs = %d daughters = %d. Returning negative label!",pdgabs,(Int_t)mother->GetNDaughters()));
299  signlabMother = -1;
300  }
301  labMother = signlabMother * labMother;
302  return labMother;
303 }
304 //-----------------------------------------------------------------------------
305 Int_t AliAODRecoDecayHF2Prong::MatchToMCB3Prong(Int_t pdgabs, Int_t pdgabs3prong, Int_t *pdgBDg, Int_t *pdgDg3prong, TClonesArray *mcArray) const
306 {
307  //
308  // Check if this candidate is matched to a MC signal
309  // If no, return -1
310  // If yes, return label (>=0) of the AliAODMCParticle
311  // If yes, but rejected by momentum conservation check, return -1*label (to allow checks at task level)
312  //
313  // NB: This function is only for non-resonant decays of the 2prong (3prong can decay resonant)
314  // NB: Loosened cut on mom. conserv. (needed because of small issue in ITS Upgrade productions)
315  //
316 
317  Int_t signlabMother = 1;
318 
319  // Check number of daughters. Candidate is AliAODRecoDecayHF2Prong, so only continue when 2 daughters
320  Int_t ndg = GetNDaughters();
321  if (!ndg) { AliWarning("HF2Prong: No daughters available"); return -1;}
322  if (ndg != 2) { AliWarning(Form("HF2Prong: %d daughters instead of 2",ndg)); return -1;}
323 
324  // loop on daughters and write the labels
325  Int_t dgLabels[2] = {-1};
326 
327  AliAODTrack *trk = (AliAODTrack*)GetDaughter(1);
328  if(!trk) return -1;
329  dgLabels[1] = trk->GetLabel();
330 
331  AliAODRecoDecayHF3Prong* daug3prong = (AliAODRecoDecayHF3Prong*)GetDaughter(0);
332  if(!daug3prong) return -1;
333  //Similar as daug3prong->MatchToMC(pdgabs3prong,mcArray,3,pdgDg3prong)
334  //but w/o mom conservation check (fix for ITSUpgrade MC's)
335  Int_t ndg3pr = daug3prong->GetNDaughters();
336  if (!ndg3pr) { AliWarning("HF3Prong: No daughters available"); return -1;}
337  if (ndg3pr != 3) { AliWarning(Form("HF3Prong: %d daughters instead of 3",ndg3pr)); return -1;}
338 
339  Int_t dgLabels3pr[3] = {-1};
340  for(Int_t i=0; i<ndg3pr; i++) {
341  AliAODTrack *trk = (AliAODTrack*)daug3prong->GetDaughter(i);
342  dgLabels3pr[i] = trk->GetLabel();
343  }
344 
345  Int_t labMom3pr[3]={0,0,0};
346  Int_t lab3pr, labMother3pr, pdgMother3pr, pdgPart3pr;
347  AliAODMCParticle *part3pr=0;
348  AliAODMCParticle *mother3pr=0;
349  Double_t pxSumDgs3pr=0., pySumDgs3pr=0., pzSumDgs3pr=0.;
350  Bool_t pdgUsed3pr[3]={kFALSE,kFALSE,kFALSE};
351 
352  // loop on daughter labels
353  for(int i=0; i<ndg3pr; i++) {
354  labMom3pr[i]=-1;
355 
356  lab3pr = TMath::Abs(dgLabels3pr[i]);
357  if(lab3pr<0) { AliWarning(Form("daughter with negative label %d",lab3pr)); return -1; }
358 
359  part3pr = (AliAODMCParticle*)mcArray->At(lab3pr);
360  if(!part3pr) { AliWarning("no MC particle"); return -1; }
361 
362  // check the PDG of the daughter
363  pdgPart3pr=TMath::Abs(part3pr->GetPdgCode());
364  for(int j=0; j<ndg3pr; j++) {
365  if(!pdgUsed3pr[j] && pdgPart3pr==pdgDg3prong[j]) {
366  pdgUsed3pr[j]=kTRUE;
367  break;
368  }
369  }
370 
371  mother3pr = part3pr;
372  while(mother3pr->GetMother()>=0) {
373  labMother3pr=mother3pr->GetMother();
374  mother3pr = (AliAODMCParticle*)mcArray->At(labMother3pr);
375  if(!mother3pr) { AliWarning("no MC mother particle"); break; }
376 
377  pdgMother3pr = TMath::Abs(mother3pr->GetPdgCode());
378  if(pdgMother3pr==pdgabs3prong) {
379  labMom3pr[i]=labMother3pr;
380  // keep sum of daughters' momenta, to check for mom conservation
381  pxSumDgs3pr += part3pr->Px();
382  pySumDgs3pr += part3pr->Py();
383  pzSumDgs3pr += part3pr->Pz();
384  break;
385  } else if(pdgMother3pr>pdgabs3prong || pdgMother3pr<10) {
386  break;
387  }
388  }
389  if(labMom3pr[i]==-1) return -1; // mother PDG not ok for this daughter
390  } // end loop on daughters
391 
392  // check if the candidate is signal
393  labMother3pr=labMom3pr[0];
394  // all labels have to be the same and !=-1
395  for(int i=0; i<ndg3pr; i++) {
396  if(labMom3pr[i]==-1) return -1;
397  if(labMom3pr[i]!=labMother3pr) return -1;
398  }
399 
400  // check that all daughter PDGs are matched
401  for (int i = 0; i < ndg3pr; i++){
402  if (pdgUsed3pr[i] == kFALSE) return -1;
403  }
404 
405  mother3pr = (AliAODMCParticle*)mcArray->At(labMother3pr);
406  Double_t pxMother3pr = mother3pr->Px();
407  Double_t pyMother3pr = mother3pr->Py();
408  Double_t pzMother3pr = mother3pr->Pz();
409  // within 2% (fix for ITSUpgrade MC's: 0.00001 -> 0.02)
410  if((TMath::Abs(pxMother3pr-pxSumDgs3pr)/(TMath::Abs(pxMother3pr)+1.e-13)) > 0.02 &&
411  (TMath::Abs(pyMother3pr-pySumDgs3pr)/(TMath::Abs(pyMother3pr)+1.e-13)) > 0.02 &&
412  (TMath::Abs(pzMother3pr-pzSumDgs3pr)/(TMath::Abs(pzMother3pr)+1.e-13)) > 0.02) {
413  AliWarning(Form("Mom. cons. not within 2.0%% perc for decay pdgabs = %d daughters = %d. Returning negative label!",pdgabs3prong,(Int_t)mother3pr->GetNDaughters()));
414  signlabMother = -1;
415  }
416  dgLabels[0] = labMother3pr;
417 
418  Int_t labMom[2] = {0, 0};
419  Int_t lab, labMother, pdgMother, pdgPart;
420  AliAODMCParticle *part = 0;
421  AliAODMCParticle *mother = 0;
422  Double_t pxSumDgs = 0., pySumDgs = 0., pzSumDgs = 0.;
423  Bool_t pdgUsed[2] = {kFALSE, kFALSE};
424 
425  // loop on daughter labels
426  for (int i = 0; i < ndg; i++){
427  labMom[i] = -1;
428 
429  lab = TMath::Abs(dgLabels[i]);
430  if(lab<0) { AliWarning(Form("daughter with negative label %d",lab)); return -1; }
431 
432  part = (AliAODMCParticle*)mcArray->At(lab);
433  if(!part) { AliWarning("no MC particle"); return -1; }
434 
435 
436  // check the PDG of the daughter, if requested
437  pdgPart = TMath::Abs(part->GetPdgCode());
438  for (int j = 0; j < ndg; j++){
439  if (!pdgUsed[j] && pdgPart == pdgBDg[j]){
440  pdgUsed[j] = kTRUE;
441  break;
442  }
443  }
444 
445  mother = part;
446  while (mother->GetMother() >= 0){
447  labMother = mother->GetMother();
448  mother = (AliAODMCParticle*)mcArray->At(labMother);
449  if(!mother) { AliWarning("no MC mother particle"); break; }
450 
451  pdgMother = TMath::Abs(mother->GetPdgCode());
452  if (pdgMother == pdgabs){
453  labMom[i] = labMother;
454  // keep sum of daughters' momenta, to check for mom conservation
455  pxSumDgs += part->Px();
456  pySumDgs += part->Py();
457  pzSumDgs += part->Pz();
458  break;
459  } else break;
460  }
461  if (labMom[i] == -1) return -1; // mother PDG not ok for this daughter
462  } // end loop on daughters
463 
464  // check if the candidate is signal
465  labMother = labMom[0];
466  // all labels have to be the same and !=-1
467  for (int i = 0; i < ndg; i++){
468  if (labMom[i] == -1) return -1;
469  if (labMom[i] != labMother) return -1;
470  }
471 
472  // check that all daughter PDGs are matched
473  for (int i = 0; i < ndg; i++){
474  if (pdgUsed[i] == kFALSE) return -1;
475  }
476 
477  // check the number of daughters (we are not looking at resonant decay)
478  if(mother->GetNDaughters() != 2) return -1;
479 
480  // Check for mom conservation
481  mother = (AliAODMCParticle*)mcArray->At(labMother);
482  Double_t pxMother = mother->Px();
483  Double_t pyMother = mother->Py();
484  Double_t pzMother = mother->Pz();
485  // within 2% (fix for ITSUpgrade MC's: 0.00001 -> 0.02)
486  if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.02 &&
487  (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.02 &&
488  (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.02) {
489  AliWarning(Form("Mom. cons. not within 2.0%% perc for decay pdgabs = %d daughters = %d. Returning negative label!",pdgabs,(Int_t)mother->GetNDaughters()));
490  signlabMother = -1;
491  }
492  labMother = signlabMother * labMother;
493 
494  return labMother;
495 }
AliAODRecoDecayHF & operator=(const AliAODRecoDecayHF &source)
double Double_t
Definition: External.C:58
Int_t MatchToMCB2Prong(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray) const
Double_t CosThetaStarJPSI() const
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Bool_t SelectD0(const Double_t *cuts, Int_t &okD0, Int_t &okD0bar) const
AliAODRecoDecayHF2Prong & operator=(const AliAODRecoDecayHF2Prong &source)
Double_t InvMassJPSIee() const
angle of e-
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Bool_t SelectBtoJPSI(const Double_t *cuts, Int_t &okB) const
Int_t MatchToMCB3Prong(Int_t pdgabs, Int_t pdgabs3prong, Int_t *pdgBDg, Int_t *pdgDg3prong, TClonesArray *mcArray) const