--- trieste/pamVMC/src/PamVMCPrimaryGenerator.cxx 2009/03/04 12:51:18 1.1.1.1 +++ trieste/pamVMC/src/PamVMCPrimaryGenerator.cxx 2009/03/24 14:04:06 1.3 @@ -43,7 +43,7 @@ ClassImp(PamVMCPrimaryGenerator) -PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) + PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed) : TObject(), fStack(stack), fevno(0), @@ -53,7 +53,7 @@ { // Standard constructor fprimColl = new TClonesArray("PamVMCPrimary"); - frnd = new TRandom3(0); + frnd = new TRandom3(seed); fprim.fPDG=kProton; fprim.fX0=1.; @@ -65,6 +65,27 @@ } +// 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), @@ -114,7 +135,6 @@ fvx=fprim.fX0; fvy=fprim.fY0; fvz=fprim.fZ0; - // Position Double_t tof = 0.; @@ -126,12 +146,12 @@ // Particle momentum Double_t px, py, pz; - px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); - 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); - //py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI); - //pz = -fprim.fP0*Cos(fprim.fTHETA); + //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; @@ -172,6 +192,9 @@ } + + + void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy) { Double_t alpha = 1.+gamma; //integral spectral index @@ -186,19 +209,99 @@ } + +// 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= "< 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 = (Pi()/180.)*frnd->Uniform(Phmin,Phmax); + phi = frnd->Uniform(Phmin,Phmax); do { y = frnd->Uniform(0.,1.); - theta = (Pi()/180.)*frnd->Uniform(Thmin,Thmax); + theta = frnd->Uniform(Thmin,Thmax); f = Sin(theta); } while (y>f); - //theta = phi = (Pi()/180.)*45.; - SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); + //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. @@ -206,7 +309,7 @@ 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)))/(Pi()/180.), (frnd->Uniform(0.,2*Pi()))/(Pi()/180.)); + SetDirection( frnd->Uniform((theta-ang),(theta+ang)) , frnd->Uniform(0.,2*Pi()) ); //SetDirection(0., 0.); } @@ -224,7 +327,8 @@ 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)); + //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 @@ -234,20 +338,22 @@ do { //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; - GenSphDist(r, 0., thmax, 0., phmax); - - Double_t th = (Pi()/180.)*fprim.fTHETA; - Double_t ph = (Pi()/180.)*fprim.fPHI; + + 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); - y1 = s1_pz*Tan(th)*Sin(ph); - x2 = s2_pz*Tan(th)*Cos(ph); - y2 = s2_pz*Tan(th)*Sin(ph); - x3 = s3_pz*Tan(th)*Cos(ph); - y3 = s3_pz*Tan(th)*Sin(ph); + 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) && @@ -255,4 +361,5 @@ //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; }while (!trkGood); //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; + }