FairRoot/PandaRoot
PndSttHitCorrector.cxx
Go to the documentation of this file.
1 /*
2  * PndSttHitCorrector.cxx
3  *
4  * Created on: May 20, 2014
5  * Author: schumann
6  */
7 
8 #include "PndSttHitCorrector.h"
9 #include "PndSttHit.h"
10 #include "TVector3.h"
11 
12 using namespace std;
13 
14 ClassImp(PndSttHitCorrector);
15 
16 void PndSttHitCorrector::CorrectHits() {
17 
18  if (fVerbose > 1)
19  cout << "PndSttHitCorrector::CorrectHits" << endl;
20 
21  //for each tube with 2 active neighbors
22  for (size_t i = 0; i < fSeparations[2].size(); ++i) {
23  int actualTubeId = fSeparations[2][i];
24 
25  //continue if straw is skewed
26  if (fStrawMap->IsSkewedStraw(actualTubeId))
27  continue;
28 
29  if (fVerbose > 2)
30  cout << "Calculation for tube: " << actualTubeId << endl;
31 
32  //vector of tangent angles for each neighbor
33  vector<vector<double> > angles;
34 
35  //calculate if both neighbors are unskewed
36  if (!fStrawMap->IsSkewedStraw(fHitNeighbors[actualTubeId][0])
37  && !fStrawMap->IsSkewedStraw(fHitNeighbors[actualTubeId][1])) {
38 
39  //for each neighbor of actual tube
40  for (size_t j = 0; j < fHitNeighbors[actualTubeId].size(); j++) {
41  int neighborTubeId = fHitNeighbors[actualTubeId][j];
42 
43  if (fVerbose > 3)
44  cout << "NeighborTube: " << neighborTubeId << " : " << endl;
45 
46  //calculate angles of tangents between the isochrones
47  angles.push_back(
48  CalculateTangentAngles(
49  (PndSttHit*) fHits[fMapTubeIdToHit[actualTubeId]],
50  (PndSttHit*) fHits[fMapTubeIdToHit[neighborTubeId]]));
51  }
52 
53  //search for the best combination of tangents and store the angles of the best one
54  set<double> bestPhi = GetBestCombinatedPhi(angles);
55 
56  if (bestPhi.size() == 2) {
57  double absDiff = abs(*bestPhi.begin() - *(++bestPhi.begin()));
58  if (absDiff < fDeltaDiff) {
59 
60  set<double> tmp;
61  tmp.insert(
62  GetAverageOfAngles(*bestPhi.begin(),
63  *(++bestPhi.begin())));
64  fTangentAngles[actualTubeId] = tmp;
65  } else {
66  //store both values
67  fTangentAngles[actualTubeId] = bestPhi;
68  }
69  } else if (bestPhi.size() == 1) {
70  fTangentAngles[actualTubeId] = bestPhi;
71 
72  } else if (fVerbose > 2)
73  cout << "Found no good combination of tangent angles for tube "
74  << actualTubeId << endl;
75 
76  }
77  } //Initialization of possible tangent angles for tubes with two active neighbors done
78 
79 //calculate unambiguous tangent angles of cells with 3 active neighbors
80  if (fVerbose > 2)
81  cout << "Search for special positioning of 3 active neighbors" << endl;
82  map<int, set<double> > tmpTangentAngles;
83 
84  int midNeighbor;
85  double phi1, phi2;
86  for (size_t i = 0; i < fSeparations[3].size(); ++i) {
87  int currentTube = fSeparations[3].at(i);
88 
89  //if currentTube is not in fTangentAngles
90  if (fTangentAngles.find(currentTube) == fTangentAngles.end()) {
91 
92  //check constellation of neighbors
93  midNeighbor = GetMiddleHitNeighbor(currentTube);
94  if (midNeighbor != 0) {
95  //calculate angle of the connecting vector between the tubes
96  phi1 = fGeometryMap->GetAngleBetweenTubes(currentTube,
97  midNeighbor);
98  fTangentAngles[currentTube].insert(phi1);
99 
100  phi2 = phi1 + TMath::Pi();
101  if (phi2 > 2 * TMath::Pi())
102  phi2 -= 2 * TMath::Pi();
103  //update the neighbor
104  fTangentAngles[midNeighbor].insert(phi2);
105  if (fVerbose > 3)
106  cout << "current tube: " << currentTube << ", midNeighbor: "
107  << midNeighbor << ", angles: "
108  << phi1 * TMath::RadToDeg() << ", "
109  << phi2 * TMath::RadToDeg() << endl;
110  }
111  }
112  }
113 
114  if (fVerbose > 2) {
115  cout << "Initialization done. Current fTangentAngles: " << endl;
116  PrintTangentAngles();
117  cout << "Adaption of angles: " << endl;
118  }
119 
120  /* Adapt angles of the unambiguous tubes. Approach: take tubes with 2, 3 or 4 neighbors into account
121  * and search for unambiguous neighbors, that form a straight line*/
122 
123  set<int> tubesToCheck;
124 
125  for (int i = 2; i < 5; ++i) {
126  for (size_t j = 0; j < fSeparations[i].size(); ++j) {
127  int actualTubeId = fSeparations[i].at(j);
128  if (fTangentAngles.find(actualTubeId) == fTangentAngles.end()
129  || fTangentAngles[actualTubeId].size() != 1) {
130  //if two neighbors and ambiguous or if more neighbors and angles were not calculated yet
131  tubesToCheck.insert(actualTubeId);
132  }
133 
134  }
135  }
136 
137  set<int> tmpTubesToCheck;
138  bool anglesChanged = true;
139 
140 // while unambiguous tangent angles were found
141  while (anglesChanged) {
142 
143  if (fVerbose > 3) {
144  cout << "start loop - adaption of angles" << endl;
145  cout << "tubes to Check: " << tubesToCheck.size() << endl;
146  }
147 
148  anglesChanged = false;
149  tmpTubesToCheck.clear();
150 
151  for (set<int>::iterator it = tubesToCheck.begin();
152  it != tubesToCheck.end(); ++it) {
153 
154  //must be checked because neighbor of actual tube could be changed in one situation
155  if (fTangentAngles.find(*it) == fTangentAngles.end()
156  || fTangentAngles[*it].size() != 1) {
157 
158  int actualTubeId = *it;
159  pair<int, int> straightNeighbors = GetStraightNeighbors(
160  actualTubeId);
161 
162  if (fVerbose > 3)
163  cout << "currentTube: " << actualTubeId << endl;
164 
165  //found straight neighbors?
166  if (straightNeighbors.first != 0) {
167  //get unambiguous neighbor
168  //search for unambiguous neighbor
169  int unambiguousNeighbor = 0;
170 
171  if (fHitNeighbors[actualTubeId].size() == 3) {
172  //get neighbor in straight line that lies next to the other neighbor (which is not in straight line)
173  int thirdNeighbor;
174  for (size_t j = 0; j < fHitNeighbors[actualTubeId].size();
175  ++j) {
176  if (fHitNeighbors[actualTubeId].at(j)
177  != straightNeighbors.first
178  && fHitNeighbors[actualTubeId].at(j)
179  != straightNeighbors.second)
180  thirdNeighbor = fHitNeighbors[actualTubeId].at(
181  j);
182  }
183 
184  //set unambiguous neighbor
185  if (AreHitNeihbors(straightNeighbors.first,
186  thirdNeighbor)
187  && HasUnambiguousAngle(
188  straightNeighbors.first)) {
189  unambiguousNeighbor = straightNeighbors.first;
190  } else if (AreHitNeihbors(straightNeighbors.second,
191  thirdNeighbor)
192  && HasUnambiguousAngle(
193  straightNeighbors.second)) {
194  unambiguousNeighbor = straightNeighbors.second;
195  }
196 
197  } else {
198  //acutalTube has 2 or 4 hit-neighbors
199 
200  //track passed the isochrones of the straight neighbors on the same side
201  if (HasUnambiguousAngle(straightNeighbors.first)) {
202 
203  unambiguousNeighbor = straightNeighbors.first;
204 
205  } else if (HasUnambiguousAngle(
206  straightNeighbors.second)) {
207 
208  unambiguousNeighbor = straightNeighbors.second;
209  }
210  }
211 
212  if (unambiguousNeighbor != 0) {
213  if (fVerbose > 3)
214  cout << "unambiguous neighbor: "
215  << unambiguousNeighbor << endl;
216 
217  //for checking the radius of the isochrones
218  PndSttHit* actualHit =
219  (PndSttHit*) fHits[fMapTubeIdToHit[actualTubeId]];
220  PndSttHit* firstNeigh =
221  (PndSttHit*) fHits[fMapTubeIdToHit[straightNeighbors.first]];
222  PndSttHit* secondNeigh =
223  (PndSttHit*) fHits[fMapTubeIdToHit[straightNeighbors.second]];
224 
225  //double keep; //[R.K. 01/2017] unused variable?
226  double similarAngle;
227 
228  if (actualHit->GetIsochrone()
229  > firstNeigh->GetIsochrone()
230  || actualHit->GetIsochrone()
231  > secondNeigh->GetIsochrone()) {
232  //isochrone of actual tube is not the smallest one
233  //assumption: track passes the isochrones along the same side --> keep same/similar angle
234  similarAngle =
235  *fTangentAngles[unambiguousNeighbor].begin();
236  if (fVerbose > 3)
237  cout
238  << "isochrone is not the smallest one --> keep similar"
239  << endl;
240 
241  } else if (actualHit->GetIsochrone()
242  / min(secondNeigh->GetIsochrone(),
243  firstNeigh->GetIsochrone())
244  >= fIsochroneEquality) {
245  //isochrones of actual tube and neighbors are nearly equal, track does not change the side
246  // --> keep similar
247 
248  similarAngle =
249  *fTangentAngles[unambiguousNeighbor].begin();
250 
251  if (fVerbose > 3)
252  cout << "nearly equal isochrones for tube "
253  << actualTubeId << ": actual "
254  << actualHit->GetIsochrone()
255  << ", second neigh "
256  << secondNeigh->GetIsochrone()
257  << ", first neigh: "
258  << firstNeigh->GetIsochrone()
259  << "--> keep similar" << endl;
260  } else {
261  //isochrone of actual tube is the smallest one
262  //assumption: track changes the side. Straight neighbor with the biggest isochrone must have
263  // the same/similar angle. The other neighbor the opposing one.
264 
265  int biggestNeigh, otherNeigh;
266  if (firstNeigh->GetIsochrone()
267  > secondNeigh->GetIsochrone()) {
268  biggestNeigh = straightNeighbors.first;
269  otherNeigh = straightNeighbors.second;
270  } else {
271  biggestNeigh = straightNeighbors.second;
272  otherNeigh = straightNeighbors.first;
273  }
274 
275  if (fVerbose > 3) {
276  cout << "Isochrones (actual, first, second): "
277  << actualHit->GetIsochrone() << ", "
278  << firstNeigh->GetIsochrone() << ", "
279  << secondNeigh->GetIsochrone() << endl;
280  cout
281  << "isochrone is the smallest one, biggest neighbor: "
282  << biggestNeigh << endl;
283  }
284 
285  if (unambiguousNeighbor == biggestNeigh) {
286  //kepp same/similar angle
287  similarAngle =
288  *fTangentAngles[unambiguousNeighbor].begin();
289 
290  //update angle of the other neighbor (if neccessary)
291  if (fTangentAngles.find(otherNeigh)
292  == fTangentAngles.end()
293  || fTangentAngles[otherNeigh].size()
294  != 1) {
295 
296  if (fVerbose > 3)
297  cout << "keep similar, change neighbor "
298  << otherNeigh << " to opposite"
299  << endl;
300 
301  double otherAngle = similarAngle
302  + TMath::Pi();
303  if (otherAngle >= 2 * TMath::Pi())
304  otherAngle -= 2 * TMath::Pi();
305  KeepBestAngle(otherNeigh, otherAngle);
306 
307  }
308 
309  } else {
310  //keep opposing angle
311  similarAngle =
312  *fTangentAngles[unambiguousNeighbor].begin()
313  + TMath::Pi();
314  if (similarAngle >= 2 * TMath::Pi())
315  similarAngle -= 2 * TMath::Pi();
316 
317  if (fTangentAngles.find(otherNeigh)
318  == fTangentAngles.end()
319  || fTangentAngles[otherNeigh].size()
320  != 1) {
321  if (fVerbose > 3)
322  cout
323  << "keep opposite, change neighbor "
324  << otherNeigh << " to similar"
325  << endl;
326  KeepBestAngle(otherNeigh,
327  *fTangentAngles[unambiguousNeighbor].begin());
328  }
329  }
330 
331  }
332 
333  //store angle for actualTube in fTangentAngles
334  KeepBestAngle(actualTubeId, similarAngle);
335 
336  //start next loop!
337  anglesChanged = true;
338  if (fVerbose > 3)
339  cout << "angles changed for " << actualTubeId
340  << endl;
341 
342  } else {
343  //try in next loop
344  tmpTubesToCheck.insert(actualTubeId);
345  if (fVerbose > 3)
346  cout << "try " << actualTubeId << " in next loop"
347  << endl;
348  }
349  } else if (fTangentAngles.find(actualTubeId)
350  != fTangentAngles.end()
351  && fHitNeighbors[actualTubeId].size() == 2
352  && !AreHitNeihbors(fHitNeighbors[actualTubeId].at(0),
353  fHitNeighbors[actualTubeId].at(1))) {
354  // tube has two neighbors and two possible angles
355 
356  if (fVerbose > 3)
357  cout << "handle zigzag for tube " << actualTubeId
358  << endl;
359 
360  int firstNeigh = fHitNeighbors[actualTubeId].at(0);
361  int secondNeigh = fHitNeighbors[actualTubeId].at(1);
362  int unambiguousNeighor = 0;
363 
364  if (HasUnambiguousAngle(firstNeigh)) {
365  unambiguousNeighor = firstNeigh;
366  } else if (HasUnambiguousAngle(secondNeigh)) {
367  unambiguousNeighor = secondNeigh;
368  }
369 
370  if (unambiguousNeighor != 0) {
371  double similarAngle =
372  *fTangentAngles[unambiguousNeighor].begin();
373 
374  // check if angles lie close to each other, compare with first
375  if (GetDiffBetweenAngles(similarAngle,
376  *fTangentAngles[actualTubeId].begin())
377  > TMath::Pi() / 2) {
378 
379  // opposing angle is similar --> change
380  similarAngle = (similarAngle + TMath::Pi());
381 
382  if (similarAngle >= 2 * TMath::Pi())
383  similarAngle -= 2 * TMath::Pi();
384 
385  }
386 
387  double best = GetBestFittingAngle(similarAngle,
388  fTangentAngles[actualTubeId]);
389  KeepBestAngle(actualTubeId, best);
390 
391  if (fVerbose > 3)
392  cout << "unambiguous neighbor: "
393  << unambiguousNeighor << ", angle: "
394  << *fTangentAngles[unambiguousNeighor].begin()
395  * TMath::RadToDeg() << ", similar: "
396  << similarAngle * TMath::RadToDeg()
397  << ", keep:" << best * TMath::RadToDeg()
398  << endl;
399 
400  } else {
401  //try in next loop
402  tmpTubesToCheck.insert(actualTubeId);
403  if (fVerbose > 3)
404  cout << "try " << actualTubeId << " in next loop"
405  << endl;
406  }
407 
408  } else if (fVerbose > 3)
409  cout << "found no straight neighbors" << endl;
410 
411  }
412  }
413 
414  tubesToCheck.clear();
415  tubesToCheck = tmpTubesToCheck;
416 
417  if (fVerbose > 3) {
418  cout << "new fTangentAngles: " << endl;
419  PrintTangentAngles();
420  }
421  }
422 
423 //If the isochrone is very small --> track passes in the middle of the tube.
424 // Set fTangentAngles to -1.
425  for (vector<FairHit*>::iterator it = fHits.begin(); it != fHits.end();
426  ++it) {
427  PndSttHit* actualHit = (PndSttHit*) (*it);
428 
429 //no tangent angles found and isochrone is smaller than fMinIsochrone?
430  if (!fStrawMap->IsSkewedStraw(actualHit->GetTubeID())
431  && !HasUnambiguousAngle(actualHit->GetTubeID())
432  && actualHit->GetIsochrone() < fMinIsochrone) {
433  set<double> tmp;
434  tmp.insert(-1.0);
435  fTangentAngles[actualHit->GetTubeID()] = tmp;
436  }
437  }
438 
439  if (fVerbose > 2) {
440  cout << "Result:" << endl;
441  PrintTangentAngles();
442  }
443 
444  CorrectIsochrones();
445 
446 }
447 
448 void PndSttHitCorrector::PrintTangentAngles() {
449  cout << "Calculated tangent angles: " << endl;
450 
451  for (map<int, set<double> >::iterator mapIt = fTangentAngles.begin();
452  mapIt != fTangentAngles.end(); ++mapIt) {
453  cout << "angles for tube " << mapIt->first << ": ";
454  for (set<double>::iterator setIt = mapIt->second.begin();
455  setIt != mapIt->second.end(); ++setIt) {
456  if (*setIt == -1.0)
457  cout << *setIt << ", ";
458  else
459  cout << TMath::RadToDeg() * (*setIt) << ", ";
460  }
461  cout << endl;
462  }
463 }
464 
465 vector<double> PndSttHitCorrector::CalculateTangentAngles(PndSttHit* tube1,
466  PndSttHit* tube2) {
467 
468  if (fVerbose > 4)
469  cout << "PndSttHitCorrector::CalculateTangentAngles" << endl;
470 
471  Double_t radius1 = tube1->GetIsochrone();
472  Double_t radius2 = tube2->GetIsochrone();
473 
474  TVector3 origin1, origin2;
475  tube1->Position(origin1);
476  tube2->Position(origin2);
477 
478  TVector3 directionBetweenTubes = origin2 - origin1;
479 
480  if (fVerbose > 4) {
481  cout << "---- Input ----" << endl;
482  cout << "Origin1: ";
483  origin1.Print();
484  cout << " Origin2: ";
485  origin2.Print();
486  cout << " Dir: ";
487  directionBetweenTubes.Print();
488  cout << endl << "Radius1: " << radius1 << " Radius2: " << radius2
489  << endl;
490  }
491 
492  Double_t R, r;
493  Bool_t tube1Larger = kTRUE;
494  if (radius1 > radius2) {
495  R = radius1;
496  r = radius2;
497  tube1Larger = kTRUE;
498  } else {
499  R = radius2;
500  r = radius1;
501  tube1Larger = kFALSE;
502  }
503 
504  Double_t phi0 = directionBetweenTubes.Phi();
505  Double_t phi1 = TMath::ASin((R - r) / directionBetweenTubes.Mag());
506  Double_t phi2 = TMath::ASin((R + r) / directionBetweenTubes.Mag());
507 
508  if (fVerbose > 4) {
509  cout << "Phi0: " << TMath::RadToDeg() * phi0 << " Phi1: "
510  << TMath::RadToDeg() * phi1 << " Phi2: "
511  << TMath::RadToDeg() * phi2 << endl << endl;
512  }
513 
514  vector<double> results;
515  if (tube1Larger) {
516  Double_t angle1 = phi0 + phi1 - TMath::Pi() / 2;
517  if (angle1 < 0)
518  angle1 += TMath::Pi() * 2;
519  results.push_back(angle1);
520 
521  Double_t angle2 = phi0 - phi1 + TMath::Pi() / 2;
522  if (angle2 < 0)
523  angle2 += TMath::Pi() * 2;
524  results.push_back(angle2);
525 
526  Double_t angle3 = phi0 + phi2 - TMath::Pi() / 2;
527  if (angle3 < 0)
528  angle3 += TMath::Pi() * 2;
529  results.push_back(angle3);
530 
531  Double_t angle4 = phi0 - phi2 + TMath::Pi() / 2;
532  if (angle4 < 0)
533  angle4 += TMath::Pi() * 2;
534  results.push_back(angle4);
535 
536  if (fVerbose > 4) {
537  cout << "---- Output ----" << endl;
538  cout << "Phi0 + Phi1 - Pi/2 = " << TMath::RadToDeg() * angle1
539  << endl;
540  cout << "Phi0 - Phi1 + Pi/2 = " << TMath::RadToDeg() * angle2
541  << endl;
542  cout << "Phi0 + Phi2 - Pi/2 = " << TMath::RadToDeg() * angle3
543  << endl;
544  cout << "Phi0 - Phi2 + Pi/2 = " << TMath::RadToDeg() * angle4
545  << endl << endl;
546  }
547  } else {
548  Double_t angle1 = phi0 + phi1 + TMath::Pi() / 2;
549  if (angle1 < 0)
550  angle1 += TMath::Pi() * 2;
551  results.push_back(angle1);
552 
553  Double_t angle2 = phi0 - phi1 - TMath::Pi() / 2;
554  if (angle2 < 0)
555  angle2 += TMath::Pi() * 2;
556  results.push_back(angle2);
557 
558  Double_t angle3 = phi0 + phi2 - TMath::Pi() / 2;
559  if (angle3 < 0)
560  angle3 += TMath::Pi() * 2;
561  results.push_back(angle3);
562 
563  Double_t angle4 = phi0 - phi2 + TMath::Pi() / 2;
564  if (angle4 < 0)
565  angle4 += TMath::Pi() * 2;
566  results.push_back(angle4);
567 
568  if (fVerbose > 4) {
569  cout << "---- Output ----" << endl;
570  cout << "Phi0 + Phi1 + Pi/2 = " << TMath::RadToDeg() * angle1
571  << endl;
572  cout << "Phi0 - Phi1 - Pi/2 = " << TMath::RadToDeg() * angle2
573  << endl;
574  cout << "Phi0 + Phi2 - Pi/2 = " << TMath::RadToDeg() * angle3
575  << endl;
576  cout << "Phi0 - Phi2 + Pi/2 = " << TMath::RadToDeg() * angle4
577  << endl << endl;
578  }
579 
580  }
581 
582  return results;
583 
584 }
585 
586 int PndSttHitCorrector::GetMiddleHitNeighbor(int tubeID) {
587  if (fVerbose > 4)
588  cout << "PndSttHitCorrctor::GetMiddleHitNeighbor(" << tubeID << ")"
589  << endl;
590 
591  int middleTubeID = 0;
592 
593  if (fHitNeighbors[tubeID].size() == 3) {
594  double phi0 = fGeometryMap->GetAngleBetweenTubes(tubeID,
595  fHitNeighbors[tubeID].at(0));
596  double phi1 = fGeometryMap->GetAngleBetweenTubes(tubeID,
597  fHitNeighbors[tubeID].at(1));
598  double phi2 = fGeometryMap->GetAngleBetweenTubes(tubeID,
599  fHitNeighbors[tubeID].at(2));
600 
601 //Calculate the relative angles between the neighboring tubes
602  double diff01 = GetDiffBetweenAngles(phi0, phi1);
603  double diff12 = GetDiffBetweenAngles(phi1, phi2);
604  double diff02 = GetDiffBetweenAngles(phi0, phi2);
605 
606 //consider rounding errors
607  double delta = pow(10, -4);
608 
609 //two angles must be 60 degree
610  if ((diff01 - (TMath::Pi() / 3)) < delta
611  && (diff12 - (TMath::Pi() / 3)) < delta)
612  middleTubeID = fHitNeighbors[tubeID].at(1);
613  else if ((diff01 - (TMath::Pi() / 3)) < delta
614  && (diff02 - (TMath::Pi() / 3)) < delta)
615  middleTubeID = fHitNeighbors[tubeID].at(0);
616  else if ((diff12 - (TMath::Pi() / 3)) < delta
617  && (diff02 - (TMath::Pi() / 3)) < delta)
618  middleTubeID = fHitNeighbors[tubeID].at(2);
619 
620  if (fVerbose > 4) {
621  cout << "diff01: " << diff01 * TMath::RadToDeg() << ", diff12: "
622  << diff12 * TMath::RadToDeg() << ", diff02: "
623  << diff02 * TMath::RadToDeg() << endl;
624  cout << "middle-tube: " << middleTubeID << endl;
625  }
626 
627  }
628 
629  return middleTubeID;
630 }
631 
632 pair<int, int> PndSttHitCorrector::GetStraightNeighbors(int tubeID) {
633 
634  if (fVerbose > 4)
635  cout << "PndSttHitCorrector::GetStraightNeighbors(" << tubeID << ")"
636  << endl;
637 
638  pair<int, int> straightNeighbors;
639  straightNeighbors.first = 0;
640  straightNeighbors.second = 0;
641 
642  for (size_t i = 0; i < fHitNeighbors[tubeID].size() - 1; ++i) {
643 
644  for (size_t j = i + 1; j < fHitNeighbors[tubeID].size(); ++j) {
645 
646  if (fGeometryMap->InStraightLine(tubeID,
647  fHitNeighbors[tubeID].at(i), fHitNeighbors[tubeID].at(j))) {
648  straightNeighbors.first = fHitNeighbors[tubeID].at(i);
649  straightNeighbors.second = fHitNeighbors[tubeID].at(j);
650 
651  if (fVerbose > 4)
652  cout << "Found straight neighbors: "
653  << straightNeighbors.first << ", "
654  << straightNeighbors.second << endl;
655  break;
656  }
657  }
658  }
659 
660  return straightNeighbors;
661 }
662 
663 set<double> PndSttHitCorrector::GetBestCombinatedPhi(
664  const vector<vector<double> > &angles) {
665 
666  if (fVerbose > 3)
667  cout << "PndSttHitCorrector::GetBestCombinatedPhi()" << endl;
668 
669  set<double> ret;
670 
671 // Calculate differences and corresponding classification for angles
672  vector<vector<vector<double> > > differences;
673  vector<vector<double> > smallestDiff;
674  vector<vector<pair<int, int> > > smallestPairs;
675 
676  CalcDifferencesBetweenAngles(angles, differences, smallestDiff,
677  smallestPairs);
678 
679  vector<vector<double> > classification;
680  classification = CalcClassification(smallestDiff);
681 
682 // get the two best classification-values
683  int biggestIndex, secondBiggestIndex;
684  double biggest = 0, secondBiggest = 0;
685 
686  vector<double> tmp = classification.at(0);
687 
688  for (size_t j = 0; j < tmp.size(); ++j) {
689  if (tmp[j] > biggest) {
690  secondBiggest = biggest;
691  secondBiggestIndex = biggestIndex;
692 
693  biggestIndex = j;
694  biggest = tmp[j];
695 
696  } else if (tmp[j] > secondBiggest) {
697  secondBiggest = tmp[j];
698  secondBiggestIndex = j;
699 
700  }
701  }
702 
703 // get corresponding angles
704  pair<int, int> bestPair = smallestPairs.at(0)[biggestIndex];
705  double firstAngle = angles[0][bestPair.first];
706  double secondAngle = angles[1][bestPair.second];
707 
708  if (fVerbose > 3)
709  cout << "Best pair of angles: " << firstAngle * TMath::RadToDeg()
710  << ", " << secondAngle * TMath::RadToDeg() << endl;
711 
712  if (abs(firstAngle - secondAngle) < fMaxDiffBetweenAngles) {
713 
714  double bestPhi = (firstAngle + secondAngle) / 2;
715 
716  ret.insert(bestPhi);
717 
718  if (fVerbose > 3) {
719  cout << "Two best values of classification: " << biggest << ", "
720  << secondBiggest << endl;
721  cout << "Best Phi: " << TMath::RadToDeg() * bestPhi << endl;
722  }
723 
724 //check similarity of the two best (classification-)values
725  if (secondBiggest / biggest > 0.90) {
726  // found 2 appropriate angles
727  pair<int, int> secondBestPair =
728  smallestPairs.at(0)[secondBiggestIndex];
729  double secondBestPhi = (angles[0][secondBestPair.first]
730  + angles[1][secondBestPair.second]) / 2;
731  ret.insert(secondBestPhi);
732 
733  if (fVerbose > 3)
734  cout << "Second best Phi: " << TMath::RadToDeg() * secondBestPhi
735  << endl;
736  }
737  }
738 
739  return ret;
740 }
741 
742 void PndSttHitCorrector::CalcDifferencesBetweenAngles(
743  const vector<vector<double> > &angles,
744  vector<vector<vector<double> > > &retDifferences,
745  vector<vector<double> > &retSmallestDiff,
746  vector<vector<pair<int, int> > > &retPairsOfSmallest) {
747 
748  if (fVerbose > 4) {
749  cout << "PndSttHitCorrector::CalcDifferencesBetweenAngles()" << endl;
750 
751  cout << "Calculate differences between:";
752  for (vector<vector<double> >::const_iterator iterFirst = angles.begin();
753  iterFirst != angles.end(); ++iterFirst) {
754  vector<double> tmp = *iterFirst;
755  cout << "(";
756  for (size_t i = 0; i < tmp.size(); ++i) {
757  cout << TMath::RadToDeg() * tmp[i] << ",";
758  }
759  cout << ") ";
760  }
761  cout << endl;
762  }
763 // for each neighbor X
764  for (vector<vector<double> >::const_iterator iterFirst = angles.begin();
765  iterFirst != angles.end() - 1; ++iterFirst) {
766 
767 // calculate distances between angles of X and angles of all the other neighbors Y
768  for (vector<vector<double> >::const_iterator iterSecond = iterFirst + 1;
769  iterSecond != angles.end(); ++iterSecond) {
770 
771  // create a vector for the pair (X,Y) that will be filled with differences of the angles
772  vector<vector<double> > tmp2D;
773  retDifferences.push_back(tmp2D);
774 
775  // create vector for smallest differences and pairs
776  vector<double> tmp1D;
777  retSmallestDiff.push_back(tmp1D);
778  vector<pair<int, int> > tmpPair;
779  retPairsOfSmallest.push_back(tmpPair);
780 
781  vector<double> first = *iterFirst;
782  vector<double> second = *iterSecond;
783 
784  double currentDiff;
785  // for all angles k of X
786  for (size_t k = 0; k < first.size(); k++) {
787 
788  // create a vector for differences between k and all angles of Y
789  vector<double> tmp;
790  retDifferences.back().push_back(tmp);
791 
792  Double_t smallestDifference = 1000;
793  pair<int, int> indexOfSmallest;
794 
795  // calculate distances between k and all angles j of Y
796  for (size_t j = 0; j < second.size(); j++) {
797 
798  if (fVerbose > 4) {
799  cout << "Phi 0: " << k << " : "
800  << TMath::RadToDeg() * first[k] << " Phi 1: "
801  << j << " : " << TMath::RadToDeg() * second[j]
802  << endl;
803  }
804 
805  // save difference
806  currentDiff = GetDiffBetweenAngles(first[k], second[j]);
807  retDifferences.back().at(k).push_back(currentDiff);
808 
809  // get smallest difference
810  if (currentDiff < smallestDifference) {
811 
812  smallestDifference = currentDiff;
813  indexOfSmallest = make_pair(k, j);
814 
815  }
816  }
817 
818  // save smallest difference and associated indizes
819  retSmallestDiff.back().push_back(smallestDifference);
820  retPairsOfSmallest.back().push_back(indexOfSmallest);
821 
822  }
823  }
824  }
825 
826  if (fVerbose > 4) {
827 
828  cout << "Result of Calculation:" << endl;
829  int first = 0, second = 1;
830  for (size_t i = 0; i < retDifferences.size(); ++i) {
831 
832  cout << "Smallest Difference for vectors (" << first << ", "
833  << second << "): " << endl;
834 
835  for (size_t k = 0; k < retSmallestDiff.at(i).size(); ++k) {
836  cout << TMath::RadToDeg() * retSmallestDiff.at(i)[k]
837  << " for Phi [" << retPairsOfSmallest.at(i)[k].first
838  << ", " << retPairsOfSmallest.at(i)[k].second << "]"
839  << endl;
840  }
841 
842  if (second == (int)retDifferences.size() - 1) {
843  ++first;
844  second = first + 1;
845  } else {
846  ++second;
847  }
848  }
849 
850  cout << endl << endl;
851 
852  }
853 }
854 
855 vector<vector<double> > PndSttHitCorrector::CalcClassification(
856  const vector<vector<double> > &smallestDiff) {
857 
858  if (fVerbose > 4)
859  cout << "PndSttHitCorrector::CalcClassification(smallestDiff)" << endl;
860 
861  vector<vector<double> > classification;
862  vector<double> sumOfSmallest;
863 
864  for (size_t i = 0; i < smallestDiff.size(); ++i) {
865  sumOfSmallest.push_back(0);
866 
867 // calculate sum of smallest differences
868  for (size_t j = 0; j < smallestDiff[i].size(); ++j) {
869  sumOfSmallest[i] += smallestDiff[i][j];
870  }
871 
872  vector<double> tmp;
873  classification.push_back(tmp);
874 
875 // calculate classification
876  for (size_t j = 0; j < smallestDiff[i].size(); ++j) {
877  classification.back().push_back(
878  1 - smallestDiff[i][j] / sumOfSmallest[i]);
879  }
880  }
881 
882  if (fVerbose > 4) {
883  cout << "Calculated classification: " << endl;
884 
885  for (size_t i = 0; i < classification.size(); ++i) {
886  cout << "Data for " << i << ". vector of smallest differences: ";
887  for (size_t j = 0; j < classification[i].size(); ++j) {
888  cout << classification[i][j] << ", ";
889  }
890  cout << endl;
891  }
892 
893  cout << endl << endl;
894  }
895 
896  return classification;
897 }
898 
899 vector<vector<double> > PndSttHitCorrector::CalcClassification(
900  const vector<vector<vector<double> > > &differences) {
901 
902  vector<vector<double> > smallestValues;
903  vector<vector<pair<int, int> > > pairsOfSmallestValues;
904  vector<vector<double> > classification;
905 
906  if (fVerbose > 4)
907  cout << "PndSttHitCorrector::CalcClassification: classification"
908  << endl;
909 
910  vector<double> sumOfSmallestValues;
911 
912  int out = 0;
913 
914  for (vector<vector<vector<double> > >::const_iterator outerIter =
915  differences.begin(); outerIter != differences.end(); outerIter++) {
916  sumOfSmallestValues.push_back(0);
917  vector<double> tmp;
918  smallestValues.push_back(tmp);
919  classification.push_back(tmp);
920 
921  vector<pair<int, int> > tmp2;
922  pairsOfSmallestValues.push_back(tmp2);
923 
924  for (size_t i = 0; i < outerIter->size(); i++) {
925  double smallestValue = 10000;
926  pair<int, int> pairOfSmallest;
927 
928  for (size_t j = 0; j < outerIter->at(i).size(); j++) {
929  if (fVerbose > 4)
930  cout << "Pairs: " << i << "/" << j << " : "
931  << outerIter->at(i)[j] << endl;
932  if (outerIter->at(i)[j] < smallestValue) {
933  smallestValue = outerIter->at(i)[j];
934  pairOfSmallest = make_pair(i, j);
935  }
936  }
937 
938  cout << endl;
939 
940  smallestValues[out].push_back(smallestValue);
941  pairsOfSmallestValues[out].push_back(pairOfSmallest);
942  sumOfSmallestValues[out] += smallestValue;
943  }
944  for (size_t i = 0; i < smallestValues[out].size(); i++) {
945  classification[out].push_back(
946  1 - smallestValues[out][i] / sumOfSmallestValues[out]);
947  if (fVerbose > 4)
948  cout << "Pair: " << pairsOfSmallestValues[out][i].first << "/"
949  << pairsOfSmallestValues[out][i].second << " Value: "
950  << TMath::RadToDeg() * smallestValues[out][i]
951  << " Classification: " << classification[out][i]
952  << endl;
953  }
954  out++;
955 
956  }
957  return classification;
958 }
959 
960 void PndSttHitCorrector::CorrectIsochrones() {
961  if (fVerbose > 2)
962  cout << "void PndSttHitCorrector::CorrectIsochrones()" << endl;
963  int actualTubeId;
964  double bestPhi;
965 
966  for (map<int, set<double> >::iterator it = fTangentAngles.begin();
967  it != fTangentAngles.end(); ++it) {
968 
969  if (it->second.size() == 1) {
970  actualTubeId = it->first;
971  bestPhi = *it->second.begin();
972 
973  PndSttHit* actualHit =
974  (PndSttHit*) fHits[fMapTubeIdToHit[actualTubeId]];
975  TVector3 origin(actualHit->GetX(), actualHit->GetY(),
976  actualHit->GetZ());
977  TVector3 hitError(0.01, 0.01, 100);
978 
979  if (bestPhi != -1) {
980 
981  TVector3 correctedPosition(
982  origin.x()
983  + actualHit->GetIsochrone()
984  * TMath::Cos(bestPhi),
985  origin.y()
986  + actualHit->GetIsochrone()
987  * TMath::Sin(bestPhi), origin.z());
988 
989  if (fVerbose > 4) {
990  cout << "Corrected Position for " << actualTubeId << ": ";
991  correctedPosition.Print();
992  cout << endl;
993  }
994  fCorrectedIsochrones[fMapTubeIdToHit[actualTubeId]] =
995  new FairHit(-1, correctedPosition, hitError, -1);
996  } else {
997  //small isochrone --> take midpoint
998  fCorrectedIsochrones[fMapTubeIdToHit[actualTubeId]] =
999  new FairHit(-1, origin, hitError, -1);
1000  }
1001 
1002  fCorrectedIsochrones[fMapTubeIdToHit[actualTubeId]]->SetLink(
1003  FairLink("STTHit", fMapTubeIdToHit[actualTubeId]));
1004  } else {
1005  if (fVerbose > 2)
1006  cout << "[warning: ambiguous tangent angles for tube "
1007  << it->first << "]" << endl;
1008  }
1009  }
1010  int indexCounter = 0;
1011  for (map<int, FairHit*>::iterator it = fCorrectedIsochrones.begin();
1012  it != fCorrectedIsochrones.end(); ++it) {
1013  it->second->SetEntryNr(FairLink("CorrectedIsochrones", indexCounter++));
1014  }
1015 }
1016 
1017 void PndSttHitCorrector::KeepBestAngle(int tubeID, double angle) {
1018 
1019  if (fVerbose > 3)
1020  cout << "void PndSttHitCorrector::KeepBestAngle" << endl;
1021 
1022  double keep;
1023 
1024 //tangent angles were already calculated for this tube?
1025  if ((fTangentAngles.find(tubeID) != fTangentAngles.end())
1026  && (fTangentAngles.size() != 1)) {
1027 //--> keep the angle that is most similar to the angle of the unambiguous neighbor
1028  keep = GetBestFittingAngle(angle, fTangentAngles[tubeID]);
1029 
1030  } else {
1031 //no tangent angle calculated
1032 //--> store the angle of the unambiguous neighbor
1033  keep = angle;
1034  }
1035 
1036  fTangentAngles[tubeID].clear();
1037  fTangentAngles[tubeID].insert(keep);
1038 
1039  if (fVerbose > 3)
1040  cout << "Result for tube: " << tubeID << ", keep "
1041  << keep * TMath::RadToDeg() << endl;
1042 }
1043 
1044 double PndSttHitCorrector::GetDiffBetweenAngles(double angle1, double angle2) {
1045  double diff = abs(angle1 - angle2);
1046  if (diff > TMath::Pi())
1047  return 2 * TMath::Pi() - diff;
1048  else
1049  return diff;
1050 }
1051 
1052 double PndSttHitCorrector::GetAverageOfAngles(double angle1, double angle2) {
1053 
1054  double diff = abs(angle1 - angle2);
1055 
1056 //return the middle of the smaller angle between both angles
1057  if (diff > TMath::Pi()) {
1058  double tmp = max(angle1, angle2) + (2 * TMath::Pi() - diff) / 2;
1059  if (tmp >= 2 * TMath::Pi())
1060  tmp -= 2 * TMath::Pi();
1061  return tmp;
1062  } else
1063  return (angle1 + angle2) / 2;
1064 }
1065 
1066 double PndSttHitCorrector::GetBestFittingAngle(double angle,
1067  set<double> anglesToTest) {
1068 
1069  double minDiff = 1000;
1070  double tmpDiff;
1071  double best = -1;
1072 
1073  for (set<double>::iterator it = anglesToTest.begin();
1074  it != anglesToTest.end(); ++it) {
1075  tmpDiff = GetDiffBetweenAngles(angle, *it);
1076  if (tmpDiff < minDiff) {
1077  minDiff = tmpDiff;
1078  best = *it;
1079  }
1080  }
1081 
1082  return best;
1083 }
1084 
1085 bool PndSttHitCorrector::AreHitNeihbors(int tubeID1, int tubeID2) {
1086 
1087 //if both signal a hit
1088  if (fHitNeighbors.find(tubeID1) == fHitNeighbors.end()
1089  || fHitNeighbors.find(tubeID2) == fHitNeighbors.end()) {
1090  return false;
1091  }
1092 
1093  for (size_t i = 0; i < fHitNeighbors[tubeID1].size(); ++i) {
1094  if (fHitNeighbors[tubeID1].at(i) == tubeID2) {
1095  return true;
1096  }
1097  }
1098 
1099  return false;
1100 }
1101 
1102 bool PndSttHitCorrector::HasUnambiguousAngle(int tubeID) {
1103 
1104  if (fTangentAngles.find(tubeID) != fTangentAngles.end()
1105  && fTangentAngles[tubeID].size() == 1)
1106  return true;
1107  else
1108  return false;
1109 }
1110 
CAMath::ASin
static T ASin(const T &x)
fVerbose
int fVerbose
Definition: poormantracks.C:24
PndSttHitCorrector::CalcDifferencesBetweenAngles
void CalcDifferencesBetweenAngles(const std::vector< std::vector< double > > &angles, std::vector< std::vector< std::vector< double > > > &retDifferences, std::vector< std::vector< double > > &retSmallestDiff, std::vector< std::vector< std::pair< int, int > > > &retPairsOfSmallest)
Definition: PndSttHitCorrector.cxx:742
r
double r
Definition: RiemannTest.C:14
i
Int_t i
Definition: run_full.C:25
PndSttHitCorrector::HasUnambiguousAngle
bool HasUnambiguousAngle(int tubeID)
Definition: PndSttHitCorrector.cxx:1102
PndSttHitCorrector::GetDiffBetweenAngles
double GetDiffBetweenAngles(double angle1, double angle2)
Definition: PndSttHitCorrector.cxx:1044
map
PndTransMap * map
Definition: sim_emc_apd.C:99
PndSttHitCorrector::GetBestFittingAngle
double GetBestFittingAngle(double angle, std::set< double > anglesToTest)
Definition: PndSttHitCorrector.cxx:1066
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
PndSttHitCorrector::GetStraightNeighbors
std::pair< int, int > GetStraightNeighbors(int tubeID)
Definition: PndSttHitCorrector.cxx:632
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
PndSttHitCorrector::GetMiddleHitNeighbor
int GetMiddleHitNeighbor(int tubeID)
Definition: PndSttHitCorrector.cxx:586
PndSttHitCorrector::KeepBestAngle
void KeepBestAngle(int tubeID, double angle)
Definition: PndSttHitCorrector.cxx:1017
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndSttHitCorrector::CalculateTangentAngles
std::vector< double > CalculateTangentAngles(PndSttHit *tube1, PndSttHit *tube2)
Definition: PndSttHitCorrector.cxx:465
PndSttHitCorrector::CalcClassification
std::vector< std::vector< double > > CalcClassification(const std::vector< std::vector< double > > &smallestDiff)
Definition: PndSttHitCorrector.cxx:855
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
PndSttHit::GetTubeID
Int_t GetTubeID() const
Definition: PndSttHit.h:75
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
out
TFile * out
Definition: reco_muo.C:20
PndSttHitCorrector::AreHitNeihbors
bool AreHitNeihbors(int tubeID1, int tubeID2)
Definition: PndSttHitCorrector.cxx:1085
ClassImp
ClassImp(PndAnaContFact)
PndSttHit::Print
virtual void Print(const Option_t *opt=0) const
Definition: PndSttHit.h:52
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
PndSttHitCorrector::GetAverageOfAngles
double GetAverageOfAngles(double angle1, double angle2)
Definition: PndSttHitCorrector.cxx:1052
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
R
Double_t R
Definition: checkhelixhit.C:61
PndSttHitCorrector::GetBestCombinatedPhi
std::set< double > GetBestCombinatedPhi(const std::vector< std::vector< double > > &angles)
Definition: PndSttHitCorrector.cxx:663
Generated on Wed Apr 3 2019 10:02:57 for FairRoot/PandaRoot by   doxygen 1.8.5