AliPhysics  63e47e1 (63e47e1)
AliHFCutOptTreeHandler.cxx
Go to the documentation of this file.
1 /* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
2  * See cxx source for full Copyright notice */
3 
4 /* $Id$ */
5 
6 //*************************************************************************
7 // Class AliHFCutOptTreeHandler
8 // helper class to handle a tree for cut optimisation and MVA analyses
9 // Authors:
10 // F. Catalano, fabio.catalano@cern.ch
11 // A. Festanti, andrea.festanti@cern.ch
12 // F. Grosa, fabrizio.grosa@cern.ch
13 // G. Innocenti, gian.michele.innocenti@cern.ch
14 // F. Prino, prino@to.infn.it
15 // L. Vermunt, luuk.vermunt@cern.ch
17 
18 #include "AliHFCutOptTreeHandler.h"
19 
20 #include <TDatabasePDG.h>
21 #include <TMath.h>
22 #include "AliAODTrack.h"
23 #include "AliAODRecoDecayHF.h"
26 #include "AliAODRecoCascadeHF.h"
27 #include "AliAODMCParticle.h"
28 #include "AliVertexingHFUtils.h"
29 
31 ClassImp(AliHFCutOptTreeHandler);
33 
34 //________________________________________________________________
36  TObject(),
37  fTreeTopolVar(0x0),
38  fDecayChannel(kD0toKpi),
39  fPidOpt(kNsigmaPID),
40  fCandType(kBkg),
41  fUseCentrality(false),
42  fCentrality(0),
43  fFillOnlySignal(false),
44  fIsMC(false),
45  fIsSignal(-999),
46  fIsPrompt(-999),
47  fIsRefl(-999),
48  fIsSelStd(false),
49  fUseSelFlag(true)
50 {
51  // Default constructor
52  SetPdgCodes();
53 
54  for(int iVar=0; iVar<knTopolVars; iVar++) fTopolVarVector[iVar] = -999.;
55  for(int iVar=0; iVar<knPidVars; iVar++) {
56  fPIDnSigmaVector[iVar]=-1;
57  fPIDnSigmaCharVector[iVar]=-1;
58  }
59 }
60 
61 //________________________________________________________________
63  TObject(),
64  fTreeTopolVar(0x0),
65  fDecayChannel(decay),
66  fPidOpt(PIDopt),
67  fCandType(kBkg),
68  fUseCentrality(false),
69  fCentrality(0),
70  fFillOnlySignal(false),
71  fIsMC(isMC),
72  fIsSignal(-999),
73  fIsPrompt(-999),
74  fIsRefl(-999),
75  fIsSelStd(false),
76  fUseSelFlag(true)
77 {
78  // Standard constructor
79  SetPdgCodes();
80  for(int iVar=0; iVar<knTopolVars; iVar++) fTopolVarVector[iVar] = -999.;
81  for(int iVar=0; iVar<knPidVars; iVar++) {
82  fPIDnSigmaVector[iVar]=-1;
83  fPIDnSigmaCharVector[iVar]=-1;
84  }
85 }
86 
87 //________________________________________________________________
89 {
90  // Destructor
91  if(fTreeTopolVar) delete fTreeTopolVar;
92 }
93 
94 //________________________________________________________________
95 bool AliHFCutOptTreeHandler::SetVariables(AliAODRecoDecayHF* d, int masshypo, AliAODPidHF* pidHF, TClonesArray* arrayMC)
96 {
97 
98  if(!d) return false;
99 
100  //candidate type
101  if(fIsMC) {
102  if(fIsSignal==-999 || fIsPrompt==-999 || fIsRefl==-999) { // if not setted, recompute
103  if(!arrayMC) fIsSignal=0;
104  else {
105  int labD=-1;
106  int orig=-1;
107  int labDau0=-1;
108  int pdgdau0=-1;
109  bool isrefl=false;
110  if(fDecayChannel!=kD0toKpi) labD = d->MatchToMC(fPdgCode,arrayMC,3,fPdgCodeProngs);
111  else labD = d->MatchToMC(fPdgCode,arrayMC,2,fPdgCodeProngs);
112  if(labD<0) fIsSignal=0;
113  else {
114  fIsSignal=1;
115 
116  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
117  if(partD) orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partD,kTRUE); //4 --> prompt, 5 --> FD
118 
119  if(orig==4) fIsPrompt=1;
120  else if(orig==5) fIsPrompt=0;
121  else fIsSignal=0;
122 
123  if(fDecayChannel==kDplustoKpipi) fIsRefl=0; // no reflected signal for D+ -> Kpipi
124  else {
125  labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
126  AliAODMCParticle* dau0=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
127  if(dau0) pdgdau0=TMath::Abs(dau0->GetPdgCode());
128  if((masshypo==0 && pdgdau0!=321) || (masshypo==1 && pdgdau0!=211)) fIsRefl=1;
129  else fIsRefl=0;
130  }
131  }
132  }
133  }
134  if(fIsSignal==0) fCandType=kBkg;
135  else {
136  if(fIsPrompt==1) {
138  else fCandType=kPromptRefl;
139  }
140  else {
141  if(fIsRefl==0) fCandType=kFDSig;
142  else fCandType=kFDRefl;
143  }
144  }
145  }
146  else { //if data not possible to tag
147  fCandType=kBkg;
148  }
149 
150  if(fCandType==kBkg && fFillOnlySignal) return true;
151 
152  //reset candidate flags
153  fIsSignal=-999;
154  fIsPrompt=-999;
155  fIsRefl=-999;
156 
157  //topological variables
158  fTopolVarVector[1]=d->Pt();
164  fTopolVarVector[7]=d->ImpParXY();
165  fTopolVarVector[8]=d->PtProng(0);
166  fTopolVarVector[9]=d->PtProng(1);
167 
168  switch(fDecayChannel) {
169  case 0: //D0 -> Kpi
170  fTopolVarVector[11]=((AliAODRecoDecayHF2Prong*)d)->Getd0Prong(0);
171  fTopolVarVector[12]=((AliAODRecoDecayHF2Prong*)d)->Getd0Prong(1);
173  if(masshypo==0) {
174  fTopolVarVector[0]=((AliAODRecoDecayHF2Prong*)d)->InvMassD0();
175  fTopolVarVector[10]=((AliAODRecoDecayHF2Prong*)d)->CosThetaStarD0();
176  }
177  else {
178  fTopolVarVector[0]=((AliAODRecoDecayHF2Prong*)d)->InvMassD0bar();
179  fTopolVarVector[10]=((AliAODRecoDecayHF2Prong*)d)->CosThetaStarD0bar();
180  }
181  break;
182  case 1: //D+ -> Kpipi
183  fTopolVarVector[0]=((AliAODRecoDecayHF3Prong*)d)->InvMassDplus();
184  fTopolVarVector[10]=((AliAODRecoDecayHF3Prong*)d)->PtProng(2);
185  fTopolVarVector[11]=((AliAODRecoDecayHF3Prong*)d)->GetSigmaVert();
186  break;
187  case 2: //Ds+ -> KKpi
188  fTopolVarVector[10]=((AliAODRecoDecayHF3Prong*)d)->PtProng(2);
189  fTopolVarVector[11]=((AliAODRecoDecayHF3Prong*)d)->GetSigmaVert();
190  float massPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
191  if(masshypo==0){ //phiKKpi
192  fTopolVarVector[0]=((AliAODRecoDecayHF3Prong*)d)->InvMassDsKKpi();
193  fTopolVarVector[12]=TMath::Abs(((AliAODRecoDecayHF3Prong*)d)->InvMass2Prongs(0,1,321,321)-massPhi);
194  fTopolVarVector[13]=((AliAODRecoDecayHF3Prong*)d)->CosPiDsLabFrameKKpi();
195  fTopolVarVector[14]=((AliAODRecoDecayHF3Prong*)d)->CosPiKPhiRFrameKKpi();
196  }
197  if(masshypo==1){ //phipiKK
198  fTopolVarVector[0]=((AliAODRecoDecayHF3Prong*)d)->InvMassDspiKK();
199  fTopolVarVector[12]=TMath::Abs(((AliAODRecoDecayHF3Prong*)d)->InvMass2Prongs(1,2,321,321)-massPhi);
200  fTopolVarVector[13]=((AliAODRecoDecayHF3Prong*)d)->CosPiDsLabFramepiKK();
201  fTopolVarVector[14]=((AliAODRecoDecayHF3Prong*)d)->CosPiKPhiRFramepiKK();
202  }
203  fTopolVarVector[14] = TMath::Abs(fTopolVarVector[14]*fTopolVarVector[14]*fTopolVarVector[14]);
204  break;
205  }
206 
207  if(fPidOpt==kNoPID || !pidHF) return true; //if no PID, return before
208 
209  SetPidVars(d,pidHF);
210  return true;
211 }
212 
213 //________________________________________________________________
215 {
216  if(fTreeTopolVar) delete fTreeTopolVar;
217  fTreeTopolVar=0x0;
218  fTreeTopolVar = new TTree(name.Data(),title.Data());
219 
220  TString topolvarnamesCommon[knTopolVarsCommon] = {"inv_mass","pt_cand","d_len","d_len_xy","norm_dl_xy","cos_p","cos_p_xy","imp_par_xy","pt_prong0","pt_prong1"};
221  TString topolvarNamesDzero[knTopolVarsDzero] = {"cos_t_star", "imp_par_prong0", "imp_par_prong1","imp_par_prod"};
222  TString topolvarNamesDs[knTopolVarsDs] = {"pt_prong2","sig_vert","delta_mass_KK","cos_PiDs","cos_PiKPhi_3"};
223  TString topolvarNamesDplus[knTopolVarsDplus] = {"pt_prong2","sig_vert"};
224  TString PIDvarnamesNsigma[knPidVars] = {"nsigTPC_Pi_0","nsigTPC_K_0","nsigTOF_Pi_0","nsigTOF_K_0","nsigTPC_Pi_1","nsigTPC_K_1","nsigTOF_Pi_1","nsigTOF_K_1","nsigTPC_Pi_2","nsigTPC_K_2","nsigTOF_Pi_2","nsigTOF_K_2"};
225  TString PIDvarnamesNsigmaComb[knPidVars] = {"nsigComb_Pi_0","nsigComb_K_0","nsigComb_Pi_1","nsigComb_K_1","nsigComb_Pi_2","nsigComb_K_2","","","","","",""};
226 
227  for(int iVar=0; iVar<knTopolVarsCommon; iVar++) {
228  fTreeTopolVar->Branch(topolvarnamesCommon[iVar].Data(),&fTopolVarVector[iVar],Form("%s/F",topolvarnamesCommon[iVar].Data()));
229  }
230  switch(fDecayChannel){
231  case 0: //D0 -> Kpi
232  for(int iVar=0; iVar<knTopolVarsDzero; iVar++){
233  fTreeTopolVar->Branch(topolvarNamesDzero[iVar].Data(),&fTopolVarVector[knTopolVarsCommon+iVar],Form("%s/F",topolvarNamesDzero[iVar].Data()));
234  }
235  break;
236  case 1: //D+ -> Kpipi
237  for(int iVar=0; iVar<knTopolVarsDplus; iVar++){
238  fTreeTopolVar->Branch(topolvarNamesDplus[iVar].Data(),&fTopolVarVector[knTopolVarsCommon+iVar],Form("%s/F",topolvarNamesDplus[iVar].Data()));
239  }
240  break;
241  case 2: //Ds -> KKpi
242  for(int iVar=0; iVar<knTopolVarsDs; iVar++){
243  fTreeTopolVar->Branch(topolvarNamesDs[iVar].Data(),&fTopolVarVector[knTopolVarsCommon+iVar],Form("%s/F",topolvarNamesDs[iVar].Data()));
244  }
245  break;
246  }
247 
248  fTreeTopolVar->Branch("cand_type",&fCandType,"cand_type/B");
249  if(fUseCentrality) fTreeTopolVar->Branch("centrality",&fCentrality,"centrality/B");
250  if(fUseSelFlag) fTreeTopolVar->Branch("isselectedstd",&fIsSelStd,"isselectedstd/O");
251 
252  switch(fPidOpt) {
253  case 0: //no PID
254  return fTreeTopolVar;
255  break;
256  case 1: // nsigmaPID
257  for(int iVar=0; iVar<knPidVars; iVar++) {
258  if(iVar>7 && fDecayChannel==kD0toKpi) continue;
259  fTreeTopolVar->Branch(PIDvarnamesNsigma[iVar].Data(),&fPIDnSigmaVector[iVar],Form("%s/F",PIDvarnamesNsigma[iVar].Data()));
260  }
261  break;
262  case 2: //nsigmaPIDchar
263  for(int iVar=0; iVar<knPidVars; iVar++) {
264  if(iVar>7 && fDecayChannel==kD0toKpi) continue;
265  fTreeTopolVar->Branch(PIDvarnamesNsigma[iVar].Data(),&fPIDnSigmaCharVector[iVar],Form("%s/B",PIDvarnamesNsigma[iVar].Data()));
266  }
267  break;
268  case 3: //nsigmaPIDfloatandchar
269  for(int iVar=0; iVar<knPidVars; iVar++) {
270  if(iVar>7 && fDecayChannel==kD0toKpi) continue;
271  fTreeTopolVar->Branch(PIDvarnamesNsigma[iVar].Data(),&fPIDnSigmaVector[iVar],Form("%s/F",PIDvarnamesNsigma[iVar].Data()));
272  fTreeTopolVar->Branch(Form("%s_char",PIDvarnamesNsigma[iVar].Data()),&fPIDnSigmaCharVector[iVar],Form("%s_char/B",PIDvarnamesNsigma[iVar].Data()));
273  }
274  break;
275  case 4: // nsigmaCombPID
276  for(int iVar=0; iVar<knPidVars/2; iVar++) {
277  if(iVar>3 && fDecayChannel==kD0toKpi) continue;
278  fTreeTopolVar->Branch(PIDvarnamesNsigmaComb[iVar].Data(),&fPIDnSigmaVector[iVar],Form("%s/F",PIDvarnamesNsigmaComb[iVar].Data()));
279  }
280  break;
281  case 5: //nsigmaCombPIDchar
282  for(int iVar=0; iVar<knPidVars/2; iVar++) {
283  if(iVar>3 && fDecayChannel==kD0toKpi) continue;
284  fTreeTopolVar->Branch(PIDvarnamesNsigmaComb[iVar].Data(),&fPIDnSigmaCharVector[iVar],Form("%s/B",PIDvarnamesNsigmaComb[iVar].Data()));
285  }
286  break;
287  case 6: //nsigmaCombPIDfloatandchar
288  for(int iVar=0; iVar<knPidVars/2; iVar++) {
289  if(iVar>3 && fDecayChannel==kD0toKpi) continue;
290  fTreeTopolVar->Branch(PIDvarnamesNsigmaComb[iVar].Data(),&fPIDnSigmaVector[iVar],Form("%s/F",PIDvarnamesNsigmaComb[iVar].Data()));
291  fTreeTopolVar->Branch(Form("%s_char",PIDvarnamesNsigmaComb[iVar].Data()),&fPIDnSigmaCharVector[iVar],Form("%s_char/B",PIDvarnamesNsigmaComb[iVar].Data()));
292  }
293  break;
294  default: //no PID
295  return fTreeTopolVar;
296  break;
297  }
298 
299  return fTreeTopolVar;
300 }
301 
302 //________________________________________________________________
304  //PID variables
305  double sigTPC_K[knMaxProngs]={-999.,-999.,-999.};
306  double sigTPC_Pi[knMaxProngs]={-999.,-999.,-999.};
307  double sigTOF_K[knMaxProngs]={-999.,-999.,-999.};
308  double sigTOF_Pi[knMaxProngs]={-999.,-999.,-999.};
309 
310  float sigComb_K[knMaxProngs] = {-1.,-1.,-1.};
311  float sigComb_Pi[knMaxProngs] = {-1.,-1.,-1.};
312 
313  AliAODTrack *track[knMaxProngs] = {0x0,0x0,0x0};
314 
315  for(int iProng=0; iProng<knMaxProngs; iProng++) {
316  track[iProng]=(AliAODTrack*)d->GetDaughter(0);
317 
318  pidHF->GetnSigmaTPC(track[iProng],3,sigTPC_K[iProng]);
319  pidHF->GetnSigmaTPC(track[iProng],2,sigTPC_Pi[iProng]);
320 
321  pidHF->GetnSigmaTOF(track[iProng],3,sigTOF_K[iProng]);
322  pidHF->GetnSigmaTOF(track[iProng],2,sigTPC_Pi[iProng]);
323 
324  if(fDecayChannel==kD0toKpi && iProng>1) continue; //D0 -> Kpi only 2 prongs
325 
327  fPIDnSigmaVector[4*iProng]=sigTPC_Pi[iProng]*10;
328  fPIDnSigmaVector[4*iProng+1]=sigTPC_K[iProng]*10;
329  fPIDnSigmaVector[4*iProng+2]=sigTPC_Pi[iProng]*10;
330  fPIDnSigmaVector[4*iProng+3]=sigTOF_K[iProng]*10;
331 
332  //to be changed below
333  if(sigTPC_Pi[iProng]>0) {
334  if(sigTPC_Pi[iProng]<12.7) fPIDnSigmaCharVector[4*iProng]=sigTPC_Pi[iProng]*10;
335  else fPIDnSigmaCharVector[4*iProng]=127;
336  }
337  else {
338  if(sigTPC_Pi[iProng]>-12.7) fPIDnSigmaCharVector[4*iProng]=sigTPC_Pi[iProng]*10;
339  else fPIDnSigmaCharVector[4*iProng]=-127;
340  }
341  if(sigTPC_K[iProng]>0) {
342  if(sigTPC_K[iProng]<12.7) fPIDnSigmaCharVector[4*iProng+1]=sigTPC_K[iProng]*10;
343  else fPIDnSigmaCharVector[4*iProng+1]=127;
344  }
345  else {
346  if(sigTPC_K[iProng]>-12.7) fPIDnSigmaCharVector[4*iProng+1]=sigTPC_K[iProng]*10;
347  else fPIDnSigmaCharVector[4*iProng+1]=-127;
348  }
349  if(sigTOF_Pi[iProng]>0) {
350  if(sigTOF_Pi[iProng]<12.7) fPIDnSigmaCharVector[4*iProng+2]=sigTOF_Pi[iProng]*10;
351  else fPIDnSigmaCharVector[4*iProng+2]=127;
352  }
353  else {
354  if(sigTOF_Pi[iProng]>-12.7) fPIDnSigmaCharVector[4*iProng+2]=sigTOF_Pi[iProng]*10;
355  else fPIDnSigmaCharVector[4*iProng+2]=-127;
356  }
357  if(sigTOF_K[iProng]>0) {
358  if(sigTOF_K[iProng]<12.7) fPIDnSigmaCharVector[4*iProng+3]=sigTOF_K[iProng]*10;
359  else fPIDnSigmaCharVector[4*iProng+3]=127;
360  }
361  else {
362  if(sigTOF_K[iProng]>-12.7) fPIDnSigmaCharVector[4*iProng+3]=sigTOF_K[iProng]*10;
363  else fPIDnSigmaCharVector[4*iProng+3]=-127;
364  }
365  }
367 
368  //PID TPC and TOF squared sum, normalized to sqrt(2)
369  if(sigTPC_K[iProng] > -998.){
370  if(sigTOF_K[iProng] > -998.)
371  sigComb_K[iProng] = TMath::Sqrt((sigTPC_K[iProng]*sigTPC_K[iProng] + sigTOF_K[iProng]*sigTOF_K[iProng]) / 2.);
372  else
373  sigComb_K[iProng] = TMath::Abs(sigTPC_K[iProng]);
374  }
375 
376  if(sigTPC_Pi[iProng] > -998.){
377  if(sigTOF_Pi[iProng] > -998.)
378  sigComb_Pi[iProng] = TMath::Sqrt((sigTPC_Pi[iProng]*sigTPC_Pi[iProng] + sigTOF_Pi[iProng]*sigTOF_Pi[iProng]) / 2.);
379  else
380  sigComb_Pi[iProng] = TMath::Abs(sigTPC_Pi[iProng]);
381  }
382 
383  fPIDnSigmaVector[2*iProng]=sigComb_Pi[iProng]*10;
384  fPIDnSigmaVector[2*iProng+1]=sigComb_K[iProng]*10;
385 
386  sigComb_Pi[iProng]<12.7 ? fPIDnSigmaCharVector[2*iProng]=sigComb_Pi[iProng]*10 : 127; //set last bin
387  sigComb_K[iProng]<12.7 ? fPIDnSigmaCharVector[2*iProng+1]=sigComb_K[iProng]*10 : 127; //set last bin
388  }
389  }
390 }
391 
392 //________________________________________________________________
394 
395  switch(fDecayChannel) {
396  case 0:
397  fPdgCode=421;
398  fPdgCodeProngs[0]=321;
399  fPdgCodeProngs[1]=211;
400  fPdgCodeProngs[2]=-1;
401  break;
402  case 1:
403  fPdgCode=411;
404  fPdgCodeProngs[0]=321;
405  fPdgCodeProngs[1]=211;
406  fPdgCodeProngs[2]=211;
407  break;
408  case 2:
409  fPdgCode=431;
410  fPdgCodeProngs[0]=321;
411  fPdgCodeProngs[1]=321;
412  fPdgCodeProngs[2]=211;
413  break;
414  default:
415  fPdgCode=-1;
416  fPdgCodeProngs[0]=-1;
417  fPdgCodeProngs[1]=-1;
418  fPdgCodeProngs[2]=-1;
419  break;
420  }
421 }
Double_t NormalizedDecayLengthXY() const
int fDecayChannel
tree for cut optimisation
const char * title
Definition: MakeQAPdf.C:27
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
bool SetVariables(AliAODRecoDecayHF *d, int masshypo, AliAODPidHF *pidHF=0x0, TClonesArray *arrayMC=0x0)
Double_t ImpParXY() const
int fPidOpt
array with nSigma PID (char)
Double_t CosPointingAngleXY() const
bool fIsMC
flag to enable only signal filling
int fIsPrompt
flag for signal=1 (including prompt, FD, reflected), bkg=0
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
char fIsSelStd
flag for reflected signal=1, non-reflected signal=0
bool fUseSelFlag
flag to tag selected candidates by "standard" cuts
void SetPidVars(AliAODRecoDecayHF *d, AliAODPidHF *pidHF)
int fIsSignal
flag to enable checks on MC truth
float fTopolVarVector[knTopolVars]
absolute values of pdg codes of the daughters
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
int fIsRefl
flag for prompt=1 (inluding reflected), FD=0
Double_t DecayLengthXY() const
Bool_t isMC
decay
Definition: HFPtSpectrum.C:41
char fCentrality
flag to enable centrality
bool fFillOnlySignal
centrality in case of p-Pb or Pb-Pb
Double_t CosPointingAngle() const
float fPIDnSigmaCharVector[knPidVars]
array with nSigma PID
char fCandType
option for PID variables
float fPIDnSigmaVector[knPidVars]
array with topological variables
bool fUseCentrality
flag for candidate type (bkg, prompt signal, FD signal, prompt refl, FD refl)
Double_t DecayLength() const
TTree * BuildTree(TString name="tree", TString title="tree")
int fPdgCodeProngs[knMaxProngs]
absolute value of pdg code of the particle of interest