{

  #include "addToChain.C"

  gStyle->SetPalette(1); // choose sensible rainbow-style color palette
  gStyle->SetOptStat("nemruoi");

  gROOT->ProcessLine(".L getPads.C");
  gROOT->ProcessLine(".L clearCutStrings.C");
  gROOT->ProcessLine(".L makeCanvas.C");
  gROOT->ProcessLine(".L makeCanvasHist1d.C");

  double degPerRad = 57.295780;

  TChain* chainPtr = new TChain("fitter_tree");
  addToChain(chainPtr);
  int xsize = 1280;
  int ysize = 1024;

  char title[128];
  char drawString[50][128];
  char cutString[50][128];

  int canvasId = 0;

  sprintf(title, "multiplicities");
  int npad = 4, nxpad = 2, nypad = 2;
  clearCutStrings(cutString);
  sprintf(drawString[0], "nCDC >> hNCDC(41, -0.5, 40.5)");
  sprintf(drawString[1], "nFDC >> hNFDC(41, -0.5, 40.5)");
  sprintf(drawString[2], "nCDC + nFDC >> hNtot(41, -0.5, 40.5)");
  sprintf(drawString[3], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "scalers");
  int npad = 4, nxpad = 2, nypad = 2;
  clearCutStrings(cutString);
  sprintf(drawString[0], "event");
  sprintf(drawString[1], "progress >> hProgress(10, -0.5, 9.5)");
  sprintf(drawString[2], "status >> hStatus(10, -4.5, 5.5)");
  sprintf(drawString[3], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "momenta, true and fit");
  clearCutStrings(cutString);
  int npad = 4, nxpad = 2, nypad = 2;
  sprintf(drawString[0], "ptinvTrue >> hPtinvTrue(100, 0.0, 5.0)");
  sprintf(drawString[1], "1.0/ptinvTrue/sin(thetaTrue) >> hPtotTrue(100, 0.0, 5.0)");
  sprintf(drawString[2], "ptinvFit >> hPtinvFit(100, 0.0, 5.0)");
  sprintf(cutString[2], "status == 0");
  sprintf(drawString[3], "1.0/ptinvFit/sin(thetaFit) >> hPtotFit(100, 0.0, 5.0)");
  sprintf(cutString[3], "status == 0");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "angles, true and fit");
  clearCutStrings(cutString);
  int npad = 4, nxpad = 2, nypad = 2;
  sprintf(drawString[0], "thetaTrue*57.295780 >> hThetaTrue(100, -10.0, 190.0)");
  sprintf(drawString[1], "phiTrue*57.295780 >> hPhiTrue(100, -250.0, 250.0)");
  sprintf(drawString[2], "thetaFit*57.295780 >> hThetaFit(100, -10.0, 190.0)");
  sprintf(drawString[3], "phiFit*57.295780 >> hPhiFit(100, -250.0, 250.0)");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "launch point, true and fit");
  sprintf(drawString[0], "xp0True >> hXp0True(100, -25.0, 25.0)");
  sprintf(drawString[1], "z0True >> hZ0True(100, 55.0, 75.0)");
  sprintf(drawString[2], "xp0Fit >> hXp0Fit(100, -5.0, 5.0)");
  sprintf(drawString[3], "z0Fit >> hZ0Fit(100, 55.0, 75.0)");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  npad = 6; nxpad = 3; nypad = 2;

  for (int ip = 0; ip < 5; ip++) {
      sprintf(cutString[ip], "status == 0");
  }

  sprintf(title, "all variables, true, good fit status");
  // OK here means that the fit was OK
  sprintf(drawString[0], "ptinvTrue >> hPtinvTrueOK(100, 0.0, 5.0)");
  sprintf(drawString[1], "thetaTrue*57.295780 >> hThetaTrueOK(100, -10.0, 190.0)");
  sprintf(drawString[2], "phiTrue*57.295780 >> hPhiTrueOK(100, -250.0, 250.0)");
  sprintf(drawString[3], "xp0True >> hXp0TrueOK(100, -25.0, 25.0)");
  sprintf(drawString[4], "z0True >> hZ0TrueOK(100, 55.0, 75.0)");
  sprintf(drawString[5], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "fit - true");
  sprintf(drawString[0], "xp0Fit - xp0True >> hXp0Diff(100, -5, 5)");
  sprintf(drawString[1], "z0Fit - z0True >> hZ0Diff(100, -25, 25)");
  sprintf(drawString[2], "thetaFit - thetaTrue >> hThetaDiff(100, -0.25, 0.25)");
  sprintf(drawString[3], "phiFit - phiTrue >> hPhiDiff(100, -0.25, 0.25)");
  sprintf(drawString[4], "ptinvFit - ptinvTrue >> hPtinvDiff(100, -0.25, 0.25)");
  sprintf(drawString[5], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "iterations, chisq, probabilities");
  clearCutStrings(cutString);
  int npad = 6, nxpad = 3, nypad = 2;
  sprintf(drawString[0], "iters >> hIters(50, -0.5, 49.5)");
  sprintf(cutString[0], "status == 0");
  sprintf(drawString[1], "chisq/(nCDC + nFDC - 5) >> hChisq(100, 0.0, 5)");
  sprintf(cutString[1], "status == 0");
  sprintf(drawString[2], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProb(100, 0.0, 1.0)");
  sprintf(cutString[2], "status == 0");
  sprintf(drawString[3], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbZ(100, 0.0, 0.1)");
  sprintf(cutString[3], "status == 0");
  sprintf(drawString[4], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbZZ(100, 0.0, 0.01)");
  sprintf(cutString[4], "status == 0");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "probabilities, CDC only, FDC only");
  clearCutStrings(cutString);
  int npad = 6, nxpad = 3, nypad = 2;
  sprintf(drawString[0], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbCDC(100, 0.0, 1.0)");
  sprintf(cutString[0], "status == 0 && nFDC == 0");
  sprintf(drawString[1], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbCDCZ(100, 0.0, 0.1)");
  sprintf(cutString[1], "status == 0 && nFDC == 0");
  sprintf(drawString[2], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbCDCZZ(100, 0.0, 0.01)");
  sprintf(cutString[2], "status == 0 && nFDC == 0");
  sprintf(drawString[3], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbFDC(100, 0.0, 1.0)");
  sprintf(cutString[3], "status == 0 && nCDC == 0");
  sprintf(drawString[4], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbFDCZ(100, 0.0, 0.1)");
  sprintf(cutString[4], "status == 0 && nCDC == 0");
  sprintf(drawString[5], "TMath::Prob(chisq, nCDC + nFDC - 5) >> hProbFDCZZ(100, 0.0, 0.01)");
  sprintf(cutString[5], "status == 0 && nCDC == 0");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  sprintf(title, "residuals");
  clearCutStrings(cutString);
  int npad = 6, nxpad = 3, nypad = 2;
  sprintf(drawString[0], "resids >> hResids(100, -.1, 0.1)");
  sprintf(cutString[0], "status == 0");
  sprintf(drawString[1], "resids >> hResidsCDC(100, -0.1, 0.1)");
  sprintf(cutString[1], "status == 0 && nFDC == 0");
  sprintf(drawString[2], "resids >> hResidsCDCGood(100, -0.1, 0.1)");
  sprintf(cutString[2], "status == 0 && nFDC == 0 && TMath::Prob(chisq, nCDC + nFDC - 5) > 0.01");
  sprintf(drawString[3], "resids >> hResidsFDC(100, -0.01, 0.19)");
  sprintf(cutString[3], "status == 0 && nCDC == 0");
  sprintf(drawString[4], "resids >> hResidsFDCGood(300, -.01, 0.19)");
  sprintf(cutString[4], "status == 0 && nCDC == 0 && TMath::Prob(chisq, nCDC + nFDC - 5) > 0.01");
  sprintf(drawString[5], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  clearCutStrings(cutString);
  int npad = 6, nxpad = 3, nypad = 2;
  sprintf(title, "pulls and param chisq");
  for (int ip = 0; ip < 5; ip++) {
      sprintf(cutString[ip], "status == 0 && TMath::Prob(chisq, nFDC + nCDC - 5) > 0.01");
  }
  sprintf(drawString[0], "(xp0Fit - xp0True)/sigmaXp0 >> hPullXp0(100, -25., 25.)");
  sprintf(drawString[1], "(z0Fit - z0True)/sigmaZ0 >> hPullZ0(100, -25., 25.)");
  sprintf(drawString[2], "(thetaFit - thetaTrue)/sigmaTheta >> hPullTheta(100, -25., 25.)");
  sprintf(drawString[3], "(phiFit - phiTrue)/sigmaPhi >> hPullPhi(100, -25., 25.)");
  sprintf(drawString[4], "(ptinvFit - ptinvTrue)/sigmaPtinv >> hPullPtinv(100, -25., 25.)");
  sprintf(drawString[5], "chisqParam >> hChisqParam(100, 0., 10000.)");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  clearCutStrings(cutString);
  int npad = 6, nxpad = 3, nypad = 2;

  for (int ip = 0; ip < 5; ip++) {
      sprintf(cutString[ip], "status == 0");
  }

  TH1F hEffXp0 = *hXp0TrueOK / *hXp0True;
  TH1F hEffZ0 = *hZ0TrueOK / *hZ0True;
  TH1F hEffTheta = *hThetaTrueOK / *hThetaTrue;
  TH1F hEffPhi = *hPhiTrueOK / *hPhiTrue;
  TH1F hEffPtinv = *hPtinvTrueOK / *hPtinvTrue;

  TH1F* hist1dPtrAry[50]; // array of pointers to 1d histograms
  int npad = 6, nxpad = 3, nypad = 2;

  sprintf(title, "efficiencies");
  hist1dPtrAry[0] = &hEffXp0;
  hist1dPtrAry[1] = &hEffZ0;
  hist1dPtrAry[2] = &hEffTheta;
  hist1dPtrAry[3] = &hEffPhi;
  hist1dPtrAry[4] = &hEffPtinv;
  hist1dPtrAry[5] = NULL;
  makeCanvasHist1d(++canvasId, title, xsize, ysize, npad, nxpad, nypad, hist1dPtrAry);

  clearCutStrings(cutString);
  int npad = 4, nxpad = 2, nypad = 2;

  sprintf(drawString[0], "1.0/ptinvTrue/sin(thetaTrue):thetaTrue*57.295 >> hPThetaAll(100, -10.0, 190.0, 50, 0.0, 5.0");
  sprintf(cutString[0], "");
  sprintf(drawString[1], "1.0/ptinvTrue/sin(thetaTrue):thetaTrue*57.295 >> hPThetaFit(100, -10.0, 190.0, 50, 0.0, 5.0");
  sprintf(cutString[1], "status == 0");
  sprintf(drawString[2], "1.0/ptinvTrue/sin(thetaTrue):thetaTrue*57.295 >> hPThetaFitBad(100, -10.0, 190.0, 50, 0.0, 5.0");
  sprintf(cutString[2], "status == 0 && TMath::Prob(chisq, nCDC + nFDC - 5) < 0.01");
  sprintf(drawString[3], "NULL");
  makeCanvas(++canvasId, title, xsize, ysize, npad, nxpad, nypad, drawString, cutString);

  gStyle->SetOptStat(kFALSE);

  sprintf(title, "fit efficiency, p vs. theta");
  TH2F hEff = *hPThetaFit / *hPThetaAll;
  TCanvas* canvasPtr = new TCanvas("canvasFitEfficiency", title, xsize, ysize);
  hEff->Draw("colZ");
  char pngFile[128];
  sprintf(pngFile, "standard%d.png", ++canvasId);
  canvasPtr->Print(pngFile,"png");

  sprintf(title, "good fit ineficiency, p vs. theta");
  TH2F hFitErrorFraction = *hPThetaFitBad / *hPThetaFit;
  TCanvas* canvasPtr2 = new TCanvas("canvasGoodFit", title, xsize, ysize);
  hFitErrorFraction->Draw("colZ");
  sprintf(pngFile, "standard%d.png", ++canvasId);
  canvasPtr2->Print(pngFile,"png");

}