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

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

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

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

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

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

// This is the deltaphi between the di-lepton system and the Higgs
// boost in the approximate Higgs rest frame, or R-FRAME
// L1 and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the 
// leptons for you
double HWWKinematics::CalcDoubleDphiRFRAME(){
  // first calculate pt-corrected MR
  float mymrnew = CalcMRNEW();
  
  TVector3 BL = L1.Vect()+L2.Vect();
  BL.SetX(0.0);
  BL.SetY(0.0);
  BL = (1./(L1.P()+L2.P()))*BL;
  L1.Boost(-BL);
  L2.Boost(-BL);
  
  //Next, calculate the transverse Lorentz transformation
  TVector3 B = L1.Vect()+L2.Vect()+MET;
  B.SetZ(0.0);
  B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
  
  L1.Boost(B);
  L2.Boost(B);
  
  // Now, calculate the delta phi
  // between di-lepton axis and boost
  // in new reference frame
  
  return B.DeltaPhi(L1.Vect()+L2.Vect());
  
}

// This is the pt corrected delta phi between the 2 leptons
// P and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the 
// leptons for you
double HWWKinematics::CalcDeltaPhiRFRAME(){
  // first calculate pt-corrected MR
  float mymrnew = CalcMRNEW();
  
  // Now, boost lepton system to rest in z
  // (approximate accounting for longitudinal boost)
  TVector3 BL = L1.Vect()+L2.Vect();
  BL.SetX(0.0);
  BL.SetY(0.0);
  BL = (1./(L1.P()+L2.P()))*BL;
  L1.Boost(-BL);
  L2.Boost(-BL);
  
  // Next, calculate the transverse Lorentz transformation
  // to go to Higgs approximate rest frame
  TVector3 B = L1.Vect()+L2.Vect()+MET;
  B.SetZ(0.0);
  B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
  
  L1.Boost(B);
  L2.Boost(B);
  
  //Now, re-calculate the delta phi
  // in the new reference frame:
  return L1.DeltaPhi(L2);
  
}

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

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

// M is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
double HWWKinematics::CalcMRNEW(TLorentzVector P, TLorentzVector Q, TVector3 M){
 TVector3 vI = M+P.Vect()+Q.Vect();
 vI.SetZ(0.0);
 double PpQ = CalcMR(P,Q); //Note - this calls the old MR function
 double vptx = (P+Q).Px();
 double vpty = (P+Q).Py();
 TVector3 vpt;
 vpt.SetXYZ(vptx,vpty,0.0);

 float MR2 = 0.5*(PpQ*PpQ-vpt.Dot(vI)+PpQ*sqrt(PpQ*PpQ+vI.Dot(vI)-2.*vI.Dot(vpt)));

 return sqrt(MR2);

}

// This is the pt corrected delta phi between the 2 mega-jets
// P and Q are the 4-vectors for the 2 hemispheres 
// M is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the 
// leptons for you
double HWWKinematics::CalcDeltaPhiNEW(TLorentzVector P, TLorentzVector Q, TVector3 M){
    // first calculate pt-corrected MR
 float mymrnew = CalcMRNEW(L1,L2,MET);

    //Next, calculate the transverse Lorentz transformation
 TVector3 B = P.Vect()+Q.Vect()+MET;
 B.SetZ(0.0);
 B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;

 P.Boost(B);
 Q.Boost(B);

    //Now, re-calculate the delta phi
    // in the new reference frame:

 return P.DeltaPhi(Q);

}

double HWWKinematics::CalcUnboostedMTR(TLorentzVector P, TLorentzVector Q, TVector3 M){
    // first calculate pt-corrected MR
 float mymrnew = CalcMRNEW(L1,L2,MET);

    //Next, calculate the transverse Lorentz transformation
 TVector3 B = P.Vect()+Q.Vect()+MET;
 B.SetZ(0.0);
 B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;

 P.Boost(B);
 Q.Boost(B);

    //Now, re-calculate MTR in the new reference frame:
 float mymtrnew = CalcMTRNEW(P, Q);

    //R is now just the ratio of mymrnew and mymtrnew;

 return mymtrnew;
}

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

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