AliPhysics  9b6b435 (9b6b435)
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  //
191  // NB: This function is only for non-resonant decays
192  // NB: No cut on mom conservation due to long decay time
193 
194  // Check number of daughters. Candidate is AliAODRecoDecayHF2Prong, so only continue when 2 daughters
195  Int_t ndg = GetNDaughters();
196  if (!ndg) { AliError("HF2Prong: No daughters available"); return -1;}
197  if (ndg != 2) { AliError(Form("HF2Prong: %d daughters instead of 2",ndg)); return -1;}
198 
199  // loop on daughters and write the labels
200  Int_t dgLabels[2] = {-1};
201  if(pdgabs2prong == 0){
202  for(Int_t i=0; i<ndg; i++) {
203  AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
204  dgLabels[i] = trk->GetLabel();
205  }
206  } else {
207  AliAODTrack *trk = (AliAODTrack*)GetDaughter(0);
208  dgLabels[0] = trk->GetLabel();
209  AliAODRecoDecayHF2Prong* daug2prong = (AliAODRecoDecayHF2Prong*)GetDaughter(1);
210  //Daughter prong (independent on their decay time) will also be far from PV, so call this function again.
211  Int_t pdgempty[2]={0,0};
212  dgLabels[1] = daug2prong->MatchToMCB2Prong(pdgabs2prong, 0, pdgDg2prong, pdgempty, mcArray);
213  }
214  if (dgLabels[0] == -1) return -1;
215  if (dgLabels[1] == -1) return -1;
216 
217  Int_t labMom[2] = {0, 0};
218  Int_t i, j, lab, labMother, pdgMother, pdgPart;
219  AliAODMCParticle *part = 0;
220  AliAODMCParticle *mother = 0;
221  Double_t pxSumDgs = 0., pySumDgs = 0., pzSumDgs = 0.;
222  Bool_t pdgUsed[2] = {kFALSE, kFALSE};
223 
224  // loop on daughter labels
225  for (i = 0; i < ndg; i++){
226  labMom[i] = -1;
227  lab = TMath::Abs(dgLabels[i]);
228  if (lab < 0){
229  printf("daughter with negative label %d\n", lab);
230  return -1;
231  }
232  part = (AliAODMCParticle*)mcArray->At(lab);
233  if (!part){
234  printf("no MC particle\n");
235  return -1;
236  }
237 
238  // check the PDG of the daughter, if requested
239  pdgPart = TMath::Abs(part->GetPdgCode());
240  for (j = 0; j < ndg; j++){
241  if (!pdgUsed[j] && pdgPart == pdgDg[j]){
242  pdgUsed[j] = kTRUE;
243  break;
244  }
245  }
246 
247  mother = part;
248  while (mother->GetMother() >= 0){
249  labMother = mother->GetMother();
250  mother = (AliAODMCParticle*)mcArray->At(labMother);
251  if (!mother){
252  printf("no MC mother particle\n");
253  break;
254  }
255  pdgMother = TMath::Abs(mother->GetPdgCode());
256  if (pdgMother == pdgabs){
257  labMom[i] = labMother;
258  // keep sum of daughters' momenta, to check for mom conservation
259  pxSumDgs += part->Px();
260  pySumDgs += part->Py();
261  pzSumDgs += part->Pz();
262  break;
263  } else break;
264  }
265  if (labMom[i] == -1) return -1; // mother PDG not ok for this daughter
266  } // end loop on daughters
267 
268  // check if the candidate is signal
269  labMother = labMom[0];
270  // all labels have to be the same and !=-1
271  for (i = 0; i < ndg; i++){
272  if (labMom[i] == -1) return -1;
273  if (labMom[i] != labMother) return -1;
274  }
275 
276  // check that all daughter PDGs are matched
277  for (i = 0; i < ndg; i++){
278  if (pdgUsed[i] == kFALSE) return -1;
279  }
280 
281  // check the number of daughters (we are not looking at resonant decay)
282  if(mother->GetNDaughters() != 2) return -1;
283 
284  // Check for mom conservation
285  mother = (AliAODMCParticle*)mcArray->At(labMother);
286  Double_t pxMother = mother->Px();
287  Double_t pyMother = mother->Py();
288  Double_t pzMother = mother->Pz();
289  if ((TMath::Abs(pxMother - pxSumDgs) / (TMath::Abs(pxMother) + 1.e-13)) > 0.005 ||
290  (TMath::Abs(pyMother - pySumDgs) / (TMath::Abs(pyMother) + 1.e-13)) > 0.005 ||
291  (TMath::Abs(pzMother - pzSumDgs) / (TMath::Abs(pzMother) + 1.e-13)) > 0.005)
292  {
293  // Only show warning if mom conservation is not within 0.5%.
294  // This can be due to large propagation distance through magnetic field.
295  AliWarning(Form("Mom. cons. not within 0.5%% perc for decay pdgabs = %d daughters = %d",pdgabs,(Int_t)mother->GetNDaughters()));
296  }
297  return labMother;
298 }
299 //-----------------------------------------------------------------------------
300 Int_t AliAODRecoDecayHF2Prong::MatchToMCB3Prong(Int_t pdgabs,Int_t pdgabs3prong, Int_t *pdgBDg,Int_t *pdgDg3prong, TClonesArray *mcArray) const
301 {
302  //std::cout<<"MCmatch 1"<<std::endl;
303  //
304  // Check if this candidate is matched to a MC signal
305  // If no, return -1
306  // If yes, return label (>=0) of the AliAODMCParticle
307  //
308  Int_t ndg=GetNDaughters();
309  if(ndg==0) {
310  AliError("No daughters available");
311  return -1;
312  }
313  if(ndg>10){
314  AliError("Only decays with <10 daughters supported");
315  return -1;
316  }
317 
318 
319 
320 
321  AliAODRecoDecayHF3Prong *DDaughter = (AliAODRecoDecayHF3Prong*)GetDaughter(0);
322  if(!DDaughter)return -1;
323  Int_t DLabel = DDaughter->MatchToMC(pdgabs3prong,mcArray,3,pdgDg3prong);
324  if(DLabel<0) return -1;
325 
326  Int_t dgLabels[10]={0};
327 
328 
329  // loop on daughters and write labels
330  for(Int_t i=0; i<ndg; i++) {
331  AliVTrack *trk = (AliVTrack*)GetDaughter(i);
332  Int_t lab = trk->GetLabel();
333  if(lab==-1) { // this daughter is the 3prong
334  lab=DLabel;
335  } else if(lab<-1) {
336  printf("daughter with negative label\n");
337  continue;
338  }
339  dgLabels[i] = lab;
340 
341  }
342 
343 
344  Int_t label = MatchToMC(pdgabs,mcArray,dgLabels,ndg,2,pdgBDg);
345 
346 
347 
348  return label;
349 
350 }
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