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

float computeR2(TLorentzVector hem1, TLorentzVector hem2, TLorentzVector pfMet) {
    double mR = computeMR(hem1, hem2);
    double term1 = pfMet.Pt()/2*(hem1.Pt() + hem2.Pt());
    double term2 = pfMet.Px()/2*(hem1.Px() + hem2.Px()) + pfMet.Py()/2*(hem1.Py() + hem2.Py()); //dot product of MET with (p1T + p2T)
    double mTR = sqrt(term1 - term2);
    return (mTR / mR) * (mTR / mR);
}

//___________________[ TAUS RECONSTRUCTION AND INVMASS CALCULATION ]_________________________________________
double SUSYLooperHistsSoftBase::DiTau_InvMass( TLorentzVector Met, TLorentzVector L1, TLorentzVector L2, float Al ) {

    TLorentzVector T1,T2;    double DiTauMass;    TMatrixF A(2,2);    TVectorF C(2),X(2);
    A(0,0)=L1.Px();
    A(0,1)=L2.Px();
    A(1,0)=L1.Py();
    A(1,1)=L2.Py();
    A=A.Invert();
    C(0)=(Met+L1+L2).Px();
    C(1)=(Met+L1+L2).Py();
    X=A*C;
    //double X0i=X(0), X1i=X(1);
    if ( (X(0)<0.||X(1)<0.) && Al>0. ) {//---[MET Alignement subsection]--------------------------------------------------------------------------------------------
        if      ( fabs(L1.DeltaPhi(Met))>Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) {}//{DO NOTHING just normaly a non-Z event!}
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
        else if ( fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))<fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
    }//-------------------------------------------------------------------------------------------------------------------------------------------------------------
    C(0)=(Met+L1+L2).Px();
    C(1)=(Met+L1+L2).Py();
    X=A*C;
    T1.SetPxPyPzE( L1.Px()*X(0), L1.Py()*X(0), L1.Pz()*X(0), sqrt( 3.1571 +L1.P()*L1.P()*X(0)*X(0) ) );
    T2.SetPxPyPzE( L2.Px()*X(1), L2.Py()*X(1), L2.Pz()*X(1), sqrt( 3.1571 +L2.P()*L2.P()*X(1)*X(1) ) );
    if (  X(0)>0.  &&  X(1)>0.   )   DiTauMass=(T1+T2).M();
    //if (  (X(0)>0. && 0.<X(1)&&X(1)<1.)  ||  (X(1)>0. && 0.<X(0)&&X(0)<1.)  )  DiTauMass=(T1+T2).M(); // B
    //if (  (X(0)>1.&& X(1)<0.) || (X(1)>1.&& X(0)<0.)    )  DiTauMass=(T1+T2).M(); // C
    //if (  X(0)<0.  &&  X(1)<0.   )  DiTauMass=(T1+T2).M(); //  D    
    else DiTauMass=-(T1+T2).M();
    return DiTauMass;
}//-----------------------------------------------------------------------------------------------------------

double SUSYLooperHists::DiTau_InvMass( TLorentzVector Met, TLorentzVector L1, TLorentzVector L2, float Al ){   
  TLorentzVector T1,T2;
  double DiTauMass; 
  TMatrixF A(2,2); 
  TVectorF C(2),X(2); 
  A(0,0)=L1.Px();
  A(0,1)=L2.Px();
  A(1,0)=L1.Py();
  A(1,1)=L2.Py();
  A=A.Invert();
  C(0)=(Met+L1+L2).Px();
  C(1)=(Met+L1+L2).Py();
  X=A*C;// double X0i=X(0), X1i=X(1);
  //---------------[ MET ReAlignement subsection ]------------------------------
  if(X(0)<0||X(1)<0){ 
    if     ( fabs(L1.DeltaPhi(Met))>Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) {/*DO NOTHING just normaly a non-Z event!*/}
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
    else if( fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))<fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
  }//---------------------------------------------------------------------------
  C(0)=(Met+L1+L2).Px(); C(1)=(Met+L1+L2).Py(); X=A*C;
  T1.SetPxPyPzE( L1.Px()*X(0), L1.Py()*X(0), L1.Pz()*X(0), sqrt( 3.1571 +L1.P()*L1.P()*X(0)*X(0) ) );
  T2.SetPxPyPzE( L2.Px()*X(1), L2.Py()*X(1), L2.Pz()*X(1), sqrt( 3.1571 +L2.P()*L2.P()*X(1)*X(1) ) );
  if( X(0)>0 && X(1)>0 ) DiTauMass=(T1+T2).M();  else DiTauMass=-(T1+T2).M();  return DiTauMass;
  //if((X(0)!=X0i||X(1)!=X1i))std::cout<<X(0)<<" "<<X(1)<<" <--"<<X0i<<" "<<X1i<<" RMETal.phi="<<(T1+T2-L1-L2).Phi()<<" RMETal.eta"<<(T1+T2-L1-L2).Eta()<<" MZ="<<DiTauMass<<endl; 
}

/*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 rochcor2012::musclefit_data( TLorentzVector& mu, TLorentzVector& mubar){

  float dpar1 = 0.0;
  float dpar2 = 0.0;
  float epar1 = 0.0;
  float epar2 = 0.0;
  
  if(fabs(mu.PseudoRapidity())<=0.9){
    dpar1 = d0par;
    epar1 = e0par;
  }else if(mu.PseudoRapidity()>0.9){
    dpar1 = d1par;
    epar1 = e1par;
  }else if(mu.PseudoRapidity()<-0.9){
    dpar1 = d2par;
    epar1 = e2par;
  }

  if(fabs(mubar.PseudoRapidity())<=0.9){
    dpar2 = d0par;
    epar2 = e0par;
  }else if(mubar.PseudoRapidity()>0.9){
    dpar2 = d1par;
    epar2 = e1par;
  }else if(mubar.PseudoRapidity()<-0.9){
    dpar2 = d2par;
    epar2 = e2par;
  }

  float corr1 = 1.0 + bpar*mu.Pt() + (-1.0)*cpar*mu.Pt()*TMath::Sign(float(1.0),float(mu.PseudoRapidity()))*TMath::Power(mu.PseudoRapidity(),2)
    + (-1.0)*dpar1*mu.Pt()*sin(mu.Phi() + epar1);
  float corr2 = 1.0 + bpar*mubar.Pt() + (1.0)*cpar*mubar.Pt()*TMath::Sign(float(1.0),float(mubar.PseudoRapidity()))*TMath::Power(mubar.PseudoRapidity(),2)
    + (1.0)*dpar2*mubar.Pt()*sin(mubar.Phi() + epar2);
  
  float px1 = mu.Px();
  float py1 = mu.Py();
  float pz1 = mu.Pz();
  float e1 = mu.E();
  
  float px2 = mubar.Px();
  float py2 = mubar.Py();
  float pz2 = mubar.Pz();
  float e2 = mubar.E();

  px1 *= corr1;
  py1 *= corr1;
  pz1 *= corr1;
  e1 *= corr1;
  
  px2 *= corr2;
  py2 *= corr2;
  pz2 *= corr2;
  e2 *= corr2;
  
  mu.SetPxPyPzE(px1,py1,pz1,e1);
  mubar.SetPxPyPzE(px2,py2,pz2,e2);

}

//these can be done better, pass in a specific struct and populate that struct
//could reduce these 6 calls into 1 or 2
void DalitzChiSq2::setap2(TLorentzVector vec){
	ap2.v.SetXYZM(vec.Px(),vec.Py(),vec.Pz(),vec.M());

	double Pt = TMath::Sqrt(vec.Px()*vec.Px() + vec.Py()*vec.Py());
	double p = TMath::Sqrt(Pt*Pt + vec.Pz()*vec.Pz());
	double thet = TMath::ACos(vec.Pz()/p);
	
	ap2.pt=Pt;
	ap2.theta=thet;
	ap2.x_m=1/Pt;
	
}

HHKinFit2::HHKinFitMasterSingleHiggs::HHKinFitMasterSingleHiggs(TLorentzVector const& tauvis1,
                                                                TLorentzVector const& tauvis2,
                                                                TVector2 const& met, 
                                                                TMatrixD const& met_cov, 
                                                                bool istruth,
                                                                TLorentzVector const& higgsgen)
:m_MET_COV(TMatrixD(4,4))
{
  
  m_tauvis1 = HHLorentzVector(tauvis1.Px(), tauvis1.Py(), tauvis1.Pz(), tauvis1.E());  
  m_tauvis2 = HHLorentzVector(tauvis2.Px(), tauvis2.Py(), tauvis2.Pz(), tauvis2.E());
 
  m_tauvis1.SetMkeepE(1.77682);
  m_tauvis2.SetMkeepE(1.77682);
   
  m_MET = met;
  m_MET_COV = met_cov;

  m_chi2_best = pow(10,10);
  m_bestHypo = 0;

  //  if (istruth){
  //    TRandom3 r(0);
  //
  //    HHLorentzVector recoil;
  //    if(heavyhiggsgen != NULL){
  //       Double_t pxRecoil = r.Gaus(-(heavyhiggsgen->Px() ), 10.0);
  //       Double_t pyRecoil = r.Gaus(-(heavyhiggsgen->Py() ), 10.0);
  //
  //       recoil = HHLorentzVector(pxRecoil, pyRecoil, 0,
  //				sqrt(pxRecoil*pxRecoil+pyRecoil*pyRecoil));
  //    }
  //    else{
  //      recoil = HHLorentzVector(0,0,0,0);
  //      std::cout << "WARNING! Truthinput mode active but no Heavy Higgs gen-information given! Setting Recoil to Zero!" << std::endl;
  //    }
  //
  //    TMatrixD recoilCov(2,2);
  //    recoilCov(0,0)=100;  recoilCov(0,1)=0;
  //    recoilCov(1,0)=0;    recoilCov(1,1)=100;
  //
  //    HHLorentzVector recoHH = m_bjet1 + m_bjet2 + m_tauvis1 + m_tauvis2 + recoil;
  //    m_MET = TVector2(-recoHH.Px(), -recoHH.Py() );
  //
  //    m_MET_COV = TMatrixD(2,2);
  //    m_MET_COV = recoilCov + bjet1Cov + bjet2Cov;
  //
  //  }
}

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;

}

double  METMHTAsys(MissingET* met,vector<Jet> jetvec,vector<Muon> muonvec,vector<Electron> electronvec,vector<Photon> photonvec){
  double Met=-99;
  double METAsys=-99;
  TVector2 PUCorMet, RawMet;
   TLorentzVector allvecsum;
  allvecsum.SetPxPyPzE(0, 0, 0, 0);
  PUCorMet.Set(0., 0.);
  RawMet.Set(0.0, 0.0);
  
  for(int i=0; i<(int)jetvec.size(); i++) {allvecsum += jetvec.at(i).P4();}
  for(int j=0; j<(int)muonvec.size(); j++) {allvecsum += muonvec.at(j).P4();}
  for(int k=0; k<(int)electronvec.size(); k++) {allvecsum += electronvec.at(k).P4();}
  for(int l=0; l<(int)photonvec.size(); l++) {allvecsum += photonvec.at(l).P4();}

  PUCorMet.Set(-allvecsum.Px(),-allvecsum.Py());
  Met= PUCorMet.Mod();
  RawMet.SetMagPhi(met->MET, met->Phi);
  
  METAsys=fabs(Met-(RawMet.Mod()))/(Met+(RawMet.Mod()));//this is funny. RawMet.Mod() must return met->MET. We didn't need to build RawMet to obtain its magnitude :):0
  //cout << "......................RawMet.Mod(): " << RawMet.Mod() << endl;
  //cout << "...................... Met: " << Met << endl; 
  //cout << "...................... METAsys: " << METAsys << endl;
  return METAsys;

}

double computeMETPhi(const double met, const double metphi, const TLorentzVector &photon)
{
  double metx = met*cos(metphi) + photon.Px();
  double mety = met*sin(metphi) + photon.Py();
  TLorentzVector mvec(metx, mety, 0, 0);
  return mvec.Phi();
}

float getMT(TLorentzVector pZ, TLorentzVector pH)
{
  //take MET from pZ
  float myMET = pZ.Pt();
  float myMETx = pZ.Px();
  float myMETy = pZ.Py();

  float Et = pH.Et();
  float px = pH.Px();
  float py = pH.Py();
  float m  = pH.M();

  float MT = -99;
  float MT2 = m*m + 2*( Et * myMET - (px*myMETx + py*myMETy) ); 
  if(MT2>0.) {MT=sqrt(MT2);}
  
  return MT;
}

void rochcor2012::momcor_mc( TLorentzVector& mu, float charge, float sysdev, int runopt,bool sync=false){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  //float mptsys = sran.Gaus(0.0,sysdev);
  
  float dm = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/mmavg[mu_phibin][mu_etabin];
  float da = mcor_ma[mu_phibin][mu_etabin] + mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin];
  
  float cor = 1.0/(1.0 + dm + charge*da*ptmu);

  //for the momentum tuning - eta,phi,Q correction
  px *= cor;
  py *= cor;
  pz *= cor;
  e *= cor;
  
  float gscler = 0.0;
  float deltaer = 0.0;
  float sfer = 0.0;
  
  gscler = TMath::Sqrt( TMath::Power(mgscl_stat,2) + TMath::Power(mgscl_syst,2) );
  deltaer = TMath::Sqrt( TMath::Power(delta_stat,2) + TMath::Power(delta_syst,2) );
  sfer = TMath::Sqrt( TMath::Power(sf_stat,2) + TMath::Power(sf_syst,2) );

  float tune;
  if (sync) tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*(1.0 + (sf + sysdev*sfer)));
  else tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*eran.Gaus(1.0,(sf + sysdev*sfer)));
  
  px *= (tune);
  py *= (tune);
  pz *= (tune);
  e *= (tune);
  
  float gscl = (genm_smr/mrecm);
  
  px *= (gscl + gscler_mc_dev*gscler);
  py *= (gscl + gscler_mc_dev*gscler);
  pz *= (gscl + gscler_mc_dev*gscler);
  e *= (gscl + gscler_mc_dev*gscler);

  mu.SetPxPyPzE(px,py,pz,e);
  
}

void rochcor2012::momcor_mc( TLorentzVector& mu, float charge, float sysdev, int runopt, float& qter){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  //float mptsys = sran.Gaus(0.0,sysdev);
  
  float Mf = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/(mpavg[mu_phibin][mu_etabin]+mmavg[mu_phibin][mu_etabin]);
  float Af = ((mcor_ma[mu_phibin][mu_etabin]+mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin]) - Mf*(mpavg[mu_phibin][mu_etabin]-mmavg[mu_phibin][mu_etabin]));     
  
  float cor = 1.0/(1.0 + 2.0*Mf + charge*Af*ptmu);
  
  //for the momentum tuning - eta,phi,Q correction
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  float gscler = mgscl_stat;
  float deltaer = delta_stat;
  float sfer = sf_stat;
    
  float gscl = (genm_smr/mrecm);
  
  px *= (gscl + gscler_mc_dev*gscler);
  py *= (gscl + gscler_mc_dev*gscler);
  pz *= (gscl + gscler_mc_dev*gscler);
  e  *= (gscl + gscler_mc_dev*gscler);
  
  float momscl = sqrt(px*px + py*py)/ptmu;
  
  float tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*eran.Gaus(1.0,(sf + sysdev*sfer)));
  
  px *= (tune); 
  py *= (tune);  
  pz *= (tune);  
  e  *= (tune);   
  
  qter *= (momscl*momscl + (1.0-tune)*(1.0-tune));
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

inline float kinematics::betaiSystem( TLorentzVector* pa, TLorentzVector* pb, int i)
{
	TLorentzVector pTmp = (*pa)+(*pb);
	float E = pTmp.E();
	if(E<=0.) return -99999999.;

	if     (i==1) return pTmp.Px()/E;
	else if(i==2) return pTmp.Py()/E;
	else if(i==3) return pTmp.Pz()/E;
	else _WARNING("i needs to be 1,2,3 (x,y,z), returning -99999999.");
	return -99999999.;
}

///------------------------------------------------------------------------
/// MEKD::computeKD - compute KD and MEs for the input processes A and B
///------------------------------------------------------------------------
int MEKD::computeKD( TString processA, TString processB,
                     TLorentzVector lept1P, int lept1Id, TLorentzVector lept2P, int lept2Id,
                     TLorentzVector lept3P, int lept3Id, TLorentzVector lept4P, int lept4Id,
                     double& kd, double& me2processA, double& me2processB )
{
	/// Prepare 4-momenta in the required format
	lept1P_i[0] = lept1P.E();
	lept1P_i[1] = lept1P.Px();
	lept1P_i[2] = lept1P.Py();
	lept1P_i[3] = lept1P.Pz();
	
	lept2P_i[0] = lept2P.E();
	lept2P_i[1] = lept2P.Px();
	lept2P_i[2] = lept2P.Py();
	lept2P_i[3] = lept2P.Pz();
	
	lept3P_i[0] = lept3P.E();
	lept3P_i[1] = lept3P.Px();
	lept3P_i[2] = lept3P.Py();
	lept3P_i[3] = lept3P.Pz();
	
	lept4P_i[0] = lept4P.E();
	lept4P_i[1] = lept4P.Px();
	lept4P_i[2] = lept4P.Py();
	lept4P_i[3] = lept4P.Pz();
	
	/// Load internal containers
	four_particle_Ps_i[0] = lept1P_i;
	four_particle_Ps_i[1] = lept2P_i;
	four_particle_Ps_i[2] = lept3P_i;
	four_particle_Ps_i[3] = lept4P_i;
	
	four_particle_IDs_i[0] = lept1Id;
	four_particle_IDs_i[1] = lept2Id;
	four_particle_IDs_i[2] = lept3Id;
	four_particle_IDs_i[3] = lept4Id;
	
	return computeKD( (string) processA.Data(), (string) processB.Data(), four_particle_Ps_i, four_particle_IDs_i, kd, me2processA, me2processB );
}