C++ TH1F::GetMeanError方法代码示例

void GetRebalanceCorrectionFactor()
{
   gROOT->ProcessLine(".L /afs/desy.de/user/c/csander/xxl-af-cms/CMSSW_7_2_3_patch1/src/macros/PlottingUtils.C+");

   SetPlotStyle();

   // --- define output file for ps ---//
   TString pt = "pt10";
   double x_min = 10.;

   TString outfile = "RebalanceCorrectionFactor/RebalanceCorrectionFactors_madgraph_phys14_withoutPUReweighting_" + pt;

   TH2F* RebCorrection_vsReco = new TH2F("RebCorrection_vsReco", "Jet pt", 1000, 0., 1000., 100, 0., 3.);
   
   string root_file;

   // madgraph
   ifstream myfile ("filelists_phys14/filelist_madgraph_phys14_v1.txt");
   if (myfile.is_open()) {
      while( myfile.good() ) {
         getline (myfile,root_file);
         cout << root_file << endl;
                     
         TH2F* RebCorrection_vsReco_temp; 
               
         TString path = root_file;

         TFile* input_file = TFile::Open(path, "READ");
         input_file->cd("QCDfromSmearing");

         gDirectory->GetObject("RebCorrection_vsReco;1", RebCorrection_vsReco_temp);
         RebCorrection_vsReco->Add(RebCorrection_vsReco_temp);

         input_file->Close();
         
      }
      myfile.close();
   }

   // ---------------------------------------------------- //
   // madgraph
   TH1F* correction_vsReco = new TH1F();
   correction_vsReco = (TH1F*) RebCorrection_vsReco->ProjectionX();
   correction_vsReco->Reset();
   for (int i = 0; i <= RebCorrection_vsReco->GetXaxis()->GetNbins(); ++i) {
      TH1F h = *((TH1F*) RebCorrection_vsReco->ProjectionY("py", i, i));
            
      double mean = h.GetMean();
      double error = h.GetMeanError();

      cout << "i: " << i << " " << "mean: " << mean << " " << "error: " << error << endl;
            
      correction_vsReco->SetBinContent(i, mean);
      correction_vsReco->SetBinError(i, error);
   }

   // ---------------------------------------------------- //

   TCanvas *c = new TCanvas("c", "", CanvasPlot[2], CanvasPlot[3]);
   correction_vsReco->SetAxisRange(x_min, correction_vsReco->GetXaxis()->GetXmax());
   correction_vsReco->SetXTitle("reco jet p_{T} [GeV]");
   correction_vsReco->SetYTitle("reb jet p_{T} / gen jet p_{T} ");
   correction_vsReco->Draw();

   TPaveText* pt1 = CMSLabelMC("4.0");
   pt1->Draw();

   c->Print(outfile + "_vsReco.ps");  
   c->Print(outfile + "_vsReco.png"); 

   // ---------------------------------------------------- //

   TFile* RebalanceCorrection = new TFile(outfile + ".root", "RECREATE");
   correction_vsReco->Write();

   RebalanceCorrection->Write();
}

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

  //d->SaveAs("tmp/d_hist.root");
  s->SaveAs("tmp/s_hist.root");
  s2->SaveAs("tmp/s2_hist.root");


  yLabel.Form("Events /(%.2f %s)", s->GetBinWidth(2), yLabelUnit.Data());
  float max = 1.1 * std::max(
			     d->GetMaximum(),///d->Integral(),
			     s->GetMaximum() ///s->Integral()
			     );
  max=1.1*d->GetMaximum();
  std::cout << "max = " << max << std::endl;
  std::cout << "nEvents data: " << d->Integral() << "\t" << d->GetEntries() << std::endl;
  std::cout << "nEvents signal: " << s->Integral() << "\t" << s->GetEntries() << std::endl;
  std::cout << "nEvents signal2: " << s2->Integral() << "\t" << s2->GetEntries() << std::endl;
  if(logy){
    max*=10;
    d->GetYaxis()->SetRangeUser(0.1,max);
    s->GetYaxis()->SetRangeUser(0.1,max);
    s2->GetYaxis()->SetRangeUser(0.1,max);
    c->SetLogy();
  } else {
    d->GetYaxis()->SetRangeUser(0,max);
    s->GetYaxis()->SetRangeUser(0,max);
    s2->GetYaxis()->SetRangeUser(0,max);
  }
  s->GetYaxis()->SetTitle(yLabel);
  s->GetXaxis()->SetTitle(xLabel);
  s2->GetYaxis()->SetTitle(yLabel);
  s2->GetXaxis()->SetTitle(xLabel);
  d->GetYaxis()->SetTitle(yLabel);
  d->GetXaxis()->SetTitle(xLabel);


  d->SetMarkerStyle(20);
  d->SetMarkerSize(1);
  if(d != s){
    s->SetMarkerStyle(20);
    s->SetMarkerSize(1);
    s->SetFillStyle(3001);
    s->SetFillColor(kRed);
  }
  if(s2 != s){
    //s2->SetMarkerStyle(1);
    //s2->SetMarkerSize(0);
    //s->SetFillStyle(0);
    //s->SetFillColor(kB);
    s2->SetLineWidth(3);
    s2->SetLineColor(kBlack);
  }


  TH1F* s_norm = (TH1F *) (s->DrawNormalized("hist", d->Integral()));
  TH1F* s2_norm = (TH1F *) (s2->DrawNormalized("hist same", d->Integral()));
  //TH1F* d_norm = s_norm;
  //if(d!=s) d_norm = (TH1F *) (d->DrawNormalized("p same", d->Integral()));
  if(d!=s) d->Draw("p same");

  if(logy){
    //d_norm->GetYaxis()->SetRangeUser(0.1,max);
    s_norm->GetYaxis()->SetRangeUser(0.1,max);
    c->SetLogy();
  } else {
    //d_norm->GetYaxis()->SetRangeUser(0,max);  
    s_norm->GetYaxis()->SetRangeUser(0,max);  
  }
  std::cout << "Variable  & Data & Simulation & Simulation2 \\" << std::endl;
  std::cout << "Mean      & " << d->GetMean() << " " << d->GetMeanError() 
	    << " & " << s_norm->GetMean() <<  " " << s_norm->GetMeanError() 
	    << " & " << s2_norm->GetMean() <<  " " << s2_norm->GetMeanError() 
	    << " \\" << std::endl;
  std::cout << "Std. dev. & " << d->GetRMS() << " " << d->GetRMSError() 
	    << " & " << s_norm->GetRMS() << " " << s_norm->GetRMSError() 
	    << " & " << s2_norm->GetRMS() << " " << s2_norm->GetRMSError() 
	    << " \\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\Chi^2$ " <<  d->Chi2Test(s_norm, "UW CHI2/NDF NORM") << std::endl;
  

  TLegend *leg = new TLegend(0.6,0.8,1,1);
  if(dataLabel !="") leg->AddEntry(d,dataLabel,"p");
  if(mcLabel   !="") leg->AddEntry(s,mcLabel, "lf");
  if(mc2Label   !="") leg->AddEntry(s2,mc2Label, "l");
  leg->SetBorderSize(1);
  leg->SetFillColor(0);
  leg->SetTextSize(0.04);
  if(dataLabel !="" && mcLabel !="") leg->Draw();
  //c->GetListOfPrimitives()->Add(leg,"");

  TPaveText *pv = new TPaveText(0.23,0.95,0.6,1,"NDC");
  pv->AddText("CMS Preliminary 2016");
  pv->SetFillColor(0);
  pv->SetBorderSize(0);
  pv->Draw();


  return c;

}

void rateStudy() {
  const int NRUNS = 25;
  const int NCH = 32;
  const int NBINS = 32;

  TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);

  TMultiGraph *mg = new TMultiGraph();

  TLegend *legend=new TLegend(0.65,0.15,0.88,0.55);
  legend->SetNColumns(4);
  legend->SetFillColor(0);

  TH1F* hRate = new TH1F("hRate", "hist", 32.0, 0, 8.0);

  //Color buffer
  const int NCOLORS = 32;
  int color[NCOLORS] = {73, 2, 3, 4, 99, 6, 7, 8, 9, 12, 28, 32, 34,
                        28, 50, 51, 56, 58, 88, 99, 1, 208, 209,
                        218, 212, 210, 221, 224, 225, 226, 227, 228 };

  ifstream fin;
  //fin.open("runlist.txt");
  fin.open("filter_runlist.txt");

  string line = "";

  TFile* out = new TFile("outtemp.root", "REACREATE");
  TH1F* h = new TH1F("h","hist", NBINS, 0, NBINS);
  TF1* pois = new TF1("pois","[0]*TMath::Poisson(x,[1])",0,50);
  TF1* ppp = new TF1("ppp","[0]*TMath::Power(0.5,x*[1])",0.01,1.0);


  for (int ch = 0; ch < NCH; ++ch) {

    //if ( ch==26 || ch==27 )
    //continue;

    //Graph points and errors
    Double_t x[NRUNS];
    Double_t y[NRUNS];
    Double_t errX[NRUNS] = {0};
    Double_t errY[NRUNS] = {0};

    int fileCounter = 0;
    while(getline(fin, line)) {
      vector<double> data = parse(line);
      stringstream filePath;
      filePath << "pmtratestudy/run" << data[0] << "*.root";
      cout << "opening file at " << filePath.str() << endl;
      cout << "file counter: " << fileCounter << " channel=" << ch << endl;
      //TFile* f = new TFile(filePath.str().c_str());
      //TTree* t = (TTree *)f->Get("eventtree");
      TChain* t = new TChain("eventtree");
      t->Add( filePath.str().c_str() );
      out->cd();

      x[fileCounter] = data[1];

      int nfires[NCH] = {0};
      int samples = 0;
      float chmax = 0.0;
      t->SetBranchAddress("nfires", &nfires);
      t->SetBranchAddress("samples", &samples);
      t->SetBranchAddress("chmax", &chmax);
      
      h->Reset();
      
      int nentries = t->GetEntries();
      for (int entry = 0; entry < nentries; ++entry) {
        t->GetEntry(entry);
        if (chmax < 100.0) {
          h->Fill(nfires[ch]);
        }
      }


      
      pois->SetParameter(0,1);
      pois->SetParameter(1, h->GetMean());

      h->Fit(pois,"RQ","",0,50);
      //TF1 *myfit = (TF1 *)h->GetFunction("pois");
      TF1 *myfit = (TF1 *)pois;
      Double_t lambda = myfit->GetParameter(1);  
      h->Draw();
      stringstream histFileName;
      histFileName << "hist/h" << data[0] << "_ch" << ch << ".png";
      c1->SaveAs(histFileName.str().c_str());
      //Graph with poisson method
#if 1
      y[fileCounter] = lambda / ((samples - 1) * 15.625E-6);
      errY[fileCounter] = myfit->GetParError(1) / ((samples - 1) * 15.625E-6);
#endif
      //Graph with mean method
#if 0
      y[fileCounter] = h->GetMean() / ((samples - 1) * 15.625E-6);
      errY[fileCounter] = h->GetMeanError() / ((samples - 1) * 15.625E-6);
#endif
      cout << x[fileCounter] << ", " << y[fileCounter] 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      while( myfile8.good() ) {
         getline (myfile8,root_file);
         cout << root_file << endl;
                     
         TH2F* RebCorrection_vsReco_NVtx21_Inf_temp; 
               
         TString path = root_file;

         TFile* input_file = TFile::Open(path, "READ");
         input_file->cd("QCDfromSmearing");

         gDirectory->GetObject("RebCorrection_vsReco;1", RebCorrection_vsReco_NVtx21_Inf_temp);
         RebCorrection_vsReco_NVtx21_Inf->Add(RebCorrection_vsReco_NVtx21_Inf_temp);

         input_file->Close();
         
      }
      myfile8.close();
   }

//    RebCorrection_vsReco_NVtx0_10->Rebin2D(2, 1);
//    RebCorrection_vsReco_NVtx11_20->Rebin2D(2, 1);
//    RebCorrection_vsReco_NVtx21_Inf->Rebin2D(2, 1);

   // ---------------------------------------------------- //
   // pythia
   TH1F* correction_vsReco = new TH1F();
   correction_vsReco = (TH1F*) RebCorrection_vsReco->ProjectionX();
   correction_vsReco->Reset();
   for (int i = 0; i <= RebCorrection_vsReco->GetXaxis()->GetNbins(); ++i) {
      TH1F h = *((TH1F*) RebCorrection_vsReco->ProjectionY("py", i, i));
            
      double mean = h.GetMean();
      double error = h.GetMeanError();

      cout << "i: " << i << " " << "mean: " << mean << " " << "error: " << error << endl;
            
      correction_vsReco->SetBinContent(i, mean);
      correction_vsReco->SetBinError(i, error);
   }

   // madgraph
   TH1F* correction_vsReco_madgraph = new TH1F();
   correction_vsReco_madgraph = (TH1F*) RebCorrection_vsReco_madgraph->ProjectionX();
   correction_vsReco_madgraph->Reset();
   for (int i = 0; i <= RebCorrection_vsReco_madgraph->GetXaxis()->GetNbins(); ++i) {
      TH1F h = *((TH1F*) RebCorrection_vsReco_madgraph->ProjectionY("py", i, i));
            
      double mean = h.GetMean();
      double error = h.GetMeanError();

      cout << "i: " << i << " " << "mean: " << mean << " " << "error: " << error << endl;
            
      correction_vsReco_madgraph->SetBinContent(i, mean);
      correction_vsReco_madgraph->SetBinError(i, error);
   }

   // pythia NVtx 0-10
   TH1F* correction_vsReco_NVtx0_10 = new TH1F();
   correction_vsReco_NVtx0_10 = (TH1F*) RebCorrection_vsReco_NVtx0_10->ProjectionX();
   correction_vsReco_NVtx0_10->Reset();
   for (int i = 0; i <= RebCorrection_vsReco_NVtx0_10->GetXaxis()->GetNbins(); ++i) {
      TH1F h = *((TH1F*) RebCorrection_vsReco_NVtx0_10->ProjectionY("py", i, i));
            
      double mean = h.GetMean();
      double error = h.GetMeanError();

//.........这里部分代码省略.........
	  if(!Histo) 
       { std::cout<<" Histo "<<thestring.Data()<<"_"<<i<<" not found."<<std::endl; i++; continue;}
 
      if(Histo->GetEntries()) hNhits->Fill(Histo->Integral());
	  
	  if(Histo->Integral()<20) //0.1
	   { //std::cout<<" Not enought entries for histo "<<thestring.Data()<<std::endl;
	    i++; continue;}
 
	  
	  detid = iter->first;

	  bool rmfit=false;

      if( rmfit || 
      // TIB modules
          // TIB - 1.4.2.5
      detid==369121605 || detid==369121606 || detid==369121614 || 
      detid==369121613 || detid==369121610 || detid==369121609 ||
          // TIB - 1.2.2.1
      detid==369121390 || detid==369121382 || detid==369121386 || 
      detid==369121385 || detid==369121389 || detid==369121381 ||
          // others in TIB  
      detid==369121437 || detid==369142077 || detid==369121722 || 
      detid==369125534 || detid==369137018 || detid==369121689 ||
      detid==369121765 || detid==369137045 || detid==369169740 ||
      detid==369121689 ||
      // TOB modules 
	      // TOB + 4.3.3.8
      detid/10==436281512 || detid/10==436281528 || detid/10==436281508 ||
      detid/10==436281524 || detid/10==436281520 || detid/10==436281516 ||
          // others in TOB  
      detid/10==436228249 || detid/10==436232694 || detid/10==436228805 ||
      detid/10==436244722 || detid/10==436245110 || detid/10==436249546 ||
      detid/10==436310808 || detid/10==436312136 || detid/10==436315600 ||
	      // without 'sensors' option 
      detid==436281512 || detid==436281528 || detid==436281508 ||
      detid==436281524 || detid==436281520 || detid==436281516 ||
      detid==436228249 || detid==436232694 || detid==436228805 ||
      detid==436244722 || detid==436245110 || detid==436249546 ||
      detid==436310808 || detid==436312136 || detid==436315600 || 
      // TID modules
      detid==402664070 || detid==402664110 ||
	  // TEC modules in small scans
      detid==470148196 || detid==470148200 || detid==470148204 ||
      detid==470148228 || detid==470148232 || detid==470148236 ||
      detid==470148240 || detid==470148261 || detid==470148262 ||
	  detid==470148265 || detid==470148266 || detid==470148292 ||
	  detid==470148296 || detid==470148300 || detid==470148304 ||
	  detid==470148324 || detid==470148328 || detid==470148332 ||
	  detid==470148336 || detid==470148340 )  { 
	    Histo->Write();
        std::cout << " Saving histo : " << thestring.Data() << std::endl;
      }  


	  if(rmfit) {nfitrm++; i++; continue;}

          int subdet = ((detid>>25)&0x7);
          int TECgeom=0;
          if(subdet==6) TECgeom = ((detid>>5)&0x7);

      // save values
	  detid = iter->first;
	  voltage  = *itVolt;
	  index = i;
	  errvoltage = 2 ;
	  Mean = Histo->GetMean();
	  errMean = Histo->GetMeanError();
	  RMS = Histo->GetRMS();
	  errRMS = Histo->GetRMSError();
	  nhits = (int) Histo->Integral();
	  tree->Fill();
	    
	  i++;

    }  

  }
  
  tree->Write();
  hNhits->Write();
  

  for(unsigned int ih=0; ih<commonHistos.size(); ih++) commonHistos[ih]->Write();

  std::map<ULong64_t, TProfile* >::iterator itMon;
  for(itMon=Monitors.begin(); itMon!=Monitors.end(); itMon++)
  {
    itMon->second->GetXaxis()->SetTimeDisplay(1);
	itMon->second->GetXaxis()->SetTimeFormat("%H:%M");
	itMon->second->GetXaxis()->SetTimeOffset(t_monitor_start);
	itMon->second->Write();
  }
  
  //// If you want to store all the individual detId histograms uncomments this line !!!!
  //myFile->Write();
  myFile->Close();

}

void plotPatDep(Int_t slugNum, TString device, TString asymType){

gROOT->Reset();
gROOT->SetStyle("Plain");

const Double_t nonLin = -0.01; //relative MD:BCM non-linearity

const Bool_t saveFile = kTRUE; //save canvases to a file
const TString plotDir = "subblocks";

Int_t histPos[7] = {1,2,3,6,8,9,12}; //positioning scheme within the canvas
TH1F *histoCharge[7]; //Aq histogram pointer array
TH1F *histoDet[7]; //MD histogram pointer array
TH1F *histoDetCorr[7]; //Aq-corrected MD histogram pointer array

//arrays containing the mean and error for the pattern-history-dependent asymmetries, the ordering is the same as in plotTitles
Double_t aqCorr[7]; Double_t aqErr[7];
Double_t detPatMean[7]; Double_t detPatErr[7];
Double_t detPatCorrMean[7]; Double_t detPatCorrErr[7];

//arrays to be used in graphing the relative sizes of the lowest-level pattern asymmetries
Double_t patType[4] = {1,2,3,4}; // 1=(-++--++-) 2=(-++-+--+) 3=(+--+-++-) 4=(+--++--+)
Double_t patTypeErr[4] = {0,0,0,0};
Double_t patAsymMean[4]; Double_t patAsymErr[4]; //temporary holders

const TString cutList[7]={
"1",
"pat_phase==2",
"pat_phase==2 && actual_previous_pattern_polarity==0",
"pat_phase==2 && actual_previous_pattern_polarity==1",
"pat_phase==1",
"pat_phase==1 && actual_previous_pattern_polarity==0",
"pat_phase==1 && actual_previous_pattern_polarity==1"
}; //cuts to impose pattern history cuts

/*
const TString cutList[7]={
"1",
"actual_pattern_polarity==0",
"actual_pattern_polarity==0 && actual_previous_pattern_polarity==0",
"actual_pattern_polarity==0 && actual_previous_pattern_polarity==1",
"actual_pattern_polarity==1",
"actual_pattern_polarity==1 && actual_previous_pattern_polarity==0",
"actual_pattern_polarity==1 && actual_previous_pattern_polarity==1"
}; //cuts to impose pattern history cuts
*/

const TString plotTitles[7]={
"All Patterns",
"-++-",
"-++- -++-",
"-++- +--+",
"+--+",
"+--+ -++-",
"+--+ +--+"
}; //titles corresponding to above cuts

//open the file
Char_t filename[100];
sprintf(filename,"/w/hallc/qweak/smacewan/scratch/slugrootfiles/QwPass3_slug%i.root",slugNum);
f = new TFile(filename);
if(!f->IsOpen()){
	std::cerr<<"Error opening ROOTFile "<<filename<<".\n"<<std::endl;
	return 0;
}
std::cout<<"Successfully opened ROOTFile "<<filename<<".\n"<<std::endl;

//grab the tree
TTree *slugTree = f->Get("slug");
gStyle->SetOptStat("neMmRr");
gStyle->SetOptFit(1);
gStyle->SetStatW(0.32);


//----------------------------------------------------------------------------------------
//CHARGE ASYMMETRY
TCanvas *cCharge = new TCanvas("cCharge","Pattern History Study: qwk_charge",1400,1100);
cCharge->Divide(3,4);

for (Int_t i=0; i<7; i++){
	histoCharge[i] = new TH1F(Form("%s",plotTitles[i].Data()),"",100,0,0);
	cCharge->cd(histPos[i]);
	histoCharge[i]->SetDirectory(0);

	slugTree->Draw(Form("1e6*asym_qwk_charge>>slug%i_q_%s",slugNum,plotTitles[i].Data()),Form("%s && ErrorFlag==0 && asym_qwk_charge_Device_Error_Code==0",cutList[i].Data()));
	
	TH1F *htemp = (TH1F*)gPad->GetPrimitive(Form("slug%i_q_%s",slugNum,plotTitles[i].Data()));
	aqCorr[i] = htemp->GetMean();
	aqErr[i] = htemp->GetMeanError();
}

cCharge->cd(10);
patAsymMean[0] = aqCorr[2]; patAsymErr[0] = aqErr[2];
patAsymMean[1] = aqCorr[3]; patAsymErr[1] = aqErr[3];
patAsymMean[2] = aqCorr[5]; patAsymErr[2] = aqErr[5];
patAsymMean[3] = aqCorr[6]; patAsymErr[3] = aqErr[6];

TGraphErrors *grp = new TGraphErrors(4,patType,patAsymMean,patTypeErr,patAsymErr);
grp->SetMarkerStyle(21);
grp->Draw("ap");
//.........这里部分代码省略.........

//--------------------------------------
//function to calculate sampling factors
std::pair<Double_t,Double_t> g4_sample(int snum, Double_t energy, bool do_pion, bool do_show, bool do_print=false, bool set_val=true){
	Sample* sp = sample_map[snum];
	if(!sp) { std::cout << "Sample " << snum << " is not loaded." << std::endl; return std::pair<Double_t,Double_t>(0.,0.); }

	//select correct file
	std::string fpre = sp->fpre;
	if(do_pion) fpre += "_pion";
	else fpre += "_elec";

	//make filenames
	std::stringstream drawname, fname, piname;
	fname << sp->dir << "/" << fpre << "_" << energy << "gev_10k.root";
	if(do_pion) piname << "#pi^{-} " << energy << " GeV";
	else piname << "e^{-} " << energy << " GeV";

	//open file and tree
	TFile* _file;
	_file = TFile::Open((fname.str()).c_str());
	TTree* totalTree = (TTree*)_file->Get("Total");

	//get histo from tree (no display)
	//define mip as sam_ecal*ecal < 1 gev = 1000 mev (for pions in HCAL)
	if(sp->det==Hcal) drawname << "(hcal+" << sp->zeroWt << "*zero)/1000>>hsam(200)";
	else drawname << "(ecal)/1000>>hsam(200)";
	
	totalTree->Draw((drawname.str()).c_str(),"","hist goff");
	TH1F* hsam = (TH1F*)gDirectory->Get("hsam");
	
	//use parameters from histo to start fit
	TSpectrum* spec = new TSpectrum(5);
	spec->Search(hsam,5,"nodraw goff");
	Float_t* xpos = spec->GetPositionX();
	Float_t* ypos = spec->GetPositionY();

	Double_t m = xpos[0];
	Double_t me = hsam->GetMeanError();
	Double_t N = hsam->GetEntries();
	std::stringstream s_mean;
	s_mean.precision(3);
	Double_t f = energy/m;
	Double_t f_err = energy*(me/(m*m));
	s_mean << f << " #pm " << f_err;

	TPolyMarker* pm = new TPolyMarker(1, xpos, ypos);
	hsam->GetListOfFunctions()->Add(pm);
	pm->SetMarkerStyle(23);
	pm->SetMarkerColor(kRed);
	pm->SetMarkerSize(1.3);

	std::cout.precision(6);
	std::cout << "f_" << (do_pion ? "pion" : "elec") << " = " << f << " +/- " << f_err << std::endl;
	
	//plotting and printing
	if (do_show){
		TCanvas* can = new TCanvas("sample","sample",700,500);
		can->cd();
		TPad* pad = new TPad("graph","",0,0,1,1);
		pad->SetMargin(0.12,0.05,0.15,0.05);
		pad->Draw();
		pad->cd();
		
		//formatting
		hsam->SetTitle("");
		hsam->GetXaxis()->SetTitle("Energy [GeV]");
		//hsam->SetStats(kTRUE);
		//gStyle->SetOptStat("mr");
		hsam->SetLineWidth(2);
		hsam->SetLineColor(kBlack);
		hsam->GetYaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetYaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetYaxis()->SetTickLength(12/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetTickLength(12/(pad->GetWh()*pad->GetAbsHNDC()));
		
		hsam->Draw();
		
		std::stringstream Nname;
		Nname << "N = " << N;
		
		//determine placing of pave
		Double_t xmin;
		if (m/((hsam->GetXaxis()->GetXmax() + hsam->GetXaxis()->GetXmin())/2) < 1) xmin = 0.65;
		else xmin = 0.2;
		
		//legend
		TPaveText *pave = new TPaveText(xmin,0.65,xmin+0.2,0.85,"NDC");
		pave->AddText((piname.str()).c_str());
		pave->AddText((Nname.str()).c_str());
		pave->AddText("Peak sampling factor:");
		pave->AddText((s_mean.str()).c_str());
		pave->SetFillColor(0);
		pave->SetBorderSize(0);
		pave->SetTextFont(42);
		pave->SetTextSize(0.05);
		pave->Draw("same");

		if(do_print) {
//.........这里部分代码省略.........