pair pt constrained QED matrix integrator

qed_cal_ptw.C

double cal_ptw(const double *x){
	TLorentzVector p1(0,0,0,0); // electron daughter 
	p1.SetPtEtaPhiM(x[0],x[1],x[2],mlepton);

	TLorentzVector p2(0,0,0,0); // positron daughter
	double pt = sqrt((par[1]*TMath::Cos(x[3])-p1.X())*(par[1]*TMath::Cos(x[3])-p1.X())+(par[1]*TMath::Sin(x[3])-p1.Y())*(par[1]*TMath::Sin(x[3])-p1.Y()));
	// express the positron in terms of the pair and electron.
	p2.SetXYZM(par[1]*TMath::Cos(x[3])-p1.X(),par[1]*TMath::Sin(x[3])-p1.Y(),pt*TMath::SinH(x[4]),mlepton);


	double q10=0.5*(p1.E()+p2.E()+beta*(p1.Pz()+p2.Pz()));

	// photon one
	TLorentzVector q1(x[5]*TMath::Cos(x[6]),x[5]*TMath::Sin(x[6]),q10/beta,q10);

	// propogator
	TLorentzVector q2(q1.X()+par[0],q1.Y(),q10/beta,q10);

	TLorentzVector pair = p1+p2;

	if(fabs(pair.Rapidity())>Ycut) return 0;
	if(pair.M()>Mmax||pair.M()<Mmin) return 0;
	
	double temp=cal_M2(q1,q2,p1,p2)*x[0]*par[1]*x[5]*p1.P()*p2.P();
	// notice here the post-factors: pT(e-) * pair_pT? * photon kt * P(e-) * P(e+)

	return temp;
}