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