Skip to content

Commit

Permalink
Create IntegralF2DCollidertprim4A.exe
Browse files Browse the repository at this point in the history
  • Loading branch information
@parkkj
parkkj committed Sep 8, 2014
1 parent e9705f3 commit 24e0c6c
Showing 1 changed file with 207 additions and 0 deletions.
207 changes: 207 additions & 0 deletions Collider/IntegralF2DCollidertprim4A.exe
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,207 @@
#! /usr/bin/perl
use Math::Trig;
use Math::Complex;
# @XBJmin = (0.02,0.04,0.06,0.08);
# @XBJmax = (0.04,0.06,0.08,0.10);
@XBJmin = (0.04,0.06,0.08,0.10);
@XBJmax = (0.06,0.08,0.10,0.12);
#
# This is the script to caluclate the cross-section and F2D from the our MC event generator
# with taking into account ISS (Initial State Smearing) as well as Cross-section angle(50mrad)
#

# wave function with Paris potential
$CC =0.3939;
# wave function with Bonn potential
# $CC =0.3930;


# @XScale = (1E-8,1E-8);
@XScale = (3.2E-5,3.2E-5);

# $Lumi = 1E34; #1/cm2/sec
# $time = 1E7; # approx. 4month
# #$time = 1.5E7; # approx. 6month

# Ch. Weiss used which is similar as (HERA) luminosity (10^-6 nb)
# $Lumi = 1E33; #1/cm2/sec
# $time = 1.2E6; # approx. 2weeks
$Lumi = 1.E33; #1/cm2/sec
$time = 1.9E6;
#2.6E6; # approx. 30days


$cm2nbarn = 1E33;
# @Q2min = (5.0,6.0,7.0,8.0,9.0,11.0,13.0);
# @Q2max = (6.0,7.0,8.0,9.0,11.0,13.0,15.0);

# @Q2min = (18.0,20.0);
@Q2min = (15.0,20.0);
@Q2max = (20.0,22.0);

$DeltaX = 0.02/2;
$DeltaQ2=5.0;
$DeltaAlpha=0.04; # should be range of half
$DeltaTprime=0.005;


$pi=3.14159265358979;
$alpha_R =1.00;


$SED = 2002.4442;


system("rm -f ./f2d_coll_stf_a/x_*.data4a");
system("rm -f ./f2d_coll_stf_a/x_*.f2d");
system("rm -f ./f2d_coll_stf_a/x_*.sf2d");
$DIR = ".";


$filename1a = "f2d_coll_stf_ascii_tprim_crs1.on_smear.dat";
$filename2a = "f2d_coll_stf_ascii_tprim_crs1.on_smear.dat";



$LT = $Lumi*$time/$cm2nbarn; # inverse nbarn


open(itf,"$DIR/$filename2a");
@itemp = <itf>;
close(itf);
$mprot =0.93827;
$mneu = 0.93955;
$Edeu = 0.00222;
$mdeu = 2*$mneu - $Edeu;

$total =0;
for($pp=0;$pp<=$#itemp;$pp++){
@vv = split(/,/,$itemp[$pp]);
$ccount = $vv[3];
$total = $total + $ccount;
}
printf ("total number of event = %11.2f\n",$total);

for($ix=0;$ix<4;$ix++){
$ixxmin = $XBJmin[$ix];
$ixxmax = $XBJmax[$ix];
$ixx = ($ixxmin+$ixxmax)/2;
$ixn = $ix +1;
$scalefactor = $XScale[$ix];

# for($iq=0;$iq<7;$iq++){
for($iq=0;$iq<=1;$iq++){
$iq2min = $Q2min[$iq];
$iq2max = $Q2max[$iq];
$iq2 = ($iq2min+$iq2max)/2;
$iqn = $iq +1;

open(out,">./f2d_coll_stf_a/x\_$ixn\_q2\_$iqn.data4a");
open(out2,">./f2d_coll_stf_a/x\_$ixn\_q2\_$iqn.sf2d");
open(f2dout,">./f2d_coll_stf_a/x\_$ixn\_q2\_$iqn.f2d");

open(ifile1,"$DIR/$filename1a");
@read1 = <ifile1>;
close(ifile1);

open(ifile2,"$DIR/$filename2a");
@read2 = <ifile2>;
close(ifile2);

for($jj=0;$jj<=$#read1;$jj++){
@items1a =split(/,/,$read1[$jj]);
@items2a =split(/,/,$read2[$jj]);
$xb = $items1a[0]/1.;
$q2 = $items1a[1]/1.;
$tprim = $items1a[2]/1.;
$nevnt1a = $items1a[3]/1.;
$nevnt2a = $items1a[4]/1.;
# $Addfactor = 1/$tprim**(1/6);
$Addfactor = 2.0;

if($xb>$ixxmin&&$xb<$ixxmax){

if($q2>$iq2min&&$q2<$iq2max){

$error1a = sqrt($nevnt1a);


if($nevnt1a!=0){

$Gamma = $DeltaX*$DeltaQ2*$DeltaAlpha*$DeltaTprime;

$TPmin0 = -(-2*$Edeu*$mneu + $Edeu**2/2.);
# print "Tprime min = $TPmin0 \n";
# $PR2 = ($tprim-$TPmin0)/2;
$PR2 = ($tprim)/2;
$ER = sqrt($PR2 + $mneu**2);
$wfactor = $pi*$mdeu/4/$ER*$DeltaAlpha*$DeltaTprime;

$xbprime = $xb*$Addfactor;
$ReNormalFactor =0.4;
$JACOBIAN = $xbprime/sqrt($mprot**2-$tprim/2)*$ReNormalFactor;
$CRS_eff = ($nevnt1a/$total)/$Gamma*$JACOBIAN;

$ncount= $CRS_eff*$DeltaX*$DeltaQ2*$LT*$wfactor;
$nerr = sqrt($ncount);

$eCRS_eff = sqrt(2)*$CRS_eff/sqrt($ncount);

# C ...FINE STRUCTURE CONSTANT
$ALEM = (1.)/(137.);
# C ...CONVERSION FACTOR FOR CROSS SECTION UNITS
# C 1/GEV**2 = 0.398 *10**6 NANOBARN
$CONVNB = 0.389*10**6;
#C ...HULTHEN WAVE FUNCTION
$ScaleFactor = 0.85;

$A = 0.045647;
$B = 0.2719;
$C = ($pi**2)*(($A-$B)**2)/$A/$B/($A+$B);
$EPOL = sqrt(-1.*($A**2) + $mneu**2);
#C
#C ...RESIDUE
$RES = 4.0*$EPOL/$C/$alpha_R;
$SPOL = $RES/($tprim)/($tprim);

$XD = $xb/2;
$YY = $q2/$XD/($SED - $mdeu**2);
$CORR = ($YY*$XD*$mdeu)**2/$q2;
$EPS = (1 - $YY - $CORR)/(1 - $YY + ($YY**2)/2 + $CORR);
$FAC = 2*$pi*($ALEM**2)*($YY**2)/($q2**2)/(1 - $EPS);
$FACD = ($FAC/2.)*$CONVNB;


$CorrFirstBin_AlphaCut = 1.0;
$F2D = $CRS_eff*$XD/$FACD*$ScaleFactor*$CorrFirstBin_AlphaCut;
$eF2D = $eCRS_eff*$XD/$FACD;

$F2DSPOL = $F2D/$SPOL;
$eF2DSPOL = $eF2D/$SPOL;

$MCcount = $nevnt2a;
$MCerr = 1/sqrt($MCcount)*100;

print "$ncount $CRS_eff $F2D \n";

# print " $F2D $eF2D $ncount \n";
##### this is only for F2 structure function with nucleon pole allowed

}
printf out ("%6.3e, %6.3e, %8.4e, %16.8e, %16.8e\n",$xb,$q2,$tprim,$F2D,$eF2D);
printf f2dout ("%6.3e, %6.3e, %8.4e, %16.8e, %16.8e\n",$xb,$q2,$tprim,$F2DSPOL,$eF2DSPOL);
if($tprim<0.05){
printf out2 ("%6.3e, %6.3e, %8.4e, %16.8e, %16.8e\n",$xb,$q2,$tprim,$F2DSPOL,$eF2DSPOL);
}

}

}
}
close(out);
close(cout);
close(nout);
close(MCout);
close(out2);
}
}

0 comments on commit 24e0c6c

Please sign in to comment.