// AvdS/AP function, 5 parameters

Double_t e_michp(Double_t *x, Double_t *par) {
  Double_t result = 0; 

  TRandom r;
  r.SetSeed();
  
  for (Int_t i=0;i<10000;i++) {
     Double_t t=x[0];
     Double_t a=r.Uniform(0,1);
     Double_t b=r.Uniform(0,1);
     Double_t c=r.Uniform(0,1);
       if (a<0.01) t=x[0]*(1-b);
 //without gaussian edges (gives errors)
  result+=fabs(TMath::Erf(t)*(t-par[1])**2*(3-2*t/par[2])*(1-TMath::Erf((t-par[3])/par[4])));
 //result=fabs(TMath::Erf(x[0])*(x[0]-par[1])**2*(3-2*x[0]/par[2])*(1-TMath::Erf((x[0]-par[3])/par[4]))); 

 //with Gaussians at edges
    // if (x[0]<=2) 
    //result=exp(x[0]);
    //if (x[0]>2 && x[0]<=par[3])
    //result=fabs((x[0]-par[1])**2*(3-2*x[0]/par[2])*(1-TMath::Erf((x[0]-par[3])/par[4])));
    //if (x[0]>par[3]) 
    //result=1000*exp(par[3]-x[0]);

 //original e_mich

    // if (x[0]<=2)
    //result=exp(x[0]);
    //if (x[0]>2 && x[0]<=par[3])
    //result=abs((x[0]-par[1])**2*(3-2*x[0]/par[2])*(1-TMath::Erf((x[0]-par[3])/par[4])));
    //if (x[0]>par[3])
    //result=1000*exp(par[3]-x[0]);

} // end for

 result=par[0]*result/10000;

 return result;
}