StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructPairLUT.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructPairLUT.cxx,v 1.4 2010/09/02 21:24:08 prindle Exp $
4  *
5  * Author: Duncan Prindle
6  *
7  **********************************************************************
8  *
9  * Description: If a pair of tracks gets too close it may be lost (due to splitting
10  * or merging). Separation distance depends on reference radius, radii
11  * of the two tracks, difference in phi, difference in eta and charge
12  * signs of the tracks. (Actually, there is an eta dependence, but it is small)
13  * Create Look Up Tables (LUT) by integrating over eta and reference
14  * radius. Use these tables to accept/reject pairs.
15  *
16  *
17  ***********************************************************************/
18 #include "TMath.h"
19 #include "StEStructPairLUT.h"
20 
21 
22 ClassImp(StEStructPairLUT)
23 
24 StEStructPairLUT::StEStructPairLUT() {
25  mR1 = 125.0;
26  mEta1 = 0.0;
27  mPhi1 = 0.0;
28  mSign1 = 1;
29  mR2 = 1250.0;
30  mEta2 = 0.0;
31  mPhi2 = 0.0;
32  mSign2 = 1;
33  mDelEta = 0.25;
34  mDelPhi = 2.5;
35  nDelEta = 25;
36  nDelPhi = 50;
37  mDelXYCut = 5;
38  mDelZCut = 5;
39  mPi = TMath::ACos(-1);
40  mRRefMin = 50;
41  mRRefMax = 200;
42 
43  for (int i=0;i<4;i++) {
44  mDists[i] = 0;
45  }
46  int iR = 0;
47  for (int ir2=0;ir2<10;ir2++) {
48  for (int ir1=0;ir1<=ir2;ir1++) {
49  mCutLS[iR] = 0;
50  mCutUS[iR] = 0;
51  iR++;
52  }
53  }
54 
55 }
56 StEStructPairLUT::~StEStructPairLUT() {
57  for (int i=0;i<4;i++) {
58  delete mDists[i];
59  }
60  int iR = 0;
61  for (int ir2=0;ir2<10;ir2++) {
62  for (int ir1=0;ir1<=ir2;ir1++) {
63  delete mCutLS[iR];
64  delete mCutUS[iR];
65  iR++;
66  }
67  }
68 }
69 double StEStructPairLUT::s (double Ref, double Rad, double eta) {
70  double l = lambda(eta);
71  return Rad*TMath::ACos(1-0.5*TMath::Power(Ref/Rad,2))/TMath::Cos(l);
72 }
73 double StEStructPairLUT::alpha (double Ref, double Rad, double eta, double sign) {
74  double l = lambda(eta);
75  double arc = s(Ref,Rad,eta);
76  return sign*arc*TMath::Cos(l)/Rad;
77 }
78 double StEStructPairLUT::lambda (double eta) {
79  double theta = 2*atan(TMath::Exp(-eta));
80  return 0.5*3.1415926 - theta;
81 }
82 double StEStructPairLUT::delZ (double Ref) {
83  double l1 = lambda(mEta1);
84  double l2 = lambda(mEta2);
85  double s1 = s(Ref,mR1,mEta1);
86  double s2 = s(Ref,mR2,mEta2);
87  return s1*TMath::Sin(l1) - s2*TMath::Sin(l2);
88 }
89 double StEStructPairLUT::delXY (double Ref) {
90  // Compared to tracks.i of Dsv, phi is direction of track (pi/2 different than in tracks.i)
91  // For s small and phi = 0 x is linear in s and y decreases quadratically (i.e. is negative).
92  // (v x B points in -Y direction if v is in X and B is in Z directions)
93  double alpha1 = alpha(Ref,mR1,mEta1,mSign1);
94  double alpha2 = alpha(Ref,mR2,mEta2,mSign2);
95  double x1 = +mSign1*mR1*(TMath::Sin(mPhi1)-TMath::Sin(mPhi1-alpha1));
96  double y1 = -mSign1*mR1*(TMath::Cos(mPhi1)-TMath::Cos(mPhi1-alpha1));
97  double x2 = +mSign2*mR2*(TMath::Sin(mPhi2)-TMath::Sin(mPhi2-alpha2));
98  double y2 = -mSign2*mR2*(TMath::Cos(mPhi2)-TMath::Cos(mPhi2-alpha2));
99  double d = TMath::Sqrt(TMath::Power(x1-x2,2)+TMath::Power(y1-y2,2));
100  mDPhi_Ref = x1*y2-x2*y1;
101 
102  alpha1 = alpha(50,mR1,mEta1,mSign1);
103  alpha2 = alpha(50,mR2,mEta2,mSign2);
104  mX1_50 = +mSign1*mR1*(TMath::Sin(mPhi1)-TMath::Sin(mPhi1-alpha1));
105  mY1_50 = -mSign1*mR1*(TMath::Cos(mPhi1)-TMath::Cos(mPhi1-alpha1));
106  mX2_50 = +mSign2*mR2*(TMath::Sin(mPhi2)-TMath::Sin(mPhi2-alpha2));
107  mY2_50 = -mSign2*mR2*(TMath::Cos(mPhi2)-TMath::Cos(mPhi2-alpha2));
108  mDPhi_50 = mX1_50*mY2_50-mX2_50*mY1_50;
109  if (mDPhi_50*mDPhi_Ref < 0) {
110  return -d;
111  } else {
112  return d;
113  }
114 }
115 void StEStructPairLUT::initHists() {
116  for (int i=0;i<4;i++) {
117  if (mDists[i]) {
118  delete mDists[i];
119  }
120  }
121  char bufferPhi[1024];
122  char bufferEta[1024];
123  sprintf(bufferPhi,"#phi_{1} (#phi_{2}=%f)",mPhi2);
124  sprintf(bufferEta,"#eta_{1} (#eta_{2}=%f)",mEta2);
125  mDists[0] = new TH2D("dists0","#delta#phi vs. Reference radius",nDelPhi,mPhi2-mDelPhi,mPhi2+mDelPhi, 100,50.0,200.0);
126  mDists[0]->GetXaxis()->SetTitle(bufferPhi);
127  mDists[0]->GetYaxis()->SetTitle("R_{ref}");
128  mDists[0]->SetMaximum(+5);
129  mDists[0]->SetMinimum(-5);
130  mDists[1] = new TH2D("dists1","#delta#eta vs. Reference radius",nDelEta,mEta2-mDelEta,mEta2+mDelEta, 100,50.0,200.0);
131  mDists[1]->GetXaxis()->SetTitle(bufferEta);
132  mDists[1]->GetYaxis()->SetTitle("R_{ref}");
133  mDists[1]->SetMaximum(+5);
134  mDists[1]->SetMinimum(-5);
135  mDists[2] = new TH2D("dists2","#delta#phi vs. #delta#eta: |#delta_{XY}| < 5cm and |#delta_{Z}| < 5cm",nDelPhi,-mDelPhi,+mDelPhi, nDelEta,-mDelEta,+mDelEta);
136  mDists[2]->GetXaxis()->SetTitle("#phi_{1}-#phi_{2}");
137  mDists[2]->GetYaxis()->SetTitle("#eta_{1}-#eta_{2}");
138  mDists[3] = new TH2D("dists3","#delta#phi vs. #delta#eta: |#delta_{Z}| < 5cm and tracks cross in #phi between 50cm and R_{ref}",nDelPhi,-mDelPhi,+mDelPhi, nDelEta,-mDelEta,+mDelEta);
139  mDists[3]->GetXaxis()->SetTitle("#phi_{1}-#phi_{2}");
140  mDists[3]->GetYaxis()->SetTitle("#eta_{1}-#eta_{2}");
141 
142  int iR = 0;
143  char buf1[1024], buf2[1024];
144  for (int ir2=0;ir2<10;ir2++) {
145  for (int ir1=0;ir1<=ir2;ir1++) {
146  sprintf(buf1,"LS_%i_%i",ir1,ir2);
147  sprintf(buf2,"cuts for LS, rad1 bin %i, rad2 bin %i",ir1,ir2);
148  if (mCutLS[iR]) delete mCutLS[iR];
149  // 2 extra bins in \delta\phi to account for high p_t tracks with wrong opening angle that stay close.
150  // extra bin in \delta\eta is the bin a 0 (not really extra).
151  mCutLS[iR] = new TH2F(buf1,buf2, 2+nDelPhi,-2*mDelPhi/nDelPhi,mDelPhi, 1+2*nDelEta,-mDelEta,+mDelEta);
152  sprintf(buf1,"US_%i_%i",ir1,ir2);
153  sprintf(buf2,"cuts for US, rad1 bin %i, rad2 bin %i",ir1,ir2);
154  if (mCutUS[iR]) delete mCutUS[iR];
155  mCutUS[iR] = new TH2F(buf1,buf2, 2+nDelPhi,-2*mDelPhi/nDelPhi,mDelPhi, 1+2*nDelEta,-mDelEta,+mDelEta);
156  iR++;
157  }
158  }
159 }
160 void StEStructPairLUT::fillDists() {
161  for (int i=0;i<4;i++) {
162  mDists[i]->Reset();
163  }
164 
165  double dZ, dXY;
166  int ix = 1;
167  double saveEta = mEta1;
168  for (mEta1=mEta2-mDelEta;mEta1<mEta2+mDelEta;mEta1+=2*mDelEta/nDelEta) {
169  int iy = 1;
170  for (double rRef=50;rRef<200;rRef+=150/100.0) {
171  dZ = delZ(rRef);
172  mDists[1]->SetBinContent(ix,iy,dZ);
173  iy++;
174  }
175  ix++;
176  }
177  mEta1 = saveEta;
178  ix = 1;
179  double savePhi = mPhi1;
180  for (mPhi1=mPhi2-mDelPhi;mPhi1<mPhi2+mDelPhi;mPhi1+=2*mDelPhi/nDelPhi) {
181  int iy = 1;
182  for (double rRef=50;rRef<200;rRef+=150/100.0) {
183  dXY = delXY(rRef);
184  mDists[0]->SetBinContent(ix,iy,dXY);
185  iy++;
186  }
187  ix++;
188  }
189  mPhi1 = savePhi;
190 
191  // For each phi1 bin scan over R_ref (from 65cm to 200cm) for dXY.
192  // If |dXY| < 5cm scan over eta1 for this R_ref. For each bin with |dZ| < 5cm increment mDists[2](phi1,eta1)
193  // If dXY < 0 scan over eta1 for this R_ref. For each bin with |dZ| < 5cm increment mDists[3](phi1,eta1)
194  ix = 1;
195  for (mPhi1=mPhi2-mDelPhi;mPhi1<mPhi2+mDelPhi;mPhi1+=2*mDelPhi/nDelPhi) {
196  int iy = 11;
197  for (double rRef=65;rRef<200;rRef+=150/100.0) {
198  dXY = mDists[0]->GetBinContent(ix,iy);
199  if (TMath::Abs(dXY) < 5) {
200  int iz = 1;
201  for (mEta1=mEta2-mDelEta;mEta1<mEta2+mDelEta;mEta1+=2*mDelEta/nDelEta) {
202  dZ = mDists[1]->GetBinContent(iz,iy);
203  if (TMath::Abs(dZ) < 5) {
204  mDists[2]->Fill(mPhi1-mPhi2,mEta1-mEta2,1.0);
205  }
206  iz++;
207  }
208  mEta1 = saveEta;
209  }
210  if (dXY < 0) {
211  int iz = 1;
212  for (mEta1=mEta2-mDelEta;mEta1<mEta2+mDelEta;mEta1+=2*mDelEta/nDelEta) {
213  dZ = mDists[1]->GetBinContent(iz,iy);
214  if (TMath::Abs(dZ) < 5) {
215  mDists[3]->Fill(mPhi1-mPhi2,mEta1-mEta2,1);
216  }
217  iz++;
218  }
219  mEta1 = saveEta;
220  }
221  iy++;
222  }
223  ix++;
224  }
225  mPhi1 = savePhi;
226 }
227 
228 void StEStructPairLUT::integrateEta(TH2F *h) {
229  h->Reset();
230 
231  // Can only depend on \delta\phi. Choose \phi_1 = 0, scan over \phi_2.
232  // For each \phi_2 bin scan over R_ref for dXY.
233  // If |dXY| < mDelXYCut and TMath::Abs(dZ) < mDelZCut we think of this as a merged pair.
234  // If dXY < 0 and TMath::Abs(dZ) < mDelZCut we think of this a a crosTMath::Sing pair.
235  // Approximately only depends on \delta\eta.
236  // Scan over \eta_2 and for each value scan over \eta_1 within some range.
237  // Need to scan over \eta_1 for R_ref where dXY < mDelXYCut to see if TMath::Abs(dZ) is ever < mDelZCut.
238  // Note:
239 
240  double savePhi1 = mPhi1;
241  double savePhi2 = mPhi2;
242  mPhi2 = 0;
243  int ix = 1;
244  for (mPhi2=-1.5*mDelPhi/nDelPhi;mPhi2<mDelPhi;mPhi2+=mDelPhi/nDelPhi) {
245  int iy = 1;
246  for (double rRef=mRRefMin;rRef<mRRefMax;rRef+=150/100.0) {
247  double dXY = delXY(rRef);
248  if (dXY < mDelXYCut) {
249  double saveEta1 = mEta1;
250  double saveEta2 = mEta2;
251  int iz2 = 1;
252  for (mEta2=0;mEta2<1.0;mEta2+=0.1) {
253  int iz1 = 1;
254  for (mEta1=mEta2-mDelEta;mEta1<mEta2+mDelEta;mEta1+=mDelEta/nDelEta) {
255  if (TMath::Abs(mEta1) <= 1) {
256  double dZ = delZ(rRef);
257  if (TMath::Abs(dZ) < mDelZCut) {
258  h->Fill(mPhi2,mEta1-mEta2,1.0);
259  }
260  }
261  iz1++;
262  }
263  iz2++;
264  }
265  mEta1 = saveEta1;
266  mEta2 = saveEta2;
267  }
268  iy++;
269  }
270  ix++;
271  }
272  mPhi1 = savePhi1;
273  mPhi2 = savePhi2;
274 }
275 
276 void StEStructPairLUT::fillLUTs() {
277 
278  // Ten equal steps in 1/R. Maximum value of 1/R is 1/75, close to 150MeV/c.
279  for (int ir=9;ir>=0;ir--) {
280  mRad[ir] = 750.0/(ir+0.5);
281  }
282 
283  int iR = 0;
284  // Choose mR1 >= mR2 => for LS need phi2 > phi1. (and can choose phi1 = 0)
285  // For US case choose mSign2 < 0 => phi2 > phi1 (and can choose phi1 = 0)
286  for (int ir2=0;ir2<10;ir2++) {
287  for (int ir1=0;ir1<=ir2;ir1++) {
288  mSign1 = +1;
289  mR1 = mRad[ir1];
290  mSign2 = +1;
291  mR2 = mRad[ir2];
292  integrateEta(mCutLS[iR]);
293  mSign1 = -1;
294  mR1 = mRad[ir1];
295  mSign2 = +1;
296  mR2 = mRad[ir2];
297  integrateEta(mCutUS[iR]);
298  iR++;
299  }
300  }
301 }
302 int StEStructPairLUT::cut(double curvature1, double curvature2, double delPhi, double delEta) {
303  // For high enough pt we can have phi1-phi2 with wrong sign but tracks still merged.
304  // We are trying to offset by 2 bins in delPhi to catch this.
305 
306  // Probably called with phi1-phi2 as argument. They might be on opposite sides of cut.
307  while (delPhi > +2*mPi) delPhi -= 2*mPi;
308  while (delPhi < -2*mPi) delPhi += 2*mPi;
309 
310  // Calculate indices into Look Up Table.
311  if (TMath::Abs(delEta) > mDelEta) return 0;
312  if (TMath::Abs(delPhi) > mDelPhi) return 0;
313  int idEta = int( delEta / (mDelEta/nDelEta) );
314  int idPhi = int( delPhi / (mDelPhi/nDelPhi) );
315  int ir1 = int(750*TMath::Abs(curvature1) -0.5);
316  if (ir1 > 9) ir1 = 9;
317  int ir2 = int(750*TMath::Abs(curvature2) -0.5);
318  if (ir2 > 9) ir2 = 9;
319 
320  // LUT has R1 >= R2.
321  // Calculate iR and adjust idEta and idPhi appropriately
322  int iR, is1, is2;
323  if (TMath::Abs(curvature1) < TMath::Abs(curvature2)) {
324  iR = ir2*(ir2+1)/2 + ir1;
325  idEta = nDelEta + idEta;
326  is1 = (curvature1 < 0) ? -1 : 1;
327  is2 = (curvature2 < 0) ? -1 : 1;
328  } else {
329  iR = ir1*(ir1+1)/2 + ir2;
330  idEta = nDelEta - idEta;
331  idPhi = -idPhi;
332  is1 = (curvature2 < 0) ? -1 : 1;
333  is2 = (curvature1 < 0) ? -1 : 1;
334  }
335 
336  if (is1 > 0) {
337  if (is2 > 0) {
338  // Need delPhi = phi_1 - phi_2 < 0
339  // Rely on GetBinContent returning 0 for out of bound index.
340  // Remember, root starts histograms with bin=1, _not_ 0.
341  return int(mCutLS[iR]->GetBinContent(-idPhi+3,idEta+1));
342  } else {
343  return int(mCutUS[iR]->GetBinContent(idPhi+3,idEta+1));
344  }
345  } else {
346  if (is2 > 0) {
347  // Need delPhi = phi_1 - phi_2 < 0
348  return int(mCutUS[iR]->GetBinContent(-idPhi+3,idEta+1));
349  } else {
350  return int(mCutLS[iR]->GetBinContent(idPhi+3,idEta+1));
351  }
352  }
353 }
Generated by   doxygen 1.8.5