本文整理汇总了C++中TF1类的典型用法代码示例。如果您正苦于以下问题：C++ TF1类的具体用法？C++ TF1怎么用？C++ TF1使用的例子？那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了TF1类的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: fitResolandScale
//.........这里部分代码省略.........
// return;
// pt spectrum
// TCanvas* c1= new TCanvas("c1", "", 500,500);
// TH1D* hptHat = new TH1D("hptHat",";pt hat;Entries",200,0,200);
// yJet->Draw3( hptHat, "yPhotonTree.ptHat"," photonEt>40 && genPhotonEt> 30 && abs(genMomId)<=22","");
//yJet->Draw3( hptHat, "yPhotonTree.pt","","");
// return;
// Energy Scale
TCanvas* c2 = new TCanvas("c2", "pt/refPt distribution", 1200, 900);
TCanvas* ccc = new TCanvas("ccc", "pt/refpt 30-40GeV", 400, 400);
makeMultiPanelCanvas(c2,5,4,0.0,0.0,0.2,0.15,0.02);
TH1D* Escale[nCentBinPa+5][nPtBin];
double mean[nCentBinPa+5][nPtBin], var[nCentBinPa+5][nPtBin], resol[nCentBinPa+5][nPtBin], resolVar[nCentBinPa+5][nPtBin];
int icent =1;
// for(int icent=1; icent <= nCentBinPa ; icent++){
for(int i=0; i < nPtBin ; i++){
// c2 -> cd((icent-1)*4+i+1);
c2 -> cd(i+1);
Escale[icent][i] = new TH1D(Form("Escale%d_%d",icent, i) , " ; p_{T}^{RECO}/p_{T}^{GEN}", 50, 0, 2);
if ( genOpt == 0 ) {
yJet -> Draw2(Escale[icent][i], "pt/refPt", centCut && partonCut && Form(" (abs(eta) < 1.6) && (dphi > 7*3.141592/8.0) && (refPt >= %d && refPt < %d)", (int)ptBin[i], (int)ptBin[i+1]),"");
}
else if (genOpt == 1) {
yJet -> Draw2(Escale[icent][i], "pt/refPt", centCut && partonCut && Form(" (abs(eta) < 1.6) && (dphi > 7*3.141592/8.0) && (pt >= %d && pt < %d) ", (int)ptBin[i], (int)ptBin[i+1]),"");
}
// Escale[icent][i] -> Draw();
TF1* ff = cleverGaus(Escale[icent][i]);
gPad->SetLogy();
mean[icent][i] = ff->GetParameter(1);
var[icent][i] = ff->GetParError(1);
resol[icent][i] = ff->GetParameter(2);
resolVar[icent][i] = ff->GetParError(2);
float dx1;
// ((icent==1)||(icent==4))? dx1=0.15 : dx1=0 ;
dx1=0;
// if ( icent == nCentBinPa )
// drawText(Form("E_{T}^{HF|#eta|>4} > %dGeV, ", (int)centBinPa[icent-1]), 0.12+dx1,0.929118,1,15);//yeonju 130805
//else
// drawText(Form("%dGeV < E_{T}^{HF|#eta|>4} < %dGeV, ", (int)centBinPa[icent-1], (int)centBinPa[icent]), 0.12+dx1,0.929118,1,15);
if ( i+1 == nPtBin )
drawText(Form("p_{T}^{GEN Jet} > %dGeV, ", (int)ptBin[i]), 0.17+dx1,0.84,1,15);//yeonju 130823
else
drawText(Form("%dGeV < p_{T}^{GEN Jet} < %dGeV, ", (int)ptBin[i], (int)ptBin[i+1]), 0.17+dx1,0.84,1,12);//yeonju 130823
// TLegend *l1 = new TLegend(0.6365615,0.6445304,0.9577623,0.846736,NULL,"brNDC");
// easyLeg(l1,"p+Pb 5.02TeV");
// l1->AddEntry(hxjgNorm[kPADATA][icent + kPADATA*50][j],"pPb ","p");
// if ( icent==1 && j==1) l1->Draw();
}
// }
c2 -> Update();
ccc -> Divide(2,1);
ccc -> cd(1);
Escale[1][0]->Draw();
ccc -> cd(2);
Escale[1][1]->Draw();示例2: checkRsnPIDqa
void checkRsnPIDqa(TString filename, TString foldername, Bool_t savePng,
TString plotTPCpi, TString plotTPCka, TString plotTPCpro,
TString plotTTOFpi, TString plotTOFka, TString plotTOFpro)
{
//Open input file
TFile * fin = TFile::Open(filename.Data());
if (!fin) return 0x0;
//Access output of specific wagon
TList * list = (TList*) fin->Get(foldername.Data());
if (!list) return 0x0;
//Set range for fit
Float_t RangeFitMomMin = 0.1; //range in momentum where to check the mean and pull
Float_t RangeFitMomMax = 2.0;
Int_t xbinFitMin = 0;
Int_t xbinFitMax = -1;
Float_t RangeFitNsigmaPIDmin = -2.0; //range in Nsigma where the fit is to be performed
Float_t RangeFitNsigmaPIDmax = 2.0;
//Set range for visualisation
Float_t RangeShowTPC[2] = {0.1, 2.0};
Float_t RangeShowTOF[2] = {0.25, 2.0};
//--------------------------
// TPC PID Nsigma
// fit with simple gaussian
//--------------------------
//Gaussian function
TF1 *fGaus = new TF1("f","gaus", -7.0, 7.0);
//--- pions
TH2F * hTPCsigmaPi = (TH2F*)list->FindObject(plotTPCpi.Data());
hTPCsigmaPi->RebinX(2);
hTPCsigmaPi->SetTitle("TPC Pions");
MakeUpHisto(hTPCsigmaPi,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
hTPCsigmaPi->GetYaxis()->SetRangeUser(-5.1,5.1);
hTPCsigmaPi->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
xbinFitMin = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMin);
xbinFitMax = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMax);
hTPCsigmaPi->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
TH1D * hTPCsigmaPi_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_mean");
TH1D * hTPCsigmaPi_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_pull");
MakeUpHisto(hTPCsigmaPi_mean, "", "", 1, kBlack, 2);
MakeUpHisto(hTPCsigmaPi_pull, "", "", 1, kRed+2, 2);
//--- kaons
TH2F * hTPCsigmaKa = (TH2F*)list->FindObject(plotTPCka.Data());
hTPCsigmaKa->RebinX(2);
hTPCsigmaKa->SetTitle("TPC Kaons");
hTPCsigmaKa->GetYaxis()->SetRangeUser(-5.1,5.1);
hTPCsigmaKa->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
hTPCsigmaKa->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
MakeUpHisto(hTPCsigmaKa,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
TH1D * hTPCsigmaKa_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCka.Data())))->Clone("hNsigmaTPCka_mean");
TH1D * hTPCsigmaKa_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCka.Data())))->Clone("hNsigmaTPCka_pull");
MakeUpHisto(hTPCsigmaKa_mean, "", "", 1, kBlack, 2);
MakeUpHisto(hTPCsigmaKa_pull, "", "", 1, kRed+2, 2);
//--- protons
TH2F * hTPCsigmaPro = (TH2F*)list->FindObject(plotTPCpro.Data());
hTPCsigmaPro->RebinX(2);
hTPCsigmaPro->SetTitle("TPC Protons");
MakeUpHisto(hTPCsigmaPro,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
hTPCsigmaPro->GetYaxis()->SetRangeUser(-5.1,5.1);
hTPCsigmaPro->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
hTPCsigmaPro->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
TH1D * hTPCsigmaPro_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_mean");
TH1D * hTPCsigmaPro_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_pull");
MakeUpHisto(hTPCsigmaPro_mean, "", "", 1, kBlack, 2);
MakeUpHisto(hTPCsigmaPro_pull, "", "", 1, kRed+2, 2);
//--- plot TPC
TLine *l11=new TLine(RangeShowTPC[0],0.,RangeShowTPC[1],0.); l11->SetLineWidth(1); l11->SetLineStyle(7);
TLine *l12=new TLine(RangeShowTPC[0],1.,RangeShowTPC[1],1.); l12->SetLineWidth(1); l12->SetLineStyle(7);
gStyle->SetOptStat(0);
TCanvas *cPidPerformance4 = new TCanvas("cPIDperformance4","TPC PID",1200,500);
cPidPerformance4->Divide(3,1);
cPidPerformance4->cd(1);
gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
hTPCsigmaPi->DrawCopy("colz");
hTPCsigmaPi_mean->DrawCopy("same");
hTPCsigmaPi_pull->DrawCopy("same");
l11->Draw("same"); l12->Draw("same");
cPidPerformance4->cd(2);
gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
hTPCsigmaKa->DrawCopy("colz");
hTPCsigmaKa_mean->DrawCopy("same");
hTPCsigmaKa_pull->DrawCopy("same");
l11->Draw("same"); l12->Draw("same");
cPidPerformance4->cd(3);
gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
hTPCsigmaPro->DrawCopy("colz");
hTPCsigmaPro_mean->DrawCopy("same");
hTPCsigmaPro_pull->DrawCopy("same");
l11->Draw("same"); l12->Draw("same");
//.........这里部分代码省略.........
示例3: checkPullTree
//.........这里部分代码省略.........
TH2D* hPullAdditionalCorr = (TH2D*)hPull->Clone("hPullAdditionalCorr");
hPullAdditionalCorr->SetTitle("Pull vs. p_{TPC} integrated over tan(#Theta) with additional dEdx correction w.r.t. tan(#Theta)");
/*
const Int_t nThetaHistos = 3;
TH2D* hPullTheta[nThetaHistos];
TH2D* hPullAdditionalCorrTheta[nThetaHistos];
Double_t tThetaLow[nThetaHistos] = { 0.0, 0.4, 0.9 };
Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.5, 1.0 };
*/
const Int_t nThetaHistos = 10;
TH2D* hPullTheta[nThetaHistos];
TH2D* hPullAdditionalCorrTheta[nThetaHistos];
Double_t tThetaLow[nThetaHistos] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
for (Int_t i = 0; i < nThetaHistos; i++) {
hPullTheta[i] = new TH2D(Form("hPullTheta_%d", i),
Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f;p_{TPC} (GeV/c);Pull", tThetaLow[i], tThetaHigh[i]),
nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
hPullAdditionalCorrTheta[i] =
new TH2D(Form("hPullAdditionalCorrTheta_%d", i),
Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f with additional dEdx correction w.r.t. tan(#Theta);p_{TPC} (GeV/c);Pull",
tThetaLow[i], tThetaHigh[i]),
nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
}
TF1 corrFuncMult("corrFuncMult", "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
0., 0.2);
TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (x -[2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
TF1 corrFuncSigmaMult("corrFuncSigmaMul", "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
// LHC13b.pass2
if (isNonPP)
printf("Using corr Parameters for 13b.pass2\n!");
corrFuncMult.SetParameter(0, -5.906e-06);
corrFuncMult.SetParameter(1, -5.064e-04);
corrFuncMult.SetParameter(2, -3.521e-02);
corrFuncMult.SetParameter(3, 2.469e-02);
corrFuncMult.SetParameter(4, 0);
corrFuncMultTanTheta.SetParameter(0, -5.32e-06);
corrFuncMultTanTheta.SetParameter(1, 1.177e-05);
corrFuncMultTanTheta.SetParameter(2, -0.5);
corrFuncSigmaMult.SetParameter(0, 0.);
corrFuncSigmaMult.SetParameter(1, 0.);
corrFuncSigmaMult.SetParameter(2, 0.);
corrFuncSigmaMult.SetParameter(3, 0.);
/* OK, but PID task was not very satisfying
corrFuncMult.SetParameter(0, -6.27187e-06);
corrFuncMult.SetParameter(1, -4.60649e-04);
corrFuncMult.SetParameter(2, -4.26450e-02);
corrFuncMult.SetParameter(3, 2.40590e-02);
corrFuncMult.SetParameter(4, 0);
示例4: fit_dscb
int fit_dscb(TH1F*& hrsp,
const double nsigma,
const double jtptmin,
const int niter,
const string alg,
const double fitmin,
const double fitmax
)
{
if (0==hrsp) {
cout<<"ERROR: Empty pointer to fit_dscb()"<<endl;return -1;
}
// first use a gaussian to constrain crystal ball gaussian core
fit_gaussian(hrsp, nsigma, jtptmin, niter, alg);
TF1* fgaus = hrsp->GetFunction("fgaus");
if (0==fgaus) {
hrsp->GetListOfFunctions()->Delete();
return -1;
}
// implementation of the low pt bias threshold
string histname = hrsp->GetName();
double ptRefMax(1.0),rspMax(0.0);
int pos1 = histname.find("RefPt");
int pos2 = histname.find("to",pos1);
string ss = histname.substr(pos1+5,pos2);
if (from_string(ptRefMax,ss,std::dec)) {
if (histname.find("RelRsp")==0)
rspMax = jtptmin/ptRefMax;
if (histname.find("AbsRsp")==0)
rspMax = jtptmin-ptRefMax;
}
double fitrange_min(fitmin);
double fitrange_max(fitmax);
fitrange_min = std::max(rspMax,fitmin);
adjust_fitrange(hrsp,fitrange_min,fitrange_max);
TF1* fdscb = new TF1("fdscb",fnc_dscb,fitrange_min,fitrange_max,7);
fdscb->SetLineWidth(2);
fdscb->SetLineStyle(2);
double norm = 2*fgaus->GetParameter(0);
double mean = fgaus->GetParameter(1);
double sigma= fgaus->GetParameter(2);
//cout << hrsp->GetName() << " fgaus "<< mean << " \t " << hrsp->GetMean() << endl;
// double norm = 2*hrsp->GetMaximum();
// double mean = hrsp->GetMean();
// //double mean = GetPeak(hrsp);
// cout << " mean : " << mean << " hist mean : "<< hrsp->GetMean() << endl;
// double sigma= hrsp->GetRMS();
double aone(2.0),atwo(2.0),pone(10.0),ptwo(10.0);
//double aone(1.0),atwo(1.0),pone(10.0),ptwo(10.0);
TVirtualFitter::SetDefaultFitter("Minuit");
int fitstatus(0);
for (int i=0;i<niter;i++) {
fdscb->SetParameter(0,norm); // N
fdscb->SetParameter(1,mean); // mean
fdscb->SetParameter(2,sigma);// sigma
fdscb->SetParameter(3,aone); // a1
fdscb->SetParameter(4,pone); // p1
fdscb->SetParameter(5,atwo); // a2
fdscb->SetParameter(6,ptwo); // p2
fdscb->FixParameter(1,mean);
fdscb->FixParameter(2,sigma);
if (i>0) fdscb->FixParameter(3,aone);
else{
fdscb->SetParLimits(3,0.,20.); //! best
}
if (i>1) fdscb->FixParameter(5,atwo);
else{
fdscb->SetParLimits(5,0.,20.); //! best
}
//! best
fdscb->SetParLimits(4,0.,60.);
fdscb->SetParLimits(6,0.,60.);
fitstatus = hrsp->Fit(fdscb,"RQ+");
if (0==fitstatus) i=999;
delete fdscb;
fdscb = hrsp->GetFunction("fdscb");
if (0==fdscb) return -1;
norm = fdscb->GetParameter(0);
aone = fdscb->GetParameter(3);
pone = fdscb->GetParameter(4);
atwo = fdscb->GetParameter(5);
ptwo = fdscb->GetParameter(6);
//.........这里部分代码省略.........
示例5: ConfidenceIntervals
void ConfidenceIntervals()
{
TCanvas *myc = new TCanvas("myc",
"Confidence intervals on the fitted function",1200, 500);
myc->Divide(3,1);
/////1. A graph
//Create and fill a graph
Int_t ngr = 100;
TGraph *gr = new TGraph(ngr);
gr->SetName("GraphNoError");
Double_t x, y;
Int_t i;
for (i=0; i<ngr; i++){
x = gRandom->Uniform(-1, 1);
y = -1 + 2*x + gRandom->Gaus(0, 1);
gr->SetPoint(i, x, y);
}
//Create the fitting function
TF1 *fpol = new TF1("fpol", "pol1", -1, 1);
fpol->SetLineWidth(2);
gr->Fit(fpol, "Q");
//Create a TGraphErrors to hold the confidence intervals
TGraphErrors *grint = new TGraphErrors(ngr);
grint->SetTitle("Fitted line with .95 conf. band");
for (i=0; i<ngr; i++)
grint->SetPoint(i, gr->GetX()[i], 0);
//Compute the confidence intervals at the x points of the created graph
(TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint);
//Now the "grint" graph contains function values as its y-coordinates
//and confidence intervals as the errors on these coordinates
//Draw the graph, the function and the confidence intervals
myc->cd(1);
grint->SetLineColor(kRed);
grint->Draw("ap");
gr->SetMarkerStyle(5);
gr->SetMarkerSize(0.7);
gr->Draw("psame");
/////2. A histogram
myc->cd(2);
//Create, fill and fit a histogram
Int_t nh=5000;
TH1D *h = new TH1D("h",
"Fitted gaussian with .95 conf.band", 100, -3, 3);
h->FillRandom("gaus", nh);
TF1 *f = new TF1("fgaus", "gaus", -3, 3);
f->SetLineWidth(2);
h->Fit(f, "Q");
h->Draw();
//Create a histogram to hold the confidence intervals
TH1D *hint = new TH1D("hint",
"Fitted gaussian with .95 conf.band", 100, -3, 3);
(TVirtualFitter::GetFitter())->GetConfidenceIntervals(hint);
//Now the "hint" histogram has the fitted function values as the
//bin contents and the confidence intervals as bin errors
hint->SetStats(kFALSE);
hint->SetFillColor(2);
hint->Draw("e3 same");
/////3. A 2d graph
//Create and fill the graph
Int_t ngr2 = 100;
Double_t z, rnd, e=0.3;
TGraph2D *gr2 = new TGraph2D(ngr2);
gr2->SetName("Graph2DNoError");
TF2 *f2 = new TF2("f2",
"1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+250",-6,6,-6,6);
f2->SetParameters(1,1);
for (i=0; i<ngr2; i++){
f2->GetRandom2(x,y);
// Generate a random number in [-e,e]
rnd = 2*gRandom->Rndm()*e-e;
z = f2->Eval(x,y)*(1+rnd);
gr2->SetPoint(i,x,y,z);
}
//Create a graph with errors to store the intervals
TGraph2DErrors *grint2 = new TGraph2DErrors(ngr2);
for (i=0; i<ngr2; i++)
grint2->SetPoint(i, gr2->GetX()[i], gr2->GetY()[i], 0);
//Fit the graph
f2->SetParameters(0.5,1.5);
gr2->Fit(f2, "Q");
//Compute the confidence intervals
(TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint2);
//Now the "grint2" graph contains function values as z-coordinates
//and confidence intervals as their errors
//draw
myc->cd(3);
f2->SetNpx(30);
f2->SetNpy(30);
f2->SetFillColor(kBlue);
f2->Draw("surf4");
grint2->SetNpx(20);
grint2->SetNpy(20);
grint2->SetMarkerStyle(24);
grint2->SetMarkerSize(0.7);
//.........这里部分代码省略.........
示例6: DrawCalibrationPlotsEE
//.........这里部分代码省略.........
sigma_vs_ring[0]->SetMarkerStyle(20);
sigma_vs_ring[0]->SetMarkerSize(1);
sigma_vs_ring[0]->SetMarkerColor(kBlue+2);
sigma_vs_ring[1] = new TGraphErrors();
sigma_vs_ring[1]->SetMarkerStyle(20);
sigma_vs_ring[1]->SetMarkerSize(1);
sigma_vs_ring[1]->SetMarkerColor(kBlue+2);
sigma_vs_ring[2] = new TGraphErrors();
sigma_vs_ring[2]->SetMarkerStyle(20);
sigma_vs_ring[2]->SetMarkerSize(1);
sigma_vs_ring[2]->SetMarkerColor(kBlue+2);
/// Graph for scale vs ring EE+, EE- and folded
TGraphErrors *scale_vs_ring[3];
scale_vs_ring[0] = new TGraphErrors();
scale_vs_ring[0]->SetMarkerStyle(20);
scale_vs_ring[0]->SetMarkerSize(1);
scale_vs_ring[0]->SetMarkerColor(kBlue+2);
scale_vs_ring[1] = new TGraphErrors();
scale_vs_ring[1]->SetMarkerStyle(20);
scale_vs_ring[1]->SetMarkerSize(1);
scale_vs_ring[1]->SetMarkerColor(kBlue+2);
scale_vs_ring[2] = new TGraphErrors();
scale_vs_ring[2]->SetMarkerStyle(20);
scale_vs_ring[2]->SetMarkerSize(1);
scale_vs_ring[2]->SetMarkerColor(kBlue+2);
TF1 *fgaus = new TF1("fgaus","gaus",-10,10);
int np[3] = {0};
/// Gaussian fit for EE+ and EE-
for (int k = 0; k < 2 ; k++){
for (int iring = 0; iring < 40; iring++){
if (hspread[k][iring]-> GetEntries() == 0) continue;
float e = 0.5*hcmap[k]-> GetYaxis()->GetBinWidth(1);
fgaus->SetParameter(1,1);
fgaus->SetParameter(2,hspread[k][iring]->GetRMS());
fgaus->SetRange(1-5*hspread[k][iring]->GetRMS(),1+5*hspread[k][iring]->GetRMS());
hspread[k][iring]->Fit("fgaus","QR");
sigma_vs_ring[k]-> SetPoint(np[k],iring,fgaus->GetParameter(2));
sigma_vs_ring[k]-> SetPointError(np[k], e ,fgaus->GetParError(2));
scale_vs_ring[k]-> SetPoint(np[k],iring,fgaus->GetParameter(1));
scale_vs_ring[k]-> SetPointError(np[k],e,fgaus->GetParError(1));
np[k]++;
}
}
for (int iring = 0; iring < 40; iring++){
if (hspreadAll[iring]-> GetEntries() == 0) continue;
float e = 0.5*hcmap[0]-> GetYaxis()->GetBinWidth(1);
fgaus->SetParameter(1,1);
fgaus->SetParameter(2,hspreadAll[iring]->GetRMS());
fgaus->SetRange(1-5*hspreadAll[iring]->GetRMS(),1+5*hspreadAll[iring]->GetRMS());
hspreadAll[iring]->Fit("fgaus","QR");
sigma_vs_ring[2]-> SetPoint(np[2],iring,fgaus->GetParameter(2));
sigma_vs_ring[2]-> SetPointError(np[2], e ,fgaus->GetParError(2));
scale_vs_ring[2]-> SetPoint(np[2],iring,fgaus->GetParameter(1));
scale_vs_ring[2]-> SetPointError(np[2],e,fgaus->GetParError(1));
np[2]++; 示例7: DoAnalysisWithTree
void DoAnalysisWithTree(TFile* file, int run, int x, int y)
{
gout << 0 << " " << run << " " << x << " " << y << " ";
TF1* flandau = new TF1("flandau","landau",0,400);
TCanvas *c1 = new TCanvas("c1","",1280,960);
c1->Divide(4,2);
TCanvas *c2 = new TCanvas("c2","",1280,960);
c2->Divide(4,2);
TCanvas *c3 = new TCanvas("c3","",1280,960);
c3->Divide(4,2);
TCanvas *c4 = new TCanvas("c4","",1280,960);
c4->Divide(4,2);
TTree* tree = (TTree*)file->Get("T");
tree->SetAlias("Average_HODO_HORIZONTAL","Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
tree->SetAlias("Average_HODO_VERTICAL","Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
tree->SetAlias("Valid_HODO_HORIZONTAL","Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) > 0");
tree->SetAlias("Valid_HODO_VERTICAL","Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) > 0");
ofstream fout((const char*)Form("DataRunByRun/mapper_run%d.dat",run));
for ( int i = 0; i < 8; ++i )
{
TString drawstring1 = "TOWER_CALIB_TILE_MAPPER[";
drawstring1 += i*2;
drawstring1 += "].energy>>hs1_";
drawstring1 += i+1;
drawstring1 += "(100,0,400)";
c1->cd(i+1);
tree->Fit("flandau",drawstring1,"","","");
TString drawstring2 = "TOWER_CALIB_TILE_MAPPER[";
drawstring2 += i*2;
drawstring2 += "].energy>>hs2_";
drawstring2 += i+1;
drawstring2 += "(100,0,400)";
c2->cd(i+1);
tree->Fit("flandau",drawstring2,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL","","");
fout << run << " " << 9999 << " " << 9999 << " " << flandau->GetParameter(1) << " " << flandau->GetParError(1) << endl;
gout << flandau->GetParameter(1) << " " << flandau->GetParError(1) << " ";
TString drawstring3 = "TOWER_CALIB_TILE_MAPPER[";
drawstring3 += i*2;
drawstring3 += "].energy>>hs3_";
drawstring3 += i+1;
drawstring3 += "(100,0,400)";
c3->cd(i+1);
tree->Fit("flandau",drawstring3,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL && abs(TOWER_CALIB_C2[3].energy)<200","","");
TString drawstring4 = "TOWER_CALIB_TILE_MAPPER[";
drawstring4 += i*2;
drawstring4 += "].energy>>hs4_";
drawstring4 += i+1;
drawstring4 += "(100,0,400)";
c4->cd(i+1);
tree->Fit("flandau",drawstring4,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL && abs(TOWER_CALIB_C2[3].energy)>200","","");
}
c1->Print(Form("FigsRunByRun/mapper_run%d_step1.png",run));
c2->Print(Form("FigsRunByRun/mapper_run%d_step2.png",run));
c3->Print(Form("FigsRunByRun/mapper_run%d_step3.png",run));
c4->Print(Form("FigsRunByRun/mapper_run%d_step4.png",run));
c1->Print(Form("FigsRunByRun/mapper_run%d_step1.pdf",run));
c2->Print(Form("FigsRunByRun/mapper_run%d_step2.pdf",run));
c3->Print(Form("FigsRunByRun/mapper_run%d_step3.pdf",run));
c4->Print(Form("FigsRunByRun/mapper_run%d_step4.pdf",run));
delete c1;
delete c2;
delete c3;
delete c4;
gout << endl;
return;
// --- let's try to have a look at the hodoscope positions
TCanvas* c5 = new TCanvas("c5","",800,800);
TH2D* th2d_hodo_fine = new TH2D("th2d_hodo_fine","",40,-0.5,7.5,40,-0.5,7.5);
TH2D* th2d_hodo_coarse = new TH2D("th2d_hodo_coarse","",8,-0.5,7.5,8,-0.5,7.5);
th2d_hodo_fine->GetXaxis()->SetTitle("Horizontal Hodoscope");
th2d_hodo_fine->GetYaxis()->SetTitle("Vertical Hodoscope");
th2d_hodo_coarse->GetXaxis()->SetTitle("Horizontal Hodoscope");
th2d_hodo_coarse->GetYaxis()->SetTitle("Vertical Hodoscope");
tree->Draw("Average_HODO_HORIZONTAL:Average_HODO_VERTICAL>>th2d_hodo_fine","Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL","colz");
c5->Print(Form("FigsRunByRun/hodoscope_fine_run%d.png",run));
//.........这里部分代码省略.........
示例8: extractLimitAtQuantile
double extractLimitAtQuantile(TString inFileName, TString plotName, double d_quantile ){
TFile *f = TFile::Open(inFileName);
TF1 *expoFit = new TF1("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
TGraphErrors *limitPlot_ = new TGraphErrors();
/* bool done = false; */
if (_debug > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;
readAllToysFromFile(limitPlot_, f, d_quantile );
f->Close();
limitPlot_->Sort();
double minDist=1e3;
int n= limitPlot_->GetN();
cout<<" Number of points in limitPlot_ : "<<n<<endl;
if(n<=0) return 0;
clsMin.first=0;
clsMin.second=0;
clsMax.first=0;
clsMax.second=0;
limit = 0; limitErr = 0;
for (int i = 0; i < n; ++i) {
double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i]; //, ey = limitPlot_->GetErrorY(i);
if (fabs(y-clsTarget) < minDist) { limit = x; minDist = fabs(y-clsTarget); }
}
int ntmp =0;
for (int j = 0; j < n; ++j) {
int i = n-j-1;
double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
if (y-3*ey >= clsTarget && ntmp<=2) {
rMin = x; clsMin = CLs_t(y,ey);
ntmp ++ ;
}
}
ntmp =0;
for (int i = 0; i < n; ++i) {
double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
if (y+3*ey <= clsTarget && ntmp<=2) {
rMax = x; clsMax = CLs_t(y,ey);
ntmp ++ ;
}
}
if((clsMin.first==0 and clsMin.second==0) )
{
rMin = limitPlot_->GetX()[0]; clsMin=CLs_t(limitPlot_->GetY()[0], limitPlot_->GetErrorY(0));
}
if((clsMax.first==0 and clsMax.second==0))
{
rMax = limitPlot_->GetX()[n-1]; clsMax=CLs_t(limitPlot_->GetY()[n-1], limitPlot_->GetErrorY(n-1));
}
if (_debug > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
limitErr = std::max(limit-rMin, rMax-limit);
expoFit->SetRange(rMin,rMax);
//expoFit.SetRange(limitPlot_->GetXaxis()->GetXmin(),limitPlot_->GetXaxis()->GetXmax());
//expoFit->SetRange(1.7,2.25);
if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
if (_debug > 1) std::cout << " reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
/* done = true; */
}
//if (!done) { // didn't reach accuracy with scan, now do fit
if (1) { // didn't reach accuracy with scan, now do fit
if (_debug) {
std::cout << "\n -- HybridNew, before fit -- \n";
std::cout << "Limit: r" << " < " << limit << " +/- " << limitErr << " [" << rMin << ", " << rMax << "]\n";
std::cout<<"rMin="<<rMin<<" clsMin="<<clsMin.first<<", rMax="<<rMax<<" clsMax="<<clsMax.first<<endl;
}
expoFit->FixParameter(0,clsTarget);
expoFit->SetParameter(1,log(clsMax.first/clsMin.first)/(rMax-rMin));
expoFit->SetParameter(2,limit);
double rMinBound, rMaxBound; expoFit->GetRange(rMinBound, rMaxBound);
limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
int npoints = 0;
for (int j = 0; j < limitPlot_->GetN(); ++j) {
if (limitPlot_->GetX()[j] >= rMinBound && limitPlot_->GetX()[j] <= rMaxBound) npoints++;
}
for (int i = 0, imax = 0; i <= imax; ++i, ++npoints) {
limitPlot_->Sort();
limitPlot_->Fit(expoFit,(_debug <= 1 ? "QNR EX0" : "NR EXO"));
if (_debug) {
std::cout << "Fit to " << npoints << " points: " << expoFit->GetParameter(2) << " +/- " << expoFit->GetParError(2) << std::endl;
}
// only when both "cls+3e<0.05 and cls-3e>0.05" are satisfied, we require below ...
// if ((rMin < expoFit->GetParameter(2)) && (expoFit->GetParameter(2) < rMax) && (expoFit->GetParError(2) < 0.5*(rMaxBound-rMinBound))) {
// sanity check fit result
limit = expoFit->GetParameter(2);
limitErr = expoFit->GetParError(2);
if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) break;
// }
}
}
if (limitPlot_) {
TCanvas *c1 = new TCanvas("c1","c1");
//.........这里部分代码省略.........
示例9: mk_sigaccjecupdown
//.........这里部分代码省略.........
// h_gauswidth->SetLineColor(kBlack);
string title;
string systematic = "pile-up";
// title="Gaussian Width vs. Mass for Light-ptcut RPV";
string tag = "Heavy", sphericity = " Sphericity Cut";
if (flavor.compare("112") == 0)
tag = "Light";
else if (ptcut == "60")
sphericity = "";
string titlepart = tag + "-flavor RPV " + flavor + " p_{T} = " + ptcut + sphericity;
title = "Width for " + titlepart;
/*
if(k==0){
title="RPV gluino #bf{" + systematic + " up} m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 60 GeV}";
if (i>=5 || p==0) title="RPV gluino #bf{" + systematic + " up} m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 110 GeV}";
}
if(k==1){
title="RPV gluino #bf{" + systematic + " down} m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 60 GeV}";
if (i>=5 || p==0) title="RPV gluino #bf{" + systematic + " down} m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 110 GeV}";
}
if(k==2){
title="RPV gluino m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 60 GeV}";
if (i>=5 || p==0) title="RPV gluino m="+to_string(masses[i])+", ptcut = "+ptcut+", #Delta = 110 GeV, #bf{6^{th} Jet p_{T} = 110 GeV}";
}
*/
TLegend *leg;
leg = new TLegend(0.6,0.2,0.8994975,0.5,NULL,"brNDC");
float linesiz = 2.0;
h_gauswidthup->SetLineColor(kBlue);
h_gauswidthup->SetLineWidth(linesiz);
h_gauswidthup->SetMarkerColor(kBlue);
TF1 *fitfunc = h_gauswidthup->GetFunction("GausWidth");
if (fitfunc != NULL)
fitfunc->SetLineColor(kBlue);
leg->AddEntry(h_gauswidthup,"JES up", "L");
h_gauswidthup->GetYaxis()->SetTitle("Width [GeV]");
h_gauswidthup->GetYaxis()->SetTitleOffset(1.3);
h_gauswidthup->GetXaxis()->SetTitle("RPV Mass [GeV]");
h_gauswidthup->SetTitle(title.c_str());
h_gauswidthup->Draw("A*");
fitfunc = h_gauswidth->GetFunction("GausWidth");
if (fitfunc != NULL)
fitfunc->SetLineColor(kBlack);
else cout << "Bad func name\n";
h_gauswidth->SetLineColor(kBlack);
h_gauswidth->SetMarkerColor(kBlack);
h_gauswidth->SetLineWidth(linesiz);
// h_gauswidth->SetTitleSize(0.01);
// h_gauswidth->Draw("AL");
leg->AddEntry(h_gauswidth,"Nominal JES", "L");
h_gauswidth->Draw("same*");
h_gauswidthdown->SetLineWidth(linesiz);
h_gauswidthdown->SetMarkerColor(kRed);
h_gauswidthdown->SetLineColor(kRed);
fitfunc = h_gauswidthdown->GetFunction("GausWidth");
if (fitfunc != NULL)
fitfunc->SetLineColor(kRed);
leg->AddEntry(h_gauswidthdown,"JES down", "L");
h_gauswidthdown->Draw("same*");
leg->Draw();
cGluinoFitsOpti->Write();
cGluinoFitsOpti->SaveAs((folder + "RPVwidthjesupdn" +flavor + ptcut+uncert+postfix).c_str());
TCanvas * cGluinoFitsOpt2 = new TCanvas(("RPVacc_"+ptcut+"_"+cuts).c_str(), ("RPV_" + ptcut+"_"+cuts).c_str(), 800, 600);
leg = new TLegend(0.6,0.2,0.8994975,0.5,NULL,"brNDC");示例10: Calibrate
void Calibrate()
{
TCanvas *mycan1 = (TCanvas*)gROOT->FindObject("mycan1");
if(!mycan1)
{
mycan1 = new TCanvas("mycan1","",1200,1000);
mycan1->Divide(4,4,0.0000001,0.0000001);
}
float energy[2][32][5] = {0.};
float peak[2][32][5] = {0.};
ifstream file("SiCalibPoints.txt");
if(!file.is_open())
{
cout << "No Si Calib Data" << endl;
return;
}
ofstream out("SiRecalib.cal");
int itele = -1;
int istrip = -1;
for(int i = 0;i<64;i++)
{
file >> itele >> istrip;
itele = itele-6;
file >>energy[itele][istrip][0];
file >> energy[itele][istrip][1] >> energy[itele][istrip][2];
file >> energy[itele][istrip][3] >> energy[itele][istrip][4];
file >> peak[itele][istrip][0] >> peak[itele][istrip][1];
file >> peak[itele][istrip][2] >> peak[itele][istrip][3];
file >> peak[itele][istrip][4];
}
float slope = 0.;
float inter = 0.;
for(int i = 0;i<16;i++)
{
mycan1->cd(i+1);
TGraph *mygraph = new TGraph(5,peak[0][i],energy[0][i]);
mygraph->SetMarkerStyle(20);
mygraph->Draw("AP");
TF1 *fit = new TF1("fit","pol1",0,500);
mygraph->Fit("fit","RQ");
slope = fit->GetParameter(1);
inter = fit->GetParameter(0);
cout << istrip <<" Slope = " << slope;
cout << " inter = " << inter << endl;
out << 0 << " " << i << " " << slope << " " << inter << endl;
}
TCanvas *mycan2 = (TCanvas*)gROOT->FindObject("mycan2");
if(!mycan2)
{
mycan2 = new TCanvas("mycan2","",1200,1000);
mycan2->Divide(4,4,0.0000001,0.0000001);
}
for(int i = 0;i<16;i++)
{
int istrip = i+16;
mycan2->cd(i+1);
TGraph *mygraph = new TGraph(5,peak[0][istrip],energy[0][istrip]);
mygraph->SetTitle(Form("Strip %i",istrip));
mygraph->SetMarkerStyle(20);
mygraph->Draw("AP");
TF1 *fit = new TF1("fit","pol1",0,500);
mygraph->Fit("fit","RQ");
cout << istrip << " Slope = " << fit->GetParameter(1);
cout << " inter = " << fit->GetParameter(0) << endl;
slope = fit->GetParameter(1);
inter = fit->GetParameter(0);
out << 0 << " " << istrip << " " << slope << " " << inter << endl;
}
TCanvas *mycan3 = (TCanvas*)gROOT->FindObject("mycan3");
if(!mycan3)
{
mycan3 = new TCanvas("mycan3","",1200,1000);
mycan3->Divide(4,4,0.0000001,0.0000001);
}
for(int i = 0;i<16;i++)
{
int istrip = i;
mycan3->cd(i+1);
TGraph *mygraph = new TGraph(5,peak[1][istrip],energy[1][istrip]);
mygraph->SetTitle(Form("Strip %i",istrip));
mygraph->SetMarkerStyle(20);
//.........这里部分代码省略.........
示例11: DrawNumLayersSiECal100GeV
// /r06/lc/sg568/HCAL_Optimisation_Studies/EnergyResolutionResults/Detector_Model_103/Reco_Stage_38/Photon
// EnergyResolution_PandoraSettingsDefault_DetectorModel_103_ReconstructionVariant_38_Photon.root
void DrawNumLayersSiECal100GeV()
{
const int recoVar(71);
const int energy(100);
std::vector<int> detectorModels;
detectorModels.push_back(96);
detectorModels.push_back(97);
detectorModels.push_back(98);
detectorModels.push_back(99);
std::map<int, int> detModelToLayerNumber;
detModelToLayerNumber[96] = 30;
detModelToLayerNumber[97] = 26;
detModelToLayerNumber[98] = 20;
detModelToLayerNumber[99] = 16;
TCanvas *pTCanvas = new TCanvas();
TGraphErrors *pTGraphErrors = new TGraphErrors("EResVsLayer","EResVsLayer");
for (std::vector<int>::iterator it = detectorModels.begin(); it != detectorModels.end(); it++)
{
const int detModel(*it);
const int layerNumber(detModelToLayerNumber.find(detModel)->second);
std::string rootFile("/r06/lc/sg568/HCAL_Optimisation_Studies/EnergyResolutionResults/Detector_Model_" + NumberToString(detModel) + "/Reco_Stage_" + NumberToString(recoVar) + "/Photon/EnergyResolution_PandoraSettingsDefault_DetectorModel_" + NumberToString(detModel) + "_ReconstructionVariant_" + NumberToString(recoVar) + "_Photon.root");
std::cout << rootFile << std::endl;
TFile *pTFile = new TFile(rootFile.c_str());
std::string histogramName("Resolution_DetectorModel_" + NumberToString(detModel) + "_ReconstructionVariant_" + NumberToString(recoVar) + "/PFOEnergyHistogram_DetectorModel_" + NumberToString(detModel) + "_ReconstructionVariant_" + NumberToString(recoVar) + ";4");
std::cout << histogramName << std::endl;
TH1F *pTH1F = (TH1F*)pTFile->Get(histogramName.c_str());
std::string fitTitle = "PFOEnergyHistogramGaussianFit_DetectorModel_" + NumberToString(detModel) + "_ReconstructionVariant_" + NumberToString(recoVar) + "_Energy" + NumberToString(energy) + "GeV";
TF1 *pGaussianFit = new TF1(fitTitle.c_str(),"gaus",0,1000);
pTH1F->Fit(fitTitle.c_str());
const float fitAmplitude(pGaussianFit->GetParameter(0));
const float fitMean(pGaussianFit->GetParameter(1));
const float fitStdDev(pGaussianFit->GetParameter(2));
const float energyResolution(fitStdDev/fitMean);
const float meanError(pGaussianFit->GetParError(1));
const float meanFracError(meanError / fitMean);
const float stdDevError(pGaussianFit->GetParError(2));
const float stdDevFracError(stdDevError / fitStdDev);
const float energyResolutionError = energyResolution * std::pow( (meanFracError*meanFracError) + (stdDevFracError*stdDevFracError) ,0.5);
pTGraphErrors->SetPoint(pTGraphErrors->GetN(),layerNumber,energyResolution*100.f);
pTGraphErrors->SetPointError(pTGraphErrors->GetN()-1,0,energyResolutionError*100.f);
std::cout << "For energy : " << energy << std::endl;
std::cout << "Amplitude : " << fitAmplitude << std::endl;
std::cout << "Mean : " << fitMean << std::endl;
std::cout << "Standard Deviation : " << fitStdDev << std::endl;
std::cout << "Det model " << detModel << std::endl;
std::cout << "Energy Resolution : " << energyResolution*100 << std::endl;
}
TH2F *pAxes = new TH2F("axesEj","",100,14,32,1000,2.4,3.4);
pAxes->SetTitle("100 GeV Photon Energy Resolution vs Number of Layers in ECal (Si)");
pAxes->GetYaxis()->SetTitle("#sigma_{Reco} / E_{Reco}");
pAxes->GetXaxis()->SetTitle("Number of ECal Layers");
pAxes->Draw();
pTGraphErrors->Draw("same PL");
const std::string name("SiECal" + NumberToString(energy) + "GeVPhotonResVsLayerNumber.C");
pTCanvas->SaveAs(name.c_str());
}示例12: fitB_extend
void fitB_extend(bool ispPb=true,TString infname="",bool doweight = 1)
{
if(ispPb==true){
inputdata="/afs/cern.ch/work/w/wangj/public/nt_20140727_PAMuon_HIRun2013_Merged_y24_Using03090319Bfinder.root";
inputmc="/afs/cern.ch/work/w/wangj/public/nt_20140801_mixed_fromQMBFinder_Kp.root";
luminosity=34.6*1e-3;
outputname="../ResultsBplus/SigmaBplus_extend.root";
cut="abs(y)<2.4&&(HLT_PAMu3_v1)&&abs(mumumass-3.096916)<0.15&&mass>5&&mass<6&& isbestchi2&&trk1Pt>0.9&&chi2cl>1.32e-02&&(d0/d0Err)>3.41&&cos(dtheta)>-3.46e01&&mu1pt>1.5&&mu2pt>1.5";
}
else{
//inputdata="/data/bmeson/data//nt_20141022_PPMuon_Run2013A_PromptReco_v1_resub20141126.root";
//inputmc="/afs/cern.ch/work/w/wangj/public/nt_20140801_mixed_fromQMBFinder_Kp.root";
inputdata="/data/bmeson/data/nt_20141022_PPMuon_Run2013A_PromptReco_v1_resub20141126_HLT_PAL1DoubleMu0_HighQ_v1.root";
inputmc="/data/bmeson/MC/nt_2015251_PYTHIA6_BuJpsiK_pp_PAL1DoubleMu0_HighQ_v1.root";
luminosity=5400*1e-3;
outputname="../ResultsBplus_pp/SigmaBplus_extend.root";
cut="abs(y)<2.4&&(HLT_PAL1DoubleMu0_HighQ_v1)&&abs(mumumass-3.096916)<0.15&&mass>5&&mass<6&& isbestchi2&&trk1Pt>0.9&&chi2cl>1.32e-02&&(d0/d0Err)>3.41&&cos(dtheta)>-3.46e01&&mu1pt>1.5&&mu2pt>1.5";
}
seldata_2y=Form("%s",cut.Data());
selmc=Form("abs(y)<2.4&&gen==23333&&%s",cut.Data());
selmcgen="abs(y)<2.4&&abs(pdgId)==521&&isSignal==1";
if (doweight==0) weight="1";
if (infname=="") infname=inputdata.Data();
TFile *inf = new TFile(infname.Data());
TTree *nt = (TTree*) inf->Get("ntKp");
TFile *infMC = new TFile(inputmc.Data());
TTree *ntGen = (TTree*)infMC->Get("ntGen");
TTree *ntGen2 = (TTree*)inf->Get("ntGen");
TTree *ntMC = (TTree*)infMC->Get("ntKp");
ntGen->AddFriend(ntMC);
ntGen2->AddFriend(ntMC);
const int nBins = 5;
double ptBins[nBins+1] = {10,15,20,25,30,60};
//const int nBins = 1;
//double ptBins[nBins+1] = {10,60};
TH1D *hPt = new TH1D("hPt","",nBins,ptBins);
TH1D *hPtRecoTruth = new TH1D("hPtRecoTruth","",nBins,ptBins);
TH1D *hGenPtSelected = new TH1D("hGenPtSelected","",nBins,ptBins);
TH1D *hPtMC = new TH1D("hPtMC","",nBins,ptBins);
TH1D *hPtGen = new TH1D("hPtGen","",nBins,ptBins);
TH1D *hPtGen2 = new TH1D("hPtGen2","",nBins,ptBins);
TFile *outf = new TFile(outputname.Data(),"recreate");
for (int i=0;i<nBins+1;i++)
{
if (i==nBins) {fit(nt,ntMC,10,60,ispPb,i);continue;}
TF1 *f = fit(nt,ntMC,ptBins[i],ptBins[i+1],ispPb,i);
double yield = f->Integral(5,6)/0.02;
double yieldErr = f->Integral(5,6)/0.02*f->GetParError(0)/f->GetParameter(0);
hPt->SetBinContent(i+1,yield/(ptBins[i+1]-ptBins[i]));
hPt->SetBinError(i+1,yieldErr/(ptBins[i+1]-ptBins[i]));
}
TCanvas *c= new TCanvas("cResult","",600,600);
hPt->SetXTitle("B^{+} p_{T} (GeV/c)");
hPt->SetYTitle("Uncorrected B^{+} dN/dp_{T}");
hPt->Sumw2();
hPt->Draw();
ntMC->Project("hPtMC","pt",TCut(weight)*(TCut(selmc.Data())&&"gen==23333"));
nt->Project("hPtRecoTruth","pt",TCut(seldata_2y.Data())&&"gen==23333");
ntGen->Project("hPtGen","pt",TCut(weight)*(TCut(selmcgen.Data())));
ntGen2->Project("hPtGen2","pt",(TCut(selmcgen.Data())));
divideBinWidth(hPtRecoTruth);
hPtRecoTruth->Draw("same hist");
divideBinWidth(hPtMC);
divideBinWidth(hPtGen);
divideBinWidth(hPtGen2);
hPtMC->Sumw2();
TH1D *hEff = (TH1D*)hPtMC->Clone("hEff");
hPtMC->Sumw2();
hEff->Divide(hPtGen);
TH1D *hPtCor = (TH1D*)hPt->Clone("hPtCor");
hPtCor->Divide(hEff);
TCanvas *cCor= new TCanvas("cCorResult","",600,600);
hPtCor->SetYTitle("Corrected B^{+} dN/dp_{T}");
hPtCor->Draw();
hPtGen->Draw("same hist");
hPtGen2->Draw("same hist");
TH1D *hPtSigma= (TH1D*)hPtCor->Clone("hPtSigma");
double BRchain=6.09604e-5;
hPtSigma->Scale(1./(2*luminosity*BRchain));
hPtSigma->SetYTitle("d#sigma (B^{+})/dp_{T}");
//.........这里部分代码省略.........
示例13: getYield
TH1D* getYield(TTree* nt, TTree* ntMC, TString triggerpass, TString triggername, TString prescale, TString variable, TString varname, TString varlatex, Int_t BIN_NUM, Double_t BIN_MIN, Double_t BIN_MAX, TString addcut="")
{
TH1D* hDistrib = new TH1D(Form("h%s_Distrib_%s",triggername.Data(),varname.Data()),"",BIN_NUM,BIN_MIN,BIN_MAX);
for(float ivar=0;ivar<BIN_NUM;ivar++)
{
TCanvas* c = new TCanvas(Form("c%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",500,500);
TH1D* h = new TH1D(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),";D^{0} mass (GeV/c^{2});Candidates",60,1.7,2.0);
TH1D* hMC = new TH1D(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
TH1D* hSW = new TH1D(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
Float_t varmin = BIN_MIN+ivar*((BIN_MAX-BIN_MIN)/BIN_NUM);
Float_t varmax = BIN_MIN+(ivar+1)*((BIN_MAX-BIN_MIN)/BIN_NUM);
nt->Project(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),Form("Dmass%s",prescale.Data()),Form("%s%s&&(%s>%f&&%s<%f)%s",prefilter.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,triggerpass.Data()));
ntMC->Project(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterMC.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
ntMC->Project(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterSW.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
h->SetMaximum(h->GetMaximum()*1.20);
h->Draw();
TF1* f = new TF1(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.14159)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.14159)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.14159)*[8]*(1+[11])))+[3]+[4]*x+[5]*x*x+[6]*x*x*x", 1.7, 2.0);
f->SetParLimits(4,-1000,1000);
f->SetParLimits(10,0.001,0.05);
f->SetParLimits(2,0.01,0.1);
f->SetParLimits(8,0.02,0.2);
f->SetParLimits(7,0,1);
f->SetParLimits(9,0,1);
f->SetParameter(0,setparam0);
f->SetParameter(1,setparam1);
f->SetParameter(2,setparam2);
f->SetParameter(10,setparam10);
f->SetParameter(9,setparam9);
f->FixParameter(8,setparam8);
f->FixParameter(7,1);
f->FixParameter(1,fixparam1);
f->FixParameter(3,0);
f->FixParameter(4,0);
f->FixParameter(5,0);
f->FixParameter(6,0);
f->FixParameter(11,0);
h->GetEntries();
hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
f->ReleaseParameter(1);
hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);
f->FixParameter(1,f->GetParameter(1));
f->FixParameter(2,f->GetParameter(2));
f->FixParameter(10,f->GetParameter(10));
f->FixParameter(9,f->GetParameter(9));
f->FixParameter(7,0);
f->ReleaseParameter(8);
f->SetParameter(8,setparam8);
hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);
f->FixParameter(7,hMC->Integral(0,1000)/(hSW->Integral(0,1000)+hMC->Integral(0,1000)));
f->FixParameter(8,f->GetParameter(8));
f->ReleaseParameter(3);
f->ReleaseParameter(4);
f->ReleaseParameter(5);
f->ReleaseParameter(6);
f->SetLineColor(kRed);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
f->ReleaseParameter(1);
f->SetParLimits(1,1.86,1.87);
f->ReleaseParameter(11);
f->SetParLimits(11,-1.,1.);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);
h->SetMarkerSize(0.8);
h->SetMarkerStyle(20);
TF1* background = new TF1(Form("background%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]+[1]*x+[2]*x*x+[3]*x*x*x");
background->SetParameter(0,f->GetParameter(3));
background->SetParameter(1,f->GetParameter(4));
background->SetParameter(2,f->GetParameter(5));
background->SetParameter(3,f->GetParameter(6));
background->SetLineColor(4);
background->SetRange(1.7,2.0);
background->SetLineStyle(2);
TF1* mass = new TF1(Form("fmass%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.14159)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.14159)*[5]*(1+[6]))))");
mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10),f->GetParameter(11));
mass->SetParError(0,f->GetParError(0));
mass->SetParError(1,f->GetParError(1));
mass->SetParError(2,f->GetParError(2));
mass->SetParError(3,f->GetParError(7));
mass->SetParError(4,f->GetParError(9));
mass->SetParError(5,f->GetParError(10));
//.........这里部分代码省略.........
示例14: TCanvas
TF1 *fit(TTree *nt,TTree *ntMC,double ptmin,double ptmax, bool ispPb, int count){
//cout<<cut.Data()<<endl;
//static int count=0;
//count++;
TCanvas *c= new TCanvas(Form("c%d",count),"",600,600);
TH1D *h = new TH1D(Form("h%d",count),"",50,5,6);
TH1D *hMC = new TH1D(Form("hMC%d",count),"",50,5,6);
TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
TF1 *f = new TF1(Form("f%d",count),"[0]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*("+iNP+")");
nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));
ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));
clean0(h);
TH1D *hraw = new TH1D(Form("hraw%d",count),"",50,5,6);
clean0(hraw);
hraw = (TH1D*)h->Clone(Form("hraw%d",count));
h->Draw();
f->SetParLimits(4,-1000,0);
f->SetParLimits(2,0.01,0.05);
f->SetParLimits(8,0.01,0.05);
f->SetParLimits(7,0,1);
f->SetParLimits(5,0,1000);
f->SetParameter(0,setparam0);
f->SetParameter(1,setparam1);
f->SetParameter(2,setparam2);
f->SetParameter(8,setparam3);
f->FixParameter(1,fixparam1);
h->GetEntries();
hMC->Fit(Form("f%d",count),"q","",5,6);
hMC->Fit(Form("f%d",count),"q","",5,6);
f->ReleaseParameter(1);
hMC->Fit(Form("f%d",count),"L q","",5,6);
hMC->Fit(Form("f%d",count),"L q","",5,6);
hMC->Fit(Form("f%d",count),"L q","",5,6);
hMC->Fit(Form("f%d",count),"L m","",5,6);
f->FixParameter(1,f->GetParameter(1));
f->FixParameter(2,f->GetParameter(2));
f->FixParameter(7,f->GetParameter(7));
f->FixParameter(8,f->GetParameter(8));
h->Fit(Form("f%d",count),"q","",5,6);
h->Fit(Form("f%d",count),"q","",5,6);
f->ReleaseParameter(1);
h->Fit(Form("f%d",count),"L q","",5,6);
h->Fit(Form("f%d",count),"L q","",5,6);
h->Fit(Form("f%d",count),"L q","",5,6);
h->Fit(Form("f%d",count),"L m","",5,6);
h->SetMarkerSize(0.8);
h->SetMarkerStyle(20);
cout <<h->GetEntries()<<endl;
// function for background shape plotting. take the fit result from f
TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
background->SetParameter(0,f->GetParameter(3));
background->SetParameter(1,f->GetParameter(4));
background->SetParameter(2,f->GetParameter(5));
background->SetParameter(3,f->GetParameter(6));
background->SetLineColor(4);
background->SetRange(5,6);
background->SetLineStyle(2);
// function for signal shape plotting. take the fit result from f
TF1 *Bkpi = new TF1(Form("fBkpi%d",count),"[0]*("+iNP+")");
Bkpi->SetParameter(0,f->GetParameter(5));
Bkpi->SetLineColor(kGreen+1);
Bkpi->SetFillColor(kGreen+1);
// Bkpi->SetRange(5.00,5.28);
Bkpi->SetRange(5.00,6.00);
Bkpi->SetLineStyle(1);
Bkpi->SetFillStyle(3004);
// function for signal shape plotting. take the fit result from f
TF1 *mass = new TF1(Form("fmass%d",count),"[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(8));
mass->SetParError(0,f->GetParError(0));
mass->SetParError(1,f->GetParError(1));
mass->SetParError(2,f->GetParError(2));
mass->SetParError(7,f->GetParError(7));
mass->SetParError(8,f->GetParError(8));
mass->SetLineColor(2);
mass->SetLineStyle(2);
// cout <<mass->Integral(0,1.2)<<" "<<mass->IntegralError(0,1.2)<<endl;
h->SetMarkerStyle(24);
h->SetStats(0);
h->Draw("e");
h->SetXTitle("M_{B} (GeV/c^{2})");
h->SetYTitle("Entries / (20 MeV/c^{2})");
h->GetXaxis()->CenterTitle();
h->GetYaxis()->CenterTitle();
h->SetTitleOffset(1.5,"Y");
h->SetAxisRange(0,h->GetMaximum()*1.2,"Y");
Bkpi->Draw("same");
background->Draw("same");
//.........这里部分代码省略.........
示例15: fit
TF1* fit(Double_t ptmin, Double_t ptmax)
{
TCanvas* c = new TCanvas(Form("c_%.0f_%.0f",ptmin,ptmax),"",600,600);
TFile* infile = new TFile(Form("%s_%s_%.0f_%.0f.root",infname.Data(),collisionsystem.Data(),ptmin,ptmax));
TH1D* h = (TH1D*)infile->Get("h"); h->SetName(Form("h_%.0f_%.0f",ptmin,ptmax));
TH1D* hMCSignal = (TH1D*)infile->Get("hMCSignal"); hMCSignal->SetName(Form("hMCSignal_%.0f_%.0f",ptmin,ptmax));
TH1D* hMCSwapped = (TH1D*)infile->Get("hMCSwapped"); hMCSwapped->SetName(Form("hMCSwapped_%.0f_%.0f",ptmin,ptmax));
TF1* f = new TF1(Form("f_%.0f_%.0f",ptmin,ptmax),"[0]*([7]*([9]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[9])*Gaus(x,[1],[10])/(sqrt(2*3.14159)*[10]))+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*x*x", 1.7, 2.0);
f->SetParLimits(4,-1000,1000);
f->SetParLimits(10,0.001,0.05);
f->SetParLimits(2,0.01,0.1);
f->SetParLimits(8,0.02,0.2);
f->SetParLimits(7,0,1);
f->SetParLimits(9,0,1);
f->SetParameter(0,setparam0);
f->SetParameter(1,setparam1);
f->SetParameter(2,setparam2);
f->SetParameter(10,setparam10);
f->SetParameter(9,setparam9);
f->FixParameter(8,setparam8);
f->FixParameter(7,1);
f->FixParameter(1,fixparam1);
f->FixParameter(3,0);
f->FixParameter(4,0);
f->FixParameter(5,0);
h->GetEntries();
hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
f->ReleaseParameter(1);
hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
f->FixParameter(1,f->GetParameter(1));
f->FixParameter(2,f->GetParameter(2));
f->FixParameter(10,f->GetParameter(10));
f->FixParameter(9,f->GetParameter(9));
f->FixParameter(7,0);
f->ReleaseParameter(8);
f->SetParameter(8,setparam8);
hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
f->FixParameter(7,hMCSignal->Integral(0,1000)/(hMCSwapped->Integral(0,1000)+hMCSignal->Integral(0,1000)));
f->FixParameter(8,f->GetParameter(8));
f->ReleaseParameter(3);
f->ReleaseParameter(4);
f->ReleaseParameter(5);
f->SetLineColor(kRed);
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
f->ReleaseParameter(1);
//f->ReleaseParameter(2); // you need to release these two parameters if you want to perform studies on the sigma shape
//f->ReleaseParameter(10); // you need to release these two parameters if you want to perform studies on the sigma shape
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
TF1* background = new TF1(Form("background_%.0f_%.0f",ptmin,ptmax),"[0]+[1]*x+[2]*x*x");
background->SetParameter(0,f->GetParameter(3));
background->SetParameter(1,f->GetParameter(4));
background->SetParameter(2,f->GetParameter(5));
background->SetLineColor(4);
background->SetRange(minhisto,maxhisto);
background->SetLineStyle(2);
TF1* mass = new TF1(Form("fmass_%.0f_%.0f",ptmin,ptmax),"[0]*([3]*([4]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[4])*Gaus(x,[1],[5])/(sqrt(2*3.14159)*[5])))");
mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10));
mass->SetParError(0,f->GetParError(0));
mass->SetParError(1,f->GetParError(1));
mass->SetParError(2,f->GetParError(2));
mass->SetParError(3,f->GetParError(7));
mass->SetParError(4,f->GetParError(9));
mass->SetParError(5,f->GetParError(10));
mass->SetFillColor(kOrange-3);
mass->SetFillStyle(3002);
mass->SetLineColor(kOrange-3);
mass->SetLineWidth(3);
mass->SetLineStyle(2);
TF1* massSwap = new TF1(Form("fmassSwap_%.0f_%.0f",ptmin,ptmax),"[0]*(1-[2])*Gaus(x,[1],[3])/(sqrt(2*3.14159)*[3])");
massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8));
massSwap->SetParError(0,f->GetParError(0));
massSwap->SetParError(1,f->GetParError(1));
massSwap->SetParError(2,f->GetParError(7));
massSwap->SetParError(3,f->GetParError(8));
massSwap->SetFillColor(kGreen+4);
massSwap->SetFillStyle(3005);
massSwap->SetLineColor(kGreen+4);
massSwap->SetLineWidth(3);
//.........这里部分代码省略.........