// $Id$
//
//    File: JEventProcessor_tagg.cc
// Created: Fri Mar 31 16:34:30 EDT 2017
// Creator: zihlmann (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64)
//

#include "JEventProcessor_tagg.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_tagg());
  }
} // "C"


//------------------
// JEventProcessor_tagg (Constructor)
//------------------
JEventProcessor_tagg::JEventProcessor_tagg()
{

}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_tagg::init(void)
{
  // This is called once at program startup. 
  
  // ADC timing offsets for Microscope counters based on run 31054
  double a[102] = { 0., 0., 1., -3., 0., 0., -3, 4., 5., -3., 
		    0., 0., 0.,  0., 0., -3., 0., 0., -2., 0., 
		    0., 1., 2., 2., 0., 3., 3., 0., 1., 3. ,
		    4., 1., 4., 0., 5., 5., -3., 1., 1., -4., 
		    -3.,-3., -4., -3., 0., 0., -4., 0., 0., 0., 
		    1., -4., -4., -4., -4., -3., -3., -4., -4., -4., 
		    -4., -4., -3., 0., 0., 0., -4., -4., 0., -4., 
		    -4.,-4.,-4.,-4.,-4.,-4.,-4.,-4.,-4.,-4., 
		    -4.,-3.,-3.,0.,0.,0.,0., -4., 0., 0., 
		    -2., 0., 0., -4., 0., 4., -5., -4., -4., -6., -4., -4.};
  
  for (int k=0 ; k<102; k++){ 
    TagMTShifts[k] = a[k];
  }

  double b[240] = { 8.,8.,8.,8.,8.,8.,8.,8.,8.,8.,
		 8.,8.,4.,8.,4.,8.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,4.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 -4.,0.,0.,4.,0.,0.,-4.,-4.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,5.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,4.,5.,0.,-4.,-2.,
		 -2.,-4.,-3.,-3.,-2.,-2.,-2.,-2.,-2.,-3.,
		 0.,0.,0.,-4.,-4.,-4.,-4.,-4.,-4.,-4.,
		 -4.,-4.,-4.,-4.,0.,-4.,0.,0.,0.,-4.,
		 -4.,-4.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		 0.,0.,0.,0.,0.,0.,0.,0.,-13.,-10.,
		 -10.,-10.,-13.,-13.,-13.,-13.,-13.,-13.,-13.,-13.,
		 -13.,-13.,-13.,-13.,-13.,-10.,-10.,-10.,-13.,-7.,
		 -9.,-16.,-9.,-16.,-9.,-13., -13., -13., -13., -13.
		 -10, -13.,-10.,-13.,-10.,-13.,-13.,-16.,-16.,-13.,
		 -13.,-16.,-16.,-16.-16.,-16.,-16.,-13.,-13.,0.,
		 19.,19.,16.,0.,0.,0.,0.,0.,0.,0.};

  
  for (int k=0 ; k<240; k++){ 
    TagHTShifts[k] = b[k];
  }


  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_tagg::brun(JEventLoop *eventLoop, int32_t runnumber)
{
  
  // This is called whenever the run number changes
  
  sprintf(RootFile,"tagger_data_run%d.root",runnumber);
  RF = new TFile(RootFile,"RECREATE");

  tagmPeak = new TH2D("tagmPeak"," TAGM ADC Peak",1000,0., 4096, 110, 0.,110.);
  taghPeak = new TH2D("taghPeak"," TAGH ADC Peak",1000,0., 4096, 350, 0.,350.);
  tagmIntegral = new TH2D("tagmIntegral"," TAGM ADC Integral",1000,0., 25000, 110, 0.,110.);
  taghIntegral = new TH2D("taghIntegral"," TAGH ADC Integral",1000,0., 25000, 350, 0.,350.);

  tagmPeakCol = new TH2D("tagmPeakCol"," TAGM ADC Peak Columns",1000,0., 4096, 24, 0., 24.);
  tagmPeakColCnt = new TH2D("tagmPeakColCnt"," TAGM Columns Counters", 6, 0., 6., 6, 0., 6.);

  tagHhits = new TH1D("tagHhits", "Number of Tagger Hodoscope Hits",100,0.,100.);
  tagMhits = new TH1D("tagMhits", "Number of Tagger Microscope Hits",100,0.,100.);

  tagHtime = new TH1D("tagHtime", "Tagger Hodoscope ADC time", 400.,0.,400.);
  tagMtime = new TH1D("tagMtime", "Tagger Microscope ADC time", 400.,0.,400.);

  tagHCvsT = new TH2D("tagHCvsT", "Tagger Hodoscope Counter vs time", 400, 0., 400., 350, 0., 350.);
  tagMCvsT = new TH2D("tagMCvsT", "Tagger Microscope Counter vs time", 400, 0., 400., 120, 0., 120.);

  peak1vspeak2H = new TH2D("peak1vspeak2H", "Peak1 vs Peak2 Hodoscope", 1000, 0., 4096, 1000, 0., 4096.);
  peak1vspeak2M = new TH2D("peak1vspeak2M", "Peak1 vs Peak2 Microscope", 1000, 0., 4096, 1000, 0., 4096.);
  peak1vspeak2MCol[0] = new TH2D("peak1vspeak2MCol0", "Peak1 vs Peak2 Microscope in Col 9", 1000, 0., 4096, 1000, 0., 4096.);
  peak1vspeak2MCol[1] = new TH2D("peak1vspeak2MCol1", "Peak1 vs Peak2 Microscope in Col 27", 1000, 0., 4096, 1000, 0., 4096.);
  peak1vspeak2MCol[2] = new TH2D("peak1vspeak2MCol2", "Peak1 vs Peak2 Microscope in Col 81", 1000, 0., 4096, 1000, 0., 4096.);
  peak1vspeak2MCol[3] = new TH2D("peak1vspeak2MCol3", "Peak1 vs Peak2 Microscope in Col 99", 1000, 0., 4096, 1000, 0., 4096.);

  char htit[128];
  char hnam[128];
  for (int k=1;k<127;k++){
    sprintf(hnam,"taggH_LvsR%d_%d",k,k+1);
    sprintf(htit,"TaggerH, ADC Counter %d vs ADC Counter %d",k,k+1);
    taggH_LvsR[k-1] = new TH2D(hnam, htit, 1000, 0., 4096, 1000, 0., 4096.);
  }
  for (int k=179;k<250;k++){
    sprintf(hnam,"taggH_LvsR%d_%d",k,k+1);
    sprintf(htit,"TaggerH, ADC Counter %d vs ADC Counter %d",k,k+1);
    taggH_LvsR[k-1] = new TH2D(hnam, htit, 1000, 0., 4096, 1000, 0., 4096.);
  }

  for (int k=1;k<102;k++){
    sprintf(hnam,"taggM_LvsR%d_%d",k,k+1);
    sprintf(htit,"TaggerHM, ADC Counter %d vs ADC Counter %d",k,k+1);
    taggM_LvsR[k-1] = new TH2D(hnam, htit, 1000, 0., 4096, 1000, 0., 4096.);
  }






  return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_tagg::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);


  vector <const DTAGHDigiHit*> TagH;
  vector <const DTAGMDigiHit*> TagM;
  
  loop->Get(TagH);
  loop->Get(TagM);

  //cout<<"Number of Tag H hits: "<<TagH.size()<<endl;
  tagHhits->Fill(TagH.size());
  vector <const DTAGHDigiHit*> HHitsInTime;

  for (unsigned int k=0; k<TagH.size(); k++){

    const DTAGHDigiHit *hit = TagH[k];
    double amp = (double)hit->pulse_peak - (double)hit->pedestal/(double)hit->nsamples_pedestal;
    taghPeak->Fill(amp, hit->counter_id);
    
    double intgr = (double)hit->pulse_integral - 
      (double)hit->pedestal/(double)hit->nsamples_pedestal*(double)hit->nsamples_integral;
    taghIntegral->Fill(intgr, hit->counter_id);

    double t = hit->pulse_time*0.0625 + TagHTShifts[hit->counter_id-1];
    tagHtime->Fill(t);
    tagHCvsT->Fill(t, hit->counter_id);

    if (TMath::Abs(t-93.)<5.){
	HHitsInTime.push_back(hit);
    }

    
  }
  
  tagMhits->Fill(TagM.size());

  vector <const DTAGMDigiHit*> MHitsInTime;
  vector <const DTAGMDigiHit*> MHitsInTimeCol[4];

  for (unsigned int k=0; k<TagM.size(); k++){

    const DTAGMDigiHit *hit = TagM[k];
    double amp = (double)hit->pulse_peak - (double)hit->pedestal/(double)hit->nsamples_pedestal;
    double tagmt = hit->pulse_time*0.0625 + TagMTShifts[hit->column-1];
    
    if (hit->row<1){
      tagmPeak->Fill(amp, hit->column);
      tagMtime->Fill(tagmt);
      tagMCvsT->Fill(tagmt, hit->column);

      double t = TMath::Abs(tagmt-74.5);
      if ((t<5.) && (hit->column>42)){
	MHitsInTime.push_back(hit);
      }
      
    } else {
      int col = 4;
      if (hit->column == 9){
	col = 0;
      } else if (hit->column == 27){
	col = 1;
      } else if (hit->column == 81){
	col = 2;
      } else if (hit->column == 99){
	col = 3;
      }
      if (col<4){
	double t = TMath::Abs(tagmt-74.5);
	if (t<5.){
	  MHitsInTimeCol[col].push_back(hit);
	}
	
      }
    }

    int col = 99;
    int row = hit->row;
    
    if (hit->column == 9) {
      col = 0;
    }else if (hit->column == 27) {
      col = 1;
    } else if (hit->column == 81) {
      col = 2;
    } else  if (hit->column == 99) {
      col = 3;
    } 

    if (col<4) {
      //cout<<hit->column<<"   "<<hit->row<<"   "<<col<<"   "<<row<<"   "<<amp<<endl;
      if (row>0){
	tagmPeakColCnt->Fill((double)col,(double)row);
	tagmPeakCol->Fill(amp,(double)(row + col*6));
      }
    }

    
    double intgr = (double)hit->pulse_integral - 
      (double)hit->pedestal/(double)hit->nsamples_pedestal*(double)hit->nsamples_integral;
    tagmIntegral->Fill(intgr, hit->column);
    
  }

  // looking at neighbouring hodo hits
  if (HHitsInTime.size()>1){
    for (unsigned int n=0; n<HHitsInTime.size()-1; n++){
      const DTAGHDigiHit *hit1 = HHitsInTime[n]; 
      for (unsigned int j=n; j<HHitsInTime.size(); j++){
	const DTAGHDigiHit *hit2 = HHitsInTime[j]; 
	int dcol = TMath::Abs(hit1->counter_id - hit2->counter_id);
	if (dcol == 1){
	  double amp1 = (double)hit1->pulse_peak - (double)hit1->pedestal/(double)hit1->nsamples_pedestal;
	  double amp2 = (double)hit2->pulse_peak - (double)hit2->pedestal/(double)hit2->nsamples_pedestal;
	  peak1vspeak2H->Fill(amp1,amp2);

	  int idx = hit1->counter_id;
	  if (hit2->counter_id < idx){
	    idx =  hit2->counter_id;
	    double save1 = amp1;
	    amp1 = amp2;
	    amp2 = save1;
	  }
	  taggH_LvsR[idx-1]->Fill(amp2, amp1); 
	}
      }
    }
  }

  // looking at neighbouring columns micro hits 
  if (MHitsInTime.size()>1){
    for (unsigned int n=0; n<MHitsInTime.size()-1; n++){
      const DTAGMDigiHit *hit1 = MHitsInTime[n]; 
      for (unsigned int j=n; j<MHitsInTime.size(); j++){
	const DTAGMDigiHit *hit2 = MHitsInTime[j]; 
	int dcol = TMath::Abs(hit1->column - hit2->column);
	if (dcol == 1){
	  double amp1 = (double)hit1->pulse_peak - (double)hit1->pedestal/(double)hit1->nsamples_pedestal;
	  double amp2 = (double)hit2->pulse_peak - (double)hit2->pedestal/(double)hit2->nsamples_pedestal;
	  peak1vspeak2M->Fill(amp1,amp2);

	  int idx = hit1->column;
	  if (hit2->column < idx){
	    idx =  hit2->column;
	    double save1 = amp1;
	    amp1 = amp2;
	    amp2 = save1;
	  }
	  taggM_LvsR[idx-1]->Fill(amp2, amp1); 
	}
      }
    }
  }

  // looking at neighbouring row micro hits 
  for (int s=0;s<4;s++){

    if (MHitsInTimeCol[s].size()>1){

      for (unsigned int n=0; n<MHitsInTimeCol[s].size()-1; n++){
	const DTAGMDigiHit *hit1 = MHitsInTimeCol[s][n]; 
	for (unsigned int j=n; j<MHitsInTimeCol[s].size(); j++){
	  const DTAGMDigiHit *hit2 = MHitsInTimeCol[s][j]; 
	  int drow = TMath::Abs(hit1->row - hit2->row);
	  if (drow == 1){
	    double amp1 = (double)hit1->pulse_peak - (double)hit1->pedestal/(double)hit1->nsamples_pedestal;
	    double amp2 = (double)hit2->pulse_peak - (double)hit2->pedestal/(double)hit2->nsamples_pedestal;
	    peak1vspeak2MCol[s]->Fill(amp1,amp2);
	  }
	}
      }
    }
  }




  
  return NOERROR;
}

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

  
  tagmPeak->Write();
  tagmIntegral->Write();

  tagmPeakColCnt->Write();
  tagmPeakCol->Write();

  taghPeak->Write();
  taghIntegral->Write();

  tagHCvsT->Write();
  tagMCvsT->Write();
  
  tagHhits->Write();
  tagHtime->Write();
  tagMhits->Write();
  tagMtime->Write();
  peak1vspeak2H->Write();
  peak1vspeak2M->Write();
  for (int k=0; k<4; k++){
    peak1vspeak2MCol[k]->Write();
  }

  for (int k=0;k<260;k++){
    if (taggH_LvsR[k]){
      taggH_LvsR[k]->Write(); 
    }
  }
  for (int k=0;k<102;k++){
    if (taggM_LvsR[k]){
      taggM_LvsR[k]->Write(); 
    }
  }
  


  RF->Close();

  return NOERROR;
}