/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //------------------------------------------------------------------------------ // Implementation of AliPerfAnalyzeInvPt class. It analyzes the output of // AliPerformancePtCalib.cxx and AliPerformancePtCalibMC.cxx: via // AliPerformancePtCalib::Analyse() or AliPerformancePtCalibMC::Analyse() // // Projection of charge/pt vs theta and vs phi resp. histoprams will be fitted with either // polynomial or gaussian fit function to extract minimum position of charge/pt. // Fit options and theta, phi bins can be set by user. // Attention: use the Set* functions of AliPerformancePtCalib.h and AliPerformancePtCalibMC.h // when running AliPerformancePtCalib*::Analyse() // The result of the analysis (histograms/graphs) are stored in the folder which is // a data member of AliPerformancePtCalib*. // // Author: S.Schuchmann 11/13/2009 //------------------------------------------------------------------------------ /* // after running the performance task of AliPerformancePtCalib* or AliPerformancePtCalibMC* // read the file, and get component: TFile f("Output.root"); AliPerformancePtCalib *compObj = (AliPerformancePtCalib*)coutput->FindObject("AliPerformancePtCalib"); //or AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC"); // set phi and theta bins for fitting and analyse comparison data compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi); compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta); compObj->SetMakeFitOption(kFALSE,exclRange,fitRange); compObj->SetDoRebin(rebin); compObj->Analyse(); //for details see functions of class AliPerformancePtCalib* // the output histograms/graphs will be stored in the folder "folderRes" compObj->GetAnalysisFolder()->ls("*"); // user can save whole comparison object (or only folder with anlysed histograms) // in the seperate output file (e.g.) TFile fout("Analysed_InvPt.root","recreate"); compObj->Write(); // compObj->GetAnalysisFolder()->Write(); fout.Close(); */ #include "TNamed.h" #include "TF1.h" #include "TH1F.h" #include "TH2F.h" #include "TH1D.h" #include "TMath.h" #include "TGraphErrors.h" #include "TCanvas.h" #include "TObjArray.h" #include "AliPerfAnalyzeInvPt.h" ClassImp(AliPerfAnalyzeInvPt) // fit functions //____________________________________________________________________________________________________________________________________________ Double_t AliPerfAnalyzeInvPt::Polynomial(Double_t *x, const Double_t *par) { // fit function for fitting of 1/pt with a polynomial of 4th order, rejecting points between +/-par[4] if (x[0] > -par[4] && x[0] < par[4]) { TF1::RejectPoint(); return 0; } return par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4); } //____________________________________________________________________________________________________________________________________________ Double_t AliPerfAnalyzeInvPt::PolynomialRejP(Double_t *x, const Double_t *par) { // fit function for fitting of 1/pt with a polynomial of 4th order to improve result (fit range and rejection zone is adjusted to first guess of minimum position) Double_t pos = par[5]; Double_t neg = - par[5]; pos += par[4]; neg += par[4]; if (x[0] > neg && x[0] < pos) { TF1::RejectPoint(); return 0; } return par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4); } //____________________________________________________________________________________________________________________________________________ Double_t AliPerfAnalyzeInvPt::InvGauss(Double_t *x, const Double_t *par) { // fit function for fitting of 1/pt with gaussian, rejecting points between +/-par[4] if (x[0] > -par[6] && x[0] < par[6]) { TF1::RejectPoint(); return 0; } 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)) ; } //____________________________________________________________________________________________________________________________________________ Double_t AliPerfAnalyzeInvPt::InvGaussRejP(Double_t *x, const Double_t *par) { // fit function for fitting of 1/pt with gaussian to improve result (fit range and rejection zone is adjusted to first guess of minimum position) Double_t pos = par[7];//0.12; Double_t neg = - par[7];//0.12; pos += par[6]; neg += par[6]; if (x[0] > neg && x[0] < pos) { TF1::RejectPoint(); return 0; } 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) ; } //_____________________________________________________________________________________________________________________________________________ AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt(const Char_t* name, const Char_t* title): TNamed(name,title), fThetaBins(), fPhiBins(), fNThetaBins(0), fNPhiBins(0), fRange(0), fExclRange(0), fFitGaus(0) , fDoRebin(0), fRebin(0), fHistFitTheta(), fHistFitPhi(), //histograms,graphs and functions fHistH2InvPtTheta(0), fHistH2InvPtPhi(0), fGrMinPosTheta(0), fGrMinPosPhi(0), fFitMinPos(0), fFitMinPosRejP(0), fFitInvGauss(0), fFitInvGaussRejP(0) { // Default constructor fFitGaus = kFALSE;// flag for gaussian fit fNThetaBins = 0; // theta bins for projections fNPhiBins = 0; //phi bins for projections fRange = 0; //fit range fExclRange = 0; //fit exclusion range fDoRebin = kFALSE; // flag for rebin fRebin = 0; // bins for rebin for(Int_t i=0;i<100;i++){ fHistFitTheta[i] = NULL; fHistFitPhi[i] = NULL; } } //______________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::InitGraphs(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi){ // initialize graphs for storing fit parameter fGrMinPosTheta = new TGraphErrors(fNThetaBins,binsXTheta,fitParamTheta,0, errFitParamTheta); fGrMinPosTheta->SetMarkerStyle(20); fGrMinPosTheta->SetMarkerColor(2); fGrMinPosTheta->SetLineColor(2); fGrMinPosTheta->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}"); fGrMinPosTheta->GetXaxis()->SetTitle("#theta bin no."); fGrMinPosTheta->GetYaxis()->SetTitleOffset(1.2); fGrMinPosTheta->SetTitle("#theta bins"); fGrMinPosPhi = new TGraphErrors(fNPhiBins,binsXPhi,fitParamPhi,0,errFitParamPhi); fGrMinPosPhi->SetMarkerStyle(20); fGrMinPosPhi->SetMarkerColor(4); fGrMinPosPhi->SetLineColor(4); fGrMinPosPhi->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}"); fGrMinPosPhi->GetXaxis()->SetTitle("#phi bin no."); fGrMinPosPhi->GetYaxis()->SetTitleOffset(1.2); fGrMinPosPhi->SetTitle("#phi bins"); } //______________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::InitFitFcn(){ // fit functions fFitMinPos = new TF1("fFitMinPos", Polynomial,-4.0,4.0,5); fFitMinPos->SetLineColor(4); fFitMinPos->SetLineWidth(1); fFitMinPos->SetParameter(0,1.0); fFitMinPos->SetParameter(1,0.0); fFitMinPos->SetParameter(2,0.0); fFitMinPosRejP = new TF1("fFitMinPosRejP",PolynomialRejP,-4.0,4.0,6); fFitMinPosRejP->SetLineColor(2); fFitMinPosRejP->SetLineWidth(1); fFitMinPosRejP->SetParameter(0,1.0); fFitMinPosRejP->SetParameter(1,0.0); fFitMinPosRejP->SetParameter(2,0.0); fFitInvGauss = new TF1("fFitInvGauss", InvGauss,-4.0,4.0,6); fFitInvGauss->SetLineColor(4); fFitInvGauss->SetLineWidth(1); fFitInvGauss->SetParameter(2,1.0); fFitInvGaussRejP = new TF1("fFitInvGaussRejP", InvGaussRejP,-4.0,4.0,7); fFitInvGaussRejP->SetLineColor(2); fFitInvGaussRejP->SetLineWidth(1); fFitInvGaussRejP->SetParameter(2,1.0); } //______________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::StartAnalysis(const TH2F *histThetaInvPt, const TH2F *histPhiInvPt, TObjArray* aFolderObj){ //start analysis: fitting charge/pt spectra if(!histThetaInvPt) { Printf("warning: no 1/pt histogram to analyse in theta bins!"); return; } if(!histPhiInvPt) { Printf("warning: no 1/pt histogram to analyse in phit bins!"); return; } Double_t thetaBins[9] = {0.77,0.97,1.17,1.37,1.57,1.77,1.97,2.17,2.37}; // default theta bins 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}; // default phi bins Int_t nThetaBins = 9; Int_t nPhiBins = 13; Double_t range = 1.0; // fit range Double_t exclRange = 0.13; // range of point rejection in fit Bool_t fitGaus=kFALSE; // fit with Gaussian or with polynomial (kFALSE) if(!fNThetaBins){ fNThetaBins = nThetaBins; for(Int_t k = 0;kClone("invPtVsTheta"); fHistH2InvPtPhi = (TH2F*)histPhiInvPt->Clone("invPtVsPhi"); // initialize fit functions InitFitFcn(); // canvas for 2D histograms (input) TCanvas *thephiCan =new TCanvas("thephiCan","theta and phi vs invPt",800,400); thephiCan->Divide(2,1); thephiCan->cd(1); fHistH2InvPtTheta->SetTitle("#theta vs charge/pt for selected #phi range"); fHistH2InvPtTheta->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}"); fHistH2InvPtTheta->GetYaxis()->SetTitle("#theta (rad)"); fHistH2InvPtTheta->DrawCopy("col1z"); thephiCan->cd(2); fHistH2InvPtPhi->SetTitle("#phi vs charge/pt for selected #theta range"); fHistH2InvPtPhi->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}"); fHistH2InvPtPhi->GetYaxis()->SetTitle("#phi (rad)"); fHistH2InvPtPhi->DrawCopy("col1z"); // canvas for 1D histograms (projections) TCanvas *theCan =new TCanvas("theCan","invPt theta bins",1200,900); theCan->Divide((fNThetaBins+2)/2,2); TCanvas *phiCan =new TCanvas("phiCan","invPt phi bins",1200,900); phiCan->Divide((fNPhiBins+2)/2,2); Int_t countPad = 1; // analyse charge/pt in bins of theta for(Int_t i=0;iSetName(name.Data()); //find bins for projection Int_t firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]); if(i>0) firstBin +=1; Int_t lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i+1]); if( i == fNThetaBins-1) { firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[0]); lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]); } //make projection fHistFitTheta[i] = (TH1D*)fHistH2InvPtTheta->ProjectionX("_px",firstBin,lastBin,"e"); Char_t titleTheta[50]; if(TMath::Abs(i-(fNThetaBins-1))<0.5) snprintf(titleTheta,50,"charge/pt (GeV/c) integrated over #theta"); else snprintf(titleTheta,50,"charge/pt (GeV/c) for #theta range: %1.3f - %1.3f",fThetaBins[i],fThetaBins[i+1]); fHistFitTheta[i]->SetTitle(titleTheta); Double_t invPtMinPos = 0; Double_t invPtMinPosErr = 0; Double_t invPtMinPosImpr = 0; Double_t invPtMinPosErrImpr = 0; if(fDoRebin) fHistFitTheta[i]->Rebin(fRebin); //start fitting if(!fFitGaus){ Printf("making polynomial fit in charge/pt in theta bins"); theCan->cd(countPad); MakeFit(fHistFitTheta[i],fFitMinPos, invPtMinPos,invPtMinPosErr, fExclRange,fRange); MakeFitBetter(fHistFitTheta[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos,fExclRange,fRange); //plot projection and fit fHistFitTheta[i]->DrawCopy(); fFitMinPos->DrawCopy("L,same"); fFitMinPosRejP->DrawCopy("L,same"); } else{ Printf("making gauss fit in charge/pt in theta bins"); theCan->cd(countPad); MakeFitInvGauss(fHistFitTheta[i],fFitInvGauss, invPtMinPos, invPtMinPosErr,fExclRange,fRange); MakeFitInvGaussBetter(fHistFitTheta[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange); //plot projection and fit fHistFitTheta[i]->DrawCopy(); fFitInvGauss->DrawCopy("L,same"); fFitInvGaussRejP->DrawCopy("L,same"); } //add objects for analysis folder aFolderObj->Add(fHistFitTheta[i]); //store fit parameter fitParamTheta[i] = invPtMinPosImpr; errFitParamTheta[i] = invPtMinPosErrImpr; //count bins and pad number for displaying binsXTheta[i] = i+1.0; countPad++; } countPad = 1; // analyse charge/pt in bins of phi for(Int_t i=0;iSetName(name.Data()); //find bins for projection Int_t firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]); if(i>0) firstBin +=1; Int_t lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i+1]); if(TMath::Abs(i-(fNPhiBins-1))<0.5){ firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[0]); lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]); } //make projection fHistFitPhi[i] = (TH1D*) fHistH2InvPtPhi->ProjectionX("_px",firstBin,lastBin,"e"); Char_t titlePhi[50]; if(TMath::Abs(i-(fNPhiBins-1))<0.5) snprintf(titlePhi,50,"charge/pt (GeV/c) integrated over #phi"); else snprintf(titlePhi,50,"charge/pt (GeV/c) for #phi range: %1.3f - %1.3f",fPhiBins[i],fPhiBins[i+1]); fHistFitPhi[i]->SetTitle(titlePhi); Double_t invPtMinPos = 0; Double_t invPtMinPosErr = 0; Double_t invPtMinPosImpr = 0; Double_t invPtMinPosErrImpr = 0; if(fDoRebin) fHistFitPhi[i]->Rebin(fRebin); //start fitting if(!fFitGaus){ Printf("making polynomial fit in charge/pt in phi bins"); phiCan->cd(countPad); MakeFit(fHistFitPhi[i],fFitMinPos, invPtMinPos, invPtMinPosErr,fExclRange,fRange); MakeFitBetter(fHistFitPhi[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange); //plot projection and fit fHistFitPhi[i]->DrawCopy(); fFitMinPos->DrawCopy("L,same"); fFitMinPosRejP->DrawCopy("L,same"); } else { Printf("making gauss fit in charge/pt in phi bins"); phiCan->cd(countPad); MakeFitInvGauss(fHistFitPhi[i],fFitInvGauss, invPtMinPos, invPtMinPosErr, exclRange,fRange); MakeFitInvGaussBetter(fHistFitPhi[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos, fExclRange,fRange); //plot projection and fit fHistFitPhi[i]->DrawCopy(); fFitInvGauss->DrawCopy("L,same"); fFitInvGaussRejP->DrawCopy("L,same"); } //add objects for analysis folder aFolderObj->Add(fHistFitPhi[i]); //store fit parameter fitParamPhi[i] = invPtMinPosImpr; errFitParamPhi[i] = invPtMinPosErrImpr; //count bins and pad number for displaying binsXPhi[i] = i+1.0; countPad++; } //initialize graphs for displaying the fit parameter InitGraphs(binsXTheta,fitParamTheta,errFitParamTheta,binsXPhi,fitParamPhi,errFitParamPhi); //plot fit values = minimum positions of charge/pt TCanvas *canFitVal = new TCanvas("canFitVal","min pos histos",800,400); canFitVal->Divide(2,1); canFitVal->cd(1); fGrMinPosTheta->Draw("AP"); canFitVal->cd(2); fGrMinPosPhi->Draw("AP"); Printf("AliPerfAnalyzeInvPt::StartAnalysis: ***** NOTE: last bin is always fit result of integral over all angle BINS which have been set by user! angle ranges of 2D projections are kept.*****"); //add objects to folder aFolderObj->Add(fGrMinPosTheta); aFolderObj->Add(fGrMinPosPhi); aFolderObj->Add(fHistH2InvPtTheta); aFolderObj->Add(fHistH2InvPtPhi); } //____________________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::MakeFit(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range) { // fit charge/pt and extract minimum position fitpb->SetRange(-range,range); fitpb->SetParLimits(1,-0.05,0.05); fitpb->SetParameter(0,1.0); fitpb->FixParameter(4,excl); fitpb->FixParameter(2,0.0); hproy->Fit(fitpb,"RM"); mean = fitpb->GetParameter(1); errMean = fitpb->GetParError(1); if(!mean) errMean = 0.0; } //____________________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::MakeFitBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range) { // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFit and fit 1/pt and extract new minimum position fitpb2->FixParameter(5,excl); fitpb2->FixParameter(4,f); // fitpb2->FixParameter(2,0.0); fitpb2->SetRange(-range+f,range+f); fitpb2->SetParLimits(1,-0.05,0.05); fitpb2->SetParameter(0,1.0); hproy->Fit(fitpb2,"RM"); mean = fitpb2->GetParameter(1); errMean = fitpb2->GetParError(1); if(!mean) errMean = 0.0; } //____________________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::MakeFitInvGauss(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range) { // fit charge/pt and extract minimum position fitpb->FixParameter(6,excl); fitpb->SetRange(-range,range); fitpb->SetParameter(0,-1.0); fitpb->SetParameter(2,1.0); fitpb->SetParameter(3,25000.0); fitpb->SetParameter(4,-1.0); fitpb->SetParLimits(1,-0.02,0.02); hproy->Fit(fitpb,"RM"); mean = fitpb->GetParameter(1); errMean = fitpb->GetParError(1); if(!mean) errMean = 0.0; } //____________________________________________________________________________________________________________________________________________ void AliPerfAnalyzeInvPt::MakeFitInvGaussBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range) { // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFitInvGauss and fit 1/pt and extract new minimum position fitpb2->FixParameter(7,excl); fitpb2->FixParameter(6,f); fitpb2->SetRange(-range+f,range+f); fitpb2->SetParameter(0,-1.0); fitpb2->SetParameter(2,1.0); fitpb2->SetParameter(3,25000.0); fitpb2->SetParameter(4,-1.0); fitpb2->SetParLimits(1,-0.02,0.02); hproy->Fit(fitpb2,"RM"); mean = fitpb2->GetParameter(1); errMean = fitpb2->GetParError(1); if(!mean) errMean = 0.0; } //____________________________________________________________________________________________________________________________________________ // set variables for analysis void AliPerfAnalyzeInvPt::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins){ // set theta bins for Analyse() //set phi bins as array and set number of this array which is equal to number of bins analysed //the last analysed bin will always be the projection from first to last bin in the array if(nphBins){ fNPhiBins = nphBins; for(Int_t k = 0;k