// $Id$
//
//    File: JEventProcessor_pip_pim.cc
// Created: Wed Jul 27 08:51:36 EDT 2011
// Creator: davidl (on Linux ifarm1102 2.6.18-128.7.1.el5 x86_64)
//

#include <TRACKING/DTrackTimeBased.h>
#include <TRACKING/DMCThrown.h>
#include <PID/DPhoton.h>

#include "JEventProcessor_pip_pim.h"
using namespace jana;


// Routine used to create our JEventProcessor
#include <JANA/JApplication.h>
extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new JEventProcessor_pip_pim());
}
} // "C"


//------------------
// JEventProcessor_pip_pim (Constructor)
//------------------
JEventProcessor_pip_pim::JEventProcessor_pip_pim()
{

}

//------------------
// ~JEventProcessor_pip_pim (Destructor)
//------------------
JEventProcessor_pip_pim::~JEventProcessor_pip_pim()
{

}

//------------------
// init
//------------------
jerror_t JEventProcessor_pip_pim::init(void)
{
	// This is called once at program startup

	pip_mom = new TH1D("pip_mom", "Momentum Magnitude (GeV)", 1200, 0.0, 12.0);
	pip_theta = new TH1D("pip_theta", "Polar angle (degrees)", 200, 0.0, 20.0);
	pip_phi = new TH1D("pip_phi", "Azimuthal angle (degrees)", 720, -180.0, 180.0);

	pim_mom = (TH1D*)pip_mom->Clone("pim_mom");
	pim_theta = (TH1D*)pip_theta->Clone("pim_theta");
	pim_phi = (TH1D*)pip_phi->Clone("pim_phi");

	pip_pim_inv_mass = new TH1D("pip_pim_inv_mass", "#pi^{+}#pi^{-} invariant mass (GeV/c^2)", 1000, 0.0, 2.0);

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_pip_pim::brun(JEventLoop *eventLoop, int runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_pip_pim::evnt(JEventLoop *loop, int eventnumber)
{
	vector<const DMCThrown *> throwns;
	vector<const DTrackTimeBased *> tracks;
	vector<const DPhoton *> photons;

	loop->Get(throwns); // generated particles
	loop->Get(tracks);  // reconstructed charged tracks
	loop->Get(photons); // reconstructed photons
	
	// Find pip and pim, (if any). We keep the number
	// of each type found in case a cut is needed
	// on there being exactly one found.
	//
	// (n.b. We could do exactly the same thing for the
	// "throwns" if the generated particles are needed.)
	int Npip = 0;
	int Npim = 0;
	TLorentzVector pip, pim;
	for(unsigned int i=0; i<tracks.size(); i++){
		const DTrackTimeBased *trk = tracks[i];
		if(fabs(trk->mass() - 0.13957)<0.001){
			if(trk->charge()>0){
				Npip++;
				pip.SetXYZT(trk->px(), trk->py(), trk->pz(), trk->energy());
			}else{
				Npim++;
				pim.SetXYZT(trk->px(), trk->py(), trk->pz(), trk->energy());
			}
		}
	}
	
	// Fill in pip histos
	if(Npip==1){

		pip_mom->Fill(pip.P());
		pip_theta->Fill(pip.Theta()*TMath::RadToDeg());
		pip_phi->Fill(pip.Phi()*TMath::RadToDeg());

	}

	// Fill in pim histos
	if(Npim==1){

		pim_mom->Fill(pim.P());
		pim_theta->Fill(pim.Theta()*TMath::RadToDeg());
		pim_phi->Fill(pim.Phi()*TMath::RadToDeg());

	}

	// Fill in invariant mass histo
	if(Npip==1 && Npim==1){
		pip_pim_inv_mass->Fill( (pip+pim).M() );
	}

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_pip_pim::erun(void)
{
	// This is called whenever the run number changes, before it is
	// changed to give you a chance to clean up before processing
	// events from the next run number.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t JEventProcessor_pip_pim::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}