StRoot  1
FtfDedx.cxx
1 /*:>-------------------------------------------------------------------
2 **: FILE: FtfDedx.cxx
3 **: HISTORY:
4 **:<------------------------------------------------------------------*/
5 #include "FtfDedx.h"
6
7 #include <algorithm>
8
9 //***********************************************************************
10 // constructor
11 //***********************************************************************
12 FtfDedx::FtfDedx(l3CoordinateTransformer *trafo, float cutLow, float cutHigh, float driftLoss)
13 {
14  fCoordTransformer = trafo;
15  fCutLow = cutLow;
16  fCutHigh = cutHigh;
17  fNTruncate = 0;
18  fDriftLoss = driftLoss;
19  // initialize unit vector perpendicular to rows
20  // for each sector
21  for (int sector=0; sector<24; sector++) {
22  // caution: sector>12 needs x->-x and y->y (east side!)
23  fUnitVec[sector].x = fCoordTransformer->GetSectorSin(sector);
24  if (sector+1>12) fUnitVec[sector].x = -1 * fUnitVec[sector].x;
25  fUnitVec[sector].y = fCoordTransformer->GetSectorCos(sector);
26  }
27 }
28
29
30 //***********************************************************************
31 // dEdx calculation using Truncated Mean algorithm
32 // averaging method: dE/dx = (sum of q/s)/(# of points)
33 //
34 // author: christof
35 //***********************************************************************
36 int FtfDedx::TruncatedMean(FtfTrack *track) {
37
39  int nPointsInArray;
40
41  // dx calculation:
42  // straight-line approx. in x-y-plane
43  // assume perfect hit position on track
44  // :=> don't calculate intersection of circle and row
45  // 1/cosl corr. factor for dip angle
46
47 // printf("Track: id %i, pt %f, nHits %i, q %i, tanl %f\n",
48 // fTrack->id, fTrack->pt, fTrack->nHits, fTrack->q, fTrack->tanl);
49
50  // calculate dip angle correction factor
51  double cosl = 1 / (sqrt( 1 + (double) track->tanl*track->tanl));
52
53 // printf(" cosl %f\n", cosl);
54  if ( cosl==0 ) return 0;
55
56  // calculate center of circle
57  double tPhi0 = track->psi + track->q * track->getPara()->bFieldPolarity * pi/2;
58  if ( tPhi0 > 2*pi ) tPhi0 = tPhi0 - 2*pi;
59  else if ( tPhi0 < 0. ) tPhi0 = tPhi0 + 2*pi;
60
61  double x0 = track->r0 * cos(track->phi0);
62  double y0 = track->r0 * sin(track->phi0);
63  double rr = track->pt / ( bFactor *track->getPara()->bField );
64  double xc = x0 - rr * cos(tPhi0);
65  double yc = y0 - rr * sin(tPhi0);
66
67 // printf(" xc %f yc %f\n", xc, yc);
68
69  nPointsInArray = 0;
70
71  //now loop over hits
72  for ( FtfHit *ihit = (FtfHit *)track->firstHit ;
73  ihit != 0 ;
74  ihit = (FtfHit *)ihit->nextTrackHit) {
75
76  // discard hits with wrong charge information
77  // i.e. hits of one-pad cluster or merged cluster
78  if (ihit->flags) continue;
79
80  // calculate direction of the tangent of circle at position
81  // of given hit:
82  // - rotate radius vector to hit by -90 degrees
83  // matrix ( 0 1 )
84  // ( -1 0 )
85  // - divide by length to get unity vector
86  struct christofs_2d_vector tangent;
87  tangent.x = (ihit->y - yc)/rr;
88  tangent.y = -(ihit->x - xc)/rr;
89  //printf(" tangent x %f y %f\n", tangent.x, tangent.y);
90
91  // get crossing angle by calculating dot product of
92  // tanget and unity vector perpendicular to row (look-up table)
93  double cosCrossing = fabs(tangent.x * fUnitVec[ihit->sector-1].x + tangent.y * fUnitVec[ihit->sector-1].y);
94
98
99  // ==> finally dx:
100  if ( cosCrossing==0 ) continue;
101  double dx = padLength / (cosCrossing * cosl);
102
103  // drift length correctrion: d_loss [m^-1], l_drift [cm]
104  // q_corr = q_meas / ( 1 - l_drift * d_loss/100 )
105  double scaleDrift = 1 - fabs(ihit->z) * fDriftLoss/100;
106
107
108  // first shot: de_scale as implemented in tph.F,
109  // i.e. gas gain, wire-to-pad coupling, etc.
110  // !! hard coded, has to come from the db somewhere sometime !!
111  double deScale;
112  if (ihit->row<14) deScale = 5.0345021e-9;
113  else deScale = 1.4234702e-8;
114
115  fDedxArray[nPointsInArray] = ihit->q * deScale / ( dx * scaleDrift );
116  nPointsInArray++;
117
118 // printf(" ihit row %i x %f y %f z %f q %d q_scaled %e scale %e flags %d\n",
119 // ihit->row, ihit->x, ihit->y, ihit->z, ihit->q, ihit->q*deScale, deScale, ihit->flags);
120 // printf(" cosCros: %f, cosl: %f ===>> dx %f\n", cosCrossing, cosl, dx);
121
122  }
123
124  // sort dEdxArray
125  sort( fDedxArray, fDedxArray+nPointsInArray );
126
127  // calculate absolute cuts
128  int cLow = (int) floor(nPointsInArray * fCutLow);
129  int cHigh = (int) floor(nPointsInArray * fCutHigh);
130
131  track->nDedx = 0;
132  track->dedx = 0;
133
134  for (int i=cLow; i<cHigh; i++) {
135  track->nDedx++;
136  track->dedx += fDedxArray[i];
137  }
138
139  if (track->nDedx>0) track->dedx = track->dedx/track->nDedx;
140
141 // printf(" %e %i\n", track->dedx, track->nDedx);
142
143  return 0;
144 }
Definition: FtfHit.h:16