// $Id$
//
//    File: JEventProcessor_AntiProtonStudies.cc
// Created: Mon Jul  8 13:08:29 EDT 2019
// Creator: staylor (on Linux ifarm1801 3.10.0-327.el7.x86_64 x86_64)
//

#include "JEventProcessor_AntiProtonStudies.h"
using namespace jana;

#define LAMBDA_CUT(mass) (fabs(mass-ParticleMass(Lambda))<0.012)

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


//------------------
// JEventProcessor_AntiProtonStudies (Constructor)
//------------------
JEventProcessor_AntiProtonStudies::JEventProcessor_AntiProtonStudies()
{

}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_AntiProtonStudies::init(void)
{
  // This is called once at program startup.
  japp->CreateLock("myAntiPlock");

  gDirectory->mkdir("LambdaLambdaBarAnalysis")->cd();
  LambdaLambdaBarConfidenceLevel=new TH1F("LambdaLambdaBarConfidenceLevel","CL",100,0,1);
  LambdaLambdaBarChiSq=new TH1F("LambdaLambdaBarChisSq","#chi^{2}/ndf",1000,0,100);
  AntiPPipVsPPim=new TH2F("AntiPPipVSPPim","M(p-bar#pi+) vs M(p#pi-)",1000,1.,3,1000,1,3);
  
  gDirectory->mkdir("PID")->cd();
 
  AntiPTOFdt=new TH2F("AntiPTOFdt","#deltat(TOF)",100,0,10,200,-5,5);	
  AntiPFCALdt=new TH2F("AntiPFCALdt","#deltat(FCAL)",100,0,10,200,-5,5);
  AntiPBCALdt=new TH2F("AntiPBCALdt","#deltat(BCAL)",100,0,10,200,-5,5);
  AntiPdEdxCDC=new TH2F("AntiPdEdxCDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0);	
  AntiPdEdxFDC=new TH2F("AntiPdEdxFDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0); 
  AntiPEOverPFCAL=new TH2F("AntiPEOverPFCAL","E/p vs p",100,0,10,150,0,1.5);
  AntiPEOverPBCAL=new TH2F("AntiPEOverPBCAL","E/p vs p",100,0,10,150,0,1.5); 

  PTOFdt=new TH2F("PTOFdt","#deltat(TOF)",100,0,10,200,-5,5);	
  PFCALdt=new TH2F("PFCALdt","#deltat(FCAL)",100,0,10,200,-5,5);
  PBCALdt=new TH2F("PBCALdt","#deltat(BCAL)",100,0,10,200,-5,5);
  PdEdxCDC=new TH2F("PdEdxCDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0);	
  PdEdxFDC=new TH2F("PdEdxFDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0); 
  PEOverPFCAL=new TH2F("PEOverPFCAL","E/p vs p",100,0,10,150,0,1.5);
  PEOverPBCAL=new TH2F("PEOverPBCAL","E/p vs p",100,0,10,150,0,1.5);

  
  gDirectory->mkdir("purity")->cd();
  
  AntiPDenom=new TH2F("AntiPDenom","P vs #theta",100,0,100,100,0.05,10.05); 
  AntiPNumer_Km=new TH2F("AntiPNumer_Km","P vs #theta",100,0,100,100,0.05,10.05);
  AntiPNumer_Pim=new TH2F("AntiPNumer_Pim","P vs #theta",100,0,100,100,0.05,10.05);
  AntiPNumer_aP=new TH2F("AntiPNumer_aP","P vs #theta",100,0,100,100,0.05,10.05);
  AntiPNumer_Em=new TH2F("AntiPNumer_Em","P vs #theta",100,0,100,100,0.05,10.05);
  
  PDenom=new TH2F("PDenom","P vs #theta",80,0,80,100,0.05,10.05); 
  PNumer_Kp=new TH2F("PNumer_Kp","P vs #theta",80,0,80,100,0.05,10.05);
  PNumer_Pip=new TH2F("PNumer_Pip","P vs #theta",80,0,80,100,0.05,10.05);
  PNumer_P=new TH2F("PNumer_P","P vs #theta",80,0,80,100,0.05,10.05);
  PNumer_Ep=new TH2F("PNumer_Ep","P vs #theta",80,0,80,100,0.05,10.05);

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

  CL_CUT=0.01;
  gPARMS->SetDefaultParameter("ANTIPROTONSTUDIES:CL_CUT",CL_CUT);

  
  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_AntiProtonStudies::brun(JEventLoop *eventLoop, int32_t runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_AntiProtonStudies::evnt(JEventLoop *loop, uint64_t eventnumber)
{
  vector<const DChargedTrack*>ctracks;
  loop->Get(ctracks); 
  if (ctracks.size()!=5) return VALUE_OUT_OF_RANGE;

  vector<const DParticleID *> pid_algorithms;
  loop->Get(pid_algorithms);
  if(pid_algorithms.size()<1){
    _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
    return RESOURCE_UNAVAILABLE;
  }
  const DParticleID *pid_algorithm = pid_algorithms[0];
  
  vector<const DBeamPhoton*>beamphotons;
  loop->Get(beamphotons);
  
  vector<const DMCThrown *>throwns;
  loop->Get(throwns);

  japp->WriteLock("myAntiPlock");

  vector<unsigned int>neglist;
  vector<unsigned int>poslist;
  for (int i=0;i<5;i++){
    if (ctracks[i]->Get_BestTrackingFOM()->charge()<0){
      neglist.push_back(i);
    }
    else poslist.push_back(i);
  }
  if (neglist.size()==2){     
    DAnalysisUtilities *dAnalysisUtilities=new DAnalysisUtilities(loop);
    DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop);
    DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils);   
    dKinFitUtils->Reset_NewEvent();
    dKinFitter->Reset_NewEvent();
    
    unsigned int posindexes[6][3];
    unsigned int negindexes[2][2];
       
    posindexes[0][0]=poslist[0];
    posindexes[0][1]=poslist[1];
    posindexes[0][2]=poslist[2]; 
    posindexes[1][0]=poslist[0];
    posindexes[1][1]=poslist[2];
    posindexes[1][2]=poslist[1];
    posindexes[2][0]=poslist[1];
    posindexes[2][1]=poslist[2];
    posindexes[2][2]=poslist[0];
    posindexes[3][0]=poslist[1];
    posindexes[3][1]=poslist[0];
    posindexes[3][2]=poslist[2]; 
    posindexes[4][0]=poslist[2];
    posindexes[4][1]=poslist[0];
    posindexes[4][2]=poslist[1];
    posindexes[5][0]=poslist[2];
    posindexes[5][1]=poslist[1];
    posindexes[5][2]=poslist[0];
    
    negindexes[0][0]=neglist[0];
    negindexes[0][1]=neglist[1];  
    negindexes[1][0]=neglist[1];
    negindexes[1][1]=neglist[0]; 
  
    for (unsigned int i=0;i<beamphotons.size();i++){    
      for (int j=0;j<2;j++){
	const DChargedTrackHypothesis *piminus_hyp=ctracks[negindexes[j][0]]->Get_Hypothesis(PiMinus);
	if (piminus_hyp==NULL) continue;
	const DChargedTrackHypothesis *antiproton_hyp=ctracks[negindexes[j][1]]->Get_Hypothesis(AntiProton);
	if (antiproton_hyp==NULL) continue;
	
	const DTrackTimeBased *piminus_track=piminus_hyp->Get_TrackTimeBased();
	if (piminus_track==NULL) continue;
	const DTrackTimeBased *antiproton_track=antiproton_hyp->Get_TrackTimeBased();
	if (antiproton_hyp==NULL) continue;
	
	for (int k=0;k<6;k++){
	  const DChargedTrackHypothesis *piplus_hyp=ctracks[posindexes[k][0]]->Get_Hypothesis(PiPlus);
	  if (piplus_hyp==NULL) continue;
	  const DChargedTrackHypothesis *proton_hyp1=ctracks[posindexes[k][1]]->Get_Hypothesis(Proton);
	  if (proton_hyp1==NULL) continue;
	  const DChargedTrackHypothesis *proton_hyp2=ctracks[posindexes[k][2]]->Get_Hypothesis(Proton);
	  if (proton_hyp2==NULL) continue;
	   
	  const DTrackTimeBased *piplus_track=piplus_hyp->Get_TrackTimeBased();
	  if (piplus_track==NULL) continue;
	  const DTrackTimeBased *proton_track1=proton_hyp1->Get_TrackTimeBased();
	  if (proton_track1==NULL) continue;
	  const DTrackTimeBased *proton_track2=proton_hyp2->Get_TrackTimeBased();
	  if (proton_track2==NULL) continue;
	  
	  double t0_rf=proton_hyp1->t0();
	  double dt_st=beamphotons[i]->time()-t0_rf; 
	  
	  double weight=1.;
	  bool got_beam_photon=false;
	  if (fabs(dt_st)>6.&&fabs(dt_st)<18.){
	    weight=-1./6.;
	    got_beam_photon=true;
	   }
	  if (fabs(dt_st)<2.){
	    got_beam_photon=true;
	  }
	  if (got_beam_photon==false) continue;

	  //--------------------------------
	  // Kinematic fit
	  //-------------------------------- 
	  dKinFitter->Reset_NewFit();
      
	  set<shared_ptr<DKinFitParticle>> InitialParticles, FinalParticles;
	  
	  shared_ptr<DKinFitParticle>myBeam=dKinFitUtils->Make_BeamParticle(beamphotons[i]);
	  shared_ptr<DKinFitParticle>myTarget=dKinFitUtils->Make_TargetParticle(Proton);
	  
	  InitialParticles.insert(myBeam);  
	  InitialParticles.insert(myTarget);
	   
	  shared_ptr<DKinFitParticle>myProton1=dKinFitUtils->Make_DetectedParticle(proton_track1);
	  shared_ptr<DKinFitParticle>myPiPlus=dKinFitUtils->Make_DetectedParticle(piplus_track);
	  shared_ptr<DKinFitParticle>myPiMinus=dKinFitUtils->Make_DetectedParticle(piminus_track);  
	  shared_ptr<DKinFitParticle>myProton2=dKinFitUtils->Make_DetectedParticle(proton_track2);
	  shared_ptr<DKinFitParticle>myAntiProton=dKinFitUtils->Make_DetectedParticle(antiproton_track); 
	  
	  FinalParticles.insert(myProton1);
	  FinalParticles.insert(myPiPlus);
	  FinalParticles.insert(myPiMinus);
	  FinalParticles.insert(myProton2);
	  FinalParticles.insert(myAntiProton);
	   
	  //MAKE LAMBDA
	  set<shared_ptr<DKinFitParticle>> locFromInitialState, locFromFinalState;
	  locFromFinalState.insert(myProton2);  
	  locFromFinalState.insert(myPiMinus);
	  shared_ptr<DKinFitParticle> myLambda = dKinFitUtils->Make_DecayingParticle(Lambda, locFromInitialState, locFromFinalState);
	  
	   // Decay Vertex constraint
	  set<shared_ptr<DKinFitParticle>> locFullConstrainParticles, locNoConstrainParticles;
	  locFullConstrainParticles.insert(myProton2);  
	  locFullConstrainParticles.insert(myPiMinus);
	  locNoConstrainParticles.insert(myLambda);
	  
	  DVector3 second_vertex;
	  dAnalysisUtilities->Calc_DOCAVertex(myProton2.get(),myPiMinus.get(),second_vertex);
		  
	  shared_ptr<DKinFitConstraint_Vertex> locDecayVertexConstraint = dKinFitUtils->Make_VertexConstraint(locFullConstrainParticles, locNoConstrainParticles,second_vertex);
	  dKinFitter->Add_Constraint(locDecayVertexConstraint);
		  
	  //MAKE ANTILAMBDA
	  set<shared_ptr<DKinFitParticle>> locFromInitialState2, locFromFinalState2;
	  locFromFinalState2.insert(myAntiProton);  
	  locFromFinalState2.insert(myPiPlus);
	  shared_ptr<DKinFitParticle> myAntiLambda = dKinFitUtils->Make_DecayingParticle(AntiLambda, locFromInitialState2, locFromFinalState2);
	  
	  // Decay Vertex constraint
	  set<shared_ptr<DKinFitParticle>> locFullConstrainParticles2, locNoConstrainParticles2;
	  locFullConstrainParticles2.insert(myAntiProton);  
	  locFullConstrainParticles2.insert(myPiPlus);
	  locNoConstrainParticles2.insert(myAntiLambda);
	  
	  DVector3 second_vertex2;
	  dAnalysisUtilities->Calc_DOCAVertex(myAntiProton.get(),myPiPlus.get(),second_vertex2);
		  
	  shared_ptr<DKinFitConstraint_Vertex> locDecayVertexConstraint2 = dKinFitUtils->Make_VertexConstraint(locFullConstrainParticles2, locNoConstrainParticles2,second_vertex2);
	  dKinFitter->Add_Constraint(locDecayVertexConstraint2);
	  
	  // make constraint
	  shared_ptr<DKinFitConstraint_P4>locP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles);
	  // set constraint
	  dKinFitter->Add_Constraint(locP4Constraint);
	  
	  // PERFORM THE KINEMATIC FIT
	  dKinFitter->Fit_Reaction();
	  
	  // GET THE FIT RESULTS
	  double locConfidenceLevel = dKinFitter->Get_ConfidenceLevel();
	  LambdaLambdaBarConfidenceLevel->Fill(locConfidenceLevel);
	  LambdaLambdaBarChiSq->Fill(dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF()));
      
	  if (locConfidenceLevel>CL_CUT){
	    if (LambdaLambdaBarAnalysis(dKinFitter,
					ctracks[posindexes[k][0]] /* pi+ guess */,
					ctracks[posindexes[k][1]] /* proton guess */,
					ctracks[posindexes[k][2]] /* proton guess */,
					ctracks[negindexes[j][0]], // pi-
					ctracks[negindexes[j][1]],weight)){
	      FillPIDHistos(t0_rf,proton_hyp1,proton_track1,proton_hyp2,
			    proton_track2,piplus_hyp,piplus_track,piminus_hyp,
			    piminus_track,antiproton_hyp,antiproton_track,
			    pid_algorithm);

	    }
	  }
	  
	}
      }
    }

    delete dKinFitter;
    delete dKinFitUtils;
    delete dAnalysisUtilities;
  }
	
  japp->Unlock("myAntiPlock");
  
  return NOERROR;
}

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


bool JEventProcessor_AntiProtonStudies::LambdaLambdaBarAnalysis(DKinFitter *dKinFitter,
								const DChargedTrack *piplus_hyp,
								const DChargedTrack *proton_hyp1, 
								const DChargedTrack *proton_hyp2,
								
								const DChargedTrack *piminus_hyp,
								const DChargedTrack *antiproton_hyp,				   
								double weight){
  DLorentzVector beam_kf;
  DLorentzVector pim_kf,pip_kf,antiproton_kf;
  vector<DLorentzVector>protons_kf;
  DVector3 piminus_pos,piplus_pos;
  vector<DVector3>proton_pos;
  
  set<shared_ptr<DKinFitParticle>>myParticles=dKinFitter->Get_KinFitParticles();
  set<shared_ptr<DKinFitParticle>>::iterator locParticleIterator=myParticles.begin();
  for(; locParticleIterator != myParticles.end(); ++locParticleIterator){
    if ((*locParticleIterator)->Get_KinFitParticleType()==d_BeamParticle){
      beam_kf=(*locParticleIterator)->Get_P4();
    }
    else if ((*locParticleIterator)->Get_KinFitParticleType()==d_DetectedParticle){
      switch ((*locParticleIterator)->Get_PID()){
      case 2212:  // Proton
	protons_kf.push_back((*locParticleIterator)->Get_P4());
	proton_pos.push_back((*locParticleIterator)->Get_Position());
	break;
      case 211: // pi+
	pip_kf=(*locParticleIterator)->Get_P4();
	piplus_pos=(*locParticleIterator)->Get_Position();
	break;
      case -211: // pi-
	pim_kf=(*locParticleIterator)->Get_P4();
	piminus_pos=(*locParticleIterator)->Get_Position();
	break;
      case -2212: // antiproton
	antiproton_kf=(*locParticleIterator)->Get_P4();
	break;
      default:
	break;
      }
    }
  }
  DLorentzVector antip_pip_kf=antiproton_kf+pip_kf;
  DLorentzVector ppim1_kf=protons_kf[0]+pim_kf; 
  DLorentzVector ppim2_kf=protons_kf[1]+pim_kf;
  
  AntiPPipVsPPim->Fill(ppim1_kf.M(),antip_pip_kf.M(),weight); 
  AntiPPipVsPPim->Fill(ppim2_kf.M(),antip_pip_kf.M(),weight);
  if (weight>0 && LAMBDA_CUT(antip_pip_kf.M())
      && (LAMBDA_CUT(ppim1_kf.M()) || LAMBDA_CUT(ppim2_kf.M()))){
    
    double antiproton_theta=antiproton_kf.Theta()*180./M_PI;
    double antiproton_momentum=antiproton_kf.P();
    AntiPDenom->Fill(antiproton_theta,antiproton_momentum);

    const DTrackTimeBased *best_fom_track=antiproton_hyp->Get_BestFOM()->Get_TrackTimeBased();
    switch(best_fom_track->PID()){
    case PiMinus:
      AntiPNumer_Pim->Fill(antiproton_theta,antiproton_momentum);
      break;
    case KMinus:
      AntiPNumer_Km->Fill(antiproton_theta,antiproton_momentum);
      break;
    case AntiProton:
      AntiPNumer_aP->Fill(antiproton_theta,antiproton_momentum);
      break;
    case Electron:
      AntiPNumer_Em->Fill(antiproton_theta,antiproton_momentum);
      break;
    default:
      break;
    }
    
    for (int i=0;i<2;i++){
      double proton_theta=protons_kf[i].Theta()*180./M_PI;
      double proton_momentum=protons_kf[i].P();
      PDenom->Fill(proton_theta,proton_momentum);
      if (i==0){
	best_fom_track=proton_hyp1->Get_BestFOM()->Get_TrackTimeBased();
      }
      else{
	best_fom_track=proton_hyp2->Get_BestFOM()->Get_TrackTimeBased();
      }
      switch(best_fom_track->PID()){
      case PiPlus:
	PNumer_Pip->Fill(proton_theta,proton_momentum);
	break;
      case KPlus:
	PNumer_Kp->Fill(proton_theta,proton_momentum);
	break;
      case Proton:
	PNumer_P->Fill(proton_theta,proton_momentum);
	break;
      case Positron:
	PNumer_Ep->Fill(proton_theta,proton_momentum);
	break;
      default:
	break;
      }
    }


    return true;
  }

  return false;
}

void JEventProcessor_AntiProtonStudies::FillPIDHistos(double t0_rf,const DChargedTrackHypothesis *proton_hyp1,
						      const DTrackTimeBased *proton_track1,
						      const DChargedTrackHypothesis *proton_hyp2,
						      const DTrackTimeBased *proton_track2,  
						      const DChargedTrackHypothesis *piplus_hyp,
						      const DTrackTimeBased *piplus_track, 
						      const DChargedTrackHypothesis *piminus_hyp,
						      const DTrackTimeBased *piminus_track, 
						      const DChargedTrackHypothesis *antiproton_hyp, 
						      const DTrackTimeBased *antiproton_track,
						      const DParticleID *pid_algorithm){  
  double p=antiproton_track->momentum().Mag();
  double dEdxAmp=1e6*antiproton_track->ddEdx_CDC_amp;	 
  double dEdx_fdc=1e6*antiproton_track->ddEdx_FDC;
  if (dEdxAmp>0){
    AntiPdEdxCDC->Fill(p/ParticleMass(Proton),dEdxAmp);
  }
  if (dEdx_fdc>0){
    AntiPdEdxFDC->Fill(p/ParticleMass(Proton),dEdx_fdc);
  }
  
  shared_ptr<const DBCALShowerMatchParams>bcalparms=antiproton_hyp->Get_BCALShowerMatchParams(); 
  shared_ptr<const DFCALShowerMatchParams>fcalparms=antiproton_hyp->Get_FCALShowerMatchParams();	  
  shared_ptr<const DTOFHitMatchParams>tofparms=antiproton_hyp->Get_TOFHitMatchParams();
	
  if (bcalparms!=NULL){
    double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf;
    AntiPBCALdt->Fill(p,dt_bcal);  
  }
  if (tofparms!=NULL){ 
    double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf;
    AntiPTOFdt->Fill(p,dt);
  } 
  if (fcalparms!=NULL){
    double dt_fcal=fcalparms->dFCALShower->getTime()
      -fcalparms->dFlightTime-t0_rf;
    AntiPFCALdt->Fill(p,dt_fcal);  

    double E_over_p=fcalparms->dFCALShower->getEnergy()/p;
    AntiPEOverPFCAL->Fill(p,E_over_p);

  }

  p=proton_track1->momentum().Mag();
  dEdxAmp=1e6*proton_track1->ddEdx_CDC_amp;	 
  dEdx_fdc=1e6*proton_track1->ddEdx_FDC;
  if (dEdxAmp>0){
    PdEdxCDC->Fill(p/ParticleMass(Proton),dEdxAmp);
  }
  if (dEdx_fdc>0){
    PdEdxFDC->Fill(p/ParticleMass(Proton),dEdx_fdc);
  }
  
  bcalparms=proton_hyp1->Get_BCALShowerMatchParams(); 
  fcalparms=proton_hyp1->Get_FCALShowerMatchParams();	  
  tofparms=proton_hyp1->Get_TOFHitMatchParams();
	
  if (bcalparms!=NULL){
    double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf;
    PBCALdt->Fill(p,dt_bcal);  
  }
  if (tofparms!=NULL){ 
    double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf;
    PTOFdt->Fill(p,dt);
  } 
  if (fcalparms!=NULL){
    double dt_fcal=fcalparms->dFCALShower->getTime()
      -fcalparms->dFlightTime-t0_rf;
    PFCALdt->Fill(p,dt_fcal);  

    double E_over_p=fcalparms->dFCALShower->getEnergy()/p;
    PEOverPFCAL->Fill(p,E_over_p);

  }
  
  
  p=proton_track2->momentum().Mag();
  dEdxAmp=1e6*proton_track2->ddEdx_CDC_amp;	 
  dEdx_fdc=1e6*proton_track2->ddEdx_FDC;
  if (dEdxAmp>0){
    PdEdxCDC->Fill(p/ParticleMass(Proton),dEdxAmp);
  }
  if (dEdx_fdc>0){
    PdEdxFDC->Fill(p/ParticleMass(Proton),dEdx_fdc);
  }
  
  bcalparms=proton_hyp2->Get_BCALShowerMatchParams(); 
  fcalparms=proton_hyp2->Get_FCALShowerMatchParams();	  
  tofparms=proton_hyp2->Get_TOFHitMatchParams();
	
  if (bcalparms!=NULL){
    double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf;
    PBCALdt->Fill(p,dt_bcal);  
  }
  if (tofparms!=NULL){ 
    double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf;
    PTOFdt->Fill(p,dt);
  } 
  if (fcalparms!=NULL){
    double dt_fcal=fcalparms->dFCALShower->getTime()
      -fcalparms->dFlightTime-t0_rf;
    PFCALdt->Fill(p,dt_fcal);  

    double E_over_p=fcalparms->dFCALShower->getEnergy()/p;
    PEOverPFCAL->Fill(p,E_over_p);

  }

}