TpcLheTrackCuts.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463
  1. #include "TpcLheTrackCuts.h"
  2. ClassImp(TpcLheTrackCuts)
  3. TpcLheTrackCuts* TpcLheTrackCuts::fInstance = 0;
  4. //________________________________________________________________
  5. TpcLheTrackCuts* TpcLheTrackCuts::Instance() {
  6. //--- Returns instance.
  7. if (fInstance == 0) {
  8. fInstance = new TpcLheTrackCuts();
  9. }
  10. return fInstance;
  11. }
  12. //________________________________________________________________
  13. TpcLheTrackCuts::TpcLheTrackCuts() {
  14. //---
  15. }
  16. //________________________________________________________________
  17. TpcLheTrackCuts::~TpcLheTrackCuts() {
  18. //---
  19. if (fInstance)
  20. delete fInstance;
  21. }
  22. //________________________________________________________________
  23. void TpcLheTrackCuts::Reset() {
  24. }
  25. //________________________________________________________________
  26. Double_t const TpcLheTrackCuts::DistFromLine(Double_t Xh,
  27. Double_t Yh,
  28. Double_t *coeff) {
  29. // Returns the distance of a point to a straight line.
  30. // straight line is given by its to coefficients: y = coeff[0]*x + coeff[1].
  31. Double_t x = (coeff[0] / (1 + coeff[0]*coeff[0])) *
  32. (1/coeff[0] * Xh + Yh - coeff[1]);
  33. return TMath::Sqrt(TMath::Power(x - Xh, 2) +
  34. TMath::Power(coeff[0]*x + coeff[1] - Yh, 2));
  35. }
  36. //________________________________________________________________
  37. void TpcLheTrackCuts::SetHighMomTrackCuts() {
  38. //--- Sets cuts of tracks
  39. // Range for prediction
  40. fTrackCuts[0].fDelX = .001; //
  41. fTrackCuts[1].fDelX = 3*.0055;
  42. // 2 -- seed
  43. fTrackCuts[0].fMaxHitDist = 3.;
  44. fTrackCuts[1].fMaxHitDist = 5.;
  45. }
  46. //________________________________________________________________
  47. void TpcLheTrackCuts::SetLowMomTrackCuts() {
  48. // Sets cuts of tracks for the given vertex constraint.
  49. // Range for prediction
  50. fTrackCuts[0].fDelX = 3*.025; //
  51. fTrackCuts[1].fDelX = 3*.034;
  52. }
  53. //________________________________________________________________
  54. Bool_t TpcLheTrackCuts::VerifyTrack(TpcLheCMTrack *track,
  55. TpcLheCMPoint *hit, Bool_t back) {
  56. // --------------------------------
  57. if (track->GetNumberOfPoints() < 3 ) {
  58. // tracklet
  59. return kTRUE;
  60. }
  61. else {
  62. // track
  63. Double_t coeff[3];
  64. coeff[0] = coeff[1] = coeff[2] = 0;
  65. CMLineFit(track, coeff);
  66. Double_t dist_cm = DistFromLine(hit->GetXprime(),
  67. hit->GetYprime(), coeff);
  68. if(dist_cm < fTrackCuts[0].fDelX) {
  69. return kTRUE;
  70. }
  71. else
  72. return kFALSE;
  73. }
  74. } ///:~
  75. //________________________________________________________________
  76. Double_t TpcLheTrackCuts::GetChi2Bend(Int_t pl) {
  77. //---
  78. return fTrackCuts[pl].fBendChi2;
  79. }
  80. //________________________________________________________________
  81. Double_t TpcLheTrackCuts::GetChi2Deep(Int_t pl) {
  82. //---
  83. return fTrackCuts[pl].fDeepChi2;
  84. }
  85. //________________________________________________________________
  86. Double_t TpcLheTrackCuts::GetDelX(Int_t pl) {
  87. //---
  88. return fTrackCuts[pl].fDelX;
  89. }
  90. //________________________________________________________________
  91. void TpcLheTrackCuts::Circle3pnts(Double_t x[],Double_t y[], Double_t r[]) {
  92. // calc center and R of circle from 3 points
  93. Double_t m_a = (y[1] - y[0]) / (x[1] - x[0]);
  94. Double_t m_b = (y[2] - y[1]) / (x[2] - x[1]);
  95. if (m_a != m_b) {
  96. Double_t x_c = .5*(m_a* m_b*(y[0] - y[2]) +
  97. m_b*(x[0] + x[1]) -
  98. m_a*(x[1] + x[2])) / (m_b - m_a);
  99. Double_t y_c = - (x_c - .5*(x[0] + x[1])) / m_a +
  100. .5*(y[0] + y[1]);
  101. Double_t dely = y_c - y[0];
  102. Double_t delx = x_c - x[0];
  103. Double_t rad =TMath::Sqrt(delx*delx + dely*dely);
  104. r[0] = x_c;
  105. r[1] = y_c;
  106. r[2] = rad;
  107. }
  108. }
  109. //________________________________________________________________
  110. int TpcLheTrackCuts::
  111. Circle_Circle_Intersection(double x0, double y0, double r0,
  112. double x1, double y1, double r1,
  113. double *xi, double *yi,
  114. double *xi_prime, double *yi_prime) {
  115. //---
  116. double a, dx, dy, d, h, rx, ry;
  117. double x2, y2;
  118. /* dx and dy are the vertical and horizontal distances between
  119. * the circle centers.
  120. */
  121. dx = x1 - x0;
  122. dy = y1 - y0;
  123. /* Determine the straight-line distance between the centers. */
  124. d = sqrt((dy*dy) + (dx*dx));
  125. /* Check for solvability. */
  126. if (d > (r0 + r1)) {
  127. /* no solution. circles do not intersect. */
  128. return 0;
  129. }
  130. if (d < TMath::Abs(r0 - r1)) {
  131. /* no solution. one circle is contained in the other */
  132. return 0;
  133. }
  134. /* 'point 2' is the point where the line through the circle
  135. * intersection points crosses the line between the circle
  136. * centers.
  137. */
  138. /* Determine the distance from point 0 to point 2. */
  139. a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
  140. /* Determine the coordinates of point 2. */
  141. x2 = x0 + (dx * a/d);
  142. y2 = y0 + (dy * a/d);
  143. /* Determine the distance from point 2 to either of the intersection
  144. * points.
  145. */
  146. h = sqrt((r0*r0) - (a*a));
  147. /* Now determine the offsets of the intersection points from point
  148. * 2.
  149. */
  150. rx = -dy * (h/d);
  151. ry = dx * (h/d);
  152. /* Determine the absolute intersection points. */
  153. *xi = x2 + rx;
  154. *xi_prime = x2 - rx;
  155. *yi = y2 + ry;
  156. *yi_prime = y2 - ry;
  157. return 1;
  158. }
  159. //________________________________________________________________
  160. Double_t TpcLheTrackCuts::
  161. GetPhiPrediction(TpcLheCMTrack *track) {
  162. //---
  163. Int_t last = (track->GetRHits())->GetLast();
  164. TpcLheHit *hit0 = (TpcLheHit *)(track->GetRHits())->At(last);
  165. TpcLheHit *hit1 = (TpcLheHit *)(track->GetRHits())->At(last-1);
  166. Double_t x[3], y[3], rez[3];
  167. x[0] = 0.; // hit0->GetX();
  168. x[1] =hit0->GetX();
  169. x[2] =hit1->GetX();
  170. y[0] =0.; // hit0->GetY();
  171. y[1] =hit0->GetY();
  172. y[2] =hit1->GetY();
  173. Circle3pnts(x, y, rez);
  174. Double_t x_c = rez[0];
  175. Double_t y_c = rez[1];
  176. Double_t Rad = rez[2];
  177. track->SetRadius(Rad);
  178. track->SetVertex(x_c, y_c, 0.);
  179. TVector3 v(x_c, y_c, 0.);
  180. Double_t phi_c = (v.Phi() >= 0.) ? phi_c = v.Phi() :
  181. phi_c = v.Phi() + 2.*TMath::Pi();
  182. return phi_c;
  183. }
  184. //________________________________________________________________
  185. Double_t TpcLheTrackCuts::
  186. GetThetPrediction(TpcLheCMTrack *track, Double_t zst, Bool_t back) {
  187. //--- return alpha angle in the next plane
  188. TpcLheHit *hit1, *hit0;
  189. if (back) {
  190. hit1 = (TpcLheHit *)(track->GetRHits())->At(1);
  191. hit0 = (TpcLheHit *)(track->GetRHits())->At(0);
  192. }
  193. else {
  194. Int_t last = (track->GetRHits())->GetLast();
  195. hit0 = (TpcLheHit *)(track->GetRHits())->At(last);
  196. hit1 = (TpcLheHit *)(track->GetRHits())->At(last-1);
  197. }
  198. Double_t a0 = (hit0->GetY() - hit1->GetY()) / (hit0->GetZ() - hit1->GetZ());
  199. Double_t a1 = (hit0->GetY()*hit1->GetZ() - hit1->GetY()*hit0->GetZ()) /
  200. (hit1->GetZ() - hit0->GetZ());
  201. Double_t y_cross = a0 * zst + a1;
  202. // return TMath::ATan2(y_cross, zst);
  203. return y_cross;
  204. }
  205. //________________________________________________________________
  206. void TpcLheTrackCuts::CMLineFit(TpcLheCMTrack *track,
  207. Double_t *a) {
  208. //---
  209. TObjArray *trackpoints = track->GetCMHits();
  210. Int_t n = trackpoints->GetEntriesFast();
  211. TpcLheCMPoint *trackpoint = NULL;
  212. TArrayD *x = new TArrayD(); x->Set(n);
  213. TArrayD *y = new TArrayD(); y->Set(n);
  214. TArrayD *delx = new TArrayD(); delx->Set(n);
  215. TArrayD *dely = new TArrayD(); dely->Set(n);
  216. Int_t ip = 0;
  217. for (Int_t is = 0; is < n; is++) {
  218. trackpoint = (TpcLheCMPoint *)trackpoints->At(is);
  219. // trackpoint->Print();
  220. if (trackpoint->GetXprime() != 0.) {
  221. x->AddAt(trackpoint->GetXprime(), ip);
  222. y->AddAt(trackpoint->GetYprime(), ip);
  223. delx->AddAt(trackpoint->GetXprimeerr(), ip);
  224. dely->AddAt(trackpoint->GetYprimeerr(), ip);
  225. ip++;
  226. }
  227. }
  228. Double_t coeff[3];
  229. LineFit(x, delx, y, dely, coeff, ip);
  230. a[0]= coeff[0];
  231. a[1]= coeff[1];
  232. a[2]= coeff[2];
  233. if (coeff[1] != 0.) {
  234. Double_t D = TMath::Sqrt((coeff[0]*coeff[0] + 1.) /
  235. (4.* coeff[1]*coeff[1]));
  236. track->SetRadius(D);
  237. }
  238. else {
  239. track->SetRadius(0.);
  240. }
  241. delete x;
  242. delete y;
  243. delete delx;
  244. delete dely;
  245. }
  246. //________________________________________________________________
  247. void TpcLheTrackCuts::LineFit(TArrayD *x, TArrayD *delx,
  248. TArrayD *y, TArrayD *dely,
  249. Double_t *coeff, Int_t np) {
  250. // fit points to a straight line y = ax + b
  251. // NRC p. 661
  252. Double_t L11 = 0.;
  253. Double_t L12 = 0.;
  254. Double_t L22 = 0.;
  255. Double_t g1 = 0.;
  256. Double_t g2 = 0.;
  257. for (Int_t i = 0; i < np; i++) {
  258. L11 += 1.; // S
  259. L12 += x->At(i); // Sx
  260. L22 += x->At(i) * x->At(i); // Sxx
  261. g1 += y->At(i); // Sy
  262. g2 += x->At(i) * y->At(i); // Sxy
  263. }
  264. Double_t D = L11*L22 - L12*L12; // S*Sxx - (Sx)^2
  265. coeff[0] = (g2*L11 - g1*L12)/D; // a = (Sxy*S - Sy*Sx)/D
  266. coeff[1] = (g1*L22 - g2*L12)/D; // b = (Sy*Sxx - Sxy*Sx)/D
  267. // Calculate chi squared
  268. Double_t chi2 = 0.;
  269. if (np > 2) {
  270. for ( Int_t i = 0; i < np; i++) {
  271. chi2 += TMath::Power(y->At(i) - coeff[0] * x->At(i) - coeff[1], 2.) /
  272. TMath::Abs(coeff[0] * x->At(i) + coeff[1]);
  273. }
  274. chi2 = chi2/float(np - 2);
  275. }
  276. coeff[2] = chi2;
  277. }
  278. //________________________________________________________________
  279. void TpcLheTrackCuts::DeepAngleFit(const TpcLheCMTrack *track,
  280. Double_t *a) {
  281. //---
  282. TObjArray *trackpoints = track->GetCMHits();
  283. Int_t n = trackpoints->GetEntriesFast();
  284. if (n > 2) {
  285. TArrayD *zv = new TArrayD(); zv->Set(n);
  286. TArrayD *yv = new TArrayD(); yv->Set(n);
  287. TArrayD *delzv = new TArrayD(); delzv->Set(n);
  288. TArrayD *delyv = new TArrayD(); delyv->Set(n);
  289. for (Int_t is = 0; is < n; is++) {
  290. TpcLheCMPoint *trackpoint = (TpcLheCMPoint *)trackpoints->At(is);
  291. // trackpoint->Print();
  292. zv->AddAt(trackpoint->GetZv(), is);
  293. yv->AddAt(trackpoint->GetYv(), is);
  294. delzv->AddAt(trackpoint->GetZverr(), is);
  295. delyv->AddAt(trackpoint->GetYverr(), is);
  296. }
  297. Double_t coeff[3];
  298. LineFit(zv, delzv, yv, delyv, coeff, n);
  299. a[0]= coeff[0];
  300. a[1]= coeff[1];
  301. a[2]= coeff[2];
  302. delete zv;
  303. delete yv;
  304. delete delzv;
  305. delete delyv;
  306. }
  307. else {
  308. TpcLheCMPoint *hit1 = (TpcLheCMPoint *)trackpoints->At(0);
  309. TpcLheCMPoint *hit2 = (TpcLheCMPoint *)trackpoints->At(1);
  310. a[0] = (hit1->GetY() - hit2->GetY()) / (hit1->GetZ() - hit2->GetZ());
  311. a[1] = (hit1->GetY()*hit2->GetZ() - hit2->GetY()*hit1->GetZ()) /
  312. (hit2->GetZ() - hit1->GetZ());
  313. }
  314. }
  315. //__________________________________________________________________
  316. Bool_t TpcLheTrackCuts::IsGoodFoundTrack(TpcLheCMTrack* track) {
  317. //---
  318. if (track->GetNumberOfPoints() < 3)
  319. return kFALSE;
  320. else
  321. return kTRUE;
  322. }
  323. //__________________________________________________________________
  324. Bool_t TpcLheTrackCuts::IsGoodGeantTrack(TpcLheTrack* track) {
  325. //--- Returns true if the given track fulfills all requirements to
  326. //be a "good" track.
  327. return kTRUE;
  328. }
  329. //////////////////////////////////////////////////////////////////////
  330. //////////////////////////////////////////////////////////////////////
  331. //////////////////////////////////////////////////////////////////////
  332. //////////////////////////////////////////////////////////////////////