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

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

// // //
double analysisClass::visPzeta( unsigned int iMuR, unsigned int iTauR ){
  TVector3 muT; TVector3 tauT; TVector3 unitmuT; TVector3 unittauT; TVector3 unitbisecT; TVector3 MET;
  muT.SetPtEtaPhi(  muPtcorr(iMuR),    0, MuonPhi->at(iMuR)    );
  tauT.SetPtEtaPhi( tauPtcorr(iTauR), 0, HPSTauPhi->at(iTauR) );
  unitmuT=muT*(1./muT.Mag());  unittauT=tauT*(1./tauT.Mag());
  unitbisecT=(unitmuT+unittauT)*(1./((unitmuT+unittauT).Mag()));
  MET.SetPtEtaPhi( METcorr("Pt"), 0, METcorr("Phi")  );
  double pZetaVis;
  pZetaVis = unitbisecT.Dot( (muT+tauT)     );
  return pZetaVis;
}

  //returns a projection onto the 2D plane 
  TVector3 Projection(TVector3 jaxis){
    //Find the projection of a jet onto this subspace

    if(v1.Mag() == 0) { scalar1 = 0; }   else { scalar1 = jaxis.Dot(v1)/(v1.Dot(v1)); } 
    if(v2.Mag() == 0) { scalar2 = 0; }   else { scalar2 = jaxis.Dot(v2)/(v2.Dot(v2)); } 
    v1 = scalar1*v1;
    v2 = scalar2*v2;
    proj(0) = v1(0) + v2(0);
    proj(1) = v1(1) + v2(1);
    proj(2) = v1(2) + v2(2); 
    
    return proj;
  }//end of projection

int get_tl_efield(TVector3 &p,TVector3 &efield)
{
  //function returns the electric field between two infinite parallel wires
  //efield normalized the jackson way with 1/cm units
  //wires extend in z-direction and are can be offset 
  double e_amp = 1. / 2.0 / M_PI * sqrt(tl_data.C / (EPS0 / M2CM));

  //calculate effective wire position for efield
  TVector3 x1(tl_data.x1, 0, 0);	//real position of first wire in cm
  TVector3 x2(tl_data.x2, 0, 0);	//real position of second wire in cm
  TVector3 l = x1 - x2;
  l *= 1./2;
  double a = sqrt(l.Mag2() - pow(tl_data.rI, 2));	//effective wire half separation in cm
  TVector3 a1 = x1 - l + a*l.Unit();	//effective position of first wire in cm
  TVector3 a2 = x1 - l - a*l.Unit();	//effective position of second wire in cm
 
  //vector from point to wires
  p.SetZ(0);
  TVector3 r1 = p + a1;
  TVector3 r2 = p + a2;
   
  //check for hitting wires 
  int status = 0;
  if ( ( (x1-r1).Mag() < tl_data.rI ) || ( (x2-r2).Mag() < tl_data.rI ) ) {
    cout << "Problem!!! Electron hit a wire! Radius " << p.Mag() << endl;
    status = 1;
  }
  //efield = e_amp * (r1 * 1/r1.Mag2() - r2 * 1/r2.Mag2());
  efield.SetX( e_amp * ( r1.X() / r1.Mag2() - r2.X() / r2.Mag2() ) );
  efield.SetY( e_amp * ( r1.Y() / r1.Mag2() - r2.Y() / r2.Mag2() ) );
  efield.SetZ( 0 );                //only true for TE or TEM modes
  return status;
}

int get_circ_wg_efield_vert(TVector3 &pos, TVector3 &efield)
{
  //returns the TE_11 efield inside an infinite circ. waveguide 
  //efield normalized the jackson way with 1/cm units
  //waveguide extends in x-direction
  // wavelength of 27 GHz radiation is 1.1 cm
  double p11 = 1.841;            //1st zero of the derivate of bessel function
  //k11 is angular wavenumber for cutoff frequency for TE11
  double k11 = p11 / tl_data.rO;
  //convert position to cylindrical
  pos.SetZ(0);
  double radius = pos.Mag();
  double phi = pos.Phi();//azimuthal position
  //double phi = acos(pos.Y()/radius);//azimuthal position, see definition of phase
  double J1 = TMath::BesselJ1(k11 * radius);//this term cancels in dot product w/ vel
  double Jp = TMath::BesselJ0(k11 * radius) - TMath::BesselJ1(k11 * radius) / k11 / radius;
  double e_amp = 1.63303/tl_data.rO;//cm
  int status = 0;

  if (radius > tl_data.rO) {
    cout << "Problem!!! Electron hit a wall! At radius " << radius << endl;
    status = 1;
  }
  efield.SetX(e_amp * (J1 / k11 / radius * cos(phi) * sin(phi) - Jp * sin(phi) * cos(phi)));
  efield.SetY(e_amp * (J1 / k11 / radius * sin(phi) * sin(phi) + Jp * cos(phi) * cos(phi)));
  efield.SetZ(0);                //only true for TE mode
  if (radius == 0) {
    Jp = 1.0/2; 
    phi = 0;
    efield.SetX(e_amp * (1. / 2 * cos(phi) * sin(phi) - Jp * sin(phi) * cos(phi)));
    efield.SetY(e_amp * (1. / 2 * sin(phi) * sin(phi) + Jp * cos(phi) * cos(phi)));
  }
  return status;
}

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

//________________________________________________________
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
}

int get_coax_efield(TVector3 &p, TVector3 &efield)
{
  //function returns the electric field inside an infinite coaxial cable 
  //efield normalized the jackson way with 1/cm units
  //cable extend in z-direction 
  double e_amp = 1 / sqrt(2.0 * M_PI * log(tl_data.rO / tl_data.rI));
  int status = 0;

  p.SetZ(0);
  if (p.Mag() < tl_data.rI || p.Mag() > tl_data.rO) {
    cout << "Problem!!! Electron hit the cable! " << endl;
    status = 1;
  }
  efield.SetX( e_amp * p.X() / p.Mag2() );
  efield.SetY( e_amp * p.Y() / p.Mag2() );
  efield.SetZ( 0 );                //only true for TE or TEM modes
  return status;
}

// // //
double analysisClass::deltaPzeta( unsigned int iMuR, unsigned int iTauR ){
  TVector3 muT; TVector3 tauT; TVector3 unitmuT; TVector3 unittauT; TVector3 unitbisecT; TVector3 MET;
  muT.SetPtEtaPhi(  muPtcorr(iMuR),    0, MuonPhi->at(iMuR)    );
  tauT.SetPtEtaPhi( tauPtcorr(iTauR), 0, HPSTauPhi->at(iTauR) );
  unitmuT=muT*(1./muT.Mag());  unittauT=tauT*(1./tauT.Mag());
  unitbisecT=(unitmuT+unittauT)*(1./((unitmuT+unittauT).Mag()));
  MET.SetPtEtaPhi( METcorr("Pt"), 0, METcorr("Phi")  );
  double pZeta;  double pZetaVis;
  pZeta    = unitbisecT.Dot( (muT+tauT+MET) );
  pZetaVis = unitbisecT.Dot( (muT+tauT)     );
  /*std::cout<<" TauMag, MuMag: "<<unitmuT.Mag()<<" "<<unittauT.Mag()<<std::endl;
  std::cout<<"  Mu-bisec: "<<muT.DeltaPhi(unitbisecT)<<std::endl;
  std::cout<<" Tau-bisec: "<<tauT.DeltaPhi(unitbisecT)<<std::endl;
  std::cout<<"    Tau-Mu: "<<tauT.DeltaPhi(muT)<<std::endl;
  std::cout<<"(pZeta-1.5*pZetaVis): "<<(pZeta-1.5*pZetaVis)<<std::endl;
  std::cout<<std::endl;*/
  return (pZeta-1.5*pZetaVis);
}

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

Double_t LineLength(Double_t alpha, TVector3 CircleCenter,
		    Double_t CircleRadius, TVector3 LineStart) {                   

  Double_t yend = CircleRadius*1.1;
  Double_t xend = yend*TMath::Tan(TMath::DegToRad()*alpha);
  TVector3 LineEnd(xend, yend, 0);
  TVector3 rend = LineCrossesCircle(CircleCenter, CircleRadius, LineStart, LineEnd);
  if (rend != LineStart) {
    TVector3 llen = rend - LineStart;
    return llen.Mag();
  }
  return 0;
} 

//find
int EUTelState::findIntersectionWithCertainID(int nextSensorID, float intersectionPoint[], TVector3& momentumAtIntersection, float& arcLength ){
	streamlog_out(DEBUG5) << "-EUTelState::findIntersectionWithCertainID---------------------------BEGIN" << std::endl;
	TVector3 pVec = computeCartesianMomentum();
	streamlog_out(DEBUG5) << "Momentum (Global): " << pVec[0]<<","<<pVec[1]<<","<<pVec[2]<<","<<" Position (local): "<<getPosition()[0]<<","<<getPosition()[1]<<","<<getPosition()[2]<< std::endl;
	if(pVec.Mag() == 0){
		throw(lcio::Exception( Utility::outputColourString("The momentum is 0","RED"))); 
	}
	double posLocal[] =  {getPosition()[0],getPosition()[1],getPosition()[2] };
	double temp[] = {0.,0.,0.};
	geo::gGeometry().local2Master(getLocation() , posLocal, temp);//IMPORTANT:For strip sensors this will make the hit strip look like a pixel at (Xstriplocal,somevalue,somevalue).
	float posGlobal[] = { static_cast<float>(temp[0]), static_cast<float>(temp[1]), static_cast<float>(temp[2]) };
	int sensorID = geo::gGeometry().findIntersectionWithCertainID(posGlobal[0] ,  posGlobal[1] , posGlobal[2], pVec[0],pVec[1],pVec[2], getBeamCharge(), nextSensorID, intersectionPoint, momentumAtIntersection, arcLength ); 
	streamlog_out(DEBUG5) << "-EUTelState::findIntersectionWithCertainID--------------------------END" << std::endl;
	return sensorID;
}

//This sets the LOCAL frame intersection. Not the curvilinear frames intersection
void EUTelState::setLocalXZAndYZIntersectionAndCurvatureUsingGlobalMomentum(TVector3 momentumIn){
	streamlog_out(DEBUG5) << "EUTelState::setLocalXZAndYZIntersectionAndCurvatureUsingGlobalMomentum--------------------------BEGIN" << std::endl;

	//set the beam energy and curvature
	setBeamEnergy(momentumIn.Mag());
	initialiseCurvature();//You must set beam charge before you call this.
	//Now calculate the momentum in LOCAL coordinates.
	const double momentum[]	= {momentumIn[0], momentumIn[1],momentumIn[2]};//Need this since geometry works with const doubles not floats 
	double localMomentum [3];
	geo::gGeometry().master2LocalVec(getLocation(), momentum, localMomentum );
	//In the LOCAL coordinates this is just dx/dz and dy/dz in the LOCAL frame
	streamlog_out(DEBUG5) << "The local momentum (x,y,z) is: "<< localMomentum[0]<<","<< localMomentum[1] <<"," <<localMomentum[2] << std::endl;
	setIntersectionLocalXZ(localMomentum[0]/localMomentum[2]);
	setIntersectionLocalYZ(localMomentum[1]/localMomentum[2]);
	streamlog_out(DEBUG5) << "The XZ tilt is: "<< getIntersectionLocalXZ()<<" The YZ tilt is: "<<  getIntersectionLocalYZ()<< std::endl;
	streamlog_out(DEBUG5) << "EUTelState::setLocalXZAndYZIntersectionAndCurvatureUsingGlobalMomentum--------------------------END" << std::endl;
}

//find
bool EUTelState::findIntersectionWithCertainID(int nextSensorID, float intersectionPoint[], TVector3& momentumAtIntersection, float& arcLength, int& newNextPlaneID )
{
	TVector3 pVec = getMomGlobal();

	//TODO: better exception
	if(pVec.Mag() == 0)
	{
		throw(lcio::Exception( "The momentum is 0")); 
	}

	double posLocal[] =  {getPosition()[0],getPosition()[1],getPosition()[2] };
	double temp[] = {0.,0.,0.};

	//IMPORTANT:For strip sensors this will make the hit strip look like a pixel at (Xstriplocal,somevalue,somevalue).
	geo::gGeometry().local2Master(getLocation() , posLocal, temp);
    int charge = -1;
	float posGlobal[] = { static_cast<float>(temp[0]), static_cast<float>(temp[1]), static_cast<float>(temp[2]) };
	return  EUTelNav::findIntersectionWithCertainID(	posGlobal[0], posGlobal[1], posGlobal[2], 
								pVec[0], pVec[1], pVec[2], charge,
								nextSensorID, intersectionPoint, 
								momentumAtIntersection, arcLength, newNextPlaneID); 
}

//_________________________________________________________________
void KVRelativeVelocity::Fill(KVNucleus* cc)
{

//
// Filling method
//
//
// Use the FillVar(val) or the FillVar(val,w) method to increment the quantity
// of the global variable.
// The weight w is evaluated from the properties of the KVNucleus passed as an
// argument.
//
// For example, to evaluate the mean charge of all fragments,
// you may proceed as follows:
//
// FillVar(c->GetZ());
//
// Another example: to evaluate the mean parallel velocity weighted by the charge
// of the nucleus:
//
// FillVar(c->GetV().Z(),c->GetZ());
//
   if (!heaviest) {
      heaviest = new TList;
      heaviest->SetOwner(kFALSE);
      heaviest->Add(cc);
   }
   else {
      for (Int_t ii = 0; ii < heaviest->GetEntries(); ii += 1) {
         TVector3 ww = cc->BoostVector() - ((KVNucleus*)heaviest->At(ii))->BoostVector();
         FillVar(ww.Mag());
      }
      heaviest->Add(cc);
   }

}