From Cecilia.Pizzolotto@ts.infn.it Fri Apr 10 09:30:53 2009
Date: Fri, 10 Apr 2009 09:30:35 +0200
From: Cecilia Pizzolotto <Cecilia.Pizzolotto@ts.infn.it>
To: Emiliano Mocchiutti <Emiliano.Mocchiutti@ts.infn.it>
Subject: spherical distribution

// Spherical distribution -- Test by Cecilia P march 2009 
-----------------------------///
-25 25 -25 25 108 108
void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, 
Double_t ymin, Double_t ymax,
                      Double_t zmin, Double_t zmax)
{
  Bool_t trkGood = kFALSE;
  Double_t theta = 999.;
  Double_t phi = 0.;
  Double_t x2,y2,x3,y3;

  //static const Double_t rad2deg = 57.2958;

  //static const Double_t s1_xmax=20.4, s1_ymax=16.5, s1_pz=102.8866;
  // calo cavity 8.07x6.57 with z2=71.6 z3=27.4 circa
  static const Double_t s2_xmax=7.8,  s2_ymax=6.0,  s2_z=73.489;
  static const Double_t s3_xmax=8.0,  s3_ymax=6.0,  s3_z=25.3159;

  // Generate random position unif. distr.
  // GenPosition(-25.,25., 25.,25., 108.0,108.0);
//   Double_t x0= fprim.fX0;
//   Double_t y0= fprim.fY0;
//   Double_t z0= fprim.fZ0;
  Double_t x0=  frnd->Uniform(xmin,xmax);
  Double_t y0=  frnd->Uniform(ymin,ymax);
  Double_t z0=  frnd->Uniform(zmin,zmax);
  Int_t posGood = 0;


  while (trkGood!=kTRUE)
    {
      // Generate phi theta
      theta=999.; // init
      while (theta>=0.82) //take only theta smaller than 30deg=0.52rad
    {
      Double_t xcos = sqrt( frnd->Uniform(0.,1.) );
      theta = acos(xcos);   //RAD
    }
      phi = frnd->Uniform(0.,2.*Pi());
 
      // Calculate xy at beginning/end of magnet cavity
      Double_t fact2 = (s2_z-z0)/cos(theta);
      x2 = x0 + fabs(fact2) *  sin(theta) * cos(phi);
      y2 = y0 + fabs(fact2) *  sin(theta) * sin(phi);
      Double_t fact3 = (s3_z-z0)/cos(theta);
      x3 = x0 + fabs(fact3) *  sin(theta) * cos(phi);
      y3 = y0 + fabs(fact3) *  sin(theta) * sin(phi);
    
      // Test condition on the direction
      if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) &&
       Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) {
    trkGood = kTRUE;
    //cout<<" x/y0= "<<x0<<" "<<y0<<"  x/y2= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<"  xy3= "<<
    //  fact3*sin(theta)*cos(phi)<<" "<<x3<<endl;
      }
    
      //if from current x0y0z0 the condition are not satysfied for any angle => def new start position
      posGood++;
      if ( posGood == 100 )
    {
      x0=  frnd->Uniform(xmin,xmax);
      y0=  frnd->Uniform(ymin,ymax);
      z0=  frnd->Uniform(zmin,zmax);
      posGood = 0;
    }
      
    }

  // Set direction and position:
  SetDirection(theta, phi);
  SetPosition(x0, y0, z0);
 
  return;
}
//--- end Test by Cecilia ------------------------------------///