本文整理汇总了C++中TLorentzVector::Rapidity方法的典型用法代码示例。如果您正苦于以下问题:C++ TLorentzVector::Rapidity方法的具体用法?C++ TLorentzVector::Rapidity怎么用?C++ TLorentzVector::Rapidity使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TLorentzVector
的用法示例。
在下文中一共展示了TLorentzVector::Rapidity方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: AnalyseParticles
void AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
{
const LHEF::HEPEUP &hepeup = reader->hepeup;
Int_t particle;
Double_t signPz, cosTheta;
TLorentzVector momentum;
TRootLHEFParticle *element;
for(particle = 0; particle < hepeup.NUP; ++particle)
{
element = (TRootLHEFParticle*) branch->NewEntry();
element->PID = hepeup.IDUP[particle];
element->Status = hepeup.ISTUP[particle];
element->Mother1 = hepeup.MOTHUP[particle].first;
element->Mother2 = hepeup.MOTHUP[particle].second;
element->ColorLine1 = hepeup.ICOLUP[particle].first;
element->ColorLine2 = hepeup.ICOLUP[particle].second;
element->Px = hepeup.PUP[particle][0];
element->Py = hepeup.PUP[particle][1];
element->Pz = hepeup.PUP[particle][2];
element->E = hepeup.PUP[particle][3];
element->M = hepeup.PUP[particle][4];
momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
cosTheta = TMath::Abs(momentum.CosTheta());
signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
element->PT = momentum.Perp();
element->Eta = (cosTheta == 1.0 ? signPz*999.9 : momentum.Eta());
element->Phi = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());
element->Rapidity = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());
element->LifeTime = hepeup.VTIMUP[particle];
element->Spin = hepeup.SPINUP[particle];
}
}
示例2: loop
//.........这里部分代码省略.........
}
Btypesize[tidx]++;
}
}
if(t!=3)
{
if(ptflag>=0) Bmaxpt[ptflag] = true;
if(probflag>=0) Bmaxprob[probflag] = true;
if(tktkflag>=0) Bbesttktkmass[tktkflag] = true;
if(ptMatchedflag>=0) BmaxptMatched[ptMatchedflag] = true;
if(probMatchedflag>=0) BmaxprobMatched[probMatchedflag] = true;
if(tktkMatchedflag>=0) BbesttktkmassMatched[tktkMatchedflag] = true;
}
if(t==0) nt0->Fill();
else if(t==1) nt1->Fill();
else if(t==2) nt2->Fill();
else if(t==4) nt3->Fill();
else if(t==5) nt5->Fill();
else if(t==6) nt6->Fill();
}
}
ntHlt->Fill();
if(isPbPb) ntHi->Fill();
if(!REAL)
{
Int_t gt=0,sigtype=0;
Int_t gsize=0;
Gsize = 0;
for(int j=0;j<GenInfo_size;j++)
{
if((TMath::Abs(GenInfo_pdgId[j])!=BPLUS_PDGID&&TMath::Abs(GenInfo_pdgId[j])!=BZERO_PDGID&&TMath::Abs(GenInfo_pdgId[j])!=BSUBS_PDGID) && gskim) continue;
Gsize = gsize+1;
Gpt[gsize] = GenInfo_pt[j];
Geta[gsize] = GenInfo_eta[j];
Gphi[gsize] = GenInfo_phi[j];
GpdgId[gsize] = GenInfo_pdgId[j];
bGen->SetPtEtaPhiM(GenInfo_pt[j],GenInfo_eta[j],GenInfo_phi[j],GenInfo_mass[j]);
Gy[gsize] = bGen->Rapidity();
sigtype=0;
for(gt=1;gt<8;gt++)
{
if(signalGen(gt,j))
{
sigtype=gt;
break;
}
}
GisSignal[gsize] = sigtype;
Gmu1pt[gsize] = -1;
Gmu1eta[gsize] = -20;
Gmu1phi[gsize] = -20;
Gmu1p[gsize] = -1;
Gmu2pt[gsize] = -1;
Gmu2eta[gsize] = -20;
Gmu2phi[gsize] = -20;
Gmu2p[gsize] = -1;
Gtk1pt[gsize] = -1;
Gtk1eta[gsize] = -20;
Gtk1phi[gsize] = -20;
Gtk2pt[gsize] = -1;
Gtk2eta[gsize] = -20;
Gtk2phi[gsize] = -20;
if(sigtype!=0)
{
Gmu1pt[gsize] = GenInfo_pt[GenInfo_da1[GenInfo_da1[j]]];
Gmu1eta[gsize] = GenInfo_eta[GenInfo_da1[GenInfo_da1[j]]];
Gmu1phi[gsize] = GenInfo_phi[GenInfo_da1[GenInfo_da1[j]]];
Gmu1p[gsize] = Gmu1pt[gsize]*cosh(Gmu1eta[gsize]);
Gmu2pt[gsize] = GenInfo_pt[GenInfo_da2[GenInfo_da1[j]]];
Gmu2eta[gsize] = GenInfo_eta[GenInfo_da2[GenInfo_da1[j]]];
Gmu2phi[gsize] = GenInfo_phi[GenInfo_da2[GenInfo_da1[j]]];
Gmu2p[gsize] = Gmu2pt[gsize]*cosh(Gmu2eta[gsize]);
if(sigtype==1||sigtype==2)
{
Gtk1pt[gsize] = GenInfo_pt[GenInfo_da2[j]];
Gtk1eta[gsize] = GenInfo_eta[GenInfo_da2[j]];
Gtk1phi[gsize] = GenInfo_phi[GenInfo_da2[j]];
}
else
{
Gtk1pt[gsize] = GenInfo_pt[GenInfo_da1[GenInfo_da2[j]]];
Gtk1eta[gsize] = GenInfo_eta[GenInfo_da1[GenInfo_da2[j]]];
Gtk1phi[gsize] = GenInfo_phi[GenInfo_da1[GenInfo_da2[j]]];
Gtk2pt[gsize] = GenInfo_pt[GenInfo_da2[GenInfo_da2[j]]];
Gtk2eta[gsize] = GenInfo_eta[GenInfo_da2[GenInfo_da2[j]]];
Gtk2phi[gsize] = GenInfo_phi[GenInfo_da2[GenInfo_da2[j]]];
}
}
gsize++;
}
ntGen->Fill();
}
}
outf->Write();
outf->Close();
return 1;
}
示例3: Loop
void oniaEff::Loop(const char* fname, bool ispbpb, bool isPsip)
{
// In a ROOT session, you can do:
// root> .L oniaEff.C
// root> oniaEff t
// root> t.GetEntry(12); // Fill t data members with entry number 12
// root> t.Show(); // Show values of entry 12
// root> t.Show(16); // Read and show values of entry 16
// root> t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
fChain->SetBranchStatus("*",0); // disable all branches
fChain->SetBranchStatus("Centrality",1);
fChain->SetBranchStatus("HLTriggers",1);
fChain->SetBranchStatus("Reco_QQ_*",1);
fChain->SetBranchStatus("Gen_QQ_*",1);
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
TFile *f = new TFile(fname, "RECREATE");
f->cd();
// define the histos
// we need 2 numerators (with and without ctau3d cut), 1 denominator. Let's make them in arrays
TH1F *hnum_centmid = new TH1F("hnum_centmid","hnum_centmid",nbins_centmid,bins_centmid);
TH1F *hnum_centfwd = new TH1F("hnum_centfwd","hnum_centfwd",nbins_centfwd,bins_centfwd);
TH1F *hnum_ptmid = new TH1F("hnum_ptmid","hnum_ptmid",nbins_ptmid,bins_ptmid);
TH1F *hnum_ptfwd = new TH1F("hnum_ptfwd","hnum_ptfwd",nbins_ptfwd,bins_ptfwd);
TH1F *hnumcut_centmid = new TH1F("hnumcut_centmid","hnumcut_centmid",nbins_centmid,bins_centmid);
TH1F *hnumcut_centfwd = new TH1F("hnumcut_centfwd","hnumcut_centfwd",nbins_centfwd,bins_centfwd);
TH1F *hnumcut_ptmid = new TH1F("hnumcut_ptmid","hnumcut_ptmid",nbins_ptmid,bins_ptmid);
TH1F *hnumcut_ptfwd = new TH1F("hnumcut_ptfwd","hnumcut_ptfwd",nbins_ptfwd,bins_ptfwd);
TH1F *hnumptdepcut_centmid = new TH1F("hnumptdepcut_centmid","hnumptdepcut_centmid",nbins_centmid,bins_centmid);
TH1F *hnumptdepcut_centfwd = new TH1F("hnumptdepcut_centfwd","hnumptdepcut_centfwd",nbins_centfwd,bins_centfwd);
TH1F *hnumptdepcut_ptmid = new TH1F("hnumptdepcut_ptmid","hnumptdepcut_ptmid",nbins_ptmid,bins_ptmid);
TH1F *hnumptdepcut_ptfwd = new TH1F("hnumptdepcut_ptfwd","hnumptdepcut_ptfwd",nbins_ptfwd,bins_ptfwd);
TH1F *hden_centmid = new TH1F("hden_centmid","hden_centmid",nbins_centmid,bins_centmid);
TH1F *hden_centfwd = new TH1F("hden_centfwd","hden_centfwd",nbins_centfwd,bins_centfwd);
TH1F *hden_ptmid = new TH1F("hden_ptmid","hden_ptmid",nbins_ptmid,bins_ptmid);
TH1F *hden_ptfwd = new TH1F("hden_ptfwd","hden_ptfwd",nbins_ptfwd,bins_ptfwd);
// also store more finely binned histos
TH1F *hcentfine = new TH1F("hcentfine","hcentfine",200,0,200);
TH2F *hnum2d = new TH2F("hnum2d","hnum2d",48,0,2.4,150,0,30);
TH2F *hden2d = new TH2F("hden2d","hden2d",48,0,2.4,150,0,30);
//nentries = 200;
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
// if (ientry>1000) break;
// skip the event if Gen_QQ_size != 1 for now
b_Gen_QQ_size->GetEntry(ientry); //read only this branch
if (Gen_QQ_size !=1) continue;
// apply gen acceptance cuts
b_Gen_QQ_mupl_4mom->GetEntry(ientry); //read only this branch
b_Gen_QQ_mumi_4mom->GetEntry(ientry); //read only this branch
TLorentzVector *tlvgenpl = (TLorentzVector*) Gen_QQ_mupl_4mom->At(0);
TLorentzVector *tlvgenmi = (TLorentzVector*) Gen_QQ_mumi_4mom->At(0);
//if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) continue;
bool isgenok = true;
if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) isgenok=false;
// fill denominators
b_Gen_QQ_4mom->GetEntry(ientry); //read only this branch
TLorentzVector *tlvgenqq = (TLorentzVector*) Gen_QQ_4mom->At(0);
double genpt = tlvgenqq->Pt();
b_Centrality->GetEntry(ientry); //read only this branch
// check that the dimuon is inside the analysis bins
bool gen_inbin = false;
if (fabs(tlvgenqq->Rapidity())<1.6 && tlvgenqq->Pt()>6.5 && tlvgenqq->Pt()<30.) gen_inbin = true;
if (fabs(tlvgenqq->Rapidity())>1.6 && fabs(tlvgenqq->Rapidity())<2.4 && tlvgenqq->Pt()>3. && tlvgenqq->Pt()<30.) gen_inbin = true;
//if (!gen_inbin) continue;
if (!gen_inbin) isgenok=false;
double weight = ispbpb ? fChain->GetWeight()*findNcoll(Centrality) : 1.;
double weight_DataMC = 0.0;
if(fabs(tlvgenqq->Rapidity())<1.6){
weight_DataMC = 1.659 - 0.0844*genpt + 0.0019*genpt*genpt;
if(ispbpb && !isPsip) weight_DataMC = 1.88 - 0.120*genpt + 0.0033*genpt*genpt;
if(!ispbpb && isPsip) weight_DataMC = 1.059 + 0.0001*genpt - 0.0005*genpt*genpt;
if(ispbpb && isPsip) weight_DataMC = 0.254 + 0.109*genpt - 0.0033*genpt*genpt;
}
//.........这里部分代码省略.........
示例4: Loop
void test_allEvt::Loop()
{
std::map<std::string, TH1D*> h_1d;
std::map<std::string, TH2D*> h_2d;
GoToHXframe = new HTransformToHelicityFrame();
GoToCSframe = new HTransformToCS();
TFile *fin_ewk_syst = new TFile("fout_invariant_mass_qed.root");
TH1D *hvpt_ewk_syst = (TH1D*) fin_ewk_syst->Get("vpt_ewk_p8std_div_ewk_p8bad");
TH1D *hmass_ewk_syst = (TH1D*) fin_ewk_syst->Get("ewk_p8std_div_ewk_p8bad");
TH1D *hmass_ewk_syst_photos = (TH1D*) fin_ewk_syst->Get("ewk_p8photos_div_ewk_p8bad");
// TH1D *hmass=new TH1D("hmass","hmass",150,50,200);
// hmass->Sumw2();
// TH1D *hvpt=new TH1D("hvpt","hvpt",300,0,300);
// hvpt->Sumw2();
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
// for (Long64_t jentry=0; jentry<nentries;jentry++) {
for (Long64_t jentry=0; jentry<10e6;jentry++) {
// for (Long64_t jentry=0; jentry<10e3;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
if(jentry%100000==0) std::cout << "jentry= "<< jentry << "/" << nentries << std::endl;
if(ZGen_PostFSR_mass>50
&& MuPosGenStatus1_pt>25
&& MuNegGenStatus1_pt>25
&& TMath::Abs(MuPosGenStatus1_eta)<2.5
&& TMath::Abs(MuNegGenStatus1_eta)<2.5
){
muPosGen_status1.SetPtEtaPhiM(MuPosGenStatus1_pt,MuPosGenStatus1_eta,MuPosGenStatus1_phi,MuPosGenStatus1_mass);
muNegGen_status1.SetPtEtaPhiM(MuNegGenStatus1_pt,MuNegGenStatus1_eta,MuNegGenStatus1_phi,MuNegGenStatus1_mass);
ZGen_status1 = muPosGen_status1 + muNegGen_status1;
// mtv = sqrt(2*ptl1*ptl2*(1d0-cos(delphi)))
double dphi = MuPosGenStatus1_phi - MuNegGenStatus1_phi;
if ( dphi > TMath::Pi() ) {
dphi -= 2.0*TMath::Pi();
} else if ( dphi <= -TMath::Pi() ) {
dphi += 2.0*TMath::Pi();
}
double ZGen_status1_mT = TMath::Sqrt(2*MuPosGenStatus1_pt*MuNegGenStatus1_pt*(1-TMath::Cos(dphi)));
// cout << "ZGen_status1_mT= " << ZGen_status1_mT << " ZGen_status1.Mt()= " << ZGen_status1.Mt() << endl;
double vpt_rew_factor = hvpt_ewk_syst->Interpolate(ZGen_status1.Pt());
double vmass_std_rew_factor = hmass_ewk_syst->Interpolate(ZGen_status1.M());
double vmass_photos_rew_factor = hmass_ewk_syst_photos->Interpolate(ZGen_status1.M());
// hmass->Fill(ZGen_PostFSR_mass);
common_stuff::plot1D("hmass", ZGen_PostFSR_mass, 1, h_1d, 150,50,200 );
common_stuff::plot1D("hNegleptonpt", muNegGen_status1.Pt(), 1, h_1d, 160,25,65 );
common_stuff::plot1D("hPosleptonpt", muPosGen_status1.Pt(), 1, h_1d, 160,25,65 );
common_stuff::plot1D("hvpt", ZGen_status1.Pt(), 1, h_1d, 300,0,300 );
common_stuff::plot1D("hrapidity", ZGen_status1.Rapidity(), 1, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hNegleptoneta", muNegGen_status1.Eta(), 1, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hPosleptoneta", muPosGen_status1.Eta(), 1, h_1d, 25,-2.5,2.5 );
// common_stuff::plot1D("hmT", ZGen_status1.Mt(), 1, h_1d, 100,50,100 );
common_stuff::plot1D("hmT", ZGen_status1_mT, 1, h_1d, 100,50,100 );
common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std", ZGen_PostFSR_mass, vpt_rew_factor, h_1d, 150,50,200 );
common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std", muNegGen_status1.Pt(), vpt_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std", muPosGen_status1.Pt(), vpt_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std", ZGen_status1.Pt(), vpt_rew_factor, h_1d, 300,0,300 );
common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std", ZGen_status1.Rapidity(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std", muNegGen_status1.Eta(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std", muPosGen_status1.Eta(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
// common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std", ZGen_status1.Mt(), vpt_rew_factor, h_1d, 100,50,100 );
common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std", ZGen_status1_mT, vpt_rew_factor, h_1d, 100,50,100 );
common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std_massStd", ZGen_PostFSR_mass, vpt_rew_factor*vmass_std_rew_factor, h_1d, 150,50,200 );
common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std_massStd", muNegGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std_massStd", muPosGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 300,0,300 );
common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Rapidity(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std_massStd", muNegGen_status1.Eta(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std_massStd", muPosGen_status1.Eta(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
// common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Mt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 100,50,100 );
common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massStd", ZGen_status1_mT, vpt_rew_factor*vmass_std_rew_factor, h_1d, 100,50,100 );
common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std_massPhotos", ZGen_PostFSR_mass, vpt_rew_factor*vmass_photos_rew_factor, h_1d, 150,50,200 );
common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std_massPhotos", muNegGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std_massPhotos", muPosGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 160,25,65 );
common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 300,0,300 );
common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Rapidity(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std_massPhotos", muNegGen_status1.Eta(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std_massPhotos", muPosGen_status1.Eta(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
// common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Mt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 100,50,100 );
common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1_mT, vpt_rew_factor*vmass_photos_rew_factor, h_1d, 100,50,100 );
}
//.........这里部分代码省略.........
示例5: fill
void fill(int const kf, TLorentzVector* b, double weight, TLorentzVector const& p1Mom, TLorentzVector const& p2Mom, TVector3 v00)
{
int const centrality = floor(nCent * gRandom->Rndm());
TVector3 const vertex = getVertex(centrality);
// smear primary vertex
// float const sigmaVertex = sigmaVertexCent[cent];
// TVector3 const vertex(gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex));
v00 += vertex;
// smear momentum
TLorentzVector const p1RMom = smearMom(0, p1Mom);
TLorentzVector const p2RMom = smearMom(0, p2Mom);
// smear position
TVector3 const p1RPos = smearPosData(0, vertex.z(), centrality, p1RMom, v00);
TVector3 const p2RPos = smearPosData(0, vertex.z(), centrality, p2RMom, v00);
// TVector3 const kRPos = smearPos(kMom, kRMom, v00);
// TVector3 const pRPos = smearPos(pMom, pRMom, v00);
// reconstruct
TLorentzVector const rMom = p1RMom + p2RMom;
float const p1Dca = dca(p1Mom.Vect(), v00, vertex);
float const p2Dca = dca(p2Mom.Vect(), v00, vertex);
float const p1RDca = dca(p1RMom.Vect(), p1RPos, vertex);
float const p2RDca = dca(p2RMom.Vect(), p2RPos, vertex);
TVector3 v0;
float const dca12 = dca1To2(p1RMom.Vect(), p1RPos, p2RMom.Vect(), p2RPos, v0);
float const decayLength = (v0 - vertex).Mag();
float const dcaD0ToPv = dca(rMom.Vect(), v0, vertex);
float const cosTheta = (v0 - vertex).Unit().Dot(rMom.Vect().Unit());
float const angle12 = p1RMom.Vect().Angle(p2RMom.Vect());
TLorentzVector p1RMomRest = p1RMom;
TVector3 beta;
beta.SetMagThetaPhi(rMom.Beta(), rMom.Theta(), rMom.Phi());
p1RMomRest.Boost(-beta);
float const cosThetaStar = rMom.Vect().Unit().Dot(p1RMomRest.Vect().Unit());
// save
float arr[100];
int iArr = 0;
arr[iArr++] = centrality;
arr[iArr++] = vertex.X();
arr[iArr++] = vertex.Y();
arr[iArr++] = vertex.Z();
arr[iArr++] = kf;
arr[iArr++] = b->M();
arr[iArr++] = b->Perp();
arr[iArr++] = b->PseudoRapidity();
arr[iArr++] = b->Rapidity();
arr[iArr++] = b->Phi();
arr[iArr++] = v00.X();
arr[iArr++] = v00.Y();
arr[iArr++] = v00.Z();
arr[iArr++] = rMom.M();
arr[iArr++] = rMom.Perp();
arr[iArr++] = rMom.PseudoRapidity();
arr[iArr++] = rMom.Rapidity();
arr[iArr++] = rMom.Phi();
arr[iArr++] = v0.X();
arr[iArr++] = v0.Y();
arr[iArr++] = v0.Z();
arr[iArr++] = dca12;
arr[iArr++] = decayLength;
arr[iArr++] = dcaD0ToPv;
arr[iArr++] = cosTheta;
arr[iArr++] = angle12;
arr[iArr++] = cosThetaStar;
arr[iArr++] = p1Mom.M();
arr[iArr++] = p1Mom.Perp();
arr[iArr++] = p1Mom.PseudoRapidity();
arr[iArr++] = p1Mom.Rapidity();
arr[iArr++] = p1Mom.Phi();
arr[iArr++] = p1Dca;
arr[iArr++] = p1RMom.M();
arr[iArr++] = p1RMom.Perp();
arr[iArr++] = p1RMom.PseudoRapidity();
arr[iArr++] = p1RMom.Rapidity();
arr[iArr++] = p1RMom.Phi();
arr[iArr++] = p1RPos.X();
arr[iArr++] = p1RPos.Y();
arr[iArr++] = p1RPos.Z();
arr[iArr++] = p1RDca;
arr[iArr++] = tpcReconstructed(0,1,centrality,p1RMom);
arr[iArr++] = p2Mom.M();
arr[iArr++] = p2Mom.Perp();
arr[iArr++] = p2Mom.PseudoRapidity();
arr[iArr++] = p2Mom.Rapidity();
arr[iArr++] = p2Mom.Phi();
arr[iArr++] = p2Dca;
//.........这里部分代码省略.........
示例6: fill
void fill(int const kf, TLorentzVector* b, double const weight, TClonesArray& daughters)
{
TLorentzVector kMom;
TLorentzVector p1Mom;
TLorentzVector p2Mom;
int nTrk = daughters.GetEntriesFast();
for (int iTrk = 0; iTrk < nTrk; ++iTrk)
{
TParticle* ptl0 = (TParticle*)daughters.At(iTrk);
switch(abs(ptl0->GetPdgCode()))
{
case 321:
ptl0->Momentum(kMom);
break;
case 211:
if(!p1Mom.P()) ptl0->Momentum(p1Mom);
else ptl0->Momentum(p2Mom);
break;
default:
break;
}
}
daughters.Clear();
// smear and get total momentum
TLorentzVector kRMom = smearMom(kMom,fKaonMomResolution);
TLorentzVector p1RMom = smearMom(p1Mom,fPionMomResolution);
TLorentzVector p2RMom = smearMom(p2Mom,fPionMomResolution);
TLorentzVector rMom = kRMom + p1RMom + p2RMom;
// save
float arr[100];
int iArr = 0;
arr[iArr++] = kf;
arr[iArr++] = weight;
arr[iArr++] = b->M();
arr[iArr++] = b->Perp();
arr[iArr++] = b->PseudoRapidity();
arr[iArr++] = b->Rapidity();
arr[iArr++] = b->Phi();
arr[iArr++] = rMom.M();
arr[iArr++] = rMom.Perp();
arr[iArr++] = rMom.PseudoRapidity();
arr[iArr++] = rMom.Rapidity();
arr[iArr++] = rMom.Phi();
arr[iArr++] = kMom.M();
arr[iArr++] = kMom.Perp();
arr[iArr++] = kMom.PseudoRapidity();
arr[iArr++] = kMom.Rapidity();
arr[iArr++] = kMom.Phi();
arr[iArr++] = kRMom.M();
arr[iArr++] = kRMom.Perp();
arr[iArr++] = kRMom.PseudoRapidity();
arr[iArr++] = kRMom.Rapidity();
arr[iArr++] = kRMom.Phi();
arr[iArr++] = p1Mom.M();
arr[iArr++] = p1Mom.Perp();
arr[iArr++] = p1Mom.PseudoRapidity();
arr[iArr++] = p1Mom.Rapidity();
arr[iArr++] = p1Mom.Phi();
arr[iArr++] = p1RMom.M();
arr[iArr++] = p1RMom.Perp();
arr[iArr++] = p1RMom.PseudoRapidity();
arr[iArr++] = p1RMom.Rapidity();
arr[iArr++] = p1RMom.Phi();
arr[iArr++] = p2Mom.M();
arr[iArr++] = p2Mom.Perp();
arr[iArr++] = p2Mom.PseudoRapidity();
arr[iArr++] = p2Mom.Rapidity();
arr[iArr++] = p2Mom.Phi();
arr[iArr++] = p2RMom.M();
arr[iArr++] = p2RMom.Perp();
arr[iArr++] = p2RMom.PseudoRapidity();
arr[iArr++] = p2RMom.Rapidity();
arr[iArr++] = p2RMom.Phi();
nt->Fill(arr);
}
示例7: RunPidGetterQAEff
void RunPidGetterQAEff()
{
TString PidFrameworkDir = "/lustre/nyx/cbm/users/klochkov/soft/PidFramework/";
gSystem->Load( PidFrameworkDir + "build/libPid");
gStyle->SetOptStat(0000);
TFile *f2 = new TFile("pid_0.root");
TTree *PidTree = (TTree*) f2->Get("PidTree");
TofPidGetter *getter = new TofPidGetter();
TBranch *PidGet = PidTree->GetBranch("TofPidGetter");
PidGet->SetAddress(&getter);
PidGet->GetEntry(0);
Float_t ret[3];
TProfile *hEffP = new TProfile ("hEffP", "", 100, 0, 10);
TProfile *hEffPi = new TProfile ("hEffPi", "", 100, 0, 6);
TProfile *hEffK = new TProfile ("hEffK", "", 100, 0, 5);
TProfile *hEffPSigma = new TProfile ("hEffPSigma", "", 100, 0, 10);
TProfile *hEffPiSigma = new TProfile ("hEffPiSigma", "", 100, 0, 6);
TProfile *hEffKSigma = new TProfile ("hEffKSigma", "", 100, 0, 5);
TProfile2D *hEffPtYP = new TProfile2D ("hEffPtYP", "", 100, -2.5, 2.5, 100, 0, 4);
TProfile2D *hEffPtYK = new TProfile2D ("hEffPtYK", "", 100, -2.5, 2.5, 100, 0, 4);
TProfile2D *hEffPtYPi = new TProfile2D ("hEffPtYPi", "", 100, -2.5, 2.5, 100, 0, 4);
TString InTreeFileName = "/lustre/nyx/cbm/users/dblau/cbm/mc/UrQMD/AuAu/10AGeV/sis100_electron/SC_ON/2016_09_01/tree/11111.root";
TFile *InFile = new TFile(InTreeFileName);
TTree *InTree = (TTree*) InFile->Get("fDataTree");
DataTreeEvent* DTEvent;
InTree -> SetBranchAddress("DTEvent",&DTEvent);
int nevents = 100000;//InTree->GetEntries();
int outputstep = 100;
std::cout << "Entries = " << nevents << std::endl;
for (int j=0;j<nevents;j++)
{
if ( (j+1) % outputstep == 0) std::cout << j+1 << "/" << nevents << "\r" << std::flush;
InTree->GetEntry(j);
Int_t Nmc[3] = {100,100,100};
Int_t Ntof[3] = {0,0,0};
Int_t PdgCode[3] = {2212, 212, 211};
Double_t sigmas [3] = {0,0,0};
for (int i=0;i<DTEvent->GetNTracks(); i++)
{
TLorentzVector v;
DataTreeTrack* track = DTEvent -> GetTrack(i);
DataTreeMCTrack* mctrack = DTEvent -> GetMCTrack(i);
Double_t p = mctrack->GetPt() * TMath::CosH( mctrack->GetEta() );
if (track->GetTOFHitId() < 0)
{
if (mctrack->GetPdgId() == 2212 )
{
v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.9386);
hEffP->Fill ( p, 0 );
hEffPSigma->Fill ( p, 0 );
hEffPtYP -> Fill( v.Rapidity() - 1.52, v.Pt(), 0 );
}
if (mctrack->GetPdgId() == 321 )
{
v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.5);
hEffK->Fill ( p, 0 );
hEffKSigma->Fill ( p, 0 );
hEffPtYK -> Fill( v.Rapidity() - 1.52, v.Pt(), 0 );
}
if (mctrack->GetPdgId() == 211 )
{
v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.14);
hEffPi->Fill ( p, 0 );
hEffPiSigma->Fill ( p, 0);
hEffPtYPi -> Fill( v.Rapidity() - 1.52, v.Pt(), 0 );
}
continue;
}
//
DataTreeTOFHit* toftrack = DTEvent -> GetTOFHit(track->GetTOFHitId());
p = toftrack->GetP();
Double_t m2 = toftrack->GetMass2 ();
Bool_t cut = toftrack->GetBeta() > 0.1 && ( track->GetChiSq(0)/track->GetNDF(0) < 3 ) && p > 1.0 ;
if ( !cut ) continue;
getter->GetBayesProbability (m2, p, ret);
sigmas[0] = getter->GetSigmaProton (m2, p);
sigmas[1] = getter->GetSigmaKaon (m2, p);
sigmas[2] = getter->GetSigmaPion (m2, p);
// std::cout << "pdg = " << mctrack->GetPdgId() << " p = " << p << " Sp = " << sigmas[0] << " Sk = " << sigmas[1]<< " Spi = " << sigmas[2] << std::endl;
//.........这里部分代码省略.........
示例8: createWorkspace
//.........这里部分代码省略.........
const double bordersPT[nPT+1] = {0., 12., 16., 20., 30., 100.};
const double bordersRAP[nRAP+1] = {0., 0.6, 2.};
const double bordersL[nL+1] = {onia::ctVarMin, -0.05, -0.03, -0.02, -0.015, -0.01, -0.005, 0., 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1, onia::ctVarMax};
TH1D* h_ctauerr_2011[nRAP+1][nPT+1][nL+1];
TH1D* h_ctauerr_2012[nRAP+1][nPT+1][nL+1];
for(int iRAP = 0; iRAP < nRAP+1; iRAP++){
for(int iPT = 0; iPT < nPT+1; iPT++){
for(int iL = 0; iL < nL+1; iL++){
h_ctauerr_2011[iRAP][iPT][iL] = (TH1D*)infile->Get(Form("h_ctauerr_2011_rap%d_pt%d_l%d",iRAP, iPT, iL));
h_ctauerr_2012[iRAP][iPT][iL] = (TH1D*)infile->Get(Form("h_ctauerr_2012_rap%d_pt%d_l%d",iRAP, iPT, iL));
}
}
}
cout<<"opened hists"<<endl;
/// Finished reading in 2011 data ctauErr-histos
*/
// loop through events in tree and save them to dataset
for (int ientries = 0; ientries < entries; ientries++) {
numEntriesTotal++;
if (ientries%10000==0) std::cout << "event " << ientries << " of " << entries << std::endl;
tree->GetEntry(ientries);
double M_jpsi =jpsi->M();
double pt_jpsi =jpsi->Pt();
double y_jpsi =jpsi->Rapidity();
double M =chic_rf->M();
//double M =chic->M()-jpsi->M()+onia::MpsiPDG;
double y=chic->Rapidity();
double pt=chic->Pt();
//if (ientries%3==0){
// M_jpsi = gRandom->Uniform(JpsiMass->getMin(), JpsiMass->getMax());
//}
//double JpsictErrRand = h_JpsictErr->GetRandom();
//double JpsictErrRand2 = h_JpsictErr->GetRandom();
//double JpsictMeanRand=0.;
////double pTcorrection=(pt-20.)*0.002;
//
////JpsictErrRand-=pTcorrection;
//if(JpsictErrRand<0) JpsictErrRand=0.001;
////JpsictErrRand2-=pTcorrection;
//if(JpsictErrRand2<0) JpsictErrRand2=0.001;
//
//if (ientries%1000000==0){
// double exponent=0.4;
// JpsictMeanRand=gRandom->Exp(exponent);
//}
//
//
//lifetime = gRandom->Gaus(JpsictMeanRand,0.8*JpsictErrRand);
//lifetimeErr = JpsictErrRand2;
//if (ientries%3==0){
// lifetime = gRandom->Gaus(JpsictMeanRand,1.5*JpsictErrRand);
//}
//
示例9: loopNonpromptBzero
//.........这里部分代码省略.........
int type,flag;
int flagEvt=0;
int offsetHltTree=0;
int testevent=0,testcand=0;
for (Long64_t i=0; i<nentries;i++) {
nbytes += root->GetEntry(i);
flagEvt=0;
while (flagEvt==0)
{
hlt->GetEntry(i+offsetHltTree);
//cout <<offsetHltTree<<" "<<Bfr_HLT_Event<<" "<<EvtInfo_EvtNo<<endl;
if (Bfr_HLT_Event==EvtInfo_EvtNo && Bfr_HLT_Run==EvtInfo_RunNo) flagEvt=1; else offsetHltTree++;
}
if (i%10000==0) cout <<i<<" / "<<nentries<<" offset HLT:"<<offsetHltTree<<endl;
int type1size=0,type2size=0,type3size=0,type4size=0,type5size=0,type6size=0,type7size=0;
float best,best2,temy;
int bestindex,best2index;
size=0;
best=-1;
bestindex=-1;
best2=10000.;
best2index=-1;
for (int j=0;j<BInfo_size;j++)
{
if(BInfo_type[j]>7) continue;
if (ifchannel[BInfo_type[j]-1]!=1) continue;
//skim{{{
b4Pout->SetXYZM(BInfo_px[j],BInfo_py[j],BInfo_pz[j],BInfo_mass[j]);
temy = b4Pout->Rapidity();
if(REAL)
{
if(!(((EvtInfo_RunNo>=210498&&EvtInfo_RunNo<=211256&&abs(temy+0.465)<1.93)||(EvtInfo_RunNo>=211313&&EvtInfo_RunNo<=211631&&abs(temy-0.465)<1.93)))) continue;
}
else
{
if(abs(temy+0.465)>=1.93) continue;
}
if(BInfo_mass[j]<5 || BInfo_mass[j]>6) continue;
if(BInfo_pt[j]<10.) continue;
//}}}
if(BInfo_type[j]==1)
{
//if(TrackInfo_pt[BInfo_rftk1_index[j]]<0.9) continue;
fillTree(bP,bVtx,b4P,j,type1size,KAON_MASS,0,REAL);
if(chi2cl[type1size]>best)
{
best = chi2cl[type1size];
bestindex = type1size;
}
type1size++;
}
}
if(size>0)
{
bestchi2 = bestindex;
isbestchi2[bestindex] = 1;
}
nt0->Fill();
size=0;
best=-1;
示例10: analyzer
//.........这里部分代码省略.........
continue;
/////////////////////////////////////////////////
// Muon Resolution Smearing to match MuScleFit
if(isSignal) // smear only signal because it has muons from higgs
{
if(reco1GenPostFSR.pt<0.)
cout << "Muon 1 Post FSR not valid!\n";
if(reco2GenPostFSR.pt<0.)
cout << "Muon 2 Post FSR not valid!\n";
float ptReco1 = -1.;
float ptReco2 = -1.;
if(runPeriod == "7TeV")
{
ptReco1 = smearPT2011 -> PTsmear(reco1GenPostFSR.pt, reco1GenPostFSR.eta, reco1GenPostFSR.charge, reco1.pt, ISMEAR);
ptReco2 = smearPT2011 -> PTsmear(reco2GenPostFSR.pt, reco2GenPostFSR.eta, reco2GenPostFSR.charge, reco2.pt, ISMEAR);
}
else
{
ptReco1 = smearPT -> PTsmear(reco1GenPostFSR.pt, reco1GenPostFSR.eta, reco1GenPostFSR.charge, reco1.pt, ISMEAR);
ptReco2 = smearPT -> PTsmear(reco2GenPostFSR.pt, reco2GenPostFSR.eta, reco2GenPostFSR.charge, reco2.pt, ISMEAR);
}
TLorentzVector reco1Vec;
TLorentzVector reco2Vec;
reco1Vec.SetPtEtaPhiM(ptReco1,reco1.eta,reco1.phi,MASS_MUON);
reco2Vec.SetPtEtaPhiM(ptReco2,reco2.eta,reco2.phi,MASS_MUON);
TLorentzVector diMuonVec = reco1Vec + reco2Vec;
reco1.pt = ptReco1;
reco2.pt = ptReco2;
recoCandMass = diMuonVec.M();
recoCandPt = diMuonVec.Pt();
recoCandY = diMuonVec.Rapidity();
recoCandPhi = diMuonVec.Phi();
reco1Vec.SetPtEtaPhiM(ptReco1,reco1.eta,reco1.phi,MASS_MUON);
reco2Vec.SetPtEtaPhiM(ptReco2,reco2.eta,reco2.phi,MASS_MUON);
diMuonVec = reco1Vec + reco2Vec;
}
//////////////////////////////////////////
// Muon-related cuts
if (recoCandMass > maxMmm || recoCandMass < minMmm)
continue;
bool muon1PassId = (*muonIdFuncPtr)(reco1);
bool muon2PassId = (*muonIdFuncPtr)(reco2);
if ( !(muon1PassId && muon2PassId))
continue;
bool muon1PassIso = (getPFRelIso(reco1) <= 0.12);
bool muon2PassIso = (getPFRelIso(reco2) <= 0.12);
if ( !(muon1PassIso && muon2PassIso))
continue;
// Order muons by pt
if (reco1.pt < reco2.pt)
{
_MuonInfo tmpMuon = reco1;
reco1 = reco2;
reco1 = tmpMuon;
}
// PU reweight
示例11: createWorkspace
void createWorkspace(const std::string &infilename, int nState, bool correctCtau, bool drawRapPt2D, bool drawPtCPM2D){
gROOT->SetStyle("Plain");
gStyle->SetTitleBorderSize(0);
// Set some strings
const std::string workspacename = "ws_masslifetime",
treename = "selectedData";
// Get the tree from the data file
TFile *f = TFile::Open(infilename.c_str());
TTree *tree = (TTree*)f->Get(treename.c_str());
// Set branch addresses in tree to be able to import tree to roofit
TLorentzVector* jpsi = new TLorentzVector;
tree->SetBranchAddress("JpsiP",&jpsi);
double CPMval = 0;
tree->SetBranchAddress("CPM",&CPMval);
double massErr = 0;
tree->SetBranchAddress("JpsiMassErr",&massErr);
double Vprob = 0;
tree->SetBranchAddress("JpsiVprob",&Vprob);
double lifetime = 0;
tree->SetBranchAddress("Jpsict",&lifetime);
double lifetimeErr = 0;
tree->SetBranchAddress("JpsictErr",&lifetimeErr);
// define variables necessary for J/Psi(Psi(2S)) mass,lifetime fit
RooRealVar* JpsiMass =
new RooRealVar("JpsiMass", "M [GeV]", onia::massMin, onia::massMax);
RooRealVar* JpsiMassErr =
new RooRealVar("JpsiMassErr", "#delta M [GeV]", 0, 5);
RooRealVar* JpsiRap =
new RooRealVar("JpsiRap", "y", -onia::rap, onia::rap);
RooRealVar* JpsiPt =
new RooRealVar("JpsiPt", "p_{T} [GeV]", 0. ,100.);
RooRealVar* JpsiCPM =
new RooRealVar("JpsiCPM", "N_{ch}", 0. ,100.);
RooRealVar* Jpsict =
new RooRealVar("Jpsict", "lifetime [mm]", -1., 2.5);
RooRealVar* JpsictErr =
new RooRealVar("JpsictErr", "Error on lifetime [mm]", 0.0001, 1);
RooRealVar* JpsiVprob =
new RooRealVar("JpsiVprob", "", 0.01, 1.);
// Set bins
Jpsict->setBins(10000,"cache");
Jpsict->setBins(100);
JpsiMass->setBins(100);
JpsictErr->setBins(100);
// The list of data variables
RooArgList dataVars(*JpsiMass,*JpsiMassErr,*JpsiRap,*JpsiPt,*JpsiCPM,*Jpsict,*JpsictErr,*JpsiVprob);
// construct dataset to contain events
RooDataSet* fullData = new RooDataSet("fullData","The Full Data From the Input ROOT Trees",dataVars);
int entries = tree->GetEntries();
cout << "entries " << entries << endl;
// loop through events in tree and save them to dataset
for (int ientries = 0; ientries < entries; ientries++) {
if (ientries%100000==0) std::cout << "event " << ientries << " of " << entries << std::endl;
tree->GetEntry(ientries);
double M =jpsi->M();
double y=jpsi->Rapidity();
double pt=jpsi->Pt();
double cpm=CPMval;
if (M > JpsiMass->getMin() && M < JpsiMass->getMax()
&& massErr > JpsiMassErr->getMin() && massErr < JpsiMassErr->getMax()
&& pt > JpsiPt->getMin() && pt < JpsiPt->getMax()
&& cpm > JpsiCPM->getMin() && cpm < JpsiCPM->getMax()
&& y > JpsiRap->getMin() && y < JpsiRap->getMax()
&& lifetime > Jpsict->getMin() && lifetime < Jpsict->getMax()
&& lifetimeErr > JpsictErr->getMin() && lifetimeErr < JpsictErr->getMax()
&& Vprob > JpsiVprob->getMin() && Vprob < JpsiVprob->getMax()
){
JpsiPt ->setVal(pt);
JpsiCPM ->setVal(cpm);
JpsiRap ->setVal(y);
JpsiMass ->setVal(M);
JpsiMassErr ->setVal(massErr);
JpsiVprob ->setVal(Vprob);
//cout<<"before lifetime correction \n"
// <<"Jpsict: "<<lifetime<<" JpsictErr: "<<lifetimeErr<<endl;
if(correctCtau){
lifetime = lifetime * onia::MpsiPDG / M ;
lifetimeErr = lifetimeErr * onia::MpsiPDG / M ;
Jpsict ->setVal(lifetime);
JpsictErr ->setVal(lifetimeErr);
//cout<<"MpsiPDG: "<<onia::MpsiPDG<<endl;
//cout<<"after lifetime correction \n"
// <<"Jpsict: "<<lifetime<<" JpsictErr: "<<lifetimeErr<<endl;
//.........这里部分代码省略.........
示例12: Loop
//==============================================
void PolMC::Loop(Int_t selDimuType)
{
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntries();
Long64_t countGenEvent = 0;
Long64_t nb = 0;
printf("number of entries = %d\n", (Int_t) nentries);
//loop over the events
for (Long64_t jentry=0; jentry<nentries;jentry++) {
if(jentry % 100000 == 0) printf("event %d\n", (Int_t) jentry);
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry);
//protection against dummy events:
if(muPosPx_Gen < -9900.)
continue;
hGen_StatEv->Fill(0.5);//count all events
//do NOT select on the dimuon type (only a RECO variable)
hGen_StatEv->Fill(1.5);//count all events
Double_t enMuPos = sqrt(muPosPx_Gen*muPosPx_Gen + muPosPy_Gen*muPosPy_Gen + muPosPz_Gen*muPosPz_Gen + jpsi::muMass*jpsi::muMass);
Double_t enMuNeg = sqrt(muNegPx_Gen*muNegPx_Gen + muNegPy_Gen*muNegPy_Gen + muNegPz_Gen*muNegPz_Gen + jpsi::muMass*jpsi::muMass);
TLorentzVector *muPos = new TLorentzVector();
TLorentzVector *muNeg = new TLorentzVector();
muPos->SetPxPyPzE(muPosPx_Gen, muPosPy_Gen, muPosPz_Gen, enMuPos);
muNeg->SetPxPyPzE(muNegPx_Gen, muNegPy_Gen, muNegPz_Gen, enMuNeg);
Double_t etaMuPos = muPos->PseudoRapidity();
Double_t etaMuNeg = muNeg->PseudoRapidity();
Double_t pTMuPos = muPos->Pt();
Double_t pTMuNeg = muNeg->Pt();
//take muons only within a certain eta range
if(TMath::Abs(etaMuPos) > jpsi::etaPS || TMath::Abs(etaMuNeg) > jpsi::etaPS){
// printf("eta(pos. muon) = %f, eta(neg. muon) = %f\n", etaMuPos, etaMuNeg);
continue;
}
hGen_StatEv->Fill(2.5);//count all events
if(pTMuPos < jpsi::pTMuMin && pTMuNeg < jpsi::pTMuMin){
// printf("pT(pos. muon) = %f, pT(neg. muon) = %f\n", pTMuPos, pTMuNeg);
continue;
}
hGen_StatEv->Fill(3.5);
//test according to Gavin's proposal:
//if any of the two muons is within 1.4 < eta < 1.6 AND
//the two muons are close in eta (deltaEta < 0.2)
//reject the dimuon (no matter whether it is "Seagull" or
//"Cowboy"):
// if(TMath::Abs(etaMuPos - etaMuNeg) < 0.2 &&
// ((etaMuPos > 1.4 && etaMuPos < 1.6) || (etaMuNeg > 1.4 && etaMuNeg < 1.6))){
// printf("rejecting the event!\n");
// continue;
// }
//build the invariant mass, pt, ... of the two muons
TLorentzVector *onia = new TLorentzVector();
*onia = *(muPos) + *(muNeg);
Double_t onia_mass = onia->M();
Double_t onia_pt = onia->Pt();
Double_t onia_P = onia->P();
Double_t onia_eta = onia->PseudoRapidity();
Double_t onia_rap = onia->Rapidity();
Double_t onia_phi = onia->Phi();
Double_t onia_mT = sqrt(onia_mass*onia_mass + onia_pt*onia_pt);
// Int_t rapIndex = -1;
// for(int iRap = 0; iRap < 2*jpsi::kNbRapBins; iRap++){
// if(onia_rap > jpsi::rapRange[iRap] && onia_rap < jpsi::rapRange[iRap+1]){
// rapIndex = iRap+1;
// break;
// }
// }
Int_t rapForPTIndex = -1;
for(int iRap = 0; iRap < jpsi::kNbRapForPTBins; iRap++){
if(TMath::Abs(onia_rap) > jpsi::rapForPTRange[iRap] &&
TMath::Abs(onia_rap) < jpsi::rapForPTRange[iRap+1]){
rapForPTIndex = iRap+1;
break;
}
}
Int_t pTIndex = -1;
for(int iPT = 0; iPT < jpsi::kNbPTBins[rapForPTIndex]; iPT++){
if(onia_pt > jpsi::pTRange[rapForPTIndex][iPT] && onia_pt < jpsi::pTRange[rapForPTIndex][iPT+1]){
pTIndex = iPT+1;
break;
}
}
Int_t rapIntegratedPTIndex = -1;
for(int iPT = 0; iPT < jpsi::kNbPTBins[0]; iPT++){
if(onia_pt > jpsi::pTRange[0][iPT] && onia_pt < jpsi::pTRange[0][iPT+1]){
rapIntegratedPTIndex = iPT+1;
//.........这里部分代码省略.........
示例13: dieleAna
//.........这里部分代码省略.........
eventmixer_eff[multbin].addVector(vep_eff,2);
eventmixer_eff[multbin].addVector(vem_eff,3);
vector<pair<HGeantKine *, HGeantKine* > >& pairsVec_eff_multbin = eventmixer_eff[multbin].getMixedVector();
for (int imix = 0; imix < 3; ++imix) {
vector<pair<HGeantKine *, HGeantKine* > > pairsVec;
TString suffix;
switch (imix) {
case 0:
pairsVec = pairsVec_acc;
suffix = TString("");
break;
case 1:
pairsVec = pairsVec_eff;
suffix = TString("_eff_multbin0");
break;
case 2:
pairsVec = pairsVec_eff_multbin;
suffix = TString("_eff_multbin")+TString::Itoa(multbin,10);
break;
}
size = pairsVec.size();
for (Int_t j = 0; j < size; j ++) {
pair<HGeantKine*,HGeantKine*>& pair = pairsVec[j];
kine1 = pair.first;
kine2 = pair.second;
TLorentzVector vec1, vec2;
HParticleTool::getTLorentzVector(kine1,vec1,kine1->getID());
HParticleTool::getTLorentzVector(kine2,vec2,kine2->getID());
Float_t mom1 = vec1.Vect().Mag();
Float_t the1 = vec1.Theta()*TMath::RadToDeg();
Float_t mom2 = vec2.Vect().Mag();
Float_t the2 = vec2.Theta()*TMath::RadToDeg();
TLorentzVector dilep = vec1 + vec2;
Float_t oAngle = vec1.Angle(vec2.Vect())*TMath::RadToDeg();
Float_t mass = dilep.M()/1000;
Float_t pt = dilep.Perp();
Float_t y = dilep.Rapidity();
Float_t mbinw = hM.get("hmassNP")->GetBinWidth(hM.get("hmassNP")->FindBin(mass));
if (oAngle < 9) continue;
TString chg = "NP";
if (kine1->getID() == kine2->getID()) {
if (kine1->getID() == 2) chg = "PP";
if (kine1->getID() == 3) chg = "NN";
}
hM.get( TString("hmass") +chg+suffix)->Fill(mass, 1./mbinw);
hM.get( TString("hoAngle") +chg+suffix)->Fill(oAngle );
hM.get( TString("hy") +chg+suffix)->Fill(y );
hM.get( TString("hpt") +chg+suffix)->Fill(pt );
hM.get2(TString("hoAnglemass")+chg+suffix)->Fill(oAngle,mass,1./mbinw);
hM.get2(TString("hoAnglept") +chg+suffix)->Fill(oAngle,pt );
hM.get2(TString("hmasspt") +chg+suffix)->Fill(mass,pt, 1./mbinw);
hM.get2(TString("hoAngley") +chg+suffix)->Fill(oAngle,y );
hM.get2(TString("hmassy") +chg+suffix)->Fill(mass,y, 1./mbinw);
hM.get2(TString("hpty") +chg+suffix)->Fill(pt,y );
hM.get2(TString("hth1th2") +chg+suffix)->Fill(the1,the2 );
hM.get2(TString("hp1p2") +chg+suffix)->Fill(mom1,mom2 );
}
}
//#define DELETE_MIX
#ifdef DELETE_MIX
vector <HGeantKine *>* toDel = eventmixer.getObjectsToDelete();
for (unsigned int ii = 0; ii < toDel->size(); ++ii) {
delete toDel->at(ii);
}
toDel->clear();
delete toDel;
vector <HGeantKine *>* toDel_eff = eventmixer_eff[0].getObjectsToDelete();
for (unsigned int ii = 0; ii < toDel_eff->size(); ++ii) {
delete toDel_eff->at(ii);
}
toDel_eff->clear();
delete toDel_eff;
vector <HGeantKine *>* toDel_eff_multbin = eventmixer_eff[multbin].getObjectsToDelete();
for (unsigned int ii = 0; ii < toDel_eff_multbin->size(); ++ii) {
delete toDel_eff_multbin->at(ii);
}
toDel_eff_multbin->clear();
delete toDel_eff_multbin;
#endif
} // end event loop
timer.Stop();
hM.getFile()->cd();
TMacro m1(__DIELEANA_FILE__);
m1.Write();
hM.writeHists("nomap");
cout<<"####################################################"<<endl;
return 0;
}
示例14: MakeTree_2015
//.........这里部分代码省略.........
TClonesArray *Reco_QQ_vtx = 0;;
TClonesArray *Reco_QQ_4mom = 0;
TClonesArray *Reco_QQ_mupl_4mom = 0;
TClonesArray *Reco_QQ_mumi_4mom = 0;
TClonesArray *Reco_mu_4mom = 0;
Int_t Reco_mu_size;
Int_t Reco_mu_type[NMAX];
Int_t Reco_mu_charge[NMAX];
int nPlusMu;
int nMinusMu;
float muPt[NMAX];
float muEta[NMAX];
float muPhi[NMAX];
//Int_t nReco_QQ;
Float_t invariantMass;
Int_t QQtrig;
Int_t QQsign;
Int_t QQTrkDr03;
Int_t QQTrkDr04;
Int_t QQTrkDr05;
Int_t QQTrkPt04;
Int_t QQTrkPt03;
Int_t QQTrkPt02;
float weight;
float weight2;
Float_t upsPt;
Float_t upsEta;
Float_t upsPhi;
Float_t upsRapidity;
Float_t vProb;
Float_t muPlusPt;
Float_t muMinusPt;
Float_t muPlusEta;
Float_t muMinusEta;
Float_t muPlusPhi;
Float_t muMinusPhi;
Float_t RmuPlusPhi = 0;
Float_t RmuMinusPhi = 0;
// centrliaty extra stuff
Int_t Npix, NpixelTracks,Ntracks;
Float_t SumET_HF, SumET_HFplus, SumET_HFminus, SumET_HFplusEta4, SumET_HFminusEta4, SumET_EB, SumET_ET, SumET_EE, SumET_EEplus, SumET_EEminus, SumET_ZDC, SumET_ZDCplus, SumET_ZDCminus;
// track extra stuff
const int NTRKMAX=100000000;
int Reco_trk_size=0;
TClonesArray *Reco_trk_vtx = 0;
TClonesArray *Reco_trk_4mom = 0;
float trkPt[NTRKMAX];
float trkEta[NTRKMAX];
float trkPhi[NTRKMAX];
//---------------------------
TBranch *b_Centrality; //!
TBranch *b_eventNb; //!
TBranch *b_runNb; //!
TBranch *b_HLTriggers; //!
TBranch *b_Reco_QQ_size; //!
TBranch *b_Reco_QQ_trig; //!
TBranch *b_Reco_QQ_type; //!
TBranch *b_Reco_QQ_sign; //!
TBranch *b_Reco_QQ_VtxProb; //!
示例15: Loop
void ljpsiresol::Loop(const char* fname, bool isdata, bool ispbpb)
{
// In a ROOT session, you can do:
// root> .L ljpsiresol.C
// root> ljpsiresol t
// root> t.GetEntry(12); // Fill t data members with entry number 12
// root> t.Show(); // Show values of entry 12
// root> t.Show(16); // Read and show values of entry 16
// root> t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
fChain->SetBranchStatus("*",0); // disable all branches
fChain->SetBranchStatus("HLTriggers",1);
fChain->SetBranchStatus("Centrality",1);
fChain->SetBranchStatus("Reco_QQ_*",1);
fChain->SetBranchStatus("Gen_QQ_*",!isdata);
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
TFile *f = new TFile(fname, "RECREATE");
f->cd();
// define the histos
// we need 2 numerators (with and without ctau3d cut), 1 denominator. Let's make them in arrays
// TH1F *hnum_centmid = new TH1F("hnum_centmid","hnum_centmid",nbins_centmid,bins_centmid);
vector<TH1F*> hists_mid, hists_fwd, hists_midcent, hists_fwdcent;
for (int i=0; i<nbins_ptmid; i++) {
double ptmin = bins_ptmid[i];
double ptmax = bins_ptmid[i+1];
TH1F *hist = new TH1F(Form("hmid_pt%.1f-%.1f",ptmin,ptmax),"",70,-0.15,0.2);
hists_mid.push_back(hist);
}
for (int i=0; i<nbins_ptfwd; i++) {
double ptmin = bins_ptfwd[i];
double ptmax = bins_ptfwd[i+1];
TH1F *hist = new TH1F(Form("hfwd_pt%.1f-%.1f",ptmin,ptmax),"",70,-0.15,0.2);
hists_fwd.push_back(hist);
}
for (int i=0; i<nbins_centmid; i++) {
double centmin = bins_centmid[i];
double centmax = bins_centmid[i+1];
TH1F *hist = new TH1F(Form("hmid_cent%.1f-%.1f",centmin,centmax),"",70,-0.15,0.2);
hists_midcent.push_back(hist);
}
for (int i=0; i<nbins_centfwd; i++) {
double centmin = bins_centfwd[i];
double centmax = bins_centfwd[i+1];
TH1F *hist = new TH1F(Form("hfwd_cent%.1f-%.1f",centmin,centmax),"",70,-0.15,0.2);
hists_fwdcent.push_back(hist);
}
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
// if (ientry>1000) break;
// skip the event if Gen_QQ_size != 1 for now
TLorentzVector *tlvgenpl=NULL, *tlvgenmi=NULL, *tlvgenqq=NULL;
if (!isdata) {
b_Gen_QQ_size->GetEntry(ientry); //read only this branch
if (Gen_QQ_size !=1) continue;
// apply gen acceptance cuts
b_Gen_QQ_mupl_4mom->GetEntry(ientry); //read only this branch
b_Gen_QQ_mumi_4mom->GetEntry(ientry); //read only this branch
tlvgenpl = (TLorentzVector*) Gen_QQ_mupl_4mom->At(0);
tlvgenmi = (TLorentzVector*) Gen_QQ_mumi_4mom->At(0);
if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) continue;
// fill denominators
b_Gen_QQ_4mom->GetEntry(ientry); //read only this branch
TLorentzVector *tlvgenqq = (TLorentzVector*) Gen_QQ_4mom->At(0);
double genpt = tlvgenqq->Pt();
b_Centrality->GetEntry(ientry); //read only this branch
// check that the dimuon is inside the analysis bins
bool gen_inbin = false;
if (fabs(tlvgenqq->Rapidity())<1.6 && tlvgenqq->Pt()>6.5 && tlvgenqq->Pt()<30.) gen_inbin = true;
if (fabs(tlvgenqq->Rapidity())>1.6 && fabs(tlvgenqq->Rapidity())<2.4 && tlvgenqq->Pt()>3. && tlvgenqq->Pt()<30.) gen_inbin = true;
if (!gen_inbin) continue;
}
double weight = ispbpb ? fChain->GetWeight()*findNcoll(Centrality) : 1.;
// is there at least one reco dimuon?
b_Reco_QQ_size->GetEntry(ientry);
//.........这里部分代码省略.........