[G8b_run] Polarization modification code

[Ken Livingston]

Hi Micheal and all,

I've made a root function out of this.
It's in a macro for testing. Could be pasted into user code.
http://nuclear.gla.ac.uk/~kl/g8b/CbremJlab/mikeTest.C

The function is:
Double_t mikeTest(Double_t pol, Double_t egamma, Double_t cohEdge, Int_t 
Plane, Int_t setting){
   //where pol is from lookup tables
   //egamma = photon energy in MeV
   //cohEdge = current coherent edge pos
   //plane = PARA (0), PERP (1)
   //setting = 1300,1500,1700,1701,1900,2100  (1701 = 1700 auto)

It can be tested at the root prompt like this:
root [0] .L sample_code/mikeTest.C

root [1] mikeTest(1.0,1370.0,1505.0,PARA,1500)
(Double_t)1.00000000000000000e+00

root [2] mikeTest(1.0,1470.0,1505.0,PARA,1500)
(Double_t)9.98414846634865683e-01

root [3] mikeTest(1.0,1470.0,1505.0,PERP,1500)
(Double_t)9.78333263695538746e-01


Seems to give the right kind of answers - Notice there's no correction 
at egamme = coh_edge - 135.0, which is what we'd expect.

I'll try to implement it with the KLambda data as soon as possible.

Cheers,
Ken


On 11/09/2011 11:29 PM, Michael Dugger wrote:
> Hi,
>
> I have some code that I call "polValNew code" that takes
> the polarization from Ken's tables and modifies it to a
> value named polValNew.
>
> Dependencies:
>
> The polValNew code calls a function named polFix.
>
> The polFix function calls a function named polFix0.
>
> The polFix0 function requires a common block vector
> named ldPar(6). Note: The vector is of length 6, but
> only 4 values are used in the code.
>
> The common block vector ldPar has elements that depend
> on the run conditions: Manual/Auto, coherent edge,
> polarization plane.
>
> Even for those that do not use FORTRAN, the code is simple
> and should be easy to read.
>
> All of the required information can be found below:
>
> c-->  polValNew code:
> ###################
>         eDiff = (edgeVec(i) - egVec(j))/1000
>
> c-->  Note: edgeVec(i) is the coherent edge in MeV
> c-->  Note: egVec(j) is the incident photon energy in MeV
>
>         polValNew = polVal
>
> c-->  Note: polVal is the polarization value from Ken's table
> c-->  Note: I get the best results for diffFix = 0.135
> c-->  Note: Use plane = 1 for PARA and plane = 2 for PERP
>
>         if (polFix(eDiff,plane,diffFix).gt.0) then
>            polValNew = polVal/polFix(eDiff,plane,diffFix)
>         endif
>
> ###################
>
> c-->  polFix.f code:
> ###################
>         REAL function polFix(xE,type,diffFix)
>         REAL xE,norm,diffFix
>         INTEGER type
>         polFix = polFix0(xE,type)
>         norm = polFix0(diffFix,type)
>         polFix = polFix/norm
>         END
> ###################
>
> c-->  polFix0.f code:
> ###################
>         REAL function polFix0(xE,type)
>         REAL ldPar(6),fun,xE,yVal,m,b,par(3)
>         REAL par1,par2,par3,tmpFun
>         INTEGER i,j,type,region
>         COMMON /ldBlock/ ldPar
>
>         region = -1
>         if (xE.le.0.025) region = 0
>         if (xE.gt.0.025.and.xE.le.0.050) region = 1
>         if (xE.gt.0.050.and.xE.le.0.075) region = 2
>         if (xE.gt.0.075.and.xE.le.0.100) region = 3
>         if (xE.gt.0.100.and.xE.le.0.125) region = 4
>         if (xE.gt.0.125.and.xE.le.0.150) region = 5
>         if (xE.gt.0.150.and.xE.le.0.175) region = 6
>         if (xE.gt.0.175.and.xE.le.0.200) region = 7
>
>         m = 1.0/0.025
>         b = -0.5
>         yVal = m*xE + b
>
>         if (type.eq.1) then
>            par(1) = ldPar(1)
>            par(2) = ldPar(2)
>         else
>            par(1) = ldPar(4)
>            par(2) = ldPar(5)
>         endif
>
>         par1 = par(1)
>         par2 = par(2)
>
>         tmpFun = 1.0
>         if (region.ge.1) then
>            do i=1,region
>               tmpFun = tmpFun*(par1 + par2*yVal)
>               yVal = yVal - 1.0
>            enddo
>         endif
>         polFix0 = tmpFun
>
>         if (tmpFun.eq.0) polFix0 = 1.0
>
>         END
> ###################
>
> ldPar values:
>
> For 1.3 Manual:
> ldPar(1) = 1.004337
> ldPar(2) = -.4986358E-02
> ldPar(4) = 1.009058
> ldPar(5) = -.5687501E-02
>
> For 1.5 Manual:
> ldPar(1) = 1.012704
> ldPar(2) = -.3850369E-02
> ldPar(4) = 1.002215
> ldPar(5) = -.2256782E-02
>
> For 1.7 Manual:
> ldPar(1) = 0.9888860
> ldPar(2) = -.2256782E-02
> ldPar(4) = 0.9924032
> ldPar(5) = -.3019766E-02
>
> For 2.1 Manual:
> ldPar(1) = 1.004469
> ldPar(2) = -.2188895E-01
> ldPar(4) = 1.067167
> ldPar(5) = -.1876592E-01
>
> For 1.7 Auto:
> ldPar(1) = 0.9844326
> ldPar(2) = 0.7438295E-03
> ldPar(4) = 1.032343
> ldPar(5) = -.1255454E-01
>
> For 1.9 Auto:
> ldPar(1) = 1.008734
> ldPar(2) = -.1557868E-02
> ldPar(4) = 0.9865551
> ldPar(5) = 0.2070320E-02
>
> ###################
>
> Please let me know if there are any questions or if it appears that I have
> left something out.
>
> Take care,
> Michael
>
> _______________________________________________
> G8b_run mailing list
> G8b_run at jlab.org
> https://mailman.jlab.org/mailman/listinfo/g8b_run