ParticleID/RICHAnalysis/dualrich_analyzer.cc

0001 /////////Code developer Alessio Del Dotto/////////
0002 ////////Compilation and usage in ROOT/////////
0003 //.L dualrich_analyzer.cpp++
0004 //eic_dual_rich *a = new eic_dual_rich()
0005 //a->some method of the class
0006
0007 #include "dualrich_analyzer.h"
0008
0009 using namespace std;
0010
0011
0012 void eic_dual_rich::set_mirror(double center_posx, double center_posy, double center_posz, double radius){
0013
0014   cx = center_posx;
0015   cy = center_posy;
0016   cz = center_posz;
0017   R_mirror = radius;
0018
0019 }
0020
0021 void eic_dual_rich::set_radiator_one(double mean_refraction_index1){
0022
0023   refidx1 = mean_refraction_index1;
0024
0025 }
0026
0027 void eic_dual_rich::set_radiator_two(double mean_refraction_index2){
0028
0029   refidx1 = mean_refraction_index2;
0030
0031 }
0032
0033 double eic_dual_rich::ind_ray(double Ex, double Ey, double Ez, double Dx, double Dy, double Dz, double vx, double vy, double vz, int select_radiator){
0034
0035   //int sel_radiator = 2;
0036
0037   //if(select_radiator == 1) sel_radiator = 1;
0038
0039   double cex,cey,cez;
0040   double cdx,cdy,cdz;
0041
0042   int i;//,iwhere;
0043
0044   double th,a,d;
0045   double x,dx;
0046
0047   double y,y1;
0048
0049   double eps = 0.00000000001;
0050   double R = 0;
0051
0052   R = R_mirror;
0053
0054   double esx, esy, esz, es;
0055   double ref_frac, theta1, theta2;
0056
0057   double theta_c = 0.;
0058
0059   cex = -cx+Ex;
0060   cey = -cy+Ey;
0061   cez = -cz+Ez;
0062
0063   cdx = -cx+Dx;
0064   cdy = -cy+Dy;
0065   cdz = -cz+Dz;
0066
0067   //a = TMath::Sqrt(cex*cex+cey*cey+cez*cez);
0068   //d = TMath::Sqrt(cdx*cdx+cdy*cdy+cdz*cdz);
0069   //th = TMath::ACos((cdx*cex+cdy*cey+cdz*cez)/a/d);
0070
0071   a = sqrt(cex*cex+cey*cey+cez*cez);
0072   d = sqrt(cdx*cdx+cdy*cdy+cdz*cdz);
0073   th = acos((cdx*cex+cdy*cey+cdz*cez)/a/d);
0074
0075   i = 0;
0076   x = th/2.;
0077   y = R*(a*sin(x)-d*sin(th-x))+a*d*sin(th-2*x);
0078   y1 = R*(a*cos(x)+d*cos(th-x))-2*a*d*cos(th-2*x);
0079   dx = -1*y/y1;
0080
0081   while(abs(dx)>eps && i<100){
0082
0083     x+=dx;
0084     y = R*(a*sin(x)-d*sin(th-x))+a*d*sin(th-2*x);
0085     y1 = R*(a*cos(x)+d*cos(th-x))-2*a*d*cos(th-2*x);
0086     dx = -1*y/y1;
0087     i++;
0088
0089   }
0090
0091   if(i>=100) cout<<"Not convergent"<<endl;
0092
0093   if(i<100){
0094     sx = cx + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cex + (R*sin(x)/d/sin(th))*cdx;
0095     sy = cy + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cey + (R*sin(x)/d/sin(th))*cdy;
0096     sz = cz + (R*cos(x)/a-R*sin(x)/tan(th)/a)*cez + (R*sin(x)/d/sin(th))*cdz;
0097   }
0098
0099   esx = sx - Ex;
0100   esy = sy - Ey;
0101   esz = sz - Ez;
0102
0103   es = sqrt(esx*esx+esy*esy+esz*esz);
0104
0105   esx = esx/es;
0106   esy = esy/es;
0107   esz = esz/es;
0108
0109   if(select_radiator == 1){
0110     ref_frac = refidx2/refidx1;
0111
0112     theta2 = acos(esz);
0113     theta1 = asin(sin(theta2)*ref_frac);
0114
0115     esx = esx*sin(theta1)/sin(theta2);
0116     esy = esy*sin(theta1)/sin(theta2);
0117     esz = cos(theta1);
0118   }
0119   theta_c = acos((esx*vx+esy*vy+esz*vz));
0120
0121   return theta_c;
0122
0123 }
0124
0125 void eic_dual_rich::fill_cherenkov_array(double angle){
0126
0127   ch_vector.push_back (angle);
0128
0129 }
0130
0131 void eic_dual_rich::cut_cherenkov_array(double theta_min, double theta_max){
0132
0133   vector<double> cut_vector;
0134
0135   for (unsigned int i=0; i<ch_vector.size(); i++){
0136
0137     if(theta_min<theta_max && ch_vector.at(i)>=theta_min && ch_vector.at(i)<=theta_max) cut_vector.push_back (ch_vector.at(i));
0138     else if (theta_min>=theta_max) cout<<"Applied cut wrong: theta_min>theta_max"<<endl;
0139
0140   }
0141
0142   ch_vector.clear();
0143
0144   for (unsigned int i=0; i<cut_vector.size(); i++){
0145
0146     ch_vector.push_back(cut_vector.at(i));
0147
0148   }
0149
0150 }
0151
0152 double eic_dual_rich::mean_cherenkov_angle(){
0153
0154   double sum = 0.;
0155
0156   for (unsigned int i=0; i<ch_vector.size(); i++){
0157
0158     sum += ch_vector.at(i);
0159
0160   }
0161
0162   sum = sum/ch_vector.size();
0163
0164   return sum;
0165
0166 }
0167
0168 double eic_dual_rich::SD_cherenkov_angle(){
0169
0170   double sum = 0.;
0171
0172   for (unsigned int i=0; i<ch_vector.size(); i++){
0173
0174     sum += (ch_vector.at(i)-mean_cherenkov_angle())*(ch_vector.at(i)-mean_cherenkov_angle());
0175
0176   }
0177
0178   sum = sum/ch_vector.size();
0179
0180   return sqrt(sum);
0181
0182 }
0183
0184 void eic_dual_rich::clear_cherenkov_array(){
0185
0186   ch_vector.clear();
0187
0188 }
0189
0190
0191
0192