PWGHF/vertexingHF/AliRDHFCutsD0toKpipipi.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, 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 for cuts on AOD reconstructed D0->Kpipipi
  21 //
  22 // Author: r.romita@gsi.de, andrea.dainese@pd.infn.it,
  23 //         fabio.colamaria@ba.infn.it
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <TDatabasePDG.h>
  27 #include <Riostream.h>
  28
  29 #include "AliRDHFCutsD0toKpipipi.h"
  30 #include "AliAODRecoDecayHF4Prong.h"
  31 #include "AliAODTrack.h"
  32 #include "AliESDtrack.h"
  33 #include "AliAODPidHF.h"
  34
  35 using std::cout;
  36 using std::endl;
  37
  38 ClassImp(AliRDHFCutsD0toKpipipi)
  39
  40 //--------------------------------------------------------------------------
  41 AliRDHFCutsD0toKpipipi::AliRDHFCutsD0toKpipipi(const char* name) :
  42 AliRDHFCuts(name)
  43 {
  44   //
  45   // Default Constructor
  46   //
  47   Int_t nvars=9;
  48   SetNVars(nvars);
  49   TString varNames[9]={"inv. mass [GeV]",
  50                        "dca [cm]",
  51                        "Dist 2-trk Vtx to PrimVtx [cm]",
  52                        "Dist 3-trk Vtx to PrimVtx [cm]",
  53                        "Dist 4-trk Vtx to PrimVtx [cm]",
  54                        "cosThetaPoint",
  55                        "pt [GeV/c]",
  56                        "rho mass [GeV]",
  57                        "PID cut"};
  58   Bool_t isUpperCut[9]={kTRUE,
  59                         kTRUE,
  60                         kFALSE,
  61                         kFALSE,
  62                         kFALSE,
  63                         kFALSE,
  64                         kFALSE,
  65                         kTRUE,
  66                         kFALSE};
  67   SetVarNames(nvars,varNames,isUpperCut);
  68   Bool_t forOpt[9]={kFALSE,
  69                     kTRUE,
  70                     kTRUE,
  71                     kTRUE,
  72                     kTRUE,
  73                     kTRUE,
  74                     kFALSE,
  75                     kFALSE,
  76                     kFALSE};
  77   SetVarsForOpt(5,forOpt);
  78   Float_t limits[2]={0,999999999.};
  79   SetPtBins(2,limits);
  80 }
  81 //--------------------------------------------------------------------------
  82 AliRDHFCutsD0toKpipipi::AliRDHFCutsD0toKpipipi(const AliRDHFCutsD0toKpipipi &source) :
  83   AliRDHFCuts(source)
  84 {
  85   //
  86   // Copy constructor
  87   //
  88
  89 }
  90 //--------------------------------------------------------------------------
  91 AliRDHFCutsD0toKpipipi &AliRDHFCutsD0toKpipipi::operator=(const AliRDHFCutsD0toKpipipi &source)
  92 {
  93   //
  94   // assignment operator
  95   //
  96   if(&source == this) return *this;
  97
  98   AliRDHFCuts::operator=(source);
  99
 100   return *this;
 101 }
 102
 103
 104 //---------------------------------------------------------------------------
 105 void AliRDHFCutsD0toKpipipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent* aod) {
 106   //
 107   // Fills in vars the values of the variables
 108   //
 109
 110   if(nvars!=fnVarsForOpt) {
 111     printf("AliRDHFCutsD0toKpipipi::GetCutsVarsForOpt: wrong number of variables\n");
 112     return;
 113   }
 114
 115   AliAODRecoDecayHF4Prong *dd = (AliAODRecoDecayHF4Prong*)d;
 116
 117   //recalculate vertex w/o daughters
 118   Bool_t cleanvtx=kFALSE;
 119   AliAODVertex *origownvtx=0x0;
 120   if(fRemoveDaughtersFromPrimary) {
 121     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
 122     cleanvtx=kTRUE;
 123     if(!RecalcOwnPrimaryVtx(dd,aod)) {
 124       CleanOwnPrimaryVtx(dd,aod,origownvtx);
 125       cleanvtx=kFALSE;
 126     }
 127   }
 128
 129   Int_t iter=-1;
 130
 131   if(fVarsForOpt[0]) {
 132     iter++;
 133     Double_t mD0[2],mD0bar[2];
 134     if(TMath::Abs(pdgdaughters[1])==321 || TMath::Abs(pdgdaughters[3])==321) {
 135       dd->InvMassD0(mD0);
 136       if(TMath::Abs(pdgdaughters[1])==321) {
 137        vars[iter]=mD0[0];
 138       }else{
 139        vars[iter]=mD0[1];
 140       }
 141     } else {
 142       dd->InvMassD0bar(mD0bar);
 143       if(TMath::Abs(pdgdaughters[0])==321) {
 144        vars[iter]=mD0bar[0];
 145       }else{
 146        vars[iter]=mD0bar[1];
 147       }
 148    }
 149   }
 150
 151   if(fVarsForOpt[1]){
 152     iter++;
 153     vars[iter]=dd->GetDCA();
 154   }
 155
 156   if(fVarsForOpt[2]){
 157     iter++;
 158     vars[iter]=dd->GetDist12toPrim();
 159   }
 160   if(fVarsForOpt[3]){
 161     iter++;
 162     vars[iter]=dd->GetDist3toPrim();
 163   }
 164   if(fVarsForOpt[4]){
 165     iter++;
 166     vars[iter]=dd->GetDist4toPrim();
 167   }
 168   if(fVarsForOpt[5]){
 169     iter++;
 170     vars[iter]=dd->CosPointingAngle();
 171   }
 172   if(fVarsForOpt[6]){
 173     iter++;
 174     vars[iter]=dd->Pt();
 175   }
 176   if(fVarsForOpt[7]){
 177     iter++;
 178     vars[iter]=999999999.;
 179     printf("ERROR: optmization for rho mass cut not implemented\n");
 180   }
 181   if(fVarsForOpt[8]){
 182     iter++;
 183     vars[iter]=999999999.;
 184     printf("ERROR: optmization for PID cut not implemented\n");
 185   }
 186
 187     if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
 188   return;
 189 }
 190 //---------------------------------------------------------------------------
 191 Int_t AliRDHFCutsD0toKpipipi::IsSelected(TObject* obj,Int_t selectionLevel) {
 192   //
 193   // Apply selection
 194   //
 195
 196   if(!fCutsRD){
 197     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 198     return 0;
 199   }
 200   //PrintAll();
 201   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;
 202
 203   if(!d){
 204     cout<<"AliAODRecoDecayHF4Prong null"<<endl;
 205     return 0;
 206   }
 207
 208   Double_t ptD=d->Pt();
 209   if(ptD<fMinPtCand) return 0;
 210   if(ptD>fMaxPtCand) return 0;
 211
 212   if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
 213
 214   // selection on daughter tracks
 215   if(selectionLevel==AliRDHFCuts::kAll ||
 216      selectionLevel==AliRDHFCuts::kTracks) {
 217     if(!AreDaughtersSelected(d)) return 0;
 218   }
 219
 220
 221   Int_t returnvalue=1;
 222
 223   // selection on candidate
 224   if(selectionLevel==AliRDHFCuts::kAll ||
 225      selectionLevel==AliRDHFCuts::kCandidate) {
 226
 227     Int_t ptbin=PtBin(d->Pt());
 228
 229     Int_t okD0=1,okD0bar=1;
 230     Double_t mD0[2],mD0bar[2];
 231     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 232
 233     d->InvMassD0(mD0);
 234     if(TMath::Abs(mD0[0]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)] &&
 235        TMath::Abs(mD0[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
 236     d->InvMassD0bar(mD0bar);
 237     if(TMath::Abs(mD0bar[0]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)] &&
 238        TMath::Abs(mD0bar[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
 239     if(!okD0 && !okD0bar) return 0;
 240
 241     if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]
 242        || d->GetDCA(3) > fCutsRD[GetGlobalIndex(1,ptbin)]
 243        || d->GetDCA(2) > fCutsRD[GetGlobalIndex(1,ptbin)]
 244        || d->GetDCA(5) > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
 245     if(d->GetDist12toPrim() < fCutsRD[GetGlobalIndex(2,ptbin)]) return 0;
 246     if(d->GetDist3toPrim() < fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;
 247     if(d->GetDist4toPrim() < fCutsRD[GetGlobalIndex(4,ptbin)]) return 0;
 248     if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;
 249     if(d->Pt() < fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;
 250     if(!d->CutRhoMass(mD0,mD0bar,fCutsRD[GetGlobalIndex(0,ptbin)],fCutsRD[GetGlobalIndex(7,ptbin)])) return 0;
 251
 252     if (okD0) returnvalue=1; //cuts passed as D0
 253     if (okD0bar) returnvalue=2; //cuts passed as D0bar
 254     if (okD0 && okD0bar) returnvalue=3; //cuts passed as D0 and D0bar
 255   }
 256
 257   // selection on PID (from AliAODPidHF)
 258   if(selectionLevel==AliRDHFCuts::kAll ||
 259      selectionLevel==AliRDHFCuts::kPID) {
 260
 261     Int_t selD01 = D01Selected(d,AliRDHFCuts::kCandidate);
 262     Int_t selD02 = D02Selected(d,AliRDHFCuts::kCandidate);
 263     Int_t selD0bar1 = D0bar1Selected(d,AliRDHFCuts::kCandidate);
 264     Int_t selD0bar2 = D0bar2Selected(d,AliRDHFCuts::kCandidate);
 265
 266     Int_t d01PID = 0, d02PID = 0, d0bar1PID = 0, d0bar2PID = 0;
 267     returnvalue = IsSelectedFromPID(d, &d01PID, &d02PID, &d0bar1PID, &d0bar2PID);  //This returnvalue is dummy! Now it's modified as it must be!
 268
 269 returnvalue = 0;
 270
 271     if((selD01 == 1 && d01PID == 1)||(selD02 == 1 && d02PID == 1)||(selD0bar1 == 1 && d0bar1PID == 1)||(selD0bar2 == 1 && d0bar2PID == 1)) returnvalue = 1;
 272   }
 273
 274   return returnvalue;
 275 }
 276
 277 //---------------------------------------------------------------------------
 278 Int_t AliRDHFCutsD0toKpipipi::IsSelectedFromPID(AliAODRecoDecayHF4Prong *d, Int_t *hyp1, Int_t *hyp2, Int_t *hyp3, Int_t *hyp4) {
 279   //
 280   // Apply selection (using AliAODPidHF methods)
 281   // Mass hypothesis true if each particle is at least compatible with specie of hypothesis
 282   //
 283
 284   Int_t output=0;
 285
 286   Int_t matchK[4], matchPi[4];
 287   Double_t ptlimit[2] = {0.6,0.8};
 288   AliAODTrack* trk[4];
 289   trk[0] = (AliAODTrack*)d->GetDaughter(0);
 290   trk[1] = (AliAODTrack*)d->GetDaughter(1);
 291   trk[2] = (AliAODTrack*)d->GetDaughter(2);
 292   trk[3] = (AliAODTrack*)d->GetDaughter(3);
 293 //  if(!trk[0] || !trk[1] || !trk[2] || !trk[3]) {          //REMOVED (needs also AliAODEvent to be passed, here and in IsSelected method)
 294 //    trk[0]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 295 //    trk[1]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
 296 //    trk[2]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(2)]);
 297 //    trk[3]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(3)]);}
 298
 299   AliAODPidHF* pidObj = new AliAODPidHF();
 300   pidObj->SetAsym(kTRUE);
 301   pidObj->SetPLimit(ptlimit);
 302   pidObj->SetSigma(0,2.);  //TPC sigma, in three pT ranges
 303   pidObj->SetSigma(1,1.);
 304   pidObj->SetSigma(2,0.);
 305   pidObj->SetSigma(3,2.);  //TOF sigma, whole pT range
 306   pidObj->SetTPC(kTRUE);
 307   pidObj->SetTOF(kTRUE);
 308   pidObj->SetMatch(1);
 309   pidObj->SetCompat(kTRUE);
 310   for(Int_t ii=0; ii<4; ii++) {
 311     pidObj->SetSigma(0,2.);
 312     matchK[ii] = pidObj->MatchTPCTOF(trk[ii],3); //arguments: track, mode, particle#, compatibility allowed
 313     pidObj->SetSigma(0,2.);
 314     matchPi[ii] = pidObj->MatchTPCTOF(trk[ii],2);
 315   }
 316
 317   //Rho invariant mass under various hypotheses (to be matched with PID infos in order to accet the candidate)
 318   Int_t d01rho03 = 0, d01rho23 = 0, d02rho01 = 0, d02rho12 = 0, d0bar1rho12 = 0, d0bar1rho23 = 0, d0bar2rho01 = 0, d0bar2rho03 = 0;
 319   if(TMath::Abs(0.775 - d->InvMassRho(0,3))<0.1)  {d01rho03 = 1; d0bar2rho03 = 1;}
 320   if(TMath::Abs(0.775 - d->InvMassRho(2,3))<0.1)  {d01rho23 = 1; d0bar1rho23 = 1;}
 321   if(TMath::Abs(0.775 - d->InvMassRho(0,1))<0.1)  {d02rho01 = 1; d0bar2rho01 = 1;}
 322   if(TMath::Abs(0.775 - d->InvMassRho(1,2))<0.1)  {d02rho12 = 1; d0bar1rho12 = 1;}
 323   Int_t d01rho = 0, d02rho = 0, d0bar1rho = 0, d0bar2rho = 0;
 324   if(d01rho03==1||d01rho23==1) d01rho = 1;
 325   if(d02rho01==1||d02rho12==1) d02rho = 1;
 326   if(d0bar1rho12==1||d0bar1rho23==1) d0bar1rho = 1;
 327   if(d0bar2rho01==1||d0bar2rho03==1) d0bar2rho = 1;
 328
 329   //This way because there could be multiple hypotheses accepted
 330   if(d01rho==1 && (matchK[1]>=0 && matchPi[0]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp1 = 1; output = 1;} //d01 hyp
 331   if(d02rho==1 && (matchK[3]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0)) {*hyp2 = 1; output = 1;} //d02 hyp
 332   if(d0bar1rho==1 && (matchK[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp3 = 1; output = 1;} //d0bar1 hyp
 333   if(d0bar2rho==1 && (matchK[2]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[3]>=0)) {*hyp4 = 1; output = 1;} //d0bar2 hyp
 334
 335   return output;
 336 }
 337 //---------------------------------------------------------------------------
 338 Int_t AliRDHFCutsD0toKpipipi::D01Selected(TObject* obj,Int_t selectionLevel) {
 339   //
 340   // Apply selection
 341   //
 342
 343   if(!fCutsRD){
 344     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 345     return 0;
 346   }
 347   //PrintAll();
 348   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;
 349
 350   if(!d){
 351     cout<<"AliAODRecoDecayHF4Prong null"<<endl;
 352     return 0;
 353   }
 354
 355   Int_t returnvalue=0;
 356
 357   // selection on candidate
 358   if(selectionLevel==AliRDHFCuts::kAll ||
 359      selectionLevel==AliRDHFCuts::kCandidate) {
 360
 361     Int_t ptbin=PtBin(d->Pt());
 362
 363     Double_t mD0[2];
 364     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 365
 366     d->InvMassD0(mD0);
 367     if(TMath::Abs(mD0[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
 368   }
 369
 370   return returnvalue;
 371 }
 372 //---------------------------------------------------------------------------
 373 Int_t AliRDHFCutsD0toKpipipi::D02Selected(TObject* obj,Int_t selectionLevel) {
 374   //
 375   // Apply selection
 376   //
 377
 378   if(!fCutsRD){
 379     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 380     return 0;
 381   }
 382   //PrintAll();
 383   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;
 384
 385   if(!d){
 386     cout<<"AliAODRecoDecayHF4Prong null"<<endl;
 387     return 0;
 388   }
 389
 390   Int_t returnvalue=0;
 391
 392   // selection on candidate
 393   if(selectionLevel==AliRDHFCuts::kAll ||
 394      selectionLevel==AliRDHFCuts::kCandidate) {
 395
 396     Int_t ptbin=PtBin(d->Pt());
 397
 398     Double_t mD0[2];
 399     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 400
 401     d->InvMassD0(mD0);
 402     if(TMath::Abs(mD0[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
 403   }
 404
 405   return returnvalue;
 406 }//---------------------------------------------------------------------------
 407 Int_t AliRDHFCutsD0toKpipipi::D0bar1Selected(TObject* obj,Int_t selectionLevel) {
 408   //
 409   // Apply selection
 410   //
 411
 412   if(!fCutsRD){
 413     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 414     return 0;
 415   }
 416   //PrintAll();
 417   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;
 418
 419   if(!d){
 420     cout<<"AliAODRecoDecayHF4Prong null"<<endl;
 421     return 0;
 422   }
 423
 424   Int_t returnvalue=0;
 425
 426   // selection on candidate
 427   if(selectionLevel==AliRDHFCuts::kAll ||
 428      selectionLevel==AliRDHFCuts::kCandidate) {
 429
 430     Int_t ptbin=PtBin(d->Pt());
 431
 432     Double_t mD0bar[2];
 433     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 434
 435     d->InvMassD0bar(mD0bar);
 436     if(TMath::Abs(mD0bar[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
 437   }
 438
 439   return returnvalue;
 440 }
 441 //---------------------------------------------------------------------------
 442 Int_t AliRDHFCutsD0toKpipipi::D0bar2Selected(TObject* obj,Int_t selectionLevel) {
 443   //
 444   // Apply selection
 445   //
 446
 447   if(!fCutsRD){
 448     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 449     return 0;
 450   }
 451   //PrintAll();
 452   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;
 453
 454   if(!d){
 455     cout<<"AliAODRecoDecayHF4Prong null"<<endl;
 456     return 0;
 457   }
 458
 459   Int_t returnvalue=0;
 460
 461   // selection on candidate
 462   if(selectionLevel==AliRDHFCuts::kAll ||
 463      selectionLevel==AliRDHFCuts::kCandidate) {
 464
 465     Int_t ptbin=PtBin(d->Pt());
 466
 467     Double_t mD0bar[2];
 468     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 469
 470     d->InvMassD0bar(mD0bar);
 471     if(TMath::Abs(mD0bar[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
 472   }
 473
 474   return returnvalue;
 475 }
 476 //---------------------------------------------------------------------------
 477 Bool_t AliRDHFCutsD0toKpipipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 478 {
 479   //
 480   // Checking if D0 is in fiducial acceptance region
 481   //
 482
 483   if(fMaxRapidityCand>-998.){
 484     if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
 485     else return kTRUE;
 486   }
 487
 488   if(pt > 5.) {
 489     // applying cut for pt > 5 GeV
 490     AliDebug(4,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
 491     if (TMath::Abs(y) > 0.8){
 492       return kFALSE;
 493     }
 494   } else {
 495     // appliying smooth cut for pt < 5 GeV
 496     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 497     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 498     AliDebug(4,Form("pt of D0 = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 499     if (y < minFiducialY || y > maxFiducialY){
 500       return kFALSE;
 501     }
 502   }
 503
 504   return kTRUE;
 505 }
 506