FairRoot/PandaRoot
PndFtsLineApproximator.cxx
Go to the documentation of this file.
1 /*
2  * PndFtsLineApproximator.cpp
3  *
4  * Created on: Jun 2, 2016
5  * Author: kibellus
6  */
7 
8 #include <PndFtsLineApproximator.h>
9 
10 Double_t PndFtsLineApproximator::compDist = 0.2;
11 Double_t PndFtsLineApproximator::compAngle = 0.2;
12 
13 PndFtsLineApproximator::~PndFtsLineApproximator() {
14  // TODO Auto-generated destructor stub
15 }
16 
17 PndFtsExpandedTrackCand PndFtsLineApproximator::createExpandedTrackCand(PndTrackCand &cand){
18  FairRootManager *ioman = FairRootManager::Instance();
19  Int_t branchID = ioman->GetBranchId("FTSHit"); //right hits
20  vector<PndFtsHit*> hits;
21  set<FairLink> linksToHits = cand.GetLinksWithType(branchID).GetLinks();
22 
23  for (std::set<FairLink>::iterator it = linksToHits.begin();it!=linksToHits.end();it++) {
24  FairLink link = *it;
25  PndFtsHit* hit = (PndFtsHit*) ioman->GetCloneOfLinkData(link);
26  transform(kTRUE,hit);
27  hit->SetY(0); //project to y=0
28  hits.push_back(hit);
29  }
30  fComparator.setZValue(hits[0]->GetZ());
31 
32  //track with only one hit
33  if(hits.size()==1){
34  vector<PndLineApproximation> approx;//empty vector
35  return PndFtsExpandedTrackCand(cand,hits,approx);
36  } else if(hits.size()==2){
37  //two hits => 4 lines
38  vector<PndLineApproximation> approx = approxLines(hits[0],hits[1]);
39  return PndFtsExpandedTrackCand(cand,hits,approx);
40  } else {
41  //more than 2 hits => search the best line
42  vector<vector<PndLineApproximation>> lines;
43  //create all combinations of possible lines
44  for(size_t i=0;i<hits.size();i++){
45  for(size_t j=i+1;j<hits.size();j++){
46  vector<PndLineApproximation> lines2 = approxLines(hits[i],hits[j]);
47  for(size_t k=0;k<lines2.size();k++)
48  addLine(lines,lines2[k]);
49 
50  }
51  }
52  //find the max
53  size_t max= 0;
54  for(size_t i=0;i<lines.size();i++){
55  if(lines[i].size()>max)
56  max = lines[i].size();
57  }
58  vector<PndLineApproximation> result;
59  //find all approximations with size=max
60  for(size_t i=0;i<lines.size();i++){
61  if(lines[i].size()==max){
62  PndLineApproximation a = lines[i][0];
63  for(size_t j=1;j<lines[i].size();j++)
64  a = a.newApproximation(lines[i][j]);
65  a.getLine().setRating(1);
66  result.push_back(a);
67  }
68  }
69  return PndFtsExpandedTrackCand(cand,hits,result);
70  }
71 }
72 
73 vector<PndLineApproximation> PndFtsLineApproximator::approxLines(PndFtsHit* h1 ,PndFtsHit* h2){
74  vector<PndLineApproximation> result;
75  vector<PndLineApproximation> r1 = approxInnerLine(h1,h2);
76  result.push_back(r1[0]);
77  result.push_back(r1[1]);
78  vector<PndLineApproximation> r2 = approxOuterLine(h1,h2);
79  result.push_back(r2[0]);
80  result.push_back(r2[1]);
81  return result;
82 }
83 
84 vector<PndLineApproximation> PndFtsLineApproximator::approxInnerLine(PndFtsHit* h1 ,PndFtsHit* h2){
85  vector<PndLineApproximation> result;
86  Double_t r0 = h1->GetIsochrone();
87  Double_t r1 = h2->GetIsochrone();
88  Double_t a = h1->GetX();
89  Double_t b = h1->GetZ();
90  Double_t c = h2->GetX();
91  Double_t d = h2->GetZ();
92 
93  Double_t xp = (c*r0-a*r1) / (r0 - r1);
94  Double_t yp = (d*r0-b*r1) / (r0 - r1);
95 
96  Double_t xt1Z = r0*r0*(xp-a) + r0*(yp-b)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
97  Double_t xt1N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
98  Double_t xt1 = xt1Z/xt1N+a;
99 
100  Double_t yt1Z = r0*r0*(yp-b) + r0*(xp-a)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
101  Double_t yt1N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
102  Double_t yt1 = yt1Z/yt1N+b;
103 
104  Double_t xt2Z = r0*r0*(xp-a) - r0*(yp-b)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
105  Double_t xt2N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
106  Double_t xt2 = xt2Z/xt2N+a;
107 
108  Double_t yt2Z = r0*r0*(yp-b) - r0*(xp-a)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
109  Double_t yt2N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
110  Double_t yt2 = yt2Z/yt2N+b;
111 
112  Double_t xt3Z = r1*r1*(xp-c) + r1*(yp-d)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
113  Double_t xt3N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
114  Double_t xt3 = xt3Z/xt3N+c;
115 
116  Double_t yt3Z = r1*r1*(yp-d) + r1*(xp-c)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
117  Double_t yt3N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
118  Double_t yt3 = yt3Z/yt3N+d;
119 
120  Double_t xt4Z = r1*r1*(xp-c) - r1*(yp-d)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
121  Double_t xt4N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
122  Double_t xt4 = xt4Z/xt4N+c;
123 
124  Double_t yt4Z = r1*r1*(yp-d) - r1*(xp-c)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
125  Double_t yt4N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
126  Double_t yt4 = yt4Z/yt4N+d;
127 
128  //outer lines
129  Double_t s = TMath::Abs((b-yt1)*(yp-yt1)/((xt1-a)*(xt1-xp)));
130  Double_t p1x = xt1;
131  Double_t p1y;
132  Double_t p2x = xt2;
133  Double_t p2y;
134  if(TMath::Abs(s-1)<0.0000000001){
135  p1y = yt1;
136  p2y = yt2;
137  } else {
138  p1y = yt2;
139  p2y = yt1;
140  }
141 
142  s = TMath::Abs((d-yt3)*(yp-yt3)/((xt3-c)*(xt3-xp)));
143 
144  Double_t p3x = xt3;
145  Double_t p3y;
146  Double_t p4x = xt4;
147  Double_t p4y;
148  if(TMath::Abs(s-1)<0.0000000001){
149  p3y = yt3;
150  p4y = yt4;
151  } else {
152  p3y = yt4;
153  p4y = yt3;
154  }
155 
156  //create the corrected hits
157  PndFtsHit* newHit1 = copyHitWithNewPosition(h1,p1x,p1y);
158  PndFtsHit* newHit2 = copyHitWithNewPosition(h1,p2x,p2y);
159  PndFtsHit* newHit3 = copyHitWithNewPosition(h2,p3x,p3y);
160  PndFtsHit* newHit4 = copyHitWithNewPosition(h2,p4x,p4y);
161 
162  //transform to normal coordinate-system
163  transform(kFALSE,newHit1);
164  transform(kFALSE,newHit2);
165  transform(kFALSE,newHit3);
166  transform(kFALSE,newHit4);
167 
168  //create 1st line
169  vector<PndFtsHit*> hits1;
170  hits1.push_back(newHit1);
171  hits1.push_back(newHit3);
172  PndLine l1(newHit1,newHit3);
173  result.push_back(PndLineApproximation(l1,hits1));
174 
175  //create 2nd line
176  vector<PndFtsHit*> hits2;
177  hits2.push_back(newHit2);
178  hits2.push_back(newHit4);
179  PndLine l2(newHit2,newHit4);
180  result.push_back(PndLineApproximation(l2,hits2));
181 
182  return result;
183 }
184 
185 vector<PndLineApproximation> PndFtsLineApproximator::approxOuterLine(PndFtsHit* h1 ,PndFtsHit* h2){
186  vector<PndLineApproximation> result;
187  Double_t r0 = h1->GetIsochrone();
188  Double_t r1 = h2->GetIsochrone();
189 
190  Double_t a = h1->GetX();
191  Double_t b = h1->GetZ();
192  Double_t c = h2->GetX();
193  Double_t d = h2->GetZ();
194 
195  Double_t xp = (c*r0+a*r1) / (r0 + r1);
196  Double_t yp = (d*r0+b*r1) / (r0 + r1);
197 
198  Double_t xt1Z = r0*r0*(xp-a) + r0*(yp-b)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
199  Double_t xt1N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
200  Double_t xt1 = xt1Z/xt1N+a;
201 
202  Double_t yt1Z = r0*r0*(yp-b) + r0*(xp-a)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
203  Double_t yt1N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
204  Double_t yt1 = yt1Z/yt1N+b;
205 
206  Double_t xt2Z = r0*r0*(xp-a) - r0*(yp-b)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
207  Double_t xt2N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
208  Double_t xt2 = xt2Z/xt2N+a;
209 
210  Double_t yt2Z = r0*r0*(yp-b) - r0*(xp-a)*TMath::Sqrt((xp-a)*(xp-a)+(yp-b)*(yp-b)-r0*r0);
211  Double_t yt2N = (xp-a)*(xp-a)+(yp-b)*(yp-b);
212  Double_t yt2 = yt2Z/yt2N+b;
213 
214  Double_t xt3Z = r1*r1*(xp-c) + r1*(yp-d)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
215  Double_t xt3N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
216  Double_t xt3 = xt3Z/xt3N+c;
217 
218  Double_t yt3Z = r1*r1*(yp-d) + r1*(xp-c)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
219  Double_t yt3N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
220  Double_t yt3 = yt3Z/yt3N+d;
221 
222  Double_t xt4Z = r1*r1*(xp-c) - r1*(yp-d)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
223  Double_t xt4N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
224  Double_t xt4 = xt4Z/xt4N+c;
225 
226  Double_t yt4Z = r1*r1*(yp-d) - r1*(xp-c)*TMath::Sqrt((xp-c)*(xp-c)+(yp-d)*(yp-d)-r1*r1);
227  Double_t yt4N = (xp-c)*(xp-c)+(yp-d)*(yp-d);
228  Double_t yt4 = yt4Z/yt4N+d;
229 
230  //outer lines
231  Double_t s = TMath::Abs((b-yt1)*(yp-yt1)/((xt1-a)*(xt1-xp)));
232  Double_t p1x = xt1;
233  Double_t p1y;
234  Double_t p2x = xt2;
235  Double_t p2y;
236  if(TMath::Abs(s-1)<0.0000000001){
237  p1y = yt1;
238  p2y = yt2;
239  } else {
240  p1y = yt2;
241  p2y = yt1;
242  }
243 
244  s = TMath::Abs((d-yt3)*(yp-yt3)/((xt3-c)*(xt3-xp)));
245 
246  Double_t p3x = xt3;
247  Double_t p3y;
248  Double_t p4x = xt4;
249  Double_t p4y;
250  if(TMath::Abs(s-1)<0.0000000001){
251  p3y = yt3;
252  p4y = yt4;
253  } else {
254  p3y = yt4;
255  p4y = yt3;
256  }
257 
258  //create the corrected hits
259  PndFtsHit* newHit1 = copyHitWithNewPosition(h1,p1x,p1y);
260  PndFtsHit* newHit2 = copyHitWithNewPosition(h1,p2x,p2y);
261  PndFtsHit* newHit3 = copyHitWithNewPosition(h2,p3x,p3y);
262  PndFtsHit* newHit4 = copyHitWithNewPosition(h2,p4x,p4y);
263 
264  //transform to normal coordinate-system
265  transform(kFALSE,newHit1);
266  transform(kFALSE,newHit2);
267  transform(kFALSE,newHit3);
268  transform(kFALSE,newHit4);
269 
270  //create 1st line
271  vector<PndFtsHit*> hits1;
272 
273  hits1.push_back(newHit1);
274  hits1.push_back(newHit3);
275  PndLine l1(newHit1,newHit3);
276  //TLine l1(newHit1->GetX(),newHit1->GetZ(),newHit3->GetX(),newHit3->GetZ());
277  result.push_back(PndLineApproximation(l1,hits1));
278 
279  //create 2nd line
280  vector<PndFtsHit*> hits2;
281  hits2.push_back(newHit2);
282  hits2.push_back(newHit4);
283  //TLine l2(newHit2->GetX(),newHit2->GetZ(),newHit4->GetX(),newHit4->GetZ());
284  PndLine l2(newHit2,newHit4);
285  result.push_back(PndLineApproximation(l2,hits2));
286 
287  return result;
288 }
289 
290 void PndFtsLineApproximator::addLine(vector<vector<PndLineApproximation>> &lines, PndLineApproximation &l){
291  Bool_t used = kFALSE;
292  for(size_t i=0;i<lines.size();i++){
293  PndLine l2 = lines[i][0].getLine();
294  if(fComparator.areEqual2D(l.getLine(),l2)){
295  lines[i].push_back(l);
296  used=kTRUE;
297  }
298  }
299  if(used)return;
300  vector<PndLineApproximation> newLine;
301  newLine.push_back(l);
302  lines.push_back(newLine);
303 }
304 /*
305 PndLineApproximation PndFtsLineApproximator::createLine(vector<PndLineApproximation> lines){
306  vector<PndFtsHit*> correctedHits;
307  for(int i=0;i<lines.size();i++){
308  vector<PndFtsHit*> hitsOfLine = lines[i].getHits();
309  for(int j=0;j<hitsOfLine.size();j++){
310  Bool_t exist = kFALSE;
311  for(int k=0;k<correctedHits.size();k++){
312  if(correctedHits[k]->GetTubeID()==hitsOfLine[j]->GetTubeID()){
313  exist = kTRUE;
314  break;
315  }
316  }
317  if(!exist)
318  correctedHits.push_back(hitsOfLine[j]);
319  }
320  }
321  PndLineApproximation newApprox(lines[0].getLine(),correctedHits);
322  return newApprox;
323 }*/
324 
325 PndFtsHit* PndFtsLineApproximator::copyHitWithNewPosition(PndFtsHit *h, Double_t x, Double_t z){
326  PndFtsHit* h2 = new PndFtsHit();
327  h2->SetDetectorID(h->GetDetectorID());
328  h2->SetTubeID(h->GetTubeID());
329  h2->SetChamberID(h->GetChamberID());
330  h2->SetLayerID(h->GetLayerID());
331  h2->SetXYZ(x,h->GetY(),z);
332  h2->SetIsochrone(h->GetIsochrone());
333  h2->SetIsochroneError(h->GetIsochroneError());
334  h2->SetRefIndex(h->GetRefIndex());
335  h2->SetEntryNr(h->GetEntryNr());
336  return h2;
337 }
338 
339 TMatrix PndFtsLineApproximator::getRotationMatrix(Double_t angle){
340  TMatrix a(3,3);
341  a[0][0] = TMath::Cos(angle);
342  a[0][1] = -TMath::Sin(angle);
343  a[0][2] = 0;
344  a[1][0] = TMath::Sin(angle);
345  a[1][1] = TMath::Cos(angle);
346  a[1][2] = 0;
347  a[2][0] = 0;
348  a[2][1] = 0;
349  a[2][2] = 1;
350  return a;
351 }
352 
353 void PndFtsLineApproximator::transform(Bool_t transToNewSystem, PndFtsHit* hit){
354  Double_t angle = 0;
355  if(((hit->GetLayerID() - 1) / 2)%4==1) angle = -5;
356  if(((hit->GetLayerID() - 1) / 2)%4==2) angle = 5;
357  if(angle==0) return;
358  if(!transToNewSystem) angle*=-1;
359  TMatrix rotMat = getRotationMatrix(angle*TMath::DegToRad());
360  TMatrix hitMat(3,1);
361  hitMat[0][0] = hit->GetX();
362  hitMat[1][0] = hit->GetY();
363  hitMat[2][0] = hit->GetZ();
364  TMatrix newPos = rotMat*hitMat;
365  hit->SetXYZ(newPos[0][0],newPos[1][0],newPos[2][0]);
366 }
PndFtsHit::SetIsochrone
void SetIsochrone(Double_t isochrone)
Definition: PndFtsHit.h:64
PndFtsLineApproximator::addLine
void addLine(vector< vector< PndLineApproximation >> &lines, PndLineApproximation &l)
Definition: PndFtsLineApproximator.cxx:290
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
i
Int_t i
Definition: run_full.C:25
PndFtsLineApproximator::approxOuterLine
vector< PndLineApproximation > approxOuterLine(PndFtsHit *h1, PndFtsHit *h2)
Definition: PndFtsLineApproximator.cxx:185
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndLineApproximation::newApproximation
PndLineApproximation newApproximation(PndLineApproximation &approx2)
Definition: PndLineApproximation.cxx:54
PndFtsHit::SetChamberID
void SetChamberID(Int_t chamberid)
Definition: PndFtsHit.h:71
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
s
TLorentzVector s
Definition: Pnd2DStar.C:50
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
PndFtsLineApproximator::copyHitWithNewPosition
PndFtsHit * copyHitWithNewPosition(PndFtsHit *h, Double_t x, Double_t z)
Definition: PndFtsLineApproximator.cxx:325
PndFtsHit::SetTubeID
void SetTubeID(Int_t tubeid)
Definition: PndFtsHit.h:69
PndFtsLineApproximator::fComparator
PndFtsLineComparator fComparator
Definition: PndFtsLineApproximator.h:41
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
r1
double r1
Definition: RiemannTrackFinderTest.C:14
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndFtsHit::GetTubeID
Int_t GetTubeID() const
Definition: PndFtsHit.h:70
PndFtsHit::GetLayerID
Int_t GetLayerID() const
Definition: PndFtsHit.h:74
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
PndFtsLineComparator::setZValue
void setZValue(Double_t z)
Definition: PndFtsLineComparator.h:28
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndFtsLineApproximator::createExpandedTrackCand
PndFtsExpandedTrackCand createExpandedTrackCand(PndTrackCand &cand)
Definition: PndFtsLineApproximator.cxx:17
a
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndLine::setRating
void setRating(Int_t r)
Definition: PndLine.h:44
PndFtsLineApproximator::getRotationMatrix
TMatrix getRotationMatrix(Double_t angle)
Definition: PndFtsLineApproximator.cxx:339
z
Double_t z
Definition: createRootFscGeometryFile.C:146
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
PndFtsHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndFtsHit.h:57
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndFtsLineComparator::areEqual2D
Bool_t areEqual2D(PndLine l1, PndLine l2)
Definition: PndFtsLineComparator.cxx:29
PndFtsHit::GetIsochroneError
Double_t GetIsochroneError() const
Definition: PndFtsHit.h:58
PndFtsHit::SetLayerID
void SetLayerID(Int_t layerid)
Definition: PndFtsHit.h:73
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
PndFtsLineApproximator::transform
void transform(Bool_t transToNewSystem, PndFtsHit *hit)
Definition: PndFtsLineApproximator.cxx:353
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
r2
double r2
Definition: RiemannTrackFinderTest.C:19
PndFtsLineApproximator::approxInnerLine
vector< PndLineApproximation > approxInnerLine(PndFtsHit *h1, PndFtsHit *h2)
Definition: PndFtsLineApproximator.cxx:84
PndFtsLineApproximator::approxLines
vector< PndLineApproximation > approxLines(PndFtsHit *h1, PndFtsHit *h2)
Definition: PndFtsLineApproximator.cxx:73
PndFtsHit::SetIsochroneError
void SetIsochroneError(Double_t isochroneError)
Definition: PndFtsHit.h:65
PndFtsHit::GetChamberID
Int_t GetChamberID() const
Definition: PndFtsHit.h:72
Generated on Wed Apr 3 2019 08:43:56 for FairRoot/PandaRoot by   doxygen 1.8.5