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 //_________________________________________
  50 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
  51   AliCFVertexingHF(mcArray, originDselection),
  52   fDecay(decay)
  53
  54 {
  55   //
  56   SetNProngs(3);
  57 }
  58
  59
  60 //_____________________________________
  61 AliCFVertexingHF3Prong& AliCFVertexingHF3Prong::operator=(const AliCFVertexingHF3Prong& c){
  62   //
  63   if  (this != &c) {
  64
  65     AliCFVertexingHF::operator=(c);
  66
  67   }
  68     return *this;
  69 }
  70
  71 //__________________________________________
  72 Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
  73   // Checks if candidate is signal and D meson is present in MC array
  74
  75   Bool_t bSignAssoc = kFALSE;
  76   fRecoCandidate = recoCand;
  77
  78   if (!fRecoCandidate) {
  79     AliError("fRecoCandidate not found, problem in assignement\n");
  80     return bSignAssoc;
  81   }
  82
  83   Int_t pdgCand = -1;
  84   Int_t pdgDaughter[3]={-1,-1,-1};
  85   if(fDecay==kDplustoKpipi){
  86     pdgCand=411;
  87     pdgDaughter[0]=321;
  88     pdgDaughter[1]=211;
  89     pdgDaughter[2]=211;
  90   }else if(fDecay==kDstoKKpi){
  91     pdgCand=431;
  92     pdgDaughter[0]=321;
  93     pdgDaughter[1]=321;
  94     pdgDaughter[2]=211;
  95   }else if(fDecay==kLctopKpi){
  96     AliError("LambdaC not yet implemented");
  97     return bSignAssoc;
  98   }else{
  99     AliError("WRONG DECAY SETTING");
 100     return bSignAssoc;
 101   }
 102
 103   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);
 104   if (mcLabel == -1) return bSignAssoc;
 105   SetMCLabel(mcLabel);
 106   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 107
 108   if (!fmcPartCandidate){
 109     AliDebug(3,"No part candidate");
 110     return bSignAssoc;
 111   }
 112
 113   bSignAssoc = kTRUE;
 114   return bSignAssoc;
 115 }
 116
 117 //______________________________________________
 118 Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
 119   //
 120   // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
 121   // pt_D
 122   // y_D
 123   // phi_D
 124   // ctau
 125   // cos point
 126   // pt_1
 127   // pt_2
 128   // pt_3
 129   // d0_1
 130   // d0_2
 131   // d0_3
 132   // zPrimVert
 133
 134   Bool_t bGenValues = kFALSE;
 135
 136   Int_t pdgCand = -1;
 137   if(fDecay==kDplustoKpipi){
 138     pdgCand=411;
 139   }else if(fDecay==kDstoKKpi){
 140     pdgCand=431;
 141   }else if(fDecay==kLctopKpi){
 142     AliError("LambdaC not yet implemented");
 143     return bGenValues;
 144   }else{
 145     AliError("WRONG DECAY SETTING");
 146     return bGenValues;
 147   }
 148
 149   Double_t vertD[3] = {0,0,0};   // D origin
 150   fmcPartCandidate->XvYvZv(vertD);  // cm
 151
 152   Int_t nprongs = 3;
 153   Int_t daughter[3];
 154   Short_t charge = fmcPartCandidate->Charge();
 155
 156   // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
 157   // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
 158   Int_t index=0;
 159   Int_t nDauLS=0;
 160   Int_t nDauOS=0;
 161
 162
 163   Int_t nDau=fmcPartCandidate->GetNDaughters();
 164   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 165   if(nDau==3){
 166     for(Int_t iDau=0; iDau<3; iDau++){
 167       Int_t ind = labelFirstDau+iDau;
 168       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 169       Short_t signDau=part->Charge();
 170       if(signDau==charge){
 171         nDauLS++;
 172         daughter[index] = ind;
 173         index=2;
 174       }else{
 175         daughter[1] = ind;
 176         nDauOS++;
 177       }
 178     }
 179   }else if(nDau==2){
 180     for(Int_t iDau=0; iDau<2; iDau++){
 181       Int_t ind = labelFirstDau+iDau;
 182       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 183       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 184       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 185         Short_t signDau=part->Charge();
 186         if(signDau==charge){
 187           nDauLS++;
 188           daughter[index] = ind;
 189           index=2;
 190         }else{
 191           daughter[1] = ind;
 192           nDauOS++;
 193         }
 194       }else{
 195         Int_t nDauRes=part->GetNDaughters();
 196         if(nDauRes!=2){
 197           AliError("Wrong resonant decay");
 198           return bGenValues;
 199         }
 200         Int_t labelFirstDauRes = part->GetDaughter(0);
 201         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 202           Int_t indDR = labelFirstDauRes+iDauRes;
 203           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 204           Short_t signDau=partDR->Charge();
 205           if(signDau==charge){
 206             nDauLS++;
 207             daughter[index] = ind;
 208             index=2;
 209           }else{
 210             daughter[1] = ind;
 211             nDauOS++;
 212           }
 213         }
 214       }
 215     }
 216   }else{
 217     AliError(Form("Wrong number of daughters %d",nDau));
 218     return bGenValues;
 219   }
 220
 221   if(nDauLS!=2 || nDauOS!=1){
 222     AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
 223     return bGenValues;
 224   }
 225   if(daughter[0]>daughter[2]){
 226     Int_t tmp=daughter[0];
 227     daughter[0]=daughter[2];
 228     daughter[2]=tmp;
 229   }
 230
 231   // getting the momentum from the daughters and decay vertex
 232   Double_t px[3],py[3],pz[3],pt[3];
 233   Double_t vertDec[3] = {0,0,0};   // decay vertex
 234   for(Int_t iDau=0; iDau<3; iDau++){
 235     AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
 236     px[iDau]=part->Px();
 237     py[iDau]=part->Py();
 238     pz[iDau]=part->Pz();
 239     pt[iDau]=part->Pt();
 240     if(iDau==0) part->XvYvZv(vertDec);
 241   }
 242
 243   Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
 244
 245   AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
 246   Double_t cT = decay->Ct(pdgCand);
 247
 248   vectorMC[0] = fmcPartCandidate->Pt();
 249   vectorMC[1] = fmcPartCandidate->Y() ;
 250   vectorMC[2] = fmcPartCandidate->Phi();
 251   vectorMC[3] = cT*1.E4 ;  // in micron
 252   vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
 253   vectorMC[5] = pt[0];
 254   vectorMC[6] = pt[1];
 255   vectorMC[7] = pt[2];
 256   vectorMC[8] = 0.;   // imppar0, dummy value, meaningless in MC
 257   vectorMC[9] = 0.;   // imppar1, dummy value, meaningless in MC, in micron
 258   vectorMC[10] = 0.;   // imppar2, dummy value, meaningless in MC, in micron
 259   vectorMC[11] = fzMCVertex;    // z of reconstructed of primary vertex
 260
 261   bGenValues = kTRUE;
 262   return bGenValues;
 263 }
 264
 265
 266 //____________________________________________
 267 Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 268 {
 269   // Fill vector (see above) with reconstructed quantities
 270   Bool_t bFillRecoValues=kFALSE;
 271
 272   Int_t pdgCand = -1;
 273   if(fDecay==kDplustoKpipi){
 274     pdgCand=411;
 275   }else if(fDecay==kDstoKKpi){
 276     pdgCand=431;
 277   }else if(fDecay==kLctopKpi){
 278     AliError("LambdaC not yet implemented");
 279     return bFillRecoValues;
 280   }else{
 281     AliError("WRONG DECAY SETTING");
 282     return bFillRecoValues;
 283   }
 284
 285   AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
 286   Short_t charge=decay3->Charge();
 287   Double_t rapidity=decay3->Y(pdgCand);
 288   Double_t cT=decay3->Ct(pdgCand);
 289   Double_t pt = decay3->Pt();
 290   Double_t cosPointingAngle = decay3->CosPointingAngle();
 291   Double_t phi = decay3->Phi();
 292
 293   Int_t daughtSorted[3];
 294   Int_t tmpIndex=0;
 295   Int_t nDauLS=0;
 296   Int_t nDauOS=0;
 297   for(Int_t iDau=0; iDau<3; iDau++){
 298     AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
 299     Int_t label = trk->GetLabel();
 300     Short_t chargedau=trk->Charge();
 301     if(chargedau==charge){
 302       daughtSorted[tmpIndex]=label;
 303       tmpIndex=2;
 304       nDauLS++;
 305     }else{
 306       daughtSorted[1]=label;
 307       nDauOS++;
 308     }
 309   }
 310
 311   if(nDauLS!=2 || nDauOS!=1){
 312     AliError("Wrong decay channel: number of OS and LS tracks not OK");
 313     return bFillRecoValues;
 314   }
 315
 316   if(daughtSorted[0]>daughtSorted[2]){
 317     Int_t tmp=daughtSorted[0];
 318     daughtSorted[0]=daughtSorted[2];
 319     daughtSorted[2]=tmp;
 320   }
 321
 322
 323   vectorReco[0] = pt;
 324   vectorReco[1] = rapidity;
 325   vectorReco[2] = phi;
 326   vectorReco[3] = cT*1.E4;  // in micron
 327   vectorReco[4] = cosPointingAngle;  // in micron
 328   vectorReco[5] = decay3->PtProng(daughtSorted[0]);
 329   vectorReco[6] = decay3->PtProng(daughtSorted[1]);
 330   vectorReco[7] = decay3->PtProng(daughtSorted[2]);
 331   vectorReco[8] = decay3->Getd0Prong(daughtSorted[0]);
 332   vectorReco[9] = decay3->Getd0Prong(daughtSorted[1]);
 333   vectorReco[10] = decay3->Getd0Prong(daughtSorted[2]);
 334   vectorReco[11] = fzPrimVertex;    // z of reconstructed of primary vertex
 335
 336   bFillRecoValues = kTRUE;
 337   return bFillRecoValues;
 338 }
 339
 340
 341 //_____________________________________________________________
 342 Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
 343 {
 344   // Check the pdg codes of the daughters
 345   Bool_t checkCD = kFALSE;
 346
 347   Int_t pdgCand = -1;
 348   Int_t pdgDaughter[3]={-1,-1,-1};
 349   if(fDecay==kDplustoKpipi){
 350     pdgCand=411;
 351     pdgDaughter[0]=321;
 352     pdgDaughter[1]=211;
 353     pdgDaughter[2]=211;
 354   }else if(fDecay==kDstoKKpi){
 355     pdgCand=431;
 356     pdgDaughter[0]=321;
 357     pdgDaughter[1]=321;
 358     pdgDaughter[2]=211;
 359   }else if(fDecay==kLctopKpi){
 360     AliError("LambdaC not yet implemented");
 361     return checkCD;
 362   }else{
 363     AliError("WRONG DECAY SETTING");
 364     return checkCD;
 365   }
 366
 367
 368   Int_t daughter[3];
 369
 370   Int_t nDau=fmcPartCandidate->GetNDaughters();
 371   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 372   if(nDau==3){
 373     for(Int_t iDau=0; iDau<3; iDau++){
 374       Int_t ind = labelFirstDau+iDau;
 375       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 376       daughter[iDau]=TMath::Abs(part->GetPdgCode());
 377     }
 378   }else if(nDau==2){
 379     Int_t nDauFound=0;
 380     for(Int_t iDau=0; iDau<2; iDau++){
 381       Int_t ind = labelFirstDau+iDau;
 382       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 383       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 384       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 385         if(nDauFound>=3) return checkCD;
 386         daughter[nDauFound]=pdgCode;
 387         nDauFound++;
 388       }else{
 389         Int_t nDauRes=part->GetNDaughters();
 390         if(nDauRes!=2) return checkCD;
 391         Int_t labelFirstDauRes = part->GetDaughter(0);
 392         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 393           Int_t indDR = labelFirstDauRes+iDauRes;
 394           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 395           Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 396           if(nDauFound>=3) return checkCD;
 397           daughter[nDauFound]=pdgCodeDR;
 398           nDauFound++;
 399         }
 400       }
 401     }
 402   }else{
 403     return checkCD;
 404   }
 405   for(Int_t iDau1=0; iDau1<3; iDau1++){
 406     for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
 407       if(daughter[iDau1]<daughter[iDau2]){
 408         Int_t tmp=daughter[iDau1];
 409         daughter[iDau1]=daughter[iDau2];
 410         daughter[iDau2]=tmp;
 411       }
 412     }
 413   }
 414   for(Int_t iDau=0; iDau<3; iDau++){
 415     if(daughter[iDau]!=pdgDaughter[iDau]){
 416       AliDebug(2, "Wrong decay channel from MC, skipping!!");
 417       return checkCD;
 418     }
 419   }
 420
 421   checkCD = kTRUE;
 422   return checkCD;
 423
 424 }
 425