C++ TGraph::Write方法代码示例

void processmarriage(const char* ifname="marry.txt",const char* ofname="marry.root"){

  vector<int> mage;
  vector<float> mpercent;

  vector<int> fage;
  vector<float> fpercent;

  mage.push_back(0);
  mpercent.push_back(0.0);
  fage.push_back(0);
  fpercent.push_back(0.0);

  ifstream marry(ifname);
  while(1){
    int gender,year;
    float percent;
    marry >> gender >> year >> percent;
    if(!marry.good())break;
    if(gender == 1){
      mage.push_back(year);
      mpercent.push_back(percent);
    }
    if(gender == 0){
      fage.push_back(year);
      fpercent.push_back(percent);
    }
  }

  TFile* outfile = new TFile(ofname,"RECREATE");

  TGraph* man = new TGraph(mage.size());
  TGraph* woman = new TGraph(fage.size());

  for(int i = 0; i < (int)mage.size(); i++){
    man->SetPoint(i,mage[i],mpercent[i]);
  }
  for(int i = 0; i < (int)fage.size(); i++){
    woman->SetPoint(i,fage[i],fpercent[i]);
  }

  man->SetName("man");
  woman->SetName("woman");
  man->Write();
  woman->Write();
  outfile->Write();
  outfile->Close();

}

void createInputs(int n = 2) 
{
   for(UInt_t i = 0; i < (UInt_t)n; ++i ) {
      TFile *file = TFile::Open(TString::Format("input%d.root",i),"RECREATE");
      TH1F * h = new TH1F("h1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
 
      Int_t nbins[5];
      Double_t xmin[5];
      Double_t xmax[5];
      for(UInt_t j = 0; j < 5; ++j) {
         nbins[j] = 10; xmin[j] = 0; xmax[j] = 10;
      }
      THnSparseF *sparse = new THnSparseF("sparse", "sparse", 5, nbins, xmin, xmax);
      Double_t coord[5] = {0.5, 1.5, 2.5, 3.5, 4.5};
      sparse->Fill(coord);
      sparse->Write();
      
      THStack *stack = new THStack("stack","");
      h = new TH1F("hs_1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
      h->SetDirectory(0);
      stack->Add(h);
      h = new TH1F("hs_2","",10,0,100);
      h->Fill(30.5); h->Fill(40.5);
      h->SetDirectory(0);
      stack->Add(h);
      stack->Write();

      TGraph *gr = new TGraph(3);
      gr->SetName("exgraph");
      gr->SetPoint(0,1,1);
      gr->SetPoint(1,2,2);
      gr->SetPoint(2,3,3);
      
      gr->Write();
      
      TTree *tree = new TTree("tree","simplistic tree");
      Int_t data = 0;
      tree->Branch("data",&data);
      for(Int_t l = 0; l < 2; ++l) {
         data = l;
         tree->Fill();
      }
      
      file->Write();
      delete file;
   }
}

//.........这里部分代码省略.........
							//hMuonCutQDC[filesScanned][k]->Fill(QDC[k]);
							hMuonManyQDC[k]->Fill(QDC[k]);
							MuonQDCtot += QDC[k];
						}
					}			
					AvgMuonQDCvalue = MuonQDCtot/((double) muonnumPanelsHit);
					hAvgMuonQDC->Fill(AvgMuonQDCvalue);
					MuonQDCtot = 0;
					
				}
				
				//=====================END ACTUAL GODDAMMED ANALYSIS===================
			
			lednumPanelsHit=0;
			muonnumPanelsHit=0;

		}	// End loop over VetoTree entries.

		// === END OF FILE Output & Plotting ===
		if (filesScanned == 0){
			TDirectory *runmultiplicity = RootFile->mkdir("RunMultiplicity");	
			TDirectory *corruptionintime = RootFile->mkdir("CorruptionInTime");
		}
		
		corruption = ((Float_t)BadTSInFile/nentries)*100;
		printf(" Corrupted scaler entries: %i of %lli, %.3f %%.\n",BadTSInFile,nentries,corruption);
		if(run>45000000) SCorruption->SetPoint(filesScanned,run-45000000,corruption);
		else SCorruption->SetPoint(filesScanned,run,corruption);

		if (PlotCorruptedEntries) {
			RootFile->cd("CorruptionInTime");
			char outfile1[200];	
			sprintf(outfile1,"CorruptionInTime_Run%i",run);
			CorruptionInTime->Write(outfile1,TObject::kOverwrite);
			RootFile->cd();
		}
		if (PlotMultiplicity) {
			RootFile->cd("RunMultiplicity");
			char outfile2[200];
			sprintf(outfile2,"Multiplicity_Run%i",run);
			OneRunMultiplicity->Write(outfile2,TObject::kOverwrite);
			RootFile->cd();
		}

		if (!IsEmpty && ledcount[filesScanned] > 2){
			ledTime[filesScanned]->Fit("pol1");
			slope[filesScanned] = ledTime[filesScanned]->GetFunction("pol1")->GetParameter(1);
		}
		
		for (int j=0; j<numFiles; ++j){
			ledFiletime->SetPoint(j, runnum[j], ledcount[j]);
		}

		// ==========================

		delete VetoTree;
		cout << "files Scanned: " << filesScanned << endl;
		if (duration >= 0){
			totdur += duration;
		}
		if (duration < 0){
			omitteddur[baddurcount]=run;
			baddurcount +=1;
		}
		cout << "totdur = " << totdur << " | duration = " << duration << endl;
		filesScanned++;

//.........这里部分代码省略.........
y[1]=1+s_200/ratio->Eval(200);
y[2]=1+s_400/ratio->Eval(400);
y[3]=1+s_600/ratio->Eval(600);
y[4]=1+s_800/ratio->Eval(800);
y[5]=1+s_1000/ratio->Eval(1000);

ym[0]=-s_100/m_100;
ym[1]=-s_200/m_200;
ym[2]=-s_400/m_400;
ym[3]=-s_600/m_600;
ym[4]=-s_800/m_800;
ym[5]=-s_1000/m_1000;


TGraph* g = new TGraph(6,x,y);
TGraph* gm = new TGraph(6,x,ym);
TGraph* gup = new TGraph(6,x,yup);
TGraph* gdown = new TGraph(6,x,ydown);


TCanvas *c3 = new TCanvas("c3","c3",1000,1000);
c3->cd();

//gup->Draw("AC*");
//gdown->Draw("C*");
g->Draw("AC*");

gPad->SetBottomMargin(0.2);
gPad->SetLeftMargin(0.2);
gStyle->SetOptStat(0);

g->GetXaxis()->SetTitle("M_{Zbb}");
g->GetXaxis()->SetRangeUser(50,1100);
g->GetYaxis()->SetLabelSize(0.06);
g->GetYaxis()->SetTitle("Uncertainty");
g->GetYaxis()->SetTitleSize(0.06);
g->GetYaxis()->SetTitleOffset(1.4);
g->GetXaxis()->SetLabelSize(0.06);
g->GetXaxis()->SetTitleSize(0.06);
g->GetXaxis()->SetTitleOffset(1);
g->GetYaxis()->SetNdivisions(5);

TFile f("syst_zxx.root","recreate");
g->Write();
f.Close();


//gm->Draw("C*");
//g->SetMaximum(1);
//g->SetMinimum(-1);
//h2c->Draw("same");


TH1D *h22=h2->Clone();

TCanvas *c5 =  new TCanvas("c5","c5",1000,1000);

gPad->SetBottomMargin(0.2);
gPad->SetLeftMargin(0.2);
gStyle->SetOptStat(0);

h22->Draw();
h22->GetXaxis()->SetRangeUser(50,1100);
h22->GetYaxis()->SetLabelSize(0.06);
h22->GetYaxis()->SetTitleSize(0.06);
h22->GetYaxis()->SetTitleOffset(1.4);
h22->GetXaxis()->SetLabelSize(0.06);
h22->GetXaxis()->SetTitleSize(0.06);
h22->GetXaxis()->SetTitleOffset(1);

ratio->SetLineColor(kRed);
ratio->Draw("same");

gup->Draw("C");
gdown->Draw("C");


TLegend *leg = new TLegend(0.6,0.7,0.89,0.89);
leg->SetLineColor(0);
leg->SetFillColor(0);
leg->AddEntry(h22,"[email protected] / MG5","lep");
leg->AddEntry(ratio,"best fit","l");
leg->AddEntry(gup,"Syst Error (#pm 1 #sigma)","l");
leg->Draw();



TCanvas *c4 = new TCanvas("c4","c4",1000,1000);
c4->Divide(2,2);
c4->cd(1);
histp0->Draw();
c4->cd(2);
histp1->Draw();
c4->cd(3);
histp2->Draw();
c4->cd(4);
histp3->Draw();


}

//.........这里部分代码省略.........
      TH1F *innereta = (TH1F *)MyFile->Get("innereta");
      TFile *rebinfile = new TFile("rebinout1.root", "Recreate"); // 80 TO 100 PT AVERAGE
    }
    if (filecounter == 3) {
      TFile *MyFile = new TFile("balanceout.root", "Read"); // 100 TO 140 PT AVERAGE
      if (MyFile->IsOpen()) printf("File Opened Successfully.\n");
      gFile = MyFile;
      MyFile->ls();
      TH1F *balance = (TH1F *)MyFile->Get("balance");
      TH1F *outereta = (TH1F *)MyFile->Get("outereta");
      TH1F *innereta = (TH1F *)MyFile->Get("innereta");
      TFile *rebinfile = new TFile("rebinout.root", "Recreate"); // 100 TO 140 PT AVERAGE
    }
    if (filecounter == 4) {
      TFile *MyFile = new TFile("balanceout12.root", "Read"); // 140 TO 200 PT AVERAGE
      if (MyFile->IsOpen()) printf("File Opened Successfully.\n");
      gFile = MyFile;
      MyFile->ls();
      TH1F *balance = (TH1F *)MyFile->Get("balance");
      TH1F *outereta = (TH1F *)MyFile->Get("outereta");
      TH1F *innereta = (TH1F *)MyFile->Get("innereta");
      TFile *rebinfile = new TFile("rebinout12.root", "Recreate"); // 140 TO 200 PT AVERAGE
    }
 
  Float_t  cms_hcal_edge_pseudorapidity[83] = {-5.191, -4.889, -4.716, -4.538, -4.363, -4.191, -4.013, -3.839, -3.664, -3.489, -3.314, -3.139, -2.964, -2.853, -2.650, -2.500, -2.322, -2.172, -2.043, -1.930, -1.830, -1.740, -1.653, -1.566, -1.479, -1.392, -1.305, -1.218, -1.131, -1.044, -0.957, -0.879, -0.783, -0.696, -0.609, -0.522, -0.435, -0.348, -0.261, -0.174, -0.087, 0.000, 0.087,  0.174,  0.261,  0.348,  0.435,  0.522,  0.609, 0.696,  0.783,  0.879,  0.957,  1.044,  1.131,  1.218, 1.305,  1.392,  1.479,  1.566,  1.653,  1.740,  1.830, 1.930,  2.043,  2.172,  2.322,  2.500,  2.650,  2.853, 2.964,  3.139,  3.314,  3.489,  3.664,  3.839,  4.013, 4.191,  4.363,  4.538,  4.716,  4.889,  5.191};
  Float_t xcoordinate[82];
  Float_t intbalance[82];
  Float_t intinnereta[82];
  Float_t intoutereta[82];
 

  for (int count = 0; count < 82; count++) {
    Float_t lowerbound = cms_hcal_edge_pseudorapidity[count];
    Float_t upperbound = cms_hcal_edge_pseudorapidity[count + 1]; 
    Int_t lowerbin = balance->FindBin(lowerbound);
    Int_t upperbin = balance->FindBin(upperbound);
    Float_t integral = balance->Integral(lowerbin, upperbin);
    Float_t integral1 = innereta->Integral(lowerbin, upperbin);
    Float_t integral2 = outereta->Integral(lowerbin, upperbin);
    xcoordinate[count] = (lowerbound + upperbound)/2; 
    intbalance[count] = integral;
    //cout << "My x-coordinate is " << xcoordinate[count] << " and my integral is " << intbalance[count] << endl;
    intinnereta[count] = integral1;
    intoutereta[count] = integral2;
    //balancerebin->Fill(integral);
    //inneretarebin->Fill(integral1);
    //outeretarebin->Fill(integral2);
  }
TGraph *balancerebin = new TGraph(82, xcoordinate, intbalance);
TGraph *inneretarebin = new TGraph(82, xcoordinate, intinnereta);
TGraph *outeretarebin = new TGraph(82, xcoordinate, intoutereta);

balancerebin->SetTitle("Rebin Full Eta Range");
balancerebin->GetXaxis()->SetTitle("Dijet Balance");
balancerebin->GetYaxis()->SetTitle("Event Fraction");
balancerebin->SetMarkerStyle(8);
balancerebin->SetLineColor(0);
inneretarebin->SetTitle("Rebin Inner Eta Range");
inneretarebin->GetXaxis()->SetTitle("Dijet Balance");
inneretarebin->GetYaxis()->SetTitle("Event Fraction");
inneretarebin->SetMarkerStyle(8);
inneretarebin->SetLineColor(0);
outeretarebin->SetTitle("Rebin Outer Eta Range");
outeretarebin->GetXaxis()->SetTitle("Dijet Balance");
outeretarebin->GetYaxis()->SetTitle("Event Fraction");
outeretarebin->SetMarkerStyle(8);
outeretarebin->SetLineColor(0);

//balancerebin->Draw("ACP");
//inneretarebin->Draw("ACP");
//outeretarebin->Draw("ACP");
  
TCanvas *brebin = new TCanvas(1);
brebin->cd();
balancerebin->Draw("ACP");

TCanvas *binnerrebin = new TCanvas(2);
binnerrebin->cd();
inneretarebin->Draw("ACP");

TCanvas *bouterrebin = new TCanvas(3);
bouterrebin->cd();
outeretarebin->Draw("ACP");
  //Double_t scalerebin = 1/balancerebin->Integral();
  //Double_t scaleinnerrebin = 1/inneretarebin->Integral();
  //Double_t scaleouterrebin = 1/outeretarebin->Integral();

  //balancerebin->Scale(scalerebin);
  //inneretarebin->Scale(scaleinnerrebin);
  //outeretarebin->Scale(scaleouterrebin);
  
  rebinfile->cd();
  balancerebin->Write();
  inneretarebin->Write();
  outeretarebin->Write();
  rebinfile->Close();
  MyFile->Close();

}
} //END FILE LOOP

//.........这里部分代码省略.........
  TCanvas *can3 = new TCanvas("can3","can3",1200,800);
  can3->Divide(3,2);

  t->SetTextSize(0.07);

  gr->SetMarkerColor(6);
  gr_obsp1->SetMarkerColor(6);
  gr_obsm1->SetMarkerColor(6);
  gr_exp->SetMarkerColor(6);
  gr_expp1->SetMarkerColor(6);
  gr_expm1->SetMarkerColor(6);

  can3->cd(1);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl->GetYaxis()->SetRangeUser(0,200);
  hexcl->Draw("colz");
  gr->Draw("lp");
  t->DrawLatex(0.3,0.8,"observed");

  can3->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl_obsp1->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl_obsp1->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl_obsp1->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl_obsp1->GetYaxis()->SetRangeUser(0,200);
  hexcl_obsp1->Draw("colz");
  gr_obsp1->Draw("lp");
  t->DrawLatex(0.3,0.8,"observed (+1#sigma)");

  can3->cd(3);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl_obsm1->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl_obsm1->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl_obsm1->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl_obsm1->GetYaxis()->SetRangeUser(0,200);
  hexcl_obsm1->Draw("colz");
  gr_obsm1->Draw("lp");
  t->DrawLatex(0.3,0.8,"observed (-1#sigma)");

  can3->cd(4);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl_exp->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl_exp->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl_exp->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl_exp->GetYaxis()->SetRangeUser(0,200);
  hexcl_exp->Draw("colz");
  gr_exp->Draw("lp");
  t->DrawLatex(0.3,0.8,"expected");

  can3->cd(5);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl_expp1->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl_expp1->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl_expp1->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl_expp1->GetYaxis()->SetRangeUser(0,200);
  hexcl_expp1->Draw("colz");
  gr_expp1->Draw("lp");
  t->DrawLatex(0.3,0.8,"expected (+1#sigma)");

  can3->cd(6);
  gPad->SetGridx();
  gPad->SetGridy();
  hexcl_expm1->GetXaxis()->SetTitle("stop mass [GeV]");
  hexcl_expm1->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
  hexcl_expm1->GetXaxis()->SetRangeUser(xaxismin,600);
  hexcl_expm1->GetYaxis()->SetRangeUser(0,200);
  hexcl_expm1->Draw("colz");
  gr_expm1->Draw("lp");
  t->DrawLatex(0.3,0.8,"expected (-1#sigma)");
*/

  if( print ){
    if     ( TString(sample).Contains("T2tt") ){
      can1->Print("../../plots/combinePlots_T2tt.pdf");
      //can2->Print("../../plots/combinePlots_T2tt_bestSignalRegion.pdf");
      //can3->Print("../../plots/combinePlots_T2tt_excludedPoints.pdf");
    }
    else if( TString(sample).Contains("T2bw") ){
      can1->Print(Form("../../plots/combinePlots_T2bw_x%i.pdf",x));
      //can2->Print(Form("../../plots/combinePlots_T2bw_x%i_bestSignalRegion.pdf",x));
      //can3->Print(Form("../../plots/combinePlots_T2bw_x%i_excludedPoints.pdf",x));
    }
  }

  TFile* fout = TFile::Open(Form("%s_x%icombinePlots.root",sample,x),"RECREATE");
  fout->cd();
  hxsec_best->Write();
  hxsec_best_exp->Write();
  gr->Write();
  fout->Close();

}

//______________________________________________________________________________
void CheckTime(const Char_t* loc)
{
    // Check time.

    Char_t t[256];

    // number of runs
    Int_t nRuns = gFiles->GetSize();

    // create arrays
    const Int_t nCh = 352;
    Double_t** pedPos = new Double_t*[nCh];
    Double_t* runNumbersD = new Double_t[nRuns];
    for (Int_t i = 0; i < nCh; i++) pedPos[i] = new Double_t[nRuns];

    // open the output files
    TFile* fROOTout = new TFile("/tmp/tagger_time.root", "RECREATE");

    // create directories
    for (Int_t i = 0; i < nCh; i++)
    {
        sprintf(t, "%03d", i);
        fROOTout->mkdir(t);
    }
    
    TF1* func = new TF1("func", "gaus(0)+pol1(3)", -1000 , 1000);

    // loop over runs
    for (Int_t i = 0; i < nRuns; i++)
    {   
        // get the file
        TFile* f = (TFile*) gFiles->At(i);

        // extract run number
        Int_t runNumber;
        sprintf(t, "%s/ARHistograms_CB_%%d.root", loc);
        sscanf(f->GetName(), t, &runNumber);
        runNumbersD[i] = (Double_t)runNumber;

        printf("Processing run %d (%d/%d)\n", runNumber, i+1, nRuns);

        fROOTout->cd();
        
        TH2* h2 = (TH2*) f->Get("CaLib_Tagger_Time");

        // loop over channels
        for (Int_t j = 0; j < nCh; j++)
        {
            // load histogram
            sprintf(t, "%03d", j);
            fROOTout->cd(t);
            sprintf(t, "%d_%d", i, j);
            TH1* h = h2->ProjectionX(t, j+1, j+1);
            h->Rebin(2);
            
            sprintf(t, "Run_%d", runNumber);
            TCanvas* c = new TCanvas(t, t);
            TLine* tline = 0;
            
            // check entries
            if (h->GetEntries())
            {
                // fit gaussian to peak
                Double_t maxPos = h->GetXaxis()->GetBinCenter(h->GetMaximumBin());
                func->SetParameters(1, maxPos, 0.5, 1, 0.1);
                func->SetRange(maxPos - 2, maxPos + 2);
                func->SetParLimits(0, 0, 1000);
                for (Int_t k = 0; k < 10; k++)
                    if (!h->Fit(func, "RBQ")) break;

                // save position in file and memory
                maxPos = func->GetParameter(1);
                pedPos[j][i] = maxPos;

                h->GetXaxis()->SetRangeUser(maxPos - 10, maxPos + 10);
                h->Draw();
                
                tline = new TLine(maxPos, 0, maxPos, 10000);
                tline->SetLineColor(kRed);
                tline->SetLineWidth(2);
                tline->Draw();
            }
            else
            {
                h->Draw();
            }

            c->Write(c->GetName(), TObject::kOverwrite);
    
            delete h;
            delete c;
            if (tline) delete tline;
        }
        
        delete h2;
    }
    
    // create pedestal evolution graphs
    fROOTout->cd();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
//         message = TString::Format("unnormalize run bgsubbed %s", pol.Data());
//          if (OkToContinue(message))
//          {
//              f = new TFile(combined_hist_run.Data(),"UPDATE");
//
//              std::cout << "\tUnnormalizing adc hists" << std::endl;
//              TDirectory* dir = f->GetDirectory("bgsubbed_shifted_ltf_corr_adc");
//              if (dir==0)  std::cout << "Cannot find folder containing hists" << std::endl;
//              else         UnnormalizeAllHists(dir,16042.0);
//
//              std::cout << "\tUnnormalizing adc_gt_thresh hists" << std::endl;
//              dir = f->GetDirectory("bgsubbed_shifted_ltf_corr_adc_gt_thresh");
//              if (dir==0)  std::cout << "Cannot find folder containing hists" << std::endl;
//              else         UnnormalizeAllHists(dir,16042.0);
//
//              f->Close();
//          }
     }

     if (OkToContinue("compute systematic unc. from overnight bgnd histograms"))
     {
         TString uncorr_dir = "ltf_corr_adc_gt_thresh";
         TString corr_dir = "bgsubbed_shifted_ltf_corr_adc_gt_thresh";
         f = new TFile(combined_hist_run.Data(),"update");
         TDirectory* du = f->GetDirectory(uncorr_dir.Data());
         TDirectory* ds = f->GetDirectory(corr_dir.Data());

         if (du==0)      std::cout << "Cannot find " << uncorr_dir << std::endl;
         else if (ds==0) std::cout << "Cannot find " << corr_dir << std::endl;
         else
         {
             TGraph* gr = GenerateSubtractionUncertainties(du, ds);
             f->cd();
             gr->Write("", TObject::kOverwrite);
         }

         f->Close();
     }

//     if (OkToContinue("Form alpha normalized spectra"))
//     {
//         f = new TFile(combined_hist_run.Data(),"update");
//
//         TString uncorr_dir = "ltf_corr_adc_gt_thresh";
//         TString corr_dir = "bgsubbed_shifted_ltf_corr_adc_gt_thresh";
//         TDirectory* du = f->GetDirectory(uncorr_dir.Data());
//
//
//         if (du==0)      std::cout << "Cannot find " << uncorr_dir << std::endl;
//         else
//         {
//            CreateAlphaNormedHists(du, "g4_sa_corr_232Th.dat");
//         }
//
//         TDirectory* dc = f->GetDirectory(corr_dir);
//
//         if (dc==0)      std::cout << "Cannot find " << corr_dir << std::endl;
//         else
//         {
//            CreateAlphaNormedHists(dc, "g4_sa_corr_232Th.dat");
//         }
//
//         f->Close();
//     }

    if (OkToContinue("integrate all \"adc_gt_thresh_tofcut\" histograms"))

//.........这里部分代码省略.........
 MyTreeMC->SetBranchAddress("p",&p);
 MyTreeMC->SetBranchAddress("eleES",&eleES);
 MyTreeMC->SetBranchAddress("eleFBrem",&eleFBrem);
 
 numEntriesMC = MyTreeMC->GetEntries();
 
 ///==== prepare minuit ====

 fitMin->SetRange(MinScanRange,MaxScanRange); 
 
 double step[1] = {0.001};
 double variable[1] = {0.0};
 minuit->SetLimitedVariable(0,"Scale" , variable[0]  , step[0] , MinScan  , MaxScan );
 

///===========================
///==== DATA Scale search ====
 ScaleTrue = -1000; ///==== default
 Data_or_MC = 1; ///=== 1 = Data;  0 = MC;
 numEvents = MyTreeDATA->GetEntries(); //==== number of events in Data sample
 outFile->cd();
 vET_data.clear();
 nIter = 1000000000; ///==== less than 1000000000 iterations at the end !!!
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
 
 MyTreeDATA->Draw(">> myList",(AdditionalCut + Form(" && ET > %f",minET)).Data(),"entrylist");
 TEntryList *mylist = (TEntryList*)gDirectory->Get("myList");
 MyTreeDATA->SetEntryList(mylist);
 
 MyTreeDATA->Draw(Form("%s >> %s",variableName.c_str(),nameDATA.Data()));
 ConvertStdVectDouble(vET_data,MyTreeDATA->GetV1(),mylist->GetN());
 
 hDATA.Write();
  
 std::cerr << "... I'm minimizing ... DATA analysis" << std::endl;
 std::cerr << ">>>>>>> numEvents = " << numEvents << " => " << vET_data.size() << " selected (=" << mylist->GetN() << ")" << std::endl;
 numSelectedData = vET_data.size();
 
 
 ///===== Chi2 ====
 std::cerr << " === Chi2 === " << std::endl;
 minuit->SetFunction(functorChi2);

 TGraph * grChi2 = new TGraph(iNoSteps);
 minuit->Scan(iPar_NoBG,iNoSteps,grChi2->GetX(),grChi2->GetY(),MinScan,MaxScan);

// TGraph * grChi2 = new TGraph();
// for (int iStep = 0; iStep < iNoSteps; iStep++){
//  double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
//  double y = Chi2F(&x);
//  grChi2->SetPoint(iStep+1,x,y);
// }
 grChi2->Draw("AL");
 outFile->cd();
 minuit->PrintResults();
 outFile->cd();
 grChi2->SetTitle("grChi2");
 grChi2->Write();
 const double *outParametersTemp = minuit->X();
 const double *errParametersTemp = minuit->Errors();
  
 double *outParameters = new double;
 double *errParameters = new double;
 outParameters[0] = outParametersTemp[0];
 errParameters[0] = errParametersTemp[0];

//.........这里部分代码省略.........
  for (int iEvt = 0; iEvt < numSelectedData; iEvt ++){
   listMCHere->Enter(myListMC->GetEntry(gRandom->Uniform(0,myListMC->GetN())));
  }
  
  MyTreeMC->SetEntryList(listMCHere);
  MyTreeMC->Draw(Form("(1+%f) * %s >> %s",ScaleTrue,variableName.c_str(),nameDATA.Data()));
  
  ConvertStdVectDouble(vET_data,MyTreeMC->GetV1(),numSelectedData);
  
  ///==== likelihood ====
 std::cerr << " === LL === " << std::endl;
 std::cerr << " === pseudo vET_data.size() = " << vET_data.size() << std::endl;
 
 
  minuit->SetFunction(functorLL); 
  TGraph * grLL_temp = new TGraph(iNoSteps);
  minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
  TGraph * grLL = new TGraph();
  int nPointLL = 0;
 for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
   double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
   double y = LLFunc(&x);
   if (y != numberDATA * numEvents) {
    grLL->SetPoint(nPointLL,x,y);
    nPointLL++;
   }
  }
  grLL->Draw("AL");
  outFile->cd();
  minuit->PrintResults();
  const double *outParametersTemp2 = minuit->X();
  const double *errParametersTemp2 = minuit->Errors();
  
  double *outParametersLL = new double;
  double *errParametersLL = new double;
  outParametersLL[0] = outParametersTemp2[0];
  errParametersLL[0] = errParametersTemp2[0];
  
  std::cerr << " nPointLL = " << nPointLL << std::endl;
  
  double minLL = grLL->Eval(outParametersLL[0]);
  ///==== end likelihood ====
  
    ///==== Save the whole shape of LL/Chi2 ====
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   
   Alpha   = X_ii;
   Chi2    = 0;
   LL      = grLL->Eval(X_ii);
   NewChi2 = 0;
   
   myTreeChi2->Fill();
  }

  ///===== Look for minima =====
  double a;
  double b;
  double c;
  
  
  double errX_low = -9999;
  double errX_up = 9999;
  int err_low = 0;
  int err_up = 0;
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   double here = grLL->Eval(X_ii);
   if (err_low == 0){
    if (here < (minLL + DELTA_LL)){
     errX_low = X_ii;
     err_low = 1;
    }
   }
   else if (err_up == 0 && here > (minLL + DELTA_LL) && X_ii > outParametersLL[0]){
    errX_up = X_ii; 
    err_up = 1;
   }
  }
  
  AlphaMean = outParametersLL[0];
  AlphaMinus = errX_low;
  AlphaPlus = errX_up;
 
  grLL->Fit("fitMin","RMQ");
  c = fitMin->GetParameter(0);
  b = fitMin->GetParameter(1);
  a = fitMin->GetParameter(2);
  AlphaMean_Fit  = -b / (2*a);
  AlphaMinus_Fit = (-b + sqrt(2*a)) / (2*a);  ///==== delta LL = 0.5
  AlphaPlus_Fit  = (-b - sqrt(2*a)) / (2*a);  ///==== delta LL = 0.5   
  
  myTreeLL_Result->Fill();
 
  grLL->Write();
  
  //delete listMCHere;
  
 }
}

//.........这里部分代码省略.........
      nBins = hDen1DvsTimeOld->GetNbinsX()+1;
      Float_t binwidth =  (Time - hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1))/(nBins-1);
      edge0 = hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1) - binwidth/2.;
      edge1 = Time + binwidth/2.;
    }
    hDen1DvsTime = new TH2F("temp","",nBins,edge0,edge1,
		       	hDen1D->GetNbinsX(),
			hDen1D->GetBinLowEdge(1),
			hDen1D->GetBinLowEdge(hDen1D->GetNbinsX()+1));
    
    for(Int_t ix=1;ix<hDen1DvsTime->GetNbinsX();ix++) {
      for(Int_t iy=1;iy<hDen1DvsTime->GetNbinsY();iy++) {
	hDen1DvsTime->SetBinContent(ix,iy,hDen1DvsTimeOld->GetBinContent(ix,iy));
      }
    }  
    delete hDen1DvsTimeOld;
  
    // Fill last bin with the newest values.
    for(Int_t iy=1;iy<=hDen1D->GetNbinsX();iy++) {
      hDen1DvsTime->SetBinContent(nBins,iy,hDen1D->GetBinContent(iy));
    }   

    hDen1DvsTime->GetZaxis()->SetTitle("n_{b}/n_{0}");
    hDen1DvsTime->GetYaxis()->SetTitle("k_{p}#zeta");
    hDen1DvsTime->GetXaxis()->SetTitle("k_{p}z");
    hDen1DvsTime->GetZaxis()->CenterTitle();
    hDen1DvsTime->GetYaxis()->CenterTitle();
    hDen1DvsTime->GetXaxis()->CenterTitle();
    hDen1DvsTime->SetName(hName);

    // Change the range of z axis 
    Float_t Denmax = hDen1DvsTime->GetMaximum();
    hDen1DvsTime->GetZaxis()->SetRangeUser(0,Denmax); 
    hDen1DvsTime->Write(hName,TObject::kOverwrite);

  }

  // RMS (vs z) of the beam's charge distribution: 
  TProfile *hDen2Dprof = NULL;
  TH1F *hRms = NULL;
  Double_t axisPos = x2Mid;
  if(hDen2D) {
    TString pname = hDen2D->GetName();
    pname += "_pfx";
    
    hDen2Dprof =  (TProfile*) gROOT->FindObject(pname.Data());
    if(hDen2Dprof) { delete hDen2Dprof; hDen2Dprof = NULL; }
    hDen2Dprof = hDen2D->ProfileX("_pfx",1,-1,"s");
    
    hRms = (TH1F*) gROOT->FindObject("hRms");
    if(hRms) delete hRms;
    
    hRms = new TH1F("hRms","",x1Nbin,x1Min,x1Max);
    
    if(CYL) axisPos = 0.0;
    
    for(Int_t j=0;j<hRms->GetNbinsX();j++) {
      Double_t rms = 0;
      Double_t total = 0;
      for(Int_t k=1;k<=x2Nbin;k++) {
	Double_t value  = hDen2D->GetBinContent(j,k);
	Double_t radius = hDen2D->GetYaxis()->GetBinCenter(k) - axisPos;
	if(CYL) {
	  rms += radius*radius*radius*value;
	  total += radius*value;
	} else {

//.........这里部分代码省略.........

  TGraph *gChi2EBP = new TGraph();
  gChi2EBP->SetMarkerStyle(20);
  gChi2EBP->SetMarkerSize(0.7);
  gChi2EBP->SetMarkerColor(3);

  TGraph *gOutOfTimeChi2EBP = new TGraph();
  gOutOfTimeChi2EBP->SetMarkerStyle(20);
  gOutOfTimeChi2EBP->SetMarkerSize(0.7);
  gOutOfTimeChi2EBP->SetMarkerColor(4);
 
  TGraph *gTimeEBM = new TGraph();
  gTimeEBM->SetMarkerStyle(24);
  gTimeEBM->SetMarkerSize(0.7);
  gTimeEBM->SetMarkerColor(2);

  TGraph *gChi2EBM = new TGraph();
  gChi2EBM->SetMarkerStyle(24);
  gChi2EBM->SetMarkerSize(0.7);
  gChi2EBM->SetMarkerColor(3);

  TGraph *gOutOfTimeChi2EBM = new TGraph();
  gOutOfTimeChi2EBM->SetMarkerStyle(24);
  gOutOfTimeChi2EBM->SetMarkerSize(0.7);
  gOutOfTimeChi2EBM->SetMarkerColor(4);


  float effS, effN;

  for (int i = 0; i < n ; i++){

    effN = (float)nNormalTime[i]/(float)nNormalTot;
    effS = (float)nSpikeTime[i]/(float)nSpikeTot;
    gTime->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2[i]/(float)nNormalTot;
    effS = (float)nSpikeChi2[i]/(float)nSpikeTot;
    gChi2->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2OutOfTime[i]/(float)nNormalTot;
    effS = (float)nSpikeChi2OutOfTime[i]/(float)nSpikeTot;
    gOutOfTimeChi2->SetPoint(i,effS,effN);

    //EB-
    effN = (float)nNormalTimeEB[i][0]/(float)nNormalTotEB[0];
    effS = (float)nSpikeTimeEB[i][0]/(float)nSpikeTotEB[0];
    gTimeEBM->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2EB[i][0]/(float)nNormalTotEB[0];
    effS = (float)nSpikeChi2EB[i][0]/(float)nSpikeTotEB[0];
    gChi2EBM->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2OutOfTimeEB[i][0]/(float)nNormalTotEB[0];
    effS = (float)nSpikeChi2OutOfTimeEB[i][0]/(float)nSpikeTotEB[0];
    gOutOfTimeChi2EBM->SetPoint(i,effS,effN);

    //EB+
    effN = (float)nNormalTimeEB[i][1]/(float)nNormalTotEB[1];
    effS = (float)nSpikeTimeEB[i][1]/(float)nSpikeTotEB[1];
    gTimeEBP->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2EB[i][1]/(float)nNormalTotEB[1];
    effS = (float)nSpikeChi2EB[i][1]/(float)nSpikeTotEB[1];
    gChi2EBP->SetPoint(i,effS,effN);

    effN = (float)nNormalChi2OutOfTimeEB[i][1]/(float)nNormalTotEB[1];
    effS = (float)nSpikeChi2OutOfTimeEB[i][1]/(float)nSpikeTotEB[1];
    gOutOfTimeChi2EBP->SetPoint(i,effS,effN);

  }




  TFile saving (outputRootName.c_str (),"recreate") ;
  saving.cd () ;  
  
  // saving distributions
  hRun->Write();

  gTime->Write("g_Time");
  gChi2->Write("g_Chi2");
  gOutOfTimeChi2->Write("g_OutOfTimeChi2");

  gTimeEBP->Write("g_Time_EBP");
  gChi2EBP->Write("g_Chi2_EBP");
  gOutOfTimeChi2EBP->Write("g_OutOfTimeChi2_EBP");

  gTimeEBM->Write("g_Time_EBM");
  gChi2EBM->Write("g_Chi2_EBM");
  gOutOfTimeChi2EBM->Write("g_OutOfTimeChi2_EBM");


  hTime[0]->Write();
  hTime[1]->Write();

  saving.Close () ;
 
  return 0 ;
}

int main(int argc, char **argv) {
    char *ttbarFileName = NULL;
    char *wp3jetsFileName = NULL;
    char *wp4jetsFileName = NULL;
    char *smwwFileName = NULL;
    char *opsFileName[10];
    int nOpsFile = 0;
    char *outputFileName = NULL;
    int c;
    while(1) {
        int option_index = 0;
        static struct option long_options[] = {
            {"ttbar", required_argument, 0, 'a'},
            {"wp3jets", required_argument, 0, 'b'},
            {"wp4jets", required_argument, 0, 'c'},
            {"smww", required_argument, 0, 'd'},
            {"opsww", required_argument, 0, 'e'},
            {"output", required_argument, 0, 'f'}
        };
        c = getopt_long(argc, argv, "abcdef:", long_options, &option_index);
        if (c==-1)
            break;
        switch(c) {
            case 'a':
                ttbarFileName = optarg;
                break;
            case 'b':
                wp3jetsFileName = optarg; 
                break;
            case 'c':
                wp4jetsFileName = optarg;
                break;
            case 'd':
                smwwFileName = optarg;
                break;
            case 'e':
                opsFileName[nOpsFile] = optarg;
                nOpsFile++;
                break;
            case 'f':
                outputFileName = optarg;
                break;
        }
    }
    /*
    double cww[4] = {4*pow(10, -6), 4.6*pow(10, -6), 5.33*pow(10, -6), 6.25*pow(10, -6)};
    double lambda[4];
    for (int i = 0; i < 4; i++) {
        lambda[i] = 1.0/sqrt(cww[i]);
    }
    */
    double lambda[4] = {500, 466, 433, 400};

    TFile *outputFile = new TFile(outputFileName, "RECREATE");
    outputFile->Close();
    double significance[10];
    for (int i = 0; i < nOpsFile; i++) {
        std::cerr << "using signal file " << opsFileName[i] << "\n";
        significance[i] = RunHypoTest(smwwFileName, ttbarFileName, wp3jetsFileName, wp4jetsFileName, opsFileName[i], outputFileName, lambda[i]);
        std::cerr << "Significance = " << significance[i] << "\n";
    }
    TGraph *graph = new TGraph(nOpsFile, lambda, significance);
    graph->SetFillColor(kRed);
    graph->SetTitle("");
    graph->GetYaxis()->SetTitle("Significance (#sigma)");
    graph->GetYaxis()->CenterTitle();
    graph->GetXaxis()->SetTitle("#Lambda (GeV)");
    graph->GetXaxis()->CenterTitle();
    TFile *outputFile2 = new TFile(outputFileName, "UPDATE");
    TCanvas *canvas = new TCanvas();
    graph->Draw("A*");
    canvas->Write();
    graph->Write();
    std::cerr << "---------------------------------------------------\n";
    for (int i = 0; i < nOpsFile; i++) {
        std::cerr << "Significance = " << significance[i] << "\n";
    }
    std::cerr << "n files: " << nOpsFile << "\n";
    outputFile2->Close();
}

//.........这里部分代码省略.........
   }

   // MIGRAD: error optimiser
   if(ierflg==0 && 0){
     arglist[0] = 500000;
     gMinuit->mnexcm("MIGRAD", arglist ,1,ierflg);
     cout << " MIGRAD flag is " << ierflg << endl;
     if(ierflg!=0){
       cout << "+++++++++++++++++++++" << endl;
       cout << "+ ups.. check again +" << endl;
       cout << "+++++++++++++++++++++" << endl;
     }
   }

   // Minos: error optimiser
   if(ierflg==0 && 0){
     arglist[0] = 500000;
     gMinuit->mnexcm("MINOS", arglist ,1,ierflg);
     cout << " MINOS flag is " << ierflg << endl;
     if(ierflg!=0){
       cout << "+++++++++++++++++++++" << endl;
       cout << "+ ups.. check again +" << endl;
       cout << "+++++++++++++++++++++" << endl;
     }
   }

// Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(3,amin);
 
   cout << "prob = " << TMath::Prob(amin,gEntriesUsed-gMinuit->GetNumFreePars())*100 << " %" << endl;
   cout << "ndf = " << gEntriesUsed-gMinuit->GetNumFreePars() << endl;

   
   cout << "output" << endl;
//   gMinuit->mnexcm("SHO COV", arglist ,0,ierflg);
//   gMinuit->mnexcm("SHO COR", arglist ,0,ierflg);
//   gMinuit->mnexcm("SHO EIG", arglist ,0,ierflg);
   
   //get fit parameters
   double par[15], epar;
   for(int i=0; i<NoPar; i++){
     gMinuit->GetParameter(i,par[i],epar);
   }

   
   TFile *output = new TFile("FitParFieldStrength_out.root","RECREATE");

   int tcount = 0;
   float tmpBAngle, tmpZenith;
   float r0x[gMaxEvents], r0y[gMaxEvents];
   float bx[gMaxEvents], by[gMaxEvents];
   
   for(int i=0; i<gEntries; i++){
     if(gFieldStrengthAnt[i]>0 && gFieldStrengthNS[i]>0 && gFieldStrengthEW[i]>0){
       tmpBAngle = gBAngle[i] * TMath::Pi() / 180.;
       tmpZenith = gZenith[i] * TMath::Pi() / 180.;
       
       r0x[tcount] = -gDistanceShowerAxis[i];
       r0y[tcount] = (log( (gFieldStrengthAnt[i]/gEffBand) / ( par[0]*cos(tmpZenith+par[3])*(par[5]+sin(tmpBAngle+par[4]))*powf(10,par[2]*(gEg[i]-8))))  );
       //if(r0y[tcount] < -4 && r0x[tcount] > 50) cout << "E = " << gEg[i] << " --  i = " << i << endl;
       
       bx[tcount] = tmpBAngle;
       by[tcount] = gFieldStrengthAnt[i]/gEffBand / (par[0]*cos(tmpZenith+par[3])*exp(-gDistanceShowerAxis[i]/par[1])*powf(10,par[2]*(gEg[i]-8))) - par[5];
//       if(by[tcount]>1.5) cout << " E = " << gEg[i] << " -- i = " << i << endl;

       tcount++;
     }
   }
   TCanvas *can1 = new TCanvas("can1");
   can1->Divide(2,2);

   TGraph *gr0 = new TGraph(tcount,r0x, r0y);
   gr0->SetName("gr0");
   gr0->SetMarkerStyle(20);
   gr0->SetTitle("lateral distribution");

   TGraph *gb = new TGraph(tcount,bx, by);
   gb->SetName("gb");
   gb->SetMarkerStyle(20);
   gb->SetTitle("geomagnetic angle");
 
   can1->cd(1);
   gr0->Draw("ap");
   
   can1->cd(2);
   gb->Draw("ap");
 
 
   gr0->Write();
   gb->Write();
   can1->Write();
   output->Close();
   
   
   

}

//.........这里部分代码省略.........
    if (cat==4 || cat==5 || cat==6 || cat==7 || cat==8) {
      modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      */
    }
    map<string,RooAbsPdf*> bkgPdfs = modelBuilder.getBkgPdfs();
    modelBuilder.setSignalPdfFromMC(sigMC);
    modelBuilder.makeSBPdfs();
    map<string,RooAbsPdf*> sbPdfs = modelBuilder.getSBPdfs();

    modelBuilder.fitToData(dataBinned,true,true);
    modelBuilder.fitToData(dataBinned,false,true);

    modelBuilder.throwToy(Form("cat%d_toy0",cat),dataBinned->sumEntries(),true,true);

    // Profile this category using ProfileMultiplePdfs
    ProfileMultiplePdfs profiler;
    for (map<string,RooAbsPdf*>::iterator pdf=sbPdfs.begin(); pdf!=sbPdfs.end(); pdf++) {
      string bkgOnlyName = pdf->first.substr(pdf->first.find("sb_")+3,string::npos);
      if (bkgPdfs.find(bkgOnlyName)==bkgPdfs.end()){
        cerr << "ERROR -- couldn't find bkg only pdf " << bkgOnlyName << " for SB pdf " << pdf->first << endl;
        pdf->second->fitTo(*dataBinned);
        exit(1);
      }
      int nParams = bkgPdfs[bkgOnlyName]->getVariables()->getSize()-1;
      profiler.addPdf(pdf->second,2*nParams);
      //profiler.addPdf(pdf->second);
      cout << pdf->second->GetName() << " nParams=" << pdf->second->getVariables()->getSize() << " nBkgParams=" << nParams << endl;
    }
    profiler.printPdfs();
    //cout << "Continue?" << endl;
    //string bus; cin >> bus;
    profiler.plotNominalFits(dataBinned,mass,80,Form("cat%d",cat));
    pair<double,map<string,TGraph*> > minNlls = profiler.profileLikelihood(dataBinned,mass,mu,mu_low,mu_high,mu_step);
    pair<double,map<string,TGraph*> > correctedNlls = profiler.computeEnvelope(minNlls,Form("cat%d",cat),2.);
    minNlltrack.push_back(make_pair(correctedNlls.first,correctedNlls.second["envelope"]));
    //minNlls.second.insert(pair<string,TGraph*>("envelope",envelopeNll.second));
    //map<string,TGraph*> minNLLs = profiler.profileLikelihoodEnvelope(dataBinned,mu,mu_low,mu_high,mu_step);
    profiler.plot(correctedNlls.second,Form("cat%d_nlls",cat));
    //profiler.print(minNLLs,mu_low,mu_high,mu_step);
    /*
    if (minNLLs.find("envelope")==minNLLs.end()){
      cerr << "ERROR -- envelope TGraph not found in minNLLs" << endl;
      exit(1);
    }
    */
    //minNlltrack.push_back(make_pair(profiler.getGlobalMinNLL(),minNLLs["envelope"]));
  }
  //exit(1);
  TGraph *comb = new TGraph();
  for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
    if (it->second->GetN()!=minNlltrack.begin()->second->GetN()){
      cerr << "ERROR -- unequal number of points for TGraphs " << it->second->GetName() << " and " << minNlltrack.begin()->second->GetName() << endl;
      exit(1);
    }
  }
  for (int p=0; p<minNlltrack.begin()->second->GetN(); p++){
    double x,y,sumy=0;
    for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
      it->second->GetPoint(p,x,y);
      sumy += (y+it->first);
    }
    comb->SetPoint(p,x,sumy);
  }
  pair<double,double> globalMin = getGraphMin(comb);
  for (int p=0; p<comb->GetN(); p++){
    double x,y;
    comb->GetPoint(p,x,y);
    comb->SetPoint(p,x,y-globalMin.second);
  }
  vector<double> fitVal = getValsFromLikelihood(comb);

  cout << "Best fit.." << endl;
  cout << "\t mu = " << Form("%4.3f",fitVal[0]) << " +/- (1sig) = " << fitVal[2]-fitVal[0] << " / " << fitVal[0]-fitVal[1] << endl;
  cout << "\t      " << "    " << " +/- (2sig) = " << fitVal[4]-fitVal[0] << " / " << fitVal[0]-fitVal[3] << endl;

  cout << comb->Eval(fitVal[0]) << " " << comb->Eval(fitVal[1]) << " " << comb->Eval(fitVal[2]) << " " << comb->Eval(fitVal[3]) << " " << comb->Eval(fitVal[4]) << endl;

  double quadInterpVal = ProfileMultiplePdfs::quadInterpMinimum(comb);
  cout << "quadInterp: mu = " << quadInterpVal << endl;
  cout << "\t " << comb->Eval(quadInterpVal) << " " << comb->Eval(quadInterpVal-0.005) << " " << comb->Eval(quadInterpVal-0.01) << " " << comb->Eval(quadInterpVal+0.005) << " " << comb->Eval(quadInterpVal+0.01) << endl;
  
  comb->SetLineWidth(2);
  TCanvas *canv = new TCanvas();
  comb->Draw("ALP");
  canv->Print("plots/comb.pdf");
  TFile *tempOut = new TFile("tempOut.root","RECREATE");
  tempOut->cd();
  comb->SetName("comb");
  comb->Write();
  tempOut->Close();
  return 0;
}