HLT/trigger/AliD0toKpi.cxx

   1 #include <iostream.h>
   2 #include <TH1.h>
   3 #include <TH2.h>
   4 #include <TCanvas.h>
   5 #include <TPaveLabel.h>
   6 #include "AliD0toKpi.h"
   7
   8 ClassImp(AliD0toKpi)
   9   /*******************************************************************
  10    *                                                                 *
  11    *  Reconstructed D^0 -> K^- pi^+   candidate class                *
  12    *                                                                 *
  13    *  origin: A. Dainese    andrea.dainese@pd.infn.it                *
  14    *******************************************************************/
  15
  16 //________________________________________________________
  17 AliD0toKpi::AliD0toKpi() {
  18
  19   // Default constructor
  20   fSignal = kFALSE;
  21
  22   fEvent = 0;
  23
  24   fV1x = 0.;
  25   fV1y = 0.;
  26   fV1z = 0.;
  27   fV2x = 0.;
  28   fV2y = 0.;
  29   fV2z = 0.;
  30   fDCA = 0.;
  31
  32   fPx[0] = 0.;
  33   fPy[0] = 0.;
  34   fPz[0] = 0.;
  35   fPx[1] = 0.;
  36   fPy[1] = 0.;
  37   fPz[1] = 0.;
  38
  39   fd0[0] = 0.;
  40   fd0[1] = 0.;
  41
  42   fPdg[0] = 0;
  43   fPdg[1] = 0;
  44   fMum[0] =0;
  45   fMum[1] =0;
  46
  47
  48   fWgtAD0 = fWgtAD0bar = fWgtBD0 = fWgtBD0bar = fWgtCD0 = fWgtCD0bar = 0;
  49
  50 }
  51
  52 //________________________________________________________
  53 AliD0toKpi::AliD0toKpi(Bool_t sgn,Int_t ev,
  54                        Double_t v1[3],Double_t v2[3],
  55                        Double_t dca,
  56                        Double_t mom[6],Double_t d0[2],
  57                        Int_t pdg[2],Int_t mum[2]) {
  58
  59   // Standard constructor
  60   fSignal = sgn;
  61
  62   fEvent = ev;
  63
  64   fV1x = v1[0];
  65   fV1y = v1[1];
  66   fV1z = v1[2];
  67   fV2x = v2[0];
  68   fV2y = v2[1];
  69   fV2z = v2[2];
  70   fDCA = dca;
  71
  72   fPx[0] = mom[0];
  73   fPy[0] = mom[1];
  74   fPz[0] = mom[2];
  75   fPx[1] = mom[3];
  76   fPy[1] = mom[4];
  77   fPz[1] = mom[5];
  78
  79   fd0[0] =  d0[0];
  80   fd0[1] = d0[1];
  81
  82   fPdg[0] = pdg[0];
  83   fPdg[1] = pdg[1];
  84   fMum[0] = mum[0];
  85   fMum[1] = mum[1];
  86
  87
  88   fWgtAD0 = fWgtAD0bar = fWgtBD0 = fWgtBD0bar = fWgtCD0 = fWgtCD0bar = 0;
  89 }
  90
  91 //_____________________________________________________________________
  92 Double_t AliD0toKpi::EtaChild(Int_t child) const {
  93
  94   Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
  95   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
  96   return eta;
  97 }
  98 //_____________________________________________________________________
  99 Double_t AliD0toKpi::Eta() const {
 100
 101   Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
 102   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
 103   return eta;
 104 }
 105
 106 //_____________________________________________________________________
 107 Double_t AliD0toKpi::CPta() const {
 108
 109   TVector3 mom(Px(),Py(),Pz());
 110   TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
 111
 112   Double_t pta = mom.Angle(flight);
 113
 114   return TMath::Cos(pta);
 115 }
 116 //_____________________________________________________________________
 117 Double_t AliD0toKpi::CPtaXY() const {
 118
 119   TVector3 momXY(Px(),Py(),0.);
 120   TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
 121
 122   Double_t ptaXY = momXY.Angle(flightXY);
 123
 124   return TMath::Cos(ptaXY);
 125 }
 126 //_____________________________________________________________________
 127 Double_t AliD0toKpi::ChildrenRelAngle() const {
 128
 129   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
 130   TVector3 mom1(fPx[1],fPy[1],fPz[1]);
 131
 132   Double_t angle = mom0.Angle(mom1);
 133
 134   return angle;
 135 }
 136 //_____________________________________________________________________
 137 void AliD0toKpi::DrawPIDinTOF() const {
 138
 139   TH2F* framePi = new TH2F("framePi","Tag probabilities for PIONS",2,0,2.5,2,0,1);
 140   framePi->SetXTitle("p [GeV/c]");
 141   framePi->SetStats(0);
 142   TH2F* frameK = new TH2F("frameK","Tag probabilities for KAONS",2,0,2.5,2,0,1);
 143   frameK->SetXTitle("p [GeV/c]");
 144   frameK->SetStats(0);
 145   TH2F* frameP = new TH2F("frameP","Tag probabilities for PROTONS",2,0,4.5,2,0,1);
 146   frameP->SetXTitle("p [GeV/c]");
 147   frameP->SetStats(0);
 148
 149   TH1F* hPiPi = new TH1F("hPiPi","Tag probabilities for PIONS",kPiBins,0,2.5);
 150   TH1F* hPiNid = new TH1F("hPiNid","Tag probabilities for PIONS",kPiBins,0,2.5);
 151
 152   TH1F* hKK = new TH1F("hKK","Tag probabilities for KAONS",kKBins,0,2.5);
 153   TH1F* hKNid = new TH1F("hKNid","Tag probabilities for KAONS",kKBins,0,2.5);
 154   TH1F* hKPi = new TH1F("hKPi","Tag probabilities for KAONS",kKBins,0,2.5);
 155
 156   TH1F* hPP = new TH1F("hPP","Tag probabilities for PROTONS",kPBins,0,4.5);
 157   TH1F* hPNid = new TH1F("hPNid","Tag probabilities for PROTONS",kPBins,0,4.5);
 158   TH1F* hPPi = new TH1F("hPPi","Tag probabilities for PROTONS",kPBins,0,4.5);
 159
 160   for(Int_t i=1; i<=kPiBins; i++) {
 161     hPiPi->SetBinContent(i,kPiTagPi[i-1]);
 162     hPiNid->SetBinContent(i,kPiTagPi[i-1]+kPiTagNid[i-1]);
 163
 164     hKK->SetBinContent(i,kKTagK[i-1]);
 165     hKPi->SetBinContent(i,kKTagK[i-1]+kKTagPi[i-1]);
 166     hKNid->SetBinContent(i,kKTagK[i-1]+kKTagPi[i-1]+kKTagNid[i-1]);
 167   }
 168   for(Int_t i=1; i<=kPBins; i++) {
 169     hPP->SetBinContent(i,kPTagP[i-1]);
 170     hPPi->SetBinContent(i,kPTagP[i-1]+kPTagPi[i-1]);
 171     hPNid->SetBinContent(i,kPTagP[i-1]+kPTagPi[i-1]+kPTagNid[i-1]);
 172   }
 173
 174   TCanvas* c = new TCanvas("c","PID",0,0,1000,400);
 175   c->Divide(3,1);
 176   c->cd(1);
 177   framePi->Draw();
 178   hPiNid->SetFillColor(18); hPiNid->Draw("same");
 179   hPiPi->SetFillColor(4); hPiPi->Draw("same");
 180   TPaveLabel* pav1 = new TPaveLabel(1,.2,1.4,.3,"#pi");
 181   pav1->SetBorderSize(0);
 182   pav1->Draw("same");
 183   TPaveLabel* pav2 = new TPaveLabel(1,.8,1.8,.9,"non-id");
 184   pav2->SetBorderSize(0);
 185   pav2->Draw("same");
 186
 187   c->cd(2);
 188   frameK->Draw();
 189   hKNid->SetFillColor(18); hKNid->Draw("same");
 190   hKPi->SetFillColor(4); hKPi->Draw("same");
 191   hKK->SetFillColor(7); hKK->Draw("same");
 192   TPaveLabel* pav3 = new TPaveLabel(1,.2,1.5,.3,"K");
 193   pav3->SetBorderSize(0);
 194   pav3->Draw("same");
 195   TPaveLabel* pav4 = new TPaveLabel(1,.8,1.8,.9,"non-id");
 196   pav4->SetBorderSize(0);
 197   pav4->Draw("same");
 198   TPaveLabel* pav5 = new TPaveLabel(.4,.5,.8,.6,"#pi");
 199   pav5->SetBorderSize(0);
 200   pav5->Draw("same");
 201
 202   c->cd(3);
 203   frameP->Draw();
 204   hPNid->SetFillColor(18); hPNid->Draw("same");
 205   hPPi->SetFillColor(4); hPPi->Draw("same");
 206   hPP->SetFillColor(3); hPP->Draw("same");
 207   TPaveLabel* pav6 = new TPaveLabel(1,.2,1.5,.3,"p");
 208   pav6->SetBorderSize(0);
 209   pav6->Draw("same");
 210   TPaveLabel* pav7 = new TPaveLabel(1,.8,2.6,.9,"non-id");
 211   pav7->SetBorderSize(0);
 212   pav7->Draw("same");
 213   TPaveLabel* pav8 = new TPaveLabel(.2,.5,1,.6,"#pi");
 214   pav8->SetBorderSize(0);
 215   pav8->Draw("same");
 216
 217
 218   return;
 219 }
 220 //___________________________________________________________________________
 221 Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,const Double_t *values) const {
 222
 223   Double_t value=0;
 224   Double_t slope;
 225
 226   if(p<0.5*Bin) {
 227     value = values[0];
 228   } else if(p>=(nBins-0.5)*Bin) {
 229     //slope = (values[nBins-1]-values[nBins-2])/Bin;
 230     slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
 231     value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
 232   } else {
 233     for(Int_t i=0; i<nBins; i++) {
 234       if(p<(i+0.5)*Bin) {
 235         slope = (values[i]-values[i-1])/Bin;
 236         value = values[i-1]+slope*(p-Bin*(i-0.5));
 237         break;
 238       }
 239     }
 240   }
 241
 242   if(value<0) value=0.;
 243   if(value>1) value=1.;
 244   return value;
 245 }
 246 //______________________________________________________________________
 247 void AliD0toKpi::ComputeWgts() {
 248
 249   Double_t probTagPi[2]  = {0.,0.};
 250   Double_t probTagK[2]   = {0.,0.};
 251   Double_t probTagNid[2] = {0.,0.};
 252   Double_t probTagP[2]   = {0.,0.};
 253
 254   // tagging of the positive track
 255   if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13
 256      || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
 257     probTagPi[0]  = LinearInterpolation(PChild(0),kPiBins,kPiBinWidth,kPiTagPi);
 258     probTagNid[0] = 1.-probTagPi[0];
 259     probTagK[0]   = 0.;
 260     probTagP[0]   = 0.;
 261   }
 262   if(TMath::Abs(fPdg[0])==321) { // kaon
 263     probTagK[0]   = LinearInterpolation(PChild(0),kKBins,kKBinWidth,kKTagK);
 264     probTagNid[0] = LinearInterpolation(PChild(0),kKBins,kKBinWidth,kKTagNid);
 265     if((probTagNid[0]+probTagK[0])>1.) probTagNid[0] = 1.-probTagK[0];
 266     probTagPi[0] = 1.-probTagNid[0]-probTagK[0];
 267     probTagP[0] = 0.;
 268   }
 269   if(TMath::Abs(fPdg[0])==2212) { // proton
 270     probTagP[0]  = LinearInterpolation(PChild(0),kPBins,kPBinWidth,kPTagP);
 271     probTagNid[0] = LinearInterpolation(PChild(0),kPBins,kPBinWidth,kPTagNid);
 272     if((probTagNid[0]+probTagP[0])>1.) probTagNid[0] = 1.-probTagP[0];
 273     probTagPi[0] = 1.-probTagNid[0]-probTagP[0];
 274     probTagK[0]   = 0.;
 275   }
 276
 277
 278   // tagging of the negative track
 279   if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13
 280      || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
 281     probTagPi[1]  = LinearInterpolation(PChild(1),kPiBins,kPiBinWidth,kPiTagPi);
 282     probTagNid[1] = 1.-probTagPi[1];
 283     probTagK[1]   = 0.;
 284     probTagP[1]   = 0.;
 285   }
 286   if(TMath::Abs(fPdg[1])==321) { // kaon
 287     probTagK[1]   = LinearInterpolation(PChild(1),kKBins,kKBinWidth,kKTagK);
 288     probTagNid[1] = LinearInterpolation(PChild(1),kKBins,kKBinWidth,kKTagNid);
 289     if((probTagNid[1]+probTagK[1])>1.) probTagNid[1] = 1.-probTagK[1];
 290     probTagPi[1] = 1.-probTagNid[1]-probTagK[1];
 291     probTagP[1] = 0.;
 292   }
 293   if(TMath::Abs(fPdg[1])==2212) { // proton
 294     probTagP[1]  = LinearInterpolation(PChild(1),kPBins,kPBinWidth,kPTagP);
 295     probTagNid[1] = LinearInterpolation(PChild(1),kPBins,kPBinWidth,kPTagNid);
 296     if((probTagNid[1]+probTagP[1])>1.) probTagNid[1] = 1.-probTagP[1];
 297     probTagPi[1] = 1.-probTagNid[1]-probTagP[1];
 298     probTagK[1]   = 0.;
 299   }
 300
 301   // assignement of the weights from PID
 302   fWgtAD0    = probTagK[1]*(probTagPi[0]+probTagNid[0]);
 303   fWgtAD0bar = probTagK[0]*(probTagPi[1]+probTagNid[1]);
 304   fWgtBD0    = probTagPi[0]*probTagNid[1];
 305   fWgtBD0bar = probTagPi[1]*probTagNid[0];
 306   fWgtCD0    = probTagNid[0]*probTagNid[1];
 307   fWgtCD0bar = probTagNid[0]*probTagNid[1];
 308
 309   // pt-dependent weight to reproduce pt distr. given by NLO QCD (MNR)
 310   // [only for candidates with at least one track coming from a D0]
 311   SetPtWgts();
 312
 313
 314
 315   /*
 316   for(Int_t j=0;j<2;j++) cerr<<" PDG = "<<GetPdgChild(j)<<" p = "<<PChild(j)<<" TagPi = "<<ProbTagPi[j]<<" TagK = "<<ProbTagK[j]<<" TagNid = "<<ProbTagNid[j]<<endl;
 317
 318   cerr<<fWgtAD0<<"  "<<fWgtAD0bar<<endl;
 319   cerr<<fWgtBD0<<"  "<<fWgtBD0bar<<endl;
 320   cerr<<fWgtCD0<<"  "<<fWgtCD0bar<<endl;
 321
 322   if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 323   if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 324   if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 325   if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 326   if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 327   if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 328   */
 329
 330   return;
 331 }
 332 //____________________________________________________________________________
 333 void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
 334
 335   fWgtAD0    *= factor;
 336   fWgtAD0bar *= factor;
 337   fWgtBD0    *= factor;
 338   fWgtBD0bar *= factor;
 339   fWgtCD0    *= factor;
 340   fWgtCD0bar *= factor;
 341
 342   return;
 343 }
 344 //____________________________________________________________________________
 345 void AliD0toKpi::SetPtWgts() {
 346
 347   if(TMath::Abs(fMum[0]) != 421 && TMath::Abs(fMum[1]) != 421 &&
 348      TMath::Abs(fMum[0]) != 411 && TMath::Abs(fMum[1]) != 411) return;
 349
 350   Double_t ptWgt = 1.;
 351   ptWgt = 2.05-0.47*Pt()+0.02*Pt()*Pt();
 352   if(Pt() >= 5.) ptWgt = 0.56*TMath::Exp(-0.12*Pt());
 353
 354   fWgtAD0    *= ptWgt;
 355   fWgtAD0bar *= ptWgt;
 356   fWgtBD0    *= ptWgt;
 357   fWgtBD0bar *= ptWgt;
 358   fWgtCD0    *= ptWgt;
 359   fWgtCD0bar *= ptWgt;
 360
 361   return;
 362 }
 363 //____________________________________________________________________________
 364 void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,Option_t *sample) const {
 365
 366   const char *sampleA = strstr(sample,"A");
 367   const char *sampleB = strstr(sample,"B");
 368   const char *sampleC = strstr(sample,"C");
 369
 370   if(sampleA) { WgtD0 = fWgtAD0;  WgtD0bar = fWgtAD0bar; }
 371   if(sampleB) { WgtD0 = fWgtBD0;  WgtD0bar = fWgtBD0bar; }
 372   if(sampleC) { WgtD0 = fWgtCD0;  WgtD0bar = fWgtCD0bar; }
 373
 374   if(fSignal) {
 375     if(fMum[0]==421)  WgtD0bar = 0.;
 376     if(fMum[0]==-421) WgtD0 = 0.;
 377   }
 378
 379   return;
 380 }
 381 //____________________________________________________________________________
 382 void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
 383
 384   Double_t energy[2];
 385
 386   // D0 -> K- Pi+
 387   energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
 388   energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));
 389
 390   mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
 391
 392
 393   // D0bar -> K+ Pi-
 394   energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
 395   energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));
 396
 397   mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
 398
 399   return;
 400
 401 }
 402 //________________________________________________________________
 403 Double_t AliD0toKpi::qT() const {
 404
 405   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
 406   TVector3 momD(Px(),Py(),Pz());
 407
 408   return mom0.Perp(momD);
 409 }
 410 //________________________________________________________________
 411 Double_t AliD0toKpi::qL(Int_t child) const {
 412
 413   Double_t qL;
 414   TVector3 mom(fPx[child],fPy[child],fPz[child]);
 415   TVector3 momD(Px(),Py(),Pz());
 416
 417   qL = mom.Dot(momD)/momD.Mag();
 418
 419   return qL ;
 420 }
 421 //________________________________________________________________
 422 void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
 423
 424   Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
 425
 426   Double_t beta = P()/Energy();
 427   Double_t gamma = Energy()/kMD0;
 428
 429   ctsD0 = (qL(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 430 // if(ctsD0 > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
 431 // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
 432
 433   ctsD0bar = (qL(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 434 // if(ctsD0bar > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
 435 // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";}
 436
 437   return;
 438 }
 439 //____________________________________________________________________
 440 Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) const {
 441 //
 442 // This function compares the D0 with a set of cuts:
 443 //
 444 // cuts[0] = inv. mass half width [GeV]
 445 // cuts[1] = dca [micron]
 446 // cuts[2] = cosThetaStar
 447 // cuts[3] = pTK [GeV/c]
 448 // cuts[4] = pTPi [GeV/c]
 449 // cuts[5] = d0K [micron]   upper limit!
 450 // cuts[6] = d0Pi [micron]  upper limit!
 451 // cuts[7] = d0d0 [micron^2]
 452 // cuts[8] = cosThetaPoint
 453 //
 454 // If the the D0/D0bar doesn't pass the cuts it sets the weights to 0
 455 // If neither D0 nor D0bar pass the cuts return kFALSE
 456 //
 457   Double_t mD0,mD0bar,ctsD0,ctsD0bar;
 458   okD0=1; okD0bar=1;
 459
 460   if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
 461   if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
 462   if(!okD0 && !okD0bar) return kFALSE;
 463
 464   if(TMath::Abs(Getd0Child(1)) > cuts[5] ||
 465      TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
 466   if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
 467      TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
 468   if(!okD0 && !okD0bar) return kFALSE;
 469
 470   if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
 471
 472   InvMass(mD0,mD0bar);
 473   if(TMath::Abs(mD0-kMD0) > cuts[0])    okD0 = 0;
 474   if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
 475   if(!okD0 && !okD0bar) return kFALSE;
 476
 477   CosThetaStar(ctsD0,ctsD0bar);
 478   if(TMath::Abs(ctsD0) > cuts[2])    okD0 = 0;
 479   if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
 480   if(!okD0 && !okD0bar) return kFALSE;
 481
 482   if(d0d0() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
 483
 484   if(CPta() < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
 485
 486   return kTRUE;
 487 }
 488
 489
 490
 491
 492
 493
 494