// MyProcessor.cc
//

#include <iostream>
#include "FRCTProcessor.h"
#include "residFDCAnode.h"

//------------------------------------------------------------------
// FRCTProcessor 
//------------------------------------------------------------------
FRCTProcessor::FRCTProcessor()
{
  cout << "FRCTProcessor constructor called\n";
}

//------------------------------------------------------------------
// ~FRCTProcessor 
//------------------------------------------------------------------
FRCTProcessor::~FRCTProcessor()
{
  cout << "FRCTProcessor destructor called\n";
}

//------------------------------------------------------------------
// init 
//------------------------------------------------------------------
jerror_t FRCTProcessor::init(void)
{
  rootfile = new TFile("fdcCathodeResidTest.root", "RECREATE");
  nTest = new TNtuple("nTest", "test results", "event:resid");
  tTest = new TTree("fdcCathodeResidTestTree", "to store the data");
  tTest->Branch("event", &ei.run, "run/I:event/I:x0/D:y0/D:z0/D");
  tTest->Branch("nResids", &nResids, "nResids/I");
  tTest->Branch("wire", wire, "wire[nResids]/I");
  tTest->Branch("layer", layer, "layer[nResids]/I");
  tTest->Branch("resids", resids, "resids[nResids]/D");

  return NOERROR;
}

//------------------------------------------------------------------
// fini
//------------------------------------------------------------------
jerror_t FRCTProcessor::fini(void)
{
  cout << "fini called\n";
  rootfile->Write();
  return NOERROR;
}

//------------------------------------------------------------------
// brun 
//------------------------------------------------------------------
jerror_t FRCTProcessor::brun(JEventLoop *eventLoop, int runnumber)
{
  return NOERROR;
}

//------------------------------------------------------------------
// erun 
//------------------------------------------------------------------
jerror_t FRCTProcessor::erun()
{
  return NOERROR;
}

//------------------------------------------------------------------
// evnt 
//------------------------------------------------------------------

jerror_t FRCTProcessor::evnt(JEventLoop *eventLoop, int eventnumber)
{

  cout << "event number " << eventnumber << endl;

  ei.run = 4321;
  ei.event = eventnumber;

  int debug_level = 3;
  // get the pseudopoints, print out info about them
  vector<const DFDCPseudo*>pseudopoints;
  eventLoop->Get(pseudopoints);
  int size_fdc =  pseudopoints.size();
  cout << "number of pseudopoints " << size_fdc << endl;
  const DFDCPseudo* thisPseudo;
  if (size_fdc > 0) {
    for (int i = 0; i < size_fdc; i++) {
      thisPseudo = pseudopoints[i];
      const DFDCWire *wire = thisPseudo->wire;
      double xw = wire->origin.X();
      double yw = wire->origin.Y();
      double zw = wire->origin.Z();
      double thetaw = acos(wire->udir.z());
      double phiw = atan2(wire->udir.y(), wire->udir.x());
      cout << setprecision(14) << "wire info: x = " << xw << " y = " << yw << " z = " << zw << " theta " << thetaw << " phi " << phiw << endl;
    }
  }

  vector<const DMCThrown*> mcthrowns;
  eventLoop->Get(mcthrowns);
  int ierror = 1;
  if (mcthrowns.size() != 1) throw ierror;
  DVector3 pos = mcthrowns[0]->position();
  cout << "X=" << pos.X() << "Y=" << pos.Y() << "Z=" << pos.Z() << endl;

  // create a trajectory, swim a straight line print results
  MyTrajectory lineTraj;
  HepVector startingVector(3);
  startingVector(1) = pos.X();
  startingVector(2) = pos.Y();
  startingVector(3) = pos.Z();
  double theta = 0.0;
  double phi = 0.0;
  lineTraj.swim(startingVector, theta, phi);
  //lineTraj.print();

  ei.vertex[0] = startingVector(1);
  ei.vertex[1] = startingVector(2);
  ei.vertex[2] = startingVector(3);

  // find the residuals from the FDC anodes
  residFDCAnode rFDCA(&pseudopoints, &lineTraj, debug_level);
  rFDCA.setInnerResidFrac(1.0);
  rFDCA.calcResids();
  vector<double> res;
  rFDCA.getResids(res);
  nResids = res.size();
  cout << "number of resids = " << nResids << endl;
  for (uint i = 0; i < res.size(); i++) {
    cout << i << "," << res[i] << "/";
    nTest->Fill(eventnumber, res[i]);
    wire[i] = pseudopoints[i]->wire->wire;
    layer[i] = pseudopoints[i]->wire->layer;
    resids[i] = res[i];
  }
  cout << endl;
  tTest->Fill();

  return NOERROR;
 
}