#include "PamVMC2PiAux.h"

void PamVMC2PiAux::gsurf(){

  Double_t th1 = DegToRad()*17.;
  Double_t th2 = DegToRad()*20.;
  Double_t z1 = 87.4389954 + 0.150935784  +  0.4*Sin(th1) + 18.15*Cos(th1)/2.;
  Double_t z2 = 87.4389954 + 0.00949157123 + 0.4*Sin(th2) + 18.15*Cos(th2)/2.;
  Double_t xcard = GPAXCARD - 0.4*Cos(th2) + 18.15*Sin(th2)/2.;
  Double_t ycard = GPAYCARD - 0.4*Cos(th1) + ( z2-(87.4389954 + 0.150935784 + 0.4*Sin(th1)) )*Sin(th1);
  Double_t th = ATan(Sqrt(Power((xcard-XS2),2.) + Power((ycard-YS2),2.))/(z2-ZS2MAX));
  Double_t ph = ATan((ycard-YS2)/(xcard-XS2));

  if (TWOPIDEBUG){
    cout<<"th1="<<th1<<" th2="<<th2<<" z1="<<z1<<" z2="<<z2<<endl;
    cout<<"xcard="<<xcard<<" ycard="<<ycard<<" th="<<th<<" ph="<<ph<<endl; 
  } 

  Double_t b = 2.*(XS2*Sin(th)*Cos(ph)+YS2*Sin(ph)*Sin(th) + ZS2MAX*Cos(th)+SHELL3RZ*Cos(th));
  Double_t c = XS2*XS2 + YS2*YS2 + (ZS2MAX+SHELL3RZ)*(ZS2MAX+SHELL3RZ) - SHELL3RR*SHELL3RR;
  Double_t l = (-b+Sqrt(b*b-4.*c))/2.;

  Double_t ZMINS = Sqrt(Power(SHELL3RR,2.) - Power(RCY,2.)) - SHELL3RZ;
  fzmins = ZMINS;
    

  //  if (TWOPIDEBUG) cout<<"ZMINS="<<ZMINS<<" RCY="<<RCY<<endl;

  Double_t h1 = (SHELL3RR-SHELL3RZ-ZMINS);
  Double_t SPHAR = TwoPi()*h1*SHELL3RR;

  
  Double_t ta = SPHAR;
  Double_t THMAX = ATan(RCY/(ZMINS+SHELL3RZ));
  fthmax = THMAX;
  if (TWOPIDEBUG){
    cout<<"Spherical surface:"<<endl;
    cout<<"  Zcenter: "<<-SHELL3RZ<<" cm, Radius: "<<SHELL3RR
	<<" cm, MaxAng: "<<THMAX*RadToDeg()<<endl;
    cout<<"    Spherical surface area: "<<SPHAR<<" cm^2"<<endl;
  }

  Double_t tgy = (ZS1-z2)/(YS1+ycard);
  Double_t ry = tgy/(RCY+YS1);
  Double_t h2 = (ZS2MAX+l*Cos(th)-ZS1);
  h2 = ry + h2;
  Double_t ZMINC = ZMINS - h2;
  fzminc = ZMINC;

  Double_t CYLAR = TwoPi()*h2*RCY;

  if (TWOPIDEBUG){
    cout<<"Cylindrical surface:"<<endl;
    cout<<"  Height of cylinder: "<<h2<<" cm, Radius: "<<RCY<<" cm"<<endl;
    cout<<"  Zmin: "<<fzminc<<" cm, Zmax: "<<fzmins<<" cm"<<endl;
    cout<<"    Cylindrical surface area: "<<CYLAR<<" cm^2"<<endl;
  }
 
  ta += CYLAR;

  if(TWOPIDEBUG) cout<<"Total generation surface area = "<<ta<<" cm^2"<<endl;

  SPHAR/=ta;
  CYLAR/=ta;

  if(TWOPIDEBUG) cout<<"SPHAR,CYLAR="<<SPHAR<<","<<CYLAR<<endl;

  fsphar = SPHAR;
  farea = ta;
}

void PamVMC2PiAux::GenPosAng(TVector3& pos, TVector2& ang){

  Double_t phi = TwoPi()*frandom->Uniform();
  Double_t the;
  Double_t ct, st, cf, sf, xbeam, ybeam, zbeam;


  //POSITIONS
 if( (frandom->Uniform()) > fsphar ){
    //generation on the lateral cylindrical surface
    the = Pi()/2.;
    ct = Cos(the);
    st = Sin(the);
    cf = Cos(phi);
    sf = Sin(phi);
    xbeam = RCY * cf;
    ybeam = RCY * sf;
    zbeam = -(fzminc+frandom->Uniform()*(fzmins-fzminc));
  } else {
    //generation from top spherical surface
    the = ASin( Sqrt(frandom->Uniform())*Sin(fthmax) );  //CORRECT!
    ct = Cos(the);
    st = Sin(the);
    cf = Cos(phi);
    sf = Sin(phi);
    xbeam = SHELL3RR * st * cf;
    ybeam = SHELL3RR * st * sf;
    zbeam = -SHELL3RR * ct + SHELL3RZ;
   }
 
  pos.SetXYZ(xbeam,ybeam,-zbeam);

  //DIRECTION COSINES
  // local angles. COS^2 generation (theta_min=0, theta_max=90)
  Double_t  lthe = ASin(Sqrt(frandom->Uniform()));
  Double_t  lphi = TwoPi() * frandom->Uniform();

  // local director cosines

  Double_t aaa = Sin(lthe) * Cos(lphi);
  Double_t bbb = Sin(lthe) * Sin(lphi);
  Double_t ccc = Cos(lthe);

  // global director cosines
  //MY CORRECT COSINES


  Double_t ubeam = aaa*cf*ct - bbb*sf - ccc*cf*st; //cos(px)
  Double_t vbeam = aaa*sf*ct + bbb*cf - ccc*sf*st; //cos(py)
  Double_t wbeam = aaa*st + ccc*ct; //cos(pz)
  
  Double_t theta_dir = ACos(wbeam);
  Double_t phi_dir = ATan2(vbeam,ubeam);
  if(phi_dir<0.) phi_dir+=TwoPi();
  
  ang.Set(theta_dir,phi_dir);
}