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

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

void FillK2pi(superCmpEvent* sevt, Hist_dir* dir, TLorentzVector pion,TLorentzVector el1,TLorentzVector el2,TLorentzVector gamma){

    TLorentzVector Pi0d_Momentum;
    TLorentzVector Kaon_Momentum;
    //int Kch =;
    Pi0d_Momentum = el1 + el2 + gamma;
    Kaon_Momentum = pion + el1 + el2 + gamma;
    dir->fh_Ntracks->Fill(sevt->Ntrack);
    dir->fh_Nclusters->Fill(sevt->Ncluster);
    dir->fh_Nvtx->Fill(sevt->Nvtx);

    dir->fh_gamma_momentum ->Fill(gamma.P());
    dir->fh_pi_momentum    ->Fill(pion.P());
    dir->fh_el1_momentum   ->Fill(el1.P());
    dir->fh_el2_momentum   ->Fill(el2.P());
    dir->fh_k2pi0d_P       ->Fill(Kaon_Momentum.P());
    dir->fh_k2pi0d_Pt      ->Fill(Kaon_Momentum.Pt());
    dir->fh_k2pi0d_M       ->Fill(Kaon_Momentum.M());
    dir->fh_pi0d_P         ->Fill(Pi0d_Momentum.P());
    dir->fh_pi0d_Pt        ->Fill(Pi0d_Momentum.Pt());
    dir->fh_pi0d_M         ->Fill(Pi0d_Momentum.M());

    //dir->fh_muee_M->Fill(Pi0d_Momentum.M());
    //dir->fh_Muee_M_3pi_assumption->Fill(Kaon_Momentum.M());
    //dir->fh_muee_Pt->Fill(Kaon_Momentum.Pt());
    //dir->fh_muee_P->Fill(Kaon_Momentum.P());
}

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

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

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

void swaps::Loop()
{
	gROOT->ProcessLine(".x /home/gcowan/lhcb/lhcbStyle.C");
	if (fChain == 0) return;

	Long64_t nentries = fChain->GetEntriesFast();

	TH1D * mKK  = new TH1D("mKK", "mKK", 100, 0.99, 1.050);
	TH1D * mKpi = new TH1D("mKpi", "mKpi", 100, 0.826, 0.966);
	TH1D * mKp  = new TH1D("mKp", "mKp", 50, 1.45, 1.55);
	TH1D * massHisto  = new TH1D("mJpsiKK_DTF", "mJpsiKK_DTF", 100, 5.20, 5.55);
	TH1D * mBs  = new TH1D("mJpsiKK", "mJpsiKK", 100, 5.20, 5.55);
	TH1D * mBsVeto = new TH1D("mJpsiKKveto", "mJpsiKKveto", 100, 5.20, 5.55);
	//TH1D * mBsKpi = new TH1D("upper_Bs_sideband", "mJpsiKpi", 40, 5.21, 5.41);
	//TH1D * mBspiK = new TH1D("lower_Bs_sideband", "mJpsipiK", 40, 5.04, 5.24);
	TH1D * mBsKpi = new TH1D("upper_Bs_sideband", "mJpsiKpi", 40, 5.51, 5.71);
	TH1D * mBspiK = new TH1D("lower_Bs_sideband", "mJpsipiK", 40, 5.34, 5.54);
	TH1D * mJpsi_ = new TH1D("mJpsi", "mJpsi", 100, 3.03, 3.15);
	TH1D * mJpsi_constr = new TH1D("mJpsi_constr", "mJpsi", 100, 3.03, 3.15);
	
	Long64_t nbytes = 0, nb = 0;

	TFile * file = TFile::Open("reflection_upper_sideband.root", "RECREATE");
	TNtuple * tuple = new TNtuple("DecayTree","DecayTree", "mass");

	for (Long64_t jentry=0; jentry<nentries;jentry++) {
	//for (Long64_t jentry=0; jentry<1000;jentry++) {
		Long64_t ientry = LoadTree(jentry);
		if (ientry < 0) break;
		nb = fChain->GetEntry(jentry);   nbytes += nb;

		if ( !(sel_cleantail==1&&sel==1&&(triggerDecisionUnbiased==1||triggerDecisionBiasedExcl==1)&&time>.3&&time<14&&sigmat>0&&sigmat<0.12
//			&&(Kplus_pidK-Kplus_pidp)>-5
//			&&(Kminus_pidK-Kminus_pidp)>-5
		   )
			) continue;


		//double mpi = 139.57018;
		double mK = 493.68;
		double mmu = 105.658;
		double mJpsi = 3096.916;
		double mpi = 938.27;
		TLorentzVector Kplus(Kplus_PX, Kplus_PY, Kplus_PZ, sqrt(Kplus_PX*Kplus_PX+Kplus_PY*Kplus_PY+Kplus_PZ*Kplus_PZ + mK*mK));
		TLorentzVector Kminus(Kminus_PX, Kminus_PY, Kminus_PZ, sqrt(Kminus_PX*Kminus_PX+Kminus_PY*Kminus_PY+Kminus_PZ*Kminus_PZ + mK*mK));
		TLorentzVector KplusWrong(Kplus_PX, Kplus_PY, Kplus_PZ, sqrt(Kplus_PX*Kplus_PX+Kplus_PY*Kplus_PY+Kplus_PZ*Kplus_PZ + mpi*mpi));
		TLorentzVector KminusWrong(Kminus_PX, Kminus_PY, Kminus_PZ, sqrt(Kminus_PX*Kminus_PX+Kminus_PY*Kminus_PY+Kminus_PZ*Kminus_PZ + mpi*mpi));
		TLorentzVector muplus(muplus_PX, muplus_PY, muplus_PZ, sqrt(muplus_P*muplus_P + mmu*mmu));
		TLorentzVector muminus(muminus_PX, muminus_PY, muminus_PZ, sqrt(muminus_P*muminus_P + mmu*mmu));
		TLorentzVector Jpsi = muplus + muminus;
		TLorentzVector Jpsi_constr(muminus_PX+muplus_PX, muminus_PY+muplus_PY, muminus_PZ+muplus_PZ, sqrt(Jpsi.P()*Jpsi.P() + mJpsi*mJpsi));
		
		TLorentzVector KK = Kplus + Kminus;
		TLorentzVector Kpi = Kplus + KminusWrong;
		TLorentzVector piK = KplusWrong + Kminus;
		TLorentzVector B = Jpsi_constr + KK;
		TLorentzVector BKpi = Jpsi_constr + Kpi;
		TLorentzVector BpiK = Jpsi_constr + piK;
		//if ( ((BpiK.M() > 5600 && BpiK.M() < 5640) || (BKpi.M() > 5600 && BKpi.M() < 5640)) ) mBsVeto->Fill(mass/1000.); 
		//if ( ( Kpi.M() > 1490 && Kpi.M() < 1550 ) || ( piK.M() > 1490 && piK.M() < 1550 ) ) continue;//mBsVeto->Fill(mass/1000.); 
		mKK->Fill(KK.M()/1000.);
		mKp->Fill(Kpi.M()/1000.);
		mKp->Fill(piK.M()/1000.);
		mBs->Fill(B.M()/1000.);
		massHisto->Fill(mass/1000.);
		mJpsi_->Fill(Jpsi.M()/1000.);
		mJpsi_constr->Fill(Jpsi_constr.M()/1000.);
		if (mass > 5420) mBsKpi->Fill(BKpi.M()/1000.);
		if (mass > 5420) mBsKpi->Fill(BpiK.M()/1000.);
		if (mass < 5330) mBspiK->Fill(BKpi.M()/1000.);
		if (mass < 5330) mBspiK->Fill(BpiK.M()/1000.);
		if (mass > 5400) tuple->Fill(BKpi.M());
		if (mass > 5400) tuple->Fill(BpiK.M());
	}

	tuple->Write();
	file->Close();

	std::cout << "Number of B candidates " << mBs->GetEntries() << std::endl;
	std::cout << "Number of B candidates after Lambda_b veto" << mBsVeto->GetEntries() << std::endl;

	gROOT->SetStyle("Plain");
	
	TCanvas * c = new TCanvas("c","c",1600,1200);
	c->Divide(3,2);
	c->cd(1);
	mKK->Draw();
	mKK->SetTitle("");
	mKK->GetXaxis()->SetTitle("m(KK) [GeV/c^{2}]");
	c->cd(2);
	mJpsi_->Draw();
	mJpsi_->SetTitle("");
	mJpsi_->GetXaxis()->SetTitle("m(#mu#mu) [GeV/c^{2}]");
	c->cd(3);
	//mKp->Draw();
	mKp->SetTitle("");
	mKp->GetXaxis()->SetTitle("m(Kp) [GeV/c^{2}]");
	c->cd(4);
	massHisto->Draw();
	//massHisto->SetMaximum(1500);
//.........这里部分代码省略.........

void NewVariables(){

  const double protonmass = 938.272013; //MeV
  const double pionmass = 139.57018; //MeV
  const double kaonmass = 493.677; //MeV
  //const double muonmass = 105.6583715; //MeV

  TStopwatch *clock = new TStopwatch();
  clock->Start(1);

  double p_PT, p_ETA, p_PHI;
  double K_PT, K_ETA, K_PHI;
  double pi_PT, pi_ETA, pi_PHI;
  double Xb_OWNPV_X, Xb_OWNPV_Y, Xb_OWNPV_Z;
  double Xb_ENDVERTEX_X, Xb_ENDVERTEX_Y, Xb_ENDVERTEX_Z;
  double Xb_PT, Xb_ETA, Xb_PHI, Xb_M;
  double Xc_PT, Xc_ETA, Xc_PHI, Xc_M;
  float Added_H_PT[200], Added_H_ETA[200], Added_H_PHI[200];
  int Added_n_Particles;

  gErrorIgnoreLevel = kError;
  TFile *fSLBS = new TFile("/auto/data/mstahl/SLBaryonSpectroscopy/SLBaryonSpectroscopyStrp21.root","read");
  TTree *Xic_tree = (TTree*)gDirectory->Get("Xib02XicMuNu/Xic2pKpi/DecayTree");
  gErrorIgnoreLevel = kPrint;
  Xic_tree->SetBranchStatus("*",0); //disable all branches
  //now switch on the ones we need (saves a lot of time)  
  Xic_tree->SetBranchStatus("Xib_M",1);
  Xic_tree->SetBranchStatus("Xib_PT",1);
  Xic_tree->SetBranchStatus("Xib_ETA",1);
  Xic_tree->SetBranchStatus("Xib_PHI",1);
  Xic_tree->SetBranchStatus("Xib_OWNPV_X",1);
  Xic_tree->SetBranchStatus("Xib_OWNPV_Y",1);
  Xic_tree->SetBranchStatus("Xib_OWNPV_Z",1);
  Xic_tree->SetBranchStatus("Xib_ENDVERTEX_X",1);
  Xic_tree->SetBranchStatus("Xib_ENDVERTEX_Y",1);
  Xic_tree->SetBranchStatus("Xib_ENDVERTEX_Z",1);

  Xic_tree->SetBranchStatus("Xic_M",1);
  Xic_tree->SetBranchStatus("Xic_PT",1);
  Xic_tree->SetBranchStatus("Xic_ETA",1);
  Xic_tree->SetBranchStatus("Xic_PHI",1);

  Xic_tree->SetBranchStatus("Added_n_Particles",1);
  Xic_tree->SetBranchStatus("Added_H_PT",1);
  Xic_tree->SetBranchStatus("Added_H_ETA",1);
  Xic_tree->SetBranchStatus("Added_H_PHI",1);

  Xic_tree->SetBranchStatus("p_PT",1);
  Xic_tree->SetBranchStatus("p_ETA",1);
  Xic_tree->SetBranchStatus("p_PHI",1);
  Xic_tree->SetBranchStatus("K_PT",1);
  Xic_tree->SetBranchStatus("K_ETA",1);
  Xic_tree->SetBranchStatus("K_PHI",1);
  Xic_tree->SetBranchStatus("pi_PT",1);
  Xic_tree->SetBranchStatus("pi_ETA",1);
  Xic_tree->SetBranchStatus("pi_PHI",1);

  //set the branch addresses
  Xic_tree->SetBranchAddress("Xib_M",&Xb_M);
  Xic_tree->SetBranchAddress("Xib_PT",&Xb_PT);
  Xic_tree->SetBranchAddress("Xib_ETA",&Xb_ETA);
  Xic_tree->SetBranchAddress("Xib_PHI",&Xb_PHI);
  Xic_tree->SetBranchAddress("Xib_OWNPV_X",&Xb_OWNPV_X);
  Xic_tree->SetBranchAddress("Xib_OWNPV_Y",&Xb_OWNPV_Y);
  Xic_tree->SetBranchAddress("Xib_OWNPV_Z",&Xb_OWNPV_Z);
  Xic_tree->SetBranchAddress("Xib_ENDVERTEX_X",&Xb_ENDVERTEX_X);
  Xic_tree->SetBranchAddress("Xib_ENDVERTEX_Y",&Xb_ENDVERTEX_Y);
  Xic_tree->SetBranchAddress("Xib_ENDVERTEX_Z",&Xb_ENDVERTEX_Z);

  Xic_tree->SetBranchAddress("Xic_M",&Xc_M);
  Xic_tree->SetBranchAddress("Xic_PT",&Xc_PT);
  Xic_tree->SetBranchAddress("Xic_ETA",&Xc_ETA);
  Xic_tree->SetBranchAddress("Xic_PHI",&Xc_PHI);

  Xic_tree->SetBranchAddress("Added_n_Particles",&Added_n_Particles);
  Xic_tree->SetBranchAddress("Added_H_PT",&Added_H_PT);
  Xic_tree->SetBranchAddress("Added_H_ETA",&Added_H_ETA);
  Xic_tree->SetBranchAddress("Added_H_PHI",&Added_H_PHI);

  Xic_tree->SetBranchAddress("p_PT",&p_PT);
  Xic_tree->SetBranchAddress("p_ETA",&p_ETA);
  Xic_tree->SetBranchAddress("p_PHI",&p_PHI);
  Xic_tree->SetBranchAddress("K_PT",&K_PT);
  Xic_tree->SetBranchAddress("K_ETA",&K_ETA);
  Xic_tree->SetBranchAddress("K_PHI",&K_PHI);
  Xic_tree->SetBranchAddress("pi_PT",&pi_PT);
  Xic_tree->SetBranchAddress("pi_ETA",&pi_ETA);
  Xic_tree->SetBranchAddress("pi_PHI",&pi_PHI);
  //SLBS_tree->AddBranchToCache("*");
  //SLBS_tree->LoadBaskets(1000000000);//Load baskets up to 1 GB to memory

  double Xb_CorrM, p_beta, K_beta, pi_beta;
  float Xcpi_CosTheta[200],XcK_CosTheta[200],Xcp_CosTheta[200];
  double p_as_piKpi_M, p_as_KKpi_M, pK_as_pipi_M, pK_as_ppi_M, pKpi_as_K_M, pKpi_as_p_M;

  TFile *f1 = new TFile("/auto/data/mstahl/SLBaryonSpectroscopy/SLBaryonSpectroscopyStrp21_friend.root","RECREATE");
  //f1->mkdir("Xib02XicMuNu/Xic2pKpi");
  //f1->cd("Xib02XicMuNu/Xic2pKpi");
  TTree added_Xic_tree("Xic2pKpi","Xic2pKpi");

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

void Bd2JpsiKst_reflections::Loop()
{
	gROOT->ProcessLine(".x /Users/gcowan/lhcb/lhcbStyle.C");
	if (fChain == 0) return;

	Long64_t nentries = fChain->GetEntriesFast();

	//TH1D * mKK  = new TH1D("mKK", "mKK", 100, 0.99, 1.050);
	TH1D * mKK  = new TH1D("mKK", "mKK", 100, 1.45, 1.55);
	//TH1D * mKK  = new TH1D("mKK", "mKK", 100, 0.4, 0.6);
	TH1D * mKpi = new TH1D("mKpi", "mKpi", 100, 0.826, 0.966);
	TH1D * massHisto  = new TH1D("mBd_DTF", "mBd_DTF", 100, 5.15, 5.4);
	TH1D * mBd  = new TH1D("mBd", "mBd", 100, 5.15, 5.4);
	TH1D * mBs  = new TH1D("mBs", "mBs", 100, 5.20, 5.55);
	//TH1D * mB_lower = new TH1D("lower_sideband", "lower_sideband", 100, 5.15, 5.4);
	//TH1D * mB_upper = new TH1D("upper_sideband", "upper_sideband", 100, 5.15, 5.4);
	//TH1D * mB_lower = new TH1D("lower_sideband", "lower_sideband", 60, 5.3, 5.45); // for looking at Bs -> JpsiKK where K -> pi misid
	//TH1D * mB_upper = new TH1D("upper_sideband", "upper_sideband", 60, 5.3, 5.45);// for looking at Bs -> JpsiKK where K -> pi misid
	TH1D * mB_lower = new TH1D("lower_sideband", "lower_sideband", 40, 5.55, 5.7);
	TH1D * mB_upper = new TH1D("upper_sideband", "upper_sideband", 40, 5.55, 5.7);
	TH1D * mJpsi_ = new TH1D("mJpsi", "mJpsi", 100, 3.03, 3.15);
	TH1D * mJpsi_constr = new TH1D("mJpsi_constr", "mJpsi", 100, 3.03, 3.15);
	
	Long64_t nbytes = 0, nb = 0;

	TFile * file = TFile::Open("reflection_upper_sideband.root", "RECREATE");
	TNtuple * tuple = new TNtuple("DecayTree","DecayTree", "mass");

	for (Long64_t jentry=0; jentry<nentries;jentry++) {
	//for (Long64_t jentry=0; jentry<1000;jentry++) {
		Long64_t ientry = LoadTree(jentry);
		if (ientry < 0) break;
		nb = fChain->GetEntry(jentry);   nbytes += nb;

		double mpi = 139.57018;
		double mK = 493.68;
		double mmu = 105.658;
		double mJpsi = 3096.916;
		double mp = 938.27;
		TLorentzVector Kplus(K1_Px, K1_Py, K1_Pz, sqrt(K1_Px*K1_Px+K1_Py*K1_Py+K1_Pz*K1_Pz + mK*mK));
		TLorentzVector Piminus(Pi1_Px, Pi1_Py, Pi1_Pz, sqrt(Pi1_Px*Pi1_Px+Pi1_Py*Pi1_Py+Pi1_Pz*Pi1_Pz + mpi*mpi));
		TLorentzVector KplusWrong(K1_Px, K1_Py, K1_Pz, sqrt(K1_Px*K1_Px+K1_Py*K1_Py+K1_Pz*K1_Pz + mpi*mpi));
		//TLorentzVector PiminusWrong(Pi1_Px, Pi1_Py, Pi1_Pz, sqrt(Pi1_Px*Pi1_Px+Pi1_Py*Pi1_Py+Pi1_Pz*Pi1_Pz + mK*mK));
		TLorentzVector PiminusWrong(Pi1_Px, Pi1_Py, Pi1_Pz, sqrt(Pi1_Px*Pi1_Px+Pi1_Py*Pi1_Py+Pi1_Pz*Pi1_Pz + mp*mp));
		TLorentzVector muplus(Mu1_Px, Mu1_Py, Mu1_Pz, sqrt(Mu1_P*Mu1_P + mmu*mmu));
		TLorentzVector muminus(Mu2_Px, Mu2_Py, Mu2_Pz, sqrt(Mu2_P*Mu2_P + mmu*mmu));
		TLorentzVector Jpsi = muplus + muminus;
		TLorentzVector Jpsi_constr(Mu2_Px+Mu1_Px, Mu2_Py+Mu1_Py, Mu2_Pz+Mu1_Pz, sqrt(Jpsi.P()*Jpsi.P() + mJpsi*mJpsi));
		
		TLorentzVector KK = Kplus + PiminusWrong; // testing Bs -> JpsiKK hypothesis
		//TLorentzVector KK = Piminus + KplusWrong; // testing Bd -> Jpsi pipi hypothesis
		TLorentzVector Kpi = Kplus + Piminus;
		TLorentzVector B = Jpsi_constr + Kpi;
		TLorentzVector BKK = Jpsi_constr + KK;
		
		mKK->Fill(KK.M()/1000.);
		mKpi->Fill(Kpi.M()/1000.);
		mBd->Fill(B.M()/1000.);
		massHisto->Fill(Bd_M/1000.);
		mJpsi_->Fill(Jpsi.M()/1000.);
		mJpsi_constr->Fill(Jpsi_constr.M()/1000.);
		if (Bd_M > 5320) mB_upper->Fill(BKK.M()/1000.);
		if (Bd_M < 5230) mB_lower->Fill(BKK.M()/1000.);
		if (Bd_M > 5320) tuple->Fill(BKK.M());
	}

	tuple->Write();
	file->Close();

	std::cout << "Number of B candidates " << mBs->GetEntries() << std::endl;

	TCanvas * c = new TCanvas("c","c",1600,1200);
	c->SetBottomMargin(0);
	c->Divide(3,2, 0.01, 0.01);
	c->cd(1);
	mKK->Draw();
	mKK->SetTitle("");
	//mKK->GetXaxis()->SetTitle("m(KK) [GeV/c^{2}]");
	mKK->GetXaxis()->SetTitle("m(Kp) [GeV/c^{2}]");
	c->cd(2);
	mJpsi_->Draw();
	mJpsi_->SetTitle("");
	mJpsi_->GetXaxis()->SetTitle("m(#mu#mu) [GeV/c^{2}]");
	c->cd(3);
	mKpi->Draw();
	mKpi->SetTitle("");
	mKpi->GetXaxis()->SetTitle("m(K#pi) [GeV/c^{2}]");
	c->cd(4);
	massHisto->Draw();
	//massHisto->SetMaximum(1500);
	mBd->SetLineColor(kRed);
	mBd->Draw("same");
	massHisto->SetTitle("");
	massHisto->GetXaxis()->SetTitle("DTF m(J/#psi K#pi) [GeV/c^{2}]");
	c->cd(5);
	mB_lower->Draw();
	mB_lower->SetTitle("");
	mB_lower->GetXaxis()->SetTitle("m(J/#psi Kp) [GeV/c^{2}]");
	//mB_lower->GetXaxis()->SetTitle("m(J/#psi KK) [GeV/c^{2}]");
	c->cd(6);
//.........这里部分代码省略.........

double KinUtils::doElasticRecoil(const TLorentzVector &chi, TLorentzVector &recoil, TLorentzVector &recoil_chi, const int &procID) {

	TVector3 v0, v1, v2;
	TVector3 p0, pr, pchi;
	double E0, Er, Echi;
	double Tr_max; //this is the maximum recoil KINETIC energy for this incoming chi
	double P0, Pr, Pchi;
	double ctheta_r, stheta_r, phi_r, sigma;
	int ii;
	E0 = chi.E();
	P0 = chi.P();
	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);
//.........这里部分代码省略.........

bool leptonic_fitter_algebraic::fit( const TLorentzVector& B, const TH1& BITF, const TF1& Beff, 
				     const TLorentzVector& lep, 
				     double MEX, double MEY, const TF1& dnuPDF )
{
  if( _dbg > 19 ) cout<<"DBG20 Entered leptonic_fitter_algebraic::fit with B mass: "<<B.M()<<", l_m:"<<lep.M()<<", MET: "<<MEX<<" "<<MEY<<endl;
  if( B.M() <= 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a b-jet with an illegal (non-positive) mass!"); 
  if( lep.M() < 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a lepton with an illegal (negative) mass!"); 
  _converged = _swapped = false;
  _obsB = B;
  _obsL = lep;

  _BITF = &BITF;
  _Beff = &Beff;
  _dnuPDF = dnuPDF;

  _b_m2 = B.M2();

  double lep_b_angle = lep.Angle( B.Vect() );
  double cos_lep_b = TMath::Cos( lep_b_angle );
  double sin_lep_b = TMath::Sin( lep_b_angle );
  double b_p = B.P();
  double b_e = B.E();
  _denom = b_e - cos_lep_b * b_p;
  
  _lep_p = lep.P();
  _x0 = - _W_m2 / ( 2 * _lep_p );
  _y1 = - sin_lep_b * _x0 * b_p / _denom;
  _x1_0 = _x0 * b_e / _denom  -  _y1*_y1 / _x0;
  _Z2_0 = _x0*_x0 - _W_m2 - _y1*_y1;
  if( _dbg > 219 ) cout<<"DBG220 lfa updated lepton with: "<<lv2str( lep )<<" -> x0:"<<_x0<<", y1: "<<_y1<<", x1_0: "<<_x1_0<<", Z2_0: "<<_Z2_0<<endl;

  static double bnums[3];
  bnums[0] = B.X();
  bnums[1] = B.Y();
  bnums[2] = B.Z();
  TMatrixD bXYZ( 3, 1, bnums );
  _R_T = rotation( 2, lep.Phi() ); // R_z^T
  _R_T *= rotation( 1, lep.Theta() - 0.5*TMath::Pi() ); // R_z^T R_y^T
  TMatrixD rotation_vect( _R_T, TMatrixD::kTransposeMult, bXYZ ); // R_y R_z
  double* rotation_array = rotation_vect.GetMatrixArray();
  double phi_x = - TMath::ATan2( rotation_array[2], rotation_array[1] );
  if( _dbg > 99 ) cout<<"DBG100 lfa x rotation vector is:"<<rotation_array[0]<<" "<<rotation_array[1]<<" "<<rotation_array[2]<<" -> phi_x:"<<phi_x<<endl;
  _R_T *= rotation( 0, - phi_x ); // R_z^T R_y^T R_x^T

  // set up _Nu's non-zero elements so that \vec{nu} = Nu \vec{t} for any \vec{t} (since only t's 3nd component is used, and its always 1).
  _Nu[0][2] = MEX;
  _Nu[1][2] = MEY;

  double iVarMET = TMath::Power( TMath::Max( 1., dnuPDF.GetHistogram()->GetRMS() ), -2 );
  _invFlatVar[0][0] = _invFlatVar[1][1] = iVarMET; // set up the chi^2 distance with the right order of magnitude (generalizes to rotated covariance matrix)
  if( _dbg > 209 ) cout<<"DBG210 lfa "<<dnuPDF.GetName()<<" --> iVarMET:"<<iVarMET<<endl;

  // (re)define fit parameter, so all fits start off on an equal footing
  _mini->SetPrintLevel( _minimizer_print_level );
  _mini->Clear();
  _mini->SetFunction( _functor );
  leptonic_fitter_algebraic_object = this; // set the function in the functor pointing back to this object. Doubtfull that all this redirection is needed...
  _mini->SetTolerance( _tolerance );
  bool OK = _mini->SetLimitedVariable( 0, "sB", 1.0, 0.4, 0.1, 6.0 );
  //bool OK = _mini->SetVariable( 0, "sB", 1.0, 0.4 );
  if( ! OK ) {cerr<<"minimizer (@lfa) failed to SetVariable."<<endl; return false;}

  // define 1 sigma in terms of the function
  _mini->SetErrorDef( 0.5 ); // since this is a likelihood fit

  // do the minimization
  OK = _mini->Minimize(); 
  if( _dbg > 19 && ( ! OK || _dbg > 59 ) ) cout<<"DBG INFO: initial fit @lfa returned OK: "<<OK<<", has status: "<<_mini->Status()<<endl;

  _converged = OK; // use status somehow? depends on fitter?

  // read parameters
  const double *xs = _mini->X();
  for( int ip = 0; ip < 1; ++ip ) _params[ ip ] = xs[ ip ];

  // return all intermediate results to the minimum, in particular, the discriminant
  calc_MLL( _params, true );
  TMatrixD nu_vec( _Emat, TMatrixD::kMult, _tvec );
  update_nu_and_decay_chain( nu_vec );
  if( _dbg > 203 ) cout<<"DBG204 lfa finalized _genN: "<<lv2str(_genN)<<", _W: "<<lv2str(_W)<<", & _t: "<<lv2str(_T)<<endl;

  _MLL = _mini->MinValue();
  return true;
} 

float computeMR(TLorentzVector hem1, TLorentzVector hem2) {
    return sqrt(pow(hem1.P() + hem2.P(), 2) - pow(hem1.Pz() + hem2.Pz(), 2));
}

map<string, TH1F *>
readSample (string sampleName, string radice, LHAPDF::PDF * pdf, float referenceScale = 0., int maxevents = -1)
{
  cout << "reading " << sampleName << endl ;
  std::ifstream ifs (sampleName.c_str ()) ;
  LHEF::Reader reader (ifs) ;
  
  map<string, TH1F *> histos ;

  TH1F * h_vbf0_eta    = addHistoToMap (histos, string ("vbf0_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_vbf0_pt     = addHistoToMap (histos, string ("vbf0_pt_")      + radice, 100, 0, 400) ;
  TH1F * h_vbf0_phi    = addHistoToMap (histos, string ("vbf0_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_vbf1_eta    = addHistoToMap (histos, string ("vbf1_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_vbf1_pt     = addHistoToMap (histos, string ("vbf1_pt_")      + radice, 100, 0, 250) ;
  TH1F * h_vbf1_phi    = addHistoToMap (histos, string ("vbf1_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_lep0_eta    = addHistoToMap (histos, string ("lep0_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_lep0_pt     = addHistoToMap (histos, string ("lep0_pt_")      + radice, 100, 0, 400) ;
  TH1F * h_lep0_phi    = addHistoToMap (histos, string ("lep0_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_lep1_eta    = addHistoToMap (histos, string ("lep1_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_lep1_pt     = addHistoToMap (histos, string ("lep1_pt_")      + radice, 100, 0, 250) ;
  TH1F * h_lep1_phi    = addHistoToMap (histos, string ("lep1_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_lep2_eta    = addHistoToMap (histos, string ("lep2_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_lep2_pt     = addHistoToMap (histos, string ("lep2_pt_")      + radice, 100, 0, 250) ;
  TH1F * h_lep2_phi    = addHistoToMap (histos, string ("lep2_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_met_eta    = addHistoToMap (histos, string ("met_eta_")     + radice, 40, -6, 6) ;
  TH1F * h_met_pt     = addHistoToMap (histos, string ("met_pt_")      + radice, 100, 0, 250) ;
  TH1F * h_met_phi    = addHistoToMap (histos, string ("met_phi_")     + radice, 30, -3.14, 3.14) ;
                                                                        
  TH1F * h_mjj_vbf     = addHistoToMap (histos, string ("mjj_vbf_")      + radice, 70, 0, 4000) ;
  TH1F * h_deta_vbf    = addHistoToMap (histos, string ("deta_vbf_")     + radice, 70, 0, 10) ;
                                                                    
  TH1F * h_NJ          = addHistoToMap (histos, string ("NJ_")           + radice, 5, 0, 5) ;
  TH1F * h_NG          = addHistoToMap (histos, string ("NG_")           + radice, 5, 0, 5) ;

  TH1F * h_scale       = addHistoToMap (histos, string ("scale_")        + radice, 100, 0., 500.) ;
  TH1F * h_weight      = addHistoToMap (histos, string ("weight_")       + radice, 100, 0., 10.) ;
 
  int ieve = 0 ;
  // loop over events
  while ( reader.readEvent () ) 
    {
      if (ieve % 10000 == 0) std::cout << "event " << ieve << "\n" ;
      if (maxevents > 0 && ieve >= maxevents) break ;
      ++ieve;
  
      vector<TLorentzVector> finalJets ;      
      vector<TLorentzVector> initialQuarks ;      
      vector<TLorentzVector> finalQuarks ;      
      vector<TLorentzVector> finalGluons ;      
      vector<TLorentzVector> leptons ;      
      vector<TLorentzVector> neutrinos ;      
      
      double x[2] = {0., 0.} ;
      int flavour[2] = {0, 0} ;

      int iquark = 0 ;
      //PG loop over particles in the event
      //PG and fill the variables of leptons and quarks
      for (int iPart = 0 ; iPart < reader.hepeup.IDUP.size (); ++iPart)
        {
//          std::cout << "\t part type [" << iPart << "] " << reader.hepeup.IDUP.at (iPart)
//                    << "\t status " << reader.hepeup.ISTUP.at (iPart)
//                    << "\n" ;

          TLorentzVector particle
                (
                  reader.hepeup.PUP.at (iPart).at (0), //PG px
                  reader.hepeup.PUP.at (iPart).at (1), //PG py
                  reader.hepeup.PUP.at (iPart).at (2), //PG pz
                  reader.hepeup.PUP.at (iPart).at (3) //PG E
                ) ;
          //PG incoming particle          
          if (reader.hepeup.ISTUP.at (iPart) == -1)
            {
               x[iquark] = particle.P () / 7000. ;
               flavour[iquark++] = reader.hepeup.IDUP.at (iPart) ;
               initialQuarks.push_back (particle) ;
            } //PG incoming particle          

          //PG outgoing particles          
          if (reader.hepeup.ISTUP.at (iPart) == 1)
            {
              // jets  
              if (abs (reader.hepeup.IDUP.at (iPart)) < 7) // quarks
                {
                  finalJets.push_back (particle) ;
                  finalQuarks.push_back (particle) ;
                } // quarks
              else if (abs (reader.hepeup.IDUP.at (iPart)) == 21 ) // gluons
                {
                  finalJets.push_back (particle) ;
                  finalGluons.push_back (particle) ;
                } // gluons
              else if (abs (reader.hepeup.IDUP.at (iPart)) == 11 ||
                       abs (reader.hepeup.IDUP.at (iPart)) == 13) // charged leptons (e,u)
//.........这里部分代码省略.........

int main(int argc, char ** argv) 
{
  if(argc < 2)
    {
      cout << "Usage:   " << argv[0] 
           << " input.lhe " << endl ;
      return -1;
    }

  bool REW = false ;
  if (argc > 2) REW = true ;
  cout << REW << endl ;
 
  const int SUBSET = 0;
  const string NAME = "cteq6ll"; //"cteq6l1"

  LHAPDF::initPDFSet(NAME, LHAPDF::LHPDF, SUBSET);
  const int NUMBER = LHAPDF::numberPDF();

  LHAPDF::initPDF (0) ;

  std::ifstream ifs (argv[1]) ;
  LHEF::Reader reader (ifs) ;

  TH1F h_phScale ("h_phScale", "h_phScale", 100, 0, 1000) ;
  TH1F h_scale ("h_scale", "h_scale", 100, 0, 1000) ;
  TH1F h_x ("h_x", "h_x", 100, 0, 1) ;
  TH1F h_weight ("h_weight", "h_weight", 100, 0, 2) ;
  TH1F h_phaseSp ("h_phaseSp", "h_phaseSp", 1000, 0, 10) ;

  TH1F h_ptj1 ("h_ptj1", "h_ptj1", 60, 0, 400) ;
  TH1F h_ptj2 ("h_ptj2", "h_ptj2", 50, 0, 300) ;
  TH1F h_etaj1 ("h_etaj1", "h_etaj1", 40, 0, 6) ;
  TH1F h_etaj2 ("h_etaj2", "h_etaj2", 40, 0, 6) ;
  TH1F h_mjj ("h_mjj", "h_mjj", 50, 0, 3500) ;
  TH1F h_detajj ("h_detajj", "h_detajj", 100, 30, 10) ;
  TH1F h_mll ("h_mll", "h_mll", 40, 0, 300) ;

  h_phScale.Sumw2 () ;
  h_scale.Sumw2 () ;
  h_x.Sumw2 () ;
  h_weight.Sumw2 () ;
  h_phaseSp.Sumw2 () ;
  h_ptj1.Sumw2 () ;
  h_ptj2.Sumw2 () ;
  h_etaj1.Sumw2 () ;
  h_etaj2.Sumw2 () ;
  h_mjj.Sumw2 () ;
  h_detajj.Sumw2 () ;
  h_mll.Sumw2 () ;

  int number_total = 0 ;
  int number_selec = 0 ;
  //PG loop over input events
  while (reader.readEvent ()) 
    {
      if ( reader.outsideBlock.length() ) std::cout << reader.outsideBlock;

      ++number_total ;
      TLorentzVector Higgs;
      int iPartHiggs;
      
      std::vector<int> leptonsFlavour ;      
      std::vector<int> finalJets ;      
      std::vector<int> initialQuarks ;      
      
      vector<TLorentzVector> v_f_jets ; //PG w/ b's
      vector<TLorentzVector> v_f_quarks ; //PG w/o b's
      vector<TLorentzVector> v_f_leptons ;
      vector<TLorentzVector> v_f_neutrinos ;
      vector<TLorentzVector> v_f_intermediate ;
      
      int nele = 0 ;
      int nmu  = 0 ;
      int ntau = 0 ;

      double x[2] = {0., 0.} ;
      int flavour[2] = {0, 0} ;
    
      // loop over particles in the event
      // and fill the variables of leptons and quarks
      for (int iPart = 0 ; iPart < reader.hepeup.IDUP.size (); ++iPart) 
        {
           TLorentzVector dummy = buildP (reader.hepeup, iPart) ;

           // incoming particle          
           if (reader.hepeup.ISTUP.at (iPart) == -1) 
             {
               initialQuarks.push_back (iPart) ;
               x[iPart] = dummy.P () / 4000. ;
               flavour[iPart] = reader.hepeup.IDUP.at (iPart) ;

             } // incoming particle          


           if (reader.hepeup.ISTUP.at (iPart) == 2) 
             {
               if (abs (reader.hepeup.IDUP.at (iPart)) == 6)
                 {
                   v_f_intermediate.push_back (dummy) ;
//.........这里部分代码省略.........