Qweak Qsqrd Error Code 2Trackers

gcc -ansi QsqrdErr.c -lm -o QsqrdErr /* Q^2 = 4*E*Eprim*sin^2(\frac{\theta}{2}) */ #include <stdio.h> #include <stdlib.h> #include <math.h> #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