--- PamVMC/src/PamVMCPrimaryGenerator.cxx	2009/02/18 17:38:17	1.2
+++ PamVMC/src/PamVMCPrimaryGenerator.cxx	2010/09/15 07:01:57	1.6
@@ -8,16 +8,11 @@
 #include <TParticlePDG.h>
 #include <TVector3.h>
 #include <TMath.h>
+#include <Riostream.h>
 
 #include "PamVMCPrimaryGenerator.h"
 
-using TMath::Sqrt;
-using TMath::Sin;
-using TMath::Cos;
-using TMath::ATan;
-using TMath::Log;
-using TMath::Power;
-using TMath::Exp;
+using namespace TMath;
 
 ClassImp(PamVMCPrimary)
 
@@ -44,12 +39,14 @@
     fevno(0),
     fmass(0.),
     fcharge(0.),
-    frnd(0)
+    frandom(0)
 {
 // Standard constructor
-  fprimColl = new TClonesArray("PamVMCPrimary");
-  frnd = new TRandom3(0);
 
+  ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
+  ftheta->SetNpx(1000);
+
+  fprimColl = new TClonesArray("PamVMCPrimary");
   fprim.fPDG=kProton;
   fprim.fX0=1.;
   fprim.fY0=1.;
@@ -67,11 +64,13 @@
     fmass(0.),
     fcharge(0.),
     fprimColl(0),
-    frnd(0)
+    frandom(0)
 {    
-  frnd = new TRandom3(0);
-// Default constructor
+ // Default constructor
   //Default primary proton
+  ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
+  ftheta->SetNpx(1000);
+
   fprim.fPDG=kProton;
   fprim.fX0=1.;
   fprim.fY0=1.;
@@ -84,13 +83,12 @@
 PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
 {
 // Destructor
-  delete frnd;
+  delete ftheta;
   delete fprimColl;
 }
 
 // private methods
 
-#include <Riostream.h>
 
 void PamVMCPrimaryGenerator::GeneratePrimary()
 {    
@@ -115,6 +113,7 @@
   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;
  
@@ -157,9 +156,9 @@
 void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
 {
   if(isEnergy) {
-    fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
+    fprim.fP0=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
   } else{
-    fprim.fP0=frnd->Uniform(PEmin,PEmax);
+    fprim.fP0=frandom->Uniform(PEmin,PEmax);
   }
 
 }
@@ -168,12 +167,165 @@
 {
   Double_t alpha = 1.+gamma; //integral spectral index
   if(alpha==0.){
-    fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
+    fprim.fP0=Exp(Log(PEmin)+frandom->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);
+    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;
+}