LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
CFAlgo3DAngle.cxx
Go to the documentation of this file.
1 #ifndef RECOTOOL_CFALGO3DANGLE_CXX
2 #define RECOTOOL_CFALGO3DANGLE_CXX
3 
4 #include "CFAlgo3DAngle.h"
5 
6 namespace cmtool {
7 
8  //-------------------------------------------------------
10  //-------------------------------------------------------
11  {
12 
13  SetDebug(false) ;
14  SetThetaCut(30) ;
15  SetPhiCut(40) ;
16  SetRatio(0.0005) ;
17 
18  }
19 
20  //-----------------------------
22  //-----------------------------
23  {
24 
25  }
26 
27  //----------------------------------------------------------------------------------------------
28  float CFAlgo3DAngle::Float(const std::vector<const cluster::ClusterParamsAlg*> &clusters)
29  //----------------------------------------------------------------------------------------------
30  {
31  // Code-block by Kazu starts
32  // This ensures the algorithm works only if # clusters is > 2 (and not =2)
33  // You may take out this block if you want to allow matching using clusters from only 2 planes.
34  if(clusters.size()==2 || clusters.size()==1) return -1;
35  // Code-block by Kazu ends
36 
37 
38  int plane_0 = clusters.at(0)->Plane();
39  int plane_1 = clusters.at(1)->Plane();
40  int plane_2 = clusters.at(2)->Plane();
41  double angle_2d_0 = clusters.at(0)->GetParams().angle_2d;
42  double angle_2d_1 = clusters.at(1)->GetParams().angle_2d;
43  double angle_2d_2 = clusters.at(2)->GetParams().angle_2d;
44  auto hits_0 = clusters.at(0)->GetParams().N_Hits ;
45  auto hits_1 = clusters.at(1)->GetParams().N_Hits ;
46  auto hits_2 = clusters.at(2)->GetParams().N_Hits ;
47 
48  auto sumCharge0 = clusters.at(0)->GetParams().sum_charge;
49  auto sumCharge1 = clusters.at(1)->GetParams().sum_charge;
50  auto sumCharge2 = clusters.at(2)->GetParams().sum_charge;
51 
52  //Calculate angles theta and phi for cluster pairs across 2 planes
53  //Theta goes from -90 to 90, phi from -180 to 180
54  double phi_01 = 0;
55  double theta_01 = 0;
56  double phi_12 = 0;
57  double theta_12 = 0;
58  double phi_02 = 0;
59  double theta_02 = 0;
60 
61  double max_phi(0), middle_phi(0), min_phi(0);
62  double max_theta(0), middle_theta(0), min_theta(0);
63  //When theta1, theta2, phi ratios are calcualted, find best 2
64  double max_ratio(0), middle_ratio(0), min_ratio(0);
65 
66  //For ordering hits and rejecting pairs which have a large difference in hit number
67  double max_hits(0), middle_hits(0), min_hits(0) ;
68 // double max_hits1(0), middle_hits1(0), min_hits1(0) ;
69 
70  //Calculate phi and theta from first 2 planes; check if third plane is consistent
72  geou.Get3DaxisN(plane_0,plane_1,angle_2d_0,angle_2d_1,phi_01,theta_01);
73  geou.Get3DaxisN(plane_1,plane_2,angle_2d_1,angle_2d_2,phi_12,theta_12);
74  geou.Get3DaxisN(plane_2,plane_0,angle_2d_2,angle_2d_0,phi_02,theta_02);
75 
76  //Adjust the range of phis/thetas that are bigger than 360 or less than 0.
77  FixPhiTheta(phi_01,theta_01);
78  FixPhiTheta(phi_12,theta_12);
79  FixPhiTheta(phi_02,theta_02);
80 
81  //Order phi and theta for ratio calculation later
82  SetMaxMiddleMin(phi_01,phi_12,phi_02,max_phi,middle_phi,min_phi);
83  SetMaxMiddleMin(theta_01,theta_12,theta_02,max_theta,middle_theta,min_theta);
84 
85  //Order hits from most to least
86  SetMaxMiddleMin(hits_0,hits_1,hits_2,max_hits,middle_hits,min_hits);
87 
88  //Ratio for hits
89  double ratio_max_min = 1;
90  double ratio_max_middle = 1;
91 
92  //Ratio for theta angles
93  double ratio_theta1 = 1;
94  double ratio_theta2 = 1;
95 
96  //Ratio for phi--only 1 ratio because 2 of the angles come from collection plane 2
97  //and are thus the same.
98  double ratio_phi = 1;
99 
100  //Total ratio
101  double ratio = 1;
102 
103  //This takes into account the fact that 0 and 360 having the same relative value (for phi; 0 and 180 for theta)
104  for(int i=0; i<2 ; i++){
105  while(min_phi + 360 < max_phi +_phi_cut && min_phi +360 > max_phi - _phi_cut)
106  {
107  min_phi +=360 ;
108  SetMaxMiddleMin(max_phi, middle_phi, min_phi,max_phi,middle_phi,min_phi);
109  }
110  }
111 
112  ratio_theta2 = min_theta / max_theta;
113  ratio_theta1 = middle_theta / max_theta;
114  ratio_phi = min_phi / max_phi ;
115 
116  //Get the biggest ratios for the thetas and phi
117  SetMaxMiddleMin(ratio_phi, ratio_theta1, ratio_theta2, max_ratio, middle_ratio, min_ratio);
118 
119  ratio_max_min = min_hits / max_hits ;
120  ratio_max_middle = middle_hits / max_hits ;
121 
122  ratio = max_ratio * sumCharge1* sumCharge2 * sumCharge0; //* ; //* ratio_phi ; //* ratio_max_middle ;
123 
124  //Test to make sure that max hits is not too much bigger than min
125  if( ratio_max_min <0.3)
126  ratio *= ratio_max_min ;
127 
128  //GeometryUtilities returns theta=-999 when 2d Angle=0--Don't know
129  //how else to deal with that. This seems to work.
130  if( theta_01 == -999 || theta_12 == -999 || theta_02 == -999)
131  return -1 ;
132 
133 
134  if(_debug && ratio > _ratio_cut ){
135  std::cout<<"\nNhits planes 0, 1, 2: " <<clusters.at(0)->GetParams().N_Hits<<", "<<clusters.at(1)->GetParams().N_Hits
136  <<", "<<clusters.at(2)->GetParams().N_Hits ;
137  std::cout<<"\nTheta1 , Theta2 : "<<ratio_theta1<<", "<<ratio_theta2;
138  std::cout<<"\nPhi Ratio: "<<ratio_phi;
139  std::cout<<"\nHits ratio mid : "<<ratio_max_middle ;
140  // std::cout<<"\nHits ratio min : "<<ratio_max_min ;
141  std::cout<<"\nTotal is: " <<ratio<<" ***************";
142 
143 
144  std::cout<<"\n\n0: 2dAngle: "<<clusters.at(0)->GetParams().cluster_angle_2d<<std::endl;
145  std::cout<<"1: 2dAngle: "<<clusters.at(1)->GetParams().cluster_angle_2d<<std::endl;
146  std::cout<<"2: 2dAngle: "<<clusters.at(2)->GetParams().cluster_angle_2d<<std::endl;
147  std::cout<<"\nTheta,Phi MaxMM : "<<max_theta<<", "<<middle_theta<<", "<<min_theta<<"\n\t\t"
148  <<max_phi<<", "<<middle_phi<<", "<<min_phi;
149 
150  std::cout<<"\nStart End Points: "<<clusters.at(0)->GetParams().start_point.t<<", "<<clusters.at(0)->GetParams().end_point.t<<"\n\t\t"
151  <<clusters.at(1)->GetParams().start_point.t<<", "<<clusters.at(1)->GetParams().end_point.t<<"\n\t\t "
152  <<clusters.at(2)->GetParams().start_point.t<<", "<<clusters.at(2)->GetParams().end_point.t;
153 
154 /* std::cout<<"\nFor planes 0 and 1: "<<std::endl ;
155  std::cout<<"\tPhi : "<<phi_01<<std::endl;
156  std::cout<<"\tTheta : "<<theta_01<<std::endl ;
157  std::cout<<"For planes 1 and 2: "<<std::endl ;
158  std::cout<<"\tPhi : "<<phi_12<<std::endl;
159  std::cout<<"\tTheta : "<<theta_12<<std::endl ;
160  std::cout<<"For plane 0 and 2: "<<std::endl ;
161  std::cout<<"\tPhi : "<<phi_02<<std::endl ;
162  std::cout<<"\tTheta : "<<theta_02<<std::endl; */
163  // std::cout<<"\nNEW CLUSTERS PAIRS NOW\n\n\n"<<std::endl<<std::endl;
164  std::cout<<"\n\n\n";
165  }
166 
167  return(ratio > _ratio_cut ? ratio : -1) ;
168  }
169 
170  //--------------------------------
171  void CFAlgo3DAngle::FixPhiTheta(double &phi, double &theta)
172  //--------------------------------
173  {
174  // while(phi <= 0)
175  // phi += 360 ;
176  // while(phi > 720)
177  // phi -= 360 ;
178 // if(phi <=0)
179 
180  phi+=360;
181 
182  if(theta != -999)
183  theta += 180 ;
184 
185  }
186 
187 
188  //------------------------------
189  void CFAlgo3DAngle::SetMaxMiddleMin(const double first, const double second, const double third, double &max, double &middle, double &min)
190  //------------------------------
191  {
192 
193  if(first > second && first > third){
194  max = first;
195  }
196  else if (first > second && first < third){
197  max = third ;
198  middle = first ;
199  min = second ;
200  }
201  else if (first > third && first < second){
202  max = second ;
203  middle = first ;
204  min = third ;
205  }
206  else if(first <second && first <third)
207  min = first ;
208 
209 
210  if (max == first && second > third){
211  middle = second ;
212  min = third ;
213  }
214  else if (max ==first && third > second){
215  middle = third ;
216  min = second ;
217  }
218 
219  if(min ==first && second > third){
220  middle = third ;
221  max = second;
222  }
223  else if(min ==first && third > second){
224  middle = second ;
225  max = third ;
226  }
227 
228 
229  //Very rarely, the angles(or hits) may be equal
230  if( first == second && first > third ){
231  max = first;
232  middle = second ;
233  min = third ;
234  }
235  else if( first == second && first < third){
236  max = third;
237  middle = first ;
238  min = second ;
239  }
240 
241  else if( first ==third && first > second){
242  max = first;
243  middle = third;
244  min = second;
245  }
246 
247  else if( first == third && first < second){
248  max = second ;
249  middle = first;
250  min = third ;
251  }
252 
253  else if( second ==third && second < first){
254  max = first;
255  middle = third;
256  min = second;
257  }
258 
259  else if( second == third && second > first){
260  max = second;
261  middle = third;
262  min = first ;
263  }
264 
265 
266 }
267 
268 
269  //------------------------------
271  //------------------------------
272  {
273 
274  }
275 
276 }
277 #endif
void SetRatio(float ratio)
Definition: CFAlgo3DAngle.h:56
Class def header for a class CFAlgo3DAngle.
virtual void SetMaxMiddleMin(const double first, const double second, const double third, double &most, double &middle, double &least)
CFAlgo3DAngle()
Default constructor.
Int_t max
Definition: plot.C:27
void SetDebug(bool debug)
Definition: CFAlgo3DAngle.h:50
Int_t min
Definition: plot.C:26
virtual float Float(const std::vector< const cluster::ClusterParamsAlg * > &clusters)
void SetThetaCut(float theta_cut)
Definition: CFAlgo3DAngle.h:52
void SetPhiCut(float phi_cut)
Definition: CFAlgo3DAngle.h:54
virtual void FixPhiTheta(double &phi, double &theta)
Int_t Get3DaxisN(Int_t iplane0, Int_t iplane1, Double_t omega0, Double_t omega1, Double_t &phi, Double_t &theta) const
virtual void Reset()
Function to reset the algorithm instance called within CMergeManager/CMatchManager&#39;s Reset() ...