MomentumScaleCalibration/test/fit2DProj.c

0001 #include "TH1.h"
0002 #include "TH2.h"
0003 #include "TF1.h"
0004 #include "TGraphErrors.h"
0005 #include "TLegend.h"
0006 #include "TCanvas.h"
0007 #include "TSystem.h"
0008 #include "TFile.h"
0009 #include <vector>
0010 #include <iostream>
0011
0012 void fit2DProj(TString name, TString path, int minEntries, int rebin, int fitType);
0013
0014 Double_t gaussian(Double_t *x, Double_t *par);
0015 Double_t lorentzian(Double_t *x, Double_t *par);
0016 Double_t lorentzianPlusLinear(Double_t *x, Double_t *par);
0017
0018 TF1* gaussianFit(TH1* histoY);
0019 TF1* lorentzianFit(TH1* histoY);
0020 TF1* linLorentzianFit(TH1* histoY);
0021
0022
0023 //
0024 // Gaussian function
0025 // -----------------------
0026 Double_t gaussian(Double_t *x, Double_t *par) {
0027   return par[0]*exp(-0.5*((x[0]-par[2])/par[1])*((x[0]-par[2])/par[1])) ; }
0028 //
0029 // Lorentzian function
0030 // -----------------------
0031 Double_t lorentzian(Double_t *x, Double_t *par) {
0032   return (0.5*par[0]*par[1]/TMath::Pi()) /
0033     TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
0034 }
0035 //
0036 // Linear + Lorentzian function
0037 // -----------------------
0038 Double_t lorentzianPlusLinear(Double_t *x, Double_t *par) {
0039   return  ((0.5*par[0]*par[1]/TMath::Pi())/
0040     TMath::Max(1.e-10,((x[0]-par[2])*(x[0]-par[2])+.25*par[1]*par[1])))+par[3]+par[4]*x[0];
0041 }
0042
0043 void fit2DProj(TString name, TString path, int minEntries, int rebin, int fitType) {
0044
0045   //Read the TH2 from file
0046   TFile *inputFile = new TFile(path);
0047   TH2 * histo = (TH2*) inputFile->Get(name);;
0048   histo->RebinX(rebin);
0049
0050   //Declare some variables
0051   TH1 * histoY;
0052   TString nameY;
0053   vector<double> Ftop;
0054   vector<double> Fwidth;
0055   vector<double> Fmass;
0056   vector<double> Etop;
0057   vector<double> Ewidth;
0058   vector<double> Emass;
0059   vector<double> Fchi2;
0060   vector<double> Xcenter;
0061   vector<double> Ex;
0062
0063   //TFile *fileOut=new TFile("fitCompare2.root","update");
0064
0065   for (int i=1; i<=(histo->GetXaxis()->GetNbins());i++) {
0066
0067     //Project on Y (set name and title)
0068     char str[100];
0069     sprintf(str,"%i",i);
0070     TString iS(str);
0071     nameY = name + "_"+ iS;
0072     std::cout << "nameY  " << nameY << std::endl;
0073     histoY = histo->ProjectionY(nameY, i, i);
0074     //histoY->Rebin(2);
0075     char title[80];
0076     double x= histo->GetXaxis()-> GetBinCenter(i);
0077     //sprintf(title,"Projection of bin=%i",i);
0078     sprintf(title,"Projection of x=%f",x);
0079     histoY->SetTitle(title);
0080
0081     if (histoY ->GetEntries() > minEntries) {
0082       //Make the dirty work!
0083       TF1 *fit;
0084       if(fitType == 1){
0085     fit = gaussianFit(histoY);
0086       }
0087       else if(fitType == 2){
0088     fit = lorentzianFit(histoY);
0089       }
0090       else if(fitType == 3)
0091     fit = linLorentzianFit(histoY);
0092       else {
0093     std::cout<<"Wrong fit type: 1=gaussian, 2=lorentzian, 3=lorentzian+linear."<<std::endl;
0094     abort();
0095       }
0096
0097       double *par = fit->GetParameters();
0098       double *err = fit->GetParErrors();
0099
0100       if(par[2]>150 || par[2]<50 || err[2]>5)
0101     continue;
0102       //Only for check
0103       TCanvas *c = new TCanvas(nameY, nameY);
0104       histoY->Draw();
0105       fit->Draw("same");
0106       //fileOut->cd();
0107       //c->Write();
0108
0109       //Store the fit results
0110       Ftop.push_back(par[0]);
0111       Fwidth.push_back(fabs(par[1]));//sometimes the gaussian has negative width (checked with Rene Brun)
0112       Fmass.push_back(par[2]);
0113       Etop.push_back(err[0]);
0114       Ewidth.push_back(err[1]);
0115       Emass.push_back(err[2]);
0116
0117       Fchi2.push_back(fit->GetChisquare()/fit->GetNDF());
0118
0119       double xx= histo->GetXaxis()->GetBinCenter(i);
0120       Xcenter.push_back(xx);
0121       double ex = 0;
0122       Ex.push_back(ex);
0123     }
0124   }
0125
0126   //fileOut->Close();
0127
0128   //Plots the fit results in  TGraphs
0129   const int nn= Ftop.size();
0130   double x[nn],ym[nn],e[nn],eym[nn];
0131   double yw[nn],eyw[nn],yc[nn];
0132
0133   for (int j=0;j<nn;j++){
0134     x[j]=Xcenter[j];
0135     ym[j]=Fmass[j];
0136     eym[j]=Emass[j];
0137     yw[j]=Fwidth[j];
0138     eyw[j]=Ewidth[j];
0139     yc[j]=Fchi2[j];
0140     e[j]=0;
0141   }
0142
0143   TGraphErrors *grM = new TGraphErrors(nn,x,ym,e,eym);
0144   grM->SetTitle(name+"_M");
0145   grM->SetName(name+"_M");
0146   TGraphErrors *grW = new TGraphErrors(nn,x,yw,e,eyw);
0147   grW->SetTitle(name+"_W");
0148   grW->SetName(name+"_W");
0149   TGraphErrors *grC = new TGraphErrors(nn,x,yc,e,e);
0150   grC->SetTitle(name+"_chi2");
0151   grC->SetName(name+"_chi2");
0152
0153   grM->SetMarkerColor(4);
0154   grM->SetMarkerStyle(20);
0155   grW->SetMarkerColor(4);
0156   grW->SetMarkerStyle(20);
0157   grC->SetMarkerColor(4);
0158   grC->SetMarkerStyle(20);
0159
0160   //Draw the graphs
0161   TCanvas * c1 = new TCanvas(name+"_W",name+"_W");
0162   grW->Draw("AP");
0163   TCanvas * c2 = new TCanvas(name+"_M",name+"_M");
0164   grM->Draw("Ap");
0165   TCanvas * c3 = new TCanvas(name+"_C",name+"_C");
0166   grC->Draw("Ap");
0167
0168   //Save the graphs in a file
0169   TFile *file = new TFile("filegraph.root","update");
0170   grW->Write();
0171   grM->Write();
0172   grC->Write();
0173   file->Close();
0174 }
0175
0176
0177 TF1* gaussianFit(TH1* histoY){
0178   TString name = histoY->GetName() + TString("Fit");
0179
0180     // Fit slices projected along Y from bins in X
0181   TF1 *fit = new TF1(name,gaussian,70,110,3);
0182   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean());
0183   fit->SetParNames("norm","width","mean");
0184   fit->SetLineWidth(2);
0185   histoY -> Fit(name,"0","",85,97);
0186   return fit;
0187 }
0188
0189 TF1* lorentzianFit(TH1* histoY){
0190   TString name = histoY->GetName() + TString("Fit");
0191
0192     // Fit slices projected along Y from bins in X
0193     TF1 *fit = new TF1(name,lorentzian,70,110,3);
0194   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean());
0195   fit->SetParNames("norm","width","mean");
0196   fit->SetLineWidth(2);
0197   histoY -> Fit(name,"0","",85,97);
0198   return fit;
0199 }
0200
0201 TF1* linLorentzianFit(TH1* histoY){
0202   TString name = histoY->GetName() + TString("Fit");
0203
0204   // Fit slices projected along Y from bins in X
0205   TF1 *fit = new TF1(name,lorentzianPlusLinear,70,110,5);
0206   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean(),10,-0.1);
0207   fit->SetParNames("norm","width","mean","offset","slope");
0208   fit->SetParLimits(1,-10,10);
0209   //fit->SetParLimits(2,85,95);
0210   fit->SetParLimits(2,90,93);
0211   fit->SetParLimits(3,0,100);
0212   fit->SetParLimits(4,-1,0);
0213   fit->SetLineWidth(2);
0214   histoY -> Fit(name,"0","",85,97);
0215   return fit;
0216 }
0217
0218