C++ TH1D::Fill方法代码示例

void GetRandomTest(){
	double k0=1.39;
	double k1 = 0.425;
	double theta0 = 3.41;
	double theta1 = 1.30;
	int iNpart=2;
	double k_=k0+k1*(iNpart-2);
        double theta_=theta0+theta1*TMath::Log(iNpart-1);
	TF1 *f = new TF1("f","TMath::GammaDist(x,[0],0,[1])",0,200);
	f->SetParameters(k1,theta_);
	cout<<"Value at 0 = "<<f->Eval(0)<<endl;
	cout<<"Value at 1e-11 = "<<f->Eval(1e-11)<<endl;
	cout<<"Integral 1= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	f->SetRange(1e-12,200);
	cout<<"Integral 2= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	cout<<"fXmin = "<<f->GetXmin()<<"\tfXmax = "<<f->GetXmax()<<"\tfNpx = "<<f->GetNpx()<<endl;
	f->SetNpx(1e5);
	TCanvas *c1 = new TCanvas();
	c1->SetLogy();
	cout<<"f mean = "<<f->Mean(0,200)<<endl;
	cout<<"math mean = "<<f->GetParameter(0)*f->GetParameter(1)<<endl;
	TH1D* h = new TH1D("h","h",1000,0,200);
	for(int i=0;i<1e6;i++){
		double para = f->GetRandom();
		h->Fill(para);
	}
	h->Scale(1.0/h->Integral()*1000/200);
	h->GetYaxis()->SetRangeUser(1e-10,1);
	h->SetMarkerStyle(24);
	h->SetMarkerColor(4);
	h->SetMarkerSize(1.1);
	TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
        leg->SetFillColor(0);
        leg->SetFillStyle(0);
        leg->SetBorderSize(0);
        leg->SetTextFont(42);
        leg->SetTextSize(0.03);
        leg->AddEntry(f,"function","lp");
        leg->AddEntry(h,"filled histogram","lp");
	h->Draw("P");
	f->Draw("same");
        leg->Draw("same");
	cout<<"h mean = "<<h->GetMean(1)<<endl;
	c1->Print("TestGetRandom.png");
	}

//Plot channels with multiple hits
void PlotMultCh(int run) {
	etf.SetPath(mypath.Data());
	etf.Open(run); 
	int maxch = MAXCH*MAXRIO;
	TH1D* h = MakeHist("hnch",maxch);
	for (int i=0;i<etf.GetNumEvents();i++) {
		etf.GetEvent(i);
		int nch = 0;
		for (int ch=0;ch<maxch;ch++) {
			if (etf.myEvent.E[ch] > 0) nch++;
		}
		if (nch > 1)
			for (int ch=0;ch<maxch;ch++) {
				if (etf.myEvent.E[ch] > 0)
					h->Fill(ch);
			}
	}
	h->Draw();
}

void SearchBinEventCount::Fill(double HT, double MHT, int NJets, int BTags, double Weight)
{
  double bin = GetBinNumber(HT,MHT,NJets,BTags);
	
  if(bin<bins_.size()+2) 
    {
      fullTH1D_->Fill(bin-0.01, Weight);
      unsigned int splitHist=0;
      // 	std::cout<<"bin before split: "<<bin<<std::endl;
      for(int ii=0;bin>splitAfter_;ii++)
	{
	  splitHist++;
	  bin = bin-splitAfter_;
	}
      // 		if(splitHist==3)std::cout<<"BinForSplit: "<<bin<<" with splitHistNumber "<<splitHist<<" and TH1DSearchBinsSplit_.size(): "<<TH1DSearchBinsSplit_.size()<<std::endl;
		
      splitTH1D_[splitHist]->Fill(bin-0.1, Weight);
    }
}

//################################################################################################################################
 // Get the 68% confindence interval of mu1/mu2
void getPoissonIntervalls(double mu1, double mu2){

  double result = mu1/mu2;
  TH1D *expHist = new TH1D("mu1/mu2","mu1/mu2",1000,0,1);
  TRandom3 rand1(0);
  TRandom3 rand2(0);

  cout<<"mu1 = "<<mu1<<endl;
  cout<<"mu2 = "<<mu2<<endl;
  cout<<"ratio = "<<result<<endl;

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

    expHist->Fill(1.*rand1.Poisson(mu1)/rand2.Poisson(mu2));
  }

  // Get now the 68% upper and lower interval
  double errUp  = 0;
  double errLow = 0;
  for(int i=1; i<=expHist->GetNbinsX();i++){

    double upperIntegral = expHist->Integral(i,expHist->GetNbinsX());
    double lowerIntegral = expHist->Integral(1,i);

    if(abs(upperIntegral/expHist->Integral()-0.32)<0.02){
      cout<<"upper bound = "<<expHist->GetBinCenter(i)<<endl;
      errUp = expHist->GetBinCenter(i)-result;
    }
    if(abs(lowerIntegral/expHist->Integral()-0.32)<0.02){
      cout<<"lower bound = "<<expHist->GetBinCenter(i)<<endl;
      errLow = result - expHist->GetBinCenter(i);
    }
  }


  cout<<"error Up = "<<errUp<<endl;
  cout<<"error Low = "<<errLow<<endl;

  TCanvas *c = new TCanvas("cExp","cExp",0,0,500,500);
  c->cd();
  expHist->Draw();
}    

void L3_pValues(double l0, double l1, double Metl, double ll, int Jets)
{
	InitExterns();
	RAND.SetSeed(5);

	Data d(l0,l1,Metl,ll,Jets);
	TH1D* h_b = Likelihood::GetBGEstimation(d);
	TH1D* h_pValues = new TH1D("pvalues","pvalues",22,0,1.1);
	int numMC = 200;

	for (int i=0; i<numMC; i++)
	{
		if (i % (numMC/10) == 0){ cout << i*(100./numMC) << "%-" << flush;}
		Data dRand;
		dRand = Toys::ToyData(d,h_b);
		TH1D* h_b_rand = Likelihood::GetBGEstimation(dRand);
		TH1D* pdf = new TH1D("pdf","pdf",100000,0,100);
		int numMC2 = 500;
		for (int j=0; j<numMC2; j++)
		{
			Data dToy;
			dToy = Toys::ToyData(dRand,h_b_rand);
			double qZero = Likelihood::qZero(dToy);
			dToy.free();
			pdf->Fill(qZero,1./numMC2);
		}

		double obs = Likelihood::qZero(dRand);
		dRand.free();
//		pdf->Draw();
		double pV = utilities::pValue(pdf,obs);
		delete pdf;
		delete h_b_rand;
		h_pValues->Fill(pV,1./numMC);
	}

	h_pValues->Draw();
	TFile* file = new TFile("f5.root","RECREATE");
	h_pValues->Write("L3_pValues");
	file->Close();
}

	double FitConfidence::choleskyUncertainty( double xx, double * fCov, TF1 * f, int nSamples ){
		int nP = f->GetNpar();
		INFO( FitConfidence::classname(), "Num Params : " << nP  );
		
		double *fCovSqrt = new double[ nP * nP ];
		calcCholesky( nP, fCov, fCovSqrt );

		double yerr = 0;

		TH1D *hDistributionAtX = new TH1D("hDistributionAtX","",200,f->Eval(xx) - .2,f->Eval(xx) + .2);
		
		for (int n = 0; n < nSamples; n++ ) {
			double val = randomSqrtCov(xx,f,nP,fCovSqrt);
			hDistributionAtX->Fill( val );
		}
		yerr = hDistributionAtX->GetRMS();
		hDistributionAtX->Delete();

		return yerr;

	}

void drawJacob(string fnm="../../output/thinRadiator/o_msc_L001um_03MeV_1e5ev.root"){
  TFile *fin=TFile::Open(fnm.c_str(),"READ");
  TTree *t=(TTree*)fin->Get("t");
  int material, pType, trackID, parentID;
  double pX, pY, pZ;
  
  t->SetBranchAddress("pType", &pType);
  t->SetBranchAddress("material", &material);
  t->SetBranchAddress("trackID", &trackID);
  t->SetBranchAddress("parentID", &parentID);
  t->SetBranchAddress("postMomX", &pX);
  t->SetBranchAddress("postMomY", &pY);
  t->SetBranchAddress("postMomZ", &pZ);

  int nev = t->GetEntries();
  TH1D *scatAng = new TH1D("scatAng","Rate/sin(theta); scattering angle [deg]", 400,0,100);
  double pi = acos(-1);
  cout<<"pi "<<pi<<endl;
  for(int i=0;i<nev;i++){
    t->GetEntry(i);

    if(pType!=11) continue;
    if(material != 1) continue;
    if(trackID != 1 || parentID!=0) continue;

    double r = sqrt( pX*pX + pY*pY + pZ*pZ);
    if(r==0) continue;

    double thetaRad = acos(pZ/r);
    double theta = thetaRad * 180/pi;
    if(theta<0.01) continue;
    
    scatAng->Fill(theta, 1./sin(thetaRad));
    
  }
  gStyle->SetOptStat("eMRou");
  scatAng->DrawCopy();
  fin->Close();
}

int main(int argc, char *argv[])
{
  int binCount = 100;
  if (argc > 2) {
    std::cout << "too many argument. exit." << std::endl;
    std::exit(1);
  } else if (argc == 2) {
    binCount = std::atoi(argv[1]);
    --argc;
  }
  TApplication app("app", &argc, argv);
  gStyle->SetOptStat(0);

  TH1D *hist = new TH1D("hist", "ns vs count", binCount, minNanoSec, maxNanoSec);

  char filename[256];
  GetFileName(filename, "../data/muon_lifetime/muon_lifetime.dat");
  std::ifstream ifs(filename);
  for (int i = 0; !ifs.eof(); ++i) {
    Double_t data;
    ifs >> data;
    if (data) {
      NanoSecWithError ns;
      ConvertTdcChannelToNanoSec(i%8, data, &ns);
      hist->Fill(ns.time);
    }
  }
  hist->Draw();

  TF1 *fit = new TF1("fit", "[0] * exp(- x / [1]) + [2]", 1000, 20000);
  fit->SetParameters(200, 2100, 40);
  hist->Fit(fit, "R+");
  std::cout << "chi^2/ndf: ";
  std::cout << fit->GetChisquare() / fit->GetNDF() << std::endl;

  app.Run();
  return 0;
}

TH1D* LoopGen(TTree* ntuple,double ptmin,double ptmax){

  int Run,size;
  bool cut_yvsRun,cut_pt;
  Float_t gen[NUM_BX];
  Float_t pdgId[NUM_BX];
  Float_t isSignal[NUM_BX];
  Float_t y[NUM_BX];
  Float_t pt[NUM_BX];
  

  ntuple->SetBranchAddress("size",&size);
  ntuple->SetBranchAddress("Run",&Run);
  ntuple->SetBranchAddress("gen",gen);
  ntuple->SetBranchAddress("pdgId",pdgId);
  ntuple->SetBranchAddress("isSignal",isSignal);
  ntuple->SetBranchAddress("y",y);
  ntuple->SetBranchAddress("pt",pt);

  TH1D *hPtMCGen = new TH1D("hPtMCGen","",nBins,ptBins);

  Int_t entries = (Int_t)ntuple->GetEntries();
  
  for (int i=0; i<entries; i++) {
    ntuple->GetEntry(i);
        
    for(int j=0;j<size;j++){
    
      cut_yvsRun=((Run<=1&&abs(y[j]+0.465)<1.93)||(Run>1&&abs(y[j]-0.465)<1.93))&&abs(pdgId[j])==511&&isSignal[j]!=0;
      cut_pt=((pt[j]>ptmin)&&(pt[j]<ptmax));
      if(cut_yvsRun&&cut_pt) hPtMCGen->Fill(pt[j]);
          	  
  	}//loop over candidates
  }// loop over events

  return hPtMCGen; 
}

void test_mc_rejection_method_lin()
{
  const int    N     = 500000;
  const int    nbins = 300;
  const double ymax  = 1.1;
  const double xmin  = 0;
  const double xmax  = 10;
  const double dx    = xmax-xmin;

  TRandom rg;

  TF1 *  func = new TF1  ("func","1/(x+1)",0,10);
  TH1D * hgen = new TH1D ("hgen","generated",nbins,xmin,xmax);

  for(int i=0; i<N; i++) {
     cout << "..................." << i << endl;
     bool selected=false;
     while(1) {
       double xg = xmin + rg.Uniform() * dx;
       double yg = ymax * rg.Uniform();
       double yc = func->Eval(xg);
       selected = (yg<yc);
       if(selected) {
           hgen->Fill(xg);
           break;
       }
     }
  }

  double IF = func->Integral(xmin,xmax);
  double IH = hgen->Integral("width");
  double sc = IF/IH;
  hgen->Scale(sc);

  hgen->Draw();
  func->Draw("same");
}

TH1D *getpuweights(TFile *file, TH1D *target) {
  
  TDirectory *dirmcpv = (TDirectory*)file->FindObjectAny("AnaFwkMod");
  TH1D *hnpu = (TH1D*)dirmcpv->Get("hNPU");
  TH1D *hpumc = (TH1D*)hnpu->Clone();
  
  hpumc->Sumw2();
  hpumc->Scale(1.0/hpumc->Integral(0,hpumc->GetNbinsX()+1));
  
  
  TH1D *htargettmp = new TH1D("htargettmp","", hpumc->GetNbinsX(), hpumc->GetXaxis()->GetXmin(), hpumc->GetXaxis()->GetXmax());
  htargettmp->Sumw2();
  for (int ibin = 0; ibin<=(htargettmp->GetNbinsX()+1); ++ibin) {
    htargettmp->Fill(htargettmp->GetBinCenter(ibin),target->GetBinContent(target->FindFixBin(htargettmp->GetBinCenter(ibin))));
  }
  htargettmp->Scale(1.0/htargettmp->Integral(0,htargettmp->GetNbinsX()+1));
  
  TH1D *puweights = new TH1D((*htargettmp)/(*hpumc));
    
  delete htargettmp;
  
  return puweights;
    
}

//*************************************************************************************************
//Main part of the macro
//*************************************************************************************************
void AnalyzeYYPhotonEvent(const string InputFilename) {

  
  TTree* YYTree = getTreeFromFile(InputFilename.c_str());
  assert(YYTree);    
  YYPhotonEvent photon(YYTree);
      
  TH1D *plot = new TH1D("hPlot", ";XAxis; Number of Events ", 6, -3,3);

  for (int n=0;n<YYTree->GetEntries();n++) { 
    photon.GetEntry(n);

    Double_t weight = photon.photon_branch_weight;    
    //cout << weight << " " << photon.photon_branch_NPhotons << endl;
    if (photon.photon_branch_NPhotons > 0 ) {
      plot->Fill(photon.photon_branch_Photon1Eta);
    }

  }

  TCanvas *c = new TCanvas("canvas" , "canvas" , 800, 600);
  plot->DrawCopy();
  
}

void DrawphiDis() {
    int ifile=0;
    TFile *f = TFile::Open(Form("/cms/store/user/qixu/flow/NewSTEG/pPbDataV205m300/vndata_50k_%d.root",ifile));
    TTree *t = (TTree*)f->Get("tree");
    Float_t eta[1000];
    Float_t phi[1000];
    TH1D* hphi = new TH1D("hphi","hphi",100,-2.4,2.4);
    int n;
    t->SetBranchAddress("etag",eta);
    t->SetBranchAddress("phig",phi);
    t->SetBranchAddress("n",&n);
    Int_t N = t->GetEntries();
    for(int ievt=0; ievt<N; ievt++) {
        //      if(ievt!=0) continue;
        // if(ievt%50000==0)        cout<<"Processing "<<ievt<<" events"<<endl;
        t->GetEntry(ievt);
        for(int imult = 0; imult<n; imult++) {
            hphi->Fill(eta[imult]);
        }
    }
//    cout<<hphi->Integral("width")/2/TMath::Pi()<<endl;
//    cout<<hphi->GetMean(2)<<endl;
    hphi->Draw();
}

//.........这里部分代码省略.........
                                       11,0.,11.);
  cutflow_preselection->GetXaxis()->SetBinLabel(1,"All Events");
  cutflow_preselection->GetXaxis()->SetBinLabel(2,"Sample based gen-selection");
  cutflow_preselection->GetXaxis()->SetBinLabel(3,"HBHEIsoNoiseFilter");
  cutflow_preselection->GetXaxis()->SetBinLabel(4,"eeBadScFilter");
  cutflow_preselection->GetXaxis()->SetBinLabel(5,"HBHENoiseFilter");
  cutflow_preselection->GetXaxis()->SetBinLabel(6,"GoodVtx");
  cutflow_preselection->GetXaxis()->SetBinLabel(7,"JetID Cleaning");

  int sampletype=-1;
  if(subSampleKey.find("TTbar_Inclusive")!=string::npos)sampletype=0; //TTbar_Inclusive
  else if(subSampleKey.find("TTbar_Tbar_SingleLep")!=string::npos || subSampleKey.find("TTbar_T_SingleLep")!=string::npos)sampletype=1;
  else if(subSampleKey.find("TTbar_DiLept")!=string::npos)sampletype=2;
  else if(subSampleKey.find("TTbar_HT_600_800")!=string::npos)sampletype=3;
  else if(subSampleKey.find("TTbar_HT_800_1200")!=string::npos)sampletype=4;
  else if(subSampleKey.find("TTbar_HT_1200_2500")!=string::npos)sampletype=5;
  else if(subSampleKey.find("TTbar_HT_2500_Inf")!=string::npos)sampletype=6;
  else if(subSampleKey.find("TTbar")!=string::npos){
    cout << " TT sample is not known. Please check the second input \n " ;
    return 2;
  }

  int TotNEve_ = utils2::TotNEve(subSampleKey);
  // Loop over the events (tree entries)
  int eventN=0;
  while( evt->loadNext() ){

  //if(eventN>1000)break;
    // Total weight
    //double totWeight = evt->weight()*1.;
    double totWeight = 10000.*evt->XS()/TotNEve_;
    //printf(" XS: %g NEvents: %d weight: %g \n ",evt->XS(),TotNEve_,totWeight);

    cutflow_preselection->Fill(0.,totWeight); // keep track of all events processed

    if(!evt->DataBool_()){

      if(sampletype==0){
        if(evt->gen_ht()>600||evt->GenElecPtVec_().size()>0||evt->GenMuPtVec_().size()>0||evt->GenTauPtVec_().size()>0)continue;
      }

      if(sampletype==1){
        if(evt->gen_ht()>600)continue;
      }

      if(sampletype==2){
        if(evt->gen_ht()>600)continue;
      }

    }

    cutflow_preselection->Fill(1.,totWeight);
    if(evt->HBHEIsoNoiseFilter_()==0)continue;
    cutflow_preselection->Fill(2.,totWeight);
    if(evt->eeBadScFilter_()==0)continue;
    cutflow_preselection->Fill(3.,totWeight);
    if(evt->HBHENoiseFilter_()==0)continue;
    cutflow_preselection->Fill(4.,totWeight);
    if(!(evt->NVtx_() >0))continue;
    cutflow_preselection->Fill(5.,totWeight);
    // Through out an event that contains HTjets with bad id
    if(evt->JetId()==0)continue;
    cutflow_preselection->Fill(6.,totWeight); // events passing JetID event cleaning


  vector<TLorentzVector> genTauJetLorVec;

//.........这里部分代码省略.........
    intree->SetBranchAddress("lumiSec",  &lumiSec);   // event lumi section
    intree->SetBranchAddress("evtNum",   &evtNum);    // event number
    intree->SetBranchAddress("npv",      &npv);       // number of primary vertices
    intree->SetBranchAddress("npu",      &npu);       // number of in-time PU events (MC)
    intree->SetBranchAddress("genVPt",   &genVPt);    // GEN W boson pT (signal MC)
    intree->SetBranchAddress("genVPhi",  &genVPhi);   // GEN W boson phi (signal MC)   
    intree->SetBranchAddress("scale1fb", &scale1fb);  // event weight per 1/fb (MC)
    intree->SetBranchAddress("met",      &met);       // MET
    intree->SetBranchAddress("metPhi",   &metPhi);    // phi(MET)
    intree->SetBranchAddress("sumEt",    &sumEt);     // Sum ET
    intree->SetBranchAddress("mt",       &mt);        // transverse mass
    intree->SetBranchAddress("u1",       &u1);        // parallel component of recoil
    intree->SetBranchAddress("u2",       &u2);        // perpendicular component of recoil
    intree->SetBranchAddress("q",        &q);	      // lepton charge
    intree->SetBranchAddress("lep",      &lep);       // lepton 4-vector
    intree->SetBranchAddress("pfChIso",  &pfChIso);
    intree->SetBranchAddress("pfGamIso", &pfGamIso);
    intree->SetBranchAddress("pfNeuIso", &pfNeuIso);
  
    //
    // loop over events
    //
    for(UInt_t ientry=0; ientry<intree->GetEntries(); ientry++) {
      intree->GetEntry(ientry);
      
      if(lep->Pt()        < PT_CUT)  continue;	
      if(fabs(lep->Eta()) > ETA_CUT) continue;
      
      if( (typev[ifile]==eAntiData || typev[ifile]==eAntiWmunu || typev[ifile]==eAntiEWK) &&
          (pfChIso+pfGamIso+pfNeuIso)>0.5*(lep->Pt()) ) 
	  continue;
      
      if(typev[ifile]==eData) {
        hDataMet->Fill(met);
	if(q>0) { hDataMetp->Fill(met); } 
	else    { hDataMetm->Fill(met); }
      
      } else if(typev[ifile]==eAntiData) {
        hAntiDataMet->Fill(met);
	if(q>0) { hAntiDataMetp->Fill(met); } 
	else    { hAntiDataMetm->Fill(met); }      
      
      } else {
        Double_t weight = 1;
        weight *= scale1fb*lumi;
	
	if(typev[ifile]==eWmunu) {
          Double_t corrMet=met, corrMetPhi=metPhi;
        
	  // apply recoil corrections to W MC
	  Double_t lepPt = lep->Pt();
	  //Double_t lepPt = gRandom->Gaus(lep->Pt(),0.5);  // (!) uncomment to apply scale/res corrections to MC
	  recoilCorr.Correct(corrMet,corrMetPhi,genVPt,genVPhi,lepPt,lep->Phi(),nsigma,q);
	
          Double_t nnlocorr=1;
          for(Int_t ibin=1; ibin<=hNNLOCorr->GetNbinsX(); ibin++) {
            if(genVPt >= hNNLOCorr->GetBinLowEdge(ibin) &&
               genVPt < (hNNLOCorr->GetBinLowEdge(ibin)+hNNLOCorr->GetBinWidth(ibin)))
              nnlocorr = hNNLOCorr->GetBinContent(ibin);
          }
	  //weight *= nnlocorr;  // (!) uncomment to apply NNLO corrections
	  
          hWmunuMet->Fill(corrMet,weight);
	  if(q>0) { hWmunuMetp->Fill(corrMet,weight); } 
	  else    { hWmunuMetm->Fill(corrMet,weight); }
        }