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

void twoscales()
{
   TCanvas *c1 = new TCanvas("c1","hists with different scales",600,400);

   //create/fill draw h1
   gStyle->SetOptStat(kFALSE);
   TH1F *h1 = new TH1F("h1","my histogram",100,-3,3);
   Int_t i;
   for (i=0;i<10000;i++) h1->Fill(gRandom->Gaus(0,1));
   h1->Draw();
   c1->Update();

   //create hint1 filled with the bins integral of h1
   TH1F *hint1 = new TH1F("hint1","h1 bins integral",100,-3,3);
   Float_t sum = 0;
   for (i=1;i<=100;i++) {
      sum += h1->GetBinContent(i);
      hint1->SetBinContent(i,sum);
   }

   //scale hint1 to the pad coordinates
   Float_t rightmax = 1.1*hint1->GetMaximum();
   Float_t scale = gPad->GetUymax()/rightmax;
   hint1->SetLineColor(kRed);
   hint1->Scale(scale);
   hint1->Draw("same");

   //draw an axis on the right side
   TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
         gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"+L");
   axis->SetLineColor(kRed);
   axis->SetLabelColor(kRed);
   axis->Draw();
}

//==========================================
//==========================================
void plotPair(TH1F *h1,TH1F *h0 ) {
  h1->Draw();

  // edit fonts/sizes
  TAxis *ax =h1->GetYaxis(); 
  float ss=ax->GetTitleSize();
  //printf("ss=%f\n",ss);
  ax->SetTitleSize(2*ss);
  ax->SetTitleOffset(0.5);
  ax =h1->GetXaxis(); 
 ax->SetTitleSize(1.5*ss);
 ax->SetLabelSize(1.5*ss);
 ax->SetTitleOffset(0.7);
 

  // edit fonts/sizes DONE



  gPad->Update();
  //scale hint1 to the pad coordinates
  Float_t rightmax = 1.1*h0->GetMaximum();
  Float_t scale = gPad->GetUymax()/rightmax;
  h0->Scale(scale);
  h0->Draw("same");
   
  //draw an axis on the right side
  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
	      gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"-R");
  int col=h0->GetLineColor();
  axis->SetLineColor(col);
  
  axis->SetTextColor(col);
  axis->SetLabelColor(col);
  axis ->SetTitle("LCP yield");
 axis->SetTitleSize(2*ss);
 axis->SetTitleOffset(.5);

  axis->Draw();
  
  
  TPaveStats *st =( TPaveStats *)gPad->GetPrimitive("stats");
  st->SetX1NDC(0.35);
  st->SetX2NDC(0.5);
  st->SetY1NDC(0.7);
  st->SetY2NDC(1.);
}

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;
}

//.........这里部分代码省略.........
        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()));

        for(int i=0; i<numberBin; i++) {
            if(i==0) {

//.........这里部分代码省略.........
            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);
        ((TH1D*)(myEfficiencyHist->At(Index)))->Draw("same,ah");
    }
    myPad->Update();
    if(debug) cout << "Y: " << myPad->GetUymax() << endl;
    double YAxisMinValue=((TH1D*)(myEfficiencyHist->At(0)))->GetYaxis()->GetXmin();

//.........这里部分代码省略.........
  tt->SetBranchAddress("lumi", &lumi, &b_lumi);

  tt->SetBranchAddress("caloMjj", &caloMjj, &b_caloMjj);
  tt->SetBranchAddress("PFMjj", &PFMjj, &b_PFMjj);

  tt->SetBranchAddress("caloJet1Pt", &caloJet1Pt_, &b_caloJet1Pt);
  tt->SetBranchAddress("caloJet2Pt", &caloJet2Pt_, &b_caloJet2Pt);
  tt->SetBranchAddress("caloJet1Eta", &caloJet1Eta_, &b_caloJet1Eta);
  tt->SetBranchAddress("caloJet2Eta", &caloJet2Eta_, &b_caloJet2Eta);
  tt->SetBranchAddress("caloDeltaEta", &caloDeltaEta_, &b_caloDeltaEta);

  tt->SetBranchAddress("PFJet1Pt", &PFJet1Pt_, &b_PFJet1Pt);
  tt->SetBranchAddress("PFJet2Pt", &PFJet2Pt_, &b_PFJet2Pt);
  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)
    {

//.........这里部分代码省略.........
   InGE11_V.close();
   string CanvasName = "RunNumber"+std::to_string(RunNumber);
   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]);

//.........这里部分代码省略.........
  tt->SetBranchAddress("PFWJet2Eta", &PFWJet2Eta_, &b_PFWJet2Eta);
  tt->SetBranchAddress("PFWJet1Phi", &PFWJet1Phi_, &b_PFWJet1Phi);
  tt->SetBranchAddress("PFWJet2Phi", &PFWJet2Phi_, &b_PFWJet2Phi);
  tt->SetBranchAddress("PFWDeltaEta", &PFWDeltaEta_, &b_PFWDeltaEta);

  tt->SetBranchAddress("l1Jet1Pt", &l1Jet1Pt_, &b_l1Jet1Pt);
  tt->SetBranchAddress("l1Jet2Pt", &l1Jet2Pt_, &b_l1Jet2Pt);
  tt->SetBranchAddress("l1Jet1Eta", &l1Jet1Eta_, &b_l1Jet1Eta);
  tt->SetBranchAddress("l1Jet2Eta", &l1Jet2Eta_, &b_l1Jet2Eta);
  tt->SetBranchAddress("l1Jet1Phi", &l1Jet1Phi_, &b_l1Jet1Phi);
  tt->SetBranchAddress("l1Jet2Phi", &l1Jet2Phi_, &b_l1Jet2Phi);
  tt->SetBranchAddress("l1DeltaEta", &l1DeltaEta_, &b_l1DeltaEta);

  tt->SetBranchAddress("l1JetPt", &l1JetPt_, &b_l1JetPt);
  tt->SetBranchAddress("l1JetEta", &l1JetEta_, &b_l1JetEta);
  tt->SetBranchAddress("l1JetPhi", &l1JetPhi_, &b_l1JetPhi);

  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);

  //##############################################
  //##############################################

  //book graphs and plots
  TGraphErrors* totL1RateVsCut = new TGraphErrors();
  totL1RateVsCut->SetMinimum(0);
  TGraphErrors* pureL1RateVsCut = new TGraphErrors();

  TGraphErrors* totHltRateVsCut = new TGraphErrors();
  totHltRateVsCut->SetMinimum(0);
  TGraphErrors* pureHltRateVsCut = new TGraphErrors();

  //loops
  int bin = 0;
  for (int cut = 200; cut < 350; cut=cut+10)
    {
      std::cout << "analyzing point at " << cut << " GeV" << std::endl;

      int mjjHltPassed = 0;
      int excl_mjjHltPassed = 0;

      int mjjL1Passed = 0;
      int excl_mjjL1Passed = 0;

      int HTT240Passed = 0;

      for (Long64_t jentry=0; jentry<nentries;++jentry) 

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

    // Extract contours
    TCanvas* c = new TCanvas("c","Contour List",0,0,600,600);
    c->cd();

    // Potential
    TH2F *hV2Dc = (TH2F*) hV2D->Clone("hV2Dc");
    const Int_t Ncontours = 25;
    Double_t contours[Ncontours];
    for(Int_t i=0; i<Ncontours; i++) {
        contours[i] = i*(trapPotential/5.0) - trapPotential;
    }
    hV2Dc->SetContour(Ncontours, contours);
    hV2Dc->Draw("cont list");

    c->Update();
    TObjArray *contsV2D = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
    TClonesArray graphsV2D("TGraph",Ncontours);
    {
        Int_t ncontours = contsV2D->GetSize();
        TList* clist = NULL;
        Int_t nGraphs = 0;
        TGraph *gr = NULL;
        for(Int_t i = 0; i < ncontours; i++) {
            if(i==0) continue;

            clist = (TList*) contsV2D->At(i);

            for(Int_t j = 0 ; j < clist->GetSize(); j++) {
                gr = (TGraph*) clist->At(j);
                if(!gr) continue;

                gr->SetLineWidth(1);
                gr->SetLineColor(kGray+1);

                if( !((i)%5) ) {
                    gr->SetLineWidth(2);
                    gr->SetLineColor(kGray+2);
                }
                new(graphsV2D[nGraphs]) TGraph(*gr) ;
                nGraphs++;
            }
        }
    }

    // Ion probability
    hIonProb2D->GetZaxis()->SetRangeUser(0.00501,80);

    TH2F *hIonProb2Dc = (TH2F*) hIonProb2D->Clone("hIonProb2Dc");
    const Int_t NcontI = 4;
    Double_t contI[NcontI] = {0.01,0.1,1.0,10.0};
    hIonProb2Dc->SetContour(NcontI, contI);
    hIonProb2Dc->Draw("cont list");

    c->Update();
    TObjArray *contsI2D = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
    TClonesArray graphsI2D("TGraph",NcontI);
    {
        Int_t ncontours = contsI2D->GetSize();
        TList* clist = NULL;
        Int_t nGraphs = 0;
        TGraph *gr = NULL;
        for(Int_t i = 0; i < ncontours; i++) {
            clist = (TList*) contsI2D->At(i);

            for(Int_t j = 0 ; j < clist->GetSize(); j++) {

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

  if(opt.Contains("units") && n0) 
    sprintf(ctext,"Z = %5.1f mm", 1e3 * skindepth * Time);
  else
    sprintf(ctext,"T = %5.1f 1/#omega_{p}",Time);
  textTime->AddText(ctext);

  
  TPaveText *textRange = new TPaveText(0.13,0.87,0.38,0.92,"NDC");
  PlasmaGlob::SetPaveTextStyle(textRange,12); 
  textRange->SetTextColor(kGray+2);
  if(opt.Contains("units") && n0)
    sprintf(ctext,"%5.3f < y < %5.3f mm",ymin,ymax);
  else
    sprintf(ctext,"%5.3f < y < %5.3f c/#omega_{p}",ymin,ymax);
  textRange->AddText(ctext);

  // Actual Plotting!
  // ------------------------------------------------------------

  // Output file
  TString fOutName = Form("./%s/Plots/Density1D/Density1D",sim.Data());
  fOutName += Form("-%s_%i",sim.Data(),time);

  // Colors
  Int_t plasmaC = kGray+1;
  Int_t beamC   = kAzure-5;
  Int_t fieldC  = kOrange+10;
  Int_t fieldCb = kGray+1;

  C->cd(0);
  gPad->SetFrameLineWidth(2);  

  hDen1D[0]->SetLineColor(plasmaC);
  hDen1D[0]->SetLineWidth(1);
  hDen1D[0]->Draw("C");
  hDen1D[0]->GetYaxis()->CenterTitle();
  hDen1D[0]->GetXaxis()->CenterTitle();
  
  C->Update();

  TLine *line0 = new TLine(hDen1D[0]->GetXaxis()->GetXmin(),
			  (gPad->GetUymin()+gPad->GetUymax())/2.,
			  hDen1D[0]->GetXaxis()->GetXmax(),
			  (gPad->GetUymin()+gPad->GetUymax())/2.);
  line0->SetLineColor(kGray+1);
  line0->SetLineStyle(2);
  line0->Draw();

  Float_t rightmax = 2.5 * hDen1D[1]->GetMaximum();
  Float_t slope = (gPad->GetUymax() - gPad->GetUymin())/rightmax;

  for(Int_t i=0;i<hDen1D[1]->GetNbinsX();i++) {
    hDen1D[1]->SetBinContent(i+1,hDen1D[1]->GetBinContent(i+1)*slope + Min);
  }

  hDen1D[1]->SetLineWidth(2);
  hDen1D[1]->SetLineColor(beamC);
  hDen1D[1]->Draw("same C");
  // hTest->Draw("same");
 
  //draw an axis on the right side
  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(),
			    gPad->GetUymax(),0,rightmax,505,"+L");
  
  axis->SetLineWidth(1);

//.........这里部分代码省略.........
  
  // Canvas setup
  TCanvas *C = new TCanvas("C","Electric wakefield on axis",850,1000);
  C->Divide(1,2);

  // Draw objects
  TPaveText *textTime = new TPaveText(0.70,0.87,0.85,0.92,"NDC");
  PlasmaGlob::SetPaveTextStyle(textTime); 
  char ctext[128];
  if(opt.Contains("units") && n0) 
    sprintf(ctext,"Z = %5.1f mm", 1e3 * skindepth * Time);
  else
    sprintf(ctext,"T = %5.1f 1/#omega_{p}",Time);
  textTime->AddText(ctext);
  
  // Colors
  Int_t fieldC  = PlasmaGlob::fieldLine;
  Int_t phaseC  = kGray+1;
 
  // Actual Plotting!
  // ------------------------------------------------------------

  // More makeup
  C->cd(1);
  gPad->SetGridy(0);
  gPad->SetGridx(0);
  gPad->SetFrameLineWidth(2);  

  hE1D[0]->SetLineWidth(1);
  hE1D[0]->GetYaxis()->CenterTitle();
  hE1D[0]->GetXaxis()->CenterTitle();
  hE1D[0]->SetLineStyle(1);
  hE1D[0]->SetLineWidth(3);
  hE1D[0]->SetLineColor(fieldC);
  hE1D[0]->SetMarkerStyle(20);

  if(Nfields>1) {
    hE1D[1]->GetYaxis()->CenterTitle();
    hE1D[1]->GetXaxis()->CenterTitle();
    hE1D[1]->SetLineStyle(1);
    hE1D[1]->SetLineWidth(1);  
    hE1D[1]->SetLineColor(fieldC);
    hE1D[1]->SetMarkerStyle(24);
  }

  Float_t factor = 1.5;
  Float_t minimum = factor * hE1D[0]->GetMinimum();
  Float_t maximum = factor * hE1D[0]->GetMaximum();
    
  if(Nfields>1) {
    if(hE1D[1]->GetMaximum() > hE1D[0]->GetMaximum()) {
      maximum = factor * hE1D[1]->GetMaximum();
    } 
    
    if(hE1D[1]->GetMinimum() < hE1D[0]->GetMinimum()) {
    minimum = factor * hE1D[1]->GetMinimum();
    } 
  }
  
  if( maximum >= TMath::Abs(minimum)) minimum = -maximum;
  else maximum = - minimum;
  
  hE1D[0]->GetYaxis()->SetRangeUser(minimum,maximum);
  hE1D[0]->Draw("C");
  
  if(Nfields>1)