ht_test.cxx 2.45 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include "GenFind/GenFindHits.h"
#include "GenFind/HoughTransform.h"

void ht_test()
{
  using namespace genfind;

  std::vector<std::shared_ptr<Hit>> hits; 
  std::vector<ConformalHit> chits; 
  std::vector<ROOT::Math::XYZTVector> hits2; 
  std::vector<ConformalHit> chits2; 
  double degree = TMath::Pi()/180.0;

  double r0  = 10.0;
  double x0  = 0.0;
  double y0  = r0;
    ROOT::Math::XYZTVector ref_hit = {0.0,0.0,0.0,0.0};
  for(int ith=0; ith < 10; ith++) {
    double x = r0*TMath::Sin((5.0 + ith*3.0)*degree)- x0 ;
    double y = TMath::Sqrt(r0*r0 - (x+x0)*(x+x0)) - y0 ;
    auto ahit = std::make_shared<Hit>(x,y,0);
    hits.push_back(ahit);
    chits.push_back(compute_conformal_hit(ahit));
    std::cout << " x " << x << ", y " << y << std::endl;

    ROOT::Math::XYZTVector pos_hit = {x,y,0.0,0.0};
    hits2.push_back(pos_hit);
    //chits2.push_back(compute_conformal_hit(pos_hit, ref_hit));
  }
  chits2 = compute_conformal_hits(hits2,ref_hit);

  HoughTransform ht;
  auto thits = ht(hits2);
34

35
36
37
38
39
40
  std::cout << "ht test\n";
  for(auto ahit : thits[0]){
    std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
  }
  std::cout << "ht test\n";

41
42
43
44
45
  auto peaks = ht.FindPhiPeaks(chits2);
  for(auto ahit : peaks){
    std::cout << " phi " << ahit << std::endl;
  }

46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
  TGraph* gr_hits = new TGraph();
  int i_point = 0;
  for(auto ahit : hits){
     gr_hits->SetPoint(i_point, ahit->X(), ahit->Y());
     i_point++;
  }

  TGraph* gr_chits = new TGraph();
  i_point = 0;
  for(auto ahit : chits){
    gr_chits->SetPoint(i_point, ahit.X(), ahit.Y());
    std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
    i_point++;
  }

  TGraph* gr_hits2 = new TGraph();
  i_point = 0;
63
  for(auto ahit : thits[0]){
64
65
66
67
68
69
70
71
72
73
74
75
76
     gr_hits2->SetPoint(i_point, ahit.X(), ahit.Y());
     i_point++;
  }

  TGraph* gr_chits2 = new TGraph();
  i_point = 0;
  for(auto ahit : chits2){
    gr_chits2->SetPoint(i_point, ahit.X(), ahit.Y());
    std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
    i_point++;
  }
  
  TCanvas* c = new TCanvas();
77
  TMultiGraph* mg = new TMultiGraph();
78
79
80
81
  c->Divide(2,2);

  c->cd(1);
  gr_hits->SetMarkerStyle(20);
82
83
84
85
86
  gr_hits2->SetMarkerStyle(21);
  gr_hits2->SetMarkerColor(2);
  mg->Add(gr_hits,"p");
  mg->Add(gr_hits2,"p");
  mg->Draw("a");
87
88
89
90
91
92
93
94
95

  c->cd(2);
  gr_chits->SetMarkerStyle(20);
  gr_chits->Draw("ap");

  c->cd(3);
  gr_hits2->Draw("ap");

  c->cd(4);
96
97
98
  ht.fPhi->DrawCopy();
  //gr_chits2->SetMarkerStyle(20);
  //gr_chits2->Draw("ap");
99
100
101


}