FairRoot/PandaRoot
PndTrkClean.cxx
Go to the documentation of this file.
1 //
2 // PndTrkClean.cxx
3 //
4 //
5 //
6 //
7 // authors: Lia Lavezzi - INFN Pavia (2012)
8 //
9 
10 #include "PndTrkClean.h"
11 
12 #include "PndTrkTools.h"
13 #include "PndTrkLegendreCluster.h"
14 #include "PndSttTube.h"
15 
16 #include "TArc.h"
17 #include <iostream>
18 
19 
20 using namespace std;
21 
22 PndTrkClean::PndTrkClean() : fTubeArray(NULL) {
23  fGeoH = PndGeoHandling::Instance();
24 }
25 
26 PndTrkClean::PndTrkClean(TClonesArray *tubearray) {
27  fGeoH = PndGeoHandling::Instance();
28  fTubeArray = tubearray;
29 }
30 
31 PndTrkClean::~PndTrkClean() {}
32 
33 // /**
34 // LAY NAME
35 
36 // barrel, layers from innermost to outermost:
37 // 1 PixeloBlo1
38 // 6 PixeloBlo2
39 // 15 StripoBl3o(Silicon)
40 // 17 StripoBl4o(Silicon)
41 
42 // forward wheels, inner, increasing z:
43 // left
44 // 2 PixeloSdkoco(Silicon)_1
45 // 3 PixeloSdkoco(Silicon)_2
46 // 7 PixeloLdkoco(Silicon)_1
47 // 8 PixeloLdkoco(Silicon)_2
48 // 9 PixeloLdkoco(Silicon)_3
49 // 10 PixeloLdkoco(Silicon)_4
50 
51 // right
52 // 4 PixeloSdkoco(Silicon)_3
53 // 5 PixeloSdkoco(Silicon)_4
54 // 11 PixeloLdkoco(Silicon)_5
55 // 12 PixeloLdkoco(Silicon)_6
56 // 13 PixeloLdkoco(Silicon)_7
57 // 14 PixeloLdkoco(Silicon)_8
58 
59 // forward wheels, last two, outer:
60 // left
61 // 16 Fwdo(Silicon)_1
62 // right
63 // 20 StripoLdko5-6oTrapSo(Silicon)_1
64 
65 
66 // **/
67 // int PndTrkClean::FindMvdLayer(int sensorID) {
68 // TString geoPath;
69 // int Layer = 0;
70 
71 // bool flag = true;
72 
73 // geoPath=fGeoH->GetPath(sensorID);
74 // cout << "path " << geoPath << endl;
75 // TVector3 o, u, v;
76 // fGeoH->GetOUVPath(geoPath, o, u, v);
77 // o.Print();
78 // u.Print();
79 // v.Print();
80 
81 
82 // if (flag && geoPath.Contains("PixeloBlo1")){Layer=1;flag=false;}
83 // if (flag && geoPath.Contains("PixeloSdko(Silicon)_1")){Layer=2;flag=false;} //naming after MVD2.2
84 // if (flag && geoPath.Contains("PixeloSdkoco(Silicon)_1")){Layer=2;flag=false;}
85 // if (flag && geoPath.Contains("PixeloSdko(Silicon)_2")){Layer=3;flag=false;} //naming after MVD2.2
86 // if (flag && geoPath.Contains("PixeloSdkoco(Silicon)_2")){Layer=3;flag=false;}
87 // if (flag && geoPath.Contains("PixeloSdko(Silicon)_3")){Layer=4;flag=false;} //naming after MVD2.2
88 // if (flag && geoPath.Contains("PixeloSdkoco(Silicon)_3")){Layer=4;flag=false;}
89 // if (flag && geoPath.Contains("PixeloSdko(Silicon)_4")){Layer=5;flag=false;} //naming after MVD2.2
90 // if (flag && geoPath.Contains("PixeloSdkoco(Silicon)_4")){Layer=5;flag=false;}
91 // if (flag && geoPath.Contains("PixeloBlo2")){Layer=6;flag=false;}
92 // if (flag && geoPath.Contains("PixeloLdkoio(Silicon)_1")){Layer=7;flag=false;} //naming after MVD2.2
93 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_1")){Layer=7;flag=false;}
94 // if (flag && geoPath.Contains("PixeloLdkoiio(Silicon)_1")){Layer=8;flag=false;} //naming after MVD2.2
95 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_2")){Layer=8;flag=false;}
96 // if (flag && geoPath.Contains("PixeloLdkoiiio(Silicon)_1")){Layer=9;flag=false;} //naming after MVD2.2
97 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_3")){Layer=9;flag=false;}
98 // if (flag && geoPath.Contains("PixeloLdkoiiio(Silicon)_2")){Layer=10;flag=false;} //naming after MVD2.2
99 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_4")){Layer=10;flag=false;}
100 // if (flag && geoPath.Contains("PixeloLdkoiio(Silicon)_2")){Layer=11;flag=false;} //naming after MVD2.2
101 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_5")){Layer=11;flag=false;}
102 // if (flag && geoPath.Contains("PixeloLdkoio(Silicon)_2")){Layer=12;flag=false;} //naming after MVD2.2
103 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_6")){Layer=12;flag=false;}
104 // if (flag && geoPath.Contains("PixeloLdkoiiio(Silicon)_4")){Layer=13;flag=false;} //naming after MVD2.2
105 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_7")){Layer=13;flag=false;}
106 // if (flag && geoPath.Contains("PixeloLdkoiiio(Silicon)_3")){Layer=14;flag=false;} //naming after MVD2.2
107 // if (flag && geoPath.Contains("PixeloLdkoco(Silicon)_8")){Layer=14;flag=false;}
108 
109 // if (flag && geoPath.Contains("StripoBl3o(Silicon)")){Layer=15;flag=false;}
110 // if (flag && geoPath.Contains("Fwdo(Silicon)_1")){Layer=16;flag=false;}
111 // if (flag && geoPath.Contains("StripoBl4o(Silicon)")){Layer=17;flag=false;}
112 // //if (flag && geoPath.Contains("Fwdo(Silicon)_2")){Layer=10;flag=false;}
113 
114 // if (flag && geoPath.Contains("StripoLdkoTrapSoRingAoSilicon_1")){Layer=18;flag=false;}
115 // if (flag && geoPath.Contains("StripoLdkoTrapSoRingAoSilicon_2")){Layer=19;flag=false;}
116 // if (flag && geoPath.Contains("StripoLdkoTrapSoRingBoSilicon_1")){Layer=18;flag=false;}
117 // if (flag && geoPath.Contains("StripoLdkoTrapSoRingBoSilicon_2")){Layer=19;flag=false;}
118 // if (flag && geoPath.Contains("StripoLdko5-6oTrapSo(Silicon)_1")){Layer=20;flag=false;}
119 // return Layer;
120 // }
121 
165 int PndTrkClean::FindMvdLayer(int sensorID) {
166  TString geoPath;
167  int Layer = 0;
168 
169  geoPath=fGeoH->GetPath(sensorID);
170  // cout << "path " << geoPath << endl;
171  TVector3 o, u, v;
172  fGeoH->GetOUVPath(geoPath, o, u, v);
173 // o.Print();
174 // u.Print();
175 // v.Print();
176 
177  // barrel
178  // left
179  if (geoPath.Contains("PixeloBlo1")) {
180  Layer=1;
181  // right
182  if(o.X() < 0) Layer += 1;
183  }
184  else if (geoPath.Contains("PixeloBlo2")) {
185  Layer=3;
186  // right
187  if(o.X() < 0) Layer += 1;
188  }
189  else if (geoPath.Contains("StripoBl3o(Silicon)")) {
190  Layer=5;
191  // right
192  if(o.X() < 0) Layer += 1;
193  }
194  else if (geoPath.Contains("StripoBl4o(Silicon)")) {
195  Layer=7;
196  // right
197  if(o.X() < 0) Layer += 1;
198  }
199  // inner wheels
200  // left
201  else if (geoPath.Contains("PixeloSdkoco(Silicon)_1")){Layer=9;}
202  else if (geoPath.Contains("PixeloSdkoco(Silicon)_2")){Layer=11;}
203  // right0
204  else if (geoPath.Contains("PixeloSdkoco(Silicon)_3")){Layer=10;}
205  else if (geoPath.Contains("PixeloSdkoco(Silicon)_4")){Layer=12;}
206  //
207  // left
208  else if (geoPath.Contains("PixeloLdkoco(Silicon)_1")){Layer=13;}
209  else if (geoPath.Contains("PixeloLdkoco(Silicon)_2")){Layer=15;}
210  else if (geoPath.Contains("PixeloLdkoco(Silicon)_3")){Layer=17;}
211  else if (geoPath.Contains("PixeloLdkoco(Silicon)_4")){Layer=19;}
212  // right
213  else if (geoPath.Contains("PixeloLdkoco(Silicon)_5")){Layer=14;}
214  else if (geoPath.Contains("PixeloLdkoco(Silicon)_6")){Layer=16;}
215  else if (geoPath.Contains("PixeloLdkoco(Silicon)_7")){Layer=18;}
216  else if (geoPath.Contains("PixeloLdkoco(Silicon)_8")){Layer=20;}
217 
218  // last two wheels
219  else if (geoPath.Contains("Fwdo(Silicon)_1")){ // left
220  // 1st
221  Layer=21;
222  // 2nd
223  if(o.Z() > 20) Layer += 2;
224  }
225  else if (geoPath.Contains("Fwdo(Silicon)_2")){ // right
226  // 1st
227  Layer=22;
228  // 2nd
229  if(o.Z() > 20) Layer += 2;
230  }
231  else cout << "I have no idea" << endl;
232 
233 
234  return Layer;
235 }
236 
237 
238 Bool_t PndTrkClean::CheckPairOfHits(PndTrkHit *hit1, PndTrkHit *hit2) {
239 
240  // cout << "hit1 " << hit1 << endl;
241  // cout << "hit2 " << hit2 << endl;
242  // both mvd
243  if(hit1->IsMvd() && hit2->IsMvd()) {
244  // cout << "BOTH MVD " << hit1->GetDistance(hit2) << endl;
245  if(hit1->GetDistance(hit2) < 7) return kTRUE;
246  }
247 
248  // one mvd - one stt
249  else if((hit1->IsMvd() && hit2->IsStt()) || (hit1->IsStt() && hit2->IsMvd())) {
250  // cout << "ONE MVD/ONE STT " << hit1->GetXYDistance(hit2) << endl;
251 
252 
253  // if the stt one is skew -> WRONG for sure!
254  if((hit2->IsStt() && hit2->IsSttSkew()) || (hit1->IsStt() && hit1->IsSttSkew())) return kFALSE;
255  // if parallel:
256  else if(hit1->GetXYDistance(hit2) < 8) return kTRUE;
257  }
258 
259  // both stt
260  else if(hit1->IsStt() && hit2->IsStt()) {
261  // one parallel one skew
262  if((hit1->IsSttParallel() && hit2->IsSttSkew()) || (hit2->IsSttParallel() && hit1->IsSttSkew())) {
263 
264 // cout << "PARALLEL/SKEW" << endl;
265 // cout << hit1->IsSttSkew() << " " << hit1->GetHitID() << " " << hit1->GetDetectorID() << endl;
266 // cout << hit2->IsSttSkew() << " " << hit2->GetHitID() << " " << hit2->GetDetectorID() << endl;
267 // cout << hit1->GetXYDistance(hit2) << endl;
268 
269  if(hit1->GetXYDistance(hit2) < 3) return kTRUE;
270  }
271  // both parallel
272  else if(hit1->IsSttParallel() && hit2->IsSttParallel()) {
273  PndSttTube *tube1 = (PndSttTube*) fTubeArray->At(hit1->GetTubeID());
274  int layerID1 = tube1->GetLayerID();
275  PndSttTube *tube2 = (PndSttTube*) fTubeArray->At(hit2->GetTubeID());
276  int layerID2 = tube2->GetLayerID();
277  //
278  // cout << "layers " << layerID1 << " " << layerID2 << endl;
279 
280  // same layer
281  if(layerID1 == tube2->GetLayerID()) return kTRUE;
282  // both inner/outer parallel layers
283  else if((layerID1 < 9 && layerID2 < 9) || (layerID1 > 8 && layerID2 > 8)) {
284  // adjacent layers
285  if(layerID1 == layerID2 - 1) return kTRUE;
286  else if(layerID1 == layerID2 + 1) return kTRUE;
287  }
288  // one inner(outer) one outer(inner)
289  else {
290  // must be on boundaries // CHECK isboundary
291  if(layerID1 == 7 && layerID2 == 16) return kTRUE;
292  else if(layerID1 == 16 && layerID2 == 7) return kTRUE;
293  }
294  }
295  // both skew
296  else if(hit1->IsSttSkew() && hit2->IsSttSkew()) {
297  PndSttTube *tube1 = (PndSttTube*) fTubeArray->At(hit1->GetTubeID());
298  int layerID1 = tube1->GetLayerID();
299  PndSttTube *tube2 = (PndSttTube*) fTubeArray->At(hit2->GetTubeID());
300  int layerID2 = tube2->GetLayerID();
301  // same layer
302  if(layerID1 == layerID2) return kTRUE;
303  // adjacent layers
304  else if(layerID1 == layerID2 - 1) return kTRUE;
305  else if(layerID1 == layerID2 + 1) return kTRUE;
306  }
307  }
308 
309  return kFALSE;
310 }
311 
312 
313 // vicinity of hits
314 PndTrkClusterList PndTrkClean::Cleanup2(PndTrkCluster *cluster) {
315 
316  // get most distance point
317  PndTrkHit *disthit = NULL;
318  double tmpdistance = -1;
319  for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) {
320  PndTrkHit *hit = cluster->GetHit(ihit);
321  double dist = hit->GetXYDistance(TVector3(0., 0., 0.));
322  if(dist > tmpdistance) {
323  disthit = hit;
324  tmpdistance = dist;
325 
326  }
327  }
328  // sort
329  for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) {
330  PndTrkHit *hit = cluster->GetHit(ihit);
331  double dist = hit->GetXYDistance(disthit);
332  hit->SetSortVariable(dist);
333  }
334  cluster->Sort();
335  //for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { //[R.K. 01/2017] unused variable
336  //PndTrkHit *hit = cluster->GetHit(ihit); //[R.K. 01/2017] unused variable
337  // cout << "sorted hit " << ihit << " " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << hit->GetSortVariable() << " " << hit->IsSortable() << endl;
338 
339  //} //[R.K. 01/2017] unused variable
340  std::vector< int > failedhits, breakpoints;
341  for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) {
342  PndTrkHit *hit = cluster->GetHit(ihit);
343 // cout << "HIT " << hit->GetHitID() << " " << hit->GetDetectorID() << endl;
344 // hit->GetPosition().Print();
345 
346  PndTrkHit *follow = cluster->GetNextHit(ihit);
347  PndTrkHit *previous = cluster->GetPreviousHit(ihit);
348  // cout << "pre " << previous << " fol " << follow << endl;
349  Bool_t pre = kTRUE, fol = kTRUE;
350  if(previous != NULL) pre = CheckPairOfHits(previous, hit);
351  if(follow != NULL) fol = CheckPairOfHits(follow, hit);
352  // cout << "pre " << pre << " fol " << fol << endl;
353 
354  // isolated hit wrong
355  if(pre == kFALSE && fol == kFALSE) {
356  // cout << "isolated " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl;
357  failedhits.push_back(ihit);
358  }
359  // first hit wrong
360  else if(previous == NULL && fol == kFALSE) {
361  PndTrkHit *follow2 = cluster->GetNextHit(ihit + 1);
362  if(CheckPairOfHits(follow, follow2) == kFALSE) {
363  failedhits.push_back(ihit);
364  // cout << "first " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl;
365 
366  }
367  }
368  // last hit wrong
369  else if(follow == NULL && pre == kFALSE) {
370  PndTrkHit *previous2 = cluster->GetPreviousHit(ihit - 1);
371  if(CheckPairOfHits(previous, previous2) == kFALSE) {
372  // cout << "last " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl;
373  failedhits.push_back(ihit);
374  }
375  }
376  // breakpoints
377  else if(pre == kTRUE && fol == kFALSE) {
378  if(find(failedhits.begin(), failedhits.end(), ihit - 1) == failedhits.end()) {
379  // cout << "BREAKPOINT " << hit->GetHitID() << " " << hit->GetDetectorID() << endl;
380  breakpoints.push_back(ihit);
381  }
382  }
383 
384  }
385 
386 
387 // cout << "CLEANUP PROCEDURE " << endl;
388 // cout << "failedhits " << failedhits.size() << " breakpoints " << breakpoints.size() << endl;
389 
390  // if(display == kTRUE)
391  {
392  for(size_t ihit = 0; ihit < failedhits.size(); ihit++) {
393 // cout << "failedhits " << failedhits.size() << " " << failedhits[ihit] << endl;
394  int hitno = failedhits[ihit];
395  PndTrkHit *hit = cluster->GetHit(hitno);
396  hit->Draw(kRed);
397  }
398 
399  for(size_t ihit = 0; ihit < breakpoints.size(); ihit++) {
400  int hitno = breakpoints[ihit];
401  PndTrkHit *hit = cluster->GetHit(hitno);
402  hit->Draw(kGreen);
403  }
404 
405  }
406 // cout << "before delete " << cluster-> GetNofHits() << endl;
407  cluster->DeleteHits(failedhits);
408 // cout << "dopo delete " << cluster-> GetNofHits() << endl;
409  //
410  // 0 1 2 3 4 5 6 7 8 9
411  // failed ihit 2 3 5
412  // 0 1 4 6 7 8 9
413  // break orig ihit 4 8
414  // break fix ihit 2 5
415  // fix the breakcounter
416  for(size_t jhit = 0; jhit < breakpoints.size(); jhit++) {
417  int counter = 0;
418  for(size_t ihit = 0; ihit < failedhits.size(); ihit++) {
419  if(failedhits[ihit] < breakpoints[jhit]) counter++;
420  else break;
421  }
422  breakpoints[jhit] -= counter;
423  }
424 
425 
426 // cout << "dopo delete " << cluster-> GetNofHits() << endl;
427 
428  // PndTrkClusterList list = Split(cluster, breakpoints);
429 
430  PndTrkClusterList list;
431  list.AddCluster(cluster);
432 // cout << "dopo delete " << cluster-> GetNofHits() << " " << list.GetNofClusters() << endl;
433  return list;
434 
435 }
436 
437 
438 // CHECK if it works
439 Bool_t PndTrkClean::SplitAtHit(PndTrkCluster *hitlist, PndTrkHit *athit, PndTrkCluster &cluster1, PndTrkCluster &cluster2) {
440 
441  //PndTrkHit *tmphit = hitlist->GetHit(0); //[R.K. 01/2017] unused variable?
442 
443  cluster1.AddHit(athit);
444  cluster2.AddHit(athit);
445  TVector3 tmpposition = athit->GetPosition();
446  int tmptubeid = athit->GetTubeID();
447  PndSttTube * tmptube = (PndSttTube *) fTubeArray->At(tmptubeid);
448  PndTrkHit *secondhit = NULL;
449 
450  // cout << "CLUS1: " << tmphit->GetHitID() << endl;
451  for(int ihit = 0; ihit < hitlist->GetNofHits(); ihit++) {
452  PndTrkHit *hit = hitlist->GetHit(ihit);
453  if(*hit == *athit) {
454  cout << "hit == ahit " << endl;
455  continue;
456  }
457  // if(hit->GetXYDistance(tmpposition) < STTPARALDISTANCE) { // CHECK
458  int tubeid = hit->GetTubeID();
459  PndSttTube * tube = (PndSttTube *) fTubeArray->At(tubeid);
460  int layerid = tube->GetLayerID();
461  tmptube = (PndSttTube *) fTubeArray->At(tmptubeid);
462  int tmplayerid = tmptube->GetLayerID();
463  cout << "*** " << hit->GetHitID() << " " << layerid << " " << tmplayerid << endl;
464  if(layerid == tmplayerid || layerid == tmplayerid + 1 || layerid == tmplayerid - 1) {
465 
466  // if(hit->GetTubeXYDistance(tmpposition) < STTPARALDISTANCE) { // CHECK
467  cluster1.AddHit(hit);
468  tmpposition = hit->GetPosition();
469  tmptubeid = tubeid;
470  if(secondhit == NULL) {
471  // cout << "secondhit " << endl;
472  secondhit = hit;
473  }
474  // else cout << "sechit " << secondhit << endl;
475  // cout << "CLUS1: " << hit->GetHitID() << endl;
476 
477  }
478  }
479 
480  // cout << "CLUS2: " << tmphit->GetHitID() << endl;
481 
482  tmpposition = athit->GetPosition();
483  for(int ihit = 0; ihit < hitlist->GetNofHits(); ihit++) {
484  PndTrkHit *hit = (PndTrkHit*) hitlist->GetHit(ihit);
485 // cout << "athit " << athit << endl;
486 // cout << "hit " << hit << endl;
487 // cout << "secondhit " << secondhit << endl;
488  if(*hit == *athit) continue;
489  if(*hit == *secondhit) continue;
490 
491  if(hit->GetXYDistance(tmpposition) < STTPARALDISTANCE) { // CHECK
492  cluster2.AddHit(hit);
493  tmpposition = hit->GetPosition();
494 // cout << "CLUS2: " << hit->GetHitID() << endl;
495  }
496  }
497 
498  return kTRUE;
499 }
500 
501 // CHECK if it works now
502 PndTrkClusterList PndTrkClean::Split(PndTrkCluster *cluster, std::vector< int > breakpoints) {
503 
504  PndTrkClusterList list;
505  if(breakpoints.size() == 0) {
506  list.AddCluster(cluster);
507 // cout << "0add cluster to list " << cluster->GetNofHits() << endl;
508  return list;
509  }
510 
511  // 0 1 2 3 4 5 6 7 8 9
512  // split @ ihit = 4 and 7
513  // 4 - 0 --> 0123 456789
514  // 7 - 4 --> 456 789
515 
516  PndTrkCluster cluster1, cluster2;
517  PndTrkCluster *tmpcluster = cluster;
518  for(size_t ihit = 0; ihit < breakpoints.size(); ihit++) {
519  int hitno = breakpoints[ihit] - cluster1.GetNofHits();
520  PndTrkHit *athit = tmpcluster->GetHit(hitno);
521  SplitAtHit(tmpcluster, athit, cluster1, cluster2);
522  tmpcluster = &cluster2;
523  list.AddCluster(&cluster1);
524 // cout << "1add cluster to list " << cluster1.GetNofHits() << endl;
525  }
526 
527 // cout << "2add cluster to list " << cluster2.GetNofHits() << endl;
528  list.AddCluster(&cluster2);
529  tmpcluster = NULL;
530  delete tmpcluster;
531  return list;
532 }
533 
534 
535 
536 
537 // SECTORS
538 
539 int PndTrkClean::CheckSectorDistribution(PndTrkCluster *cluster) {
540  int sectorcounter[6] = {0, 0, 0, 0, 0, 0}; // CHECK hard coded nof sectors
541  for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) {
542  PndTrkHit *hit = cluster->GetHit(ihit);
543  PndSttTube *tube = (PndSttTube*) fTubeArray->At(hit->GetTubeID());
544  int sectorID = tube->GetSectorID();
545  sectorcounter[sectorID]++;
546  }
547 
548  int tmpsector = -1, tmpcounter = -1;
549  for(int isec = 0; isec < 6; isec++) { // CHECK hard coded nof sectors
550  if(sectorcounter[isec] > tmpcounter) {
551  tmpcounter = sectorcounter[isec];
552  tmpsector = isec;
553  }
554  }
555  return tmpsector;
556 }
557 
558 PndTrkCluster PndTrkClean::CleanSectors(PndTrkCluster *cluster, int sector) {
559  PndTrkCluster newcluster;
560  for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) {
561  PndTrkHit *hit = cluster->GetHit(ihit);
562  PndSttTube *tube = (PndSttTube*) fTubeArray->At(hit->GetTubeID());
563  int sectorID = tube->GetSectorID();
564  if(fabs(sectorID - sector) < 1.2) newcluster.AddHit(hit);
565  }
566  return newcluster;
567 }
568 
569 // MERGING
570 
571 PndTrkClusterList PndTrkClean::MergeClusters(PndTrkClusterList *clusterlist) {
572  int nmax = clusterlist->GetNofClusters();
573 
574  cout << "MERGING on " << nmax << endl;
575  std::vector< int> merged;
576  PndTrkClusterList newclusterlist;
577  for(int imax = 0; imax < nmax; imax++) {
578 
579  if(find(merged.begin(), merged.end(), imax) != merged.end()) continue;
580 
581  PndTrkLegendreCluster *iclus = (PndTrkLegendreCluster*) clusterlist->GetCluster(imax);
582 
583  // dont try to merge 1 cluster alone
584  if(nmax < 2) {
585  newclusterlist.AddCluster(iclus);
586  continue;
587  }
588  PndTrkLegendreCluster *merge = iclus;
589 
590  for(int jmax = imax + 1; jmax < nmax; jmax++) {
591  cout << "COMPARING " << imax << " " << jmax << endl;
592  PndTrkLegendreCluster *jclus = (PndTrkLegendreCluster*) clusterlist->GetCluster(jmax);
593  bool isSimilar = iclus->IsSimilarTo(jclus);
594  if(isSimilar) {
595  cout << imax << " is similar to " << jmax << endl;
596  merge->MergeTo(jclus); // CHECK //int nhits = //[R.K.03/2017] unused variable
597  merged.push_back(jmax);
598  cout << "MERGING " << iclus->GetNofHits() << " + " << jclus->GetNofHits() << " = " << merge->GetNofHits() << endl;
599  // iclus->Print();
600  // jclus->Print();
601  // merge->Print();
602  cout << "theta/r " << merge->GetTheta() << " " << merge->GetR() << endl;
603  }
604  }
605 
606 
607  newclusterlist.AddCluster(merge);
608  cout << "ADD CLUSTER " << merge->GetNofHits() << endl;
609  }
610  cout << "FINAL CHECK "<< newclusterlist.GetNofClusters() << endl;
611  for(int iclus = 0; iclus < newclusterlist.GetNofClusters(); iclus++) {
612  PndTrkLegendreCluster *cluster = (PndTrkLegendreCluster* ) newclusterlist.GetCluster(iclus);
613  cluster->Print();
614  cout << "theta-r " << cluster->GetTheta() << " " << cluster->GetR() << endl;
615 
616  }
617 
618  return newclusterlist;
619 }
620 
621 ClassImp(PndTrkClean)
622 
PndTrkClean::FindMvdLayer
int FindMvdLayer(int sensorID)
Definition: PndTrkClean.cxx:165
PndTrkClean::MergeClusters
PndTrkClusterList MergeClusters(PndTrkClusterList *clusterlist)
Definition: PndTrkClean.cxx:571
PndTrkClean::fTubeArray
TClonesArray * fTubeArray
Definition: PndTrkClean.h:35
PndTrkHit::SetSortVariable
void SetSortVariable(Double_t sortvar)
Definition: PndTrkHit.h:44
PndTrkCluster::AddHit
void AddHit(PndTrkHit *hit)
Definition: PndTrkCluster.cxx:70
PndTrkClean::CheckPairOfHits
Bool_t CheckPairOfHits(PndTrkHit *hit1, PndTrkHit *hit2)
Definition: PndTrkClean.cxx:238
PndGeoHandling::GetPath
TString GetPath(Int_t shortID)
for a given shortID the path is returned
Definition: PndGeoHandling.cxx:222
v
__m128 v
Definition: P4_F32vec4.h:4
PndTrkCluster::GetNextHit
PndTrkHit * GetNextHit(int index)
Definition: PndTrkCluster.cxx:657
PndTrkCluster::GetHit
PndTrkHit * GetHit(int index)
Definition: PndTrkCluster.cxx:643
PndTrkLegendreCluster::MergeTo
int MergeTo(PndTrkLegendreCluster *cluster2)
Definition: PndTrkLegendreCluster.cxx:51
PndTrkClusterList::AddCluster
void AddCluster(PndTrkCluster *cluster)
Definition: PndTrkClusterList.cxx:35
PndTrkHit::GetHitID
Int_t GetHitID()
Definition: PndTrkHit.h:56
PndTrkHit::IsSttParallel
Bool_t IsSttParallel()
Definition: PndTrkHit.h:70
PndTrkHit::IsMvd
Bool_t IsMvd()
Definition: PndTrkHit.h:77
PndSttTube::GetLayerID
int GetLayerID()
Definition: PndSttTube.cxx:128
isec
int isec
Definition: f_Init.h:19
counter
int counter
Definition: ZeeAnalysis.C:59
PndTrkCluster::DeleteHits
void DeleteHits(std::vector< int > todelete)
Definition: PndTrkCluster.cxx:124
merge
int merge(TString ntp="", TString fout="", TString f1="", TString f2="", TString f3="", TString f4="", TString f5="")
Definition: examples/analysistools/merge.C:1
PndGeoHandling::GetOUVPath
void GetOUVPath(TString path, TVector3 &o, TVector3 &u, TVector3 &v)
for a volume given by its path the o, u, v vectors for the plane are returned
Definition: PndGeoHandling.cxx:320
PndTrkCluster::IsSimilarTo
Bool_t IsSimilarTo(PndTrkCluster *cluster2)
Definition: PndTrkCluster.cxx:665
PndTrkHit::GetDistance
Double_t GetDistance(PndTrkHit *fromhit)
Definition: PndTrkHit.cxx:83
STTPARALDISTANCE
#define STTPARALDISTANCE
Definition: PndTrkParameters.h:5
PndTrkClean::SplitAtHit
Bool_t SplitAtHit(PndTrkCluster *hitlist, PndTrkHit *athit, PndTrkCluster &cluster1, PndTrkCluster &cluster2)
Definition: PndTrkClean.cxx:439
PndTrkClean::CheckSectorDistribution
int CheckSectorDistribution(PndTrkCluster *cluster)
Definition: PndTrkClean.cxx:539
PndTrkHit::GetPosition
TVector3 GetPosition()
Definition: PndTrkHit.h:62
PndTrkClean::Split
PndTrkClusterList Split(PndTrkCluster *cluster, std::vector< int > breakpoints)
Definition: PndTrkClean.cxx:502
PndTrkClean::CleanSectors
PndTrkCluster CleanSectors(PndTrkCluster *cluster, int sector)
Definition: PndTrkClean.cxx:558
PndTrkClean::fGeoH
PndGeoHandling * fGeoH
Definition: PndTrkClean.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndGeoHandling::Instance
static PndGeoHandling * Instance()
Definition: PndGeoHandling.cxx:36
PndTrkCluster::GetPreviousHit
PndTrkHit * GetPreviousHit(int index)
Definition: PndTrkCluster.cxx:651
PndTrkHit::Draw
void Draw(Color_t color)
Definition: PndTrkHit.cxx:109
PndTrkHit::GetXYDistance
Double_t GetXYDistance(PndTrkHit *fromhit)
Definition: PndTrkHit.cxx:94
PndTrkHit::IsStt
Bool_t IsStt()
Definition: PndTrkHit.h:72
ClassImp
ClassImp(PndAnaContFact)
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndTrkClusterList::GetCluster
PndTrkCluster * GetCluster(Int_t index)
Definition: PndTrkClusterList.h:34
PndTrkCluster::GetNofHits
Int_t GetNofHits()
Definition: PndTrkCluster.h:52
PndSttTube::GetSectorID
int GetSectorID()
Definition: PndSttTube.cxx:124
PndTrkHit::GetTubeID
Int_t GetTubeID()
Definition: PndTrkHit.h:61
PndTrkHit::IsSttSkew
Bool_t IsSttSkew()
Definition: PndTrkHit.h:71
PndTrkClean::Cleanup2
PndTrkClusterList Cleanup2(PndTrkCluster *cluster)
Definition: PndTrkClean.cxx:314
Generated on Wed Apr 3 2019 08:44:03 for FairRoot/PandaRoot by   doxygen 1.8.5