PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx

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