// $Id$
//
//    File: JEventProcessor_TOFmon.cc
// Created: Mon Jul 18 13:20:12 EDT 2016
// Creator: zihlmann (on Linux gluon47.jlab.org 2.6.32-431.20.3.el6.x86_64 x86_64)
//

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


//------------------
// JEventProcessor_TOFmon (Constructor)
//------------------
JEventProcessor_TOFmon::JEventProcessor_TOFmon()
{

  
}

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

  
}

//------------------
// init
//------------------
jerror_t JEventProcessor_TOFmon::init(void)
{
  // This is called once at program startup. If you are creating
  // and filling historgrams in this plugin, you should lock the
  // ROOT mutex like this:
  //
  // japp->RootWriteLock();  //  ... fill historgrams or trees ...
  // japp->RootUnLock();
  //
  TDirectory *top = gDirectory;

  sprintf(RootFile,"tofmon.root");
  //  OUTF = new TFile(RootFile,"RECREATE");
  OUTF = new TFile(RootFile,"RECREATE");
  OUTF->cd();
  OUTF->mkdir("TOFmon")->cd();

  pedestals = new TH2F("pedestals", "TOF pedestals", 176, 0., 176., 100, 50., 150.);
  amplitude = new TH2F("amplitude", "TOF signal amplitude", 176, 0., 176., 1000, 0., 4096.);
  amplitudeMatch = new TH2F("amplitudeMatch", "TOF signal amplitude with TDC Hit", 176, 0., 176., 1000, 0., 4096.);
  integrals = new TH2F("integrals", "TOF signal integrals", 176, 0., 176., 1000, 0., 30000.);

  ADCtime = new TH2F("ADCtime", "Signal timing ADC",1000, 0., 800., 176, 0., 176.);
  TDCtime = new TH2F("TDCtime", "Signal timing TDC",1000, 0., 5000., 176, 0., 176.);
  TDCtimeMatch = new TH2F("TDCtimeMatch", "Signal timing TDC with ADC Match", 1000, 0., 5000., 176, 0., 176.);
  TDChits = new TH1F("TDChits", "TDC hits",176, 0., 176);

  dt_21_p1_p2 = new TH1F("dt_21_p1_p2","Time difference dt paddle 21 plane 1 plane 2", 300, -10., 10.);
  dt_21_p1_p2corr = new TH1F("dt_21_p1_p2corr","Time difference dt paddle 21 plane 1 plane 2", 300, -10., 10.);

  char hnam[128];
  char htit[128];
  for (int k=0;k<44;k++){  
    sprintf(hnam,"MTDiff%d",k);
    sprintf(htit,"MTDiff paddle %d of plane 1 with all others of plane2",k);
    MTDiff[k] = new TH2F(hnam, htit, 600, -15., 15., 44, 0., 44.);
  }

  ADC2time = 0.0625;
  TDC2time = 0.0234375;

  //TDCtOffset = 1.24043e+03;
  //ADCtOffset = 3.01637e+02;
  TDCtOffset = 405.;
  ADCtOffset = 175.;

  /*
  TDCtOffset = 3.55270e+02;
  ADCtOffset = 1.37790e+02;
  */

  for (int k=0;k<176;k++){
    sprintf(hnam,"TWalk%d",k);
    sprintf(htit,"t_TDC - T_ADC vs amplitude");
    TWalk[k] = new TH2F(hnam, htit, 1000, 0., 4096., 400, -30., 50.);    
  }
  
  for (int k=0;k<88;k++){
    sprintf(hnam,"TDiffRaw%d",k);
    sprintf(htit,"Time difference (TDC) RAW");
    DTRaw[k] = new TH1F(hnam,htit, 100, -30., 30.);    

    sprintf(hnam,"AMPvsDTRaw0%d",k);
    sprintf(htit,"Amplitude vs DTRaw plane 0, PMT %d",k);
    AMPvsDTRaw[0][k] = new TH2F(hnam, htit,100, -20,20., 1000, 0., 4096.);
    sprintf(hnam,"AMPvsDTRaw1%d",k);
    sprintf(htit,"Amplitude vs DTRaw plane 1, PMT %d",k);
    AMPvsDTRaw[1][k] = new TH2F(hnam, htit,100, -20,20., 1000, 0., 4096.);

  }

  OUTF->cd();
  OUTF->mkdir("TCORR")->cd();
  
  for (int k=0;k<6;k++){
    int SLOT = k+3;
    for (int n=0; n<32; n++){
      sprintf(hnam,"Intensities%d%02d",SLOT,n);
      sprintf(htit,"Intensity Distribution SLOT %d, Channel: %02d",SLOT,n);
      Intensities[k][n] = new TH1F(hnam,htit,1024,0.,1024.);

      sprintf(hnam,"INL%d%02d",SLOT,n);
      sprintf(htit,"Integral Non Linearity SLOT %d, Channel: %02d",SLOT,n);
      INL[k][n] = new TH1F(hnam,htit,1024,0.,1024.);

    }   
 
  }


  top->cd();

  NOBEAM = 0;
  BEAM = 0;
  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_TOFmon::brun(JEventLoop *eventLoop, int32_t runnumber)
{
  // This is called whenever the run number changes
  map<string,double> tdcshift;
  if (!eventLoop->GetCalib("/TOF/tdc_shift", tdcshift)){
    TOF_TDC_SHIFT = tdcshift["TOF_TDC_SHIFT"];
  }
  cout<<endl<<endl<<"*******************"<<endl;
  cout<<"the run number is "<<runnumber<<endl;
  cout<<"*******************"<<endl;

  TIMINGCUT = 60.;

  cout<<"Timing offsets: "<<TDCtOffset <<" / "<<ADCtOffset<<endl;

  char fnam[128];
  sprintf(fnam,"walk_fitresults_run%06d.dat",runnumber);
  struct stat buffer;
  if ((stat(fnam, &buffer) == 0)) {
    ifstream inf;
    inf.open(fnam);
    for (int k=0;k<176;k++){
      inf>>WalkPar[k][0]>>WalkPar[k][1]>>WalkPar[k][2];
    }
    inf.close();
  } else {
    for (int k=0;k<176;k++){
      WalkPar[k][0] = 0;
      WalkPar[k][1] = 0;
      WalkPar[k][2] = 0;
    }
  }

  return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_TOFmon::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->RootWriteLock();
  //  ... fill historgrams or trees ...
  // japp->RootUnLock();

  vector<const DEPICSvalue*> epicsvalues;
  loop->Get(epicsvalues); 
 
  if (0) {
    
    bool isEPICS = loop->GetJEvent().GetStatusBit(kSTATUS_EPICS_EVENT);
    if (isEPICS) {
      for(vector<const DEPICSvalue*>::const_iterator val_itr = epicsvalues.begin();
	  val_itr != epicsvalues.end(); val_itr++) {
	const DEPICSvalue* epics_val = *val_itr;
	if (epics_val->name == "IBCAD00CRCUR6") {
	  float fconv = atof(epics_val->sval.c_str());
	  if (fconv<1.) {
	    NOBEAM = 1;
	    BEAM = 0;
	  } else {
	    NOBEAM = 0;
	    BEAM = 1;
	  }
	}
      }
      
      return NOERROR;
    }
  }

  if (NOBEAM){
    return NOERROR; // only look a nobeam data	
  }


  vector< const DCAEN1290TDCHit*> CAENHits;
  loop->Get(CAENHits);
  if (CAENHits.size()<=0){
    return NOERROR;
  }

  uint32_t locROCID = CAENHits[0]->rocid;
  vector <const DCODAROCInfo*> ROCS;
  loop->Get(ROCS);
  int indx = -1;
  for ( unsigned int n=0; n<ROCS.size(); n++) {
    if (ROCS[n]->rocid == locROCID){
      indx = n;
      break;
    }
  }

  if (indx<0){
    cout<<"ERROR: NO ROCID 1 FOUND!!!"<<endl;
    return NOERROR;
  }

  // 
  // count hit intensity for integral non linearity test
  //
  for (unsigned int j=0; j<CAENHits.size(); j++){
    const DCAEN1290TDCHit *hit = CAENHits[j];
    int S = hit->slot-3;
    int C = hit->channel;
    int bin = hit->time & 0x3FF;

    float t = ((float)hit->time)*TDC2time - TDCtOffset;
    if (fabs(t)<TIMINGCUT+20.) {
      Intensities[S][C]->Fill((float)bin);
    }    
  }

  //
  // get timing shift for 6 fold ambiguity
  //
  uint64_t TriggerTime = ROCS[indx]->timestamp;
  int TriggerBIT = TriggerTime % ((uint64_t)6);
  // float TimingShift = TOF_TDC_SHIFT - (float)TriggerBIT;
  int TShift = 1 - TriggerBIT;
  if(TShift <= 0) {
    TShift += 6;
  }
  float TimingShift = 4. * (float)TShift ;

  vector <const DTOFTDCDigiHit*> TDCHits;
  vector <const DTOFDigiHit*> ADCHits;

  float MeanTime[2][44];
  int MTHit[2][44];

  memset(MTHit,0,4*2*44);

  loop->Get(TDCHits);
  loop->Get(ADCHits);
  
  for (unsigned int k=0; k<ADCHits.size(); k++){
  
    const DTOFDigiHit *hit = ADCHits[k];

    int plane = hit->plane;
    int bar = hit->bar-1;
    int side = hit->end;
    float idx = plane*88 + bar + side*44;

    int nsamples_integral = 3+25;
    int nsamples_pedestal = 4;
    float ped = hit->pedestal/(float)nsamples_pedestal;
      
    //cout<< idx<<"    "<< ped<<"  "<<hit->pulse_peak<<"     "<<hit->pulse_integral<<endl;

    float integr = (float)hit->pulse_integral - (float)nsamples_integral*ped;
    float max = 0;
    max = (float)hit->pulse_peak;// - ped;

    amplitude->Fill(idx, max);
    integrals->Fill(idx, integr);
    pedestals->Fill(idx, ped);

    ADCtime->Fill((float)hit->pulse_time*ADC2time, idx);

    for (unsigned int s=0; s<TDCHits.size(); s++){

      const DTOFTDCDigiHit *thit = TDCHits[s];
      if ((thit->plane == plane) && (thit->bar-1 == bar) && (thit->end == side)){
	// found matching TDC Hit

	float tADC = (float)hit->pulse_time*ADC2time - ADCtOffset;
	float tTDC = (float)thit->time*TDC2time - TDCtOffset;
	
	if ((fabs(tADC)<TIMINGCUT) && (fabs(tTDC)<TIMINGCUT)){
	  float dt = tTDC-tADC+TimingShift; 
	  TWalk[(int)idx]->Fill(max, dt); 

	  // found matching TDC Hit
	  amplitudeMatch->Fill(idx, max);
	  TDCtimeMatch->Fill((float)thit->time*TDC2time, idx);

	}

      }

    }

  }


  float DT21_p1 = 999.;
  float DT21_p2 = 999.;
  float DT21_p1corr = 999.;
  float DT21_p2corr = 999.;
  int DT_found = 0;

  float TDC_Time[2];
  //float ADC_Time[2];
  float ADC_Amp[2];

  for (unsigned int j=0; j<TDCHits.size(); j++){  // first loop over TDC hits
    const DTOFTDCDigiHit *hit = TDCHits[j];
    int plane = hit->plane;
    int bar = hit->bar-1;
    int side = hit->end;
    float idx = plane*88 + bar + side*44;
    float tTDC = (float)hit->time*TDC2time - TDCtOffset;	
    TDCtime->Fill((float)hit->time*TDC2time, idx);




    if (fabs(tTDC)<TIMINGCUT) { // TDC in time
      TDC_Time[side] = tTDC;
      TDChits->Fill(idx);
     
      for (unsigned int n=j+1; n<TDCHits.size(); n++){ // loop over TDC again to find maching hit on other side
	const DTOFTDCDigiHit *hit2 = TDCHits[n];
	int p2 = hit2->plane;
	int b2 = hit2->bar-1;
	int side2 = hit2->end;
	
	if ((plane == p2) && (bar == b2) && (side != side2)){
	  float tTDC2 = (float)hit2->time*TDC2time - TDCtOffset;	
	  
	  if (fabs(tTDC2)<TIMINGCUT) { // TDC in time
	    TDC_Time[side2] = tTDC2;
	    
	    float dt = (tTDC - tTDC2)/2.; // left minus right
	    if (side){
	      dt = (tTDC2 - tTDC)/2.;
	    }

	    int loc = plane*44 + bar;
	    DTRaw[loc]->Fill(dt); // histogram is called TDiffRaw%d
	    
	    for (unsigned int k=0; k<ADCHits.size(); k++){ // now loop over ADC hits to find maching ADC for both TDC hits
	      const DTOFDigiHit *hita = ADCHits[k];
	      int pa = hita->plane;
	      int ba = hita->bar-1;
	      int ea = hita->end;
	      if ((pa == plane) && (ba == bar)) { // found matching ADC

		//const Df250PulsePedestal *PT = NULL;
		//hita->GetSingle(PT);
		float max = (float)hita->pulse_peak;
		float tADC = (float)hita->pulse_time*ADC2time - ADCtOffset;

		if (fabs(tADC)<TIMINGCUT) { // make sure ADC hit has good timing
		  //ADC_Time[ea] = tADC;
		  ADC_Amp[ea] = max;

		  for (unsigned int s=k+1; s<ADCHits.size(); s++){ // now find matching ADC hit for the other side  
		    const DTOFDigiHit *hita2 = ADCHits[s];
		    int pa2 = hita2->plane;
		    int ba2 = hita2->bar-1;
		    int ea2 = hita2->end;
		    if ((pa2 == plane) && (ba2 == bar)) {
		      //const Df250PulsePedestal *PT2 = NULL;
		      //hita2->GetSingle(PT2);

		      float max2 = (float)hita2->pulse_peak;
		      float tADC2 = (float)hita2->pulse_time*ADC2time - ADCtOffset;
		    
		      if (fabs(tADC2)<TIMINGCUT) { // make sure ADC hit has good timing

			//ADC_Time[ea2] = tADC2;
			ADC_Amp[ea2] = max2;

			int idx = pa*88+ba + 44*ea;
			float amp = ADC_Amp[ea];
			float corr[2];
			corr[ea]= WalkPar[idx][1]*TMath::Power(amp,WalkPar[idx][2]) -  
			  WalkPar[idx][1]*TMath::Power(4095,WalkPar[idx][2]);

			idx =  pa2*88+ba2+ea2*44;
			amp = ADC_Amp[ea2];
			corr[ea2]= WalkPar[idx][1]*TMath::Power(amp,WalkPar[idx][2]) -  
			  WalkPar[idx][1]*TMath::Power(4095,WalkPar[idx][2]);

			float dt_walk_corr = corr[0] - corr[1];
			AMPvsDTRaw[plane][44*ea + bar]->Fill(dt-dt_walk_corr, ADC_Amp[ea]);
			AMPvsDTRaw[plane][44*ea2+ bar]->Fill(dt-dt_walk_corr, ADC_Amp[ea2]);

			//cout<<corr1<<" / "<<corr2<<endl;

			MeanTime[plane][bar] = (TDC_Time[0]-corr[0] + TDC_Time[1]-corr[1])/2.;
			MTHit[plane][bar] = 1; // this is a good hit to keep MeanTime
			
			if ((bar == 20)){
			  
			  if ( ( TMath::Abs(((TDC_Time[0]-corr[0]) - (TDC_Time[1] - corr[1]))/2.)<1.  ) && 
			       (ADC_Amp[0]>200) && (ADC_Amp[1]>200) ) {			    
			    if (plane){
			      DT21_p2 = (TDC_Time[0] + TDC_Time[1])/2.;//((TDC_Time[0]-corr1) - (TDC_Time[1] - corr2))/2.;
			      DT21_p2corr = ((TDC_Time[0]-corr[0]) + (TDC_Time[1] - corr[1]))/2.;
			      DT_found += 1;
			   
			    } else {
			      DT21_p1 = (TDC_Time[0] + TDC_Time[1])/2.;//((TDC_Time[0]-corr1) - (TDC_Time[1] - corr2))/2.;
			      DT21_p1corr = ((TDC_Time[0]-corr[0]) + (TDC_Time[1] - corr[1]))/2.;
			      DT_found += 1;
			    }
			  }
			}
		      }
		    }
		  }
		}
	      }
	    }
	  } // tdc in time

	} // other tdc in time too 

      }
    }
  }

  if (DT_found>1){
    dt_21_p1_p2->Fill(DT21_p1-DT21_p2);
    dt_21_p1_p2corr->Fill(DT21_p1corr-DT21_p2corr);
  }

  // find half paddles and match with overlaying double ended paddle
  for (unsigned int j=0; j<TDCHits.size(); j++){
    const DTOFTDCDigiHit *hit = TDCHits[j];
    int bar = hit->bar-1;
    if ( (bar<21) || (bar>22)){
      continue;
    }
    // allright we have a single ended bar
    int plane = hit->plane;
    int side = hit->end;
    int idx = plane*88 + bar + side*44;
    float tTDC = (float)hit->time*TDC2time - TDCtOffset;	

    if (fabs(tTDC)>TIMINGCUT) { // TDC in time
      continue;
    }

    int OtherPlane = 1;
    if (plane){
      OtherPlane = 0;
    }
    int PadOff = 23;
    if (side){
      PadOff = 20;
    }

    if (MTHit[OtherPlane][PadOff]){
      float TD = MeanTime[OtherPlane][PadOff] - tTDC;
      
      for (unsigned int k=0; k<ADCHits.size(); k++){ // now loop over ADC hits to find maching ADC
	const DTOFDigiHit *hita = ADCHits[k];
	int pa = hita->plane;
	int ba = hita->bar-1;
	int ea = hita->end;
	if ((pa == plane) && (ba == bar) && (ea == side)) { // found matching ADC	
	  //const Df250PulsePedestal *PT = NULL;
	  //hita->GetSingle(PT);
	  float max = (float)hita->pulse_peak;
	  float tADC = (float)hita->pulse_time*ADC2time - ADCtOffset;	
	  if (fabs(tADC)<TIMINGCUT) { // make sure ADC hit has good timing

	    float walkcorr = WalkPar[idx][1]*TMath::Power(max,WalkPar[idx][2]) -  
	      WalkPar[idx][1]*TMath::Power(4095,WalkPar[idx][2]);
	    TD += walkcorr;
	    
	    AMPvsDTRaw[plane][44*side + bar]->Fill(TD, max);
	  }
	}
      }
    }

  }

  for (int k=0; k<44; k++){

    if (MTHit[0][k]){
      for (int n=0; n<44; n++){

	if (MTHit[1][n]){
	  float mtdiff =  MeanTime[0][k] - MeanTime[1][n];
	  MTDiff[k]->Fill(mtdiff,n);
	}
      }
    }
  }
  

  return NOERROR;
}

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

  if (0) {
    for (int k=0;k<176;k++){
      if (TWalk[k]->GetEntries()>100){
	WalkFit( TWalk[k], k);
      } else {
	for (int n=0;n<6;n++) {
	  walkpar[k][n] = 0.;
	}
      }
    }
    ofstream of;
    of.open("walk_parameters.dat");
    for (int k=0;k<176;k++){
      of<<k<<" "<<walkpar[k][0]<<" "<<walkpar[k][1]<<" "<<walkpar[k][2]<<" "<<walkpar[k][3]<<" "<<walkpar[k][4]
	<<" "<<walkpar[k][5]<<endl;
    }
    of.close();
  }

  TDirectory *top = gDirectory;

  OUTF->cd();
  OUTF->cd("TOFmon");
 
  dt_21_p1_p2->Write();
  dt_21_p1_p2corr->Write();

  pedestals->Write();
  amplitude->Write();
  amplitudeMatch->Write();
  integrals->Write();
  ADCtime->Write();
  TDCtime->Write();
  TDCtimeMatch->Write();
  TDChits->Write();
  for (int k=0;k<176;k++){
    TWalk[k]->Write();
  }

  for (int k=0;k<44;k++){
    DTRaw[k]->Write();
    AMPvsDTRaw[0][k]->Write();
    AMPvsDTRaw[0][k+44]->Write();
    DTRaw[k+44]->Write();
    AMPvsDTRaw[1][k]->Write();
    AMPvsDTRaw[1][k+44]->Write();
    MTDiff[k]->Write();
  }

  OUTF->cd();
  OUTF->cd("TCORR");

  for (int k=0;k<6;k++){
    for (int n=0; n<32; n++){
      Intensities[k][n]->Write();

      double I = Intensities[k][n]->Integral(1,1024);
      if (I>0.){

	float sum = 0.;
	float RunningInt[1024];
	for (int j=0; j<1024; j++){
	  RunningInt[j] = sum; // this makes sure the first entry is zero!
	  sum += Intensities[k][n]->GetBinContent(j+1);
	}

	for (int j=0; j<1024; j++){
	  float kp = RunningInt[j]/I*1024.;
	  float cor = ((float)j) - kp;
	  INL[k][n]->Fill( ((float)j), cor);
	}
      }      
      INL[k][n]->Write();
    }   
  }

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


  return NOERROR;

}

jerror_t JEventProcessor_TOFmon::WalkFit( TH2F *h, int idx)
{

  TF1 *f1 = new TF1("f1", "[0]+[1] * ( pow(x,[2]))",20.,2000.);

  f1->SetLineColor(2);
  f1->SetParameter(0,248.);
  f1->SetParameter(1, 1.);
  f1->SetParameter(2, -0.7);
  
  f1->SetParLimits(0,-10.,20.);
  f1->SetParLimits(1, 1., 100000.);
  f1->SetParLimits(2, -1.99, -0.25);

  h->Fit(f1,"","RQ",180.,4095.);
  TF1 *thefit = h->GetFunction("f1");
  for (int k=0;k<3;k++){
    walkpar[idx][2*k] = thefit->GetParameter(k);
    walkpar[idx][2*k+1] = thefit->GetParError(k);
  }

  return NOERROR;

}