/[PAMELA software]/PamVMC/src/PamVMCPrimaryGenerator.cxx
ViewVC logotype

Diff of /PamVMC/src/PamVMCPrimaryGenerator.cxx

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

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