#include <iostream>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>
#include "hddm_fitter.h"
using namespace std;
using namespace TMath;

const int debug_level = 2;

int main()
{


  TFile *file = new TFile("resids.root", "RECREATE");
  Int_t bins = 50;
  Double_t xlow = -0.1, xhigh = 0.1;
  TH1F *h1 = new TH1F("residsH1", "residuals", bins, xlow, xhigh);
  TH1F *h1Resids2 = new TH1F("resids2", "squared residuals", 100, 0.0, 0.001);
  TH1F *h1ChiSq = new TH1F("chisq", "chisq from fit", 100, 0.0, 1.0);
  TH1F *h1Prob = new TH1F("prob", "tracking chi-squared probability", 100, 0.0, 1.0);
  char hddmFileName[128] = "fitter.hddm";

  int count = 0;
  fitter_HDDM_t *thisInputEvent = 0;
  fitter_iostream_t *thisInputFile = 0;
  if (! (thisInputFile = open_fitter_HDDM(hddmFileName))) {
    printf("Error - could not open input file\n");
    exit(1);
  }

  int number = -1;
  float resid = -1.0;
  double chisq;
  int nCDC;
  double residFit = 4.62908e-03;
  double residFit2 = residFit*residFit;
  double prob;
  while (thisInputEvent = read_fitter_HDDM(thisInputFile)) {
    count++;
    if (debug_level > 1) printf("count = %d\n", count);
    number = thisInputEvent->events->in[0].number;
    if (debug_level > 1) printf("event number = %d\n", number);
    int ntracks = thisInputEvent->events->in[0].tracks->mult;
    if (debug_level > 1) cout << "ntracks = " << ntracks << endl;
    if (ntracks == 1) {
      fitter_Track_t* trackPtr = &(thisInputEvent->events->in[0].tracks->in[0]);
      if (trackPtr->fitStatus == 0) {
	chisq = trackPtr->chisq;
	nCDC = trackPtr->nCDC;
	int ntrajs = trackPtr->trajectorys->mult;
	if (debug_level > 1) cout << "ntrajs = " << ntrajs << endl;
	h1ChiSq->Fill(chisq);
	prob = Prob(chisq/residFit2, nCDC);
	h1Prob->Fill(prob);
	if (prob < 0.01) {cout << "bad event=" << number << " prob=" << prob << " chisq=" << chisq << " nCDC=" << nCDC << " count=" << count << endl;}
	fitter_Trajectory_t* trajPtr = &(trackPtr->trajectorys->in[3]); 
	int nresids = trajPtr->residuals->mult;
	if (debug_level > 1) cout << "nresids = " << nresids << endl;
	for (int i = 0; i < nresids; i++) {
	  resid = trajPtr->residuals->in[i].value;
	  if (debug_level > 1) cout << i << " resid = " << resid << endl;
	  h1->Fill(resid);
	  h1Resids2->Fill(resid*resid);
	}
      }
    }
    flush_fitter_HDDM(thisInputEvent,0);
  }
  close_fitter_HDDM(thisInputFile);

  file->Write();

  return 0;

}