Qweak Qsqrd Error Code 2Trackers

 gcc -ansi QsqrdErr.c -lm -o QsqrdErr

/*

Q^2 = 4*E*Eprim*sin^2(\frac{\theta}{2})

• /

1. include <stdio.h>
2. include <stdlib.h>
3. include <math.h>
4. define data_prec float

data_prec QsqrdErr(data_prec Qsqrd, data_prec Ebeam, data_prec Eprime, data_prec Theta, data_prec EbeamErr, data_prec E primeErr, data_prec ThetaErr); main(int argc, char *argv[]) {

 int i,j;
 data_prec Qsqrd,Ebeam,Eprime,Theta,Zdc,y1,y2,ThetaRad;
 data_prec EbeamErr,EprimeErr,ThetaErr,ZdcErr;
 data_prec EbeamRelativeErr,QsqrdRelativeErr,ZdcRelativeErr;
 data_prec TrackerErr,TrackerErrSqrd;
 data_prec Mnucleon=938.27231;
 data_prec PI=3.1415927;

 if(argc<7)
   {
     fprintf(stderr,"usage:Ebeam(MeV) Theta(Degrees) y1(cm) y2(cm) EbeamErr(MeV) ZdcErr(cm)\n");
     exit(0);
   }
 else
   {
     Ebeam=atof(*(argv+1));
     ++argv;
     Theta=atof(*(argv+1));
     ++argv;
     y1=atof(*(argv+1));
     ++argv;
     y2=atof(*(argv+1));
     ++argv;
     EbeamErr=atof(*(argv+1));
     ++argv;
     ZdcErr=atof(*(argv+1));
     ++argv;
   }

 ThetaRad=Theta*PI/180;

 EbeamRelativeErr= 4*(1+Ebeam*sin(ThetaRad/2)*sin(ThetaRad/2)/Mnucleon)*(1+Ebeam*sin(ThetaRad/2)*sin(ThetaRad/2)/Mnucl

eon)*EbeamErr*EbeamErr/Ebeam/Ebeam;

 Zdc=40.0;
 fprintf(stderr,"TrackerErr (um),Zdc (m),QsqrdRelativeErr(\%)\n");
 for(i=0;i<500;i++)
   {
     Zdc=Zdc+ 1.0;
     TrackerErr=0;
     for(j=0;j<10;j++)
       {
         TrackerErr=TrackerErr+100;
         TrackerErrSqrd=2*TrackerErr*TrackerErr/10000/10000; /* in units of cm*/
         ZdcRelativeErr=(y2-y1)*(y2-y1)*ZdcErr*ZdcErr/Zdc/Zdc;

         QsqrdRelativeErr=Zdc*Zdc/(Zdc*Zdc+y2*y2-2*y2*y1+y1*y1)/(Zdc*Zdc+y2*y2-2*y2*y1+y1*y1)/tan(ThetaRad/2)/tan(Thet

aRad/2);

         if(Zdc==50 && j==1)
             fprintf(stderr,"QsqrdRelativeErr=%g\n",Zdc*Zdc/(Zdc*Zdc+y2*y2-2*y2*y1+y1*y1)/(Zdc*Zdc+y2*y2-2*y2*y1+y1*y1

)/tan(ThetaRad/2)/tan(ThetaRad/2));

         QsqrdRelativeErr=QsqrdRelativeErr*(TrackerErrSqrd+ZdcRelativeErr);
         if(Zdc==50 && j==1)
           fprintf(stderr,"QsqrdRelativeErr=%g\n",QsqrdRelativeErr);
         QsqrdRelativeErr=QsqrdRelativeErr+EbeamRelativeErr;

         QsqrdRelativeErr=QsqrdRelativeErr/(1+2*Ebeam*sin(ThetaRad/2)*sin(ThetaRad/2)/Mnucleon)*(1+2*Ebeam*sin(ThetaRa

d/2)*sin(ThetaRad/2)/Mnucleon);

         if(Zdc==50 && j==1)
           {
             fprintf(stderr,"QsqrdRelativeErr/%g=%g\n",(1+2*Ebeam*sin(ThetaRad/2)*sin(ThetaRad/2)/Mnucleon)*(1+2*Ebeam

• sin(ThetaRad/2)*sin(ThetaRad/2)/Mnucleon),QsqrdRelativeErr);

             fprintf(stderr,"ZdcRelativeErr=%g\n",ZdcRelativeErr);
             fprintf(stderr,"EbeamRelativeErr=%g\n",EbeamRelativeErr);
             fprintf(stderr,"%g,%g,%g\n",TrackerErr,Zdc/100,sqrt(QsqrdRelativeErr)*100);
           }
         if(QsqrdRelativeErr < 0.00012)
                fprintf(stdout,"%g,%g,%g\n",TrackerErr/1000,Zdc/100,sqrt(QsqrdRelativeErr)*100);
       }
   }
 

}

data_prec QsqrdErr(data_prec Qsqrd, data_prec Ebeam, data_prec Eprime, data_prec Theta, data_prec EbeamErr, data_prec E primeErr, data_prec ThetaErr) {

 int phase;
 data_prec Pi,Angle,DeltaX,Err;
 data_prec TaylorCosine,TaylorSine,Number;
 data_prec Cbeam,Ceprime,Ctheta;

 Pi = 3.1415926535898; 
 Angle=Theta*Pi/180;  /* The angle in units of radians*/

 TaylorCosine=1.0;
 phase=(-1);
 Number=2;
 DeltaX=Angle*Angle;
 /*
 while(DeltaX>1E-29 && Number<1E29)
   {
     TaylorCosine=TaylorCosine+phase*DeltaX/Number;
     phase=(-1)*phase;
     Number=Number*(Number+1);
     DeltaX=DeltaX*Angle*Angle;
     fprintf(stderr,"TaylorCosine=%g+%g\n",TaylorCosine,phase*Angle/Number);
   }
 fprintf(stderr,"Cos(%g)=%g\n",Theta,cos(Angle));
 fprintf(stderr,"Sin(%g)=%g\n",Theta,sin(Angle));
 */

 Cbeam=EbeamErr/Ebeam;
 Ceprime=EprimeErr/Eprime;
 Ctheta=cos(Angle/2)*ThetaErr*Pi/sin(Angle/2)/180;

 fprintf(stderr,"Relative Errors\n\t E\tEprime\tTheta\n");
 fprintf(stderr,"\t%g\t%g\t%g\n",Cbeam,Ceprime,Ctheta);
 fprintf(stderr,"contributions (percent)\n");

 Err=sqrt(Cbeam*Cbeam+Ceprime*Ceprime+Ctheta*Ctheta);
 fprintf(stderr,"\t%g\t%g\t%g\n",Cbeam/Err,Ceprime/Err,Ctheta/Err);
 fprintf(stderr,"Err*Qsqrd=%g * %g=%g\n",Err,Qsqrd,Err*Qsqrd);
 return(Err);

}

Go Back Qweak_Qsqrd_Tracking