PWGPP/pid/doeffPr.C

   1 TObject* fContPid1;
   2 TObject* fContPid2;
   3 const Int_t nBinPid = 14; // pt,eta, ptPip, ptPin, PPip, PPin, TOF3sigmaPip, TOF3sigmaPin, isPhiTrue, nsigmaPip, nsigmaPin
   4 // 0.985 < mass < 1.045 (60) and 0 < centrality < 100 (10)
   5 Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaPr*/,20/*pt+*/,20/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
   6 Float_t xmin[nBinPid] = {1,-0.8,0.3,0.3,0,0,-0.5,-0.5,-0.5,0,0,-TMath::Pi(),-TMath::Pi(),-TMath::Pi()/2};
   7 Float_t xmax[nBinPid] = {5,0.8,4.3,4.3,1,1,1.5,1.5,1.5,7.5,7.5,TMath::Pi(),TMath::Pi(),TMath::Pi()/2};
   8
   9 TF1 *fsign;
  10 TF1 *fall;
  11 TF1 *fback;
  12
  13 Int_t ifunc=0;
  14
  15 Float_t fitmin = 1.08;
  16 Float_t fitmax = 1.155;
  17
  18 Int_t cmin = 1;
  19 Int_t cmax = 8;
  20
  21 Float_t weightS = -1.;
  22
  23 Int_t rebinsize = 2;
  24
  25 Int_t parplotted = 2;
  26
  27 Bool_t isMC = kFALSE; // don't change this (is set automatically)
  28 Bool_t selectTrue = kTRUE; // put it to true to remove background (only for MC)
  29 Bool_t keepTrue = kFALSE; // put it to false to fit only background (only for MC)
  30
  31 Bool_t kGoodMatch = kFALSE; // to check good matching
  32
  33 Bool_t kSigma2vs3 = kFALSE; // to check good matching
  34
  35 Bool_t require5sigma = kFALSE; // don't touch this flag
  36
  37 Bool_t bayesVsigma = kFALSE; // only to do checks
  38
  39 Bool_t kTOFmatch = kFALSE; // for combined PID requires TOF matching
  40
  41 Bool_t kOverAll = kTRUE;
  42 Bool_t kOverAllTOFmatch = kFALSE;
  43 Bool_t kOverAll2Sigma = kTRUE;
  44
  45 Bool_t kLoaded=kFALSE;
  46 LoadLib(){
  47   weightS = -1.;
  48
  49   if(! kLoaded){
  50     gSystem->Load("libVMC.so");
  51     gSystem->Load("libPhysics.so");
  52     gSystem->Load("libTree.so");
  53     gSystem->Load("libMinuit.so");
  54     gSystem->Load("libSTEERBase.so");
  55     gSystem->Load("libANALYSIS.so");
  56     gSystem->Load("libAOD.so");
  57     gSystem->Load("libESD.so");
  58     gSystem->Load("libANALYSIS.so");
  59     gSystem->Load("libANALYSISalice.so");
  60     gSystem->Load("libCORRFW.so");
  61     gSystem->Load("libNetx.so");
  62     gSystem->Load("libPWGPPpid.so");
  63
  64     TFile *f = new TFile("AnalysisResults.root");
  65     f->ls();
  66     TList *l = (TList *) f->Get("contLambdaBayes1");
  67     l->ls();
  68     fContPid1 = (AliPIDperfContainer *) l->At(0);
  69     fContPid2 = (AliPIDperfContainer *) l->At(1);
  70   }
  71   kLoaded = kTRUE;
  72
  73   // check if MC
  74   Float_t x[] = {xmin[0]+0.001,xmin[1]+0.001,xmin[2]+0.001,xmin[3]+0.001,xmin[4]+0.001,xmin[5]+0.001,xmin[6]+0.001,xmin[7]+0.001,1/*trueMC*/,xmin[9],xmin[10]};
  75   Float_t x2[] = {xmax[0],xmax[1],xmax[2],xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],xmax[8],xmax[9],xmax[10]};
  76
  77   AliPIDperfContainer *tmp = fContPid1;
  78   TH1D *h = tmp->GetQA(0, x, x2)->ProjectionX("checkMC");
  79
  80   if(h->GetEntries()) isMC = kTRUE;
  81   else isMC=kFALSE;
  82
  83   if(!isMC){
  84     selectTrue = kFALSE;
  85     keepTrue = kTRUE;
  86   }
  87   else{
  88     printf("MC truth found!!!!!!\nIt is MC!!!!!!");
  89   }
  90
  91   fsign = new TF1("fsign","[0]*TMath::Voigt(x-[1],[3],[2])",fitmin,fitmax);
  92   fback = new TF1("fback","pol2",fitmin,fitmax);
  93   fall = new TF1("fall","[0]*TMath::Voigt(x-[1],[3],[2]) + pol2(4)",fitmin,fitmax);
  94
  95   fsign->SetLineColor(2);
  96   fback->SetLineColor(4);
  97
  98   if(kSigma2vs3){
  99     kGoodMatch=kFALSE;
 100     kOverAll = 0;
 101   }
 102
 103   if(bayesVsigma){
 104     kOverAll = 0;
 105     kGoodMatch=kFALSE;
 106     kSigma2vs3=kFALSE;
 107     kTOFmatch=kTRUE;
 108     weightS = -0.7;
 109   }
 110 }
 111
 112 doeffPr(Int_t pos=1,Float_t prob=0.1,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
 113   LoadLib();
 114
 115   Int_t nptbin = binPid[2];
 116   Float_t minptbin = xmin[2];
 117   Float_t maxptbin = xmax[2];
 118
 119   if(pos == 0){
 120     nptbin = binPid[3];
 121     minptbin = xmin[3];
 122     maxptbin = xmax[3];
 123   }
 124
 125   if(prob > 0.1999){
 126     kGoodMatch = kFALSE;
 127     kSigma2vs3 = kFALSE;
 128     if(! kOverAll) require5sigma = kTRUE;
 129     if(!isMC) weightS = -0.95;
 130   }
 131
 132   TCanvas *c = new TCanvas();
 133   c->Divide((nptbin+1)/2,2);
 134   TH2F *hh.*hh2;
 135   TH1D *h,*h2;
 136   char name[100];
 137   Float_t b[50][3];
 138
 139   Double_t xx[50],yy[50];
 140   Double_t exx[50],eyy[50];
 141
 142   for(Int_t i=0;i < nptbin;i++){
 143     c->cd(i+1)->SetLogy();
 144     Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
 145     Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
 146
 147     xx[i] = (ptmin+ptmax)/2;
 148     exx[i] = (-ptmin+ptmax)/2;
 149
 150     Float_t pp=0.1;
 151     if(prob < 0.2) pp = 0.;
 152     if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
 153     else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
 154     sprintf(name,"TOF matched: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
 155     hh->SetTitle(name);
 156     sprintf(name,"hNoPid%i",i);
 157
 158     pp=prob;
 159     if(prob < 0.2) pp = 0.1;
 160     if(pos) hh2=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
 161     else hh2=GetHistoPin(ptmin,ptmax,pp,0.0);
 162     AddHisto(hh,hh2,weightS);
 163
 164     h = hh->ProjectionX(name,cmin,cmax);
 165     h->RebinX(rebinsize);
 166     h->Draw("ERR");
 167     h->SetMarkerStyle(24);
 168     b[i][0]=-1;
 169     Int_t ntrial = 0;
 170     Float_t chi2 = 10000;
 171     while(ntrial < 10 && (chi2 > 20 + 1000*selectTrue)){
 172       fit(h,b[i],"WW","",xx[i]);
 173       c1->Update();
 174 //       getchar();
 175       fit(h,b[i],"","",xx[i]);
 176       ntrial++;
 177       chi2 = b[i][2];
 178       printf("chi2 = %f\n",chi2);
 179       c1->Update();
 180 //       getchar();
 181
 182     }
 183
 184     yy[i] = fall->GetParameter(parplotted);
 185     eyy[i] = fall->GetParError(parplotted);
 186   }
 187
 188   TGraphErrors *gpar = new TGraphErrors(nptbin,xx,yy,exx,eyy);
 189   c->cd(8);
 190   gpar->Draw("AP");
 191   gpar->SetMarkerStyle(20);
 192
 193   TCanvas *c2 = new TCanvas();
 194   c2->Divide((nptbin+1)/2,2);
 195   Float_t b2[50][3];
 196
 197   for(Int_t i=0;i < nptbin;i++){
 198     c2->cd(i+1);
 199     Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
 200     Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
 201
 202     Float_t pp=prob;
 203     if(prob < 0.2) pp = 0.1;
 204     if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
 205     else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
 206     sprintf(name,"P_{TOF} > 0.8: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
 207     hh->SetTitle(name);
 208     sprintf(name,"hPid60_%i",i);
 209     h = hh->ProjectionX(name,cmin,cmax);
 210     h->RebinX(rebinsize);
 211     h->Draw("ERR");
 212     h->SetMarkerStyle(24);
 213     b2[i][0]=-1;
 214     Int_t ntrial = 0;
 215     Float_t chi2 = 10000;
 216     while(ntrial < 40 && (chi2 > 20 + 1000*selectTrue)){
 217       fit(h,b2[i],"WW","");
 218       fit(h,b2[i],"","");
 219       ntrial++;
 220       chi2 = b2[i][2];
 221       printf("chi2 = %f\n",chi2);
 222     }
 223     yy[i] = fall->GetParameter(parplotted);
 224     eyy[i] = fall->GetParError(parplotted);
 225
 226   }
 227
 228   TGraphErrors *gpar2 = new TGraphErrors(nptbin,xx,yy,exx,eyy);
 229   c2->cd(8);
 230   gpar2->Draw("AP");
 231   gpar2->SetMarkerStyle(20);
 232
 233   Double_t xpt[50],expt[50],eff[50],efferr[50];
 234   for(Int_t i=0;i<nptbin;i++){
 235     printf("%f +/- %f -  %f +/- %f\n",b[i][0],b[i][1],b2[i][0],b2[i][1]);
 236
 237     Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
 238     Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
 239
 240     xpt[i] = (ptmin+ptmax)/2;
 241     expt[i] = (-ptmin+ptmax)/2;
 242     eff[i] = b2[i][0]/(b[i][0]-b2[i][0]*weightS);
 243
 244     b[i][0] = b[i][0]-b2[i][0]*weightS;
 245
 246     efferr[i] = TMath::Sqrt(b[i][1]*b[i][1]/b[i][0]/b[i][0] + b2[i][1]*b2[i][1]/b2[i][0]/b2[i][0])*(b2[i][0]+b2[i][1])*(1+weightS*(b2[i][0]-b2[i][1])/b[i][0])/b[i][0];//*(1-eff[i]);//der2*der2*(b[i][1]*b[i][1] - b2[i][1]*b2[i][1]));
 247
 248     if(TMath::Abs(efferr[i]) > 1)efferr[i]=1;
 249   }
 250   new TCanvas();
 251   TGraphErrors *geff = new TGraphErrors(nptbin,xpt,eff,expt,efferr);
 252   geff->Draw("AP");
 253
 254   char flag[100];
 255   sprintf(flag,"");
 256
 257   if(isMC){
 258     if(selectTrue) sprintf(flag,"true");
 259     else if(!keepTrue) sprintf(flag,"back");
 260   }
 261
 262   char flag2[100];
 263   sprintf(flag2,"");
 264
 265   char etarange[100];
 266   sprintf(etarange,"_%.1f-%.1f_",etaminkp,etamaxkp);
 267
 268   if(kGoodMatch)
 269     sprintf(flag2,"GM");
 270
 271   if(bayesVsigma)
 272     sprintf(flag2,"BayesVsSigma");
 273
 274   if(kSigma2vs3)
 275     sprintf(flag2,"Sigma2vs3");
 276
 277   if(kOverAll)
 278     sprintf(flag2,"OverAll");
 279   if(kOverAllTOFmatch)
 280     sprintf(flag2,"OverAllTOF");
 281   if(kOverAll2Sigma)
 282     sprintf(flag2,"OverAll2sigma");
 283
 284   if(pos){
 285     if(prob >=0.2) sprintf(name,"protonPos%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
 286     else sprintf(name,"protonPos%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
 287   }
 288   else{
 289     if(prob >=0.2) sprintf(name,"protonNeg%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
 290     else sprintf(name,"protonNeg%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
 291   }
 292
 293   TFile *fout = new TFile(name,"RECREATE");
 294   geff->Write();
 295   fout->Close();
 296
 297
 298   TF1 *ff = new TF1("ff","[0] - [1]*TMath::Exp([2]*x)",0,3);
 299   ff->SetParameter(0,0.67);
 300   ff->SetParameter(1,1.14383e+00);
 301   ff->SetParameter(2,-2.29910);
 302   ff->SetLineColor(4);
 303   ff->SetLineColor(2);
 304   ff->Draw("SAME");
 305
 306   TF1 *ff2 = new TF1("ff2","[0] - [1]*TMath::Exp([2]*x)",0,3);
 307   ff2->SetParameter(0,0.67);
 308   ff2->SetParameter(1,9.23126e-01);
 309   ff2->SetParameter(2,-1.851);
 310   ff2->SetLineColor(4);
 311   ff2->Draw("SAME");
 312 }
 313
 314 TH2F *GetHistoPip(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkp=0,Float_t pMinkn=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
 315
 316   Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,pt+0.001,xmin[3]+0.001,pMinkp+0.001,pMinkn+0.001,(pMinkp>0.09)+0.001,kTOFmatch+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
 317   Float_t x2[] = {xmax[0],etamaxkp-0.001,ptM-0.001,xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
 318
 319   if(kOverAll){
 320     x[6] = 0.0001;
 321     x2[9] = 5.9;
 322     if(pMinkp > 0.19) x2[9] = 4.9;
 323   }
 324
 325   if(kOverAllTOFmatch && pMinkp > 0.19){
 326     x[6] = 1.0001;
 327     x2[9] = 4.9;
 328   }
 329
 330   if(kOverAll2Sigma && pMinkp > 0.09){
 331     x2[9] = 2;
 332     x[6] = 1.0001;
 333   }
 334
 335   if(kGoodMatch){
 336     x[6] = 1.0001;
 337     if(pMinkp > 0)
 338       x2[9] = 4.9;
 339
 340   }
 341
 342   if(kTOFmatch){
 343     x[6] = 1.0001;
 344   }
 345
 346   if(kSigma2vs3){
 347     x[6] = 1.0001;
 348     x2[9] = 3;
 349     if(pMinkp > 0)
 350       x2[9] = 2;
 351   }
 352
 353   if(bayesVsigma){
 354     if(pMinkp > 0){
 355       x[4] = 0.2001;
 356       x2[9] = 5;
 357     }
 358     else{
 359       x2[9] = 3;
 360     }
 361
 362
 363   }
 364
 365   if(require5sigma) x2[9] = 4.9;
 366
 367   AliPIDperfContainer *tmp = fContPid1;
 368
 369   TH2F *h = tmp->GetQA(0, x, x2);
 370
 371   h->GetXaxis()->SetTitle("M_{#Lambda} (GeV/#it{c}^{2})");
 372   h->GetYaxis()->SetTitle("centrality [%]");
 373
 374   return h;
 375 }
 376
 377 TH2F *GetHistoPin(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkn=0,Float_t pMinkp=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8){
 378
 379   Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,xmin[2]+0.001,pt+0.001,pMinkp+0.001,pMinkn+0.001,kTOFmatch+0.001,(pMinkn>0.09)+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
 380   Float_t x2[] = {xmax[0],etamaxkp-0.001,xmax[2],ptM-0.001,xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
 381
 382   if(kOverAll){
 383     x[7] = 0.0001;
 384     x2[10] = 5.9;
 385     if(pMinkn > 0.19) x2[10] = 4.9;
 386   }
 387
 388   if(kOverAllTOFmatch && pMinkn > 0.19){
 389     x[7] = 1.0001;
 390     x2[10] = 4.9;
 391   }
 392
 393   if(kOverAll2Sigma && pMinkn > 0.09){
 394     x2[10] = 2;
 395     x[7] = 1.0001;
 396   }
 397
 398   if(kGoodMatch){
 399     x[7] = 1.0001;
 400     if(pMinkn > 0)
 401       x2[10] = 4.9;
 402
 403   }
 404
 405   if(kTOFmatch){
 406     x[7] = 1.0001;
 407   }
 408
 409   if(kSigma2vs3){
 410     x[7] = 1.0001;
 411     x2[10] = 3;
 412     if(pMinkn > 0)
 413       x2[10] = 2;
 414   }
 415
 416  if(bayesVsigma){
 417     if(pMinkn > 0){
 418       x[5] = 0.2001;
 419       x2[10] = 5;
 420     }
 421     else{
 422       x2[10] = 3;
 423     }
 424   }
 425
 426   if(require5sigma) x2[10] = 4.9;
 427
 428   AliPIDperfContainer *tmp = fContPid2;
 429
 430   TH2F *h = tmp->GetQA(0, x, x2);
 431
 432   h->GetXaxis()->SetTitle("M_{#Lambda} (GeV/#it{c}^{2})");
 433   h->GetYaxis()->SetTitle("centrality [%]");
 434
 435   return h;
 436 }
 437
 438 fit(TH1D *h,Float_t *a=NULL,char *opt="",char *opt2="",Float_t pt=1.5){
 439   if(h->Integral(1,h->GetNbinsX()) < 1){
 440     if(a){
 441       a[0]=0.01;
 442       a[1]=1;
 443     }
 444     return;
 445   }
 446
 447
 448  fall->SetParameter(0,100);
 449  fall->SetParameter(1,1.115);
 450  fall->SetParameter(2,0.0003);
 451  fall->SetParameter(3,0.001);
 452
 453  fall->SetParLimits(0,0.00001,10000);
 454  fall->SetParLimits(1,1.105,1.125);//1.01898 + 2.4e-04*pt-1e-03,1.01898 + 2.4e-04*pt+1e-03);
 455  fall->SetParLimits(2,0.00001,0.001);
 456  fall->SetParLimits(3,0.0005,0.01);
 457
 458  // fall->FixParameter(2,5E-4);
 459
 460 //  fall->FixParameter(1,1.01884 + 2.9891e-04*pt);
 461 //  fall->FixParameter(2,0.0044);
 462 //  fall->FixParameter(3,7.57574e-04 + 3.85408e-04*pt);
 463
 464  fall->ReleaseParameter(4);
 465  fall->ReleaseParameter(5);
 466  fall->ReleaseParameter(6);
 467
 468  if(selectTrue){
 469    fall->FixParameter(4,0);
 470    fall->FixParameter(5,0);
 471    fall->FixParameter(6,0);
 472  }
 473
 474  char name[100];
 475  TF1 *ftmp=fall;
 476
 477  TF1 *ftmp2=new TF1(*fsign);
 478  sprintf(name,"fsign%i",ifunc);
 479  ftmp2->SetName(name);
 480
 481  TF1 *ftmp3=new TF1(*fback);
 482  sprintf(name,"ftmp3%i",ifunc);
 483  ftmp3->SetName(name);
 484
 485  ifunc++;
 486
 487  h->Fit(ftmp,opt,opt2,fitmin,fitmax);
 488  h->Draw("ERR");
 489
 490  ftmp2->SetParameter(0,ftmp->GetParameter(0));
 491  ftmp2->SetParameter(1,ftmp->GetParameter(1));
 492  ftmp2->SetParameter(2,ftmp->GetParameter(2));
 493  ftmp2->SetParameter(3,ftmp->GetParameter(3));
 494  ftmp2->Draw("SAME");
 495  ftmp3->SetParameter(0,ftmp->GetParameter(4));
 496  ftmp3->SetParameter(1,ftmp->GetParameter(5));
 497  ftmp3->SetParameter(2,ftmp->GetParameter(6));
 498  ftmp3->Draw("SAME");
 499
 500  Float_t mean = ftmp->GetParameter(1);
 501  Float_t sigma = TMath::Abs(ftmp->GetParameter(3)) + TMath::Abs(ftmp->GetParameter(2));
 502
 503  Float_t signI = ftmp2->Integral(mean-10*sigma,mean+10*sigma)/h->GetBinWidth(1);
 504  Float_t backI = ftmp3->Integral(mean-3*sigma,mean+3*sigma)/h->GetBinWidth(1);
 505
 506  Float_t errI = TMath::Sqrt(ftmp->GetParError(0)*ftmp->GetParError(0)/(0.001+ftmp->GetParameter(0))/(0.001+ftmp->GetParameter(0)));
 507
 508  printf("signal(5 sigma) = %f +/- %f(fit) +/- %f(stat)\n",signI,errI*signI,TMath::Sqrt(signI));
 509  printf("backgr(3sigma) = %f\n",backI);
 510  printf("significance(3 sigma) = %f\n",signI/sqrt(signI+backI));
 511
 512  if(a){
 513    a[0]=signI;
 514    a[1]=signI*errI*signI*errI + signI;
 515    a[1] = TMath::Sqrt(a[1]);
 516    if(ftmp->GetNDF()) a[2] = ftmp->GetChisquare()/ftmp->GetNDF();
 517
 518
 519     if(selectTrue){
 520       a[0] = h->Integral(1,h->GetNbinsX());
 521       a[1] = TMath::Sqrt(a[0]);
 522     }
 523  }
 524 }
 525
 526 AddHisto(TH2F *h1,TH2F *h2,Float_t w){
 527   Int_t nbinx = h1->GetNbinsX();
 528   Int_t nbiny = h1->GetNbinsY();
 529
 530   for(Int_t i=1;i<=nbinx;i++){
 531     for(Int_t j=1;j<=nbiny;j++){
 532       Float_t val = h1->GetBinContent(i,j) + h2->GetBinContent(i,j)*w;
 533       Float_t err = TMath::Min(TMath::Sqrt(val),val);
 534       h1->SetBinContent(i,j,val);
 535       h1->SetBinError(i,j,err);
 536     }
 537   }
 538 }