本文整理汇总了C++中TH1F::Integral方法的典型用法代码示例。如果您正苦于以下问题：C++ TH1F::Integral方法的具体用法？C++ TH1F::Integral怎么用？C++ TH1F::Integral使用的例子？那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TH1F的用法示例。
在下文中一共展示了TH1F::Integral方法的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: createNNLOplot
//.........这里部分代码省略.........
// perform a linear fit wrt previous point
// get the two points (this bin and the previous one)
int binPrev= (bin==1&&variable=="topPt") ? 0 : hist->FindBin(Xvalues[bin-2]+0.5*binwidth);
double x2=-1;
double x1= 0;
double y2=-1;
double y1= 0.;
rawHist->GetPoint(bin-1, x2, y2);
x2+=0.5*binwidth;
if(bin==1&&variable=="topPt"){
y1=0;
x1=0;
}
else{
rawHist->GetPoint(bin-2, x1, y1);
x1+=0.5*binwidth;
}
// calculate linear funtion
double a=(y2-y1)/(x2-x1);
double b=y1-a*x1;
TF1* linInterpol=new TF1("linInterpol"+getTStringFromInt(bin), "[0]*x+[1]", x1, x2);
linInterpol->SetParameter(0,a);
linInterpol->SetParameter(1,b);
double xlow = (bin==1&&variable=="topPt") ? 0. : hist->GetBinLowEdge(binPrev+1);
double xhigh = hist->GetBinLowEdge(bin2+1);
std::cout << " lin. interpolation [ (" << x1 << "," << y1 << ") .. (" << x2 << "," << y2 << ") ]: " << a << "*x+" << b << std::endl;
linInterpol->SetRange(xlow, xhigh);
hist->Add(linInterpol);
}
}
}
// check theory curve
std::cout << std::endl << "analyze theory curve:" << std::endl;
double integralRawTheory=hist->Integral(0.,hist->GetNbinsX()+1);
std::cout << "Integral: " << integralRawTheory << std::endl;
if(integralRawTheory==0.){
std::cout << "ERROR: Integral can not be 0!" << std::endl;
exit(0);
}
std::cout << "binwidth: " << binwidth << std::endl;
std::cout << "Integral*binwidth: " << integralRawTheory*binwidth << std::endl;
std::cout << "binRange: " << xmin << ".." << xmax << std::endl;
std::cout << " -> reco range: " << low << ".." << high << std::endl;
std::vector<double> recoBins_=bins_[variable];
for(unsigned int i=0; i<recoBins_.size(); ++i){
i==0 ? std::cout << "(" : std::cout << ", ";
std::cout << recoBins_[i];
if(i==recoBins_.size()-1) std::cout << ")" << std::endl;
}
// check if you need to divide by binwidth
if(std::abs(1.-integralRawTheory)<0.03){
std::cout << "Integral is approx. 1 -> need to divide by binwidth!" << std::endl;
hist->Scale(1./binwidth);
}
else if(std::abs(1-integralRawTheory*binwidth)<0.03){
std::cout << "Integral*binwidth is approx. 1 -> no scaling needed!" << std::endl;
}
else{
std::cout << "need to normalize and divide by binwidth";
hist->Scale(1./(binwidth*integralRawTheory));
}
// styling
hist->GetXaxis()->SetRangeUser(low,high);
hist->SetMarkerColor(kMagenta);
hist->SetLineColor( kMagenta);
hist->SetMarkerStyle(24);示例2: fit_h_Xm_extended
void fit_h_Xm_extended() {
TCanvas* c1;
c1 = new TCanvas("c1","",1440,900);
TFile *f;
int mass[13]={600,800,1000,1200,1400,1600,1800,2000,2500,3000,3500,4000,4500};
float width [4]={0.05,0.1,0.2,0.3};
string widthdata[4]={"0.05","0.1","0.2","0.3"};
//for(int i=0 ;i<4;i++){
//for(int j=0)
//}
for(int i=0;i<4;i++){
for(int j=0;j<13;j++){
f = TFile::Open(Form("offshell_root_files/Zprime_Zh_Zlephbb_w%s_M%d.root",widthdata[i].data(),mass[j]));
TH1F* h =(TH1F*) f->FindObjectAny("h_Xm_extended");
int bmin=0,bmax=0;
for (int k=1;k<25001;k++){
bmin=k;
if (h->GetBinContent(k)!=0) break;
}
for (int k=25000;k>0;k--){
bmax=k;
if (h->GetBinContent(k)!=0) break;
}
if((bmin-300)<0)bmin=0;
else bmin=bmin-300;
bmax=bmax+bmax/10;
cout <<bmin<<","<<bmax<<endl;
//gStyle->SetOptFit(111111);
//gStyle->SetStatW (0.1);
//gStyle->SetStatH (0.1);
//gStyle->SetStatX (0.941349);
//gStyle->SetStatY (0.912491);
fiveWidthLow= mass[j] -5* (width [i])*mass[j];
fiveWidthUp= mass[j] +5* (width [i])*mass[j];
float inte=h->Integral(fiveWidthLow,fiveWidthUp);
inte=inte*100/(bmax-bmin);
//h->Rebin((bmax-bmin)/100);
int massfix=200;
if(i>1&&j>8)massfix=400;
if(i>2&&j>9)massfix=500;
float widthfix=0;
if(i>2&&j>9)widthfix=0.1;
TF1* f1 = new TF1("f1","[2]*TMath::BreitWigner(x,[0],[1])");
f1->SetNpx(2500);
f1->SetParameters(mass[j],width [i]*mass[j],h->Integral());
//f1->SetLineColor(3);
h->GetXaxis()->SetRangeUser(bmin,bmax);
h->SetTitle(Form("offshell_root_files/Zprime_Zh_Zlephbb_w%s_M%d.root",widthdata[i].data(),mass[j]));
TH1F* h2=(TH1F*) h->Clone("h2");
//h->Fit(f1,"","",fiveWidthLow,fiveWidthUp);
//h->Fit(f1,"L");
TF1* f2 = new TF1("f2",fline,bmin,bmax,3);
if(j>2)inte=h->Integral();
f2->SetParameters(mass[j],width [i]*mass[j],h->Integral());
f2->SetNpx(2500);
//f2->SetLineColor(3);
//h2->Fit(f2);
//TF1* fleft = new TF1("fleft","[2]*TMath::BreitWigner(x,[0],[1])");
//TF1 *fleft = new TF1("fleft","[2]*TMath::BreitWigner(x,[0],[1])",bmin,bmax,3);
//fleft->SetParameters(f2->GetParameters());
//fleft->SetLineColor(3);
h->SetMinimum(0);
h->Draw();
//fleft->Draw();
//h2->Draw("same");
//f1->Draw();
if (i==0 && j==0)c1->Print("fit_h_Xm_extended.pdf(");
else if (i==3 && j==12)c1->Print("fit_h_Xm_extended.pdf)");
else c1->Print("fit_h_Xm_extended.pdf");
c1->SaveAs(Form("png/Zprime_Zh_Zlephbb_w%s_M%d.png",widthdata[i].data(),mass[j]));
//.........这里部分代码省略.........
示例3: RunMakeRazorPlots
//.........这里部分代码省略.........
if (!(MR > 400 && Rsq > 0.25)) continue;
// if (!(MR < 500 )) continue;
// if (!(Rsq < 0.3)) continue;
if (!(fabs(dPhiRazor) > 2.8)) continue;
if (!hasSignal || i>1) {
histMRAllBkg->Fill(MR, intLumi*weight);
histRsqAllBkg->Fill(Rsq, intLumi*weight);
histMRRsq[i]->Fill(MR, Rsq, intLumi*weight);
}
if(i==1){
if (intLumi*weight > 30) continue;
float qcdweight = 1.56841;
histMRQCD->Fill(MR, intLumi*weight*qcdweight);
histRsqQCD->Fill(Rsq, intLumi*weight*qcdweight);
histMRRsq[i]->Fill(MR, Rsq, intLumi*weight*qcdweight);
}
if (hasSignal && i==0) {
histMRData->Fill(MR);
histRsqData->Fill(Rsq);
histMRRsq[i]->Fill(MR, Rsq);
}
}
inputFile->Close();
delete inputFile;
}
std::cout<<"Data: "<<histMRData->Integral()<<", All Backgrounds: "<<histMRAllBkg->Integral()<<" , QCD: "<<histMRQCD->Integral()<<", Data2 "<< histMRRsq[0]->Integral() <<" QCD2 "<<histMRRsq[1]->Integral() <<std::endl;
//*******************************************************************************************
//Draw Plots
//*******************************************************************************************
// fill out the unrolled histograms
std::cout<<"Rsq bins: "<<nRsqBins<<" "<<nMRBins<<std::endl;
for (uint i=0; i < histMRRsq.size(); ++i) {
int binN = 0;
for(int ii = 0; ii<nMRBins; ii++)
for (int jj = 0; jj<nRsqBins; jj++)
{
float value = (histMRRsq[i]->GetBinContent(ii+1, jj+1) > 0) ? histMRRsq[i]->GetBinContent(ii+1, jj+1) : 0. ;
float Xrange = histMRRsq[i]->GetXaxis()->GetBinLowEdge(jj+2) - histMRRsq[i]->GetXaxis()->GetBinLowEdge(jj+1);
float Yrange = histMRRsq[i]->GetYaxis()->GetBinLowEdge(ii+2) - histMRRsq[i]->GetYaxis()->GetBinLowEdge(ii+1);
float area =1.;
if(density) area = Xrange*Yrange; //normalize each bin by its area
histUnrolled[i]->SetBinContent(binN+1, value/area);
binN++;
}
histUnrolled[i]->SetMinimum(0.00001);
if ( histUnrolled[i]->Integral() > 0) {
if( !hasSignal || i > 0 )
stackUnrolled->Add(histUnrolled[i]);示例4: FitHistos
void ClusterWidthAnalysisTreeMaker::FitHistos(std::map<ULong64_t , std::vector<TH1F*> > &HistSoN, string output_file,
std::vector< TH1F* > commonHistos, std::map<ULong64_t, TProfile* > Monitors){
TFile * myFile = new TFile(output_file.c_str(), "recreate");
ULong64_t detid;
double voltage;
double errvoltage;
double Mean;
double errMean;
double RMS;
double errRMS;
int index;
int nhits;
TTree *tree = new TTree("T", "summary information");
tree->Branch("DetID",&detid, "DetID/l");
tree->Branch("Voltage",&voltage,"Voltage/D");
tree->Branch("Index",&index,"Index/I");
tree->Branch("errVoltage",&errvoltage,"errVoltage/D");
tree->Branch("Mean",&Mean,"Mean/D");
tree->Branch("errMean",&errMean,"errMean/D");
tree->Branch("RMS",&RMS,"RMS/D");
tree->Branch("errRMS",&errRMS,"errRMS/D");
tree->Branch("Nhits",&nhits,"Nhits/I");
//TCanvas* c1 = new TCanvas();
TH1F* hNhits = new TH1F("hNhits", "hNhits", 1000, 0,1000); // N hits per module
unsigned int nfitrm=0;
for(std::map<ULong64_t , std::vector<TH1F*> >::iterator iter = HistSoN.begin(); iter != HistSoN.end(); ++iter){
unsigned int i=0; // voltage index
std::set< int >::iterator itVolt;
std::set< int > Voltage = VSmaker.getVoltageList();
for( itVolt=Voltage.begin(); itVolt!=Voltage.end(); itVolt++){
//std::cout<<"going through the measurement: " << i << std::endl;
TString thestring;
thestring.Form("DetID_%llu_%u",iter->first,i);
//std::cout << "searching for " << thestring.Data() << std::endl;
//TH1F* SoNHisto= (TH1F*)gROOT->FindObject( thestring.Data() );
if(i>=iter->second.size())
{ std::cout<<" Wrong number of voltage steps. "<<std::endl; i++; continue;}
TH1F* Histo = iter->second[i];
if(!Histo)
{ std::cout<<" Histo "<<thestring.Data()<<"_"<<i<<" not found."<<std::endl; i++; continue;}
if(Histo->GetEntries()) hNhits->Fill(Histo->Integral());
if(Histo->Integral()<20) //0.1
{ //std::cout<<" Not enought entries for histo "<<thestring.Data()<<std::endl;
i++; continue;}
detid = iter->first;
bool rmfit=false;
if( rmfit ||
// TIB modules
// TIB - 1.4.2.5
detid==369121605 || detid==369121606 || detid==369121614 ||
detid==369121613 || detid==369121610 || detid==369121609 ||
// TIB - 1.2.2.1
detid==369121390 || detid==369121382 || detid==369121386 ||
detid==369121385 || detid==369121389 || detid==369121381 ||
// others in TIB
detid==369121437 || detid==369142077 || detid==369121722 ||
detid==369125534 || detid==369137018 || detid==369121689 ||
detid==369121765 || detid==369137045 || detid==369169740 ||
detid==369121689 ||
// TOB modules
// TOB + 4.3.3.8
detid/10==436281512 || detid/10==436281528 || detid/10==436281508 ||
detid/10==436281524 || detid/10==436281520 || detid/10==436281516 ||
// others in TOB
detid/10==436228249 || detid/10==436232694 || detid/10==436228805 ||
detid/10==436244722 || detid/10==436245110 || detid/10==436249546 ||
detid/10==436310808 || detid/10==436312136 || detid/10==436315600 ||
// without 'sensors' option
detid==436281512 || detid==436281528 || detid==436281508 ||
detid==436281524 || detid==436281520 || detid==436281516 ||
detid==436228249 || detid==436232694 || detid==436228805 ||
detid==436244722 || detid==436245110 || detid==436249546 ||
detid==436310808 || detid==436312136 || detid==436315600 ||
// TID modules
detid==402664070 || detid==402664110 ||
// TEC modules in small scans
detid==470148196 || detid==470148200 || detid==470148204 ||
detid==470148228 || detid==470148232 || detid==470148236 ||
detid==470148240 || detid==470148261 || detid==470148262 ||
detid==470148265 || detid==470148266 || detid==470148292 ||
//.........这里部分代码省略.........
示例5: phi2
//.........这里部分代码省略.........
leg->AddEntry(hbkg2,"Drell Yann","f");
leg->AddEntry(hbkg3,"#gamma + Jets","f");
leg->AddEntry(h2,"m_{Z'} = 800 GeV","l");
leg->AddEntry(hbkg5,"QCD","f");
leg->AddEntry(hbkg4,"ggH","f");
leg->AddEntry(h3,"m_{Z'} = 1000 GeV","l");
leg->AddEntry(hbkg6,"VH","f");
leg->AddEntry(hbkg7,"ttH","f");
leg->AddEntry(h4,"m_{Z'} = 1200 GeV","l");
leg->AddEntry(hbkg8,"VBF H","f");
leg->AddEntry(hbkg9,"t + #gamma + Jets","f");
leg->AddEntry(h5,"m_{Z'} = 1400 GeV","l");
leg->AddEntry(hbkg10,"tt + #gamma +Jets","f");
leg->AddEntry(hbkg11,"#gamma+W","f");
leg->AddEntry(h6,"m_{Z'} = 1700 GeV","l");
leg->AddEntry(hbkg12,"#gamma+Z","f");
leg->AddEntry(hsum,"Bkg uncertainty","f");
leg->AddEntry(h7,"m_{Z'} = 2500 GeV","l");
leg->Draw("same");
c1->cd();
smallPad->Draw();
smallPad->cd();
TGraphErrors *gr = new TGraphErrors(0);
double integralData=hdata->Integral();
double integralBKG=hsum->Integral();
double error, ratio;
for(int w=1; w<20; w++){
if((hdata->GetBinContent(w)!=0) && (hsum->GetBinContent(w)!=0)){
gr->SetPoint(w, hdata->GetBinCenter(w),(hdata->GetBinContent(w))/(hsum->GetBinContent(w)));
ratio= (hdata->GetBinContent(w))/(hsum->GetBinContent(w));
error= (hdata->GetBinContent(w)*sqrt(hsum->GetBinContent(w))/(hsum->GetBinContent(w)*hsum->GetBinContent(w)) + sqrt(hdata->GetBinContent(w))/hsum->GetBinContent(w));
std::cout<<"VALUE: "<<ratio<<" ERROR: "<<error<<std::endl;
gr->SetPointError(w, hdata->GetBinWidth(w)/2,error);
}else{
gr->SetPoint(w, hdata->GetBinCenter(w),10);
}
}
gStyle->SetPadTickY(1);
gStyle->SetPadTickX(1);
gr->GetHistogram()->SetMaximum(2);
gr->GetHistogram()->SetMinimum(0.1);
gStyle->SetTextSize(14);
gROOT->ForceStyle();
gr->GetXaxis()->SetLabelFont(43);
gr->GetXaxis()->SetLabelSize(15);
gr->GetYaxis()->SetLabelFont(43);
gr->GetYaxis()->SetLabelSize(15);
gr->GetXaxis()->SetLimits(-4,4);
示例6: spectraSysTrig_dNdetashift
void spectraSysTrig_dNdetashift() {
// === General Settings ===
// Note centrality is by 2.5% bins, low bin is inclusive, high bin is exclusive
// i.e. for 0-10% use minCent=0, maxCent=4.
int minCent[12] = {0,2,4,6,8,10,12,14,16,20,24,28};
int maxCent[12] = {2,4,6,8,10,12,14,16,20,24,28,32};
// absolute values of eta are used, negative eta values are included
// for full eta range set etaMin=0 etaMax=2.4
Double_t etaMin = 0.0;
Double_t etaMax = 0.4;
double dndeta[12]= {1612,1313,1082,893.9,731.6,596.8, 482.3,383.9, 266.5, 153.2,79.5, 36.3 };
double dndetaerr[12] = {55,45, 38, 31.4, 26.3, 23.1, 18.7, 16.2, 14.2, 9.7, 6.57, 3.99};
double ptbins[20]={0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.6,2.0,2.5,3.0,3.5,4.0,5.0,6.0,8.0,10.0,12.0,16.,22.};
// Input File
TFile *f = new TFile("LowPtSpectrum_fine_Full_d.root");
TFile *fcorr = new TFile("validation3D_HydjetNoWedge_100k_flowSQ_vertexZ10.root");
TFile *fcorrAMPT;
TFile *fcorrDataMix = new TFile("validation3D_20pionBadWedge_SQ12_vertexZ10.root");
// =========================
gStyle->SetErrorX(0);
gStyle->SetOptStat(0);
set_plot_style();
char dirname[256] = "spectrumGoodTightdz10chi40";
TH1D * spec[10][12];
TH1D * eff[10][12];
TH1D * fak[10][12];
TH1D * sim[10][12];
TH1D * rec[10][12];
TGraphErrors * gSpectrum[10][12];
TF1 * f1[10][12];
Double_t meanpt[10][12];
Double_t meanpterr[10][12];
double integral[4][12];
Double_t px[100],py[100],pxe[100],pye[100];
for( int c=0; c<12; c++)
{
for( int j=0; j<4; j++)
{
if ( j==0) spec[j][c] = getSpectrum( f, Form("%s/tracks3D",dirname), minCent[c], maxCent[c], etaMin, etaMax );
if ( j==1) spec[j][c] = getSpectrum( f, Form("%s/tracks3D_c94",dirname), minCent[c], maxCent[c], etaMin, etaMax );
if ( j==2) spec[j][c] = getSpectrum( f, Form("%s/tracks3D_c97",dirname), minCent[c], maxCent[c], etaMin, etaMax );
if ( j==3) spec[j][c] = getSpectrum( f, Form("%s/tracks3D_c100",dirname), minCent[c], maxCent[c], etaMin, etaMax );
eff[j][c] = getSpectrum( fcorr, "hitrkPixelEffAnalyzer/heff3D", minCent[c], maxCent[c], etaMin, etaMax );
sim[j][c] = getSpectrum( fcorr, "hitrkPixelEffAnalyzer/hsim3D", minCent[c], maxCent[c], etaMin, etaMax );
fak[j][c] = getSpectrum( fcorr, "hitrkPixelEffAnalyzer/hfak3D", minCent[c], maxCent[c], etaMin, etaMax );
rec[j][c] = getSpectrum( fcorr, "hitrkPixelEffAnalyzer/hrec3D", minCent[c], maxCent[c], etaMin, etaMax );
// determine correction factors
TH1F * nevts;
if( j==0) nevts = (TH1F *) f->Get(Form("%s/nevts",dirname));
if( j==1) nevts = (TH1F *) f->Get(Form("%s/nevts_c94",dirname));
if( j==2) nevts = (TH1F *) f->Get(Form("%s/nevts_c97",dirname));
if( j==3) nevts = (TH1F *) f->Get(Form("%s/nevts_c100",dirname));
double Nevt = nevts->Integral( nevts->GetXaxis()->FindBin(minCent[c]+0.001),
nevts->GetXaxis()->FindBin(maxCent[c]-0.001) );
Double_t maxy = 0.;
Double_t miny = 100000.;
integral[j][c] = 0;
for( int i=0; i<=19;i++)
{
double ptmin = ptbins[i]; double ptmax = ptbins[i+1];
double iptmin = spec[j][c]->FindBin(ptmin+1e-3);
double iptmax = spec[j][c]->FindBin(ptmax-1e-3);
double icptmin = eff[j][c]->FindBin(ptmin+1e-3);
double icptmax = eff[j][c]->FindBin(ptmax-1e-3);
double pt = 0.;
for( int k=iptmin;k<=iptmax;k++) pt += spec[j][c]->GetBinCenter(k) * spec[j][c]->GetBinContent(k);
pt /= spec[j][c]->Integral(iptmin,iptmax);
double yielderr;
double yield = spec[j][c]->IntegralAndError(iptmin,iptmax,yielderr);
yield /= (ptmax-ptmin) * Nevt * 2 * (etaMax-etaMin);
yielderr /= (ptmax-ptmin) * Nevt * 2 * (etaMax-etaMin);
double efmin = eff[j][c]->GetBinContent(icptmin)/sim[j][c]->GetBinContent(icptmin);
double efmax = eff[j][c]->GetBinContent(icptmax)/sim[j][c]->GetBinContent(icptmax);
double famin = fak[j][c]->GetBinContent(icptmin)/rec[j][c]->GetBinContent(icptmin);
//.........这里部分代码省略.........
示例7: SetStyle
//.........这里部分代码省略.........
TH1F* ttbar = refill((TH1F*)input->Get(TString::Format("%s/TT" , directory)), "TT" ); InitHist(ttbar, "", "", TColor::GetColor(155,152,204), 1001);
TH1F* Ztt = refill((TH1F*)input->Get(TString::Format("%s/ZTT" , directory)), "ZTT"); InitHist(Ztt , "", "", TColor::GetColor(248,206,104), 1001);
#ifdef MSSM
TH1F* ggH = refill((TH1F*)input2->Get(TString::Format("%s/ggH$MA" , directory)), "ggH"); InitSignal(ggH); ggH->Scale($TANB);
TH1F* bbH = refill((TH1F*)input2->Get(TString::Format("%s/bbH$MA" , directory)), "bbH"); InitSignal(bbH); bbH->Scale($TANB);
TH1F* ggH_SM125= refill((TH1F*)input->Get(TString::Format("%s/ggH_SM125" , directory)), "ggH_SM125"); InitHist(ggH_SM125, "", "", kGreen+2, 1001);
TH1F* qqH_SM125= refill((TH1F*)input->Get(TString::Format("%s/qqH_SM125" , directory)), "qqH_SM125"); InitHist(qqH_SM125, "", "", kGreen+2, 1001);
TH1F* VH_SM125 = refill((TH1F*)input->Get(TString::Format("%s/VH_SM125" , directory)), "VH_SM125" ); InitHist(VH_SM125, "", "", kGreen+2, 1001);
#else
#ifndef DROP_SIGNAL
TH1F* ggH = refill((TH1F*)input->Get(TString::Format("%s/ggH125" , directory)), "ggH"); InitSignal(ggH); ggH->Scale(SIGNAL_SCALE);
TH1F* qqH = refill((TH1F*)input->Get(TString::Format("%s/qqH125" , directory)), "qqH"); InitSignal(qqH); qqH->Scale(SIGNAL_SCALE);
TH1F* VH = refill((TH1F*)input->Get(TString::Format("%s/VH125" , directory)), "VH" ); InitSignal(VH ); VH ->Scale(SIGNAL_SCALE);
#endif
#endif
#ifdef ASIMOV
TH1F* data = refill((TH1F*)input->Get(TString::Format("%s/data_obs_asimov", directory)), "data", true);
#else
TH1F* data = refill((TH1F*)input->Get(TString::Format("%s/data_obs", directory)), "data", true);
#endif
InitHist(data, "#bf{m_{#tau#tau} [GeV]}", "#bf{dN/dm_{#tau#tau} [1/GeV]}"); InitData(data);
TH1F* ref=(TH1F*)Fakes->Clone("ref");
ref->Add(EWK0 );
ref->Add(EWK1 );
#ifdef EXTRA_SAMPLES
ref->Add(EWK2 );
#endif
ref->Add(EWK );
ref->Add(ttbar);
ref->Add(Ztt );
double unscaled[7];
unscaled[0] = Fakes->Integral();
unscaled[1] = EWK ->Integral();
unscaled[1]+= EWK0 ->Integral();
unscaled[1]+= EWK1 ->Integral();
#ifdef EXTRA_SAMPLES
unscaled[1]+= EWK2 ->Integral();
#endif
unscaled[2] = ttbar->Integral();
unscaled[3] = Ztt ->Integral();
#ifdef MSSM
unscaled[4] = ggH ->Integral();
unscaled[5] = bbH ->Integral();
unscaled[6] = 0;
#else
#ifndef DROP_SIGNAL
unscaled[4] = ggH ->Integral();
unscaled[5] = qqH ->Integral();
unscaled[6] = VH ->Integral();
#endif
#endif
if(scaled){
/* Fakes = refill(shape_histos(Fakes, datacard, "QCD"), "QCD");
EWK0 = refill(shape_histos(EWK0, datacard, "VV"), "VV");
EWK1 = refill(shape_histos(EWK1, datacard, "W"), "W");
#ifdef EXTRA_SAMPLES
EWK2 = refill(shape_histos(EWK2, datacard, "ZJ"), "ZJ");
EWK = refill(shape_histos(EWK, datacard, "ZL"), "ZL");
#else
// EWK = refill(shape_histos(EWK, datacard, "ZLL"), "ZLL");
#endif
ttbar = refill(shape_histos(ttbar, datacard, "TT"), "TT");示例8: SetStyle
//.........这里部分代码省略.........
TH1F* ggH = refill((TH1F*)input ->Get(TString::Format("%s/ggH125" , directory)), "ggH" );
InitSignal(ggH);
ggH->Scale(SIGNAL_SCALE);
TH1F* qqH = refill((TH1F*)input ->Get(TString::Format("%s/qqH125" , directory)), "qqH" );
InitSignal(qqH);
qqH->Scale(SIGNAL_SCALE);
TH1F* VH = refill((TH1F*)input ->Get(TString::Format("%s/VH125" , directory)), "VH" );
InitSignal(VH );
VH ->Scale(SIGNAL_SCALE);
#endif
#endif
#ifdef ASIMOV
TH1F* data = refill((TH1F*)input->Get(TString::Format("%s/data_obs_asimov", directory)), "data", true);
#else
TH1F* data = refill((TH1F*)input->Get(TString::Format("%s/data_obs", directory)), "data", true);
#endif
#ifdef MSSM
InitHist(data, "#bf{m_{#tau#tau} [GeV]}" , "#bf{dN/dm_{#tau#tau} [1/GeV]}");
InitData(data);
#else
InitHist(data, "#bf{D}", "#bf{dN/dD}" );
InitData(data);
#endif
TH1F* ref=(TH1F*)ZTT->Clone("ref");
ref->Add(ZMM);
ref->Add(TTJ);
ref->Add(QCD);
ref->Add(Dibosons);
if(WJets) {
ref->Add(WJets);
}
double unscaled[9];
unscaled[0] = ZTT->Integral();
unscaled[1] = ZMM->Integral();
unscaled[2] = TTJ->Integral();
unscaled[3] = QCD->Integral();
unscaled[4] = Dibosons->Integral();
unscaled[5] = 0;
if(WJets) {
unscaled[5] = WJets->Integral();
}
#ifdef MSSM
unscaled[6] = ggH->Integral();
unscaled[7] = bbH->Integral();
unscaled[8] = 0;
#else
#ifndef DROP_SIGNAL
unscaled[6] = ggH->Integral();
unscaled[7] = qqH->Integral();
unscaled[8] = VH ->Integral();
#endif
#endif
if(scaled) {
rescale(ZTT, 1);
rescale(ZMM, 2);
rescale(TTJ, 3);
rescale(QCD, 4);
rescale(Dibosons, 5);
if(WJets) {
rescale(WJets, 6);
}
#ifdef MSSM
rescale(ggH, 7);
rescale(bbH, 8);示例9: main
int main(int argc, char **argv) {
TFile *fshape = new TFile("mlb.root");
TH1F *shapemlb = (TH1F *) fshape->Get("mlb");
MassReconstructor theMass(100, shapemlb);
TFile *f = new TFile(argv[1]);
f->cd();
TTree *tree = (TTree *)f->Get("events");
setSetBranchAddresses(tree);
Long64_t nentries = tree->GetEntriesFast();
Long64_t nb = 0;
//nentries = 5000;
char w_char[50];
sprintf(w_char, "%s_w", argv[2]);
char pt_char[50];
sprintf(pt_char, "%s_pt", argv[2]);
char comp_char[50];
sprintf(comp_char, "%s_comp", argv[2]);
TH1F *weights = new TH1F(w_char, "", 50, 0, 0.01);
TH1F *pt_med = new TH1F(pt_char, "", 50, 0, 1000);
TH1F *pt_med_meas = new TH1F("pt_med_meas", "", 50, 0, 1000);
TH1F *pt_comp = new TH1F(comp_char, "", 50, -1, 1);
TH2F *med_comp = new TH2F("med_res", "", 100, -1000, 1000, 100, -1000, 1000);
nentries = 5000;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
std::cout << "Starting event" << jentry << std::endl;
nb = tree->GetEntry(jentry);
//Loading into TLorentzVector
TLorentzVector med; med.SetPtEtaPhiM(ev.ptMED, ev.etaMED, ev.phiMED, ev.mMED);
TLorentzVector t1; t1.SetPtEtaPhiM(ev.ptTOP1, ev.etaTOP1, ev.phiTOP1, ev.mTOP1);
TLorentzVector t2; t2.SetPtEtaPhiM(ev.ptTOP2, ev.etaTOP2, ev.phiTOP2, ev.mTOP2);
TLorentzVector l1; l1.SetPtEtaPhiM(ev.ptlep1, ev.etalep1, ev.philep1, ev.mlep1);
TLorentzVector l2; l2.SetPtEtaPhiM(ev.ptlep2, ev.etalep2, ev.philep2, ev.mlep2);
TLorentzVector j1; j1.SetPtEtaPhiM(ev.ptb1, ev.etab1, ev.phib1, ev.mb1);
TLorentzVector j2; j2.SetPtEtaPhiM(ev.ptb2, ev.etab2, ev.phib2, ev.mb2);
TLorentzVector n1; n1.SetPtEtaPhiM(ev.ptnu1, ev.etanu1, ev.phinu1, ev.mnu1);
TLorentzVector n2; n2.SetPtEtaPhiM(ev.ptnu2, ev.etanu2, ev.phinu2, ev.mnu2);
//Collections of jets
std::vector<TLorentzVector> jets;
std::vector<TLorentzVector> bjets; bjets.push_back(j1); bjets.push_back(j2);
//MET
TVector2 origMET(ev.ptMET*cos(ev.phiMET), ev.ptMET*sin(ev.phiMET));
TVector2 MET(ev.ptMET*cos(ev.phiMET), ev.ptMET*sin(ev.phiMET));
//The tops
TVector2 top1, top2;
//And the weight
Float_t w = 0;
Float_t step = 0.1;
Float_t lambda = 1;
for(unsigned int jtrial = 0; jtrial < 1000; ++jtrial) {
theMass.startVariations(l1, l2, bjets, jets, MET, top1, top2, w);
if(w > 0) {
std::cout << "Solution was found" << std::endl;
if(jtrial == 0) break;
weights->Fill(w);
TVector2 ptmed(med.Pt()*cos(med.Phi()), med.Pt()*sin(med.Phi()));
TVector2 ptmed_meas = -1.0 * (top1 + top2);
pt_med->Fill(ptmed.Mod());
pt_med_meas->Fill(ptmed_meas.Mod());
pt_comp->Fill((top1.Mod() - ev.ptTOP1)/ev.ptTOP1);
TVector2 med2 = origMET - MET;
med_comp->Fill(ptmed.X(), med2.X());
med_comp->Fill(ptmed.Y(), med2.Y());
break;
}
vector<double> grad;
theMass.massSolver->gradient(MET, j1 , j2 , l1, l2, 80.385, 80.385, 173.34, 173.34, grad, step, lambda);
TVector2 theGrad(grad[0], grad[1]);
MET = MET + theGrad;
}
}
weights->Scale(1.0/weights->Integral());
pt_med->Scale(1.0/pt_med->Integral());
pt_comp->Scale(1.0/pt_comp->Integral());
TFile *foutput = new TFile("output.root", "UPDATE");
foutput->cd();
weights->Write();
pt_med_meas->Write();
pt_med->Write();
pt_comp->Write();
med_comp->Write();
foutput->Close();
f->Close();
fshape->Close();
return 1;
//.........这里部分代码省略.........
示例10: syst_MCeta
//.........这里部分代码省略.........
cout<<"ERROR !! The syst2 plot is void. Exiting now ..."<<endl;
return;
}
TH1F* factor = new TH1F("factor","factor",hbase->GetNbinsX(),hbase->GetXaxis()->GetXbins()->GetArray());
if(E==0.9){
for(int i=1;i<=hbase->GetNbinsX();++i){
if(hsyst1->GetBinContent(i)!=0)
factor->SetBinContent(i,hsyst2->GetBinContent(i) / hsyst1->GetBinContent(i));
else
factor->SetBinContent(i,1.);
}
}
else if(E==2.36){
double start = -0.5;
double ratio = 1;
if((getEdgeLastFilledBin(hsyst1) - start)!=0)
ratio = (getEdgeLastFilledBin(hbase) - start) / (getEdgeLastFilledBin(hsyst1) - start);
//part with no stretching
for( int i = 1 ; i < hbase->GetXaxis()->FindFixBin(start) ; ++i ){
if(hsyst1->GetBinContent(i)!=0)
factor->SetBinContent(i,hsyst2->GetBinContent(i) / hsyst1->GetBinContent(i));
else
factor->SetBinContent(i,1.);
}
//part with stretching
for( int i = hbase->GetXaxis()->FindFixBin(start) ; i <= hbase->GetNbinsX() ; ++i ){
int bin = hsyst1->GetXaxis()->FindFixBin( (hbase->GetBinCenter(i)-start) / ratio + start);
if(hsyst1->GetBinContent(bin)!=0)
factor->SetBinContent(i,hsyst2->GetBinContent(bin) / hsyst1->GetBinContent(bin));
else
factor->SetBinContent(i,1.);
}
}
TH1F* hsyst = (TH1F*) hbase->Clone("eta_syst");
hsyst->Multiply(factor);
TMoments* msyst = new TMoments(hsyst);
msyst->ComputeMoments();
msyst->print();
//Making output file
outstr.str("");
outstr << "hyp" << 1 << "_niter" << 0 << "_cut" << cut << "_DataType" << 0;
TString toutput = fileManager(3,typeMC,E,1,500,-1,outstr.str());
TFile* foutput = TFile::Open(toutput,"RECREATE");
cout<<"Output file : " << toutput << endl;
//Making kno plot
double knomean = hsyst->GetMean();
TString tkno = toutput;
tkno.ReplaceAll("plots/","plots/current_b1_2/");
cout<<"Opening for kno mean the file : "<<tkno<<endl;
TFile* fkno = TFile::Open(tkno,"READ");
if(fkno!=0){
TMoments* mom_kno = (TMoments*) fkno->Get("unfolding/moments/moments");
knomean = mom_kno->mean->GetMean();
}
else{
cout<<"WARNING !! The file does not exist, taking the mean of the hist instead"<<endl;
}
TH1F* kno = new TH1F("kno_corrected","kno_corrected;z = n_{ch} / < n_{ch} >;#psi(z)",hsyst->GetNbinsX(),Divide(hsyst->GetXaxis()->GetXbins(),knomean));
kno->Sumw2();
/*for( int k = 60 ; k <= nch_corrected->GetNbinsX() ; ++k)
nch_corrected->SetBinContent(k,0);*/
for( int k = 1 ; k <= hsyst->GetNbinsX() ; ++k){
kno->SetBinContent(k , knomean * hsyst->GetBinContent(k) / hsyst->Integral());
kno->SetBinError(k , knomean * hsyst->GetBinError(k) / hsyst->Integral());
}
foutput->cd();
foutput->mkdir("unfolding");
foutput->cd("unfolding");
factor->Write();
hsyst->Write("nch_data_corrected");
kno->Write();
gDirectory->mkdir("moments");
gDirectory->cd("moments");
msyst->Write("moments");
foutput->Close();
}示例11: fitTF1
void fitTF1(TCanvas *canvas, TH1F h, double XMIN_, double XMAX_, double dX_, double params[], Color_t LC_=kBlack) { //double& FWHM_, double& x0_, double& x1_, double& x2_, double& y0_, double& y1_, double& INT_, double& YIELD_) {
//TCanvas* fitTF1(TH1F h, double HSCALE_, double XMIN_, double XMAX_) {
//TF1* fitTF1(TH1F h, double HSCALE_, double XMIN_, double XMAX_) {
gROOT->ForceStyle();
RooMsgService::instance().setSilentMode(kTRUE);
for(int i=0;i<2;i++) RooMsgService::instance().setStreamStatus(i,kFALSE);
//TCanvas *canvas = new TCanvas(TString::Format("canvas_%s",h.GetName()),TString::Format("canvas_%s",h.GetName()),1800,1200);
RooDataHist *RooHistFit = 0;
RooAddPdf *model = 0;
RooWorkspace w = RooWorkspace("w","workspace");
RooRealVar x("mbbReg","mbbReg",XMIN_,XMAX_);
RooRealVar kJES("CMS_scale_j","CMS_scale_j",1,0.9,1.1);
RooRealVar kJER("CMS_res_j","CMS_res_j",1,0.8,1.2);
kJES.setConstant(kTRUE);
kJER.setConstant(kTRUE);
TString hname = h.GetName();
RooHistFit = new RooDataHist("fit_"+hname,"fit_"+hname,x,&h);
RooRealVar YieldVBF = RooRealVar("yield_"+hname,"yield_"+hname,h.Integral());
RooRealVar m("mean_"+hname,"mean_"+hname,125,100,150);
RooRealVar s("sigma_"+hname,"sigma_"+hname,12,3,30);
RooFormulaVar mShift("mShift_"+hname,"@0*@1",RooArgList(m,kJES));
RooFormulaVar sShift("sShift_"+hname,"@0*@1",RooArgList(m,kJER));
RooRealVar a("alpha_"+hname,"alpha_"+hname,1,-10,10);
RooRealVar n("exp_"+hname,"exp_"+hname,1,0,100);
RooRealVar b0("b0_"+hname,"b0_"+hname,0.5,0.,1.);
RooRealVar b1("b1_"+hname,"b1_"+hname,0.5,0.,1.);
RooRealVar b2("b2_"+hname,"b2_"+hname,0.5,0.,1.);
RooRealVar b3("b3_"+hname,"b3_"+hname,0.5,0.,1.);
RooBernstein bkg("signal_bkg_"+hname,"signal_bkg_"+hname,x,RooArgSet(b0,b1,b2));
RooRealVar fsig("fsig_"+hname,"fsig_"+hname,0.7,0.0,1.0);
RooCBShape sig("signal_gauss_"+hname,"signal_gauss_"+hname,x,mShift,sShift,a,n);
model = new RooAddPdf("signal_model_"+hname,"signal_model_"+hname,RooArgList(sig,bkg),fsig);
//RooFitResult *res = model->fitTo(*RooHistFit,RooFit::Save(),RooFit::SumW2Error(kFALSE),"q");
model->fitTo(*RooHistFit,RooFit::Save(),RooFit::SumW2Error(kFALSE),"q");
//res->Print();
//model->Print();
canvas->cd();
canvas->SetTopMargin(0.1);
RooPlot *frame = x.frame();
// no scale
RooHistFit->plotOn(frame);
model->plotOn(frame,RooFit::LineColor(LC_),RooFit::LineWidth(2));//,RooFit::LineStyle(kDotted));
model->plotOn(frame,RooFit::Components(bkg),RooFit::LineColor(LC_),RooFit::LineWidth(2),RooFit::LineStyle(kDashed));
// with scale
// RooHistFit->plotOn(frame,RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent));
// model->plotOn(frame,RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent),RooFit::LineWidth(1));
// model->plotOn(frame,RooFit::Components(bkg),RooFit::LineColor(kBlue),RooFit::LineWidth(1),RooFit::LineStyle(kDashed),RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent));
frame->GetXaxis()->SetLimits(50,200);
frame->GetXaxis()->SetNdivisions(505);
frame->GetXaxis()->SetTitle("M_{b#bar{b}} (GeV)");
frame->GetYaxis()->SetTitle("Events");
frame->Draw();
h.SetFillColor(kGray);
h.Draw("hist,same");
frame->Draw("same");
gPad->RedrawAxis();
TF1 *tmp = model->asTF(x,fsig,x);
//tmp->Print();
double y0_ = tmp->GetMaximum();
double x0_ = tmp->GetMaximumX();
double x1_ = tmp->GetX(y0_/2.,XMIN_,x0_);
double x2_ = tmp->GetX(y0_/2.,x0_,XMAX_);
double FWHM_ = x2_-x1_;
double INT_ = tmp->Integral(XMIN_,XMAX_);
double YIELD_= YieldVBF.getVal();
double y1_ = dX_*0.5*y0_*(YieldVBF.getVal()/tmp->Integral(XMIN_,XMAX_));
params[0] = x0_;
params[1] = x1_;
params[2] = x2_;
params[3] = y0_;
params[4] = y1_;
params[5] = FWHM_;
params[6] = INT_;
params[7] = YIELD_;
//cout<<"Int = "<<tmp->Integral(XMIN_,XMAX_)<<", Yield = "<<YieldVBF.getVal()<<", y0 = "<<y0_<<", y1 = "<<y1_ <<", x0 = "<< x0_ << ", x1 = "<<x1_<<", x2 = "<<x2_<<", FWHM = "<<FWHM_<<endl;
TLine ln = TLine(x1_,y1_,x2_,y1_);
ln.SetLineColor(kMagenta+3);
ln.SetLineStyle(7);
ln.SetLineWidth(2);
ln.Draw();
//.........这里部分代码省略.........
示例12: GenerateSusyFile
//.........这里部分代码省略.........
sprintf( hname, "h_susy_sl_%db", bi+1 ) ;
h_susy_sl[bi] = new TH2F( hname, hname, nBinsVar1, Mbins, nBinsVar2, Hbins ) ;
h_susy_sl[bi] -> Sumw2() ;
sprintf( hname, "h_susy_ldp_%db", bi+1 ) ;
h_susy_ldp[bi] = new TH2F( hname, hname, nBinsVar1, Mbins, nBinsVar2, Hbins ) ;
h_susy_ldp[bi] -> Sumw2() ;
}
// the following is just an example, to run on a bunch of points
// in the T2tt plane
int mGls[9] = {350,450,500,550,600,650,700,750,800};
int mLsps[9] = {0,0,0,0,0,0,0,0,0};
for ( int iGl = 0 ; iGl < 9 ; iGl++ ) {
int mGl = mGls[iGl] ;
int mLsp = mLsps[iGl] ;
TString cutSMS = "mgluino>";
cutSMS += mGl-1;
cutSMS += "&&mgluino<";
cutSMS += mGl+1;
cutSMS += "&&mlsp>";
cutSMS += mLsp-1;
cutSMS += "&&mlsp<";
cutSMS += mLsp+1;
TH1F *dummyHist = new TH1F("dummyHist","",100,0.,10000.);
chainT2tt.Project("dummyHist","MET",cutSMS);
inFile << mGl << " " << mLsp << " " << dummyHist->Integral() << " " ;
printf(" mGl=%4d, mLsp=%4d\n", mGl, mLsp ) ; cout << flush ;
cutSMS += "&&";
printf("\n\n") ;
for (int k = 0 ; k < nBinsBjets ; k++) {
TString cut = cBbins[k] ;
TString allSigCuts = cutSMS+cutsSig+cut ;
TString allSLSigCuts = cutSMS+cutsSLSig+cut ;
TString allSLCuts = cutSMS+cutsSL+cut ;
TString allLDPCuts = cutSMS+cutsLDP+cut ;
char hname[1000] ;
sprintf( hname, "h_susy_sig_%db", k+1 ) ;
chainT2tt.Project( hname,s2DVars,allSigCuts);
printf(" mGl=%d, mLsp=%d, nBjets = %d, SIG selection %9.1f events.\n", mGl, mLsp, k+1, h_susy_sig[k]->Integral() ) ; cout << flush ;
sprintf( hname, "h_susy_slsig_%db", k+1 ) ;
chainT2tt.Project( hname,s2DVars,allSLSigCuts);
printf(" mGl=%d, mLsp=%d, nBjets = %d, SL Sig selection %6.1f events.\n", mGl, mLsp, k+1, h_susy_slsig[k]->Integral() ) ; cout << flush ;
sprintf( hname, "h_susy_sl_%db", k+1 ) ;
chainT2tt.Project( hname,s2DVars,allSLCuts);
printf(" mGl=%d, mLsp=%d, nBjets = %d, SL selection %9.1f events.\n", mGl, mLsp, k+1, h_susy_sl[k]->Integral() ) ; cout << flush ;
sprintf( hname, "h_susy_ldp_%db", k+1 ) ;
chainT2tt.Project( hname,s2DVars,allLDPCuts);
printf(" mGl=%d, mLsp=%d, nBjets = %d, LDP selection %9.1f events.\n", mGl, mLsp, k+1, h_susy_ldp[k]->Integral() ) ; cout << flush ;示例13: calccorr_sandv
//.........这里部分代码省略.........
for(int icent_b=0; icent_b<ncent+1; icent_b++)
if(selcentbin[icent_a] == centbin[icent_b]) break;
int xcentmin = icent_b;
for(int icent_b=0; icent_b<ncent+1; icent_b++)
if(selcentbin[icent_a+1] == centbin[icent_b]) break;
int xcentmax = icent_b;
if((xcentmin == ncent+1) || (xcentmax == ncent+1)) exit(0);
for(int ipt_a=0;ipt_a<npt_a;ipt_a++){
for(int ipt_b=0; ipt_b<npt+1; ipt_b++)
if(selptbin[ipt_a] == ptbin[ipt_b]) break;
int xptmin = ipt_b;
for(int ipt_b=0; ipt_b<npt+1; ipt_b++)
if(selptbin[ipt_a+1] == ptbin[ipt_b]) break;
int xptmax = ipt_b;
if((xptmin == npt+1) || (xptmax == npt+1)) exit(0);
cout<<xcentmin<<"\t"<<xcentmax<<"\t"<<endl;
cout<<xptmin<<"\t"<<xptmax<<"\t"<<endl;
kforebbcw_In = (TH1F*)kforebbcw[xcentmin][xptmin]->Clone();
hforebbcw_In = (TH1F*)hforebbcw[xcentmin][xptmin]->Clone();
kbackbbcw2_In = (TH1F*)kbackbbcw2[xcentmin][xptmin]->Clone();
kforebbcw_In->Reset();
hforebbcw_In->Reset();
kbackbbcw2_In->Reset();
for(int icent=xcentmin; icent<xcentmax; icent++){
for(int ipt=xptmin; ipt<xptmax; ipt++){
kforebbcw_In->Add(kforebbcw[icent][ipt]);
hforebbcw_In->Add(hforebbcw[icent][ipt]);
kbackbbcw2_In->Add(kbackbbcw2[icent][ipt]);
}
}
kforebbcw_In->Rebin(2);
hforebbcw_In->Rebin(2);
kbackbbcw2_In->Rebin(2);
TH1F* hpp;
TH1F* hbackpp;
hpp = (TH1F*)kforebbcw_In->Clone();
hbackpp = (TH1F*)kbackbbcw2_In->Clone();
float nbackpp = hbackpp->Integral()/2.0/PI;
float nforepp = hpp->Integral()/2.0/PI;
//float ntrig0 = ptforedis_0->Integral(11,30);
for(int i=0; i<20; i++){
float pp_cont = 1.0*hpp->GetBinContent(i+1);
//float pp0_err = 1.0*hpp->GetBinError(i+1);
float weight2 = sqrt(1.0*hforebbcw_In->GetBinContent(i+1));
float backpp_cont = 1.0*hbackpp->GetBinContent(i+1);
float con = pp_cont/backpp_cont*nbackpp/nforepp;
float err = weight2/backpp_cont*nbackpp/nforepp;
hpp->SetBinContent(i+1, con);
hpp->SetBinError(i+1, err);
if(i==0) fsout<<pp_cont<<" "<<weight2<<" "<<con<<" "<<err<<endl;
}
TF1 *fun0 = new TF1("fun0","[0]*(1+2*[1]*cos(x)+2*[2]*cos(2*x)+2*[3]*cos(3*x)+2*[4]*cos(4*x))", -0.5*PI, 1.5*PI);
fun0->SetLineColor(1);
hpp->Fit("fun0","NORQ");
TF1 *fun1 = new TF1("fun1","[0]*(1+2*[1]*cos(x))", -0.5*PI, 1.5*PI);
TF1 *fun2 = new TF1("fun2","[0]*(1+2*[1]*cos(2*x))", -0.5*PI, 1.5*PI);
TF1 *fun3 = new TF1("fun3","[0]*(1+2*[1]*cos(3*x))", -0.5*PI, 1.5*PI);
TF1 *fun4 = new TF1("fun4","[0]*(1+2*[1]*cos(4*x))", -0.5*PI, 1.5*PI);
fout<<fun0->GetParameter(1)<<" "<<fun0->GetParError(1)<<" "
// <<-fun0->GetParameter(2)/fun0->GetParameter(1)<<" "<<fun0->GetParameter(2)/fun0->GetParameter(1)*(sqrt((fun0->GetParError(2)/fun0->GetParameter(2))**2+(fun0->GetParError(1)/fun0->GetParameter(1))**2))<<" "
<<fun0->GetParameter(2)<<" "<<fun0->GetParError(2)<<" "
<<fun0->GetParameter(3)<<" "<<fun0->GetParError(3)<<endl;
}
}
fout.close();
fsout.close();
/*
cout<<"*************** v2 ***********"<<endl;
float v2_0 = funvn0->GetParameter(1)/(zym_pp0 + funvn0->GetParameter(0));
float v2_1 = funvn1->GetParameter(1)/(zym_pp1 + funvn1->GetParameter(0));
float v2_2 = funvn2->GetParameter(1)/(zym_pp2 + funvn2->GetParameter(0));
float v2_3 = funvn3->GetParameter(1)/(zym_pp3 + funvn3->GetParameter(0));
cout<<v2_0<<" "<<v2_1<<" "<<v2_2<<" "<<v2_3<<endl;
cout<<funvn0->GetParameter(0)*funvn0->GetParameter(1)<<" "
<<funvn1->GetParameter(0)*funvn1->GetParameter(1)<<" "
<<funvn2->GetParameter(0)*funvn2->GetParameter(1)<<" "
<<funvn3->GetParameter(0)*funvn3->GetParameter(1)<<endl;
cout<<"*************** v2 ***********"<<endl;
cout<<funvn0->GetParameter(2)<<" "<<funvn1->GetParameter(2)<<" "<<funvn2->GetParameter(2)<<" "<<funvn3->GetParameter(2)<<endl;
cout<<"*************** v3 ***********"<<endl;
cout<<funvn0->GetParameter(3)<<" "<<funvn1->GetParameter(3)<<" "<<funvn2->GetParameter(3)<<" "<<funvn3->GetParameter(3)<<endl;
*/
}示例14: dataCardSM
//.........这里部分代码省略.........
TH1F* VV = 0;
if(sm==0 || sm==2 || sm==1 || sm==3) VV=analysis->getDiBoson();
if(sm==4) VV=analysis->getDiBosonVBFHCP();
VV->SetName("VV");
TH1F* ZLL=(TH1F*)ZL->Clone("ZLL");
ZLL->SetName("ZLL");
ZLL->Add(ZJ);
//blind
TString tmpsel=analysis->extrasel_;
if(Blind)analysis->extrasel_ += "*(svfitmass<100||160<svfitmass)";
TH1F* data_obs = analysis->getTotalData();
data_obs->SetName("data_obs");
analysis->extrasel_ =tmpsel;
TH1F* MC=(TH1F*)ZTT->Clone("MC");//needed below
MC->Add(ZL);
MC->Add(ZJ);
MC->Add(W);
MC->Add(TT);
MC->Add(VV);
MC->Add(QCD);
dir->cd();
fix0Bins(ZTT); ZTT->Write();
ZTT->SetName("ZTT125"); ZTT->Write();//needed for ZTT fits
fix0Bins(ZL); ZL->Write();
fix0Bins(ZJ); ZJ->Write();
fix0Bins(ZLL); ZLL->Write();
fix0Bins(W); W->Write();
fix0Bins(TT); TT->Write();
fix0Bins(VV); VV->Write();
fix0Bins(QCD); QCD->Write();
data_obs->Write();
gROOT->cd();
delete ZTT ;
delete ZL;
delete ZJ;
delete ZLL;
delete W;
delete TT;
delete VV;
delete QCD;
for(Int_t m=0;m<NMASS;m++){
long ma=massValues[m];
//Nominal h
TH1F* SM = analysis->getSample(TString("HiggsGG")+ma);
SM->SetName(TString("ggH")+ma);
TH1F* VBF = analysis->getSample(TString("HiggsVBF")+ma);
VBF->SetName(TString("qqH")+ma);
TH1F* VH = analysis->getSample(TString("HiggsVH")+ma);
VH->SetName(TString("VH")+ma);
SM->Scale(1./analysis->findSample(TString("HiggsGG")+ma)->getCrossection());
VBF->Scale(1./analysis->findSample(TString("HiggsVBF")+ma)->getCrossection());
VH->Scale(1./analysis->findSample(TString("HiggsVH")+ma)->getCrossection());
//check for empty histos
if( SM->Integral()<=0.){ SM->SetBinContent(SM->GetNbinsX()/2,1e-4); SM->SetBinError(SM->GetNbinsX()/2,1e-4); }
if( VBF->Integral()<=0.){ VBF->SetBinContent(VBF->GetNbinsX()/2,1e-4); VBF->SetBinError(VBF->GetNbinsX()/2,1e-4); }
if( VH->Integral()<=0.){ VH->SetBinContent(VH->GetNbinsX()/2,1e-4); VH->SetBinError(VH->GetNbinsX()/2,1e-4); }
dir->cd();
fixSignal(data_obs,MC,VH); fix0Bins(VH); VH->Write();
fixSignal(data_obs,MC,SM); fix0Bins(SM); SM->Write();
fixSignal(data_obs,MC,VBF); fix0Bins(VBF); VBF->Write();
gROOT->cd();
delete VH;
delete SM;
delete VBF;
}
delete MC;
delete data_obs;
}
output.ls();
output.Close();
gROOT->ProcessLine(".q");
}示例15: IndResponse
//.........这里部分代码省略.........
//! 2D correlation between refpt and recopt
hcorrptrefpt[nj][centb]->Fill(refpt,recopt,wxs*wcen*wvz);
hrawptrefpt [nj][centb]->Fill(refpt,rawpt,wxs*wcen*wvz);
hcorrptrawpt[nj][centb]->Fill(rawpt,recopt,wxs*wcen*wvz);
//! Very fine bin in ref pt
hrescrpt[nj][centb]->Fill(refpt,recopt/refpt,wxs*wcen*wvz);
hresrrpt[nj][centb]->Fill(refpt,rawpt/refpt,wxs*wcen*wvz);
hresrcrpt[nj][centb]->Fill(recopt,recopt/rawpt,wxs*wcen*wvz);
int ieta=-1;
if(fabs(recoeta)<1.3)ieta=0; //! barrel region
else ieta=1; //! HCAL region
if(ieta>=0){
hratiocorrrefpt_eta[nj][centb][ieta]->Fill(refpt,recopt/refpt,wxs*wcen*wvz);
hratiorawrefpt_eta [nj][centb][ieta]->Fill(refpt,rawpt/refpt,wxs*wcen*wvz);
}
//! pileup study
hjetptpu[nj][centb]->Fill(recopt,iJet->jtpu[ij],wxs*wcen*wvz);
//! Response in different eta and phi bins
int etabin = GetEtaBin(fabs(refeta));
int phibin = GetPhiBin(refphi);
if(etabin < 0 || etabin>=maxe || phibin < 0 || phibin>=maxph)continue;
//! Response in eta and phi bins
hpteta[nj][centb][etabin]->Fill(refpt,recopt/refpt,wxs*wcen*wvz);
hptphi[nj][centb][phibin]->Fill(refpt,recopt/refpt,wxs*wcen*wvz);
//! Gen:Reco correlation plots
hgenjrecoj [nj][centb]->Fill(refpt,recopt,wxs*wcen*wvz);
if(fabs(refparton_flavor)<=21)hgenjrecoj_pflavor [nj][centb]->Fill(refpt,recopt,wxs*wcen*wvz);
}//! ij loop
hNevt[nj][centb]->Fill(vz);
if(njets>0){
hAvjets[nj][centb]->Fill(pthat,njets,wxs*wcen*wvz);
hMeanPtjets[nj][centb]->Fill(pthat,meanpt/njets,wxs*wcen*wvz);
}
delete [] ljet;
}//! nj jet loop
if(istat){
hBin->Fill(hiBin,wxs*wcen*wvz);
hVz->Fill(vz,wxs*wcen*wvz);
hHF->Fill(hiHF,wxs*wcen*wvz);
iEvent++;
}
//std::cout<<"Completed event # "<<ievt<<std::endl;
}//! event loop ends
std::cout<<std::endl;
std::cout<<std::endl;
std::cout<<std::endl;
std::cout<<"Events which passed the pT hat cut : "<<hTotEve->Integral()<<" out of : "<<hEvt->Integral()
<<" efficiency of all the cuts : " <<hTotEve->Integral()/hEvt->Integral()<<std::endl;
std::cout<<std::endl;
if(strcmp(ksp,"pp")==0){
for(int nj=0;nj<knj;nj++){
//std::cout<<"# of Events for : "<<c->GetCjets[Nj](nj)<<"\t"<<hNevt[nj][0]->Integral()<<"\t # of Jets : "<<hNjets[nj][0]->Integral()<<std::endl;
std::cout<<"# of Events for : "<<cjets[nj]<<"\t"<<hNevt[nj][0]->Integral()<<"\t # of Jets : "<<hNjets[nj][0]->Integral()<<std::endl;
}
}else{
for(int nj=0;nj<knj;nj++){
for(int icen=0;icen<ncen;icen++){
std::cout<<"# of Events for : "<<cjets[nj]<<"\t icen : "<<icen<<"\t"<<hNevt[nj][icen]->Integral()<<"\t # of Jets : "<<hNjets[nj][icen]->Integral()<<std::endl;
}
std::cout<<std::endl;
}
}
//! Write to output file
fout->cd();
fout->Write();
fout->Close();
//! Check
timer.Stop();
float mbytes = 0.000001*nb;
double rtime = timer.RealTime();
double ctime = timer.CpuTime();
std::cout<<std::endl;
std::cout<<Form("RealTime=%f seconds, CpuTime=%f seconds",rtime,ctime)<<std::endl;
std::cout<<Form("You read %f Mbytes/RealTime seconds",mbytes/rtime)<<std::endl;
std::cout<<Form("You read %f Mbytes/CpuTime seconds",mbytes/ctime)<<std::endl;
std::cout<<std::endl;
std::cout<<"Good bye : " <<"\t"<<std::endl;
return 1;
}