Changeset 1357 for trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

Timestamp:

06/12/02 14:34:27 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (modified) (4 diffs)

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

@@ -71 +71 @@
 
 // --------------------------------------------------------------------------
-Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list)
+Bool_t MPairProduction::Process(MParticle *gamma, const Double_t Ep, const Double_t theta, TList *list)
 {
     //
     //  gamma: primary particle from source
     //  phot:  infrared photon from background. (angle = interaction angle)
-    //
+    //  Ep: Energy photon [GeV]
+    //  theta: interaction angle [rad]
+    //
+
+
     const Double_t E0     = 511e-6;              // [GeV]
-
-    //const Double_t theta  = phot->GetAngle();    // [2pi]
-    const Double_t Ep     = phot->GetEnergy();   // [GeV]
     const Double_t Eg     = gamma->GetEnergy();  // [GeV]
 
     const Double_t ctheta = cos(theta);
 
-    const Double_t s      = (Ep*Eg/(2*E0*E0))*(1-ctheta); //[1]
-    if (s<1)
+    const Double_t s      = Ep*Eg*2*(1-ctheta); //[1]
+    const Double_t S      = s/(E0*E0*4); //[1]
+    if (S<1)
         return kFALSE;
 
+    const Double_t stheta = sin(theta);
+
+    const Double_t sqrbetah = s/((Eg+Ep)*(Eg+Ep)) + 1;
+    const Double_t sqrbetae = 1.-1./S;
+
+    const Double_t GammaH = (Eg+Ep)/sqrt(s);
+
+    Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
+
+    fAngle->SetParameter(0, sqrt(sqrbetae));
+
+    const Double_t alpha = psi-acos(fAngle->GetRandom());
+
+    Double_t salpha = sin(alpha);
+    Double_t calpha = cos(alpha);
+
+    const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi)
+
+    Double_t bb = sqrt(sqrbetah/sqrbetae);
+
+    Double_t s1 = calpha/GammaH;
+    Double_t s2 = tphi*s1 - salpha - bb;
+
+    Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
+    Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
+
+    const Double_t E = (Eg+Ep)/2;;
+    const Double_t f = sqrt(sqrbetah*sqrbetae)*salpha;
+
+    cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;
+
+    MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ());
+    MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ());
+    p0 = *gamma; // Set Position and direction
+    p1 = *gamma; // Set Position and direction
+
+    static TRandom rand(0);
+    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
+    p0.SetNewDirection(atan(+tan1), rnd);
+    p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi());
+
+    list->Add(&p0);
+    list->Add(&p1);
+
+    /*
     const Double_t Epg    = Ep/Eg; // 1+Epg ~ 1
     const Double_t Egp    = Eg/Ep; // 1+Egp ~ Egp
@@ -120 +167 @@
     const Double_t E = sqrt(s)*E0/gamma1;
     const Double_t dE = E*Bcosd;
+
+    const Double_t E1 = E0/(E+dE);
+    const Double_t E2 = E0/(E-dE);
+
+    const Double_t beta1 = sqrt(1.-E1*E1);
+    const Double_t beta2 = sqrt(1.-E2*E2);
+
+    const Double_t Bscp = Bsind*cphi;
+    const Double_t spg  = sphi/gamma1;
+    const Double_t cpg  = cphi/gamma1;
+
+    const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp));
+    const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
 
     MElectron &p0 = *new MElectron(E+dE, gamma->GetZ());
@@ -126 +186 @@
     p1 = *gamma; // Set Position and direction
 
-    const Double_t E1 = E0/(E+dE);
-    const Double_t E2 = E0/(E-dE);
-
-    const Double_t beta1 = sqrt(1.-E1*E1);
-    const Double_t beta2 = sqrt(1.-E2*E2);
-
     p0.SetBeta(beta1);
     p1.SetBeta(beta2);
-
-    const Double_t Bscp = Bsind*cphi;
-    const Double_t spg  = sphi/gamma1;
-    const Double_t cpg  = cphi/gamma1;
-
-    const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp));
-    const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
 
     static TRandom rand(0);
@@ -149 +196 @@
     list->Add(&p0);
     list->Add(&p1);
-
+    */
     return kTRUE;
 }