C++ TPaveText::Draw方法代码示例

//.........这里部分代码省略.........
   SegtoLCTXY_1__85->SetBinContent(83,0.003510582);
   SegtoLCTXY_1__85->SetBinContent(84,0.00228951);
   SegtoLCTXY_1__85->SetBinContent(85,0.003663216);
   SegtoLCTXY_1__85->SetBinContent(86,0.00305268);
   SegtoLCTXY_1__85->SetBinContent(87,0.00305268);
   SegtoLCTXY_1__85->SetBinContent(88,0.003968484);
   SegtoLCTXY_1__85->SetBinContent(89,0.002594778);
   SegtoLCTXY_1__85->SetBinContent(90,0.002747412);
   SegtoLCTXY_1__85->SetBinContent(91,0.003510582);
   SegtoLCTXY_1__85->SetBinContent(92,0.00457902);
   SegtoLCTXY_1__85->SetBinContent(93,0.003968484);
   SegtoLCTXY_1__85->SetBinContent(94,0.003205314);
   SegtoLCTXY_1__85->SetBinContent(95,0.002442144);
   SegtoLCTXY_1__85->SetBinContent(96,0.00152634);
   SegtoLCTXY_1__85->SetBinContent(97,0.003510582);
   SegtoLCTXY_1__85->SetBinContent(98,0.00228951);
   SegtoLCTXY_1__85->SetBinContent(99,0.002900046);
   SegtoLCTXY_1__85->SetBinContent(100,0.002594778);
   SegtoLCTXY_1__85->SetBinContent(101,3.150366);
   SegtoLCTXY_1__85->SetMinimum(0.02);
   SegtoLCTXY_1__85->SetMaximum(50);
   SegtoLCTXY_1__85->SetEntries(675802);
   SegtoLCTXY_1__85->SetStats(0);

   Int_t ci;      // for color index setting
   TColor *color; // for color definition with alpha
   ci = TColor::GetColor("#ff00ff");
   SegtoLCTXY_1__85->SetLineColor(ci);

   ci = TColor::GetColor("#ff00ff");
   SegtoLCTXY_1__85->SetMarkerColor(ci);
   SegtoLCTXY_1__85->GetXaxis()->SetTitle("cm");
   SegtoLCTXY_1__85->GetYaxis()->SetTitle("scaled number of entries");
   SegtoLCTXY_1__85->Draw("H");
   
   TLegend *leg = new TLegend(0.5,0.7,0.9,0.9,NULL,"brNDC");
   leg->SetBorderSize(1);
   leg->SetTextFont(62);
   leg->SetTextSize(0.02);
   leg->SetLineColor(1);
   leg->SetLineStyle(1);
   leg->SetLineWidth(1);
   leg->SetFillColor(0);
   leg->SetFillStyle(1001);
   TLegendEntry *entry=leg->AddEntry("SegtoLCTXY_1","ME11A: mean:0.4cm;RMS:1.9cm","l");

   ci = TColor::GetColor("#ff00ff");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
   entry->SetMarkerColor(1);
   entry->SetMarkerStyle(21);
   entry->SetMarkerSize(1);
   entry->SetTextFont(62);
   entry=leg->AddEntry("SegtoLCTXY_2","ME11B: mean:0.4cm;RMS:1.8cm","l");

   ci = TColor::GetColor("#ff9999");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
   entry->SetMarkerColor(1);
   entry->SetMarkerStyle(21);
   entry->SetMarkerSize(1);
   entry->SetTextFont(62);
   entry=leg->AddEntry("SegtoLCTXY_3","ME12+13: mean:0.8cm;RMS:1.8cm","l");
   entry->SetLineColor(1);

void SegX_2016B_June22all_sameYrange_fullIntegral_10k9k()
{
//=========Macro generated from canvas: SegX/SegX
//=========  (Fri Aug  5 09:22:37 2016) by ROOT version6.06/01
   TCanvas *SegX = new TCanvas("SegX", "SegX",0,0,500,500);
   gStyle->SetOptStat(0);
   SegX->SetHighLightColor(2);
   SegX->Range(-125,-2.848556,125,2.626708);
   SegX->SetFillColor(0);
   SegX->SetBorderMode(0);
   SegX->SetBorderSize(2);
   SegX->SetLogy();
   SegX->SetFrameBorderMode(0);
   SegX->SetFrameBorderMode(0);
   
   TH1D *SegX_1__19 = new TH1D("SegX_1__19","SegX",100,-100,100);
   SegX_1__19->SetBinContent(0,399.7427);
   SegX_1__19->SetBinContent(39,0.0007394809);
   SegX_1__19->SetBinContent(40,0.01109221);
   SegX_1__19->SetBinContent(41,0.02366339);
   SegX_1__19->SetBinContent(42,0.04880574);
   SegX_1__19->SetBinContent(43,0.09613251);
   SegX_1__19->SetBinContent(44,0.5671818);
   SegX_1__19->SetBinContent(45,1.913037);
   SegX_1__19->SetBinContent(46,3.895585);
   SegX_1__19->SetBinContent(47,8.348);
   SegX_1__19->SetBinContent(48,11.60985);
   SegX_1__19->SetBinContent(49,11.7363);
   SegX_1__19->SetBinContent(50,12.04171);
   SegX_1__19->SetBinContent(51,12.05058);
   SegX_1__19->SetBinContent(52,11.89159);
   SegX_1__19->SetBinContent(53,11.26081);
   SegX_1__19->SetBinContent(54,8.149079);
   SegX_1__19->SetBinContent(55,3.796495);
   SegX_1__19->SetBinContent(56,1.839089);
   SegX_1__19->SetBinContent(57,0.5302078);
   SegX_1__19->SetBinContent(58,0.09835096);
   SegX_1__19->SetBinContent(59,0.05472159);
   SegX_1__19->SetBinContent(60,0.02292391);
   SegX_1__19->SetBinContent(61,0.008873771);
   SegX_1__19->SetBinContent(62,0.005176366);
   SegX_1__19->SetMinimum(0.005);
   SegX_1__19->SetMaximum(120);
   SegX_1__19->SetEntries(675802);
   SegX_1__19->SetStats(0);

   Int_t ci;      // for color index setting
   TColor *color; // for color definition with alpha
   ci = TColor::GetColor("#ff00ff");
   SegX_1__19->SetLineColor(ci);

   ci = TColor::GetColor("#ff00ff");
   SegX_1__19->SetMarkerColor(ci);
   SegX_1__19->GetXaxis()->SetTitle("cm");
   SegX_1__19->GetYaxis()->SetTitle("scaled number of entries");
   SegX_1__19->Draw("H");
   
   TLegend *leg = new TLegend(0.5,0.7,0.9,0.9,NULL,"brNDC");
   leg->SetBorderSize(1);
   leg->SetTextFont(62);
   leg->SetTextSize(0.02);
   leg->SetLineColor(1);
   leg->SetLineStyle(1);
   leg->SetLineWidth(1);
   leg->SetFillColor(0);
   leg->SetFillStyle(1001);
   TLegendEntry *entry=leg->AddEntry("SegX_1","ME11A: mean:-0.0cm;RMS:5.4cm","l");

   ci = TColor::GetColor("#ff00ff");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
   entry->SetMarkerColor(1);
   entry->SetMarkerStyle(21);
   entry->SetMarkerSize(1);
   entry->SetTextFont(62);
   entry=leg->AddEntry("SegX_2","ME11B: mean:-0.1cm;RMS:8.4cm","l");

   ci = TColor::GetColor("#ff9999");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
   entry->SetMarkerColor(1);
   entry->SetMarkerStyle(21);
   entry->SetMarkerSize(1);
   entry->SetTextFont(62);
   entry=leg->AddEntry("SegX_3","ME12+13: mean:-0.1cm;RMS:17.6cm","l");
   entry->SetLineColor(1);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
   entry->SetMarkerColor(1);
   entry->SetMarkerStyle(21);
   entry->SetMarkerSize(1);
   entry->SetTextFont(62);
   entry=leg->AddEntry("SegX_4","ME2: mean:-0.1cm;RMS:23.2cm","l");

   ci = TColor::GetColor("#ff0000");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);
//.........这里部分代码省略.........

void RDK2AnalysisPlotter::makeEPPlot(CoDet detType)
{
    gROOT->cd();
    TCanvas* theCanvas;
    TPad* mainPad;
    TPad* titlePad;
    TPad* pTPad;
    TPad* eEPad;
    TPad* pEPad;
    TPad* pTPadSubs[3];
    TPad* eEPadSubs[3];
    TPad* pEPadSubs[3];


    int numTitleLines=getTitleBoxLines( detType);
    int canvasNumPixelsYPlot=600;
    int canvasNumPixelsYPlots=3*canvasNumPixelsYPlot;
    int canvasNumPixelsYTitle=40*numTitleLines;
    int canvasNumPixelsY=canvasNumPixelsYPlots+canvasNumPixelsYTitle;

    theCanvas = new TCanvas("EPExpMCAnalysisComparisonPlot","EPExpMCAnalysisComparisonPlot",10,10,1200,canvasNumPixelsY);
     mainPad=(TPad*) theCanvas->GetPad(0);


    double ylow,yhigh;

    yhigh=1;
    ylow=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY;
    titlePad=new TPad("titlePad", "titlePad", 0., ylow, 1., yhigh, -1, 1, 1);

    yhigh=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY;
    ylow=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY-canvasNumPixelsYPlot/(double)canvasNumPixelsY;
    pTPad=new TPad("pTPad", "pTPad", 0., ylow, 1., yhigh, -1, 1, 1);

    yhigh=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY-canvasNumPixelsYPlot/(double)canvasNumPixelsY;
    ylow=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY- 2*canvasNumPixelsYPlot/(double)canvasNumPixelsY;
    eEPad=new TPad("eEPad", "eEPad", 0., ylow, 1., yhigh, -1, 1, 1);

    yhigh=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY-2*canvasNumPixelsYPlot/(double)canvasNumPixelsY;
    ylow=0;
    pEPad=new TPad("pEPad", "pEPad", 0., ylow, 1., yhigh, -1, 1, 1);

    titlePad->Draw();
    pTPad->Draw();
    eEPad->Draw();
    pEPad->Draw();



//    titlePad=(TPad*) theCanvas->GetPad(1);
//    pTPad=(TPad*) theCanvas->GetPad(2);
//    eEPad=(TPad*) theCanvas->GetPad(3);
//    pEPad=(TPad*) theCanvas->GetPad(4);

    titlePad->SetMargin(0.05,0.0,0.0,0.0);
    pTPad->SetMargin(0.05,0.0,0.0,0.05);
    pEPad->SetMargin(0.05,0.0,0.0,0.05);
    eEPad->SetMargin(0.05,0.0,0.0,0.05);


    if(numExp+numMC >1)
    {
        pTPad->Divide(2);
        eEPad->Divide(2);
        pEPad->Divide(2);

        pTPadSubs[0]=(TPad*) pTPad->GetPad(1);
        eEPadSubs[0]=(TPad*) eEPad->GetPad(1);
        pEPadSubs[0]=(TPad*) pEPad->GetPad(1);

        TPad* tempPad;
        tempPad=(TPad*) pTPad->GetPad(2);
        tempPad->Divide(1,2);
        pTPadSubs[1]=(TPad*) tempPad->GetPad(1);
        pTPadSubs[2]=(TPad*) tempPad->GetPad(2);

        tempPad=(TPad*) pEPad->GetPad(2);
        tempPad->Divide(1,2);
        pEPadSubs[1]=(TPad*) tempPad->GetPad(1);
        pEPadSubs[2]=(TPad*) tempPad->GetPad(2);

        tempPad=(TPad*) eEPad->GetPad(2);
        tempPad->Divide(1,2);
        eEPadSubs[1]=(TPad*) tempPad->GetPad(1);
        eEPadSubs[2]=(TPad*) tempPad->GetPad(2);

        TPad* allBasePads[10]={titlePad,pTPadSubs[0],pTPadSubs[1],pTPadSubs[2],eEPadSubs[0],eEPadSubs[1],eEPadSubs[2],pEPadSubs[0],pEPadSubs[1],pEPadSubs[2]};
        for (int i = 1;i< 10;i++)
        {
            allBasePads[i]->SetGrid(1,1);
            allBasePads[i]->SetTickx(1);
            allBasePads[i]->SetTicky(1);
        }
        pTPadSubs[0]->SetMargin(0.13,0.1,.1,0.1);
        pEPadSubs[0]->SetMargin(0.13,0.1,.1,0.1);
        eEPadSubs[0]->SetMargin(0.13,0.1,.1,0.1);

        pTPadSubs[1]->SetMargin(0.1,0.1,.1,0.1);
        pEPadSubs[1]->SetMargin(0.1,0.1,.1,0.1);
        eEPadSubs[1]->SetMargin(0.1,0.1,.1,0.1);
//.........这里部分代码省略.........

void proc_pileup(const AliVdMMetaData& vdmMetaData,
                 AliVdMScanData& allData,
                 const char* classAC,
                 const char* classAnotC,
                 const char* classCnotA,
                 const std::vector<Double_t>& par0)
{
  typedef std::map<Short_t, TGraphErrors> map_t; // BCID -> TGraphErrors
  map_t gAnotC, gCnotA;  // one-arm/two-arm ratios

  // (1) fill one-arm/two-arm ratio graphs for all BCIDs
  for (Int_t iScan=0; iScan<4; ++iScan) {
    AliVdMTree& vtAND   = allData.GetMap(iScan)[classAC];
    AliVdMTree& vtAnotC = allData.GetMap(iScan)[classAnotC];
    AliVdMTree& vtCnotA = allData.GetMap(iScan)[classCnotA];
    vtAND.Zip3([&gAnotC,&gCnotA](const AliVdMTree::DefaultBranchData& d,
                                 AliVdMTree::branchMap_t& mapAC,
                                 AliVdMTree::branchMap_t& mapAnotC,
                                 AliVdMTree::branchMap_t& mapCnotA)
               {
                 AliVdMTree::ValErr v1 = mapAnotC["rate"];
                 v1 /= mapAC["rate"];
                 // v1 /= mapAC["relBkgd"];
                 v1 *= mapAnotC["relBkgd"];
                 if (!v1.isInf() && v1.val()) {
                   const Int_t m1 = gAnotC[d.BCID()].GetN();
                   gAnotC[d.BCID()].SetPoint     (m1, mapAC["mu"].val(), v1.val());
                   gAnotC[d.BCID()].SetPointError(m1, mapAC["mu"].err(), v1.err());
                 }
                 AliVdMTree::ValErr v2 = mapCnotA["rate"];
                 v2 /= mapAC["rate"];
                 // v2 /= mapAC["relBkgd"];
                 v2 *= mapCnotA["relBkgd"];
                 if (!v2.isInf() && v2.val()) {
                   const Int_t m2 = gCnotA[d.BCID()].GetN();
                   gCnotA[d.BCID()].SetPoint     (m2, mapAC["mu"].val(), v2.val());
                   gCnotA[d.BCID()].SetPointError(m2, mapAC["mu"].err(), v2.err());
                 }
               },
               vtAnotC, vtCnotA);
  }

  // (2) fit model
  AliVdMPileup pileupModel;

  TString pn = TString::Format("pileup_%s.pdf", classAC);
  TCanvas *c1 = new TCanvas;
  c1->SaveAs(pn+"[");

  const AliTriggerBCMask& bcMask = vdmMetaData.GetTriggerBCMask();
  const Int_t             nBCs   = bcMask.GetNUnmaskedBCs();

  TH1 *hPar[5] = {
    SetAttr(new TH1D("hrA",      ";BCID;r_{A}",           nBCs,0,nBCs), kRed),
    SetAttr(new TH1D("hrC",      ";BCID;r_{C}",           nBCs,0,nBCs), kBlue),
    SetAttr(new TH1D("hbkgdA",   ";BCID;bkgd_{A}",        nBCs,0,nBCs), kRed),
    SetAttr(new TH1D("hbkgdC",   ";BCID;bkgd_{C}",        nBCs,0,nBCs), kBlue),
    SetAttr(new TH1D("hChi2NDF", ";BCID;#chi^{2}/n.d.f.", nBCs,0,nBCs), kBlue)
  };

  for (Int_t bc=0, counter=0; bc<3564; ++bc) {
    if (bcMask.GetMask(bc))
      continue;

    const TString binLabel = TString::Format("%d", bc);
    for (Int_t i=0; i<5; ++i)
      hPar[i]->GetXaxis()->SetBinLabel(1+counter, binLabel);

    c1->Clear();
    c1->SetLogx();
    c1->SetLogy();
    TH1 *hf = c1->DrawFrame(1e-6, 0.01, 0.5, 20);
    hf->SetTitle(TString::Format("BCID=%d %s;two-arm #mu;one-arm/two-arm", bc, classAC));
    pileupModel.DoFit(&gAnotC[bc], &gCnotA[bc], &par0[0]);

    SetAttr(&gAnotC[bc], kBlue);
    SetAttr(&gCnotA[bc], kRed);
    gAnotC[bc].Draw("PE");
    gCnotA[bc].Draw("PE");

    TF1 *fAnotC = SetAttr(new TF1("fAnotC", &pileupModel, &AliVdMPileup::fcnAnotC, 1e-6, 0.5, 5), kBlue);
    fAnotC->SetParameters(pileupModel.GetPar());
    fAnotC->SetNpx(1000);
    fAnotC->Draw("same");

    TF1 *fCnotA = SetAttr(new TF1("fCnotA", &pileupModel, &AliVdMPileup::fcnCnotA, 1e-6, 0.5, 5), kRed);
    fCnotA->SetParameters(pileupModel.GetPar());
    fCnotA->SetNpx(1000);
    fCnotA->Draw("same");

    TLegend *leg = new TLegend(0.6, 0.75, 0.9, 0.9);
    leg->AddEntry(&gCnotA[bc], "AnotC/AandC", "PEL");
    leg->AddEntry(&gAnotC[bc], "CnotA/AandC", "PEL");
    leg->Draw();

    TPaveText *pt = new TPaveText(0.6, 0.4, 0.9, 0.7, "NDC NB");
    pt->SetFillStyle(0);
    pt->AddText(TString::Format("#chi^{2}/n.d.f = %.0f/%.0f = %.2f", pileupModel.GetChi2(), pileupModel.GetNDF(), pileupModel.GetChi2()/pileupModel.GetNDF()));
    {
      double curval,err, lowlim, uplim;
//.........这里部分代码省略.........

int main(int argc, char* argv[]) {
  AmFitter amfitter(argv[1]);
  const double* results;

  int count = 0;

  double x_wdeg[300];
  double y_fval[300];
  double y_X[300];

  int n_blacklist;

  for(double curr_wdeg = 2.7; curr_wdeg <= 2.7; curr_wdeg += 0.01) {
    amfitter.set_wdeg(curr_wdeg);

    for(int i=0;i<N_HSK;i++) {
      amfitter.set_file_id(i);
      results = amfitter.findmin();
      //amfitter._blacklist[i] = 1;
    }
    
    //amfitter._blacklist[19] = 0;

    n_blacklist = 0;
    for(int i=0;i<N_HSK;i++) {
      n_blacklist += amfitter._blacklist[i];
    }

    cout << n_blacklist << " files blacklisted" << endl;

    // do final fit accounting for blacklisted files
    amfitter.set_file_id(-1);
    results = amfitter.findmin();
  
    x_wdeg[count] = curr_wdeg;
    y_fval[count] = amfitter.chi_sq_sum(results[0],results[1],results[2],results[3]);
    y_X[count] = results[1];
    cout << ": wdeg = " << curr_wdeg << ": wfval = " << y_fval[count] << ": X = " << y_X[count] << endl;

    count++;

  }

  // make some plot to see if things are well behaved
  /*TCanvas cfval;
  TGraph grfval = TGraph(count,x_wdeg,y_fval);
  grfval.Draw("A*");
  cfval.SaveAs("fval.gif");

  TCanvas cX;
  TGraph grX = TGraph(count,x_wdeg,y_X);
  grX.Draw("A*");
  cX.SaveAs("X.gif");*/


  cout << endl << "AMBER GALACTIC CALIBRATION FIT RESULTS" << endl << endl;

  cout << "const B: " << results[0] << endl;
  cout << "const X: " << results[1] << endl;
  cout << "noise: " << results[2] << endl;
  cout << "time shift: " << results[3] << endl;
  cout << endl;

  // stuff to write the resulting fit to a .gif
  //double wdeg = results[0];
  //double A = results[1];
  double B = results[0];
  double X = results[1];
  //double C = results[3];
  double n = results[2];
  double t_shift = results[3];
  
  double val = amfitter.chi_sq_sum(B, X, n, t_shift);

  cout << "VALUE: " << val << endl;
  cout << endl;

  // if we're fitting by hand, this is in effect
  //wdeg = 2.2;
  //A = 66.2721;
  //B = 0.02714;
  //X = -133.019;
  //n = 6.89e-10;
  //t_shift = 1.846;

  double y1[215];
  double y2[215];
  double x[215];

  for(int i=0;i<215;i++) {
    y1[i-0] = 10*log10((5.0E8)*(1.38065E-23)*(amfitter.temps_func(amfitter.time_arr[0][i]))/.001);
    //y2[i] = 10*log((A*exp(amfitter.hsk_func(amfitter.time_arr[0][i])*B+C) - n)/1);
    y2[i-0] = 10*log10( (pow(10,0.1*(amfitter.hsk_func(amfitter.time_arr[0][i]+t_shift)*B+X)) - n)/1);
    x[i-0] = amfitter.time_arr[0][i] - amfitter.time_arr[0][0];
    if( (i>100) && (i<120) ) {
      cout << x[i-0] << " " << y1[i-0] << " " << y2[i-0] << " " << amfitter.hsk_func(amfitter.time_arr[0][i]) << " " << amfitter.temps_func(amfitter.time_arr[0][i]) << " " << amfitter.time_arr[0][i] << endl;
    }
  }

  TGraph* gr1 = new TGraph(215,x,y1);
//.........这里部分代码省略.........

void RooToyMCFit1Bin(const int _bin){
    if(_bin<1 || _bin>8){
        cout << "Bin number should be between 1 and 8" << endl;
        return;
    }
    stringstream out;
    out.str("");
    out << "toyMC_" << _sigma_over_tau << "_" << _purity << "_" << mistag_rate << ".root";
    TFile* file = TFile::Open(out.str().c_str());
    TTree* tree = (TTree*)file->Get("ToyTree");
    cout << "The tree has been readed from the file " << out.str().c_str() << endl;

    RooRealVar tau("tau","tau",_tau,"ps"); tau.setConstant(kTRUE);
    RooRealVar dm("dm","dm",_dm,"ps^{-1}"); dm.setConstant(kTRUE);
    RooRealVar sin2beta("sin2beta","sin2beta",_sin2beta,-5.,5.); if(constBeta) sin2beta.setConstant(kTRUE);
    RooRealVar cos2beta("cos2beta","cos2beta",_cos2beta,-5.,5.); if(constBeta) cos2beta.setConstant(kTRUE);
    RooRealVar dt("dt","#Deltat",-5.,5.,"ps");
    RooRealVar avgMisgat("avgMisgat","avgMisgat",mistag_rate,0.0,0.5); if(constMistag) avgMisgat.setConstant(kTRUE);
    RooRealVar delMisgat("delMisgat","delMisgat",0); delMisgat.setConstant(kTRUE);
    RooRealVar mu("mu","mu",0); mu.setConstant(kTRUE);
    RooRealVar moment("moment","moment",0.);  moment.setConstant(kTRUE);
    RooRealVar parity("parity","parity",-1.); parity.setConstant(kTRUE);

    cout << "Preparing coefficients..." << endl;
    RooRealVar*    K  = new RooRealVar("K","K",K8[_bin-1],0.,1.); if(constK) K->setConstant(kTRUE);
    RooFormulaVar* Kb = new RooFormulaVar("Kb","Kb","[email protected]",RooArgList(*K));

    RooRealVar* C = new RooRealVar("C","C",_C[_bin-1]); C->setConstant(kTRUE);
    RooRealVar* S = new RooRealVar("S","S",_S[_bin-1]); S->setConstant(kTRUE);

    RooFormulaVar* a1  = new RooFormulaVar("a1","a1","-(@[email protected])/(@[email protected])",RooArgList(*K,*Kb));
    RooFormulaVar* a1b = new RooFormulaVar("a1b","a1b","(@[email protected])/(@[email protected])",RooArgList(*K,*Kb));
    RooFormulaVar* a2  = new RooFormulaVar("a2","a2","(@[email protected])/(@[email protected])",RooArgList(*K,*Kb));
    RooFormulaVar* a2b = new RooFormulaVar("a2b","a2b","-(@[email protected])/(@[email protected])",RooArgList(*K,*Kb));

    RooFormulaVar* b1  = new RooFormulaVar("b1","b1","2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected]);",RooArgList(*K,*Kb,*C,*S,sin2beta,cos2beta));
    RooFormulaVar* b1b = new RooFormulaVar("b1b","b1b","2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected]);",RooArgList(*K,*Kb,*C,*S,sin2beta,cos2beta));
    RooFormulaVar* b2  = new RooFormulaVar("b2","b2","-2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected]);",RooArgList(*K,*Kb,*C,*S,sin2beta,cos2beta));
    RooFormulaVar* b2b = new RooFormulaVar("b2b","b2b","-2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected]);",RooArgList(*K,*Kb,*C,*S,sin2beta,cos2beta));

    RooRealVar* dgamma = new RooRealVar("dgamma","dgamma",0.); dgamma->setConstant(kTRUE);
    RooRealVar* f0 = new RooRealVar("f0","f0",1.); f0->setConstant(kTRUE);
    RooRealVar* f1 = new RooRealVar("f1","f1",0.); f1->setConstant(kTRUE);

    RooCategory tag("tag","tag");
    tag.defineType("B0",1);
    tag.defineType("anti-B0",-1);

    RooCategory bin("bin","bin");
    bin.defineType("bin",_bin);
    bin.defineType("binb",-_bin);

    RooSuperCategory bintag("bintag","bintag",RooArgSet(bin,tag));

    RooDataSet d("data","data",tree,RooArgSet(dt,bin,tag));
    cout << "DataSet is ready." << endl;
    d.Print();

    RooRealVar mean("mean","mean",0.,"ps"); mean.setConstant(kTRUE);
    RooRealVar sigma("sigma","sigma",_sigma_over_tau*_tau,0.,_tau,"ps"); if(constSigma) sigma.setConstant(kTRUE);
    RooGaussModel rf("rf","rf",dt,mean,sigma);
//    RooTruthModel rf("rf","rf",dt);
    RooGaussian rfpdf("rfpdf","rfpdf",dt,mean,sigma);

    cout << "Preparing PDFs..." << endl;
    RooBDecay* sigpdf1  = new RooBDecay("sigpdf1","sigpdf1",dt,tau,*dgamma,*f0,*f1,*a1,*b1,dm,rf,RooBDecay::DoubleSided);
    RooBDecay* sigpdf2  = new RooBDecay("sigpdf2","sigpdf2",dt,tau,*dgamma,*f0,*f1,*a2,*b2,dm,rf,RooBDecay::DoubleSided);
    RooBDecay* sigpdf1b  = new RooBDecay("sigpdf1b","sigpdf1b",dt,tau,*dgamma,*f0,*f1,*a1b,*b1b,dm,rf,RooBDecay::DoubleSided);
    RooBDecay* sigpdf2b  = new RooBDecay("sigpdf2b","sigpdf2b",dt,tau,*dgamma,*f0,*f1,*a2b,*b2b,dm,rf,RooBDecay::DoubleSided);

    RooRealVar fsig("fsig","fsigs",_purity,0.,1.);  if(constFSig) fsig.setConstant(kTRUE);
    RooAddPdf* PDF1 = new RooAddPdf("PDF1","PDF1",RooArgList(*sigpdf1,rfpdf),RooArgList(fsig));
    RooAddPdf* PDF2 = new RooAddPdf("PDF2","PDF2",RooArgList(*sigpdf2,rfpdf),RooArgList(fsig));
    RooAddPdf* PDF1b= new RooAddPdf("PDF1b","PDF1b",RooArgList(*sigpdf1b,rfpdf),RooArgList(fsig));
    RooAddPdf* PDF2b= new RooAddPdf("PDF2b","PDF2b",RooArgList(*sigpdf2b,rfpdf),RooArgList(fsig));

    //Adding mistaging
    RooAddPdf* pdf1 = new RooAddPdf("pdf1","pdf1",RooArgList(*PDF2,*PDF1),RooArgList(avgMisgat));
    RooAddPdf* pdf2 = new RooAddPdf("pdf2","pdf2",RooArgList(*PDF1,*PDF2),RooArgList(avgMisgat));
    RooAddPdf* pdf1b= new RooAddPdf("pdf1b","pdf1b",RooArgList(*PDF2b,*PDF1b),RooArgList(avgMisgat));
    RooAddPdf* pdf2b= new RooAddPdf("pdf2b","pdf2b",RooArgList(*PDF1b,*PDF2b),RooArgList(avgMisgat));

    RooSimultaneous pdf("pdf","pdf",bintag);
    pdf.addPdf(*pdf1,"{bin;B0}");
    pdf.addPdf(*pdf2,"{bin;anti-B0}");
    pdf.addPdf(*pdf1b,"{binb;B0}");
    pdf.addPdf(*pdf2b,"{binb;anti-B0}");

    cout << "Fitting..." << endl;
    pdf.fitTo(d,Verbose(),Timer());

    cout << "Drawing plots." << endl;
    // Plus bin
    RooPlot* dtFrame = dt.frame();

    // B0
    out.str("");
    out << "tag == 1 && bin == " << _bin;
    cout << out.str() << endl;
    RooDataSet* ds = d.reduce(out.str().c_str());
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   phimuphimet90->SetBinError(4,0.000274051);
   phimuphimet90->SetBinError(5,0.0002641202);
   phimuphimet90->SetBinError(6,0.0002554581);
   phimuphimet90->SetBinError(7,0.0002714728);
   phimuphimet90->SetBinError(8,0.0002560081);
   phimuphimet90->SetBinError(9,0.0002614441);
   phimuphimet90->SetBinError(10,0.0002430443);
   phimuphimet90->SetEntries(2551);

   ci = TColor::GetColor("#0000ff");
   phimuphimet90->SetFillColor(ci);

   ci = TColor::GetColor("#0000ff");
   phimuphimet90->SetLineColor(ci);

   ci = TColor::GetColor("#0000ff");
   phimuphimet90->SetMarkerColor(ci);
   phimuphimet90->SetMarkerStyle(21);
   phimuphimet90->GetXaxis()->SetTitle("phimuphimet");
   phimuphimet90->GetXaxis()->SetLabelFont(42);
   phimuphimet90->GetXaxis()->SetLabelSize(0.035);
   phimuphimet90->GetXaxis()->SetTitleSize(0.035);
   phimuphimet90->GetXaxis()->SetTitleFont(42);
   phimuphimet90->GetYaxis()->SetTitle("Events/pb");
   phimuphimet90->GetYaxis()->SetLabelFont(42);
   phimuphimet90->GetYaxis()->SetLabelSize(0.035);
   phimuphimet90->GetYaxis()->SetTitleSize(0.035);
   phimuphimet90->GetYaxis()->SetTitleFont(42);
   phimuphimet90->GetZaxis()->SetLabelFont(42);
   phimuphimet90->GetZaxis()->SetLabelSize(0.035);
   phimuphimet90->GetZaxis()->SetTitleSize(0.035);
   phimuphimet90->GetZaxis()->SetTitleFont(42);
   phimuphimet->Add(phimuphimet,"");
   phimuphimet->Draw("nostack");
   
   TPaveText *pt = new TPaveText(0.3326724,0.9342857,0.6673276,0.995,"blNDC");
   pt->SetName("title");
   pt->SetBorderSize(0);
   pt->SetFillColor(0);
   pt->SetFillStyle(0);
   pt->SetTextFont(42);
   TText *AText = pt->AddText("phimuphimet");
   pt->Draw();
   
   TLegend *leg = new TLegend(0.54023,0.639881,0.938218,0.924107,NULL,"brNDC");
   leg->SetBorderSize(1);
   leg->SetTextSize(0.034965);
   leg->SetLineColor(1);
   leg->SetLineStyle(1);
   leg->SetLineWidth(1);
   leg->SetFillColor(0);
   leg->SetFillStyle(1001);
   TLegendEntry *entry=leg->AddEntry("phimuphimet","phimuphimet_QCD","lp");

   ci = TColor::GetColor("#00cc00");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);

   ci = TColor::GetColor("#00cc00");
   entry->SetMarkerColor(ci);
   entry->SetMarkerStyle(22);
   entry->SetMarkerSize(1);
   entry->SetTextFont(42);
   entry=leg->AddEntry("phimuphimet","phimuphimet_WJetsToLNu","lp");

//.........这里部分代码省略.........
   hmt60->SetBinContent(2,0.0002683464);
   hmt60->SetBinContent(3,0.0001341732);
   hmt60->SetBinContent(4,1.677165e-05);
   hmt60->SetBinError(1,0.0001006299);
   hmt60->SetBinError(2,6.708659e-05);
   hmt60->SetBinError(3,4.743738e-05);
   hmt60->SetBinError(4,1.677165e-05);
   hmt60->SetEntries(61);

   ci = TColor::GetColor("#0000ff");
   hmt60->SetFillColor(ci);

   ci = TColor::GetColor("#0000ff");
   hmt60->SetLineColor(ci);

   ci = TColor::GetColor("#0000ff");
   hmt60->SetMarkerColor(ci);
   hmt60->SetMarkerStyle(21);
   hmt60->GetXaxis()->SetTitle("hmt");
   hmt60->GetXaxis()->SetLabelFont(42);
   hmt60->GetXaxis()->SetLabelSize(0.035);
   hmt60->GetXaxis()->SetTitleSize(0.035);
   hmt60->GetXaxis()->SetTitleFont(42);
   hmt60->GetYaxis()->SetTitle("Events/pb");
   hmt60->GetYaxis()->SetLabelFont(42);
   hmt60->GetYaxis()->SetLabelSize(0.035);
   hmt60->GetYaxis()->SetTitleSize(0.035);
   hmt60->GetYaxis()->SetTitleFont(42);
   hmt60->GetZaxis()->SetLabelFont(42);
   hmt60->GetZaxis()->SetLabelSize(0.035);
   hmt60->GetZaxis()->SetTitleSize(0.035);
   hmt60->GetZaxis()->SetTitleFont(42);
   hmt->Add(hmt,"");
   hmt->Draw("nostack");
   
   TPaveText *pt = new TPaveText(0.444023,0.94,0.555977,0.995,"blNDC");
   pt->SetName("title");
   pt->SetBorderSize(0);
   pt->SetFillColor(0);
   pt->SetFillStyle(0);
   pt->SetTextFont(42);
   TText *AText = pt->AddText("hmt");
   pt->Draw();
   
   TLegend *leg = new TLegend(0.54023,0.639881,0.938218,0.924107,NULL,"brNDC");
   leg->SetBorderSize(1);
   leg->SetTextSize(0.034965);
   leg->SetLineColor(1);
   leg->SetLineStyle(1);
   leg->SetLineWidth(1);
   leg->SetFillColor(0);
   leg->SetFillStyle(1001);
   TLegendEntry *entry=leg->AddEntry("hmt","hmt_QCD","lp");

   ci = TColor::GetColor("#00cc00");
   entry->SetLineColor(ci);
   entry->SetLineStyle(1);
   entry->SetLineWidth(1);

   ci = TColor::GetColor("#00cc00");
   entry->SetMarkerColor(ci);
   entry->SetMarkerStyle(22);
   entry->SetMarkerSize(1);
   entry->SetTextFont(42);
   entry=leg->AddEntry("hmt","hmt_WJetsToLNu","lp");

//.........这里部分代码省略.........
  }
  else{
#ifdef MSSM
    bbH  ->Add(WJets);
    ggH  ->Add(bbH);
#else
#ifndef DROP_SIGNAL
    VH   ->Add(WJets);
    qqH  ->Add(VH );
    ggH  ->Add(qqH);
#endif
#endif
  }


  /*
    mass plot before and after fit
  */
  TCanvas* canv = MakeCanvas("canv", "histograms", 600, 600);
  canv->cd();
  if(log){ canv->SetLogy(1); }
#if defined MSSM
  if(!log){ data->GetXaxis()->SetRange(0, data->FindBin(345)); } else{ data->GetXaxis()->SetRange(0, data->FindBin(UPPER_EDGE)); };
#else
  data->GetXaxis()->SetRange(0, data->FindBin(345));
#endif
  data->SetNdivisions(505);
  data->SetMinimum(min);
#ifndef DROP_SIGNAL
  data->SetMaximum(max>0 ? max : std::max(std::max(maximum(data, log), maximum(ZTT, log)), maximum(ggH, log)));
#else
  data->SetMaximum(max>0 ? max : std::max(maximum(data, log), maximum(ZTT, log)));
#endif
  data->Draw("e");

  TH1F* errorBand = (TH1F*)ZEE ->Clone("errorBand");
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(1);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
    if(errorBand->GetBinContent(idx)>0){
      std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
      break;
    }
  }
  if(log){
    ZEE->Draw("histsame");
    ZTT->Draw("histsame");
    TTJ->Draw("histsame");
    QCD->Draw("histsame");
    WJets->Draw("histsame");
    //Dibosons->Draw("histsame");
    $DRAW_ERROR
#ifndef DROP_SIGNAL
    ggH->Draw("histsame");
#endif

  }
  else{
#ifndef DROP_SIGNAL
    ggH  ->Draw("histsame");
#endif
    ZEE->Draw("histsame");
    ZTT->Draw("histsame");
    TTJ->Draw("histsame");

void drawComparisonResponse(const string& canvasName, TMultiGraph *mgResponse, TGraphErrors *gResponseMC, TGraphErrors *gResponseData, TGraph *gratio,const string& mcSample, const string& path, bool doFit = false, bool setRangeUser = true) { // Mikko's tdrstyle_mod14
  TCanvas *c2 = tdrDiCanvas(canvasName.c_str(),mgResponse,gratio,2,0);
  c2->cd(1);

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  if (setRangeUser) {
    mgResponse->SetMaximum(1.05);
    mgResponse->SetMinimum(0.9);
  } else {
    mgResponse->SetMaximum(0.5);
    mgResponse->SetMinimum(0.3);    
  }
  mgResponse->Draw("P");
  gPad->RedrawAxis();
  c2->SetLogx(1);
  if(doFit) {
    TF1* myLinFit = new TF1("myLinFit", "[0]*x + [1]", mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());
    myLinFit->SetParameter(0.5, 1.);
    myLinFit->SetLineColor(1);
    myLinFit->SetLineWidth(2);


    c2->Update();
    gResponseData->Fit(myLinFit, "RQ");
    double linfitValue = myLinFit->GetParameter(0);
    double linfitError = myLinFit->GetParError(0);

    TPaveText* linfitlabel = new TPaveText(0.55, 0.77, 0.88, 0.83, "brNDC");
    linfitlabel->SetTextSize(0.035);
    linfitlabel->SetFillColor(0);
    linfitlabel->SetTextFont(43);
    TString linfitLabelText = TString::Format("Data slope: %.5f #pm %.5f", linfitValue, linfitError);
    linfitlabel->AddText(linfitLabelText);
    linfitlabel->Draw("same");

    gPad->RedrawAxis();	
  }

  //TLegend* legend = new TLegend(0.55, 0.15, 0.92, 0.38);
  TLegend* legend = new TLegend(0.65, 0.15, 0.95, 0.30);
  legend->SetFillColor(kWhite);
  legend->SetFillStyle(kNone);
  legend->SetTextFont(43);
  legend->SetBorderSize(0);
  //legend->SetTextSize(0.045);
  legend->AddEntry(gResponseMC,"MC","p");
  legend->AddEntry(gResponseData,"Data 2012","p");
  legend->Draw("same");

  legend->Draw("same");

  c2->cd(2);

  TF1* ratioFit = new TF1("ratioFit", "[0]", mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());
  ratioFit->SetParameter(0, 0.);
  ratioFit->SetLineColor(46);
  ratioFit->SetLineWidth(2);
  gratio->Draw("APE1");
  gratio->GetYaxis()->SetTitle("Data / MC");
  // 	//gratio->GetXaxis()->SetLabelOffset(0.1);
  // 	gratio->GetXaxis()->SetLabelFont(42);
  // 	gratio->GetXaxis()->SetLabelSize(0.06);
  // 	//gratio->GetYaxis()->SetLabelOffset(0.1);
  // 	gratio->GetYaxis()->SetLabelFont(42);
  // 	gratio->GetYaxis()->SetLabelSize(0.06);
  // 	gratio->GetXaxis()->SetTitleOffset(-0.6);
  // 	gratio->GetXaxis()->SetTitleFont(42);
  // 	gratio->GetXaxis()->SetTitleSize(0.06);
  //gratio->GetYaxis()->SetTitleOffset(1.3);
  // 	gratio->GetYaxis()->SetTitleFont(42);
  // 	gratio->GetYaxis()->SetTitleSize(0.06);

  gratio->SetMarkerSize(1.0);
  gratio->SetMaximum(1.05);
  gratio->SetMinimum(0.95);
  gratio->GetXaxis()->SetLimits(mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());

  c2->Update();
  gratio->Fit(ratioFit, "RQ");
  //gratio->GetYaxis()->SetRangeUser(-1,3);
  double fitValue = ratioFit->GetParameter(0);
  double fitError = ratioFit->GetParError(0);

  TPaveText* fitlabel = new TPaveText(0.55, 0.77, 0.88, 0.83, "brNDC");
  fitlabel->SetTextSize(0.08);
  fitlabel->SetFillColor(0);
  fitlabel->SetTextFont(42);
  TString fitLabelText = TString::Format("Fit: %.4f #pm %.4f", fitValue, fitError);
  fitlabel->AddText(fitLabelText);
  fitlabel->Draw("same");

  gPad->RedrawAxis();
  c2->SaveAs(path.c_str());
}

void mergePlotResponse(const string& canvasName, TMultiGraph *mgResponse, TGraph *gResponseMC_MJB, TGraph *gResponseData_MJB, TGraph *gratio_MJB, TGraph *gResponseMC_MPF, TGraph *gResponseData_MPF, TGraph *gratio_MPF, double y_min, double y_max, const string& label_MJB, const string& label_MPF, const string& path) { // Mikko's tdrstyle_mod14
  TCanvas *c3 = tdrDiCanvas(canvasName.c_str(),mgResponse,gratio_MJB,2,0);
  c3->cd(1);

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  mgResponse->SetMaximum(y_max);
  mgResponse->SetMinimum(y_min);
  mgResponse->Draw("P");
  gPad->RedrawAxis();
  //c3->SetLogx(1);

  //TLegend* legend = new TLegend(0.55, 0.15, 0.92, 0.38);
  TLegend* legend = new TLegend(0.63, 0.15, 0.93, 0.30);
  legend->SetFillColor(kWhite);
  legend->SetFillStyle(kNone);
  legend->SetTextFont(43);
  legend->SetBorderSize(0);
  //legend->SetTextSize(0.045);
  string toLeg = label_MJB + ", MC";
  legend->AddEntry(gResponseMC_MJB, toLeg.c_str(),"p");
  toLeg = label_MJB + ", data";
  legend->AddEntry(gResponseData_MJB,toLeg.c_str(),"p");
  toLeg = label_MPF + ", MC";
  legend->AddEntry(gResponseMC_MPF,toLeg.c_str(),"p");
  toLeg = label_MPF + ", data";
  legend->AddEntry(gResponseData_MPF,toLeg.c_str(),"p");
  legend->Draw("same");


  c3->cd(2);
  TF1* ratioFit_MJB = new TF1("ratioFit_MJB", "[0]", mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());
  ratioFit_MJB->SetParameter(0, 0.);
  ratioFit_MJB->SetLineColor(gratio_MJB->GetLineColor());
  ratioFit_MJB->SetLineWidth(0);

  gratio_MJB->Draw("APE1");
  gratio_MJB->GetYaxis()->SetTitle("Data / MC");
  gratio_MJB->SetMarkerSize(1.0);
  gratio_MJB->SetMaximum(1.05);
  gratio_MJB->SetMinimum(0.95);
  gratio_MJB->GetXaxis()->SetLimits(mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());

  gratio_MPF->Draw("PE1same");

  c3->Update();
  gratio_MJB->Fit(ratioFit_MJB, "RQ");
  //gratio_MJB->GetYaxis()->SetRangeUser(-1,3);
  double fitValue_MJB = ratioFit_MJB->GetParameter(0);
  double fitError_MJB = ratioFit_MJB->GetParError(0);

  TF1* ratioFit_MPF = new TF1("ratioFit_MPF", "[0]", mgResponse->GetXaxis()->GetXmin(),mgResponse->GetXaxis()->GetXmax());
  ratioFit_MPF->SetParameter(0, 0.);
  ratioFit_MPF->SetLineColor(gratio_MPF->GetLineColor());
  ratioFit_MPF->SetLineWidth(0);
  gratio_MPF->Fit(ratioFit_MPF, "RQ");
  //gratio_MPF->GetYaxis()->SetRangeUser(-1,3);
  double fitValue_MPF = ratioFit_MPF->GetParameter(0);
  double fitError_MPF = ratioFit_MPF->GetParError(0);

  TPaveText* fitlabel = new TPaveText(0.57, 0.82, 0.85, 0.84, "brNDC");
  fitlabel->SetTextSize(0.08);
  fitlabel->SetFillColor(0);
  fitlabel->SetTextFont(42);
  fitlabel->SetTextColor(gratio_MJB->GetLineColor());
  TString fitLabelTextMJB = TString::Format("Fit %s: %.4f #pm %.4f", label_MJB.c_str(), fitValue_MJB, fitError_MJB);
  fitlabel->AddText(fitLabelTextMJB);
  fitlabel->Draw("same");

  TPaveText* fitlabel_MPF = new TPaveText(0.57, 0.72, 0.85, 0.74, "brNDC");
  fitlabel_MPF->SetTextSize(0.08);
  fitlabel_MPF->SetFillColor(0);
  fitlabel_MPF->SetTextFont(42);
  fitlabel_MPF->SetTextColor(gratio_MPF->GetLineColor());
  TString fitLabelTextMPF = TString::Format("Fit %s: %.4f #pm %.4f", label_MPF.c_str(), fitValue_MPF, fitError_MPF);
  fitlabel_MPF->AddText(fitLabelTextMPF);
  fitlabel_MPF->Draw("same");

  gPad->RedrawAxis();
  c3->SaveAs(path.c_str());
}

void gluinoMass(double lumi=-1., double maxInstLumi=-1.) {
	
  if (lumi<0)
    lumi=877.;
  if (maxInstLumi<0)
    maxInstLumi=1300.;
  LimitPlots plots(lumi);
  
  plots.calculateCrossSections(7,4,39,9);
	
  // expected limit (1 and 2 sigma bands)
  TGraph* g_exp = plots.getExpMassLimitGluino();
  TGraphAsymmErrors* g_exp1 = plots.getExpMassLimitGluino1Sig();
  TGraphAsymmErrors* g_exp2 = plots.getExpMassLimitGluino2Sig();
  
  // three points on counting expt curve
  TGraph* g_gluino = plots.getMassLimitGluino();
  TGraph* g_stop = plots.getMassLimitStop();
  
  // one point from lifetime fit
  TGraph* g_tp = plots.getMassLimitGluinoTP();
  
  // theory prediction
  TGraph* g_thGluino = plots.getGluinoTheory();
  TGraph* g_thStop = plots.getStopTheory();
  
  TCanvas* canvas = new TCanvas("canvas");
  
  //canvas->SetGrid();
  canvas->SetLogy();
  
  TH1 * h;
  h = canvas->DrawFrame(300., .02, 1000., 1e2);
  //h->SetTitle("Beamgap Expt;m_{#tilde{g}} [GeV/c^{2}]; Stopped HSCP Cross Section #times BR [pb]");
  h->SetTitle("Beamgap Expt;m_{#tilde{g}} [GeV/c^{2}]; #sigma(pp #rightarrow #tilde{g}#tilde{g}) #times BR(#tilde{g} #rightarrow g#tilde{#chi}^{0}) [pb]");
  
  // not covered region
  TBox* nc = new TBox(100., .1, 150., 5e2);
  nc->SetFillStyle(3354);
  nc->SetFillColor(kRed-4);
  //nc->Draw();
  
  // details
  TPaveText* blurb = new TPaveText(300., 2, 550., 1e2);
  blurb->AddText("CMS Preliminary 2012");

  std::stringstream label;
  label<<"#int L dt = "<<lumi<<" fb^{-1}";
  blurb->AddText(label.str().c_str());
  label.str("");
  double peakInstLumi=maxInstLumi;
  int exponent=30;
  while (peakInstLumi>10) {
    peakInstLumi/=10.;
    ++exponent;
  }
  label<<"L^{max}_{inst} = "<<peakInstLumi<<" x 10^{"<<exponent<<"} cm^{-2}s^{-1}";
  blurb->AddText(label.str().c_str());
  label.str("");
  label << "#sqrt{s} = " << ENERGY << " TeV";
  blurb->AddText(label.str().c_str());
  blurb->AddText("m_{#tilde{g}} - m_{#tilde{#chi}^{0}} = 100 GeV/c^{2}");
  //blurb->AddText("m_{#tilde{t}} - m_{#tilde{#chi}^{0}} = 200 GeV/c^{2}");	
  blurb->SetTextFont(42);
  blurb->SetBorderSize(0);
  blurb->SetFillColor(0);
  blurb->SetShadowColor(0);
  blurb->SetTextAlign(12);
  blurb->SetTextSize(0.032);

  
  // legend
  TBox *legbg = new TBox(600., 2., 900., 1e2);
  legbg->Draw();
  TLegend *leg = new TLegend(600., 2., 900., 1e2,"95% C.L. Limits","");
  leg->SetTextSize(0.028);
  leg->SetBorderSize(0);
  leg->SetTextFont(42);
  leg->SetFillColor(0);
  leg->AddEntry(g_exp, "Expected: 10 #mus - 1000 s Counting Exp. ", "l");
  leg->AddEntry(g_exp1, "Expected #pm1#sigma: 10 #mus - 1000 s Counting Exp. ", "f");
  leg->AddEntry(g_exp2, "Expected #pm2#sigma: 10 #mus - 1000 s Counting Exp. ", "f");
  //	  leg->AddEntry(graph3, "Obs.: 10^{6} s Counting Exp.", "l");
  leg->AddEntry(g_gluino, "Obs.: 10 #mus - 1000 s Counting Exp. ", "l");
  leg->AddEntry(g_tp, "Obs.: 10 #mus Timing Profile ", "l");
  //leg->AddEntry(g_stop, "Obs.: 10 #mus - 1000 s Counting Exp. (#tilde{t})", "l");
  //leg->AddEntry(graph_em, "Obs.: 10 #mus - 1000 s Counting Exp. (EM only)", "l");
  //  leg->AddEntry(graph1, "Obs.: 570 ns Counting Exp.", "l");
  leg->Draw();
  
  
  
  // 2 sigma expected band
  g_exp2->SetLineColor(0);
  g_exp2->SetLineStyle(0);
  g_exp2->SetLineWidth(0);
  g_exp2->SetFillColor(5);
  g_exp2->SetFillStyle(1001);
  g_exp2->Draw("3");
  
//.........这里部分代码省略.........

void getCalibrationResults(TString url="TopMassCalibration.root",int nCategs=4)
{
  TString mpts[]={"161.5","163.5","166.5","169.5","172.5","175.5","178.5","181.5","184.5"};
  const size_t nmpts=sizeof(mpts)/sizeof(TString);

  std::vector<TString> labels;
  if(nCategs==2)
    {
      labels.push_back("1 b-tags");
      labels.push_back("#geq 2 b-tags");
    }
  else if(nCategs==4)
    {
      labels.push_back("1 b-tags (ee+#mu#mu)");
      labels.push_back("#geq 2 b-tags (ee+#mu#mu)");
      labels.push_back("1 b-tags (e#mu)");
      labels.push_back("#geq 2 b-tags (e#mu)");
    }
  else 
    {
      labels.push_back("1 b-tags (ee)");
      labels.push_back("#geq 2 b-tags (ee)");
      labels.push_back("1 b-tags (#mu#mu)");
      labels.push_back("#geq 2 b-tags (#mu#mu)");
      labels.push_back("1 b-tags (e#mu)");
      labels.push_back("#geq 2 b-tags (e#mu)");
    }
  
  TString outDir=gSystem->DirName(url);
  report << "Plots will be stored at @ " << outDir << endl;

  TGraphErrors *lingr = new TGraphErrors;  
  lingr->SetMarkerStyle(20);
  TGraphErrors *biasgr = new TGraphErrors; 
  biasgr->SetMarkerStyle(20);
  std::vector<TGraphErrors *> categbiasgr;
  for(int icat=0; icat<nCategs; icat++) categbiasgr.push_back( new TGraphErrors );
  TObjArray statUncResults;
  TObjArray pullResults;

  //get results from file
  TFile *fin=TFile::Open(url);
  for(size_t i=0; i<nmpts; i++)
    {
      Float_t trueMass=mpts[i].Atof();
      
      TH1D *massFitH=(TH1D *)fin->Get(mpts[i]+"/massfit");
      if(massFitH==0) continue;

      TF1 *ffunc=massFitH->GetFunction("gaus");
      float avgMass=ffunc->GetParameter(1);
      float avgMassErr=ffunc->GetParError(1);
      int ipt=lingr->GetN();
      lingr->SetPoint(ipt,trueMass,avgMass);
      lingr->SetPointError(ipt,0,avgMassErr);

      for(int icat=0; icat<=nCategs;icat++)
	{
	  TString postfix("");
	  TGraphErrors *gr=biasgr;
	  if(icat<nCategs){ postfix ="_"; postfix +=( icat+1); gr=categbiasgr[icat]; }
	  TH1D *biasH=(TH1D *)fin->Get(mpts[i]+"/bias"+postfix);
	  ffunc=biasH->GetFunction("gaus");
	  float avgBias=ffunc->GetParameter(1);
	  float avgBiasErr=ffunc->GetParError(1);   
	  ipt=gr->GetN();
	  gr->SetPoint(ipt,trueMass,avgBias);
	  gr->SetPointError(ipt,0,avgBiasErr);
	
	  if(trueMass==172.5) 
	    { 
	      TH1D *h=(TH1D *)fin->Get(mpts[i]+"/statunc"+postfix); 
	      h->SetDirectory(0);
	      h->SetLineColor(icat==nCategs? 1 : icat+2); 
	      h->SetMarkerColor(icat==nCategs? 1 : icat+2); 
	      h->SetMarkerStyle(1);
	      h->SetLineWidth(2);
	      statUncResults.Add(h);
	    }  
	}



      TH1D *pullH=(TH1D *)fin->Get(mpts[i]+"/pull");
      pullH->SetDirectory(0);
      pullResults.Add(pullH);
    }
  fin->Close();


  //plot results
  setStyle();
  gStyle->SetOptFit(0);

  TCanvas *c=new TCanvas("linc","linc",600,600);
  c->SetWindowSize(600,600);
  c->Divide(1,2);
  TPad *p=(TPad *)c->cd(1);
  p->SetPad(0,0.3,1.0,1.0);
  lingr->Draw("ap");
//.........这里部分代码省略.........

void sphere()
{
   /********************/
   /*Mpi initialization*/
   /********************/
   ROOT::Mpi::TEnvironment env(gApplication->Argc(), gApplication->Argv());
   ROOT::Mpi::TIntraCommunicator world;

   /*********************************/
   /*Error control, need 4 processes*/
   /*********************************/   
   if(world.Size()!=4)
   {
      std::cout<<"Error: You need run with 4 processes exactly, the sphere was spplitted in four regions. \nExample: mpirun -np 4 root -l  sphere.C\n";
      std::cout.flush();
      world.Abort(ROOT::Mpi::ERR_SIZE);  
   }
   
   /*********************************************************************/
   /*Section to calculate volume with parallel multimensional integrator*/
   /*********************************************************************/   
   //integrator creation root's style
   ROOT::Math::Functor1D function(&sphere_function);

   ROOT::Math::Integrator integrator;
   integrator.SetFunction(function);

   //Mpi's integrator wrapper, integrator from ROOT in one dimension
   ROOT::Mpi::Math::TIntegratorWrapper MpiIntegrator(&integrator);

   //integrate in parallel and send the result to rank zero, intervar [0,radius]
   MpiIntegrator.Integrate(0, 0, radius);

   
   /*****************************************************************************************/
   /*this is the root process,calculate the surface z=+qsrt(1−(x2+b2)) for x∈[−1,0],y∈[−1,1]*/
   /*and recv the other surface regions in TPolyMarker3D object from others ranks           */
   /*****************************************************************************************/
   if (world.Rank() == 0) {
      gBenchmark->Start("sphere");
      // create and open a canvas
      TCanvas *canvas = new TCanvas("sphere", "sphere", 300, 10, 700, 500);
      canvas->SetFillColor(14);
      // creating view
      TView *view = TView::CreateView(1, 0, 0);
      view->RotateView(90,-10);
      view->SetRange(-radius*2, -radius*2, -2*radius, 2*radius, 2*radius, 2*radius);
      // create a PolyMarker3D for region/segment 1 calculate in rank 0(this rank)
      TPolyMarker3D *surfaceSeg13d = new TPolyMarker3D();
      
      for (theta = 0; theta <0.5*M_PI ; theta+= delta) {
         for (phi = 0; phi <M_PI; phi += delta) {
	   
	   x=radius*cos(theta)*sin(phi);
	   y=radius*sin(theta)*sin(phi);
	   z=radius*cos(phi);
            // set points
            surfaceSeg13d->SetPoint(point, x, y, z);
            point++;
         }
      }
      // set marker size, color & style
      surfaceSeg13d->SetMarkerSize(1);
      surfaceSeg13d->SetMarkerColor(world.Rank() + world.Size());
      surfaceSeg13d->SetMarkerStyle(1);
      //drawing region 1
      surfaceSeg13d->Draw();

      //getting region 2 from rank 1
      TPolyMarker3D *surfaceSeg23d = new TPolyMarker3D;
      world.RecvObject(1, 1, *surfaceSeg23d);
      surfaceSeg23d->Draw();

      //getting region 3 from rank 2
      TPolyMarker3D *surfaceSeg33d = new TPolyMarker3D;
      world.RecvObject(2, 1, *surfaceSeg33d);
      surfaceSeg33d->Draw();

      //getting region 4 from rank 3
      TPolyMarker3D *surfaceSeg43d = new TPolyMarker3D;
      world.RecvObject(3, 1, *surfaceSeg43d);
      surfaceSeg43d->Draw();

      Char_t timeStr[60];
      Char_t pointsStr[100];
      Char_t radiusStr[100];
      Char_t deltaStr[100];
      gBenchmark->Show("sphere");

      Float_t ct = gBenchmark->GetCpuTime("sphere");
      sprintf(timeStr, "Execution time: %g sec.", ct);
      sprintf(radiusStr, "Radius : %g", radius);
      sprintf(deltaStr, "Delta : %g", delta);
      sprintf(pointsStr, "Total surface points : %g", point*4);
      
      TPaveText *text = new TPaveText(0.1, 0.6, 0.9, 0.99);
      text->SetFillColor(42);
      text->AddText("ROOTMpi example: sphere.C");
      text->AddText(timeStr);
      text->AddText(radiusStr);
//.........这里部分代码省略.........

//*************************************************************
void arrangeBiasCanvas(TCanvas *canv,TH1F* dxyPhiMeanTrend[100],TH1F* dzPhiMeanTrend[100],TH1F* dxyEtaMeanTrend[100],TH1F* dzEtaMeanTrend[100],Int_t nFiles, TString LegLabels[10]){
//*************************************************************

  TLegend *lego = new TLegend(0.19,0.70,0.79,0.92);
  lego-> SetNColumns(2);
  lego->SetFillColor(10);
  lego->SetTextSize(0.042);
  lego->SetTextFont(42);
  lego->SetFillColor(10);
  lego->SetLineColor(10);
  lego->SetShadowColor(10);

  TPaveText *pt = new TPaveText(0.13,0.95,0.89,0.97,"NDC");
  pt->SetFillColor(10);
  pt->SetTextColor(1);
  pt->SetTextFont(42);
  pt->SetTextAlign(11);
  TText *text1 = pt->AddText("CMS Preliminary 2015 - 3.8T collision data");
  text1->SetTextSize(0.05);
 
  TPaveText *pt2 = new TPaveText(0.70,0.75,0.89,0.92,"NDC");
  pt2->SetFillColor(10);
  pt2->SetTextColor(kBlue);
  pt2->SetTextFont(62);
  pt2->SetTextAlign(11);
  TText *text2 = pt2->AddText("run 247078");
  text2->SetTextSize(0.05);

  canv->SetFillColor(10);  
  canv->Divide(2,2);
 
  canv->cd(1)->SetBottomMargin(0.14);
  canv->cd(1)->SetLeftMargin(0.18);
  canv->cd(1)->SetRightMargin(0.01);
  canv->cd(1)->SetTopMargin(0.06);  

  canv->cd(2)->SetBottomMargin(0.14);
  canv->cd(2)->SetLeftMargin(0.18);
  canv->cd(2)->SetRightMargin(0.01);
  canv->cd(2)->SetTopMargin(0.06);  
  
  canv->cd(3)->SetBottomMargin(0.14);
  canv->cd(3)->SetLeftMargin(0.18);
  canv->cd(3)->SetRightMargin(0.01);
  canv->cd(3)->SetTopMargin(0.06);  

  canv->cd(4)->SetBottomMargin(0.14);
  canv->cd(4)->SetLeftMargin(0.18);
  canv->cd(4)->SetRightMargin(0.01);
  canv->cd(4)->SetTopMargin(0.06); 

  TH1F *dBiasTrend[4][nFiles]; 
  
  for(Int_t i=0;i<nFiles;i++){
    dBiasTrend[0][i] = dxyPhiMeanTrend[i];
    dBiasTrend[1][i] = dzPhiMeanTrend[i];
    dBiasTrend[2][i] = dxyEtaMeanTrend[i];
    dBiasTrend[3][i] = dzEtaMeanTrend[i];
  }

  Double_t absmin[4]={999.,999.,999.,999.};
  Double_t absmax[4]={-999.,-999.-999.,-999.};

  for(Int_t k=0; k<4; k++){

    canv->cd(k+1);
    
    for(Int_t i=0; i<nFiles; i++){
      if(dBiasTrend[k][i]->GetMaximum()>absmax[k]) absmax[k] = dBiasTrend[k][i]->GetMaximum();
      if(dBiasTrend[k][i]->GetMinimum()<absmin[k]) absmin[k] = dBiasTrend[k][i]->GetMinimum();
    }

   
    Double_t safeDelta=(absmax[k]-absmin[k])/8.;
    Double_t theExtreme=std::max(absmax[k],TMath::Abs(absmin[k]));

    for(Int_t i=0; i<nFiles; i++){
      if(i==0){
	//dBiasTrend[i]->GetYaxis()->SetRangeUser(absmin-safeDelta/2.,absmax+safeDelta);
	//std::cout<<"name is: "<< dBiasTrend[k][i]->GetName() <<std::endl;
	TString theTitle = dBiasTrend[k][i]->GetName();
	if( theTitle.Contains("Norm")){
	  dBiasTrend[k][i]->GetYaxis()->SetRangeUser(std::min(-0.48,absmin[k]-safeDelta/2.),std::max(0.48,absmax[k]+safeDelta/2.));
	} else {
	  //dBiasTrend[k][i]->GetYaxis()->SetRangeUser(std::min(-8.8,absmin[k]-safeDelta/2.),std::max(8.8,absmax[k]+safeDelta/2.));
	  dBiasTrend[k][i]->GetYaxis()->SetRangeUser(-theExtreme-(safeDelta/2.),theExtreme+(safeDelta/2.));
	  //dBiasTrend[k][i]->GetYaxis()->SetRangeUser(-theExtreme,theExtreme);
	} 
	dBiasTrend[k][i]->Draw("Le1");
      }
      else dBiasTrend[k][i]->Draw("Le1sames");
      if(k==0){
	lego->AddEntry(dBiasTrend[k][i],LegLabels[i]);
      } 
    }  
  
    lego->Draw();
    pt->Draw("same");
  }
//.........这里部分代码省略.........