本文整理汇总了C++中TFile::Get方法的典型用法代码示例。如果您正苦于以下问题:C++ TFile::Get方法的具体用法?C++ TFile::Get怎么用?C++ TFile::Get使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TFile
的用法示例。
在下文中一共展示了TFile::Get方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: PP_MC_normResponseMatrix
void PP_MC_normResponseMatrix(int radius = 3){
gStyle->SetOptStat(0);
bool printDebug = true;
TFile * fin = TFile::Open(Form("pp_MC_new_ak%d_20_eta_20.root",radius));
TH2F * mPP_Matrix;
TH2F * mPP_Response;
TH2F * mPP_ResponseNorm;
TH2F * mPP_ResponseNorm_recoTrunc;
// TH2F * mPbPb_Matrix_recoTrunc[nbins_cent];
// TH2F * mPbPb_Response_recoTrunc[nbins_cent];
// TH2F * mPbPb_ResponseNorm_recoTrunc[nbins_cent];
TFile fout(Form("PP_kineticEfficiency_correctionFactors_R%d.root",radius),"RECREATE");
fout.cd();
TH1F * hGen_Projection;
TH1F * hGen_Projection_recoTrunc;
TH1F * hKineticEfficiency;
TCanvas * cGenProj = new TCanvas("cGenProj","",800,600);
TCanvas * cKineticEfficiency = new TCanvas("cKineticEfficiency","",800,600);
TCanvas * cMatrix = new TCanvas("cMatrix","",800,600);
TCanvas * cMatrix_recoTrunc = new TCanvas("cMatrix_recoTrunc","",800,600);
mPP_Matrix = (TH2F*)fin->Get(Form("hpp_matrix_HLT_R%d_20_eta_20",radius));
mPP_Matrix->Rebin2D(10, 10);
for (int y=1;y<=mPP_Matrix->GetNbinsY();y++) {
double sum=0;
for (int x=1;x<=mPP_Matrix->GetNbinsX();x++) {
if (mPP_Matrix->GetBinContent(x,y)<=1*mPP_Matrix->GetBinError(x,y)) {
mPP_Matrix->SetBinContent(x,y,0);
mPP_Matrix->SetBinError(x,y,0);
}
sum+=mPP_Matrix->GetBinContent(x,y);
}
for (int x=1;x<=mPP_Matrix->GetNbinsX();x++) {
double ratio = 1;
// if (hGenSpectraCorrPP->GetBinContent(x)!=0) ratio = 1e5/hGenSpectraCorrPP->GetBinContent(x);
mPP_Matrix->SetBinContent(x,y,mPP_Matrix->GetBinContent(x,y)*ratio);
mPP_Matrix->SetBinError(x,y,mPP_Matrix->GetBinError(x,y)*ratio);
}
}
// mPbPb_Matrix[i]->Smooth(0);
// Ok major differences here between my code and Kurt in b-jet Tools under Unfold - lines 469 and above.
if(printDebug) cout<<"getting the response matrix"<<endl;
mPP_Response = (TH2F*)mPP_Matrix->Clone("mPP_Response");
TH1F *hProjPP = (TH1F*)mPP_Response->ProjectionY()->Clone("hProjPP");
for (int y=1;y<=mPP_Response->GetNbinsY();y++) {
double sum=0;
for (int x=1;x<=mPP_Response->GetNbinsX();x++) {
if (mPP_Response->GetBinContent(x,y)<=1*mPP_Response->GetBinError(x,y)) {
// in the above if statement, kurt has 1*error and my old has 0*error
mPP_Response->SetBinContent(x,y,0);
mPP_Response->SetBinError(x,y,0);
}
sum+=mPP_Response->GetBinContent(x,y);
}
for (int x=1;x<=mPP_Response->GetNbinsX();x++) {
if (sum==0) continue;
double ratio = 1;
//if(dPbPb_TrgComb[i]->GetBinContent(y)==0) ratio = 1e-100/sum;
// else ratio = dPbPb_TrgComb[i]->GetBinContent(y)/sum
ratio = 1./sum;
if (hProjPP->GetBinContent(y)==0) ratio = 1e-100/sum;
else ratio = hProjPP->GetBinContent(y)/sum;
mPP_Response->SetBinContent(x,y,mPP_Response->GetBinContent(x,y)*ratio);
mPP_Response->SetBinError(x,y,mPP_Response->GetBinError(x,y)*ratio);
}
}
if(printDebug) cout<<"getting the normalized response matrix"<<endl;
mPP_ResponseNorm = (TH2F*)mPP_Matrix->Clone("mPP_ResponseNorm");
for (int x=1;x<=mPP_ResponseNorm->GetNbinsX();x++) {
double sum=0;
for (int y=1;y<=mPP_ResponseNorm->GetNbinsY();y++) {
if (mPP_ResponseNorm->GetBinContent(x,y)<=1*mPP_ResponseNorm->GetBinError(x,y)) {
mPP_ResponseNorm->SetBinContent(x,y,0);
mPP_ResponseNorm->SetBinError(x,y,0);
}
sum+=mPP_ResponseNorm->GetBinContent(x,y);
}
for (int y=1;y<=mPP_ResponseNorm->GetNbinsY();y++) {
if (sum==0) continue;
double ratio = 1./sum;
mPP_ResponseNorm->SetBinContent(x,y,mPP_ResponseNorm->GetBinContent(x,y)*ratio);
mPP_ResponseNorm->SetBinError(x,y,mPP_ResponseNorm->GetBinError(x,y)*ratio);
//.........这里部分代码省略.........
示例2: main
int main() {
float min_logL1 = 5986.94;
float min_logL0 = 5987.16;
string filepath="FINAL_RESULT_AB.root_RESULT__RESULT";
filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_175622__RNSG_46692.4.root__RESULT__RESULT"; // final MCwS
filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_175622__RNSG_46692.4.root__RESULT__RESULT";
filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_141126__RNSG_97048.1.root__RESULT__RESULT";
// *************************************************************************************
setlumi(PlottingSetup::luminosity);
setessentialcut(PlottingSetup::essential); // this sets the essential cut; this one is used in the draw command so it is AUTOMATICALLY applied everywhere. IMPORTANT: Do NOT store weights here!
stringstream resultsummary;
// write_analysis_type(PlottingSetup::RestrictToMassPeak,PlottingSetup::DoBTag);
do_png(true);
do_pdf(true);
do_eps(false);
do_C(true);
do_root(false);
PlottingSetup::directoryname = "pValuePlot";
gROOT->SetStyle("Plain");
bool do_fat_line = false; // if you want to have HistLineWidth=1 and FuncWidth=1 as it was before instead of 2
setTDRStyle(do_fat_line);
gStyle->SetTextFont(42);
bool showList = true;
set_directory(PlottingSetup::directoryname); // Indicate the directory name where you'd like to save the output files in Setup.C
set_treename("events"); // you can set the treename here to be used; options are "events" (for reco) for "PFevents" (for particle flow)
TFile *f = new TFile(filepath.c_str());
if(f->IsZombie()) {
cout << "Seems to be a zombie. goodbye." << endl;
return -1;
}
RooWorkspace *wa = (RooWorkspace*)f->Get("transferSpace");
RooPlot *plot = (RooPlot*) wa->obj("frame_mlledge_109fde50");
// cout << plot << endl;
wa->Print("v");
TCanvas *can = new TCanvas("can","can");
cout << "Address of plot : " << plot << endl;
// plot->Draw();
float pVal_mllmin=35;
float pVal_mllmax=90;
int is_data=PlottingSetup::data;
vector < std::pair < float, float> > loglikelihoods;
string function="";
for(int i=0; i< plot->numItems();i++){
string name = plot->getObject(i)->GetName();
if (plot->getObject(i)->IsA()->InheritsFrom( "RooCurve" ))function=name;
}
RooCurve* curve = (RooCurve*) plot->findObject(function.c_str(),RooCurve::Class()) ;
if (!curve) {
dout << "RooPlot::residHist(" << plot->GetName() << ") cannot find curve" << endl ;
return 0 ;
}
int iMinimum=0;
float min=1e7;
for(int i=0;i<curve->GetN();i++) {
double x,y;
curve->GetPoint(i,x,y);
if(y<min & y>=0) {
min=y;
iMinimum=i;
}
}
double x,y;
curve->GetPoint(iMinimum,x,y);
cout << "Minimum is at " << x << " : " << y << endl;
loglikelihoods.push_back(make_pair(x,y+min_logL1));
//move right starting from the minimum
for(int i=iMinimum+1;i<curve->GetN();i++) {
float yold=y;
curve->GetPoint(i,x,y);
//if(abs((y-yold)/yold)>0.5) continue;
loglikelihoods.push_back(make_pair(x,y+min_logL1));
}
/*
//.........这里部分代码省略.........
示例3: PhotonIsolation
void PhotonIsolation()
{
const float chargeHadronIsoCut_LWP_EB = 1.3;
const float chargeHadronIsoCut_MWP_EB = 0.44;
const float chargeHadronIsoCut_TWP_EB = 0.2;
const float chargeHadronIsoCut_LWP_EE = 1.36;
const float chargeHadronIsoCut_MWP_EE = 0.82;
const float chargeHadronIsoCut_TWP_EE = 0.27;
const float HoverECut_LWP_EB = 0.06;
const float HoverECut_MWP_EB = 0.04;
const float HoverECut_TWP_EB = 0.027;
const float HoverECut_LWP_EE = 0.027;
const float HoverECut_MWP_EE = 0.021;
const float HoverECut_TWP_EE = 0.020;
const float SigIEtaIEtaCut_LWP_EB = 0.01031;
const float SigIEtaIEtaCut_MWP_EB = 0.01023;
const float SigIEtaIEtaCut_TWP_EB = 0.00994;
const float SigIEtaIEtaCut_LWP_EE = 0.0269;
const float SigIEtaIEtaCut_MWP_EE = 0.0259;
const float SigIEtaIEtaCut_TWP_EE = 0.0258;
const float beamSpotZ = 0.282329;
const float beamSpotSigmaZ = 43.131299;
const float NPU = 140.0;
const int n_density = 5;
float IsoEff_NoTiming[n_density] = {0.0};// for n_density pu density bins, ranging from 0 to 2.0, bin width 0.1
float Err_IsoEff_NoTiming[n_density] = {0.0};// for n_density pu density bins, ranging from 0 to 2.0, bin width 0.1
float IsoEff_Timing50_TrkVtx[n_density] = {0.0};
float IsoEff_Timing80_TrkVtx[n_density] = {0.0};
float Err_IsoEff_Timing80_TrkVtx[n_density] = {0.0};
float IsoEff_Timing120_TrkVtx[n_density] = {0.0};
float IsoEff_Timing50_TrkPho[n_density] = {0.0};
float IsoEff_Timing80_TrkPho[n_density] = {0.0};
float IsoEff_Timing120_TrkPho[n_density] = {0.0};
float Err_IsoEff_Timing120_TrkPho[n_density] = {0.0};
int N_Pho_Total_NoBin = 0;
int N_Pho_Total_NoBin200 = 0;
int N_Pho_Total[n_density] = {0};
int N_Pho_Total200[n_density] = {0};
int N_Pho_PassIso_NoTiming[n_density] = {0};
int N_Pho_PassIso_Timing50_TrkVtx[n_density] = {0};
int N_Pho_PassIso_Timing80_TrkVtx[n_density] = {0};
int N_Pho_PassIso_Timing120_TrkVtx[n_density] = {0};
int N_Pho_PassIso_Timing50_TrkPho[n_density] = {0};
int N_Pho_PassIso_Timing80_TrkPho[n_density] = {0};
int N_Pho_PassIso_Timing120_TrkPho[n_density] = {0};
TTree * tree = 0;
TFile *f =new TFile("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU140_Timing_Iso.root");
TDirectory * dir = (TDirectory*)f->Get("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU140_Timing_Iso.root:/ntuples");
dir->GetObject("RazorEvents",tree);
//cout<<tree->GetEntries()<<endl;
int NEntries = tree->GetEntries();
TTree * tree200 = 0;
TFile *f200 =new TFile("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU200_Timing_Iso.root");
TDirectory * dir200 = (TDirectory*)f200->Get("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU200_Timing_Iso.root:/ntuples");
dir200->GetObject("RazorEvents",tree200);
//cout<<tree->GetEntries()<<endl;
int NEntries200 = tree200->GetEntries();
Int_t nPhotons;
Float_t pvZ;
Float_t phoSigmaIetaIeta[700];
Float_t pho_superClusterEta[700];
Float_t pho_superClusterEnergy[700];
Float_t pho_HoverE[700];
Bool_t pho_passEleVeto[700];
Float_t pho_sumChargedHadronPt_NewPV_NoTiming[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing50_TrkVtx[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing80_TrkVtx[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing120_TrkVtx[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing50_TrkPho[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing80_TrkPho[700];
Float_t pho_sumChargedHadronPt_NewPV_Timing120_TrkPho[700];
tree->SetBranchAddress("nPhotons", &nPhotons);
tree->SetBranchAddress("pvZ_New", &pvZ);
tree->SetBranchAddress("phoSigmaIetaIeta", phoSigmaIetaIeta);
tree->SetBranchAddress("pho_superClusterEta", pho_superClusterEta);
tree->SetBranchAddress("pho_superClusterEnergy", pho_superClusterEnergy);
tree->SetBranchAddress("pho_HoverE", pho_HoverE);
tree->SetBranchAddress("pho_passEleVeto", pho_passEleVeto);
//.........这里部分代码省略.........
示例4: mumucross
int mumucross()
{
TCanvas *c1 = new TCanvas("c1","di-mu cross section");
TF1 *f1 = new TF1("f1",mumuxsdang,2,3.1,0);
f1->GetXaxis()->SetTitle("Ecm (GeV)");
f1->GetYaxis()->SetTitle("#sigma (nbar)");
f1->Draw();
TCanvas *c2 = new TCanvas("c2","BhaBha cross section");
TF1 *f2 = new TF1("f2",bhabhaxsdang,2,3.1,0);
f2->GetXaxis()->SetTitle("Ecm (GeV)");
f2->GetYaxis()->SetTitle("#sigma (nbar)");
f2->Draw();
//return 0;
double energys[22] = {
2.0000, 2.0500, 2.1000, 2.1500, 2.1750,
2.2000, 2.2324, 2.3094, 2.3864, 2.3960,
2.5000, 2.6444, 2.6464, 2.7000, 2.8000,
2.9000, 2.9500, 2.9810, 3.0000, 3.0200,
3.0800, 2.125
};
double luminosity[22] = {
10074, 3343, 12167, 2841, 10625,
13699, 11856, 21089, 22549, 66869,
1098, 33722, 34003, 1034, 1008,
105253, 15942, 16071, 15881, 17290,
126188, 108490
};
double mccross[22] = {
22.3763923, 21.2939032, 20.4018357, 19.4157221, 18.9842334,
18.5074456, 18.0584889, 16.9121088, 15.8285505, 15.6795155,
14.3957219, 12.9117968, 12.9257851, 12.3962790, 11.5661025,
10.7698318, 10.4156161, 10.2122681, 10.0576468, 9.9391192,
9.5904453, 21.1564523
};
double nevt[22] = {
5.0e5, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
5.0e5, 5.0e5, 5.0e5, 5.0e5, 1.07e6,
5.0e5, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
1.13e6, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
1.16e6, 1.93e6
};
TFile *file = new TFile("output/output_bck.root");
TH1D *hdata;
TH1D *hmcmu;
TH1D *hmcKK;
gStyle->SetOptStat(0);
for (int i=0 ;i<22;i++) {
//double nmumu = luminosity[i]*f1->Eval(energys[i]);
//cout<<nmumu<<"\t events at "<< energys[i] <<" GeV."<<endl;
//double nee = luminosity[i]*f2->Eval(energys[i]);
//cout<<nee<<"\t events at "<< energys[i] <<" GeV."<<endl;
cout<<"ene "<<energys[i]<<", scale is "<< luminosity[i]/(nevt[i]/mccross[i])<<endl;
//continue;
char hname[100];
// get data spectrum
//sprintf(hname,"hp_KK_%d",(int)(energys[i]*10000+0.5));
sprintf(hname,"hp_mcKK_%.4f",energys[i]);
hdata = (TH1D*)file->Get(hname);
hdata->SetFillColor(0);
//hdata->SetTitle("momentum spectrum of data and di-#mu MC");
hdata->SetTitle("");
hdata->GetXaxis()->SetTitle("p (GeV/c)");
hdata->GetYaxis()->SetTitle("Counts");
// get dimu spectrum
sprintf(hname,"hp_mcmumu_%.4f",energys[i]);
hmcmu = (TH1D*)file->Get(hname);
// get mc KK spectrum
sprintf(hname,"hp_mcKK2_%.4f",energys[i]);
hmcKK = (TH1D*)file->Get(hname);
int binid;
double maxdata = getMax(hdata,1,50,binid);
double maxmcKK = getMax(hmcKK,1,50,binid);
cout<< hdata <<"\t"<<hmcmu<<" ndata,nmcKK "<< maxdata<<" "<<maxmcKK<<endl;
sprintf(hname,"Canvas_%02d",i+1);
TCanvas *chp = new TCanvas(hname,"Background analysis");
chp->SetMargin(0.15,0.1,0.15,0.1);
hdata->SetLineColor(1);
hdata->SetTitle("");
hdata->GetXaxis()->SetNdivisions(505);
hdata->GetXaxis()->SetTitleSize(0.05);
hdata->GetXaxis()->SetLabelSize(0.05);
hdata->GetYaxis()->SetNdivisions(505);
hdata->GetYaxis()->SetTitleSize(0.05);
hdata->GetYaxis()->SetLabelSize(0.05);
hdata->GetYaxis()->SetTitleOffset(1.1);
hdata->Draw("E");
// scale dimu to data according luminosity
hmcmu->Scale(luminosity[i]/(nevt[i]/mccross[i]));
hmcmu->SetLineColor(2);
hmcmu->SetMarkerColor(2);
hmcmu->SetFillColor(0);
hmcmu->Draw("same");
//.........这里部分代码省略.........
示例5: rootelectron
//注意,主函数名称和文件名称要一样
int rootelectron(void)
{
const double PI = 3.1415926535897932384626433832795;
///////////////////////%%%%%%输入参数%%%%%%%%%%%/////////////////////////////////////////////////////////////
string filename="spectrum.root"; //数据文件名称
//root文件中tree的名称,需要和Geant4中的一致
string treename="selectree";
string branchname = "particledata"; //root文件中branch的名称,需要和Geant4中的一致
// string treename[4]={"pelectree","pgamatree","selectree","sgamatree"};
double Ebins=200; //统计能谱的直方图中bin的个数,
double Elow=-75.0; //能谱的低能截止,单位keV
double Eup=95.5555673; //能谱的高能能截止,单位keV
double cbins[3] = {180,180,180};//统计角分布的三个直方图中bin的个数【x,y,z】
double chlow[3] = {0,-90,0};//统计角分布的三个直方图的左边界【x,y,z】
double chup[3] = {180,90,90};//统计角分布的三个直方图的左边界【x,y,z】
double eneconv=1;//1000->MeV,1->keV将输入的能量转换单位
///////////////////////%%%%%%读入ROOT文件和Tree%%%%%%%%%%%/////////////////////////////////////////////////////////////
gROOT->Reset();
cout<<endl;
cout << "\nStart\n" <<endl;
TFile *myfile = new TFile(&filename[0]); //新建对象,指向数据文件
if (!myfile) {
cout <<"no input file!" <<endl;
return;
}
cout<<"Open file: "<<filename<<endl;
TTree *mytree = (TTree *)myfile->Get(&treename[0]);//新建tree,用于处理数据
if (!mytree) {
cout <<"no tree named: "<<treename<<endl;
return;
}
cout << "input tree: "<<treename <<endl;
//////////////////////////////%%设置新tree,用于数据处理%%%%%%/////////////////////////////////////%%%%%%%%%%%%%%%
double data[4]={0};//保存从root文件中取出的数据。
mytree->SetBranchAddress(&branchname[0],data);//设置新建的tree中branch的指向。
//create two histograms
double EEup=Eup/eneconv;
double EElow=Elow/eneconv;
double sintheta=0;
TH1D * h1 = new TH1D("E","Energy Spectrum of source",Ebins,EElow,EEup); //统计能谱
// TH2D *h2Ex = new TH2D("h2Ex","h2 Energy vs x",cbins[0],chlow[0],chup[0],Ebins,EElow,EEup);//统计角分布,出射方向与X方向的夹角
// TH2D *h2Ey = new TH2D("h2Ey","h2 Energy vs y",cbins[1],chlow[1],chup[1],Ebins,EElow,EEup);//统计角分布,出射方向与Y方向的夹角
TH2D *h2Ez = new TH2D("h2Ez","h2 Energy vs z",cbins[2],chlow[2],chup[2],Ebins,EElow,EEup);//统计角分布,出射方向与Z方向的夹角(横坐标名,纵坐标名,横坐标bins,横坐标低截止,横坐标高截止,纵坐标bins,纵坐标低截止,纵坐标高截止)
//read all entries and fill the histograms
Long64_t nentries = mytree->GetEntries(); //注意:因为事例entry很多,所以需要Long64_t
for (Long64_t i=0;i<nentries;i++) {
mytree->GetEntry(i);
//将处理数据的tree指向数据文件中的第i个事例,由于上面设置了的tree中branch的指针为变量data,所以处理data变量就处理了数据。
h1->Fill(data[3]); //统计能谱
// sintheta= -sqrt(1-data[3]*data[3]);//data[3]=cos(theta)
// h2Ex->Fill(acos(data[1]/sintheta)/PI*180,data[0]);
// h2Ey->Fill(asin(data[2]/sintheta)/PI*180,data[0]);
h2Ez->Fill(acos(data[1])/PI*180,data[0]);
}
///////////////////////%%%%%%%%%Profile Histogram%%%%%%%/////////////////////////////////////
TH1D* az;
// ax=h2Ex->ProjectionX(); //将二维散点图投影到1维就得到了一维的角分布,
// ay=h2Ey->ProjectionX(); //将二维散点图投影到1维就得到了一维的角分布,
az=h2Ez->ProjectionX(); //将二维散点图投影到1维就得到了一维的角分布,
///////////////////////%%%%%%%%%Plot Histogram%%%%%%%/////////////////////////////////////
//gROOT->SetStyle("Plain"); // uncomment to set a different style
gStyle->SetPalette(1); // use precomputed color palette 1
gStyle->SetOptStat("ner");
/* TCanvas *c1 = new TCanvas();
c1->cd();
h2Ex->Draw("C"); //energy
h2Ex->SetTitle("Energy Spectrum"); //散点图
h2Ex->GetXaxis()->SetTitle("alpha(x)");
h2Ex->GetYaxis()->SetTitle("Energy (keV)");
h2Ex->GetXaxis()->CenterTitle();
h2Ex->GetYaxis()->CenterTitle();
gPad->Update();
// eneh->GetXaxis()->SetRangeUser(0, 1200);
// eneh->GetYaxis()->SetRangeUser(1e-8,1e-1);
// gPad->SetLogx();
// gPad->SetLogy();
TCanvas *c2 = new TCanvas();
c2->cd();
h2Ey->Draw("C"); //energy
h2Ey->SetTitle("Energy Spectrum");
h2Ey->GetXaxis()->SetTitle("beta(y)");
h2Ey->GetYaxis()->SetTitle("Energy (MeV)");
h2Ey->GetXaxis()->CenterTitle();
h2Ey->GetYaxis()->CenterTitle();
gPad->Update(); */
/*
TCanvas *c3 = new TCanvas();
c3->cd();
h2Ez->Draw("C"); //energy
h2Ez->SetTitle("Energy Spectrum");
h2Ez->GetXaxis()->SetTitle("gamma(z)");
//.........这里部分代码省略.........
示例6: hackedbob
void hackedbob(bool eff = false)
{
using namespace std;
double integral, preInt = 1.0;
vector<pair<string,string> > fnames;
vector<pair<string,string> > hnames;
//TFile * tfile = new TFile("/home/ugrad/pastika/cms/HeavyNu/CMSSW_4_2_7/src/HeavyNu/AnalysisModules/analysis.root");
TFile * tfile = new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_DYToLL_M-50_7TeV-sherpa_nov3.root");
hnames.push_back(pair<string,string>("none","hNu%s/cut0_none/mWR"));
hnames.push_back(pair<string,string>("LLJJ Pt","hNu%s/cut1_LLJJpt/mWR"));
hnames.push_back(pair<string,string>("trig","hNu%s/cut2_TrigMatches/mWR"));
hnames.push_back(pair<string,string>("vertex","hNu%s/cut3_Vertex/mWR"));
hnames.push_back(pair<string,string>("mu1 pt","hNu%s/cut4_Mu1HighPt/mWR"));
hnames.push_back(pair<string,string>("Mll","hNu%s/cut5_diLmass/mWR"));
hnames.push_back(pair<string,string>("MWR","hNu%s/cut6_mWRmass/mWR"));
fnames.push_back(pair<string,string>("nom Mu24", "Mu24"));
fnames.push_back(pair<string,string>("nom Mu40", "Mu40"));
//fnames.push_back(pair<string,string>("JES Hi Mu24", "Mu24jesHi"));
//fnames.push_back(pair<string,string>("JES Lo Mu24", "Mu24jesLo"));
//fnames.push_back(pair<string,string>("JES Hi Mu40", "Mu40jesHi"));
//fnames.push_back(pair<string,string>("JES Lo Mu40", "Mu40jesLo"));
//fnames.push_back(pair<string,string>("MES Hi Mu24", "Mu24mesHi"));
//fnames.push_back(pair<string,string>("MES Lo Mu24", "Mu24mesLo"));
//fnames.push_back(pair<string,string>("MES Hi Mu40", "Mu40mesHi"));
//fnames.push_back(pair<string,string>("MES Lo Mu40", "Mu40mesLo"));
//fnames.push_back(pair<string,string>("MuID Hi Mu24", "Mu24midHi"));
//fnames.push_back(pair<string,string>("MuID Lo Mu24", "Mu24midLo"));
//fnames.push_back(pair<string,string>("MuID Hi Mu40", "Mu40midHi"));
//fnames.push_back(pair<string,string>("MuID Lo Mu40", "Mu40midLo"));
//fnames.push_back(pair<string,string>("Trig Hi Mu24", "Mu24trigHi"));
//fnames.push_back(pair<string,string>("Trig Lo Mu24", "Mu24trigLo"));
//fnames.push_back(pair<string,string>("Trig Hi Mu40", "Mu40trigHi"));
//fnames.push_back(pair<string,string>("Trig Lo Mu40", "Mu40trigLo"));
fnames.push_back(pair<string,string>("Pileup Hi Mu24", "Mu24puHi"));
fnames.push_back(pair<string,string>("Pileup Lo Mu24", "Mu24puLo"));
fnames.push_back(pair<string,string>("Pileup Hi Mu40", "Mu40puHi"));
fnames.push_back(pair<string,string>("Pileup Lo Mu40", "Mu40puLo"));
printf("%8s", "cut");
for(vector<pair<string,string> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
{
printf(eff?" & %15s ":" & %15s", i->first.c_str());
}
printf(" \\\\ \\hline\n");
for(vector<pair<string,string> >::const_iterator j = hnames.begin(); j != hnames.end(); j++)
{
printf("%8s", j->first.c_str());
for(vector<pair<string,string> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
{
char histname[256];
sprintf(histname, j->second.c_str(), i->second.c_str());
TH1* h = (TH1*)tfile->Get(histname);
integral = h->Integral(0, h->GetNbinsX()+1);
if(integral > 10) printf((j == hnames.begin() && eff)?" & %15.0f ":" & %15.0f", integral);
else if(integral > 1) printf((j == hnames.begin() && eff)?" & %15.2f ":" & %15.2f", integral);
else printf((j == hnames.begin() && eff)?" & %15.3f ":" & %15.3f", integral);
//if(j != hnames.begin() && eff)
//{
// TH1* hpre = (TH1*)tfile->Get((j-1)->second.c_str());
// preInt = hpre->Integral(0, h->GetNbinsX()+1);
// printf(" (%4.2f)", integral/preInt);
//}
}
printf(" \\\\ \\hline\n");
}
}
示例7: plotMVAOutput
void plotMVAOutput( bool printgif = false ){
//gROOT->ProcessLine(".x selection.h");
//char* path = "Trainings/H130_WWTo2L2Nu/output";
//char* path = "Trainings/H130_WWTo2L2Nu_WJetsToLNu/output";
char* path = "Trainings/H130_allbkg_4vars/output";
//char* mvaname = "MVA_PDERS";
//char* mvaname = "MVA_MLPBNN";
//char* mvaname = "MVA_BDT";
//char* mvaname = "LikelihoodPCA";
vector<char*> mvanames;
mvanames.push_back("BDT");
mvanames.push_back("MLPBNN");
const unsigned int nmva = mvanames.size();
int rebin = 10;
int colors[] = { 5 , 2 , 4 , 3 , 7 , 8 , 6 , 9 , 10};
vector<char*> samples;
samples.push_back("WWTo2L2Nu");
samples.push_back("GluGluToWWTo4L");
samples.push_back("WZ");
samples.push_back("ZZ");
samples.push_back("TTJets");
samples.push_back("tW");
samples.push_back("WJetsToLNu");
samples.push_back("DY");
const unsigned int nsamples = samples.size();
char* higgssample = "Higgs130";
//char* higgssample = "Higgs160";
//char* higgssample = "Higgs200";
TCanvas *can[nmva];
for( unsigned int imva = 0 ; imva < nmva ; ++imva ){
TFile* file = new TFile();
TH1F* hist = new TH1F();
TH1F* bkghist = new TH1F();
THStack* bkgstack = new THStack("bkgstack","bkgstack");
TLegend *leg = new TLegend(0.3,0.7,0.5,0.9);
leg->SetBorderSize(1);
leg->SetFillColor(0);
//loop over backgrounds
for( unsigned int i = 0 ; i < nsamples ; ++i ){
file = TFile::Open(Form("%s/%s.root",path,samples.at(i)));
hist = cloneHist( (TH1F*) file->Get( Form("MVA_%s",mvanames.at(imva) ) ) );
hist->SetFillColor(colors[i]);
leg->AddEntry(hist,samples.at(i),"f");
if( i == 0 ) bkghist = (TH1F*) hist->Clone();
else bkghist -> Add(hist);
hist->Rebin( rebin );
bkgstack->Add(hist);
}
//higgs sample
file = TFile::Open(Form("%s/%s.root",path,higgssample));
TH1F* higgshist = cloneHist( (TH1F*) file->Get( Form("MVA_%s",mvanames.at(imva) ) ) );
higgshist->SetLineWidth(2);
leg->AddEntry(higgshist,higgssample,"l");
float bkg = 0;
float sig = 0;
float minbkg = 1.48;
//float minbkg = 1.10;
float cut = 0.;
for( int ibin = 1 ; ibin < bkghist->GetNbinsX() ; ibin++ ){
bkg = bkghist->Integral( ibin , 10000 );
sig = higgshist->Integral( ibin , 10000 );
if( bkg < minbkg ){
cut = bkghist->GetBinCenter(ibin);
cout << endl;
cout << "S/B " << sig/bkg << endl;
cout << "Sig " << sig << endl;
cout << "Bkg " << bkg << endl;
cout << "cut value " << cut << endl;
break;
}
}
float cutsig = sig;
//.........这里部分代码省略.........
示例8: ComputeTestStat
float ComputeTestStat(TString wsfile, double mu_susy_sig_val) {
gROOT->Reset();
TFile* wstf = new TFile( wsfile ) ;
RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
ws->Print() ;
ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
modelConfig->Print() ;
RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;
rds->Print() ;
rds->printMultiline(cout, 1, kTRUE, "") ;
RooAbsPdf* likelihood = modelConfig->GetPdf() ;
RooRealVar* rrv_mu_susy_sig = ws->var("mu_susy_all0lep") ;
if ( rrv_mu_susy_sig == 0x0 ) {
printf("\n\n\n *** can't find mu_susy_all0lep in workspace. Quitting.\n\n\n") ;
return ;
} else {
printf(" current value is : %8.3f\n", rrv_mu_susy_sig->getVal() ) ; cout << flush ;
rrv_mu_susy_sig->setConstant(kFALSE) ;
}
/*
// check the impact of varying the qcd normalization:
RooRealVar *rrv_qcd_0lepLDP_ratioH1 = ws->var("qcd_0lepLDP_ratio_H1");
RooRealVar *rrv_qcd_0lepLDP_ratioH2 = ws->var("qcd_0lepLDP_ratio_H2");
RooRealVar *rrv_qcd_0lepLDP_ratioH3 = ws->var("qcd_0lepLDP_ratio_H3");
rrv_qcd_0lepLDP_ratioH1->setVal(0.3);
rrv_qcd_0lepLDP_ratioH2->setVal(0.3);
rrv_qcd_0lepLDP_ratioH3->setVal(0.3);
rrv_qcd_0lepLDP_ratioH1->setConstant(kTRUE);
rrv_qcd_0lepLDP_ratioH2->setConstant(kTRUE);
rrv_qcd_0lepLDP_ratioH3->setConstant(kTRUE);
*/
printf("\n\n\n ===== Doing a fit with SUSY component floating ====================\n\n") ;
RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
double logLikelihoodSusyFloat = fitResult->minNll() ;
double logLikelihoodSusyFixed(0.) ;
double testStatVal(-1.) ;
if ( mu_susy_sig_val >= 0. ) {
printf("\n\n\n ===== Doing a fit with SUSY fixed ====================\n\n") ;
printf(" fixing mu_susy_sig to %8.2f.\n", mu_susy_sig_val ) ;
rrv_mu_susy_sig->setVal( mu_susy_sig_val ) ;
rrv_mu_susy_sig->setConstant(kTRUE) ;
fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
logLikelihoodSusyFixed = fitResult->minNll() ;
testStatVal = 2.*(logLikelihoodSusyFixed - logLikelihoodSusyFloat) ;
printf("\n\n\n ======= test statistic : -2 * ln (L_fixed / ln L_max) = %8.3f\n\n\n", testStatVal ) ;
}
return testStatVal ;
}
示例9: makeStackedPlot_Jet1pT
void makeStackedPlot_Jet1pT(){
gROOT->SetStyle("Plain");
gStyle->SetErrorX(0);
TCanvas c1("c1","Stacked Histogram",10,10,700,800);
TPad *p_2=new TPad("p_2", "p_2", 0, 0, 1, 0.35);
TPad *p_1=new TPad("p_1", "p_1", 0, 0.35, 1, 1);
p_1->SetBottomMargin(0.005);
p_1->SetFillStyle(4000);
p_1->SetFrameFillColor(0);
p_2->SetFillStyle(4000);
p_2->SetFrameFillColor(0);
p_1->SetLogy();
p_1->Draw();
p_2->Draw();
p_1->cd();
double qcd_100_200 = 27540000.0*209.0/80093092.0;
double qcd_200_300 = 1735000.0*209.0/18717349.0;
double qcd_300_500 = 366800.0*209.0/20086103.0;
double qcd_500_700 = 29370.0*209.0/19542847.0;
double qcd_700_1000 = 6524.0*209.0/15011016.0;
double qcd_1000_1500 = 1064.0*209.0/4963895.0;
double qcd_1500_2000 = 121.5*209.0/3848411.0;
double qcd_2000_inf = 25.4*209.0/1961774.0;
THStack hs("hs","Jet pT > 50 GeV");
TFile* file0 = TFile::Open("output_QCD_HT_100_200.root");
TH1F *h0 = (TH1F*)file0->Get("general/jet_pT50");
h0->Rebin(1);
h0->GetXaxis()->SetRangeUser(0, 300);
h0->Scale(qcd_100_200);
h0->SetFillColor(kYellow);
h0->Draw();
hs.Add(h0);
TFile* file1 = TFile::Open("output_QCD_HT_200_300.root");
TH1F *h1 = (TH1F*)file1->Get("general/jet_pT50");
h1->Rebin(1);
h1->GetXaxis()->SetRangeUser(0, 300);
h1->Scale(qcd_200_300);
h1->SetFillColor(kYellow);
h1->Draw();
hs.Add(h1);
TFile* file2 = TFile::Open("output_QCD_HT_300_500.root");
TH1F *h2 = (TH1F*)file2->Get("general/jet_pT50");
h2->Rebin(1);
h2->GetXaxis()->SetRangeUser(0, 300);
h2->Scale(qcd_300_500);
h2->SetFillColor(kYellow);
h2->Draw();
hs.Add(h2);
TFile* file3 = TFile::Open("output_QCD_HT_500_700.root");
TH1F *h3 = (TH1F*)file3->Get("general/jet_pT50");
h3->Rebin(1);
h3->GetXaxis()->SetRangeUser(0, 300);
h3->Scale(qcd_500_700);
h3->SetFillColor(kYellow);
h3->Draw();
hs.Add(h3);
TFile* file4 = TFile::Open("output_QCD_HT_700_1000.root");
TH1F *h4 = (TH1F*)file4->Get("general/jet_pT50");
h4->Rebin(1);
h4->GetXaxis()->SetRangeUser(0, 300);
h4->Scale(qcd_700_1000);
h4->SetFillColor(kYellow);
h4->Draw();
hs.Add(h4);
TFile* file5 = TFile::Open("output_QCD_HT_1000_1500.root");
TH1F *h5 = (TH1F*)file5->Get("general/jet_pT50");
h5->Rebin(1);
h5->GetXaxis()->SetRangeUser(0, 300);
h5->Scale(qcd_1000_1500);
h5->SetFillColor(kYellow);
h5->Draw();
hs.Add(h5);
TFile* file6 = TFile::Open("output_QCD_HT_1500_2000.root");
TH1F *h6 = (TH1F*)file6->Get("general/jet_pT50");
h6->Rebin(1);
h6->GetXaxis()->SetRangeUser(0, 300);
h6->Scale(qcd_1500_2000);
h6->SetFillColor(kYellow);
h6->Draw();
hs.Add(h6);
TFile* file7 = TFile::Open("output_QCD_HT_2000_inf.root");
TH1F *h7 = (TH1F*)file7->Get("general/jet_pT50");
h7->Rebin(1);
h7->GetXaxis()->SetRangeUser(0, 300);
h7->Scale(qcd_2000_inf);
h7->SetFillColor(kYellow);
h7->Draw();
hs.Add(h7);
//.........这里部分代码省略.........
示例10: Plot_AM_events_07Sep_Susy_1_auto
void Plot_AM_events_07Sep_Susy_1_auto() {
TString cutNameBefore = Form("Data/");
// cut_variable
TString cutNameAfter = Form("_1_6_Tot_temp");
gROOT->LoadMacro("PlotVHqqHggH.C+");
gInterpreter->ExecuteMacro("LatinoStyle2.C");
TCanvas* c1 = new TCanvas("events","events",500,600);
TFile* f = new TFile("~/Cern/Code/VBF/qqHWW/AnalysisPackage_qqHWWlnulnu/out_test_Latinos_07Sep2012_2000_Run2012AB_8TeV_SUSY.root");
PlotVHqqHggH* hs = new PlotVHqqHggH();
hs->setLumi(5.063);
hs->setLabel("event");
hs->addLabel(" #sqrt{s} = 8 TeV");
TString name;
std::vector<int> vectColourBkg;
std::vector<double> vectSystBkg;
std::vector<std::string> vectNameBkg;
std::vector<TH1F*> vectTHBkg;
std::vector<int> vectColourSig;
std::vector<double> vectSystSig;
std::vector<std::string> vectNameSig;
std::vector<TH1F*> vectTHSig;
///==== signal (begin) ====
name = Form("%sT2tt-350-70%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop350-lsp70");
vectColourSig.push_back(6);
name = Form("%sT2tt-500-70%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop500-lsp70");
vectColourSig.push_back(97);
name = Form("%sT2tt-350-100%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop350-lsp100");
vectColourSig.push_back(70);
name = Form("%sT2tt-500-100%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop500-lsp100");
vectColourSig.push_back(65);
name = Form("%sT2tt-500-200%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop500-lsp200");
vectColourSig.push_back(5);
name = Form("%sT2tt-200-0%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHSig.push_back ( (TH1F*) f->Get(name) );
vectNameSig.push_back ("T2tt-stop200-lsp0");
vectColourSig.push_back(7);
///==== signal (end) ====
name = Form("%sDATA%s",cutNameBefore.Data(),cutNameAfter.Data());
hs->setDataHist ((TH1F*)f->Get(name));
hs->setBlindBinSx(0);
hs->setBlindBinDx(0);
hs->setCutSx(-999,">");
hs->setCutDx(-999,"<");
///==== background (begin) ====
name = Form("%sH-125%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHBkg.push_back ( (TH1F*) f->Get(name) );
vectNameBkg.push_back ("H-125");
vectColourBkg.push_back(633);
vectSystBkg.push_back(0.00);
name = Form("%sVV%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHBkg.push_back ( (TH1F*) f->Get(name) );
vectNameBkg.push_back ("WZ/ZZ");
vectColourBkg.push_back(858);
vectSystBkg.push_back(0.00);
name = Form("%sWJets%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHBkg.push_back ( (TH1F*) f->Get(name) );
vectNameBkg.push_back ("W+jets");
vectColourBkg.push_back(921);
vectSystBkg.push_back(0.36);
name = Form("%sTop%s",cutNameBefore.Data(),cutNameAfter.Data());
vectTHBkg.push_back ( (TH1F*) f->Get(name) );
//.........这里部分代码省略.........
示例11: DrawStatUncertainty
void DrawStatUncertainty()
{
TFile *f;
TString ss("ExpStat_PYTHIA_canvasOutFile.root");
f = new TFile(ss);
TCanvas *c
c = (TCanvas*)f->Get("_R92777");
TString algoTitle;
algoTitle = "k_{T} D = 0.6";
TList *li = (TList*)gPad->GetListOfPrimitives();
TObject *obj;
TIter next(li);
TH1D *histList[100], *htmp;
int i,j,N;
while ((obj = (TObject*)next()))
{
TString cname = obj->ClassName();
TString name = obj->GetName();
cout << cname <<" "<<name<<endl;
if (cname=="TH1D")
{
histList[N] = (TH1D*)gPad->FindObject(obj);
histList[N]->SetFillColor(0);
histList[N]->SetFillStyle(0);
N++;
}
}
TCanvas *c = new TCanvas("c","c");
gPad->SetLogx();
gPad->SetLogy();
//histList[0]->SetMaximum(1e+10);
histList[0]->SetTitle("");
histList[0]->GetXaxis()->SetTitle("jet p_{T} (GeV)");
histList[0]->GetXaxis()->SetTitleFont(42);
histList[0]->GetXaxis()->SetLabelFont(42);
histList[0]->GetXaxis()->SetTitleSize(0.05);
histList[0]->GetXaxis()->SetLabelSize(0.05);
histList[0]->GetYaxis()->SetTitle("Relative Statistical Uncertainty");
histList[0]->GetYaxis()->SetTitleFont(42);
histList[0]->GetYaxis()->SetLabelFont(42);
histList[0]->GetYaxis()->SetTitleSize(0.05);
histList[0]->GetYaxis()->SetLabelSize(0.05);
histList[0]->GetYaxis()->SetNdivisions(505);
histList[0]->SetLineWidth(1);
histList[0]->SetLineStyle(1);
histList[0]->SetLineColor(1);
histList[0]->SetMarkerColor(1);
histList[0]->SetMarkerStyle(20);
histList[0]->SetMarkerSize(1.2);
histList[1]->SetLineColor(2);
histList[1]->SetMarkerColor(2);
histList[1]->SetLineWidth(1);
histList[1]->SetMarkerStyle(25);
histList[1]->SetLineStyle(1);
histList[1]->SetMarkerSize(1.2);
histList[2]->SetLineWidth(1);
histList[2]->SetLineStyle(1);
histList[2]->SetLineColor(4);
histList[2]->SetMarkerColor(4);
histList[2]->SetMarkerStyle(22);
histList[2]->SetMarkerSize(1.2);
histList[0]->Draw("P");
histList[1]->Draw("sameP");
histList[2]->Draw("sameP");
TLegend *leg = new TLegend(0.47,0.2,0.92,0.35);
leg->SetTextFont(42);
leg->AddEntry(histList[0],"|y| < 0.55","P");
leg->AddEntry(histList[1],"1.10 < |y| < 1.70","P");
leg->AddEntry(histList[2],"1.70 < |y| < 2.50","P");
leg->SetFillColor(0);
leg->SetBorderSize(0);
leg->Draw();
TPaveText *pave3 = new TPaveText(0.2,0.65,0.45,0.9,"NDC");
pave3->SetTextFont(42);
pave3->AddText("CMS preliminary");
pave3->AddText(algoTitle);
pave3->AddText("#sqrt{s} = 10 TeV");
pave3->AddText("");
pave3->AddText("#int L dt = 10 pb^{-1}");
pave3->SetFillColor(0);
pave3->SetBorderSize(0);
pave3->Draw();
c->Print("KT6_StatUncertainty.eps");
}
示例12: Ana
int Ana(const char *filename, const char* outdir, TFile *fileout)
{
//TFile *file = new TFile("mc_KPI_22324.root");
TFile *file;
if (filename==0) file= new TFile("KK_22324.root");
else file = new TFile(filename);
std::cout<<"File name is "<<file->GetName()<<std::endl;
if (file==0) return -1;
TTree *tree = (TTree*)file->Get("TwoProng");
if (tree==0) return -2;
double ene[22];
double pcut[22];
double epcut[22];
double thecut[22];
double m_pcut;
double m_epcut;
double m_thecut;
ifstream cutin("cutpar");
char line[1000];
if (cutin.is_open()) cutin.getline(line,1000);
int iene=0;
// set E/p and p cut, mark: set cut
while (!cutin.eof()){
cutin.getline(line,1000);
istringstream iss;
iss.str(line);
iss>>ene[iene]>>pcut[iene]>>epcut[iene];
iene++;
if (iene==22) break;
}
for (int i=0;i<22;i++) epcut[i] = 1.64295 - 0.629622*ene[i] + 0.104755 *pow(ene[i],2);
//for (int i=0;i<21;i++) thecut[i] = 173.946 + 1.74736*ene[i];
//for (int i=0;i<22;i++) thecut[i] = 179;
double kappx,kappy,kappz,kampx,kampy,kampz;
int nneu;
int run;
int idxmc;
int pdgid[100];
int motheridx[100];
int emcstatusInt;
short emcstatusShort;
double emctrk1;
double emctrk2;
double epratio1;
double epratio2;
int ntof1;
int ntof2;
int tofl1[5];
int tofl2[5];
double tof1[5];
double tof2[5];
tree->SetBranchAddress("run",&run);
tree->SetBranchAddress("indexmc",&idxmc);
tree->SetBranchAddress("pdgid",pdgid);
tree->SetBranchAddress("motheridx",motheridx);
tree->SetBranchAddress("kappx",&kappx);
tree->SetBranchAddress("kappy",&kappy);
tree->SetBranchAddress("kappz",&kappz);
tree->SetBranchAddress("kampx",&kampx);
tree->SetBranchAddress("kampy",&kampy);
tree->SetBranchAddress("kampz",&kampz);
//tree->SetBranchAddress("costheta",&costheta);
tree->SetBranchAddress("nneu",&nneu);
if (strncmp(tree->GetBranch("emcstatus")->GetLeaf("emcstatus")->GetTypeName(),"Short_t",7)==0)
tree->SetBranchAddress("emcstatus",&emcstatusShort);
else
tree->SetBranchAddress("emcstatus",&emcstatusInt);
tree->SetBranchAddress("epratio1",&epratio1);
tree->SetBranchAddress("epratio2",&epratio2);
tree->SetBranchAddress("emctrk1",&emctrk1);
tree->SetBranchAddress("emctrk2",&emctrk2);
tree->SetBranchAddress("ntof1",&ntof1);
tree->SetBranchAddress("toflayer1",tofl1);
tree->SetBranchAddress("tof1",tof1);
tree->SetBranchAddress("ntof2",&ntof2);
tree->SetBranchAddress("toflayer2",tofl2);
tree->SetBranchAddress("tof2",tof2);
double pi = TMath::Pi();
double mka = 0.493677;
//double mka = 0.13957;
//double mka = 0.1057;
tree->GetEntry(0);
double Ebeam = GetEnergy(run);
if (Ebeam<1.0) Ebeam = getEne(filename);
//Ebeam = 3.08;
for (int iene=0;iene<22;iene++){
if (fabs(Ebeam-ene[iene])<1e-4) {
m_pcut = pcut[iene];
m_epcut = epcut[iene];
m_thecut = thecut[iene];
break;
}
}
//m_epcut = m_epcut +0.05; // change E/p to determine uncertainty from this cut
//m_thecut = 179+0.2; // uncertainty from theta cut
//m_pcut = 0.956409;
//.........这里部分代码省略.........
示例13: Draweff
void Draweff(){
int sth=0, Gth=0;
TFile *f = TFile::Open(outG);
if(sth==0){TString dirname = "std";}
else if(sth==1){TString dirname ="Gri055";}
else {TString dirname ="Gri101";}
gStyle->SetErrorX(0);
TString name;
TObjString* dataname = (TObjString*)f->Get(Form("dataname"));
TObjString* histoname = (TObjString*)f->Get(Form("histoname"));
if(Gth==0)
name = "G0";
else if(Gth<nGlau)
name = Form("Glau_%d",Gth);
else
name = Form("bin_%d",Gth-nGlau+1);
TObjString* Glaubername = (TObjString*)f->Get(Form("%s/%s/Glaubername",dirname.Data(),name.Data()));
TVectorD* xmin = (TVectorD*)f->Get(Form("%s/%s/xmin",dirname.Data(),name.Data()));
TVectorD* xmax = (TVectorD*)f->Get(Form("%s/%s/xmax",dirname.Data(),name.Data()));
TVectorD* mubest = (TVectorD*)f->Get(Form("%s/%s/mubest",dirname.Data(),name.Data()));
TVectorD* kbest = (TVectorD*)f->Get(Form("%s/%s/kbest",dirname.Data(),name.Data()));
TVectorD* Ndf = (TVectorD*)f->Get(Form("%s/%s/Ndf",dirname.Data(),name.Data()));
TVectorD* chis = (TVectorD*)f->Get(Form("%s/%s/chis",dirname.Data(),name.Data()));
TVectorD *kpoint = (TVectorD*)f->Get(Form("%s/%s/kpoint",dirname.Data(),name.Data()));
TFile *fdata = TFile::Open(dataname->GetString());
TH1D *histo_obs = (TH1D*)fdata->Get(histoname->GetString());
histo_obs->Sumw2();
TFile *fGlauber = TFile::Open(Glaubername->GetString());
int binnum = histo_obs->GetNbinsX();
double Minx = histo_obs->GetXaxis()->GetXmin();
double Maxx = histo_obs->GetXaxis()->GetXmax();
double binsize = (Double_t)(Maxx-Minx)/binnum;
int xbinmin=(int)(((*xmin)[0]-Minx)/binsize);
int xbinmax=(int)(((*xmax)[0]-Minx)/binsize);
TH1D *histo_exp = new TH1D("histo_exp","Simulated distribution;Multiplicity;Event Fraction",binnum,Minx,Maxx);
histo_exp->Sumw2();
Int_t ibin;
TH1D *histo_obs_norm = (TH1D*)histo_obs->Clone();
histo_obs_norm->Scale(1/histo_obs->Integral(xbinmin,xbinmax));
TF1 *NBD_fun = new
TF1("NBD_fun","[0]*TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1]))*TMath::Power([2]/[1],x)/TMath::Power([2]/[1]+1,x+[1])",0,100);
NBD_fun->SetParameter(0,1); //[0]: Normalized constant
NBD_fun->SetParameter(1,(*kbest)[0]); //[1]: k value
NBD_fun->SetParameter(2,(*mubest)[0]); //[2]: mu value
TTree *t = (TTree*) fGlauber->Get("nt_Pb_Pb");
Long_t Nevent;
Nevent = (Long_t) t->GetEntries();
Long_t Ev; Int_t Bino; Double_t Para, Bi_Para, Mult;
Float_t Ncoll;
t->SetBranchAddress("Ncoll",&Ncoll);
for(Ev=0; Ev<Nevent; Ev++){
if(Ev%100000==0) cout<<"Have run "<<Ev<<" events"<<endl;
t->GetEntry(Ev);
Para = 0; //make sure that Para doesn't accumulate through loops
for(Bino=0; Bino<Ncoll; Bino++){
Bi_Para = NBD_fun->GetRandom();
Para += Bi_Para;
}
histo_exp->Fill(Para);
}
Double_t SumEvent, scale;
SumEvent = histo_exp->Integral(xbinmin,xbinmax);
scale = 1/SumEvent;
TH1D *histo_exp_norm = (TH1D*) histo_exp->Clone();
histo_exp_norm->Scale(scale);
TCanvas *c1 = new TCanvas();
gStyle->SetOptStat(kFALSE);
double hfbin[]={0,1,2,3,4,6,8,10,13,16,20,25,30,40,55,70,90};
int nhfbin = 16;
rehisto_obs_norm = (TH1D*)histo_obs_norm->Rebin(nhfbin,"rehisto_obs_norm",hfbin);
normalizeByBinWidth(rehisto_obs_norm);
rehisto_exp_norm = (TH1D*)histo_exp_norm->Rebin(nhfbin,"rehisto_exp_norm",hfbin);
normalizeByBinWidth(rehisto_exp_norm);
TH1D* ratio = (TH1D*)rehisto_obs_norm->Clone("ratio");
ratio->Divide(rehisto_exp_norm);
ratio->SetMaximum(1.2);
ratio->SetMinimum(0);
ratio->GetXaxis()->SetTitle("HF #Sigma E_{T}");
ratio->GetYaxis()->SetTitle("ratio");
/*
示例14: make_rValues
void make_rValues(std::string model, int m1, int m2){
//setup file to hold r-values if it doesn't already exist
TFile *f_temp = new TFile(Form("r-values_%s.root", model.c_str()), "NEW");
if(f_temp){
// default case: T5ZZ binning
int m1_max = 2125;
int m1_min = 775;
int m2_max = 2025;
int m2_min = 75;
int m1_div = 50;
int m2_div = 50;
if (model.find("TChiWZ") != std::string::npos) {
m1_max = 712.5;
m1_min = 87.5;
m2_max = 305;
m2_min = -5;
m1_div = 25;
m2_div = 10;
}
int nbinsx = (m1_max - m1_min)/m1_div;
int nbinsy = (m2_max - m2_min)/m2_div;
TH2F *hExp = new TH2F("hExp", "hExp" , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
TH2F *hObs = new TH2F("hObs", "hObs" , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
TH2F *hExp1m = new TH2F("hExp1m", "hExp1m", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
TH2F *hExp2m = new TH2F("hExp2m", "hExp2m", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
TH2F *hExp1p = new TH2F("hExp1p", "hExp1p", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
TH2F *hExp2p = new TH2F("hExp2p", "hExp2p", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
//TH2F *hSig = new TH2F("hSig", "hSig" , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
f_temp->Write();
f_temp->Close();
}
delete f_temp;
//This file is created earlier by running combine
TFile *limit_file = new TFile(Form("limit_%s_%d_%d.root", model.c_str(), m1, m2), "READ");
TTree *limit_tree = (TTree*)limit_file->Get("limit");
//This file is created earlier by running combine
//TFile *significance_file = new TFile(Form("significance_%s_%d_%d.root", model.c_str(), m1, m2), "READ");
//TTree *significance_tree = (TTree*)significance_file->Get("limit");
//This file will hold 2d histograms with limit r-values and discovery significance
TFile *f = new TFile(Form("r-values_%s.root", model.c_str()), "UPDATE");
double value = -1.0;
limit_tree->SetBranchAddress("limit", &value);
double limit = -1.0; //observed limit
double rm2s = -1.0; //expected - 2 sigma
double rm1s = -1.0; //expected - 1 sigma
double rmedian = -1.0; //expected limit
double rp1s = -1.0; //expected + 1 sigma
double rp2s = -1.0; //expected + 2 sigma
for(int i=0; i< 6; i++){
limit_tree->GetEntry(i);
if(i==0) rm2s = value;
else if(i==1) rm1s = value;
else if(i==2) rmedian = value;
else if(i==3) rp1s = value;
else if(i==4) rp2s = value;
else if(i==5) limit = value;
}
//double sig = -1.0;
//significance_tree->SetBranchAddress("limit", &sig);
//significance_tree->GetEntry(0);
delete limit_tree;
limit_file->Close();
delete limit_file;
//delete significance_tree;
//significance_file->Close();
//delete significance_file;
// // not sure what this is protecting against..
// if( value > 21263 ) exit(0);
f->cd();
TH2F *hExp = (TH2F*)f->Get("hExp");
TH2F *hObs = (TH2F*)f->Get("hObs");
TH2F *hExp1m = (TH2F*)f->Get("hExp1m");
TH2F *hExp2m = (TH2F*)f->Get("hExp2m");
TH2F *hExp1p = (TH2F*)f->Get("hExp1p");
TH2F *hExp2p = (TH2F*)f->Get("hExp2p");
//TH2F *hSig = (TH2F*)f->Get("hSig");
Fill2d(hExp , rmedian, m1, m2);
Fill2d(hObs , limit , m1, m2);
Fill2d(hExp1m, rm1s , m1, m2);
Fill2d(hExp2m, rm2s , m1, m2);
Fill2d(hExp1p, rp1s , m1, m2);
Fill2d(hExp2p, rp2s , m1, m2);
//.........这里部分代码省略.........
示例15: M2gHist
void M2gHist( UInt_t chan)
{
UInt_t cbin;
Double_t factor;
TString name;
TFile* file;
cbin = chan+1;
// Target FULL
// Prompt
TH2D *hP = (TH3D*)full.Get( "THR_M2gTaggP_v_TChanM2gP");
hP->GetXaxis()->SetRange( cbin, cbin);
TH1D *hPf_y = hP->ProjectionY( "M2gP");
// Random
TH3D *hR = (TH3D*)full.Get( "THR_M2gTaggR_v_TChanM2gR");
hR->GetXaxis()->SetRange( cbin, cbin);
TH1D *hRf_y = hR->ProjectionY( "M2gR");
hRf_y->Scale( 0.5);
// Subtracted
hSf = (TH1D*) hPf_y->Clone( "FullSubt");
hSf->Sumw2();
hSf->Add( hRf_y, -1.0);
// Target EMPTY
// Prompt
TH2D *hP = (TH3D*)empty.Get( "THR_M2gTaggP_v_TChanM2gP");
hP->GetXaxis()->SetRange( cbin, cbin);
TH1D *hPf_y = hP->ProjectionY( "M2gP");
// Random
TH3D *hR = (TH3D*)empty.Get( "THR_M2gTaggR_v_TChanM2gR");
hR->GetXaxis()->SetRange( cbin, cbin);
TH1D *hRe_y = hR->ProjectionY( "M2gR");
hRe_y->Scale( 0.5);
// Subtracted
hSe = (TH1D*) hPe_y->Clone( "EmptySubt");
hSe->Sumw2();
hSe->Add( hRe_y, -1.0);
// COMPLETELY SUBTRACTED
// Dead-time correction
file = full;
f_dead_f = DeadTimeSF( file);
file = empty;
f_dead_e = DeadTimeSF( file);
// Scale full to empty
factor = HistSF( hesc, hfsc, cbin, cbin)*f_dead_f/f_dead_e;
hSf->Scale( factor);
hS = (TH1D*) hSf->Clone( "Subtracted");
hS->Add( hSe, -1.0);
fhist = (TH1D*) hSf->Clone( "Full");
ehist = (TH1D*) hSe->Clone( "Empty");
shist = (TH1D*) hS->Clone( "Subt");
}