STARLIGHT update:
authorcmayer <Christoph.Mayer@cern.ch>
Thu, 8 Jan 2015 09:08:24 +0000 (10:08 +0100)
committercmayer <Christoph.Mayer@cern.ch>
Thu, 8 Jan 2015 09:08:24 +0000 (10:08 +0100)
 - update to r193
 - propagation of kinematical parameters for gamma-gamma processes
 - added a cut on the rapidity of the mother partice

16 files changed:
STARLIGHT/starlight/include/gammaavm.h
STARLIGHT/starlight/include/gammagammaleptonpair.h
STARLIGHT/starlight/include/gammagammasingle.h
STARLIGHT/starlight/include/photonNucleusCrossSection.h
STARLIGHT/starlight/include/readinluminosity.h
STARLIGHT/starlight/src/beambeamsystem.cpp
STARLIGHT/starlight/src/gammaaluminosity.cpp
STARLIGHT/starlight/src/gammaavm.cpp
STARLIGHT/starlight/src/gammagammaleptonpair.cpp
STARLIGHT/starlight/src/gammagammasingle.cpp
STARLIGHT/starlight/src/incoherentPhotonNucleusLuminosity.cpp
STARLIGHT/starlight/src/incoherentVMCrossSection.cpp
STARLIGHT/starlight/src/narrowResonanceCrossSection.cpp
STARLIGHT/starlight/src/photonNucleusCrossSection.cpp
STARLIGHT/starlight/src/readinluminosity.cpp
STARLIGHT/starlight/src/wideResonanceCrossSection.cpp

index 42f9425..47886e0 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+// $Rev:: 192                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:25 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -77,6 +77,9 @@ class Gammaavectormeson : public eventChannel
   int _VMinterferencemode;
   int _VMCoherence;
   int _ProductionMode;
+  int _TargetBeam; 
+  int N1;
+  int N2; 
   double _VMCoherenceFactor;
   double _VMgamma_em;
   double _VMNPT;
index de7b6f7..5a82db7 100644 (file)
@@ -73,7 +73,7 @@ class Gammagammaleptonpair : public eventChannel
   void pairMomentum(double w,double y,double &E,double &px,double &py,double&pz);
   double pp_1(double E);
   double pp_2(double E);
-  void twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double E,double W,double px0,double py0,double pz0,double &px1,double &py1,double&pz1,double &px2,double &py2,/*double &py2,*/double &pz2,int &iFbadevent);
+  void twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double E,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &pt1,double &px2,double &py2,double &pz2,double &pt2,double &mdec,int &iFbadevent);
   double thetalep(double W,double theta);
   void tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2);
   
index 344663d..0a8354e 100644 (file)
@@ -70,7 +70,12 @@ class Gammagammasingle : public eventChannel
   
   void parentMomentum(double w,double y,double &E,double &px,double &py,double&pz);
   double pp(double E);
-  void twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double E,double W,double px0,double py0,double pz0,double &px1,double &py1,double&pz1,double &px2,double &py2,/*double &py2,*/double &pz2,int &iFbadevent);
+  void twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double E,double W,
+                   double px0,double py0,double pz0,
+                   double &px1,double &py1,double &pz1, double &pt1,
+                   double &px2,double &py2,double &pz2, double &pt2,
+                   double &mdec,
+                   int &iFbadevent);
   // void transform(double betax,double betay,double betaz,double &E,double &px,double &py,double &pz,int &iFbadevent);
   void thephi(double W,double px,double py,double pz,double E,double &theta,double &phi);
   
index ca33b83..9fb0df2 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+// $Rev:: 192                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:25 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -60,10 +60,13 @@ public:
        // So just use the wide or narrow constructor to calculate it
        // wide/narrow will inherit this.
        double getcsgA(const double Egamma,
-                      const double W);
-       double photonFlux(const double Egamma);
+                      const double W,
+                       const int beam);
+       double photonFlux(const double Egamma,
+                       const int beam);
        double sigmagp(const double Wgp);
-       double sigma_A(const double sig_N);
+       double sigma_A(const double sig_N, 
+                       const int beam);
         double sigma_N(const double Wgp);
        double breitWigner(const double W,
                           const double C);
index 8755817..69381c1 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+// $Rev:: 192                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:25 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -52,11 +52,16 @@ class readLuminosity
   */
   double *_Warray;
   double *_Yarray;
-  double **_Farray;
+  double **_Farray; 
+  double **_Farray1;
+  double **_Farray2;
   
   double _f_max;
+  double _f_max1;
+  double _f_max2;
+
   double _fptarray[500][500];
-  //           inputParameters inputread;
+
   double _bwnormsave;
 
  protected:
index 0b48e71..ac483d1 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 176                         $: revision of last commit
-// $Author:: jseger                   $: author of last commit
-// $Date:: 2014-06-20 22:15:20 +0200 #$: date of last commit
+// $Rev:: 193                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -188,7 +188,6 @@ beamBeamSystem::generateBreakupProbabilities()
         if (_beamBreakupMode == 7)
             printInfo << "Requiring breakup of one nucleus (Xn,0n). " << endl;
 
-        //pp may cause segmentation fault in here and it does not use this...
        double pOfB = 0;
        double b = bMin;
         double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
@@ -212,7 +211,7 @@ beamBeamSystem::generateBreakupProbabilities()
             {
                 if((1-pOfB)<_breakupCutOff)
                 {
- //                         std::cout << "Break off b: " << b << std::endl;
+//                         std::cout << "Break off b: " << b << std::endl;
 //                         std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
                         break;
                 }
@@ -239,7 +238,7 @@ beamBeamSystem::generateBreakupProbabilities()
             b *= _breakupImpactParameterStep;
            pOfB = exp(-1 * _pHadronBreakup) * _pPhotonBreakup;
             _breakupProbabilities.push_back(pOfB);
-        }
+        } // End while(1)
     }
     else if (((_beam1.Z() == 1) && (_beam1.A() == 1)) || ((_beam2.Z() == 1) && (_beam2.A() == 1))) {  
       
@@ -282,21 +281,20 @@ beamBeamSystem::generateBreakupProbabilities()
 double
 beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
 {
-       //      double pbreakup =0.; 
        //probability of hadron breakup, 
        //this is what is returned when the function is called
        double gamma = _beamLorentzGamma; 
        //input for gamma_em
        //this will need to be StarlightInputParameters::gamma or whatever
        double b = impactparameter;
-       int ap = _beam1.A();  
-       //Notice this is  taking from nucleus 1.Still assuming symmetric system?
+       int a1 = _beam1.A();  
+       int a2 = _beam2.A();  
 
        static int IFIRSTH = 0;
-       static double DELL=0., DELR=0., SIGNN=0., R1=0., A1=0., R2=0., RHO1=0.;
-       static double RHO2=0., NZ1=0., NZ2=0., NR1=0., NR2=0.,RR1=0., NY=0., NX=0.;
+       static double DELL=0., DELR=0., SIGNN=0., R1=0., A1=0., A2=0., R2=0., RHO1=0.;
+       static double RHO2=0., NZ1=0., NZ2=0., NR1=0., NR2=0.,RR1=0., RR2=0., NY=0., NX=0.;
        static double AN1=0., AN2=0.;
-       double delo=0.,RSQ=0.,Z1=0.,Y=0.,X=0.,XB=0.,RPU=0.,IRUP=0.,RTU=0.;
+       double delo=0.,RSQ=0.,Z1=0.,Z2=0.,Y=0.,X=0.,XB=0.,RPU=0.,IRUP=0.,RTU=0.;
        double IRUT=0.,T1=0.,T2=0.;
        static double DEN1[20002], DEN2[20002];
        if (IFIRSTH != 0) goto L100;
@@ -308,21 +306,25 @@ beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
 
        //use two sigma_NN's. 52mb at rhic 100gev/beam, 88mb at LHC 2.9tev/beam, gamma is in cm system
        SIGNN = 5.2;
-       if ( gamma > 500. ) SIGNN = 8.8;
+       if ( gamma > 500. ) SIGNN = 8.8; 
        //use parameter from Constants
-       R1 = ( _beam1.nuclearRadius());  //remember _beam2? better way to do this generically
-       A1 = 0.535; //This is woodsaxonskindepth?
+       R1 = ( _beam1.nuclearRadius());  
+        R2 = ( _beam2.nuclearRadius());
+       A1 = 0.535; //This is woodsaxonskindepth
+        A2 = 0.535; 
        //write(6,12)r1,a1,signn  Here is where we could probably set this up asymmetrically R2=_beam2.nuclearRadius() and RHO2=ap2=_beam2.A()
-       R2 = R1;
-       RHO1 = ap;
-       RHO2 = RHO1;
+       // R2 = R1;
+       RHO1 = a1;
+       RHO2 = a2;
        NZ1  = ((R1+5.)/DELR);
        NR1  = NZ1;
        NZ2  = ((R2+5.)/DELR);
        NR2  = NZ2;
        RR1  = -DELR;
+        RR2  = -DELR; 
        NY   = ((R1+5.)/DELL);
        NX   = 2*NY;
+       // This calculates T_A(b) for beam 1 and stores it in DEN1[IR1] 
        for ( int IR1 = 1; IR1 <= NR1; IR1++) {
                DEN1[IR1] = 0.;
                RR1       = RR1+DELR;
@@ -332,21 +334,40 @@ beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
                        Z1  = Z1+DELR;
                        RSQ = RR1*RR1+Z1*Z1;
                        DEN1[IR1] = DEN1[IR1]+1./(1.+exp((sqrt(RSQ)-R1)/A1));
-               }//for(IZ)
+               }
 
                DEN1[IR1] = DEN1[IR1]*2.*DELR;
-               DEN2[IR1] = DEN1[IR1];
-       }//for(IR)
+       }
+
+       // This calculates T_A(b) for beam 2 and stores it in DEN2[IR2] 
+       for ( int IR2 = 1; IR2 <= NR2; IR2++) {
+               DEN2[IR2] = 0.;
+               RR2       = RR2+DELR;
+               Z2        = -DELR/2;
+
+               for ( int IZ2 = 1; IZ2 <= NZ2; IZ2++) {
+                       Z2  = Z2+DELR;
+                       RSQ = RR2*RR2+Z2*Z2;
+                       DEN2[IR2] = DEN2[IR2]+1./(1.+exp((sqrt(RSQ)-R2)/A2));
+               }
+
+               DEN2[IR2] = DEN2[IR2]*2.*DELR;
+       }
 
        AN1 = 0.;
        RR1 = 0.;
+        RR2 = 0.;
     
        for ( int IR1 =1; IR1 <= NR1; IR1++) {
                RR1 = RR1+DELR;
                AN1 = AN1+RR1*DEN1[IR1]*DELR*2.*3.141592654;
-       }//for(IR)
+       }
+       for ( int IR2 =1; IR2 <= NR2; IR2++) {
+               RR2 = RR2+DELR;
+               AN2 = AN2+RR2*DEN2[IR2]*DELR*2.*3.141592654;
+       }
         
-       AN2 = AN1; //This will also probably need to be changed?
+       // AN2 = AN1; //This will also probably need to be changed?
 
        delo = .05;
        //.1 to turn mb into fm^2
index 86fcb4e..029e1fd 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 176                         $: revision of last commit
-// $Author:: jseger                   $: author of last commit
-// $Date:: 2014-06-20 22:15:20 +0200 #$: date of last commit
+// $Rev:: 193                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -69,14 +69,12 @@ photonNucleusLuminosity::~photonNucleusLuminosity()
 //______________________________________________________________________________
 void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
 {
-       // double Av,Wgp,cs,cvma;
   double W,dW,dY;
   double Egamma,Y;
-  // double t,tmin,tmax;
   double testint,dndWdY;
   double csgA;
-  // double ax,bx;
   double C;  
+  int beam; 
 
   ofstream wylumfile;
   wylumfile.precision(15);
@@ -129,34 +127,41 @@ void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
   Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
                                                            +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
           (inputParametersInstance.protonEnergy()+sqrt(inputParametersInstance.protonEnergy()*
-                                              inputParametersInstance.protonEnergy()-starlightConstants::protonMass*starlightConstants::protonMass)));
-  
+                                              inputParametersInstance.protonEnergy()-starlightConstants::protonMass*starlightConstants::protonMass))); 
+
+  int A_1 = getbbs().beam1().A(); 
+  int A_2 = getbbs().beam2().A();
+
+  // Do this first for the case when the first beam is the photon emitter 
+  // Treat pA separately with defined beams 
+  // The variable beam (=1,2) defines which nucleus is the target 
   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 
     W = _wMin + double(i)*dW + 0.5*dW;
     
-    for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
+    for(unsigned int j = 0; j <= _nYbins - 1; ++j) { 
 
       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
 
-      int A_1 = getbbs().beam1().A(); 
-      int A_2 = getbbs().beam2().A();
       if( A_2 == 1 && A_1 != 1 ){
-        // pA, first beam is the nucleus 
-        Egamma = 0.5*W*exp(-Y);
+        // pA, first beam is the nucleus and is in this case the target  
+        Egamma = 0.5*W*exp(-Y); 
+        beam = 1; 
       } else if( A_1 ==1 && A_2 != 1){
-        // pA, second beam is the nucleus 
+        // pA, second beam is the nucleus and is in this case the target 
         Egamma = 0.5*W*exp(Y); 
+        beam = 2; 
       } else {
         Egamma = 0.5*W*exp(Y);        
+        beam = 2; 
       }
 
       dndWdY = 0.; 
 
       if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
 
-       csgA=getcsgA(Egamma,W);
-       dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);
+       csgA=getcsgA(Egamma,W,beam);
+        dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
 
       }
 
@@ -165,16 +170,49 @@ void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
     }
   }
 
+  // Repeat the loop for the case when the second beam is the photon emitter. 
+  // Don't repeat for pA
+  if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){ 
+    
+    for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
+
+      W = _wMin + double(i)*dW + 0.5*dW;
+    
+      for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
+
+        Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
+
+        beam = 1; 
+        Egamma = 0.5*W*exp(-Y);        
+
+        dndWdY = 0.; 
+        if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
+
+         csgA=getcsgA(Egamma,W,beam);
+          dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
+
+        }
+
+        wylumfile << dndWdY << endl;
+
+      }
+    }
+  }
+
+  wylumfile << bwnorm << endl;
+  wylumfile << inputParametersInstance.parameterValueKey() << endl;
   wylumfile.close();
   
   if(inputParametersInstance.interferenceEnabled()==1) 
     pttablegen();
  
-  wylumfile.open("slight.txt",ios::app);
-  cout << "bwnorm: "<< bwnorm <<endl;
-  wylumfile << bwnorm << endl;
-  wylumfile << inputParametersInstance.parameterValueKey() << endl;
-  wylumfile.close();
+  //  wylumfile.open("slight.txt",ios::app);
+  // cout << "bwnorm: "<< bwnorm <<endl;
+  //wylumfile << bwnorm << endl;
+  //wylumfile << inputParametersInstance.parameterValueKey() << endl;
+  //wylumfile.close();
+
 }
 
 
@@ -240,6 +278,7 @@ void photonNucleusLuminosity::pttablegen()
   double NPT = inputParametersInstance.nmbPtBinsInterference();
   double gamma_em=inputParametersInstance.beamLorentzGamma();
   double mass= getChannelMass();
+  int beam; 
   
   //  loop over y from 0 (not -ymax) to yma
   
@@ -256,6 +295,7 @@ void photonNucleusLuminosity::pttablegen()
     //  Find the sigma(gammaA) for the two directions
     //  Photonuclear Cross Section 1
     //  Gamma-proton CM energy
+    beam=2; 
     
     Wgp=sqrt(2.*Egamma1*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
                                 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
@@ -267,8 +307,7 @@ void photonNucleusLuminosity::pttablegen()
            /starlightConstants::alpha);
     
     // Calculate V.M.+nucleus cross section
-    
-    cvma=sigma_A(cs);
+    cvma=sigma_A(cs,beam);
     
     // Calculate Av = dsigma/dt(t=0) Note Units: fm**2/Gev**2
     
@@ -293,6 +332,7 @@ void photonNucleusLuminosity::pttablegen()
     sig_ga_1 = csgA;
           
     // Photonuclear Cross Section 2
+    beam=1; 
     
     Wgp=sqrt(2.*Egamma2*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
                                 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
@@ -300,7 +340,7 @@ void photonNucleusLuminosity::pttablegen()
     cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
            starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)/starlightConstants::alpha);
     
-    cvma=sigma_A(cs);
+    cvma=sigma_A(cs,beam);
     
     Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
                                              *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
index b95284d..576733a 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 179                         $: revision of last commit
+// $Rev:: 193                         $: revision of last commit
 // $Author:: jnystrand                $: author of last commit
-// $Date:: 2014-09-11 20:58:05 +0200 #$: date of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //    Added incoherent t2-> pt2 selection.  Following pp selection scheme
@@ -65,6 +65,8 @@ Gammaavectormeson::Gammaavectormeson(beamBeamSystem& bbsystem):eventChannel(bbsy
        _VMCoherenceFactor=inputParametersInstance.coherentProduction();
         _ProductionMode=inputParametersInstance.productionMode();
 
+        N1 = 0; N2 = 0; 
+
        switch(_VMpidtest){
        case starlightConstants::RHO:
        case starlightConstants::RHOZEUS:
@@ -131,7 +133,7 @@ Gammaavectormeson::~Gammaavectormeson()
 //______________________________________________________________________________
 void Gammaavectormeson::pickwy(double &W, double &Y)
 {
-       double dW, dY, xw,xy,xtest;
+        double dW, dY, xw,xy,xtest,btest;
        int  IW,IY;
   
        dW = (_VMWmax-_VMWmin)/double(_VMnumw);
@@ -153,7 +155,32 @@ void Gammaavectormeson::pickwy(double &W, double &Y)
 
        if( xtest > _Farray[IW][IY] )
                goto L201pwy;
-  
+
+       // Determine the target nucleus 
+       // For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2 
+        if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
+           if( _ProductionMode == 2 ){
+            _TargetBeam = 2;
+          } else {
+             _TargetBeam = 1;
+           }
+        } else if(  _bbs.beam1().A() != 1 && _bbs.beam2().A()==1 ){
+           if( _ProductionMode == 2 ){
+            _TargetBeam = 1;
+          } else {
+             _TargetBeam = 2;
+           }
+        } else {
+          btest = randyInstance.Rndom();
+          if ( btest < _Farray1[IW][IY]/_Farray[IW][IY] ){
+            _TargetBeam = 2;
+            N2++;
+          }  else {
+            _TargetBeam = 1;
+            N1++; 
+          }
+        }
+       // cout<<"N1: "<<N1<<" N2; "<<N2<<endl;
 }         
 
 
@@ -432,17 +459,18 @@ void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &p
            Egam = 0.5*W*exp(Y);
            Epom = 0.5*W*exp(-Y);
           }  
-        } else {
-         Egam = 0.5*W*exp(Y);
+        } else if ( _TargetBeam == 1 ) {
+         Egam = 0.5*W*exp(-Y);
+         Epom = 0.5*W*exp(Y);
+        } else { 
+          Egam = 0.5*W*exp(Y);
          Epom = 0.5*W*exp(-Y);
-        }
+       }
 
-        // cout<<" Y: "<<Y<<" W: "<<W<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl; 
+        // cout<<" Y: "<<Y<<" W: "<<W<<" TargetBeam; "<<_TargetBeam<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl; 
         pt1 = pTgamma(Egam);  
        phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
 
-       //      cout<<" pt1: "<<pt1<<endl; 
-
        if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) || 
             (_ProductionMode == 4) ) {
            if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
@@ -487,6 +515,8 @@ void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &p
                   }else{
                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();  
                   }  
+                } else if (_TargetBeam==1) {
+                    pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();  
                 } else {
                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();  
                 }
@@ -513,27 +543,18 @@ void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &p
                     }else{
                       comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
                     }  
+                  } else if (_TargetBeam==1) {
+                    comp = _bbs.beam2().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
                   } else {
-                    comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
+                    comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2; 
                   }
              
                   if( xtest > comp ) goto L203vm;
                        }
 
-               /*
-               if(_VMCoherence==0 && (!(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2))){
-                    //dsig/dt= exp(-_VMbslope*t)
-                    xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-                    t2 = (-1./_VMbslope)*log(xtest);
-                    pt2 = sqrt(1.*t2);
-               }
-               */
-
        }//else end from pp
        phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 
-       //        cout<<" pt2: "<<pt2<<endl;  
-
        px1 = pt1*cos(phi1);
        py1 = pt1*sin(phi1);
        px2 = pt2*cos(phi2);
@@ -549,16 +570,19 @@ void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &p
  
         // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl; 
 
-       // Randomly choose to make pz negative 50% of the time
+       // Keep this special case for d+A 
        if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
                pz = -pz;
-       }else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
+       }
+
+        /*
+        else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
          // Don't switch      
         }
        else{
                if (randyInstance.Rndom() >= 0.5) pz = -pz;
        }
-
+        */
 }
 
 //______________________________________________________________________________
@@ -591,6 +615,8 @@ double Gammaavectormeson::pTgamma(double E)
       }else{
          singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
       }  
+    } else if (_TargetBeam == 1) {
+      singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
     } else {
       singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
     }
@@ -616,6 +642,9 @@ double Gammaavectormeson::pTgamma(double E)
         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
       }  
+    } else if (_TargetBeam == 1) {
+        pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
+        singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
     } else {
         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
@@ -647,6 +676,9 @@ double Gammaavectormeson::pTgamma(double E)
                 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
                 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
               }  
+            } else if (_TargetBeam == 1) {
+              pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
+              singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
             } else {
               pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
               singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
@@ -664,8 +696,8 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
                              int&) // tcheck (unused)
 {
        //    This function calculates momentum and energy of vector meson
-       //     given W and Y, including interference.
-       //     It gets the pt distribution from a lookup table.
+       //    given W and Y, including interference.
+       //    It gets the pt distribution from a lookup table.
        double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
        int IY=0,j=0;
   
@@ -674,7 +706,7 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
   
        //  Y is already fixed; choose a pt
        //  Follow the approavh in pickwy.f
-       // in   _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
+       //  in  _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
   
        IY=int(fabs(Y)/dY);//+1;
        if (IY > (_VMnumy/2)-1){
@@ -683,17 +715,16 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
   
        yleft=fabs(Y)-(IY)*dY;
        yfract=yleft*dY;
-                                                                                                                                  
+  
        xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
-                                                                                                                                  
+        
        for(j=0;j<_VMNPT+1;j++){
                if (xpt < _fptarray[IY][j]) goto L60;
        }
  L60:
   
        //  now do linear interpolation - start with extremes
-  
-       if (j == 0){
+       if (j == 0){
                pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
                goto L80;
        }
@@ -703,8 +734,7 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
        }
   
        //  we're in the middle
-  
-       ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
+       ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
        pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
   
        //  at an extreme in y?
@@ -713,15 +743,14 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
                goto L120;
        }
  L80:
-       //  interpolate in y repeat for next fractional y bin
-                                                                                                                                  
+
+       //  interpolate in y repeat for next fractional y bin      
        for(j=0;j<_VMNPT+1;j++){
                if (xpt < _fptarray[IY+1][j]) goto L90;
        }
  L90:
   
        //  now do linear interpolation - start with extremes
-                                                                                                                                  
        if (j == 0){
                pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
                goto L100;
@@ -732,27 +761,20 @@ void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py,
        }
   
        //  we're in the middle
-                                                                                                                                  
        ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
        pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
-                                                                                                                                  
  L100:
-                                                                                                                                  
-       //  now interpolate in y
-                                                                                                                                  
+
+       //  now interpolate in y  
        pt=yfract*pt2+(1-yfract)*pt1;
-                                                                                                                                  
+
  L120:
-                                                                                                                                  
-       //  we have a pt
-                                                                                                                                  
+
+       //  we have a pt 
        theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
        px=pt*cos(theta);
        py=pt*sin(theta);
-                                                                                                                                  
-       //      I guess W is the mass of the vector meson (not necessarily
-       //      on-mass-shell), and E is the energy
-                                                                                                                                  
+
        E  = sqrt(W*W+pt*pt)*cosh(Y);
        pz = sqrt(W*W+pt*pt)*sinh(Y);
        //      randomly choose to make pz negative 50% of the time
index cdcb458..e4a751b 100644 (file)
@@ -385,17 +385,17 @@ double Gammagammaleptonpair::pp_2(double E)
 
 //______________________________________________________________________________
 void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
-                                    double  ,  // E (unused)
-                                    double  W,
-                                    double  px0, double  py0, double  pz0,
-                                    double& px1, double& py1, double& pz1,
-                                    double& px2, double& py2, double& pz2,
-                                    int&    iFbadevent)
+                                       double  ,  // E (unused)
+                                       double  W,
+                                       double  px0, double  py0, double  pz0,
+                                       double& px1, double& py1, double& pz1, double& E1,
+                                       double& px2, double& py2, double& pz2, double& E2,
+                                       double& mdec,
+                                       int&    iFbadevent)
 {
     //     This routine decays a particle into two particles of mass mdec,
     //     taking spin into account
 
-    double mdec=0.,E1=0.,E2=0.;
     double pmag, anglelep[20001];
     // double ytest=0.,dndtheta;
     double phi,theta,xtest,Ecm;
@@ -536,7 +536,8 @@ starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
     int iFbadevent=0;
     starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
        
-    double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
+    double px2=0.,px1=0.,py2=0.,pt2=0.,py1=0.,pz2=0.,pz1=0.,pt1=0.;
+    double E1=0.,E2=0.,mdec=0.;
 //this function decays particles and writes events to a file
     //zero out the event structure
     leptonpair._numberOfTracks=0;
@@ -554,32 +555,25 @@ starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
     picky(rapidity);
 
     pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
-    twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+    twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px1,px2,py2,E2,mdec,iFbadevent);
     if (iFbadevent==0){
-       int q1=0,q2=0; 
-
-       double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-       if (xtest<0.5)
-       {
-           q1=1;
-           q2=-1;
-       }
-       else{
-           q1=-1;
-           q2=1;
-       }       
+      const bool xtest(randyInstance.Rndom() < 0.5);
+      const bool charges[2] = {
+       ( xtest) ? +1 : -1,
+       (!xtest) ? -1 : +1
+      };
        leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
        leptonpair.px[0]=px1;
        leptonpair.py[0]=py1;
        leptonpair.pz[0]=pz1;
        leptonpair._fsParticle[0]=ipid; 
-       leptonpair._charge[0]=q1;
+       leptonpair._charge[0]=charges[0];
 
        leptonpair.px[1]=px2;
        leptonpair.py[1]=py2;
        leptonpair.pz[1]=pz2;
        leptonpair._fsParticle[1]=ipid;
-       leptonpair._charge[1]=q2;
+       leptonpair._charge[1]=charges[1];
 
        ievent=ievent+1;
     }
@@ -602,7 +596,7 @@ upcEvent Gammagammaleptonpair::produceEvent()
    int iFbadevent=0;
    starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
    
-   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.;
+   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.,mdec=0.;
    bool accepted = false;
    do{ 
      //this function decays particles and writes events to a file
@@ -615,7 +609,7 @@ upcEvent Gammagammaleptonpair::produceEvent()
    
   
      _nmbAttempts++;
-     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px2,py2,pz2,E2,mdec,iFbadevent);
      double pt1chk = sqrt(px1*px1+py1*py1);
      double pt2chk = sqrt(px2*px2+py2*py2);
     
@@ -649,30 +643,18 @@ upcEvent Gammagammaleptonpair::produceEvent()
    //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
    
    if (iFbadevent==0){
-     int q1=0,q2=0; 
-     
-     double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-     if (xtest<0.5)
-       {
-        q1=1;
-        q2=-1;
-       }
-     else{
-       q1=-1;
-       q2=1;
-     }
+     const bool xtest(randyInstance.Rndom() < 0.5);
+     const int charges[2] = {
+       ( xtest) ? +1 : -1,
+       (!xtest) ? -1 : +1
+     };
 
      // The new stuff
-     double mlepton = getMass(); 
-     E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 ); 
-     E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 ); 
-
-     starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1);
+     starlightParticle particle1(px1, py1, pz1, E1, mdec, -charges[0]*ipid, charges[0]);
      event.addParticle(particle1);
      
-     starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2);
-     event.addParticle(particle2);
-     
+     starlightParticle particle2(px2, py2, pz2, E2, mdec, -charges[1]*ipid, charges[1]);
+     event.addParticle(particle2);     
     }
    return event;
 }
index 45875ff..fa507db 100644 (file)
@@ -327,12 +327,11 @@ double Gammagammasingle::pp(double E)
 
 
 //______________________________________________________________________________
-void Gammagammasingle::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double /*E*/,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &px2,double &py2,double &pz2,int &iFbadevent)
+void Gammagammasingle::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double /*E*/,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2, double &mdec,int &iFbadevent)
 {
   //     This routine decays a particle into two particles of mass mdec,
   //     taking spin into account
   
-  double mdec=0.,E1=0.,E2=0.;
   double pmag,ytest=0.;
   double phi,theta,xtest,dndtheta,Ecm;
   double  betax,betay,betaz;
@@ -450,34 +449,19 @@ starlightConstants::event Gammagammasingle::produceEvent(int &/*ievent*/)
 //starlightConstants::event Gammagammasingle::produceEvent(int &ievent)
 upcEvent Gammagammasingle::produceEvent()
 {
-  //    cout << "NOT IMPLEMENTED!" << endl;
-        
-  //    return upcEvent();
-
+  upcEvent event;
+  
   //    returns the vector with the decay particles inside.
   //   onedecayparticle single;
-  starlightConstants::event single;
   double comenergy = 0.;
   double rapidity = 0.;
   double parentE = 0.;
   double parentmomx=0.,parentmomy=0.,parentmomz=0.;
 
-  //this function decays particles and writes events to a file
-  //zeroing out the event structure
-  single._numberOfTracks=0;
-  for(int i=0;i<4;i++){
-    single.px[i]=0.;
-    single.py[i]=0.;
-    single.pz[i]=0.;
-    single._fsParticle[i]=starlightConstants::UNKNOWN;
-    single._charge[i]=0;
-  }
-  
   pickw(comenergy);
   picky(rapidity);
   parentMomentum(comenergy,rapidity,parentE,parentmomx,parentmomy,parentmomz);
-  
-  
+    
   if(_GGsingInputpidtest != starlightConstants::F2 && _GGsingInputpidtest != starlightConstants::F2PRIME)
   {
 #ifdef ENABLE_PYTHIA
@@ -490,104 +474,83 @@ upcEvent Gammagammasingle::produceEvent()
   }
 
 
-  int ievent = 0;
   int iFbadevent=0;
-  starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
-  double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
-  double px3=0.,px4=0.,py3=0.,py4=0.,pz3=0.,pz4=0.;
   //  double theta=0.,phi=0.;//angles from jetset
-  double xtest=0.,ztest=0.;
   switch(_GGsingInputpidtest){
-  case starlightConstants::ZOVERZ03:
+  case starlightConstants::ZOVERZ03: {
     //Decays into two pairs.
-    parentmomx=parentmomx/2.;
-    parentmomy=parentmomy/2.;
-    parentmomz=parentmomz/2.;
+    parentmomx /= 2.;
+    parentmomy /= 2.;
+    parentmomz /= 2.;
+    starlightConstants::particleTypeEnum ipid[2] = { starlightConstants::UNKNOWN };
+    double mdec[2] = { 0 };
+    double px[4] = { 0 };
+    double py[4] = { 0 };
+    double pz[4] = { 0 };
+    double pt[4] = { 0 };
     //Pair #1  
-    twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+    twoBodyDecay(ipid[0],parentE,comenergy,
+                parentmomx,parentmomy,parentmomz,
+                px[0],py[0],pz[0],pt[0],
+                px[1],py[1],pz[1],pt[1],
+                mdec[0],iFbadevent);
     //Pair #2
-    twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px3,py3,pz3,px4,py4,pz4,iFbadevent);
+    twoBodyDecay(ipid[1],parentE,comenergy,
+                parentmomx,parentmomy,parentmomz,
+                px[2],py[2],pz[2],pt[2],
+                px[3],py[3],pz[3],pt[3],
+                mdec[1],iFbadevent);
     //Now add them to vectors to be written out later.
                
-    single._numberOfTracks=4;//number of tracks per event
     if (iFbadevent==0){
-      xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-      ztest = randyInstance.Rndom();
       //Assigning charges randomly.
-      if (xtest<0.5){
-       single._charge[0]=1;//q1=1;
-       single._charge[1]=-1;//q2=-1;
-      }
-      else{
-       single._charge[0]=1;//q1=-1;
-       single._charge[1]=-1;//q2=1;
-      }
-      if (ztest<0.5){
-       single._charge[2]=1;//q3=1;
-       single._charge[3]=-1;//q4=-1;
-      }
-      else{
-       single._charge[2]=-1;//q3=-1;
-       single._charge[3]=1;//q4=1;
-      }
-      //Track #1
-      single.px[0]=px1;
-      single.py[0]=py1;
-      single.pz[0]=pz1;
-      single._fsParticle[0]=ipid;
-      //Track #2                                                                                                                      
-      single.px[1]=px2;
-      single.py[1]=py2;
-      single.pz[1]=pz2;
-      single._fsParticle[1]=ipid;
-      //Track #3
-      single.px[2]=px3;
-      single.py[2]=py3;
-      single.pz[2]=pz3;
-      single._fsParticle[2]=ipid;
-      //Track #4
-      single.px[3]=px4;
-      single.py[3]=py4;
-      single.pz[3]=pz4;
-      single._fsParticle[3]=ipid;
-      
-      ievent=ievent+1;
-    }  
-    
+      const bool xtest(randyInstance.Rndom() < 0.5);
+      const bool ztest(randyInstance.Rndom() < 0.5);
+      const int charges[4] = {
+       ( xtest) ? +1 : -1,
+       (!xtest) ? -1 : +1,
+       ( ztest) ? +1 : -1,
+       (!ztest) ? -1 : +1,
+      };
+      for (int i=0; i<4; ++i) {
+       starlightParticle particle(px[i], py[i], pz[i], pt[i], mdec[i/2], charges[i]*ipid[i/2], charges[i]);    
+       event.addParticle(particle);
+      } 
+    }
     break;
+  }
   case starlightConstants::F2:
-  case starlightConstants::F2PRIME:
-    twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+  case starlightConstants::F2PRIME: {
+    double px[2] = { 0 };
+    double py[2] = { 0 };
+    double pz[2] = { 0 };
+    double pt[2] = { 0 };
+    double mdec  = 0;
+    starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
+    twoBodyDecay(ipid,parentE,comenergy,
+                parentmomx,parentmomy,parentmomz,
+                px[0],py[0],pz[0],pt[0],
+                px[1],py[1],pz[1],pt[1],
+                mdec,iFbadevent);
     
-    single._numberOfTracks=2;
-    if (iFbadevent==0){
-      xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-      if (xtest<0.5){
-       single._charge[0]=1;//q1=1;
-       single._charge[1]=-1;//q2=-1;
+    if (iFbadevent==0) {
+      const double xtest = randyInstance.Rndom();
+      const int charges[2] = {
+       (xtest < 0.5) ? +1 : -1,
+       (xtest < 0.5) ? -1 : +1,
+      };
+      for (int i=0; i<2; ++i) {
+       starlightParticle particle(px[i], py[i], pz[i], pt[i], mdec, charges[i]*ipid, charges[i]);
+       event.addParticle(particle);
       }
-      else{
-       single._charge[0]=-1;//q1=-1;
-       single._charge[1]=1;//q2=1;
-      }        
-      //Track #1
-      single.px[0]=px1;
-      single.py[0]=py1;
-      single.pz[0]=pz1;
-      single._fsParticle[0]=ipid*single._charge[0]; 
-      //Track #2
-      single.px[1]=px2;
-      single.py[1]=py2;
-      single.pz[1]=pz2;
-      single._fsParticle[1]=ipid*single._charge[1];
-      ievent=ievent+1;
     }
     break;
+  }
   default:
     break;
   }
-  
-  return upcEvent(single);
+
+  return event;
 }
 
 
index d092734..6f72c7c 100644 (file)
@@ -66,13 +66,11 @@ incoherentPhotonNucleusLuminosity::~incoherentPhotonNucleusLuminosity()
 //______________________________________________________________________________
 void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLuminosity()
 {
-       // double Av,Wgp,cs,cvma;
   double W,dW,dY;
   double Egamma,Y;
-  // double t,tmin,tmax;
   double testint,dndWdY;
-  // double ax,bx;
   double C;  
+  int beam; 
 
   ofstream wylumfile;
   wylumfile.precision(15);
@@ -121,11 +119,13 @@ void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLumin
     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
     wylumfile << Y << endl;
   }
-  //  Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
-  //                                                       +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
-  //      (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
-  
+
+  int A_1 = getbbs().beam1().A(); 
+  int A_2 = getbbs().beam2().A();
+
+  // Do this first for the case when the first beam is the photon emitter 
+  // Treat pA separately with defined beams 
+  // The variable beam (=1,2) defines which nucleus is the target 
   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 
     W = _wMin + double(i)*dW + 0.5*dW;
@@ -139,16 +139,17 @@ void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLumin
 
       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
 
-      int A_1 = getbbs().beam1().A(); 
-      int A_2 = getbbs().beam2().A();
       if( A_2 == 1 && A_1 != 1 ){
-        // pA, first beam is the nucleus 
+        // pA, first beam is the nucleus and photon emitter
         Egamma = 0.5*W*exp(Y);
+        beam = 2; 
       } else if( A_1 ==1 && A_2 != 1){
-        // pA, second beam is the nucleus 
+        // pA, second beam is the nucleus and photon emitter
         Egamma = 0.5*W*exp(-Y); 
+        beam = 1; 
       } else {
         Egamma = 0.5*W*exp(Y);        
+        beam = 2; 
       }
       
       dndWdY = 0.; 
@@ -165,17 +166,17 @@ void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLumin
           // localbmin = getbbs().beam2().nuclearRadius() + 0.7; 
           // localz = getbbs().beam2().Z(); 
          //   dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm); 
-          dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm); 
+          dndWdY = Egamma*photonFlux(Egamma,beam)*localsig*breitWigner(W,bwnorm); 
         }else if (A_2 ==1 && A_1 !=1){
           // localbmin = getbbs().beam1().nuclearRadius() + 0.7; 
           // localz = getbbs().beam1().Z(); 
          //   dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm); 
-          dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm); 
+          dndWdY = Egamma*photonFlux(Egamma,beam)*localsig*breitWigner(W,bwnorm); 
         }else{ 
           double csVN = sigma_N(Wgp);         
-          double csVA = sigma_A(csVN); 
+          double csVA = sigma_A(csVN,beam); 
           double csgA= (csVA/csVN)*sigmagp(Wgp); 
-          dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm); 
+          dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm); 
         }
       }
 
@@ -184,12 +185,49 @@ void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLumin
     }
   }
 
+  // Repeat the loop for the case when the second beam is the photon emitter. 
+  // Don't repeat for pA
+  if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){ 
+    for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
+
+      W = _wMin + double(i)*dW + 0.5*dW;
+
+      double Ep = inputParametersInstance.protonEnergy();
+
+      Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
+          (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
+    
+      for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
+
+        Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
+
+        beam = 1; 
+        Egamma = 0.5*W*exp(-Y);        
+      
+        dndWdY = 0.; 
+
+        if(Egamma > Eth){
+         if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy();
+          double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
+                                 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
+
+          double csVN = sigma_N(Wgp);         
+          double csVA = sigma_A(csVN,beam); 
+          double csgA= (csVA/csVN)*sigmagp(Wgp); 
+          dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm); 
+        
+        }
+
+        wylumfile << dndWdY << endl;
+
+      }
+    }
+  }
+
   wylumfile << bwnorm << endl;
   wylumfile << inputParametersInstance.parameterValueKey() << endl;
   wylumfile.close();
   
-  // cout << "bwnorm: "<< bwnorm <<endl;
-
 }
 
 
index 8c0c57d..f119b51 100644 (file)
@@ -66,19 +66,15 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
        // This subroutine calculates the vector meson cross section assuming
        // a narrow resonance.  For reference, see STAR Note 386.
   
-       // double Av,Wgp,cs,cvma;
        double W,dY;
        double y1,y2,y12,ega1,ega2,ega12;
         // double y1lab,y2lab,y12lab;
-       // double t,tmin,tmax;
-       double csgA1,csgA2,csgA12,int_r,dR,rate;
+       double csgA1,csgA2,csgA12,int_r,dR;
         double Wgp,csVN,csVA; 
         double Wgpm; 
        double tmp;
-       // double ax,bx;
        double Eth;
-       int    J,NY;
-       // int    K,NGAUSS;
+       int    J,NY,beam;
   
        NY   =  _narrowNumY;
        dY   = (_narrowYmax-_narrowYmin)/double(NY);
@@ -94,6 +90,12 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
   
        tmp = 0.0;
   
+        int A_1 = getbbs().beam1().A(); 
+        int A_2 = getbbs().beam2().A();
+
+        // Do this first for the case when the first beam is the photon emitter 
+        // Treat pA separately with defined beams 
+        // The variable beam (=1,2) defines which nucleus is the target 
        for(J=0;J<=(NY-1);J++){
     
                // This is the fdefault
@@ -101,9 +103,31 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
                y2  = _narrowYmin + double(J+1)*dY;
                y12 = 0.5*(y1+y2);
     
-               ega1  = 0.5*W*exp(y1);
-               ega2  = 0.5*W*exp(y2);
-               ega12 = 0.5*W*exp(y12);
+                if( A_2 == 1 && A_1 != 1 ){
+                  // pA, first beam is the nucleus and photon emitter
+                  // Egamma = 0.5*W*exp(Y);
+                 ega1  = 0.5*W*exp(y1);
+                 ega2  = 0.5*W*exp(y2);
+                 ega12 = 0.5*W*exp(y12);
+                  beam = 2; 
+                } else if( A_1 ==1 && A_2 != 1){
+                  // pA, second beam is the nucleus and photon emitter
+                  // Egamma = 0.5*W*exp(-Y); 
+                 ega1  = 0.5*W*exp(-y1);
+                 ega2  = 0.5*W*exp(-y2);
+                 ega12 = 0.5*W*exp(-y12);
+                  beam = 1; 
+                } else {
+                  // Egamma = 0.5*W*exp(Y);        
+                 ega1  = 0.5*W*exp(y1);
+                 ega2  = 0.5*W*exp(y2);
+                 ega12 = 0.5*W*exp(y12);
+                  beam = 2; 
+                }
+
+               // ega1  = 0.5*W*exp(y1);
+               // ega2  = 0.5*W*exp(y2);
+               // ega12 = 0.5*W*exp(y12);
 
                 // This is for checking things in the lab frame 
                 // y1lab  = _narrowYmin + double(J)*dY;
@@ -136,7 +160,7 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
                 Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
                                  +starlightConstants::protonMass*starlightConstants::protonMass);
                 csVN = sigma_N(Wgp);            
-                csVA = sigma_A(csVN); 
+                csVA = sigma_A(csVN,beam); 
                 csgA1 = (csVA/csVN)*sigmagp(Wgp); 
                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
                   csgA1 = sigmagp(Wgp);
@@ -147,7 +171,7 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
                                  +starlightConstants::protonMass*starlightConstants::protonMass);
                 Wgpm = Wgp; 
                 csVN = sigma_N(Wgp);            
-                csVA = sigma_A(csVN); 
+                csVA = sigma_A(csVN,beam); 
                 csgA12 = (csVA/csVN)*sigmagp(Wgp); 
                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
                   csgA12 = sigmagp(Wgp);
@@ -157,28 +181,118 @@ incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave
                 Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
                                  +starlightConstants::protonMass*starlightConstants::protonMass);
                 csVN = sigma_N(Wgp);            
-                csVA = sigma_A(csVN); 
+                csVA = sigma_A(csVN,beam); 
                 csgA2 = (csVA/csVN)*sigmagp(Wgp); 
                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
                   csgA2 = sigmagp(Wgp);
                 }
 
-                dR = ega1*photonFlux(ega1)*csgA1;  
-                dR = dR + 4*ega12*photonFlux(ega12)*csgA12;
-                dR = dR + ega2*photonFlux(ega2)*csgA2; 
+                dR = ega1*photonFlux(ega1,beam)*csgA1;  
+                dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
+                dR = dR + ega2*photonFlux(ega2,beam)*csgA2; 
                 dR = dR*(dY/6.); 
 
                // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
                 // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
 
                // The 2 accounts for the 2 beams    
-               if(getbbs().beam1().A()==getbbs().beam2().A()){
-                       dR  = 2.*dR;
-               }
+               //              if(getbbs().beam1().A()==getbbs().beam2().A()){
+               //      dR  = 2.*dR;
+               //}
 
                int_r = int_r+dR;
 
        }
-       rate=luminosity()*int_r;
+
+        // Repeat the loop for the case when the second beam is the photon emitter. 
+        // Don't repeat for pA
+        if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){ 
+         for(J=0;J<=(NY-1);J++){
+    
+               // This is the fdefault
+               y1  = _narrowYmin + double(J)*dY;
+               y2  = _narrowYmin + double(J+1)*dY;
+               y12 = 0.5*(y1+y2);
+    
+                beam = 1; 
+               ega1  = 0.5*W*exp(-y1);
+               ega2  = 0.5*W*exp(-y2);
+               ega12 = 0.5*W*exp(-y12);
+
+                // This is for checking things in the lab frame 
+                // y1lab  = _narrowYmin + double(J)*dY;
+                // y2lab  = _narrowYmin + double(J+1)*dY;
+                // y12lab = 0.5*(y1lab+y2lab);
+
+                // p+Pb
+                // y1 = y1lab + 0.465;
+                // y2 = y2lab + 0.465;
+                // y12 = y12lab + 0.465; 
+                // ega1  = 0.5*W*exp(y1);
+                // ega2  = 0.5*W*exp(y2);
+                // ega12 = 0.5*W*exp(y12);
+
+                // Pb+p
+                // y1 = y1lab - 0.465;
+                // y2 = y2lab - 0.465;
+                // y12 = y12lab - 0.465; 
+                // ega1  = 0.5*W*exp(-y1);
+                // ega2  = 0.5*W*exp(-y2);
+                // ega12 = 0.5*W*exp(-y12);
+
+    
+               if(ega2 < Eth)   
+                       continue;
+               if(ega1 > maxPhotonEnergy()) 
+                       continue;
+
+               // First point 
+                Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
+                                 +starlightConstants::protonMass*starlightConstants::protonMass);
+                csVN = sigma_N(Wgp);            
+                csVA = sigma_A(csVN,beam); 
+                csgA1 = (csVA/csVN)*sigmagp(Wgp); 
+                if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
+                  csgA1 = sigmagp(Wgp);
+                }
+
+               // Middle point 
+                Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
+                                 +starlightConstants::protonMass*starlightConstants::protonMass);
+                Wgpm = Wgp; 
+                csVN = sigma_N(Wgp);            
+                csVA = sigma_A(csVN,beam); 
+                csgA12 = (csVA/csVN)*sigmagp(Wgp); 
+                if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
+                  csgA12 = sigmagp(Wgp);
+                }
+
+               // Last point 
+                Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
+                                 +starlightConstants::protonMass*starlightConstants::protonMass);
+                csVN = sigma_N(Wgp);            
+                csVA = sigma_A(csVN,beam); 
+                csgA2 = (csVA/csVN)*sigmagp(Wgp); 
+                if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
+                  csgA2 = sigmagp(Wgp);
+                }
+
+                dR = ega1*photonFlux(ega1,beam)*csgA1;  
+                dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
+                dR = dR + ega2*photonFlux(ega2,beam)*csgA2; 
+                dR = dR*(dY/6.); 
+
+               // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
+                // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
+
+               // The 2 accounts for the 2 beams    
+               // if(getbbs().beam1().A()==getbbs().beam2().A()){
+               //      dR  = 2.*dR;
+               //}
+
+               int_r = int_r+dR;
+
+         }
+        }
        cout<<" Cross section (mb): " <<10.*int_r<<endl;
 }
index f238c5c..8d4529e 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+// $Rev:: 193                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -66,16 +66,11 @@ narrowResonanceCrossSection::crossSectionCalculation(const double)  // _bwnormsa
        // This subroutine calculates the vector meson cross section assuming
        // a narrow resonance.  For reference, see STAR Note 386.
   
-       // double Av,Wgp,cs,cvma;
        double W,dY;
        double y1,y2,y12,ega1,ega2,ega12;
-       // double t,tmin,tmax;
-       double csgA1,csgA2,csgA12,int_r,dR,rate;
-       double tmp;
-       // double ax,bx;
+       double csgA1,csgA2,csgA12,int_r,dR;
        double Eth;
-       int    J,NY;
-       // int    K,NGAUSS;
+       int    J,NY,beam;
   
        NY   =  _narrowNumY;
        dY   = (_narrowYmax-_narrowYmin)/double(NY);
@@ -88,47 +83,106 @@ narrowResonanceCrossSection::crossSectionCalculation(const double)  // _bwnormsa
   
        cout<<" gamma+nucleon  Threshold: "<<Eth<<endl;
        int_r=0.;
+
+        int A_1 = getbbs().beam1().A(); 
+        int A_2 = getbbs().beam2().A();
   
-       tmp = 0.0;
-  
+        // Do this first for the case when the first beam is the photon emitter 
+        // Treat pA separately with defined beams 
+        // The variable beam (=1,2) defines which nucleus is the target 
        for(J=0;J<=(NY-1);J++){
     
                y1  = _narrowYmin + double(J)*dY;
                y2  = _narrowYmin + double(J+1)*dY;
                y12 = 0.5*(y1+y2);
     
-               ega1  = 0.5*W*exp(y1);
-               ega2  = 0.5*W*exp(y2);
-               ega12 = 0.5*W*exp(y12);
+                if( A_2 == 1 && A_1 != 1 ){
+                  // pA, first beam is the nucleus and is in this case the target  
+                  // Egamma = 0.5*W*exp(-Y); 
+                  ega1  = 0.5*W*exp(-y1);
+                  ega2  = 0.5*W*exp(-y2);
+                  ega12 = 0.5*W*exp(-y12);
+                  beam = 1; 
+                } else if( A_1 ==1 && A_2 != 1){
+                  // pA, second beam is the nucleus and is in this case the target 
+                  // Egamma = 0.5*W*exp(Y); 
+                  ega1  = 0.5*W*exp(y1);
+                  ega2  = 0.5*W*exp(y2);
+                  ega12 = 0.5*W*exp(y12);
+                  beam = 2; 
+                } else {
+                  // Egamma = 0.5*W*exp(Y);        
+                  ega1  = 0.5*W*exp(y1);
+                  ega2  = 0.5*W*exp(y2);
+                  ega12 = 0.5*W*exp(y12);
+                  beam = 2; 
+                }
+
+               //              ega1  = 0.5*W*exp(y1);
+               //              ega2  = 0.5*W*exp(y2);
+               //              ega12 = 0.5*W*exp(y12);
     
                if(ega1 < Eth)   
                        continue;
                if(ega2 > maxPhotonEnergy()) 
                        continue;
 
-               csgA1=getcsgA(ega1,W);
+               csgA1=getcsgA(ega1,W,beam);
     
                // Middle Point                      =====>>>
-               csgA12=getcsgA(ega12,W); 
+               csgA12=getcsgA(ega12,W,beam); 
 
                // Second Point                      =====>>>
-               csgA2=getcsgA(ega2,W);
+               csgA2=getcsgA(ega2,W,beam);
     
                // Sum the contribution for this W,Y.
-               dR  = ega1*photonFlux(ega1)*csgA1;
-               dR  = dR + 4.*ega12*photonFlux(ega12)*csgA12;
-               dR  = dR + ega2*photonFlux(ega2)*csgA2;
-               tmp = tmp+2.*dR*(dY/6.);
+               dR  = ega1*photonFlux(ega1,beam)*csgA1;
+               dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
+               dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
                dR  = dR*(dY/6.);
 
-               // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
+               // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12,beam)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
 
-               // The 2 accounts for the 2 beams
-               if(getbbs().beam1().A()==getbbs().beam2().A()){
-                       dR  = 2.*dR;
-               }
                int_r = int_r+dR;
        }
-       rate=luminosity()*int_r;
+
+        // Repeat the loop for the case when the second beam is the photon emitter. 
+        // Don't repeat for pA
+        if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){ 
+         for(J=0;J<=(NY-1);J++){
+    
+               y1  = _narrowYmin + double(J)*dY;
+               y2  = _narrowYmin + double(J+1)*dY;
+               y12 = 0.5*(y1+y2);
+    
+                beam = 1; 
+               ega1  = 0.5*W*exp(-y1);
+               ega2  = 0.5*W*exp(-y2);
+               ega12 = 0.5*W*exp(-y12);
+      
+               if(ega2 < Eth)   
+                       continue;
+               if(ega1 > maxPhotonEnergy()) 
+                       continue;
+
+               csgA1=getcsgA(ega1,W,beam);
+    
+               // Middle Point                      =====>>>
+               csgA12=getcsgA(ega12,W,beam); 
+
+               // Second Point                      =====>>>
+               csgA2=getcsgA(ega2,W,beam);
+    
+               // Sum the contribution for this W,Y.
+               dR  = ega1*photonFlux(ega1,beam)*csgA1;
+               dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
+               dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
+               dR  = dR*(dY/6.);
+
+               // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12,beam)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
+
+               int_r = int_r+dR;
+         }
+        }
        cout<<" Cross section (mb): " <<10.*int_r<<endl;
 }
index 8c281cf..044407b 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 189                         $: revision of last commit
-// $Author:: srklein                  $: author of last commit
-// $Date:: 2014-10-27 04:29:14 +0100 #$: date of last commit
+// $Rev:: 193                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -61,39 +61,8 @@ photonNucleusCrossSection::photonNucleusCrossSection(const beamBeamSystem& bbsys
           _productionMode    (inputParametersInstance.productionMode()    ),
          _sigmaNucleus      (_bbs.beam2().A()          )
 {
-       // define luminosity for various beam particles in units of 10^{26} cm^{-2} sec^{-1}
-       switch(_bbs.beam1().Z()) {
-       case 1:   // proton
-               _luminosity = 1.E8;
-               break;
-       case 8:   // O
-               _luminosity = 980.;
-               break;
-       case 14:  // Si
-               _luminosity = 440.;
-               break;
-       case 20:  // Ca
-               _luminosity = 2000.;
-               break;
-       case 29:  // Cu
-               _luminosity = 95.;
-               break;
-       case 49:  // Indium, uses same as Iodine
-               _luminosity = 27.;
-               break;
-       case 53:  // I
-               _luminosity = 27.;
-               break;
-       case 79:  // Au
-               _luminosity = 2.0;
-               break;
-       case 82:  // Pb
-               _luminosity = 1.;
-               break;
-       default:
-               printWarn << "luminosity is not defined for beam with Z = " << _bbs.beam1().Z()
-                         << ". using " << _luminosity << " 10^{26} cm^{-2} sec^{-1}" << endl;
-       }
+       // The _luminosity variable should note be used. Set it to 1.0 here. 
+        _luminosity = 1.0; 
 
        switch(_particleType) {
        case RHO:
@@ -201,8 +170,8 @@ photonNucleusCrossSection::photonNucleusCrossSection(const beamBeamSystem& bbsys
                     <<" GammaAcrosssection"<<endl;
        }
 
-        double maxradius = (_bbs.beam1().nuclearRadius()<_bbs.beam2().nuclearRadius())?_bbs.beam2().nuclearRadius():_bbs.beam1().nuclearRadius();
-       _maxPhotonEnergy = 5. * _beamLorentzGamma * hbarc/maxradius; 
+        // double maxradius = (_bbs.beam1().nuclearRadius()<_bbs.beam2().nuclearRadius())?_bbs.beam2().nuclearRadius():_bbs.beam1().nuclearRadius();
+       _maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius()); 
 }
 
 
@@ -222,7 +191,8 @@ photonNucleusCrossSection::crossSectionCalculation(const double)
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::getcsgA(const double Egamma,
-                                   const double W)
+                                   const double W, 
+                                   const int beam)
 {
        //This function returns the cross-section for photon-nucleus interaction 
        //producing vectormesons
@@ -230,7 +200,7 @@ photonNucleusCrossSection::getcsgA(const double Egamma,
        double Av,Wgp,cs,cvma;
        double t,tmin,tmax;
        double csgA,ax,bx;
-       int NGAUSS;                                                                                                                                       
+       int NGAUSS; 
   
        //     DATA FOR GAUSS INTEGRATION
        double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
@@ -266,20 +236,19 @@ photonNucleusCrossSection::getcsgA(const double Egamma,
                }
                csgA = 0.5 * (tmax - tmin) * csgA;
                csgA = Av * csgA;
-       }       else if (!_coherentProduction &&
-                        (!((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)))) {  // incoherent AA interactions
+               //      }       else if (!_coherentProduction &&
+               //                       (!((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)))) {  // incoherent AA interactions
                // For incoherent AA interactions, since incoherent treating it as gamma-p
                // Calculate the differential V.M.+proton cross section
-               csgA = 1.E-4 * _incoherentFactor * _sigmaNucleus * _slopeParameter * sigmagp(Wgp);
-
+               //              csgA = 1.E-4 * _incoherentFactor * _sigmaNucleus * _slopeParameter * sigmagp(Wgp);
        }       else {  // coherent AA interactions
                // For typical AA interactions.
                // Calculate V.M.+proton cross section
-               cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
+               cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha); 
     
                // Calculate V.M.+nucleus cross section
                // cvma = _bbs.beam1().A()*cs; 
-               cvma = sigma_A(cs); 
+               cvma = sigma_A(cs,beam); 
 
                // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
                Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
@@ -300,7 +269,13 @@ photonNucleusCrossSection::getcsgA(const double Egamma,
                         }else if(A_2 ==1 && A_1 != 1){
                          csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                         }else{     
-                         csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
+                          if( beam==1 ){
+                           csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
+                          }else if(beam==2){
+                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);      
+                         }else{
+                           cout<<"Something went wrong here, beam= "<<beam<<endl; 
+                          }
                         }
 
                        t    = ax * (-xg[k]) + bx;
@@ -309,7 +284,13 @@ photonNucleusCrossSection::getcsgA(const double Egamma,
                         }else if(A_2 ==1 && A_1 != 1){
                          csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                         }else{     
-                         csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
+                          if( beam==1 ){
+                           csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
+                          }else if(beam==2){
+                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);      
+                         }else{
+                           cout<<"Something went wrong here, beam= "<<beam<<endl; 
+                          }
                         }
                }
                csgA = 0.5 * (tmax - tmin) * csgA;
@@ -322,7 +303,7 @@ photonNucleusCrossSection::getcsgA(const double Egamma,
 
 //______________________________________________________________________________
 double
-photonNucleusCrossSection::photonFlux(const double Egamma)
+photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
 {
        // This routine gives the photon flux as a function of energy Egamma
        // It works for arbitrary nuclei and gamma; the first time it is
@@ -331,27 +312,36 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
        // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
        // energies are in GeV, in the lab frame
        // rewritten 4/25/2001 by SRK
-  
+
+        // NOTE: beam (=1,2) defines the photon TARGET
+
        double lEgamma,Emin,Emax;
        static double lnEmax, lnEmin, dlnE;
-       double stepmult,energy,rZ,rA;
+       // double stepmult,energy,rZ,rA;
+       double stepmult,energy,rZ;
        int nbstep,nrstep,nphistep,nstep;
        double bmin,bmax,bmult,biter,bold,integratedflux;
        double fluxelement,deltar,riter;
        double deltaphi,phiiter,dist;
        static double dide[401];
        double lnElt;
-       double rA2, rZ2; 
+       // double rA2, rZ2; 
        double flux_r; 
        double Xvar;
        int Ilt;
-       double RNuc=0.,RNuc2=0.,maxradius=0.;
+       double RNuc=0.,RSum=0.;
+
+       RSum=_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius();
+        if( beam == 1){
+          rZ=double(_bbs.beam2().Z());
+          RNuc = _bbs.beam1().nuclearRadius();
+        } else { 
+         rZ=double(_bbs.beam1().Z());
+          RNuc = _bbs.beam2().nuclearRadius();
+        }
 
-       RNuc=_bbs.beam1().nuclearRadius();
-       RNuc2=_bbs.beam2().nuclearRadius();
-        maxradius =  (RNuc<RNuc2)?RNuc2:RNuc;
-       // static ->>> dide,lnEMax,lnEmin,dlnE
        static int  Icheck = 0;
+       static int  Ibeam  = 0; 
   
        //Check first to see if pp 
        if( _bbs.beam1().A()==1 && _bbs.beam2().A()==1 ){
@@ -359,7 +349,7 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
                double bmin = 0.5;
                double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
                double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
-
                double local_sum=0.0;
 
                // Impact parameter loop 
@@ -402,19 +392,13 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
                double flux_r=local_sum*pi; 
                return flux_r;
 
-               //    bmin = nuclearRadius+nuclearRadius;
-               //    flux_r = nepoint(Egamma,bmin);
-               //    return flux_r;
        }
 
-       //   first call?  - initialize - calculate photon flux
-       Icheck=Icheck+1;
-       if(Icheck > 1) goto L1000f;
+       //   first call or new beam?  - initialize - calculate photon flux
+       Icheck=Icheck+1; 
+       if(Icheck > 1 && beam == Ibeam ) goto L1000f; 
+        Ibeam = beam; 
   
-       rZ=double(_bbs.beam1().Z());
-       rA=double(_bbs.beam1().A());
-       rZ2=double(_bbs.beam2().Z());  //Sergey--dAu
-       rA2=double(_bbs.beam2().A());  //Sergey
   
        //  Nuclear breakup is done by PofB
        //  collect number of integration steps here, in one place
@@ -431,19 +415,18 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
        if (_beamLorentzGamma < 500) 
                Emin=1.E-3;
   
-       //  maximum energy is 6 times the cutoff
-       //      Emax=12.*hbarc*_beamLorentzGamma/RNuc;
-       Emax=6.*hbarc*_beamLorentzGamma/maxradius;
-  
+       //  maximum energy is 12 times the cutoff
+        Emax=_maxPhotonEnergy; 
+       // Emax=12.*hbarc*_beamLorentzGamma/RSum;
        //     >> lnEmin <-> ln(Egamma) for the 0th bin
        //     >> lnEmax <-> ln(Egamma) for the last bin
-  
-       lnEmin=log(Emin);
+       lnEmin=log(Emin);
        lnEmax=log(Emax);
-       dlnE=(lnEmax-lnEmin)/nstep;                                                                                                                  
+       dlnE=(lnEmax-lnEmin)/nstep; 
 
        cout<<" Calculating flux for photon energies from E= "<<Emin 
-           <<" to  "<<Emax<<"  GeV (CM frame) "<<endl;
+           <<" to  "<<Emax<<"  GeV (CM frame) for source nucleus with Z = "<<rZ<<endl;
 
 
        stepmult= exp(log(Emax/Emin)/double(nstep));
@@ -455,7 +438,7 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
                //  integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
                //  use exponential steps
     
-               bmin=0.8*(RNuc+RNuc2); //Start slightly below 2*Radius 
+               bmin=0.8*RSum; //Start slightly below 2*Radius 
                bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
     
                bmult=exp(log(bmax/bmin)/double(nbstep));
@@ -464,7 +447,7 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
     
                if (_bbs.beam2().Z()==1&&_bbs.beam1().A()==2){
                     //This is for deuteron-gold
-                    Xvar = (RNuc+RNuc2)*energy/(hbarc*(_beamLorentzGamma));
+                    Xvar = RSum*energy/(hbarc*(_beamLorentzGamma));
       
                     fluxelement = (2.0/pi)*rZ*rZ*alpha/
                         energy*(Xvar*bessel::dbesk0(Xvar)*bessel::dbesk1(Xvar)-(1/2)*Xvar*Xvar*
@@ -484,8 +467,8 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
                       }
 
                      int nbsteps = 400;
-                     double bmin = 0.7*(RNuc+RNuc2);
-                     double bmax = 2.0*(RNuc+RNuc2) + (8.0*_beamLorentzGamma*hbarc/energy);
+                     double bmin = 0.7*RSum;
+                     double bmax = 2.0*RSum + (8.0*_beamLorentzGamma*hbarc/energy);
                      double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
 
                      double local_sum=0.0;
@@ -572,7 +555,8 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
             
                                                for( int jphi=1;jphi<= nphistep;jphi++){
                                                        phiiter=(double(jphi)-0.5)*deltaphi;
-                                                       //  dist is the distance from the center of the emitting nucleus to the point in question
+                                                       // dist is the distance from the center of the emitting nucleus 
+                                                       // to the point in question
                                                        dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
                                                                   cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
                                                        Xvar=energy*dist/(hbarc*_beamLorentzGamma);  
@@ -593,9 +577,11 @@ photonNucleusCrossSection::photonFlux(const double Egamma)
                                //  multiply by volume element to get total flux in the volume element
                                fluxelement=fluxelement*2.*pi*biter*(biter-bold);
                                //  modulate by the probability of nuclear breakup as f(biter)
+                                // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
                                if (_beamBreakupMode > 1){
                                        fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
                                }
+                                // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
                                integratedflux=integratedflux+fluxelement;
        
                        } //end loop over impact parameter 
@@ -670,9 +656,9 @@ photonNucleusCrossSection::nepoint(const double Egamma,
 double
 photonNucleusCrossSection::sigmagp(const double Wgp)
 {
-       //     >> Function for the gamma-proton --> VectorMeson
-       //     >> cross section. Wgp is the gamma-proton CM energy.
-       //     >> Unit for cross section: fm**2
+       // Function for the gamma-proton --> VectorMeson
+       // cross section. Wgp is the gamma-proton CM energy.
+       // Unit for cross section: fm**2
   
        double sigmagp_r=0.;
   
@@ -730,7 +716,7 @@ photonNucleusCrossSection::sigmagp(const double Wgp)
 
 //______________________________________________________________________________
 double
-photonNucleusCrossSection::sigma_A(const double sig_N)
+photonNucleusCrossSection::sigma_A(const double sig_N, const int beam)
 {                                                         
        // Nuclear Cross Section
        // sig_N,sigma_A in (fm**2) 
@@ -782,13 +768,19 @@ photonNucleusCrossSection::sigma_A(const double sig_N)
                   arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
                 }else if(A_2 == 1 && A_1 != 1){
                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
-                }else{ 
-                  arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
+                }else{
+                 // Check which beam is target 
+                  if( beam == 1 ){ 
+                    arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
+                  }else if( beam==2 ){
+                    arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
+                  }else{
+                    cout<<" Something went wrong here, beam= "<<beam<<endl; 
+                  } 
                 }
     
                Pint=1.0-exp(arg);
                sum=sum+2.*pi*b*Pint*ag[IB];
-
     
                b = 0.5*bmax*(-xg[IB])+0.5*bmax;
 
@@ -797,7 +789,14 @@ photonNucleusCrossSection::sigma_A(const double sig_N)
                 }else if(A_2 == 1 && A_1 != 1){
                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
                 }else{ 
-                  arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
+                 // Check which beam is target 
+                  if( beam == 1 ){ 
+                    arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
+                  }else if(beam==2){
+                    arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
+                  }else{
+                    cout<<" Something went wrong here, beam= "<<beam<<endl; 
+                  } 
                 }
 
                Pint=1.0-exp(arg);
index ccfe5c2..72186c3 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+// $Rev:: 193                         $: revision of last commit
+// $Author:: jnystrand                $: author of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //    Added 18->19 for reading in the luminosity table
 #include "starlightconstants.h"
 #include "inputParameters.h"
 
-
 using namespace std;
 
 
 //______________________________________________________________________________
 readLuminosity::readLuminosity()//:inputread(input)
-: _Warray(0), _Yarray(0), _Farray(0)
+  : _Warray(0), _Yarray(0), _Farray(0), _Farray1(0), _Farray2(0)
 {
   //storing inputparameters into protected variables for the object to use them
   _ReadInputNPT=inputParametersInstance.nmbPtBinsInterference();
@@ -64,6 +63,8 @@ readLuminosity::~readLuminosity()
   if(_Warray) delete [] _Warray;
   if(_Yarray) delete [] _Yarray;
   if(_Farray) delete [] _Farray;
+  if(_Farray1) delete [] _Farray1;
+  if(_Farray2) delete [] _Farray2; 
 }
 
 
@@ -81,40 +82,87 @@ void readLuminosity::read()
       _Farray[i] = new double[_ReadInputnumy];
     }
   }
+  if(!_Farray1) 
+  {
+    _Farray1 = new double*[_ReadInputnumw];
+    for(int i = 0; i < _ReadInputnumw; i++)
+    {
+      _Farray1[i] = new double[_ReadInputnumy];
+    }
+  }
+  if(!_Farray2) 
+  {
+    _Farray2 = new double*[_ReadInputnumw];
+    for(int i = 0; i < _ReadInputnumw; i++)
+    {
+      _Farray2[i] = new double[_ReadInputnumy];
+    }
+  }
+
   double dummy[19]; //14//18
-//  double (*finterm)[starlightLimits::MAXWBINS]=new double[starlightLimits::MAXWBINS][starlightLimits::MAXYBINS];  
-  
+
+  //  double (*finterm)[starlightLimits::MAXWBINS]=new double[starlightLimits::MAXWBINS][starlightLimits::MAXYBINS];  
   //decreased from 1000*1000; too big! causes fault!
+
   double fpart =0.;
   double fptsum=0.;
   ifstream wylumfile;
 
   _f_max=0.0;
+  _f_max1=0.0;
+  _f_max2=0.0;
 
   wylumfile.open("slight.txt");
   for(int i=0;i < 19;i++){ // was 14; this is to account for sergei's additional parameters ie d-Au//was19
     wylumfile >> dummy[i];
   }
+  int A_1 = int(dummy[1]);
+  int A_2 = int(dummy[3]);
+
   for(int i=0;i<_ReadInputnumw;i++){
     wylumfile >> _Warray[i];
   }
   for(int i=0;i<_ReadInputnumy;i++){
     wylumfile >> _Yarray[i];
   }
-  for(int i=0;i<_ReadInputnumw;i++){
-    for(int j=0;j<_ReadInputnumy;j++){
-      wylumfile >> _Farray[i][j];
-      if( _Farray[i][j] > _f_max ) _f_max=_Farray[i][j];
+
+  if( (_ReadInputgg_or_gP == 1) || (A_2 == 1 && A_1 != 1) || (A_1 ==1 && A_2 != 1) ){ 
+    for(int i=0;i<_ReadInputnumw;i++){
+      for(int j=0;j<_ReadInputnumy;j++){
+        wylumfile >> _Farray[i][j];
+        if( _Farray[i][j] > _f_max ) _f_max=_Farray[i][j];
+      }
     }
-  }
-  //Normalize farray (JN 010402)
-  for(int i=0;i<_ReadInputnumw;i++){
-    for(int j=0;j<_ReadInputnumy;j++){
-      _Farray[i][j]=_Farray[i][j]/_f_max;
+    //Normalize farray 
+    for(int i=0;i<_ReadInputnumw;i++){
+      for(int j=0;j<_ReadInputnumy;j++){
+        _Farray[i][j]=_Farray[i][j]/_f_max;
+      }
+    }
+  } else {
+    for(int i=0;i<_ReadInputnumw;i++){
+      for(int j=0;j<_ReadInputnumy;j++){
+        wylumfile >> _Farray1[i][j];
+       //        if( _Farray1[i][j] > _f_max ) _f_max=_Farray1[i][j];
+      }
+    }
+    for(int i=0;i<_ReadInputnumw;i++){
+      for(int j=0;j<_ReadInputnumy;j++){
+        wylumfile >> _Farray2[i][j];
+        if( _Farray1[i][j] + _Farray2[i][j] > _f_max ) _f_max=(_Farray1[i][j] + _Farray2[i][j]);
+      }
+    }
+    //Normalize farray, farray1, farray2 
+    for(int i=0;i<_ReadInputnumw;i++){
+      for(int j=0;j<_ReadInputnumy;j++){
+        _Farray1[i][j]=_Farray1[i][j]/_f_max;
+        _Farray2[i][j]=_Farray2[i][j]/_f_max;
+        _Farray[i][j]=_Farray1[i][j]+_Farray2[i][j];
+      }
     }
   }
 
- if (_ReadInputgg_or_gP != 1 && _ReadInputinterferencemode != 0) {
+  if (_ReadInputgg_or_gP != 1 && _ReadInputinterferencemode != 0) {
         // only numy/2 y bins here, from 0 (not -ymax) to ymax
         double **finterm  = new double*[starlightLimits::MAXWBINS];
         for (int i = 0; i < starlightLimits::MAXWBINS; i++) finterm[i] = new double[starlightLimits::MAXYBINS];
index 328f320..4af4ecc 100644 (file)
@@ -20,9 +20,9 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
-// $Rev:: 191                         $: revision of last commit
+// $Rev:: 193                         $: revision of last commit
 // $Author:: jnystrand                $: author of last commit
-// $Date:: 2014-11-11 11:51:05 +0100 #$: date of last commit
+// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
 //
 // Description:
 //
@@ -69,28 +69,21 @@ wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
        //     This subroutine calculates the cross-section assuming a wide
        //     (Breit-Wigner) resonance.
 
-       // double Av,Wgp,cs,cvma;
        double W,dW,dY;
        double y1,y2,y12,ega1,ega2,ega12;
-       // double t,tmin,tmax;
-       double csgA1,csgA2,csgA12,int_r,dR,rate;
-       double dsigdW,dsigdWalt,dndW,tmp;
+       double csgA1,csgA2,csgA12,int_r,dR;
+       double dsigdW,dsigdWalt,dndW;
        double dsigdW2;
-       // double ax,bx;
        double Eth;
-       int    I,J,NW,NY;
-       // int    K,NGAUSS;
-                                                                                                                                                      
-       // ----------------- !!!!!!!!!!!!!!!!!!!! -----------------------------
-                                                                                                                                                      
-       double bwnorm =bwnormsave;//used to transfer the bwnorm from the luminosity tables
-
-       // --------------------------------------------------------------------
-       //gamma+nucleon threshold.
+       int    I,J,NW,NY,beam;
+
+       double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
 
+        cout<<" in wideResonanceCrossSection, bwnorm: "<<bwnorm<<endl; 
+       //gamma+nucleon threshold.
        Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
                  -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
-                                                                                                                                                      
+                   
        NW   = 100;
        dW   = (_wideWmax-_wideWmin)/double(NW);
   
@@ -105,14 +98,19 @@ wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
        }
   
        cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
-                                                                                                                                                      
+
+        int A_1 = getbbs().beam1().A(); 
+        int A_2 = getbbs().beam2().A();
+
        int_r=0.;
  
+        // Do this first for the case when the first beam is the photon emitter 
+        // Treat pA separately with defined beams 
+        // The variable beam (=1,2) defines which nucleus is the target 
        for(I=0;I<=NW-1;I++){
     
                W = _wideWmin + double(I)*dW + 0.5*dW;
     
-               tmp = 0.0;
                dsigdW=0.0;
                dsigdW2=0.0;
                dsigdWalt=0.0;
@@ -123,41 +121,111 @@ wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
                        y1  = _wideYmin + double(J)*dY;
                        y2  = _wideYmin + double(J+1)*dY;
                        y12 = 0.5*(y1+y2);
+
+                        if( A_2 == 1 && A_1 != 1 ){
+                          // pA, first beam is the nucleus and is in this case the target  
+                          // Egamma = 0.5*W*exp(-Y); 
+                          ega1  = 0.5*W*exp(-y1);
+                          ega2  = 0.5*W*exp(-y2);
+                          ega12 = 0.5*W*exp(-y12);
+                          beam = 1; 
+                        } else if( A_1 ==1 && A_2 != 1){
+                          // pA, second beam is the nucleus and is in this case the target 
+                          // Egamma = 0.5*W*exp(Y); 
+                          ega1  = 0.5*W*exp(y1);
+                          ega2  = 0.5*W*exp(y2);
+                          ega12 = 0.5*W*exp(y12);
+                          beam = 2; 
+                        } else {
+                          // Egamma = 0.5*W*exp(Y);        
+                          ega1  = 0.5*W*exp(y1);
+                          ega2  = 0.5*W*exp(y2);
+                          ega12 = 0.5*W*exp(y12);
+                          beam = 2; 
+                        }
       
-                       ega1  = 0.5*W*exp(y1);
-                       ega2  = 0.5*W*exp(y2);
-                       ega12 = 0.5*W*exp(y12);
+                       //                      ega1  = 0.5*W*exp(y1);
+                       //      ega2  = 0.5*W*exp(y2);
+                       //      ega12 = 0.5*W*exp(y12);
       
                        if(ega1 < Eth) continue;
                        if(ega2 > maxPhotonEnergy()) continue;
           
-                       csgA1=getcsgA(ega1,W);
+                       csgA1=getcsgA(ega1,W,beam);
 
                        //         >> Middle Point                      =====>>>
-                       csgA12=getcsgA(ega12,W);         
+                       csgA12=getcsgA(ega12,W,beam);         
 
                        //         >> Second Point                      =====>>>
-                       csgA2=getcsgA(ega2,W);
+                       csgA2=getcsgA(ega2,W,beam);
       
                        //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
-                       dR  = ega1*photonFlux(ega1)*csgA1;
-                       dR  = dR + 4.*ega12*photonFlux(ega12)*csgA12;
-                       dR  = dR + ega2*photonFlux(ega2)*csgA2;
-                       tmp = tmp+2.*dR*(dY/6.);
+                       dR  = ega1*photonFlux(ega1,beam)*csgA1;
+                       dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
+                       dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
                        dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
       
+                        // cout<<" W: "<<W<<" Y: "<<y12<<" dR: "<<dR<<endl;
                        //For identical beams, we double.  Either may emit photon/scatter
                        //For large differences in Z, we approx, that only beam1 emits photon
                        //and beam2 scatters, eg d-Au where beam1=au and beam2=d
-                       if(getbbs().beam1().A()==getbbs().beam2().A()){
-                               dR  = 2.*dR;
-                       }
+                       //if(getbbs().beam1().A()==getbbs().beam2().A()){
+                       //      dR  = 2.*dR;
+                       //}
                        int_r = int_r+dR;  
                }
        }
-                                                                                                                                                      
-       rate=luminosity()*int_r;
-  
+
+        // Repeat the loop for the case when the second beam is the photon emitter. 
+        // Don't repeat for pA
+        if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){ 
+         for(I=0;I<=NW-1;I++){
+    
+               W = _wideWmin + double(I)*dW + 0.5*dW;
+    
+               dsigdW=0.0;
+               dsigdW2=0.0;
+               dsigdWalt=0.0;
+               dndW=0.0;
+    
+               for(J=0;J<=NY-1;J++){
+      
+                       y1  = _wideYmin + double(J)*dY;
+                       y2  = _wideYmin + double(J+1)*dY;
+                       y12 = 0.5*(y1+y2);
+      
+                        beam = 1; 
+                       ega1  = 0.5*W*exp(-y1);
+                       ega2  = 0.5*W*exp(-y2);
+                       ega12 = 0.5*W*exp(-y12);
+      
+                       if(ega2 < Eth) continue;
+                       if(ega1 > maxPhotonEnergy()) continue;
+          
+                       csgA1=getcsgA(ega1,W,beam);
+
+                       //         >> Middle Point                      =====>>>
+                       csgA12=getcsgA(ega12,W,beam);         
+
+                       //         >> Second Point                      =====>>>
+                       csgA2=getcsgA(ega2,W,beam);
+      
+                       //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
+                       dR  = ega1*photonFlux(ega1,beam)*csgA1;
+                       dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
+                       dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
+                       dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
+      
+                       //For identical beams, we double.  Either may emit photon/scatter
+                       //For large differences in Z, we approx, that only beam1 emits photon
+                       //and beam2 scatters, eg d-Au where beam1=au and beam2=d
+                       // if(getbbs().beam1().A()==getbbs().beam2().A()){
+                       //      dR  = 2.*dR;
+                       // }
+                       int_r = int_r+dR;  
+               }
+         }
+        }
        cout<<" Cross section (mb): "<<10.*int_r<<endl;
-       cout<<" Production rate   : "<<rate<<" Hz"<<endl;
+
 }