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

//------------------------------------------------------------------------
void DrawQTC()
{
  
  TCanvas *c1 = new TCanvas("c1", "QTC C side",0,48,1280,951);
   c1->Divide(4,3);
   // c1->Divide(2,3);
  
  Char_t buf1[10];
  for (Int_t i=0; i<12; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"QTC%i",i+1);
      
      TH1F *qtc = (TH1F*) gFile->Get(buf1);
      
      Float_t mean = qtc->GetMean();
      Float_t rms = qtc->GetRMS();
      Float_t hminR=mean - 0.1*mean;
      Float_t hmaxR =mean + 0.1*mean;
      qtc->GetXaxis()->SetRange(hminR,hmaxR);
      Float_t hmin=mean - 3*rms;
      Float_t hmax =mean + 3*rms;
      qtc->GetXaxis()->SetRange(hmin,hmax);

     // TF1 *g2 = new TF1("g2", "gaus", hmin, hmax);
       //             qtc->Fit("g2","RQ");
	qtc->GetXaxis()->SetLabelSize(0.03);
       
       qtc->Draw();
    }
  
  TCanvas *c2 = new TCanvas("c2", "QTC A side",0,48,1280,951);
   c2->Divide(4,3);
   // c1->Divide(2,3);
  
  Char_t buf1[10];
  for (Int_t i=12; i<24; i++)
    {
      c2->cd(i+1-12);
      sprintf(buf1,"QTC%i",i+1);
      
      TH1F *qtc = (TH1F*) gFile->Get(buf1);
      
      Float_t mean = qtc->GetMean();
      Float_t rms = qtc->GetRMS();
      Float_t hminR=mean - 0.1*mean;
      Float_t hmaxR =mean + 0.1*mean;
      qtc->GetXaxis()->SetRange(hminR,hmaxR);
      Float_t hmin=mean - 3*rms;
      Float_t hmax =mean + 3*rms;
      qtc->GetXaxis()->SetRange(hmin,hmax);
      // TF1 *g2 = new TF1("g2", "gaus", hmin, hmax);
       //             qtc->Fit("g2","RQ");
 	qtc->GetXaxis()->SetLabelSize(0.03);
      
       qtc->Draw();
    }
  
  
}

void GenerateStatErrorMC(double Nt, double Nl, double Nt_err, double Nl_err, double &NtS_err, double &NtZ_err, double Eff, double Fake)
{
  double NtS, NtZ;
  TH1F* htS = new TH1F("htS", "htS", 11000, -100, 1000);
  TH1F* htZ = new TH1F("htZ", "htZ", 11000, -100, 1000);
  
  for(int itest=0; itest<1000; itest++)
  {
    double N_err = 0;
    if(Nt_err<=Nl_err) N_err = sqrt(Nl_err*Nl_err-Nt_err*Nt_err);
    else N_err = sqrt(Nl_err*Nl_err+Nt_err*Nt_err);
    double N_rand=gRandom->Gaus(Nl-Nt, N_err);
    double Nt_rand=gRandom->Gaus(Nt, Nt_err);
    double Nl_rand=N_rand+Nt_rand;
    
    ComputeEstimation(Nt_rand, Nl_rand, NtS, NtZ, Eff, Fake);
    htS->Fill(NtS);
    htZ->Fill(NtZ);
  }
  
  NtS_err = htS->GetRMS();
  NtZ_err = htZ->GetRMS();
  
  cout<<"S MC rms : "<<htS->GetRMS()<<endl;
  cout<<"Z MC rms : "<<htZ->GetRMS()<<endl;

}

vector<double> one_fit(bool do_fit)
{
    TH1F* h = new TH1F("h","h", 20, -10, 10);
    h->FillRandom("gaus");
    vector<double> ret;
    if (do_fit) {
        h->Fit("gaus","");
        TF1* g = h->GetFunction("gaus");
        h->Draw();
        canvas->Modified();
        canvas->Update();
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(g->GetParameter(ind));
        }
    }
    else {
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(0.0);
        }
    }
    
    ret.push_back(h->GetMean());
    ret.push_back(h->GetRMS());
    delete h;
    return ret;
}

void explore_zmass_boost()
{
  // Tell root not to draw everything to the screen.
  gROOT->SetBatch();

  //    TString dyfilename         = TString("/home/acarnes/h2mumu/samples/stage1/monte_carlo/bg/stage_1_dy_jetsToLL_asympt50_ALL.root");
  TString dyfilename         = TString("/home/acarnes/h2mumu/samples/stage1/monte_carlo/bg/stage_1_dy_ZToMuMu_asympt50_ALL.root");
  TString datafilename         = TString("/home/acarnes/h2mumu/samples/stage1/data_from_json/Cert_246908-251883_13TeV_PromptReco_Collisions15_JSON_v2/stage_1_doubleMuon_RunBPrompt_MINIAOD.root");

  TString savedir = TString("../png/dy_vs_data/run1cuts_v2_golden_json/dyzmumu/");

  // Initialize the DiMuPlottingSystems for MC and data
  DiMuPlottingSystem* dpsdata = new DiMuPlottingSystem(datafilename);
  DiMuPlottingSystem* dpsdy = new DiMuPlottingSystem(dyfilename);
  addDiMuMassPrimeBranch(dpsdata);
  dpsdata->applyRun1Cuts();
  dpsdy->applyRun1Cuts();

  // Get the 2D histos
  TH2F* hdataZPt = ZMassVsZPtHist2D("data_zpt", "recoCandMassPrime", "(14,0,60,30,87,95)", dpsdata);
  TH2F* hdyZPt = ZMassVsZPtHist2D("dy_zpt", "recoCandMass", "(14,0,60,30,87,95)", dpsdy);

  //overlayTProfiles(TH2F* hdata, TH2F* hdy, TString bininfo, TString savename)
  TProfile* p = overlayTProfiles(hdataZPt, hdyZPt, 90, 92, savedir+"z_mass_vs_z_pt.png");
  gStyle->SetOptFit(0011);
  fitTProfileCustom(p);
  TCanvas* c = new TCanvas();
  c->cd();
  p->Draw("hist c");
  p->Draw("E0 X0 same");
  dpsdata->arrangeStatBox(c);
  c->Draw();
  c->Print("blah2.png");

  // Look at RMS
  c->Clear();
  c->cd();
  c->SetGridx(kTRUE);
  c->SetGridy(kTRUE);
  c->cd();
  TH1F* h = dpsdata->hist1D("recoCandMass", "(100,87,95)", "");
  TH1F* h2 = dpsdata->hist1D("recoCandMassPrime", "(100,87,95)", "");
  h2->SetLineColor(2);
  h->Draw("");
  h2->Draw("same");

  std::cout << "mass  rms: " << h->GetRMS() << std::endl;
  std::cout << "mass' rms: " << h2->GetRMS() << std::endl;
  std::cout << "mass  mean: " << h->GetMean() << std::endl;
  std::cout << "mass' mean: " << h2->GetMean() << std::endl;
  std::cout << "mass  num: " << h->Integral() << std::endl;
  std::cout << "mass' num: " << h2->Integral() << std::endl;

  DiMuPlottingSystem* dps = new DiMuPlottingSystem();
  dps->arrangeStatBox(c);
  c->Draw();
  c->Print("newfit.png");
}

void PlotHistsNhitsPerModule(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
{
  TString canvName="c_";
  canvName+=strMillepedeRes;
  canvName+="_";
  canvName+=StrPlotType(NHITS);
  canvName.ReplaceAll(".res","");


  //enum {PXB,PXF,TIB,TID,TOB,TEC};
  int colors[6]={1,2,3,4,6,7};
//  TString labels[6]={"PXB","PXF","TIB","TID","TOB","TEC"};

  f->cd();
  TCanvas* canv = new TCanvas(canvName,canvName,600,600);
  canv->SetLogx();
  canv->SetLogy();

  for (int ind=1; ind<=1; ind++){
    TString strHist = "hNhits_";
    strHist+=StrPar(ind);
    TString strCut="label<700000 && ((label%20-1)%9+1)==";
    strCut+=ind;
    TStyle style; 
    style.SetTitleFontSize(0.2);
    THStack *hSt = new THStack("hNhits","# of derivatives (~tracks or hits) per module");
    TLegend *leg = new TLegend(0.75,0.65,0.95,0.95);
    for (int inv=0; inv<6; inv++){
      std::cout<<"- - - - - -"<<std::endl;
      std::cout<<subdLabels[inv]<<":"<<std::endl;
      std::cout<<StrCutSubd(inv)<<": "<<tr->GetEntries(StrCutSubd(inv))<<" parameters"<<std::endl;
      TString strHist1=strHist;
      strHist1+=ind;
      strHist1+=inv;
      TH1F* hValInt = new TH1F(strHist1,strHist1,300,10,15000);  
      TString strCut1 = strCut+TString(" && ")+StrCutSubd(inv);
      tr->Draw(TString("Nhits>>")+strHist1,strCut1,"goff");
      std::cout<<"# hits = "<<(int)hValInt->GetMean()<<"+-"<<(int)hValInt->GetRMS()<<std::endl;
      hValInt->SetLineColor(1);
      hValInt->SetFillColor(colors[inv]);
      hValInt->SetLineWidth(2);
      hSt->Add(hValInt);
      leg->AddEntry(hValInt,subdLabels[inv],"f");
      leg->SetFillColor(0);
    }
    hSt->Draw();
    leg->Draw("same");
    
  }//end of loop over ind

  canvName+=".png";
  TString saveName=strOutdir+canvName;
  canv->SaveAs(saveName);
  saveName.ReplaceAll(".png",".pdf");
  canv->SaveAs(saveName);
}//end of PlotHistsNhitsPerModule

void GenerateStatError(double Nt, double Nl, double &NtS_err, double &NtZ_err, double Eff, double Fake)
{
  double NtS, NtZ;
  TH1F* htS = new TH1F("htS", "htS", 11000, -100, 1000);
  TH1F* htZ = new TH1F("htZ", "htZ", 11000, -100, 1000);
  
  for(int itest=0; itest<1000; itest++)
  {
    double N_rand=gRandom->Poisson(Nl-Nt);
    double Nt_rand=gRandom->Poisson(Nt);
    double Nl_rand=N_rand+Nt_rand;
    
    ComputeEstimation(Nt_rand, Nl_rand, NtS, NtZ, Eff, Fake);
    htS->Fill(NtS);
    htZ->Fill(NtZ);
  }
  
  NtS_err = htS->GetRMS();
  NtZ_err = htZ->GetRMS();
  
  cout<<"S rms : "<<htS->GetRMS()<<endl;
  cout<<"Z rms : "<<htZ->GetRMS()<<endl;

}

void makePlot(vector<TTree*> sigTree,vector<double> sigWeight,
	      vector<int> ptHatLo, vector<int> ptHatHi,
              std::string var,TCut cut,TH1F* h,bool norm)
{
   TH1F *hRes = (TH1F*)h->Clone();
   hRes->SetName("hRes");
   hRes->Sumw2();
  
   char tmp[300];
   for (unsigned int i=0; i<sigTree.size(); i++)
     {
       // first determine the pthat cut
       
       sprintf(tmp, "ptHat >= %d && ptHat <= %d",ptHatLo[i], ptHatHi[i]);
       TCut ptHatCut = tmp;
       TCut allCut = cut + ptHatCut;
       cout << "Current cut = " << allCut.GetTitle() << endl;
       TH1F *htmp = (TH1F*)h->Clone();
       htmp->SetName("htmp");
       sigTree[i]->Draw(Form("%s>>htmp",var.data()),allCut);
       htmp->Sumw2();
       htmp->Scale(sigWeight[i]);
       cout << "scale = " << sigWeight[i] << endl;
       cout << "After scaling htmp -> entries() " << htmp->GetEntries() << endl;
       cout << "After scaling htmp -> Integral() " << htmp->Integral() << endl;
       cout << "After scaling htmp -> GetMean()  " << htmp->GetMean() << endl;
       cout << "After scaling htmp -> GetRMS()  " << htmp->GetRMS() << endl;
       hRes->Add(htmp);
       delete htmp;
   }
   h->Sumw2();
   h->Add(hRes);
   if(norm)h->Scale(1.0/(double)h->Integral(0,1000));
   cout << "After scaling h-> entries() " << h->GetEntries() << endl;
   cout << "After scaling h-> Integral() " << h->Integral() << endl;
   cout << "After scaling h -> GetMean()  " << h->GetMean() << endl;
   cout << "After scaling h -> GetRMS()  " << h->GetRMS() << endl;
   
   delete hRes;
}

void DumpRecenterParv2(){
    TString pro = "Pro104";
    int taxi = 8583;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit("Empty run list file!");

    TFile *fout = new TFile(Form("Recentering.root"),"Recreate");
    TVectorD vecmean;
    TVectorD vecrms;
    vecmean.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    vecrms.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    for(int irun=0;irun<nrun;irun++){
     cout<<irun<<" of total "<<nrun<<" runs"<<endl;
     int RunNumber=GetRun(irun);
     fin = TFile::Open(Form("Run15pAu200MinBias/output_fvtxwithcntrecenter_%d.root",RunNumber));
     //ofstream fout(Form("Calibration/%s/Recentering_%d.dat",dataset.Data(),GetRun(irun)));
     for(int icent=0;icent<ncent;icent++){
      for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
       for(int ihar=0;ihar<nhar;ihar++){
        for(int isub=0;isub<nsub;isub++){
         for(int ixy=0;ixy<nxy;ixy++){
             TH1F* q = (TH1F*)fin->Get(Form("q_%d_%d_%d_%d_%d",icent,ibbcz,ihar,isub,ixy));
             float mean = q->GetMean();
             float rms = q->GetRMS();
             vecmean[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = mean;
             vecrms[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = rms;
         }
        }
       }
      }
     }
    fout->cd();
    vecmean.Write(Form("mean_%d",RunNumber));
    vecrms.Write(Form("rms_%d",RunNumber));
    fin->Close();
    }
}

void  PlotAlignmentValidation::plotDMR(const std::string variable, Int_t minHits )
{
 setNiceStyle(); 
 gStyle->SetOptStat(0);

 // TList treeList=getTreeList();

 TCanvas *c = new TCanvas("canv", "canv", 600, 600);
 setCanvasStyle( *c );
 //loop over sub-detectors 
 for (int i=1;i<7;++i){
 
   int histo_Counter=1;
   TLegend *leg_hist = new TLegend(0.17,0.8,0.85,0.88);
   setLegendStyle(*leg_hist);
   //loop over file list
   //TTree *tree= (TTree*)treeList.First();
   //binning
   int nbinsX=100;
   double xmin=0;
   double xmax=0;
   float maxY=0;
   bool isHisto = false;
   std::string plotVar=variable;
   THStack *hstack=new THStack("hstack","hstack");

   for(std::vector<TkOfflineVariables*>::iterator it = sourceList.begin();
       it != sourceList.end(); ++it){
      
     //while ( tree ){
     plotVar=variable;
     TString subdet = "entries>=";
     subdet+=minHits; 
     subdet+=" && subDetId==";
     subdet+=i;
      
     char binning [50]="";
     sprintf (binning, ">>myhisto(%d,  %f , %f)", nbinsX, xmin, xmax);
     TH1F *h = 0;
     
     if (histo_Counter==1&&plotVar=="meanX")(*it)->getTree()->Draw( (plotVar+=">>myhisto(50,-0.001,0.001)").c_str(),subdet,"goff");
     else if (histo_Counter==1&&plotVar=="meanY")(*it)->getTree()->Draw( (plotVar+=">>myhisto(50,-0.005,0.005)").c_str(),subdet,"goff");
     else if (histo_Counter==1&&plotVar=="medianX")(*it)->getTree()->Draw( (plotVar+=">>myhisto(50,-0.005,0.005)").c_str(),subdet,"goff");
     else if (histo_Counter==1&&plotVar=="medianY")(*it)->getTree()->Draw( (plotVar+=">>myhisto(50,-0.005,0.005)").c_str(),subdet,"goff");
     else if (histo_Counter==1&&plotVar=="meanNormX")(*it)->getTree()->Draw( (plotVar+=">>myhisto(100,-2,2)").c_str(),subdet,"goff");
     else if (histo_Counter==1&&plotVar=="rmsX")(*it)->getTree()->Draw( (plotVar+=">>myhisto(100,0.,0.1)").c_str(),subdet,"goff");
     else if (histo_Counter!=1)(*it)->getTree()->Draw( (plotVar+=binning).c_str(),subdet,"goff");

     if (gDirectory) gDirectory->GetObject("myhisto", h);
     std::pair<float,float> fitResults(9999., 9999.);
     if (h){
       if (h->GetEntries()>0) {
	 isHisto = true;
	 h->SetDirectory(0);
	 //general draw options
	 h->SetLineWidth(2);
	 //first histo only, setting optStat...
	 if (histo_Counter==1)
	   setHistStyle(*h,plotVar.c_str() ,"#modules", 1 ); //set color later

	 h->SetLineColor( (*it)->getLineColor() );
	 h->SetLineStyle( (*it)->getLineStyle() );
	   //h->SetMarkerStyle(20+file_Counter);
      
	 //draw options
	 
	 if (maxY<h->GetMaximum()){
	   maxY=h->GetMaximum();
	 }
      
	 //fit histogram for median and mean
	 if (variable=="medianX"||variable =="meanX")fitResults=fitGauss(h, (*it)->getLineColor() );

	 if (histo_Counter==1){
	   //get mean and sigma from fit: gauss for 2sigma range
	   hstack->Add(h);
	   nbinsX=h->GetXaxis()->GetNbins();
	   xmin=h->GetXaxis()->GetXmin();
	   xmax=h->GetXaxis()->GetXmax();
	  
	 }else if (histo_Counter!=1 &&  h->GetEntries()>0)hstack->Add(h);
     
	 char legend [50]="";
	 std::string legEntry = (*it)->getName();
	 if ( (variable=="medianX"||variable =="meanX") && useFit_)
	   sprintf (legend, "%s: #mu = %4.2f#mum, #sigma = %4.2f#mum ",legEntry.c_str(),fitResults.first ,fitResults.second);
	 if ( (variable=="medianX"||variable =="meanX"||variable=="medianY"||variable =="meanY" )&&useFit_ == false)
	   sprintf (legend, "%s: #mu = %4.2f#mum, rms = %4.2f#mum ",legEntry.c_str(),h->GetMean(1)*10000 ,h->GetRMS(1)*10000);
	 else sprintf (legend, "%s ",legEntry.c_str());
	 if(h)
	   leg_hist->AddEntry(h,legend,"l");
	 else
	   std::cerr<< "histogram did not exist!";
       }
     }
     //     tree= (TTree*)treeList.After( tree );
     //     file_Counter++;
     histo_Counter++;
        
   }
//.........这里部分代码省略.........

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

  //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 drawNoisePlotsMC(TString InputFileName, TString suffix)
{
  TFile* f = TFile::Open(InputFileName.Data(),"READ");
  f -> cd();
  
  
  TH1F* h_occupancy_vsNvtx   = (TH1F*)( f->Get("h_occupancy_vsNvtx_EB") );
  TH1F* h_occupancy_vsNavgPU = (TH1F*)( f->Get("h_occupancy_vsNavgPU_EB") );
  TH1F* h_occupancy_vsIeta   = (TH1F*)( f->Get("h_occupancy_vsIeta_EB") );
  TH1F* h_occupancy_vsIring  = (TH1F*)( f->Get("h_occupancy_vsIring_EE") );
  
  
    
  TGraphErrors* g_RMSNoise_vsNvtx_EB = new TGraphErrors();
  TGraphErrors* g_RMSNoise_vsNvtx_EE = new TGraphErrors();

  TGraphErrors* g_RMSNoise_vsNavgPU_EB = new TGraphErrors();
  TGraphErrors* g_RMSNoise_vsNavgPU_EE = new TGraphErrors();
  
  int point = 0;
  for(int bin = 1; bin <= nBins_nVtx; ++bin)
  {  
    TH1F* histo;
    char histoName[50];
    
    float binCenter  = h_occupancy_vsNvtx -> GetBinCenter(bin);
    float binLowEdge = h_occupancy_vsNvtx -> GetBinLowEdge(bin);
    float binHigEdge = h_occupancy_vsNvtx -> GetBinLowEdge(bin) + h_occupancy_vsNvtx_EB->GetBinWidth(bin);
    
    if( binCenter < 5. ) continue;
    if( binCenter > 30. ) continue;
        
    sprintf(histoName,"hEB_recHitE_nVtx%02.1f-%02.1f",binLowEdge,binHigEdge);
    histo = (TH1F*)( f->Get(histoName) );
    g_RMSNoise_vsNvtx_EB -> SetPoint(point,binCenter,histo->GetRMS());
    g_RMSNoise_vsNvtx_EB -> SetPointError(point,0.,histo->GetRMSError());
    
    sprintf(histoName,"hEE_recHitE_nVtx%02.1f-%02.1f",binLowEdge,binHigEdge);
    histo = (TH1F*)( f->Get(histoName) );
    g_RMSNoise_vsNvtx_EE -> SetPoint(point,binCenter,histo->GetRMS());
    g_RMSNoise_vsNvtx_EE -> SetPointError(point,0.,histo->GetRMSError());    
    
    ++point;
  }
  
  g_RMSNoise_vsNvtx_EB -> SetPoint(g_RMSNoise_vsNvtx_EB->GetN(),-1.,-1.);
  g_RMSNoise_vsNvtx_EE -> SetPoint(g_RMSNoise_vsNvtx_EE->GetN(),-1.,-1.);
  
  
  
  point = 0;
  for(int bin = 1; bin <= nBins_nAvgPU; ++bin)
  {  
    TH1F* histo;
    char histoName[50];
    
    float binCenter  = h_occupancy_vsNavgPU -> GetBinCenter(bin);
    float binLowEdge = h_occupancy_vsNavgPU -> GetBinLowEdge(bin);
    float binHigEdge = h_occupancy_vsNavgPU -> GetBinLowEdge(bin) + h_occupancy_vsNavgPU_EB->GetBinWidth(bin);
    
    if( binCenter < 8. ) continue;
    if( binCenter > 30. ) continue;
        
    sprintf(histoName,"hEB_recHitE_nAvgPU%02.1f-%02.1f",binLowEdge,binHigEdge);
    histo = (TH1F*)( f->Get(histoName) );
    g_RMSNoise_vsNavgPU_EB -> SetPoint(point,binCenter,histo->GetRMS());
    g_RMSNoise_vsNavgPU_EB -> SetPointError(bin-1,0.,histo->GetRMSError());
    
    sprintf(histoName,"hEE_recHitE_nAvgPU%02.1f-%02.1f",binLowEdge,binHigEdge);
    histo = (TH1F*)( f->Get(histoName) );
    g_RMSNoise_vsNavgPU_EE -> SetPoint(point,binCenter,histo->GetRMS());
    g_RMSNoise_vsNavgPU_EE -> SetPointError(bin-1,0.,histo->GetRMSError());    
    
    ++point;
  }
  
  
  
  
  
  
  TCanvas* c1_RMSNoise_EB = new TCanvas("c1_RMSNoise_EB","EB RMS noise");
  c1_RMSNoise_EB -> cd();
  c1_RMSNoise_EB -> SetGridx();
  c1_RMSNoise_EB -> SetGridy();
  
  g_RMSNoise_vsNvtx_EB -> SetMarkerColor(kRed+2);
  g_RMSNoise_vsNavgPU_EB -> SetMarkerColor(kGreen+2);
  
  g_RMSNoise_vsNvtx_EB -> GetXaxis() -> SetTitle("N_{vtx} or Num^{Int}_{True}");
  g_RMSNoise_vsNvtx_EB -> GetYaxis() -> SetTitle("noise RMS (GeV)");
  g_RMSNoise_vsNvtx_EB -> GetXaxis() -> SetRangeUser(0.,35.);
  g_RMSNoise_vsNvtx_EB -> GetYaxis() -> SetRangeUser(0.040,0.080);
  
  g_RMSNoise_vsNvtx_EB -> Draw("AP");
  g_RMSNoise_vsNavgPU_EB -> Draw("P,same");
  
  TF1* f_pol1_vsNvtx_EB   = new TF1("f_pol1_vsNvtx_EB","[0]+[1]*x",0.,35.);
  f_pol1_vsNvtx_EB -> SetLineColor(kRed+2);
  TF1* f_pol1_vsNavgPU_EB = new TF1("f_pol1_vsNavgPU_EB","[0]+[1]*x",0.,35.);
//.........这里部分代码省略.........

double runGraphicHistos(TString ifile = "tmpfitdirc.root", bool verbose_out=true, double ienergy = 5, int iupdown = 0)
{
	TCanvas *c1 = new TCanvas("myc1","myc1",1000,1200);

	if (verbose_out == false)
	{
		//supresses canvas message 
		//gErrorIgnoreLevl=kInfo;
		gROOT->ProcessLine( "gErrorIgnoreLevel = kWarning;");
	}

	double hmin = -100;
	double hmax = 100;

	hmin = -25;
	hmax = 25;

	double energy = ienergy;

	double pi_mass = .13957;
	double k_mass = .49367;

	double pi_beta = sqrt(1-pi_mass*pi_mass/(energy*energy));
	double k_beta = sqrt(1-k_mass*k_mass/(energy*energy));

	double quartz_index = 1.47;

	double pi_mrad = 1000*acos(1/(pi_beta*quartz_index));
	double k_mrad = 1000*acos(1/(k_beta*quartz_index));

	//double seperation = 5.2;
	double seperation = pi_mrad - k_mrad;
	double mean_pion = 0;
	double mean_kaon = mean_pion + seperation;
	if (verbose_out == true)
	{
		printf("Energy: %4.02f\n",energy);
		printf("Mrad Seperation: %8.03f\n",seperation);
	}
	//denominator pf spread

	TRandom3* randgen = new TRandom3();

	TFile *f1 = new TFile(ifile);
	TH1F *hpion;
	TH1F *hkaon;
	TH1F *phots_pion;
	if (iupdown == 0)
	{
		hpion = (TH1F*) f1->Get("ll_diff_pion");
		hkaon = (TH1F*) f1->Get("ll_diff_kaon");
		phots_pion = (TH1F*) f1->Get("phot_found_pion");
	}
	else if (iupdown == 1)
	{
		hpion = (TH1F*) f1->Get("ll_diff_pion_up");
		hkaon = (TH1F*) f1->Get("ll_diff_kaon_up");
		phots_pion = (TH1F*) f1->Get("phot_found_pion_up");
	}
	else if (iupdown == -1)
	{
		hpion = (TH1F*) f1->Get("ll_diff_pion_down");
		hkaon = (TH1F*) f1->Get("ll_diff_kaon_down");
		phots_pion = (TH1F*) f1->Get("phot_found_pion_down");
	}
	else
	{
		printf("Unrecognize updown arguement: %d \nFailing....\n",iupdown);
		return -1;
	}
	if (verbose_out == true)
	{
		printf("pion_ll mean, spread: %12.04f, %12.04f\n",hpion->GetMean(),hpion->GetRMS());
		printf("kaon_ll mean, spread: %12.04f, %12.04f\n",hkaon->GetMean(),hkaon->GetRMS());
	}
	double spread = seperation/2;
	double spreadsq2 = 2*spread*spread;
	//Swap kaon and pion numbers

	for (int i = 1; i < hpion->GetNbinsX()/2; i++)
	{
		double t_swap = hpion->GetBinContent(i);
		hpion->SetBinContent(i,hpion->GetBinContent(hpion->GetNbinsX() - i + 1));
		hpion->SetBinContent(hpion->GetNbinsX() - i + 1, t_swap);
	}
	for (int i = 1; i < hkaon->GetNbinsX()/2; i++)
	{
		double t_swap = hkaon->GetBinContent(i);
		hkaon->SetBinContent(i,hkaon->GetBinContent(hkaon->GetNbinsX() - i + 1));
		hkaon->SetBinContent(hkaon->GetNbinsX() - i + 1, t_swap);
	}

	double titlesize=1.2*.04;

	int rebin = 20;
	rebin = 160;
	hpion->Rebin(rebin);
	hkaon->Rebin(rebin);

	hpion->SetAxisRange(hmin,hmax);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    
    //g -> SetMinimum(min-0.05*fabs(max-min));
    //g -> SetMaximum(max+0.05*fabs(max-min));
    g -> SetLineWidth(2);
    g -> SetLineColor(fibIt+1);
    g -> SetMarkerSize(0.2);
    g -> GetXaxis() -> SetTitle("sample time (ns)");
    g -> Draw("APL");
    
    c_waveform_cut_fib_all -> cd();
    
    if( plotIt == 0 ) g -> Draw("APL");
    else              g -> Draw("PL,same");
    
    ++plotIt;
  }
  
  c_waveform_cut_fib_all -> Print(Form("/afs/cern.ch/user/a/abenagli/www/TBatFNAL/%s/waveform_cut_fib_all.png",inFileName.c_str()),"png");
  
  
  
  for(int fibIt = 0; fibIt < nFib; ++fibIt)
  {
    TH1F* h = h_ped_fib[fibIt];
    
    TCanvas* c_ped_fib = new TCanvas();
    c_ped_fib -> SetLogy();
    
    h -> SetLineWidth(2);
    h -> GetXaxis() -> SetTitle("max sample");
    h -> Draw();
    h -> Fit("gaus","Q");
    
    TLatex* latex1 = new TLatex(0.60,0.90,Form("RMS = %.1f",h->GetRMS()));
    latex1 -> SetNDC();
    latex1 -> SetTextFont(42);
    latex1 -> SetTextSize(0.04);
    latex1 -> Draw("same");
    
    TLatex* latex2 = new TLatex(0.60,0.85,Form("#sigma = %.1f",h->GetFunction("gaus")->GetParameter(2)));
    latex2 -> SetNDC();
    latex2 -> SetTextFont(42);
    latex2 -> SetTextSize(0.04);
    latex2 -> Draw("same");
    
    c_ped_fib -> Print(Form("/afs/cern.ch/user/a/abenagli/www/TBatFNAL/%s/plotsPerFib/ped_fib%02d.png",inFileName.c_str(),fibIt),"png");
    
    h = h_maximum_fib[fibIt];
    
    TCanvas* c_maximum_fib = new TCanvas();
    c_maximum_fib -> SetLogy();
    
    h -> SetLineWidth(2);
    h -> GetXaxis() -> SetTitle("max sample");
    h -> Draw();
    
    c_maximum_fib -> Print(Form("/afs/cern.ch/user/a/abenagli/www/TBatFNAL/%s/plotsPerFib/maximum_fib%02d.png",inFileName.c_str(),fibIt),"png");
  }
  
  TCanvas* c_ped_fib_all = new TCanvas();
  c_ped_fib_all -> SetLogy();
  
  h_ped_fib_all -> SetLineWidth(2);
  h_ped_fib_all -> GetXaxis() -> SetTitle("pedestal");
  h_ped_fib_all -> Draw();
  h_ped_fib_all -> Fit("gaus","Q");

//.........这里部分代码省略.........
      sprintf(ctit,"phaseHF/espec/E_+%d_%d_1",ieta,iphi);
      hLdatPx[ii][ll/2]  =   new TH1F(*((TH1F*)fila->Get(ctit)));
      hLdatPx[ii][ll/2]->SetAxisRange(Ethr1,Ethr2);
      rLP = hLdatPx[ii][ll/2]->Integral()*hLdatPx[ii][ll/2]->GetMean();
      hLmapP0->SetBinContent(ii+1,ll/2+1,rLP);
      sprintf(ctit,"phaseHF/espec/E_+%d_%d_2",ieta,iphi);
      hSdatPx[ii][ll/2]  =   new TH1F(*((TH1F*)fila->Get(ctit)));
      hSdatPx[ii][ll/2]->SetAxisRange(Ethr1,Ethr2);
      rSP = hSdatPx[ii][ll/2]->Integral()*hSdatPx[ii][ll/2]->GetMean();
      hSmapP0->SetBinContent(ii+1,ll/2+1,rSP);
      if (ieta<=32 && iphi==67) continue;
      if (rLP>0) {
	htL->Fill(rLP);
	mLE += rLP;
	nmLE++;
      }
      if (rSP>0) {
	htS->Fill(rSP);
	mSE += rSP;
	nmSE++;
      }
    }
    if (nmLE>0) mLE /= nmLE; 
    else continue;
    if (nmSE>0) mSE /= nmSE; 
    else continue;
    ccxx->cd(1); htL->Draw("hist");
    ccxx->cd(2); htS->Draw("hist");
    ccxx->Update();
    //histspec(htL,mLE,ermean,rms,4,3);
    //histspec(htS,mSE,ermean,rms,4,3);
    mLEphi[ii]=mLE;
    mSEphi[ii]=mSE;
    dmLEphi[ii]=htL->GetRMS();
    dmSEphi[ii]=htS->GetRMS();
    printf("ieta %2d :  <E>L= %8.1f (%6.1f) x %d    <E>S= %8.1f (%6.1f) x %d \n",
	   ieta,mLE,dmLEphi[ii],nmLE,mSE,dmSEphi[ii],nmSE);
    
    for (int jj=1;jj<=72;jj+=2) {
      int iphi=jj;
      if (abs(ieta)>39 && (iphi-1)%4==0) continue;
      if (ieta<=32 && iphi==67) {
	hLmapP->SetBinContent(ii+1,jj/2+1,hLmapP0->GetBinContent(ii+1,jj/2+1));
	hSmapP->SetBinContent(ii+1,jj/2+1,hSmapP0->GetBinContent(ii+1,jj/2+1));
	continue;
      }

      for (nIter=1;nIter<nIterN;nIter++) { //cout<<nIter<<" |  ";
	corrL=hLmapPc->GetBinContent(ii+1,jj/2+1);
	hLdatP[ii][jj/2]->Reset();

	for (int kk=1;kk<=hLdatPx[ii][jj/2]->GetNbinsX();kk++) {
	  xxL[kk-1]=hLdatPx[ii][jj/2]->GetBinCenter(kk);
	  yyL[kk-1]=hLdatPx[ii][jj/2]->GetBinContent(kk);
	}
	ttL = new TSpline5("tt",xxL,yyL,1000,"",10,20);

	for (int kk=1;kk<=hLdatP[ii][jj/2]->GetNbinsX();kk++) {
	  x=hLdatP[ii][jj/2]->GetBinCenter(kk);
	  y=hLdatP[ii][jj/2]->GetBinContent(kk);
	  hLdatP[ii][jj/2]->Fill(x*corrL,ttL->Eval(x)/8.0);
	}
	ttL->Delete();

	hLdatP[ii][jj/2]->SetAxisRange(Ethr1,Ethr2);
	rLP = hLdatP[ii][jj/2]->Integral()*hLdatP[ii][jj/2]->GetMean();

//.........这里部分代码省略.........
	      fHistCosp[index]->SetLineWidth(2);
	      fHistCosp[index]->GetXaxis()->SetTitle("pT (GeV/c)");
	      fHistCosp[index]->GetYaxis()->SetTitle("Cos Point");
	      fHistCosp[index]->SetTitle("Prompt CosPoint vs pT");

	      if (index % 2 == 1)
		fHistDecLen[index]->SetTitle("FeedDown Dec Len vs pT");

	      htemp = (TH1F *)fHistPtYMultGenDauInAcc[index]->ProjectionX(Form("hPtDen%s%s", type[i].Data(), names[j].Data()));
	      hEffPt[index] = (TH1F *)fHistPtYMultRecoFilt[index]->ProjectionX(Form("hPtNum%s%s", type[i].Data(), names[j].Data()));
	      fHistPtYMultGenDauInAcc[index]->Sumw2();
	      fHistPtYMultRecoFilt[index]->Sumw2();
	      hEffPt[index]->Sumw2();
	      hEffPt[index]->Divide(htemp);
	      hEffPt[index]->SetLineColor(color[j]);
	      hEffPt[index]->SetLineWidth(2);
	      hEffPt[index]->GetXaxis()->SetTitle("pT (GeV/c)");
	      hEffPt[index]->GetYaxis()->SetTitle("Prompt Efficiency");
	      hEffPt[index]->SetTitle("Prompt Efficiency");
	      hEffPt[index]->SetStats(0);
	      fHistCosp[index]->SetStats(0);
	      fHistDecLen[index]->SetStats(0);

	      if (index % 2 == 1)
		{
		  hEffPt[index]->GetYaxis()->SetTitle("Feeddown Efficiency");
		  hEffPt[index]->SetTitle("Feeddown Efficiency");
		}

	      fHistInvMassRes[index] = new TH1F(*hEffPt[index]);
	      for (Int_t jj = 1; jj < hEffPt[index]->GetNbinsX() + 1; jj++)
		{
		  TH1F *hTemp = (TH1F *)fHistInvMassVsPt[index]->ProjectionY("htemp", jj, jj);
		  fHistInvMassRes[index]->SetBinContent(jj, hTemp->GetRMS());
		  fHistInvMassRes[index]->SetBinError(jj, hTemp->GetRMSError());
		  fHistInvMassRes[index]->SetLineColor(color[j]);
		  fHistInvMassRes[index]->SetLineWidth(2);
		  fHistInvMassRes[index]->GetXaxis()->SetTitle("pT (GeV/c)");
		  fHistInvMassRes[index]->GetYaxis()->SetTitle("Inv Mass RMS (GeV/c2)");
		  fHistInvMassRes[index]->SetTitle("Inv Mass RMS vs pT");
		  if (index == 0)
		    printf("D0: pt=%f, res=%f \n", fHistInvMassRes[index]->GetBinCenter(jj), fHistInvMassRes[index]->GetBinContent(jj));
		  TH1F *hTempX = (TH1F *)fHistXvtxResVsPt[index]->ProjectionY("htempX", jj, jj);
		  TH1F *hTempY = (TH1F *)fHistYvtxResVsPt[index]->ProjectionY("htempY", jj, jj);
		  TH1F *hTempZ = (TH1F *)fHistZvtxResVsPt[index]->ProjectionY("htempZ", jj, jj);

		  fHistXvtxRes2[index]->SetBinContent(jj, hTempX->GetRMS());
		  fHistXvtxRes2[index]->SetBinError(jj, hTempX->GetRMSError());

		  fHistYvtxRes2[index]->SetBinContent(jj, hTempY->GetRMS());
		  fHistYvtxRes2[index]->SetBinError(jj, hTempY->GetRMSError());

		  fHistZvtxRes2[index]->SetBinContent(jj, hTempZ->GetRMS());
		  fHistZvtxRes2[index]->SetBinError(jj, hTempZ->GetRMSError());
		}
	    }
	}

      fHistNCandDplus->SetLineColor(2);
      fHistNCandDstar->SetLineColor(3);
      fHistNCandDs->SetLineColor(4);
      fHistNCandLc->SetLineColor(kOrange);
      fHistNCandDplus->SetLineWidth(2);
      fHistNCandDstar->SetLineWidth(2);
      fHistNCandDs->SetLineWidth(2);
      fHistNCandLc->SetLineWidth(2);