// $Id$
//
//    File: JEventProcessor_pi0.cc
// Created: Thu Apr  6 11:51:33 EDT 2017
// Creator: zihlmann (on Linux gluon49.jlab.org 2.6.32-431.20.3.el6.x86_64 x86_64)
//

#include "JEventProcessor_pi0.h"
using namespace jana;


// 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_pi0());
  }
} // "C"


//------------------
// JEventProcessor_pi0 (Constructor)
//------------------
JEventProcessor_pi0::JEventProcessor_pi0()
{
  
}

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

//------------------
// init
//------------------
jerror_t JEventProcessor_pi0::init(void)
{
  // This is called once at program startup. 
  ECounter = 0;

  japp->RootFillLock(this);

  CreateHistogramsAndTree();

  japp->RootFillUnLock(this);

  return NOERROR;
}


//------------------
// brun
//------------------
jerror_t JEventProcessor_pi0::brun(JEventLoop *eventLoop, int32_t runnumber)
{
  // This is called whenever the run number changes


  return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_pi0::evnt(JEventLoop *loop, uint64_t eventnumber)
{
  // This is called for every event. Use of common resources like writing
  // to a file or filling a histogram should be mutex protected. Using
  // loop->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;
  // loop->Get(mydataclasses);
  //
  // japp->RootFillLock(this);
  //  ... fill historgrams or trees ...
  // japp->RootFillUnLock(this);

  int CheckTrigger = 1;
  if (CheckTrigger) {

    vector <const DL1Trigger*> trig;
    loop->Get(trig);
    
    if (trig.size() == 0){
      //cout<<"no trigger found"<<endl;
      return NOERROR;     
    }
    
    if (trig[0]->trig_mask > 7){ // do not look at PS triggers
      return NOERROR;
    }
    if (trig[0]->trig_mask < 1){ // require physics trigger
      return NOERROR;
    }
  }

  vector <const DTrackTimeBased*> TrkCand;
  loop->Get(TrkCand);
  if (TrkCand.size()>0){ // only look at no charged track events
    return NOERROR;
  }

  vector <const DNeutralParticle*> Neutrals;
  loop->Get(Neutrals);
  int NNeutrals = Neutrals.size();
  if (NNeutrals<2){  // require at least two photons/neutrals
    return NOERROR;
  }

  vector <const DSCHit*> SCHits;
  loop->Get(SCHits);
  int NSCHits = SCHits.size();


  vector <const DBeamPhoton*> Beam;
  loop->Get(Beam);
  int NBeamPhotons = Beam.size();
  if (NBeamPhotons<1){
    return NOERROR;
  }

  vector <const DRFTDCDigiTime*> RFTimes;
  loop->Get(RFTimes);

  vector <uint32_t> T;
  for (unsigned int k=0; k<RFTimes.size(); k++){
    const DRFTDCDigiTime* RFT = RFTimes[k];
    if (RFT->dIsCAENTDCFlag){
      T.push_back(RFT->time);
    }
  }

  vector <const DEventRFBunch*> RFBunch;
  loop->Get(RFBunch);
  int NRFBunches = RFBunch.size();

  vector <DLorentzVector> Photons;
  vector <DLorentzVector> ShowerSpaceTime;
  vector <int> DetectorSystem;

  for (int k=0; k<NNeutrals; k++){
    
    //const DNeutralParticleHypothesis *gamHyp = Neutrals[k]->Get_Hypothesis(Gamma);
    //DLorentzVector gamma = gamHyp->lorentzMomentum();

    const DNeutralShower *shower = Neutrals[k]->dNeutralShower;
    DLorentzVector vec1 = shower->dSpacetimeVertex;
    float E = shower->dEnergy;
    // assume VERTEX at center of the target 
    // cout<<vec1.X()<<"  "<<vec1.Y()<<"  "<<vec1.Z()<<endl; 
    double R2 = vec1.X()*vec1.X() + vec1.Y()*vec1.Y() + (vec1.Z()-65.)*(vec1.Z()-65.);
    double R = TMath::Sqrt(R2);
    DLorentzVector gamma(vec1.X()/R*E, vec1.Y()/R*E, (vec1.Z()-65.)/R*E, E);
    
    if (gamma.E()>0.15){           // Photon Energy Cut 100MeV
      Photons.push_back(gamma);
      ShowerSpaceTime.push_back(shower->dSpacetimeVertex);
      int CAL = 1;
      if (shower->dDetectorSystem==SYS_BCAL){
	CAL = 0;
      }
      DetectorSystem.push_back(CAL);
    }
    
  }

  if (Photons.size()<2){
    return NOERROR;
  }

  for (int k=0; k<NRFBunches; k++){
    rfbunch->Fill(RFBunch[k]->dTime);
  }

  // look at raw SC hits
  for (int n=0;n<NSCHits;n++){
    if ((SCHits[n]->has_fADC) && (SCHits[n]->has_TDC)){
      sctimeraw->Fill(SCHits[n]->t);
      scenergyraw->Fill(SCHits[n]->dE);
    }
  }


  int DSYS = 2;
  for (unsigned int s=0; s<DetectorSystem.size(); s++){
    if (DetectorSystem[s] == 0){
      DSYS = 1;
    }
  }
  eventtype->Fill(DSYS);
   
  //rfTOF->Fill(TOF_RF_Time);
  
  vector <int> FoundBeam;
  int RFCUT = 0;
  for (int k=0; k<NBeamPhotons; k++){
    BeamPhotonE->Fill(Beam[k]->momentum().Mag());
    BeamPhotonT->Fill(Beam[k]->time());
    float dt = Beam[k]->time() - RFBunch[0]->dTime;
    float dt1 = Beam[k]->time() - RFBunch[0]->dTime;
    rfbeam->Fill(dt);
    rfbeamNoOffset->Fill(dt1);
    
    if (Beam[k]->dSystem == SYS_TAGM){
      rfbeamM->Fill(dt1);
      mictime->Fill(dt1, Beam[k]->dCounter);  
      if (DSYS==2){
	rfbeamMNOBCAL->Fill(dt1);
      }
    } else {
      rfbeamH->Fill(dt1);
      hodtime->Fill(dt1, Beam[k]->dCounter);  
      if (DSYS==2){
	rfbeamHNOBCAL->Fill(dt1);
      }
    }
    
    
    if (TMath::Abs(dt1)<1.){
      FoundBeam.push_back(k);
      FoundBeamPhotonE->Fill(Beam[k]->momentum().Mag());
      FoundBeamPhotonT->Fill(Beam[k]->time());
    }

    if (TMath::Abs(dt)<22.){
      RFCUT = 1;
    }
   

  }
  FoundBeamPhotons->Fill(FoundBeam.size());

  for (unsigned int k=0; k<Photons.size()-1; k++){
    if (DetectorSystem[k]){
      PhotonTFcal->Fill(ShowerSpaceTime[k].T());
    } else {
      PhotonTBcal->Fill(ShowerSpaceTime[k].T());
    }
  }
  int FoundSCHit = 0;
  for (int k=0; k<NSCHits; k++){
    const DSCHit* hit = SCHits[k];
    
    if ((hit->dE > 0.0015) && (hit->has_fADC) && (hit->has_TDC)){
      FoundSCHit = 1;
    }
  }

  vector <const DSCHit*> GoodSCHits;

  int FoundSCHitT = 0;
  for (int k=0; k<NSCHits; k++){
    const DSCHit* hit = SCHits[k];
    if ((hit->t > -20.) && (hit->t<30.) && (hit->has_fADC) && (hit->has_TDC)){
      FoundSCHitT = 1;
      if ((hit->dE > 0.0015) && (hit->has_fADC) && (hit->has_TDC)){
	GoodSCHits.push_back(hit);
      }
    }
  }
  
  if (FoundBeam.size()>0) {
    
    for (int k=0; k<NSCHits; k++){
      const DSCHit* hit = SCHits[k];
      if (( hit->has_fADC) && (hit->has_TDC)) {
        if (FoundSCHitT){
          SCE[hit->sector-1]->Fill(hit->dE);
        }
        if (FoundSCHit){
          SCT[hit->sector-1]->Fill(hit->t);      
        }
      }
    }
  }

  if ((FoundBeam.size()>0) && (FoundSCHit) && (FoundSCHitT) && (RFCUT) ) {
    
    //rfTOFAfter->Fill(TOF_RF_Time);
    
    for (unsigned int k=0; k<Photons.size()-1; k++){
      if (DetectorSystem[k]){
        FoundPhotonTFcal->Fill(ShowerSpaceTime[k].T());
        FoundPhotonEFcal->Fill(Photons[k].E());
      } else {
        FoundPhotonTBcal->Fill(ShowerSpaceTime[k].T());
        FoundPhotonEBcal->Fill(Photons[k].E());
      }

      DLorentzVector vec1 = ShowerSpaceTime[k];
      double R2 = vec1.X()*vec1.X() + vec1.Y()*vec1.Y() + (vec1.Z()-65.)*(vec1.Z()-65.);
      double R = TMath::Sqrt(R2);
      double TargetTime1 = vec1.T() - R/29.97;

      for (unsigned int j=k+1; j<Photons.size(); j++){

	DLorentzVector vec2 = ShowerSpaceTime[j];
	double R2_2 = vec2.X()*vec2.X() + vec2.Y()*vec2.Y() + (vec2.Z()-65.)*(vec2.Z()-65.);
	double R_2 = TMath::Sqrt(R2_2);
	double TargetTime2 = vec2.T() - R_2/29.97;
	

	// match time of the two photons
	if (TMath::Abs(TargetTime2-TargetTime1)<2.) {

	  DLorentzVector Ptot = Photons[k] + Photons[j];
	  if (Photons[k].E()>Photons[j].E()){
	    PE1vsPE2->Fill(Photons[j].E(),Photons[k].E());
	  } else {
	    PE1vsPE2->Fill(Photons[k].E(),Photons[j].E());
	  }
	  
	  double phi = 0;
	  if (Ptot.X()>0){
	    if (Ptot.Y()>0){
	      phi = TMath::ATan(Ptot.Y()/Ptot.X());
	    } else {
	      phi = 3.1415926*2. + TMath::ATan(Ptot.Y()/Ptot.X()); 
	    }
	  } else if (Ptot.X()<0){
            phi = 3.1415926 + TMath::ATan(Ptot.Y()/Ptot.X());     
	  }
	  PHI->Fill(phi);
	  double theta = Ptot.Theta();
	  THETA->Fill(theta);
	  
	  for (unsigned int s=0; s<GoodSCHits.size(); s++){
	    double phiSC = (GoodSCHits[s]->sector+1.) * 3.1415926*2./30. + 2.1415926/2.;
	    
	    double dphi = (phi+3.14159)-phiSC;
	    DPHI->Fill(dphi);
	  }
	  
	  for (unsigned int n=0; n<FoundBeam.size(); n++){
	    DLorentzVector v4 = Beam[FoundBeam[n]]->lorentzMomentum();
	    DLorentzVector d4 = v4-Ptot;
	    MissingP->Fill(d4.P());
	    MissingE->Fill(d4.E());
	    MissingPperp->Fill(d4.Pt());
	    MissingM->Fill(d4.M());
	  }
	  
	  double InvariantMass = Ptot.M();
	  twophotonm->Fill(InvariantMass);
	  if ( Photons.size() == 2) {
	    twophotonmextended->Fill(InvariantMass);
	  }
	  //cout<<InvariantMass<<endl;


	  for (unsigned int s=0; s< GoodSCHits.size(); s++){
	    const DSCHit* hit = GoodSCHits[s];
	    tSC_vs_tgamma->Fill(TargetTime1, hit->t);
	    ESC_vs_Mgamma->Fill(InvariantMass, hit->dE);
	  }

	}
      }
    }
  }


  // Do fill tree here:
  //if ( (Neutrals.size()>1)  && (NSCHits>0) && (RFCUT) && (Beam.size()>0) ){
  if ( (Neutrals.size()>1)  && (RFCUT) ){
 
    japp->RootFillLock(this);

    EventNUM = ECounter;
    RFT=999;
    if (NRFBunches){
      RFT = RFBunch[0]->dTime;
    }

    Nbeam = 0;
    int NMax = Beam.size();

    for (int k=0; k<NMax; k++) {
      if (Nbeam<400) {
	BeamE[Nbeam] = Beam[k]->momentum().Mag();
	BeamT[Nbeam] = Beam[k]->time();
	Nbeam++;
      }
    }
    

    Nneut = 0;
    NMax = Neutrals.size();
    if (NMax>20){
      NMax = 20;
    }
    for (int k=0; k<NMax; k++) {
      const DNeutralShower *shower = Neutrals[k]->dNeutralShower;
      DLorentzVector vec1 = shower->dSpacetimeVertex;
      float E = shower->dEnergy;
      if (E>0.15){
	NeutralX[Nneut] =  vec1.X();
	NeutralY[Nneut] =  vec1.Y();
	NeutralZ[Nneut] =  vec1.Z();
	NeutralT[Nneut] =  vec1.T();
	NeutralE[Nneut] =  E;
	int sys = 1;
	if (shower->dDetectorSystem==SYS_BCAL){
	  sys = 0;
	}
	DetSys[Nneut] = sys;	
	Nneut++;
      }
    }

    Nstart = 0;
    NMax = SCHits.size();
    if (NMax>30){
      NMax = 30;
    }
    for (int k=0; k<NMax; k++){
      const DSCHit* hit = SCHits[k];
      SC_E[Nstart] = hit->dE;
      SC_T[Nstart] = hit->t;
      SC_S[Nstart] = hit->sector;
      Nstart++;
    }


    t3->Fill();
    ECounter++;
    
    japp->RootFillUnLock(this);

  }

  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_pi0::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_pi0::fini(void)
{
  // Called before program exit after event processing is finished.
  
  WriteHistograms();

  return NOERROR;
}
jerror_t JEventProcessor_pi0::WriteHistograms(void) {
  
  TDirectory *top = gDirectory;

  F1->cd();
  F1->cd("TreeDir");

  t3->Write();

  rfbunch->Write();
  rfbeam->Write();
  rfbeamNoOffset->Write();
  rfTOF->Write();
  rfTOFAfter->Write();
  rfbeamH->Write();
  rfbeamM->Write();
  rfbeamHNOBCAL->Write();
  rfbeamMNOBCAL->Write();
  eventtype->Write();
  mictime->Write();
  hodtime->Write();
  BeamPhotonE->Write();
  BeamPhotonT->Write();
  FoundBeamPhotonE->Write();
  FoundBeamPhotonT->Write();
  FoundBeamPhotons->Write();
  sctimeraw->Write();
  scenergyraw->Write();
  twophotonm->Write(); 
  twophotonmextended->Write(); 
  PhotonTBcal->Write(); 
  PhotonTFcal->Write(); 
  FoundPhotonTBcal->Write(); 
  FoundPhotonTFcal->Write(); 
  FoundPhotonEBcal->Write(); 
  FoundPhotonEFcal->Write(); 
  PE1vsPE2->Write();
  
  MissingP->Write();
  MissingE->Write();
  MissingM->Write();
  MissingPperp->Write();
  
  PHI->Write();
  THETA->Write();
  DPHI->Write();
  
  tSC_vs_tgamma->Write();
  ESC_vs_Mgamma->Write();
  
  
  for (int k=0;k<30;k++){
    SCE[k]->Write();
    SCT[k]->Write();
  }

  F1->Close();
  top->cd();

  return NOERROR;
 
}

jerror_t JEventProcessor_pi0::CreateHistogramsAndTree(void) {

  TDirectory *top = gDirectory;

  sprintf(RootFile,"pi0notracks.root");

  F1 = new TFile(RootFile,"RECREATE");
  F1->cd();
  F1->mkdir("TreeDir");
  F1->cd("TreeDir");

  rfbunch = new TH1D("rfbunch","RF time",1000, -100., 100.);
  rfbeam = new TH1D("rfbeam","BeamPhoton time minus RF time",1000, -100., 100.);
  rfbeamNoOffset = new TH1D("rfbeamNoOffset","BeamPhoton NoOffset time minus RF time",1000, -100., 100.);
  
  rfTOF = new TH1D("rfTOF","RF time from TOF",5000, 0., 1000.);
  rfTOFAfter = new TH1D("rfTOFAfter","RF time from TOF After cut",5000, 0., 1000.);
  
  rfbeamH = new TH1D("rfbeamH","BeamPhoton time TAGH minus RF time",4000, -150., 150.);
  rfbeamM = new TH1D("rfbeamM","BeamPhoton time TAGM minus RF time",4000, -150., 150.);

  rfbeamHNOBCAL = new TH1D("rfbeamHNOBCAL","BeamPhoton time TAGH minus RF time NO BCAL HITS",4000, -150., 150.);
  rfbeamMNOBCAL = new TH1D("rfbeamMNOBCAL","BeamPhoton time TAGM minus RF time NO BCAL HITS",4000, -150., 150.);

  eventtype = new TH1D("eventtype", "1=BcalHits, 2=No BcalHits",4,0,4);
  
  mictime = new TH2D("mictime","BeamPhoton time TAGM minus RF time for all sectors",2000, -100., 100., 160., 0.,160.);
  hodtime = new TH2D("hodtime","BeamPhoton time TAGH minus RF time for all detectors",2000, -100., 100., 400., 0.,400.);

  sctimeraw = new TH1D("sctimeraw", "Start Counter timing raw for all multi photon events",800, -100., 100.);
  scenergyraw = new TH1D("scenergyraw", "Start Counter Energy raw for all multi photon events",100, 0., 0.01);

  BeamPhotonE = new TH1D("BeamPhotonE", "Photon Beam Energy", 300,0., 12.);
  BeamPhotonT = new TH1D("BeamPhotonT", "Photon Beam Time", 1000,-100., 100.);
  FoundBeamPhotonE = new TH1D("FoundBeamPhotonE", "Photon Beam Energy", 300,0., 12.);
  FoundBeamPhotonT = new TH1D("FoundBeamPhotonT", "Photon Beam Time", 1000,-100., 100.);
  FoundBeamPhotons = new TH1D("FoundBeamPhotons", "Number of found beam photons", 11, -0.5, 10.5);
  
  PhotonTBcal =  new TH1D("PhotonTBCAL", "BCAL Photon Time", 1000,-100., 100.);
  PhotonTFcal =  new TH1D("PhotonTFCAL", "FCAL Photon Time", 1000,-100., 100.);
  FoundPhotonTBcal =  new TH1D("FoundPhotonTBCAL", "BCAL Photon Time good RF", 1000,-100., 100.);
  FoundPhotonTFcal =  new TH1D("FoundPhotonTFCAL", "FCAL Photon Time good RF", 1000,-100., 100.);
  FoundPhotonEBcal =  new TH1D("FoundPhotonEBCAL", "BCAL Photon Energy good RF", 300, 0., 12.);
  FoundPhotonEFcal =  new TH1D("FoundPhotonEFCAL", "FCAL Photon Energy good RF", 300, 0., 12.);
  
  MissingP = new TH1D("MissingP", "Missing Momentum", 200, -3., 6.);
  MissingE = new TH1D("MissingE", "Missing Energy", 200, -3., 6.);
  MissingM = new TH1D("MissingM", "Missing Mass", 100, -3., 3.);
  MissingPperp = new TH1D("MissingPperp", "Missing Pt", 100, 0., 1.);
  
  PHI = new TH1D("PHI","phi angle of 2photon system",100,0.,3.1415926*2);
  DPHI = new TH1D("DPHI","dphi angle of 2photon minus recoil SC HIT",1000,0.,3.1415926*2);
  THETA = new TH1D("THETA","theta angle of 2photon system", 100, 0., 0.25);
  
  PE1vsPE2 = new TH2D("PE1vsPE2", "Photon1 Energy vs Photon2 energy (1>2)", 500, 0., 5., 500, 0., 12.);
  
  char htit[128];
  char hnam[128];
  for (int k=0;k<30;k++){
    sprintf(hnam,"SCE%d",k);
    sprintf(htit,"Start Counter Energy");
    SCE[k] = new TH1D(hnam,htit, 100., 0., 0.010);
    
    sprintf(hnam,"SCT%d",k);
    sprintf(htit,"Start Counter Time");
    SCT[k] = new TH1D(hnam,htit, 1000., -100., 100.);
  }
  twophotonm = new TH1D("twophotonm", "two photon invariant Mass", 1000, 0. , 4.);
  twophotonmextended = new TH1D("twophotonmextended", "two photon invariant Mass", 1000, 0. , 4.);

  tSC_vs_tgamma = new TH2D("tSC_vs_tgamma", "SC time vs pi0 time", 200, -40., 40., 200, -40., 40.);
  ESC_vs_Mgamma = new TH2D("ESC_vs_Mgamma", "SC Energy vs 2gamma invariant Mass", 800, 0., 4., 100, 0., 0.01 );
  
  // prepare the tree
  t3 = new TTree("t3", "All Neutral Events");
  
  t3->Branch("EventNUM",&EventNUM,"EventNUM/I"); // lower case l for ULong64_t
  t3->Branch("RFT", &RFT, "RFT/F");
  t3->Branch("Nbeam", &Nbeam, "Nbeam/I");
  t3->Branch("BeamE",BeamE,"BeamE[Nbeam]/F");
  t3->Branch("BeamT",BeamT,"BeamT[Nbeam]/F");
  t3->Branch("Nneut", &Nneut, "Nneut/I");
  t3->Branch("NeutralX", NeutralX, "NeutralX[Nneut]/F");
  t3->Branch("NeutralY", NeutralY, "NeutralY[Nneut]/F");
  t3->Branch("NeutralZ", NeutralZ, "NeutralZ[Nneut]/F");
  t3->Branch("NeutralT", NeutralT, "NeutralT[Nneut]/F");
  t3->Branch("NeutralE", NeutralE, "NeutralE[Nneut]/F");
  t3->Branch("DetSys", DetSys, "DetSys[Nneut]/I");
  t3->Branch("Nstart", &Nstart, "Nstart/I");
  t3->Branch("SC_E", SC_E, "SC_E[Nstart]/F");
  t3->Branch("SC_T", SC_T, "SC_T[Nstart]/F");
  t3->Branch("SC_S", SC_S, "SC_S[Nstart]/I");

  top->cd();

  return NOERROR;

}