// calculate cross-section of total inelastic
// and 2(pi+pi-)pi0
void CS_DPM(double Plab=5.61){

  //1) total c.s 2(pi+pi-)pi0, based on parametisation from A.Dbeyssi note
  double a=26.14; double b=0.599; double c=0.26;
  double cs_pions = a*TMath::Exp(-b*Plab)+c/Plab;

  //2) inelastic = total - elastic, based on parametrisation from DPM (init.f)
  double cs_tot = 38.4 + 77.6*TMath::Power(Plab,-0.64)
    +0.26*TMath::Power(TMath::Log(Plab),2)-1.2*TMath::Log(Plab);
  double cs_el = 10.2 + 52.7*TMath::Power(Plab,-1.16)
    +0.125*TMath::Power(TMath::Log(Plab),2)-1.28*TMath::Log(Plab);

  double cs_inel=cs_tot-cs_el;
  //  cout<<"cs_tot = "<<cs_tot<<" cs_el = "<<cs_el<<endl;
  cout<<"cs_inel = "<<cs_inel<<" cs_pions = "<<cs_pions<<endl;
  //      <<" cs_pions/cs_inel = "<<cs_pions/cs_inel<<endl;

  //3) sig/inelastic
  double ratio_sig =  1e-5*(cs_pions/cs_inel);
  cout<<"ratio_sig = "<<ratio_sig<<endl;

  //4) Number of DPM events for bkg, if signal=1e4 events
  double Nbkg = 1e4/ratio_sig;
  cout<<"Number of bkg events (for 1e4 sig.events) = "<<Nbkg<<endl;
}