PWG1/TPC/AliPerfAnalyzeInvPt.cxx

   1 #include "TNamed.h"
   2 #include "TList.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TH3F.h"
   6 #include "TCanvas.h"
   7 #include "TMath.h"
   8 #include "TGraphErrors.h"
   9 #include "TF1.h"
  10
  11 //#include "AliPerformancePtCalibMC.h"
  12 //#include "AliPerformancePtCalib.h"
  13 #include "AliPerfAnalyzeInvPt.h"
  14
  15 ClassImp(AliPerfAnalyzeInvPt)
  16
  17
  18 // fit functions
  19 //____________________________________________________________________________________________________________________________________________
  20    Double_t AliPerfAnalyzeInvPt::Polynomial(Double_t *x, Double_t *par)
  21 {
  22
  23
  24    if (x[0] > -par[4] && x[0] < par[4]) {
  25       TF1::RejectPoint();
  26       return 0;
  27    }
  28    return  par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
  29
  30 }
  31 //____________________________________________________________________________________________________________________________________________
  32 Double_t AliPerfAnalyzeInvPt::PolynomialRejP(Double_t *x, Double_t *par)
  33 {
  34
  35    Double_t pos  = par[5];
  36    Double_t neg = - par[5];
  37    pos += par[4];
  38    neg += par[4];
  39
  40    if (x[0] > neg && x[0] < pos) {
  41       TF1::RejectPoint();
  42       return 0;
  43    }
  44    return   par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
  45 }
  46 //____________________________________________________________________________________________________________________________________________
  47 Double_t AliPerfAnalyzeInvPt::InvGauss(Double_t *x, Double_t *par)
  48 {
  49    if (x[0] > -par[6] && x[0] < par[6]) {
  50       TF1::RejectPoint();
  51       return 0;
  52    }
  53
  54    return par[3]+par[0]*TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4) ;
  55 }
  56 //____________________________________________________________________________________________________________________________________________
  57 Double_t AliPerfAnalyzeInvPt::InvGaussRejP(Double_t *x, Double_t *par)
  58 {
  59    Double_t pos  = par[7];//0.12;
  60    Double_t neg = - par[7];//0.12;
  61    pos += par[6];
  62    neg += par[6];
  63
  64    if (x[0] > neg && x[0] < pos) {
  65       TF1::RejectPoint();
  66       return 0;
  67    }
  68
  69
  70    return par[3]+par[0]*TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4) ;
  71 }
  72
  73
  74 //_____________________________________________________________________________________________________________________________________________
  75 AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt():
  76    TNamed("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt"),
  77    fNThetaBins(0),
  78    fNPhiBins(0),
  79    fRange(0),
  80    fExclRange(0),
  81    fFitGaus(0) ,
  82    fHistH2InvPtTheta(0),
  83    fHistH2InvPtPhi(0),
  84    fGrMinPosTheta(0),
  85    fGrMinPosPhi(0),
  86    fFitMinPos(0),
  87    fFitMinPosRejP(0),
  88    fFitInvGauss(0),
  89    fFitInvGaussRejP(0)
  90 {
  91    // Default constructor
  92
  93    fFitGaus = kFALSE;
  94    fNThetaBins = 0;
  95    fNPhiBins = 0;
  96    fRange = 0;
  97    fExclRange = 0;
  98    fFitGaus = 0;
  99
 100    // projection histos
 101    TH1D *fHistFitTheta[100];
 102    TH1D *fHistFitPhi[100];
 103
 104    for(Int_t i=0;i<100;i++){
 105
 106       fHistFitTheta[i] = NULL;
 107       fHistFitPhi[i] = NULL;
 108    }
 109
 110 }
 111 //_____________________________________________________________________________________________________________________________________________
 112 AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt(Char_t* name="AliAnalyzeInvPt",Char_t* title="AliAnalyzeInvPt"): TNamed(name, title),
 113    fNThetaBins(0),
 114    fNPhiBins(0),
 115    fRange(0),
 116    fExclRange(0),
 117    fFitGaus(0) ,
 118    fHistH2InvPtTheta(0),
 119    fHistH2InvPtPhi(0),
 120    fGrMinPosTheta(0),
 121    fGrMinPosPhi(0),
 122    fFitMinPos(0),
 123    fFitMinPosRejP(0),
 124    fFitInvGauss(0),
 125    fFitInvGaussRejP(0)
 126    {
 127       //  Double_t fThetaBins[100] = {0};
 128       //   Double_t fPhiBins[100] = {0};
 129
 130       fFitGaus = kFALSE;
 131       fNThetaBins = 0;
 132       fNPhiBins =0;
 133       fRange = 0;
 134       fExclRange = 0;
 135       fFitGaus = 0;
 136       // projection histos
 137       TH1D *fHistFitTheta[100];
 138       TH1D *fHistFitPhi[100];
 139
 140       for(Int_t i=0;i<100;i++){
 141
 142          fHistFitTheta[i] = NULL;
 143          fHistFitPhi[i] = NULL;
 144       }
 145       fHistH2InvPtTheta = NULL;
 146       fHistH2InvPtPhi = NULL;
 147       fGrMinPosTheta= NULL;
 148       fGrMinPosPhi= NULL;
 149    }
 150
 151
 152    //______________________________________________________________________________________________________________________________________
 153 void AliPerfAnalyzeInvPt::InitHistos(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi){// init Histos
 154
 155
 156
 157
 158       fGrMinPosTheta = new TGraphErrors(fNThetaBins,binsXTheta,fitParamTheta,0, errFitParamTheta);
 159       fGrMinPosTheta->SetMarkerStyle(20);
 160       fGrMinPosTheta->SetMarkerColor(2);
 161       fGrMinPosTheta->SetLineColor(2);
 162       fGrMinPosTheta->GetYaxis()->SetTitle("min pos (Gev/c)^{-1}");
 163       fGrMinPosTheta->GetXaxis()->SetTitle("#theta bin no.");
 164       fGrMinPosTheta->GetYaxis()->SetTitleOffset(1.2);
 165       fGrMinPosTheta->SetTitle("#theta bins ");
 166
 167       fGrMinPosPhi = new TGraphErrors(fNPhiBins,binsXPhi,fitParamPhi,0,errFitParamPhi);
 168       fGrMinPosPhi->SetMarkerStyle(20);
 169       fGrMinPosPhi->SetMarkerColor(4);
 170       fGrMinPosPhi->SetLineColor(4);
 171       fGrMinPosPhi->GetYaxis()->SetTitle("min pos (Gev/c)^{-1}");
 172       fGrMinPosPhi->GetXaxis()->SetTitle("#phi bin no.");
 173       fGrMinPosPhi->GetYaxis()->SetTitleOffset(1.2);
 174       fGrMinPosPhi->SetTitle("#phi bins ");
 175 }
 176
 177 //______________________________________________________________________________________________________________________________________
 178
 179 void AliPerfAnalyzeInvPt::InitFitFcn(){
 180       // fit functions
 181       fFitMinPos = new TF1("fFitMinPos", Polynomial,-4.0,4.0,5);
 182       fFitMinPos->SetLineColor(4);
 183       fFitMinPos->SetLineWidth(1);
 184       fFitMinPos->SetParameter(0,1.0);
 185       fFitMinPos->SetParameter(1,0.0);
 186       fFitMinPos->SetParameter(2,0.0);
 187
 188       fFitMinPosRejP = new TF1("fFitMinPosRejP",PolynomialRejP,-4.0,4.0,6);
 189       fFitMinPosRejP->SetLineColor(2);
 190       fFitMinPosRejP->SetLineWidth(1);
 191       fFitMinPosRejP->SetParameter(0,1.0);
 192       fFitMinPosRejP->SetParameter(1,0.0);
 193       fFitMinPosRejP->SetParameter(2,0.0);
 194
 195
 196       fFitInvGauss = new TF1("fFitInvGauss", InvGauss,-4.0,4.0,6);
 197       fFitInvGauss->SetLineColor(4);
 198       fFitInvGauss->SetLineWidth(1);
 199       fFitInvGauss->SetParameter(2,1.0);
 200
 201       fFitInvGaussRejP = new TF1("fFitInvGaussRejP", InvGaussRejP,-4.0,4.0,7);
 202       fFitInvGaussRejP->SetLineColor(2);
 203       fFitInvGaussRejP->SetLineWidth(1);
 204       fFitInvGaussRejP->SetParameter(2,1.0);
 205
 206    }
 207 //______________________________________________________________________________________________________________________________________
 208 void AliPerfAnalyzeInvPt::StartAnalysis(TH2F *histThetaInvPt, TH2F *histPhiInvPt, TObjArray* aFolderObj){//start Ana
 209
 210
 211
 212    if(!histThetaInvPt) {
 213       Printf("warning: no 1/pt histogram to analyse in theta bins!");
 214    }
 215    if(!histPhiInvPt) {
 216       Printf("warning: no 1/pt histogram to analyse in phit bins!");
 217    }
 218
 219    Double_t thetaBins[9] = {0.77,0.97,1.17,1.37,1.57,1.77,1.97,2.17,2.37};                 // theta bins
 220    Double_t phiBins[13] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.5};           // phi bins
 221
 222    Int_t nThetaBins = 9;
 223    Int_t nPhiBins = 13;
 224    Double_t range = 1.0;        // fit range
 225    Double_t  exclRange  = 0.13; // range of point rejection in fit
 226    Bool_t fitGaus=kFALSE;       // fit with Gaussian or with polynomial (kFALSE)
 227
 228    if(fNThetaBins == 0){
 229       fNThetaBins = nThetaBins;
 230       for(Int_t k = 0;k<fNThetaBins;k++){
 231          fThetaBins[k] = thetaBins[k];
 232       }
 233       Printf("warning: for theta bins no user intput! bins are set to default values.");
 234    }
 235
 236    if(fNPhiBins == 0){
 237       fNPhiBins = nPhiBins;
 238       for(Int_t k = 0;k<fNPhiBins;k++){
 239          fPhiBins[k] = phiBins[k];
 240       }
 241
 242       Printf("warning: for phi bins no user intput! bins are set to default values.");
 243
 244    }
 245
 246    if(fRange==0){
 247       fRange = range;
 248       fExclRange = exclRange;
 249       fFitGaus = fitGaus;
 250       Printf("warning: no fit range is set. default fitting conditions are used.");
 251    }
 252
 253
 254    //param arrays
 255    Double_t fitParamTheta[100],fitParamPhi[100];
 256    Double_t errFitParamTheta[100], errFitParamPhi[100];
 257    Double_t binsXTheta[100],binsXPhi[100];
 258
 259    fHistH2InvPtTheta  = histThetaInvPt;
 260    fHistH2InvPtPhi    = histPhiInvPt;
 261
 262    InitFitFcn();
 263
 264    TCanvas *theCan =new TCanvas("theCan","invPt theta bins",1200,900);
 265    theCan->Divide((fNThetaBins+2)/2,2);
 266    TCanvas *phiCan =new TCanvas("phiCan","invPt phi bins",1200,900);
 267    phiCan->Divide((fNPhiBins+2)/2,2);
 268    Int_t countPad = 1;
 269
 270    // analyse 1/pt in bins of theta
 271    for(Int_t i=0;i<fNThetaBins;i++){
 272       TString name = "fit_theta_";
 273       name +=i;
 274       Int_t firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
 275       if(i>0) firstBin +=1;
 276       Int_t lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i+1]);
 277       if( i == fNThetaBins-1) {
 278          firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[0]);
 279          lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
 280       }
 281
 282       fHistH2InvPtTheta->SetName(name.Data());
 283       fHistFitTheta[i] =  (TH1F*)fHistH2InvPtTheta->ProjectionX("_px",firstBin,lastBin,"e");
 284
 285       Char_t titleTheta[50];
 286       if(i == fNThetaBins-1) sprintf(titleTheta,"1/pt (GeV/c) integrated over #theta");
 287       else  sprintf(titleTheta,"1/pt (GeV/c) for #theta range: %1.3f - %1.3f",fThetaBins[i],fThetaBins[i+1]);
 288
 289       fHistFitTheta[i]->SetTitle(titleTheta);
 290
 291       Double_t invPtMinPos  = 0;
 292       Double_t invPtMinPosErr = 0;
 293       Double_t invPtMinPosImpr  = 0;
 294       Double_t invPtMinPosErrImpr = 0;
 295
 296
 297
 298       if(fFitGaus==kFALSE){
 299          Printf("making polynomial fit in 1/pt in theta bins");
 300          theCan->cd(countPad);
 301          MakeFit(fHistFitTheta[i],fFitMinPos, invPtMinPos,invPtMinPosErr, fExclRange,fRange);
 302          MakeFitBetter(fHistFitTheta[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos,fExclRange,fRange);
 303
 304
 305          fHistFitTheta[i]->DrawCopy();
 306          fFitMinPos->DrawCopy("L,same");
 307          fFitMinPosRejP->DrawCopy("L,same");
 308       }
 309       else{
 310          Printf("making gauss fit in 1/pt in theta bins");
 311          theCan->cd(countPad);
 312          MakeFitInvGauss(fHistFitTheta[i],fFitInvGauss, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
 313          MakeFitInvGaussBetter(fHistFitTheta[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
 314
 315          fHistFitTheta[i]->DrawCopy();
 316          fFitInvGauss->DrawCopy("L,same");
 317          fFitInvGaussRejP->DrawCopy("L,same");
 318       }
 319
 320       aFolderObj->Add(fHistFitTheta[i]);
 321
 322       fitParamTheta[i] = invPtMinPosImpr;
 323       errFitParamTheta[i] = invPtMinPosErrImpr;
 324
 325       binsXTheta[i] = i+1.0;
 326       countPad++;
 327    }
 328
 329
 330    countPad = 1;
 331
 332
 333    // analyse 1/pt in bins of phi
 334
 335    for(Int_t i=0;i<fNPhiBins;i++){
 336       TString name = "fit_phi_";
 337       name +=i;
 338
 339       fHistH2InvPtPhi->SetName(name.Data());
 340       Int_t  firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
 341       if(i>0) firstBin +=1;
 342       Int_t   lastBin =  fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i+1]);
 343       if(i == fNPhiBins-1){
 344          firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[0]);
 345          lastBin =  fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
 346       }
 347       fHistFitPhi[i] =  (TH1F*) fHistH2InvPtPhi->ProjectionX("_px",firstBin,lastBin,"e");
 348
 349       Char_t titlePhi[50];
 350       if(i == fNPhiBins-1) sprintf(titlePhi,"1/pt (GeV/c) integrated over #phi");
 351           else  sprintf(titlePhi,"1/pt (GeV/c) for #phi range: %1.3f - %1.3f",fPhiBins[i],fPhiBins[i+1]);
 352
 353       fHistFitPhi[i]->SetTitle(titlePhi);
 354
 355       Double_t invPtMinPos  = 0;
 356       Double_t invPtMinPosErr = 0;
 357       Double_t invPtMinPosImpr  = 0;
 358       Double_t invPtMinPosErrImpr = 0;
 359
 360       if(fFitGaus==kFALSE){
 361          Printf("making polynomial fit in 1/pt in phi bins");
 362          phiCan->cd(countPad);
 363          MakeFit(fHistFitPhi[i],fFitMinPos, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
 364          MakeFitBetter(fHistFitPhi[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
 365
 366          fHistFitPhi[i]->DrawCopy();
 367          fFitMinPos->DrawCopy("L,same");
 368          fFitMinPosRejP->DrawCopy("L,same");
 369
 370       }
 371       else {
 372          Printf("making gauss fit in 1/pt in phi bins");
 373          phiCan->cd(countPad);
 374          MakeFitInvGauss(fHistFitPhi[i],fFitInvGauss, invPtMinPos, invPtMinPosErr, exclRange,fRange);
 375          MakeFitInvGaussBetter(fHistFitPhi[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos, fExclRange,fRange);
 376
 377          fHistFitPhi[i]->DrawCopy();
 378          fFitInvGauss->DrawCopy("L,same");
 379          fFitInvGaussRejP->DrawCopy("L,same");
 380       }
 381
 382       aFolderObj->Add(fHistFitPhi[i]);
 383
 384       fitParamPhi[i] = invPtMinPosImpr;
 385       errFitParamPhi[i] = invPtMinPosErrImpr;
 386
 387       binsXPhi[i] = i+1.0;
 388       countPad++;
 389    }
 390
 391    InitHistos(binsXTheta,fitParamTheta,errFitParamTheta,binsXPhi,fitParamPhi,errFitParamPhi);
 392
 393    TCanvas *canFitVal = new TCanvas("canFitVal","min pos histos",800,400);
 394    canFitVal->Divide(2,1);
 395
 396    canFitVal->cd(1);
 397    fGrMinPosTheta->Draw("ALP");
 398    canFitVal->cd(2);
 399    fGrMinPosPhi->Draw("ALP");
 400
 401    Printf("AliPerfAnalyzeInvPt: NOTE: last bin is always fit result  of integral over all angle ranges which have been set by user!");
 402
 403    aFolderObj->Add(fGrMinPosTheta);
 404    aFolderObj->Add(fGrMinPosPhi);
 405    aFolderObj->Add(fHistH2InvPtTheta);
 406    aFolderObj->Add(fHistH2InvPtPhi);
 407
 408
 409 }
 410
 411
 412 //____________________________________________________________________________________________________________________________________________
 413
 414 void AliPerfAnalyzeInvPt::MakeFit(TH1 *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
 415 {
 416
 417    fitpb->SetRange(-range,range);
 418    fitpb->SetParLimits(1,-0.05,0.05);
 419    fitpb->SetParameter(0,1.0);
 420    fitpb->FixParameter(4,excl);
 421    fitpb->FixParameter(2,0.0);
 422
 423    hproy->Fit(fitpb,"RQM");
 424    mean = fitpb->GetParameter(1);
 425    errMean = fitpb->GetParError(1);
 426    if(mean == 0)  errMean = 0;
 427 }
 428 //____________________________________________________________________________________________________________________________________________
 429 void AliPerfAnalyzeInvPt::MakeFitBetter(TH1 *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
 430 {
 431    fitpb2->FixParameter(5,excl);
 432    fitpb2->FixParameter(4,f);
 433    // fitpb2->FixParameter(2,0.0);
 434    fitpb2->SetRange(-range+f,range+f);// was 0.25 0.45
 435    fitpb2->SetParLimits(1,-0.05,0.05);
 436    fitpb2->SetParameter(0,1.0);
 437    hproy->Fit(fitpb2,"RQM");
 438    mean = fitpb2->GetParameter(1);
 439    errMean = fitpb2->GetParError(1);
 440    if(mean == 0)   errMean = 0;
 441
 442 }
 443 //____________________________________________________________________________________________________________________________________________
 444 void AliPerfAnalyzeInvPt::MakeFitInvGauss(TH1 *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
 445 {
 446    fitpb->FixParameter(6,excl);
 447    fitpb->SetRange(-range,range);
 448    fitpb->SetParameter(0,-1.0);
 449    fitpb->SetParameter(2,1.0);
 450    fitpb->SetParameter(3,25000.0);
 451    fitpb->SetParameter(4,-1.0);
 452    fitpb->SetParLimits(1,-0.02,0.02);
 453    hproy->Fit(fitpb,"RQM");
 454    mean = fitpb->GetParameter(1);
 455    errMean = fitpb->GetParError(1);
 456    if(mean == 0)   errMean = 0;
 457
 458 }
 459 //____________________________________________________________________________________________________________________________________________
 460 void AliPerfAnalyzeInvPt::MakeFitInvGaussBetter(TH1 *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
 461 {
 462    fitpb2->FixParameter(7,excl);
 463    fitpb2->FixParameter(6,f);
 464    fitpb2->SetRange(-range+f,range+f);
 465    fitpb2->SetParameter(0,-1.0);
 466    fitpb2->SetParameter(2,1.0);
 467    fitpb2->SetParameter(3,25000.0);
 468    fitpb2->SetParameter(4,-1.0);
 469    fitpb2->SetParLimits(1,-0.02,0.02);
 470    hproy->Fit(fitpb2,"RQM");
 471    mean = fitpb2->GetParameter(1);
 472    errMean = fitpb2->GetParError(1);
 473    if(mean == 0)   errMean = 0;
 474
 475 }
 476
 477
 478 //____________________________________________________________________________________________________________________________________________
 479 // set variables
 480
 481
 482 void AliPerfAnalyzeInvPt::SetProjBinsPhi(const Double_t *phiBinArray,Int_t nphBins){
 483
 484    fNPhiBins = nphBins;
 485
 486    for(Int_t k = 0;k<fNPhiBins;k++){
 487       fPhiBins[k] = phiBinArray[k];
 488    }
 489    Printf("number of bins in phi set to %i",fNPhiBins);
 490
 491 }
 492 //____________________________________________________________________________________________________________________________________________
 493 void AliPerfAnalyzeInvPt::SetProjBinsTheta(const Double_t *thetaBinArray, Int_t nthBins){
 494
 495    fNThetaBins = nthBins;
 496    for(Int_t k = 0;k<fNThetaBins;k++){
 497       fThetaBins[k] = thetaBinArray[k];
 498    }
 499    Printf("number of bins in theta set to %i",fNThetaBins);
 500
 501 }
 502 //____________________________________________________________________________________________________________________________________________
 503 void AliPerfAnalyzeInvPt::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
 504
 505
 506
 507    fFitGaus = setGausFit;
 508    fExclRange  = exclusionR;
 509    fRange = fitR;
 510
 511    if(fFitGaus) Printf("set MakeGausFit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
 512    else  Printf("set standard polynomial fit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
 513
 514 }