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