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

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

/*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 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);

		  }
	  }

void THisHandler::FillHists()
{
   fTree->SetBranchStatus("*", 0);

   TVector3 position;
   fTree->SetBranchStatus("x", 1);
   fTree->SetBranchStatus("y", 1);
   fTree->SetBranchStatus("z", 1);
   fTree->SetBranchAddress("x", &position[0]);
   fTree->SetBranchAddress("y", &position[1]);
   fTree->SetBranchAddress("z", &position[2]);
   
   Char_t volumeName[256]; // Array size is enough?
   fTree->SetBranchStatus("VolumeName", 1);
   fTree->SetBranchAddress("VolumeName", volumeName);

   Double_t ene;
   fTree->SetBranchStatus("DepositEnergy", 1);
   fTree->SetBranchAddress("DepositEnergy", &ene);

   const Int_t kEntries = fTree->GetEntries();
   for(Int_t iEntry = 0; iEntry < kEntries; iEntry++){
      //if(iEntry%100000 == 0) cout << iEntry <<" / "<< kEntries << endl;
      fTree->GetEntry(iEntry);
      if(!(ene > 0.)) continue;
      
      fHisAll->Fill(position.X(), position.Y(), position.Z(), ene);

      if(TString(volumeName) == "Sealing")
         fHisSealing->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Window")
         fHisWindow->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Foil")
         fHisFoil->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Holder")
         fHisHolder->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Cassette")
         fHisCassette->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Air")
         fHisAir->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Plate" ||
              TString(volumeName) == "Well" ||
              TString(volumeName) == "Outer")
         fHisPlate->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Film")
         fHisFilm->Fill(position.X(), position.Y(), position.Z(), ene);
      else if(TString(volumeName) == "Stuff")
         fHisWell->Fill(position.X(), position.Y(), position.Z(), ene);

      else if(TString(volumeName) == "Cell"){
         Int_t wellID = XYtoIndex(position.X(), position.Y());
         fHisEachCell20[wellID]->Fill(position.X(), position.Y(), ene);
         fHisEachCell50[wellID]->Fill(position.X(), position.Y(), ene);
      }

   }

}

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

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

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

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;

}

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

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

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

void drawFocalPlane(TString infile="../build/hits.root", Double_t r1 = 48.8, Double_t r2 = 29.1, Int_t it1=0, Int_t it2=0, Double_t energy=3){
  gSystem->Load("../build/libprtdirclib.so");
  PrtInit(infile,0);
  gStyle->SetOptStat(0);

  Double_t radiatorL(1250); //bar
  // Double_t radiatorL(1224.9); //plate
  TVector3 res;
  TH2F *hFp1 = new TH2F("hFp1",Form("r_{1}=%2.2f    r_{2}=%2.2f;x, [cm];y, [cm]",r1,r2),500,0,50,500,0,30 );
  TH2F *hFp2 = new TH2F("hFp2",Form("r_{1}=%2.2f    r_{2}=%2.2f;z, [cm];y, [cm]",r1,r2),500,-30,30,500,-30,50 );

  PrtHit hit[2];
  Int_t nevents = fCh->GetEntries();
  Int_t count(0);
  for (Int_t ievent=0; ievent<nevents; ievent++){
    PrtNextEvent(ievent,1000);
    Int_t nhits = prt_event->GetHitSize(); //fEvent->GetHitSize();
    if(nhits!=2) continue;
    for(Int_t h=0; h<nhits; h++) hit[h] = prt_event->GetHit(h); //fEvent->GetHit(h);
    Double_t d = findVertex(hit[0].GetGlobalPos(),hit[0].GetMomentum().Unit(),hit[1].GetGlobalPos(),hit[1].GetMomentum().Unit(), &res);
    if(d<1){
      Double_t x = -(res.X()+radiatorL/2.)/10.;
      if(x>6){
	hFp1->Fill(x,res.Z()/10.);
	hFp2->Fill(res.Y()/10.,res.Z()/10.);
      }
      count++;
    }
  } 
  Double_t eff = 100*count/(Double_t)nevents;

  TString senergy = "";
    if(energy>-1){
      Double_t m = 1, cm = 0.01, nanometer = 0.000000001, GeV = 1000000000;
      Double_t LambdaE = 2.0 * 3.14159265358979323846 * 1.973269602e-16 * m * GeV;
      Double_t lam = LambdaE/(energy*nanometer);

      senergy = Form("   E=%2.2f [eV] / #lambda=%2.0f [nm]",energy,lam);
    }
  hFp1->SetTitle(Form("r_{1}=%2.2f    r_{2}=%2.2f   #varepsilon=%2.0f",r1,r2,eff)+senergy);
  canvasAdd(Form("fp_%d_%d",it1,it2),800,500);
  hFp1->Draw("colz");
  canvasAdd(Form("fp2_%d_%d",it1,it2),600,800);
  hFp2->Draw("colz");
  //canvasSave(0,"drawFocalPlane.C",1,"fp");
  
}

  static Bool_t is_in_eta_fiducial(const TVector3 &pos) {
    // converted from emid::is_in_eta_fiducial

    // NOTE: this routine assumes pos is in local coordinates!
    //convert pos to local coordinate system
    //TVector3 v = toLocal(pos);

    // cryostat edges
    bool is_at_edge;
    if (TMath::Abs(pos.Z())>150.) {
      is_at_edge=pos.Perp()>EC_R_MAX;
    } else {
      is_at_edge=fabs(pos.Z())>CC_Z_MAX;
    }
    return (!is_at_edge);
    
  }

//_________________________________________________________________
Double_t KVTenseur3::GetPhiPlan(void)
{
   // Obtention du Phi du plan de reaction (en degrès).
   // Cet angle suit les conventions de KaliVeda i.e. il est entre 0 et 360 degrees
   if (is_diago == 0)
      Diago();

   TVector3 vp = GetVep(3);

   if (vp.Z() < 0)
      vp = -vp;
   Double_t phi = vp.Phi() * TMath::RadToDeg();
   return (phi < 0 ? 360. + phi : phi);
}