PID/hpDIRC/phi_code_dirc.C

0001 // This code calculates the input phi angle needed to shoot the charged particle such that it hits the middle of a given DIRC bar in the magnetic field
0002 // This prints the phi values for theta range 30-150 deg in 5 deg steps
0003 // Author: Nilanga Wickramaarachchi
0004 // Based on code by William Llope
0005
0006 #include "TMath.h"
0007
0008 void phi_code_dirc()
0009 {
0010   int kTargetBar = -1;
0011   double fTargetBarY =  0;
0012   double fBarsGap =  0.15; // gap between bars in y-direction (mm)
0013   double momentum = 6.0; // momentum of charged particle (GeV/c)
0014   double fRadiatorW = 34.865; // width of bar (mm)
0015
0016   double theta_value[25];
0017   double phi_value[25];
0018
0019   for(int i=0; i < 25; i++)
0020       {
0021     theta_value[i] = 30 + 5*i;
0022         double theta = theta_value[i];
0023
0024     kTargetBar = 4; // use the bar number from 0-9
0025     fTargetBarY = ((double)(kTargetBar-5))*(fRadiatorW + fBarsGap) + fRadiatorW/2. + fBarsGap/2.;
0026
0027     //std::cout << "kTargetBar = " << kTargetBar << std::endl;
0028     //std::cout << "fTargetBarY = "<< fTargetBarY <<std::endl;
0029
0030     //double partheta = (180. - fRun->getTheta());   // deg
0031     double partheta = 180. - theta; // deg
0032     double distheta = fabs(partheta-90.);// deg
0033     //std::cout << "distheta = " << distheta << std::endl;
0034     double Bval = (-1.75 + (1.75-1.6)*(distheta/60.)); // MARCO 1.7T, theta scalings to handle radial component of field
0035     //std::cout << "Bval = " << Bval << std::endl;
0036
0037     double pt = momentum*TMath::Sin(theta*TMath::DegToRad());
0038     //std::cout << "pt = " << pt << std::endl;
0039
0040     double R = 1000.*pt/0.3/fabs(Bval);
0041     double x1 = 0;// primary vertex
0042     double y1 = 0;// primary vertex
0043     double x2 = 708.5;
0044     //double x2 = 770.5; // dirc radius (mm)
0045     double y2 = fTargetBarY;// target point !!!! HERE ASSUMES TEN BARS PER BOX
0046     double xa = (x2 - x1)/2.;
0047     double ya = (y2 - y1)/2.;
0048     double x0 =  x1 + xa;
0049     double y0 =  y1 + ya;
0050     double a  = sqrt(xa*xa + ya*ya);
0051     double b  = sqrt(R*R - a*a);
0052     double x3 = x0 + b*ya/a;
0053     double y3 = y0 - b*xa/a;
0054     double x4 = x0 - b*ya/a;
0055     double y4 = y0 + b*xa/a;
0056     double ang3 = atan2(x3,y3);
0057     // reverse order to get phi angle of perpendicular to ray from PV to (x3,y3), NEG
0058     double ang4 = atan2(x4,y4);
0059     // reverse order to get phi angle of perpendicular to ray from PV to (x4,y4), POS
0060     //if (acharge>0.){ phi = -ang4; } else // pos
0061     //if (acharge<0.){ phi = -ang3; } // neg
0062
0063     phi_value[i] = 360 - ang4*TMath::RadToDeg();
0064
0065     std::cout << phi_value[i] << "\t";
0066     if(i==24) std::cout << "\n" << std::endl;
0067
0068     //phi = ang3;
0069     //std::cout<<"ang3 deg = "<< -ang3*TMath::RadToDeg() <<std::endl;
0070       }
0071
0072   for(int i=0; i < 25; i++)
0073     {
0074       std::cout << "theta = " << theta_value[i] << "\t" << "phi = " << phi_value[i] << std::endl;
0075     }
0076
0077 }
0078
0079
0080