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

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

//.........这里部分代码省略.........
	p0 = chi.Vect();

	/*1: extract the recoil total energy from the cross-section*/
	ii = 0;
	ii = (int) (E0 / (Ebeam / nFunctionsElastic));
	if (ii >= nFunctionsElastic) ii = nFunctionsElastic - 1; //should not happen!!!

	if (procID == Proc_Pelastic) {
		Tr_max = (2 * Mn * (E0 * E0 - Mchi * Mchi)) / (2 * E0 * Mn + Mn * Mn + Mchi * Mchi); //maximum energy transfer
		if (Tr_max < (Pthr + Pbinding)) return 0; //this event is not compatible with the threshold, it is useless to proceed further
	} else if (procID == Proc_Eelastic) {
		Tr_max = (2 * Me * (E0 * E0 - Mchi * Mchi)) / (2 * E0 * Me + Me * Me + Mchi * Mchi); //maximum energy transfer
		if (Tr_max < Ethr) {
			//cout<<"THR IS: "<<Ethr<<" "<<Tr_max<<" "<<E0<<endl;
			return 0; //this event is not compatible with the threshold, it is useless to proceed further
		}
	} else if (procID == Proc_Nuclelastic) {
		Tr_max = (2 * Mnucl * (E0 * E0 - Mchi * Mchi)) / (2 * E0 * Mnucl + Mnucl * Mnucl + Mchi * Mchi); //maximum energy transfer
		if (Tr_max < Nuclthr) {
			//cout<<"THR IS: "<<Ethr<<" "<<Tr_max<<" "<<E0<<endl;
			return 0; //this event is not compatible with the threshold, it is useless to proceed further
		}
	}
	if (procID == Proc_Pelastic) {
		Er = f_chipXsection[ii]->GetRandom(Pthr + Pbinding + Mn, Tr_max + Mn);
	} else if (procID == Proc_Eelastic) {
		Er = f_chieXsection[ii]->GetRandom(Ethr + Me, Tr_max + Me);
	} else if (procID == Proc_Nuclelastic) {
		Er = f_chinuclXsection[ii]->GetRandom(Nuclthr, Tr_max); //Here the variable is the KINETIC energy of the recoiling nucleus
		Er = Er + Mnucl;
	}
	/*1a: correct the proton energy for binding effects*/
	if (procID == Proc_Pelastic) {
		Er = Er - Pbinding; /*Effective binding energy correction*/
	}

	/*1b: compute x-section total . No time consuming, since integrals are cached!*/
	if (procID == Proc_Pelastic) {
		sigma = f_chipXsection[ii]->Integral(Pthr + Pbinding + Mn, Tr_max + Mn);
	} else if (procID == Proc_Eelastic) {
		sigma = f_chieXsection[ii]->Integral(Ethr + Me, Tr_max + Me);
	} else if (procID == Proc_Nuclelastic) {
		sigma = f_chinuclXsection[ii]->Integral(Nuclthr, Tr_max); //Has to be integrated in this range since the variable is the KINETIC energy here (and not the total as before)
	}

	/*2: compute recoil chi TOTAL energy*/
	if (procID == Proc_Pelastic) {
		Echi = E0 + Mn - Er;
	} else if (procID == Proc_Eelastic) {
		Echi = E0 + Me - Er;
	} else if (procID == Proc_Nuclelastic) {
		Echi = E0 + Mnucl - Er;
	}

	/*3: compute the momenta*/
	Pchi = sqrt(Echi * Echi - Mchi * Mchi);
	if (procID == Proc_Pelastic) {
		Pr = sqrt(Er * Er - Mn * Mn);
	} else if (procID == Proc_Eelastic) {
		Pr = sqrt(Er * Er - Me * Me);
	} else if (procID == Proc_Nuclelastic) {
		Pr = sqrt(Er * Er - Mnucl * Mnucl);
	}
	/*4: compute the angle of the recoil nucleon wrt the initial chi momentum direction*/
	if (procID == Proc_Pelastic) {
		ctheta_r = E0 * E0 - Echi * Echi + Er * Er - Mn * Mn;
		ctheta_r /= 2 * P0 * Pr;
	} else if (procID == Proc_Eelastic) {
		ctheta_r = E0 * E0 - Echi * Echi + Er * Er - Me * Me;
		ctheta_r /= 2 * P0 * Pr;
	} else if (procID == Proc_Nuclelastic) {
		ctheta_r = E0 * E0 - Echi * Echi + Er * Er - Mnucl * Mnucl;
		ctheta_r /= 2 * P0 * Pr;
	}
	if (ctheta_r > 1) ctheta_r = 1;
	if (ctheta_r < -1) ctheta_r = -1;
	stheta_r = sqrt(1 - ctheta_r * ctheta_r);
	/*5: The azimuthal angle (around the incoming chi momentum direction) is flat*/
	phi_r = Rand.Uniform(-PI, PI);

	/*6: Now set the 4-vectors*/
	/*6a: build an orthogonal coordinate system, with v0 along the initial chi momentum direction*/
	v0 = chi.Vect().Unit();
	v1 = v0.Orthogonal();
	v1 = v1.Unit();
	v2 = v0.Cross(v1); //v2 = v0 x v1

	/*write the 3-momenta*/
	pr = v0 * Pr * ctheta_r + v1 * Pr * stheta_r * sin(phi_r) + v2 * Pr * stheta_r * cos(phi_r);
	pchi = p0 - pr;

	/*6b:Set them */
	recoil.SetVect(pr);
	recoil.SetE(Er);
	recoil_chi.SetVect(pchi);
	recoil_chi.SetE(Echi);

	return sigma;

}