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

// ----------------------------------------------------------------------
void plotRatioRaw(string era = "2016BF") {
  gStyle->SetOptStat(0);

  if (era == "all") {
    plotRatioRaw("2011");
    plotRatioRaw("2012");
    plotRatioRaw("2016BF");
    plotRatioRaw("2016GH");
    return;
  }


  double eps(0.00001);
  TH1D *hresult = new TH1D("hresult", Form("%s", era.c_str()), 80, 0., 0.4);
  hresult->GetXaxis()->SetTitle("BDT > ");
  hresult->SetMarkerStyle(24); hresult->SetLineColor(kRed+2); hresult->SetMarkerColor(kRed+2);
  TH1D *hsresult = new TH1D("hsresult", Form("%s", era.c_str()), 80, 0.+eps, 0.4+eps);
  hsresult->SetMarkerStyle(25); hsresult->SetLineColor(kBlue);  hsresult->SetMarkerColor(kBlue);

  TH1D *hbdtcut = new TH1D("hbdtcut", Form("%s", era.c_str()), 80, 0., 0.4);
  hbdtcut->SetMarkerColor(kRed+2);  hbdtcut->SetMarkerStyle(20);
  TH1D *hsbdtcut = new TH1D("hsbdtcut", Form("%s", era.c_str()), 80, 0.+eps, 0.4+eps);
  hsbdtcut->SetMarkerColor(kBlue);  hsbdtcut->SetMarkerStyle(21);

  invertedMuonIDRaw(era, 0.15, hresult, hsresult);
  invertedMuonIDRaw(era, 0.20, hresult, hsresult);
  invertedMuonIDRaw(era, 0.25, hresult, hsresult);
  invertedMuonIDRaw(era, 0.30, hresult, hsresult);
  invertedMuonIDRaw(era, -99., hbdtcut, hsbdtcut);

  c0->Clear();
  double ymax(50.);
  double smax = hresult->GetMaximum();
  double dmax = hsresult->GetMaximum();
  double themax = (smax>dmax? smax: dmax);
  if (themax < 50.)  ymax = 50.;
  if (themax < 20.)  ymax = 30.;
  if (themax < 5.)  ymax = 10.;
  if (themax < 0.5) ymax = 1.0;
  hresult->SetMinimum(-ymax);
  hresult->SetMaximum(ymax);
  hresult->Draw();
  hsresult->Draw("same");
  hbdtcut->Draw("same");
  hsbdtcut->Draw("same");
  pl->DrawLine(0., 0., 0.4, 0.);

  TLegend *tle = new TLegend(0.25, 0.2, 0.50, 0.5);
  tle->SetFillStyle(0);
  tle->SetBorderSize(0);
  tle->SetHeader(era.c_str());
  tle->AddEntry(hresult, "yield at B0", "p");
  tle->AddEntry(hsresult, "yield at Bs", "p");
  tle->Draw();

  c0->SaveAs(Form("invertedMuonIDRaw-%s-result.pdf", era.c_str()));
}

void testCorr3D(Int_t corrLevel=0,
    TString mod="B2InConeGluon",
    TString modref="B2InConeQuark",
    TString app="_jtv3_2_cv4qg",
    TString appref="_jtv3_2_cv4qg",
    Int_t isample=0, // -1 for all samples
    Int_t etaPM=2 // +/- 2 for |eta|<1
    )
{
  Int_t mode=1; // 0 for write, 1 for read
  Int_t cbin=0;
  Float_t ptHatMin=110;

  Corrector3D trkCorr("trkCorrHisAna_djuq",app,mod);
  trkCorr.ptHatMin_.clear();
  trkCorr.ptHatMin_.push_back(ptHatMin);
  trkCorr.sampleMode_ = 0; // 0 for choosing individual sample, 1 for merge samples
  trkCorr.Init();
  Double_t corr[4];
  cout << trkCorr.GetCorr(10,0,110,0,corr) << endl;

  Corrector3D trkCorrRef("trkCorrHisAna_djuq",appref,modref);
  trkCorrRef.ptHatMin_.clear();
  trkCorrRef.ptHatMin_.push_back(ptHatMin);
  trkCorrRef.sampleMode_ = 0; // 0 for choosing individual sample, 1 for merge samples
  trkCorrRef.Init();

  /*
  // Plot 2D Corr.
  gStyle->SetPadRightMargin(0.15);

  TCanvas * c2 = new TCanvas("c2","c2",500,500);
  TH2D *hCorr2D = (TH2D*)trkCorr.InspectCorr(corrLevel,isample,cbin,4,30,0);
  hCorr2D->Draw("colz");
  */

  TCanvas * cEff = new TCanvas("cEff","cEff",500,500);
  TH1D * hCorrPtRef = (TH1D*) trkCorrRef.InspectCorr(corrLevel,isample,cbin,3,15,2,7-etaPM,7+etaPM);
  hCorrPtRef->SetAxisRange(0,120,"X");
  hCorrPtRef->SetMarkerStyle(0);
  hCorrPtRef->SetMarkerColor(kBlack);
  hCorrPtRef->SetLineColor(kBlack);
  hCorrPtRef->Draw("histE");
  TH1D * hCorrPt = (TH1D*) trkCorr.InspectCorr(corrLevel,isample,cbin,3,15,2,7-etaPM,7+etaPM);
  hCorrPt->SetMarkerStyle(0);
  hCorrPt->SetMarkerColor(kRed);
  hCorrPt->SetLineColor(kRed);
  hCorrPt->Draw("hist same E");

  /*
  TCanvas * cEffEta = new TCanvas("cEffEta","cEffEta",500,500);
  TH1D * hCorrEta = (TH1D*)trkCorr.InspectCorr(corrLevel,isample,cbin,4,30,1);
  hCorrEta->Draw("histE");
  */
}

void anaFragShape(TString infile="dj_HCPR-GoodTrkAndPixel_CleanEvt1130.root")
{
  TChain * djcalo = new TChain("djcalo/djTree");
  djcalo->Add(infile);
  aliases_dijet(djcalo);
  TChain * djcalopix = new TChain("djcalo_pxl/djTree");
  djcalopix->Add(infile);
  aliases_dijet(djcalopix);
  TString evtSel("(cent<10 && nljet>120 && abs(nljeta)<2 && aljet>50 && abs(aljeta)<2 && jdphi>2.5 && Aj>0.24)");
  TH1::SetDefaultSumw2();

  TCanvas * c0 = new TCanvas("c0","c0",500,500);
  djcalo->Draw("Aj>>hAj(20,0,1)",evtSel);
  djcalopix->Draw("Aj",evtSel,"Esame");

  TCanvas * c2 = new TCanvas("c2","c2",500,500);
  TH1D * hJDPhi = new TH1D("hJDPhi","hJDPhi",50,0,3.1416);
  TH1D * hJDPhi2 = new TH1D("hJDPhi2","hJDPhi",50,0,3.1416);
  Float_t numDJ = djcalo->Draw("jdphi>>hJDPhi",evtSel);
  Float_t numDJ2 = djcalopix->Draw("jdphi>>hJDPhi2",evtSel);
  cout << "num djs: " << numDJ << "  djs(pix)" << numDJ2 << endl;

  TH1D * hPNDRTrk = new TH1D("hPNDRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPNDRPix = new TH1D("hPNDRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRTrk = new TH1D("hPADRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRPix = new TH1D("hPADRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  djcalo->Draw("pndr>>hPNDRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalopix->Draw("pndr>>hPNDRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalo->Draw("padr>>hPADRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  djcalopix->Draw("padr>>hPADRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  hPNDRTrk->Scale(1./(numDJ*hPNDRTrk->GetBinWidth(1)));
  hPNDRPix->Scale(1./(numDJ*hPNDRPix->GetBinWidth(1)));
  hPADRTrk->Scale(1./(numDJ*hPADRTrk->GetBinWidth(1)));
  hPADRPix->Scale(1./(numDJ*hPADRPix->GetBinWidth(1)));
  hPNDRPix->SetMinimum(50);
  hPNDRPix->SetMaximum(100);
  hPNDRTrk->SetMinimum(10);
  hPNDRTrk->SetMaximum(2000);
  hPNDRPix->SetMarkerStyle(kFullCircle);
  hPNDRTrk->SetMarkerStyle(kFullSquare);
  hPNDRPix->SetMarkerColor(kRed);
  hPNDRTrk->SetMarkerColor(kRed);
  hPADRPix->SetMarkerStyle(kOpenCircle);
  hPADRTrk->SetMarkerStyle(kOpenSquare);
  hPADRPix->SetMarkerColor(kBlue);
  hPADRTrk->SetMarkerColor(kBlue);
  TCanvas * c3 = new TCanvas("c3","c3",500,500);
  hPNDRPix->Draw("E");
  hPADRPix->Draw("Esame");
  TCanvas * c3_2 = new TCanvas("c3_2","c3_2",500,500);
  c3_2->SetLogy();
  hPNDRTrk->Draw("E");
  hPADRTrk->Draw("Esame");
}

void ptDependence()
{
   TFile *inf = new TFile("histos/ppMC.root");
   TTree *t = (TTree*) inf->Get("nt");

   const int nBins = 4;
   double ptBin[nBins+1] = {100,120,140,160,200};
//   const int nBins = 1;
//   double ptBin[nBins+1] = {100,400};
   
   TH1D *hProb = new TH1D("hProb","",nBins,ptBin);
   TH1D *hCSV = new TH1D("hCSV","",nBins,ptBin);
   TH1D *hSVTXM = new TH1D("hSVTXM","",nBins,ptBin);
   TProfile *pGen = new TProfile("pGen","",nBins,ptBin);
   
   for (int n=0; n<nBins;n++)
   {
      RooRealVar f1 = bfractionFit("discr_prob",0,3.5,ptBin[n],ptBin[n+1]);
      RooRealVar f2 = bfractionFit("discr_csvSimple",0,1,ptBin[n],ptBin[n+1]);
      RooRealVar f3 = bfractionFit("svtxm",0,6,ptBin[n],ptBin[n+1]);
      hProb->SetBinContent(n+1,f1.getVal());    
      hProb->SetBinError(n+1,f1.getError());    
      hCSV->SetBinContent(n+1,f2.getVal());    
      hCSV->SetBinError(n+1,f2.getError());    
      hSVTXM->SetBinContent(n+1,f3.getVal());    
      hSVTXM->SetBinError(n+1,f3.getError());    
   }
   
   TCanvas *c2 = new TCanvas("c2","",600,600);
   hProb->SetAxisRange(0,0.05,"Y");
   hProb->SetXTitle("Jet p_{T} (GeV/c)");
   hProb->SetYTitle("b-jet fraction");
   hProb->SetTitleOffset(1.5,"Y");
   hProb->Draw();
   hCSV->SetLineColor(2);
   hCSV->SetMarkerColor(2);
   hCSV->SetMarkerStyle(24);
   hCSV->Draw("same");
   hSVTXM->SetLineColor(4);
   hSVTXM->SetMarkerColor(4);
   hSVTXM->SetMarkerStyle(24);
//   hSVTXM->Draw("same");
   t->Draw("abs(refparton_flavorForB)==5:jtpt","","prof same");
   
   TLegend *leg = new TLegend(0.2,0.7,0.5,0.9);
   leg->SetBorderSize(0);
   leg->SetFillStyle(0);
   leg->SetFillColor(0);
   leg->AddEntry(hProb,"Jet Probability","pl");
   leg->AddEntry(hCSV,"CSV","pl");
//   leg->AddEntry(hSVTXM,"SV mass","pl");
   leg->Draw();
}

void Drawhisto() {
    gStyle->SetOptStat(kFALSE);
    gStyle->SetErrorX(0);
    TFile *fin = TFile::Open("PythiaMBhistos.root");
    TH1D* hratioplus = (TH1D*)fin->Get("NjetratioPlus");
    TCanvas *c1 = new TCanvas();
    c1->SetTickx(1);
    c1->SetTicky(1);
    hratioplus->GetXaxis()->SetTitle("#Sigma E_{T} (2.9<#eta<5.2) (GeV)");
    hratioplus->GetXaxis()->SetRangeUser(0,30);
    hratioplus->GetYaxis()->SetRangeUser(0,0.07);
    hratioplus->GetYaxis()->SetTitle("N_{jet}/N_{evt}");
    hratioplus->SetMarkerSize(1);
    hratioplus->SetMarkerStyle(20);
    hratioplus->SetMarkerColor(1);
    hratioplus->SetLineColor(1);
    hratioplus->Draw("P");
    TLatex *tl = new TLatex();
    tl->SetTextSize(0.045);
    tl->SetNDC();
    TLegend *tg = new TLegend(0.15,0.7,0.6,0.75);
    tg->SetTextSize(0.045);
    tg->SetFillColor(0);
    tg->SetBorderSize(0);
    tg->AddEntry(hratioplus,"5.02TeV, Z2Star Tune");
    tg->Draw("same");
    tl->DrawLatex(0.2,0.85,"Pythia 6");
    tl->DrawLatex(0.2,0.8,"anti-k_{T} R=0.3 PF jets, E_{T} > 20 GeV, |#eta| < 2.8");
    c1->Print("jetYratioHFPlus.png");
    c1->Print("jetYratioHFPlus.pdf");
}

void mult()
{
  //TFile* f = new TFile("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_UCC2011_pixeltracks_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin500_centmax1000.root");
  TFile* f = new TFile("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_UCC2011_pixeltracks_INCLMULTUCC010_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root");

  TNtuple* h = (TNtuple*)f->Get("higeninfontuple");
 
  TH1D* hMultAll = new TH1D("hMultAll",";nMult",300,0,6000); 
  TH1D* hMult = new TH1D("hMult",";nMult",300,0,6000);
  TH1D* hMult010 = new TH1D("hMult010",";nMult",300,0,6000);
  TH1D* hMult001 = new TH1D("hMult001",";nMult",300,0,6000);

  higeninfontuple->Draw("nmult>>hMultAll","","GOFF");
  higeninfontuple->Draw("nmult>>hMult","(7.0*hft+zdc)<36000","GOFF");
  higeninfontuple->Draw("nmult>>hMult010","(7.0*hft+zdc)<36000 && npixel>51400 && hft>3260","GOFF");
  higeninfontuple->Draw("nmult>>hMult001","(7.0*hft+zdc)<36000 && npixel>53450 && hft>3393","GOFF");

  hMult->SetMarkerColor(kRed);
  hMult010->SetMarkerColor(kBlue);
  hMult001->SetMarkerColor(kGreen);

  TCanvas* c = new TCanvas("c","c",550,500);
  c->SetLogy();
  hMultAll->Draw("PE");
  hMult->Draw("PESAME");
  hMult010->Draw("PESAME");
  hMult001->Draw("PESAME");
}

TH1D* Prediction::GetRatio(TH1D* h1, TH1D* h2, int mergenbins, bool binomial){
	// compute Zll/Photon ratio
	*fLogStream << "------------------------------------------------------------------------" << endl;
	*fLogStream << "Compute ratio between " << h1->GetName() << " and " << h2->GetName()      << endl;           

	
	if(h1==0 || h2==0 ){
		cout << "GetMCZnunuToPhotonRatio: received 0 pointer!" << endl;
		exit(-1);
	}

	// GetScaled histos with only one bin and probagated errors: stat error and error on scale factor
	TH1D *currh1   = GetScaledHisto(h1, 1, 0, mergenbins); 
	TH1D *currh2   = GetScaledHisto(h2, 1, 0, mergenbins); 

	TString name = h1->GetName();
	name         += "_";
	name         += h2->GetName();
	name         += "_Ratio";
	TH1D *ratio = (TH1D*) currh1->Clone(name);
	ratio->Divide(currh1, currh2, 1,1, binomial?"B":""); // binomial errors
	
	ratio->SetLineColor(kBlack);
	ratio->SetMarkerColor(kBlack);
	ratio->SetMarkerStyle(4);

	if(fWriteToFile){
		ratio ->Write();
	}
	delete currh1;
	delete currh2;
	return ratio;
}

// ----------------------------------------------------------------------
void plotRatio(string era = "2016BF") {
  gStyle->SetOptStat(0);

  if (era == "all") {
    plotRatio("2011");
    plotRatio("2012");
    plotRatio("2016BF");
    plotRatio("2016GH");
    return;
  }


  TH1D *hresult = new TH1D("hresult", Form("n(gauss(5.28))/n(gauss(5.37)) %s", era.c_str()), 40, 0., 0.4);
  hresult->GetXaxis()->SetTitle("BDT > ");
  TH1D *hbdtcut = new TH1D("hbdtcut", Form("n(gauss(5.28))/n(gauss(5.37)) %s", era.c_str()), 40, 0., 0.4);
  hbdtcut->SetMarkerColor(kBlue);
  invertedMuonID(era, 0.15, hresult);
  invertedMuonID(era, 0.20, hresult);
  invertedMuonID(era, 0.25, hresult);
  invertedMuonID(era, 0.30, hresult);
  invertedMuonID(era, -99., hbdtcut);

  c0->Clear();
  double ymax(50.);
  if (hresult->GetMaximum() < 5.)  ymax = 10.;
  if (hresult->GetMaximum() < 0.5) ymax = 1.0;
  hresult->SetMinimum(-ymax);
  hresult->SetMaximum(ymax);
  hresult->Draw();
  hbdtcut->Draw("same");
  pl->DrawLine(0., 0., 0.4, 0.);
  c0->SaveAs(Form("invertedMuonID-%s-result.pdf", era.c_str()));
}

// ----------------------------------------------------------------------
void plotRatio2Comb(string era = "2016BF") {
  gStyle->SetOptStat(0);

  if (era == "all") {
    plotRatio2Comb("2011");
    plotRatio2Comb("2012");
    plotRatio2Comb("2016BF");
    plotRatio2Comb("2016GH");
    return;
  }

  double eps(0.00001);
  TH1D *hresultBd = new TH1D("hresultBd", era.c_str(), 80, 0., 0.4);
  TH1D *hresultBs = new TH1D("hresultBs", era.c_str(), 80, 0.+eps, 0.4+eps);
  hresultBd->GetYaxis()->SetTitle("a.u.");
  hresultBd->GetXaxis()->SetTitle("BDT > "); hresultBd->SetMarkerStyle(24); hresultBd->SetLineColor(kRed+2); hresultBd->SetMarkerColor(kRed+2);
  hresultBs->GetXaxis()->SetTitle("BDT > "); hresultBs->SetMarkerStyle(25); hresultBs->SetLineColor(kBlue);  hresultBs->SetMarkerColor(kBlue);
  TH1D *hbdtcutBd = new TH1D("hbdtcutbd", "", 80, 0., 0.4);
  TH1D *hbdtcutBs = new TH1D("hbdtcutBs", "", 80, 0.+eps, 0.4+eps);
  hbdtcutBs->SetMarkerColor(kBlue);  hbdtcutBd->SetMarkerStyle(24); hbdtcutBs->SetLineColor(kBlue);
  hbdtcutBd->SetMarkerColor(kRed+2); hbdtcutBd->SetMarkerStyle(25); hbdtcutBd->SetLineColor(kRed+2);

  invertedMuonID2Comb(era, 0.15, hresultBd, hresultBs);
  invertedMuonID2Comb(era, 0.20, hresultBd, hresultBs);
  invertedMuonID2Comb(era, 0.25, hresultBd, hresultBs);
  invertedMuonID2Comb(era, 0.30, hresultBd, hresultBs);
  invertedMuonID2Comb(era, -99., hbdtcutBd, hbdtcutBs);

  c0->Clear();
  double ymax(5.);
  hresultBd->SetMinimum(-ymax);
  hresultBd->SetMaximum(ymax);
  hresultBd->Draw();
  hbdtcutBd->Draw("same");
  pl->DrawLine(0., 0., 0.4, 0.);

  tl->SetTextSize(0.035);
  tl->SetTextColor(kRed+2); tl->DrawLatexNDC(0.2, 0.85, "n(gauss(5.28))/combinatorial");
  tl->SetTextColor(kBlue);  tl->DrawLatexNDC(0.2, 0.81, "n(gauss(5.37))/combinatorial");
  tl->SetTextColor(kBlack);

  hresultBs->Draw("same");
  hbdtcutBs->Draw("same");

  c0->SaveAs(Form("invertedMuonID2Comb-%s-result.pdf", era.c_str()));

}

void plotCorrectionGenNested(){

	for(int deta = 0; deta < NdEtaBins+1; deta++){
		dEtaBins[deta] = dEtaBins[deta] - 0.00001;//fix bin boundary
	}


	TFile* file = new TFile("../rootfiles/CME_QvsdEta_pPb_EPOS_GEN_NestedLoop_v3.root");

	TH1D* QvsdEta[48];

	for(int deta = 0; deta < NdEtaBins; deta++){
		  
	  QvsdEta[deta] = (TH1D*) file->Get( Form("ana/QvsdEta_%d",deta) );
			
	}


	TH1D* hist1[3][2];
	TH1D* hist2[3][2];
	for(int sign = 0; sign < 1; sign++){
		for(int HF = 0; HF < 1; HF++){
			hist1[sign][HF] = new TH1D(Form("hist1_%d_%d",sign,HF),"test", NdEtaBins, dEtaBins);
			hist2[sign][HF] = new TH1D(Form("hist2_%d_%d",sign,HF),"test", NdEtaBins, dEtaBins);
		}
	}

	for(int deta = 0; deta < NdEtaBins; deta++){

		double Q_total_real_dEta = QvsdEta[deta]->GetMean();
		double Q_total_real_dEta_error = QvsdEta[deta]->GetMeanError();

		hist1[0][0]->SetBinContent(deta+1, Q_total_real_dEta );
		hist1[0][0]->SetBinError(deta+1,  Q_total_real_dEta_error);
	
	}

	TH1D* base3 = makeHist("base3","like sign(++)","#Delta#eta", "cos(#phi_{1}+#phi_{2}-2#phi_{3})/v2_{3}", 48,0,4.8);
    base3->GetXaxis()->SetTitleColor(kBlack);
    base3->GetYaxis()->SetRangeUser(-0.0015,0.001);
    base3->GetYaxis()->SetTitleOffset(1.9);

    TH1D* base4 = (TH1D*) base3->Clone("base4");

    TCanvas* c4 = makeCanvas("c4","c4");
	gPad->SetTicks();
	gPad->SetLeftMargin(0.20);
	gPad->SetBottomMargin(0.16);
	base3->Draw();
	TH1D* temp = (TH1D*)hist1[0][0]->Clone("temp");
	temp->Scale(1.0/0.20);
	temp->SetMarkerColor(kRed);
	temp->SetLineColor(kRed);
	temp->SetMarkerStyle(20);
	temp->Draw("Psame");


}

void dmesondecaylength()
{


    //gStyle->SetOptStat("nemruoi");
    gStyle->SetTitleSize(.04,"S");
    gStyle->SetOptTitle(1);
    gStyle->SetTitleOffset(1.0,"X");
    gStyle->SetTitleOffset(.88,"Y");
    gStyle->SetTitleSize(.04,"X");
    gStyle->SetTitleSize(.04,"Y");
    gStyle->SetLabelSize(.035,"X");
    gStyle->SetLabelSize(.035,"Y");
    gStyle->SetHistLineWidth(2);
    gStyle->SetOptFit(1);
    gStyle->SetOptStat(0);
// -----   Timer   --------------------------------------------------------
    TStopwatch timer;
    timer.Start();
    // ------------------------------------------------------------------------
    double c= 3* std::pow(10.,8.);
    double mt= 1040*std::pow(10.,-15.); // mean decay time
    double mass= 1.869; // rest mass in GeV/c^2
//TCanvas* can = new TCanvas("can","Radiation Length for start detector",0,0,100,100);
    TCanvas *c1 = new TCanvas("c1", "c1",0,52,1191,692);
    TH1D* h = new TH1D("hist","D-meson, D-meson decay length",24,0,24);
    h->SetTitle("D^{+} meson decay length = c#tau#sqrt{(#gamma_{D^{+}}^{2}-1)};Momentum (GeV/c); Decay length (mm)");
    Int_t ci;   // for color index setting
    ci = TColor::GetColor("#000099");
    for (int p=1; p<=24; p++)
    {
        double E = p*p + mass*mass;
        E= std:: sqrt(E);

        double gamma = E/mass;
        double decaylength = c* std::sqrt(gamma*gamma -1)*mt*1000;
        std:: cout<<" Decay length=" <<decaylength<<std::endl;
        h->SetLineColor(ci);
        h->GetXaxis()->CenterTitle(true);
        h->GetYaxis()->CenterTitle(true);
        h->SetMarkerColor(2);
        h->SetMarkerStyle(20);
        h->SetBinContent(p,decaylength);
        h->Draw("E2-text");
    }
// -----  Finish   -------------------------------------------------------
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished succesfully." << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
    cout << endl;
    // ------------------------------------------------------------------------

}

TH1D * DrawHistogram(TList *list, char *name,char *option, Int_t rebin, Int_t color, Int_t marker, Float_t scale, Bool_t draw = kTRUE){

  TH1D *histo = list->FindObject(name);
  histo->SetMarkerColor(color);
  histo->SetMarkerStyle(marker);
  histo->SetLineColor(color);
  histo->Rebin(rebin);
  histo->Scale(1.0/((Float_t)rebin));
  histo->Scale(scale);
  histo->SetMarkerSize(1.2);
  if(draw) histo->Draw(option);
  return histo;
}

TH1D *
GetITSsaSpectrum(TFile *file, Int_t part, Int_t charge, Int_t cent, Bool_t cutSpectrum = kTRUE, Bool_t addSystematicError = kTRUE)
{
  /* pt limits for combined spectra */
  Double_t ptMin[AliPID::kSPECIES] = {0., 0., 0.1, 0.2, 0.3};
  Double_t ptMax[AliPID::kSPECIES] = {0., 0., 0.6, 0.5, 0.6};

  TList *list = (TList *)file->Get("output");
  TH1D *hin = (TH1D *)list->FindObject(Form("h_%s_%s_cen_%d", ITSsaPartName[part], ITSsaChargeName[charge], cent));
  if (!hin) return NULL;  

  /* get systematics */
  TFile *fsys = TFile::Open("SPECTRASYS_ITSsa.root");
  TH1 *hsys = fsys->Get(Form("hSystTot%s%s", ITSsaChargeName[charge], ITSsaPartName[part]));
			
  TH1D *h = new TH1D(Form("hITSsa_cent%d_%s_%s", cent, AliPID::ParticleName(part), chargeName[charge]), "ITSsa", NptBins, ptBin);
  Double_t pt, width, value, error, sys;
  Int_t bin;
  for (Int_t ipt = 0; ipt < NptBins; ipt++) {
    /* get input bin */
    pt = h->GetBinCenter(ipt + 1);
    width = h->GetBinWidth(ipt + 1);
    bin = hin->FindBin(pt);
    /* sanity check */
    if (TMath::Abs(hin->GetBinCenter(bin) - pt) > 0.001 ||
	TMath::Abs(hin->GetBinWidth(bin) - width) > 0.001)
      continue;
    /* check pt limits */
    if (cutSpectrum && (pt < ptMin[part] || pt > ptMax[part])) continue;
    /* copy bin */
    value = hin->GetBinContent(bin);
    error = hin->GetBinError(bin);
    /*** TEMP ADD SYS ***/
    if (addSystematicError) {
      sys = hsys->GetBinContent(bin) * value;
      error = TMath::Sqrt(error * error + sys * sys);
    }
    h->SetBinContent(ipt + 1, value);
    h->SetBinError(ipt + 1, error);
  }

  h->SetTitle("ITSsa");
  h->SetLineWidth(1);
  h->SetLineColor(1);
  h->SetMarkerStyle(20);
  h->SetMarkerColor(1);
  h->SetFillStyle(0);
  h->SetFillColor(0);

  return h;
}

TH1D* DrawTH1(int isOld, TString var, int pthat, int hibinmin, int hibinmax, int color, int marker){
	if(pthat!=0){
		if(isOld==0)
			TFile *f = TFile::Open(Form("/cms/store/user/ymao/pA5TEV/Mixing/STARTHI53V27/merged/QCD_Pt_%d_TuneZ2_Mix_5p02TeV_pythia.root",pthat));
		else if(isOld==1)
			TFile *f = TFile::Open(Form("/cms/store/user/qixu/jetRpA/OldForest/pPb_Full_BTagForest%d_Fix3output.root",pthat));
		TTree *t = (TTree*)f->Get("hiEvtAnalyzer/HiTree");
	}
	else{
		if(isOld==0)
			TFile *f = TFile::Open(Form("/store/user/qixu/jetRpA/skimTree/MC%sakPu3PFskimfile0_10.root",coll.Data()));
		else if(isOld==1)
			TFile *f = TFile::Open(Form("/store/user/qixu/jetRpA/skimTree/MC%sakPu3PFskimfile0_8.root",coll.Data()));
		else if(isOld==2)
			TFile *f = TFile::Open(Form("/store/user/qixu/jetRpA/skimTree/DATA%sakPu3PFskimfile.root",coll.Data()));
		TTree *t = (TTree*)f->Get("nt");
		
	}
	Float_t hiHFplusEta4, hiHFminusEta4, vz;
	Int_t hiBin;
	TH1D* h = new TH1D("h","h",100,0,100);
	h->Sumw2();
	t->SetBranchAddress("hiHFplusEta4",&hiHFplusEta4);
	t->SetBranchAddress("hiHFminusEta4",&hiHFminusEta4);
	t->SetBranchAddress("vz",&vz);
	t->SetBranchAddress("hiBin",&hiBin);
	cout<<"Total Events="<<t->GetEntries()<<endl;
	for(int Ev=0;Ev<t->GetEntries();Ev++){
		if(Ev%1000000==0) cout<<"has processed "<<Ev<<endl;
		t->GetEntry(Ev);
		if(hiBin > hibinmax || hiBin < hibinmin) continue;
		if(TMath::Abs(vz)>=15) continue;
		if(var == "hiBin")
		h->Fill(hiBin);
		else if(var == "hiHFplusEta4")
		h->Fill(hiHFplusEta4);
		else if(var == "hiHFEta4")
		h->Fill(hiHFplusEta4+hiHFminusEta4);
	}
//	t->Draw(Form("%s>>h",var.Data()),"TMath::Abs(vz)<15");
//	h = (TH1D*)gDirectory->Get("h");
//	h->Scale(1.0/h->Integral());
	h->SetMarkerColor(color);
	h->SetMarkerStyle(marker);
	h->SetLineColor(color);
	h->SetMarkerSize(1.3);
	return h;
}

void plotter::draw_bias(TH1* output_, TH1D* truth_, TH1* bias_, TString file_name){
  TH1D* bias = (TH1D*) bias_->Clone("bias");
  TH1D* truth = (TH1D*) truth_->Clone("truth");
  TH1D* output = (TH1D*) output_->Clone("output");

  double ymax_temp = 0;
  if(truth->GetMaximum() > ymax_temp) ymax_temp = truth->GetMaximum();
  if(bias->GetMaximum() > ymax_temp) ymax_temp = bias->GetMaximum();
  if(output->GetMaximum() > ymax_temp) ymax_temp = output->GetMaximum();
  double ymax = 1.5 * ymax_temp;

  TCanvas *c= new TCanvas("Bias Distribution","",600,600);
  TGaxis::SetMaxDigits(3);
  gPad->SetLeftMargin(0.15);
  truth->SetTitle(file_name);
  truth->GetYaxis()->SetRangeUser(0., ymax);
  truth->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  truth->GetYaxis()->SetTitle("events");
  truth->GetYaxis()->SetTitleOffset(1.1);
  truth->GetXaxis()->SetTitleOffset(0.9);
  truth->GetYaxis()->SetTitleSize(0.05);
  truth->GetXaxis()->SetTitleSize(0.05);
  truth->GetYaxis()->SetNdivisions(505);

  truth->SetLineWidth(4);
  truth->SetLineColor(kRed);
  bias->SetLineWidth(4);
  bias->SetLineColor(kAzure+7);
  output->SetLineColor(kBlack);
  output->SetMarkerColor(kBlack);
  output->SetMarkerStyle(8);
  output->SetMarkerSize(1);

  truth->Draw("HIST SAME");
  bias->Draw("HIST SAME");
  output->Draw("E1 SAME");
  TLegend *l=new TLegend(0.55,0.6,0.85,0.85);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  l->AddEntry(output,"unfolded result","pl");
  l->AddEntry(truth,"truth","pl");
  l->AddEntry(bias,"bias distribution","pl");
  l->SetTextSize(0.04);
  l->Draw();
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}