本文整理汇总了C++中TNtuple::Write方法的典型用法代码示例。如果您正苦于以下问题：C++ TNtuple::Write方法的具体用法？C++ TNtuple::Write怎么用？C++ TNtuple::Write使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TNtuple的用法示例。
在下文中一共展示了TNtuple::Write方法的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: calculateTstatistic
double pValueFromPoint2PointDissimilarity(IDataSet * data, IDataSet * mcData)
{
std::cout << "GC: calculating T stat for data" << std::endl;
double T = calculateTstatistic( data, mcData );
char buffer[20];
sprintf( buffer, "Tdata = %f", T );
cout << buffer << endl;
int nPerm = 25;
vector<double> Tvalues = permutation( data, mcData, nPerm );
int count = 0;
vector<double>::iterator Titer;
for ( Titer = Tvalues.begin(); Titer != Tvalues.end(); ++Titer ){
if ( T < *Titer ) count++;
}
double pvalue = count/double(nPerm);
string fileName = ResultFormatter::GetOutputFolder();
fileName.append("/tvalues.root");
TFile * outputFile = new TFile(fileName.c_str(), "RECREATE");
TNtuple * ntuple = new TNtuple("tvalues", "tvalues", "T:Tdata:pvalue");
for ( int i = 0; i < nPerm; i++ ) ntuple->Fill(Tvalues[i], T, pvalue);
ntuple->Write();
outputFile->Close();
delete outputFile;
return pvalue;
}示例2: main
int main(){
std::vector<double> initparams;
initparams.push_back(5.0);
FFF::FitSimple fit(&myNLL,&myMC);
FFF::DataWrapper* data = fit.GenerateData(initparams);
std::cout << "The correct answer is 5.0, with a likelihood of 0.0" << std::endl;
std::vector<double> results = fit.MigradFit(initparams,data);
std::cout << "Migrad: " << results[0] << " L=" << results[1] << std::endl;
TNtuple *lspace = fit.MCMCMapLikelihood(initparams,data);
TFile f("lspace.root","RECREATE");
lspace->Write();
f.Close();
results = fit.MCMCFit(initparams,data);
std::cout << "MCMC: " << results[0] << " L=" << results[1] << std::endl;
results = fit.SAnnealFit(initparams,data);
std::cout << "SAnneal: " << results[0] << " L=" << results[1] << std::endl;
// You can also do normal minuit2 stuff
ROOT::Minuit2::FCNBase* migradfunc = fit.GetMinuit2FCNBase(data);
std::vector<double> verrors;
for (size_t i=0;i<initparams.size();i++)
verrors.push_back(0.0); //FIXME
ROOT::Minuit2::MnUserParameters mnParams(initparams,verrors);
ROOT::Minuit2::MnMigrad migrad(*migradfunc,mnParams);
ROOT::Minuit2::FunctionMinimum theMin = migrad(); //FIXME maxfcn, tol
ROOT::Minuit2::MnUserParameters migradresult = theMin.UserParameters();
ROOT::Minuit2::MnMinos minos(*migradfunc,theMin);
std::pair<double, double> errors = minos(0);
std::cout << "Minos errors:" << std::endl;
std::cout << migradresult.Params()[0] << " +" << errors.second << " -" << -1*errors.first << std::endl;
return 0;
}示例3: testPicoD0EventRead
void testPicoD0EventRead(TString filename)
{
gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
loadSharedLibraries();
gSystem->Load("StPicoDstMaker");
gSystem->Load("StPicoD0Maker");
TFile* f = new TFile(filename.Data());
TTree* T = (TTree*)f->Get("T");
StPicoD0Event* event = new StPicoD0Event();
T->SetBranchAddress("dEvent",&event);
TFile ff("read_test.root","RECREATE");
TNtuple* nt = new TNtuple("nt","","m:pt:eta:phi:theta:"
"decayL:kDca:pDca:dca12:cosThetaStar");
StKaonPion* kp = 0;
for(Int_t i=0;i<100000;++i)
{
T->GetEntry(i);
TClonesArray* arrKPi = event->kaonPionArray();
for(int idx=0;idx<event->nKaonPion();++idx)
{
kp = (StKaonPion*)arrKPi->At(idx);
nt->Fill(kp->m(),kp->pt(),kp->eta(),kp->phi(),kp->pointingAngle(),
kp->decayLength(),kp->kaonDca(),kp->pionDca(),kp->dcaDaughters(),kp->cosThetaStar());
}
}
nt->Write();
ff.Close();
}示例4: main
//.........这里部分代码省略.........
TVector3 vLje; vLje.SetPtEtaPhi(sortedJets[0].pt(), sortedJets[0].eta(), sortedJets[0].phi());
TVector3 vJet;
int kJrl = -1; double dJrl = -1.;
int kTrl = -1; double dTrl = -1.;
for (int j=0; j<sortedJets.size(); j++) {
double dJet = sortedJets[j].pt();
sortedJets[j].set_user_index(-1);
vJet.SetPtEtaPhi(dJet, sortedJets[j].eta(), sortedJets[j].phi());
if (TMath::Abs(vJet.DeltaPhi(vLje))>dCut) { sortedJets[j].set_user_index(1); if (dJet>dJrl) { dJrl = dJet; kJrl = j; } }
if (TMath::Abs(vJet.DeltaPhi(vLtk))>dCut) { sortedJets[j].set_user_index(2); if (dJet>dTrl) { dTrl = dJet; kTrl = j; } }
}
//=============================================================================
TVector3 v1sj, v2sj;
for (int j=0; j<sortedJets.size(); j++) {
Float_t dVar[kVar]; for (Int_t i=0; i<kVar; i++) dVar[i] = -1.;
dVar[kWgt] = 1.; dVar[kXsc] = 1.;
vJet.SetPtEtaPhi(sortedJets[j].pt(), sortedJets[j].eta(), sortedJets[j].phi());
//=============================================================================
dVar[kLje] = vLje.Pt(); if (sortedJets[j].user_index()==1) { dVar[kLjr] = ((kJrl==j) ? 1.5 : 0.5); }
dVar[kLtk] = vLtk.Pt(); if (sortedJets[j].user_index()==2) { dVar[kLtr] = ((kTrl==j) ? 1.5 : 0.5); }
//=============================================================================
dVar[kJet] = sortedJets[j].pt();
dVar[kAje] = sortedJets[j].area();
dVar[kMje] = sortedJets[j].m();
//=============================================================================
fastjet::Filter trimmer(subjDef, fastjet::SelectorPtFractionMin(0.));
fastjet::PseudoJet trimmdJet = trimmer(sortedJets[j]);
std::vector<fastjet::PseudoJet> trimmdSj = trimmdJet.pieces();
double d1sj = -1.; int k1sj = -1;
double d2sj = -1.; int k2sj = -1;
for (int i=0; i<trimmdSj.size(); i++) {
double dIsj = trimmdSj[i].pt(); if (dIsj<0.001) continue;
if (dIsj>d1sj) {
d2sj = d1sj; k2sj = k1sj;
d1sj = dIsj; k1sj = i;
} else if (dIsj>d2sj) {
d2sj = dIsj; k2sj = i;
}
}
//=============================================================================
if (d1sj>0.) {
v1sj.SetPtEtaPhi(d1sj, trimmdSj[k1sj].eta(), trimmdSj[k1sj].phi());
dVar[k1sz] = d1sj;
dVar[k1sA] = trimmdSj[k1sj].area();
dVar[k1sm] = trimmdSj[k1sj].m();
dVar[k1sr] = v1sj.DeltaR(vJet);
}
if (d2sj>0.) {
v2sj.SetPtEtaPhi(d2sj, trimmdSj[k2sj].eta(), trimmdSj[k2sj].phi());
dVar[k2sz] = d2sj;
dVar[k2sA] = trimmdSj[k2sj].area();
dVar[k2sm] = trimmdSj[k2sj].m();
dVar[k2sr] = v2sj.DeltaR(vJet);
}
if ((d1sj>0.) && (d2sj>0.)) dVar[kDsr] = v2sj.DeltaR(v1sj);
nt->Fill(dVar);
}
}
//=============================================================================
delete evt;
ascii_in >> evt;
}
//=============================================================================
file->cd(); nt->Write(); file->Close();
//=============================================================================
TString sXsec = sFile; sXsec.ReplaceAll("out", "xsecs");
file = TFile::Open(Form("%s/xsecs/%s.root",sPath.Data(),sXsec.Data()), "READ");
TH1D *hPtHat = (TH1D*)file->Get("hPtHat"); hPtHat->SetDirectory(0);
TH1D *hWeightSum = (TH1D*)file->Get("hWeightSum"); hWeightSum->SetDirectory(0);
TProfile *hSigmaGen = (TProfile*)file->Get("hSigmaGen"); hSigmaGen->SetDirectory(0);
file->Close();
//=============================================================================
sFile.ReplaceAll("out", "wgt");
file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
hPtHat->Write();
hWeightSum->Write();
hSigmaGen->Write();
file->Close();
//=============================================================================
cout << "DONE" << endl;
//=============================================================================
return 0;
}示例5: buildtupledata
//.........这里部分代码省略.........
foundJ1 ? rawpt[ind1] : NaN,
foundJ1 ? jtpt[ind1] : NaN,
foundJ1 ? jtphi[ind1] : NaN,
foundJ1 ? jteta[ind1] : NaN,
foundJ1 ? discr_csvV1[ind1] : NaN,
foundJ1 ? svtxm[ind1] : NaN,
foundJ1 ? discr_prob[ind1] : NaN,
foundJ1 ? svtxdls[ind1] : NaN,
foundJ1 ? svtxpt[ind1] : NaN,
foundJ1 ? (float)svtxntrk[ind1] : NaN,
foundJ1 ? (float)nsvtx[ind1] : NaN,
foundJ1 ? (float)nselIPtrk[ind1] : NaN,
foundJ2 ? rawpt[ind2] : NaN,
foundJ2 ? jtpt[ind2] : NaN,
foundJ2 ? jtphi[ind2] : NaN,
foundJ2 ? jteta[ind2] : NaN,
foundJ2 ? discr_csvV1[ind2] : NaN,
foundJ2 ? svtxm[ind2] : NaN,
foundJ2 ? discr_prob[ind2] : NaN,
foundJ2 ? svtxdls[ind2] : NaN,
foundJ2 ? svtxpt[ind2] : NaN,
foundJ2 ? (float)svtxntrk[ind2] : NaN,
foundJ2 ? (float)nsvtx[ind2] : NaN,
foundJ2 ? (float)nselIPtrk[ind2] : NaN,
foundJ2 && foundJ1 ? acos(cos(jtphi[ind2]-jtphi[ind1])) : NaN,
foundJ3 ? rawpt[ind3] : NaN,
foundJ3 ? jtpt[ind3] : NaN,
foundJ3 ? jtphi[ind3] : NaN,
foundJ3 ? jteta[ind3] : NaN,
foundJ3 ? discr_csvV1[ind3] : NaN,
foundJ3 ? svtxm[ind3] : NaN,
foundJ3 ? discr_prob[ind3] : NaN,
foundJ3 ? svtxdls[ind3] : NaN,
foundJ3 ? svtxpt[ind3] : NaN,
foundJ3 ? (float)svtxntrk[ind3] : NaN,
foundJ3 ? (float)nsvtx[ind3] : NaN,
foundJ3 ? (float)nselIPtrk[ind3] : NaN,
foundJ3 && foundJ1 ? acos(cos(jtphi[ind3]-jtphi[ind1])) : NaN,
foundJ3 && foundJ2 ? acos(cos(jtphi[ind3]-jtphi[ind2])) : NaN,
foundSL ? (float)SLord : NaN,
foundSL ? rawpt[indSL] : NaN,
foundSL ? jtpt[indSL] : NaN,
foundSL ? jtphi[indSL] : NaN,
foundSL ? jteta[indSL] : NaN,
foundSL ? discr_csvV1[indSL] : NaN,
foundSL ? svtxm[indSL] : NaN,
foundSL ? discr_prob[indSL] : NaN,
foundSL ? svtxdls[indSL] : NaN,
foundSL ? svtxpt[indSL] : NaN,
foundSL ? (float)svtxntrk[indSL] : NaN,
foundSL ? (float)nsvtx[indSL] : NaN,
foundSL ? (float)nselIPtrk[indSL] : NaN,
foundSL && foundJ1 ? acos(cos(jtphi[indSL]-jtphi[ind1])) : NaN};
ntdj->Fill(&vdj[0]);
}
f->Close();
}
}
foutevt->cd();
ntevt->Write();
foutevt->Close();
foutdj->cd();
ntdj->Write();
foutdj->Close();
foutinc->cd();
ntinc->Write();
foutinc->Close();
cout<<endl;
cout<<"Total input entries "<<totentries<<endl;
//making centrality-dependent ntuples
//PutInCbins(outputfolder, code, {{0,40}, {80,200}});
if (PbPb && sample=="bjt"){
auto w = calculateWeightsBjet(outputfilenamedj);
updatePbPbBtriggerweight(outputfilenamedj,w);
updatePbPbBtriggerweight(outputfilenameinc,w);
updatePbPbBtriggerweight(outputfilenameevt,w);
} else {
updateweight(outputfilenamedj);
updateweight(outputfilenameinc);
updateweight(outputfilenameevt);
}
}示例6: TFile
//.........这里部分代码省略.........
hSideband->Fit("fSideband","m q","",fitRangeL,fitRangeH);
hSideband->Fit("fSideband","m q","",fitRangeL,fitRangeH);
hSideband->Fit("fSideband"," q","",fitRangeL,fitRangeH);
double normSideband = hSideband->Integral(0,0.12);
cSideband->cd(2);
hMCPSignal->Fit("fMCPSignal","LL q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","LL q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal"," q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal"," q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","LL q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","LL q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal"," q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal"," q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","m q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","m q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal","m q","",fitRangeL,fitRangeH);
hMCPSignal->Fit("fMCPSignal"," m q","",fitRangeL,fitRangeH);
double normMCPSignal = hMCPSignal->Integral(0,0.12);
cSideband->cd(3);
hMCNPSignal->Fit("fMCNPSignal","LL q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","LL q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal"," q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal"," q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","LL q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","LL q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal"," q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal"," q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","m q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","m q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","m q","",fitRangeL,fitRangeH);
hMCNPSignal->Fit("fMCNPSignal","m q","",fitRangeL,fitRangeH);
double normMCNPSignal = hMCNPSignal->Integral(0,0.12);
string sSideband = Form("(%e)*(exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e)))",fSideband->GetParameter(0),fSideband->GetParameter(1),fSideband->GetParameter(2),fSideband->GetParameter(3),fSideband->GetParameter(4),fSideband->GetParameter(5),fSideband->GetParameter(6),fSideband->GetParameter(7),fSideband->GetParameter(8));
string sMCPSignal = Form("(%e)*(exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e)))",fMCPSignal->GetParameter(0),fMCPSignal->GetParameter(1),fMCPSignal->GetParameter(2),fMCPSignal->GetParameter(3),fMCPSignal->GetParameter(4),fMCPSignal->GetParameter(5),fMCPSignal->GetParameter(6),fMCPSignal->GetParameter(7),fMCPSignal->GetParameter(8));
string sMCNPSignal = Form("(%e)*(exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e))+(%e)*exp((%e)*x+(%e)))",fMCNPSignal->GetParameter(0),fMCNPSignal->GetParameter(1),fMCNPSignal->GetParameter(2),fMCNPSignal->GetParameter(3),fMCNPSignal->GetParameter(4),fMCNPSignal->GetParameter(5),fMCNPSignal->GetParameter(6),fMCNPSignal->GetParameter(7),fMCNPSignal->GetParameter(8));
string function;
function="0*("+sSideband+")+[1]*([2]*"+sMCPSignal+"+(1-[2])*("+sMCNPSignal+"))";
string functionNP;
functionNP="0*("+sSideband+")+[1]*(0*"+sMCPSignal+"+(1-[2])*("+sMCNPSignal+"))";
TCanvas *cFit = new TCanvas("cFit","Fit",600,600);
TF1 *f = new TF1("f",function.c_str());
f->SetParameters(1,0.1,0.9,0,0.1);
f->SetParLimits(0,0,1000);
f->SetParLimits(2,-1,2);
f->SetLineColor(2);
f->SetFillColor(2);
f->SetFillStyle(3001);
f->Draw("same");
hData->SetAxisRange(fitRangeL,fitRangeH,"X");
hData->Fit("f","LL");
hData->Fit("f","LL");
hData->Fit("f","LL");
hData->Fit("f","LL");
hData->Fit("f","m");
hData->SetXTitle("Flight Distance Significance");
hData->SetYTitle("Event Fraction");
TF1 *fNP = new TF1("fNP",functionNP.c_str());
fNP->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2));
fNP->SetRange(fitRangeL,fitRangeH*2);
fNP->SetLineColor(4);
fNP->SetFillStyle(3001);
fNP->SetFillColor(4);
hData->SetStats(0);
hData->Draw("same");
fNP->Draw("same");
cout <<fNP->Integral(3.5,60)/f->Integral(3.5,60)<<endl;
TLegend *leg = new TLegend(0.5,0.7,0.9,0.9);
leg->SetBorderSize(0);
leg->SetFillStyle(0);
leg->AddEntry(hData,"pp data","pl");
leg->AddEntry(f,Form("Prompt Fraction %.1f #pm %.1f %%",100*f->GetParameter(2),100*f->GetParError(2)),"");
leg->AddEntry(f,"Prompt D^{0}","f");
leg->AddEntry(fNP,"Non-Prompt D^{0}","f");
leg->Draw();
cSideband->cd(4);
hMCPSignal->Draw("hist");
hSideband->Draw("hist same");
hMCNPSignal->Draw("hist same");
hData->Draw("same");
nt->Fill(ptMin,ptMax,f->GetParameter(2),f->GetParError(2));
nt->Write();
outf->Write();
// outf->Close();
return f;
}示例7: 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);
//.........这里部分代码省略.........
示例8: main
int main(int argc, char **argv)
{
TFile *f;
TFile *newf;
TNtuple *ntuple;
TNtuple *newntuple;
TH1F *h1;
TF1 *func;
TString Metal;
TString dataLoc;
Int_t dataSL;
if(argc < 5 || argc > 6){
std::cout << "The number of arguments is incorrect" << std::endl;
return 0;
}
dataLoc = argv[1];
PHI_MIN = (Double_t) std::stod(argv[2]);
PHI_MAX = (Double_t) std::stod(argv[3]);
N_PHI = (Int_t) std::stoi(argv[4]);
if(argc >= 6) Metal = argv[5];
dataSL = dataLoc.Length();
if(dataLoc(dataSL-1, 1) != "/"){
dataLoc = dataLoc + "/";
}
std::cout << "The following settings are being used" << std::endl;
std::cout << "Data Directory = " << dataLoc << std::endl;
std::cout << PHI_MIN << " < PHI < " << PHI_MAX << ", N_PHI = " << N_PHI << std::endl;
if(argc == 6) std::cout << "Running only for Metal "<< Metal << std::endl;
else std::cout << "Running all Metals" << std::endl;
Float_t Q2, Xb, Zh, Pt, Phi, Val, Err;
Float_t A, Ac, Acc;
Float_t AErr, AcErr, AccErr;
Float_t ChiSQ;
Int_t nentries, empty;
if(Metal == "") N_METAL = 6;
else N_METAL = 1;
for(Int_t met = 0; met < N_METAL; met++){
if(Metal == ""){
if(met == 0) Metal = "C";
else if(met == 1) Metal = "Fe";
else if(met == 2) Metal = "Pb";
else if(met == 3) Metal = "D_C";
else if(met == 4) Metal = "D_Fe";
else if(met == 5) Metal = "D_Pb";
}
if(dataLoc == "")
f = new TFile("../Gen5DimData/" + Metal + "_5_dim_dist.root", "READ");
else
f = new TFile(dataLoc + Metal + "_5_dim_dist.root", "READ");
ntuple = (TNtuple*) f->Get("fit_data");
nentries = ntuple->GetEntries();
ntuple->SetBranchAddress("Xb", &Xb);
ntuple->SetBranchAddress("Q2", &Q2);
ntuple->SetBranchAddress("Xb", &Xb);
ntuple->SetBranchAddress("Zh", &Zh);
ntuple->SetBranchAddress("Pt", &Pt);
ntuple->SetBranchAddress("Phi", &Phi);
ntuple->SetBranchAddress("Val", &Val);
ntuple->SetBranchAddress("Err", &Err);
newntuple = new TNtuple("AAcAcc_data", "AAcAcc_data", "Q2:Xb:Zh:Pt:A:AErr:Ac:AcErr:Acc:AccErr:ChiSQ");
func = new TF1("fit", "[0]+TMath::Cos(x*TMath::Pi()/180)*[1]+TMath::Cos(2*x*TMath::Pi()/180)*[2]");
newf = new TFile(Metal + "newphihist.root", "RECREATE");
newf->cd();
for(Int_t i = 0; i < nentries; i = i + N_PHI){
ntuple->GetEntry(i);
h1 = new TH1F((const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt),
(const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt), N_PHI, PHI_MIN, PHI_MAX);
empty = 0;
for(Int_t j = 1; j <= N_PHI; j++){
ntuple->GetEntry(i+j-1);
h1->SetBinContent(j, Val);
if(h1->GetBinContent(j) == 0)
empty++;
h1->SetBinError(j, Err*1.04);
}
if(empty <= (N_PHI - MIN_BIN)){
h1->Fit(func, "q");
h1->Write();
if(func->GetNDF() != 0){
ChiSQ = func->GetChisquare();
A = func->GetParameter(0);
AErr = func->GetParError(0);
Ac = func->GetParameter(1);
AcErr = func->GetParError(1);
Acc = func->GetParameter(2);
AccErr = func->GetParError(2);
newntuple->Fill(Q2, Xb, Zh, Pt, A, AErr, Ac, AcErr, Acc, AccErr, ChiSQ);
}
//.........这里部分代码省略.........
示例9: createGlauberTree
//.........这里部分代码省略.........
MyResponse *myresp = new MyResponse;
TFile *outFile = TFile::Open(outFileName, "RECREATE");
outFile->SetCompressionLevel(5);
TDirectory::TContext context(outFile);
TTree *tree = new TTree("glaubertree", "Glauber tree");
tree->Branch("header",&myheader, 32*1024, 99);
tree->Branch("response",&myresp, 32*1024, 99);
TNtuple *ntuple = new TNtuple("gnt", "Glauber ntuple", "npart:ncoll:b");
Double_t etas[] = {-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10};
TH1D *hNEta = new TH1D("hNeta","",12,etas);
TH1D *hEtEta = new TH1D("hEteta","",12,etas);
// create events and fill them
for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
stack->Reset();
hNEta->Reset();
hEtEta->Reset();
genHi->Generate();
AliStack *s = genHi->GetStack();
const TObjArray *parts = s->Particles();
Int_t nents = parts->GetEntries();
for (Int_t i = 0; i<nents; ++i) {
TParticle *p = (TParticle*)parts->At(i);
//p->Print();
TParticlePDG *pdg = p->GetPDG(1);
Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
if (c!=0) {
hNEta->Fill(p->Eta());
hEtEta->Fill(p->Eta(),p->Pt());
}
}
AliGenHijingEventHeader *h = (AliGenHijingEventHeader*)genHi->CollisionGeometry();
myheader->fNATT = nents;
myheader->fEATT = h->TotalEnergy();
myheader->fJATT = h->HardScatters();
myheader->fNT = h->TargetParticipants();
myheader->fNP = h->ProjectileParticipants();
myheader->fN00 = h->NwNw();
myheader->fN01 = h->NwN();
myheader->fN10 = h->NNw();
myheader->fN11 = h->NN();
myheader->fBB = h->ImpactParameter();
myheader->fRP = h->ReactionPlaneAngle();
myheader->fPSn = h->ProjSpectatorsn();
myheader->fPSp = h->ProjSpectatorsp();
myheader->fTSn = h->TargSpectatorsn();
myheader->fTSp = h->TargSpectatorsn();
myresp->fEtch0p = hEtEta->GetBinContent(hEtEta->FindBin(0.5));
myresp->fEtch1p = hEtEta->GetBinContent(hEtEta->FindBin(1.5));
myresp->fEtch2p = hEtEta->GetBinContent(hEtEta->FindBin(2.5));
myresp->fEtch3p = hEtEta->GetBinContent(hEtEta->FindBin(3.5));
myresp->fEtch4p = hEtEta->GetBinContent(hEtEta->FindBin(4.5));
myresp->fEtch5p = hEtEta->GetBinContent(hEtEta->FindBin(5.5));
myresp->fEtchrp = hEtEta->GetBinContent(hEtEta->FindBin(10.5));
myresp->fEtch0n = hEtEta->GetBinContent(hEtEta->FindBin(-0.5));
myresp->fEtch1n = hEtEta->GetBinContent(hEtEta->FindBin(-1.5));
myresp->fEtch2n = hEtEta->GetBinContent(hEtEta->FindBin(-2.5));
myresp->fEtch3n = hEtEta->GetBinContent(hEtEta->FindBin(-3.5));
myresp->fEtch4n = hEtEta->GetBinContent(hEtEta->FindBin(-4.5));
myresp->fEtch5n = hEtEta->GetBinContent(hEtEta->FindBin(-5.5));
myresp->fEtchrn = hEtEta->GetBinContent(hEtEta->FindBin(-10.5));
myresp->fNch0p = hNEta->GetBinContent(hNEta->FindBin(0.5));
myresp->fNch1p = hNEta->GetBinContent(hNEta->FindBin(1.5));
myresp->fNch2p = hNEta->GetBinContent(hNEta->FindBin(2.5));
myresp->fNch3p = hNEta->GetBinContent(hNEta->FindBin(3.5));
myresp->fNch4p = hNEta->GetBinContent(hNEta->FindBin(4.5));
myresp->fNch5p = hNEta->GetBinContent(hNEta->FindBin(5.5));
myresp->fNchrp = hNEta->GetBinContent(hNEta->FindBin(10.5));
myresp->fNch0n = hNEta->GetBinContent(hNEta->FindBin(-0.5));
myresp->fNch1n = hNEta->GetBinContent(hNEta->FindBin(-1.5));
myresp->fNch2n = hNEta->GetBinContent(hNEta->FindBin(-2.5));
myresp->fNch3n = hNEta->GetBinContent(hNEta->FindBin(-3.5));
myresp->fNch4n = hNEta->GetBinContent(hNEta->FindBin(-4.5));
myresp->fNch5n = hNEta->GetBinContent(hNEta->FindBin(-5.5));
myresp->fNchrn = hNEta->GetBinContent(hNEta->FindBin(-10.5));
tree->Fill();
if (ntuple) {
Int_t np = h->TargetParticipants() + h->ProjectileParticipants();
Int_t nc = h->NwNw() + h->NwN() + h->NNw() + h->NN();
Double_t b = h->ImpactParameter();
ntuple->Fill(np,nc,b);
}
} // end of event loop
tree->Write();
ntuple->Write();
outFile->Close();
}示例10: makeSmearedTable
//.........这里部分代码省略.........
if (binHFminusTrunc) parameter = getHFminusByPt(npart);
if (binTracks) parameter = getTracksByGen(npart);
values.push_back(parameter);
}
if(label.compare("b") == 0) sort(values.begin(),values.end(),descend);
else sort(values.begin(),values.end());
//Finding the bin boundaries
txtfile << "-------------------------------------" << endl;
txtfile << label.data() << " based cuts are: " << endl;
txtfile << "(";
int size = values.size();
binboundaries[nbins] = values[size-1];
for(int i = 0; i < nbins; i++) {
int entry = (int)(i*(size/nbins));
if(entry < 0 || i == 0) binboundaries[i] = 0;
else binboundaries[i] = values[entry];
txtfile << binboundaries[i] << ", ";
}
txtfile << binboundaries[nbins] << ")" << endl << "-------------------------------------" << endl;
// Determining Glauber results in various bins
TH2D* hNpart = new TH2D("hNpart","",nbins,binboundaries,40,0,40);
TH2D* hNcoll = new TH2D("hNcoll","",nbins,binboundaries,40,0,40);
TH2D* hNhard = new TH2D("hNhard","",nbins,binboundaries,50,0,50);
TH2D* hb = new TH2D("hb","",nbins,binboundaries,600,0,30);
for(unsigned int iev = 0; iev < Nevents; iev++) {
if( iev % 20000 == 0 ) cout<<"Processing event : " << iev << endl;
t->GetEntry(iev);
if (binHFplusTrunc) parameter = getHFplusByPt(npart);
if (binHFminusTrunc) parameter = getHFminusByPt(npart);
if (binTracks) parameter = getTracksByGen(npart);
hNpart->Fill(parameter,npart);
hNcoll->Fill(parameter,ncoll);
hNhard->Fill(parameter,nhard);
hb->Fill(parameter,b);
int bin = hNpart->GetXaxis()->FindBin(parameter) - 1;
if(bin < 0) bin = 0;
if(bin >= nbins) bin = nbins - 1;
nt->Fill(parameter, et, bin, b, npart, ncoll, nhard);
}
TCanvas *cf = new TCanvas();
TF1* fGaus = new TF1("fb","gaus(0)",0,2);
fitSlices(hNpart,fGaus);
fitSlices(hNcoll,fGaus);
fitSlices(hNhard,fGaus);
fitSlices(hb,fGaus);
TH1D* hNpartMean = (TH1D*)gDirectory->Get("hNpart_1");
TH1D* hNpartSigma = (TH1D*)gDirectory->Get("hNpart_2");
TH1D* hNcollMean = (TH1D*)gDirectory->Get("hNcoll_1");
TH1D* hNcollSigma = (TH1D*)gDirectory->Get("hNcoll_2");
TH1D* hNhardMean = (TH1D*)gDirectory->Get("hNhard_1");
TH1D* hNhardSigma = (TH1D*)gDirectory->Get("hNhard_2");
TH1D* hbMean = (TH1D*)gDirectory->Get("hb_1");
TH1D* hbSigma = (TH1D*)gDirectory->Get("hb_2");
txtfile<<"-------------------------------------"<<endl;
txtfile<<"# Bin NpartMean NpartSigma NcollMean NcollSigma bMean bSigma BinEdge"<<endl;
for(int i = 0; i < nbins; i++){
int ii = nbins-i;
outputTable->table_[i].n_part_mean = hNpartMean->GetBinContent(ii);
outputTable->table_[i].n_part_var = hNpartSigma->GetBinContent(ii);
outputTable->table_[i].n_coll_mean = hNcollMean->GetBinContent(ii);
outputTable->table_[i].n_coll_var = hNcollSigma->GetBinContent(ii);
outputTable->table_[i].b_mean = hbMean->GetBinContent(ii);
outputTable->table_[i].b_var = hbSigma->GetBinContent(ii);
outputTable->table_[i].n_hard_mean = hNhardMean->GetBinContent(ii);
outputTable->table_[i].n_hard_var = hNhardSigma->GetBinContent(ii);
outputTable->table_[i].bin_edge = binboundaries[ii-1];
txtfile << i << " " << hNpartMean->GetBinContent(ii) << " " << hNpartSigma->GetBinContent(ii) << " " << hNcollMean->GetBinContent(ii) << " " << hNcollSigma->GetBinContent(ii) << " " << hbMean->GetBinContent(ii) << " " <<hbSigma->GetBinContent(ii) << " " << binboundaries[ii-1] << " " <<endl;
}
txtfile<<"-------------------------------------"<<endl;
outf->cd();
dir->cd();
outputTable->Write();
nt->Write();
for(int ih = 0; ih < nhist; ih++) {
Npart_vs_PtGenPlusEta4[ih]->Write();
PtGenPlusEta4_vs_HFplusEta4[ih]->Write();
Npart_vs_PtGenMinusEta4[ih]->Write();
PtGenMinusEta4_vs_HFminusEta4[ih]->Write();
Npart_vs_Ngentrk[ih]->Write();
Ngentrk_vs_Ntracks[ih]->Write();
}
outf->Write();
txtfile.close();
}示例11: jettrigger
void jettrigger()
{
TFile *f = new TFile("jettrig.root");
TNtuple *nt = (TNtuple *)f->Get("nt");
//Declaration of leaves types
Float_t pt;
Float_t eta;
Float_t phi;
Float_t mass;
Float_t jt40;
Float_t jt60;
Float_t jt80;
Float_t jt100;
Float_t pscl40;
Float_t pscl60;
Float_t pscl80;
Float_t pscl100;
// Set branch addresses.
nt->SetBranchAddress("pt",&pt);
nt->SetBranchAddress("eta",&eta);
nt->SetBranchAddress("phi",&phi);
nt->SetBranchAddress("mass",&mass);
nt->SetBranchAddress("jt40",&jt40);
nt->SetBranchAddress("jt60",&jt60);
nt->SetBranchAddress("jt80",&jt80);
nt->SetBranchAddress("jt100",&jt100);
nt->SetBranchAddress("pscl40",&pscl40);
nt->SetBranchAddress("pscl60",&pscl60);
nt->SetBranchAddress("pscl80",&pscl80);
nt->SetBranchAddress("pscl100",&pscl100);
TFile *fout = new TFile("jettrig_withweight.root","recreate");
TNtuple *ntout = new TNtuple("ntweight","ntweight","pt:eta:phi:mass:jt40:jt60:jt80:jt100:pscl40:pscl60:pscl80:weightJet:weight12003");
//TTree *ntout=new TTree("ntweight","ntweight");
ntout->SetDirectory(fout);
/* ntout->Branch("pt",&pt);
ntout->Branch("eta",&eta);
ntout->Branch("phi",&phi);
ntout->Branch("mass",&mass);
ntout->Branch("jt40",&jt40);
ntout->Branch("jt60",&jt60);
ntout->Branch("jt80",&jt80);
ntout->Branch("jt100",&jt100);
ntout->Branch("pscl40",&pscl40);
ntout->Branch("pscl60",&pscl60);
ntout->Branch("pscl80",&pscl80);
ntout->Branch("pscl100",&pscl100);
ntout->Branch("weight",&weight);
*/
Long64_t nentries = nt->GetEntries();
cout<<nentries<<endl;
int oneperc = nentries/100;
Long64_t nbytes = 0;
//#ifndef __CINT__
//#pragma omp parallel for ordered schedule(dynamic)
//#endif
for (Long64_t i=0; i<nentries;i++) {
nbytes += nt->GetEntry(i);
if (i % oneperc == 0) cout<<"\r"<<i/oneperc<<"% "<<flush;
float weightJet = 0, weight12003 = 0;
if (jt40 && !jt60 && !jt80 && !jt100) weight12003 = 1/pscl40;
if (jt60 && !jt80 && !jt100) weight12003 = 1;
if (jt80 && !jt100) weight12003 = 1;
if (jt100) weight12003 = 1;
if (jt40 && pt>40 && pt<60) weightJet = pscl40;
if (jt60 && pt>60 && pt<80) weightJet = pscl60;
if (jt80 && pt>80 && pt<100) weightJet = pscl80;
if (jt100 && pt>100) weightJet = pscl100;
//#ifndef __CINT__
//#pragma omp ordered
//#endif
ntout->Fill(pt,eta,phi,mass,jt40,jt60,jt80,jt100,pscl40,pscl60,pscl80,weightJet,weight12003);
}
cout<<endl;
cout<<ntout->GetEntries()<<endl;
ntout->Write();
fout->Close();
f->Close();
f = new TFile("jettrig_withweight.root");
nt = (TNtuple *)f->Get("ntweight");
cout<<nt->GetEntries()<<endl;
f->Close();
}示例12: 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);
//.........这里部分代码省略.........
示例13: plot_moments
void plot_moments(char * filename, char * outputfilename)
{
TFile * f = TFile::Open(filename);
TTree * t = (TTree*)f->Get("partial");
//TF1 * f1 = new TF1("f1", myfunction1, -1, 1, 3);
//TF1 * f2 = new TF1("f2", myfunc, -1, 1, 2);
TF1 * f1 = new TF1("f1","([0] + x*[1] +(3.*x*x-1.)/2.*[2])/5.", -1, 1);
TF1 * f2 = new TF1("f2", "-1./sqrt(2.)/5.*([0]*sqrt(1.-x*x) + [1]*sqrt(3./4.)*2.*sqrt(1.-x*x)*x)", -1, 1);
TF1 * f3 = new TF1("f3","[0]*(1-x*x)/5*sqrt(3./8.)", -1, 1);
f1->SetParameter(0, 10000*0.1); // G020 = 0.1
f1->SetParameter(1, 10000*0.2); // G021 = 0.2
f1->SetParameter(2, 10000*0.1); // G022 = 0.1
f2->SetParameter(0, 10000*0.5); //G121 = 0.5
f2->SetParameter(1, 10000*0.2); //G122 = 0.2
f3->SetParameter(0, 10000*0.3); //G222 = 0.3
f1->SetLineColor(kRed);
f2->SetLineColor(kRed);
f3->SetLineColor(kRed);
f1->SetLineWidth(3);
f2->SetLineWidth(3);
f3->SetLineWidth(3);
TFile * outfile = TFile::Open(outputfilename, "RECREATE");
TNtuple * tup = new TNtuple("coeffs","coeff", "f1_G020:f1_G021:f1_G022:f2_G121:f2_G122:f3_G222:f1_status:f2_status:f3_status");
TH1D * h1 = new TH1D("h1", "h1", 100, -1., 1.);
TH1D * h2 = new TH1D("h2", "h1", 100, -1., 1.);
TH1D * h3 = new TH1D("h3", "h1", 100, -1., 1.);
h1->Sumw2();
h2->Sumw2();
h3->Sumw2();
TCanvas * c = new TCanvas("c","c",1200, 1000);
c->Divide(3,2);
c->cd(1);
t->Draw("cosThetaL","","goff");
c->cd(2);
t->Draw("cosThetaK","","goff");
c->cd(3);
t->Draw("phi","","goff");
c->cd(4);
cout << "Doing Fit 1" << endl;
t->Draw("cosThetaL>>h1", "D002","goff");
TFitResultPtr f1_result = h1->Fit(f1, "RS");
h1->GetXaxis()->SetTitle("cosThetaL (weighted with D_{0,0}^{2})");
f1->Draw("samegoff");
int f1_status = gMinuit->fStatus;
c->cd(5);
cout << "Doing Fit 2" << endl;
t->Draw("cosThetaL>>h2", "D102","goff");
TFitResultPtr f2_result = h2->Fit(f2, "RS");
h2->GetXaxis()->SetTitle("cosThetaL (weighted with D_{1,0}^{2})");
f2->Draw("samegoff");
int f2_status = gMinuit->fStatus;
c->cd(6);
cout << "Doing Fit 3" << endl;
t->Draw("cosThetaL>>h3", "D202","goff");
h3->GetXaxis()->SetTitle("cosThetaL (weighted with D_{2,0}^{2})");
c->Update();
TFitResultPtr f3_result = h3->Fit(f3, "RS");
f3->Draw("samegoff");
int f3_status = gMinuit->fStatus;
tup->Fill(f1->GetParameter(0), f1->GetParameter(1), f1->GetParameter(2)
, f2->GetParameter(0), f2->GetParameter(1)
, f3->GetParameter(0)
, f1_status, f2_status, f3_status
);
tup->Write();
outfile->Close();
//c->SaveAs("B2Kstll_partial_moments_l_2.pdf");
}示例14: main
int main(int argc, char* argv[])
{
int outtype = 1;
if (argc == 2) {
outtype = atoi(argv[1]);
cout << "output type will be " << outtype << endl;
} else if (argc > 2) {
cout << "fakedata takes 0 or 1 arguments" << endl;
exit(1);
}
int status = 0;
TFile *tf = new TFile("sin.root", "recreate");
TTree *tt = new TTree("interF_1", "Just a continuous sine wave");
INTERINFO fi;
tt->Branch("fi", &fi, "i/L:t_ant/D:t_ret/D:x/D:y/D:z/D:ycen/D:zcen/D:rad/D:vx/D:vy/D:vz/D/:theta/D:Ef/D:fW_t/D:phase/D:dphdt/D:omega/D");
//set antenna:
double xend = 45;
TNtuple *runcard = new TNtuple("runcard", "run parameters", "i:repeat:xi:yi:zi:ekin:thetai:phii:mass:charge:phasei:nscatters:n_temp:imp:x_ant:t_del:atten");
float ntarray[18];
double aphase = 0;
ntarray[14] = xend;
runcard->Fill(ntarray);
runcard->Write();
int imax = 440000;
switch (outtype) {
default: //pure sine
cout << "doing sine" << endl;
for (int i = 0; i<imax; i++) {
fi.i = i;
fi.x = 50. * TMath::Sin(2 * TMath::Pi() * i * 0.00007);
fi.Ef = TMath::Sin(2 * TMath::Pi() * i * 0.01);
aphase = fmod((2 * TMath::Pi() * i * 0.00007), (TMath::Pi() * 2));
if ((aphase >= TMath::Pi() / 2.0) && (aphase < 3.0 * TMath::Pi() / 2.0)) {
fi.vx = -1;
} else {
fi.vx = 1;
}
tt->Fill();
}
break;
case 1: //sine with gaps that cause a phase shift
cout << "doing sine with gaps and phase shift per gap" << endl;
int ngaps = 0;
double phase = TMath::Pi() / 3.0;
int counted = 0;
for (int i = 0; i<imax; i++) {
fi.i = i;
fi.x = 50. * TMath::Sin(2 * TMath::Pi() * i * 0.00007);
if ((fi.x > xend) || (fi.x < -xend)){
//fi.Ef = 0;
fi.Ef = TMath::Sin(2 * TMath::Pi() * i * 0.01 + (ngaps * phase));
if (counted == 0) {
ngaps++;
counted = 1;
}
} else {
fi.Ef = TMath::Sin(2 * TMath::Pi() * i * 0.01 + (ngaps * phase));
if (counted == 1) {
counted = 0;
}
}
aphase = fmod((2 * TMath::Pi() * i * 0.00007), (TMath::Pi() * 2));
if ((aphase >= TMath::Pi() / 2.0) && (aphase < 3.0 * TMath::Pi() / 2.0)) {
fi.vx = -1;
} else {
fi.vx = 1;
}
tt->Fill();
}
break;
}
tt->Write();
tf->Close();
return status;
}示例15: main
//.........这里部分代码省略.........
TFile *file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
TNtuple *nt = new TNtuple("nt", "", "fWgt:fJet:fAje:fMje:f1sj:f1sA:f1sr:f1sm:f2sj:f2sA:f2sr:f2sm:fDsr:fIsm");
//=============================================================================
HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
HepMC::GenEvent *evt = ascii_in.read_next_event();
while (evt) {
fjInput.resize(0);
TLorentzVector vPar;
for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
vPar.SetPtEtaPhiM((*p)->momentum().perp(), (*p)->momentum().eta(), (*p)->momentum().phi(), dMass);
if ((TMath::Abs(vPar.Eta())<dCutEtaMax)) {
fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.E()));
}
}
//=============================================================================
fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
TLorentzVector vJet, v1sj, v2sj, vIsj;
for (int j=0; j<selectJets.size(); j++) {
double dJet = selectJets[j].pt();
vJet.SetPtEtaPhiM(dJet, selectJets[j].eta(), selectJets[j].phi(), selectJets[j].m());
fastjet::Filter trimmer(subjDef, fastjet::SelectorPtFractionMin(0.));
fastjet::PseudoJet trimmdJet = trimmer(selectJets[j]);
std::vector<fastjet::PseudoJet> trimmdSj = trimmdJet.pieces();
double d1sj = -1.; int k1sj = -1;
double d2sj = -1.; int k2sj = -1;
for (int i=0; i<trimmdSj.size(); i++) {
double dIsj = trimmdSj[i].pt(); if (dIsj<0.001) continue;
if (dIsj>d1sj) {
d2sj = d1sj; k2sj = k1sj;
d1sj = dIsj; k1sj = i;
} else if (dIsj>d2sj) {
d2sj = dIsj; k2sj = i;
}
}
//=============================================================================
Float_t dVar[] = { 1., dJet, selectJets[j].area(), selectJets[j].m(), d1sj, -1., -1., -1., d2sj, -1., -1., -1., -1., -1. };
if (d1sj>0.) {
v1sj.SetPtEtaPhiM(d1sj, trimmdSj[k1sj].eta(), trimmdSj[k1sj].phi(), trimmdSj[k1sj].m());
dVar[k1sr] = v1sj.DeltaR(vJet);
dVar[k1sm] = trimmdSj[k1sj].m();
dVar[k1sA] = trimmdSj[k1sj].area();
}
if (d2sj>0.) {
v2sj.SetPtEtaPhiM(d2sj, trimmdSj[k2sj].eta(), trimmdSj[k2sj].phi(), trimmdSj[k2sj].m());
dVar[k2sr] = v2sj.DeltaR(vJet);
dVar[k2sm] = trimmdSj[k2sj].m();
dVar[k2sA] = trimmdSj[k2sj].area();
}
if ((d1sj>0.) && (d2sj>0.)) {
vIsj = v1sj + v2sj;
dVar[kIsm] = vIsj.M();
dVar[kDsr] = v2sj.DeltaR(v1sj);
}
nt->Fill(dVar);
}
//=============================================================================
delete evt;
ascii_in >> evt;
}
//=============================================================================
file->cd(); nt->Write(); file->Close();
//=============================================================================
TString sXsec = sFile; sXsec.ReplaceAll("out", "xsecs");
file = TFile::Open(Form("%s/xsecs/%s.root",sPath.Data(),sXsec.Data()), "READ");
TH1D *hPtHat = (TH1D*)file->Get("hPtHat"); hPtHat->SetDirectory(0);
TH1D *hWeightSum = (TH1D*)file->Get("hWeightSum"); hWeightSum->SetDirectory(0);
TProfile *hSigmaGen = (TProfile*)file->Get("hSigmaGen"); hSigmaGen->SetDirectory(0);
file->Close();
//=============================================================================
sFile.ReplaceAll("out", "wgt");
file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
hPtHat->Write();
hWeightSum->Write();
hSigmaGen->Write();
file->Close();
//=============================================================================
cout << "DONE" << endl;
return 0;
}