// $Id$
//
//    File: DEventProcessor_primex_recon.cc
// Created: Mon Apr  3 11:38:03 EDT 2006
// Creator: davidl (on Darwin swire-b241.jlab.org 8.4.0 powerpc)
//

#include <map>
using namespace std;

#include <PID/DBeamPhoton.h>

#include <PID/DEventRFBunch.h>

#include <RF/DRFTime_factory.h>
#include <RF/DRFTime.h>



#include "DEventProcessor_primex_recon.h"

#include <TRACKING/DTrackTimeBased.h>

#include <TLorentzVector.h>

#include <DANA/DApplication.h>
#include <FCAL/DFCALCluster.h>
#include <TRACKING/DMCThrown.h>

#include <FCAL/DFCALHit.h>
#include <CCAL/DCCALHit.h>

#include <FCAL/DFCALGeometry.h>
#include <FCAL/DFCALShower.h>


#include <FCAL/DFCALDigiHit.h>
#include <CCAL/DCCALDigiHit.h>
#include <CCAL/DCCALShower.h>

#include <BCAL/DBCALDigiHit.h>

#include <DAQ/Df250WindowRawData.h>
#include <DAQ/Df250PulseIntegral.h>
#include <DAQ/Df250PulsePedestal.h>
#include <DAQ/Df250PulseData.h>

#include <TRIGGER/DL1Trigger.h>

#include <TAGGER/DTAGHDigiHit.h>
#include <TAGGER/DTAGMDigiHit.h>

#include <START_COUNTER/DSCHit.h>

#include "DAQ/DL1Info.h"


// Routine used to create our DEventProcessor
extern "C"{
void InitPlugin(JApplication *app){


	InitJANAPlugin(app);

	app->AddProcessor(new DEventProcessor_primex_recon());

}
} // "C"

// TEst
//------------------
// init
//------------------
jerror_t DEventProcessor_primex_recon::init(void)
{
  ccal_dt = new TH2D("ccal_dt","ccal_dt", 20, 0., 10., 1000, -10, 10);

  htrig_en     = new TH1F("trig_en","trig_en", 100, 0.,10.);

  htrig_bit      = new TH1D("trig_bit","trig_bit", 32, 0.5, 32.5);
  htrig_bit_fp   = new TH1D("trig_bit_fp","trig_bit_fp", 32, 0.5, 32.5);

  hfcal_time   = new TH1F("fcal_time","fcal_time", 1000, -200, 200);

  hfcal_dt        = new TH1F("fcal_dt","fcal_dt", 1000, -200, 200);
  hfcal_st_dt     = new TH1F("fcal_st_dt","fcal_st_dt", 1000, -200, 200);

  hfcal_beam_tagm      = new TH1F("fcal_beam_tagm","fcal_beam_tagm", 2000, -100, 100);
  hfcal_beam_tagh      = new TH1F("fcal_beam_tagh","fcal_beam_tagh", 2000, -100, 100);

  htagh_fcal_beam      = new TH2D("tagh_fcal_beam","tagh_fcal_beam", 300, 0.5, 300.5, 1000, -200, 200);
  htagm_fcal_beam      = new TH2D("tagm_fcal_beam","tagm_fcal_beam", 300, 0.5, 300.5, 1000, -200, 200);

  hccal_fcal_dt  = new TH1F("ccal_fcal_dt","ccal_fcal_dt", 1000, -200, 200);

  hn_pair  = new TH1D("pair","pair", 10, -0.5, 9.5);
  hn_fcal  = new TH1D("nfcal","nfcal", 10, -0.5, 9.5);
  hn_ccal  = new TH1D("nccal","nccal", 10, -0.5, 9.5);

  comp_dphi     = new TH1D("comp_dphi","comp_dphi",  720, -360., 360);

  int nbin     =  800;
  double x_min = -5;
  double x_max =  5;

  //  int nbin     =  500;

  comp_e3    = new TH1D("comp_e3","comp_e3",    nbin, x_min, x_max);
  comp_e4    = new TH1D("comp_e4","comp_e4",    nbin, x_min, x_max);
  comp_e5    = new TH1D("comp_e5","comp_e5",    nbin, x_min, x_max);
  comp_e6    = new TH1D("comp_e6","comp_e6",    nbin, x_min, x_max);
  comp_e7    = new TH1D("comp_e7","comp_e7",    nbin, x_min, x_max);
  comp_e8    = new TH1D("comp_e8","comp_e8",    nbin, x_min, x_max);
  comp_e9    = new TH1D("comp_e9","comp_e9",    nbin, x_min, x_max);
  comp_e10   = new TH1D("comp_e10","comp_e10",  nbin, x_min, x_max);
  comp_e3p   = new TH1D("comp_e3p","comp_e3p",  nbin, x_min, x_max);
  comp_e6p   = new TH1D("comp_e6p","comp_e6p",  nbin, x_min, x_max);
  comp_e9p   = new TH1D("comp_e9p","comp_e9p",  nbin, x_min, x_max);

  comp_e6p_tmp   = new TH1D("comp_e6p_tmp","comp_e6p_tmp",   nbin, x_min, x_max);

  comp_sb_e3    = new TH1D("comp_sb_e3","comp_sb_e3",     nbin, x_min, x_max);
  comp_sb_e4    = new TH1D("comp_sb_e4","comp_sb_e4",     nbin, x_min, x_max);
  comp_sb_e5    = new TH1D("comp_sb_e5","comp_sb_e5",     nbin, x_min, x_max);
  comp_sb_e6    = new TH1D("comp_sb_e6","comp_sb_e6",     nbin, x_min, x_max);
  comp_sb_e7    = new TH1D("comp_sb_e7","comp_sb_e7",     nbin, x_min, x_max);
  comp_sb_e8    = new TH1D("comp_sb_e8","comp_sb_e8",     nbin, x_min, x_max);
  comp_sb_e9    = new TH1D("comp_sb_e9","comp_sb_e9",     nbin, x_min, x_max);
  comp_sb_e10   = new TH1D("comp_sb_e10","comp_sb_e10",   nbin, x_min, x_max);
  comp_sb_e3p   = new TH1D("comp_sb_e3p","comp_sb_e3p",   nbin, x_min, x_max);
  comp_sb_e6p   = new TH1D("comp_sb_e6p","comp_sb_e6p",   nbin, x_min, x_max);
  comp_sb_e9p   = new TH1D("comp_sb_e9p","comp_sb_e9p",   nbin, x_min, x_max);

  comp_sb_e6p_tmp   = new TH1D("comp_sb_e6p_tmp","comp_sb_e6p_tmp",   nbin, x_min, x_max);

  nbin  =  500;
  x_min = -1;
  x_max =  1;

  comp_r3    = new TH1D("comp_r3","comp_r3",    nbin, x_min, x_max);
  comp_r4    = new TH1D("comp_r4","comp_r4",    nbin, x_min, x_max);
  comp_r5    = new TH1D("comp_r5","comp_r5",    nbin, x_min, x_max);
  comp_r6    = new TH1D("comp_r6","comp_r6",    nbin, x_min, x_max);
  comp_r7    = new TH1D("comp_r7","comp_r7",    nbin, x_min, x_max);
  comp_r8    = new TH1D("comp_r8","comp_r8",    nbin, x_min, x_max);
  comp_r9    = new TH1D("comp_r9","comp_r9",    nbin, x_min, x_max);
  comp_r10   = new TH1D("comp_r10","comp_r10",  nbin, x_min, x_max);
  comp_r3p   = new TH1D("comp_r3p","comp_r3p",  nbin, x_min, x_max);
  comp_r6p   = new TH1D("comp_r6p","comp_r6p",  nbin, x_min, x_max);
  comp_r9p   = new TH1D("comp_r9p","comp_r9p",  nbin, x_min, x_max);

  comp_sb_r3    = new TH1D("comp_sb_r3","comp_sb_r3",     nbin, x_min, x_max);
  comp_sb_r4    = new TH1D("comp_sb_r4","comp_sb_r4",     nbin, x_min, x_max);
  comp_sb_r5    = new TH1D("comp_sb_r5","comp_sb_r5",     nbin, x_min, x_max);
  comp_sb_r6    = new TH1D("comp_sb_r6","comp_sb_r6",     nbin, x_min, x_max);
  comp_sb_r7    = new TH1D("comp_sb_r7","comp_sb_r7",     nbin, x_min, x_max);
  comp_sb_r8    = new TH1D("comp_sb_r8","comp_sb_r8",     nbin, x_min, x_max);
  comp_sb_r9    = new TH1D("comp_sb_r9","comp_sb_r9",     nbin, x_min, x_max);
  comp_sb_r10   = new TH1D("comp_sb_r10","comp_sb_r10",   nbin, x_min, x_max);
  comp_sb_r3p   = new TH1D("comp_sb_r3p","comp_sb_r3p",   nbin, x_min, x_max);
  comp_sb_r6p   = new TH1D("comp_sb_r6p","comp_sb_r6p",   nbin, x_min, x_max);
  comp_sb_r9p   = new TH1D("comp_sb_r9p","comp_sb_r9p",   nbin, x_min, x_max);


  ccal_xy     = new TH2F("ccal_xy","ccal_xy", 12, -12., 12., 12, -12., 12.);
  fcal_xy     = new TH2F("fcal_xy","fcal_xy", 21, -42., 42., 21, -42., 42.);


  comp_tagh     = new TH2D("comp_tagh","comp_tagh", 300, 0.5, 300.5, 200, -5, 5.);
  comp_tagh_bg  = new TH2D("comp_tagh_bg","comp_tagh_bg", 300, 0.5, 300.5, 200, -5, 5.);

  comp_tagm     = new TH2D("comp_tagm","comp_tagm", 300, 0.5, 300.5, 200, -5, 5.);
  comp_tagm_bg  = new TH2D("comp_tagm_bg","comp_tagm_bg", 300, 0.5, 300.5, 200, -5, 5.);


  TDirectory *dir = new TDirectory("MIP","MIP");
  dir->cd();
  
  // Go back up to the parent directory
  dir->cd("../");

  tree1 = new TTree("tree1","Primex");

  tree1->Branch("bit1",&bit1,"bit1/I");
  tree1->Branch("bit2",&bit2,"bit2/I");
  tree1->Branch("bit3",&bit3,"bit3/I");
  tree1->Branch("bit4",&bit4,"bit4/I");
  tree1->Branch("bit11",&bit11,"bit11/I");

  tree1->Branch("nph",    &nph,    "nph/I");
  tree1->Branch("beam_en",&beam_en,"beam_en[nph]/F");

  tree1->Branch("nccal",     &nccal,     "nccal/I");
  tree1->Branch("ccal_time", &ccal_time, "ccal_time[nccal]/F");
  tree1->Branch("ccal_en",   &ccal_en,   "ccal_en[nccal]/F");

  tree1->Branch("ccal_x",   &ccal_x,   "ccal_x[nccal]/F");
  tree1->Branch("ccal_y",   &ccal_y,   "ccal_y[nccal]/F");
  tree1->Branch("ccal_z",   &ccal_z,   "ccal_z[nccal]/F");

  tree1->Branch("ccal_dime",   &ccal_dime,   "ccal_dime[nccal]/I");

  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_primex_recon::brun(JEventLoop *eventLoop, int32_t runnumber)
{

  // This is called whenever the run number changes
  dRFTimeFactory = static_cast<DRFTime_factory*>(eventLoop->GetFactory("DRFTime"));
  vector<const DRFTime*> locRFTimes;
  eventLoop->Get(locRFTimes);



  // FCAL constants - will be retrieved from the RCDB
  

  
  return NOERROR;



}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_primex_recon::evnt(JEventLoop *loop, uint64_t eventnumber){
  
         vector<const DL1Trigger*> l1trig;	 

	 vector<const DBeamPhoton*> beam_ph;    
	
	 vector<const DFCALShower  *>fcal;   
	 vector<const DCCALShower  *>ccal;
	 
	 vector<const DFCALHit  *>fcal_hit;   
	 vector<const DCCALHit  *>ccal_hit;

	 unsigned int trig_bit[33], trig_bit_fp[33];
	

	 vector <const DRFTime*> rfTimes;
	 const DRFTime* rfTime = nullptr;
	 loop->Get(rfTimes,"TAGH");
	 if (rfTimes.size() > 0)
	   rfTime = rfTimes[0];
	 else
	   return NOERROR;
	     
	 //	 cout << " SASCHA " << rfTimes.size() << endl;

	 double t_RF = rfTime->dTime;

	 loop->Get(l1trig);

	 loop->Get(beam_ph);

	 loop->Get(fcal);
	 loop->Get(ccal);

	 loop->Get(fcal_hit);
	 loop->Get(ccal_hit);

	 japp->RootWriteLock();



	 memset(trig_bit,0,sizeof(trig_bit));
	 memset(trig_bit_fp,0,sizeof(trig_bit_fp));
	 
         memset(trig_bit,0,sizeof(trig_bit));
	 memset(trig_bit_fp,0,sizeof(trig_bit_fp));
	 memset(beam_en,0,sizeof(beam_en));

	 memset(ccal_time,0,sizeof(ccal_time));
	 memset(ccal_en,0,sizeof(ccal_en));
	 memset(ccal_x,0,sizeof(ccal_x));
	 memset(ccal_y,0,sizeof(ccal_y));
	 memset(ccal_z,0,sizeof(ccal_z));
	 memset(ccal_dime,0,sizeof(ccal_dime));

	 bit1 = 0;  bit2 = 0;  bit3 = 0;  bit4 = 0;  bit11 = 0;  

	 nph = 0;

	 nccal = 0;


	 DLorentzVector beam;





	 if( l1trig.size() > 0){
	   
	   for(unsigned int bit = 0; bit < 32; bit++){
	     trig_bit[bit + 1] = (l1trig[0]->trig_mask & (1 << bit)) ? 1 : 0;	    
	     if(trig_bit[bit + 1] == 1){
	       htrig_bit->Fill(Float_t(bit+1)); 
	     }
	   } 
	   
	   for(unsigned int bit = 0; bit < 32; bit++){
	     trig_bit_fp[bit + 1] = (l1trig[0]->fp_trig_mask & (1 << bit)) ? 1 : 0;	    
	     if(trig_bit_fp[bit + 1] == 1) htrig_bit_fp->Fill(Float_t(bit+1)); 
	   } 
	   
	   if(trig_bit[1] == 1)  bit1   = 1;
           if(trig_bit[2] == 1)  bit2   = 1;
           if(trig_bit[3] == 1)  bit3   = 1;
           if(trig_bit[4] == 1)  bit4   = 1;
           if(trig_bit[11] == 1) bit11  = 1;


	 }


	 if(trig_bit[1] == 1) {

	   double trig_en = 0;
	   for (unsigned int ii = 0; ii < fcal_hit.size(); ii++){
	     const  DTAGHDigiHit* tagh_digi_hit = NULL;
          
	     //	     tagh_hits[ii]->GetSingle(tagh_digi_hit);  
          
	     int adc_time = 0;
          
	     if(tagh_digi_hit)
	       adc_time    =  (tagh_digi_hit->pulse_time & 0x7FC0) >> 6;
	     

	     trig_en += fcal_hit[ii]->E;
	   }

	   for (unsigned int ii = 0; ii < ccal_hit.size(); ii++)
	     trig_en += ccal_hit[ii]->E;	     
	   
	   htrig_en->Fill(trig_en);
	     

	 }

	 
#if 1
	 if(beam_ph.size() > 0 ){
	   for(unsigned int ii = 0; ii < beam_ph.size(); ii++){
	     if(nph > 100) continue;
	     beam_en[nph] =  beam_ph[ii]->momentum().Mag();
	     nph++;
	   }
	 }
#endif


	 int fcal_ccal_pair = 0;

	 int fcal_id = 0;
	 int ccal_id = 0;
	 int fcal_cl = 0;
	 int ccal_cl = 0;

	 int index_fcal = -10;
	 int index_ccal = -10;


	 for (unsigned int ii = 0; ii < fcal.size(); ii++){  
	   const DFCALShower* FCALShower = fcal[ii];	     
	   
	   //	   if(FCALShower->getEnergy() < 0.1) continue;

	   if(FCALShower->getEnergy() < 0.5) continue;

	   double fcal_time    =  FCALShower->getTime();


	   for (unsigned int jj = 0; jj < ccal.size(); jj++){ 

	     if(ccal[jj]->E < 1.) continue;
	     

	     ccal_time[nccal] = ccal[jj]->time;
	     ccal_dime[nccal] = ccal[jj]->dime;
	     
	     double dt = fcal_time - ccal[jj]->time;

	     hccal_fcal_dt->Fill(dt);

	     
	     // 61914
	     //	     if( (dt > -25) && (dt < -21) && (ccal[jj]->dime > 3)){ 
	     // 61950 - 61327

	     if( (dt > -24) && (dt < -19) && (ccal[jj]->dime > 3)){
 
	       if(fcal_id != ii) fcal_cl ++;
	       if(ccal_id != jj) ccal_cl ++;
   
	       fcal_id = ii;
	       ccal_id = jj;

	       index_fcal = ii;
	       index_ccal = jj;

	       fcal_ccal_pair++;
	     }
	     
	     
	   }

	 }



	 hn_pair->Fill(float(fcal_ccal_pair));
	 hn_fcal->Fill(float(fcal_cl));
	 hn_ccal->Fill(float(ccal_cl));



	 
	 //	 if(fcal_ccal_pair == 1  && fcal.size() == 1 && ccal.size() == 1){

	 if(fcal_ccal_pair == 1 ){

	   if(beam_ph.size() > 0 ){
	     
	     double fcal_time    =  fcal[index_fcal]->getTime();
	     double fcal_energy  =  fcal[index_fcal]->getEnergy();

	     double ccal_x    =  ccal[index_ccal]->x;
	     double ccal_y    =  ccal[index_ccal]->y;
	     double phic       = atan2(ccal_y,ccal_x);
	     
	     
	     double fcal_x = fcal[index_fcal]->getPosition().X();
	     double fcal_y = fcal[index_fcal]->getPosition().Y();
	     double phif   = atan2(fcal_y,fcal_x);
	     
	     
	     int bad_accept = 0;
	     
	     if( (fcal_x > -10) && (fcal_x < 10) &&  (fcal_y > -10) && (fcal_y < 10)) bad_accept = 1;
	     if( (ccal_x > -4) && (ccal_x < 4) &&  (ccal_y > -4) && (ccal_y < 4)) bad_accept = 1;
	     if( (ccal_x < -7) ||  (ccal_x > 9)) bad_accept = 1;
	     if( (ccal_y < -9) ||  (ccal_y > 9)) bad_accept = 1;
	     
	     if(bad_accept == 0){
	       
	       
	       double dphi = (phif-phic)*180./3.1415;
	       comp_dphi->Fill(dphi);
	       
	       for(unsigned int ii = 0; ii < beam_ph.size(); ii++){
		 
		 double beam_en_tmp  = beam_ph[ii]->momentum().Mag();
		 double fcal_ccal_en = fcal_energy + ccal[index_ccal]->E - beam_en_tmp;
		 
		 int detector_id  =  beam_ph[ii]->dSystem;
		 int counter_id   =  beam_ph[ii]->dCounter;
		 
		 //	       cout << counter_id << endl;
		 
		 //		 double dt_fcal_beam =  fcal_time - beam_ph[ii]->time();
		 
		 //		 double t_ccal = dRFTimeFactory->Step_TimeToNearInputTime( ccal[index_ccal]->time, t_RF);

		 double t_beam = dRFTimeFactory->Step_TimeToNearInputTime( beam_ph[ii]->time(), t_RF);

		 //		 double dt_fcal_beam =  ccal[index_ccal]->time - beam_ph[ii]->time();

		 //		 double dt_fcal_beam =  t_ccal - t_RF;
		 
		 //		 double dt_fcal_beam =  t_beam - ccal[index_ccal]->time;

		 //		 double dt_fcal_beam =   ccal[index_ccal]->time - t_RF;

		 double dt_fcal_beam = beam_ph[ii]->time() - (t_beam - t_RF) - ccal[index_ccal]->time + 40.345;

		 //		 ccal_dt->Fill(ccal[index_ccal]->E,dt_fcal_beam);





		 double t_min = 1000.;
		 double t_max = 1000.;
		 
		 double t_min_side = 1000.;
		 double t_max_side = 1000.;
		 
		 
		 if(detector_id == 2048){
		   if(fabs(fcal_ccal_en) < 1. && beam_en_tmp > 6.) 
		     hfcal_beam_tagh->Fill(dt_fcal_beam);
		   htagh_fcal_beam->Fill(float(counter_id), dt_fcal_beam);
		   t_min = 17.;
		   t_max = 21.;
		   t_min_side = 1.;
		   t_max_side = 13.;
		   
		   
		 }
		 
		 if(detector_id == 128){
		   hfcal_beam_tagm->Fill(dt_fcal_beam);
		   htagm_fcal_beam->Fill(float(counter_id), dt_fcal_beam);
		   t_min = 17.;
		   t_max = 21.;
		   t_min_side = 1.;
		   t_max_side = 13.;
		 }
		 
		 //	       if(dt_fcal_beam > 18. && dt_fcal_beam < 21.){ 
		 if(dt_fcal_beam > t_min && dt_fcal_beam < t_max){ 
		   
		   fcal_xy->Fill(fcal_x, fcal_y);
		   ccal_xy->Fill(ccal_x, ccal_y);
		   
		   if(beam_ph[ii]->momentum().Mag() > 3){
		     comp_e3p->Fill(fcal_ccal_en);
		     comp_r3p->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if(beam_ph[ii]->momentum().Mag() > 6){
		     comp_e6p->Fill(fcal_ccal_en);
		     comp_r6p->Fill(fcal_ccal_en/beam_en_tmp);
		     if( fabs(dphi) > 130 &&  fabs(dphi) < 230) comp_e6p_tmp->Fill(fcal_ccal_en);
		   }
		   if(beam_ph[ii]->momentum().Mag() > 9){
		     comp_e9p->Fill(fcal_ccal_en);
		     comp_r9p->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   
		   if((beam_ph[ii]->momentum().Mag() > 3) && (beam_ph[ii]->momentum().Mag() <= 4)){
		     comp_e3->Fill(fcal_ccal_en);
		     comp_r3->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 4) && (beam_ph[ii]->momentum().Mag() <= 5)){
		     comp_e4->Fill(fcal_ccal_en);
		     comp_r4->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 5) && (beam_ph[ii]->momentum().Mag() <= 6)){
		     comp_e5->Fill(fcal_ccal_en);
		     comp_r5->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 6) && (beam_ph[ii]->momentum().Mag() <= 7)){
		     comp_e6->Fill(fcal_ccal_en);
		     comp_r6->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 7) && (beam_ph[ii]->momentum().Mag() <= 8)){
		     comp_e7->Fill(fcal_ccal_en);
		     comp_r7->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 8) && (beam_ph[ii]->momentum().Mag() <= 9)){
		     comp_e8->Fill(fcal_ccal_en);
		     comp_r8->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 9) && (beam_ph[ii]->momentum().Mag() <= 10)){
		     comp_e9->Fill(fcal_ccal_en);
		     comp_r9->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 10) && (beam_ph[ii]->momentum().Mag() <= 12)){
		     comp_e10->Fill(fcal_ccal_en);
		     comp_r10->Fill(fcal_ccal_en/beam_en_tmp);
		   }
		   


		   if(detector_id == 2048){
		     comp_tagh->Fill(float(counter_id),fcal_ccal_en);
		   } else if(detector_id == 128){
		     comp_tagm->Fill(float(counter_id),fcal_ccal_en);
		   }
		   
		 }
		 
		 if(dt_fcal_beam > t_min_side && dt_fcal_beam < t_max_side){ 
		   
		   double scale = 1./3.;
		   
		   if(beam_ph[ii]->momentum().Mag() > 3){
		     comp_sb_e3p->Fill(fcal_ccal_en, scale);
		     comp_sb_r3p->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if(beam_ph[ii]->momentum().Mag() > 6){
		     comp_sb_e6p->Fill(fcal_ccal_en, scale);
		     comp_sb_r6p->Fill(fcal_ccal_en/beam_en_tmp, scale);
		     if( fabs(dphi) > 130 &&  fabs(dphi) < 230) comp_sb_e6p_tmp->Fill(fcal_ccal_en,scale);
		   }
		   if(beam_ph[ii]->momentum().Mag() > 9){
		     comp_sb_e9p->Fill(fcal_ccal_en,scale);
		     comp_sb_r9p->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   
		   if((beam_ph[ii]->momentum().Mag() > 3) && (beam_ph[ii]->momentum().Mag() <= 4)){
		     comp_sb_e3->Fill(fcal_ccal_en,scale);
		     comp_sb_r3->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 4) && (beam_ph[ii]->momentum().Mag() <= 5)){
		     comp_sb_e4->Fill(fcal_ccal_en,scale);
		     comp_sb_r4->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 5) && (beam_ph[ii]->momentum().Mag() <= 6)){
		     comp_sb_e5->Fill(fcal_ccal_en,scale);
		     comp_sb_r5->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 6) && (beam_ph[ii]->momentum().Mag() <= 7)){
		     comp_sb_e6->Fill(fcal_ccal_en,scale);
		     comp_sb_r6->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 7) && (beam_ph[ii]->momentum().Mag() <= 8)){
		     comp_sb_e7->Fill(fcal_ccal_en,scale);
		     comp_sb_r7->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 8) && (beam_ph[ii]->momentum().Mag() <= 9)){
		     comp_sb_e8->Fill(fcal_ccal_en,scale);
		     comp_sb_r8->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 9) && (beam_ph[ii]->momentum().Mag() <= 10)){
		     comp_sb_e9->Fill(fcal_ccal_en,scale);
		     comp_sb_r9->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   if((beam_ph[ii]->momentum().Mag() > 10) && (beam_ph[ii]->momentum().Mag() <= 12)){
		     comp_sb_e10->Fill(fcal_ccal_en,scale);
		     comp_sb_r10->Fill(fcal_ccal_en/beam_en_tmp,scale);
		   }
		   
		   
		   if(detector_id == 2048){
		     comp_tagh_bg->Fill(float(counter_id),fcal_ccal_en);
		   } else if(detector_id == 128){
		     comp_tagm_bg->Fill(float(counter_id),fcal_ccal_en);
		   }
		   
		 }
		 

	       }
	       
	     }
	   }	   	   
	 }
	   
	 

	 double fcal_time = -1000;


#if 0
	 // FCAL reconstruction
	 
	 //	 cout << fcal.size() << endl;
	 
	 for (unsigned int jj = 0; jj < fcal.size(); jj++){  
	   const DFCALShower* dFCALShower = fcal[jj];
	   
	   double fcal_x       =  dFCALShower->getPosition().X();                 
	   double fcal_y       =  dFCALShower->getPosition().Y();
	   double fcal_z       =  dFCALShower->getPosition().Z();
	   
	   //	   cout << fcal_z << endl;
	   
	   double fcal_time    =  dFCALShower->getTime();
	   double fcal_energy  =  dFCALShower->getEnergy();
	   
	   hfcal_time->Fill(fcal_time);
	   
	   if(beam_ph.size() > 0 ){
	     for(unsigned int kk = 0; kk < beam_ph.size(); kk++){
	       if(nph > 100) continue;
	       double dt_fcal_beam =  fcal_time - beam_ph[kk]->time();
	       //	       hfcal_beam->Fill(dt_fcal_beam);
	     }
	   }
	   
	   
	 }
	 
	 
	 
#endif	 
	 

	   for (unsigned int jj = 0; jj < ccal.size(); jj++){ 


	     double ccal_time   =  ccal[jj]->time;
	     double ccal_en     =  ccal[jj]->E ;

	     double t_ccal = dRFTimeFactory->Step_TimeToNearInputTime(ccal_time, t_RF);

	     ccal_dt->Fill(ccal_en, t_ccal - t_RF);

	   }




	 
	 
	 
	 //	 tree1->Fill();
	 
	 japp->RootUnLock();  
	 
	 return NOERROR;

}
	 
//------------------
// erun
//------------------
jerror_t DEventProcessor_primex_recon::erun(void)
{
	// Any final calculations on histograms (like dividing them)
	// should be done here. This may get called more than once.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DEventProcessor_primex_recon::fini(void)
{
	return NOERROR;
}