C++ TVector3::Y方法代码示例

	  void qaVtxDiff(TString pre="", RhoCandidate * c, RhoTuple * n){

		  if(n==0) return;
		  if(c==0) return;

		  RhoCandidate * mct = c->GetMcTruth();

		  if(mct){
			  TVector3 v = c->DecayVtx();
			  TVector3 mcv = mct->Daughter(0)->Pos();
			  TVector3 vdiff = v-mcv;
			  TMatrixD cov7 = c->Cov7();

			  n->Column(pre + "diffvx", (Float_t) vdiff.X(), 0.0f );
			  n->Column(pre + "diffvy", (Float_t) vdiff.Y(), 0.0f );
			  n->Column(pre + "diffvz", (Float_t) vdiff.Z(), 0.0f );

			  n->Column(pre + "pullvx", (Float_t) (vdiff.X()/TMath::Sqrt(cov7(0,0))), 0.0f);
			  n->Column(pre + "pullvy", (Float_t) (vdiff.Y()/TMath::Sqrt(cov7(1,1))), 0.0f);
			  n->Column(pre + "pullvz", (Float_t) (vdiff.Z()/TMath::Sqrt(cov7(2,2))), 0.0f);

		  }
		  else{
			  n->Column(pre + "diffvx", (Float_t) -999.0, 0.0f );
			  n->Column(pre + "diffvy", (Float_t) -999.0, 0.0f );
			  n->Column(pre + "diffvz", (Float_t) -999.0, 0.0f );

			  n->Column(pre + "pullvx", (Float_t) -999.0, 0.0f);
			  n->Column(pre + "pullvy", (Float_t) -999.0, 0.0f);
			  n->Column(pre + "pullvz", (Float_t) -999.0, 0.0f);

		  }
	  }

/*Given the impact point on the front face (vin) and the incoming particle LorentzVector (chi for invisible decay, A' for visible),
 *  determine the interaction point within the fiducial volume and save it in vhit.
 Use a random distribution along the chi flight path, with uniform probability
 This function returns the length (in m) of the trajectory within the fiducial volume.
 Displacement is the lateral displacement (in m) of the detector, along x
 */
double KinUtils::findInteractionPoint(const TLorentzVector &chi, const TVector3 &fiducialV, const TVector3 &vin, TVector3 &vout, TVector3 &vhit) {

	double tz, tx, ty, tout, L;
	int sigPx, sigPy;

	sigPx = (chi.Px() > 0 ? 1 : -1);
	sigPy = (chi.Py() > 0 ? 1 : -1);

	tz = fiducialV.Z() / chi.Pz();
	tx = (sigPx * fiducialV.X() / 2 - vin.X()) / chi.Px();
	ty = (sigPy * fiducialV.Y() / 2 - vin.Y()) / chi.Py();
	tout = 0;

	if ((tz < tx) && (tz < ty)) {
		tout = tz;
	} else if ((tx < ty) && (tx < tz)) {
		tout = tx;
	} else if ((ty < tx) && (ty < tz)) {
		tout = ty;
	}
	vout.SetXYZ(tout * chi.Px() + vin.X(), tout * chi.Py() + vin.Y(), tout * chi.Pz() + vin.Z());
	vhit.SetXYZ(Rand.Uniform(vin.X(), vout.X()), Rand.Uniform(vin.Y(), vout.Y()), Rand.Uniform(vin.Z(), vout.Z()));
	L = (vout - vin).Mag();

	return L;
}

//_____________________________________________________________________________
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;
}

TVector3 LineCrossesCircle(TVector3 CircleCenter, Double_t CircleRadius,
			   TVector3 LineStart, TVector3 LineEnd) {

  TVector3 LineVec = LineEnd - LineStart;
  TVector3 VecToLine = LineStart - CircleCenter;
  Double_t a, b, c;
  Double_t sqrtterm, res1, res2, res;
  a = LineVec.X()*LineVec.X() + LineVec.Y()*LineVec.Y();
  b = 2 * ( VecToLine.X()*LineVec.X() + VecToLine.Y()*LineVec.Y() );
  c = ( VecToLine.X()*VecToLine.X() + VecToLine.Y()*VecToLine.Y()) -
    CircleRadius*CircleRadius;
  sqrtterm = b*b - 4*a*c;
  if(sqrtterm < 0) return LineStart;
  sqrtterm = TMath::Sqrt(sqrtterm);
  res1 = ( -b - sqrtterm ) / (2 * a);
  res2 = ( -b + sqrtterm ) / (2 * a);
  if(res1 >= 0 && res1 <= 1) {
    res = res1;
  } else {
    res = res2;
  }
  TVector3 result(LineStart.X() + res * (LineEnd.X() - LineStart.X()),
		  LineStart.Y() + res * (LineEnd.Y() - LineStart.Y()), 0);            
  return result;
} 

double NonCollisionBG::tDCA(const TVector3& X, const TVector3& P){
 
  // get position, and momentum components
  double x0 = X.X();
  double y0 = X.Y();
  double px = P.X();
  double py = P.Y();

  // solve for t of closest approach to z-axis:
  double t_closestZ = -(y0*py + x0*px)/(px*px + py*py);
	
  return t_closestZ;
}

int getSurfacePoint(TVector3 intPos, TVector3 &intDir, TVector3 &surfPos) {
    Double_t b=intDir.X()*intPos.X() + intDir.Y()*intPos.Y() + intDir.Z()*intPos.Z();
    Double_t c = intPos.X()*intPos.X() + intPos.Y()*intPos.Y() + intPos.Z()*intPos.Z() - rEarth*rEarth;
    if(b*b < c)
        return 0;


    Double_t l1=-1*b + TMath::Sqrt(b*b - c);
    Double_t l2=-1*b - TMath::Sqrt(b*b - c);


    if(intPos.Mag2()> rEarth*rEarth) {
        //Start outside take l1
        //return 0;
        //    cout << "Outside:\t" << l1 << "\t" << l2 << endl;
        intDir*=-1;
        l1*=-1;
        surfPos.SetX(intPos.X() + l1*intDir.X());
        surfPos.SetY(intPos.Y() + l1*intDir.Y());
        surfPos.SetZ(intPos.Z() + l1*intDir.Z());
        //    cout << surfPos.Mag() << endl;
    }
    else {
        //Start inside take l2
        //    cout << "Inside:\t" << l1 << "\t" << l2 <<  endl;
        //    cout << l2 << endl;
        //    return 0;
        surfPos.SetX(intPos.X() + l2*intDir.X());
        surfPos.SetY(intPos.Y() + l2*intDir.Y());
        surfPos.SetZ(intPos.Z() + l2*intDir.Z());
        //    cout << surfPos.Mag() << endl;
    }

    return 1;
}

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;

}

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);
}

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;
}

vector <double> NonCollisionBG::innerIntersectionTs(double R, const TVector3& X, const TVector3& P){  
  double x0 = X.X();
  double y0 = X.Y();
  double px = P.X();
  double py = P.Y();
  
  // this b_over_2 = b/2 from eq 3.
  double b_over_2 = x0 * px + y0 * py ;
  
  // this thingUnderSqrt = (b/2)^2 - a c from eq 3.
  double thingUnderSqrt = b_over_2 * b_over_2 - (px*px + py*py) * (x0*x0 + y0*y0 - R*R) ;
  vector <double> t;
  t.push_back(   (- x0*px - y0*py - sqrt(thingUnderSqrt)) / (px*px + py*py)  );
  t.push_back(   (- x0*px - y0*py + sqrt(thingUnderSqrt)) / (px*px + py*py)  );
  return t;
  
}

int get_sq_wg_efield(TVector3 &pos, TVector3 &efield)
{
  //function returns the TE_10 efield inside an infinite square waveguide 
  //centered at x=0, y=0 (not with corner at origin as is standard)
  //efield normalized the jackson way with 1/cm units
  //waveguide extends in z-direction
  // wavelength of 27 GHz radiation is 1.1 cm
  double e_amp = sqrt(2 / tl_data.x1 / tl_data.y1);
  int status = 0;

  if ((pos.X() > tl_data.x1 / 2.) || (pos.X() < -tl_data.x1 / 2.)) {
    cout << "Problem!!! Electron hit a wall in x-dir! " << endl;
    status = 1;
  }
  if ((pos.Y() > tl_data.y1 / 2.) || (pos.Y() < -tl_data.y1 / 2.)) {
    cout << "Problem!!! Electron hit a wall in y-dir! " << endl;
    status = 1;
  }
  efield.SetX( 0 );
  efield.SetY( e_amp * sin(M_PI * (pos.X() + tl_data.x1 / 2.) / tl_data.x1) );
  efield.SetZ( 0 );                //is true for TE mode
  return status;
}

  static TVector3 toGlobal(const TVector3 &local){

    TVector3 global = local;

    if(TMath::Abs(local.Z()) < EC_DIV){ //CC
      global.SetX(global.X()+CC_XSHIFT);
      global.SetY(global.Y()+CC_YSHIFT);
      global.SetZ(global.Z()+CC_ZSHIFT);
    }else{
      if (local.Z() < 0) global.SetX(local.X()+EC_XSHIFT_NORTH);
      global.SetZ( (local.Z()>0) ? EC_3Z_SOUTH : EC_3Z_NORTH );
    }

    return global;
  }

  static TVector3 toLocal(const TVector3 &global){

    TVector3 local = global;

    if(TMath::Abs(global.Z()) < EC_DIV){ //CC
      local.SetX(local.X()-CC_XSHIFT);
      local.SetY(local.Y()-CC_YSHIFT);
      local.SetZ(local.Z()-CC_ZSHIFT);
    }else{
      if (global.Z() < 0) local.SetX(local.X()-EC_XSHIFT_NORTH);
      local.SetZ( (global.Z()>0) ? EC_3Z : -EC_3Z );
    }

    return local;
  }

int get_pp_efield(TVector3 &p, TVector3 &efield)
{
  //function returns the electric field between two parallel plates 
  //amplitude corrected for capacitance of finite width strips
  //not valid outside plates
  //efield normalized the jackson way with 1/cm units
  //strips extend in y- and z- directions, so the field is in the x-direction
  double e_amp = sqrt(1 / (tl_data.n * tl_data.l * (tl_data.x2 - tl_data.x1)));
  //factor of sqrt(n=1.8) lower than ideal capacitor 
  int status = 0;

  if ((p.X() < tl_data.x1) || (p.X() > tl_data.x2)) {
    cout << "Problem!!! Electron hit a plate! " << endl;
    status = 1;
  }
  if ((p.Y() < -tl_data.l / 2) || (p.Y() > tl_data.l / 2)) {
    cout << "Problem!!! Electron outside plates! " << endl;
    status = 1;
  }
  efield.SetX( e_amp );
  efield.SetY( 0 );
  efield.SetZ( 0 );                //only true for TE or TEM modes
  return status;
}