void KVParticle::SetFrame(const Char_t* frame, const TVector3& boost,
                          Bool_t beta)
{
   //Define a Lorentz-boosted frame in which to calculate the particle's momentum and energy.
   //
   //The new frame will have the name given in the string "frame", which can then be used to
   //access the kinematics of the particle in different frames using GetFrame().
   //
   //The boost velocity vector is that of the boosted frame with respect to the original frame of the particles in the event.
   //The velocity vector can be given either in cm/ns units (default) or in units of 'c' (beta=kTRUE).
   //
   //E.g. to define a frame moving at 0.1c in the +ve z-direction with respect to the original
   //event frame:
   //
   //      (...supposing a valid pointer KVParticle* my_part...)
   //      TVector3 vframe(0,0,0.1);
   //      my_part->SetFrame("my_frame", vframe, kTRUE);
   //
   //In order to access the kinematics of the particle in the new frame:
   //
   //      my_part->GetFrame("my_frame")->GetTransverseEnergy();// transverse energy in "my_frame"

   //set up TLorentzRotation corresponding to boosted frame
   TLorentzRotation tmp;
   if (beta) {
      tmp.Boost(boost);
   } else {
      tmp.Boost(boost.X() / kSpeedOfLight, boost.Y() / kSpeedOfLight,
                boost.Z() / kSpeedOfLight);
   }
   SetFrame(frame, tmp);
}

template <typename T> TVector4<T>::TVector4( const TVector3<T> & in, T w )
	: m_x(in.GetX())
	, m_y(in.GetY())
	, m_z(in.GetZ())
	, m_w(w)
{
}

TVector3 find_vertex(line_vec track_a, line_vec track_b) {
  
  TVector3 start_diff = track_a.start() - track_b.start();
  TVector3 unit_a = track_a.dir_unit();
  TVector3 unit_b = track_b.dir_unit();

  float ta = ( -(start_diff*unit_a) + (start_diff *  unit_b) * (unit_a * unit_b) )/ 
    ( 1.0 - ( (unit_a * unit_b) *  (unit_a * unit_b)) );
  float tb =   ( (start_diff*unit_b) - (start_diff *  unit_a) * (unit_a * unit_b) )/ 
    ( 1.0 - ( (unit_a * unit_b) *  (unit_a * unit_b)) );

  TVector3 close_a = track_a.start()  + (ta *  unit_a);
  TVector3 close_b = track_b.start()  + (tb * unit_b);
  
  // now what I really want to store as a vertex is the middle of the Vector going from close_a to close_b
  
  cout << "find_vertex::Closest point on Vector A = " << close_a.X() << ", " << close_a.Y() << " ," << close_a.Z() << endl;
  cout << "find_vertex::Closest point on Vector B = " << close_b.X() << ", " << close_b.Y() << " ," << close_b.Z() << endl;

  TVector3 from_a_to_b = close_b - close_a;
  // conceptually that's more of a point than a vector
  TVector3 vertex = close_a + 0.5 * from_a_to_b;
  cout << "find_vertex::The vertex is here:  " << vertex.X() << ", " <<  vertex.Y() << ", " <<  vertex.Z() << endl;

  return vertex;
}

vector<TParticle*> LMCphysGen::GenCerPhotons(TParticle *part, LMCstep *step, Int_t N) {

  const Double_t hbarc = 0.197326960277E-6; //GeV.nm 
  vector <TParticle*> v;
  //rotation matrix to transform vectors from geom -> part
  TRotation rm;
  TVector3 vt(part->Px(), part->Py(), part->Pz());
  rm.RotateX(vt.Theta()); //rotate on X
  rm.RotateZ(vt.Phi()); //rotate on Z
  TVector3 pdir = vt.Unit();
  // rotation matrix from part->geom
  TRotation im = rm.Inverse();
  //generate photons
  Double_t Emin = 2.*TMath::Pi()* hbarc / 200.; //GeV
  Double_t Emax = 2.*TMath::Pi()* hbarc / 700.; //GeV

  TVector3 InitPoint = step->GetInitPoint(); 
  Double_t stepLength = step->GetStepLength();
  for (int i=0; i< N; i++) {
    Double_t thetaCer = step->GetCerAngle();
    Double_t phi = (RandGen->Rndm())*2.*TMath::Pi();  
    TVector3 photonDir(sin(thetaCer)*cos(phi), sin(thetaCer)*sin(phi), cos(thetaCer));
    TVector3 photonDirDet = im*photonDir;                         //transform photon direction from particle to detector frame
    Double_t photonE      = Emin + (RandGen->Rndm())*(Emax-Emin); //cerenkov photon flat on energy (GeV)
    TVector3 photonMomDet = photonE*photonDirDet;                 //on detector frame
    TVector3 GenPoint     = InitPoint + stepLength*(RandGen->Rndm())*pdir;
    v.push_back(new TParticle(22,0,0,0,0,0,photonMomDet.X(),photonMomDet.Y(),photonMomDet.Z(),photonE,GenPoint.X(),GenPoint.Y(),GenPoint.Z(),0));
  }

  return v;

}

Double_t getShiftChi2(const Double_t* thetaPhi) {
   Double_t chi2(0), c(0);

   TVector3 norm; norm.SetMagThetaPhi(1.0, thetaPhi[0], thetaPhi[1]);
   for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
      for (Int_t xc=0; xc<ch; ++xc) {
         Double_t dtcor=0;
         for (Int_t i=(ch-xc); i>0; --i) {
            dtcor += dtCorrs[ch-i];
         }
         const TVector3& posCh = getStnPos(ch);
         const TVector3& posXc = getStnPos(xc);
         const Double_t  disCh = -(posCh.Dot(norm));
         const Double_t  disXc = -(posXc.Dot(norm));
         // !!! check sign of delta(distance) and dtcor!
         const Double_t  dt = kSmpRate * (             // from ns to samples
            ( (disCh-disXc) * kNgTopFern / kC_m_ns )   // from m to ns
            + dtcor );                                 // correct dt offset (ns)
         // FIXME: dt in samples, maxdt in ns
         if (TMath::Abs(dt)>maxdt) {
            // really don't like being out of bounds
            c = TMath::Exp(dt*dt);
            if (c>kReallyBig) {
               c = kReallyBig;
            }
         } else {
            // get the correlation coefficient for this dt (in num samples)
            c = gspc[ch][xc]->Eval(dt) - 1.0;
         }
         chi2 += c*c;
      }
   }
   return chi2;
}

double coeff_of_t(TVector3 &efield, TVector3 &vel, int dir)
{
  //returns A C(t), the coefficient of each waveguide mode in +x-direction
  // in units of C*Ohm/s or volts
  efield.SetZ( dir * efield.Z() );//only effect TM modes
  return Echarge * tl_data.Zw / 2 * ( efield * vel ) / US2S;
}

void UEAnalysisOnRootple::MPIAnalysisRECO(Float_t weight,string tkpt)
{
  vector<TVector3*> JetRECO;
  JetRECO.clear();
  
  for(int j=0;j<NumberTracksJet;j++){
    if(fabs(EtaCJ[j])<etaRegion){
      TVector3* jetvector = new TVector3;
      //jetvector->SetPtEtaPhi(CalibrationPt(TrasverseMomentumTJ[j],tkpt)*TrasverseMomentumTJ[j], EtaTJ[j], PhiTJ[j]);
      jetvector->SetPtEtaPhi(TrasverseMomentumTJ[j], EtaTJ[j], PhiTJ[j]);
      JetRECO.push_back(jetvector);
    }
  }
  
  vector<AssociatedObject> assoJetRECO;
  assoJetRECO.clear();

  while(JetRECO.size()>1){
    int oldSize = JetRECO.size();
    vector<TVector3*>::iterator itH = JetRECO.begin();
    if((*itH)->Pt()>=ptThreshold){
      for(vector<TVector3*>::iterator it=JetRECO.begin();it!=JetRECO.end();it++){
	float azimuthDistanceJet = fabs( (*itH)->Phi() - (*it)->Phi() );
	if((*it)->Pt()/(*itH)->Pt()>=0.3){
	  if( (piG - rangePhi) <  azimuthDistanceJet && azimuthDistanceJet < (piG + rangePhi)) {
	    AssociatedObject tmpPair((*itH),(*it));
	    assoJetRECO.push_back(tmpPair);
	    JetRECO.erase(it);
	    int newSize = oldSize -1;
	    oldSize = newSize;
	    JetRECO.resize(newSize);
	    break;
	  }
	}
      }
    }
    JetRECO.erase(itH);
    int newSize = oldSize -1;
    JetRECO.resize(newSize);
  }
  
  if(assoJetRECO.size()){
    fNumbMPIRECO->Fill(assoJetRECO.size());
    vector<AssociatedObject>::iterator at= assoJetRECO.begin();
    
    const TVector3* leadingJet((*at).first);
    const TVector3* secondJet((*at).second);

    pPtRatio_vs_PtJleadRECO->Fill(leadingJet->Pt(),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_EtaJleadRECO->Fill(fabs(leadingJet->Eta()),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_PhiJleadRECO->Fill(leadingJet->Phi(),(secondJet->Pt()/leadingJet->Pt()));
    
    fdEtaLeadingPairRECO->Fill(leadingJet->Eta()-secondJet->Eta());
    float dPhiJet = fabs(leadingJet->Phi()-secondJet->Phi());
    if(dPhiJet> piG) dPhiJet = 2*piG -dPhiJet;
    dPhiJet = (180*dPhiJet)/piG;
    fdPhiLeadingPairRECO->Fill(dPhiJet);
    fptRatioLeadingPairRECO->Fill(secondJet->Pt()/leadingJet->Pt());
  }
}

Path::Path(const TLorentzVector& p4, const TVector3& origin, double field)
    : m_unitDirection(p4.Vect().Unit()),
      m_speed(p4.Beta() * gconstc),
      m_origin(origin.X(), origin.Y(), origin.Z()),
      m_field(field) {
  m_points[papas::Position::kVertex] = m_origin;
}

TGraph GetCurve(int Points,const double & hi_ex_set)
{
  TGraph curve;
  if(!gPrimaryReaction.IsSet()){
    std::cout<<"Reaction Masses have not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }
  if(!gPrimaryReaction.BeamEnergy()){
    std::cout<<"Beam Energy has not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }

  sim::RN_SimEvent evt1(gPrimaryReaction.BeamEnergy(),gPrimaryReaction.M_Beam(),gPrimaryReaction.M_Target(),gPrimaryReaction.M_Recoil(),gPrimaryReaction.M_Fragment()); 
  // Fill the points of the kinematic curve
  int p=0;
  while(p<Points){
    double theta_deg = 180.0*p/Points;
    double phi=2.*M_PI*global::myRnd.Rndm();

    TVector3 nv; nv.SetMagThetaPhi(1.,theta_deg*M_PI/180.0,phi);
    if(!evt1.radiate_in_CM(nv,hi_ex_set))
      continue;
    else 
      curve.SetPoint(p, evt1.getLVrad().Theta()*180/3.14,(double)(evt1.getLVrad().E()-evt1.getLVrad().M()));      
    p++;
  }
  // end for(p)    

  return curve;
}

//Emulates missile travel. Returns distance travelled.
fixed_t FCajunMaster::FakeFire (AActor *source, AActor *dest, ticcmd_t *cmd)
{
	AActor *th = Spawn ("CajunTrace", source->PosPlusZ(4*8*FRACUNIT), NO_REPLACE);
	
	th->target = source;		// where it came from

	float speed = (float)th->Speed;

	TVector3<double> velocity = source->Vec3To(dest);
	velocity.MakeUnit();
	th->velx = FLOAT2FIXED(velocity[0] * speed);
	th->vely = FLOAT2FIXED(velocity[1] * speed);
	th->velz = FLOAT2FIXED(velocity[2] * speed);

	fixed_t dist = 0;

	while (dist < SAFE_SELF_MISDIST)
	{
		dist += th->Speed;
		th->Move(th->velx, th->vely, th->velz);
		if (!CleanAhead (th, th->X(), th->Y(), cmd))
			break;
	}
	th->Destroy ();
	return dist;
}

//________________________________________________________
void KVParticle::SetRandomMomentum(Double_t T, Double_t thmin,
                                   Double_t thmax, Double_t phmin,
                                   Double_t phmax, Option_t* opt)
{
   //Give randomly directed momentum to particle with kinetic energy T
   //Direction will be between (thmin,thmax) [degrees] limits in polar angle,
   //and (phmin,phmax) [degrees] limits in azimuthal angle.
   //
   //If opt = "" or "isotropic" (default) : direction is isotropically distributed over the solid angle
   //If opt = "random"                    : direction is randomly distributed over solid angle
   //
   //Based on KVPosition::GetRandomDirection().

   Double_t p = (T + M()) * (T + M()) - M2();
   if (p > 0.)
      p = (TMath::Sqrt(p));     // calculate momentum
   else
      p = 0.;

   TVector3 dir;
   KVPosition pos(thmin, thmax, phmin, phmax);
   dir = pos.GetRandomDirection(opt);   // get isotropic unit vector dir
   if (p && dir.Mag())
      dir.SetMag(p);            // set magnitude of vector to momentum required
   SetMomentum(dir);            // set momentum 4-vector
}

TVectorD EUTelState::getStateVec() const { 
	streamlog_out( DEBUG1 ) << "EUTelState::getTrackStateVec()------------------------BEGIN" << std::endl;
	TVector3 momentum =	computeCartesianMomentum();
	TVectorD stateVec(5);
	const float lambda = asin(momentum[2]/(momentum.Mag()));//This will be in radians.
	const float phi = asin(momentum[1]/(momentum.Mag())*cos(lambda));


	stateVec[0] = getOmega();
	stateVec[1] = getIntersectionLocalXZ();
	stateVec[2] = getIntersectionLocalYZ(); 
	stateVec[3] = getPosition()[0]; 
	stateVec[4] = getPosition()[1];
			
//	if ( streamlog_level(DEBUG0) ){
//		streamlog_out( DEBUG0 ) << "Track state:" << std::endl;
//		stateVec.Print();
//	}
	if(stateVec[0] == INFINITY or stateVec[0] == NAN ){
		throw(lcio::Exception( Utility::outputColourString("Passing a state vector where curvature is not defined","RED"))); 
	}

	streamlog_out( DEBUG1 ) << "EUTelState::getTrackStateVec()------------------------END" << std::endl;
 	return stateVec;
}

void Mpdshape::SetPosition(TVector3 vec) {
  ostringstream o;
  o.precision(6);
  o.setf(ios::showpoint); 
  o.setf(ios::fixed);
  o << vec.x() << " " << vec.y() << " " << vec.z();
  fPosition.append(o.str().c_str());
}

double CosThetaStar(TLorentzVector p1, TLorentzVector p2){
	TLorentzVector p = p1 + p2;
	TVector3 theBoost = p.BoostVector();
	TVector3 bostDir;
	if ( theBoost.Mag() != 0 ) bostDir = theBoost.Unit(); // / theBoost.Mag());
	else return -1;
	p1.Boost(-theBoost);
	if (p1.Vect().Mag()!=0) return p1.Vect().Dot(bostDir) / p1.Vect().Mag();
	else return -1;	
}

//_____________________________________________________________________________
void THaSpectrometer::LabToTransport( const TVector3& vertex, 
				      const TVector3& pvect,
				      TVector3& tvertex, Double_t* ray ) const
{
  // Convert lab coordinates to TRANSPORT coordinates in the spectrometer
  // coordinate system.
  // Inputs:
  //  vertex:  Reaction point in lab system
  //  pvect:   Momentum vector in lab
  // Outputs:
  //  tvertex: The vertex point in the TRANSPORT system, without any
  //           coordinate projections applied
  //  ray:     The TRANSPORT ray according to TRANSPORT conventions.
  //           This is an array of size 6 with elements x, tan(theta),
  //           y, tan(y), z, and delta.
  //           z is set to 0, and accordingly x and y are the TRANSPORT 
  //           coordinates in the z=0 plane. delta is computed with respect 
  //           to the present spectrometer's central momentum.
  //           Units are the same as of the input vectors.

  tvertex = fToTraRot * ( vertex - fPointingOffset );
  TVector3 pt = fToTraRot * pvect;
  if( pt.Z() != 0.0 ) {
    ray[1] = pt.X() / pt.Z();
    ray[3] = pt.Y() / pt.Z();
    // In the "ray", project the vertex to z=0
    ray[0] = tvertex.X() - tvertex.Z() * ray[1];
    ray[2] = tvertex.Y() - tvertex.Z() * ray[3];
  } else
    ray[0] = ray[1] = ray[2] = ray[3] = 0.0;
  
  ray[4] = 0.0;   // By definition for this ray, TRANSPORT z=0
  ray[5] = pt.Mag() / fPcentral - 1.0;
}

  static TVector3 floorXYZ(double zvtx, double physeta, double physphi) {
    // Function to calculate detector position when particle is extrapolated
    // to 3rd layer of em calorimeter.

    // Note that we assume an x,y vertex position of zero.
    // Also, the x,y cal shifts are not implemented...yet.

    TVector3 v;
    
    double phi = physphi;  
    double eta = physeta;
    double theta = PTools::etaToTheta(eta);
    double ztmp = zvtx + CC_3R/TMath::Tan(theta);
    
    if( TMath::Abs(ztmp) >= EC_DIV ){ // in EC
      (eta>=0) ? v.SetZ(EC_3Z_SOUTH) : v.SetZ(EC_3Z_NORTH);
      v.SetX( (v.Z()-zvtx) * TMath::Tan(theta) * TMath::Cos(phi));
      v.SetY( (v.Z()-zvtx) * TMath::Tan(theta) * TMath::Sin(phi));
      
    }
    else{ // in CC
      v.SetZ(ztmp);
      v.SetX( CC_3R * TMath::Cos(phi));
      v.SetY( CC_3R * TMath::Sin(phi));
      
    }
    
    return v;
  }

//plane class constructor
Plane::Plane(TVector3 nT){
  
  //Use TVector3 to find an orthogonal vector and a second vector orthogonal to the first and nT
  v1 = nT.Orthogonal();  v2 = nT.Cross(v1);

  //Normalize, checking for 0 length axes
  if ((v1(0) == 0) && (v1(1) == 0) && (v1(2) == 0)){  v1(0) = 0;    v1(1) = 0;    v1(2) = 0; }
  else { mag1 = v1.Mag();   v1(0) = v1(0)/mag1;    v1(1) = v1(1)/mag1;    v1(2) = v1(2)/mag1; } 
  if ((v2(0) == 0) && (v2(1) == 0) && (v2(2) == 0)){  v2(0) = 0;    v2(1) = 0;    v2(2) = 0; } 
  else { mag2 = v2.Mag();   v2(0) = v2(0)/mag2;    v2(1) = v2(1)/mag2;    v2(2) = v2(2)/mag2; }
}//end plane constructor

TVector3 KVParticle::GetVelocity() const
{
   //returns velocity vector in cm/ns units
   TVector3 beta;
   if (E()) {
      beta = GetMomentum() * (1. / E());
   } else {
      beta.SetXYZ(0, 0, 0);
   }
   return (kSpeedOfLight * beta);
}

Double_t getShiftLL(const Double_t* thetaPhi) {
   
   Double_t chi2(1), c(0), oo(0);

   TVector3 norm; norm.SetMagThetaPhi(1.0, thetaPhi[0], thetaPhi[1]);
   for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
      for (Int_t xc=0; xc<ch; ++xc) {
         Double_t dtcor=0;
         for (Int_t i=(ch-xc); i>0; --i) {
            dtcor += dtCorrs[ch-i];
         }
         const TVector3& posCh = getStnPos(ch);
         const TVector3& posXc = getStnPos(xc);
         const Double_t  disCh = -(posCh.Dot(norm));
         const Double_t  disXc = -(posXc.Dot(norm));
         // !!! check sign of delta(distance) and dtcor!
         Double_t  dt = 
            ( (disCh-disXc) * kNgTopFern / kC_m_ns ) // from m to ns
            + dtcor;                                 // correct dt offset (ns)
         const Double_t odt = dt;
         Bool_t oob=kFALSE;
         if (dt<-maxdt) {
            dt  = -maxdt;
            oob = kTRUE;
         } else if (dt>maxdt) {
            dt  = maxdt;
            oob = kTRUE;
         }
         dt *= kSmpRate;
         c = getProbFromCorrCoef(gspl[ch][xc]->Eval(dt));
         if (oob) {
            const Double_t wa = TMath::Abs(odt) - maxdt;
            oo += wa*wa;
         }
         /*
         if (TMath::Abs(dt)>maxdt) {
            // really don't like being out of bounds
            c = TMath::Exp(-dt*dt);
         } else {
            // get the correlation coefficient for this dt (in num samples)
            const Double_t  corco = gspl[ch][xc]->Eval(dt);
            c = getProbFromCorrCoef(corco);
         }
         if (c>1.0) {
            Fatal("getShiftLL","Got ll term > 1 (%g)",c);
         }
         */
         chi2 *= c;
      }
   }
   chi2 = -TMath::Log(chi2);
   chi2 += oo;
   return chi2;
}

Double_t KVMaterial::GetEffectiveAreaDensity(TVector3& norm,
      TVector3& direction)
{
   // Calculate effective area density of absorber (in g/cm**2) as 'seen' in 'direction', taking into
   // account the arbitrary orientation of the 'norm' normal to the material's surface

   TVector3 n = norm.Unit();
   TVector3 d = direction.Unit();
   //absolute value of scalar product, in case direction is opposite to normal
   Double_t prod = TMath::Abs(n * d);
   return GetAreaDensity() / TMath::Max(prod, 1.e-03);
}