test/Macros/fit2DProj.C

0001 #include "TROOT.h"
0002 #include "TH1.h"
0003 #include "TH2.h"
0004 #include "TF1.h"
0005 #include "TGraphErrors.h"
0006 #include "TLegend.h"
0007 #include "TCanvas.h"
0008 #include "TSystem.h"
0009 #include "TStyle.h"
0010 #include "TFile.h"
0011 #include "TPaveText.h"
0012 #include <vector>
0013 #include <iostream>
0014 #include <sstream>
0015 #include <cmath>
0016 #include <iomanip>
0017
0018
0019 /**
0020  * Small function to simplify the creation of text in the TPaveText. <br>
0021  * It automatically writes the corret number of figures.
0022  */
0023 TString setText(const char * text, const double & num1, const char * divider = "", const double & num2 = 0)
0024 {
0025   // std::cout << "text = " << text << ", num1 = " << num1 << ", divider = " << divider << ", num2 = " << num2 << std::endl;
0026
0027   // Counter gives the precision
0028   int precision = 1;
0029   int k=1;
0030   while( int(num2*k) == 0 ) {
0031     // std::cout << "int(num2*"<<k<<")/int("<<k<<") = " << int(num2*k)/int(k) << std::endl;
0032     k*=10;
0033     ++precision;
0034   }
0035
0036   std::stringstream numString;
0037   TString textString(text);
0038   numString << std::setprecision(precision) << std::fixed << num1;
0039   textString += numString.str();
0040   if( num2 != 0 ) {
0041     textString += divider;
0042     numString.str("");
0043     if( std::string(text).find("ndf") != std::string::npos ) precision = 0;
0044     numString << std::setprecision(precision) << std::fixed << num2;
0045     textString += numString.str();
0046   }
0047   return textString;
0048 }
0049
0050 /// This function sets up the TPaveText with chi2/ndf and parameters values +- errors
0051 void setTPaveText(const TF1 * fit, TPaveText * paveText) {
0052   Double_t chi2 = fit->GetChisquare();
0053   Int_t ndf = fit->GetNDF();
0054   paveText->AddText(setText("#chi^2 / ndf = ", chi2,  " / ", ndf));
0055
0056   for( int iPar=0; iPar<fit->GetNpar(); ++iPar ) {
0057     TString parName(fit->GetParName(iPar));
0058     Double_t parValue = fit->GetParameter(iPar); // value of Nth parameter
0059     Double_t parError = fit->GetParError(iPar);  // error on Nth parameter
0060     paveText->AddText(setText(parName + " = ", parValue, " #pm ", parError));
0061   }
0062 }
0063
0064 TGraphErrors* fit2DProj(TString name, TString path, int minEntries, int rebinX, int rebinY, int fitType,
0065                         TFile * outputFile, const TString & resonanceType = "Upsilon", const double & xDisplace = 0., const TString & append = "");
0066 void macroPlot( TString name, TString nameGen, const TString & nameFile1 = "0_MuScleFit.root", const TString & nameFile2 = "4_MuScleFit.root",
0067                 const TString & title = "", const TString & resonanceType = "Upsilon", const int rebinX = 0, const int rebinY = 0, const int fitType = 1, const TString & outputFileName = "filegraph.root", const bool genMass = false );
0068
0069 Double_t gaussian(Double_t *x, Double_t *par);
0070 Double_t lorentzian(Double_t *x, Double_t *par);
0071 Double_t lorentzianPlusLinear(Double_t *x, Double_t *par);
0072
0073 Double_t sinusoidal(Double_t *x, Double_t *par);
0074 Double_t parabolic(Double_t *x, Double_t *par);
0075 Double_t parabolic2(Double_t *x, Double_t *par);
0076 Double_t onlyParabolic(Double_t *x, Double_t *par);
0077 Double_t linear(Double_t *x, Double_t *par);
0078 Double_t overX(Double_t *x, Double_t *par);
0079
0080 TF1* gaussianFit(TH1* histoY, const TString & resonanceType = "Upsilon");
0081 TF1* lorentzianFit(TH1* histoY, const TString & resonanceType = "Upsilon");
0082 TF1* linLorentzianFit(TH1* histoY);
0083
0084 void setTDRStyle();
0085 // Gaussian function
0086 Double_t gaussian(Double_t *x, Double_t *par) {
0087   return par[0]*exp(-0.5*((x[0]-par[2])/par[1])*((x[0]-par[2])/par[1]));
0088 }
0089 // Lorentzian function
0090 Double_t lorentzian(Double_t *x, Double_t *par) {
0091   return (0.5*par[0]*par[1]/3.14) /
0092     TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
0093 }
0094 // Linear + Lorentzian function
0095 Double_t lorentzianPlusLinear(Double_t *x, Double_t *par) {
0096   return  ((0.5*par[0]*par[1]/3.14)/
0097     TMath::Max(1.e-10,((x[0]-par[2])*(x[0]-par[2])+.25*par[1]*par[1])))+par[3]+par[4]*x[0];
0098 }
0099
0100 /**
0101  * This function fits TH2F slices with a selected fitType function (gaussian, lorentz, ...).
0102  */
0103 TGraphErrors* fit2DProj(TString name, TString path, int minEntries, int rebinX, int rebinY, int fitType,
0104                         TFile * outputFile, const TString & resonanceType, const double & xDisplace, const TString & append) {
0105
0106   //Read the TH2 from file
0107   TFile *inputFile = new TFile(path);
0108   TH2 * histo = (TH2*) inputFile->Get(name);
0109   if( rebinX > 0 ) histo->RebinX(rebinX);
0110   if( rebinY > 0 ) histo->RebinY(rebinY);
0111
0112   //Declare some variables
0113   TH1 * histoY;
0114   TString nameY;
0115   std::vector<double> Ftop;
0116   std::vector<double> Fwidth;
0117   std::vector<double> Fmass;
0118   std::vector<double> Etop;
0119   std::vector<double> Ewidth;
0120   std::vector<double> Emass;
0121   std::vector<double> Fchi2;
0122   std::vector<double> Xcenter;
0123   std::vector<double> Ex;
0124
0125   TString fileOutName("fitCompare2"+name);
0126   fileOutName += append;
0127   fileOutName += ".root";
0128   TFile *fileOut=new TFile(fileOutName,"RECREATE");
0129
0130   for (int i=1; i<=(histo->GetXaxis()->GetNbins());i++) {
0131
0132     //Project on Y (set name and title)
0133     std::stringstream number;
0134     number << i;
0135     TString numberString(number.str());
0136     nameY = name + "_" + numberString;
0137     // std::cout << "nameY  " << nameY << std::endl;
0138
0139     histoY = histo->ProjectionY(nameY, i, i);
0140
0141     double xBin = histo->GetXaxis()->GetBinCenter(i);
0142     std::stringstream xBinString;
0143     xBinString << xBin;
0144     TString title("Projection of x = ");
0145     title += xBinString.str();
0146
0147     histoY->SetName(title);
0148
0149     if (histoY->GetEntries() > minEntries) {
0150
0151       //Make the dirty work!
0152       TF1 *fit;
0153       std::cout << "fitType = " << fitType << std::endl;
0154       if(fitType == 1) fit = gaussianFit(histoY, resonanceType);
0155       else if(fitType == 2) fit = lorentzianFit(histoY, resonanceType);
0156       else if(fitType == 3) fit = linLorentzianFit(histoY);
0157       else {
0158     std::cout<<"Wrong fit type: 1=gaussian, 2=lorentzian, 3=lorentzian+linear."<<std::endl;
0159     abort();
0160       }
0161
0162       double *par = fit->GetParameters();
0163       double *err = fit->GetParErrors();
0164
0165       // Check the histogram alone
0166       TCanvas *canvas = new TCanvas(nameY+"alone", nameY+" alone");
0167       histoY->Draw();
0168       canvas->Write();
0169
0170       // Only for check
0171       TCanvas *c = new TCanvas(nameY, nameY);
0172
0173       histoY->Draw();
0174       fit->Draw("same");
0175       fileOut->cd();
0176       c->Write();
0177
0178       if( par[0] == par[0] ) {
0179         //Store the fit results
0180         Ftop.push_back(par[0]);
0181         Fwidth.push_back(fabs(par[1]));//sometimes the gaussian has negative width (checked with Rene Brun)
0182         Fmass.push_back(par[2]);
0183         Etop.push_back(err[0]);
0184         Ewidth.push_back(err[1]);
0185         Emass.push_back(err[2]);
0186
0187         Fchi2.push_back(fit->GetChisquare()/fit->GetNDF());
0188
0189         double xx= histo->GetXaxis()->GetBinCenter(i);
0190         Xcenter.push_back(xx);
0191         double ex = 0;
0192         Ex.push_back(ex);
0193       }
0194       else {
0195         // Skip nan
0196         std::cout << "Skipping nan" << std::endl;
0197         Ftop.push_back(0);
0198         Fwidth.push_back(0);
0199         Fmass.push_back(0);
0200         Etop.push_back(1);
0201         Ewidth.push_back(1);
0202         Emass.push_back(1);
0203
0204         Fchi2.push_back(100000);
0205
0206         Xcenter.push_back(0);
0207         Ex.push_back(1);
0208       }
0209     }
0210   }
0211
0212   fileOut->Close();
0213
0214   //Plots the fit results in  TGraphs
0215   const int nn= Ftop.size();
0216   double x[nn],ym[nn],e[nn],eym[nn];
0217   double yw[nn],eyw[nn],yc[nn];
0218
0219   // std::cout << "number of bins = " << nn << std::endl;
0220   // std::cout << "Values:" << std::endl;
0221
0222   for (int j=0;j<nn;j++){
0223     // std::cout << "xCenter["<<j<<"] = " << Xcenter[j] << std::endl;
0224     x[j]=Xcenter[j]+xDisplace;
0225     // std::cout << "Fmass["<<j<<"] = " << Fmass[j] << std::endl;
0226     ym[j]=Fmass[j];
0227     // std::cout << "Emass["<<j<<"] = " << Emass[j] << std::endl;
0228     eym[j]=Emass[j];
0229     // std::cout << "Fwidth["<<j<<"] = " << Fwidth[j] << std::endl;
0230     yw[j]=Fwidth[j];
0231     // std::cout << "Ewidth["<<j<<"] = " << Ewidth[j] << std::endl;
0232     eyw[j]=Ewidth[j];
0233     // std::cout << "Fchi2["<<j<<"] = " << Fchi2[j] << std::endl;
0234     yc[j]=Fchi2[j];
0235     e[j]=0;
0236   }
0237
0238   TGraphErrors *grM = new TGraphErrors(nn,x,ym,e,eym);
0239   grM->SetTitle(name+"_M");
0240   grM->SetName(name+"_M");
0241   TGraphErrors *grW = new TGraphErrors(nn,x,yw,e,eyw);
0242   grW->SetTitle(name+"_W");
0243   grW->SetName(name+"_W");
0244   TGraphErrors *grC = new TGraphErrors(nn,x,yc,e,e);
0245   grC->SetTitle(name+"_chi2");
0246   grC->SetName(name+"_chi2");
0247
0248   grM->SetMarkerColor(4);
0249   grM->SetMarkerStyle(20);
0250   grW->SetMarkerColor(4);
0251   grW->SetMarkerStyle(20);
0252   grC->SetMarkerColor(4);
0253   grC->SetMarkerStyle(20);
0254
0255   //Draw and save the graphs
0256   outputFile->cd();
0257   TCanvas * c1 = new TCanvas(name+"_W",name+"_W");
0258   c1->cd();
0259   grW->Draw("AP");
0260   c1->Write();
0261   TCanvas * c2 = new TCanvas(name+"_M",name+"_M");
0262   c2->cd();
0263   grM->Draw("AP");
0264   c2->Write();
0265   TCanvas * c3 = new TCanvas(name+"_C",name+"_C");
0266   c3->cd();
0267   grC->Draw("AP");
0268   c3->Write();
0269
0270   return grM;
0271 }
0272
0273 TF1* gaussianFit(TH1* histoY, const TString & resonanceType){
0274
0275   TString name = histoY->GetName() + TString("Fit");
0276
0277   // Fit slices projected along Y from bins in X
0278   // -------------------------------------------
0279   TF1 *fit = 0;
0280   // Set parameters according to the selected resonance
0281   if( resonanceType == "JPsi" ) {
0282     fit = new TF1(name,gaussian,2,4,3);
0283     fit->SetParLimits(2, 3.09, 3.15);
0284   }
0285   if( resonanceType.Contains("Upsilon") ) {
0286     fit = new TF1(name,gaussian,8.5,10.5,3);
0287     // fit = new TF1(name,gaussian,9,11,3);
0288     fit->SetParLimits(2, 9.2, 9.6);
0289     fit->SetParLimits(1, 0.09, 0.1);
0290   }
0291   if( resonanceType == "Z" ) {
0292     fit = new TF1(name,gaussian,80,100,3);
0293     fit->SetParLimits(2, 80, 100);
0294     fit->SetParLimits(1, 0.01, 1);
0295   }
0296   if( resonanceType == "reso" ) {
0297     fit = new TF1(name,gaussian,-0.05,0.05,3);
0298     fit->SetParLimits(2, -0.5, 0.5);
0299     fit->SetParLimits(1, 0, 0.5);
0300   }
0301   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean());
0302   // fit->SetParLimits(1, 0.01, 1);
0303   //   fit->SetParLimits(1, 40, 60);
0304   fit->SetParNames("norm","width","mean");
0305   fit->SetLineWidth(2);
0306
0307   if( histoY->GetNbinsX() > 1000 ) histoY->Rebin(10);
0308   histoY->Fit(name,"R0");
0309
0310   return fit;
0311 }
0312
0313 TF1* lorentzianFit(TH1* histoY, const TString & resonanceType){
0314   TString name = histoY->GetName() + TString("Fit");
0315
0316   // Fit slices projected along Y from bins in X
0317   TF1 *fit = 0;
0318   if( resonanceType == "JPsi" || resonanceType == "Psi2S" ) {
0319     fit = new TF1(name, lorentzian, 3.09, 3.15, 3);
0320   }
0321   if( resonanceType.Contains("Upsilon") ) {
0322     fit = new TF1(name, lorentzian, 9, 10, 3);
0323   }
0324   if( resonanceType == "Z" ) {
0325     fit = new TF1(name, lorentzian, 80, 105, 3);
0326   }
0327   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean());
0328   fit->SetParNames("norm","width","mean");
0329   fit->SetLineWidth(2);
0330   histoY->Fit( name,"R0" );
0331   return fit;
0332 }
0333
0334 TF1* linLorentzianFit(TH1* histoY){
0335   TString name = histoY->GetName() + TString("Fit");
0336
0337   // Fit slices projected along Y from bins in X
0338   TF1 *fit = new TF1(name,lorentzianPlusLinear,70,110,5);
0339   fit->SetParameters(histoY->GetMaximum(),histoY->GetRMS(),histoY->GetMean(),10,-0.1);
0340   fit->SetParNames("norm","width","mean","offset","slope");
0341   fit->SetParLimits(1,-10,10);
0342   //fit->SetParLimits(2,85,95);
0343   //fit->SetParLimits(2,90,93);
0344   fit->SetParLimits(2,2,4);
0345   fit->SetParLimits(3,0,100);
0346   fit->SetParLimits(4,-1,0);
0347   fit->SetLineWidth(2);
0348   //histoY -> Fit(name,"0","",85,97);
0349   histoY -> Fit(name,"0","",2,4);
0350   return fit;
0351 }
0352
0353 /****************************************************************************************/
0354 void macroPlot( TString name, TString nameGen, const TString & nameFile1, const TString & nameFile2, const TString & title,
0355                 const TString & resonanceType, const int rebinX, const int rebinY, const int fitType,
0356                 const TString & outputFileName, const bool genMass) {
0357
0358   gROOT->SetBatch(true);
0359
0360   //Save the graphs in a file
0361   TFile *outputFile = new TFile(outputFileName,"RECREATE");
0362
0363   setTDRStyle();
0364
0365   TGraphErrors *grM_1 = fit2DProj(name, nameFile1, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "_1");
0366   TGraphErrors *grM_2 = fit2DProj(name, nameFile2, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "_2");
0367   TGraphErrors *grM_Gen = fit2DProj(nameGen, nameFile2, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "");
0368
0369   TCanvas *c = new TCanvas(name+"_Z",name+"_Z");
0370   c->SetGridx();
0371   c->SetGridy();
0372
0373   grM_1->SetMarkerColor(1);
0374   grM_2->SetMarkerColor(2);
0375   if(genMass)
0376     grM_Gen->SetMarkerColor(4);
0377
0378   TString xAxisTitle;
0379
0380   double x[2],y[2];
0381
0382   if( name.Contains("Eta") ) {
0383     x[0]=-3; x[1]=3;       //<------useful for reso VS eta
0384     xAxisTitle = "muon #eta";
0385   }
0386   else if( name.Contains("PhiPlus") || name.Contains("PhiMinus") ) {
0387     x[0] = -3.2; x[1] = 3.2;
0388     xAxisTitle = "muon #phi(rad)";
0389   }
0390   else {
0391     x[0] = 0.; x[1] = 200;
0392     xAxisTitle = "muon pt(GeV)";
0393   }
0394   if( resonanceType == "JPsi" || resonanceType == "Psi2S" ) { y[0]=0.; y[1]=6.; }
0395   else if( resonanceType.Contains("Upsilon") ) { y[0]=8.; y[1]=12.; }
0396   else if( resonanceType == "Z" ) { y[0]=80; y[1]=100; }
0397
0398   // This is used to have a canvas containing both histogram points
0399   TGraph *gr = new TGraph(2,x,y);
0400   gr->SetMarkerStyle(8);
0401   gr->SetMarkerColor(108);
0402   gr->GetYaxis()->SetTitle("Res Mass(GeV)   ");
0403   gr->GetYaxis()->SetTitleOffset(1);
0404   gr->GetXaxis()->SetTitle(xAxisTitle);
0405   gr->SetTitle(title);
0406
0407   // Text for the fits
0408   TPaveText * paveText1 = new TPaveText(0.20,0.15,0.49,0.28,"NDC");
0409   paveText1->SetFillColor(0);
0410   paveText1->SetTextColor(1);
0411   paveText1->SetTextSize(0.02);
0412   paveText1->SetBorderSize(1);
0413   TPaveText * paveText2 = new TPaveText(0.59,0.15,0.88,0.28,"NDC");
0414   paveText2->SetFillColor(0);
0415   paveText2->SetTextColor(2);
0416   paveText2->SetTextSize(0.02);
0417   paveText2->SetBorderSize(1);
0418
0419   TPaveText * paveText3 = new TPaveText(0.59,0.32,0.88,0.45,"NDC");
0420   paveText3->SetFillColor(0);
0421   paveText3->SetTextColor(4);
0422   paveText3->SetTextSize(0.02);
0423   paveText3->SetBorderSize(1);
0424
0425   /*****************PARABOLIC FIT (ETA)********************/
0426   if( name.Contains("Eta") ) {
0427     std::cout << "Fitting eta" << std::endl;
0428     TF1 *fit1 = new TF1("fit1",onlyParabolic,-3.2,3.2,2);
0429     fit1->SetLineWidth(2);
0430     fit1->SetLineColor(1);
0431     if( grM_1->GetN() > 0 ) grM_1->Fit("fit1","", "", -3, 3);
0432     setTPaveText(fit1, paveText1);
0433
0434     TF1 *fit2 = new TF1("fit2","pol0",-3.2,3.2);
0435     // TF1 *fit2 = new TF1("fit2",onlyParabolic,-3.2,3.2,2);
0436     fit2->SetLineWidth(2);
0437     fit2->SetLineColor(2);
0438     if( grM_2->GetN() > 0 ) grM_2->Fit("fit2","", "", -3, 3);
0439     // grM_2->Fit("fit2","R");
0440     setTPaveText(fit2, paveText2);
0441     if(genMass){
0442       TF1 *fit3 = new TF1("fit3",onlyParabolic,-3.2,3.2,2);
0443       fit3->SetLineWidth(2);
0444       fit3->SetLineColor(4);
0445       if( grM_Gen->GetN() > 0 ) grM_Gen->Fit("fit3","", "", -3, 3);
0446       // grM_2->Fit("fit2","R");
0447       setTPaveText(fit3, paveText3);
0448     }
0449   }
0450   /*****************SINUSOIDAL FIT (PHI)********************/
0451 //   if( name.Contains("Phi") ) {
0452 //     std::cout << "Fitting phi" << std::endl;
0453 //     TF1 *fit1 = new TF1("fit1",sinusoidal,-3.2,3.2,3);
0454 //     fit1->SetParameters(9.45,1,1);
0455 //     fit1->SetParNames("offset","amplitude","phase");
0456 //     fit1->SetLineWidth(2);
0457 //     fit1->SetLineColor(1);
0458 //     fit1->SetParLimits(2,-3.14,3.14);
0459 //     if( grM_1->GetN() > 0 ) grM_1->Fit("fit1","","",-3,3);
0460 //     setTPaveText(fit1, paveText1);
0461
0462 //     TF1 *fit2 = new TF1("fit2",sinusoidal,-3.2,3.2,3);
0463 //     fit2->SetParameters(9.45,1,1);
0464 //     fit2->SetParNames("offset","amplitude","phase");
0465 //     fit2->SetLineWidth(2);
0466 //     fit2->SetLineColor(2);
0467 //     fit2->SetParLimits(2,-3.14,3.14);
0468 //     if( grM_2->GetN() > 0 ) grM_2->Fit("fit2","","",-3,3);
0469 //     setTPaveText(fit2, paveText2);
0470
0471 //     if(genMass){
0472 //       TF1 *fit3 = new TF1("fit3",sinusoidal,-3.2,3.2,3);
0473 //       fit3->SetParameters(9.45,1,1);
0474 //       fit3->SetParNames("offset","amplitude","phase");
0475 //       fit3->SetLineWidth(2);
0476 //       fit3->SetLineColor(4);
0477 //       fit3->SetParLimits(2,-3.14,3.14);
0478 //       if( grM_Gen->GetN() > 0 ) grM_Gen->Fit("fit3","","",-3,3);
0479 //       setTPaveText(fit3, paveText3);
0480 //     }
0481 //   }
0482   /*****************************************/
0483   if( name.Contains("Pt") ) {
0484     std::cout << "Fitting pt" << std::endl;
0485     // TF1 * fit1 = new TF1("fit1", linear, 0., 150., 2);
0486     TF1 * fit1 = new TF1("fit1", "[0]", 0., 150.);
0487     fit1->SetParameters(0., 1.);
0488     fit1->SetParNames("scale","pt coefficient");
0489     fit1->SetLineWidth(2);
0490     fit1->SetLineColor(1);
0491     if( grM_1->GetN() > 0 ) {
0492       if( name.Contains("Z") ) grM_1->Fit("fit1","","",0.,150.);
0493       else grM_1->Fit("fit1","","",0.,27.);
0494     }
0495     setTPaveText(fit1, paveText1);
0496
0497     // TF1 * fit2 = new TF1("fit2", linear, 0., 150., 2);
0498     TF1 * fit2 = new TF1("fit2", "[0]", 0., 150.);
0499     fit2->SetParameters(0., 1.);
0500     fit2->SetParNames("scale","pt coefficient");
0501     fit2->SetLineWidth(2);
0502     fit2->SetLineColor(2);
0503     if( grM_2->GetN() > 0 ) {
0504       if( name.Contains("Z") ) grM_2->Fit("fit2","","",0.,150.);
0505       grM_2->Fit("fit2","","",0.,27.);
0506     }
0507     setTPaveText(fit2, paveText2);
0508
0509     if(genMass){
0510       TF1 * fit3 = new TF1("fit3", "[0]", 0., 150.);
0511       fit3->SetParameters(0., 1.);
0512       fit3->SetParNames("scale","pt coefficient");
0513       fit3->SetLineWidth(2);
0514       fit3->SetLineColor(4);
0515       if( grM_Gen->GetN() > 0 ) {
0516     if( name.Contains("Z") ) grM_Gen->Fit("fit3","","",0.,150.);
0517     grM_Gen->Fit("fit3","","",0.,27.);
0518       }
0519       setTPaveText(fit3, paveText3);
0520     }
0521   }
0522
0523   c->cd();
0524   gr->Draw("AP");
0525   grM_1->Draw("P");
0526   grM_2->Draw("P");
0527   if(genMass)
0528     grM_Gen->Draw("P");
0529
0530   paveText1->Draw("same");
0531   paveText2->Draw("same");
0532   if(genMass)
0533     paveText3->Draw("same");
0534
0535   TLegend *leg = new TLegend(0.65,0.85,1,1);
0536   leg->SetFillColor(0);
0537   leg->AddEntry(grM_1,"before calibration","P");
0538   leg->AddEntry(grM_2,"after calibration","P");
0539   if(genMass)
0540     leg->AddEntry(grM_Gen, "generated mass", "P");
0541   leg->Draw("same");
0542
0543   outputFile->cd();
0544   c->Write();
0545   outputFile->Close();
0546 }
0547
0548 Double_t sinusoidal(Double_t *x, Double_t *par) {
0549   return (par[0] + par[1]*sin(x[0]+par[2]));
0550 }
0551
0552 Double_t parabolic(Double_t *x, Double_t *par) {
0553   return (par[0] + par[1]*fabs(x[0]) + par[2]*x[0]*x[0]) ;
0554 }
0555
0556 Double_t overX(Double_t *x, Double_t *par) {
0557   return (par[0] + par[1]/x[0]) ;
0558 }
0559
0560 Double_t parabolic2(Double_t *x, Double_t *par) {
0561   if(x>0)
0562     return (par[0] + par[1]*fabs(x[0]) + par[2]*x[0]*x[0]) ;
0563   else
0564     return (par[0] - par[1]*fabs(x[0]) + par[2]*x[0]*x[0]) ;
0565 }
0566
0567 Double_t onlyParabolic(Double_t *x, Double_t *par) {
0568   return (par[0] + par[1]*x[0]*x[0]) ;
0569 }
0570
0571 Double_t linear(Double_t *x, Double_t *par) {
0572   return (par[0] + par[1]*fabs(x[0])) ;
0573 }
0574
0575
0576 /****************************************************************************************/
0577
0578 void setTDRStyle() {
0579   TStyle *tdrStyle = new TStyle("tdrStyle","Style for P-TDR");
0580
0581   // For the canvas:
0582   tdrStyle->SetCanvasBorderMode(0);
0583   tdrStyle->SetCanvasColor(kWhite);
0584   tdrStyle->SetCanvasDefH(600); //Height of canvas
0585   tdrStyle->SetCanvasDefW(600); //Width of canvas
0586   tdrStyle->SetCanvasDefX(0);   //POsition on screen
0587   tdrStyle->SetCanvasDefY(0);
0588
0589   // For the Pad:
0590   tdrStyle->SetPadBorderMode(0);
0591   tdrStyle->SetPadColor(kWhite);
0592   tdrStyle->SetPadGridX(false);
0593   tdrStyle->SetPadGridY(false);
0594   tdrStyle->SetGridColor(0);
0595   tdrStyle->SetGridStyle(3);
0596   tdrStyle->SetGridWidth(1);
0597
0598   // For the frame:
0599   tdrStyle->SetFrameBorderMode(0);
0600   tdrStyle->SetFrameBorderSize(1);
0601   tdrStyle->SetFrameFillColor(0);
0602   tdrStyle->SetFrameFillStyle(0);
0603   tdrStyle->SetFrameLineColor(1);
0604   tdrStyle->SetFrameLineStyle(1);
0605   tdrStyle->SetFrameLineWidth(1);
0606
0607   // For the histo:
0608   tdrStyle->SetHistLineColor(1);
0609   tdrStyle->SetHistLineStyle(0);
0610   tdrStyle->SetHistLineWidth(1);
0611
0612   tdrStyle->SetEndErrorSize(2);
0613   tdrStyle->SetErrorX(0.);
0614
0615   tdrStyle->SetMarkerStyle(20);
0616
0617   //For the fit/function:
0618   tdrStyle->SetOptFit(1);
0619   tdrStyle->SetFitFormat("5.4g");
0620   tdrStyle->SetFuncColor(2);
0621   tdrStyle->SetFuncStyle(1);
0622   tdrStyle->SetFuncWidth(1);
0623
0624   //For the date:
0625   tdrStyle->SetOptDate(0);
0626
0627   // For the statistics box:
0628   tdrStyle->SetOptFile(0);
0629   tdrStyle->SetOptStat(0); // To display the mean and RMS:   SetOptStat("mr");
0630   tdrStyle->SetStatColor(kWhite);
0631   tdrStyle->SetStatFont(42);
0632   tdrStyle->SetStatFontSize(0.025);
0633   tdrStyle->SetStatTextColor(1);
0634   tdrStyle->SetStatFormat("6.4g");
0635   tdrStyle->SetStatBorderSize(1);
0636   tdrStyle->SetStatH(0.1);
0637   tdrStyle->SetStatW(0.15);
0638
0639   // Margins:
0640   tdrStyle->SetPadTopMargin(0.05);
0641   tdrStyle->SetPadBottomMargin(0.13);
0642   tdrStyle->SetPadLeftMargin(0.13);
0643   tdrStyle->SetPadRightMargin(0.05);
0644
0645   // For the Global title:
0646   tdrStyle->SetOptTitle(0);
0647
0648   // For the axis titles:
0649   tdrStyle->SetTitleColor(1, "XYZ");
0650   tdrStyle->SetTitleFont(42, "XYZ");
0651   tdrStyle->SetTitleSize(0.06, "XYZ");
0652   tdrStyle->SetTitleXOffset(0.9);
0653   tdrStyle->SetTitleYOffset(1.05);
0654
0655   // For the axis labels:
0656   tdrStyle->SetLabelColor(1, "XYZ");
0657   tdrStyle->SetLabelFont(42, "XYZ");
0658   tdrStyle->SetLabelOffset(0.007, "XYZ");
0659   tdrStyle->SetLabelSize(0.05, "XYZ");
0660
0661   // For the axis:
0662   tdrStyle->SetAxisColor(1, "XYZ");
0663   tdrStyle->SetStripDecimals(kTRUE);
0664   tdrStyle->SetTickLength(0.03, "XYZ");
0665   tdrStyle->SetNdivisions(510, "XYZ");
0666   tdrStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
0667   tdrStyle->SetPadTickY(1);
0668
0669   // Change for log plots:
0670   tdrStyle->SetOptLogx(0);
0671   tdrStyle->SetOptLogy(0);
0672   tdrStyle->SetOptLogz(0);
0673
0674   //tdrStyle->SetOptFit(00000);
0675   tdrStyle->cd();
0676 }