PamVMC_update/src/PamVMCPrimaryGenerator.cxx

// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 


#include <TVirtualMC.h>
#include <TVirtualMCStack.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TVector3.h>
#include <TMath.h>
#include <Riostream.h>


#include "PamVMCPrimaryGenerator.h"

using namespace TMath;

ClassImp(PamVMCPrimary)


PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
{
  a.fID=b.fID;
  a.fPDG=b.fPDG;
  a.fX0=b.fX0;
  a.fY0=b.fY0;
  a.fZ0=b.fZ0;
  a.fTHETA=b.fTHETA;
  a.fPHI=b.fPHI;
  a.fP0=b.fP0;
  a.fGOOD=b.fGOOD;
  a.fXTOL=b.fXTOL;
  a.fYTOL=b.fYTOL;
  for (Int_t j = 0; j<14; j++){
    a.fCR_X[j] = b.fCR_X[j];    
    a.fCR_Y[j] = b.fCR_Y[j];            
  }         

  return a;
}


ClassImp(PamVMCPrimaryGenerator)

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) 
  : TObject(),
    fStack(stack),
    fevno(0),
    fngood(0),
    fmass(0.),
    fcharge(0.),
    fxprev(0.),
    fyprev(0.),
    fzprev(200.),
    frandom(0)
{
  //gObjectTable->Add(this);

// Standard constructor

  ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
  ftheta->SetNpx(1000);

  fprimColl = new TClonesArray("PamVMCPrimary");
  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV

  f2PiAux = new PamVMC2PiAux();
  fpgf = new PamVMCPrimaryGF();
} 

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
  : TObject(),
    fStack(0),
    fevno(0),
    fngood(0),
    fmass(0.),
    fcharge(0.),
    fxprev(0.),
    fyprev(0.),
    fzprev(200.),
    fprimColl(0),
    frandom(0)
{    
 // Default constructor
  //Default primary proton

  //gObjectTable->Add(this);

  ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
  ftheta->SetNpx(1000);

  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV

  f2PiAux = new PamVMC2PiAux();
  fpgf = new PamVMCPrimaryGF();
}

PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
{
// Destructor
  delete fpgf;
  delete f2PiAux;
  delete ftheta;
  delete fprimColl;
  //gObjectTable->Remove(this);
}


// private methods


void PamVMCPrimaryGenerator::GeneratePrimary()
{    
// Add one primary particle to the user stack (derived from TVirtualMCStack).
  
  // Track ID (filled by stack)
  Int_t ntr;
 
  // Option: to be tracked
  Int_t toBeDone = 1; 
 
  // Particle type
  Int_t pdg  = fprim.fPDG;
  
  Double_t fvx, fvy, fvz;
  fvx=fprim.fX0;
  fvy=fprim.fY0;
  fvz=fprim.fZ0;

  // Position
  
  Double_t tof = 0.;

  // Energy (in GeV)
  //printf("generateprimary check fprimP0 = %f\n",fprim.fP0);
  Double_t kinEnergy = MomentumToKinE(fprim.fP0);     
  Double_t e  = fmass + kinEnergy;
 
  // Particle momentum
  Double_t  px, py, pz;
   
  px = TMath::Abs(fprim.fP0)*Sin(fprim.fTHETA)*Cos(fprim.fPHI);
  py = TMath::Abs(fprim.fP0)*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
  pz = -TMath::Abs(fprim.fP0)*Cos(fprim.fTHETA);
 
  // Polarization
  TVector3 polar;

  // Add particle to stack 
  fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof, 
                   polar.X(), polar.Y(), polar.Z(), 
                   kPPrimary, ntr, 1., 0);

  PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
 
  *pc = fprim;
}


void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
  fprim.fPDG=pdg;
  //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
  fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
  fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
}

void PamVMCPrimaryGenerator::SetMomentum(
                              Double_t px, Double_t py, Double_t pz) 
{
  fprim.fP0= Sqrt(px*px+py*py+pz*pz);
  fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz);
  fprim.fPHI=ATan(py/px);
}
                                     
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
{
  if(isEnergy) {
    fprim.fP0=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
  } else{
    fprim.fP0=frandom->Uniform(PEmin,PEmax);
  }

}

void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
{
  Double_t alpha = 1.+gamma; //integral spectral index
  if(alpha==0.){
    fprim.fP0=Exp(Log(PEmin)+frandom->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
  } else {
    if(PEmin==0.) PEmin=1.E-10;
    fprim.fP0=Power((frandom->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
  }
  //  cout<<"GenSpe fprim.fP0= "<<fprim.fP0<<endl;
  if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);

}


//Cecilia Pizzolotto: powerlaw spectrum 3 with the shape
// J(E) = 0.5*(E + b * exp(-c * sqrt(E)))^-a
// between PEmin and PEmax and with the input parameters a,b,c.
// Valeria di Felice fits parameter values are:
// protons:   a,b,c= 2.70, 2.15, 0.21
// electrons: a,b,c= 0.0638, 1.248e-16, -38.248  
void PamVMCPrimaryGenerator::GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c)
{

  Bool_t found=0;
  Double_t funct_min, funct_max;
  funct_max = function3par(PEmin,a,b,c);
  funct_min = function3par(PEmax,a,b,c);
  //
  Double_t wurfP;
  Double_t wurfy ;
  //printf("in genspe3par^^^^%f ^^%f ^^^^^^^^%f ^^^^^%f^^^^^^^\n",PEmin,PEmax,funct_min,funct_max);
  //printf("in par^^ %f  %f  %f \n",a,b,c);
  while( found==0 )
    {
      wurfP = frandom->Uniform(PEmin,PEmax);
      wurfy = frandom->Uniform(funct_min,funct_max);
      if( wurfy<(function3par(wurfP,a,b,c) )) 
        {
          // this is ok!
          fprim.fP0=wurfP;
          found=1;
        }
    }
  //printf("exit+++++++++++++++++++  %f   %f \n",wurfP,fprim.fP0);
}


// cecilia pizzolotto
void PamVMCPrimaryGenerator::GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
{
  // Generates a flat spectrum from PEmin to PElim. Then a power law
  Double_t PElim = 1.;
  //Double_t alpha = 1.+gamma; //integral spectral index

  Bool_t   okflag=0.;
  Double_t throw_x =0.;
  Double_t throw_y =0.;
  
  while(okflag==0)
    {
      throw_x=frandom->Uniform(PEmin,PEmax);
      //        cout<<"        x "<<throw_x<<endl;
      if(throw_x<=PElim)
        {
          okflag=1.;
        }
      else
        {
          throw_y=frandom->Uniform(0.,1.);
          if( throw_y<(1*pow(throw_x,gamma)))
            {
              okflag=1.;
            }
        }
    }
  fprim.fP0=throw_x;
  //h->Fill(fprimf.P0);
  okflag=0.; // reset
  
  if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);

}

// Spherical distribution -- Test by Cecilia P july 2009 ----
// flusso isotropo su 2pi
void PamVMCPrimaryGenerator::GenSphericalPhiThe()
{
  // Generate phi theta
  Double_t theta=0.;
  Double_t phi=0.;
  
  Double_t xcos = sqrt( frandom->Uniform(0.,1.) );
  theta = acos(xcos);   //RAD
  
  phi = frandom->Uniform(0.,2.*Pi());
  
  SetDirection(theta, phi);
  return;
}


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;
  Double_t x0,y0,z0;

  //static const Double_t rad2deg = 57.2958;
  // S21 and S31 position/size taken as reference (z on top of det)
  // constraint: must pass throuth these planes
  static const Double_t s2_xmax=9.05, s2_ymax=7.55, s2_z=73.439; // z on top of det
  static const Double_t s3_xmax=9.05, s3_ymax=7.55, s3_z=26.093; // z on top of det

  //Double_t thetamax=3.14;
  //thetamax = atan((xmax+s3_xmax)/(zmax-s3_z));
  //cout<<" Quanto è il theta max? "<<thetamax<<" in deg "<<thetamax*(90./Pi())<<endl;

  while (trkGood!=kTRUE)
    {
       x0= frandom->Uniform(xmin,xmax);
       y0= frandom->Uniform(ymin,ymax);
       z0= frandom->Uniform(zmin,zmax);
  
      // Generate phi theta
      theta=999.; // init
      while (theta>=0.65) // take only theta smaller than 37deg=0.65rad
        {
          Double_t xcos = sqrt( frandom->Uniform(0.,1.) );
          theta = acos(xcos);   //RAD
        }
      phi = frandom->Uniform(0.,2.*Pi());
 
      // Calculate xy at the constraint
      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);
         
      //cout<<" x/y0= "<<x0<<" "<<y0<<"  x/y2= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<"  xy3= "<<
      //  fact3*sin(theta)*cos(phi)<<" "<<x3<<"   phi/the "<<phi*(90./Pi())<<" "<<theta*(90./Pi())<<endl;

      // 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;
      }      
    } 

  // Set direction and position:
  SetDirection(theta, phi);
  SetPosition(x0, y0, z0);
  
  return;
}
1	formato	1.1	// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03
2
3
4			#include <TVirtualMC.h>
5			#include <TVirtualMCStack.h>
6			#include <TPDGCode.h>
7			#include <TDatabasePDG.h>
8			#include <TParticlePDG.h>
9			#include <TVector3.h>
10			#include <TMath.h>
11			#include <Riostream.h>
12
13
14
15			#include "PamVMCPrimaryGenerator.h"
16
17			using namespace TMath;
18
19			ClassImp(PamVMCPrimary)
20
21
22			PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
23			{
24			a.fID=b.fID;
25			a.fPDG=b.fPDG;
26			a.fX0=b.fX0;
27			a.fY0=b.fY0;
28			a.fZ0=b.fZ0;
29			a.fTHETA=b.fTHETA;
30			a.fPHI=b.fPHI;
31			a.fP0=b.fP0;
32			a.fGOOD=b.fGOOD;
33			a.fXTOL=b.fXTOL;
34			a.fYTOL=b.fYTOL;
35			for (Int_t j = 0; j<14; j++){
36			a.fCR_X[j] = b.fCR_X[j];
37			a.fCR_Y[j] = b.fCR_Y[j];
38			}
39
40			return a;
41			}
42
43
44
45			ClassImp(PamVMCPrimaryGenerator)
46
47			PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)
48			: TObject(),
49			fStack(stack),
50			fevno(0),
51			fngood(0),
52			fmass(0.),
53			fcharge(0.),
54			fxprev(0.),
55			fyprev(0.),
56			fzprev(200.),
57			frandom(0)
58			{
59			//gObjectTable->Add(this);
60
61			// Standard constructor
62
63			ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
64			ftheta->SetNpx(1000);
65
66			fprimColl = new TClonesArray("PamVMCPrimary");
67			fprim.fPDG=kProton;
68			fprim.fX0=1.;
69			fprim.fY0=1.;
70			fprim.fZ0=130.;
71			fprim.fTHETA=0.;
72			fprim.fPHI=0.;
73			fprim.fP0=1.; //1GV
74
75			f2PiAux = new PamVMC2PiAux();
76			fpgf = new PamVMCPrimaryGF();
77			}
78
79			PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
80			: TObject(),
81			fStack(0),
82			fevno(0),
83			fngood(0),
84			fmass(0.),
85			fcharge(0.),
86			fxprev(0.),
87			fyprev(0.),
88			fzprev(200.),
89			fprimColl(0),
90			frandom(0)
91			{
92			// Default constructor
93			//Default primary proton
94
95			//gObjectTable->Add(this);
96
97			ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
98			ftheta->SetNpx(1000);
99
100			fprim.fPDG=kProton;
101			fprim.fX0=1.;
102			fprim.fY0=1.;
103			fprim.fZ0=130.;
104			fprim.fTHETA=0.;
105			fprim.fPHI=0.;
106			fprim.fP0=1.; //1GV
107
108			f2PiAux = new PamVMC2PiAux();
109			fpgf = new PamVMCPrimaryGF();
110			}
111
112			PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
113			{
114			// Destructor
115			delete fpgf;
116			delete f2PiAux;
117			delete ftheta;
118			delete fprimColl;
119			//gObjectTable->Remove(this);
120			}
121
122
123			// private methods
124
125
126			void PamVMCPrimaryGenerator::GeneratePrimary()
127			{
128			// Add one primary particle to the user stack (derived from TVirtualMCStack).
129
130			// Track ID (filled by stack)
131			Int_t ntr;
132
133			// Option: to be tracked
134			Int_t toBeDone = 1;
135
136			// Particle type
137			Int_t pdg = fprim.fPDG;
138
139			Double_t fvx, fvy, fvz;
140			fvx=fprim.fX0;
141			fvy=fprim.fY0;
142			fvz=fprim.fZ0;
143
144			// Position
145
146			Double_t tof = 0.;
147
148			// Energy (in GeV)
149			//printf("generateprimary check fprimP0 = %f\n",fprim.fP0);
150			Double_t kinEnergy = MomentumToKinE(fprim.fP0);
151			Double_t e = fmass + kinEnergy;
152
153			// Particle momentum
154			Double_t px, py, pz;
155
156			px = TMath::Abs(fprim.fP0)Sin(fprim.fTHETA)Cos(fprim.fPHI);
157			py = TMath::Abs(fprim.fP0)Sin(fprim.fTHETA)Sin(fprim.fPHI);
158			pz = -TMath::Abs(fprim.fP0)*Cos(fprim.fTHETA);
159
160			// Polarization
161			TVector3 polar;
162
163			// Add particle to stack
164			fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof,
165			polar.X(), polar.Y(), polar.Z(),
166			kPPrimary, ntr, 1., 0);
167
168			PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
169
170			*pc = fprim;
171			}
172
173
174			void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
175			fprim.fPDG=pdg;
176			//TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
177			fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
178			fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
179			}
180
181			void PamVMCPrimaryGenerator::SetMomentum(
182			Double_t px, Double_t py, Double_t pz)
183			{
184			fprim.fP0= Sqrt(pxpx+pypy+pz*pz);
185			fprim.fTHETA=ATan(Sqrt(pxpx+pypy)/pz);
186			fprim.fPHI=ATan(py/px);
187			}
188
189			void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
190			{
191			if(isEnergy) {
192			fprim.fP0=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
193			} else{
194			fprim.fP0=frandom->Uniform(PEmin,PEmax);
195			}
196
197			}
198
199			void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
200			{
201			Double_t alpha = 1.+gamma; //integral spectral index
202			if(alpha==0.){
203			fprim.fP0=Exp(Log(PEmin)+frandom->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
204			} else {
205			if(PEmin==0.) PEmin=1.E-10;
206			fprim.fP0=Power((frandom->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
207			}
208			// cout<<"GenSpe fprim.fP0= "<<fprim.fP0<<endl;
209			if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
210
211			}
212
213
214			//Cecilia Pizzolotto: powerlaw spectrum 3 with the shape
215			// J(E) = 0.5(E + b exp(-c * sqrt(E)))^-a
216			// between PEmin and PEmax and with the input parameters a,b,c.
217			// Valeria di Felice fits parameter values are:
218			// protons: a,b,c= 2.70, 2.15, 0.21
219			// electrons: a,b,c= 0.0638, 1.248e-16, -38.248
220			void PamVMCPrimaryGenerator::GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c)
221			{
222
223			Bool_t found=0;
224			Double_t funct_min, funct_max;
225			funct_max = function3par(PEmin,a,b,c);
226			funct_min = function3par(PEmax,a,b,c);
227			//
228			Double_t wurfP;
229			Double_t wurfy ;
230			//printf("in genspe3par^^^^%f ^^%f ^^^^^^^^%f ^^^^^%f^^^^^^^\n",PEmin,PEmax,funct_min,funct_max);
231			//printf("in par^^ %f %f %f \n",a,b,c);
232			while( found==0 )
233			{
234			wurfP = frandom->Uniform(PEmin,PEmax);
235			wurfy = frandom->Uniform(funct_min,funct_max);
236			if( wurfy<(function3par(wurfP,a,b,c) ))
237			{
238			// this is ok!
239			fprim.fP0=wurfP;
240			found=1;
241			}
242			}
243			//printf("exit+++++++++++++++++++ %f %f \n",wurfP,fprim.fP0);
244			}
245
246
247
248			// cecilia pizzolotto
249			void PamVMCPrimaryGenerator::GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
250			{
251			// Generates a flat spectrum from PEmin to PElim. Then a power law
252			Double_t PElim = 1.;
253			//Double_t alpha = 1.+gamma; //integral spectral index
254
255			Bool_t okflag=0.;
256			Double_t throw_x =0.;
257			Double_t throw_y =0.;
258
259			while(okflag==0)
260			{
261			throw_x=frandom->Uniform(PEmin,PEmax);
262			// cout<<" x "<<throw_x<<endl;
263			if(throw_x<=PElim)
264			{
265			okflag=1.;
266			}
267			else
268			{
269			throw_y=frandom->Uniform(0.,1.);
270			if( throw_y<(1*pow(throw_x,gamma)))
271			{
272			okflag=1.;
273			}
274			}
275			}
276			fprim.fP0=throw_x;
277			//h->Fill(fprimf.P0);
278			okflag=0.; // reset
279
280			if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
281
282			}
283
284			// Spherical distribution -- Test by Cecilia P july 2009 ----
285			// flusso isotropo su 2pi
286			void PamVMCPrimaryGenerator::GenSphericalPhiThe()
287			{
288			// Generate phi theta
289			Double_t theta=0.;
290			Double_t phi=0.;
291
292			Double_t xcos = sqrt( frandom->Uniform(0.,1.) );
293			theta = acos(xcos); //RAD
294
295			phi = frandom->Uniform(0.,2.*Pi());
296
297			SetDirection(theta, phi);
298			return;
299			}
300
301
302
303
304
305			void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax,
306			Double_t zmin, Double_t zmax)
307			{
308			Bool_t trkGood = kFALSE;
309			Double_t theta = 999.;
310			Double_t phi = 0.;
311			Double_t x2,y2,x3,y3;
312			Double_t x0,y0,z0;
313
314			//static const Double_t rad2deg = 57.2958;
315			// S21 and S31 position/size taken as reference (z on top of det)
316			// constraint: must pass throuth these planes
317			static const Double_t s2_xmax=9.05, s2_ymax=7.55, s2_z=73.439; // z on top of det
318			static const Double_t s3_xmax=9.05, s3_ymax=7.55, s3_z=26.093; // z on top of det
319
320			//Double_t thetamax=3.14;
321			//thetamax = atan((xmax+s3_xmax)/(zmax-s3_z));
322			//cout<<" Quanto è il theta max? "<<thetamax<<" in deg "<<thetamax*(90./Pi())<<endl;
323
324			while (trkGood!=kTRUE)
325			{
326			x0= frandom->Uniform(xmin,xmax);
327			y0= frandom->Uniform(ymin,ymax);
328			z0= frandom->Uniform(zmin,zmax);
329
330			// Generate phi theta
331			theta=999.; // init
332			while (theta>=0.65) // take only theta smaller than 37deg=0.65rad
333			{
334			Double_t xcos = sqrt( frandom->Uniform(0.,1.) );
335			theta = acos(xcos); //RAD
336			}
337			phi = frandom->Uniform(0.,2.*Pi());
338
339			// Calculate xy at the constraint
340			Double_t fact2 = (s2_z-z0)/cos(theta);
341			x2 = x0 + fabs(fact2) * sin(theta) * cos(phi);
342			y2 = y0 + fabs(fact2) * sin(theta) * sin(phi);
343			Double_t fact3 = (s3_z-z0)/cos(theta);
344			x3 = x0 + fabs(fact3) * sin(theta) * cos(phi);
345			y3 = y0 + fabs(fact3) * sin(theta) * sin(phi);
346
347			//cout<<" x/y0= "<<x0<<" "<<y0<<" x/y2= "<<fact2sin(theta)cos(phi)<<" "<<x2<<" xy3= "<<
348			// fact3sin(theta)cos(phi)<<" "<<x3<<" phi/the "<<phi(90./Pi())<<" "<<theta(90./Pi())<<endl;
349
350			// Test condition on the direction
351			if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) &&
352			Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) {
353			trkGood = kTRUE;
354			//cout<<" x/y0= "<<x0<<" "<<y0<<" x/y2= "<<fact2sin(theta)cos(phi)<<" "<<x2<<" xy3= "<<
355			// fact3sin(theta)cos(phi)<<" "<<x3<<endl;
356			}
357			}
358
359			// Set direction and position:
360			SetDirection(theta, phi);
361			SetPosition(x0, y0, z0);
362
363			return;
364			}