File indexing completed on 2023-03-17 11:20:49
0001 #include "SeedParaFit.h"
0002
0003 SeedParaFit::SeedParaFit() {
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 fName = "mix31x_v3.root";
0015
0016 dname = "Pt31x_v3";
0017
0018 suffixps = ".gif";
0019
0020 nsigmas = 1.5;
0021
0022 debug = false;
0023
0024 debugAll = false ;
0025
0026
0027
0028
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
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
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
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
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
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
0130
0131
0132
0133
0134 if ( debug ) gSystem->mkdir("BugInfo");
0135
0136
0137
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
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
0223 if ( nu >= 10 ) {
0224
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);
0238 if ( par0 < 1. || par0 > nu*0.8 ) heta_dphiPt_pjy->Fit("gausf","N0Q" );
0239 par1 = gausf->GetParameter(1);
0240 par2 = gausf->GetParameter(2);
0241
0242 L1 = par1 - (nsigmas * par2);
0243 H1 = par1 + (nsigmas * par2);
0244
0245
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);
0255 par1 = gausf->GetParameter(1);
0256 par2 = gausf->GetParameter(2);
0257 chi2 = gausf->GetChisquare();
0258 ndf = gausf->GetNDF();
0259 }
0260
0261 if (debug) {
0262 cv->Update();
0263 cv->Print("BugInfo/"+plot_id);
0264 }
0265
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
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
0317
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
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
0396
0397
0398 Double_t rErr2a[sz2] ;
0399 for ( int i=0; i< sz2; i++ ) {
0400
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
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
0470
0471
0472 gStyle->SetOptStat("ne");
0473
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
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
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 }