// $Id$
//
//    File: JEventProcessor_genmumu.cc
// Created: Mon Jan  4 08:47:50 EST 2016
// Creator: davidl (on Linux gluon02.jlab.org 2.6.32-358.18.1.el6.x86_64 x86_64)
//

#include <stdint.h>

#include <iostream>
using namespace std;

#include "JEventProcessor_genmumu.h"
using namespace jana;

#include <PID/DBeamPhoton.h>
#include <TRACKING/DMCThrown.h>


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


//------------------
// JEventProcessor_genmumu (Constructor)
//------------------
JEventProcessor_genmumu::JEventProcessor_genmumu()
{

}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_genmumu::init(void)
{
	japp->RootWriteLock();

	h_Ebeam = new TH1D("Ebeam", "Photon beam energy (all)", 100, 1.0, 12.1);
	h_Ecoherent = new TH1D("Ecoherent", "Photon beam energy for coherently produced #gammas", 100, 1.0, 12.1);
	h_Eincoherent = new TH1D("Eincoherent", "Photon beam energy for incoherently produced #gammas", 100, 1.0, 12.1);

	h_thetabeam = new TH1D("thetabeam", "#theta_{beam}", 200, 0.0, 0.005);
	h_thetabeam_coherent = new TH1D("thetabeam_coherent", "#theta_{beam} for coherently produced #gammas", 200, 0.0, 0.005);
	h_thetabeam_incoherent = new TH1D("thetabeam_incoherent", "#theta_{beam} for incoherently produced #gammas", 200, 0.0, 0.005);
	h_phibeam = new TH1D("phibeam", "#phi_{beam}", 200, -180.0, 180.0);
	h_phibeam_coherent = new TH1D("phibeam_coherent", "#phi_{beam} for coherently produced #gammas", 200, -180.0, 180.0);
	h_phibeam_incoherent = new TH1D("phibeam_incoherent", "#phi_{beam} for incoherently produced #gammas", 200, -180.0, 180.0);

	h_Pmup = new TH1D("Pmup", "P_{#mu+}", 100, 0.0, 10.0);
	h_Pmup_coherent = new TH1D("Pmup_coherent", "P_{#mu+} for coherently produced #gammas", 100, 0.0, 10.0);
	h_Pmup_incoherent = new TH1D("Pmup_incoherent", "P_{#mu+} for incoherently produced #gammas", 100, 0.0, 10.0);
	h_Pmum = new TH1D("Pmum", "P_{#mu-}", 100, 0.0, 10.0);
	h_Pmum_coherent = new TH1D("Pmum_coherent", "P_{#mu-} for coherently produced #gammas", 100, 0.0, 10.0);
	h_Pmum_incoherent = new TH1D("Pmum_incoherent", "P_{#mu-} for incoherently produced #gammas", 100, 0.0, 10.0);
	
	h_thetamup = new TH1D("thetamup", "#theta_{#mu+}", 100, 0.0, 10.0);
	h_thetamup_coherent = new TH1D("thetamup_coherent", "#theta_{#mu+} for coherently produced #gammas", 100, 0.0, 10.0);
	h_thetamup_incoherent = new TH1D("thetamup_incoherent", "#theta_{#mu+} for incoherently produced #gammas", 100, 0.0, 10.0);
	h_thetamum = new TH1D("thetamum", "#theta_{#mu-}", 100, 0.0, 10.0);
	h_thetamum_coherent = new TH1D("thetamum_coherent", "#theta_{#mu-} for coherently produced #gammas", 100, 0.0, 10.0);
	h_thetamum_incoherent = new TH1D("thetamum_incoherent", "#theta_{#mu-} for incoherently produced #gammas", 100, 0.0, 10.0);

	h_phimup = new TH1D("phimup", "#phi_{#mu+}", 100, -180.0, 180.0);
	h_phimup_coherent = new TH1D("phimup_coherent", "#phi_{#mu+} for coherently produced #gammas", 100, -180.0, 180.0);
	h_phimup_incoherent = new TH1D("phimup_incoherent", "#phi_{#mu+} for incoherently produced #gammas", 100, -180.0, 180.0);
	h_phimum = new TH1D("phimum", "#phi_{#mu-}", 100, -180.0, 180.0);
	h_phimum_coherent = new TH1D("phimum_coherent", "#phi_{#mu-} for coherently produced #gammas", 100, -180.0, 180.0);
	h_phimum_incoherent = new TH1D("phimum_incoherent", "#phi_{#mu-} for incoherently produced #gammas", 100, -180.0, 180.0);

	h_Pmumu = new TH1D("Pmumu", "P_{#mu#mu}", 100, 0.0, 12.0);
	h_Pmumu_coherent = new TH1D("Pmumu_coherent", "P_{#mu#mu} for coherently produced #gammas", 100, 0.0, 12.0);
	h_Pmumu_incoherent = new TH1D("Pmumu_incoherent", "P_{#mu#mu} for incoherently produced #gammas", 100, 0.0, 12.0);
	h_thetamumu = new TH1D("thetamumu", "#theta_{#mu#mu}", 100, 0.0, 1.0);
	h_thetamumu_coherent = new TH1D("thetamumu_coherent", "#theta_{#mu#mu} for coherently produced #gammas", 100, 0.0, 1.0);
	h_thetamumu_incoherent = new TH1D("thetamumu_incoherent", "#theta_{#mu#mu} for incoherently produced #gammas", 100, 0.0, 1.0);
	h_phimumu = new TH1D("phimumu", "#phi_{#mu#mu}", 100, -180.0, 180.0);
	h_phimumu_coherent = new TH1D("phimumu_coherent", "#phi_{#mu#mu} for coherently produced #gammas", 100, -180.0, 180.0);
	h_phimumu_incoherent = new TH1D("phimumu_incoherent", "#phi_{#mu#mu} for incoherently produced #gammas", 100, -180.0, 180.0);

	japp->RootUnLock();

	return NOERROR;
}

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

//------------------
// evnt
//------------------
jerror_t JEventProcessor_genmumu::evnt(JEventLoop *loop, uint64_t eventnumber)
{
	vector<const DBeamPhoton*> beamphotons;
	vector<const DMCThrown*> mcthrowns;
	loop->Get(beamphotons, "MCGEN");
	loop->Get(mcthrowns);
	
	// Get DBeamPhoton object and guarantee it exists
	const DBeamPhoton *beam = (beamphotons.empty() ? NULL:beamphotons[0]);
	if(!beam){
		_DBG_<<"NO DBeamPhoton:MCGEN !!" << endl;
		exit(-1);
	}
	double Ebeam = beam->energy();
	
	// Get mu+/mu- objects and guarantee both exist
	const DMCThrown *muplus = NULL;
	const DMCThrown *muminus = NULL;
	for(uint32_t i=0; i<mcthrowns.size(); i++){
		const DMCThrown *thrown = mcthrowns[i];
		if(thrown->pdgtype == -13)muplus = thrown;
		if(thrown->pdgtype == +13)muminus = thrown;
	}
	if( (muplus==NULL) || (muminus==NULL) ){
		_DBG_<<"Missing mu+ or mu- !" << endl;
		exit(-1);
	}
	TLorentzVector vplus = muplus->lorentzMomentum();
	TLorentzVector vminus = muminus->lorentzMomentum();
	TLorentzVector vmumu = vplus + vminus;

	// Check if this was produced coherently or incoherently
	bool is_coherent = false;
	if(muplus && (muplus->mech==1) ) is_coherent = true;

	//--------- Lock ROOT mutex ---------
	japp->RootWriteLock();
	
	h_Ebeam->Fill(Ebeam);
	h_thetabeam->Fill(beam->momentum().Theta()*TMath::RadToDeg());
	h_phibeam->Fill(beam->momentum().Phi()*TMath::RadToDeg());
	h_Pmup->Fill(vplus.P());
	h_Pmum->Fill(vminus.P());
	h_thetamup->Fill(vplus.Theta()*TMath::RadToDeg());
	h_thetamum->Fill(vminus.Theta()*TMath::RadToDeg());
	h_phimup->Fill(vplus.Phi()*TMath::RadToDeg());
	h_phimum->Fill(vminus.Phi()*TMath::RadToDeg());
	h_Pmumu->Fill(vmumu.P());
	h_thetamumu->Fill(vmumu.Theta()*TMath::RadToDeg());
	h_phimumu->Fill(vmumu.Phi()*TMath::RadToDeg());
	if(is_coherent){
		h_Ecoherent->Fill(Ebeam);
		h_thetabeam_coherent->Fill(beam->momentum().Theta()*TMath::RadToDeg());
		h_phibeam_coherent->Fill(beam->momentum().Phi()*TMath::RadToDeg());
		h_Pmup_coherent->Fill(vplus.P());
		h_Pmum_coherent->Fill(vminus.P());
		h_thetamup_coherent->Fill(vplus.Theta()*TMath::RadToDeg());
		h_thetamum_coherent->Fill(vminus.Theta()*TMath::RadToDeg());
		h_phimup_coherent->Fill(vplus.Phi()*TMath::RadToDeg());
		h_phimum_coherent->Fill(vminus.Phi()*TMath::RadToDeg());
		h_Pmumu_coherent->Fill(vmumu.P());
		h_thetamumu_coherent->Fill(vmumu.Theta()*TMath::RadToDeg());
		h_phimumu_coherent->Fill(vmumu.Phi()*TMath::RadToDeg());
	}else{
		h_Eincoherent->Fill(Ebeam);
		h_thetabeam_incoherent->Fill(beam->momentum().Theta()*TMath::RadToDeg());
		h_phibeam_incoherent->Fill(beam->momentum().Phi()*TMath::RadToDeg());
		h_Pmup_incoherent->Fill(vplus.P());
		h_Pmum_incoherent->Fill(vminus.P());
		h_thetamup_incoherent->Fill(vplus.Theta()*TMath::RadToDeg());
		h_thetamum_incoherent->Fill(vminus.Theta()*TMath::RadToDeg());
		h_phimup_incoherent->Fill(vplus.Phi()*TMath::RadToDeg());
		h_phimum_incoherent->Fill(vminus.Phi()*TMath::RadToDeg());
		h_Pmumu_incoherent->Fill(vmumu.P());
		h_thetamumu_incoherent->Fill(vmumu.Theta()*TMath::RadToDeg());
		h_phimumu_incoherent->Fill(vmumu.Phi()*TMath::RadToDeg());
	}

	//-------- Unock ROOT mutex --------
	japp->RootUnLock();


	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_genmumu::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_genmumu::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}