C++ TH2D::Add方法代码示例

void AnaCombineEvH()
{
  TFile *file = TFile::Open("AnalysisResults_EmcalEvH.root", "READ");
  TList *list = (TList*)file->Get("listEmcalEvH");
  file->Close();
//=============================================================================

  TH2D *hist2 = 0;
  TH2D *hCentMtCh = 0;
  TH2D *hMtChMtEm = 0;
  for (Int_t i=0; i<4; i++) {
    hist2 = (TH2D*)list->FindObject(Form("hCentMtCh_%d",i));
    if (i==0) { hCentMtCh = (TH2D*)hist2->Clone("hCentMtCh"); hCentMtCh->Reset(); } hCentMtCh->Add(hist2); hist2 = 0;

    hist2 = (TH2D*)list->FindObject(Form("hMtChMtEm_%d",i));
    if (i==0) { hMtChMtEm = (TH2D*)hist2->Clone("hMtChMtEm"); hMtChMtEm->Reset(); } hMtChMtEm->Add(hist2); hist2 = 0;
  }
//=============================================================================

/*file = TFile::Open("AnalysisResults_EmcalEvH.root", "UPDATE");
  hCentMtCh->Write();
  hMtChMtEm->Write();
  file->Close();*/
//=============================================================================

  return;
}

void drawDist(const char* infilename, const char* system, Int_t rWrite, Int_t rPerformance) {

        myOptions(0);

        gROOT->ForceStyle();
        gStyle->SetPalette(1.0);

        TDatime now;
        int iDate = now.GetDate();
        int iYear=iDate/10000;
        int iMonth=(iDate%10000)/100;
        int iDay=iDate%100;
        char* cMonth[12]={"Jan","Feb","Mar","Apr","May","Jun",
                          "Jul","Aug","Sep","Oct","Nov","Dec"};
        char cStamp1[25],cStamp2[25];
        sprintf(cStamp1,"%i %s %i",iDay, cMonth[iMonth-1], iYear);
        sprintf(cStamp2,"%i/%.2d/%i",iDay, iMonth, iYear);

        TFile *f = new TFile(infilename, "read");
        // TList *list = (TList*)f->Get("femtolist");

        // pseudorapidity vs pt
        TH2D* ypt =(TH2D*)f->Get(Form("EtaPtcutPass1%stpcM%i",system,0));

        int minMultBin = 0;
        int maxMultBin = 6;
        double EvMultall = 0;

        for(int i = minMultBin; i < maxMultBin; i++) {
                TH1D* yptN =(TH1D*)f->Get(Form("EtaPtcutPass1%stpcM%i",system,i));
                ypt->Add(yptN);
                //delete hEvMult;
        }

        TCanvas *c2 = new TCanvas("pseudorapidity vs pt", "pseudorapidity vs pt");
        c2->SetGridx();
        c2->SetGridy();
        c2->SetFillColor(10);
        ypt->GetXaxis()->SetTitle("#eta");
        ypt->GetYaxis()->SetTitle("p_{T}");
        ypt->GetXaxis()->SetTitleOffset(1.3);
        ypt->GetYaxis()->SetTitleOffset(1.3);
        ypt->GetXaxis()->SetRangeUser(-0.8,0.8);
        ypt->GetYaxis()->SetRangeUser(0.1,8.);
        ypt->Draw("colz");

        // https://wiki.bnl.gov/eic/index.php/ROOT#Moving_and_resizing_the_palette_axis_of_a_2D_histogram
        gPad->SetRightMargin( 0.12 ); // The default right margin is 0.1 i.e. 10% of the image width
//   TPaletteAxis* palette
//       = dynamic_cast<TPaletteAxis*>( myHistogram.GetListOfFunctions()->FindObject( "palette" ) );
//   if( palette ) {
//     palette->SetX1NDC( 0.86 ); // Start the palette 86 % of the way across the image
//     palette->SetX1NDC( 0.91 ); // End the palette 91% of the way across the image
//     gPad->Modified(); // Update with the new position
//   } // if

        postprocess(c2,Form("ypt%s",system),rWrite,rPerformance);

}

TH2D* combinePlots(std::vector<TH2D*> plots, std::vector<double> scalingFactors) {
	TH1::SetDefaultSumw2();

	TH2D* h = (TH2D*)plots[0]->Clone(plots[0]->GetName());
	h->Scale(scalingFactors[0]);
	for (unsigned i = 1; i < plots.size(); i++) {
		h->Add(plots[i], scalingFactors[i]);
	}
	return h;
}

/* Get a fake estimate from the given histogram, i.e. return Data - MC.
 *
 * Subtract all background MCs from data but the one specified with
 * "notremove". Under- and overflows are ignored in computation.
 *
 * @param hname     Histogram name from which to compute estimate
 * @return          Histogram of fakes
 */
TH2D * get_fakes_2d(const char * hname)
{
  // get number of single fakes from data histogram
  plot2(hname);
  TH2D * hData = dataHisto2();
  if (hData == 0) {
    THROW("get_fakes() needs a data histogram");
  }
  // data
  double N = hData->Integral();
  INFO("Data events: " << N);

  TH2D * hBack = 0;
  for (unsigned int i = 0; i < sizeof(removeNames)/sizeof(void *); i++) {
    TH2D * hSub = backgroundHisto2(removeNames[i]);
    if (hSub == 0) {
      THROW(string("get_fakes() problem getting background ")+removeNames[i]);
    }
    // add backgrounds together
    if (hBack == 0) {
      hBack = hSub;
    }
    else {
      hBack->Add(hSub);
      delete hSub;
    }
  }
  N = hBack->Integral();
  INFO("Background events : " << N);

  // subtract
  hData->Add(hBack, -1.);
  // temporarily needed for function call
  delete hBack;
  return hData;
}

TH2D* combineScale(std::vector<TH2D*> hists, std::vector<double> scalingFactors) {
    TH1::SetDefaultSumw2();

    // get total
    double total = 0;
    for (int i = 0; i < scalingFactors.size(); i++) {
        total += scalingFactors[i];
    }
    TH2D* h = (TH2D*) hists[0]->Clone(hists[0]->GetName());
    h->Scale(scalingFactors[0]/total);

    for (i = 1; i < scalingFactors.size(); i++) {
        h->Add(hists[i], scalingFactors[i]/total);
    }

    return h;
}

TH2D* Plot2D (std::string var,std::vector<TFile*>& tfiles,std::vector<double>& weights) {

    double weightTot = 0.0;
    for ( int i = 0; i<weights.size(); i++) {
        weightTot+=weights[i];
    }

    std::string name = var;
    name+="NEW";
    TH2D* hVar = (TH2D*)tfiles[0]->Get(var.c_str());
    TH2D* HVar = (TH2D*)hVar->Clone(name.c_str());
    //HVar->Sumw2();

    for (int i=1; i<tfiles.size(); i++) {
        TH2D* htempVar = (TH2D*)tfiles[i]->Get(var.c_str());
        HVar->Add(htempVar,weights[i]/weightTot);
    }

    return HVar;

}

TH2D* Plots2D::readHistogram(Sample sample, Variable variable) {

	cout << "plot: " << "Binning/"+folder+variable.name << endl;
	TH2D* plot = (TH2D*) sample.file->Get("Binning/"+folder+variable.name+"_2btags");
	TH2D* plot2 = (TH2D*) sample.file->Get("Binning/"+folder+variable.name+"_3btags");
	TH2D* plot3 = (TH2D*) sample.file->Get("Binning/"+folder+variable.name+"_4orMoreBtags");

//	TH1D* plot = (TH1D*) sample.file->Get(selection+objName+"/"+variable.name+"_2btags");
//	TH1D* plot2 = (TH1D*) sample.file->Get(selection+objName+"/"+variable.name+"_3btags");
//	TH1D* plot3 = (TH1D*) sample.file->Get(selection+objName+"/"+variable.name+"_4orMoreBtags");

	plot->Add(plot2);
	plot->Add(plot3);

	plot->SetFillColor(sample.fillColor);
	plot->SetLineColor(sample.lineColor);

	plot->Rebin2D(variable.rebinFact, variable.rebinFact);

	return plot;
}

void TSVDUnfoldExample()
{
   gROOT->Reset();
   gROOT->SetStyle("Plain");
   gStyle->SetOptStat(0);

   TRandom3 R;

   const Double_t cutdummy= -99999.0;

   // --- Data/MC toy generation -----------------------------------

   // The MC input
   Int_t nbins = 40;
   TH1D *xini = new TH1D("xini", "MC truth", nbins, -10.0, 10.0);
   TH1D *bini = new TH1D("bini", "MC reco", nbins, -10.0, 10.0);
   TH2D *Adet = new TH2D("Adet", "detector response", nbins, -10.0, 10.0, nbins, -10.0, 10.0);

   // Data
   TH1D *data = new TH1D("data", "data", nbins, -10.0, 10.0);
   // Data "truth" distribution to test the unfolding
   TH1D *datatrue = new TH1D("datatrue", "data truth", nbins, -10.0, 10.0);
   // Statistical covariance matrix
   TH2D *statcov = new TH2D("statcov", "covariance matrix", nbins, -10.0, 10.0, nbins, -10.0, 10.0);

   // Fill the MC using a Breit-Wigner, mean 0.3 and width 2.5.
   for (Int_t i= 0; i<100000; i++) {
      Double_t xt = R.BreitWigner(0.3, 2.5);
      xini->Fill(xt);
      Double_t x = Reconstruct( xt, R );
      if (x != cutdummy) {
         Adet->Fill(x, xt);
         bini->Fill(x);
      }
   }

   // Fill the "data" with a Gaussian, mean 0 and width 2.
   for (Int_t i=0; i<10000; i++) {
      Double_t xt = R.Gaus(0.0, 2.0);
      datatrue->Fill(xt);
      Double_t x = Reconstruct( xt, R );
      if (x != cutdummy)
      data->Fill(x);
   }

   cout << "Created toy distributions and errors for: " << endl;
   cout << "... \"true MC\"   and \"reconstructed (smeared) MC\"" << endl;
   cout << "... \"true data\" and \"reconstructed (smeared) data\"" << endl;
   cout << "... the \"detector response matrix\"" << endl;

   // Fill the data covariance matrix
   for (int i=1; i<=data->GetNbinsX(); i++) {
       statcov->SetBinContent(i,i,data->GetBinError(i)*data->GetBinError(i));
   }

   // --- Here starts the actual unfolding -------------------------

   // Create TSVDUnfold object and initialise
   TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );

   // It is possible to normalise unfolded spectrum to unit area
   tsvdunf->SetNormalize( kFALSE ); // no normalisation here

   // Perform the unfolding with regularisation parameter kreg = 13
   // - the larger kreg, the finer grained the unfolding, but the more fluctuations occur
   // - the smaller kreg, the stronger is the regularisation and the bias
   TH1D* unfres = tsvdunf->Unfold( 13 );

   // Get the distribution of the d to cross check the regularization
   // - choose kreg to be the point where |d_i| stop being statistically significantly >>1
   TH1D* ddist = tsvdunf->GetD();

   // Get the distribution of the singular values
   TH1D* svdist = tsvdunf->GetSV();

   // Compute the error matrix for the unfolded spectrum using toy MC
   // using the measured covariance matrix as input to generate the toys
   // 100 toys should usually be enough
   // The same method can be used for different covariance matrices separately.
   TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );

   // Now compute the error matrix on the unfolded distribution originating
   // from the finite detector matrix statistics
   TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );

   // Sum up the two (they are uncorrelated)
   ustatcov->Add( uadetcov );

   //Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics.
   TH2D* utaucov = tsvdunf->GetXtau();
   utaucov->Add( uadetcov );

   //Get the computed inverse of the covariance matrix
   TH2D* uinvcov = tsvdunf->GetXinv();


   // --- Only plotting stuff below ------------------------------

   for (int i=1; i<=unfres->GetNbinsX(); i++) {
      unfres->SetBinError(i, TMath::Sqrt(utaucov->GetBinContent(i,i)));
//.........这里部分代码省略.........

void pc3matching() {

  TFile *f = TFile::Open("output_perform.root");

  ofstream fout("run16dAupc3matching.h");
  ofstream fout2("run16dAupc3matchingfirst.h");
  fout << "float pc3dphimean[2][2][5][50];" << endl;
  fout << "float pc3dphisigma[2][2][5][50];" << endl;
  fout << "float pc3dzmean[2][2][5][50];" << endl;
  fout << "float pc3dzsigma[2][2][5][50];" << endl;
  fout << "float pc3dphimeanerr[2][2][5][50];" << endl;
  fout << "float pc3dphisigmaerr[2][2][5][50];" << endl;
  fout << "float pc3dzmeanerr[2][2][5][50];" << endl;
  fout << "float pc3dzsigmaerr[2][2][5][50];" << endl;
  fout << " " << endl;
  fout << "void fetchpc3dphidz();" << endl;
  fout << " " << endl;
  fout << "void fetchpc3dphidz() {" << endl;


  float max = 0.0;
  float sigma = 0.0;
  float mean = 0.0;

  for (Int_t iarm = 0; iarm < 2; iarm++) {
    for (Int_t ich = 0; ich < 2; ich++) {
      for (Int_t ipt = 0; ipt < 50; ipt++) {
        for(Int_t ivz = 0; ivz < 5; ivz++) {
           // if(ipt!=3 || iarm!=0 || ich!=0 || ivz!=0)continue;
            cout<<iarm<<" "<<ich<<" "<<ipt<<" "<<ivz<<endl;
          double sigmaerr=0.0;
	  double meanerr=0.0;
          TString ch = "";
            if(ich == 0) ch = "pos";
            if(ich == 1) ch = "neg";
	  TString histname = Form("pc3dphidz_arm%d_%s_z%d_%d",iarm,ch.Data(),ivz*2,ipt);
	  TString histname1 = Form("pc3dphidz_arm%d_%s_z%d_%d",iarm,ch.Data(),ivz*2+1,ipt);
          TH2D *hist = (TH2D*) f->Get(histname);
          TH2D *hist1 = (TH2D*) f->Get(histname1);
          hist->Add(hist1);

          TH1D *dphi = (TH1D*) hist->ProjectionX(Form("pc3dphi_%d_%d_%d_%d",iarm,ich,ipt,ivz));
          TH1D *dz = (TH1D*) hist->ProjectionY(Form("pc3dz_%d_%d_%d_%d",iarm,ich,ipt,ivz));

	  dphi->GetXaxis()->SetRangeUser(-0.1,0.1);

          gStyle->SetOptFit(1101);

          TF1 *fphi1 = new TF1("fphi1","gaus",-0.1,0.1);
          TF1 *fz1 = new TF1("fz1","gaus",-10,10);
          TF1 *fphi2 = new TF1("fphi2","gaus(0)+gaus(3)",-0.1,0.1);
          TF1 *fz2 = new TF1("fz2","gaus(0)+gaus(3)",-10,10);
          TF1 *phi_gaus1 = new TF1("phi_gaus1","gaus",-0.1,0.1);
          TF1 *phi_gaus2 = new TF1("phi_gaus2","gaus",-0.1,0.1);
          TF1 *z_gaus1 = new TF1("z_gaus1","gaus",-10,10);
          TF1 *z_gaus2 = new TF1("z_gaus2","gaus",-10,10);


          Float_t Xbins = dphi->GetNbinsX();
          Float_t Xmin = dphi->GetXaxis()->GetXmin();
          Float_t Xmax = dphi->GetXaxis()->GetXmax();
          mean = (dphi->GetMaximumBin() * (Xmax-Xmin))/Xbins + Xmin;
          max = dphi->GetMaximum();
          fphi1->SetRange(mean-0.01,mean+0.01);
	  fphi1->SetParameters(max,mean);

          TCanvas *c = new TCanvas("c","c",500,500);
	  //dphi->Scale(1./dphi->Integral());
	  dphi->SetTitle("dphi matching");
	  dphi->GetXaxis()->SetTitle("dphi");
	  dphi->GetYaxis()->SetTitle("# tracks");
	  if (iarm == 1 && ich ==0) dphi->Rebin(2);
	//  dphi->Rebin(2);
	  dphi->SetMarkerSize(1);
	  dphi->SetMarkerStyle(kFullCircle);
          dphi->Draw("P");
          dphi->Fit("fphi1","RQ0");
          double dphi_par[6];
          fphi1->GetParameters(dphi_par);
          dphi_par[3] = 0.1*dphi_par[0];
          dphi_par[4] = dphi_par[1];
          dphi_par[5] = 8*dphi_par[2];
          fphi2->SetParameters(dphi_par);
	  fphi2->SetParLimits(3,0,5*dphi_par[3]);
	  fphi2->SetParLimits(4,-1,1);
	  fphi2->SetParLimits(5,0,100*dphi_par[2]);
          dphi->Fit("fphi2","RQ0");
	  fphi2->Draw("same");
          fphi2->GetParameters(dphi_par);
	  meanerr = fphi2->GetParError(1);
	  sigmaerr = fphi2->GetParError(2);

          fout << "pc3dphimean[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "]=" << dphi_par[1] << ";" << endl;
          fout << "pc3dphisigma[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "]=" <<  dphi_par[2] << ";" << endl;
          fout << "pc3dphimeanerr[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "]=" << meanerr << ";" << endl;
          fout << "pc3dphisigmaerr[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "]=" <<  sigmaerr << ";" << endl;
          if(ipt == 2){
            fout2 << "PC3_dphifirst_mean[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "] = " << dphi_par[1] << ";" << endl;
            fout2 << "PC3_dphifirst_sigma[" << iarm << "][" << ich << "][" << ivz << "][" << ipt << "] = " <<  dphi_par[2] << ";" << endl;
          }
//.........这里部分代码省略.........

// infilename - root file with relevant histograms
// system - PP,APAP,PP
// status - Pass,Fail
// rWrite - 0-no,1-png,2-eps
// rPerformance - 0-no,1-yes (ALICE logo etc.)
// bin: 0 - all, 1- 0:5, 2- 5:10, etc
void drawPID(const char* infilename, const char* system, const char* status, Int_t rWrite, Int_t rPerformance, Int_t bin)
{


    TFile *f = new TFile(infilename, "read");

    // TPC dEdx
    TH2D* TPCdEdx =(TH2D*)f->Get(Form("TPCdEdxcut%s1%stpcM%i",status, system,0));

    if (!bin) {
        int minMultBin = 0;
        int maxMultBin = 6; // 8
    }
    else {
        int minMultBin = bin-1;
        int maxMultBin = bin; // 8
    }

    double EvMultall = 0;

    for (int i = minMultBin; i < maxMultBin; i++) {

        TH2D* TPCdEdxN =(TH2D*)f->Get(Form("TPCdEdxcut%s1%stpcM%i",status,system,i));
        TPCdEdx->Add(TPCdEdxN);

        cout << i << " " << TPCdEdxN->GetEntries() << endl;

        //delete hEvMult;
    }

    TCanvas *c2 = new TCanvas("TPC dEdx", "TPC dEdx");
    c2->SetGridx();
    c2->SetGridy();
    c2->SetFillColor(10);
    c2->SetRightMargin(1.9);
    c2->SetLogz();

    TPCdEdx->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
    TPCdEdx->GetXaxis()->SetRangeUser(0.0,6.0);
    TPCdEdx->GetYaxis()->SetTitle("dE/dx");
    TPCdEdx->GetZaxis()->SetLabelSize(0.03);
    TPCdEdx->Draw("colz");

// 	double a1 = -3000.0; double b1 =  1280.0;
// 	double a2 = -312.5;  double b2 =  312.5;
// 	double a3 = -200.0;  double b3 =  240.0;

//    TF1 *fa1 = new TF1("fa1","-1800*x+940",0.3,0.4);
//   fa1->Draw("same");
//   TF1 *fa2 = new TF1("fa2","-500.0*x+420.0",0.4,0.6);
//   fa2->Draw("same");
// TF1 *fa3 = new TF1("fa3","-216.7*x+250.0",0.6,0.9);
//    fa3->Draw("same");
// TF1 *fa4 = new TF1("fa4","-566.7*x+570.0",0.6,0.75);
//    fa4->Draw("same");
// TF1 *fa5 = new TF1("fa5","-2076.92*x+1476.15",0.47,0.6);
//    fa5->Draw("same");
//   cout<<TPCdEdx->GetNbinsX()<<endl;
//   cout<<TPCdEdx->GetNbinsY()<<endl;
//   for (int ii=0;ii<TPCdEdx->GetNbinsX();ii++){

//     for (int jj=0;jj<TPCdEdx->GetNbinsY();jj++){

//       cout<<"binX: "<<ii<<endl;
//       cout<<"binY: "<<jj<<endl;
//       cout<<"val: "<<TPCdEdx->GetBinContent(ii,jj)<<endl;
//     }
//   }

// 	TH1D *py = TPCdEdx->ProjectionY("py", 230, 232); // where firstYbin = 0 and lastYbin = 9
// 	TCanvas *c22 = new TCanvas("TPC2", "TPC2");
// 	py->Draw();

    postprocess(c2,Form("TPCdEdx%s",status),rWrite,rPerformance,system);

    // TPC Nsigma
    TH2D* TPCNsigma =(TH2D*)f->Get(Form("TPCNSigmacut%s1%stpcM%i",status,system,0));

    // int minMultBin = 0;
    // int maxMultBin = 2;
    double EvMultall = 0;

    for(int i = minMultBin; i<maxMultBin; i++) {
        //all
        TH2D* TPCNsigmaN =(TH2D*)f->Get(Form("TPCNSigmacut%s1%stpcM%i",status,system,i));
        TPCNsigma->Add(TPCNsigmaN);
        //delete hEvMult;
    }

    TCanvas *c3 = new TCanvas("TPC Nsigma", "TPC Nsigma");
    c3->SetGridx();
    c3->SetGridy();
    c3->SetFillColor(10);
    c3->SetRightMargin(1.7);
//.........这里部分代码省略.........

/*
 * this script takes 2 TStrings as root filenames as a parameters
 * basic functionality:
 * loop through all directories (the mass bins) in the root file
 * -> create difference plots
 * -> create global plots
 * -> create 2D diff vs mass plots
 * -> etc...
 */
void plotGlobalWeightedEvts_Kpipi(TString input_filename, TString output_filename) {
  setupBookies();

  gROOT->SetStyle("Plain");
  gStyle->SetTitleFont(10*13+2,"xyz");
  gStyle->SetTitleSize(0.06, "xyz");
  gStyle->SetTitleOffset(1.3,"y");
  gStyle->SetTitleOffset(1.3,"z");
  gStyle->SetLabelFont(10*13+2,"xyz");
  gStyle->SetLabelSize(0.06,"xyz");
  gStyle->SetLabelOffset(0.009,"xyz");
  gStyle->SetPadBottomMargin(0.16);
  gStyle->SetPadTopMargin(0.16);
  gStyle->SetPadLeftMargin(0.16);
  gStyle->SetPadRightMargin(0.16);
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gROOT->ForceStyle();
  gStyle->SetFrameFillColor(0);
  gStyle->SetFrameFillStyle(0);
  TGaxis::SetMaxDigits(3);
  //IsPhDStyle = true;

  int massbins =0;
  double mass= 0.0, massstart =1000.0, massend=0.0;
  std::map<std::string, std::pair<double, std::pair<double, double> > > diffbounds;

  TFile* infile = TFile::Open(input_filename, "READ");
  TFile* outfile = new TFile(output_filename, "RECREATE");
  outfile->mkdir("global");

  TList *dirlist = infile->GetListOfKeys();
  massbins = dirlist->GetSize();
  infile->cd();
  TIter diriter(dirlist);
  TDirectory *dir;

  std::cout<< "scanning directories and creating overview canvases..." <<std::endl;
  while ((dir = (TDirectory *)diriter())) {
    std::string dirname = dir->GetName();
    // check if directory is mass bin dir
    unsigned int pointpos = dirname.find(".");
    if(pointpos == 0 || pointpos == dirname.size()) continue;
    std::string masslow = dirname.substr(0, pointpos+1);
    std::string masshigh = dirname.substr(pointpos+1);
    double massstarttemp = atof(masslow.c_str())/1000;
    double massendtemp = atof(masshigh.c_str())/1000;
    if((int)(massendtemp - massstarttemp) != massbinwidth)
      massbinwidth = (int)(massendtemp - massstarttemp);
    mass = (massstarttemp + massendtemp)/2;
    if(massstart > massstarttemp) massstart = massstarttemp;
    if(massend < massendtemp) massend = massendtemp;

    outfile->cd();
    outfile->mkdir(dir->GetName());
    infile->cd(dir->GetName());

    // make list of MC Histograms
    TList mclist;
    TList *histlist = gDirectory->GetListOfKeys();
    TIter histiter(histlist);
    TObject *obj;
    while ((obj = histiter())) {
      TString s(obj->GetName());
      if(s.EndsWith("MC"))
        mclist.Add(obj);
      else if(s.Contains("MC_"))
        mclist.Add(obj);
    }
    make1DOverviewCanvas(infile, outfile, &mclist, dirname);
    histiter = TIter(&mclist);
    TH1D *diffhist, *mchist;
    while ((mchist = (TH1D*)histiter())) {
      // generate difference histograms
      std::string hnamemc(mchist->GetName());
      // create new string with MC exchanged for Diff
      std::string hnamediff(hnamemc);
      int pos = hnamemc.find("MC");
      hnamediff.erase(pos, 2);
      hnamediff.insert(pos, "Diff");

      infile->GetObject((std::string(dir->GetName())+"/"+hnamediff).c_str(), diffhist);

      if (diffhist) {
        // get diff min max values
        std::pair<double, std::pair<double, double> > p;

        bool change =false;
        double maxdiff = diffhist->GetMaximum();
        double maxdifftemp = diffhist->GetMinimum();
        if(abs(maxdifftemp) > maxdiff) maxdiff = maxdifftemp;
//.........这里部分代码省略.........

void sysError(
    TString inFileName="jfh_HCPR_J50U_Cent30to100_Aj0to100_SubEtaRefl.root",
    Int_t compMode = 0, // Compare mode: 0 reco-genSig, 1 reco-genAll, 3 genAll-genSig, 4 calo_genp-allGen
    Int_t sysMode = 0, // Plot mode: 0 for simple plot, 1 for difference
    TString outdir = ".",
    TString title = "test"
    ) {
  // ===============================================
  // Inputs
  // ===============================================
  TFile *f = new TFile(inFileName);
  TString inFileNameStrip(inFileName); inFileNameStrip.ReplaceAll(".root","");
  TString inFileNameGen(inFileName);
  if (compMode==0) { 
    inFileNameGen.ReplaceAll("djcalo","djcalo_genp");
    inFileNameGen.ReplaceAll("HydjetAll","HydjetSig");
  }
  if (compMode==1) inFileNameGen.ReplaceAll("djcalo","djcalo_genp");
  else if (compMode==3) inFileNameGen.ReplaceAll("HydjetAll","HydjetSig");
  else if (compMode==5) inFileNameGen.ReplaceAll("djcalo_genp","djgen");
  TFile *fgen = new TFile(inFileNameGen);
  TString inFileNameStripGen(inFileNameGen); inFileNameStripGen.ReplaceAll(".root","");
  cout << "==========================================================" << endl;
  cout << "Compare: " << inFileName << endl
       << "         vs " << endl
       << "         " << inFileNameGen << endl;
  cout << "==========================================================" << endl;

  // ===============================================
  // Setup
  // ===============================================
  TString tag=Form("sysError_%s_%s_%d_%d",inFileNameStrip.Data(),title.Data(),compMode,sysMode);

  // ===============================================
  // Analyze
  // ===============================================
  TH2D * hPtPNDR = (TH2D*) f->Get("hPtPNDR");
  TH2D * hPtPADR = (TH2D*) f->Get("hPtPADR");
  TH2D * hPtPNDRBg = (TH2D*) f->Get("hPtPNDRBg");
  TH2D * hPtPADRBg = (TH2D*) f->Get("hPtPADRBg");
  TH2D * hPtPNDRSub = (TH2D*)hPtPNDR->Clone(tag+"hPtPNDRSub");
  TH2D * hPtPADRSub = (TH2D*)hPtPADR->Clone(tag+"hPtPADRSub");
  hPtPNDRSub->Add(hPtPNDR,hPtPNDRBg,1,-1);
  hPtPADRSub->Add(hPtPADR,hPtPADRBg,1,-1);

  TH2D * hPtPNDRGen = (TH2D*) fgen->Get("hPtPNDR");
  TH2D * hPtPADRGen = (TH2D*) fgen->Get("hPtPADR");
  TH2D * hPtPNDRBgGen = (TH2D*) fgen->Get("hPtPNDRBg");
  TH2D * hPtPADRBgGen = (TH2D*) fgen->Get("hPtPADRBg");
  TH2D * hPtPNDRSubGen = (TH2D*)hPtPNDRGen->Clone(tag+"hPtPNDRSub");
  TH2D * hPtPADRSubGen = (TH2D*)hPtPADRGen->Clone(tag+"hPtPADRSub");
  hPtPNDRSubGen->Add(hPtPNDRGen,hPtPNDRBgGen,1,-1);
  hPtPADRSubGen->Add(hPtPADRGen,hPtPADRBgGen,1,-1);

  // ===============================================
  // Draw
  // ===============================================
  // Get Pt info
  Int_t numPtBins=hPtPNDR->GetNbinsX();
  TH1D * hPt = (TH1D*)hPtPNDR->ProjectionX("hPt");
  /*
  cout << "Pt bins: " << numPtBins << endl;
  for (Int_t i=0; i<numPtBins+2; ++i) {
    cout << "Pt Bin " << i << " Low Edge: " << hPt->GetBinLowEdge(i) << endl;
  }
  */

  // What pt bins to draw
  const Int_t numPtBinsDraw=3;

  TCanvas * c6 = new TCanvas("c"+tag,"c"+tag,1400,500);
  c6->Divide(3,1);
  for (Int_t i=0; i<numPtBinsDraw; ++i) {
    Int_t iBeg,iEnd;
    if (i==0) { iBeg=2; iEnd=3;}
    if (i==1) { iBeg=4; iEnd=4;}
    if (i==2) { iBeg=5; iEnd=numPtBins;}
    cout << "Bin: " << iBeg <<  " to " << iEnd << endl;
    TH1D * hNr = (TH1D*)hPtPNDRSub->ProjectionY(tag+Form("hPNDRSub_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hAw = (TH1D*)hPtPADRSub->ProjectionY(tag+Form("hPADRSub_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hNrGen = (TH1D*)hPtPNDRSubGen->ProjectionY(tag+Form("hPNDRSubGen_%d_%d",iBeg,iEnd),iBeg,iEnd);
    TH1D * hAwGen = (TH1D*)hPtPADRSubGen->ProjectionY(tag+Form("hPADRSubGen_%d_%d",iBeg,iEnd),iBeg,iEnd);
    if (sysMode==1) {
      hNr->Add(hNrGen,-1);
      hAw->Add(hAwGen,-1);
    }
    if (sysMode==2) {
      hNr->Divide(hNrGen);
      hAw->Divide(hAwGen);
    }
    // Print
    cout << Form("%.1f < P_{T} < %.1f GeV: ",hPt->GetBinLowEdge(iBeg),hPt->GetBinLowEdge(iEnd+1))
      << " SigSubBkg Integral - Nr: " << hNr->Integral() << " Aw: " << hAw->Integral() << endl
      << " Gen - Nr: " << hNrGen->Integral() << " Aw: " << hAwGen->Integral() << endl;
    // Styles
    hNr->SetMarkerColor(kRed);
    hNr->SetLineColor(kRed);
    hAw->SetMarkerColor(kBlue);
    hAw->SetLineColor(kBlue);
    hAwGen->SetLineStyle(2);
//.........这里部分代码省略.........

// infilename - root file with relevant histograms
// system - PP,APAP,PP
// status - Pass,Fail
// rWrite - 0-no,1-png,2-eps
// rPerformance - 0-no,1-yes (ALICE logo etc.)
// bin: 0 - all, 1- 0:5, 2- 5:10, etc
void drawDCA(const char* infilename, const char* system, const char* status, Int_t rWrite, Int_t rPerformance, int isMC, Int_t bin, Int_t ptrange)
{

    myOptions(0);

    gROOT->ForceStyle();
    gStyle->SetPalette(1.0);

    TDatime now;
    int iDate = now.GetDate();
    int iYear=iDate/10000;
    int iMonth=(iDate%10000)/100;
    int iDay=iDate%100;
    char* cMonth[12]={"Jan","Feb","Mar","Apr","May","Jun",
                      "Jul","Aug","Sep","Oct","Nov","Dec"};
    char cStamp1[25],cStamp2[25];
    sprintf(cStamp1,"%i %s %i",iDay, cMonth[iMonth-1], iYear);
    sprintf(cStamp2,"%i/%.2d/%i",iDay, iMonth, iYear);

    TFile *f = new TFile(infilename, "read");

    // DCA xy
    TH2D* DCAxy =(TH2D*)f->Get(Form("DCARPtcut%s1%stpcM%d","Pass", system,0));


    if (!bin) {
        int minMultBin = 0;
        int maxMultBin = 6; // 8
    }
    else {
        int minMultBin = bin-1;
        int maxMultBin = bin; // 8
    }

    // int minMultBin = 0;
    // int maxMultBin = 0; // 8

    double EvMultall = 0;

    for(int i = minMultBin; i<maxMultBin; i++) {

        TH2D* DCAxyN = (TH2D*)f->Get(Form("DCARPtcut%s1%stpcM%d",status, system,i));
        DCAxy->Add(DCAxyN);

        cout<<i<<" "<<DCAxyN->GetEntries()<<endl;

        //delete hEvMult;
    }


    if (!isMC) {
        TCanvas *c2 = new TCanvas("DCA xy prim", "DCA xy prim");
        c2->SetGridx();
        c2->SetGridy();
        c2->SetFillColor(10);
        c2->SetRightMargin(1.9);
        c2->SetLogz();

        DCAxy->GetXaxis()->SetTitle("DCA_{XY} (cm)");
        DCAxy->GetXaxis()->SetRangeUser(-5.0,5.0);
        DCAxy->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
        // DCAxy->GetZaxis()->SetLabelSize(0.05);
        DCAxy->Draw("colz");
        postprocess(c2,Form("DCAxy%s",status),rWrite,rPerformance,system);

        // TCanvas *c4 = new TCanvas("DCA xy Projection X", "DCA xy Projection X");
        // c4->SetGridx();
        // c4->SetGridy();
        // c4->SetFillColor(10);
        // c4->SetRightMargin(1.9);
        // c4->SetLogy();

        gStyle->SetOptTitle(1);
        TCanvas *myCan = new TCanvas("myCan",cStamp1,600,400);
        myCan->Draw();
        myCan->cd();

        TPad *myPad = new TPad("myPad", "The pad",0,0,1,1);
        myPadSetUp(myPad,0.15,0.04,0.04,0.15);
        myPad->Draw();
        myPad->SetLogy();

        myPad->cd();

        if (ptrange == 0) {
            TH1D* pripp = (TH1D*)DCAxy->ProjectionX("zxc1",1,100);
            pripp->SetTitle("0.5 < #it{p}_{T} < 3 GeV/#it{c}");
        }
        else if (ptrange == 1) {
            TH1D* pripp = (TH1D*)DCAxy->ProjectionX("zxc1",15,33);
            pripp->SetTitle("0.5 < #it{p}_{T} < 1 GeV/#it{c}");
        }
        else if (ptrange == 2) {
            TH1D* pripp = (TH1D*)DCAxy->ProjectionX("zxc1",33,100);
//.........这里部分代码省略.........

void drawJetFragBalance_DRDiff(
			       TString inFileName=    "plot/jfhCorrEtaPtBin4RBin20v2_HCPR_J50U_djcalo_Cent0to30_Aj0to100_SubEtaRefl.root",
			       TString inFileNameHyPy="plot/jfhCorrEtaPtBin4RBin20v2_Hydjet_djcalo_Cent0to30_Aj0to100_SubEtaRefl.root",
    TString title = "test",
    Int_t drawMode=1,
			       Int_t doLeg=1,
bool cumulative = 1
    ) {
  TFile *f = new TFile(inFileName);
  TString inFileNameStrip(inFileName); inFileNameStrip.ReplaceAll(".root","");

  TH2D * hPtPNDR = (TH2D*) f->Get("hPtPNDR");
  TH2D * hPtPADR = (TH2D*) f->Get("hPtPADR");
  TH2D * hPtPNDRBg = (TH2D*) f->Get("hPtPNDRBg");
  TH2D * hPtPADRBg = (TH2D*) f->Get("hPtPADRBg");
  TH2D * hPtPNDRBgSub = (TH2D*)hPtPNDR->Clone(inFileNameStrip+"hPtPNDRBgSub");
  TH2D * hPtPADRBgSub = (TH2D*)hPtPADR->Clone(inFileNameStrip+"hPtPADRBgSub");
  hPtPNDRBgSub->Add(hPtPNDR,hPtPNDRBg,1,-1);
  hPtPADRBgSub->Add(hPtPADR,hPtPADRBg,1,-1);

  TFile *fhypy = new TFile(inFileNameHyPy);
  TH2D * hPtPNDRHyPy = (TH2D*) fhypy->Get("hPtPNDR");
  TH2D * hPtPADRHyPy = (TH2D*) fhypy->Get("hPtPADR");
  TH2D * hPtPNDRBgHyPy = (TH2D*) fhypy->Get("hPtPNDRBg");
  TH2D * hPtPADRBgHyPy = (TH2D*) fhypy->Get("hPtPADRBg");
  TH2D * hPtPNDRBgSubHyPy = (TH2D*)hPtPNDR->Clone(inFileNameStrip+"hPtPNDRBgSubHyPy");
  TH2D * hPtPADRBgSubHyPy = (TH2D*)hPtPADR->Clone(inFileNameStrip+"hPtPADRBgSubHyPy");
  hPtPNDRBgSubHyPy->Add(hPtPNDRHyPy,hPtPNDRBgHyPy,1,-1);
  hPtPADRBgSubHyPy->Add(hPtPADRHyPy,hPtPADRBgHyPy,1,-1);

  // Get Pt info
  Int_t numBinsPt=hPtPNDR->GetNbinsX();
  Int_t numBinsDR=hPtPNDR->GetNbinsY();
  TH1D * hPt = (TH1D*)hPtPNDR->ProjectionX("hPt");
  cout << "Pt bins: " << numBinsPt << endl;

  Double_t totPtBgSubNr=hPtPNDRBgSub->Integral();
  Double_t totPtBgSubAw=hPtPADRBgSub->Integral();

  Double_t totPtBgSubNrHyPy=hPtPNDRBgSubHyPy->Integral();
  Double_t totPtBgSubAwHyPy=hPtPADRBgSubHyPy->Integral();

  int ptUp = 2;
  TH1D * hDRBgSubNr = (TH1D*)hPtPNDRBgSub->ProjectionY(inFileNameStrip+"hDRBgSubNr",1,ptUp);
  TH1D * hDRBgSubAw = (TH1D*)hPtPADRBgSub->ProjectionY(inFileNameStrip+"hDRBgSubAw",1,ptUp);
  TH1D * hDRBgSubNrHyPy = (TH1D*)hPtPNDRBgSubHyPy->ProjectionY(inFileNameStrip+"hDRBgSubNrHyPy",1,ptUp);
  TH1D * hDRBgSubAwHyPy = (TH1D*)hPtPADRBgSubHyPy->ProjectionY(inFileNameStrip+"hDRBgSubAwHyPy",1,ptUp);
  // Print
  cout << Form("%.1f < p_{T} < %.1f GeV/c: ",hPt->GetBinLowEdge(1),hPt->GetBinLowEdge(ptUp+1)) << " SigSubBkg Integral - Nr: " << endl;
  cout << " Data - Nr: " << hDRBgSubNr->Integral() << " Aw: " << hDRBgSubAw->Integral() << endl;
  cout << " Pythia+Hydjet - Nr: " << hDRBgSubNrHyPy->Integral() << " Aw: " << hDRBgSubAwHyPy->Integral() << endl;

  if (drawMode==1) TCanvas * c6 = new TCanvas("c6","c6",500,500);
  hDRBgSubNr->Scale(1./totPtBgSubNr);
  hDRBgSubAw->Scale(1./totPtBgSubAw);
  hDRBgSubNrHyPy->Scale(1./totPtBgSubNrHyPy);
  hDRBgSubAwHyPy->Scale(1./totPtBgSubAwHyPy);

  if(cumulative){
     hDRBgSubNr = IntegrateFromLeft(hDRBgSubNr);
     hDRBgSubAw = IntegrateFromLeft(hDRBgSubAw);
     hDRBgSubNrHyPy = IntegrateFromLeft(hDRBgSubNrHyPy);
     hDRBgSubAwHyPy = IntegrateFromLeft(hDRBgSubAwHyPy);
  }

  // Set Styles
  hDRBgSubNr->SetMarkerStyle(kOpenSquare);
  mcStyle1(hDRBgSubNrHyPy);
  mcStyle2(hDRBgSubAwHyPy);
  hDRBgSubNrHyPy->SetMarkerStyle(0);
  hDRBgSubAwHyPy->SetMarkerStyle(0);
  // Draw
  hDRBgSubNrHyPy->SetTitle(";#DeltaR_{max};F(#DeltaR<#DeltaR_{max})");
  hDRBgSubNrHyPy->SetAxisRange(0,0.784,"X");
  hDRBgSubNrHyPy->SetAxisRange(0,0.7,"Y");
  if(!cumulative){
     hDRBgSubNrHyPy->SetAxisRange(0,0.1,"Y");
  }
  fixedFontHist(hDRBgSubNrHyPy);
  hDRBgSubNrHyPy->DrawCopy("Ehist");
  hDRBgSubAwHyPy->DrawCopy("Ehistsame");
  hDRBgSubNr->DrawCopy("Esame");
  hDRBgSubAw->DrawCopy("Esame");

  if (doLeg==1) {
     TLegend *leg = new TLegend(0.302407,0.67,0.7536548,0.9324599);
    leg->SetFillStyle(0);
    leg->SetBorderSize(0);
    leg->SetTextFont(63);
    leg->SetTextSize(16);
    leg->AddEntry(hDRBgSubNr,Form("%.1f < p_{T} < %.1f GeV/c",hPt->GetBinLowEdge(1),hPt->GetBinLowEdge(ptUp+1)),"");
    leg->AddEntry(hDRBgSubNr,"Data Leading Jet","pl");
    leg->AddEntry(hDRBgSubAw,"Data SubLeading Jet","pl");
    leg->AddEntry(hDRBgSubNrHyPy,"MC Leading Jet","l");
    leg->AddEntry(hDRBgSubAwHyPy,"MC SubLeading Jet","l");
    leg->Draw();
  }
}