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

void getKinematics(TLorentzVector& b,double const mass)
{
   float const pt = gRandom->Uniform(momentumRange.first,momentumRange.second);
   float const eta = gRandom->Uniform(-acceptanceEta,acceptanceEta);
   float const phi = TMath::TwoPi() * gRandom->Rndm();

   b.SetXYZM(pt * cos(phi),pt * sin(phi),pt * sinh(eta), mass);
}

TLorentzVector leptonic_fitter_algebraic::lv_with_mass( const TMatrixD& mat, double mass )
{
  assert( mat.GetNcols() == 1 || mat.GetNrows() == 1 );
  TLorentzVector lv;
  const double *array = mat.GetMatrixArray();
  lv.SetXYZM( array[0], array[1], array[2], mass );
  return lv;
}

TLorentzVector smearMom(TLorentzVector const& b,TF1 const * const fMomResolution)
{
  float const pt = b.Perp();
  float const sPt = gRandom->Gaus(pt,pt*fMomResolution->Eval(pt));

  TLorentzVector sMom;
  sMom.SetXYZM(sPt*cos(b.Phi()),sPt*sin(b.Phi()),sPt*sinh(b.PseudoRapidity()),b.M());
  return sMom;
}

    // Reset all private to false/zwwDefault/null vector/etc.
void ZWW::reset(){
        // For checker
    isAllOk_ = false;
    isLepOk_ = false;
    isJetOk_ = false;
    isMETOk_ = false;

        // For setter
    lepPt_. clear();
    lepEta_.clear();
    lepPhi_.clear();
    lepFl_. clear();
    lepCh_. clear();
    lepId_. clear();
    jetPt_. clear();
    jetEta_.clear();
    jetPhi_.clear();
    jetM_.  clear();
    jetTag_.clear();
    met_          = zwwDefault;
    metPhi_       = zwwDefault;
    njet_         = zwwDefault;
    nbjet_        = zwwDefault;
    st_           = zwwDefault;
    minMt_        = zwwDefault;
    minDeltaPhi_  = zwwDefault;

        // Buffer
    lepVec_.clear();
    jetVec_.clear();
//    TLorentzVector metVec_.SetXYZM(0.,0.,0.,0.);
    metVec_.SetXYZM(0.,0.,0.,0.);
    z0LepIdx_[0] = zwwDefault;
    z0LepIdx_[1] = zwwDefault;
    z1LepIdx_[0] = zwwDefault;
    z1LepIdx_[1] = zwwDefault;
    zaLepIdx_[0] = zwwDefault;
    zaLepIdx_[1] = zwwDefault;
    zbLepIdx_[0] = zwwDefault;
    zbLepIdx_[1] = zwwDefault;
    flagZ1SF_    = zwwDefault;

        // All other variables
    outPreSelection_ = false;

}

//.........这里部分代码省略.........

      otree->byCombinedIsolationDeltaBetaCorrRaw3Hits_1 = 0;
      otree->againstElectronLooseMVA5_1 = 0;
      otree->againstElectronMediumMVA5_1 = 0;
      otree->againstElectronTightMVA5_1 = 0;
      otree->againstElectronVLooseMVA5_1 = 0;
      otree->againstElectronVTightMVA5_1 = 0;
      otree->againstMuonLoose3_1 = 0;
      otree->againstMuonTight3_1 = 0;
      
      // filling tau variables
      otree->pt_2 = analysisTree.tau_pt[tauIndex];
      otree->eta_2 = analysisTree.tau_eta[tauIndex];
      otree->phi_2 = analysisTree.tau_phi[tauIndex];
      otree->q_2 = -1;
      if (analysisTree.tau_charge[tauIndex]>0)
	otree->q_2 = 1;
      otree->mva_2 = log(0);
      otree->d0_2 = analysisTree.tau_dxy[tauIndex];
      otree->dZ_2 = analysisTree.tau_dz[tauIndex];
      otree->iso_2 = analysisTree.tau_byCombinedIsolationDeltaBetaCorrRaw3Hits[tauIndex];
      otree->m_2 = analysisTree.tau_mass[tauIndex];

      otree->byCombinedIsolationDeltaBetaCorrRaw3Hits_2 = analysisTree.tau_byCombinedIsolationDeltaBetaCorrRaw3Hits[tauIndex];
      otree->againstElectronLooseMVA5_2 = analysisTree.tau_againstElectronLooseMVA5[tauIndex];
      otree->againstElectronMediumMVA5_2 = analysisTree.tau_againstElectronMediumMVA5[tauIndex];
      otree->againstElectronTightMVA5_2 = analysisTree.tau_againstElectronTightMVA5[tauIndex];
      otree->againstElectronVLooseMVA5_2 = analysisTree.tau_againstElectronVLooseMVA5[tauIndex];
      otree->againstElectronVTightMVA5_2 = analysisTree.tau_againstElectronVTightMVA5[tauIndex];
      otree->againstMuonLoose3_2 = analysisTree.tau_againstMuonLoose3[tauIndex];
      otree->againstMuonTight3_2 = analysisTree.tau_againstMuonTight3[tauIndex];

      // ditau system
      TLorentzVector electronLV; electronLV.SetXYZM(analysisTree.electron_px[electronIndex],
					    analysisTree.electron_py[electronIndex],
					    analysisTree.electron_pz[electronIndex],
					    electronMass);

      TLorentzVector tauLV; tauLV.SetXYZM(analysisTree.tau_px[tauIndex],
					  analysisTree.tau_py[tauIndex],
					  analysisTree.tau_pz[tauIndex],
					  tauMass);

      TLorentzVector dileptonLV = electronLV + tauLV;

      otree->m_vis = dileptonLV.M();
      otree->pt_tt = dileptonLV.Pt();

      // opposite charge
      otree->os = (otree->q_1 * otree->q_2) < 0.;

      // dimuon veto
      otree->dilepton_veto = 0;

      for (unsigned int ie = 0; ie<analysisTree.electron_count; ++ie) {
	if (analysisTree.electron_pt[ie]<ptDiElectronVeto) continue;
	if (fabs(analysisTree.electron_eta[ie])>etaDiElectronVeto) continue;	
	
	if (fabs(analysisTree.electron_dxy[ie])>dxyElectronCut) continue;
	if (fabs(analysisTree.electron_dz[ie])>dzElectronCut) continue;

	if (otree->q_1 * analysisTree.electron_charge[ie] > 0.) continue;

	float neutralHadIsoEle = analysisTree.electron_neutralHadIso[ie];
	float photonIsoEle = analysisTree.electron_photonIso[ie];
	float chargedHadIsoEle = analysisTree.electron_chargedHadIso[ie];

//.........这里部分代码省略.........

	    if (i_event%n_step == 0) {cout << n_step_*20 <<" % "<<endl; n_step_++;}

	    // generation of dilepton in the pp event in the pp CM

    // mass

    isBG = 0;

    if ( gRandom->Uniform() < f_BG ) { mass = gRandom->Uniform(mass_min, mass_max); isBG = 1; }
                                        else { do { mass = gRandom->Gaus(mass_signal_peak, mass_signal_sigma); }
                                               while ( mass < mass_min || mass > mass_max ); }

    // pT:

    pT = pT_distr->GetRandom();

    // pL:

    double rap_sign = gRandom->Uniform(-1., 1.); rap_sign /= TMath::Abs(rap_sign);
    rap = rap_distr->GetRandom() * rap_sign;
    double mT = sqrt( mass*mass + pT*pT );
    double pL1 = 0.5 *mT * exp(rap);
    double pL2 = - 0.5 *mT * exp(-rap);
    double pL = pL1 + pL2;

    // Phi:

    double Phi = 2. * gPI * gRandom->Uniform(1.);

    // 4-vector:

    TLorentzVector dilepton;
    dilepton.SetXYZM( pT * cos(Phi) , pT * sin(Phi), pL, mass );


    // generation of polarization (generic reference frame)

    double lambda_theta    = lambda_theta_sig;
    double lambda_phi      = lambda_phi_sig;
    double lambda_thetaphi = lambda_thetaphi_sig;
    bool HX_is_natural = HX_is_natural_sig;
    bool PX_is_natural = PX_is_natural_sig;

    if ( isBG ) {
       lambda_theta    = lambda_theta_bkg;
       lambda_phi      = lambda_phi_bkg;
       lambda_thetaphi = lambda_thetaphi_bkg;
       HX_is_natural = HX_is_natural_bkg;
       PX_is_natural = PX_is_natural_bkg;
    }


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

void leptonic_fitter_algebraic::update_nu_and_decay_chain( double nu_px, double nu_py, double nu_pz )
{
  TLorentzVector nu;
  nu.SetXYZM( nu_px, nu_py, nu_pz, 0. );
  update_nu_and_decay_chain( nu );
}

//.........这里部分代码省略.........
      iso_1 = isoMuMin;
      m_1 = muonMass;
      float dPhiMETMuon = dPhiFrom2P(analysisTree.muon_px[muonIndex],analysisTree.muon_py[muonIndex],
				     analysisTree.pfmet_ex,analysisTree.pfmet_ey);
      mt_1 = TMath::Sqrt(2*met*analysisTree.muon_pt[muonIndex]*(1-TMath::Cos(dPhiMETMuon)));


      // filling tau variables
      pt_2 = analysisTree.tau_pt[tauIndex];
      eta_2 = analysisTree.tau_eta[tauIndex];
      phi_2 = analysisTree.tau_phi[tauIndex];
      q_2 = -1;
      if (analysisTree.tau_charge[tauIndex]>0)
        q_2 = 1;
      mva_2 = -9999;
      d0_2 = analysisTree.tau_leadchargedhadrcand_dxy[tauIndex];
      dZ_2 = analysisTree.tau_leadchargedhadrcand_dz[tauIndex];
      iso_2 = isoTauMin;
      m_2 = analysisTree.tau_mass[tauIndex];
      float dPhiMETTau = dPhiFrom2P(analysisTree.tau_px[tauIndex],analysisTree.tau_py[tauIndex],
				    analysisTree.pfmet_ex,analysisTree.pfmet_ey);
      mt_2 = TMath::Sqrt(2*met*analysisTree.tau_pt[tauIndex]*(1-TMath::Cos(dPhiMETTau)));


      byCombinedIsolationDeltaBetaCorrRaw3Hits_2 = analysisTree.tau_byCombinedIsolationDeltaBetaCorrRaw3Hits[tauIndex];
      againstElectronLooseMVA5_2 = analysisTree.tau_againstElectronLooseMVA5[tauIndex];
      againstElectronMediumMVA5_2 = analysisTree.tau_againstElectronMediumMVA5[tauIndex];
      againstElectronTightMVA5_2 = analysisTree.tau_againstElectronTightMVA5[tauIndex];
      againstElectronVLooseMVA5_2 = analysisTree.tau_againstElectronVLooseMVA5[tauIndex];
      againstElectronVTightMVA5_2 = analysisTree.tau_againstElectronVTightMVA5[tauIndex];
      againstMuonLoose3_2 = analysisTree.tau_againstMuonLoose3[tauIndex];
      againstMuonTight3_2 = analysisTree.tau_againstMuonTight3[tauIndex];

      TLorentzVector muonLV; muonLV.SetXYZM(analysisTree.muon_px[muonIndex],
					    analysisTree.muon_py[muonIndex],
					    analysisTree.muon_pz[muonIndex],
					    muonMass);

      TLorentzVector tauLV; tauLV.SetXYZT(analysisTree.tau_px[tauIndex],
					  analysisTree.tau_py[tauIndex],
					  analysisTree.tau_pz[tauIndex],
					  analysisTree.tau_e[tauIndex]);

      TLorentzVector metLV; metLV.SetXYZT(analysisTree.pfmet_ex,
					  analysisTree.pfmet_ey,
					  0,
					  TMath::Sqrt(analysisTree.pfmet_ex*analysisTree.pfmet_ex + analysisTree.pfmet_ey*analysisTree.pfmet_ey));
      
      TLorentzVector dileptonLV = muonLV + tauLV;

      // visible mass
      m_vis = dileptonLV.M();
      // visible ditau pt 
      pt_tt = (dileptonLV+metLV).Pt();

      // bisector of electron and muon transverse momenta
      float tauUnitX = tauLV.Px()/tauLV.Pt();
      float tauUnitY = tauLV.Py()/tauLV.Pt();
	
      float muonUnitX = muonLV.Px()/muonLV.Pt();
      float muonUnitY = muonLV.Py()/muonLV.Pt();

      float zetaX = tauUnitX + muonUnitX;
      float zetaY = tauUnitY + muonUnitY;
      
      float normZeta = TMath::Sqrt(zetaX*zetaX+zetaY*zetaY);

//.........这里部分代码省略.........
  YieldTree -> SetBranchAddress ("mass",            &mass);
  YieldTree -> SetBranchAddress ("xF",              &xF);
  //YieldTree -> SetBranchAddress ("xB",              &xB);
  YieldTree -> SetBranchAddress ("xT",              &xT);
  //YieldTree -> SetBranchAddress ("costh",              &costh);
  //YieldTree -> SetBranchAddress ("phi",              &phi);
  //YieldTree -> SetBranchAddress ("trackSeparation", &trackSeparation);
  //YieldTree -> SetBranchAddress ("chisq_dimuon",    &chisq_dimuon);
  //YieldTree -> SetBranchAddress ("px1",             &px1);
  //YieldTree -> SetBranchAddress ("py1",             &py1);
  //YieldTree -> SetBranchAddress ("pz1",             &pz1);
  //YieldTree -> SetBranchAddress ("px2",             &px2);
  //YieldTree -> SetBranchAddress ("py2",             &py2);
  //YieldTree -> SetBranchAddress ("pz2",             &pz2);
  YieldTree -> SetBranchAddress ("target",          &target);
  YieldTree -> SetBranchAddress ("dump",            &dump);
  //YieldTree -> SetBranchAddress ("targetPos",       &targetPos);

  //here's the for loop
  int nentries = YieldTree -> GetEntries();
  int spillID0;
  int daycount = 0;
  int DspillID;
  cout << "Number of entries: " << nentries << endl;
  for (int i = 0; i < nentries; i++) {
    //for (int i = 0; i < 10000; i++) {
    YieldTree -> GetEntry (i);
    if(i%100000 == 0) cout << "entry " << i << endl;
    if(i == 0 ) spillID0 = spillID;
    if(runID > 13800 && runID < 14799)continue;
    if(xT < 0.1 || xT > 0.4)continue;
    //Now do the kinematic cuts!
    xf = xF;
    lordm.SetXYZM(dpx, dpy, dpz, mass);
    pt = lordm.Pt();
    if(pt < 0.5) continue;
    //if(fabs(lordm.Py()/lordm.Px()) < 0.176 || fabs(lordm.Py()/lordm.Px()) > 5.67) continue;
    //if(fabs(lordm.Px()/lordm.Py()) < 0.176) continue;
    phi = lordm.Phi();

    //make a left right phi!
    phi2 = phi - PI/2;
    if(phi2 < -PI){phi2 = phi2 + 2*PI;}
    phi3 = phi - PI/2;
    if(phi3 < -PI){phi3 = phi3 + 2*PI;}

    if(nbins == 1){k = 0;}
    else{
      if(Trig == 0 || Trig == 2){int k = int (10 * xf);}
    }

    if(nbins == 1){
      if(xf < 0.2 || xf > 0.6)continue;}
    else if(nbins > 1){
      if(xf < 0.2 || xf > 0.7)continue;
    }
    if(nbins == 1){k = 0;}
    else{
      int k = int (10 * (xf - 0.2));
    }

    
    if (k < 0) continue;
    if (k > nbins-1) k = nbins - 1;
    //if (k > nbins-1) continue;