C++ TGaxis::SetTitle方法代码示例

ReverseXAxis (TGraphAsymmErrors *g, double  size)
{
// Remove the current axis
   g->GetXaxis()->SetLabelOffset(999);
   g->GetXaxis()->SetTickLength(0);

//  Redraw the new axis 
   gPad->Update();

	g->GetYaxis()->SetTitle("Event Selection Efficiency");
   TGaxis *newaxis = new TGaxis(gPad->GetUxmax(),
                                gPad->GetUymin(),
                                gPad->GetUxmin(),
                                gPad->GetUymin(),
                                g->GetXaxis()->GetXmin(),
                                g->GetXaxis()->GetXmax(),
                                510,"-");
   newaxis->SetLabelOffset(-0.04);
   newaxis->SetTitle("Centrality");
   newaxis->SetTitleFont(42);
   newaxis->SetTitleSize(size);
   newaxis->SetLabelSize(size);
   newaxis->CenterTitle();
   newaxis->SetLabelFont(42);
   newaxis->Draw();


}

void ReverseXAxis (TH1 *h)
{
   h->GetXaxis()->SetLabelOffset(999);
   h->GetXaxis()->SetTickLength(0);

   gPad->Update();
   TGaxis *newaxis = new TGaxis(gPad->GetUxmax(),
                                gPad->GetUymin(),
                                gPad->GetUxmin(),
                                gPad->GetUymin(),
                                h->GetXaxis()->GetXmin(),
                                h->GetXaxis()->GetXmax(),
                                510,"-");
   newaxis->SetLabelOffset(-0.03);
   newaxis->SetTitle("X[cm]");
   newaxis->CenterTitle(1);
   newaxis->Draw();
}

ExtraAxis anotherScale (const TH1* refHist, double scale, int color, const char* title, double offset) {
  ExtraAxis result;
  double x0 = refHist->GetXaxis()->GetXmin();
  double x1 = refHist->GetXaxis()->GetXmax();
  double y0 = refHist->GetMinimum();
  double y1 = refHist->GetMaximum();
  // double y0 = refHist->GetYaxis()->GetXmin();
  // double y1 = refHist->GetYaxis()->GetXmax();
  double xoffset = exp (log(x0) - (log(x1) - log(x0))*offset);

  TGaxis* axis = new TGaxis(xoffset, y0, xoffset, y1, y0*scale,y1*scale,510,"-GS");
  axis->ImportAxisAttributes (refHist->GetXaxis());
  axis->SetTitle(title);
  axis->SetTextColor (color);
  axis->SetLineColor (color);
  axis->SetLineWidth (1);
  axis->SetTextColor (color);
  axis->SetLabelColor (color);
  axis->SetLabelOffset (0.);
  axis->SetTitleOffset (0.65);
  axis->SetTickSize(0.015);	
  result.Add (axis);

  TLine* line = new TLine (xoffset, y0, xoffset, y1);
  line->SetLineColor (color);
  line->SetLineWidth (2);
  result.Add (line);

  line = new TLine (x0, y0, xoffset, y0);
  line->SetLineColor (kGray);
  line->SetLineWidth (2);
  result.Add (line);
  
  line = new TLine (x0, y1, xoffset, y1);
  line->SetLineColor (kGray);
  line->SetLineWidth (2);
  result.Add (line);

  return result;
}

//.........这里部分代码省略.........
        TH1F * th1SignFull=new TH1F("SignFull","SignFull",numberBinFull,0,numberBinFull);
        TH1F * th1EffFull=new TH1F("th1EffFull","th1EffFull",numberBinFull,0,numberBinFull);
        TH1F * th1Eff2Full=new TH1F("th1Eff2Full","th1Eff2Full",numberBinFull,0,numberBinFull);
        for(int i=0; i<numberBin; i++) {
            th1Sign->SetBinContent(i+1,sign[signI[i]][signJ[i]]);
            th1Sign->SetBinError(i+1,signErr[signI[i]][signJ[i]]);
            th1Eff->SetBinContent(i+1,eff[signI[i]][signJ[i]]);
            th1Eff2->SetBinContent(i+1,eff2[signI[i]][signJ[i]]);
            th1Eff->SetBinError(i+1,err[signI[i]][signJ[i]]);
            th1Eff2->SetBinError(i+1,err2[signI[i]][signJ[i]]);
            th1Sign->GetXaxis()->SetBinLabel(i+1,Form("%d-%d",bmin[signJ[i]],bmin[signJ[i]]+width[signI[i]]));
            //cout<<"range=["<<bmin[signJ[i]]<<","<<bmin[signJ[i]]+width[signI[i]]<<"] and significant="<<sign[signI[i]][signJ[i]]<<endl;
        }

        //cout<<"signI.size="<<signI.size()<<endl;
        for(unsigned int i=0; i<signI.size(); i++) {
            twikiSign[massP][signI[i]][signJ[i]]=i+1;

            th1SignFull->SetBinContent(i+1,sign[signI[i]][signJ[i]]);
            th1SignFull->SetBinError(i+1,signErr[signI[i]][signJ[i]]);
            th1EffFull->SetBinContent(i+1,eff[signI[i]][signJ[i]]);
            th1Eff2Full->SetBinContent(i+1,eff2[signI[i]][signJ[i]]);
            th1EffFull->SetBinError(i+1,err[signI[i]][signJ[i]]);
            th1Eff2Full->SetBinError(i+1,err2[signI[i]][signJ[i]]);
            //th1SignFull->GetXaxis()->SetBinLabel(i+1,Form("%d-%d",bmin[signJ[i]],bmin[signJ[i]]+width[signI[i]]));
            th1SignFull->GetXaxis()->SetBinLabel(i+1,"");
            twikiSignValue[massP][signI[i]][signJ[i]]=sign[signI[i]][signJ[i]];
            twikiEffValue[massP][signI[i]][signJ[i]]=eff[signI[i]][signJ[i]];
            //cout<<"range=["<<bmin[signJ[i]]<<","<<bmin[signJ[i]]+width[signI[i]]<<"] and significant="<<sign[signI[i]][signJ[i]]<<endl;
        }
        gStyle->SetOptStat(0000000000);


        th1Sign->SetTitle(Form("largest 15 Sign.(eff),%sGev",masspoint[massP].data()));
        th1Sign->SetMinimum(0);
        th1Sign->SetMaximum(1);
        th1Sign->SetLineColor(4);
        th1Sign->Draw();
        th1Eff->SetLineColor(2);
        th1Eff2->SetLineColor(3);
        th1Eff->Draw("same");
        th1Eff2->Draw("same");
        TLegend* leg ;
        leg=new TLegend(x1NDC,y1NDC,x2NDC,y2NDC);
        setLeg(leg);
        leg->AddEntry(th1Sign,"significance");
        leg->AddEntry(th1Eff,"signal efficiency");
        leg->AddEntry(th1Eff2,"Bkg. efficiency");
        leg->Draw("same");

        if(massP==0)c1->Print(Form("pdfPR/%s.pdf(",output.data()));
        //else if(massP==12)c1->Print("pdf/signv2.pdf)");
        else c1->Print(Form("pdfPR/%s.pdf",output.data()));

        th1SignFull->SetTitle(Form("all windows(eff),%sGev",masspoint[massP].data()));
        th1SignFull->SetMinimum(0);
        th1SignFull->SetMaximum(1);
        th1SignFull->SetLineColor(4);
        th1SignFull->Draw("Hist");
        th1EffFull->SetLineColor(2);
        th1Eff2Full->SetLineColor(3);
        th1EffFull->Draw("same");
        th1Eff2Full->Draw("same");
        leg->Draw("same");

        c1->Print(Form("pdfPR/%s.pdf",output.data()));

//.........这里部分代码省略.........
            DeltaPt2simPtHist->SetBinContent(PtIndex, DeltaPtMean);
            DeltaPt2simPtHist->SetBinError(PtIndex, DeltaPtRMS);

            TH1D* DeltaPhiHist = DeltaPhimaxPurity2simPtHist->ProjectionY("DeltaPhi", PtIndex, PtIndex, "o");
            double DeltaPhiMean = DeltaPhiHist->GetMean();
            double DeltaPhiRMS = DeltaPhiHist->GetRMS();
            DeltaPhi2simPtHist->SetBinContent(PtIndex, DeltaPhiMean);
            DeltaPhi2simPtHist->SetBinError(PtIndex, DeltaPhiRMS);

            TH1D* DeltaEtaHist = DeltaEtamaxPurity2simPtHist->ProjectionY("DeltaEta", PtIndex, PtIndex, "o");
            double DeltaEtaMean = DeltaEtaHist->GetMean();
            double DeltaEtaRMS = DeltaEtaHist->GetRMS();
            DeltaEta2simPtHist->SetBinContent(PtIndex, DeltaEtaMean);
            DeltaEta2simPtHist->SetBinError(PtIndex, DeltaEtaRMS);
        }
        myEfficiencyHist->AddLast(Efficiency2simPtHist);
        myParticleHist->AddLast(Particle2simPtHist);
        mySTAHist->AddLast(STA2simPtHist);
        myChargeCheckHist->AddLast(InverseChargeRato2simPtHist);
        myDeltaPtHist->AddLast(DeltaPt2simPtHist);
    }
    double minX = 0;
    double minY = 0;
    double maxX = 110;
    double maxY = 40;

    TCanvas* myCanvas = new TCanvas("Canvas", "Canvas", 800, 600);
    myCanvas->cd();
    TPad* myPad = new TPad("Pad", "Pad", 0, 0, 1, 1);
    myPad->Draw();
    myPad->cd();

    ((TH1D*)(myParticleHist->At(0)))->SetStats(0);
    ((TH1D*)(myParticleHist->At(0)))->GetXaxis()->SetTitle("simPt/Gev");
    ((TH1D*)(myParticleHist->At(0)))->GetXaxis()->CenterTitle(1);
    ((TH1D*)(myParticleHist->At(0)))->Draw();
    for(int Index = 0; Index < FileNumber; Index++) {
	((TH1D*)(mySTAHist->At(Index)))->SetStats(0);
        ((TH1D*)(mySTAHist->At(Index)))->SetLineColor(kRed+Index);
        ((TH1D*)(mySTAHist->At(Index)))->Draw("same");
    }
    TLegend *STALeg = new TLegend(0.6,0.1,0.9,0.3);
    STALeg->SetBorderSize(1);
    TString LegKey = "ParticleTrack";
    STALeg->AddEntry(myParticleHist->At(0), LegKey, "lpf");
    for(int Index = 0; Index < FileNumber; Index++) {
        LegKey = TypeName[Index];
        STALeg->AddEntry(mySTAHist->At(Index), LegKey, "lpf");
    }
    STALeg->Draw();
    string SaveName = OutputPlotNamepreFix + "_STA2simPt" + OutputPlotNameFix;
    myCanvas->SaveAs(SaveName.c_str());

    myPad->Clear();
    myPad->Update();
    double YScale = myPad->GetUymax() / 110.;
    ((TH1D*)(myEfficiencyHist->At(0)))->GetXaxis()->SetTitle("simPt/Gev");
    ((TH1D*)(myEfficiencyHist->At(0)))->GetXaxis()->CenterTitle(1);
    ((TH1D*)(myEfficiencyHist->At(0)))->SetStats(0);
    ((TH1D*)(myEfficiencyHist->At(0)))->Scale(YScale);
    ((TH1D*)(myEfficiencyHist->At(0)))->SetLineColor(kRed);
    ((TH1D*)(myEfficiencyHist->At(0)))->Draw("same,ah");
    for(int Index = 1; Index < FileNumber; Index++) {
	((TH1D*)(myEfficiencyHist->At(Index)))->SetStats(0);
        ((TH1D*)(myEfficiencyHist->At(Index)))->Scale(YScale);
        ((TH1D*)(myEfficiencyHist->At(Index)))->SetLineColor(kRed+Index);

//.........这里部分代码省略.........
  tt->SetBranchAddress("PFJet1Eta", &PFJet1Eta_, &b_PFJet1Eta);
  tt->SetBranchAddress("PFJet2Eta", &PFJet2Eta_, &b_PFJet2Eta);
  tt->SetBranchAddress("PFDeltaEta", &PFDeltaEta_, &b_PFDeltaEta);

  tt->SetBranchAddress("hltAccept", &hltAccept, &b_hltAccept);
  tt->SetBranchAddress("l1Accept", &l1Accept, &b_l1Accept);
  tt->SetBranchAddress("l1Names", &l1Names, &b_l1Names);

  int nentries = tt->GetEntries();
  std::cout << "Number of entries: " << nentries << std::endl;

  //book graphs and plots
  float min = 0.;
  float max = 1000.;
  int nBins = 20;

  TF1* f1 = new TF1("f1","[0]*TMath::Erf((x-[1])/[2])-[0]*TMath::Erf((-x-[1])/[2])",min,max);
  f1->SetParameters(0.5,350,40);  
  f1->FixParameter(0,0.5);
  f1->SetLineWidth(2.);
  f1->SetLineColor(kRed);

  TF1* f2 = (TF1*)f1->Clone("f2");
  f2->SetParameters(0.5,150,10);
  f2->SetLineColor(kBlack);

  TH1F* caloMjjSpectrum = new TH1F("caloMjjSpectrum","caloMjjSpectrum",nBins,min,max);
  TH1F* PFMjjSpectrum = new TH1F("PFMjjSpectrum","PFMjjSpectrum",nBins,min,max);
  
  TEfficiency* mjj450_eff = new TEfficiency("mjj450_eff","mjj450_eff",nBins,min,max);
  mjj450_eff->SetMarkerColor(kRed);
  mjj450_eff->SetLineColor(kRed);
  mjj450_eff->SetLineWidth(2);
  mjj450_eff->SetTitle("turnOn;Mjj [GeV]");
  TEfficiency* mjj200_eff = new TEfficiency("mjj200_eff","mjj200_eff",nBins,min,max);
  mjj200_eff->SetLineWidth(2);
  mjj200_eff->SetTitle("turnOn;Mjj [GeV]");
  TEfficiency* pf410_eff = new TEfficiency("pf410_eff","pf410_eff",nBins,min,max);
  pf410_eff->SetMarkerColor(kOrange+1);
  pf410_eff->SetLineColor(kOrange+1);
  TEfficiency* calo250_eff = new TEfficiency("calo250_eff","calo250_eff",nBins,min,max);
  calo250_eff->SetMarkerColor(kBlue);
  calo250_eff->SetLineColor(kBlue);
  TEfficiency* HTT240_eff = new TEfficiency("HTT240_eff","HTT240_eff",nBins,min,max);
  HTT240_eff->SetMarkerColor(kGreen+2);
  HTT240_eff->SetLineColor(kGreen+2);

  TH1F* l1 = new TH1F("l1","l1",14,0.,14.);
  TH1F* l2 = new TH1F("l2","l2",14,0.,14.);
  
  //loop
  for (Long64_t jentry=0; jentry<nentries;++jentry)
    {
      tt->GetEntry(jentry);

      //remove low rate lumis.
      //see: https://cmswbm2.web.cern.ch/cmswbm2/cmsdb/servlet/ChartHLTTriggerRates?RUNID=274200&PATHID=2043408&LSLENGTH=23.31040958&TRIGGER_PATH=DST_HT250_CaloScouting_v2
      //if(lumi > 539 && lumi < 553) continue;

      //l1 and hlt rates
      for(unsigned int ii=0; ii<l1Names->size(); ++ii)
	if (l1Accept->at(ii)==1)
	  l1->Fill(ii);
      
      
      //analysis cuts needed to compare to the analysis

//.........这里部分代码省略.........
   const char * CharCanvasName = CanvasName.c_str();
   TCanvas* c1 = new TCanvas(CharCanvasName,"Efficiency Scan Plot",200,10,700,500);
   TPad *pad = new TPad("pad","",0,0,1,1);
   //pad->SetFillColor(42);
   pad->SetGrid();
   pad->Draw();
   pad->cd();

      // draw a frame to define the range
   TH1F *hr = pad->DrawFrame(0,-0.5,100,1.1);
   hr->SetXTitle("Detector Position (mm)");
   hr->SetYTitle("Efficiency");
   //pad->GetFrame()->SetFillColor(21);
   pad->GetFrame()->SetBorderSize(12);

      // create first graph
   TGraphErrors *gr_GIF = new TGraphErrors(V_Efficiency.size());
   TGraphErrors *gr_IV = new TGraphErrors(V_Efficiency.size());
   TGraphErrors *gr_V = new TGraphErrors(V_Efficiency.size());
   for(unsigned int i=0;i<V_Efficiency.size();i++)
   {
       gr_GIF->SetPoint(i, GIF_MeanPosOfSector[i], GIF_Efficiency[i]);
       gr_GIF->SetPointError(i,0, GIF_EfficiencyError[i]);
       
       gr_IV->SetPoint(i, IV_MeanPosOfSector[i], IV_Efficiency[i]);
       gr_IV->SetPointError(i,0, IV_EfficiencyError[i]);

       gr_V->SetPoint(i, V_MeanPosOfSector[i], V_Efficiency[i]);
       gr_V->SetPointError(i,0, V_EfficiencyError[i]);
   }
   gr_GIF->SetMarkerColor(kBlue);
   gr_GIF->SetLineColor(kBlue);
   gr_GIF->SetMarkerStyle(21);
   gr_GIF->GetXaxis()->SetTitle("dist (mm)");
   gr_GIF->GetYaxis()->SetTitle("Efficiency");
   gr_GIF->GetYaxis()->SetRangeUser(0,1.2);
   gr_GIF->SetTitle("Efficiency Scan");
   gr_GIF->Draw("ACP");
   gr_IV->SetMarkerColor(kGreen-6);
   gr_IV->SetLineColor(kGreen);
   gr_IV->SetMarkerStyle(21);
   gr_IV->Draw("sameCP");
   gr_V->SetMarkerColor(kBlack);
   gr_V->SetLineColor(kBlack);
   gr_V->SetMarkerStyle(21);
   gr_V->Draw("sameCP");

   //create a transparent pad drawn on top of the main pad
   c1->cd();
   TPad *overlay = new TPad("overlay","",0,0,1,1);
   overlay->SetFillStyle(4000);
   overlay->SetFillColor(0);
   overlay->SetFrameFillStyle(4000);
   overlay->Draw();
   overlay->cd();
   // create second graph
   //TGraphErrors* gr_GIF_Num = new TGraphErrors("data_noerror.dat","%lg %lg");
   TGraphErrors *gr_GIF_Num = new TGraphErrors(GIF_Nevents.size());
   TGraphErrors *gr_IV_Num = new TGraphErrors(IV_Nevents.size());
   TGraphErrors *gr_V_Num = new TGraphErrors(V_Nevents.size());

   for(unsigned int i=0;i<V_MeanPosOfSector.size();i++)
   {
       gr_GIF_Num->SetPoint(i,GIF_MeanPosOfSector[i],GIF_Nevents[i]);
       gr_IV_Num->SetPoint(i,IV_MeanPosOfSector[i],IV_Nevents[i]);
       gr_V_Num->SetPoint(i,V_MeanPosOfSector[i],V_Nevents[i]);

//.........这里部分代码省略.........
	    ++excl_mjjHltPassed;


	}
      float HTT240rate = (float)HTT240Passed/(float)nentries*PDRate;
      //std::cout << "HTT240rate = " << HTT240rate << std::endl;
      
      float sigmaNentries = sqrt((float)nentries);

      float sigmaMjjHltPassed = sqrt((float)mjjHltPassed);
      float excl_sigmaMjjHltPassed = sqrt((float)excl_mjjHltPassed);

      float sigmaMjjL1Passed = sqrt((float)mjjL1Passed);
      float excl_sigmaMjjL1Passed = sqrt((float)excl_mjjL1Passed);

      float mjjHltRate = (float)mjjHltPassed/(float)nentries*PDRate;
      float mjjHltRateE = PDRate*sqrt(pow((sigmaMjjHltPassed/nentries),2)+pow((sigmaNentries*mjjHltPassed/nentries/nentries),2));
      float excl_mjjHltRate = (float)excl_mjjHltPassed/(float)nentries*PDRate;
      float excl_mjjHltRateE = PDRate*sqrt(pow((excl_sigmaMjjHltPassed/nentries),2)+pow((sigmaNentries*excl_mjjHltPassed/nentries/nentries),2));

      float mjjL1Rate = (float)mjjL1Passed/(float)nentries*PDRate;
      float mjjL1RateE = PDRate*sqrt(pow((sigmaMjjL1Passed/nentries),2)+pow((sigmaNentries*mjjL1Passed/nentries/nentries),2));
      float excl_mjjL1Rate = (float)excl_mjjL1Passed/(float)nentries*PDRate;
      float excl_mjjL1RateE = PDRate*sqrt(pow((excl_sigmaMjjL1Passed/nentries),2)+pow((sigmaNentries*excl_mjjL1Passed/nentries/nentries),2));

      totHltRateVsCut->SetPoint(bin,cut,mjjHltRate);
      totHltRateVsCut->SetPointError(bin,0.,mjjHltRateE);
      pureHltRateVsCut->SetPoint(bin,cut,excl_mjjHltRate);
      pureHltRateVsCut->SetPointError(bin,0.,excl_mjjHltRateE);

      totL1RateVsCut->SetPoint(bin,cut,mjjL1Rate);
      totL1RateVsCut->SetPointError(bin,0.,mjjL1RateE);
      pureL1RateVsCut->SetPoint(bin,cut,excl_mjjL1Rate);
      pureL1RateVsCut->SetPointError(bin,0.,excl_mjjL1RateE);

      ++bin;
    }

  //plotting and styling
  TLegend* legHlt = new TLegend(0.62, 0.78, 0.83, 0.89);
  legHlt->AddEntry(totHltRateVsCut,"hlt total rate","P");
  legHlt->AddEntry(pureHltRateVsCut,"hlt pure rate wrt HT250","P");

  TLegend* legL1 = new TLegend(0.62, 0.78, 0.83, 0.89);
  legL1->AddEntry(totL1RateVsCut,"l1 total rate","P");
  legL1->AddEntry(pureL1RateVsCut,"l1 pure rate wrt HTT240","P");

  totHltRateVsCut->GetXaxis()->SetTitle("Mjj cut threshold [GeV]");
  totHltRateVsCut->GetYaxis()->SetTitle("Rate @7E33 [Hz]");
  totHltRateVsCut->SetMarkerColor(kRed);
  pureHltRateVsCut->SetMarkerColor(kOrange+1);

  totL1RateVsCut->GetXaxis()->SetTitle("Mjj cut threshold [GeV]");
  totL1RateVsCut->GetYaxis()->SetTitle("Rate @7E33 [Hz]");
  totL1RateVsCut->SetMarkerColor(kBlue);
  totL1RateVsCut->SetMarkerColor(kAzure);

  TCanvas* c4 = new TCanvas();
  c4->cd();
  totHltRateVsCut->Draw("AP");
  pureHltRateVsCut->Draw("P,sames");
  legHlt->Draw("sames");
  c4->Update();
  TGaxis *axisHlt = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(), gPad->GetUymax(),
    			    (totHltRateVsCut->GetYaxis()->GetBinLowEdge(1))*lumiScaleFactor,
			    (totHltRateVsCut->GetYaxis()->GetBinLowEdge(totHltRateVsCut->GetYaxis()->GetNbins())+totHltRateVsCut->GetYaxis()->GetBinWidth(1))*lumiScaleFactor,510,"+L");

  c4->SetTicky(0);
  axisHlt->SetLineColor(kRed);
  axisHlt->SetLabelColor(kRed);
  axisHlt->SetTextColor(kRed);
  axisHlt->SetTitleOffset(1.3);
  axisHlt->SetLabelSize(0.03);
  axisHlt->SetTitle("Rate @1E34 [Hz]");
  axisHlt->Draw();
  c4->Print("rates/hltRate.pdf","pdf");


  TCanvas* c5 = new TCanvas();
  c5->cd();
  totL1RateVsCut->Draw("AP");
  pureL1RateVsCut->Draw("P,sames");
  legL1->Draw("sames");
  c5->Update();
  TGaxis *axisL1 = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(), gPad->GetUymax(),
			      (totL1RateVsCut->GetYaxis()->GetBinLowEdge(1))*lumiScaleFactor,
			      (totL1RateVsCut->GetYaxis()->GetBinLowEdge(totL1RateVsCut->GetYaxis()->GetNbins())+totL1RateVsCut->GetYaxis()->GetBinWidth(1))*lumiScaleFactor,510,"+L");

  c5->SetTicky(0);
  axisL1->SetLineColor(kRed);
  axisL1->SetLabelColor(kRed);
  axisL1->SetTextColor(kRed);
  axisL1->SetTitleOffset(1.3);
  axisL1->SetLabelSize(0.03);
  axisL1->SetTitle("Rate @1E34 [Hz]");
  axisL1->Draw();
  c5->Print("rates/l1Rate.pdf","pdf");

  return 0;
}

//.........这里部分代码省略.........
    // Get charge density histos
    Int_t Nspecies = pData->NSpecies();
    TH2F **hDen2D = new TH2F*[Nspecies];
    // Get charge density on-axis
    TH1F **hDen1D = new TH1F*[Nspecies];
    // And electric current (integrated)
    TH1F **hCur1D = new TH1F*[Nspecies];
    for(Int_t i=0; i<Nspecies; i++) {

        hDen2D[i] = NULL;

        if(!pData->GetChargeFileName(i))
            continue;

        cout << Form(" Getting charge density of specie: ") << i << endl;


        char hName[24];
        sprintf(hName,"hDen2D_%i",i);
        hDen2D[i] = (TH2F*) gROOT->FindObject(hName);
        if(hDen2D[i]) delete hDen2D[i];

        if(!pData->Is3D())
            hDen2D[i] = pData->GetCharge(i,opt);
        else
            hDen2D[i] = pData->GetCharge2DSliceZY(i,-1,Nbins,opt+"avg");

        hDen2D[i]->SetName(hName);
        hDen2D[i]->GetXaxis()->CenterTitle();
        hDen2D[i]->GetYaxis()->CenterTitle();
        hDen2D[i]->GetZaxis()->CenterTitle();

        if(opt.Contains("comov"))
            hDen2D[i]->GetXaxis()->SetTitle("k_{p} #zeta");
        else
            hDen2D[i]->GetXaxis()->SetTitle("k_{p} z");

        if(pData->IsCyl())
            hDen2D[i]->GetYaxis()->SetTitle("k_{p} r");
        else
            hDen2D[i]->GetYaxis()->SetTitle("k_{p} y");

        hDen2D[i]->GetZaxis()->SetTitle("n [n_{0}]");


        hDen1D[i] = NULL;
        hCur1D[i] = NULL;

        if(!pData->GetEfieldFileName(i))
            continue;

        sprintf(hName,"hDen1D_%i",i);
        hDen1D[i] = (TH1F*) gROOT->FindObject(hName);
        if(hDen1D[i]) delete hDen1D[i];

        // 1D histograms
        if(pData->Is3D()) {
            hDen1D[i] = pData->GetH1SliceZ3D(pData->GetChargeFileName(i)->c_str(),"charge",-1,Nbins,-1,Nbins,opt+"avg");
        } else if(pData->IsCyl()) { // Cylindrical: The first bin with r>0 is actually the number 1 (not the 0).
            hDen1D[i] = pData->GetH1SliceZ(pData->GetChargeFileName(i)->c_str(),"charge",1,Nbins,opt+"avg");
        } else { // 2D cartesian
            hDen1D[i] = pData->GetH1SliceZ(pData->GetChargeFileName(i)->c_str(),"charge",-1,Nbins,opt+"avg");
        }
        hDen1D[i]->SetName(hName);

        // if(hDen1D[i]) delete hDen1D[i];

void PlotDensity1D( const TString &sim, Int_t time, Int_t Nbins=1, const TString &options="") {
  
#ifdef __CINT__  
  gSystem->Load("libplasma.so");
#endif

  PlasmaGlob::Initialize();
  
  // Init Units table
  PUnits::UnitsTable::Get();

  // Load PData
  PData *pData = PData::Get(sim.Data());
  pData->LoadFileNames(time);
  if(!pData->IsInit()) return; 

  TString opt = options;
 
  gStyle->SetPadRightMargin(0.20); // Margin right axis 
  if(opt.Contains("grid")) {
    gStyle->SetPadGridX(1);
    gStyle->SetPadGridY(1);
  }

  Bool_t CYL = kFALSE;
  if(sim.Contains("cyl")) { CYL = kTRUE; opt += "cyl"; } 
  
  Bool_t ThreeD = kFALSE;
  if(sim.Contains("3D")) ThreeD = kTRUE; 

  Bool_t INT = kTRUE; // Integrate instead of averaging.
  
  // Some plasma constants
  Float_t n0 = pData->GetPlasmaDensity(); 
  Float_t kp = pData->GetPlasmaK();       
  Float_t skindepth = (1/kp);               
  Float_t E0 = pData->GetPlasmaE0();

  // Some beam properties:
  Double_t Ebeam = pData->GetBeamEnergy();
  Double_t gamma = pData->GetBeamGamma();
  Double_t vbeam = pData->GetBeamVelocity();
  Double_t kbeta = PFunc::BeamBetatronWavenumber(gamma,n0);
  Double_t rms0  = pData->GetBeamRmsY() * kp;
  if(CYL)  rms0  = pData->GetBeamRmsR() * kp;
  
  cout << Form(" - Bunch gamma      = %8.4f", gamma ) << endl;
  cout << Form(" - Bunch velocity   = %8.4f c", vbeam ) << endl;
  cout << Form(" - Bunch betatron k = %8.4f mm-1", kbeta * PUnits::mm) << endl;
  cout << Form(" - Bunch RMS_0      = %8.4f um", rms0 * skindepth / PUnits::um) << endl;
  cout << endl;

  // Time in OU
  Float_t Time = pData->GetRealTime();
  // z start of the plasma in normalized units.
  Float_t zStartPlasma = pData->GetPlasmaStart()*kp;
  // z start of the beam in normalized units.
  Float_t zStartBeam = pData->GetBeamStart()*kp;
  
  if(opt.Contains("center")) {
    Time -= zStartPlasma;
    if(opt.Contains("comov"))      // Centers on the head of the beam.
      Time += zStartBeam;
  }
  
  // Get charge density histos
  Int_t Nspecies = pData->NSpecies();
  TH1F **hDen1D = new TH1F*[Nspecies];
  TH2F *hDen2D = NULL;
  for(Int_t i=0;i<Nspecies;i++) {
    hDen1D[i] = NULL;
    
    if(!pData->GetChargeFileName(i)) 
      continue;
    
    if(i==0) {
      if(ThreeD)
	hDen2D = pData->GetCharge2DSliceZY(i,-1,Nbins);
      else
	hDen2D = pData->GetCharge(i,opt);
      
      char hName[24];
      sprintf(hName,"hDen_%i",i);
      hDen2D->SetName(hName);
      hDen2D->GetXaxis()->CenterTitle();
      hDen2D->GetYaxis()->CenterTitle();
      hDen2D->GetZaxis()->CenterTitle();
      hDen2D->GetXaxis()->SetTitle("z [c/#omega_{p}]");
      hDen2D->GetYaxis()->SetTitle("y [c/#omega_{p}]");
      if(i==0)
	hDen2D->GetZaxis()->SetTitle("#LTn_{e}#GT [n_{0}]");
      else
	hDen2D->GetZaxis()->SetTitle("#LTn_{b}#GT [n_{0}]");
    }
    
    if(Nbins==0) {
      Nbins = TMath::Nint(rms0 / hDen2D->GetYaxis()->GetBinWidth(1)) ;
      // cout << Form(" Rms0 = %6.2f  Dx = %6.2f  Nbins = %4i .", 
      // 	   rms0, hDen2D[i]->GetYaxis()->GetBinWidth(1), Nbins) << endl;
    }
//.........这里部分代码省略.........

void PlotField1D( const TString &sim, Int_t time, Int_t Nbins=1, const TString &options="") { 
  
#ifdef __CINT__  
  gSystem->Load("libplasma.so");
#endif

  PlasmaGlob::Initialize();
    
  TString opt = options;

  gStyle->SetPadLeftMargin(0.10);  // Margin left axis  
  gStyle->SetPadRightMargin(0.12); // Margin right axis 
  if(opt.Contains("grid")) {
    gStyle->SetPadGridX(1);
    gStyle->SetPadGridY(1);
  }

  // Load PData
  PData *pData = PData::Get(sim.Data());
  pData->LoadFileNames(time);
  if(!pData->IsInit()) return;

  Bool_t CYL = kFALSE;
  if(sim.Contains("cyl")) { CYL = kTRUE; opt += "cyl"; } 
    
  Bool_t ThreeD = kFALSE;
  if(sim.Contains("3D")) ThreeD = kTRUE; 

  // Some plasma constants
  Double_t n0 = pData->GetPlasmaDensity();
  Double_t kp = pData->GetPlasmaK();
  Double_t skindepth = 1/kp;
  Double_t E0 = pData->GetPlasmaE0();

  // Some beam properties:
  Float_t Ebeam = pData->GetBeamEnergy() * PUnits::MeV;
  Float_t gamma = Ebeam / PConst::ElectronMassE;
  Float_t vbeam = TMath::Sqrt(1 - 1/(gamma*gamma));
  // cout << Form(" - Bunch gamma      = %8.4f", gamma ) << endl;
  // cout << Form(" - Bunch velocity   = %8.4f c", vbeam ) << endl;

  Float_t Time = pData->GetRealTime();
  // z start of the plasma in normalized units.
  Float_t zStartPlasma = pData->GetPlasmaStart()*kp;
  // z start of the beam in normalized units.
  Float_t zStartBeam = pData->GetBeamStart()*kp;
  Time -= zStartPlasma - zStartBeam;

  // 1D histograms
  TString opth1 = opt;
  opth1 += "avg";
  // Get electric fields
  const Int_t Nfields = 1;
  TH1F **hE1D = new TH1F*[Nfields];
  for(Int_t i=0;i<Nfields;i++) {
    hE1D[i] = NULL;
    if(!pData->GetEfieldFileName(i))
      continue;
    
    char nam[3]; sprintf(nam,"e%i",i+1);
    if(ThreeD) {
      if(i==0) 
	hE1D[i] = pData->GetH1SliceZ3D(pData->GetEfieldFileName(i)->c_str(),nam,-1,Nbins,-1,Nbins,opth1.Data());
      else  
	hE1D[i] = pData->GetH1SliceZ3D(pData->GetEfieldFileName(i)->c_str(),nam,-Nbins,Nbins,-Nbins,Nbins,opth1.Data());
    } else if(CYL) { // Cylindrical: The firt bin with r>0 is actually the number 1 (not the 0).
      if(i==0) 
	hE1D[i] = pData->GetH1SliceZ(pData->GetEfieldFileName(i)->c_str(),nam,1,Nbins,opth1.Data());
      else
	hE1D[i] = pData->GetH1SliceZ(pData->GetEfieldFileName(i)->c_str(),nam,1,Nbins,opth1.Data());
    } else { // 2D cartesian
      if(i==0) 
	hE1D[i] = pData->GetH1SliceZ(pData->GetEfieldFileName(i)->c_str(),nam,-1,Nbins,opth1.Data());
      else 
	hE1D[i] = pData->GetH1SliceZ(pData->GetEfieldFileName(i)->c_str(),nam,-Nbins,Nbins,opth1.Data());
    }
    
    char hName[24];
    sprintf(hName,"hE_%i",i);
    hE1D[i]->SetName(hName);
    if(opt.Contains("comov"))
      hE1D[i]->GetXaxis()->SetTitle("#zeta [c/#omega_{p}]");
    else
      hE1D[i]->GetXaxis()->SetTitle("z [c/#omega_{p}]");
   
    if(i==0)
      hE1D[i]->GetYaxis()->SetTitle("E_{z} [E_{0}]");
    else if(i==1)
      hE1D[i]->GetYaxis()->SetTitle("E_{y} [E_{0}}]");
    else if(i==2)
      hE1D[i]->GetYaxis()->SetTitle("E_{x} [E_{0}}]");
    
    
  }  
  
  // Chaning to user units: 
  if(opt.Contains("units") && n0) {
    for(Int_t i=0;i<Nfields;i++) {
      Int_t NbinsX = hE1D[i]->GetNbinsX();
      Float_t xMin = (skindepth/PUnits::mm) * hE1D[i]->GetXaxis()->GetXmin();
//.........这里部分代码省略.........