double myfunc(double xin[],double par[]){
  double d_2=par[4]/2.;
  double x=xin[0],y=xin[1];
  double b=par[1];
  double dxm=x-par[2]-d_2;
  double dxp=x-par[2]+d_2;
  double dym=y-par[3]-d_2;
  double dyp=y-par[3]+d_2;
  double f=atan(dxp*dyp/b/sqrt(b*b+dxp*dxp+dyp*dyp))
    -atan(dxm*dyp/b/sqrt(b*b+dxm*dxm+dyp*dyp))
    -atan(dxp*dym/b/sqrt(b*b+dxp*dxp+dym*dym))
    +atan(dxm*dym/b/sqrt(b*b+dxm*dxm+dym*dym));
  f*=par[0]/2./M_PI;
  return f;
}

void fit_shower(TH2F *Exy,bool insert=false){
  TF2 *f1=new TF2("f1",myfunc,-20,20,-20,20,5);
  f1->SetParameters(100,0.7,0,0);  
  Exy->GetXaxis()->SetRangeUser(-6,6);
  Exy->GetYaxis()->SetRangeUser(-6,6);
  f1->FixParameter(4,4.0);
  //f1->Draw();
  if (insert){ 
    f1->SetParameters(100,0.4,0,0);
    f1->FixParameter(4,2.0);
    Exy->GetXaxis()->SetRangeUser(-2.5,2.5);
    Exy->GetYaxis()->SetRangeUser(-2.5,2.5);
  }
  Exy->Fit("f1");
  TH2F *h=new TH2F("h",";#deltax[cm];#deltay[cm]",400,-20,20,400,-20,20);
  h->FillRandom("f1",100000);
  //h->Draw("colz");
}