test/MuonSeedPTAnalysis/SeedParaFit.C

0001 #include "SeedParaFit.h"
0002
0003 SeedParaFit::SeedParaFit() {
0004   /*
0005    * Get the profile histogram for Pt*dphi vs eta for segment pair case
0006    * Slice each eta bin and do the gaussian fit
0007    * plot the mean and the sigma against the middle of the Y
0008    *
0009    * Author:  Shih-Chuan Kao  --- UCR
0010    *
0011    */
0012
0013   // Name of root file
0014   fName = "mix31x_v3.root";
0015   // folder for all plots
0016   dname = "Pt31x_v3";
0017   // plot formate
0018   suffixps = ".gif";
0019   // fitting range
0020   nsigmas = 1.5;
0021   // debug for each slice and fitting => fitting for each point
0022   debug = false;
0023   // show the pre-fitting of PtxdPhi and Rel_Error
0024   debugAll = false ;
0025   // bin size for eta, unit bin size is 0.01
0026   // bin range => CSC: 0.9 - 2.7 ,  DT: 0.0 - 1.2 , OL: 0.7 - 1.3
0027   // # of bin        : 180            : 120           : 60
0028   // possible rbin value = 1, 2, 3, 4, 5, 10 ;
0029   rbin = 3 ;
0030   bsz = 0.01*rbin ;
0031
0032   ptFunc = new SeedPtFunction();
0033
0034   gSystem->mkdir(dname);
0035   dfile   = fopen(dname+"/ptSeedParameterization_cfi.py","a");
0036   if (debug) logfile = fopen(dname+"/parafitting.log","a");
0037
0038
0039 }
0040
0041 SeedParaFit::~SeedParaFit() {
0042
0043   if ( debug ) delete cv;
0044   delete c1;
0045   delete c2;
0046   delete c3;
0047   //delete c4;
0048
0049   delete ptFunc ;
0050
0051   fclose(dfile);
0052   if (debug) fclose(logfile);
0053
0054 }
0055
0056 void SeedParaFit::ParaFit(int type, int s1, double h1, double h2, int np1){
0057 //void SeedParaFit::ParaFit(int type, int s1, int b1, int b2, int np1){
0058
0059   /*
0060    *  type  : seed type, 1: CSC(CSC), 2:Overlap(OL), 3:DT(DT), 4:DT_Single(SMB) 5:CSC_Sigle(SME)
0061    *  s1    : station combination ex: st1 and st2 => 12
0062    *  h1,h2 : eta range
0063    *  np1   : number of fitting parameters p0 + p1*x + p2*x^2 + .... => np1 = 3
0064    *
0065    */
0066
0067   // name the file title by detector type
0068   if( type == 1 ) {
0069     sprintf(det, "CSC_%d", s1);
0070     det_id = det;
0071     dphi_type = det_id+"_eta_rdphiPt";
0072     if (s1 == 11) {
0073        sprintf(det,"CSC_0%d",1);
0074        det_id = det;
0075        dphi_type = det_id+"_eta_rdphiPt";
0076     }
0077     if (s1 == 12 && h1 > 1.7) {
0078        sprintf(det,"CSC_0%d",2);
0079        det_id = det;
0080        dphi_type = "CSC_12_eta_rdphiPt";
0081     }
0082     if (s1 == 13 && h1 > 1.75) {
0083        sprintf(det,"CSC_0%d",3);
0084        det_id = det;
0085        dphi_type = "CSC_13_eta_rdphiPt";
0086     }
0087   }
0088   if( type == 2 ) {
0089     sprintf(det, "OL_%d", s1);
0090     det_id = det;
0091     dphi_type = det_id+"_eta_rdphiPtA";
0092   }
0093   if( type == 3 ) {
0094     sprintf(det, "DT_%d", s1);
0095     det_id = det;
0096     dphi_type = det_id+"_eta_rdphiPt";
0097   }
0098   if( type == 4 ) {
0099     sprintf(det, "SMB_%d", s1);
0100     det_id = det;
0101     dphi_type = det_id+"_eta_rdphiPt";
0102   }
0103   if( type == 5 ) {
0104     sprintf(det, "SME_%d", s1);
0105     det_id = det;
0106     dphi_type = det_id+"_eta_rdphiPt";
0107   }
0108
0109   // Name the plots
0110   plot01 = dname+"/"+det_id+"_eta_pTxdphi"+suffixps;
0111   plot02 = dname+"/"+det_id+"_eta_RelError"+suffixps;
0112   plot03 = dname+"/"+det_id+"_eta_pTxdphi_scalar"+suffixps;
0113   plot04 = dname+"/"+det_id+"_prefit_test"+suffixps;
0114
0115 // ********************************************************************
0116 // Pointers to histograms
0117 // ********************************************************************
0118
0119   file = TFile::Open(fName);
0120   if (type ==1 )  heta_dphiPt  = (TH2F *) file->Get("CSC_All/"+dphi_type);
0121   if (type ==2 )  heta_dphiPt  = (TH2F *) file->Get("OL_All/"+dphi_type);
0122   if (type ==3 )  heta_dphiPt  = (TH2F *) file->Get("DT_All/"+dphi_type);
0123   if (type ==4 )  heta_dphiPt  = (TH2F *) file->Get("MB_All/"+dphi_type);
0124   if (type ==5 )  heta_dphiPt  = (TH2F *) file->Get("ME_All/"+dphi_type);
0125   heta_dphiPt->RebinX(rbin);
0126
0127
0128 // ********************************************************************
0129 // create a folder to store all histograms
0130 // ********************************************************************
0131
0132   //gSystem->cd(dname);
0133
0134   if ( debug ) gSystem->mkdir("BugInfo");
0135
0136  // *****************************************
0137  // ***** 1 hdeta vs. Pt*dphi   Low eta *****
0138  // *****************************************
0139
0140  Double_t r  = 0.;
0141  Double_t ini =0.;
0142  Double_t fnl =2.5;
0143  Double_t f1 = 0.;
0144  Double_t f2 = 0.;
0145  int b1 = 0;
0146  int b2 = 0;
0147
0148  if (type == 1 || type == 5) {
0149     b1 = (h1 - 0.9)/bsz ;
0150     b2 = (h2 - 0.9)/bsz ;
0151     ini = 0.9 - bsz/2.;
0152     fnl = 2.5 + bsz/2. ;
0153     r = 0.9   + (b1-1.0)*bsz + (bsz/2.);
0154     f1= 0.9   + (b1-1.0)*bsz;
0155     f2= 0.9   + (b2*bsz);
0156  }
0157  if (type ==3 || type ==4) {
0158     b1 = (h1 - 0.)/bsz ;
0159     b2 = (h2 - 0.)/bsz ;
0160     if (b1 == 0) b1 = 1 ;
0161     ini = 0.  - bsz/2. ;
0162     fnl = 1.2 + bsz/2;
0163     r  = 0.0   + (b1-1.0)*bsz + (bsz/2.);
0164     f1 = 0.0   + (b1-1.0)*bsz;
0165     f2 = 0.0   + (b2*bsz);
0166  }
0167  if (type ==2 ) {
0168     b1 = (h1 - 0.7)/bsz ;
0169     b2 = (h2 - 0.7)/bsz ;
0170     ini = 0.7 - bsz/2.;
0171     fnl = 1.3 + bsz/2.;
0172     r  = 0.7   + (b1-1.0)*bsz + (bsz/2.);
0173     f1 = 0.7   + (b1-1.0)*bsz;
0174     f2 = 0.7   + (b2*bsz);
0175     cout<<" b1:"<<b1<<"   b2:"<<b2 <<endl;;
0176     cout<<" f1:"<<f1<<"  f2:"<<f2 <<endl;
0177  }
0178
0179  vector<double> yv1;
0180  vector<double> xv1;
0181  vector<double> yv1e;
0182
0183  yv1.clear();
0184  xv1.clear();
0185  yv1e.clear();
0186
0187  // ***************************************
0188  //   Main Loop for full eta range
0189  // ***************************************
0190  if (debug) fprintf (logfile,"            ");
0191  if (debug) fprintf (logfile,"DetId     p0        p1        p2        chi2        rms       nu \n");
0192  double par0 = 0.;
0193  double par1 = 0.;
0194  double par2 = 0.;
0195  double chi2 = 0.;
0196  double  ndf = 0.;
0197
0198  TF1 *gausf = new TF1("gausf", SeedPtFunction::fgaus ,-2., 2., 3  );
0199
0200  if (debug) {
0201     cv = new TCanvas("cv");
0202     gStyle->SetOptStat(kTRUE);
0203     gStyle->SetOptFit(0111);
0204  }
0205
0206  for (int i=b1; i<b2; i++) {
0207
0208      if ( debug ) {
0209         sprintf(dbplot, "dbug_%d.gif",i);
0210         plot_id = dbplot;
0211      }
0212
0213      TH1D* heta_dphiPt_pjy = heta_dphiPt->ProjectionY("heta_dphiPt_pjy",i,i);
0214      if ( debug )  heta_dphiPt_pjy->Draw();
0215      double  nu = heta_dphiPt_pjy->Integral();
0216      double ave = heta_dphiPt_pjy->GetMean();
0217      double rms = heta_dphiPt_pjy->GetRMS();
0218
0219      double L1 = ave - (3. * rms);
0220      double H1 = ave + (3. * rms);
0221
0222      //cout<<" fk1 nu:"<< nu <<" ave:"<< ave <<"  rms:"<< rms <<endl;
0223      if ( nu >= 10  ) {
0224         //gausf->SetParameter(0, amp );
0225         gausf->SetParameter(1, ave );
0226         gausf->SetParameter(2, rms );
0227         gausf->SetParLimits(1, L1, H1 );
0228         gausf->SetParLimits(2, 0., 5. );
0229         if (!debug) heta_dphiPt_pjy->Fit("gausf", "N0RQL", "", L1, H1 );
0230         if ( debug) {
0231                     gausf->SetLineStyle(2);
0232                     gausf->SetLineColor(2);
0233                     heta_dphiPt_pjy->Fit("gausf", "RQL", "sames", L1, H1 );
0234                     cv->Update();
0235                     cv->Print("BugInfo/"+plot_id);
0236         }
0237     par0 = gausf->GetParameter(0);  // hieght
0238         if ( par0 < 1. || par0 > nu*0.8 ) heta_dphiPt_pjy->Fit("gausf","N0Q" );
0239     par1 = gausf->GetParameter(1);  // mean value
0240     par2 = gausf->GetParameter(2);  // sigma
0241
0242         L1 = par1 - (nsigmas * par2);
0243         H1 = par1 + (nsigmas * par2);
0244
0245         //gausf->SetParameter(0, par0 );
0246         gausf->SetParameter(1, par1 );
0247         gausf->SetParameter(2, par2 );
0248     heta_dphiPt_pjy->Fit("gausf","N0RQ","", L1, H1);
0249         if ( debug) {
0250                     gausf->SetLineStyle(1);
0251                     gausf->SetLineColor(4);
0252                     gausf->Draw("sames");
0253         }
0254     par0 = gausf->GetParameter(0);  // hieght
0255     par1 = gausf->GetParameter(1);  // mean value
0256     par2 = gausf->GetParameter(2);  // sigma
0257     chi2 = gausf->GetChisquare();   // chi2
0258     ndf  = gausf->GetNDF();         // ndf
0259      }
0260
0261      if (debug) {
0262         cv->Update();
0263         cv->Print("BugInfo/"+plot_id);
0264      }
0265      // exclude bad fitting!!
0266      if (ndf == 0)  ndf = 1.0 ;
0267      bool badfit1 =  ( (chi2/ndf) > 10. ) ? true : false ;
0268      bool badfit2 =  ( fabs(par1 - ave) > rms ) ? true : false ;
0269
0270      if ( badfit1 ) {
0271         if (debug) fprintf (logfile," Bad1={ %s, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f} \n"
0272                                 ,det,   par0,   par1,   par2,   chi2,    rms,    nu );
0273         par2 =0.0;
0274         par1 =0.0;
0275      }
0276      if ( badfit2 ) {
0277         if (debug) fprintf (logfile," Bad2={ %s, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f} \n"
0278                                 ,det,   par0,   par1,   par2,   chi2,    rms,    nu );
0279         par2 =0.0;
0280         par1 =0.0;
0281      }
0282
0283      if ( (nu < 15.) || (fabs(par1) > 1.5) ) {
0284         if (debug) fprintf (logfile," Bad3={ %s, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f, %-6.6f} \n"
0285                                 ,det,   par0,   par1,   par2,   chi2,    rms,    nu );
0286         par2 =0.0;
0287         par1 =0.0;
0288      }
0289
0290      if ( (par1 != 0.0) && (par2 > 0.0001) ){
0291         if ( (par2/par1) < 100. && (par2/par1) > 0.01) {
0292            yv1.push_back(par1);
0293        yv1e.push_back(par2);
0294        xv1.push_back(r);
0295         }
0296      }
0297      r = r + bsz;
0298
0299      delete heta_dphiPt_pjy;
0300  }
0301  cout<<" sliced and fit done!!! "<<endl;
0302
0303  // *************************************************
0304  //   Fill the array for histogram from vector !!!
0305  // *************************************************
0306
0307  const Int_t sz = xv1.size();
0308  Double_t ya1[sz], xa1[sz], ya1e[sz] ;
0309
0310  for(int j = 0; j < sz ; j++) {
0311     ya1[j]=yv1[j];
0312     xa1[j]=xv1[j];
0313     ya1e[j]=yv1e[j];
0314  }
0315
0316  // preliminary fitting - rejecting bad points
0317  // Chauvenet's Criterion to reject data from fitting
0318  int np2 = ( np1 == 1) ? np1 : np1 -1 ;
0319  TF1 *func0 = new TF1("func0", SeedPtFunction::linear, f1, f2, np2 );
0320  vector<int> goodYs = preFitter( func0, "func0", f1, f2, sz, xa1, ya1 );
0321  const Int_t sz2 = ( sz >=3 && goodYs.size() > 1 ) ? goodYs.size() + 2  : sz+2;
0322  cout<<" sz2 size "<<sz2<<endl;
0323  Double_t ya2[sz2] , xa2[sz2], ya2e[sz2],  xa2e[sz2] ;
0324
0325  xa2[0]= ini;
0326  ya2[0]= 0.;
0327  ya2e[0]= 0.;
0328  for(size_t j = 0; j < goodYs.size() ; j++) {
0329     int k = goodYs[j] ;
0330     xa2[j+1]  = xa1[k]  ;
0331     ya2[j+1]  = ya1[k]  ;
0332     ya2e[j+1] = ya1e[k] ;
0333     xa2e[j+1] = 0.0 ;
0334  }
0335  xa2[sz2-1] = fnl;
0336  ya2[sz2-1]  = 0.;
0337  ya2e[sz2-1] = 0.;
0338  xa2e[sz2-1] = 0.;
0339
0340  if ( sz2 < 5 ) {
0341     for (int j=1; j< sz2-1; j++ ) {
0342         xa2[j]  = xa1[j-1]  ;
0343         ya2[j]  = ya1[j-1]  ;
0344         ya2e[j] = ya1e[j-1] ;
0345         xa2e[j] = 0.0 ;
0346     }
0347  }
0348
0349
0350  gStyle->SetOptStat(kTRUE);
0351  gStyle->SetOptFit(0111);
0352  c1 = new TCanvas("c1");
0353  c1->SetFillColor(10);
0354  c1->SetGrid();
0355  c1->GetFrame()->SetFillColor(21);
0356  c1->GetFrame()->SetBorderSize(12);
0357  c1->cd();
0358
0359  double qq[2] = { 0., 0. } ;
0360  if ( sz2 < 5 ) {
0361      for (int j=1; j< sz2-1; j++ ) {
0362          qq[0] += ya2[j];
0363          qq[1] += ya2e[j]*ya2e[j] ;
0364      }
0365      qq[0] = qq[0]/(sz2-2) ;
0366      //qq[1] = sqrt( qq[1] ) /qq[0] ;
0367      qq[1] = sqrt( qq[1] /(sz2-2) ) ;
0368  }
0369
0370  TGraphErrors* eta_dphiPt_prf = new TGraphErrors(sz2,xa2,ya2,xa2e,ya2e);
0371  eta_dphiPt_prf->SetTitle(det_id+"  #eta vs. pT*#Delta#phi");
0372  eta_dphiPt_prf->SetMarkerColor(4);
0373  eta_dphiPt_prf->SetMarkerStyle(21);
0374  eta_dphiPt_prf->GetXaxis()->SetTitle(" #eta  ");
0375  eta_dphiPt_prf->GetYaxis()->SetTitle(" pT*#Delta#phi  ");
0376
0377  TF1 *func = new TF1("func", SeedPtFunction::linear, f1, f2, np1 );
0378
0379  if ( sz2 < 5 ) {
0380     func->FixParameter(0,qq[0]);
0381     func->FixParameter(1, 0.);
0382     func->FixParameter(2, 0.);
0383  }
0384  eta_dphiPt_prf->Draw("AP");
0385  eta_dphiPt_prf->Fit("func","RQ","sames",f1,f2);
0386
0387  double q0 = func->GetParameter(0);
0388  double q1 = func->GetParameter(1);
0389  double q2 = func->GetParameter(2);
0390
0391  c1->Update();
0392  c1->Print(plot01);
0393
0394  // ***********************************
0395  // the relative error parameterizaton
0396  // ***********************************
0397
0398  Double_t rErr2a[sz2] ;
0399  for ( int i=0; i< sz2; i++ ) {
0400      //rErr2a[i] = ya2e[i] / func->Eval( xa2[i] ) ;
0401      rErr2a[i] = ya2e[i] ;
0402  }
0403
0404  int np3 = ( np2 == 1) ? np2 : np2 -1 ;
0405  TF1 *func1 = new TF1("func1", SeedPtFunction::linear, f1, f2, np3 );
0406  vector<int> goodErrs = preFitter( func1, "func1", f1, f2, sz2, xa2, rErr2a );
0407
0408  const Int_t sz3 = ( sz2 >= 5 && goodErrs.size() > 1 ) ? goodErrs.size() + 2 : sz2;
0409  Double_t rErr3a[sz3] ;
0410  Double_t xa3[sz3]    ;
0411
0412  xa3[0] = ini ;
0413  rErr3a[0] = 0.0 ;
0414  for (int j=1; j< sz3-1 ; j++ ) {
0415      int k = goodErrs[j-1] ;
0416      xa3[j]    = xa2[k] ;
0417      rErr3a[j] = rErr2a[k] ;
0418  }
0419  xa3[sz3-1] = fnl ;
0420  rErr3a[sz3-1] = 0.0 ;
0421
0422  if ( sz3 < 5 ) {
0423     for (int j=1; j< sz2-1; j++ ) {
0424         xa3[j]    = xa2[j]  ;
0425         rErr3a[j] = rErr2a[j] ;
0426     }
0427  }
0428
0429  // get better fitting
0430  gStyle->SetOptStat(kTRUE);
0431  gStyle->SetOptFit(0111);
0432  c2 = new TCanvas("c2");
0433  c2->SetFillColor(10);
0434  c2->SetGrid();
0435  c2->GetFrame()->SetFillColor(21);
0436  c2->GetFrame()->SetBorderSize(12);
0437  c2->cd();
0438
0439  TGraph* eta_rErr = new TGraph(sz3,xa3,rErr3a);
0440  eta_rErr ->SetTitle(det_id+" Relative Error vs eta");
0441  eta_rErr ->SetMarkerColor(4);
0442  eta_rErr ->SetMarkerStyle(21);
0443  eta_rErr ->GetXaxis()->SetTitle(" #eta  ");
0444  eta_rErr ->GetYaxis()->SetTitle(" #sigma of pT*#Delta#phi  ");
0445
0446  int np4 = ( b2-b1 < 4) ? 1 : np2 ;
0447  TF1 *func2 = new TF1("func2", SeedPtFunction::linear, f1, f2, np4 );
0448  eta_rErr->Draw("AP");
0449  eta_rErr->Fit("func2","RQ","sames",f1,f2);
0450
0451  if ( sz2 < 5 ) {
0452     func2->FixParameter(0,qq[1]);
0453     func2->FixParameter(1, 0.);
0454     func2->FixParameter(2, 0.);
0455  }
0456
0457  double k0 = func2->GetParameter(0);
0458  double k1 = func2->GetParameter(1);
0459  double k2 = func2->GetParameter(2);
0460
0461  c2->Update();
0462  c2->Print(plot02);
0463
0464  fprintf (dfile," %s = cms.vdouble( %-6.3f, %-6.3f, %-6.3f, %-6.3f, %-6.3f, %-6.3f ), \n",
0465                  det,                   q0,     q1,     q2,     k0,     k1,     k2 );
0466
0467
0468 // ********************************************************************
0469 // Draw the origin scalar plots
0470 // ********************************************************************
0471
0472  gStyle->SetOptStat("ne");
0473  //gStyle->SetPalette(10);
0474  TCanvas *c3 = new TCanvas("c3","");
0475  c3->SetFillColor(10);
0476  c3->SetFillColor(10);
0477  heta_dphiPt->SetTitle("pT x #Delta#phi 2D histogram");
0478  heta_dphiPt->RebinY(2);
0479  heta_dphiPt->Draw("COLZ");
0480  heta_dphiPt->GetXaxis()->SetTitle(" #eta  ");
0481  heta_dphiPt->GetYaxis()->SetTitle("pT x #Delta#phi  ");
0482  c3->Update();
0483  c3->Print(plot03);
0484
0485   cout<<" Finished  !! "<<endl;
0486
0487   file->Close();
0488
0489   delete c1;
0490   delete c2;
0491   delete c3;
0492   delete gausf;
0493   delete func;
0494   delete func0;
0495   delete func1;
0496   delete func2;
0497   delete eta_dphiPt_prf;
0498   delete eta_rErr;
0499
0500 }
0501
0502 vector<int> SeedParaFit::preFitter( TF1* ffunc, TString falias, double h1, double h2, int asize, Double_t* xa, Double_t* ya ){
0503
0504     if ( debugAll ) {
0505        gStyle->SetOptStat(kFALSE);
0506        gStyle->SetOptFit(0111);
0507        c4 = new TCanvas("c4");
0508        c4->SetFillColor(10);
0509        c4->SetGrid();
0510        c4->GetFrame()->SetFillColor(21);
0511        c4->GetFrame()->SetBorderSize(12);
0512        c4->cd();
0513     }
0514
0515     vector<int> goodone;
0516     goodone.clear();
0517     if ( asize < 3 ) return goodone ;
0518
0519     TGraph* preh = new TGraph(asize,xa,ya);
0520     preh->SetTitle(" Test Fitting ");
0521     preh->SetMarkerColor(4);
0522     preh->SetMarkerStyle(21);
0523     preh->GetXaxis()->SetTitle(" #eta  ");
0524
0525     preh->Fit( falias ,"N0RQ","",h1,h2 );
0526
0527     // get the sigma of the pre-fitting
0528     double dv = 0.0;
0529     for (int j = 0; j < asize; j++) {
0530        double fite = ffunc->Eval( xa[j] ) ;
0531        dv += ( ya[j] - fite )*( ya[j] - fite ) ;
0532     }
0533     double sigma = sqrt( dv/( asize - 1.) );
0534
0535     for (int j = 0; j < asize; j++) {
0536
0537         double fitV = ffunc->Eval( xa[j] ) ;
0538         double width = fabs( ya[j] - fitV ) ;
0539
0540         /// gaussian probability for number of data expected
0541         bool reject = ptFunc->DataRejection( sigma , width, asize );
0542
0543         if ( !reject ) {
0544            goodone.push_back( j );
0545         }
0546     }
0547
0548     if ( debugAll ) {
0549        preh->Draw("AP");
0550        ffunc->Draw("sames");
0551        c4->Update();
0552        c4->Print(plot04);
0553     }
0554
0555     delete preh;
0556     if ( debugAll ) delete c4;
0557
0558     return goodone;
0559 }
0560
0561
0562 void SeedParaFit::PrintTitle(){
0563
0564   fprintf (dfile,"import FWCore.ParameterSet.Config as cms \n");
0565   fprintf (dfile,"ptSeedParameterization = cms.PSet( \n");
0566   fprintf (dfile," \n");
0567 }
0568
0569 void SeedParaFit::PrintEnd(){
0570
0571   fprintf (dfile," \n");
0572   fprintf (dfile,") ");
0573 }