TPC/SCDCorrections/sHelix.cxx

0001 #include <iostream>
0002 #include <fstream>
0003 #include "sHelix.h"
0004 #include "TMath.h"
0005
0006 //=================
0007 sHelix::sHelix() :
0008   fR(0),
0009   fW(0),
0010   fC(0),
0011   fPhi(0),
0012   fX0(0),
0013   fY0(0),
0014   fZ0(0),
0015   fDebug(false)
0016 {
0017 }
0018 //=================
0019 sHelix::sHelix(float x0, float y0, float z0, float px, float py, float pz, float q, float b) { // [cm] [GeV] [e] [Tesla]
0020   fX0 = x0;
0021   fY0 = y0;
0022   fZ0 = z0;
0023   float qB = q*b; //[e]x[Tesla]
0024   fW = qB; //[angular frequency]
0025   fPhi = TMath::Pi()+TMath::ATan2(-py,-px); //[rad]
0026   float pt = TMath::Sqrt(px*px+py*py);
0027   float p = TMath::Sqrt(px*px+py*py+pz*pz);;
0028   fR = pt/(0.3*qB)*100; //[cm]
0029   fC = pz/pt*fR*fW; //[pitch]
0030   //std::cout << "sHelix:: p " << p << " | pt " << pt << " | phi " << fPhi*180/TMath::Pi() << std::endl;
0031   fDebug = false;
0032 }
0033 //=================
0034 void sHelix::breakIntoPieces(float t1, float t2, float xx[100][3]) {
0035   // returns coordinates of 100 equidistant pieces of this track
0036   float dt = (t2-t1)/101;
0037   for(int n=0;n!=100; ++n) {
0038     float ti = t1+(n+1)*dt;
0039     xx[n][0] = x(ti);
0040     xx[n][1] = y(ti);
0041     xx[n][2] = z(ti);
0042   }
0043   return;
0044 }
0045 //=================
0046 float sHelix::findFirstInterceptTo(float rd, float hz) {
0047   if(fDebug)
0048     std::cout << "***** findFirstInterceptTo r:" << rd << " z:" << hz << std::endl;
0049   // returns t of first intercept of the helix to a cylinder(r,hz);
0050   float t_1, t_2, t_3, t_4;
0051
0052   // step 1: solving X^2 + Y^2 = r^2
0053   // equation: a cosx + b sinx + c = 0
0054   // a = - 2 R Y0 - 2 R^2 Cos(Phi)
0055   // b = 2 R X0 - 2 R^2 Sin(Phi)
0056   // c = 2 R^2 + X0^2 + Y0^2 - r^2 + 2 R ( Y0 Cos(Phi) - X0 Sin(Phi) )
0057   // x = Phi - W T
0058   // solution: (a^2 + b^2) sinx = -bc +- sqrt( a^2 (a^2+b^2-c^2) )
0059   float a = 2*fR*( fY0 - fR*TMath::Cos(fPhi) );
0060   float b = -2*fR*( fX0 + fR*TMath::Sin(fPhi) );
0061   float c = 2*fR*fR + fX0*fX0 + fY0*fY0 - rd*rd + 2*fR*(fX0*TMath::Sin(fPhi) - fY0*TMath::Cos(fPhi));
0062   //std::cout << "a =" << a << std::endl;
0063   //std::cout << "b =" << b << std::endl;
0064   //std::cout << "c =" << c << std::endl;
0065   if( a*a + b*b < c*c ) { // warranties a^2 + b^2 is more than zero
0066     if(fDebug) {
0067       std::cout << "FAIL: solution is not a real number\n";
0068       std::cout << "a*a =" << a*a << std::endl;
0069       std::cout << "b*b =" << b*b << std::endl;
0070       std::cout << "c*c =" << c*c << std::endl;
0071     }
0072     return 0;
0073   }
0074   float sinx_1 = ( -b*c + TMath::Sqrt(a*a*(a*a+b*b-c*c)) ) / ( a*a + b*b );
0075   float sinx_2 = ( -b*c - TMath::Sqrt(a*a*(a*a+b*b-c*c)) ) / ( a*a + b*b );
0076   if(fDebug) {
0077     std::cout << "Sinx_1 " << sinx_1 << std::endl;
0078     std::cout << "Sinx_2 " << sinx_2 << std::endl;
0079     //std::cout << "r = x_1 / CosPhi " << (-fR*sinx_1+fX0+fR*TMath::Sin(fPhi)) / TMath::Cos(fPhi) << std::endl;
0080     //std::cout << "r = x_2 / CosPhi " << (-fR*sinx_2+fX0+fR*TMath::Sin(fPhi)) / TMath::Cos(fPhi) << std::endl;
0081   }
0082   t_1=t_2=t_3=t_4=99999;
0083   // arc sin returns an angle from -Pi/4 to +Pi/4 so I created two more variables to extend to full 2pi
0084   float asinx_1, asinx_2, asinx_3, asinx_4;
0085   asinx_1 = TMath::ASin(sinx_1);
0086   asinx_2 = TMath::ASin(sinx_2);
0087   if(asinx_1<0) {
0088     asinx_3 = TMath::Pi() - asinx_1;
0089     asinx_1+=TMath::TwoPi();
0090   } else {
0091     asinx_3 = TMath::Pi() - asinx_1;
0092   }
0093   if(asinx_2<0) {
0094     asinx_4 = TMath::Pi() - asinx_2;
0095     asinx_2+=TMath::TwoPi();
0096   } else {
0097     asinx_4 = TMath::Pi() - asinx_2;
0098   }
0099
0100   //phi-wt = asin
0101   if(fDebug) {
0102     std::cout << "PHI0 " << fPhi << std::endl;
0103     std::cout << "ASinx_1 " << asinx_1 << std::endl;
0104     std::cout << "ASinx_2 " << asinx_2 << std::endl;
0105     std::cout << "ASinx_3 " << asinx_3 << std::endl;
0106     std::cout << "ASinx_4 " << asinx_4 << std::endl;
0107   }
0108
0109   t_1 = (fPhi - asinx_1)/fW;
0110   t_2 = (fPhi - asinx_2)/fW;
0111   t_3 = (fPhi - asinx_3)/fW;
0112   t_4 = (fPhi - asinx_4)/fW;
0113
0114   if(fDebug) {
0115     std::cout << "IN R: | x(t) | y(t) || r(t) z(t)" << std::endl;
0116     std::cout << "t_1 " << t_1 << " | " << x(t_1) << " " << y(t_1) << " || " << r(t_1) << " " << z(t_1) << std::endl;
0117     std::cout << "t_2 " << t_2 << " | " << x(t_2) << " " << y(t_2) << " || " << r(t_2) << " " << z(t_2) << std::endl;
0118     std::cout << "t_3 " << t_3 << " | " << x(t_3) << " " << y(t_3) << " || " << r(t_3) << " " << z(t_3) << std::endl;
0119     std::cout << "t_4 " << t_4 << " | " << x(t_4) << " " << y(t_4) << " || " << r(t_4) << " " << z(t_4) << std::endl;
0120   }
0121   if(t_1<0||!TMath::AreEqualAbs(rd,r(t_1),0.1)) t_1 = 99999;
0122   if(t_2<0||!TMath::AreEqualAbs(rd,r(t_2),0.1)) t_2 = 99999;
0123   if(t_3<0||!TMath::AreEqualAbs(rd,r(t_3),0.1)) t_3 = 99999;
0124   if(t_4<0||!TMath::AreEqualAbs(rd,r(t_4),0.1)) t_4 = 99999;
0125
0126   float t1 = TMath::Min( t_1, t_2 );
0127   t1 = TMath::Min( t1, t_3 );
0128   t1 = TMath::Min( t1, t_4 );
0129
0130   // step 2: solving |z| = hz
0131   t_1 = ( hz - fZ0)/fC;
0132   t_2 = (-hz - fZ0)/fC;
0133   if(fDebug) {
0134     std::cout << "IN Z: | x(t) | y(t) || r(t) z(t)" << std::endl;
0135     std::cout << "t_1 " << t_1 << " | " << x(t_1) << " " << y(t_1) << " || " << r(t_1) << " " << z(t_1) << std::endl;
0136     std::cout << "t_2 " << t_2 << " | " << x(t_2) << " " << y(t_2) << " || " << r(t_2) << " " << z(t_2) << std::endl;
0137   }
0138   if(t_1<0||!TMath::AreEqualAbs(hz,z(t_1),0.1)) t_1 = 99999;
0139   if(t_2<0||!TMath::AreEqualAbs(hz,z(t_2),0.1)) t_2 = 99999;
0140   float t2 = TMath::Min( t_1, t_2 );
0141   float tt = TMath::Min(t1,t2);
0142   if(fDebug) std::cout << "  chosen time " << tt << std::endl;
0143   if(tt>4) tt = 99999;
0144   return tt;
0145 }
0146 //=================
0147 void sHelix::SaveTracktoRootScript(float ri, float ro, float hz, char *filec) {
0148   float t1 = findFirstInterceptTo(ri,hz);
0149   float t2 = findFirstInterceptTo(ro,hz);
0150   float xx[100],yy[100],zz[100], xt[100],yt[100],zt[100];
0151   for(int t=0; t!=100; ++t) {
0152     float ti = t/float(100);
0153     xx[t] = x(ti);
0154     yy[t] = y(ti);
0155     zz[t] = z(ti);
0156     float tt = t1 + t/100.0*(t2-t1);
0157     xt[t] = x(tt);
0158     yt[t] = y(tt);
0159     zt[t] = z(tt);
0160   }
0161   std::ofstream out(filec);
0162   out << "{\n";
0163   out << "float x[100] = {";
0164   for(int t=0; t!=100; ++t) out << (t==0?"":",") << xx[t];
0165   out << "};\n";
0166   out << "float y[100] = {";
0167   for(int t=0; t!=100; ++t) out << (t==0?"":",") << yy[t];
0168   out << "};\n";
0169   out << "float z[100] = {";
0170   for(int t=0; t!=100; ++t) out << (t==0?"":",") << zz[t];
0171   out << "};\n";
0172   out << "float xt[100] = {";
0173   for(int t=0; t!=100; ++t) out << (t==0?"":",") << xt[t];
0174   out << "};\n";
0175   out << "float yt[100] = {";
0176   for(int t=0; t!=100; ++t) out << (t==0?"":",") << yt[t];
0177   out << "};\n";
0178   out << "float zt[100] = {";
0179   for(int t=0; t!=100; ++t) out << (t==0?"":",") << zt[t];
0180   out << "};\n";
0181   out << "TPolyLine3D *pol1 = new TPolyLine3D(100,x,y,z);\n";
0182   out << "TPolyLine3D *pol2 = new TPolyLine3D(100,xt,yt,zt);\n";
0183   out << "pol2->SetLineColor(kRed);\n";
0184   out << "pol2->SetLineWidth(3);\n";
0185   out << "pol1->Draw();\n";
0186   out << "pol2->Draw();\n";
0187   out << "}\n";
0188   out.close();
0189   return;
0190 }