// -----------------------------------------------------------------------------
// -----                      R3BNeutronTracker2D                          -----
// -----                Created 29-02-2012 by D.Kresan                     -----
// -----            Based on algorithm developed by M.Heil                 -----
// -----------------------------------------------------------------------------

#include <iostream>
#include <math.h>
#include <algorithm>
#include <fstream>

#include "TClonesArray.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TFile.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "R3BLandDigiPar.h"
#include "R3BLandDigi.h"
#include "R3Bveto_segDigiPar.h"
#include "R3Bveto_segDigi.h"
#include "R3BLandFirstHits.h"
#include "R3BMCTrack.h"		
#include "R3BLandPoint.h"
#include "R3BNeutHit.h"
#include "R3BNeutronTracker2D.h"

using std::cout;
using std::endl;
using std::ifstream;

ClassImp(R3BNeutronTracker2D)

Double_t fTarget_Xpos_global = 0.0;
Double_t fTarget_Ypos_global = 0.0;
Double_t fTarget_Zpos_global = 0.0;
Double_t fBeam_Xpos_global = 0.0;
Double_t fBeam_Ypos_global = 0.0;
Double_t fBeam_Zpos_global = 0.0;

Double_t c = 29.9792458;
Double_t gBeamBeta;
bool AuxSortClustersBeta(R3BNeuLandCluster*, R3BNeuLandCluster*);

// -----------------------------------------------------------------------------
R3BNeutronTracker2D::R3BNeutronTracker2D()
    : FairTask("R3B NeuLAND Neutron Tracker")
    , f2DCutEnabled(kTRUE)
    , fNNeutrons(0)
    
    // initialize global positioning functions:
    , fTarget_Xpos(0.0)
    , fTarget_Ypos(0.0)
    , fTarget_Zpos(0.0)
    , fBeam_Xpos(0.0)
    , fBeam_Ypos(0.0)
    , fBeam_Zpos(0.0)
{
    dio = 10.6; // 3 times half the diogonal of a paddle
}
// -----------------------------------------------------------------------------

// -----------------------------------------------------------------------------
R3BNeutronTracker2D::~R3BNeutronTracker2D()
{
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::SetParContainers()
{
  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  if ( ! rtdb ) Fatal("SetParContainers", "No runtime database");

  fLandDigiPar = (R3BLandDigiPar*)(rtdb->getContainer("R3BLandDigiPar"));

  if ( fLandDigiPar ) {
      cout << "-I- R3BLandDigitizer::SetParContainers() "<< endl;
      cout << "-I- Container R3BLandDigiPar  loaded " << endl;
  }
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
InitStatus R3BNeutronTracker2D::Init()
{
  // Get input array 
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No FairRootManager");
  fLandPoints = (TClonesArray*) ioman->GetObject("LandPoint");
  fLandMCTrack = (TClonesArray*) ioman->GetObject("MCTrack");
  fLandDigi = (TClonesArray*) ioman->GetObject("LandDigi");
fveto_segDigi = (TClonesArray*) ioman->GetObject("veto_segDigi");
  fLandFirstHits = (TClonesArray*) ioman->GetObject("LandFirstHits");
  fArrayCluster = (TClonesArray*) ioman->GetObject("NeuLandCluster");

  // New structure created by the Neutron Tracker
  fNeutHits = new TClonesArray("R3BNeutHit");
  ioman->Register("LandNeutHits", "Neutron Hits", fNeutHits, kTRUE);

  npaddles = fLandDigiPar->GetMaxPaddle()+1;
  nplanes = fLandDigiPar->GetMaxPlane();
  cout<<"# paddles: "<<npaddles-1<<"  # planes: "<<nplanes<<endl;
  amu = 931.494028;
  mNeutron = 1.0086649156*amu;
  cMedia = 14.;// speed of light in material in [cm/ns]
  eventNo=0;
  printing=0;
  
  PRIM_part=new R3BPrimPart*[10];
  PRIM_prot=new R3BPrimPart*[10];
  PRIM_frag=new R3BPrimPart*[10];
  PRIM_gamma=new R3BPrimPart*[10];  
     
  CreateHistograms();
  
  fCounter = -1;
  
  // Pass values to global variables:
  fTarget_Xpos_global = fTarget_Xpos;
  fTarget_Ypos_global = fTarget_Ypos;
  fTarget_Zpos_global = fTarget_Zpos;
  fBeam_Xpos_global = fBeam_Xpos;
  fBeam_Ypos_global = fBeam_Ypos;
  fBeam_Zpos_global = fBeam_Zpos;
  // NOTE: The class-variables are set with member functions AFTER the constructor, but before
  // the Init()-function is called, hence dureing Exec(), the global variables and the class-variables
  // have the same values! As it should be...

  return kSUCCESS;
}
// -----------------------------------------------------------------------------




void R3BNeutronTracker2D::ReadCalibrFiles()
{
  // Open calibration files:
  TFile* inputs = new TFile("../../Data/Inputs/Inputs.root","read");
  TH1I* Integers = (TH1I*) inputs->Get("Integers");
  fMaxProtons = Integers->GetBinContent(246);

  // Weonly need one file at a time, even that we read multiple of them.
  ifstream ThisFile;
  
  // Now read the files one-by-one:
  TString Name = "";
  TString kstr = "";
  TString st = "";
  Int_t NLines = 0;
  fNumberOfCuts = new Int_t[fMaxProtons+1];
  fCuts = new Double_t*[fMaxProtons+1];
  fKappa = new Double_t[fMaxProtons+1];
  
  for (Int_t k = 0; k<(fMaxProtons+1); ++k)
  {
      // Define the filenames:
      kstr = st.Itoa(k,10);
      Name = "../../UI/Calibration/Energy_Cuts_" + kstr + "p.txt";
  
      // Open the files:
      ThisFile.open(Name);
      // Check if opened succesfully
      if(! ThisFile.is_open()) 
      {
        Fatal("ReadCalibrFile", "Calibration file does not exist! Aborting...");
      }

      // Now we must read the calibration file. First determine its size and process:
      NLines = 0;
      std::string ThisLine;
      while (std::getline(ThisFile,ThisLine)) {NLines = NLines + 1;}
      fNumberOfCuts[k] = NLines - 1;
      fCuts[k] = new Double_t[fNumberOfCuts[k]];
      ThisFile.close();
      
      // Then re-open and read:
      ThisFile.open(Name);
      ThisFile >> fKappa[k];
      for (Int_t kE = 0; kE<fNumberOfCuts[k]; ++kE) {ThisFile >> fCuts[k][kE];}
      ThisFile.close();  // some stupid text...
      // some stupid text...
      // Print output:
      cout << "-I- R3BNeutronTracker2D : calibration for " << k << " protons: " << endl;
      cout << "                                kappa : " << fKappa[k] << endl;
      for(Int_t i = 0; i < fNumberOfCuts[k]; i++) 
      {  // some stupid text...
         cout << "                                      : " << fCuts[k][i] << endl;
      }  // some stupid text...
  } // some stupid text...
} // some stupid text...


// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::Exec(Option_t* opt)
{
    ++fCounter;
  if(0 == (eventNo%100)) {
    cout << "-I- R3BNeutronTracker2D : Event #: " << eventNo << "." << endl;
  }
  eventNo += 1;


  Reset();


  R3BMCTrack *aTrack1;
  Int_t prim;
  Int_t particleID;
  for(Int_t it = 0; it < fLandMCTrack->GetEntriesFast(); it++) {
    aTrack1 = (R3BMCTrack*) fLandMCTrack->At(it);

    prim = aTrack1->GetMotherId();
    if(-1 != prim) {
      break;
    }
    
    particleID = aTrack1->GetPdgCode();
    
    // NOTE: The following tracks do not need to be coordinate-transformed!
    // The are in global coordinates and the only thing the coordinates do is to
    // compute the length, wihich is invariant under global-to-local.
    if(2112 == particleID) {
      //neutron
      PRIM_part[fNPrimNeut] = new R3BPrimPart(aTrack1->GetPdgCode(),
					      aTrack1->GetPx()*1000.,
					      aTrack1->GetPy()*1000.,
					      aTrack1->GetPz()*1000.,
					      aTrack1->GetStartX(),
					      aTrack1->GetStartY(),
					      aTrack1->GetStartZ(),
					      aTrack1->GetStartT(),
					      1.0086649, mNeutron);
      fNPrimNeut += 1;
    } else if(2212 == particleID) {
      //proton
      PRIM_prot[fNPrimProt] = new R3BPrimPart(aTrack1->GetPdgCode(),
					      aTrack1->GetPx()*1000.,
					      aTrack1->GetPy()*1000.,
					      aTrack1->GetPz()*1000.,
					      aTrack1->GetStartX(),
					      aTrack1->GetStartY(),
					      aTrack1->GetStartZ(),
					      aTrack1->GetStartT(),
					      1.00782503207, 1.00782503207*amu);
      fNPrimProt += 1;
    } else if(22 == particleID) {
      //gamma
      PRIM_gamma[fNPrimGamma] = new R3BPrimPart(aTrack1->GetPdgCode(),
						aTrack1->GetPx()*1000.,
						aTrack1->GetPy()*1000.,
						aTrack1->GetPz()*1000.,
						aTrack1->GetStartX(),
						aTrack1->GetStartY(),
						aTrack1->GetStartZ(),
						aTrack1->GetStartT(),
						0., 0.);
      fNPrimGamma += 1;
    } else{
      //fragment
      Double_t A = aTrack1->GetMass();
      PRIM_frag[fNofFrag] = new R3BPrimPart(aTrack1->GetPdgCode(),
					    aTrack1->GetPx()*1000.,
					    aTrack1->GetPy()*1000.,
					    aTrack1->GetPz()*1000.,
					    aTrack1->GetStartX(),
					    aTrack1->GetStartY(),
					    aTrack1->GetStartZ(),
					    aTrack1->GetStartT(),
					    A, A*amu);
      fNofFrag += 1;
    }
  }
      





  fNofClusters = fArrayCluster->GetEntries();

  // ------------------------------------------------------------------
  // Determine amount of neutrons
  // Loop over digis (also compute NeuLAND TOF):
  fNeuLAND_TOF = 1e99;
  sumTotalEnergy=0;   
  R3BLandDigi *digi;
  for(Int_t id = 0; id < fLandDigi->GetEntries(); id++) {
    digi = (R3BLandDigi*) fLandDigi->At(id);
    fEtot += digi->GetQdc();
    sumTotalEnergy += digi->GetQdc();
    fTOF = digi->GetTdc();
    if (fTOF<fNeuLAND_TOF) {fNeuLAND_TOF = fTOF;}
  }
  
  // Determine amount of protons:
  // Loop over veto digis:
  fvetoTOFhits = 0;
  R3Bveto_segDigi *digi_veto;
  fvetoHits = fveto_segDigi->GetEntries();
  for(Int_t id = 0; id < fvetoHits; id++) {
    digi_veto = (R3Bveto_segDigi*) fveto_segDigi->At(id);
    fTOF = digi_veto->GetTdc();
    if (fTOF<fNeuLAND_TOF) {fvetoTOFhits = fvetoTOFhits + 1;}    
  }
  fCurrentProtons = fvetoTOFhits; // FIXME: we treat every hit as a proton!
  
  // NOTE: we cannot access moe protons then e have calibration files for:
  if (fCurrentProtons>fMaxProtons) {fCurrentProtons = fMaxProtons;}
  
  // Double_t neutmult=0.0098379928527*fEtot/beamEnergy*600.;
  Double_t neutmult=0.00960093705145622332*fEtot/beamEnergy*600.;
  nNeut = (Int_t)(neutmult+0.5);
  hNeutmult->Fill(nNeut);

  h_nofclusters->Fill(fNofClusters);
  h_etot->Fill(fEtot);
  h_ncl_etot->Fill(fEtot, fNofClusters);
  h_ndigi_etot->Fill(fEtot, nentries);

  // Apply cuts:
  if (f2DCutEnabled)
  {
    // We have fNofClusters and fEtot which hold the point in the histogram. now we need to decide
    // between which cuts we are. The cuts go from 0 MeV (first) to 1e6 MeV (last). Then the lines
    // are defined as y = fKappa*(fCuts[k] - fEtot) So compute Y for 2 successive cuts and see
    // if fNofClusters is between them. 
    
    fdClusters = (Int_t) fNofClusters;

    for (Int_t k = 0; k<(fNumberOfCuts[fCurrentProtons]-1); ++k)
    {
      Ylow = fKappa[fCurrentProtons]*(fCuts[fCurrentProtons][k] - fEtot);
      Yhigh = fKappa[fCurrentProtons]*(fCuts[fCurrentProtons][k+1] - fEtot);
      if ((Ylow<fdClusters)&&(Yhigh>fdClusters)) {nNeut = k+1;}
    }
    if (fEtot<1.0) {nNeut = 0;}
  }
  else 
  {
      nNeut = fNNeutrons;
  }
  
  h_ntracks->Fill(nNeut);
  
     

  if(nNeut > 0) {
    // Find neutrons
    nAboveThresh=0;
    AdvancedMethod();
    hClusters1->Fill(fNofClusters1);
    h_ncl1_etot->Fill(fEtot, fNofClusters1);
  }
  
  // Calculate invariant mass of neutrons
  CalculateMassInv();
//   CalculateExce();
  
  hMinv0->Fill(fMinvTrue);
  hExce0->Fill(fExceTrue);
  if(nNeut != 0 && fNofTracks == fNPrimNeut) {
    hMinv->Fill(fMinv);
    hExce->Fill(fExce);
  }

}
// ----------------------------------------------------------------------------- 


    
// -----------------------------------------------------------------------------
Int_t R3BNeutronTracker2D::AdvancedMethod()
{
  // Loop over clusters
  R3BNeuLandCluster *cluster;
  for(Int_t ic = 0; ic < fNofClusters; ic++) {
    // Get pointer to the cluster
    cluster = (R3BNeuLandCluster*) fArrayCluster->At(ic);
    // Process starting from this cluster
    NextIteration(ic, cluster);
  }// clusters

  // Sort clusters
  SortClustersBeta();

  Int_t nOutput = 0;
  Bool_t mapUsed[3000];
  for(Int_t i = 0; i < 3000; i++) {
    mapUsed[i] = kFALSE;
  }
  for(Int_t i = 0; i < nNeut; i++) {
    for(Int_t ic = 0; ic < fNofClusters1; ic++) {
      cluster = fVectorCluster.at(ic);

      // Skip found clusters
      if(mapUsed[ic]) {
  	continue;
      }

      // Check energy of the cluster
      if(cluster->GetE() < 2.5 && ic > 0) {
  	continue;
      }

      // Check if in the beta-window
      TVector3 pos;
      cluster->StartPosition(pos);
      
      // NOTE: This start position is in global coordinates:
      Double_t Travel_Dist = TMath::Sqrt((pos.X() - fTarget_Xpos)*(pos.X() - fTarget_Xpos) + (pos.Y() - fTarget_Ypos)*(pos.Y() - fTarget_Ypos) + (pos.Z() - fTarget_Zpos)*(pos.Z() - fTarget_Zpos));
      Double_t t0 = TMath::Sqrt((fTarget_Xpos - fBeam_Xpos)*(fTarget_Xpos - fBeam_Xpos) + (fTarget_Ypos - fBeam_Ypos)*(fTarget_Ypos - fBeam_Ypos) + (fTarget_Zpos - fBeam_Zpos)*(fTarget_Zpos - fBeam_Zpos));
      t0 = t0/(gBeamBeta*c);
      Double_t time = cluster->GetStartT() - t0;
      Double_t beta = Travel_Dist/(time*c);
      
      // NOTE: Now we have overall beta (coordinate-invariant) and we can test it:
      if(TMath::Abs(beta-beamBeta) > (0.05*600./beamEnergy) && ic > 0) {
  	continue;
      }
     
      // Create neutron track
      fTracks[fNofTracks][0] = ic;
      mapUsed[ic] = kTRUE;
      fNofHits[fNofTracks] += 1;
      fMapTracks[fNofTracks] = kTRUE;
      new ((*fNeutHits)[fNofTracks]) R3BNeutHit(pos.X(),
						pos.Y(),
						pos.Z(),
						cluster->GetStartT());
      nOutput += 1;
      fNofTracks += 1;
      break;
    }// clusters
  }// nNeut

  return nOutput;
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::SortClustersBeta()
{
  if(0 == fNofClusters) {
    return;
  }

  R3BNeuLandCluster *cluster;
  for(Int_t ic = 1; ic < fNofClusters; ic++) {
    cluster = (R3BNeuLandCluster*) fArrayCluster->At(ic);
    if(cluster->IsMarked()) {
      continue;
    }
    fVectorCluster.push_back(cluster);
  }// Clusters

  fNofClusters1 = fVectorCluster.size();
  if(fNofClusters1 > 1) {
    std::sort(fVectorCluster.begin(), fVectorCluster.end(), AuxSortClustersBeta);
  }

  cluster = (R3BNeuLandCluster*) fArrayCluster->At(0);
  fVectorCluster.insert(fVectorCluster.begin(), cluster);
  fNofClusters1 = fVectorCluster.size();
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::NextIteration(Int_t curIndex, R3BNeuLandCluster *curClus)
{
  if(curClus->GetSize() < 2) {
    return;
  }

  // Loop over clusters
  R3BNeuLandCluster *cluster;
  for(Int_t ic = (curIndex+1); ic < fNofClusters; ic++) {
    cluster = (R3BNeuLandCluster*) fArrayCluster->At(ic);
    // Check if elastic scattering
    if(IsElastic(curClus, cluster)) {
      // Elastic, so mark the cluster
      cluster->Mark();
      // Stop the loop
      break;
    }
  }// clusters
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
Bool_t R3BNeutronTracker2D::IsElastic(R3BNeuLandCluster *c1, R3BNeuLandCluster *c2)
{
    // Set global target coordinates:
  Double_t x0=fTarget_Xpos;
  Double_t y0=fTarget_Ypos;
  Double_t z0=fTarget_Zpos;
  
  // NOTE: t must be interaction: The beta of the beam is stored in gBeamBeta.
  // vacuum speed of light is stored in c.
  Double_t t0 = TMath::Sqrt((fTarget_Xpos - fBeam_Xpos)*(fTarget_Xpos - fBeam_Xpos) + (fTarget_Ypos - fBeam_Ypos)*(fTarget_Ypos - fBeam_Ypos) + (fTarget_Zpos - fBeam_Zpos)*(fTarget_Zpos - fBeam_Zpos));
  t0 = t0/(gBeamBeta*c);

//  Double_t v1xmin,v1ymin,v1zmin,v1xmax,v1ymax,v1zmax;
  Double_t v4xmin,v4ymin,v4zmin;//,v4xmax,v4ymax,v4zmax;
  Double_t v3xmin,v3ymin,v3zmin;//,v3xmax,v3ymax,v3zmax;
//  Double_t v6xmin,v6ymin,v6zmin,v6xmax,v6ymax,v6zmax;

  //    search for scattering
  //    elastic scattering of particle 1 (neutron) on a particle 2 (proton) at rest.
  //    Outgoing particles are 3 (scattered neutron) and 4 (recoil proton).
  //    angle of particles after scattering is theta3 and theta4
  //    K is kinetic energy, p is momentum

  TVector3 pos1;
  TVector3 pos2;

  c1->StartPosition(pos1);

  //    incoming particle
  // vector from previous interaction to present interaction      
  // NOTE: We compute differences between global coordinates ==> Correct!
  Double_t v1x=pos1.X()-x0;
  Double_t v1y=pos1.Y()-y0;
  Double_t v1z=pos1.Z()-z0;
  // time difference and distance between previous interaction to present interaction
  Double_t dt=c1->GetStartT()-t0;	 
  Double_t dr=sqrt(v1x*v1x + v1y*v1y + v1z*v1z);

  Double_t beta1=dr/(dt*c);
  Double_t beta1max=(dr+dio)/dt/c;
  Double_t beta1min=(dr-dio)/dt/c;
  if (beta1max>0.99) beta1max=0.99;
  if (beta1min<0.) beta1min=0.;
  
  Double_t gamma1=1./sqrt(1.-beta1*beta1);  
  Double_t gamma1min=1./sqrt(1.-beta1min*beta1min);  
  Double_t gamma1max=1./sqrt(1.-beta1max*beta1max);  
  Double_t p1=beta1*gamma1*1.*amu;
  Double_t p1min=beta1min*gamma1min*1.*amu;
  Double_t p1max=beta1max*gamma1max*1.*amu;
  Double_t En1=sqrt(p1*p1+amu*amu);
  Double_t En1min=sqrt(p1min*p1min+amu*amu);
  Double_t En1max=sqrt(p1max*p1max+amu*amu);
  Double_t K1=En1-amu;	    
  Double_t K1min=En1min-amu;	    
  Double_t K1max=En1max-amu;	    
  
  // particle 4 is proton
  c1->StopPosition(pos2);
  // NOTE: We compute differences between global coordinates ==> Correct!
  Double_t v4x=(pos2-pos1).X();
  Double_t v4y=(pos2-pos1).Y();
  Double_t v4z=(pos2-pos1).Z();
  //vector perpendicular to scattering plane   
  Double_t v5x=v1y*v4z-v1z*v4y;
  Double_t v5y=v1z*v4x-v1x*v4z;
  Double_t v5z=v1x*v4y-v1y*v4x;
  
  Double_t tempAngle;
  dt=c1->GetStopT()-c1->GetStartT();
  dr=sqrt(v4x*v4x + v4y*v4y + v4z*v4z);
      	    
  Double_t theta4Measured=acos((v1x*v4x + v1y*v4y + v1z*v4z )
			      /sqrt(v1x*v1x + v1y*v1y + v1z*v1z)/sqrt(v4x*v4x + v4y*v4y + v4z*v4z));		
  Double_t theta4Measuredmin=theta4Measured;
  Double_t theta4Measuredmax=theta4Measured;
  if((v4x*v4x + v4y*v4y + v4z*v4z)==0.){
    theta4Measured=1.55;
    theta4Measuredmin=0.;
    theta4Measuredmax=1.55;
  }
  for(Int_t k=1;k<3;k++){
    if (k==1){
      v4xmin=v4x-dio;
    } 
    else{
      v4xmin=v4x+dio;
    }
    for(Int_t l=1;l<3;l++){
      if(l==1){
	v4ymin=v4y-dio;
      }
      else{
	v4ymin=v4y+dio;
      }
      for(Int_t m=1;m<3;m++){
	if(m==1){
	  v4zmin=v4z-dio;
	}
	else{
	  v4zmin=v4z+dio;
	}
	tempAngle=acos((v1x*v4xmin + v1y*v4ymin + v1z*v4zmin) 
		       /sqrt(v1x*v1x + v1y*v1y + v1z*v1z)
		       /sqrt(v4xmin*v4xmin + v4ymin*v4ymin + v4zmin*v4zmin));
	if (tempAngle>theta4Measuredmax) theta4Measuredmax=tempAngle;		     
	if (tempAngle<theta4Measuredmin) theta4Measuredmin=tempAngle;		     
      }
    }
  }
  if (theta4Measuredmax>1.55) theta4Measuredmax=1.55;
  
  // calculate velocity of neutron after scattering
  TVector3 pos3;
	 
  // Vector of scattered neutron    
  // NOTE: We compute differences between global coordinates ==> Correct!
  c2->StartPosition(pos3);
  Double_t v3x=(pos3-pos1).X();
  Double_t v3y=(pos3-pos1).Y();
  Double_t v3z=(pos3-pos1).Z();
  
  Double_t v6x=v1y*v3z-v1z*v3y;
  Double_t v6y=v1z*v3x-v1x*v3z;
  Double_t v6z=v1x*v3y-v1y*v3x;
  
  Double_t theta56=acos((v5x*v6x + v5y*v6y + v5z*v6z )
		       /sqrt(v5x*v5x + v5y*v5y + v5z*v5z)/sqrt(v6x*v6x + v6y*v6y + v6z*v6z));		     
  if((v5x*v5x + v5y*v5y + v5z*v5z)==0. || (v6x*v6x + v6y*v6y + v6z*v6z)==0.){
    theta56=3.14;
  }
  
  Double_t theta3=acos((v1x*v3x + v1y*v3y + v1z*v3z )
		      /sqrt(v1x*v1x + v1y*v1y + v1z*v1z)/sqrt(v3x*v3x + v3y*v3y + v3z*v3z));		     
  Double_t theta3min=theta3;
  Double_t theta3max=theta3;
  
  for(Int_t k=1;k<3;k++){
    if (k==1){
      v3xmin=v3x-dio;
    } 
    else{
      v3xmin=v3x+dio;
    }
    for(Int_t l=1;l<3;l++){
      if(l==1){
	v3ymin=v3y-dio;
      }
      else{
	v3ymin=v3y+dio;
      }
      for(Int_t m=1;m<3;m++){
	if(m==1){
	  v3zmin=v3z-dio;
	}
	else{
	  v3zmin=v3z+dio;
	}
	tempAngle=acos((v1x*v3xmin + v1y*v3ymin + v1z*v3zmin) 
		       /sqrt(v1x*v1x + v1y*v1y + v1z*v1z)
		       /sqrt(v3xmin*v3xmin + v3ymin*v3ymin + v3zmin*v3zmin));
	if (tempAngle>theta3max) theta3max=tempAngle;		     
	if (tempAngle<theta3min) theta3min=tempAngle;		     
      }
      
    }
  }
  
  if (theta3max>1.55) theta3max=1.55;
  
  dt=c2->GetStartT()-c1->GetStartT();
  dr=sqrt(v3x*v3x + v3y*v3y + v3z*v3z);

  Double_t beta3=dr/dt/c;
  Double_t beta3max=(dr+dio)/dt/c;
  Double_t beta3min=(dr-dio)/dt/c;
  
  if (beta3max>0.99) beta3max=0.99;
  if (beta3min<0.) beta3min=0.;
  
  Double_t gamma3=1./sqrt(1.-beta3*beta3);  
  Double_t gamma3min=1./sqrt(1.-beta3min*beta3min);  
  Double_t gamma3max=1./sqrt(1.-beta3max*beta3max);  
  Double_t p3=beta3*gamma3*1.*amu;
  Double_t p3min=beta3min*gamma3min*1.*amu;
  Double_t p3max=beta3max*gamma3max*1.*amu;
  if (p3>p1) p3=p1;
//  Double_t En3=sqrt(p3*p3+amu*amu);
//  Double_t En3min=sqrt(p3min*p3min+amu*amu);
//  Double_t En3max=sqrt(p3max*p3max+amu*amu);
//  Double_t K3=En3-amu;
//  Double_t K3min=En3min-amu;	    
//  Double_t K3max=En3max-amu;	    
  
  Double_t Ma,Mb,Mc,Md,Ka,Thc,Ei,Pa,AA,BB,a,b,cc,Pc1,Pc2;
  Double_t Pd1,Thd1,Ec1,Ed1,Kc1,Kd1,Qsqr1;
//  Double_t Pd2,Thd2,Ec2,Ed2,Kc2,Kd2,Qsqr2;
  Ma = 1.0087*amu;
  Mb = 1.0073*amu;
  Mc = Ma;
  Md = Mb;
  // calculate momentum of scattered proton and neutron
  Ka = K1;
  Thc = theta3;
  Ei=Ma +Mb +Ka;
  Pa = sqrt( Ka * Ka +2. * Ka * Ma );
  AA = Ei * Ei - Md * Md +Mc * Mc - Pa * Pa;
  BB = AA * AA - 4. * Ei * Ei * Mc * Mc;	
  a = 4. * Pa * Pa * cos(Thc) * cos(Thc) - 4. * Ei * Ei;
  b = 4. * AA * Pa * cos(Thc);
  cc = BB;	
  Pc1 = -b/(2. * a)+sqrt( (b * b)/( 4. * a * a) - ( cc / a ) );
  if(Pc1<0.) Pc1=0.;
  Pc2 = -b/(2. * a)-sqrt( (b * b)/( 4. * a * a) - ( cc / a ) );	
  Pd1 = sqrt( Pc1 * Pc1 +Pa * Pa - 2. * Pc1 * Pa * cos(Thc) );
  Thd1 = acos( ( Pc1 * Pc1 - Pd1 * Pd1 - Pa * Pa) / ( -2. * Pd1 * Pa ) ) ;
  Ec1 = sqrt( Pc1 * Pc1 + Ma * Ma );
  Ed1 = sqrt( Pd1 * Pd1 + Mb * Mb );
  Kc1 = Ec1 - Ma;
  Kd1 = Ed1 - Mb;
  Qsqr1 = (- ( Ka - Kc1 ) * (Ka - Kc1 ) +  ( Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc) ))/197.327/197.327;
  Double_t p3b=Pc1;
//  Double_t K3b=Kc1;
//  Double_t E3b=Ec1;	    
  Double_t theta4=Thd1;
  Double_t p4b=Pd1;
//  Double_t K4b=Kd1;
//  Double_t E4b=Ed1;
  
  Double_t p3bmin=p3b;
  Double_t p3bmax=p3b;
  Double_t p4bmin=p4b;
  Double_t p4bmax=p4b;
  Double_t theta4min=theta4;
  Double_t theta4max=theta4;
  
  // calculate minimum and maximum(within errors) momentum of scattered proton and neutron
  for (Int_t m=1;m<5;m++){
    if(m==1){
      Ka = K1min;
      Thc = theta3min;
    }
    if(m==2){
      Ka = K1min;
      Thc = theta3max;
    }
    if(m==3){
      Ka = K1max;
      Thc = theta3min;
    }
    if(m==4){
      Ka = K1max;
      Thc = theta3max;
    }
    
    Ei=Ma +Mb +Ka;
    Pa = sqrt( Ka * Ka +2. * Ka * Ma );
    AA = Ei * Ei - Md * Md +Mc * Mc - Pa * Pa;
    BB = AA * AA - 4. * Ei * Ei * Mc * Mc;	
    a = 4. * Pa * Pa * cos(Thc) * cos(Thc) - 4. * Ei * Ei;
    b = 4. * AA * Pa * cos(Thc);
    cc = BB;	
    Pc1 = -b/(2. * a)+sqrt( (b * b)/( 4. * a * a) - ( cc / a ) );
    if (((b * b)/( 4. * a * a) - ( cc / a ))<0.) Pc1=0.;
    if(Pc1<0.) Pc1=0.;
    Pc2 = -b/(2. * a)-sqrt( (b * b)/( 4. * a * a) - ( cc / a ) );	
    Pd1 = sqrt( Pc1 * Pc1 +Pa * Pa - 2. * Pc1 * Pa * cos(Thc) );
    Thd1 = acos( ( Pc1 * Pc1 - Pd1 * Pd1 - Pa * Pa) / ( -2. * Pd1 * Pa ) ) ;
    
    if(Pc1<p3bmin) p3bmin=Pc1;
    if(Pc1>p3bmax) p3bmax=Pc1;
    if(Pd1<p4bmin) p4bmin=Pd1;
    if(Pd1>p4bmax) p4bmax=Pd1;
    if(Thd1<theta4min) theta4min=Thd1;
    if(Thd1>theta4max) theta4max=Thd1;
  }
  
  Ec1 = sqrt( p3bmin * p3bmin + Ma * Ma );
  Kc1 = Ec1 - Ma;
  Ed1 = sqrt( p4bmin * p4bmin + Mb * Mb );
  Kd1 = Ed1 - Mb;
  Qsqr1 = (- ( Ka - Kc1 ) * (Ka - Kc1 ) +  ( Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc) ))/197.327/197.327;
//  Double_t K3bmin=Kc1;
//  Double_t E3bmin=Ec1;	    
  Double_t K4bmin=Kd1;
//  Double_t E4bmin=Ed1;
  Ec1 = sqrt( p3bmax * p3bmax + Ma * Ma );
  Kc1 = Ec1 - Ma;
  Ed1 = sqrt( p4bmax * p4bmax + Mb * Mb );
  Kd1 = Ed1 - Mb;
  Qsqr1 = (- ( Ka - Kc1 ) * (Ka - Kc1 ) +  ( Pa * Pa + Pc1 * Pc1 - 2. * Pa * Pc1 * cos(Thc) ))/197.327/197.327;
  
//  Double_t K3bmax=Kc1;
//  Double_t E3bmax=Ec1;	    
  Double_t K4bmax=Kd1;
//  Double_t E4bmax=Ed1;
  
  
  // decide if Cluster comes from scattered neutron or another neutron is needed!
  Double_t protonEnergy=8.76839+4.13858*c1->GetE()-0.00337368*c1->GetE()*c1->GetE();
  if(p3bmax>p3min && p3bmin<p3max && theta4max>theta4Measuredmin && theta4min<theta4Measuredmax &&
     K4bmax>protonEnergy && K4bmin<protonEnergy && theta56*180./3.14>120.) {
    return kTRUE;
  } else {
    return kFALSE;
  }
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::CalculateMassInv()
{
  R3BLandFirstHits *fh = (R3BLandFirstHits*) fLandFirstHits->At(0);
  TVector3 fhv[6];
  fhv[0].SetXYZ(fh->GetX0(), fh->GetY0(), fh->GetZ0());
  fhv[1].SetXYZ(fh->GetX1(), fh->GetY1(), fh->GetZ1());
  fhv[2].SetXYZ(fh->GetX2(), fh->GetY2(), fh->GetZ2());
  fhv[3].SetXYZ(fh->GetX3(), fh->GetY3(), fh->GetZ3());
  fhv[4].SetXYZ(fh->GetX4(), fh->GetY4(), fh->GetZ4());
  fhv[5].SetXYZ(fh->GetX5(), fh->GetY5(), fh->GetZ5());


  Double_t sum_energy_tr = 0.;
  Double_t sum_momentumX_tr = 0.;
  Double_t sum_momentumY_tr = 0.;
  Double_t sum_momentumZ_tr = 0.;
  Double_t sum_masses_tr = 0.;
  Double_t sum_energy = 0.;
  Double_t sum_momentumX = 0.;
  Double_t sum_momentumY = 0.;
  Double_t sum_momentumZ = 0.;
  Double_t sum_masses = 0.;


  for (Int_t i = 0; i < fNPrimNeut; i++) {
    TVector3 mom;
    PRIM_part[i]->Momentum(mom);
    // NOTE: We compute differences between global coordinates ==> Correct!
    Double_t momentumZ = mom.Z();
    Double_t momentumX = mom.X();
    Double_t momentumY = mom.Y();
    Double_t energy = sqrt(mom.Mag2()+mNeutron*mNeutron);

    sum_energy_tr += energy;
    sum_momentumX_tr += momentumX;
    sum_momentumY_tr += momentumY;
    sum_momentumZ_tr += momentumZ;
    sum_masses_tr += mNeutron;
  }

  R3BNeutHit *hit;
  for (Int_t i = 0; i < fNeutHits->GetEntriesFast(); i++) {
    // add up momentum of neutrons
    // neutron: beta is calculated with position and time of first hit

    hit = (R3BNeutHit*) fNeutHits->At(i);
    
    // Now transform the coordinates of the hit. Original code:
    Double_t momentumT = hit->GetP(); // total spatial momentum: frame-invariant ==> OK!
    Double_t momentumZ = momentumT*
      cos(TMath::ATan2(TMath::Sqrt(hit->GetX()*hit->GetX() +
				   hit->GetY()*hit->GetY()),
		       hit->GetZ()));
    Double_t momentumX = momentumZ*hit->GetX()/hit->GetZ();
    Double_t momentumY = momentumZ*hit->GetY()/hit->GetZ();
    // NOTE that momentum X,Y,Z are all in global coordinates ==> code still OK so far.
    
    Double_t energy = TMath::Sqrt(momentumT*momentumT + mNeutron*mNeutron);

    sum_energy += energy;
    sum_momentumX += momentumX;
    sum_momentumY += momentumY;
    sum_momentumZ += momentumZ;
    sum_masses += mNeutron;
      
    R3BLandPoint *point;
    TVector3 posPoint;
    Double_t d;
    Int_t index;
    Double_t dmin = 1e6;
    for(Int_t ip = 0; ip < fLandPoints->GetEntriesFast(); ip++) {
      point = (R3BLandPoint*) fLandPoints->At(ip);
      point->PositionIn(posPoint);
      // NOTE: We compute differences between global coordinates ==> Correct!
      d = TMath::Sqrt(TMath::Power(hit->GetX() - posPoint.X(), 2) +
		      TMath::Power(hit->GetY() - posPoint.Y(), 2) +
		      TMath::Power(hit->GetZ() - posPoint.Z(), 2));
      if(d < dmin) {
    	dmin = d;
    	index = ip;
      }
    }
    
    if (index<0) {index = 0;}
    if (index>=(fLandPoints->GetEntries())) {index = fLandPoints->GetEntries()-1;}
    
    point = (R3BLandPoint*) fLandPoints->At(index);
    Int_t trackID = point->GetTrackID();
    
    if (trackID<0) {trackID = 0;}
    if (trackID>=fLandMCTrack->GetEntries()) {trackID = fLandMCTrack->GetEntries()-1;}
    
    R3BMCTrack *track = (R3BMCTrack*) fLandMCTrack->At(trackID);
    Int_t mothID = track->GetMotherId();
    Int_t ind;
    while(-1 != mothID) {
        
      if (mothID<0) {mothID=0;}
      if (mothID>=fLandMCTrack->GetEntries()) {mothID = fLandMCTrack->GetEntries()-1;}
        
      track = (R3BMCTrack*) fLandMCTrack->At(mothID);
      ind = mothID;
      mothID = track->GetMotherId();
    }
    // NOTE: We compute differences between global coordinates ==> Correct!
    hDeltaPx2->Fill(momentumX-track->GetPx()*1000.);
    hDeltaPy2->Fill(momentumY-track->GetPy()*1000.);
    hDeltaPz2->Fill(momentumZ-track->GetPz()*1000.);

    Double_t deltaP = momentumT - track->GetP()*1000.;
    if(0 == i) {
      hDeltaP1->Fill(deltaP);
    } else if(1 == i) {
      hDeltaP2->Fill(deltaP);
    } else if(2 == i) {
      hDeltaP3->Fill(deltaP);
    } else if(3 == i) {
      hDeltaP4->Fill(deltaP);
    }
      
    dmin = 1e6;
    for(Int_t ifh = 0; ifh < 4; ifh++) {
        // NOTE: We compute differences between global coordinates ==> Correct!
      d = TMath::Sqrt(TMath::Power(hit->GetX() - fhv[ifh].X(), 2) +
		      TMath::Power(hit->GetY() - fhv[ifh].Y(), 2) +
		      TMath::Power(hit->GetZ() - fhv[ifh].Z(), 2));
      if(d < dmin) {
    	dmin = d;
    	index = ifh;
      }
    }// first hits
// NOTE: We compute differences between global coordinates ==> Correct!
    hDeltaX->Fill(hit->GetX()-fhv[index].X());
    hDeltaY->Fill(hit->GetY()-fhv[index].Y());
    hDeltaZ->Fill(hit->GetZ()-fhv[index].Z());
  }// tracks

      	
  // fragment: This information have to come later from the Tracker
  if(fNofFrag > 0) {
    TVector3 momF;
    PRIM_frag[0]->Momentum(momF);
    Double_t beta_frag=momF.Mag()/sqrt(momF.Mag2()+
				       PRIM_frag[0]->GetM()*PRIM_frag[0]->GetM());
    Double_t gamma_frag=1./sqrt(1.-beta_frag*beta_frag);
    Double_t energy_frag = gamma_frag*PRIM_frag[0]->GetM();
    sum_energy += energy_frag;
    sum_momentumX += momF.X();
    sum_momentumY += momF.Y();
    sum_momentumZ += momF.Z();
    sum_masses += PRIM_frag[0]->GetM();
    sum_energy_tr += energy_frag;
    sum_momentumX_tr += momF.X();
    sum_momentumY_tr += momF.Y();
    sum_momentumZ_tr += momF.Z();
    sum_masses_tr += PRIM_frag[0]->GetM();
  }


  // Gamma
  TVector3 momGamma;
  for(Int_t ig = 0; ig < fNPrimGamma; ig++) {
    PRIM_gamma[ig]->Momentum(momGamma);
    sum_energy += momGamma.Mag();      
    sum_momentumX += momGamma.X();
    sum_momentumY += momGamma.Y();
    sum_momentumZ += momGamma.Z();
    sum_energy_tr += momGamma.Mag();      
    sum_momentumX_tr += momGamma.X();
    sum_momentumY_tr += momGamma.Y();
    sum_momentumZ_tr += momGamma.Z();
  }

      
  // calculate invariant mass
  // calculate relative energy
  Double_t totMom2_tr = sum_momentumX_tr*sum_momentumX_tr + sum_momentumY_tr*sum_momentumY_tr +
    sum_momentumZ_tr*sum_momentumZ_tr;
  Double_t totMom2 = sum_momentumX*sum_momentumX + sum_momentumY*sum_momentumY +
    sum_momentumZ*sum_momentumZ;
  fMinvTrue = (sqrt(sum_energy_tr*sum_energy_tr-totMom2_tr) - sum_masses_tr);
  fMinv = (sqrt(sum_energy*sum_energy-totMom2) - sum_masses);
  // cout << fMinvTrue << "   " << fMinv << endl;
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::CalculateExce()
{
  TVector3 mom1;
  TVector3 mom2;
  Double_t en_g = 0;
  for(Int_t i = 0; i < fNPrimGamma; i++) {
    PRIM_gamma[i]->Momentum(mom1);
    en_g += mom1.Mag();
  }

  Int_t N = fNofFrag + fNPrimProt + fNPrimNeut - 1;   // !!!!!!!!!!!!! FIXME
  R3BPrimPart* p[30];

  Int_t ii = 0;
  for(Int_t i = 0; i < fNofFrag; i++) {
    p[ii] = PRIM_frag[i];
    ii += 1;
  }
  for(Int_t i = 0; i < 1/*fNPrimProt*/; i++) {    // !!!!!!!!!!!!! FIXME
    p[ii] = PRIM_prot[i];
    ii += 1;
  }
  for(Int_t i = 0; i < fNPrimNeut; i++) {
    p[ii] = PRIM_part[i];
    ii += 1;
  }

  // MC truth - ideal neutron reconstruction -------------------------
  Double_t sum_mc = 0.;
  for(Int_t i = 0; i < N; i++) {
    p[i]->Momentum(mom1);
    for(Int_t j = 0; j < N; j++) {
      p[j]->Momentum(mom2);
      sum_mc += p[i]->GetGamma()*p[j]->GetGamma()*
	p[i]->GetM()*p[j]->GetM()*
	(1. - p[i]->GetBeta()*p[j]->GetBeta()*cos(mom1.Angle(mom2)));
    }
  }
  // Double_t m_proj = 135.939300*amu;
  Double_t m_proj = 2.0141017778*amu;
  fExceTrue = sqrt(sum_mc) + en_g - m_proj;
  // -----------------------------------------------------------------


  // RECO spectrum ---------------------------------------------------
  if(! fNofTracks) {
    fExce = -1.;
    return;
  }

  N = fNofFrag + fNPrimProt + fNofTracks - 1;   // !!!!!!!!!!!!! FIXME
  ii = 0;
  for(Int_t i = 0; i < fNofFrag; i++) {
    p[ii] = PRIM_frag[i];
    ii += 1;
  }
  for(Int_t i = 0; i < 1/*fNPrimProt*/; i++) {   // !!!!!!!!!!!!! FIXME
    p[ii] = PRIM_prot[i];
    ii += 1;
  }
  Int_t k = ii;
  for(Int_t i = 0; i < fNofTracks; i++) {
    R3BNeuLandCluster *c1 = fVectorCluster.at(fTracks[i][0]);
    TVector3 pos;
    c1->StartPosition(pos);
    // NOTE: Shift pos over the target position to get the right numbers:
    pos.SetX(pos.X() - fTarget_Xpos);
    pos.SetY(pos.Y() - fTarget_Ypos);
    pos.SetZ(pos.Z() - fTarget_Zpos);
    
    Double_t s = pos.Mag();
    Double_t rr = pos.Perp();
    
    // To compute beta: shift time also to target time:
    Double_t t0 = TMath::Sqrt((fTarget_Xpos - fBeam_Xpos)*(fTarget_Xpos - fBeam_Xpos) + (fTarget_Ypos - fBeam_Ypos)*(fTarget_Ypos - fBeam_Ypos) + (fTarget_Zpos - fBeam_Zpos)*(fTarget_Zpos - fBeam_Zpos));
    t0 = t0/(gBeamBeta*c);
    Double_t time = c1->GetStartT() - t0;

    Double_t beta = s/(c*time);
    Double_t gamma = 1./sqrt(1.-beta*beta);
    Double_t momT = beta*gamma*mNeutron;
    Double_t momZ = momT*cos(TMath::ATan2(rr, pos.Z()));
    Double_t momX = momZ*pos.X()/pos.Z();
    Double_t momY = momZ*pos.Y()/pos.Z();
    p[k] = new R3BPrimPart(2112, momX, momY, momZ,
			   t0, fTarget_Xpos, fTarget_Ypos, fTarget_Zpos, // NOTE: This is where the particle comes from!
			   1.0086649, mNeutron);
    k += 1;
  }
  Double_t sum = 0.;
  for(Int_t i = 0; i < N; i++) {
    p[i]->Momentum(mom1);
    for(Int_t j = 0; j < N; j++) {
      p[j]->Momentum(mom2);
      sum += p[i]->GetGamma()*p[j]->GetGamma()*
	p[i]->GetM()*p[j]->GetM()*
	(1. - p[i]->GetBeta()*p[j]->GetBeta()*cos(mom1.Angle(mom2)));
    }
  }
  fExce = sqrt(sum) + en_g - m_proj;
  for(Int_t i = ii; i < N; i++) {
    delete p[i];
  }
  // -----------------------------------------------------------------
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::Reset()
{
  for(Int_t i = 0; i < fNPrimNeut; i++) {
    delete PRIM_part[i];
  }
  for(Int_t i = 0; i < fNPrimProt; i++) {
    delete PRIM_prot[i];
  }
  for(Int_t i = 0; i < fNPrimGamma; i++) {
    delete PRIM_gamma[i];
  }
  for(Int_t i = 0; i < fNofFrag; i++) {
    delete PRIM_frag[i];
  }

  fMinvTrue = 0.;
  fMinv = 0.;
  fExceTrue = 0.;
  fExce = 0.;
  fEtot = 0.;

  fNPrimNeut = 0;
  fNPrimProt = 0;
  fNPrimGamma = 0;
  fNofFrag = 0;
  fNofTracks = 0;
  fNofClusters = 0;
  fNofClusters1 = 0;

  if(fVectorCluster.size() > 0) {
    fVectorCluster.clear();
  }
  
  // Reset the neutron tracks
  for(Int_t i = 0; i < 1000; i++) {
    fMapTracks[i] = kTRUE;
    for(Int_t k = 0; k < 300; k++) {
      fTracks[i][k] = -1;
    }
    fNofHits[i] = 0;
  }

  fNeutHits->Clear();
}   
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::Finish()
{
  hClusters1->Write();

  hNeutmult->Write();

  hMinv->Write();
  hMinv0->Write();
  hMinv1->Write();
  hMinv2->Write();

  hExce->Write();
  hExce0->Write();

  hErel1->Write();
  hErel2->Write();
  hErel3->Write();
  hErel4->Write();

  hDeltaPx1->Write();
  hDeltaPy1->Write();
  hDeltaPz1->Write();
  hDeltaPx2->Write();
  hDeltaPy2->Write();
  hDeltaPz2->Write();

  hDeltaX->Write();
  hDeltaY->Write();
  hDeltaZ->Write();
  hDeltaT->Write();

  hDeltaP1->Write();
  hDeltaP2->Write();
  hDeltaP3->Write();
  hDeltaP4->Write();
  hDeltaP5->Write();
  hDeltaP6->Write();
   
  hDelta->Write();
   
  hFirstHitZ->Write();

  h_theta_cand->Write();
  h_beta_min->Write();
  h_beta_max->Write();
  h_ntracks->Write();
  h_nhits->Write();
  h_theta->Write();

  h_nofclusters->Write();
  h_etot->Write();
  h_ncl_etot->Write();
  h_ncl_etot_1->Write();
  h_ndigi_etot->Write();
  h_ncl1_etot->Write();
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::UseBeam(Double_t beam_energy, Double_t beam_beta)
{
  beamEnergy = beam_energy;
  beamBeta = beam_beta;
  gBeamBeta = beam_beta;
}
// -----------------------------------------------------------------------------




// -----------------------------------------------------------------------------
void R3BNeutronTracker2D::CreateHistograms()
{
  hNeutmult = new TH1F("Neutmult","Neutron multiplicity from energy considerations",10,-0.5,9.5);
  hNeutmult->GetXaxis()->SetTitle("Number of Neutrons");
  hNeutmult->GetYaxis()->SetTitle("Counts");

  hErel1 = new TH1F("Erel1","Erel for 1 neutron",400,0.,40.);
  hErel1->GetXaxis()->SetTitle("Erel (MeV)");
  hErel1->GetYaxis()->SetTitle("Counts");

  hErel2 = new TH1F("Erel2","Erel for 2 neutron",400,0.,40.);
  hErel2->GetXaxis()->SetTitle("Erel (MeV)");
  hErel2->GetYaxis()->SetTitle("Counts");

  hErel3 = new TH1F("Erel3","Erel for 3 neutron",400,0.,40.);
  hErel3->GetXaxis()->SetTitle("Erel (MeV)");
  hErel3->GetYaxis()->SetTitle("Counts");

  hErel4 = new TH1F("Erel4","Erel for 4 neutron",400,0.,40.);
  hErel4->GetXaxis()->SetTitle("Erel (MeV)");
  hErel4->GetYaxis()->SetTitle("Counts");

  hMinv = new TH1F("Minv","Minv for reconstructed hits",4000,-0.01,39.99);
  // hMinv = new TH1F("Minv","Minv for reconstructed hits",1000,0.,1000.);
  hMinv->GetXaxis()->SetTitle("Erel (MeV)");
  hMinv->GetYaxis()->SetTitle("Counts");

  hExce = new TH1F("Exce","Reconstructed Exce",100,0.,600.);

  hMinv1 = new TH1F("Minv1","Minv for first hits",4000,0.,40.);
  hMinv1->GetXaxis()->SetTitle("Erel (MeV)");
  hMinv1->GetYaxis()->SetTitle("Counts");

  hMinv2 = new TH1F("Minv2","Minv for selected hits",4000,0.,40.);
  hMinv2->GetXaxis()->SetTitle("Erel (MeV)");
  hMinv2->GetYaxis()->SetTitle("Counts");

  hMinv0 = new TH1F("Minv0","Minv for ideal hits",2000,-0.01,39.99);
  // hMinv0 = new TH1F("Minv0","Minv for ideal hits",1000,0.,1000.);
  hMinv0->GetXaxis()->SetTitle("Erel (MeV)");
  hMinv0->GetYaxis()->SetTitle("Counts");

  hExce0 = new TH1F("Exce0","Input Exce",100,0.,600.);

  hDeltaX = new TH1F("DeltaX","error in x determination",300,-150.,150.);
  hDeltaX->GetXaxis()->SetTitle("x position (cm)");
  hDeltaX->GetYaxis()->SetTitle("Counts");
     
  hDeltaY = new TH1F("DeltaY","error in y determination",300,-150.,150.);
  hDeltaY->GetXaxis()->SetTitle("y position (cm)");
  hDeltaY->GetYaxis()->SetTitle("Counts");
     
  hDeltaZ = new TH1F("DeltaZ","error in z determination",500,-250.,250.);
  hDeltaZ->GetXaxis()->SetTitle("z position (cm)");
  hDeltaZ->GetYaxis()->SetTitle("Counts");
     
  hDeltaT = new TH1F("DeltaT","error in time determination",2000,-10.,10.);
  hDeltaT->GetXaxis()->SetTitle("time (ns)");
  hDeltaT->GetYaxis()->SetTitle("Counts");

  hDeltaP1 = new TH1F("DeltaP1","difference in reconstucted momenta for 1st neutron)",1000,-150.,150.);
  hDeltaP1->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP1->GetYaxis()->SetTitle("Counts");
  
  hDeltaP2 = new TH1F("DeltaP2","difference in reconstucted momenta for 2nd neutron)",1000,-150.,150.);
  hDeltaP2->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP2->GetYaxis()->SetTitle("Counts");
  
  hDeltaP3 = new TH1F("DeltaP3","difference in reconstucted momenta for 3rd neutron)",1000,-150.,150.);
  hDeltaP3->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP3->GetYaxis()->SetTitle("Counts");
  
  hDeltaP4 = new TH1F("DeltaP4","difference in reconstucted momenta for 4th neutron)",1000,-150.,150.);
  hDeltaP4->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP4->GetYaxis()->SetTitle("Counts");
  
  hDeltaP5 = new TH1F("DeltaP5","difference in reconstucted momenta for 5th neutron)",1000,-150.,150.);
  hDeltaP5->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP5->GetYaxis()->SetTitle("Counts");
  
  hDeltaP6 = new TH1F("DeltaP6","difference in reconstucted momenta for 6th neutron)",1000,-150.,150.);
  hDeltaP6->GetXaxis()->SetTitle("Delta P (MeV/c)");
  hDeltaP6->GetYaxis()->SetTitle("Counts");  

  hDeltaPx1 = new TH1F("DeltaPx1","difference in reconstucted momenta px (ideal case)",1000,-50.,50.);
  hDeltaPx1->GetXaxis()->SetTitle("Delta Px (MeV/c)");
  hDeltaPx1->GetYaxis()->SetTitle("Counts");

  hDeltaPy1 = new TH1F("DeltaPy1","difference in reconstucted momenta py (ideal case)",1000,-50.,50.);
  hDeltaPy1->GetXaxis()->SetTitle("Delta Py (MeV/c)");
  hDeltaPy1->GetYaxis()->SetTitle("Counts");

  hDeltaPz1 = new TH1F("DeltaPz1","difference in reconstucted momenta pz (ideal case)",1000,-150.,150.);
  hDeltaPz1->GetXaxis()->SetTitle("Delta Pz (MeV/c)");
  hDeltaPz1->GetYaxis()->SetTitle("Counts");

  hDeltaPx2 = new TH1F("DeltaPx2","difference in reconstucted momenta px (exp case)",1000,-50.,50.);
  hDeltaPx2->GetXaxis()->SetTitle("Delta Px (MeV/c)");
  hDeltaPx2->GetYaxis()->SetTitle("Counts");

  hDeltaPy2 = new TH1F("DeltaPy2","difference in reconstucted momenta py (exp case)",1000,-50.,50.);
  hDeltaPy2->GetXaxis()->SetTitle("Delta Py (MeV/c)");
  hDeltaPy2->GetYaxis()->SetTitle("Counts");
 
  hDeltaPz2 = new TH1F("DeltaPz2","difference in reconstucted momenta pz (exp case)",1000,-150.,150.);
  hDeltaPz2->GetXaxis()->SetTitle("Delta Pz (MeV/c)");
  hDeltaPz2->GetYaxis()->SetTitle("Counts");

  hClusters1 = new TH1F("Cluster1","Number of clusters after eliminating elastic scattering",100,-0.5,99.5);
  hClusters1->GetXaxis()->SetTitle("number of clusters");
  hClusters1->GetYaxis()->SetTitle("Counts");

  hDelta = new TH1F("Delta","distance between two primary interactions",300,-150.,150.);
  hDelta->GetXaxis()->SetTitle("distance (cm)");
  hDelta->GetYaxis()->SetTitle("Counts");

  hFirstHitZ = new TH1F("FirstHitZ","z positions of first hits",200,1000.,2000.);
  hFirstHitZ->GetXaxis()->SetTitle("z position (cm)");
  hFirstHitZ->GetYaxis()->SetTitle("Counts");

  h_theta_cand = new TH1F("h_theta_cand", "Polar angle of a track candidate", 180, 0., 180.);
  h_theta_cand->GetXaxis()->SetTitle("#theta (deg)");
  h_theta_cand->GetYaxis()->SetTitle("counts");

  h_beta_min = new TH1F("h_beta_min", "velocity of the potential neutron", 400, -2., 2.);
  h_beta_min->GetXaxis()->SetTitle("#beta_{min}");
  h_beta_min->GetYaxis()->SetTitle("counts");

  h_beta_max = new TH1F("h_beta_max", "velocity of the potential neutron", 1200, -2., 10.);
  h_beta_max->GetXaxis()->SetTitle("#beta_{max}");
  h_beta_max->GetYaxis()->SetTitle("counts");

  h_ntracks = new TH1F("h_ntracks", "Number of reconstructed neutron tracks", 50, -0.5, 49.5);
  h_ntracks->GetXaxis()->SetTitle("N_{tracks}");
  h_ntracks->GetYaxis()->SetTitle("events");

  h_nhits = new TH1F("h_nhits", "Number of hits on a neutron track", 50, -0.5, 49.5);
  h_nhits->GetXaxis()->SetTitle("N_{hits} / track");
  h_nhits->GetYaxis()->SetTitle("all tracks");

  h_theta = new TH1F("h_theta", "Polar angle", 180, 0., 180.);
  h_theta->GetXaxis()->SetTitle("#theta (deg)");
  h_theta->GetYaxis()->SetTitle("counts");

  h_nofclusters = new TH1F("h_nofclusters", "Number of clusters", 100, -0.5, 99.5);

  h_etot = new TH1F("h_etot", "Total light", 100, 0., 1000.);

  if(beamEnergy < 300) {
    h_ncl_etot = new TH2F("h_ncl_etot", "Number of clusters vs. total light", 1500, 0., 1500., 150, -0.5, 149.5);
  } else {
    h_ncl_etot = new TH2F("h_ncl_etot", "Number of clusters vs. total light", 150, 0., 1500., 150, -0.5, 149.5);
  }

  h_ncl_etot_1 = new TH2F("h_ncl_etot_1", "Number of clusters vs. total light after cut", 100, 0., 1000., 100, -0.5, 99.5);

  h_ndigi_etot = new TH2F("h_ndigi_etot", "", 100, 0., 1000., 500, -0.5, 499.5);

  h_ncl1_etot = new TH2F("h_ncl1_etot", "", 100, 0., 1000., 100, -0.5, 99.5);
}
// -----------------------------------------------------------------------------




bool AuxSortClustersBeta(R3BNeuLandCluster *c1, R3BNeuLandCluster *c2)
{
  TVector3 pos1;
  TVector3 pos2;
  c1->StartPosition(pos1);
  c2->StartPosition(pos2);
  
  // NOTE: These start positions are in global coordinates.
  Double_t t0 = TMath::Sqrt((fTarget_Xpos_global - fBeam_Xpos_global)*(fTarget_Xpos_global - fBeam_Xpos_global) + (fTarget_Ypos_global - fBeam_Ypos_global)*(fTarget_Ypos_global - fBeam_Ypos_global) + (fTarget_Zpos_global - fBeam_Zpos_global)*(fTarget_Zpos_global - fBeam_Zpos_global));
  t0 = t0/(gBeamBeta*c);
  
  Double_t Travel_Dist1 = TMath::Sqrt((pos1.X() - fTarget_Xpos_global)*(pos1.X() - fTarget_Xpos_global) + (pos1.Y() - fTarget_Ypos_global)*(pos1.Y() - fTarget_Ypos_global) + (pos1.Z() - fTarget_Zpos_global)*(pos1.Z() - fTarget_Zpos_global));
  Double_t time1 = c1->GetStartT() - t0;
  Double_t beta1 = Travel_Dist1/(time1*c);
  
  Double_t Travel_Dist2 = TMath::Sqrt((pos2.X() - fTarget_Xpos_global)*(pos2.X() - fTarget_Xpos_global) + (pos2.Y() - fTarget_Ypos_global)*(pos2.Y() - fTarget_Ypos_global) + (pos2.Z() - fTarget_Zpos_global)*(pos2.Z() - fTarget_Zpos_global));
  Double_t time2 = c2->GetStartT() - t0;
  Double_t beta2 = Travel_Dist2/(time2*c);
  
  // Now continue with the Original code:
  
  Double_t val1 = TMath::Abs(beta1-gBeamBeta)/c1->GetE();
  Double_t val2 = TMath::Abs(beta2-gBeamBeta)/c2->GetE();
  if(val1 < val2) {
    return true;
  } else {
    return false;
  }
}