BeamSpotProducer/test/BeamFit.cc

0001 #include <iostream>
0002 #include <stdlib.h>
0003 #include <cmath>
0004 #include "TROOT.h"
0005 #include "TFile.h"
0006 #include "TRint.h"
0007 #include "TH1.h"
0008 #include "TF1.h"
0009 #include "TF2.h"
0010 #include "TH2.h"
0011 #include "TCanvas.h"
0012 #include "TMinuit.h"
0013 #include "TSystem.h"
0014 #include "TGraphErrors.h"
0015 #include "TMatrixD.h"
0016 #include "TMatrixDSym.h"
0017 #include "NtupleHelper.h"
0018 #include "global.h"
0019 //#include "beamfit_fcts.h"
0020
0021 //const char par_name[dim][20]={"z0  ","sigma ","emmitance","beta*"};
0022 const char par_name[dim][20] = {
0023     "z0  ", "sigma ", "x0 ", "y0 ", "dxdz ", "dydz ", "sigmaBeam ", "c0  ", "c1    ", "emittance ", "betastar "};
0024
0025 Double_t params[dim], errparams[dim];
0026 Double_t sfpar[dim], errsfpar[dim];
0027
0028 constexpr Double_t step[dim] = {1.e-5, 1.e-5, 1.e-3, 1.e-3, 1.e-3, 1.e-3, 1.e-5, 1.e-5, 1.e-5, 1.e-5, 1.e-5};
0029 zData zdata;  //!
0030 int tmpNtrks_ = 0;
0031 int fnthite = 0;
0032 TMatrixD ftmp;
0033
0034 Double_t fd0cut = 4.0;
0035
0036 int sequence[11] = {500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000};
0037 Double_t sequenceD[11] = {500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000};
0038 Double_t zeros[11] = {0.};
0039
0040 int iloop_ = 0;
0041
0042 int Nsequence_ = 11;
0043 Double_t res_d0phi_X[11] = {0.};
0044 Double_t res_d0phi_Xerr[11] = {0.};
0045 Double_t res_d0phi_Y[11] = {0.};
0046 Double_t res_d0phi_Yerr[11] = {0.};
0047 Double_t res_d0phi_dXdz[11] = {0.};
0048 Double_t res_d0phi_dXdzerr[11] = {0.};
0049 Double_t res_d0phi_dYdz[11] = {0.};
0050 Double_t res_d0phi_dYdzerr[11] = {0.};
0051
0052 Double_t res_Z[11] = {0.};
0053 Double_t res_sigmaZ[11] = {0.};
0054 Double_t res_Zerr[11] = {0.};
0055 Double_t res_sigmaZerr[11] = {0.};
0056 Double_t res_Z_lh[11] = {0.};
0057 Double_t res_sigmaZ_lh[11] = {0.};
0058 Double_t res_Zerr_lh[11] = {0.};
0059 Double_t res_sigmaZerr_lh[11] = {0.};
0060
0061 TGraphErrors *g_X;
0062 TGraphErrors *g_Y;
0063 TGraphErrors *g_dXdz;
0064 TGraphErrors *g_dYdz;
0065 TGraphErrors *g_Z;
0066 TGraphErrors *g_sigmaZ;
0067 TGraphErrors *g_Z_lh;
0068 TGraphErrors *g_sigmaZ_lh;
0069
0070 Double_t zdis(Double_t z, Double_t sigma, Double_t *parms) {
0071   //---------------------------------------------------------------------------
0072   //  This is a simple function to parameterize the z-vertex distribution. This
0073   // is parametrized by a simple normalized gaussian distribution.
0074   //---------------------------------------------------------------------------
0075
0076   Double_t sig = sqrt(sigma * sigma + parms[Par_Sigma] * parms[Par_Sigma]);
0077   Double_t result = (exp(-((z - parms[Par_Z0]) * (z - parms[Par_Z0])) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0078   return result;
0079 }
0080
0081 void zfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0082   //----------------------------------------------------------------------------------
0083   // this is the function used by minuit to do the unbinned fit to the z distribution
0084   //----------------------------------------------------------------------------------
0085   f = 0.0;
0086   for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0087     f = log(zdis(iter->Z, iter->SigZ, params)) + f;
0088   }
0089   f = -2.0 * f;
0090   return;
0091 }
0092
0093 float betaf(float betastar, float emmitance, float z, float z0) {
0094   float x = sqrt(emmitance * (betastar + (((z - z0) * (z - z0)) / betastar)));
0095   return x;
0096 }
0097
0098 Double_t ddis(Double_t z, Double_t sigma, Double_t d, Double_t sigmad, Double_t *parms) {
0099   //---------------------------------------------------------------------------
0100   // This is a simple function to parameterize the sigma of the beam at a given z.
0101   // This is parametrized by a simple normalized gaussian distribution.
0102   //---------------------------------------------------------------------------
0103   Double_t sig = betaf(parms[Par_beta], parms[Par_eps], z, parms[Par_Z0]);
0104   sig = sqrt(sig * sig + sigmad * sigmad);
0105   Double_t result = (exp(-(d * d) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0106   return result;
0107 }
0108
0109 Double_t allddis(Double_t z, Double_t d, Double_t sigmad, Double_t phi0, Double_t *parms) {
0110   //---------------------------------------------------------------------------
0111   // This is a simple function to parameterize the beam parameters
0112   // This is parametrized by a simple normalized gaussian distribution.
0113   //---------------------------------------------------------------------------
0114   Double_t sig =
0115       sqrt(parms[Par_Sigbeam] * parms[Par_Sigbeam] + 10 * sigmad * sigmad);  //9 //2.5 factor for full simu data
0116   Double_t dprime =
0117       d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0118   Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0119   return result;
0120 }
0121
0122 Double_t allddisBeta(Double_t z, Double_t sigma, Double_t d, Double_t sigmad, Double_t phi0, Double_t *parms) {
0123   //---------------------------------------------------------------------------
0124   // This is a simple function to parameterize the beam parameters
0125   // This is parametrized by a simple normalized gaussian distribution.
0126   //---------------------------------------------------------------------------
0127   Double_t sigmabeam =
0128       sqrt(parms[Par_eps] * (parms[Par_beta] + (((z - parms[Par_Z0]) * (z - parms[Par_Z0])) / parms[Par_beta])));
0129
0130   Double_t sig = sqrt(sigmabeam * sigmabeam + sigmad * sigmad);
0131
0132   Double_t dprime =
0133       d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0134   Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0135   return result;
0136 }
0137
0138 Double_t allddis2(Double_t z, Double_t sigma, Double_t d, Double_t pt, Double_t phi0, Double_t *parms) {
0139   //---------------------------------------------------------------------------
0140   //
0141   //---------------------------------------------------------------------------
0142   Double_t sigmad = parms[Par_c0] + parms[Par_c1] / pt;
0143   Double_t sig = sqrt(parms[Par_Sigbeam] * parms[Par_Sigbeam] + sigmad * sigmad);
0144   Double_t dprime =
0145       d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0146   Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0147   return result;
0148 }
0149
0150 void dfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0151   //----------------------------------------------------------------------------------
0152   // this is the function used by minuit to do the unbinned fit to the IP distribution
0153   //----------------------------------------------------------------------------------
0154   f = 0.0;
0155   for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0156     //if(iter->weight2 == 0) continue;
0157
0158     f = log(allddis(iter->Z, iter->D, iter->SigD, iter->Phi, params)) + f;
0159   }
0160   f = -2.0 * f;
0161   return;
0162 }
0163
0164 void dfcnbeta(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0165   //----------------------------------------------------------------------------------
0166   // this is the function used by minuit to do the unbinned fit to the IP distribution
0167   //----------------------------------------------------------------------------------
0168   f = 0.0;
0169   for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0170     f = log(allddisBeta(iter->Z, iter->SigZ, iter->D, iter->SigD, iter->Phi, params)) + f;
0171   }
0172   f = -2.0 * f;
0173   return;
0174 }
0175
0176 void dfcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0177   //----------------------------------------------------------------------------------
0178   // this is the function used by minuit to do the unbinned fit to the IP distribution
0179   //----------------------------------------------------------------------------------
0180   f = 0.0;
0181   for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0182     f = log(allddis2(iter->Z, iter->SigZ, iter->D, iter->Pt, iter->Phi, params) * zdis(iter->Z, iter->SigZ, params)) +
0183         f;
0184   }
0185   f = -2.0 * f;
0186   return;
0187 }
0188
0189 void cfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0190   //-----------------------------------------------------------------------------------
0191   // this is the function used by minuit to do the unbinned finned combined in z an IP.
0192   //-----------------------------------------------------------------------------------
0193   f = 0.0;
0194   for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0195     f = log(allddis(iter->Z, iter->D, iter->SigD, iter->Phi, params) * zdis(iter->Z, iter->SigZ, params)) + f;
0196   }
0197   f = -2.0 * f;
0198   return;
0199 }
0200
0201 void fit(TMatrixD &x, TMatrixDSym &V) {
0202   TMatrixDSym Vint(4);
0203   TMatrixD b(4, 1);
0204   Double_t weightsum = 0;
0205   tmpNtrks_ = 0;
0206
0207   Vint.Zero();
0208   b.Zero();
0209   TMatrixD g(4, 1);
0210   TMatrixDSym temp(4);
0211   for (zDataConstIter i = zdata.begin(); i != zdata.end(); ++i) {
0212     //std::cout << "weight  " << sqrt(i->weight2) << "\n";
0213     if (i->weight2 == 0)
0214       continue;
0215
0216     g(0, 0) = sin(i->Phi);
0217     g(1, 0) = -cos(i->Phi);
0218     g(2, 0) = i->Z * g(0, 0);
0219     g(3, 0) = i->Z * g(1, 0);
0220
0221     temp.Zero();
0222     for (int j = 0; j < 4; ++j) {
0223       for (int k = j; k < 4; ++k) {
0224         temp(j, k) += g(j, 0) * g(k, 0);
0225       }
0226     }
0227     double sigmabeam2 = 0.002 * 0.002;
0228     //double sigma2 = sigmabeam2 +  (i->SigD)* (i->SigD) / i->weight2;
0229     double sigma2 = sigmabeam2 + (i->SigD) * (i->SigD);
0230
0231     TMatrixD ftmptrans(1, 4);
0232     ftmptrans = ftmptrans.Transpose(ftmp);
0233     TMatrixD dcor = ftmptrans * g;
0234
0235     bool pass = true;
0236     if ((fnthite > 0) && (std::abs(i->D - dcor(0, 0)) > fd0cut))
0237       pass = false;
0238     //if ( tmpNtrks_ < 10 && fnthite>0 ) {
0239     //std::cout << " d0 = " << i->D << " dcor(0,0)= " << dcor(0,0) << " fd0cut=" << fd0cut << std::endl;
0240     //}
0241
0242     if (pass) {
0243       Vint += (temp * (1 / sigma2));
0244       b += ((i->D / sigma2) * g);
0245       weightsum += sqrt(i->weight2);
0246       tmpNtrks_++;
0247     }
0248   }
0249   Double_t determinant;
0250   V = Vint.InvertFast(&determinant);
0251   x = V * b;
0252   ftmp = x;
0253
0254   std::cout << "number of tracks used in this iteration = " << tmpNtrks_ << std::endl;
0255   std::cout << "Sum of all weights:" << weightsum << "\n";
0256   std::cout << "x0                :" << x(0, 0) << " +- " << sqrt(V(0, 0)) << " cm \n";
0257   std::cout << "y0                :" << x(1, 0) << " +- " << sqrt(V(1, 1)) << " cm \n";
0258   std::cout << "x slope           :" << x(2, 0) << " +- " << sqrt(V(2, 2)) << " cm \n";
0259   std::cout << "y slope           :" << x(3, 0) << " +- " << sqrt(V(3, 3)) << " cm \n";
0260
0261   res_d0phi_X[iloop_] = x(0, 0);
0262   res_d0phi_Xerr[iloop_] = sqrt(V(0, 0));
0263   res_d0phi_Y[iloop_] = x(1, 0);
0264   res_d0phi_Yerr[iloop_] = sqrt(V(1, 1));
0265   res_d0phi_dXdz[iloop_] = x(2, 0);
0266   res_d0phi_dXdzerr[iloop_] = sqrt(V(2, 2));
0267   res_d0phi_dYdz[iloop_] = x(3, 0);
0268   res_d0phi_dYdzerr[iloop_] = sqrt(V(3, 3));
0269 }
0270
0271 //Double_t cutOnD(const TMatrixD& x, Double_t fd0cut)
0272 //{
0273 //  TMatrixD g(1,4);
0274 //  Double_t weightsum = 0;
0275 //  for(zDataConstIter i = zdata.begin() ; i != zdata.end() ; ++i) {
0276 //    g(0,0) = - sin(i->Phi);
0277 //    g(0,1) =   cos(i->Phi);
0278 //    g(0,2) = i->Z * g(0,0);
0279 //    g(0,3) = i->Z * g(0,1);
0280 //    TMatrixD dcor = g * x;
0281 //    //std::cout << dcor.GetNrows() << " , " << dcor.GetNcols() << "\n";
0282 //    if(std::abs(i->D -  dcor(0,0)) > fd0cut) {
0283 //      i->weight2 = 0;
0284 //    } else {
0285 //      i->weight2 = 1;
0286 //      weightsum += 1;
0287 //    }
0288 //  }
0289 //  return weightsum;
0290 //}
0291
0292 Double_t cutOnChi2(const TMatrixD &x, Double_t chi2cut) {
0293   double sigmabeam2 = 0.002 * 0.002;
0294   TMatrixD g(1, 4);
0295   Double_t weightsum = 0;
0296   for (zDataIter i = zdata.begin(); i != zdata.end(); ++i) {
0297     g(0, 0) = sin(i->Phi);
0298     g(0, 1) = -cos(i->Phi);
0299     g(0, 2) = i->Z * g(0, 0);
0300     g(0, 3) = i->Z * g(0, 1);
0301     TMatrixD dcor = g * x;
0302     //std::cout << dcor.GetNrows() << " , " << dcor.GetNcols() << "\n";
0303     Double_t chi2 = (i->D - dcor(0, 0)) * (i->D - dcor(0, 0)) / (sigmabeam2 + (i->SigD) * (i->SigD));
0304     if (chi2 > chi2cut) {
0305       i->weight2 = 0;
0306     } else {
0307       i->weight2 = 1;
0308       weightsum += 1;
0309     }
0310   }
0311   //std::cout << weightsum << "\n";
0312   return weightsum;
0313 }
0314
0315 void fitAll() {
0316   std::cout << "starting fits" << std::endl;
0317
0318   // chi2 fit for Z
0319   TH1F *h1z = new TH1F("h1z", "z distribution", 100, -50., 50.);
0320   for (zDataConstIter i = zdata.begin(); i != zdata.end(); ++i) {
0321     h1z->Fill(i->Z, i->SigZ);
0322   }
0323   h1z->Sumw2();
0324   h1z->Fit("gaus");
0325   //std::cout << "fitted "<< std::endl;
0326
0327   TF1 *fgaus = h1z->GetFunction("gaus");
0328   //std::cout << "got function" << std::endl;
0329   res_Z[iloop_] = fgaus->GetParameter(1);
0330   res_sigmaZ[iloop_] = fgaus->GetParameter(2);
0331   res_Zerr[iloop_] = fgaus->GetParError(1);
0332   res_sigmaZerr[iloop_] = fgaus->GetParError(2);
0333
0334   std::cout << "======== End of Chi2 Fit for Z ========\n" << std::endl;
0335
0336   // LH fit for Z
0337   TMinuit *gmMinuit = new TMinuit(2);
0338   gmMinuit->SetFCN(zfcn);
0339   std::cout << "SetFCN done" << std::endl;
0340   int ierflg = 0;
0341   sfpar[Par_Sigma] = 7.55;   //mygaus->GetParameter(2);
0342   sfpar[Par_Z0] = 0.;        //r/mygaus->GetParameter(1);
0343   errsfpar[Par_Sigma] = 0.;  //r/mygaus->GetParError(2);
0344   errsfpar[Par_Z0] = 0.;     //r/mygaus->GetParError(1);
0345
0346   // UML Fit in z
0347   //Double_t zlimitUp[2] = { 0., 8. };
0348   //Double_t zlimitDown[2] = { 0., 7. };
0349
0350   gmMinuit->SetFCN(zfcn);
0351   for (int i = 0; i < 2; i++) {
0352     gmMinuit->mnparm(i, par_name[i], sfpar[i], step[i], 0, 0, ierflg);
0353     //gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],zlimitDown[i],zlimitUp[i],ierflg);
0354   }
0355
0356   //std::cout << "fit..." << std::endl;
0357
0358   gmMinuit->Migrad();
0359   // gmMinuit->mncuve();
0360   // gmMinuit->mnmnos();
0361   //gmMinuit->Migrad();
0362   for (int i = 0; i < 2; i++) {
0363     gmMinuit->GetParameter(i, sfpar[i], errsfpar[i]);
0364   }
0365
0366   res_Z_lh[iloop_] = sfpar[Par_Z0];
0367   res_sigmaZ_lh[iloop_] = sfpar[Par_Sigma];
0368   res_Zerr_lh[iloop_] = errsfpar[Par_Z0];
0369   res_sigmaZerr_lh[iloop_] = errsfpar[Par_Sigma];
0370
0371   // get 1 sigma contour of param 0 vs 1
0372   ///gmMinuit->SetErrorDef(1);
0373   ///TGraph *gra1_sig1 = (TGraph*)gMinuit->Contour(60,0,1);
0374   //gra1_sig1->Draw("alp");
0375   //gra1_sig1->SetName("gra1_sig1");
0376   ///gmMinuit->SetErrorDef(4); // two sigma
0377   ///TGraph *gra1_sig2 = (TGraph*)gMinuit->Contour(60,0,1);
0378   //gra1_sig2->Draw("alp");
0379   //gra1_sig2->SetName("gra1_sig2");
0380   // gr12->Draw("alp");
0381
0382   std::cout << "======== End of Maximum LH Fit for Z ========\n" << std::endl;
0383
0384   //_______________ d0 phi fit _________________
0385
0386   tmpNtrks_ = 0;  // reset track counter
0387
0388   TMatrixD xmat(4, 1);
0389   TMatrixDSym Verr(4);
0390   fit(xmat, Verr);  // get initial values
0391   fnthite++;
0392
0393   //Double_t chi2cut = 20;
0394
0395   int fminNtrks = 100;
0396   double fconvergence = 0.9;
0397
0398   while (tmpNtrks_ > fconvergence * zdata.size()) {
0399     // below is a very artificial cut requesting that 50 % of the sample survive
0400     // we hould investigate if there are better criteria than that.
0401     //
0402
0403     // while( cutOnChi2(xmat,chi2cut) > 0.5 * zdata.size() ) {
0404     //   tmpNtrks_ = 0;
0405
0406     fit(xmat, Verr);
0407     fd0cut /= 1.5;
0408     //chi2cut /= 1.5;
0409     if (tmpNtrks_ > fconvergence * zdata.size() && tmpNtrks_ > fminNtrks)
0410       fnthite++;
0411   }
0412   std::cout << " number of iterations: " << fnthite << std::endl;
0413
0414   std::cout << "======== End of d0-phi Fit ========\n" << std::endl;
0415
0416   //__________________ LH Fit for beam __________
0417   /*
0418
0419  for (int i = 0; i<2; i++) {
0420      gmMinuit->GetParameter(i,sfpar[i],errsfpar[i]);
0421  }
0422  sfpar[Par_x0] = xmat(0,0);
0423  sfpar[Par_y0] = xmat(1,0);
0424  sfpar[Par_dxdz]= xmat(2,0);
0425  sfpar[Par_dydz]= xmat(3,0);
0426  sfpar[Par_Sigbeam]=0.0020; // for LHC
0427  //sfpar[Par_Sigbeam]=0.0022; // for CDF
0428  //sfpar[Par_eps] = 3.7e-8;
0429  //sfpar[Par_beta]= 55.;
0430
0431  gmMinuit->mncler();
0432  gmMinuit->SetFCN(dfcn);
0433
0434  for (int i = 2; i<7; i++) {
0435    gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],0,0,ierflg);
0436  }
0437  //gmMinuit->FixParameter(Par_x0);
0438
0439  //gmMinuit->mnparm(9,par_name[9],sfpar[9],step[9],0,0,ierflg);
0440  //gmMinuit->mnparm(10,par_name[10],sfpar[10],step[10],0,0,ierflg);
0441
0442  gmMinuit->Migrad();
0443
0444  // fix x0 and fit again
0445  //for (int i = 2; i<7; i++) {
0446 //   gmMinuit->GetParameter(i,sfpar[i],errsfpar[i]);
0447 // }
0448  //gmMinuit->mncler();
0449  //gmMinuit->SetFCN(dfcn);
0450  //for (int i = 2; i<7; i++) {
0451   // gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],0,0,ierflg);
0452 // }
0453  //gmMinuit->FixParameter(Par_x0);
0454  //gmMinuit->Migrad();
0455
0456  std::cout << "======== End of Maximum LH Fit for Beam ========\n" << std::endl;
0457
0458  //__________________ LH Fit for Product __________
0459
0460  for (int i = 2; i<7; i++) {
0461      gmMinuit->GetParameter(i,sfpar[i],errsfpar[i]);
0462  }
0463
0464  //gmMinuit->GetParameter(5,sfpar[Par_Sigbeam],errsfpar[Par_Sigbeam]);
0465
0466  gmMinuit->mncler();
0467  gmMinuit->SetFCN(cfcn);
0468  for (int i = 0; i<7; i++) {
0469    gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],0,0,ierflg);
0470  }
0471  gmMinuit->Migrad();
0472  gmMinuit->mnhess();
0473  gmMinuit->Migrad();
0474
0475  //   for (int i=0; i<2; i++) {
0476  //  gmMinuit->GetParameter(i,sfpar[i],errsfpar[i]);
0477  // }
0478
0479
0480  gmMinuit->SetErrorDef(4); // two sigma
0481  TGraph *gra1_sig2 = (TGraph*)gMinuit->Contour(60,2,3);
0482  gra1_sig2->Draw("alp");
0483  gra1_sig2->SetName("gra1_sig2");
0484
0485  std::cout << "======== End of Maximum LH Fit Product for Beam ========\n" << std::endl;
0486
0487  //__________________ LH Fit including resolution  __________
0488
0489  for (int i = 0; i<7; i++) {
0490      gmMinuit->GetParameter(i,sfpar[i],errsfpar[i]);
0491  }
0492
0493  // only for pixeless resolution use the generated value:
0494  sfpar[Par_Sigbeam] = 14.e-4;
0495  errsfpar[Par_Sigbeam] = 0.;
0496  //
0497
0498  sfpar[Par_c0] = 0.0100;
0499  sfpar[Par_c1] = 0.0100;
0500
0501  gmMinuit->mncler();
0502  gmMinuit->SetFCN(dfcn2);
0503  for (int i = 0; i<9; i++) {
0504    if (i==7) {
0505      gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],1.e-5,1.e-1,ierflg);
0506    } else if (i==8) {
0507      gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],1.e-5,1.e-1,ierflg);
0508    } else {
0509      gmMinuit->mnparm(i,par_name[i],sfpar[i],step[i],0,0,ierflg);
0510    }
0511  }
0512  // fix beam width
0513  gmMinuit->FixParameter(Par_Sigbeam);
0514
0515  gmMinuit->Migrad();
0516
0517
0518  std::cout << "======== End of Maximum LH Fit Product for Beam ========\n" << std::endl;
0519
0520 */
0521
0522   //outroot->cd();
0523   //r/h_z->Write(); //to draw only markers use "HISTPE1"
0524   //gra1_sig1->Write();
0525
0526   //gra1_sig2->Write();
0527
0528   //app.Run();
0529   // return 0;
0530 }
0531
0532 int main(int argc, char **argv) {
0533   //gROOT->SetBatch(true);
0534
0535   //tnt *t = new tnt("PythiaGenerator/lhc_2008_aa_1001.root");
0536   std::cout << "name: " << argv[1] << std::endl;
0537
0538   ftmp.ResizeTo(4, 1);
0539   ftmp.Zero();
0540
0541   int cut_on_events = 0;
0542   bool run_sequence = false;
0543
0544   if (argc >= 3) {
0545     TString tmpst(argv[2]);
0546     cut_on_events = tmpst.Atoi();
0547     if (cut_on_events == -1) {
0548       run_sequence = true;
0549       std::cout << " run a sequence of fits for 0.5k, 1k, 2k, 5k, 10k" << std::endl;
0550     } else {
0551       std::cout << " maximum number of tracks = " << cut_on_events << std::endl;
0552     }
0553   }
0554
0555   TString in_filename(argv[1]);
0556   TString tmp_string = in_filename;
0557   TString out_filename = tmp_string.Replace(in_filename.Length() - 5, in_filename.Length(), "_results.root");
0558   if (argc == 4) {
0559     out_filename = TString(argv[3]);
0560   }
0561
0562   //static TROOT beamfit("beamfit","Beam fitting");
0563   //static TRint app("app",&argc,argv,NULL,0);
0564
0565   TFile *outroot = new TFile(out_filename, "RECREATE");
0566
0567   NtupleHelper *t = new NtupleHelper(in_filename);
0568   std::cout << " ntuple initialized" << std::endl;
0569
0570   //TH1F *h_z = (TH1F*) gDirectory->Get("z0");
0571   //h_z->SetDirectory(0);
0572
0573   t->Book();
0574   std::cout << " ntuple booked" << std::endl;
0575
0576   if (run_sequence) {
0577     for (int iloop = 0; iloop < Nsequence_; iloop++) {
0578       std::cout << " loop over # " << sequence[iloop] << " tracks." << std::endl;
0579       zdata = t->Loop(sequence[iloop]);
0580       fitAll();
0581       std::cout << " loop has finished." << std::endl;
0582       iloop_++;
0583     }
0584
0585     outroot->cd();
0586     g_X = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_X, zeros, res_d0phi_Xerr);
0587     g_Y = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_Y, zeros, res_d0phi_Yerr);
0588     g_dXdz = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_dXdz, zeros, res_d0phi_dXdzerr);
0589     g_dYdz = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_dYdz, zeros, res_d0phi_dYdzerr);
0590
0591     g_Z = new TGraphErrors(Nsequence_, sequenceD, res_Z, zeros, res_Zerr);
0592     g_sigmaZ = new TGraphErrors(Nsequence_, sequenceD, res_sigmaZ, zeros, res_sigmaZerr);
0593
0594     g_Z_lh = new TGraphErrors(Nsequence_, sequenceD, res_Z_lh, zeros, res_Zerr_lh);
0595     g_sigmaZ_lh = new TGraphErrors(Nsequence_, sequenceD, res_sigmaZ_lh, zeros, res_sigmaZerr_lh);
0596
0597     g_X->SetName("beamX");
0598     g_Y->SetName("beamY");
0599     g_dXdz->SetName("beamdXdz");
0600     g_dYdz->SetName("beamdYdz");
0601     g_Z->SetName("beamZ");
0602     g_sigmaZ->SetName("beamSigmaZ");
0603     g_Z_lh->SetName("beamZ_likelihood");
0604     g_sigmaZ_lh->SetName("beamSigmaZ_likelihood");
0605
0606     //TCanvas *cv_x = new TCanvas("cv_x","beam_X",700,700);
0607     //g_X->Draw("AP");
0608
0609     g_X->Write();
0610     g_Y->Write();
0611     g_dXdz->Write();
0612     g_dYdz->Write();
0613     g_Z->Write();
0614     g_sigmaZ->Write();
0615     g_Z_lh->Write();
0616     g_sigmaZ_lh->Write();
0617
0618     outroot->Close();
0619
0620   } else {
0621     zdata = t->Loop(cut_on_events);
0622     fitAll();
0623     std::cout << " finished loop over events" << std::endl;
0624   }
0625
0626   //cv_x->Print("beam_X.png");
0627
0628   return 0;
0629 }