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 //                                                               //
  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     AliError("LambdaC not yet implemented");
 111     return bSignAssoc;
 112   }else{
 113     AliError("WRONG DECAY SETTING");
 114     return bSignAssoc;
 115   }
 116
 117   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);
 118   if (mcLabel == -1) return bSignAssoc;
 119
 120   if (fRecoCandidate->NumberOfFakeDaughters()>0){
 121           fFake = 0;    // fake candidate
 122           if (fFakeSelection==1) return bSignAssoc;
 123   }
 124   if (fRecoCandidate->NumberOfFakeDaughters()==0){
 125           fFake = 2;    // non-fake candidate
 126           if (fFakeSelection==2) return bSignAssoc;
 127   }
 128
 129   SetMCLabel(mcLabel);
 130   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 131
 132   if (!fmcPartCandidate){
 133     AliDebug(3,"No part candidate");
 134     return bSignAssoc;
 135   }
 136
 137   bSignAssoc = kTRUE;
 138   return bSignAssoc;
 139 }
 140
 141 //______________________________________________
 142 Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
 143   //
 144   // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
 145   // pt_D
 146   // y_D
 147   // phi_D
 148   // ctau
 149   // cos point
 150   // pt_1
 151   // pt_2
 152   // pt_3
 153   // d0_1
 154   // d0_2
 155   // d0_3
 156   // zPrimVert
 157   // centrality
 158
 159   Bool_t bGenValues = kFALSE;
 160
 161   Int_t pdgCand = -1;
 162   if(fDecay==kDplustoKpipi){
 163     pdgCand=411;
 164   }else if(fDecay==kDstoKKpi){
 165     pdgCand=431;
 166   }else if(fDecay==kLctopKpi){
 167     AliError("LambdaC not yet implemented");
 168     return bGenValues;
 169   }else{
 170     AliError("WRONG DECAY SETTING");
 171     return bGenValues;
 172   }
 173
 174   Double_t vertD[3] = {0,0,0};   // D origin
 175   fmcPartCandidate->XvYvZv(vertD);  // cm
 176
 177   Int_t nprongs = 3;
 178   Int_t daughter[3];
 179   Short_t charge = fmcPartCandidate->Charge();
 180
 181   // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
 182   // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
 183   Int_t index=0;
 184   Int_t nDauLS=0;
 185   Int_t nDauOS=0;
 186
 187
 188   Int_t nDau=fmcPartCandidate->GetNDaughters();
 189   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 190   if(nDau==3){
 191     for(Int_t iDau=0; iDau<3; iDau++){
 192       Int_t ind = labelFirstDau+iDau;
 193       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 194       if(!part){
 195         AliError("Daughter particle not found in MC array");
 196         return bGenValues;
 197       }
 198       Short_t signDau=part->Charge();
 199       if(signDau==charge){
 200         nDauLS++;
 201         daughter[index] = ind;
 202         index=2;
 203       }else{
 204         daughter[1] = ind;
 205         nDauOS++;
 206       }
 207     }
 208   }else if(nDau==2){
 209     for(Int_t iDau=0; iDau<2; iDau++){
 210       Int_t ind = labelFirstDau+iDau;
 211       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 212       if(!part){
 213         AliError("Daughter particle not found in MC array");
 214         return bGenValues;
 215       }
 216       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 217       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 218         Short_t signDau=part->Charge();
 219         if(signDau==charge){
 220           nDauLS++;
 221           daughter[index] = ind;
 222           index=2;
 223         }else{
 224           daughter[1] = ind;
 225           nDauOS++;
 226         }
 227       }else{
 228         Int_t nDauRes=part->GetNDaughters();
 229         if(nDauRes!=2){
 230           AliError("Wrong resonant decay");
 231           return bGenValues;
 232         }
 233         Int_t labelFirstDauRes = part->GetDaughter(0);
 234         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 235           Int_t indDR = labelFirstDauRes+iDauRes;
 236           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 237           if(!partDR){
 238             AliError("Daughter particle not found in MC array");
 239             return bGenValues;
 240           }
 241           Short_t signDau=partDR->Charge();
 242           if(signDau==charge){
 243             nDauLS++;
 244             daughter[index] = ind;
 245             index=2;
 246           }else{
 247             daughter[1] = ind;
 248             nDauOS++;
 249           }
 250         }
 251       }
 252     }
 253   }else{
 254     AliError(Form("Wrong number of daughters %d",nDau));
 255     return bGenValues;
 256   }
 257
 258   if(nDauLS!=2 || nDauOS!=1){
 259     AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
 260     return bGenValues;
 261   }
 262   if(daughter[0]>daughter[2]){
 263     Int_t tmp=daughter[0];
 264     daughter[0]=daughter[2];
 265     daughter[2]=tmp;
 266   }
 267
 268   // getting the momentum from the daughters and decay vertex
 269   Double_t px[3],py[3],pz[3],pt[3];
 270   Double_t vertDec[3] = {0,0,0};   // decay vertex
 271   for(Int_t iDau=0; iDau<3; iDau++){
 272     AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
 273     if(!part){
 274       AliError("Daughter particle not found in MC array");
 275       return bGenValues;
 276     }
 277     px[iDau]=part->Px();
 278     py[iDau]=part->Py();
 279     pz[iDau]=part->Pz();
 280     pt[iDau]=part->Pt();
 281     if(iDau==0) part->XvYvZv(vertDec);
 282   }
 283
 284   Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
 285
 286   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
 287   Double_t cT = decay->Ct(pdgCand);
 288
 289   vectorMC[0] = fmcPartCandidate->Pt();
 290   vectorMC[1] = fmcPartCandidate->Y() ;
 291   vectorMC[2] = fmcPartCandidate->Phi();
 292   vectorMC[3] = cT*1.E4 ;  // in micron
 293   vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
 294   vectorMC[5] = pt[0];
 295   vectorMC[6] = pt[1];
 296   vectorMC[7] = pt[2];
 297   vectorMC[8] = 0.;   // imppar0, dummy value, meaningless in MC
 298   vectorMC[9] = 0.;   // imppar1, dummy value, meaningless in MC, in micron
 299   vectorMC[10] = 0.;   // imppar2, dummy value, meaningless in MC, in micron
 300   vectorMC[11] = fzMCVertex;    // z of reconstructed of primary vertex
 301   vectorMC[12] = fCentValue; // reconstructed centrality value
 302   vectorMC[13] = 1.;           // always filling with 1 at MC level
 303
 304
 305   bGenValues = kTRUE;
 306   return bGenValues;
 307 }
 308
 309
 310 //____________________________________________
 311 Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 312 {
 313   // Fill vector (see above) with reconstructed quantities
 314   Bool_t bFillRecoValues=kFALSE;
 315
 316   Int_t pdgCand = -1;
 317   if(fDecay==kDplustoKpipi){
 318     pdgCand=411;
 319   }else if(fDecay==kDstoKKpi){
 320     pdgCand=431;
 321   }else if(fDecay==kLctopKpi){
 322     AliError("LambdaC not yet implemented");
 323     return bFillRecoValues;
 324   }else{
 325     AliError("WRONG DECAY SETTING");
 326     return bFillRecoValues;
 327   }
 328
 329   AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
 330   Short_t charge=decay3->Charge();
 331   Double_t rapidity=decay3->Y(pdgCand);
 332   Double_t cT=decay3->Ct(pdgCand);
 333   Double_t pt = decay3->Pt();
 334   Double_t cosPointingAngle = decay3->CosPointingAngle();
 335   Double_t phi = decay3->Phi();
 336
 337   Int_t daughtSorted[3];
 338   Int_t tmpIndex=0;
 339   Int_t nDauLS=0;
 340   Int_t nDauOS=0;
 341   for(Int_t iDau=0; iDau<3; iDau++){
 342     AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
 343     Int_t label = TMath::Abs(trk->GetLabel());
 344     Short_t chargedau=trk->Charge();
 345     if(chargedau==charge){
 346       daughtSorted[tmpIndex]=label;
 347       tmpIndex=2;
 348       nDauLS++;
 349     }else{
 350       daughtSorted[1]=label;
 351       nDauOS++;
 352     }
 353   }
 354
 355   if(nDauLS!=2 || nDauOS!=1){
 356     AliError("Wrong decay channel: number of OS and LS tracks not OK");
 357     return bFillRecoValues;
 358   }
 359
 360   if(daughtSorted[0]>daughtSorted[2]){
 361     Int_t tmp=daughtSorted[0];
 362     daughtSorted[0]=daughtSorted[2];
 363     daughtSorted[2]=tmp;
 364   }
 365
 366
 367   vectorReco[0] = pt;
 368   vectorReco[1] = rapidity;
 369   vectorReco[2] = phi;
 370   vectorReco[3] = cT*1.E4;  // in micron
 371   vectorReco[4] = cosPointingAngle;  // in micron
 372   vectorReco[5] = decay3->PtProng(daughtSorted[0]);
 373   vectorReco[6] = decay3->PtProng(daughtSorted[1]);
 374   vectorReco[7] = decay3->PtProng(daughtSorted[2]);
 375   vectorReco[8] = decay3->Getd0Prong(daughtSorted[0]);
 376   vectorReco[9] = decay3->Getd0Prong(daughtSorted[1]);
 377   vectorReco[10] = decay3->Getd0Prong(daughtSorted[2]);
 378   vectorReco[11] = fzPrimVertex;    // z of reconstructed of primary vertex
 379   vectorReco[12] = fCentValue; //reconstructed centrality value
 380   vectorReco[13] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
 381
 382
 383   bFillRecoValues = kTRUE;
 384   return bFillRecoValues;
 385 }
 386
 387
 388 //_____________________________________________________________
 389 Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
 390 {
 391   // Check the pdg codes of the daughters
 392   Bool_t checkCD = kFALSE;
 393
 394   Int_t pdgCand = -1;
 395   Int_t pdgDaughter[3]={-1,-1,-1};
 396   if(fDecay==kDplustoKpipi){
 397     pdgCand=411;
 398     pdgDaughter[0]=321;
 399     pdgDaughter[1]=211;
 400     pdgDaughter[2]=211;
 401   }else if(fDecay==kDstoKKpi){
 402     pdgCand=431;
 403     pdgDaughter[0]=321;
 404     pdgDaughter[1]=321;
 405     pdgDaughter[2]=211;
 406   }else if(fDecay==kLctopKpi){
 407     AliError("LambdaC not yet implemented");
 408     return checkCD;
 409   }else{
 410     AliError("WRONG DECAY SETTING");
 411     return checkCD;
 412   }
 413
 414
 415   Int_t daughter[3];
 416
 417   Int_t nDau=fmcPartCandidate->GetNDaughters();
 418   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 419   if(nDau==3){
 420     for(Int_t iDau=0; iDau<3; iDau++){
 421       Int_t ind = labelFirstDau+iDau;
 422       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 423       if(!part){
 424         AliError("Daughter particle not found in MC array");
 425         return checkCD;
 426       }
 427       daughter[iDau]=TMath::Abs(part->GetPdgCode());
 428     }
 429   }else if(nDau==2){
 430     Int_t nDauFound=0;
 431     for(Int_t iDau=0; iDau<2; iDau++){
 432       Int_t ind = labelFirstDau+iDau;
 433       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 434       if(!part){
 435         AliError("Daughter particle not found in MC array");
 436         return checkCD;
 437       }
 438       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 439       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 440         if(nDauFound>=3) return checkCD;
 441         daughter[nDauFound]=pdgCode;
 442         nDauFound++;
 443       }else{
 444         Int_t nDauRes=part->GetNDaughters();
 445         if(nDauRes!=2) return checkCD;
 446         Int_t labelFirstDauRes = part->GetDaughter(0);
 447         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 448           Int_t indDR = labelFirstDauRes+iDauRes;
 449           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 450           if(!partDR){
 451             AliError("Daughter particle not found in MC array");
 452             return checkCD;
 453           }
 454           Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 455           if(nDauFound>=3) return checkCD;
 456           daughter[nDauFound]=pdgCodeDR;
 457           nDauFound++;
 458         }
 459       }
 460     }
 461   }else{
 462     return checkCD;
 463   }
 464   for(Int_t iDau1=0; iDau1<3; iDau1++){
 465     for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
 466       if(daughter[iDau1]<daughter[iDau2]){
 467         Int_t tmp=daughter[iDau1];
 468         daughter[iDau1]=daughter[iDau2];
 469         daughter[iDau2]=tmp;
 470       }
 471     }
 472   }
 473   for(Int_t iDau=0; iDau<3; iDau++){
 474     if(daughter[iDau]!=pdgDaughter[iDau]){
 475       AliDebug(2, "Wrong decay channel from MC, skipping!!");
 476       return checkCD;
 477     }
 478   }
 479
 480   checkCD = kTRUE;
 481   return checkCD;
 482
 483 }