// $Id$
//
//    File: DEventProcessor_npiPlus_singlePi.cc
// Created: Fri Jul 15 11:23:47 EDT 2016
// Creator: mpatsyuk (on Linux ifarm1102 2.6.32-431.el6.x86_64 x86_64)
//

#include "DEventProcessor_npiPlus_singlePi.h"

// Routine used to create our DEventProcessor

extern "C"
{
	void InitPlugin(JApplication *locApplication)
	{
		InitJANAPlugin(locApplication);
		locApplication->AddProcessor(new DEventProcessor_npiPlus_singlePi()); //register this plugin
		locApplication->AddFactoryGenerator(new DFactoryGenerator_npiPlus_singlePi()); //register the factory generator
	}
} // "C"

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

	/*
	//OPTIONAL: Create an EventStore skim.  
	string locSkimFileName = "npiPlus_singlePi.idxa";
	dEventStoreSkimStream.open(locSkimFileName.c_str());
	dEventStoreSkimStream << "IDXA" << endl;
	*/
  /*
  japp->RootWriteLock();
  {
    TDirectory* dir = new TDirectoryFile("MCTruthStudy","MCTruthStudy");
    dir->cd();

    dHist_nNeutrons = new TH1I("dHist_nNeutrons","number of neutrons in event",10,0.,10.);
    dHist_nPartners = new TH1I("dHist_nPartners","what particles are generated together with pion, PID",52,0.,52.);
    dHist_nPions = new TH1F("dHist_nPions", "n of thrown pi+ with the same parent as neutron; PID", 10, 0., 10.);
    dHist_nExtraTracks = new TH1F("dHist_nExtraTracks", "number of extra tracks in signal events", 20, 0., 20.);
    dHist_momExtraTracks = new TH1F("dHist_momExtraTracks", "p of extra tracks; p[GeV/c]", 400, 0., 10.);
    dHist_thetaExtraTracks = new TH1F("dHist_thetaExtraTracks", "#theta of extra tracks; #theta [degrees]", 180, 0., 180.);
    dHist_mom_theta_ET = new TH2F("dHist_mom_theta_ET","p vs #theta for extra tracks; #theta [degrees]",400,0.,10.,180,0.,180.);
    dHist_mom_theta_Pi = new TH2F("dHist_mom_theta_Pi","p vs #theta for pi+; #theta [degrees]",400,0.,10.,180,0.,180.);
    dHist_EShowerR = new TH2F("dHist_EShowerR","E of the shower vs radius; E [GeV]; R[cm]", 60,0.,6.,100,0.,100.);
    dHist_EShowerR_FCAL = new TH2F("dHist_EShowerR_FCAL","E of the FCAL shower vs radius; E [GeV]; R[cm]", 60,0.,6.,100,0.,100.);
    dHist_NShowers = new TH1I("dHist_NShowers","n of showers",10,0.,10.);
    
    dir->cd("../"); //return to the action directory

    dHist_IsEvent = new TH1D("IsEvent", "Is the event an event?", 2, -0.5, 1.5);
    dHist_IsEvent->GetXaxis()->SetBinLabel(1, "False");
    dHist_IsEvent->GetXaxis()->SetBinLabel(2, "True");
  }
  japp->RootUnLock();*/
	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_npiPlus_singlePi::brun(jana::JEventLoop* locEventLoop, int32_t locRunNumber)
{
	// This is called whenever the run number changes

	
  /*	//Recommended: Create output ROOT TTrees (nothing is done if already created)
	const DEventWriterROOT* locEventWriterROOT = NULL;
	locEventLoop->GetSingle(locEventWriterROOT);
	locEventWriterROOT->Create_DataTrees(locEventLoop);
  */	

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_npiPlus_singlePi::evnt(jana::JEventLoop* locEventLoop, uint64_t locEventNumber)
{
	// This is called for every event. Use of common resources like writing
	// to a file or filling a histogram should be mutex protected. Using
	// locEventLoop->Get(...) to get reconstructed objects (and thereby activating the
	// reconstruction algorithm) should be done outside of any mutex lock
	// since multiple threads may call this method at the same time.
	//
	// Here's an example:
	//
	// vector<const MyDataClass*> mydataclasses;
	// locEventLoop->Get(mydataclasses);
	//
	// japp->RootFillLock(this);
	//  ... fill historgrams or trees ...
	// japp->RootFillUnLock(this);

	// DOCUMENTATION:
	// ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software

	/*********************************************************** REQUIRED ***********************************************************/

	//REQUIRED: To run an analysis, You MUST call one at least of the below code fragments. 
		//JANA is on-demand, so if you don't call one of these, then your analysis won't run. 

	
	//Recommended: Write surviving particle combinations (if any) to output ROOT TTree
		//If no cuts are performed by the analysis actions added to a DReaction, then this saves all of its particle combinations. 
		//The event writer gets the DAnalysisResults objects from JANA, performing the analysis. 
	// string is DReaction factory tag: will fill trees for all DReactions that are defined in the specified factory
      	const DEventWriterROOT* locEventWriterROOT = NULL;
	locEventLoop->GetSingle(locEventWriterROOT);
       	locEventWriterROOT->Fill_DataTrees(locEventLoop, "npiPlus_singlePi");
  
	
	//Optional: Get the analysis results for all DReactions. 
		//Getting these objects triggers the analysis, if it wasn't performed already. 
		//These objects contain the DParticleCombo objects that survived the DAnalysisAction cuts that were added to the DReactions
	vector<const DAnalysisResults*> locAnalysisResultsVector;
	locEventLoop->Get(locAnalysisResultsVector);
	

	/************************************************** OPTIONAL: FURTHER ANALYSIS **************************************************/
	/*
	//Optional: Get all particle combinations for all DReactions. 
		//If kinematic fits were requested, these contain both the measured and kinematic-fit track parameters
		//No cuts from DAnalysisActions are placed on these combos
	vector<const DParticleCombo*> locParticleCombos;
	locEventLoop->Get(locParticleCombos);
	for(size_t loc_i = 0; loc_i < locParticleCombos.size(); ++loc_i)
	{
		const DParticleCombo* locParticleCombo = locParticleCombos[loc_i];
		if(locParticleCombo->Get_Reaction()->Get_ReactionName() != "npiPlus_singlePi")
			continue; // particle combination was for a different reaction
		//perform further analysis steps here...
	}
	*/
	/*
	vector<const DMCThrown*> locMCThrowns;                                                                                                                                                    
        locEventLoop->Get(locMCThrowns);
	int neutron_parentID = -100.;
	const DMCThrown* mcthrown = NULL;
	int nNeu = 0;

	vector<const DMCThrownMatching*> locMCThrownMatchingVector;
	locEventLoop->Get(locMCThrownMatchingVector);
	double locMatchingFOM;

	vector<const DChargedTrackHypothesis*> locChargedTrackHypotheses;
	locEventLoop->Get(locChargedTrackHypotheses);
	deque<const DChargedTrackHypothesis*> pions;
	int i_pi = 0;
       	bool PionIsMatched = 0;

	vector<const DNeutralShower*> locNeutralShowers;
	locEventLoop->Get(locNeutralShowers);
	//const DNeutralShower* piShower = NULL;
	//double locMatchingFOM_sho;

	vector<const DFCALTruthShower*> locFCALTruthShowers;
	locEventLoop->Get(locFCALTruthShowers);

	vector<const DFCALShower*> locFCALShowers;
	locEventLoop->Get(locFCALShowers);

	if(!locMCThrowns.empty()){

	  for(unsigned int jj=0; jj<locMCThrowns.size(); jj++){
	    mcthrown = locMCThrowns[jj];
	    if(mcthrown->type == 13){ // neutron, check the id of its parent
	      for(unsigned int pp=0; pp<locMCThrowns.size(); pp++){
		const DMCThrown* mcthrown2 = locMCThrowns[pp];
		if(mcthrown2->myid == mcthrown->parentid){
		  neutron_parentID = mcthrown2->type;
		  break;
		}
		//neutron_parentID = mcthrown->parentid;
	      }
	      break;
	    }
	  }
	  
	  if(neutron_parentID > -99.){// n is found
	    for(unsigned kk=0; kk<locMCThrowns.size(); kk++){
	      const DMCThrown* mcthrown1 = locMCThrowns[kk];
	      if(mcthrown1->type == 8 && mcthrown1->parentid == mcthrown->parentid){
		PionIsMatched = locMCThrownMatchingVector[0]->Get_MatchingChargedHypotheses(mcthrown1, pions, locMatchingFOM);
		//piShower = locMCThrownMatchingVector[0]->Get_MatchingNeutralShower(mcthrown1,locMatchingFOM_sho);
		i_pi++;
	      }
	    }
	  }

	  //count and plot tracks, which are not signal pions
	  japp->RootWriteLock();
	  {
	    if(i_pi == 1){// n is primary particle, number of pions born together with neutron == 1
	      dHist_mom_theta_Pi->Fill(pions[0]->pmag(), acos(pions[0]->pz()/pions[0]->pmag())/3.1415*180.);
	      dHist_nPions->Fill(i_pi);
	      dHist_nExtraTracks->Fill(locChargedTrackHypotheses.size()-1);
	      for(unsigned int tt=0; tt<locChargedTrackHypotheses.size(); tt++){//plot properties of the extra tracks
		const DChargedTrackHypothesis* trackhyp = locChargedTrackHypotheses[tt];
		if(trackhyp->pmag() == pions[0]->pmag()) continue; // this is the matched pion, skip it
		dHist_momExtraTracks->Fill(trackhyp->pmag());
		dHist_thetaExtraTracks->Fill(acos(trackhyp->pz()/trackhyp->pmag())/3.1415*180.);
		dHist_mom_theta_ET->Fill(trackhyp->pmag(), acos(trackhyp->pz()/trackhyp->pmag())/3.1415*180.);
	      }
	    }
	    if(neutron_parentID > -99.){
	      for(unsigned int gg=0; gg<locMCThrowns.size(); gg++){
		if(locMCThrowns[gg]->parentid == mcthrown->parentid){
		  dHist_nPartners->Fill(locMCThrowns[gg]->type);
		}
		if(locMCThrowns[gg]->type == 13){nNeu++;}
	      }
	      dHist_nNeutrons->Fill(nNeu);
	    }
	    //plot info about extra showers:
	    if(neutron_parentID > -99.){
	      dHist_NShowers->Fill(locFCALShowers.size());
	      for(unsigned int ss=0; ss<locNeutralShowers.size(); ss++){
		const DLorentzVector vectSho = locNeutralShowers[ss]->dSpacetimeVertex;
		double r = sqrt(pow(vectSho.X(),2) + pow(vectSho.Y(),2));
		dHist_EShowerR->Fill(locNeutralShowers[ss]->dEnergy, r);
	      }
	      
	      for(unsigned int mm=0; mm<locFCALShowers.size(); mm++){
		const DVector3 vectFCALSho = locFCALShowers[mm]->getPosition();
		double r1 = sqrt(pow(vectFCALSho.X(),2) + pow(vectFCALSho.Y(),2));
		dHist_EShowerR_FCAL->Fill(locFCALShowers[mm]->getEnergy(), r1);
	      }
	    }
	    
	  }
	  japp->RootUnLock();

	}//if mcthrowns not empty
	*/
	/*
	//Optional: Perform further cuts on the particle combos in the analysis results. 
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "npiPlus_singlePi")
			continue; // analysis results were for a different reaction
		//get the DParticleCombo objects for this DReaction that survived all of the DAnalysisAction cuts
		deque<const DParticleCombo*> locPassedParticleCombos;
		locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);
		for(size_t loc_j = 0; loc_j < locPassedParticleCombos.size(); ++loc_j)
		{
			const DParticleCombo* locPassedParticleCombo = locPassedParticleCombos[loc_j];
			//perform further analysis steps here...
		}
	}
	*/

	/******************************************************** OPTIONAL: SKIMS *******************************************************/

	/*
	//Optional: Save event to output REST file. Use this to create physical skims.
       	const DEventWriterREST* locEventWriterREST = NULL;
	locEventLoop->GetSingle(locEventWriterREST);
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "npiPlus_singlePi")
			continue; // analysis results were for a different reaction

		//get the DParticleCombo objects for this DReaction that survived all of the DAnalysisAction cuts
		deque<const DParticleCombo*> locPassedParticleCombos;
		locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos);

		if(!locPassedParticleCombos.empty())
			locEventWriterREST->Write_RESTEvent(locEventLoop, "npiPlus_singlePi"); //string is part of output file name
	}
	*/

	/*
	//Optional: Create an EventStore skim. 

	// See whether this is MC data or real data
	vector<const DMCThrown*> locMCThrowns;
	locEventLoop->Get(locMCThrowns);

	unsigned int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
	unsigned int locUniqueID = locMCThrowns.empty() ? 1 : Get_FileNumber(locEventLoop);

	// If a particle combo passed the cuts, save the event info in the output file
	for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
	{
		const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i];
		if(locAnalysisResults->Get_Reaction()->Get_ReactionName() != "npiPlus_singlePi")
			continue; // analysis results were for a different reaction
		if(locAnalysisResults->Get_NumPassedParticleCombos() == 0)
			continue; // no combos passed

		//MUST LOCK AROUND MODIFICATION OF MEMBER VARIABLES IN brun() or evnt().
		japp->WriteLock("npiPlus_singlePi.idxa"); //Lock is unique to this output file
		{
			dEventStoreSkimStream << locRunNumber << " " << locEventNumber << " " << locUniqueID << endl;
		}
		japp->Unlock("npiPlus_singlePi.idxa");
	}
	*/

	return NOERROR;
}

int DEventProcessor_npiPlus_singlePi::Get_FileNumber(JEventLoop* locEventLoop) const
{
	//Assume that the file name is in the format: *_X.ext, where:
		//X is the file number (a string of numbers of any length)
		//ext is the file extension (probably .evio or .hddm)

	//get the event source
	JEventSource* locEventSource = locEventLoop->GetJEvent().GetJEventSource();
	if(locEventSource == NULL)
		return -1;

	//get the source file name (strip the path)
	string locSourceFileName = locEventSource->GetSourceName();

	//find the last "_" & "." indices
	size_t locUnderscoreIndex = locSourceFileName.rfind("_");
	size_t locDotIndex = locSourceFileName.rfind(".");
	if((locUnderscoreIndex == string::npos) || (locDotIndex == string::npos))
		return -1;

	size_t locNumberLength = locDotIndex - locUnderscoreIndex - 1;
	string locFileNumberString = locSourceFileName.substr(locUnderscoreIndex + 1, locNumberLength);

	int locFileNumber = -1;
	istringstream locFileNumberStream(locFileNumberString);
	locFileNumberStream >> locFileNumber;

	return locFileNumber;
}

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