本文整理汇总了C++中TLorentzVector::P方法的典型用法代码示例。如果您正苦于以下问题:C++ TLorentzVector::P方法的具体用法?C++ TLorentzVector::P怎么用?C++ TLorentzVector::P使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TLorentzVector
的用法示例。
在下文中一共展示了TLorentzVector::P方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: CalcDoubleDphiRFRAME
// 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());
}
示例2: FillK2pi
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());
}
示例3: DiTau_InvMass
//___________________[ 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;
}//-----------------------------------------------------------------------------------------------------------
示例4: DiTau_InvMass
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;
}
示例5: CalcDeltaPhiRFRAME
// 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);
}
示例6: cosThetaBoost
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;
}
示例7: Boost_To_Stop_Rest_Frame
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());
}
示例8: Loop
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);
//.........这里部分代码省略.........
示例9: NewVariables
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");
//.........这里部分代码省略.........
示例10: Kplus
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);
//.........这里部分代码省略.........
示例11: doElasticRecoil
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);
//.........这里部分代码省略.........
示例12: fit
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;
}
示例13: computeMR
float computeMR(TLorentzVector hem1, TLorentzVector hem2) {
return sqrt(pow(hem1.P() + hem2.P(), 2) - pow(hem1.Pz() + hem2.Pz(), 2));
}
示例14: ifs
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)
//.........这里部分代码省略.........
示例15: main
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) ;
//.........这里部分代码省略.........