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

Diff of /PamVMC/src/PamVMCPrimaryGenerator.cxx

Parent Directory | Revision Log | View Patch Patch

-revision 1.5 by pam-rm2,
Fri Jun 12 18:39:32 2009 UTC
+revision 1.6 by pizzolot,
Wed Sep 15 07:01:57 2010 UTC
 Line 113 
 void PamVMCPrimaryGenerator::GeneratePri
    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;
-Line 171 
 void PamVMCPrimaryGenerator::GenSpe(Doub
+Line 172 
 void PamVMCPrimaryGenerator::GenSpe(Doub
      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;
+ }

 Legend:



Removed from v.1.5
 


changed lines


 
Added in v.1.6
 Legend:



Removed from v.1.5
 


changed lines


 
Added in v.1.6
-Removed from v.1.5
+Added in v.1.6

	ViewVC Help
Powered by ViewVC 1.1.23