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