test/MuonSeedPTAnalysis/SeedPtScale.C

0001 #include "SeedPtScale.h"
0002 //ClassImp(SeedPtScale);

0003
0004 SeedPtScale::SeedPtScale() {
0005
0006   dname = "SeedPtScale3" ;
0007   gSystem->mkdir(dname);
0008   suffixps = ".gif";
0009
0010   debug = false;
0011   // bin size for eta, unit bin size is 0.01

0012   // bin range => CSC: 0.9 - 2.7 ,  DT: 0.0 - 1.2 , OL: 0.7 - 1.3

0013   // # of bin        : 180            : 120           : 60

0014   // possible rbin value = 1, 2, 3, 4, 5, 10 ;

0015   xbsz = 0.01 ;
0016   //ybsz = 0.0001;

0017
0018 }
0019
0020 SeedPtScale::~SeedPtScale() {
0021
0022   delete c1;
0023   delete c2;
0024   delete c6;
0025   if (debug) {
0026      delete c3;
0027      delete c4;
0028      delete c5;
0029      delete cv;
0030   }
0031
0032   delete plot01;
0033   delete plot02;
0034   delete plot06;
0035   if (debug) {
0036      delete plot03;
0037      delete plot04;
0038      delete plot05;
0039   }
0040
0041 }
0042
0043 void SeedPtScale::PtScale( int type, int st, double h1, double h2, int idx, int np){
0044 //void SeedPtScale::PtScale( int type, int st, int b1, int b2, int idx, int np){

0045
0046   /*

0047    * Get the profile histogram for Pt*dphi vs pt or 1/phi

0048    *

0049    * Author:  Shih-Chuan Kao  --- UCR

0050    *

0051    * type = 1:CSC,  2:Overlap,  3:DT  4:DTSingle   5:CSCSingle

0052    * st   = 12 => station1 + station2 ; for type 4 & 5 => 12:ME12 or MB12

0053    * h1,h2= eta range

0054    * b1   = lst bin for lower fitting bound

0055    * b2   = 2nd bin for upper fitting bound

0056    * idx  = index for different eta range

0057    *

0058    */
0059
0060   int b1 = 0;
0061   int b2 = 0;
0062   if (type == 1 || type == 5) {
0063      b1 = (h1 - 0.9)/xbsz ;
0064      b2 = (h2 - 0.9)/xbsz ;
0065   }
0066   if (type ==3 || type == 4) {
0067      b1 = (h1 - 0.)/xbsz ;
0068      b2 = (h2 - 0.)/xbsz ;
0069      if (b1 == 0) b1 = 1 ;
0070   }
0071   if (type ==2 ) {
0072      b1 = (h1 - 0.7)/xbsz ;
0073      b2 = (h2 - 0.7)/xbsz ;
0074   }
0075   cout<<" Bin1:"<< b1 <<"  Bin2:"<<b2 <<endl;
0076
0077 // ********************************************************************

0078 // create a folder to store all histograms

0079 // ********************************************************************

0080
0081   float nsigmas = 1.5;
0082   char det[8];
0083   char plotname[9];
0084
0085   if (type == 1)   {
0086      if (st == 01) {
0087        sprintf(det,"CSC_01");
0088        sprintf(plotname,"CSC_01_%d",idx);
0089      } else {
0090        sprintf(det,"CSC_%d",st);
0091        sprintf(plotname,"CSC_%d_%d",st,idx);
0092      }
0093      Dir = "CSC_All";
0094      det_id = det;
0095   }
0096   if (type == 2)   {
0097      sprintf(det,"OL_%d",st);
0098      sprintf(plotname,"OL_%d",st);
0099      Dir = "OL_All";
0100      det_id = det;
0101   }
0102   if (type == 3)   {
0103      sprintf(det,"DT_%d",st);
0104      sprintf(plotname,"DT_%d_%d",st,idx);
0105      Dir = "DT_All";
0106      det_id = det;
0107   }
0108   if (type == 4)   {
0109      sprintf(det,"SMB_%d",st);
0110      sprintf(plotname,"SMB_%d",st);
0111      Dir = "MB_All";
0112      det_id = det;
0113   }
0114   if (type == 5)   {
0115      sprintf(det,"SME_%d",st);
0116      sprintf(plotname,"SME_%d",st);
0117      Dir = "ME_All";
0118      det_id = det;
0119   }
0120
0121   pname = plotname;
0122
0123   cout<<" Set Path" <<endl;
0124   TString thePath  = Dir+"/"+det_id+"_eta_rdphiPt";
0125   TString thePath1 = Dir+"/"+det_id+"_eta_rdphi";
0126   if (type ==2 ) {
0127       thePath  = Dir+"/"+det_id+"_eta_rdphiPtA";
0128       thePath1 = Dir+"/"+det_id+"_eta_rdphiA";
0129   }
0130
0131   plot01 = pname+"_pt_ptxdPhi"+suffixps;
0132   plot02 = pname+"_pt_odphi"+suffixps;
0133   plot03 = pname+"_pt_dphi"+suffixps;
0134   plot04 = pname+"_corr"+suffixps;
0135   plot05 = pname+"_normal"+suffixps;
0136   plot06 = pname+"_odphi_ptxdPhi"+suffixps;
0137
0138   int ptarr[] ={10,20,50,100,150,200,350,500};
0139   vector<int> ptlist( ptarr, ptarr+8 );
0140   /*

0141   vector<int> ptlist;

0142   for (int ii = 0; ii < 7; ii++) {

0143       ptlist.push_back( ptarr[ii] );

0144   }*/
0145   const Int_t sz = ptlist.size();
0146
0147   if ( debug ) gSystem->mkdir("BugInfo");
0148
0149 // *****************************************************************

0150 // main program -- looping pt from 10 -> 1000

0151 // *****************************************************************

0152
0153    float L1,H1 = 0.0;
0154    float mean,rms,ent =0.0;
0155    bool fitfail = false ;
0156
0157    Double_t mean1[sz], sigma1[sz] ;
0158    Double_t pt[sz],     ptErr[sz] ;
0159    Double_t dphi[sz], dphiErr[sz] ;
0160    Double_t dphi_real[sz], dphiErr_real[sz] ;
0161    Double_t ophi[sz], ophiErr[sz] ;
0162    Double_t ratio[sz],   corr[sz] ;
0163
0164    //SeedPtFunction* fptr = new SeedPtFunction() ;

0165    //TF1 *gausf = new TF1("gausf", SeedFitFunc, &SeedPtFunction::fgaus, -1., 3., 3, "SeedPtFunction", "fgaus" );

0166    TF1 *gausf = new TF1("gausf", SeedPtFunction::fgaus, -1., 3., 3  );
0167
0168    if (debug) {
0169       cv = new TCanvas("cv");
0170       gStyle->SetOptStat(kTRUE);
0171       gStyle->SetOptFit(0111);
0172    }
0173
0174    for ( int i=0; i<sz; i++) {
0175
0176        if ( debug ) {
0177           sprintf(dbplot, "dbug_%d.gif",i);
0178           plot_id = dbplot;
0179           cv->cd();
0180        }
0181
0182        char filename[18];
0183        sprintf(filename,"para%d.root",ptlist[i]);
0184        TString fname = filename;
0185        pt[i] = ptlist[i] ;
0186
0187        // pt*dphi

0188        TFile *file      = TFile::Open(fname);
0189        TH2F* hdphiPt    = (TH2F *) file->Get(thePath);
0190        TH1D* hdphiPt_py = hdphiPt->ProjectionY("hdphiPt_py",b1,b2);
0191        if ( debug ) hdphiPt_py->Draw();
0192
0193        mean = hdphiPt_py->GetMean();
0194        rms  = hdphiPt_py->GetRMS();
0195        ent  = hdphiPt_py->GetEntries();
0196        L1 = mean - (1.*rms);
0197        H1 = mean + (1.*rms);
0198        fitfail = false ;
0199
0200        gausf->SetParLimits(0, 1., 0.8*ent);
0201        gausf->SetParLimits(1, mean-rms, mean+rms);
0202        gausf->SetParLimits(2, 0.2*rms, 2.*rms);
0203        gausf->SetParameter(1,    mean);
0204        gausf->SetParameter(2,     rms);
0205
0206        if ( !debug ) hdphiPt_py->Fit( "gausf" ,"N0RQC","", L1, H1 );
0207        if (  debug ) deBugger(gausf, "gausf", hdphiPt_py, L1, H1 , 2 );
0208        fitfail = BadFitting(gausf, hdphiPt_py) ;
0209        if ( fitfail ) hdphiPt_py->Fit("gausf","N0Q" );
0210        mean1[i]  = gausf->GetParameter(1);
0211        sigma1[i] = gausf->GetParameter(2);
0212        L1 = mean1[i] - (nsigmas * sigma1[i]);
0213        H1 = mean1[i] + (nsigmas * sigma1[i]);
0214
0215        if ( !debug ) hdphiPt_py->Fit( "gausf" ,"N0RQC","",L1, H1);
0216        if (  debug ) deBugger(gausf, "gausf", hdphiPt_py, L1, H1 , 4 );
0217        mean1[i]  = gausf->GetParameter(1);
0218        sigma1[i] = gausf->GetParameter(2);
0219
0220        fitfail = false;
0221
0222        // dphi information , no fit just use the mean and rms

0223        TH2F* hdphi   = (TH2F *) file->Get(thePath1);
0224        TH1D* dphi_py = hdphi->ProjectionY("dphi_py",b1,b2);
0225
0226        // calculated dphi

0227        dphi[i]    = mean1[i]/pt[i];
0228        dphiErr[i] = dphi_py->GetRMS();
0229        // real dphi

0230        dphi_real[i]    = 1000*dphi_py->GetMean();
0231        dphiErr_real[i] = 1000*dphi_py->GetRMS();
0232        ophi[i]    = 1./dphi[i] ;
0233        if ( i > 0 ) {
0234           if ( ophi[i] < ophi[i-1] ) {
0235               dphi[i] = mean1[i-1] / pt[i] ;
0236               ophi[i]    = 1./dphi[i] ;
0237           }
0238        }
0239        ophiErr[i] = dphiErr[i] / (dphi[i]*dphi[i]);
0240
0241        cout<<" pt:"<< pt[i] <<" mean = "<< mean1[i] <<" df :"<< dphi[i] <<endl;
0242
0243        delete hdphiPt ;
0244        delete hdphiPt_py ;
0245        delete hdphi ;
0246        delete dphi_py ;
0247        file->Close();
0248    }
0249
0250  // *************************************************

0251  //   Fill the array for histogram from vector !!!

0252  // *************************************************

0253
0254  gSystem->cd(dname);
0255  FILE *dfile = fopen("MuonSeedPtScale.py","a");
0256  FILE *dfile2 = fopen("MuonSeeddPhiScale.py","a");
0257
0258  // 1. pt vs pt*dphi

0259  gStyle->SetOptStat(kFALSE);
0260  gStyle->SetOptFit(0111);
0261  c1 = new TCanvas("c1");
0262  c1->SetFillColor(10);
0263  c1->SetGrid();
0264  c1->GetFrame()->SetFillColor(21);
0265  c1->GetFrame()->SetBorderSize(12);
0266  c1->cd();
0267
0268  TGraphErrors* pt_dphiPt = new TGraphErrors(sz,pt,mean1,ptErr,sigma1);
0269  pt_dphiPt->SetMarkerColor(4);
0270  pt_dphiPt->SetMarkerStyle(21);
0271  pt_dphiPt->SetTitle(plot01);
0272  pt_dphiPt->GetXaxis()->SetTitle(" pt ");
0273  pt_dphiPt->GetYaxis()->SetTitle(" pT*#Delta#phi  ");
0274
0275  TF1 *func0 = new TF1("func0", SeedPtFunction::fitf, 5, 1100, np );
0276  pt_dphiPt->Fit( func0 ,"R","",5,1100);
0277
0278  double q0 = func0->GetParameter(0);
0279  double q1 = func0->GetParameter(1);
0280  double q2 = func0->GetParameter(2);
0281  double t1 = q1/q0;
0282  double t2 = q2/q0;
0283
0284  //fprintf (dfile,"  %s_%d_scale = cms.vdouble( %f, %f, %f, %f, %f ), \n"

0285  //                  ,det ,idx                 ,q0, q1, q2, t1, t2 );

0286  fprintf (dfile,"  %s_%d_scale = cms.vdouble(  %f, %f ), \n",
0287                    det, idx,                   t1, t2 );
0288
0289  pt_dphiPt->Draw("AP");
0290  c1->Update();
0291  c1->Print(plot01);
0292
0293  // get the ratio and correction

0294  for (int j=0; j<sz; j++) {
0295      double theX = 1./ (10. + ophi[j]) ;
0296      ratio[j] = 1. + (t1 * theX) + ( t2 * theX * theX ) ;
0297      corr[j] = 1./ratio[j] ;
0298  }
0299
0300
0301  if ( debug ) {
0302     // 2. pt vs 1/dphi

0303     gStyle->SetOptStat(kFALSE);
0304     gStyle->SetOptFit(0111);
0305     c2 = new TCanvas("c2");
0306     c2->SetFillColor(10);
0307     c2->SetGrid();
0308     c2->GetFrame()->SetFillColor(21);
0309     c2->GetFrame()->SetBorderSize(12);
0310     c2->cd();
0311
0312     TGraph* ophi_Pt = new TGraph(sz, pt, ophi);
0313     ophi_Pt->SetMarkerColor(4);
0314     ophi_Pt->SetMarkerStyle(21);
0315     ophi_Pt->SetTitle(plot02);
0316     ophi_Pt->GetXaxis()->SetTitle(" pt ");
0317     ophi_Pt->GetYaxis()->SetTitle(" 1/*#Delta#phi  ");
0318
0319     TF1 *func2 = new TF1("func2", SeedPtFunction::linear, 5, 1000, 2 );
0320     ophi_Pt->Fit( func2 ,"R","",0,1000);
0321
0322     ophi_Pt->Draw("AP");
0323     c2->Update();
0324     c2->Print(plot02);
0325
0326
0327     // 4. pt vs correction factor

0328     gStyle->SetOptStat(kFALSE);
0329     gStyle->SetOptFit(0111);
0330     c4 = new TCanvas("c4");
0331     c4->SetFillColor(10);
0332     c4->SetGrid();
0333     c4->GetFrame()->SetFillColor(21);
0334     c4->GetFrame()->SetBorderSize(12);
0335     c4->cd();
0336
0337     //TGraph* cor_ophi = new TGraph(sz,ophi,corr);

0338     TGraph* cor_pt = new TGraph(sz,pt,corr);
0339     cor_pt->SetMarkerColor(4);
0340     cor_pt->SetMarkerStyle(21);
0341     cor_pt->SetTitle(plot04);
0342     cor_pt->GetXaxis()->SetTitle(" Pt  ");
0343     cor_pt->GetYaxis()->SetTitle(" correction factor ");
0344
0345     cor_pt->Draw("AP");
0346     c4->Update();
0347     c4->Print(plot04);
0348
0349     //5. 1/dphi vs normalized factor

0350     gStyle->SetOptStat(kFALSE);
0351     gStyle->SetOptFit(0111);
0352     c5 = new TCanvas("c5");
0353     c5->SetFillColor(10);
0354     c5->SetGrid();
0355     c5->GetFrame()->SetFillColor(21);
0356     c5->GetFrame()->SetBorderSize(12);
0357     c5->cd();
0358
0359     TGraph* nom_ophi = new TGraph(sz,ophi,ratio);
0360     nom_ophi->SetMarkerColor(4);
0361     nom_ophi->SetMarkerStyle(21);
0362     nom_ophi->SetTitle(plot05);
0363     nom_ophi->GetXaxis()->SetTitle(" 1/#Delta#phi  ");
0364     nom_ophi->GetYaxis()->SetTitle(" normalized factor ");
0365
0366     nom_ophi->Draw("AP");
0367     c5->Update();
0368     c5->Print(plot05);
0369
0370     delete nom_ophi;
0371     delete cor_pt;
0372     delete ophi_Pt;
0373     delete func2;
0374  }
0375
0376  // 3. pt vs dphi

0377  gStyle->SetOptStat(kFALSE);
0378  gStyle->SetOptFit(0111);
0379  c3 = new TCanvas("c3");
0380  c3->SetFillColor(10);
0381  c3->SetGrid();
0382  c3->GetFrame()->SetFillColor(21);
0383  c3->GetFrame()->SetBorderSize(12);
0384  c3->cd();
0385
0386  TGraphErrors* dphi_Pt = new TGraphErrors(sz,pt,dphi_real,ptErr,dphiErr_real);
0387  dphi_Pt->SetMarkerColor(4);
0388  dphi_Pt->SetMarkerStyle(21);
0389  dphi_Pt->SetTitle(plot03);
0390  dphi_Pt->GetXaxis()->SetTitle(" pt ");
0391  dphi_Pt->GetYaxis()->SetTitle(" #Delta#phi (mrad) ");
0392
0393  dphi_Pt->Draw("AP");
0394  c3->Update();
0395  c3->Print(plot03);
0396
0397  // 6. ptxdPhi vs 1/dphi

0398  gStyle->SetOptStat(kFALSE);
0399  gStyle->SetOptFit(0111);
0400  c6 = new TCanvas("c6");
0401  c6->SetFillColor(10);
0402  c6->SetGrid();
0403  c6->GetFrame()->SetFillColor(21);
0404  c6->GetFrame()->SetBorderSize(12);
0405  c6->cd();
0406
0407  // get the 1/dphi from ptxdPhi fitting

0408  Double_t odf_fit[sz] ;
0409  for (int j=0; j<sz; j++) {
0410      double theX = 1./ (10. + pt[j]) ;
0411      double BFit = q0 + (q1 * theX) + ( q2 * theX * theX ) ;
0412      odf_fit[j] =  pt[j] / BFit ;
0413  }
0414
0415  TGraphErrors* ophi_dphiPt = new TGraphErrors(sz, odf_fit, mean1, ptErr, sigma1);
0416  ophi_dphiPt->SetMarkerColor(4);
0417  ophi_dphiPt->SetMarkerStyle(21);
0418  ophi_dphiPt->SetTitle(plot02);
0419  ophi_dphiPt->GetXaxis()->SetTitle(" 1/#Delta#phi ");
0420  ophi_dphiPt->GetYaxis()->SetTitle(" ptx*#Delta#phi  ");
0421
0422  TF1 *func6 = new TF1("func6", SeedPtFunction::fitf, 5, 25000, np );
0423  if ( type < 4 ) {
0424     func6->SetParLimits(1, -10., -0.0001);
0425     func6->SetParameter(1, q1);
0426  }
0427  ophi_dphiPt->Fit( func6 ,"R","", 0, 25000);
0428
0429  double r0 = func6->GetParameter(0);
0430  double r1 = func6->GetParameter(1);
0431  double r2 = func6->GetParameter(2);
0432  double s1 = r1/r0;
0433  double s2 = r2/r0;
0434
0435  fprintf (dfile2,"  %s_%d_scale = cms.vdouble( %f, %f, %f, %f, %f ), \n"
0436                   ,det ,idx                   ,r0, r1, r2, s1, s2 );
0437
0438  ophi_dphiPt->Draw("AP");
0439  c6->Update();
0440  c6->Print(plot06);
0441
0442  fclose(dfile);
0443  fclose(dfile2);
0444
0445  gSystem->cd("../");
0446
0447
0448  delete gausf;
0449  delete func0;
0450  delete func6;
0451
0452  delete dphi_Pt;
0453  delete pt_dphiPt;
0454  delete ophi_dphiPt;
0455
0456  //gROOT->Reset();

0457  //gROOT->ProcessLine(".q");

0458
0459 }
0460
0461 void SeedPtScale::deBugger( TF1* fitfunc, TString funcName ,TH1D* histo, double L2, double H2, int cr ) {
0462
0463       fitfunc->SetLineStyle(cr);
0464       fitfunc->SetLineColor(cr);
0465       double D2 = H2 - L2 ;
0466
0467       histo->SetAxisRange(L2-D2,H2+D2,"X");
0468       histo->Fit( funcName , "RQ", "sames", L2, H2 );
0469
0470       cv->Update();
0471       cv->Print("BugInfo/"+plot_id);
0472 }
0473
0474 bool SeedPtScale::BadFitting( TF1* fitfunc, TH1D* histo ) {
0475
0476     double p0 = fitfunc->GetParameter(0);
0477     double p1 = fitfunc->GetParameter(1);
0478     double p2 = fitfunc->GetParameter(2);
0479     double mean_h = histo->GetMean();
0480     double rms_h  = histo->GetRMS();
0481     double ent_h  = histo->GetEntries();
0482
0483     bool badfit = false ;
0484     if ( p0 > 0.9* ent_h )               badfit = true;
0485     if ( p0 < 1. && ent_h > 10. )        badfit = true;
0486     if ( p1 < 0. )                       badfit = true;
0487     if ( fabs(p1 - mean_h) >  rms_h*3. ) badfit = true;
0488     if ( p2 < 0. )                       badfit = true;
0489     if ( p2 > rms_h*3. )                 badfit = true;
0490
0491     if (badfit) cout<<" Fit Fails !!!" <<endl;
0492     return badfit;
0493 }
0494