C++ TH1I::GetBinContent方法代码示例

//______________________________________________________________________
void EventInfo(Int_t event, Int_t px, Int_t , TObject *selected)
{
   //draw the tooltip showing the backtrace for the bin at px
   if (!gTip) return;
   gTip->Hide();
   if (event == kMouseLeave)
      return;
   Double_t xpx  = gPad->AbsPixeltoX(px);
   Int_t bin = hleaks->GetXaxis()->FindBin(xpx);
   if (bin <=0 || bin > hleaks->GetXaxis()->GetNbins()) return;
   Int_t nbytes = (Int_t)hleaks->GetBinContent(bin);
   Int_t entry  = (Int_t)hentry->GetBinContent(bin);
   Int_t btid   = (Int_t)V4[entry];
   Double_t  time = 0.0001*V3[entry];
   TH1I *hbtids = (TH1I*)T->GetUserInfo()->FindObject("btids");
   if (!hbtids) return;
   if (!btidlist) btidlist = (TObjArray*)T->GetUserInfo()->FindObject("FAddrsList");
   if (!btidlist) btidlist = (TObjArray*)f->Get("FAddrsList"); //old memstat files
   if (!btidlist) return;
   Int_t nbt = (Int_t)hbtids->GetBinContent(btid-1);
   TString ttip;
   for (Int_t i=0;i<nbt;i++) {
      Int_t j = (Int_t)hbtids->GetBinContent(btid+i);
      TNamed *nm = (TNamed*)btidlist->At(j);
      if (nm==0) break;
      char *title = (char*)nm->GetTitle();
      Int_t nch = strlen(title);
      if (nch < 20) continue;
      if (nch > 100) title[100] =0;
      const char *bar = strstr(title,"| ");
      if (!bar) continue;
      if (strstr(bar,"operator new")) continue;
      if (strstr(bar,"libMemStat")) continue;
      if (strstr(bar,"G__Exception")) continue;
      ttip += TString::Format("%2d %s\n",i,bar+1);
   }
   
   if (selected) {
      TString form1 = TString::Format("  Leak number=%d, leaking %d bytes at entry=%d    time=%gseconds\n\n",bin,nbytes,entry,time);
      gTip->SetText(TString::Format("%s%s",form1.Data(),ttip.Data() ));
      gTip->SetPosition(px+15, 100);
      gTip->Reset();
   }
}

int main()
{

  float const conS =.3/2.35;
  ofstream fout("cal/backN.cal");
  ofstream fwhm("cal/fwhmback.dat");
  TFile f("sort.root");
  TCanvas* canvas[14];
  int Ntele = 14;
  int Nstrip = 32;
  ostringstream outstring;
  TH2I frame("frame","",10,4.5,9.5,10,0,130);
  frame.SetStats(kFALSE);

  double xx[14*32];
  double yy[14*32];

  TF1 *func = new TF1("fit",ThPeaks,3,9,4);
  double para[5];
  ifstream file("cal/back.cal");
  float intercept, slope;
  int i1,i2;
  string name;

  TH1F con("con","",500,0,10);
  for (int itele=0;itele<Ntele;itele++)
    {
      outstring.str("");
      outstring << "B"<<itele;
      name = outstring.str();
      canvas[itele] = new TCanvas(name.c_str());
      canvas[itele]->Divide(6,6);
      for (int istrip =0;istrip<Nstrip;istrip++)
        {
      
          canvas[itele]->cd(istrip+1);
          file >> i1 >> i2  >> slope >> intercept;


          outstring.str("");
          outstring << "back/cal/EBC"<<itele<<"_"<<istrip;
          string name = outstring.str();
          cout <<  name << endl;
          TH1I * hist = (TH1I*) f.Get(name.c_str());
	  frame.Draw();


          hist->SetStats(kFALSE);
          hist->GetXaxis()->SetRangeUser(4.5,9.5);
          con.GetXaxis()->SetRangeUser(4.5,9.5);
	  for (int i=1;i<=500;i++)
	    for (int j=1;j<500;j++)
	    {
              float deltax = hist->GetBinCenter(i)-con.GetBinCenter(j);
	      if (fabs(deltax) > 10.*conS)continue;
              float fact = gauss(deltax,0.,conS);
	      float y = fact*hist->GetBinContent(i)*hist->GetBinWidth(i);
              con.SetBinContent(j,y+con.GetBinContent(j));
	    }



	  for (int i=1;i<=500;i++) 
	    {
	      hist->SetBinContent(i,con.GetBinContent(i));
	      con.SetBinContent(i,0.);
	    }

          hist->Draw("same");



          func->SetParameter(0,0);
          func->SetParameter(1,1.);
	  func->FixParameter(2,conS);
          //func->SetParameter(2,0.1);
          func->SetParameter(3,8.);
          func->SetLineColor(2);
          //func->Draw("same");



	  hist->Fit(func);
          func->GetParameters(para);
          cout << "chisq=" << func->GetChisquare() << endl;
           if (fabs(para[1]-1.) < .2) 
	     { 

              slope *= para[1];
              intercept = intercept*para[1] + para[0];
	     }
            fout << itele << " " << istrip << " " 
                 << slope << " " << intercept << endl;
            fwhm << itele << " " << istrip << " " 
                  << para[2]*2.35 << endl;
            int ii = itele*32+istrip;
            xx[ii] = (float)ii;
            yy[ii] = para[2]*2.35;
            cout << para[0] << " " << para[1] << " " << para[2] << endl;

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

void makehisto(TString varname, TString process, TString ofile, TString file)
{
	TFile *fo = new TFile("histos/"+ofile+".root","RECREATE");
	TFile *f = new TFile("trees_8TeV/Nominal/"+file+".root");

	TTree *tree = (TTree*)f->Get("trees/Events");

	cout << "Filling histogram " << varname << "__" << process << " with " << file << endl;

	fo->cd();

	Float_t lumi = 19739;

	// weights
	Float_t weight = 1.0;
	Float_t pu_weight = 1.0;
	Float_t btag_weight = 1.0;
	Float_t muonID_weight = 1.0;
	Float_t muonIso_weight = 1.0;
	Float_t muonTrigger_weight = 1.0;
	tree->SetBranchAddress("pu_weight", &pu_weight);
	tree->SetBranchAddress("b_weight_nominal", &btag_weight);
	tree->SetBranchAddress("muon_IDWeight", &muonID_weight);
	tree->SetBranchAddress("muon_IsoWeight", &muonIso_weight);
	tree->SetBranchAddress("muon_TriggerWeight", &muonTrigger_weight);
	// cut variables
	Float_t muonIso = 0;
	Float_t jetrms = 0;
	Float_t mtw = 0;
	tree->SetBranchAddress("mu_iso", &muonIso);
	tree->SetBranchAddress("rms_lj", &jetrms);
	tree->SetBranchAddress("mt_mu", &mtw);
	//
	Int_t wjets_class = -1;
	tree->SetBranchAddress("wjets_classification", &wjets_class);

	Float_t var = 0;
	tree->SetBranchAddress(varname, &var);
	
	TH1::SetDefaultSumw2(true);
	TH1I *hcount = (TH1I*)f->Get("trees/count_hist");
	Int_t count = hcount->GetBinContent(1); 
	cout << "sample size: " << count << endl;

	Float_t xsec = get_xsec(file);
	cout << "xsec: " << xsec << endl;
	
	// FIXME adapt range for each variable histogram
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,30,-4.5,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,15,-4.5,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,20,-4.5,4.5);
	// FIXME
	TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,15,0,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,20,0,4.5); // too much
	histo->Sumw2();
	
	Long_t nentries = tree->GetEntries();
	for(Long_t i = 0; i < nentries; i++) {
		tree->GetEntry(i);

		// wjets
		if(process == "wjets_light") {
			if(wjets_class == 0 || wjets_class == 1) continue;
		} else if (process == "wjets_heavy") {
			if(wjets_class != 0 && wjets_class != 1) continue;
		}
		
		if(jetrms >= 0.025) continue;
		if(mtw <= 50) continue;
		if(process == "qcd") {
			if (muonIso < 0.3 || muonIso > 0.5) continue;
		}
		
		if(process != "DATA"  && process != "qcd") {
			weight = pu_weight*btag_weight*muonID_weight*muonIso_weight*muonTrigger_weight;
			weight *= lumi;
			weight *= xsec;
			weight /= (Float_t)count;
		}
		
		// FIXME
		//histo->Fill(var,weight);
		histo->Fill(TMath::Abs(var),weight);
	}
	if(process == "qcd") histo->Scale(xsec);


	cout << "integral: " << histo->Integral() << endl;
	// write histos
	histo->Write();
	fo->Close();
}

void Zlumi::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L Zlumi.C
//      Root > Zlumi t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
 
 gROOT->ForceStyle();
        tdrStyle();

  if (fChain == 0) return;

int minRun=0;
int maxRun=0;
int maxLS=0;
bool forminRun=true;


//TH2I * LumiSRun = new TH2I("LumiSRun", "LS vs Run", 3000, 0., 3000., 20000, 160000., 180000.);
//TH1I * test2 = new TH1I("test2","test2", 3000,0, 3000);


Long64_t nentries = fChain->GetEntriesFast();

Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
	Long64_t ientry = LoadTree(jentry);
	if (ientry < 0) break;
	nb = fChain->GetEntry(jentry);   nbytes += nb;


	if(forminRun && (Run!=0)){minRun=Run; forminRun=false;}
	if((Run!=0) && (Run>maxRun)){maxRun=Run;}
	if((LS!=0) && (LS>maxLS)){maxLS=LS;}
	//printf("run %i ls %i \n",Run,LS);
}
cout << nentries << " nentries \n";

TH2I *LumiSRun = new TH2I("LumiSRun", "LS vs Run", maxLS, 0, maxLS, maxRun-minRun+2, minRun-1, maxRun+1);
TH1I *Runs = new TH1I("Runs","Run", maxRun-minRun+2, minRun-1, maxRun+1);
 Runs->Sumw2();
nbytes = 0; nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
	Long64_t ientry = LoadTree(jentry);
	if (ientry < 0) break;
	nb = fChain->GetEntry(jentry);   nbytes += nb;

	LumiSRun->Fill(LS,Run);
	//printf("run %i ls %i \n",Run,LS);
	Runs->Fill(Run);
	//test2->Fill(LS);
}
	printf("minRun %i maxRun %i \n",minRun,maxRun);
LumiSRun->Draw();
 for (int h=0;h<Runs->GetNbinsX();h++){
   Runs->SetBinError(h+1,sqrt(Runs->GetBinContent(h+1)) );
 }
Runs->Draw();
//test->Draw();
//test2->Draw();

TH1F *FileRuns = new TH1F("FileRuns","Run from Lumicalc", maxRun-minRun+2, minRun-1, maxRun+1);
TH1D *XsecDistro = new TH1D("XsecDistro","X sec distribution", 60, 0., 0.6);
//-------------
   const Int_t mpt = maxRun-minRun;
   int fileRun[mpt];
   double Lumi[mpt];

   int npt = 0;
   // read data file
   ifstream file;
   //file.open("./2011-run-lumi.txt");
   file.open("./LumiAeB-dav.txt");
   while (1) {

	   file >> fileRun[npt] >> Lumi[npt];
	   if ( ! file.good() ) break;
	   cout << "x = " << fileRun[npt] << " y = " << Lumi[npt] << endl;
	   
	   FileRuns->SetBinContent((fileRun[npt]-minRun+2),0.2);
	   npt++;
//.........这里部分代码省略.........