// $Id$
//
//    File: JEventProcessor_fcal2.cc
// Created: Fri Apr 26 18:14:10 EDT 2024
// Creator: staylor (on Linux ifarm1901.jlab.org 3.10.0-1160.108.1.el7.x86_64 x86_64)
//

#include "JEventProcessor_fcal2.h"
using namespace jana;

#include <FCAL/DFCALShower.h>
#include <FCAL/DFCALCluster.h>
#include <FCAL/DFCALHit.h>
#include <FCAL/DFCALGeometry.h>
#include <FCAL/DFCALTruthShower.h>
#include <ECAL/DECALHit.h>
#include <ECAL/DECALTruthShower.h>
#include <TRACKING/DMCThrown.h>
#include <PID/DChargedTrack.h>
#include <TDirectory.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_fcal2());
  }
} // "C"


//------------------
// JEventProcessor_fcal2 (Constructor)
//------------------
JEventProcessor_fcal2::JEventProcessor_fcal2()
{

}

//------------------
// ~JEventProcessor_fcal2 (Destructor)
//------------------
JEventProcessor_fcal2::~JEventProcessor_fcal2()
{
 
}

//------------------
// init
//------------------
jerror_t JEventProcessor_fcal2::init(void)
{
  gPARMS->SetDefaultParameter("FCAL2:THRESHOLD",THRESHOLD);
  gPARMS->SetDefaultParameter("FCAL2:P_MIN",P_MIN); 
  gPARMS->SetDefaultParameter("FCAL2:P_MAX",P_MAX); 
  gPARMS->SetDefaultParameter("FCAL2:R_MIN",R_MIN); 
  gPARMS->SetDefaultParameter("FCAL2:R_MAX",R_MAX);

  gDirectory->mkdir("ECAL")->cd();
  {
    ECALOccupancy = new TH2F("ECALOccupancy",";column;row;",40,-0.5,39.5,40,-0.5,39.5);
    ECAL_Etrue_vs_Ecluster=new TH2F("ECAL_Etrue_vs_Ecluster",";E(cluster) [GeV]; E(thrown) [GeV]",100,0,10,500,0,10);
    ECAL_Normalized_Eresidual_vs_E=new TH2F("ECAL_Normalized_Eresidual_vs_E",
					    ";E(true)[GeV];#DeltaE/E",100,0,10,
					  120,-0.3,0.3);
    ECAL_dx_vs_E=new TH2F("ECAL_dx_vs_E",";E(true)[GeV];#Deltax [cm]",100,0,10,
			  100,-5,5);
    ECAL_dy_vs_E=new TH2F("ECAL_dy_vs_E",";E(true)[GeV];#Deltay [cm]",100,0,10,
			  100,-5,5);
    ECAL_dx_vs_x=new TH2F("ECAL_dx_vs_x",";x[cm];#Deltax [cm]",1800,-45,45,
			  100,-5,5);
    ECAL_dy_vs_y=new TH2F("ECAL_dy_vs_y",";y[cm];#Deltay [cm]",1800,-45,45,
			  100,-5,5);

    gDirectory->mkdir("TwoShowers")->cd();
    {
      ECAL_ShowerSeparation=new TH2F("ECAL_ShowerSeparation",";#Deltar(true)[cm];#Deltar [cm]",1000,0,100,1000,0,100);
      ECAL_TrueSeparation=new TH1F("ECAL_TrueSeparation",";#Deltar(true)[cm]",
				   1000,0,100);  
      ECAL_ShowerSeparation_NoSec=new TH2F("ECAL_ShowerSeparation_NoSec",
					   ";#Deltar(true)[cm];#Deltar [cm]",
					   1000,0,100,1000,0,100);
      ECAL_TrueSeparation_NoSec=new TH1F("ECAL_TrueSeparation_NoSec",
					 ";#Deltar(true)[cm]",
					 1000,0,100);
    }
  }
  gDirectory->cd("../../");

  gDirectory->mkdir("FCAL")->cd();
  {
    FCALOccupancy = new TH2F("FCALOccupancy",";column;row;",59,-0.5,58.5,59,-0.5,58.5);
    FCAL_Etrue_vs_Ecluster=new TH2F("FCAL_Etrue_vs_Ecluster",";E(cluster) [GeV]; E(thrown) [GeV]",100,0,10,500,0,10);
    FCAL_Normalized_Eresidual_vs_E=new TH2F("FCAL_Normalized_Eresidual_vs_E",
					    ";E(true)[GeV];#DeltaE/E",100,0,10,
					    120,-0.3,0.3);
    FCAL_dx_vs_E=new TH2F("FCAL_dx_vs_E",";E(true)[GeV];#Deltax [cm]",100,0,10,
			  100,-5,5);
    FCAL_dy_vs_E=new TH2F("FCAL_dy_vs_E",";E(true)[GeV];#Deltay [cm]",100,0,10,
			  100,-5,5);
    FCAL_dx_vs_x=new TH2F("FCAL_dx_vs_x",";x[cm];#Deltax [cm]",4800,-120,120,
			  100,-5,5);
    FCAL_dy_vs_y=new TH2F("FCAL_dy_vs_y",";y[cm];#Deltay [cm]",4800,-120,120,
			  100,-5,5);

    gDirectory->mkdir("TwoShowers")->cd();
    {
      FCAL_ShowerSeparation=new TH2F("FCAL_ShowerSeparation",";#Deltar(true)[cm];#Deltar [cm]",1000,0,100,1000,0,100);
      FCAL_TrueSeparation=new TH1F("FCAL_TrueSeparation",";#Deltar(true)[cm]",
				   1000,0,100);   
      FCAL_ShowerSeparation_NoSec=new TH2F("FCAL_ShowerSeparation_NoSec",";#Deltar(true)[cm];#Deltar [cm]",1000,0,100,1000,0,100);
      FCAL_TrueSeparation_NoSec=new TH1F("FCAL_TrueSeparation_NoSec",
					 ";#Deltar(true)[cm]",
					 1000,0,100);
    }
  }
  gDirectory->cd("../../");

  gDirectory->mkdir("Combined")->cd();
  {
    NumShowers=new TH1F("NumShowers",";Number of neutral showers",20,-0.5,19.5);
    ECALFCAL_ShowerSeparation=new TH2F("ECALFCAL_ShowerSeparation",";#Deltar(true)[cm];#Deltar [cm]",1000,0,100,1000,0,100);
    ECALFCAL_TrueSeparation=new TH1F("ECALFCAL_TrueSeparation",
				  ";#Deltar(true)[cm]",1000,0,100);
  }
  gDirectory->cd("../");
  gDirectory->mkdir("Tracks")->cd();
  gDirectory->mkdir("Electrons")->cd();
  {
    Exy_electron=new TH2F("Exy_electron","E_{i}/p vs x vs y;x[cm];y[cm]",400,-20,20,400,-20,20);
    HEoverP_electron=new TH2F("HEoverP_electron",";p [GeV/c];E/p",100,0,10,100,0,1.5);   
    Hdx_vs_x_electron=new TH2F("Hdx_vs_x_electron",";x [cm];#deltax [cm]",100,-2.5,2.5,100,-2.5,2.5);
    Hdy_vs_y_electron=new TH2F("Hdy_vs_y_electron",";y [cm];#deltay [cm]",100,-2.5,2.5,100,-2.5,2.5);
  }
  gDirectory->cd("../");
  gDirectory->mkdir("Positrons")->cd();
  {
    Exy_positron=new TH2F("Exy_positron","E_{i}/p vs x vs y;x[cm];y[cm]",400,-20,20,400,-20,20);
    HEoverP_positron=new TH2F("HEoverP_positron",";p [GeV/c];E/p",100,0,10,100,0,1.5);   
    Hdx_vs_x_positron=new TH2F("Hdx_vs_x_positron",";x [cm];#deltax [cm]",100,-2.5,2.5,100,-2.5,2.5);
    Hdy_vs_y_positron=new TH2F("Hdy_vs_y_positron",";y [cm];#deltay [cm]",100,-2.5,2.5,100,-2.5,2.5);

  }
  gDirectory->cd("../../");

  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_fcal2::brun(JEventLoop *eventLoop, int32_t runnumber)
{
  eventLoop->GetSingle(fcalGeom);

  // This is called whenever the run number changes
  return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_fcal2::evnt(JEventLoop *loop, uint64_t eventnumber)
{
  vector<const DMCThrown *>throwns;
  loop->Get(throwns);

  vector<const DChargedTrack *>tracks;
  loop->Get(tracks);

  vector<const DFCALCluster *>fcalclusters;
  loop->Get(fcalclusters);

  vector<const DFCALShower *>fcalshowers;
  loop->Get(fcalshowers);

  vector<const DFCALHit*>fcalhits;
  loop->Get(fcalhits);

  vector<const DECALHit*>ecalhits;
  loop->Get(ecalhits);

  vector<const DFCALTruthShower *>fcaltruths;
  loop->Get(fcaltruths);
  vector<const DFCALTruthShower *>myfcaltruths;
  for (unsigned int i=0;i<fcaltruths.size();i++){
    const DFCALTruthShower *myfcaltruth=fcaltruths[i];
    if (myfcaltruth->primary()>0 && myfcaltruth->E()>THRESHOLD){
      myfcaltruths.push_back(myfcaltruth);
    }
  }

  vector<const DECALTruthShower *>ecaltruths;
  loop->Get(ecaltruths);
  vector<const DECALTruthShower *>myecaltruths;
  for (unsigned int i=0;i<ecaltruths.size();i++){
    const DECALTruthShower *myecaltruth=ecaltruths[i];
    if (myecaltruth->primary()>0 && myecaltruth->E()>THRESHOLD){
      myecaltruths.push_back(myecaltruth);
    }
  }

  japp->RootWriteLock();

  NumShowers->Fill(fcalshowers.size());

  for (unsigned int i=0;i<ecalhits.size();i++){ 
    ECALOccupancy->Fill(ecalhits[i]->column,ecalhits[i]->row);
  }
  for (unsigned int i=0;i<fcalhits.size();i++){ 
    FCALOccupancy->Fill(fcalhits[i]->column,fcalhits[i]->row);
  }

  double dr=1e6;
  if (fcalshowers.size()==2){
    DVector3 diff=fcalshowers[0]->getPosition()-fcalshowers[1]->getPosition();
    dr=diff.Perp();
  }
  if (ecaltruths.size()==2&&myecaltruths.size()==2){
    double dxtrue=ecaltruths[0]->x()-ecaltruths[1]->x();
    double dytrue=ecaltruths[0]->y()-ecaltruths[1]->y();
    double drtrue=sqrt(dxtrue*dxtrue+dytrue*dytrue);
    ECAL_TrueSeparation_NoSec->Fill(drtrue);
    if (dr<1e6){
      ECAL_ShowerSeparation_NoSec->Fill(drtrue,dr);
    }
  }
  if (myecaltruths.size()==2){
    double dxtrue=myecaltruths[0]->x()-myecaltruths[1]->x();
    double dytrue=myecaltruths[0]->y()-myecaltruths[1]->y();
    double drtrue=sqrt(dxtrue*dxtrue+dytrue*dytrue);
    ECAL_TrueSeparation->Fill(drtrue);
    if (dr<1e6){
      ECAL_ShowerSeparation->Fill(drtrue,dr);
    }
  }
  if (fcaltruths.size()==2&&myfcaltruths.size()==2){
    double dxtrue=fcaltruths[0]->x()-fcaltruths[1]->x();
    double dytrue=fcaltruths[0]->y()-fcaltruths[1]->y();
    double drtrue=sqrt(dxtrue*dxtrue+dytrue*dytrue);
    FCAL_TrueSeparation_NoSec->Fill(drtrue);
    if (dr<1e6){
      FCAL_ShowerSeparation_NoSec->Fill(drtrue,dr);
    }
  }
  if (myfcaltruths.size()==2){
    double dxtrue=myfcaltruths[0]->x()-myfcaltruths[1]->x();
    double dytrue=myfcaltruths[0]->y()-myfcaltruths[1]->y();
    double drtrue=sqrt(dxtrue*dxtrue+dytrue*dytrue);
    FCAL_TrueSeparation->Fill(drtrue);
    if (dr<1e6){
      FCAL_ShowerSeparation->Fill(drtrue,dr);
    }
  }
  if (myfcaltruths.size()==1 && myecaltruths.size()==1){
    double dxtrue=myfcaltruths[0]->x()-myecaltruths[0]->x();
    double dytrue=myfcaltruths[0]->y()-myecaltruths[0]->y();
    double drtrue=sqrt(dxtrue*dxtrue+dytrue*dytrue);
    ECALFCAL_TrueSeparation->Fill(drtrue);
    if (dr<1e6){
      ECALFCAL_ShowerSeparation->Fill(drtrue,dr);
    }
  }

  if (fcalclusters.size()==1&&fcalshowers.size()==1){
    double E=fcalclusters[0]->getEnergy();
    double Eshower=fcalshowers[0]->getEnergy();
    DVector3 pos=fcalshowers[0]->getPosition();
    
    if (myecaltruths.size()==1){
      double Etrue=myecaltruths[0]->E();
      ECAL_Etrue_vs_Ecluster->Fill(E,Etrue);

      double dz=pos.z()-myecaltruths[0]->z();
      double pz=myecaltruths[0]->pz();
      double xtrue=myecaltruths[0]->x()+(myecaltruths[0]->px()/pz)*dz;  
      double ytrue=myecaltruths[0]->y()+(myecaltruths[0]->py()/pz)*dz;
  
      double dE_over_E=(Eshower-Etrue)/Etrue;
      double dx=pos.x()-xtrue;
      double dy=pos.y()-ytrue;

      ECAL_Normalized_Eresidual_vs_E->Fill(Etrue,dE_over_E);
      ECAL_dx_vs_E->Fill(Etrue,dx);
      ECAL_dy_vs_E->Fill(Etrue,dy);
      ECAL_dy_vs_y->Fill(pos.y(),dy);
      ECAL_dx_vs_x->Fill(pos.x(),dx);
    }
    if (myfcaltruths.size()==1){
      double Etrue=myfcaltruths[0]->E();
      FCAL_Etrue_vs_Ecluster->Fill(E,Etrue);

      double dz=pos.z()-myfcaltruths[0]->z();
      double pz=myfcaltruths[0]->pz();
      double xtrue=myfcaltruths[0]->x()+(myfcaltruths[0]->px()/pz)*dz;  
      double ytrue=myfcaltruths[0]->y()+(myfcaltruths[0]->py()/pz)*dz;

      double dE_over_E=(Eshower-Etrue)/Etrue;
      double dx=pos.x()-xtrue;
      double dy=pos.y()-ytrue;

      FCAL_Normalized_Eresidual_vs_E->Fill(Etrue,dE_over_E);
      FCAL_dx_vs_E->Fill(Etrue,dx);
      FCAL_dy_vs_E->Fill(Etrue,dy);
      FCAL_dy_vs_y->Fill(pos.y(),dy);
      FCAL_dx_vs_x->Fill(pos.x(),dx);
    }
  } 

  for (unsigned int i=0;i<tracks.size();i++){
    const DChargedTrackHypothesis *hyp=tracks[i]->Get_Hypothesis(Electron);
    if (hyp){
      FillLeptonHistos(hyp,HEoverP_electron,Exy_electron,Hdx_vs_x_electron,
		       Hdy_vs_y_electron);
    }   
    else {
      hyp=tracks[i]->Get_Hypothesis(Positron);
      if (hyp){
	FillLeptonHistos(hyp,HEoverP_positron,Exy_positron,Hdx_vs_x_positron,
			 Hdy_vs_y_positron);
      }
    }
  }

  japp->RootUnLock();
  
  return NOERROR;
}

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

void JEventProcessor_fcal2::FillLeptonHistos(const DChargedTrackHypothesis *hyp,
					     TH2F *HEoverP,TH2F *Exy,TH2F *Hdx,
					     TH2F *Hdy){
  const DTrackTimeBased *track=hyp->Get_TrackTimeBased();	
  //double p_at_target=track->momentum().Mag();
  vector<DTrackFitter::Extrapolation_t>fcal_extraps=track->extrapolations.at(SYS_FCAL); 
  shared_ptr<const DFCALShowerMatchParams>fcalparms=hyp->Get_FCALShowerMatchParams();
  if (fcalparms!=NULL&& fcal_extraps.size()>0){
    DVector3 track_mom=fcal_extraps[0].momentum;
    double p=fcal_extraps[0].momentum.Mag();
    double Eshower=fcalparms->dFCALShower->getEnergy();
    double E_over_p=Eshower/p;
    HEoverP->Fill(p,E_over_p); 
    
    if (E_over_p>0.8 && p>P_MIN && p<P_MAX){
      const DFCALCluster *cluster=NULL;
      fcalparms->dFCALShower->GetSingle(cluster);
      if (cluster!=NULL){
	vector<DFCALCluster::DFCALClusterHit_t>hits=cluster->GetHits();
	if (hits.size()>2){
	  DVector3 track_pos=fcal_extraps[0].position;
	  DVector3 pos=fcalparms->dFCALShower->getPosition();
	  track_pos+=((pos.z()-track_pos.z())/track_mom.z())*track_mom;
	  double R=track_pos.Perp();
	  if (R>R_MIN && R<R_MAX){
	    DVector3 diff=pos-track_pos;
	    
	    for (unsigned int j=0;j<hits.size();j++){
	      double dxt=track_pos.x()-hits[j].x;
	      double dyt=track_pos.y()-hits[j].y;
	      Exy->Fill(dxt,dyt,hits[j].E/p);
	    } 
	    
	    int channel=cluster->getChannelEmax();
	    DVector2 xyblock_Emax=fcalGeom->positionOnFace(channel);   
	    double dx_Emax=track_pos.x()-xyblock_Emax.X();
	    double dy_Emax=track_pos.y()-xyblock_Emax.Y();
	    Hdx->Fill(dx_Emax,diff.x());
	    Hdy->Fill(dy_Emax,diff.y());
	  }
	}
      } 
    }
  }
}