simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
TOFExtrapolation.h
1 using namespace std;
2 using namespace io;
3 using namespace utl;
4 using namespace evt;
5 using namespace evt::rec;
6 using namespace det;
7 
8 //__________________________________________________________________________________________________________________________
9 vector<double> GetTOF_scint_fitExtrapolated_xy(const RecEvent& recEvent,
10  const Track& track, const int& TOF_scint_id,
11  const det::TOFWall& TOFwall,
12  const det::MagneticFieldTracker& tracker,
13  bool debug = false)
14 {
16 
17  if (debug == true)
18  {
19  cerr<<"\nReached till GetTOF_scint_fitExtrapolated_xy().";
20  cerr<<"\nNumber of points on track in the MTPC = "<<track.GetNumberOfClusters(TrackConst::eMTPC)<<endl;
21  }
22 
23  const utl::Point& last_point_on_track = track.GetLastPointOnTrack();
24  vector <utl::Point> clusters;
25 
26 
27  // for (ClusterIndexIterator ClusterIter = track.ClustersBegin(); ClusterIter != track.ClustersEnd(); ++ ClusterIter)
28  // for (ClusterIndexIterator ClusterIter = track.ClustersEnd() ; iPoint < 10 ; ++iPoint)
29  int iPoint = 0;
30  for (ClusterIndexIterator ClusterIter = track.ClustersEnd(); ClusterIter != track.ClustersBegin(); --ClusterIter)
31  {
32  if (clusters.size() > 9)
33  break;
34 
35  if (iPoint == 0) --ClusterIter;
36  ++iPoint;
37 
38  if (debug == true)
39  cerr<<"\nFor iPoint = "<<iPoint;
40 
41  const Cluster& cluster = recEvent.Get(*ClusterIter);
42  clusters.push_back(cluster.GetPosition());
43 
44  if (debug == true)
45  {
46  cerr<<": "<<cluster.GetPosition().GetX()<<", "<<cluster.GetPosition().GetY()<<", "<<cluster.GetPosition().GetZ();
47 
48  if (last_point_on_track == cluster.GetPosition())
49  cerr<<". Found last point on track, corresponding to iPoint = "<<iPoint;
50  }
51  }
52 
53 
54  TGraph * graph_track_XZ = new TGraph(clusters.size());
55  TGraph * graph_track_YZ = new TGraph(clusters.size());
56  double largest_gap_in_MTPCtrk = 0;
57  for (unsigned int iclust = 0; iclust < clusters.size(); iclust++)
58  {
59  graph_track_XZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetX());
60  graph_track_YZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetY());
61 
62  if (iclust < clusters.size() - 1)
63  {
64  double cluster_diff_x = clusters[iclust].GetX() - clusters[iclust+1].GetX();
65  double cluster_diff_y = clusters[iclust].GetY() - clusters[iclust+1].GetY();
66  double cluster_diff_z = clusters[iclust].GetZ() - clusters[iclust+1].GetZ();
67  double cluster_separation = sqrt(cluster_diff_x*cluster_diff_x + cluster_diff_y*cluster_diff_y + cluster_diff_z*cluster_diff_z);
68 
69  if (debug == true)
70  cout<<"\ncluster_separation = "<<cluster_separation;
71 
72  if (largest_gap_in_MTPCtrk < cluster_separation)
73  largest_gap_in_MTPCtrk = cluster_separation;
74  }
75  }
76  if (debug == true)
77  cout<<"\nlargest_gap_in_MTPCtrk = "<<largest_gap_in_MTPCtrk<<endl;
78 
79 
80  TF1 * lin_fit_XZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);
81  TF1 * lin_fit_YZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);
82 
83  graph_track_XZ->Fit(lin_fit_XZ, "Q"); //"V"
84  graph_track_YZ->Fit(lin_fit_YZ, "Q");
85  // TFitResultPtr r1 = graph_track_XZ->Fit(lin_fit_XZ, "SQ");
86  // TFitResultPtr r2 = graph_track_YZ->Fit(lin_fit_YZ, "SQ");
87  // cout<<"Chi-sq/ndf: Fit XZ = "<<r1->Chi2()/r1->Ndf()<<", Fit YZ = "<<r2->Chi2()/r2->Ndf()<<endl;
88  // delete r1;
89  // delete r2;
90 
91  const det::TOFScintillator& Iscint = TOFwall.GetScintillator(TOF_scint_id);
92  double extrapolationPosition = Iscint.GetCenterPosition().GetZ() - 2.3/2;
93  double fitX = lin_fit_XZ->Eval(extrapolationPosition);
94  double fitY = lin_fit_YZ->Eval(extrapolationPosition);
95 
96 
98  utl::Vector stopMomentum; // this variable will keep track momentum when the track reaches ToF wall
99  utl::Point stopPosition; // this variable will keep track position when the track reaches ToF wall
100  tracker.TrackToZ(extrapolationPosition, track.GetCharge(), track.GetFirstPointOnTrack(), track.GetMomentum(), stopPosition, stopMomentum);
101  double TOFx = stopPosition.GetX();
102  double TOFy = stopPosition.GetY();
103 
104  delete graph_track_XZ;
105  delete graph_track_YZ;
106  delete lin_fit_XZ;
107  delete lin_fit_YZ;
108 
109  vector<double> XY = {fitX, fitY, TOFx, TOFy, largest_gap_in_MTPCtrk};
110  return XY;
111 }
112 
113 
114 
115 
116 //__________________________________________________________________________________________________________________________
117 int GetClusterFitExtrapolatedTOFPixelID(const det::TOFWall& TOFwall,
118  const RecEvent& recEvent,
119  const Track& track,
120  double& TOF_pixel_Z,
121  const bool& debug = false)
122 {
123  if (debug == true)
124  {
125  cerr<<"\n\nNumber of points on track in the MTPC = "<<track.GetNumberOfClusters(TrackConst::eMTPC)<<endl;
126  }
127 
128  const utl::Point& last_point_on_track = track.GetLastPointOnTrack();
129  vector<Point> clusters;
130 
131 
133  int iPoint = 0;
134  for (ClusterIndexIterator ClusterIter = track.ClustersEnd(); ClusterIter != track.ClustersBegin(); --ClusterIter)
135  {
136  if (clusters.size() > 9)
137  break;
138 
139  if (debug == true)
140  cerr<<"\nClusterIter = "<<*ClusterIter;
141 
142  if (iPoint == 0) --ClusterIter;
143  ++iPoint;
144 
145  if (debug == true)
146  cerr<<"\nFor iPoint = "<<iPoint;
147 
148  const Cluster& cluster = recEvent.Get(*ClusterIter);
149  clusters.push_back(cluster.GetPosition());
150 
151  if (debug == true)
152  {
153  cerr<<": "<<cluster.GetPosition().GetX()<<", "<<cluster.GetPosition().GetY()<<", "<<cluster.GetPosition().GetZ();
154 
155  if (last_point_on_track == cluster.GetPosition())
156  cerr<<". Found last point on track, corresponding to iPoint = "<<iPoint;
157  }
158  }
159 
160  if (debug == true) cerr<<"\n";
161 
162 
164  /*
165  int pointsIter = 0;
166  for (ClusterIndexIterator ClusterIter = track.ClustersEnd(); pointsIter < 10; ++pointsIter)
167  {
168  if (debug == true)
169  cerr<<"\nClusterIter = "<<*ClusterIter;
170 
171  const Cluster& cluster = recEvent.Get(*(--ClusterIter));
172  clusters.push_back(cluster.GetPosition());
173 
174  if (debug == true)
175  {
176  cerr<<"\nFor pointsIter = "<<pointsIter;
177  cerr<<": "<<cluster.GetPosition().GetX()<<", "<<cluster.GetPosition().GetY()<<", "<<cluster.GetPosition().GetZ();
178 
179  if (last_point_on_track == cluster.GetPosition())
180  cerr<<". Found last point on track, corresponding to pointsIter = "<<pointsIter;
181  }
182  }
183 
184  if (debug == true) cerr<<"\n";
185  */
186 
187 
188  TGraph *graph_track_XZ = new TGraph(clusters.size());
189  TGraph *graph_track_YZ = new TGraph(clusters.size());
190  for(unsigned int iclust = 0; iclust < clusters.size(); iclust++)
191  {
192  graph_track_XZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetX());
193  graph_track_YZ->SetPoint((int)iclust, clusters.at(iclust).GetZ(), clusters.at(iclust).GetY());
194  }
195  TF1 *lin_fit_XZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);
196  TF1 *lin_fit_YZ = new TF1("track_XZ", "[0]*x+[1]", 650, 830);
197 
198  graph_track_XZ->Fit(lin_fit_XZ, "Q"); //"V"
199  graph_track_YZ->Fit(lin_fit_YZ, "Q");
200  // TFitResultPtr r1 = graph_track_XZ->Fit(lin_fit_XZ, "SQ");
201  // TFitResultPtr r2 = graph_track_YZ->Fit(lin_fit_YZ, "SQ");
202  // cout<<"Chi-sq/ndf: Fit XZ = "<<r1->Chi2()/r1->Ndf()<<", Fit YZ = "<<r2->Chi2()/r2->Ndf()<<endl;
203  // delete r1;
204  // delete r2;
205 
206 
207  int nHasTOF = 0;
208  vector<int> vTOFHitIndex;
209  vector<double> vTOF_pixel_Z;
210 
211  for(int hitindex = 1; hitindex <= 891; ++hitindex)
212  {
213  const det::TOFScintillator& Iscint = TOFwall.GetScintillator(hitindex);
214  double extrapolationPositionZ = Iscint.GetCenterPosition().GetZ() - 2.3/2;
215  double extrapolationPositionX = lin_fit_XZ->Eval(extrapolationPositionZ);
216  double extrapolationPositionY = lin_fit_YZ->Eval(extrapolationPositionZ);
217 
218  if (extrapolationPositionX < Iscint.GetCenterPosition().GetX()-Iscint.GetWidth()/2.) continue;
219  if (extrapolationPositionX > Iscint.GetCenterPosition().GetX()+Iscint.GetWidth()/2.) continue;
220  if (extrapolationPositionY < Iscint.GetCenterPosition().GetY()-3.4/2) continue;
221  if (extrapolationPositionY > Iscint.GetCenterPosition().GetY()+3.4/2) continue;
222 
223  // cout<<hitindex<<endl;
224  ++nHasTOF;
225  vTOFHitIndex.push_back(hitindex);
226  vTOF_pixel_Z.push_back(extrapolationPositionZ);
227  }
228 
229 
230  delete graph_track_XZ;
231  delete graph_track_YZ;
232  delete lin_fit_XZ;
233  delete lin_fit_YZ;
234 
235  if (nHasTOF == 0) // no hits
236  {
237  TOF_pixel_Z = -1;
238  return -1;
239  }
240 
241  else if (nHasTOF == 1)
242  {
243  TOF_pixel_Z = vTOF_pixel_Z[0];
244  return vTOFHitIndex[0];
245  }
246 
247  else //(nHasTOF > 1)
248  {
249  TOF_pixel_Z = -9;
250  return -9;
251  }
252 }
253 
254 
255 
256 
257 //__________________________________________________________________________________________________________________________
258 int GetBFieldExtrapolatedTOFPixelID(const det::TOFWall& TOFwall,
259  const det::MagneticFieldTracker& tracker,
260  const utl::Point& firstPointOnTrack,
261  const utl::Vector& TrkMomentum,
262  const int& track_q)
263 {
264  int nHasTOF = 0;
265  vector<int> vTOFHitIndex;
266 
267  for(int hitindex = 1; hitindex <= 891; ++hitindex)
268  {
269  const det::TOFScintillator& Iscint = TOFwall.GetScintillator(hitindex);
270  double extrapolationPositionZ = Iscint.GetCenterPosition().GetZ() - 2.3/2;
271  // cout<<hitindex<<", "<<extrapolationPositionZ<<endl;
272 
273  utl::Vector stopMomentum;
274  utl::Point stopPosition;
275 
276  tracker.TrackToZ(extrapolationPositionZ, track_q, firstPointOnTrack, TrkMomentum, stopPosition, stopMomentum); //needs track momentum, not vtxTrack momentum.
277  double extrapolationPositionX = stopPosition.GetX();
278  double extrapolationPositionY = stopPosition.GetY();
279 
280  if (extrapolationPositionX < Iscint.GetCenterPosition().GetX()-Iscint.GetWidth()/2.) continue;
281  if (extrapolationPositionX > Iscint.GetCenterPosition().GetX()+Iscint.GetWidth()/2.) continue;
282  if (extrapolationPositionY < Iscint.GetCenterPosition().GetY()-3.4/2) continue;
283  if (extrapolationPositionY > Iscint.GetCenterPosition().GetY()+3.4/2) continue;
284 
285  // cout<<hitindex<<endl;
286  ++nHasTOF;
287  vTOFHitIndex.push_back(hitindex);
288  }
289 
290 
291  if (nHasTOF == 0) // no hits
292  return -1;
293 
294  else if (nHasTOF == 1)
295  return vTOFHitIndex[0];
296 
297  else //(nHasTOF > 1)
298  return -9;
299 }
Generated by   doxygen 1.8.5