C++ TH2D类代码示例

double GetMultTrgMean(int itrg)
{
                   TH2D* hmulttrg = (TH2D*) gInputfile->Get(Form("%smult_trg_%d",subdir.Data(),itrg));
                   TH1D* hmulttrg_proj = (TH1D*)hmulttrg->ProjectionY(Form("mult_trg_proj_%d",itrg),-1,-1);
                   double  multtrg_mean = hmulttrg_proj->GetMean();
                   return multtrg_mean;
}

void plotTurnOn(TTree* inttree, TString triggerpass, TString variable, TString varname, TString varlatex, Int_t BIN_NUM, Double_t BIN_MIN, Double_t BIN_MAX, TString addcut="")
{
  TLatex* tex = new TLatex(0.18,0.96,triggerpass);
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);

  TH1D* hinclusive = new TH1D(Form("hinclusive_%s_%s",triggerpass.Data(),varname.Data()),Form(";Matched reco D^{0} %s;Candidates",varlatex.Data()),BIN_NUM,BIN_MIN,BIN_MAX);
  inttree->Project(Form("hinclusive_%s_%s",triggerpass.Data(),varname.Data()),variable,Form("%s%s",prefilter.Data(),addcut.Data()));
  hinclusive->Sumw2();
  if(triggerpass=="HLT_DmesonTrackingGlobalPt8_Dpt20_v1"||triggerpass=="HLT_DmesonPPTrackingGlobal_Dpt20_v1")
    {
      TCanvas* chinclusive = new TCanvas(Form("chinclusive_%s",varname.Data()),"",500,500);
      hinclusive->Draw();
      hinclusive->SetStats(0);
      tex->Draw();
      if(isPbPb) chinclusive->SaveAs(Form("triggerturnonPlots/pthat%.0f/MBseed/pbpb/chinclusive_%s.pdf",pthat,varname.Data()));
      else chinclusive->SaveAs(Form("triggerturnonPlots/pthat%.0f/MBseed/pp/chinclusive_%s.pdf",pthat,varname.Data()));
    }
  TH2D* hempty = new TH2D(Form("hempty_%s_%s",triggerpass.Data(),varname.Data()),Form(";Matched reco D^{0} %s;Pass efficiency (ZB seed)",varlatex.Data()),BIN_NUM,BIN_MIN,BIN_MAX,10,0,1.2);
  hempty->SetStats(0);
  TH1D* hMBseed = new TH1D(Form("h%s_MBseed_%s",triggerpass.Data(),varname.Data()),"",BIN_NUM,BIN_MIN,BIN_MAX);
  inttree->Project(Form("h%s_MBseed_%s",triggerpass.Data(),varname.Data()),variable,Form("%s%s&&%s",prefilter.Data(),addcut.Data(),triggerpass.Data()));
  hMBseed->Sumw2();
  TEfficiency* pEffMBseed = new TEfficiency(*hMBseed,*hinclusive);
  TCanvas* cMBseed = new TCanvas(Form("c%s_MBseed_%s",triggerpass.Data(),varname.Data()),"",500,500);
  hempty->Draw();
  pEffMBseed->Draw("PSAME");
  tex->Draw();
  /*
  if(isPbPb) cMBseed->SaveAs(Form("triggerturnonPlots/pthat%.0f/MBseed/pbpb/c%s_MBseed_%s.pdf",pthat,triggerpass.Data(),varname.Data()));
  else cMBseed->SaveAs(Form("triggerturnonPlots/pthat%.0f/MBseed/pp/c%s_MBseed_%s.pdf",pthat,triggerpass.Data(),varname.Data()));
  */
  cMBseed->SaveAs(Form("triggerturnonPlots/pthat%.0f/ZBseed/pp/c%s_MBseed_%s.pdf",pthat,triggerpass.Data(),varname.Data()));
}

double GetMultAssMean(int jass)
{
                   TH2D* hmultass = (TH2D*) gInputfile->Get(Form("%smult_ass_%d",subdir.Data(),jass));
                   TH1D* hmultass_proj = (TH1D*)hmultass->ProjectionY(Form("mult_ass_proj_%d",jass),-1,-1);
                   double  multass_mean = hmultass_proj->GetMean();
                   return multass_mean;
}

Double_t MakeHistsForSlides(TString goodFile, TString badFile, TString rootFile)
{
  //Extracting the histogram
  TFile* f = TFile::Open(rootFile);
  TDirectoryFile* dirFile = (TDirectoryFile*)(f->Get("AliAnalysisTaskCalibEmcal"));
  TList* list = (TList*)(dirFile->Get("histogram"));
  TH2D* cellVAmplitude = (TH2D*)(list->FindObject("_histogram_cell_id_amplitude"));

    
  //Getting a Good Cell
  TH1D* goodC = new TH1D();
  ifstream myFile(goodFile);
  TString title = "";

  const int num = 500;

  if(myFile.is_open())
    {
      string line = "";
      std::getline(myFile,line);

      int index = 0;
      double chi2 = 0.0;
      
      while(myFile >> index >> chi2 && index < num)
	{
	  //cout << index << "\t" << chi2 << endl;
	  title = Form("Energy distribution of Cell %d",index);
	  TH1D* goodCell = cellVAmplitude->ProjectionY("goodCell", index+1,index+1);
	  goodC = (TH1D*)goodCell->Clone();
	}
    }

void drawArrowPlot(TChain *mc, TString varU, TString deltaU, TString varV, TString deltaV, Int_t nbinU, Double_t minU, Double_t maxU, Int_t nbinV, Double_t minV, Double_t maxV){

  TString uvname = generateRandomName();
  TString uname = generateRandomName();
  TString vname = generateRandomName();
  TH2D * UV = new TH2D(uvname,uvname,10*nbinU,minU,maxU,10*nbinV,minV,maxV);
  TProfile2D * dU = new TProfile2D(uname,uname,nbinU,minU,maxU,nbinV,minV,maxV);
  TProfile2D * dV = new TProfile2D(vname,vname,nbinU,minU,maxU,nbinV,minV,maxV);

  mc->Draw(varV+":"+varU+">>"+uvname,"(isAssoc==1)","goff");
  mc->Draw(deltaU+":"+varV+":"+varU+">>"+uname,"(isAssoc==1)","goffprof");
  mc->Draw(deltaV+":"+varV+":"+varU+">>"+vname,"(isAssoc==1)","goffprof");
  
  SetFancyGrayscalePalette();
  UV->Draw("colsame");


  for (Int_t iU=1; iU<nbinU+1; iU++){
    for (Int_t iV=1; iV<nbinV+1; iV++){
      
      Double_t uu = dU->GetXaxis()->GetBinCenter(iU);
      Double_t vv = dU->GetYaxis()->GetBinCenter(iV);
      Double_t du = dU->GetBinContent(iU,iV);
      Double_t dv = dV->GetBinContent(iU,iV);
      Double_t due = dU->GetBinError(iU,iV);
      Double_t dve = dV->GetBinError(iU,iV);
      
      drawArrow(uu,vv,du,due,dv,dve);

    }
  }
}

TH1* getHisto(char * filename, char* histoName, char * dirName, int nBin, double lumi)
{
  TH1 * hpt_=0;
  TFile *file0 = TFile::Open(filename);
  if(!file0) return hpt_;
  TDirectory *dir;
  TH2D * hMuPt;

  if(dirName == "0") {
    hMuPt = (TH2D*) file0->Get(histoName);
  } else {
    dir = (TDirectory*) file0->Get(dirName);
    if(!dir) return hpt_;
    hMuPt = (TH2D*) dir->Get(histoName);
  }

  if(hMuPt) {
  hpt_ = (TH1*) hMuPt->Clone();
  hpt_->Sumw2();
  hpt_->Scale(1./lumi); // this take into into account the luminosity
  hpt_->SetLineWidth(2);
  hpt_->Rebin(nBin);
  double nBinX=hpt_->GetNbinsX();

  // overFlow
  hpt_->SetBinContent((int)nBinX,hpt_->GetBinContent((int)nBinX)+hpt_->GetBinContent((int)nBinX+1));
  hpt_->SetDirectory(0);
  file0->Close();
  hpt_->SetLineWidth(3);
  }
  return hpt_;
}

TProfile* GetMultTrgProfileX(int itrg)
{
                   TH2D* hmulttrg = (TH2D*) gInputfile->Get(Form("%smult_trg_%d",subdir.Data(),itrg));
                   hmulttrg->Rebin2D(8,1);
                   TProfile* hmulttrg_profx = (TProfile*)hmulttrg->ProfileX(Form("mult_trg_%d_profx",itrg),-1,-1);
                   return hmulttrg_profx;
}

void paint(TString a, TString b)
{
	TChain *t = new TChain(a);
	t->Add(a + "_" + b +".root");

	std::cout << "painting the histograms for: " << a + "_" + b +".root" << std::endl;

	TDirectory *dir = new TDirectory("dir", "dir");
	dir->cd();

//	histo = new TH2D("jets", "jets", 50, -2.5, 2.5, 60, 10, 610);
//	histo = new TH2D("jets", "jets", 50, -2.5, 2.5, 40, 4.7004803657924166, 6.404803657924166);
//	TH2D * histo = new TH2D("jets", "jets", 25, -2.5, 2.5, 20, 4.0943445622221004, 6.1943445622221004); 
	TH2D * histo = new TH2D("jets", "jets", 50, -2.5, 2.5, 40, 4.0943445622221004, 6.1943445622221004);//original
//	TH2D * histo = new TH2D("jets", "jets", 50, -2.5, 2.5, 40, 4.0943445622221004, 7.8);
	histo->SetDirectory(dir);

	//the varexp part of the draw synatx means: draw log(jetPt+50) versus jetEta and append the existing ("+" -> avoid recreation) histogram called "jets"
	//selection is an expression with a combination of the Tree variables -> no selection applied in this case ""
	//option is the drawing option -> if option contains the string "goff", no graphics is generated.
	//fourth and fifth arguments are: Int_t nevents, Int_t firstevent
	t->Draw("log(jetPt+50):jetEta >> +jets", "", "Lego goff");
	//std::cout <<"jetPt " << log(jetPt+50) << " and jetEta " << jetEta << std::endl;


	std::cout << "saving the histograms: " << a + "_" + b +"_histo.root" << std::endl;
	TFile g(a + "_" + b +"_histo.root", "RECREATE");
	histo->SetDirectory(&g);
	delete dir;

	g.cd();
	histo->Write();
	g.Close();
}

TestProblem
MonteCarloConvolution(const int m,
                      const int n,
                      const double xm1,
                      const double xm2,
                      const double xt1,
                      const double xt2,
                      TF1 *truthFn,
                      TF1 *kernelFn,
                      const int nEvents)
{
  static int id = 0; id++;
  TestProblem t;
  double dm = (xm2-xm1)/m;
  double dt = (xt2-xt1)/n;
  TMatrixD R(m,n);

  // Discretize the kernel to fill R
  for (int i=0; i<m; i++)
    for (int j=0; j<n; j++)
      R(i,j) = kernelFn->Eval(xt1+j*dt-i*dm);
  t.Response = MatrixUtils::Matrix2Hist(R, Form("R%d",id), xm1, xm2, xt1, xt2);
  t.Response->SetTitle(Form("%d x %d convolution matrix;"
                            "s (observed);t (true)", m, n));

  TH2D *RMC = new TH2D(Form("R_MC%d",id), "A_{MC}",
                       m,xm1,xm2,n,xt1,xt2);

  // Model a true and a measured distribution
  // There is no bIdeal for this problem.
  t.xTruth    = new TH1D("hTrue",    "",         n, xt1, xt2);
  t.xTruthEst = new TH1D("hTrueEst", "hTrueEst", n, xt1, xt2);
  t.bNoisy    = new TH1D("hMeas",    "hMeas",    m, xm1, xm2);
  t.xTruthEst->Sumw2();
  t.bNoisy->Sumw2();


  for (Int_t i=0; i<nEvents; i++)
  {
    Double_t xt = truthFn->GetRandom();
    t.xTruthEst->Fill(xt);
    Double_t xm = xt + kernelFn->GetRandom();
    t.bNoisy->Fill(xm);
    RMC->Fill(xm,xt);
  }

  // MC efficiency histogram
  t.eff = RMC->ProjectionY("eff",1,m);
  t.eff->Divide(t.xTruthEst);

  // Exact truth histogram
  for (int j=1; j<=n; j++)
  {
    double val = truthFn->Eval(t.xTruth->GetBinCenter(j));
    t.xTruth->SetBinContent(j, val);
  }
  t.xTruth->Scale(nEvents/t.xTruth->Integral());

  return t;
}

/**
 * Loads a histogram from a ROOT-file, scales it, and either returns it as \c
 * returnedHistogram or adds it to \c returnedHistogram . The first if
 * \c returnedHistogram is \c NULL otherwise the latter.
 *
 * @param histogramName
 *                              name of the histogram.
 * @param inputFilename
 *                              name of the ROOT-file.
 * @param scale
 *                              scale factor for histogram.
 * @param returnedHistogram
 *                              the returned histogram (input, pass-by-reference does not work).
 * @param debug
 *                              switch on debug output, defaults to \c false
 *
 * @return
 *                              the returned histogram (output, pass-by-reference does not work).
 */
void
LoadHistogramTH2D(const TString& histogramName, const TString& inputFilename, double scale, TH2D*& returnedHistogram, bool debug = false)
{
// {{{
    TFile inputFile(inputFilename);
    if (!inputFile.IsOpen())
    {
        cerr << "Could not open '" << inputFilename << "' for reading." << endl;
    }
    else
    {
        TH2D* histogram = dynamic_cast<TH2D*>( inputFile.Get(histogramName) );
        if (!histogram)
        {
            cerr << "No histogram named '" << histogramName << "' in file '" << inputFilename << "'" << endl;
        }
        else
        {
            if (debug) cerr << inputFilename << " " << histogramName << " entries=" << histogram->GetEntries() << " integral=" << histogram->Integral() << " scale*integral=" << scale*histogram->Integral() << endl;
            histogram->Scale(scale);
            if (!returnedHistogram)
            {
                returnedHistogram = new TH2D(*histogram);
                returnedHistogram->SetDirectory(0); // otherwise "TFile inputFile" owns this returnedHistogram and deletes it when "TFile inputFile" goes out of scope
            }
            else
            {
                returnedHistogram->Add(histogram);
            }
        }
        inputFile.Close();
    }
// }}}
};

void paint(TString dir, TString a, TString b)
{
	TChain *t = new TChain(a);
	t->Add(dir + a + "_" + b +".root");

	std::cout << "painting the histograms for: " << dir + a + "_" + b +".root" << std::endl;

	TDirectory *direc = new TDirectory("dir", "dir");
	direc->cd();

	TH2D * histo = new TH2D("jets", "jets", 50, -2.5, 2.5, 40, 4.17438727, 6.95654544315);//pt starting from 15 and until 1000
	histo->SetDirectory(direc);

	t->Draw("log(jetPt+50):jetEta >> +jets", "", "Lego goff");

  TH2D * histo_lin = new TH2D("jets_lin", "jets_lin", 50, -2.5, 2.5, 40, 15, 1000);//pt starting from 15 and until 1000
	t->Draw("jetPt:jetEta >> +jets_lin", "", "Lego goff");

	std::cout << "saving the histograms: " << a + "_" + b +"_histo.root" << std::endl;
	TFile g(a + "_" + b +"_histo.root", "RECREATE");
	histo->SetDirectory(&g);
	histo_lin->SetDirectory(&g);
	delete direc;

	g.cd();
	histo->Write();
	histo_lin->Write();
	g.Close();

}

void plotTowers_dj(
    TString inFile0Name="~/scratch01/ana/Hydj_BSC_HF_L1Emul/oh0928_v2/openhlt_hiReco_RAW2DIGI_RECO_*_djana.root",
    TString outdir = "out/Hydj/BSC_HF_L1Emul/devHLT381")
{
  TChain * djTree = new TChain("dijetAna_data_calojet_tower/djTree");
  djTree->Add(inFile0Name);
  aliases_dj(djTree);
  cout << " # entries: " << djTree->GetEntries() << endl;

  gSystem->mkdir(outdir.Data(),kTRUE);
  Int_t NTowerMax=70;


  pair<TString,TString> t = make_pair("name","cut");
  vector<pair<TString,TCut> > trigs;
  trigs.push_back(make_pair("All","1==1"));
  trigs.push_back(make_pair("HfAnyHit","hlt[2]"));
  trigs.push_back(make_pair("HfCoic1","hlt[3]"));
  trigs.push_back(make_pair("HfCoic2","hlt[4]"));

  for (unsigned int i=0; i<trigs.size(); ++i) {
    TCanvas * cHfNvsP = new TCanvas("HfNvsP_"+trigs[i].first,"HfNvsP_"+trigs[i].first,500,500);
    //cHfNvsP->SetLogx();
    //cHfNvsP->SetLogy();
    TH2D * hHfNvsP = new TH2D("hHfNvsP"+trigs[i].first,"hHfNvsP_",NTowerMax,0,NTowerMax,NTowerMax,0,NTowerMax);
    hHfNvsP->SetAxisRange(0.1,400,"X");
    hHfNvsP->SetAxisRange(0.1,400,"Y");
    djTree->Draw("NHfN:NHfP>>hHfNvsP"+trigs[i].first,trigs[i].second,"colz");
    hHfNvsP->SetTitle("HF Towers;# HF+ Towers (E>3GeV);# HF- Towers (E>3GeV)");
    cHfNvsP->Print(Form("%s/trHfNvsP_%s.gif",outdir.Data(),trigs[i].first.Data()));
  }
}

void CaloHitAna::DrawPID(bool cut)
{
   if (fChain == 0) return;

   TH2D *hpid = new TH2D("qpid", "QPID", 250, 0, 4E5, 250, 0, 4E5);
   TCutG *cut_proton = NULL;

   if(cut)
   {
      TFile *f = TFile::Open("qpid_cuts.root", "read");
      cut_proton = dynamic_cast<TCutG*>(f->Get("qpid_proton"));
   }

   Long64_t nentries = fChain->GetEntriesFast();

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

      for(i = 0; i < CaloHit_; i++)
      {
         if(!cut_proton || cut_proton->IsInside(CaloHit_fNf[i], CaloHit_fNs[i]))
            hpid->Fill(CaloHit_fNf[i], CaloHit_fNs[i]);
      }

   }

   hpid->Draw("colz");
}

TH1* processFile(TString fname, double scaleFactor){

  static int hcounter=0;

  cout<<"Doing fname="<<fname<<endl;
  TFile *_file2 = TFile::Open(fname);
  PS->cd();
  TString plot =  "lLV[0].lpfIso:lLV[0].fP.Eta()";
  TString hname = Form("h%i",hcounter);
  jets2p->Draw(plot +" >> " + hname+"(24,-3,3,100,0,2)","","qn");
  TH2D *aux = (TH2D*)gDirectory->Get(hname);
  //aux->Scale(scaleFactor);
  TH1D * ss = aux->ProjectionX(hname+"ss",1,10); // signal region 
  TH1D * ni = aux->ProjectionX(hname+"ni",16,100);// nonIso region (use 11
  //cout<<"h="<<hcounter<<" ss="<<ss->GetEntries()<<" ni="<<ni->GetEntries()<<endl;

  TString statis = Form("( %i / %i )",(int)ss->GetEntries(),(int)ni->GetEntries());
  ss->Divide(ni);
  ss->GetYaxis()->SetTitle(TString("N_{signal}/N_{nonIso}")+statis);
  ss->GetXaxis()->SetTitle("#eta");
  
  hcounter++;

  return ss;

}

TGraph* getContourFilledX(TH2D* inputHisto, TCanvas* goodCanvas, int Width, int Style, int FillStyle, double X){
       TCanvas* c1 = new TCanvas("temp", "temp",600,600);

       TH2D* histo = (TH2D*)inputHisto->Clone("temp");

       double levels[] = {X};
       histo->SetContour(1, levels);
       histo->Draw("CONT LIST");
       c1->Update();
       TObjArray* contours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
       Int_t ncontours     = contours->GetSize();
       TList *list         = (TList*)contours->At(0);
       delete c1;

       goodCanvas->cd();
       printf("list size = %i\n", (int)list->GetSize());
       if(list->GetSize()<=0)return new TGraph(0);

       for(unsigned int i=0;i<list->GetSize();i++){
       TGraph* EXCLUSION   = (TGraph*)(list->At(i)->Clone("copy"));
       EXCLUSION->SetLineColor(1);
       EXCLUSION->SetLineWidth(Width);
       EXCLUSION->SetLineStyle(Style);
       EXCLUSION->SetFillColor(kBlack);
       EXCLUSION->SetFillStyle(FillStyle);
       //EXCLUSION->Draw("CL F same");
       }

       return EXCLUSION;
}