{

  const double Dist = 2.46;

  const double P1 = 0.;               // m      TOF1 position

  const double P2 = P1 + Dist;        // m      TOF2 position

  const double tau_pi = 2.6e-8;       // s      pion lifetime

  const double E_pi_beam = 270.;      // MeV    beam energy

  const double m_pi = 139.57;         // MeV    pion mass

  const double m_mu = 105.66;         // MeV    muon mass

  //const double Ek_mu_cm = 4.12;     // MeV

  const double E_mu_cm = 109.78;      // MeV    muon energy in the c.m. system

  //const double P_mu_cm = 29.79;     // MeV/c

  const double P_mu_cm = sqrt(E_mu_cm*E_mu_cm - m_mu*m_mu);  // MeV/c

  TCanvas *c1 = new TCanvas("MC","MC",1200,800);

  TH1F *h1 = new TH1F("pi_time","time",200,7.,12.);   // time of flight for pions (no decay)

  TH1F *h2 = new TH1F("pimu_time","time",200,7.,12.); // time of flight for pions->muons (decay between TOF1 and TOF2)

  TH1F *h2a = new TH1F("dtimea","dtimea",200,7.,12.); 

  TH1F *h3 = new TH1F("tof","tof  MC",200,7.,12.);    // time of flight for all events

  h3->GetXaxis()->SetTitle("tof ns");

  TH1F *h4 = new TH1F("betaMu","betaMu",100,.6,1.);   // beta of the muons

  TH1F *h5 = new TH1F("betaPi","betaPi",100,.6,1.);   // beta of the pions

  TH1F *h6 = new TH1F("ThetaMu","ThetaMu",200,0,.6);//TMath::Pi());

  TH1F *h7 = new TH1F("dbeta","dbeta",100,-.5,.5);

  double P_pi_beam = sqrt(E_pi_beam*E_pi_beam - m_pi*m_pi);

  cout<<"E_pi = "<<E_pi_beam<<"   P_pi = "<<P_pi_beam<<"  beta_pi = "<< P_pi_beam/E_pi_beam <<endl;

  TRandom* Rand = new TRandom();

  Rand->SetSeed(0);

  int i = 0;

  while( i<1e4 )

  {

    double E_pi = Rand->Gaus(E_pi_beam, E_pi_beam*0.03);  // initial energy of the pion

    double Pz_pi = sqrt(E_pi*E_pi - m_pi*m_pi);

    TLorentzVector pi(0, 0, Pz_pi, E_pi);

    double beta_pi = pi.Beta();                           // pion beta

    double gamma_pi = pi.Gamma();                         // pion gamma

    // pion decay simulation

    double u = Rand->Uniform(0,1);

    double t_pi = -gamma_pi*tau_pi*log(u);

    double l_pi = t_pi*beta_pi*TMath::C();

    double l_mu = P2 - l_pi;

    double CosThCM_mu = Rand->Uniform(-1,1);

    double SinThCM_mu = sqrt(1 - CosThCM_mu*CosThCM_mu);

    double PhiCM_mu = Rand->Uniform(-TMath::Pi(),TMath::Pi());

    double px = P_mu_cm*SinThCM_mu*cos(PhiCM_mu);

    double py = P_mu_cm*SinThCM_mu*sin(PhiCM_mu);

    double pz = P_mu_cm*CosThCM_mu;

    TLorentzVector mu(px, py, pz, E_mu_cm);               // muon in the c.m. system

    mu.Boost(0,0,beta_pi);                                // muon in the lab. system

    double beta_mu = mu.Beta();

    double theta_mu = mu.Theta();

    h4->Fill(beta_mu);

    h5->Fill(beta_pi);

    h6->Fill(theta_mu);

    h7->Fill(beta_mu - beta_pi);

    if( l_pi < P1)

    {

      h2a->Fill( Dist/( beta_mu*TMath::C() )*1e9 );

      //h3->Fill( Dist/( beta_mu*TMath::C() )*1e9 );	  

    }

    else if(l_pi > P1 && l_pi < P2)  // decay between TOF1 and TOF2

    {

      double t1_pi = (l_pi - P1)/( beta_pi*TMath::C() )*1e9;

      double t2_mu = (P2 - l_pi)/( beta_mu*TMath::C() )*1e9;

      double tof = Rand->Gaus(t1_pi + t2_mu, 7e-2);  // measured time of flight [ns]

      h2->Fill(tof);

      h3->Fill(tof);

    }

    else if(l_pi > P2)   //No decay

    {

      double tof = Rand->Gaus( Dist/( beta_pi*TMath::C() )*1e9, 7e-2);

      h1->Fill(tof);

      h3->Fill(tof);

    }

    i++;

  }

  c1->Divide(2,0);

  /*

  c1->cd(4);

  h7->GetXaxis()->SetTitle("#beta_{#mu} -#beta_{#pi} ");

  h7->Draw();

  */

  c1->cd(2);

  h3->SetFillColor(kGray);

  h3->Draw();

  h2->SetLineColor(kBlue);

  h2->Draw("same");

  h2a->SetLineColor(kGreen);

  //h2a->Draw("same"); 

  h1->SetLineColor(kRed);

  h1->Draw("same"); 

  c1->cd(1);

  h5->SetLineColor(kBlue);

  h5->GetXaxis()->SetTitle("#beta");

  h5->Draw();

  h5->SetLineColor(kRed);

  h4->Draw("sames");

/*

  c1->cd(3);

  h6->Draw();

*/

}