ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FourHitSeedFinder_find_3_4.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FourHitSeedFinder_find_3_4.cpp
1 #include "FourHitSeedFinder.h"
2 #include <algorithm>
3 
4 
5 using namespace std;
6 
7 
8 void FourHitSeedFinder::findTracks_3_4(vector<SimpleHit3D>& hits, vector<SimpleTrack3D>& tracks, const HelixRange& range)
9 {
10  vector<vector<SimpleHit3D> > layer_sorted;
11  vector<SimpleHit3D> one_vec;
12  layer_sorted.assign(4, one_vec);
13  for(unsigned int i=0;i<hits.size();++i)
14  {
15  layer_sorted[hits[i].layer].push_back(hits[i]);
16  }
17 
18  float k = 0.5*(range.max_k + range.min_k);
19  float dzdl = 0.5*(range.max_dzdl + range.min_dzdl);
20  float p_inv = 3.33333333333333314e+02*k*Bfield_inv*sqrt(1. - dzdl*dzdl);
21  float scatter_weight = 2.;
22 
23  unsigned int l0_size = layer_sorted[0].size();
24  unsigned int l1_size = layer_sorted[1].size();
25  unsigned int l2_size = layer_sorted[2].size();
26  unsigned int l3_size = layer_sorted[3].size();
27 
28  vector<vector<unsigned int> > three_layers;
29  vector<unsigned int> one_hit;one_hit.assign(3,0);
30 
31  vector<double> chi2_hits;
32  vector<unsigned int> tempcomb;
33  tempcomb.assign(3,0);
34  vector<float> kappa_vals;
35  vector<float> dzdl_vals;
36  for(unsigned int i0=0;i0<l0_size;++i0)
37  {
38  for(unsigned int i1=0;i1<l1_size;++i1)
39  {
40  for(unsigned int i2=0;i2<l2_size;++i2)
41  {
42  // have we seen this combo before?
43  tempcomb[0] = layer_sorted[0][i0].index;
44  tempcomb[1] = layer_sorted[1][i1].index;
45  tempcomb[2] = layer_sorted[2][i2].index;
46  sort(tempcomb.begin(), tempcomb.end());
47  set<vector<unsigned int> >::iterator it = combos_3_fail.find(tempcomb);
48  // if we know this combo sucks, ignore it
49  if(it != combos_3_fail.end()){continue;}
50 // it = combos_3_pass.find(tempcomb);
51 // // if we know this combo is good, attach it to the list
52 // if(it != combos_3_pass.end())
53 // {
54 // one_hit[0] = i0;
55 // one_hit[1] = i1;
56 // one_hit[2] = i2;
57 // three_layers.push_back(one_hit);
58 // continue;
59 // }
60 
61  // we haven't seen this combo before ... let's fit it!
62 
63  SimpleTrack3D temp_track;
64  temp_track.hits.push_back(layer_sorted[0][i0]);
65  temp_track.hits.push_back(layer_sorted[1][i1]);
66  temp_track.hits.push_back(layer_sorted[2][i2]);
67 
68  // calculate how much this thing should scatter
69  float rad_diff = detector_radii[2] - detector_radii[1];
70  float scatter = scatter_weight*p_inv*detector_scatter[1]*7.07106781186547462e-01;
71  float dr = scatter*rad_diff;
72  layer_xy_resolution[2] += dr;
73  layer_z_resolution[2] += dr;
74  double chi2 = fitTrack(temp_track, chi2_hits);
75  layer_xy_resolution[2] -= dr;
76  layer_z_resolution[2] -= dr;
77  if(fabs(chi2) > chi2_cut){combos_3_fail.insert(tempcomb);continue;}
78  one_hit[0] = i0;
79  one_hit[1] = i1;
80  one_hit[2] = i2;
81  three_layers.push_back(one_hit);
82  combos_3_pass.insert(tempcomb);
83  kappa_vals.push_back(temp_track.kappa);
84  dzdl_vals.push_back(temp_track.dzdl);
85  }
86  }
87  }
88 
89  tempcomb.resize(4,0);
90  // now three_layers contains a list of three-hit combinations from the inner 3 layers which pass the chi^2 cut
91  // next, add in the hits from the 4th layer
92  for(unsigned int j=0;j<three_layers.size();++j)
93  {
94  for(unsigned int i3=0;i3<l3_size;++i3)
95  {
96  tempcomb[0] = layer_sorted[0][three_layers[j][0]].index;
97  tempcomb[1] = layer_sorted[1][three_layers[j][1]].index;
98  tempcomb[2] = layer_sorted[2][three_layers[j][2]].index;
99  tempcomb[3] = layer_sorted[3][i3].index;
100  sort(tempcomb.begin(), tempcomb.end());
101  set<vector<unsigned int> >::iterator it = combos.find(tempcomb);
102  if(it != combos.end()){continue;}
103  combos.insert(tempcomb);
104 
105  SimpleTrack3D temp_track;
106  temp_track.hits.push_back(layer_sorted[1][three_layers[j][1]]);
107  temp_track.hits.push_back(layer_sorted[2][three_layers[j][2]]);
108  temp_track.hits.push_back(layer_sorted[3][i3]);
109 
110  // calculate how much this thing should scatter
111  float rad_diff = detector_radii[3] - detector_radii[2];
112  float scatter = scatter_weight*p_inv*detector_scatter[2]*7.07106781186547462e-01;
113  float dr = scatter*rad_diff;
114  layer_xy_resolution[2] += dr;
115  layer_z_resolution[2] += dr;
116  double chi2 = fitTrack(temp_track, chi2_hits);
117  layer_xy_resolution[2] -= dr;
118  layer_z_resolution[2] -= dr;
119  if(fabs(chi2) > chi2_cut){continue;}
120 
121 // float sm_k = temp_track.kappa;
122 // float big_k = kappa_vals[j];
123 // if(sm_k > big_k)
124 // {
125 // big_k = temp_track.kappa;
126 // sm_k = kappa_vals[j];
127 // }
128 // float rel_diff = (big_k - sm_k)/big_k;
129 // if(rel_diff > (0.006/k)*0.5){continue;}
130 //
131 // float sm_z = fabs(temp_track.dzdl);
132 // float big_z = fabs(dzdl);
133 // if(sm_z > big_z)
134 // {
135 // big_z = fabs(temp_track.dzdl);
136 // sm_z = fabs(dzdl);
137 // }
138 // rel_diff = (big_z - sm_z)/big_z;
139 // if(rel_diff > (0.006/k)*0.5){continue;}
140 
141  SimpleTrack3D good_track;
142  good_track.hits.push_back(layer_sorted[0][three_layers[j][0]]);
143  good_track.hits.push_back(layer_sorted[1][three_layers[j][1]]);
144  good_track.hits.push_back(layer_sorted[2][three_layers[j][2]]);
145  good_track.hits.push_back(layer_sorted[3][i3]);
146 
147  tracks.push_back(good_track);
148  }
149  }
150 }
151 
152 
153 
154