{
TCanvas *c1 = new TCanvas("c1","",200,10,700,500);
    double myEne = 0.6;//T2K peak energy in GeV
    double myDmsq23 = 2.5e-3;//in eV^2
    double myDmsq21 = 7.6e-5;//in eV^2
    double myBaseL = 295;//T2K km
    double myTheta12 = asin(sqrt(0.8704))/2; //in radian 45*Pi/180
    double myTheta23 = asin(sqrt(0.5)); //in radian 45*Pi/180
    double myTheta13 = asin(sqrt(0.085))/2; //in radian 45*Pi/180
    double myDelta = 0;//-TMath::Pi()/2; //in radian
    double myDensity = 2.05e-4;
  
// the leading term
	double leading(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta);
{
    	double probmu2e = 4*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*pow(sin(1.27*dmsq_23*baseL/Ene),2);
        return probmu2e;

}  	    
  
    double myProbmu2e = leading(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta);
    cout<<"probably at the peak"<<myProbmu2e<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2e;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = leading(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2e = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2e->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
    pGraphMu2e->SetTitle("");    //set title for the graph
    pGraphMu2e->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
    pGraphMu2e->GetYaxis()->SetTitle("Oscillation Probability P(#nu_{#mu}#rightarrow #nu_{e})");
    pGraphMu2e->GetXaxis()->SetRangeUser(0.0, 2.0);
    pGraphMu2e->GetYaxis()->SetRangeUser(0, 0.06);
    pGraphMu2e->GetXaxis()->CenterTitle();
    pGraphMu2e->GetYaxis()->CenterTitle();
    pGraphMu2e->GetYaxis()->SetTitleOffset(1.2);
    pGraphMu2e->SetLineColor(2);
  //  pGraphMu2e->SetLineStyle(7);
    pGraphMu2e->Draw("AL");//draw graph
   // gPad->Print("prob_mat_mu2e_t2k.pdf");//pdf format
   

  
// the matter term
	double matter(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta, double density);
{
    	double probmu2eb = -8*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*(density*0.6/dmsq_23)*(2*pow(sin(theta_13), 2) - 1)*pow(sin(1.27*dmsq_23*baseL/Ene),2)
 + 8*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*(density*0.6*baseL/(4*Ene))*(2*pow(sin(theta_13), 2) - 1)*sin(1.27*dmsq_23*baseL/Ene)*cos(1.27*dmsq_23*baseL/Ene) ;
        return probmu2eb;

}  	    

    double myProbmu2eb = matter(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
    cout<<"probably at the peak"<<myProbmu2eb<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2eb;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = matter(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2eb = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2eb->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
 //   pGraphMu2eb->SetTitle("");    //set title for the graph
 //   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
 //   pGraphMu2eb->GetYaxis()->SetTitle("Osillation Probability #nu_{#mu}#rightarrow #nu_e");
    pGraphMu2eb->SetLineColor(3);
  //  pGraphMu2eb->SetLineStyle(7);
    pGraphMu2eb->Draw("same");//draw graph



// the CP-violating term
	double violation(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta, double density)
{
    	double probmu2ec = -8*sin(theta_12)*sin(theta_13)*sin(theta_23)*cos(theta_12)*pow(cos(theta_13),2)*cos(theta_23)*sin(delta)*sin(1.27*dmsq_21*baseL/Ene)*pow(sin(1.27*dmsq_23*baseL/Ene),2);
        return probmu2ec;

}  	    

    double myProbmu2ec = violation(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
    cout<<"probably at the peak"<<myProbmu2ec<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2ec;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = violation(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2ec = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2ec->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
 //   pGraphMu2eb->SetTitle("");    //set title for the graph
 //   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
 //   pGraphMu2eb->GetYaxis()->SetTitle("Osillation Probability #nu_{#mu}#rightarrow #nu_e");
    pGraphMu2ec->SetLineColor(4);
  //  pGraphMu2ec->SetLineStyle(7);
    pGraphMu2ec->Draw("same");//draw graph
   //save the graph
   // gPad->Print("prob_mat_mu2e_t2k.pdf");//pdf format



// the CP-conserving term
	double conservation(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta, double density)
{
    	double probmu2ed = 8*sin(theta_12)*sin(theta_13)*sin(theta_23)*pow(cos(theta_13),2)(cos(theta_12)*cos(theta_23)*cos(delta) - sin(theta_12)*sin(theta_13)*sin(theta_23))*sin(1.27*dmsq_21*baseL/Ene)*sin(1.27*dmsq_23*baseL/Ene)*cos(1.27*dmsq_23*baseL/Ene);
        return probmu2ed;

}  	    

    double myProbmu2ed = conservation(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
    cout<<"probably at the peak"<<myProbmu2ed<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2ed;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = conservation(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2ed = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2ed->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
 //   pGraphMu2eb->SetTitle("");    //set title for the graph
 //   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
 //   pGraphMu2eb->GetYaxis()->SetTitle("Osillation Probability #nu_{#mu}#rightarrow #nu_e");
    pGraphMu2ed->SetLineColor(5);
  //  pGraphMu2ed->SetLineStyle(7);
    pGraphMu2ed->Draw("same");//draw graph
   //save the graph
    //gPad->Print("prob_mat_mu2e_t2k.bmp");//pdf format




// the solar term
	double solar(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta, double density)
{
    	double probmu2ee = 4*pow(sin(theta_12),2)*pow(cos(theta_13),2)*(pow(cos(theta_12),2)*pow(cos(theta_23),2) + pow(sin(theta_12),2)*pow(sin(theta_13),2)*pow(sin(theta_23),2) - 2*sin(theta_12)*sin(theta_13)*sin(theta_23)*cos(theta_12)*cos(theta_23)*cos(delta))*pow(sin(1.27*dmsq_21*baseL/Ene),2);
        return probmu2ee;

}  	    

    double myProbmu2ee = solar(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
    cout<<"probably at the peak"<<myProbmu2ee<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2ee;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = solar(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2ee = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2ee->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
 //   pGraphMu2eb->SetTitle("");    //set title for the graph
 //   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
 //   pGraphMu2eb->GetYaxis()->SetTitle("Osillation Probability #nu_{#mu}#rightarrow #nu_e");
    pGraphMu2ee->SetLineColor(6);
 //   pGraphMu2ee->SetLineStyle(7);
    pGraphMu2ee->Draw("same");//draw graph
   //save the graph
    //gPad->Print("prob_mat_mu2e_t2k.bmp");//pdf format
  

// the total term
	double total(double Ene, double dmsq_23, double dmsq_21, double baseL, double theta_12, double theta_23, double theta_13, double delta, double density)
{
    	double probmu2ef = 4*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*pow(sin(1.27*dmsq_23*baseL/Ene),2) 

- 8*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*(density*0.6/dmsq_23)*(2*pow(sin(theta_13), 2) - 1)*pow(sin(1.27*dmsq_23*baseL/Ene),2)

+ 8*pow(sin(theta_13),2)*pow(sin(theta_23),2)*pow(cos(theta_13),2)*(density*0.6*baseL/(4*Ene))*(2*pow(sin(theta_13), 2) - 1)*sin(1.27*dmsq_23*baseL/Ene)*cos(1.27*dmsq_23*baseL/Ene) 

- 8*sin(theta_12)*sin(theta_13)*sin(theta_23)*cos(theta_12)*pow(cos(theta_13),2)*cos(theta_23)*sin(delta)*pow(sin(1.27*dmsq_23*baseL/Ene),2)*sin(1.27*dmsq_21*baseL/Ene)

+ 8*sin(theta_12)*sin(theta_13)*sin(theta_23)*pow(cos(theta_13),2)(cos(theta_12)*cos(theta_23)*cos(delta) - sin(theta_12)*sin(theta_13)*sin(theta_23))*sin(1.27*dmsq_21*baseL/Ene)*sin(1.27*dmsq_23*baseL/Ene)*cos(1.27*dmsq_23*baseL/Ene)

+ 4*pow(sin(theta_12),2) *pow(cos(theta_13),2)*(pow(cos(theta_12),2)*pow(cos(theta_23),2) + pow(sin(theta_12),2)*pow(sin(theta_13),2)*pow(sin(theta_23),2) - 2*sin(theta_12)*sin(theta_13)*sin(theta_23)*cos(theta_12)*cos(theta_23)*cos(delta))*pow(sin(1.27*dmsq_21*baseL/Ene),2);
        return probmu2ef;

}  	    

    double myProbmu2ef = total(myEne, myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
    cout<<"probably at the peak"<<myProbmu2ef<<endl;
    //set number of steps you want to calculate
    const Int_t nEnStep = 200;
    // set maximal energy to plot
    double EnergyMax = 2.0;
    
    //Graph of oscillation prob. as function of energy
    //Note:formula is not work at zero energy
    TGraph *pGraphMu2ef;
    Double_t pEne[nEnStep], pProb[nEnStep];

    for (Int_t istep=1; istep<=nEnStep; ++istep)
 {
        //energy value at each step
        pEne[istep-1] = EnergyMax*istep/nEnStep;
        //calculated osc. prob. at each step
        pProb[istep-1] = total(pEne[istep-1],myDmsq23,myDmsq21, myBaseL,myTheta12,myTheta23,myTheta13,myDelta, myDensity);
        cout<<"step "<<istep<<" Energy "<<pEne[istep-1]<<" Prob. "<<pProb[istep-1]<<endl;
  }
    //create new Graph based on the array
    pGraphMu2ef = new TGraph(nEnStep,pEne,pProb);
    //set the width for graph line
    pGraphMu2ef->SetLineWidth(2);
    
    //create new Canvas
   // new TCanvas;
    
 //   pGraphMu2eb->SetTitle("");    //set title for the graph
 //   pGraphMu2eb->GetXaxis()->SetTitle("Neutrino Energy [GeV]");
 //   pGraphMu2eb->GetYaxis()->SetTitle("Osillation Probability #nu_{#mu}#rightarrow #nu_e");
    pGraphMu2ef->SetLineColor(1);
    pGraphMu2ef->Draw("same");//draw graph

 /* Draw horizontal lines   */
    c1->Update();
    TLine *l1 = new TLine(c1->GetUxmin(), 0.0, c1->GetUxmax(), 0.0);
    l1->SetLineColor(1);
    l1->Draw("L SAME");

    TLegend * pleg = new TLegend(0.7,0.65,0.85,0.87);
    pleg -> AddEntry( pGraphMu2ef , "Total" ,"l");
    pleg -> AddEntry( pGraphMu2e , "Leading" ,"l");
    pleg -> AddEntry( pGraphMu2eb, "Matter effect" ,"l");
    pleg -> AddEntry( pGraphMu2ec , "CP-violating" ,"l");
    pleg -> AddEntry( pGraphMu2ed , "CP-conserving" ,"l");
    pleg -> AddEntry( pGraphMu2ee , "Solar" ,"l");

    pleg -> SetFillColor(0);
    pleg -> SetBorderSize(0);
    pleg ->SetTextSize(14);
    pleg ->SetTextFont(43);
    pleg -> Draw();


 gStyle->SetOptStat(0);
leg = new TLegend(0.7,0.25,0.85,0.35);
leg->SetTextSize(0.03);
leg->SetTextColor(1);
leg->SetFillColor(0);
leg->SetBorderSize(0);
leg ->SetTextSize(14);
leg ->SetTextFont(43);
leg->AddEntry(pGraphMu2e,"L = 295km","");
leg->AddEntry(pGraphMu2e,"#delta_{CP} = -#pi/2","");
leg->Draw();
   //save the graph
    gPad->Print("prob_vac_mu2e_t2k.pdf");//pdf format

 
}