void NeutrinoMass()
{
  // *** inputs ***
  // particle massees and their uncertainties in MeV
  double mass_pi = 139.57039,   merr_pi = 0.00018;
  double mass_mu = 105.6583755, merr_mu = 0.0000023;
  // neutrino mass is set manually. You can change this.
  double mass_nu = 0.2;
  
  double p_pi = 100.0;   // magnitude of the sapace momentum of pion
  double sigmaP = 60e-3; // momentum resolution of pion and muon
  int nEvent = 500;      // number of events required

  // *** calculation ***
  // 4-vectors 
  TLorentzVector pi, mu, nu;
  
  // book histograms
  TCanvas *c1 = new TCanvas("c1","c1",800,600);
  c1->Divide(2,2);
 
  // set random number generator whose seed is taken from timer
  TRandom3 rnd;
  rnd.SetSeed(time(NULL));

  // variables
  double p,px,py,pz,ene,theta,phi,s;
  double const twoPi = 2.0*3.1415926;

  while(1){
   TH1F *hp_pi = new TH1F("hp_pi","pion momentum;p [MeV];Entries",100,99,101);
   TH1F *hp_mu = new TH1F("hp_mu","muon momentum;p [MeV];Entries",100,0,150);
   TH1F *hp_nu = new TH1F("hp_nu","neutrino momentum; p [MeV];Entries",100,0,150);
   TH1F *hm_nu = new TH1F("hm_nu","neutrino mass-aquared;#m^2 [#MeV^2]",50,-10,10);

   for(int i=0; i<nEvent; i++){
     // set lab 4-vector of pi
     pi.SetXYZM(0.0, 0.0, p_pi, mass_pi);

     // cm frame of pion
     ene = (mass_pi*mass_pi+mass_mu*mass_mu-mass_nu*mass_nu)/(2*mass_pi);
     p   = sqrt(ene*ene-mass_mu*mass_mu);
     theta = acos(2.0*rnd.Uniform()-1.0);
     phi= twoPi*rnd.Uniform();
     px = p*sin(theta)*cos(phi);
     py = p*sin(theta)*sin(phi);
     pz = p*cos(theta);
     mu.SetXYZM( px, py, pz,mass_mu);
     nu.SetXYZM(-px,-py,-pz,mass_nu);

     // boost mu and  nu to lab frame 
     mu.Boost(pi.Px()/pi.E(), pi.Py()/pi.E(), pi.Pz()/pi.E());
     nu.Boost(pi.Px()/pi.E(), pi.Py()/pi.E(), pi.Pz()/pi.E());

     // now simulate detector effects by smearing momenta
      s = rnd.Gaus(mu.P(),sigmaP) / mu.P();
     px = mu.Px()*s;
     py = mu.Py()*s;
     pz = mu.Pz()*s;
     mu.SetXYZM(px,py,pz,mass_mu);

      s = rnd.Gaus(pi.P(),sigmaP) / pi.P();
     px = pi.Px()*s;
     py = pi.Py()*s;
     pz = pi.Pz()*s;
     pi.SetXYZM(px,py,pz,mass_pi);
	
     // construct 4-vector of neutrino
     nu = pi - mu;
	
	 // fill histograms
     hm_nu->Fill(nu.M2());
	 hp_pi->Fill(pi.P());
	 hp_mu->Fill(mu.P());
	 hp_nu->Fill(nu.P());
  }
  if(hm_nu->GetMean() > 0.0) break;
 }
 c1->cd(1); hm_nu->Draw();
 c1->cd(2); hp_pi->Draw();
 c1->cd(3); hp_mu->Draw();
 c1->cd(4); hp_nu->Draw();

 // setting upper limit
 
 double CL = 0.90;
 double m2 = hm_nu->GetMean();
 double s2 = hm_nu->GetRMS();
 double m  = sqrt(m2);
 double s  = s2/(2*m);
 double f  = ROOT::Math::normal_quantile(CL,1);
 double e  = f*s/sqrt(nEvent);
 double mup= m + e;
 cout << "Mean of neutrino mass-square = " << m2 << endl;
 cout << "RMS  of neutrino mass-square = " << s2 << endl;
 cout << "Mass neutrino < " << setprecision(2) 
      << mup <<" MeV " << " CL = " << int(CL*100) <<"%"<< endl;
}