Double_t e_p2ep(Double_t *x, Double_t *par) { // par[0]=norm, par[2]=peak, par[2]=width, par[3]=tail, par[4]=exponent // good starting parameters: peak=66.2, width=3.49, tail=1.19, exp=3.1 (no cuts), exp=7.3 (CsIOut<5) Double_t yield=0; TRandom r; r.SetSeed(); for (Int_t i=0;i<10000;i++) { Double_t t = (x[0]-par[1])/par[2]; if (r.Uniform(0,1)<0.01) t=(x[0]*(1-r.Uniform(0,1))-par[1])/par[2]; if (par[3] < 0) t = -t; Double_t absAlpha = fabs(par[3]); if (t >= -absAlpha) { yield+=par[0]*exp(-0.5*t*t); } else { Double_t a = TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha); Double_t b= par[4]/absAlpha - absAlpha; yield+=par[0]*a/TMath::Power(b - t, par[4]); } } // for return yield/10000; }