/[PAMELA software]/trieste/pamVMC/src/PamVMCPrimaryGenerator.cxx
ViewVC logotype

Diff of /trieste/pamVMC/src/PamVMCPrimaryGenerator.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.1 by pamelats, Wed Mar 4 12:51:18 2009 UTC revision 1.3 by pizzolot, Tue Mar 24 14:04:06 2009 UTC
# Line 43  PamVMCPrimary & operator+=(PamVMCPrimary Line 43  PamVMCPrimary & operator+=(PamVMCPrimary
43    
44  ClassImp(PamVMCPrimaryGenerator)  ClassImp(PamVMCPrimaryGenerator)
45    
46  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)    PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed)
47    : TObject(),    : TObject(),
48      fStack(stack),      fStack(stack),
49      fevno(0),      fevno(0),
# Line 53  PamVMCPrimaryGenerator::PamVMCPrimaryGen Line 53  PamVMCPrimaryGenerator::PamVMCPrimaryGen
53  {  {
54  // Standard constructor  // Standard constructor
55    fprimColl = new TClonesArray("PamVMCPrimary");    fprimColl = new TClonesArray("PamVMCPrimary");
56    frnd = new TRandom3(0);    frnd = new TRandom3(seed);
57    
58    fprim.fPDG=kProton;    fprim.fPDG=kProton;
59    fprim.fX0=1.;    fprim.fX0=1.;
# Line 65  PamVMCPrimaryGenerator::PamVMCPrimaryGen Line 65  PamVMCPrimaryGenerator::PamVMCPrimaryGen
65    
66  }  }
67    
68    // 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  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
90    : TObject(),    : TObject(),
91      fStack(0),      fStack(0),
# Line 114  void PamVMCPrimaryGenerator::GeneratePri Line 135  void PamVMCPrimaryGenerator::GeneratePri
135    fvx=fprim.fX0;    fvx=fprim.fX0;
136    fvy=fprim.fY0;    fvy=fprim.fY0;
137    fvz=fprim.fZ0;    fvz=fprim.fZ0;
   
138    // Position    // Position
139        
140    Double_t tof = 0.;    Double_t tof = 0.;
# Line 126  void PamVMCPrimaryGenerator::GeneratePri Line 146  void PamVMCPrimaryGenerator::GeneratePri
146    // Particle momentum    // Particle momentum
147    Double_t  px, py, pz;    Double_t  px, py, pz;
148        
149    px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI));    //px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); // ritabrata
150    py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI));    //py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI));
151    pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian    //pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian
152    //px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI);    px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI); // RADIANTS
153    //py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);    py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
154    //pz = -fprim.fP0*Cos(fprim.fTHETA);    pz = -fprim.fP0*Cos(fprim.fTHETA);
155    
156    // Polarization    // Polarization
157    TVector3 polar;    TVector3 polar;
# Line 172  void PamVMCPrimaryGenerator::GenSpe(Doub Line 192  void PamVMCPrimaryGenerator::GenSpe(Doub
192    
193  }  }
194    
195    
196    
197    
198  void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)  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    Double_t alpha = 1.+gamma; //integral spectral index
# Line 186  void PamVMCPrimaryGenerator::GenSpe(Doub Line 209  void PamVMCPrimaryGenerator::GenSpe(Doub
209    
210  }  }
211    
212    
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= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<"  xy3= "<<
264            //  fact3*sin(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  void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)  void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax)
291  {  {
292      // all angles in RADIANTS
293          Double_t theta, phi, y, f;          Double_t theta, phi, y, f;
294          phi = (Pi()/180.)*frnd->Uniform(Phmin,Phmax);          phi = frnd->Uniform(Phmin,Phmax);
295                    
296          do          do
297          {       y = frnd->Uniform(0.,1.);          {       y = frnd->Uniform(0.,1.);
298                  theta = (Pi()/180.)*frnd->Uniform(Thmin,Thmax);                  theta = frnd->Uniform(Thmin,Thmax);
299                  f = Sin(theta);                  f = Sin(theta);
300          } while (y>f);          } while (y>f);
301                    
302          //theta = phi = (Pi()/180.)*45.;          //SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); // ritabrata
303          SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta)));          SetPosition(1.,1.,130.); // cecilia
304    
305    
306          //random distribution of theta phi in the angle at the vertex at (0,0,r)          //random distribution of theta phi in the angle at the vertex at (0,0,r)
307          //by the S3 max distant corners.          //by the S3 max distant corners.
# Line 206  void PamVMCPrimaryGenerator::GenSphDist( Line 309  void PamVMCPrimaryGenerator::GenSphDist(
309          Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz));          Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz));
310                    
311          //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((phi-ang),(phi+ang)))/(Pi()/180.));
312          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()) );
313          //SetDirection(0., 0.);          //SetDirection(0., 0.);
314  }  }
315    
# Line 224  void PamVMCPrimaryGenerator::GenSphDist( Line 327  void PamVMCPrimaryGenerator::GenSphDist(
327          Double_t rz = s1_pz-s3_pz;          Double_t rz = s1_pz-s3_pz;
328    
329          Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz)));          Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz)));
330          Double_t phmax = (180./Pi())*(ATan2(ry,rx));          //Double_t phmax = (180./Pi())*(ATan2(ry,rx));
331    
332          //cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl;          //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          //generate a track check and let it go only if it is passing through all
# Line 234  void PamVMCPrimaryGenerator::GenSphDist( Line 338  void PamVMCPrimaryGenerator::GenSphDist(
338          do          do
339          {          {
340          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
341                  GenSphDist(r, 0., thmax, 0., phmax);  
342                                    GenSphDist(r, 0., thmax, 0., 2*Pi() );
343                  Double_t th = (Pi()/180.)*fprim.fTHETA;          
344                  Double_t ph = (Pi()/180.)*fprim.fPHI;                  Double_t th = fprim.fTHETA;
345                    Double_t ph = fprim.fPHI;
346    
347          //cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl;          //cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl;
348    
349                  Double_t x1, y1, x2, y2, x3, y3;                  Double_t x1, y1, x2, y2, x3, y3;
350                                    
351                  x1 = s1_pz*Tan(th)*Cos(ph);                  x1 = s1_pz*Tan(th)*Cos(ph) - fprim.fX0;
352                  y1 = s1_pz*Tan(th)*Sin(ph);                  y1 = s1_pz*Tan(th)*Sin(ph) - fprim.fY0;
353                  x2 = s2_pz*Tan(th)*Cos(ph);                  x2 = s2_pz*Tan(th)*Cos(ph) - fprim.fX0;
354                  y2 = s2_pz*Tan(th)*Sin(ph);                  y2 = s2_pz*Tan(th)*Sin(ph) - fprim.fY0;
355                  x3 = s3_pz*Tan(th)*Cos(ph);                  x3 = s3_pz*Tan(th)*Cos(ph) - fprim.fX0;
356                  y3 = s3_pz*Tan(th)*Sin(ph);                  y3 = s3_pz*Tan(th)*Sin(ph) - fprim.fY0;
357    
358                  if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) &&                  if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) &&
359                       Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) &&                       Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) &&
# Line 255  void PamVMCPrimaryGenerator::GenSphDist( Line 361  void PamVMCPrimaryGenerator::GenSphDist(
361          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
362          }while (!trkGood);          }while (!trkGood);
363          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;          //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl;
364    
365  }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23