pamVMC/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 "PamVMCPrimaryGenerator.h"

using TMath::Sqrt;
using TMath::Sin;
using TMath::Cos;
using TMath::ACos;
using TMath::Tan;
using TMath::ATan;
using TMath::ATan2;
using TMath::Log;
using TMath::Power;
using TMath::Exp;
using TMath::Pi;
using TMath::Abs;

ClassImp(PamVMCPrimary)

PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
{
  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;

  return a;
}


ClassImp(PamVMCPrimaryGenerator)

  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed) 
  : TObject(),
    fStack(stack),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    frnd(0)
{
// Standard constructor
  fprimColl = new TClonesArray("PamVMCPrimary");
  frnd = new TRandom3(seed);

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

} 

// PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(UInt_t seed)
//   : TObject(),
//     fStack(0),
//     fevno(0),
//     fmass(0.),
//     fcharge(0.),
//     fprimColl(0),
//     frnd(0)
// {    
//   frnd = new TRandom3(seed);
// // Default constructor
//   //Default primary proton
//   fprim.fPDG=kProton;
//   fprim.fX0=1.; 
//   fprim.fY0=1.;
//   fprim.fZ0=130.;
//   fprim.fTHETA=0.;
//   fprim.fPHI=0.;
//   fprim.fP0=1.; //1GV
// }

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
  : TObject(),
    fStack(0),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    fprimColl(0),
    frnd(0)
{    
  frnd = new TRandom3(0);
// Default constructor
  //Default primary proton
  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV
}

PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
{
// Destructor
  delete frnd;
  delete fprimColl;
}

// private methods

#include <Riostream.h>

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)
  Double_t kinEnergy = MomentumToKinE(fprim.fP0);     
  Double_t e  = fmass + kinEnergy;
 
  // Particle momentum
  Double_t  px, py, pz;
   
  //px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); // ritabrata
  //py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI));
  //pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian
  px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI); // RADIANTS
  py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
  pz = -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=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
  } else{
    fprim.fP0=frnd->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)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
  } else {
    if(PEmin==0.) PEmin=1.E-10;
    fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
  }

  if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);

}


// Spherical distribution -- Test by Cecilia P march 2009 -----------------------------///
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 ------------------------------------///


void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)
{
  // all angles in RADIANTS
        Double_t theta, phi, y, f;
        phi = frnd->Uniform(Phmin,Phmax);
        
        do
        {       y = frnd->Uniform(0.,1.);
                theta = frnd->Uniform(Thmin,Thmax);
                f = Sin(theta);
        } while (y>f);
        
        //SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); // ritabrata
        SetPosition(1.,1.,130.); // cecilia


        //random distribution of theta phi in the angle at the vertex at (0,0,r)
        //by the S3 max distant corners.
        static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;
        Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz));
        
        //SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform((phi-ang),(phi+ang)))/(Pi()/180.));
        SetDirection( frnd->Uniform((theta-ang),(theta+ang)) ,  frnd->Uniform(0.,2*Pi()) );
        //SetDirection(0., 0.);
}

void PamVMCPrimaryGenerator::GenSphDist(Double_t r)
{
        //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
        static const Double_t s1_x=20.4, s1_y=16.5, s1_pz=102.8866;
        static const Double_t s2_x=9.0, s2_y=7.5, s2_pz=73.489;
        static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;

        //calculate max theta and phi angles to be allowed (calculating wrt one 
        //corner of s3 to opposite corner of s1)
        Double_t rx = s1_x+s3_x;
        Double_t ry = s1_y+s3_y;
        Double_t rz = s1_pz-s3_pz;

        Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz)));
        //Double_t phmax = (180./Pi())*(ATan2(ry,rx));

        //cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl;

        //generate a track check and let it go only if it is passing through all
        //the 3 TOF
        Bool_t trkGood = kFALSE;

        do
        {
        //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;

                GenSphDist(r, 0., thmax, 0., 2*Pi() );
        
                Double_t th = fprim.fTHETA;
                Double_t ph = fprim.fPHI;

        //cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl;

                Double_t x1, y1, x2, y2, x3, y3;
                
                x1 = s1_pz*Tan(th)*Cos(ph) - fprim.fX0;
                y1 = s1_pz*Tan(th)*Sin(ph) - fprim.fY0;
                x2 = s2_pz*Tan(th)*Cos(ph) - fprim.fX0;
                y2 = s2_pz*Tan(th)*Sin(ph) - fprim.fY0;
                x3 = s3_pz*Tan(th)*Cos(ph) - fprim.fX0;
                y3 = s3_pz*Tan(th)*Sin(ph) - fprim.fY0;

                if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) &&
                     Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) &&
                     Abs(x3) <= Abs(s3_x) && Abs(y3) <= Abs(s3_y) ) trkGood = kTRUE;
        //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
        }while (!trkGood);
        //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;

}
1	pamelats	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
12			#include "PamVMCPrimaryGenerator.h"
13
14			using TMath::Sqrt;
15			using TMath::Sin;
16			using TMath::Cos;
17			using TMath::ACos;
18			using TMath::Tan;
19			using TMath::ATan;
20			using TMath::ATan2;
21			using TMath::Log;
22			using TMath::Power;
23			using TMath::Exp;
24			using TMath::Pi;
25			using TMath::Abs;
26
27			ClassImp(PamVMCPrimary)
28
29			PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
30			{
31			a.fPDG=b.fPDG;
32			a.fX0=b.fX0;
33			a.fY0=b.fY0;
34			a.fZ0=b.fZ0;
35			a.fTHETA=b.fTHETA;
36			a.fPHI=b.fPHI;
37			a.fP0=b.fP0;
38			a.fGOOD=b.fGOOD;
39
40			return a;
41			}
42
43
44			ClassImp(PamVMCPrimaryGenerator)
45
46	pizzolot	1.3	PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed)
47	pamelats	1.1	: TObject(),
48			fStack(stack),
49			fevno(0),
50			fmass(0.),
51			fcharge(0.),
52			frnd(0)
53			{
54			// Standard constructor
55			fprimColl = new TClonesArray("PamVMCPrimary");
56	pizzolot	1.3	frnd = new TRandom3(seed);
57	pamelats	1.1
58			fprim.fPDG=kProton;
59			fprim.fX0=1.;
60			fprim.fY0=1.;
61			fprim.fZ0=130.;
62			fprim.fTHETA=0.;
63			fprim.fPHI=0.;
64			fprim.fP0=1.; //1GV
65
66			}
67
68	pizzolot	1.3	// PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(UInt_t seed)
69			// : TObject(),
70			// fStack(0),
71			// fevno(0),
72			// fmass(0.),
73			// fcharge(0.),
74			// fprimColl(0),
75			// frnd(0)
76			// {
77			// frnd = new TRandom3(seed);
78			// // Default constructor
79			// //Default primary proton
80			// fprim.fPDG=kProton;
81			// fprim.fX0=1.;
82			// fprim.fY0=1.;
83			// fprim.fZ0=130.;
84			// fprim.fTHETA=0.;
85			// fprim.fPHI=0.;
86			// fprim.fP0=1.; //1GV
87			// }
88
89	pamelats	1.1	PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
90			: TObject(),
91			fStack(0),
92			fevno(0),
93			fmass(0.),
94			fcharge(0.),
95			fprimColl(0),
96			frnd(0)
97			{
98			frnd = new TRandom3(0);
99			// Default constructor
100			//Default primary proton
101			fprim.fPDG=kProton;
102			fprim.fX0=1.;
103			fprim.fY0=1.;
104			fprim.fZ0=130.;
105			fprim.fTHETA=0.;
106			fprim.fPHI=0.;
107			fprim.fP0=1.; //1GV
108			}
109
110			PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
111			{
112			// Destructor
113			delete frnd;
114			delete fprimColl;
115			}
116
117			// private methods
118
119			#include <Riostream.h>
120
121			void PamVMCPrimaryGenerator::GeneratePrimary()
122			{
123			// Add one primary particle to the user stack (derived from TVirtualMCStack).
124
125			// Track ID (filled by stack)
126			Int_t ntr;
127
128			// Option: to be tracked
129			Int_t toBeDone = 1;
130
131			// Particle type
132			Int_t pdg = fprim.fPDG;
133
134			Double_t fvx, fvy, fvz;
135			fvx=fprim.fX0;
136			fvy=fprim.fY0;
137			fvz=fprim.fZ0;
138			// Position
139
140			Double_t tof = 0.;
141
142			// Energy (in GeV)
143			Double_t kinEnergy = MomentumToKinE(fprim.fP0);
144			Double_t e = fmass + kinEnergy;
145
146			// Particle momentum
147			Double_t px, py, pz;
148
149	pizzolot	1.3	//px = -fprim.fP0Sin((Pi()/180.)(fprim.fTHETA))Cos((Pi()/180.)(fprim.fPHI)); // ritabrata
150			//py = -fprim.fP0Sin((Pi()/180.)(fprim.fTHETA))Sin((Pi()/180.)(fprim.fPHI));
151			//pz = -fprim.fP0Cos((Pi()/180.)(fprim.fTHETA)); // converting in radian
152			px = fprim.fP0Sin(fprim.fTHETA)Cos(fprim.fPHI); // RADIANTS
153			py = fprim.fP0Sin(fprim.fTHETA)Sin(fprim.fPHI);
154			pz = -fprim.fP0*Cos(fprim.fTHETA);
155	pamelats	1.1
156			// Polarization
157			TVector3 polar;
158
159			// Add particle to stack
160			fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof,
161			polar.X(), polar.Y(), polar.Z(),
162			kPPrimary, ntr, 1., 0);
163
164			PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
165
166			*pc = fprim;
167			}
168
169
170			void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
171			fprim.fPDG=pdg;
172			//TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
173			fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
174			fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
175			}
176
177			void PamVMCPrimaryGenerator::SetMomentum(
178			Double_t px, Double_t py, Double_t pz)
179			{
180			fprim.fP0= Sqrt(pxpx+pypy+pz*pz);
181			fprim.fTHETA=ATan(Sqrt(pxpx+pypy)/pz);
182			fprim.fPHI=ATan(py/px);
183			}
184
185			void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
186			{
187			if(isEnergy) {
188			fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
189			} else{
190			fprim.fP0=frnd->Uniform(PEmin,PEmax);
191			}
192
193			}
194
195	pizzolot	1.3
196
197
198	pamelats	1.1	void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
199			{
200			Double_t alpha = 1.+gamma; //integral spectral index
201			if(alpha==0.){
202			fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
203			} else {
204			if(PEmin==0.) PEmin=1.E-10;
205			fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
206			}
207
208			if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
209
210			}
211
212	pizzolot	1.3
213			// Spherical distribution -- Test by Cecilia P march 2009 -----------------------------///
214			void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax,
215			Double_t zmin, Double_t zmax)
216			{
217			Bool_t trkGood = kFALSE;
218			Double_t theta = 999.;
219			Double_t phi = 0.;
220			Double_t x2,y2,x3,y3;
221
222			//static const Double_t rad2deg = 57.2958;
223
224			//static const Double_t s1_xmax=20.4, s1_ymax=16.5, s1_pz=102.8866;
225			// calo cavity 8.07x6.57 with z2=71.6 z3=27.4 circa
226			static const Double_t s2_xmax=7.8, s2_ymax=6.0, s2_z=73.489;
227			static const Double_t s3_xmax=8.0, s3_ymax=6.0, s3_z=25.3159;
228
229			// Generate random position unif. distr.
230			// GenPosition(-25.,25., 25.,25., 108.0,108.0);
231			// Double_t x0= fprim.fX0;
232			// Double_t y0= fprim.fY0;
233			// Double_t z0= fprim.fZ0;
234			Double_t x0= frnd->Uniform(xmin,xmax);
235			Double_t y0= frnd->Uniform(ymin,ymax);
236			Double_t z0= frnd->Uniform(zmin,zmax);
237			Int_t posGood = 0;
238
239
240			while (trkGood!=kTRUE)
241			{
242			// Generate phi theta
243			theta=999.; // init
244			while (theta>=0.82) //take only theta smaller than 30deg=0.52rad
245			{
246			Double_t xcos = sqrt( frnd->Uniform(0.,1.) );
247			theta = acos(xcos); //RAD
248			}
249			phi = frnd->Uniform(0.,2.*Pi());
250
251			// Calculate xy at beginning/end of magnet cavity
252			Double_t fact2 = (s2_z-z0)/cos(theta);
253			x2 = x0 + fabs(fact2) * sin(theta) * cos(phi);
254			y2 = y0 + fabs(fact2) * sin(theta) * sin(phi);
255			Double_t fact3 = (s3_z-z0)/cos(theta);
256			x3 = x0 + fabs(fact3) * sin(theta) * cos(phi);
257			y3 = y0 + fabs(fact3) * sin(theta) * sin(phi);
258
259			// Test condition on the direction
260			if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) &&
261			Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) {
262			trkGood = kTRUE;
263			//cout<<" x/y0= "<<x0<<" "<<y0<<" x/y2= "<<fact2sin(theta)cos(phi)<<" "<<x2<<" xy3= "<<
264			// fact3sin(theta)cos(phi)<<" "<<x3<<endl;
265			}
266
267			//if from current x0y0z0 the condition are not satysfied for any angle => def new start position
268			posGood++;
269			if ( posGood == 100 )
270			{
271			x0= frnd->Uniform(xmin,xmax);
272			y0= frnd->Uniform(ymin,ymax);
273			z0= frnd->Uniform(zmin,zmax);
274			posGood = 0;
275			}
276
277			}
278
279			// Set direction and position:
280			SetDirection(theta, phi);
281			SetPosition(x0, y0, z0);
282
283			return;
284			}
285			//--- end Test by Cecilia ------------------------------------///
286
287
288
289
290	pamelats	1.1	void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)
291			{
292	pizzolot	1.3	// all angles in RADIANTS
293	pamelats	1.1	Double_t theta, phi, y, f;
294	pizzolot	1.3	phi = frnd->Uniform(Phmin,Phmax);
295	pamelats	1.1
296			do
297			{ y = frnd->Uniform(0.,1.);
298	pizzolot	1.3	theta = frnd->Uniform(Thmin,Thmax);
299	pamelats	1.1	f = Sin(theta);
300			} while (y>f);
301
302	pizzolot	1.3	//SetPosition((rSin(theta)Cos(phi)), (rSin(theta)Sin(phi)), (r*Cos(theta))); // ritabrata
303			SetPosition(1.,1.,130.); // cecilia
304
305	pamelats	1.1
306			//random distribution of theta phi in the angle at the vertex at (0,0,r)
307			//by the S3 max distant corners.
308			static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;
309			Double_t ang = ATan((Sqrt(s3_xs3_x+s3_ys3_y))/(r-s3_pz));
310
311			//SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform((phi-ang),(phi+ang)))/(Pi()/180.));
312	pizzolot	1.3	SetDirection( frnd->Uniform((theta-ang),(theta+ang)) , frnd->Uniform(0.,2*Pi()) );
313	pamelats	1.1	//SetDirection(0., 0.);
314			}
315
316			void PamVMCPrimaryGenerator::GenSphDist(Double_t r)
317			{
318			//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
319			static const Double_t s1_x=20.4, s1_y=16.5, s1_pz=102.8866;
320			static const Double_t s2_x=9.0, s2_y=7.5, s2_pz=73.489;
321			static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159;
322
323			//calculate max theta and phi angles to be allowed (calculating wrt one
324			//corner of s3 to opposite corner of s1)
325			Double_t rx = s1_x+s3_x;
326			Double_t ry = s1_y+s3_y;
327			Double_t rz = s1_pz-s3_pz;
328
329			Double_t thmax = (180./Pi())(ACos(rz/Sqrt(rxrx+ryry+rzrz)));
330	pamelats	1.2	//Double_t phmax = (180./Pi())*(ATan2(ry,rx));
331	pizzolot	1.3
332	pamelats	1.1	//cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl;
333
334			//generate a track check and let it go only if it is passing through all
335			//the 3 TOF
336			Bool_t trkGood = kFALSE;
337
338			do
339			{
340			//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
341	pizzolot	1.3
342			GenSphDist(r, 0., thmax, 0., 2*Pi() );
343
344			Double_t th = fprim.fTHETA;
345			Double_t ph = fprim.fPHI;
346
347	pamelats	1.1	//cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl;
348
349			Double_t x1, y1, x2, y2, x3, y3;
350
351	pamelats	1.2	x1 = s1_pzTan(th)Cos(ph) - fprim.fX0;
352			y1 = s1_pzTan(th)Sin(ph) - fprim.fY0;
353			x2 = s2_pzTan(th)Cos(ph) - fprim.fX0;
354			y2 = s2_pzTan(th)Sin(ph) - fprim.fY0;
355			x3 = s3_pzTan(th)Cos(ph) - fprim.fX0;
356			y3 = s3_pzTan(th)Sin(ph) - fprim.fY0;
357	pamelats	1.1
358			if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) &&
359			Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) &&
360			Abs(x3) <= Abs(s3_x) && Abs(y3) <= Abs(s3_y) ) trkGood = kTRUE;
361			//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
362			}while (!trkGood);
363			//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
364	pizzolot	1.3
365	pamelats	1.1	}