test/Macros/ResolFit.C

0001 #include <iostream>
0002 #include <string>
0003 #include <sstream>
0004 #include <fstream>
0005 #include <list>
0006
0007 #include "TROOT.h"
0008 #include "TH1D.h"
0009 #include "TF1.h"
0010 #include "TCanvas.h"
0011 #include "TPaveText.h"
0012 #include "TFile.h"
0013
0014 /**
0015  * This is the pt vs eta resolution by points. It uses fabs(eta) assuming symmetry.
0016  */
0017 Double_t etaByPoints(Double_t * inEta, Double_t * par) {
0018   Double_t sigmaPtVsEta = 0.;
0019   Double_t eta = fabs(inEta[0]);
0020   if( 0. <= eta && eta <= 0.2 ) sigmaPtVsEta = 0.0120913;
0021   else if( 0.2 < eta && eta <= 0.4 ) sigmaPtVsEta = 0.0122204;
0022   else if( 0.4 < eta && eta <= 0.6 ) sigmaPtVsEta = 0.0136937;
0023   else if( 0.6 < eta && eta <= 0.8 ) sigmaPtVsEta = 0.0142069;
0024   else if( 0.8 < eta && eta <= 1.0 ) sigmaPtVsEta = 0.0177526;
0025   else if( 1.0 < eta && eta <= 1.2 ) sigmaPtVsEta = 0.0243587;
0026   else if( 1.2 < eta && eta <= 1.4 ) sigmaPtVsEta = 0.019994;
0027   else if( 1.4 < eta && eta <= 1.6 ) sigmaPtVsEta = 0.0185132;
0028   else if( 1.6 < eta && eta <= 1.8 ) sigmaPtVsEta = 0.0177141;
0029   else if( 1.8 < eta && eta <= 2.0 ) sigmaPtVsEta = 0.0211577;
0030   else if( 2.0 < eta && eta <= 2.2 ) sigmaPtVsEta = 0.0255051;
0031   else if( 2.2 < eta && eta <= 2.4 ) sigmaPtVsEta = 0.0338104;
0032   // ATTENTION: This point has a big error and it is very displaced from the rest of the distribution.
0033   else if( 2.4 < eta && eta <= 2.6 ) sigmaPtVsEta = 0.31;
0034   return ( par[0]*sigmaPtVsEta );
0035 }
0036
0037 void myFunc()
0038 {
0039   TF1 *f1 = new TF1("myFunc",etaByPoints,-2.5,0,1);
0040   f1->SetParameter(0, 1);
0041   f1->SetParNames("constant");
0042   // f1->Draw();
0043 }
0044
0045 /**
0046  * This function draws a histogram and a function on a canvas and adds a text box with the results of the fit.
0047  */
0048 void draw(TH1D * h, TF1 * f) {
0049   Width_t lineWidth = Width_t(0.4);
0050   Color_t lineColor = kRed;
0051   f->SetLineColor(lineColor);
0052   f->SetLineWidth(lineWidth);
0053   TCanvas c(h->GetName(),h->GetTitle(),1000,800);
0054   c.cd();
0055   h->Draw();
0056   f->Draw("same");
0057
0058   // Text box with fit results
0059   TPaveText * fitLabel = new TPaveText(0.78,0.71,0.98,0.91,"NDC");
0060   fitLabel->SetBorderSize(1);
0061   fitLabel->SetTextAlign(12);
0062   fitLabel->SetTextSize(0.02);
0063   fitLabel->SetFillColor(0);
0064   fitLabel->AddText("Function: "+f->GetExpFormula());
0065   for( int i=0; i<f->GetNpar(); ++i ) {
0066     char name[50];
0067     std::cout << "par["<<i<<"] = " << f->GetParameter(i) << std::endl;
0068     sprintf(name, "par[%i] = %4.2g #pm %4.2g",i, f->GetParameter(i), f->GetParError(i));
0069     fitLabel->AddText(name);
0070   }
0071   fitLabel->Draw("same");
0072
0073   c.Write();
0074   h->Write();
0075 }
0076
0077 /**
0078  * This function reads parameters from the FitParameters.txt file and returns them in a vector.
0079  */
0080 pair<list<double>, list<double> > readParameters(int fitFile) {
0081   list<double> parameters;
0082   list<double> parameterErrors;
0083
0084   std::ifstream a("FitParameters.txt");
0085   std::string line;
0086   bool indexFound = false;
0087   std::string iteration("Iteration ");
0088   while (a) {
0089     getline(a,line);
0090     unsigned int lineInt = line.find("value");
0091
0092     // Take only the values from the matching iteration
0093     if( line.find(iteration) != std::string::npos ) {
0094       std::stringstream iterationNum;
0095       iterationNum << fitFile;
0096       if( line.find(iteration+iterationNum.str()) != std::string::npos ) {
0097         indexFound = true;
0098         std::cout << "In: " << line << std::endl;
0099         std::cout << "Found Index = " << iteration+iterationNum.str() << std::endl;
0100       }
0101       else indexFound = false;
0102     }
0103
0104     if ( (lineInt != std::string::npos) && indexFound ) {
0105       int subStr1 = line.find("value");
0106       int subStr2 = line.find("+-");
0107       std::stringstream paramStr;
0108       double param = 0;
0109       paramStr << line.substr(subStr1+5,subStr2);
0110       paramStr >> param;
0111       parameters.push_back(param);
0112       // std::cout << "paramStr = " << line.substr(subStr1+5,subStr2) << std::endl;
0113       std::stringstream parErrorStr;
0114       double parError = 0;
0115       parErrorStr << line.substr(subStr2+1);
0116       parErrorStr >> parError;
0117       parameterErrors.push_back(parError);
0118
0119       // std::cout << "param = " << param << std::endl;
0120       // std::cout << "parError = " << parError << std::endl;
0121     }
0122   }
0123
0124   //     std::cout << "Reading function from file" << std::endl;
0125   //     TString param = "aaa a a a  a value = 193.4+-12";
0126   //     int id = param.Index("value");
0127   //     int length = param.Length();
0128   //     std::cout << "param(id,-1)" << param(id, length) << std::endl;
0129   return make_pair(parameters, parameterErrors);
0130 }
0131
0132 /**
0133  * Set parameters from the list to the function. It empties the list while using it.
0134  */
0135 void setParameters(TF1 * f, pair<list<double>, list<double> > & parameters) {
0136   int parNum = f->GetNpar();
0137   for( int iPar=0; iPar<parNum; ++iPar ) {
0138     // Read the first element and remove it from the list
0139     f->SetParameter(iPar, parameters.first.front());
0140     f->SetParError(iPar, parameters.second.front());
0141     parameters.first.pop_front();
0142     parameters.second.pop_front();
0143   }
0144 }
0145
0146 /**
0147  * This macro fits and draws the histograms in the file written by ResolDraw.cc.
0148  * If the false is passed as argument, it searches for the file FitParameters, reads
0149  * the parameters of the functions to draw over the histograms instead of fitting.
0150  * ATTENTION: the functions used in this case must be the same of those which parameters
0151  * have been written in FitParameters.txt.
0152  */
0153 int ResolFit( int fitFile = -1 ) {
0154
0155   TString mainPtName("PtResolution");
0156   TString mainCotgThetaName("CotgThetaResolution");
0157   TString mainPhiName("PhiResolution");
0158
0159   // Read the values from the FitParameters.txt file if required
0160   pair<list<double>, list<double> > parameters;
0161   if( fitFile != -1 ) {
0162     parameters = readParameters( fitFile );
0163   }
0164
0165   TFile inputFile("redrawed.root","READ");
0166   TFile outputFile("fitted.root","RECREATE");
0167
0168   outputFile.cd();
0169
0170   // Pt resolution
0171   // -------------
0172   // VS pt
0173   TDirectory * tempDir = (TDirectory*) inputFile.Get(mainPtName+"GenVSMu");
0174   TH1D * h = (TH1D*) tempDir->Get(mainPtName+"GenVSMu_ResoVSPt_resol");
0175   TF1 *f = new TF1("f","pol1",0,100);
0176
0177   if( fitFile == -1 ) {
0178     std::cout << "Fitting Pt resolution vs Pt" << std::endl;
0179     h->Fit("f","R0");
0180   }
0181   else {
0182     setParameters(f, parameters);
0183     // Put back the constant so that it can be used also by the following function
0184     parameters.first.push_front(f->GetParameter(0));
0185     parameters.second.push_front(f->GetParError(0));
0186   }
0187
0188   h->SetMinimum(0);
0189   h->SetMaximum(0.045);
0190   h->GetXaxis()->SetTitle("pt(GeV)");
0191   h->GetYaxis()->SetTitleOffset(1.4);
0192   h->GetYaxis()->SetTitle("#sigma pt");
0193   draw(h,f);
0194
0195   std::cout << "sigmapt vs eta" << std::endl;
0196   // VS eta
0197   tempDir = (TDirectory*) inputFile.Get(mainPtName+"GenVSMu");
0198   h = (TH1D*) tempDir->Get(mainPtName+"GenVSMu_ResoVSEta_resol");
0199
0200   // f = new TF1("f","pol2",-2.5,2.5);
0201
0202   // This call is needed in order to use myFunc in the fit.
0203   myFunc();
0204   f = (TF1*)gROOT->GetFunction("myFunc");
0205   f->SetParameter(0,1.);
0206
0207   if( fitFile == -1 ) {
0208     std::cout << "Fitting Pt resolution vs Eta" << std::endl;
0209     // h->Fit("f","R0");
0210     h->Fit("myFunc","R0");
0211   }
0212   else {
0213     setParameters(f, parameters);
0214   }
0215
0216   h->SetMinimum(0);
0217   h->SetMaximum(0.045);
0218   h->GetXaxis()->SetTitle("#eta");
0219   h->GetYaxis()->SetTitleOffset(1.4);
0220   h->GetYaxis()->SetTitle("#sigma pt");
0221   draw(h,f);
0222
0223   // CotgTheta resolution
0224   // --------------------
0225   // VS pt
0226   tempDir = (TDirectory*) inputFile.Get(mainCotgThetaName+"GenVSMu");
0227   h = (TH1D*) tempDir->Get(mainCotgThetaName+"GenVSMu_ResoVSPt_resol");
0228
0229   f = new TF1("f","[0]+[1]/x",0,100);
0230   if( fitFile == -1 ) {
0231     h->Fit("f","R0");
0232     std::cout << "Fitting CotgTheta resolution vs Pt" << std::endl;
0233   }
0234   else {
0235     setParameters(f, parameters);
0236     parameters.first.push_front(f->GetParameter(0));
0237     parameters.second.push_front(f->GetParError(0));
0238   }
0239
0240   h->SetMinimum(0);
0241   h->SetMaximum(0.0014);
0242   h->GetXaxis()->SetTitle("pt(GeV)");
0243   h->GetYaxis()->SetTitleOffset(1.4);
0244   h->GetYaxis()->SetTitle("#sigma cotg#theta");
0245   draw(h,f);
0246   // VS eta
0247   tempDir = (TDirectory*) inputFile.Get(mainCotgThetaName+"GenVSMu");
0248   h = (TH1D*) tempDir->Get(mainCotgThetaName+"GenVSMu_ResoVSEta_resol");
0249
0250   f = new TF1("f","pol2",-2.5,2.5);
0251   if( fitFile == -1 ) {
0252     std::cout << "Fitting CotgTheta resolution vs Eta" << std::endl;
0253     h->Fit("f","R0");
0254   }
0255   else {
0256     setParameters(f, parameters);
0257   }
0258
0259   h->SetMinimum(0);
0260   h->SetMaximum(0.0035);
0261   h->GetXaxis()->SetTitle("eta");
0262   h->GetYaxis()->SetTitleOffset(1.4);
0263   h->GetYaxis()->SetTitle("#sigma cotg#theta");
0264   draw(h,f);
0265
0266   // Phi resolution
0267   // --------------
0268   // VS pt
0269   tempDir = (TDirectory*) inputFile.Get(mainPhiName+"GenVSMu");
0270   h = (TH1D*) tempDir->Get(mainPhiName+"GenVSMu_ResoVSPt_resol");
0271
0272   f = new TF1("f","[0]+[1]/x",0,100);
0273   if( fitFile == -1 ) {
0274     std::cout << "Fitting Phi resolution vs Pt" << std::endl;
0275     h->Fit("f","R0");
0276   }
0277   else {
0278     setParameters(f, parameters);
0279     parameters.first.push_front(f->GetParameter(0));
0280     parameters.second.push_front(f->GetParError(0));
0281   }
0282
0283   h->SetMinimum(0);
0284   h->SetMaximum(0.001);
0285   h->GetXaxis()->SetTitle("pt(GeV)");
0286   h->GetYaxis()->SetTitleOffset(1.4);
0287   h->GetYaxis()->SetTitle("#sigma #phi");
0288   draw(h,f);
0289   // VS eta
0290   tempDir = (TDirectory*) inputFile.Get(mainPhiName+"GenVSMu");
0291   h = (TH1D*) tempDir->Get(mainPhiName+"GenVSMu_ResoVSEta_resol");
0292
0293   f = new TF1("f","pol2",-2.4,2.4);
0294   if( fitFile == -1 ) {
0295     std::cout << "Fitting Phi resolution vs Eta" << std::endl;
0296     h->Fit("f","R0");
0297   }
0298   else {
0299     setParameters(f, parameters);
0300   }
0301
0302   h->SetMinimum(0);
0303   h->SetMaximum(0.005);
0304   h->GetXaxis()->SetTitle("eta");
0305   h->GetYaxis()->SetTitleOffset(1.4);
0306   h->GetYaxis()->SetTitle("#sigma #phi");
0307   draw(h,f);
0308   return 0;
0309 }