//--------------------------------------------------------------------------
// File and Version Information: MyAnaModule.cc version 1.0
//
// Description: A simple analysis module which runs over tracks, jets,
//						electrons, obsp particles, and dumps the info into a
//						simple root tree.
//
// Environment: Software developed for CDF.
//
// Created: 03/23/01 Richard E. Hughes
//                   The Ohio State University
//                   CDF Group
//
//------------------------------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "MyAnaModule.hh"

//-------------
// C Headers --
//-------------
#include <assert.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream>

//---------------
// C++ Headers --
//---------------
   
   
//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "evt/Event.hh"
#include "AbsEnv/AbsEnv.hh"
#include "Trybos/TRY_Record_Iter_Same.hh"
#include "Trybos/TRY_Bank_Type.hh"
#include "Framework/AbsParmGeneral.hh"

#include "TFile.h"
#include "TRandom.h"
#include "TTree.h"
#include "TH1.h"
#include "TNtuple.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TStopwatch.h"


#include "RootUtils/Utils/TCdfRint.hh"
#include "RootUtils/Utils/TCdfRoot.hh"

#include "MyAnaModule.hh"

#include "Edm/EventRecord.hh"
#include "Edm/GenericConstHandle.hh"
#include "Edm/ConstHandle.hh"
#include "Edm/Handle.hh"

#include "SimulationObjects/OBSP_StorableBank.hh"
#include "SimulationObjects/OBSV_StorableBank.hh"
#include "ParticleDB/CdfParticleDatabase.hh"
//
#include "TrackingUtils/GeneratorLevel/CDFPDGParticle.hh"
#include "TrackingUtils/GeneratorLevel/PDG1994Particle.hh"
#include "TrackingUtils/GeneratorLevel/GenPrimaryVertex.hh"
#include "TrackingUtils/GeneratorLevel/GenPrimaryVertexSet.hh"
#include "TrackingUtils/GeneratorLevel/ParentedParticle.hh"
#include "TrackingObjects/Tracks/SmearedFourVectorTrack.hh"
#include "TrackingObjects/Tracks/track_cut.hh"
#include "TrackingUtils/Trybos/TRYGenPrimaryVertexSet.hh"

#include "CLHEP/config/CLHEP.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/Vector.h"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "Trajectory/Helix.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"


//kludge
#ifdef USE_CDFEDM2
#define USE_BFIELD_MAP
#endif

#ifdef USE_BFIELD_MAP
#include "GeometryBase/Bfield.hh"
#include "inc/jobsta.h"
extern "C" void gufld_ (const float* x,float* h);
extern "C" void bfinit_();
extern "C" void* jobstr_address_(void);
#endif // USE_BFIELD_MAP

//#define CT_DEDX

using std::ostream;
using std::cout;
using std::endl;

TFile*         File;
TTree*         myanaTree;
TBranch*       myBranch;
TH1F*          h_ptCOT;
TH1F*          h_ptOBSP;
TNtuple*       n_COT;
TNtuple*       n_OBSP;
//  TRYGenPrimaryVertexSet myGenSet;

  GenPrimaryVertexSet *theParticleSet;

Int_t 	num_resid;
Int_t		trknum_resid[10000];
Int_t		obsp_resid[10000];
Int_t		sl_resid[10000];
Int_t		cell_resid[10000];
Int_t		wire_resid[10000];
Float_t 	d_resid[10000];

Int_t 	nTrack_t;
Int_t 	nRun_t;
Int_t 	nEvt_t;
Int_t		trknum_t[500];
Float_t 	cu_t[500];
Float_t 	cot_t[500];
Float_t 	phi_t[500];
Float_t 	d0_t[500];
Float_t 	z0_t[500];
Float_t 	cuErr_t[500];
Float_t 	cotErr_t[500];
Float_t 	phiErr_t[500];
Float_t 	d0Err_t[500];
Float_t 	z0Err_t[500];
Int_t 	alg_t[500];

Int_t 	nObsp_t;
Float_t 	pvx_t;
Float_t 	pvy_t;
Float_t 	pvz_t;
Float_t 	obsp_vx_t[500];
Float_t 	obsp_vy_t[500];
Float_t 	obsp_vz_t[500];
Int_t 	obsp_type_t[500];
Float_t 	obsp_cu_t[500];
Float_t 	obsp_cot_t[500];
Float_t 	obsp_phi_t[500];
Float_t 	obsp_d0_t[500];
Float_t 	obsp_z0_t[500];
Int_t		obsp_id_t[500];


//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

//----------------
// Constructors --
//----------------
MyAnaModule::MyAnaModule( 
    const char* const theName, 
    const char* const theDescription )
    : RootHistModule( theName, theDescription )
{
}
//--------------
// Destructor --
//--------------
MyAnaModule::~MyAnaModule( )
{
}

//--------------
// Operations --
//--------------
AppResult	 MyAnaModule::beginJob(AbsEvent* event )
{
//
   cout << " MyAnaModule::beginJob " << endl;
//
// Let's define some simple histograms
	h_ptOBSP =  new TH1F("ptOBPS","Pt of Obsp Particles",100,0.,20.);
	HistogramList()->Add(h_ptOBSP);
	h_ptCOT =  new TH1F("ptCOT","Pt of COT Trackss",100,0.,20.);
   HistogramList()->Add(h_ptCOT);
//
// Let's define some simple ntuples
	n_COT = new TNtuple("COTNtuple","COT Tracks","pt:phi:cot:d0:z0");
   HistogramList()->Add(n_COT);
	n_OBSP = new TNtuple("OBSPNtuple","OBSP Particles","pt:phi:cot:d0:z0:ID");
   HistogramList()->Add(n_OBSP);
//
// Define a more complicated object: a root tree (a column-wise ntuple)
	myanaTree = new TTree("myana","Tracking Analysis Tree");
//
// The run and event number of this event	
	myanaTree->Branch("nRun_t",&nRun_t,"nRun_t/I");
	myanaTree->Branch("nEvt_t",&nEvt_t,"nEvt_t/I");
//
// The primary vertex for this event
	myanaTree->Branch("pvx_t",&pvx_t,"pvx_t/F");
	myanaTree->Branch("pvy_t",&pvy_t,"pvy_t/F");
	myanaTree->Branch("pvz_t",&pvz_t,"pvz_t/F");
//
// An array of tracking info for this event
	myanaTree->Branch("nTrack_t",&nTrack_t,"nTrack_t/I");
	myanaTree->Branch("trknum_t",trknum_t,"trknum_t[nTrack_t]/F");
	myanaTree->Branch("cu_t",cu_t,"cu_t[nTrack_t]/F");
	myanaTree->Branch("cot_t",cot_t,"cot_t[nTrack_t]/F");
	myanaTree->Branch("phi_t",phi_t,"phi_t[nTrack_t]/F");
	myanaTree->Branch("d0_t",d0_t,"d0_t[nTrack_t]/F");
	myanaTree->Branch("z0_t",z0_t,"z0_t[nTrack_t]/F");
	myanaTree->Branch("cuErr_t",cuErr_t,"cuErr_t[nTrack_t]/F");
	myanaTree->Branch("cotErr_t",cotErr_t,"cotErr_t[nTrack_t]/F");
	myanaTree->Branch("phiErr_t",phiErr_t,"phiErr_t[nTrack_t]/F");
	myanaTree->Branch("d0Err_t",d0Err_t,"d0Err_t[nTrack_t]/F");
	myanaTree->Branch("z0Err_t",z0Err_t,"z0Err_t[nTrack_t]/F");
	myanaTree->Branch("alg_t",alg_t,"alg_t[nTrack_t]/I");
//
// An array of OBSP info for this event
	myanaTree->Branch("nObsp_t",&nObsp_t,"nObsp_t/I");
	myanaTree->Branch("obsp_cu_t",obsp_cu_t,"obsp_cu_t[nObsp_t]/F");
	myanaTree->Branch("obsp_cot_t",obsp_cot_t,"obsp_cot_t[nObsp_t]/F");
	myanaTree->Branch("obsp_phi_t",obsp_phi_t,"obsp_phi_t[nObsp_t]/F");
	myanaTree->Branch("obsp_d0_t",obsp_d0_t,"obsp_d0_t[nObsp_t]/F");
	myanaTree->Branch("obsp_z0_t",obsp_z0_t,"obsp_z0_t[nObsp_t]/F");
	myanaTree->Branch("obsp_id_t",obsp_id_t,"obsp_id_t[nObsp_t]/I");
	myanaTree->Branch("obsp_type_t",obsp_type_t,"obsp_type_t[nObsp_t]/I");
	myanaTree->Branch("obsp_vx_t",obsp_vx_t,"obsp_vx_t[nObsp_t]/F");
	myanaTree->Branch("obsp_vy_t",obsp_vy_t,"obsp_vy_t[nObsp_t]/F");
	myanaTree->Branch("obsp_vz_t",obsp_vz_t,"obsp_vz_t[nObsp_t]/F");
//
// Add the defined tree to the list of kept histograms
	HistogramList()->Add(myanaTree);

//

    return AppResult::OK;
}
//
// The event routine
AppResult MyAnaModule::event( AbsEvent* anEvent )
{
//
// Get the run and event number
	int _nrun = AbsEnv::instance()->runNumber();
   int _nevt = AbsEnv::instance()->trigNumber();
	printf("new event run number %d event number % d\n");
//
// Initialise the track set from MC truth.
	TRYGenPrimaryVertexSet *myEventStorage = new TRYGenPrimaryVertexSet(*anEvent);

//By default, this event does not pass the filter.
//
// Try OBSP first
   bool useHEPG=false;
   if (myEventStorage->canInitializeFromObs()) {
	   myEventStorage->initializeFromObs();
	   useHEPG=false;
   } else { // Either or both banks don't exist -- use HEPG
   	cerr << name() << ": either or both of OBSP and OBSV banks are missing --"
      	  << endl << "using HEPG instead." << endl;
   	useHEPG=true;
   }
//
// If OBSP was not there
   if (useHEPG) {     
		myEventStorage->initializeFromHepevt();
   }
//  
// Iterate over vertices
	nObsp_t = 0;
	Hep3Vector obs3vector[500];
   HepPoint3D obsProdVert[500];
    
	GenPrimaryVertexSet::ConstIterator currentVertex;
	for (currentVertex=myEventStorage->begin();
       currentVertex!=myEventStorage->end();
       currentVertex++) {
		cout << "New Primary Vertex at " <<  (*currentVertex).getPosition() << endl;
		pvx_t = (*currentVertex).getPosition().x();
		pvy_t = (*currentVertex).getPosition().y();
		pvz_t = (*currentVertex).getPosition().z();
      		
// Iterate over particles in each vertex
      GenPrimaryVertex::ConstParticleIterator currentParticle;
      for (currentParticle=(*currentVertex).begin();
         	currentParticle!=(*currentVertex).end();
         	currentParticle++) {
//
// Here we iterate over all known, charged, stable
// particles.
			if ((*currentParticle).isKnown() // sanity-check
          	&& !(*currentParticle).hasDecayed() // final-state
         	 && (*currentParticle).iCharge()) // PYTHIA-style integer charge (=q*3)
			{
	

// Make a helix from it

// This slightly circuitous route to get the momentum uses a
// () cast in the definition of HepLorentzVector to convert
// from the 4-momentum that ParentedParticle::momentum()
// returns
				Hep3Vector p3momentum((*currentParticle).momentum());
//
// Create a `smeared track', although in this case the
// smearing model is the null pointer
				SmearedFourVectorTrack gTrack(p3momentum,
                                        (*currentParticle).productionVertex(),
                                        (*currentParticle).charge(),
                                        NULL);
// Create the helix
				Helix gHelix=gTrack.getHelixFit().getMeasureable();
            float pt = fabs(gTrack.pt());
//
// Fill the histogram
            h_ptOBSP->Fill(pt);
//
// Keep tracks with Pt>1.0
				if (pt > 1.0) {
//
// Get the decay vertex
				const Hep3Vector vertparticle = (*currentParticle).productionVertex();
// Print out some information
//                cout <<"b daugh " << bdaughter 
//          		      << " k0s " << k0sdaughter 
//                     << " type " << (*currentParticle).type()
//                     << " x,y " 
//          			      << vertparticle.getX() << 
//                        " " << vertparticle.getY() << 
//                        " " << vertparticle.getZ() << endl;

// Now you can get all the information you like from the
// helix. See
// TrackingUtils/Trajectory/Helix.hh for
// details. 
               double phi = gHelix.getPhi0();
               double gencot = gHelix.getCotTheta();
               double gencrv = gHelix.getCurvature();
               double genZ0 = gHelix.getZ0();
               double genD0 = gHelix.getD0();
               int numCh = (*currentParticle).numChildren();
               int myid = (*currentParticle).id();
               double mass = (*currentParticle).mass();
//
// Fill the ntuple
					float nt[6] = {pt,phi,gencot,genD0,genZ0,
               							(*currentParticle).type()};
					n_OBSP->Fill(nt);
//
// Fill the tree vars
               obsp_cu_t[nObsp_t] = gencrv;
               obsp_cot_t[nObsp_t] = gencot;
               obsp_phi_t[nObsp_t] = phi;
               obsp_d0_t[nObsp_t] = genD0;
               obsp_z0_t[nObsp_t] = genZ0;
               obsp_id_t[nObsp_t] = myid;
               obsp_vx_t[nObsp_t] = (float)vertparticle.getX();
               obsp_vy_t[nObsp_t] = (float)vertparticle.getY();
               obsp_vz_t[nObsp_t] = (float)vertparticle.getZ();
               obsp_type_t[nObsp_t] = (*currentParticle).type();;

               nObsp_t++;

				} // pt > cut
			} // endif known, charged, stable.
        
		} // loop over particles
      
	} // loop over vertices
//
// Now get ready to look at COT track info
	nTrack_t = 0;
//
   CdfTrackView_h hView; // This is the handle for the "view."
	printf("Looping over tracks \n\n");
   if (CdfTrackView::defTracks(hView) == CdfTrackView::OK) { 

// The view is now filled with the default tracks, so extract contents.
		const CdfTrackView::CollType & tracks = hView->contents();

// Now loop over the tracks, doing a double-dereference to get at each.
      for (CdfTrackView::const_iterator it = tracks.begin();
      it != tracks.end(); ++it) {

         const CdfTrack & trk = **it;
         CdfTrack_lnk cdfTrack=(*it);

//Extract the helix parameters (or equivalent) for each track
         float crv = trk.curvature();
         float phi0 = trk.phi0();
         float d0 = trk.d0();
         float z0 = trk.z0();
         float pt = trk.pt();
         float ctg = trk.lambda();
         int alg = trk.algorithm().value();
//
// Fill the histogram
         h_ptCOT->Fill(pt);
//
// Get the error matrix
//      	HepSymMatrix E=trk.getHelixFit().getErrorMatrix();
// Is this the same as the previous?
      	HepSymMatrix E=(*cdfTrack).getHelixFit().getErrorMatrix();
         float ctgErr = sqrt(E[0][0]);
         float crvErr = sqrt(E[1][1]);
         float z0Err = sqrt(E[2][2]);
         float d0Err = sqrt(E[3][3]);
         float phi0Err = sqrt(E[4][4]);
//
// Do the following for high pt tracks which are of one tracking algorithm         
//	      if ((fabs(pt) > 1.0) && (alg == CdfTrack::OutsideInAlg)) { 
	      if ((fabs(pt) > 1.0) && (alg == CdfTrack::CotStandAloneAlg)) { 
//
// Fill the ntuple
				float nt[5] = {pt,phi0,ctg,d0,z0};
				n_COT->Fill(nt);
//
// Fill the tree vars
			   trknum_t[nTrack_t] = nTrack_t;
			   cu_t[nTrack_t] = crv;
			   cot_t[nTrack_t] = ctg;
			   phi_t[nTrack_t] = phi0;
			   d0_t[nTrack_t] = d0;
			   z0_t[nTrack_t] = z0;
			   cuErr_t[nTrack_t] = crvErr;
			   cotErr_t[nTrack_t] = ctgErr;
			   phiErr_t[nTrack_t] = phi0Err;
			   d0Err_t[nTrack_t] = d0Err;
			   z0Err_t[nTrack_t] = z0Err;
   		   nTrack_t++;

			} 		// end check on pt/algorithm

		}	// end loop over tracks
	}
 
   myanaTree->Fill();
//
// We are done
	return AppResult::OK;
}