C++ TLorentzVector::Boost方法代码示例

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

inline TLorentzVector* kinematics::Boost( TLorentzVector* pa, TLorentzVector* pb, TLorentzVector* p )
{
	TLorentzVector pTmp = (*pa)+(*pb);
	TLorentzVector* pBoosted = (TLorentzVector*)p->Clone("");
	pBoosted->Boost(-1.*pTmp.BoostVector());
	return pBoosted;
}

void Tran_to_fixed_target (TLorentzVector e_vec, TLorentzVector d_vec,
                            TLorentzVector &e_ft, TLorentzVector  &d_ft){
   //input: 
   //   e_vec - electron 4 vector in some frame
   //   d_vec - ion 4 vector in some frame
   //output:
   //   e_ft - electron 4 vector in rest frame of ion (fixed target frame)
   //   d_ft - ion 4 vector in its rest frame 
   e_ft = e_vec;
   d_ft = d_vec;
   
   d_ft.Boost(-d_vec.BoostVector()); //should be (0,0,0, mass)
   e_ft.Boost(-d_vec.BoostVector());
   return;


}

void Tran_to_cm(TLorentzVector e_vec, TLorentzVector d_vec, 
                TLorentzVector &e_cm, TLorentzVector &d_cm){
     //input: 
     //   e_vec - electron 4 vector in some frame
     //   d_vec - ion 4 vector in some frame
     //output:
     //   e_ft - electron 4 vector in center of mass frame
     //   d_ft - ion 4 vector in center of mass frame
     e_cm = e_vec;
     d_cm = d_vec;

     TLorentzVector cm_vec = e_vec + d_vec;

     e_cm.Boost(-cm_vec.BoostVector());
     d_cm.Boost(-cm_vec.BoostVector());
     return;

}

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

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 Lorentz()
{
  // gamma + p -> gamma ' + p '
  TLorentzVector proton(0,0,0,0.938);
  TLorentzVector photon(1,0,0,2);
  proton.Print();
  photon.Print();
  TLorentzVector CM=proton+photon;
  CM.Print();
  CM.Vect().Print();
  v= CM.Vect();
  TLorentzVector protonCM;
  TLorentzVector photonCM;
  protonCM=proton;
  protonCM.Boost(-v);
  photonCM=photon;
  photonCM.Boost(-v);
  protonCM.Print();
  photonCM.Print();
}

void Boost_To_Stop_Rest_Frame(TLorentzVector& stop4, TLorentzVector& chargino4, TLorentzVector& b4, TLorentzVector& neutralino4, TLorentzVector& W4, TLorentzVector& up4, TLorentzVector& down4, TLorentzVector& s4)
{
    TVector3 betaV(-stop4.Px()/stop4.Energy(),-stop4.Py()/stop4.Energy(),-stop4.Pz()/stop4.Energy());
    stop4.Boost(betaV);
    chargino4.Boost(betaV);
    b4.Boost(betaV);
    neutralino4.Boost(betaV);
    W4.Boost(betaV);
    up4.Boost(betaV);
    down4.Boost(betaV);
    s4.SetE(chargino4.P()/chargino4.M());
    s4.SetVect(chargino4.Vect().Unit()*chargino4.Gamma());
}

double cosRestFrame(TLorentzVector boost, TLorentzVector vect) {

  double bx = -boost.Px()/boost.E();
  double by = -boost.Py()/boost.E();
  double bz = -boost.Pz()/boost.E();

  vect.Boost(bx,by,bz);
  double prod = -vect.Px()*bx-vect.Py()*by-vect.Pz()*bz;
  double modBeta = TMath::Sqrt(bx*bx+by*by+bz*bz); 
  double modVect = TMath::Sqrt(vect.Px()*vect.Px()+vect.Py()*vect.Py()+vect.Pz()*vect.Pz());
  
  double cosinus = prod/(modBeta*modVect);

  return cosinus;

}

inline float kinematics::cosThetaBoost( TLorentzVector* pa, float ca,
										TLorentzVector* pb, float cb )
{
	// http://xrootd.slac.stanford.edu/BFROOT/www/doc/workbook_backup_010108/analysis/analysis.html
	// A useful quantity in many analyses is the helicity angle.
	// In the reaction Y -> X -> a + b, the helicity angle of 
	// particle a is the angle measured in the rest frame of the
	//decaying parent particle, X, between the direction of the
	// decay daughter a and the direction of the grandparent particle Y.

	TLorentzVector pTmp = (*pa)+(*pb); // this is the mumu system (Z) 4vector
	TVector3 ZboostVector = pTmp.BoostVector(); // this is the 3vector of the Z
	TLorentzVector p; // this is the muon 4vector
	
	if(ca<0)      p.SetPxPyPzE(pa->Px(),pa->Py(),pa->Pz(),pa->E());
	else if(cb<0) p.SetPxPyPzE(pb->Px(),pb->Py(),pb->Pz(),pb->E());
	p.Boost( -ZboostVector ); // boost p to the dimuon CM (rest) frame
	float cosThetaB = p.Vect()*pTmp.Vect()/(p.P()*pTmp.P());
	//if (ySystem(pa,pb) < 0) cosThetaB *= -1.; // reclassify ???
	return cosThetaB;
}

////////////////////////////////////////////////////////////////////////
// Calculates theta and phi in HX frame
////////////////////////////////////////////////////////////////////////
pair<double, double> GetAngles_HX( TLorentzVector a,  TLorentzVector b) {
	TLorentzVector c = a+b;                   // JPsi momentum in lab frame
	TVector3 bv = c.BoostVector();
	TLorentzVector p1(0., 0.,  1380., 1380.); // beam momentum in lab frame
	TLorentzVector p2(0., 0., -1380., 1380.); // beam momentum in lab frame
	p1.Boost(-bv);
	p2.Boost(-bv);
	TVector3 beam1 = p1.Vect().Unit();        // beam direction in JPsi rest frame
	TVector3 beam2 = p2.Vect().Unit();        // beam direction in JPsi rest frame

	TVector3 Z = c.Vect().Unit();             // JPsi direction in lab frame
	TVector3 Y = beam1.Cross( beam2 ).Unit(); // the production plane normal
	TVector3 X = Y.Cross(Z).Unit();           // completes the right-handed coordinate

	a.Boost(-bv);                             // muon+ momentum in JPsi rest frame
	TVector3 mu(a.Vect().Dot(X), a.Vect().Dot(Y), a.Vect().Dot(Z)); // transform to new coordinate
	pair<double, double> angles;
	angles.first = mu.Theta();
	angles.second = mu.Phi()>0. ? mu.Phi() : mu.Phi()+2.*TMath::Pi();
	return angles;
}

//.........这里部分代码省略.........
    double costhphidistr_max = 1. + TMath::Abs(lambda_phi) + TMath::Abs(lambda_thetaphi);
    double costhphidistr_rnd;
    double costhphidistr;
    double costh_gen;
    double sinth_gen;
    double phi_gen;

    if ( lambda_theta > 0. ) costhphidistr_max += lambda_theta;

    do { costh_gen = -1. + 2. * gRandom->Uniform(1.);
         phi_gen   = 2. * gPI * gRandom->Uniform(1.);
         sinth_gen = sqrt( 1. - costh_gen*costh_gen );
         costhphidistr_rnd = costhphidistr_max * gRandom->Uniform(1.);
         costhphidistr = 1. + lambda_theta    * costh_gen*costh_gen
                            + lambda_phi      * sinth_gen*sinth_gen * cos(2.*phi_gen)
                            + lambda_thetaphi * 2.* sinth_gen*costh_gen * cos(phi_gen);
       } while ( costhphidistr_rnd > costhphidistr );


    // lepton momentum in the dilepton rest frame:

    double p_lepton_DILEP = sqrt( 0.25*mass*mass - Mlepton*Mlepton );

    TLorentzVector lepton_DILEP;

    lepton_DILEP.SetXYZM( p_lepton_DILEP * sinth_gen * cos(phi_gen),
                          p_lepton_DILEP * sinth_gen * sin(phi_gen),
                          p_lepton_DILEP * costh_gen,
                          Mlepton );


    // reference directions to calculate angles:

    TVector3 lab_to_dilep = -dilepton.BoostVector();

    TLorentzVector beam1_DILEP = beam1_LAB;
    beam1_DILEP.Boost(lab_to_dilep);         // beam1 in the dilepton rest frame
    TLorentzVector beam2_DILEP = beam2_LAB;
    beam2_DILEP.Boost(lab_to_dilep);         // beam2 in the dilepton rest frame

    TVector3 beam1_direction     = beam1_DILEP.Vect().Unit();
    TVector3 beam2_direction     = beam2_DILEP.Vect().Unit();
    TVector3 dilep_direction     = dilepton.Vect().Unit();
    TVector3 beam1_beam2_bisect  = ( beam1_direction - beam2_direction ).Unit();


    deltaHXCS = dilep_direction.Angle(beam1_beam2_bisect) * 180./gPI;

    // all polarization frames have the same Y axis = the normal to the plane formed by
    // the directions of the colliding hadrons

    TVector3 Yaxis = ( beam1_direction.Cross( beam2_direction ) ).Unit();

    // flip of y axis with rapidity

    if ( rap < 0 ) Yaxis = - Yaxis;

    TVector3 perpendicular_to_beam = ( beam1_beam2_bisect.Cross( Yaxis ) ).Unit();


    // step 1: transform (rotation) lepton momentum components from generation frame
    // to the frame with x,y,z axes as in the laboratory

    TVector3 oldZaxis = beam1_beam2_bisect;
    if ( HX_is_natural ) oldZaxis = dilep_direction;
    if ( PX_is_natural ) oldZaxis = perpendicular_to_beam;

//.........这里部分代码省略.........
			h_lep_min2_E	-> Fill(lep_min2.E());
			h_lep_min2_Pt	-> Fill(ptlepm2);
			h_lep_min2_eta	-> Fill(lep_min2.Eta());
			h_lep_min2_phi	-> Fill(lep_min2.Phi());	

		//reconstructing the two lepton pairs Lorentz vectors
		lpair1 = lep_plus1 + lep_min1;
		lpair2 = lep_plus2 + lep_min2;

		//constructing 3-vectors in the lab frame
		lep_plus1_lab 	= lep_plus1.Vect();
		lep_plus2_lab 	= lep_plus2.Vect();	//.Vect() gives 3 vector from 4vector
		lep_min1_lab 	= lep_min1.Vect();	
		lep_min2_lab 	= lep_min2.Vect();
		lpair1_lab 	= lep_plus1_lab.Cross(lep_min1_lab);	
		lpair2_lab 	= lep_plus2_lab.Cross(lep_min2_lab);

// 		cout << " pt of lepton pair1 on rest frame is: "<< lpair1.Perp()<<endl;


	   	//Filling up Histograms with angles defined in the lab frame
		h_angle_lab_pair1 -> Fill(lep_plus1_lab.Angle(lep_min1_lab));
		h_angle_lab_pair2 -> Fill(lep_plus2_lab.Angle(lep_min2_lab));


       		//Filling up histograms with variables from articles
       		h_angle_lab_deleta1	-> Fill(fabs(lep_min1.Eta() - lep_plus1.Eta()));
       		h_angle_lab_delphi1	-> Fill(fabs(lep_min1.DeltaPhi(lep_plus1)));
       		h_angle_lab_deleta2	-> Fill(fabs(lep_min2.Eta() - lep_plus2.Eta()));
       		h_angle_lab_delphi2	-> Fill(fabs(lep_min2.DeltaPhi(lep_plus2)));


	   	//Looking at the Higgs rest frame
	   	TVector3 boost_rH 	= -rec_H.BoostVector(); //NOTE the minus sign! WHY - sign???
	   	TVector3 boost_rZ1	= -Z1.BoostVector();
	   	TVector3 boost_rZ2	= -Z2.BoostVector();	
	   	Higgs_rest	= rec_H;
	   	Z1_rH		= Z1;
	   	Z2_rH		= Z2;
	   	lep_p1_rH	= lep_plus1;	//
	   	lep_m1_rH	= lep_min1;
	   	lep_p2_rH	= lep_plus2;
	   	lep_m2_rH	= lep_min2;
	   	lep_p1_rZ1	= lep_plus1;
	   	lep_m2_rZ2	= lep_min2;
	   	lep_p2_rZ2	= lep_plus2;
	   	lep_m1_rZ1	= lep_min1;

	   	//Boosting vectors to the Higgs rest frame
	   	Higgs_rest.Boost(boost_rH);
	   	Z1_rH.Boost(boost_rH);
	   	Z2_rH.Boost(boost_rH);
	   	lep_p1_rH.Boost(boost_rH);
	   	lep_m1_rH.Boost(boost_rH);
	   	lep_p2_rH.Boost(boost_rH);
	   	lep_m2_rH.Boost(boost_rH);

	   	//Boosting leptons to Z rest frames
	   	lep_p1_rZ1.Boost(boost_rZ1);
	   	lep_m1_rZ1.Boost(boost_rZ1);
	   	lep_p2_rZ2.Boost(boost_rZ2);
	   	lep_m2_rZ2.Boost(boost_rZ2);

	   	//Setting 3Vectors in Higgs rest frame
	   	Z3_1_rH		= Z1_rH.Vect();
	   	Z3_2_rH		= Z2_rH.Vect();