LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
TrackMomentumCalculator.cxx
Go to the documentation of this file.
1 // \file TrackMomentumCalculator.cxx
2 //
3 // \author sowjanyag@phys.ksu.edu
4 
5 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
6 
7 double my_mcs_chi2( const double *x )
8 {
9  Double_t result = 0.0;
10 
11  Double_t p = x[0];
12 
13  Double_t theta0 = x[1];
14 
15  for ( Int_t i=0; i<nmeas; i++ )
16  {
17  Double_t xx = xmeas[i];
18 
19  Double_t yy = ymeas[i];
20 
21  Double_t ey = eymeas[i];
22 
23  Double_t rad_length = 14.0;
24 
25  Double_t l0 = xx/rad_length;
26 
27  Double_t res1 = 0.0;
28 
29  if ( xx>0 && p>0 ) res1 = ( 13.6/p )*sqrt( l0 )*( 1.0+0.038*TMath::Log( l0 ) );
30 
31  res1 = sqrt( res1*res1+theta0*theta0 );
32 
33  Double_t diff = yy-res1;
34 
35  if ( ey==0 ) { cout << " Zero denominator in my_mcs_chi2 ! " << endl; return -1; }
36 
37  Double_t incre = pow( diff/ey, 2.0 );
38 
39  result += incre;
40 
41  }
42 
43  result += 2.0/( 4.6 )*theta0; // *TMath::Log( 1.0/14.0 );
44 
45  if ( isnan( float(result) ) || isinf( float(result) ) ) {
46  LOG_DEBUG("TrackMomentumCalculator")<<" Is nan in my_mcs_chi2 ! ";
47  return -1; }
48 
49  return result;
50 
51 }
52 
53 namespace trkf{
54 
55  TrackMomentumCalculator::TrackMomentumCalculator()
56  {
57  n = 0;
58 
59  n_reco = 0;
60 
61  seg_stop = -1.0; n_seg = 0;
62 
63  gr_reco_xyz = new TPolyLine3D(); gr_reco_xy = new TGraph(); gr_reco_yz = new TGraph(); gr_reco_xz = new TGraph();
64 
65  gr_seg_xyz = new TPolyLine3D(); gr_seg_xy = new TGraph(); gr_seg_yz = new TGraph(); gr_seg_xz = new TGraph();
66 
67  steps_size = 10.0; n_steps = 6; for ( Int_t i=1; i<=n_steps; i++ ) { steps.push_back( steps_size*i ); }
68 
69  basex.SetXYZ( 1.0, 0.0, 0.0 ); basey.SetXYZ( 0.0, 1.0, 0.0 ); basez.SetXYZ( 0.0, 0.0, 1.0 );
70 
71  nmeas = 0; p_mcs = -1.0; p_mcs_e = -1.0; chi2 = -1.0;
72 
73  steps_size2 = 10.0;
74 
75  p_mcs_2 = -1.0; LLbf = -1.0;
76 
77  kcal = 0.0024;
78 
79  minLength = 100;
80 
81  maxLength = 1350.0;
82 
83  Float_t Range_grampercm[29] = {9.833E-1, 1.786E0, 3.321E0, 6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2, 2.385E2, 4.934E2,
84  6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
85  4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4, 4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5};
86  for (int i = 0; i<29; ++i){
87  Range_grampercm[i] /= 1.396; //convert to cm
88  }
89  Float_t KE_MeV[29] = {10, 14, 20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000, 4000, 8000, 10000, 14000, 20000, 30000,
90  40000, 80000, 100000, 140000, 200000, 300000, 400000};
91  KEvsR = new TGraph(29, Range_grampercm, KE_MeV);
92  KEvsR_spline3 = new TSpline3("KEvsRS",KEvsR);
93  //KEvsRFromData = tfs->make<TH2D>("KEvsRFromData","KE vs R from Data",1000,0,1000,1000,0,5000);
94 
95  }
96 
97  TrackMomentumCalculator::~TrackMomentumCalculator(){
98 
99  delete gr_reco_xyz;
100  delete gr_reco_xy;
101  delete gr_reco_yz;
102  delete gr_reco_xz;
103  delete gr_seg_xyz;
104  delete gr_seg_xy;
105  delete gr_seg_yz;
106  delete gr_seg_xz;
107  delete KEvsR;
108  delete KEvsR_spline3;
109  delete gr_meas;
110  delete gr_xyz;
111  delete gr_yz;
112  delete gr_xz;
113  delete gr_xy;
114 
115  }
116 
117  double TrackMomentumCalculator::GetTrackMomentum(double trkrange, int pdg)
118  {
119 
120  /* Muon range-momentum tables from CSDA (Argon density = 1.4 g/cm^3)
121  website: http://pdg.lbl.gov/2012/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_289.pdf
122 
123  CSDA table values:
124  Float_t Range_grampercm[30] = {9.833E-1, 1.786E0, 3.321E0, 6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2, 2.385E2, 4.934E2,
125  6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
126  4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4, 4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5};
127  Float_t KE_MeV[30] = {10, 14, 20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000, 4000, 8000, 10000, 14000, 20000, 30000,
128  40000, 80000, 100000, 140000, 200000, 300000, 400000};
129 
130  Functions below are obtained by fitting polynomial fits to KE_MeV vs Range (cm) graph. A better fit was obtained
131  by splitting the graph into two: Below range<=200cm,a polynomial of power 4 was a good fit; above 200cm, a
132  polynomial of power 6 was a good fit
133 
134  Fit errors for future purposes:
135  Below 200cm, Forpoly4 fit: p0 err=1.38533;p1 err=0.209626; p2 err=0.00650077; p3 err=6.42207E-5; p4 err=1.94893E-7;
136  Above 200cm, Forpoly6 fit: p0 err=5.24743;p1 err=0.0176229; p2 err=1.6263E-5; p3 err=5.9155E-9; p4 err=9.71709E-13;
137  p5 err=7.22381E-17;p6 err=1.9709E-21;*/
138 
140  //*********For muon, the calculations are valid up to 1.91E4 cm range corresponding to a Muon KE of 40 GeV**********//
142 
143  /*Proton range-momentum tables from CSDA (Argon density = 1.4 g/cm^3):
144  website: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html
145 
146  CSDA values:
147  Double_t KE_MeV_P_Nist[31]={10, 15, 20, 30, 40, 80, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
148  1500, 2000, 2500, 3000, 4000, 5000};
149 
150  Double_t Range_gpercm_P_Nist[31]={1.887E-1,3.823E-1, 6.335E-1, 1.296, 2.159, 7.375, 1.092E1, 2.215E1, 3.627E1, 5.282E1, 7.144E1, 9.184E1, 1.138E2,
151  1.370E2, 1.614E2, 1.869E2, 2.132E2, 2.403E2, 2.681E2, 2.965E2, 3.254E2, 3.548E2, 3.846E2, 4.148E2, 4.454E2,
152  7.626E2, 1.090E3, 1.418E3, 1.745E3, 2.391E3, 3.022E3};
153 
154  Functions below are obtained by fitting power and polynomial fits to KE_MeV vs Range (cm) graph. A better fit was obtained
155  by splitting the graph into two: Below range<=80cm,a a*(x^b) was a good fit; above 80cm, a
156  polynomial of power 6 was a good fit
157 
158  Fit errors for future purposes:
159  For power function fit: a=0.388873; and b=0.00347075
160  Forpoly6 fit: p0 err=3.49729;p1 err=0.0487859; p2 err=0.000225834; p3 err=4.45542E-7; p4 err=4.16428E-10;
161  p5 err=1.81679E-13;p6 err=2.96958E-17;*/
162 
164  //*********For proton, the calculations are valid up to 3.022E3 cm range corresponding to a Muon KE of 5 GeV**********//
166 
167  if (trkrange<0 || std::isnan(trkrange)){
168  mf::LogError("TrackMomentumCalculator")<<"Invalid track range "<<trkrange<<" return -1";
169  return -1.;
170  }
171 
172  double KE, Momentum, Muon_M = 105.7, Proton_M = 938.272, M;
173 
174  if (abs(pdg) == 13){
175  M = Muon_M;
176  /*if (trkrange>0 && trkrange<=200)
177  KE = 10.658+ (3.71181*trkrange) + (-0.0302898*trkrange*trkrange) +
178  (0.000228501*trkrange*trkrange*trkrange)+
179  (-5.86201E-7*trkrange*trkrange*trkrange*trkrange);
180  else if (trkrange>200 && trkrange<=1.91E4)
181  KE = 30.6157+ (2.11089*trkrange) + (0.000255267*trkrange*trkrange)+
182  (-5.19254E-8*trkrange*trkrange*trkrange)+
183  (6.39454E-12*trkrange*trkrange*trkrange*trkrange)+
184  (-3.99392E-16*trkrange*trkrange*trkrange*trkrange*trkrange)+
185  (9.73174E-21*trkrange*trkrange*trkrange*trkrange*trkrange*trkrange);
186  else
187  KE = -999;*/
188  KE = KEvsR_spline3->Eval(trkrange);
189  }
190  else if (abs(pdg) == 2212){
191  M = Proton_M;
192  if (trkrange>0 && trkrange<=80)
193  KE = 29.9317*TMath::Power(trkrange,0.586304);
194  else if (trkrange>80 && trkrange<=3.022E3)
195  KE = 149.904+ (3.34146*trkrange) + (-0.00318856*trkrange*trkrange)+
196  (4.34587E-6*trkrange*trkrange*trkrange)+
197  (-3.18146E-9*trkrange*trkrange*trkrange*trkrange)+
198  (1.17854E-12*trkrange*trkrange*trkrange*trkrange*trkrange)+
199  (-1.71763E-16*trkrange*trkrange*trkrange*trkrange*trkrange*trkrange);
200  else
201  KE = -999;
202  }
203  else
204  KE = -999;
205 
206  if (KE<0)
207  Momentum = -999;
208  else
209  Momentum = TMath::Sqrt((KE*KE)+(2*M*KE));
210 
211  Momentum = Momentum/1000;
212 
213  return Momentum;
214  }
215 
216  /*double TrackMomentumCalculator::GetProtonTrackMomentum(double trkrange)
217  {
218  double ProtonKE;
219  if (trkrange<=80)
220  ProtonKE = 29.9317*TMath::Power(trkrange,0.586304);
221  else
222  ProtonKE = 149.904+ (3.34146*trkrange) + (-0.00318856*trkrange*trkrange)+
223  (4.34587E-6*trkrange*trkrange*trkrange)+
224  (-3.18146E-9*trkrange*trkrange*trkrange*trkrange)+
225  (1.17854E-12*trkrange*trkrange*trkrange*trkrange*trkrange)+
226  (-1.71763E-16*trkrange*trkrange*trkrange*trkrange*trkrange*trkrange);
227 
228  return ProtonKE;
229  }*/
230 
231  // Momentum measurement via Multiple Coulomb Scattering (MCS)
232 
233  // author: Leonidas N. Kalousis (July 2015)
234 
235  // email: kalousis@vt.edu
236 
237  Double_t TrackMomentumCalculator::GetMomentumMultiScatterLLHD( const art::Ptr<recob::Track> &trk )
238  {
239  Double_t p = -1.0;
240 
241  std::vector<Float_t> recoX; std::vector<Float_t> recoY; std::vector<Float_t> recoZ;
242 
243  recoX.clear(); recoY.clear(); recoZ.clear();
244 
245  Int_t n_points = trk->NumberTrajectoryPoints();
246 
247  for ( Int_t i=0; i<n_points; i++ )
248  {
249  const TVector3 &pos = trk->LocationAtPoint( i );
250 
251  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
252 
253  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl; getchar();
254 
255  }
256 
257  Int_t my_steps = recoX.size();
258 
259  if ( my_steps<2 ) return -1.0;
260 
261  Int_t check0 = GetRecoTracks( recoX, recoY, recoZ );
262 
263  if ( check0!=0 ) return -1.0;
264 
265  seg_size = steps_size2;
266 
267  Int_t check1 = GetSegTracks2( recoX, recoY, recoZ );
268 
269  if ( check1!=0 ) return -1.0;
270 
271  Int_t seg_steps = segx.size();
272 
273  if ( seg_steps<2 ) return -1;
274 
275  Int_t seg_steps0 = seg_steps-1;
276 
277  Double_t recoL = segL.at(seg_steps0);
278 
279  if ( recoL<minLength || recoL>maxLength ) return -1;
280 
281  Int_t check2 = GetDeltaThetaij( dEi, dEj, dthij, seg_size, ind );
282 
283  if ( check2!=0 ) return -1.0;
284 
285  Double_t logL = 1e+16;
286 
287  Double_t bf = -666.0; // Double_t errs = -666.0;
288 
289  Double_t start1 = 0.0; Double_t end1 = 750.0;
290 
291  Double_t start2 = 0.0; Int_t end2 = 0.0; // 800.0;
292 
293  for ( Int_t k=start1; k<=end1; k++ )
294  {
295  Double_t p_test = 0.001+k*0.01;
296 
297  for ( Int_t l=start2; l<=end2; l++ )
298  {
299  Double_t res_test = 2.0; // 0.001+l*1.0;
300 
301  Double_t fv = my_mcs_llhd( p_test, res_test );
302 
303  if ( fv<logL )
304  {
305  bf = p_test;
306 
307  logL = fv;
308 
309  // errs = res_test;
310 
311  }
312 
313  }
314 
315  }
316 
317  p_mcs_2 = bf; LLbf = logL;
318 
319  p = p_mcs_2;
320 
321  return p;
322 
323  }
324 
325  TVector3 TrackMomentumCalculator::GetMultiScatterStartingPoint( const art::Ptr<recob::Track> &trk )
326  {
327  TVector3 start;
328 
329  Double_t LLHDp = GetMuMultiScatterLLHD3( trk, true );
330 
331  Double_t LLHDm = GetMuMultiScatterLLHD3( trk, false );
332 
333  if ( LLHDp!=-1 && LLHDm!=-1 )
334  {
335  if ( LLHDp>LLHDm )
336  {
337  Int_t n_points = trk->NumberTrajectoryPoints();
338 
339  const TVector3 &pos = trk->LocationAtPoint( n_points-1 );
340 
341  start.SetXYZ( pos.X(), pos.Y(), pos.Z() );
342 
343  }
344  else
345  {
346  const TVector3 &pos = trk->LocationAtPoint( 0 );
347 
348  start.SetXYZ( pos.X(), pos.Y(), pos.Z() );
349 
350  }
351 
352  }
353 
354  else
355  {
356  const TVector3 &pos = trk->LocationAtPoint( 0 );
357 
358  start.SetXYZ( pos.X(), pos.Y(), pos.Z() );
359 
360  }
361 
362  return start;
363 
364  }
365  /*
366  Double_t TrackMomentumCalculator::GetMuMultiScatterLLHD( const art::Ptr<recob::Track> &trk )
367  {
368  Double_t LLHD = -1.0;
369 
370  std::vector<Float_t> recoX; std::vector<Float_t> recoY; std::vector<Float_t> recoZ;
371 
372  recoX.clear(); recoY.clear(); recoZ.clear();
373 
374  Int_t n_points = trk->NumberTrajectoryPoints();
375 
376  for ( Int_t i=0; i<n_points; i++ )
377  {
378  const TVector3 &pos = trk->LocationAtPoint( i );
379 
380  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
381 
382  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl;
383 
384  }
385 
386  Int_t my_steps = recoX.size();
387 
388  if ( my_steps<2 ) return -1.0;
389 
390  Int_t check0 = GetRecoTracks( recoX, recoY, recoZ );
391 
392  if ( check0!=0 ) return -1.0;
393 
394  seg_size = 5.0;
395 
396  Int_t check1 = GetSegTracks2( recoX, recoY, recoZ );
397 
398  if ( check1!=0 ) return -1.0;
399 
400  Int_t seg_steps = segx.size();
401 
402  if ( seg_steps<2 ) return -1;
403 
404  Int_t seg_steps0 = seg_steps-1;
405 
406  Double_t recoL = segL.at(seg_steps0);
407 
408  if ( recoL<15.0 || recoL>1350.0 ) return -1;
409 
410  Int_t check2 = GetDeltaThetaij( dEi, dEj, dthij, seg_size, ind );
411 
412  if ( check2!=0 ) return -1.0;
413 
414  Double_t p_range = recoL*kcal;
415 
416  Double_t logL = my_mcs_llhd( p_range, 0.5 );
417 
418  LLHD = logL;
419 
420  return LLHD;
421 
422  }
423 
424  Double_t TrackMomentumCalculator::GetMuMultiScatterLLHD2( const recob::Track &trk )
425  {
426  Double_t LLHD = -1.0;
427 
428  std::vector<Float_t> recoX; std::vector<Float_t> recoY; std::vector<Float_t> recoZ;
429 
430  recoX.clear(); recoY.clear(); recoZ.clear();
431 
432  Int_t n_points = trk.NumberTrajectoryPoints();
433 
434  for ( Int_t i=0; i<n_points; i++ )
435  {
436  const TVector3 &pos = trk.LocationAtPoint( i );
437 
438  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
439 
440  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl;
441 
442  }
443 
444  Int_t my_steps = recoX.size();
445 
446  if ( my_steps<2 ) return -1.0;
447 
448  Int_t check0 = GetRecoTracks( recoX, recoY, recoZ );
449 
450  if ( check0!=0 ) return -1.0;
451 
452  seg_size = 5.0;
453 
454  Int_t check1 = GetSegTracks2( recoX, recoY, recoZ );
455 
456  if ( check1!=0 ) return -1.0;
457 
458  Int_t seg_steps = segx.size();
459 
460  if ( seg_steps<2 ) return -1;
461 
462  Int_t seg_steps0 = seg_steps-1;
463 
464  Double_t recoL = segL.at(seg_steps0);
465 
466  if ( recoL<15.0 || recoL>1350.0 ) return -1;
467 
468  Int_t check2 = GetDeltaThetaij( dEi, dEj, dthij, seg_size, ind );
469 
470  if ( check2!=0 ) return -1.0;
471 
472  Double_t p_range = recoL*kcal;
473 
474  Double_t logL = my_mcs_llhd( p_range, 0.5 );
475 
476  LLHD = logL;
477 
478  return LLHD;
479 
480  }
481  */
482  Double_t TrackMomentumCalculator::GetMuMultiScatterLLHD3( const art::Ptr<recob::Track> &trk , bool dir )
483  {
484  Double_t LLHD = -1.0;
485 
486  std::vector<Float_t> recoX; std::vector<Float_t> recoY; std::vector<Float_t> recoZ;
487 
488  recoX.clear(); recoY.clear(); recoZ.clear();
489 
490  Int_t n_points = trk->NumberTrajectoryPoints();
491 
492  if ( dir )
493  {
494  for ( Int_t i=0; i<n_points; i++ )
495  {
496  const TVector3 &pos = trk->LocationAtPoint( i );
497 
498  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
499 
500  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl;
501 
502  }
503 
504  }
505  else
506  {
507  for ( Int_t i=0; i<n_points; i++ )
508  {
509  const TVector3 &pos = trk->LocationAtPoint( n_points-1-i );
510 
511  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
512 
513  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl;
514 
515  }
516 
517  }
518 
519  Int_t my_steps = recoX.size();
520 
521  if ( my_steps<2 ) return -1.0;
522 
523  Int_t check0 = GetRecoTracks( recoX, recoY, recoZ );
524 
525  if ( check0!=0 ) return -1.0;
526 
527  seg_size = 5.0;
528 
529  Int_t check1 = GetSegTracks2( recoX, recoY, recoZ );
530 
531  if ( check1!=0 ) return -1.0;
532 
533  Int_t seg_steps = segx.size();
534 
535  if ( seg_steps<2 ) return -1;
536 
537  Int_t seg_steps0 = seg_steps-1;
538 
539  Double_t recoL = segL.at(seg_steps0);
540 
541  if ( recoL<15.0 || recoL>maxLength ) return -1;
542 
543  Int_t check2 = GetDeltaThetaij( dEi, dEj, dthij, seg_size, ind );
544 
545  if ( check2!=0 ) return -1.0;
546 
547  Double_t p_range = recoL*kcal;
548 
549  Double_t logL = my_mcs_llhd( p_range, 5.65 );
550 
551  LLHD = logL;
552 
553  return LLHD;
554 
555  }
556 
557  Int_t TrackMomentumCalculator::GetDeltaThetaij( std::vector<Float_t> &ei, std::vector<Float_t> &ej, std::vector<Float_t> &th, Double_t thick, std::vector<Float_t> &ind )
558  {
559  Int_t a1 = segx.size(); Int_t a2 = segy.size(); Int_t a3 = segz.size();
560 
561  if ( ( a1!=a2 ) || ( a1!=a3 ) ) { cout << " ( Get thij ) Error ! " << endl; return -1.0; }
562 
563  Int_t tot = a1-1; Double_t thick1 = thick+0.13;
564 
565  ei.clear(); ej.clear(); th.clear(); ind.clear();
566 
567  azx0.clear(); azy0.clear();
568 
569  for ( Int_t i=0; i<tot; i++ )
570  {
571  Double_t dx = segnx.at( i ); Double_t dy = segny.at( i ); Double_t dz = segnz.at( i );
572 
573  TVector3 vec_z( dx, dy, dz );
574 
575  TVector3 vec_x;
576 
577  TVector3 vec_y;
578 
579  Double_t switcher = basex.Dot( vec_z );
580 
581  if ( TMath::Abs( switcher )<=0.995 )
582  {
583  vec_y = vec_z.Cross( basex ); vec_y = vec_y.Unit();
584 
585  vec_x = vec_y.Cross( vec_z );
586 
587  }
588 
589  else
590  {
591  // cout << " It switched ! Isn't this lovely !!! " << endl; // getchar();
592 
593  vec_y = basez.Cross( vec_z ); vec_y = vec_y.Unit();
594 
595  vec_x = vec_y.Cross( vec_z );
596 
597  }
598 
599  TVector3 Rx;
600 
601  Double_t ex = vec_x.Dot( basex ); Rx.SetX( ex );
602 
603  ex = vec_x.Dot( basey ); Rx.SetY( ex );
604 
605  ex = vec_x.Dot( basez ); Rx.SetZ( ex );
606 
607  TVector3 Ry;
608 
609  Double_t ey = vec_y.Dot( basex ); Ry.SetX( ey );
610 
611  ey = vec_y.Dot( basey ); Ry.SetY( ey );
612 
613  ey = vec_y.Dot( basez ); Ry.SetZ( ey );
614 
615  TVector3 Rz;
616 
617  Double_t ez = vec_z.Dot( basex ); Rz.SetX( ez );
618 
619  ez = vec_z.Dot( basey ); Rz.SetY( ez );
620 
621  ez = vec_z.Dot( basez ); Rz.SetZ( ez );
622 
623  Double_t refL = segL.at( i );
624 
625  for ( Int_t j=i; j<tot; j++ )
626  {
627  Double_t L1 = segL.at( j );
628 
629  Double_t L2 = segL.at( j+1 );
630 
631  Double_t dz1 = L1 - refL;
632 
633  Double_t dz2 = L2 - refL;
634 
635  if ( dz1<=thick1 && dz2>thick1 )
636  {
637  Double_t here_dx = segnx.at( j );
638 
639  Double_t here_dy = segny.at( j );
640 
641  Double_t here_dz = segnz.at( j );
642 
643  Double_t ac1 = 0.001*180.0/3.14159*find_angle( here_dz, here_dx );
644 
645  Double_t ac2 = 0.001*180.0/3.14159*find_angle( here_dz, here_dy );
646 
647  TVector3 here_vec; here_vec.SetXYZ( here_dx, here_dy, here_dz );
648 
649  TVector3 rot_here; rot_here.SetXYZ( Rx.Dot( here_vec ), Ry.Dot( here_vec ), Rz.Dot( here_vec ) );
650 
651  Double_t scx = rot_here.X();
652 
653  Double_t scy = rot_here.Y();
654 
655  Double_t scz = rot_here.Z();
656 
657  Double_t azy = find_angle( scz, scy );
658 
659  Double_t azx = find_angle( scz, scx );
660 
661  Double_t ULim = 10000.0;
662 
663  Double_t LLim = -10000.0;
664 
665  Double_t cL = kcal;
666 
667  Double_t Li = segL.at( i );
668 
669  Double_t Lj = segL.at( j );
670 
671  if ( azy<=ULim && azy>=LLim )
672  {
673  ei.push_back( Li*cL );
674 
675  ej.push_back( Lj*cL );
676 
677  th.push_back( azy );
678 
679  ind.push_back( 2 );
680 
681  azx0.push_back( ac1 );
682 
683  azy0.push_back( ac2 );
684 
685  }
686 
687  if ( azx<=ULim && azx>=LLim )
688  {
689  ei.push_back( Li*cL );
690 
691  ej.push_back( Lj*cL );
692 
693  th.push_back( azx );
694 
695  ind.push_back( 1 );
696 
697  azx0.push_back( ac1 );
698 
699  azy0.push_back( ac2 );
700 
701  }
702 
703  break; // of course !
704 
705  }
706 
707  }
708 
709  }
710 
711  return 0;
712 
713  }
714 
715  Double_t TrackMomentumCalculator::GetMomentumMultiScatterChi2( const art::Ptr<recob::Track> &trk )
716  {
717  Double_t p = -1.0;
718 
719  std::vector<Float_t> recoX; std::vector<Float_t> recoY; std::vector<Float_t> recoZ;
720 
721  recoX.clear(); recoY.clear(); recoZ.clear();
722 
723  Int_t n_points = trk->NumberTrajectoryPoints();
724 
725  for ( Int_t i=0; i<n_points; i++ )
726  {
727  const TVector3 &pos = trk->LocationAtPoint( i );
728 
729  recoX.push_back( pos.X() ); recoY.push_back( pos.Y() ); recoZ.push_back( pos.Z() );
730 
731  // cout << " posX, Y, Z : " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << endl;
732 
733  }
734 
735  Int_t my_steps = recoX.size();
736 
737  if ( my_steps<2 ) return -1.0;
738 
739  Int_t check0 = GetRecoTracks( recoX, recoY, recoZ );
740 
741  if ( check0!=0 ) return -1.0;
742 
743  seg_size = steps_size;
744 
745  Int_t check1 = GetSegTracks2( recoX, recoY, recoZ );
746 
747  if ( check1!=0 ) return -1.0;
748 
749  Int_t seg_steps = segx.size();
750 
751  if ( seg_steps<2 ) return -1;
752 
753  Int_t seg_steps0 = seg_steps-1;
754 
755  Double_t recoL = segL.at(seg_steps0);
756 
757  if ( seg_steps<2 || recoL<minLength || recoL>maxLength ) return -1;
758 
759  Double_t mean = 666.0; Double_t rms = 666.0; Double_t rmse = 666.0;
760 
761  nmeas = 0; Double_t max1=-999.0; Double_t min1=+999.0;
762 
763  for ( Int_t j=0; j<n_steps; j++ )
764  {
765  Double_t trial = steps.at( j );
766 
767  GetDeltaThetaRMS( mean, rms, rmse, trial );
768 
769  // cout << mean << ", " << rms << ", " << rmse << ", " << trial << endl;
770 
771  xmeas[ nmeas ] = trial;
772 
773  ymeas[ nmeas ] = rms;
774 
775  eymeas[ nmeas ] = sqrt( pow( rmse, 2.0 ) + pow( 0.05*rms, 2.0 ) ); // <--- conservative syst. error to fix chi^{2} behaviour !!!
776 
777  // ymeas[ nmeas ] = 10.0; eymeas[ nmeas ] = 1.0; // <--- for debugging !
778 
779  if ( min1>ymeas[ nmeas ] ) min1=ymeas[ nmeas ];
780 
781  if ( max1<ymeas[ nmeas ] ) max1=ymeas[ nmeas ];
782 
783  nmeas++;
784 
785  }
786 
787  gr_meas = new TGraphErrors( nmeas, xmeas, ymeas, 0, eymeas );
788 
789  gr_meas->SetTitle( "(#Delta#theta)_{rms} versus material thickness; Material thickness in cm; (#Delta#theta)_{rms} in mrad" );
790 
791  gr_meas->SetLineColor( kBlack ); gr_meas->SetMarkerColor( kBlack ); gr_meas->SetMarkerStyle( 20 ); gr_meas->SetMarkerSize( 1.2 );
792 
793  gr_meas->GetXaxis()->SetLimits( ( steps.at( 0 )-steps.at( 0 ) ), ( steps.at( n_steps-1 )+steps.at( 0 ) ) );
794 
795  gr_meas->SetMinimum( 0.0 );
796 
797  gr_meas->SetMaximum( 1.80*max1 );
798 
799  // c1->cd();
800 
801  // gr_meas->Draw( "APEZ" );
802 
803  // c1->Update();
804 
805  // c1->WaitPrimitive();
806 
807  ROOT::Minuit2::Minuit2Minimizer *mP = new ROOT::Minuit2::Minuit2Minimizer( );
808 
809  ROOT::Math::Functor FCA( &my_mcs_chi2, 2 );
810 
811  mP->SetFunction( FCA );
812 
813  mP->SetLimitedVariable( 0, "p_{MCS}", 1.0, 0.01, 0.001, 7.5 );
814 
815  mP->SetLimitedVariable( 1, "#delta#theta", 0.0, 1.0, 0.0, 45.0 );
816 
817  mP->SetMaxFunctionCalls( 1.E9 );
818 
819  mP->SetMaxIterations( 1.E9 );
820 
821  mP->SetTolerance( 0.01 );
822 
823  mP->SetStrategy( 2 );
824 
825  mP->SetErrorDef( 1.0 );
826 
827  bool mstatus = mP->Minimize();
828 
829  mP->Hesse();
830 
831  const double *pars = mP->X();
832 
833  const double *erpars = mP->Errors();
834 
835  Double_t deltap = ( recoL*kcal )/2.0;
836 
837  p_mcs = pars[0]+deltap;
838 
839  p_mcs_e = erpars[0];
840 
841  if ( mstatus ) p = p_mcs;
842 
843  else p = -1.0;
844 
845  chi2 = mP->MinValue();
846 
847  delete mP;
848 
849  return p;
850 
851  }
852 
853  Int_t TrackMomentumCalculator::GetTracks( const std::vector<Float_t> &xxx, const std::vector<Float_t> &yyy, const std::vector<Float_t> &zzz )
854  {
855  Int_t a1 = xxx.size(); Int_t a2 = yyy.size(); Int_t a3 = zzz.size();
856 
857  if ( ( a1!=a2 ) || ( a1!=a3 ) || ( a2!=a3 ) ) { cout << " ( Get tracks ) Error ! " << endl; return -1; }
858 
859  n = 0;
860 
861  for ( Int_t i=0; i<a1; i++ )
862  {
863  Double_t nowx = xxx.at( i );
864 
865  Double_t nowy = yyy.at( i );
866 
867  Double_t nowz = zzz.at( i );
868 
869  x[n] = nowx;
870 
871  y[n] = nowy;
872 
873  z[n] = nowz;
874 
875  n++;
876 
877  }
878 
879  gr_xyz = new TPolyLine3D( n, z, x, y );
880 
881  gr_yz = new TGraph( n, z, y ); gr_xz = new TGraph( n, z, x ); gr_xy = new TGraph( n, x, y );
882 
883  return 0;
884 
885  }
886 
887  Int_t TrackMomentumCalculator::GetRecoTracks( const std::vector<Float_t> &xxx, const std::vector<Float_t> &yyy, const std::vector<Float_t> &zzz )
888  {
889  Int_t a1 = xxx.size(); Int_t a2 = yyy.size(); Int_t a3 = zzz.size();
890 
891  if ( ( a1!=a2 ) || ( a1!=a3 ) || ( a2!=a3 ) ) { cout << " ( Get reco tacks ) Error ! " << endl; return -1; }
892 
893  n_reco = 0;
894 
895  for ( Int_t i=0; i<a1; i++ )
896  {
897  Double_t nowx = xxx.at( i );
898 
899  Double_t nowy = yyy.at( i );
900 
901  Double_t nowz = zzz.at( i );
902 
903  x_reco[n_reco] = nowx;
904 
905  y_reco[n_reco] = nowy;
906 
907  z_reco[n_reco] = nowz;
908 
909  n_reco++;
910 
911  }
912 
913  delete gr_reco_xyz;
914  delete gr_reco_xy;
915  delete gr_reco_yz;
916  delete gr_reco_xz;
917 
918  gr_reco_xyz = new TPolyLine3D( n_reco, z_reco, x_reco, y_reco );
919 
920  gr_reco_yz = new TGraph( n_reco, z_reco, y_reco ); gr_reco_xz = new TGraph( n_reco, z_reco, x_reco ); gr_reco_xy = new TGraph( n_reco, x_reco, y_reco );
921 
922  return 0;
923 
924  }
925 
926  Int_t TrackMomentumCalculator::GetSegTracks( const std::vector<Float_t> &xxx, const std::vector<Float_t> &yyy, const std::vector<Float_t> &zzz )
927  {
928  Int_t a1 = xxx.size(); Int_t a2 = yyy.size(); Int_t a3 = zzz.size();
929 
930  if ( ( a1!=a2 ) || ( a1!=a3 ) || ( a2!=a3 ) ) { cout << " ( Digitize reco tacks ) Error ! " << endl; return -1; }
931 
932  Int_t stopper = seg_stop / seg_size;
933 
934  Int_t a4 = a1-1;
935 
936  segx.clear(); segy.clear(); segz.clear();
937 
938  segnx.clear(); segny.clear(); segnz.clear();
939 
940  segL.clear();
941 
942  Double_t x0 = xxx.at( 0 ); Double_t y0 = yyy.at( 0 ); Double_t z0 = zzz.at( 0 );
943 
944  segx.push_back( x0 ); segy.push_back( y0 ); segz.push_back( z0 ); segL.push_back( 0.0 );
945 
946  n_seg = 0; x_seg[ n_seg ] = x0; y_seg[ n_seg ] = y0; z_seg[ n_seg ] = z0;
947 
948  n_seg++;
949 
950  Int_t ntot = 0;
951 
952  Double_t Sz = 0; Double_t Szz = 0;
953 
954  Double_t Szx = 0; Double_t Sx = 0;
955 
956  Double_t Szy = 0; Double_t Sy = 0;
957 
958  Sz += z0; Szz += z0*z0;
959 
960  Szx += z0*x0; Sx += x0;
961 
962  Szy += z0*y0; Sy += y0;
963 
964  ntot++;
965 
966  for ( Int_t i=1; i<a4; i++ )
967  {
968  Double_t x1 = xxx.at( i ); Double_t y1 = yyy.at( i ); Double_t z1 = zzz.at( i );
969 
970  Double_t dr1 = sqrt( pow( x1-x0, 2 )+pow( y1-y0, 2)+pow( z1-z0, 2 ) );
971 
972  Double_t x2 = xxx.at( i+1 ); Double_t y2 = yyy.at( i+1 ); Double_t z2 = zzz.at( i+1 );
973 
974  Double_t dr2 = sqrt( pow( x2-x0, 2 )+pow( y2-y0, 2)+pow( z2-z0, 2 ) );
975 
976  if ( dr1<seg_size )
977  {
978  Sz += z1; Szz += z1*z1;
979 
980  Szx += z1*x1; Sx += x1;
981 
982  Szy += z1*y1; Sy += y1;
983 
984  ntot++;
985 
986  }
987 
988  if ( dr1<=seg_size && dr2>seg_size )
989  {
990  Double_t t = 0.0;
991 
992  Double_t dx = x2-x1; Double_t dy = y2-y1; Double_t dz = z2-z1;
993 
994  Double_t dr = sqrt( dx*dx+dy*dy+dz*dz );
995 
996  if ( dr==0 ) { cout << " ( Zero ) Error ! " << endl; return -1; }
997 
998  Double_t beta = 2.0*( (x1-x0)*dx+(y1-y0)*dy+(z1-z0)*dz )/( dr*dr );
999 
1000  Double_t gamma = ( dr1*dr1 - seg_size*seg_size )/( dr*dr );
1001 
1002  Double_t delta = beta*beta - 4.0*gamma;
1003 
1004  if ( delta<0.0 ) { cout << " ( Discriminant ) Error ! " << endl; return -1; }
1005 
1006  Double_t lysi1 = ( -beta+sqrt( delta ) )/2.0;
1007 
1008  t=lysi1;
1009 
1010  Double_t xp = x1+t*dx;
1011 
1012  Double_t yp = y1+t*dy;
1013 
1014  Double_t zp = z1+t*dz;
1015 
1016  segx.push_back( xp ); segy.push_back( yp ); segz.push_back( zp );
1017 
1018  segL.push_back( 1.0*n_seg*1.0*seg_size );
1019 
1020  x_seg[ n_seg ] = xp; y_seg[ n_seg ] = yp; z_seg[ n_seg ] = zp; n_seg++;
1021 
1022  x0 = xp; y0 = yp; z0 = zp;
1023 
1024  Sz += z0; Szz += z0*z0;
1025 
1026  Szx += z0*x0; Sx += x0;
1027 
1028  Szy += z0*y0; Sy += y0;
1029 
1030  ntot++;
1031 
1032  Double_t ax = ( Szx*ntot-Sz*Sx )/( Szz*ntot-Sz*Sz );
1033 
1034  Double_t ay = ( Szy*ntot-Sz*Sy )/( Szz*ntot-Sz*Sz );
1035 
1036  segnx.push_back( ax ); segny.push_back( ay ); segnz.push_back( 1.0 );
1037 
1038  ntot = 0;
1039 
1040  Sz = 0; Szz = 0;
1041 
1042  Szx = 0; Sx = 0;
1043 
1044  Szy = 0; Sy = 0;
1045 
1046  Sz += z0; Szz += z0*z0;
1047 
1048  Szx += z0*x0; Sx += x0;
1049 
1050  Szy += z0*y0; Sy += y0;
1051 
1052  ntot++;
1053 
1054  }
1055 
1056  else if ( dr1>seg_size )
1057  {
1058  Double_t t = 0.0;
1059 
1060  Double_t dx = x1-x0; Double_t dy = y1-y0; Double_t dz = z1-z0;
1061 
1062  Double_t dr = sqrt( dx*dx+dy*dy+dz*dz );
1063 
1064  if ( dr==0 ) { cout << " ( Zero ) Error ! " << endl; getchar(); }
1065 
1066  if ( dr!=0 ) t = seg_size/dr;
1067 
1068  Double_t xp = x0+t*dx;
1069 
1070  Double_t yp = y0+t*dy;
1071 
1072  Double_t zp = z0+t*dz;
1073 
1074  segx.push_back( xp ); segy.push_back( yp ); segz.push_back( zp );
1075 
1076  segL.push_back( 1.0*n_seg*1.0*seg_size );
1077 
1078  x_seg[ n_seg ] = xp; y_seg[ n_seg ] = yp; z_seg[ n_seg ] = zp; n_seg++;
1079 
1080  x0 = xp; y0 = yp; z0 = zp;
1081 
1082  i=( i-1 );
1083 
1084  Sz += z0; Szz += z0*z0;
1085 
1086  Szx += z0*x0; Sx += x0;
1087 
1088  Szy += z0*y0; Sy += y0;
1089 
1090  ntot++;
1091 
1092  Double_t ax = ( Szx*ntot-Sz*Sx )/( Szz*ntot-Sz*Sz );
1093 
1094  Double_t ay = ( Szy*ntot-Sz*Sy )/( Szz*ntot-Sz*Sz );
1095 
1096  segnx.push_back( ax ); segny.push_back( ay ); segnz.push_back( 1.0 );
1097 
1098  // Double_t angx = find_angle( 1.0, ax ); Double_t angy = find_angle( 1.0, ay );
1099 
1100  // cout << angx*0.001*180.0/3.14 << endl;
1101 
1102  ntot = 0;
1103 
1104  Sz = 0; Szz = 0;
1105 
1106  Szx = 0; Sx = 0;
1107 
1108  Szy = 0; Sy = 0;
1109 
1110  Sz += z0; Szz += z0*z0;
1111 
1112  Szx += z0*x0; Sx += x0;
1113 
1114  Szy += z0*y0; Sy += y0;
1115 
1116  ntot++;
1117 
1118  }
1119 
1120  if ( n_seg>=( stopper+1.0 ) && seg_stop!=-1 ) break;
1121 
1122  }
1123 
1124  delete gr_seg_xyz;
1125  delete gr_seg_xy;
1126  delete gr_seg_yz;
1127  delete gr_seg_xz;
1128 
1129  gr_seg_xyz = new TPolyLine3D( n_seg, z_seg, x_seg, y_seg );
1130 
1131  gr_seg_yz = new TGraph( n_seg, z_seg, y_seg ); gr_seg_xz = new TGraph( n_seg, z_seg, x_seg ); gr_seg_xy = new TGraph( n_seg, x_seg, y_seg );
1132 
1133  return 0;
1134 
1135  }
1136 
1137  Int_t TrackMomentumCalculator::GetSegTracks2( const std::vector<Float_t> &xxx, const std::vector<Float_t> &yyy, const std::vector<Float_t> &zzz )
1138  {
1139  Double_t stag = 0.0;
1140 
1141  Int_t a1 = xxx.size(); Int_t a2 = yyy.size(); Int_t a3 = zzz.size();
1142 
1143  if ( ( a1!=a2 ) || ( a1!=a3 ) || ( a2!=a3 ) ) { cout << " ( Digitize reco tacks ) Error ! " << endl; return -1; }
1144 
1145  Int_t stopper = seg_stop / seg_size;
1146 
1147  Int_t a4 = a1-1;
1148 
1149  segx.clear(); segy.clear(); segz.clear();
1150 
1151  segnx.clear(); segny.clear(); segnz.clear();
1152 
1153  segL.clear();
1154 
1155  Int_t ntot = 0;
1156 
1157  n_seg = 0;
1158 
1159  Double_t x0=0.0; Double_t y0=0.0; Double_t z0=0.0;
1160 
1161  Double_t x00 = xxx.at( 0 ); Double_t y00 = yyy.at( 0 ); Double_t z00 = zzz.at( 0 );
1162 
1163  Int_t indC=0;
1164 
1165  std::vector<Float_t> vx; std::vector<Float_t> vy; std::vector<Float_t> vz;
1166 
1167  vx.clear(); vy.clear(); vz.clear();
1168 
1169  for ( Int_t i=0; i<=a4; i++ )
1170  {
1171  x0 = xxx.at( i ); y0 = yyy.at( i ); z0 = zzz.at( i );
1172 
1173  Double_t RR0 = sqrt( pow(x00-x0,2.0)+pow(y00-y0,2.0)+pow(z00-z0,2.0) );
1174 
1175  if ( RR0>=stag )
1176  {
1177  segx.push_back( x0 ); segy.push_back( y0 ); segz.push_back( z0 );
1178 
1179  segL.push_back( stag );
1180 
1181  x_seg[ n_seg ] = x0; y_seg[ n_seg ] = y0; z_seg[ n_seg ] = z0;
1182 
1183  n_seg++;
1184 
1185  vx.push_back( x0 ); vy.push_back( y0 ); vz.push_back( z0 );
1186 
1187  ntot++;
1188 
1189  indC = i+1;
1190 
1191  break;
1192 
1193  }
1194 
1195  }
1196 
1197  for ( Int_t i=indC; i<a4; i++ )
1198  {
1199  Double_t x1 = xxx.at( i ); Double_t y1 = yyy.at( i ); Double_t z1 = zzz.at( i );
1200 
1201  Double_t dr1 = sqrt( pow( x1-x0, 2 )+pow( y1-y0, 2)+pow( z1-z0, 2 ) );
1202 
1203  Double_t x2 = xxx.at( i+1 ); Double_t y2 = yyy.at( i+1 ); Double_t z2 = zzz.at( i+1 );
1204 
1205  Double_t dr2 = sqrt( pow( x2-x0, 2 )+pow( y2-y0, 2)+pow( z2-z0, 2 ) );
1206 
1207  if ( dr1<seg_size )
1208  {
1209  vx.push_back( x1 ); vy.push_back( y1 ); vz.push_back( z1 );
1210 
1211  ntot++;
1212 
1213  }
1214 
1215  if ( dr1<=seg_size && dr2>seg_size )
1216  {
1217  // ..
1218 
1219  // cout << " 1 " << endl;
1220 
1221  Double_t t = 0.0;
1222 
1223  Double_t dx = x2-x1; Double_t dy = y2-y1; Double_t dz = z2-z1;
1224 
1225  Double_t dr = sqrt( dx*dx+dy*dy+dz*dz );
1226 
1227  if ( dr==0 ) { cout << " ( Zero ) Error ! " << endl; return -1; }
1228 
1229  Double_t beta = 2.0*( (x1-x0)*dx+(y1-y0)*dy+(z1-z0)*dz )/( dr*dr );
1230 
1231  Double_t gamma = ( dr1*dr1 - seg_size*seg_size )/( dr*dr );
1232 
1233  Double_t delta = beta*beta - 4.0*gamma;
1234 
1235  if ( delta<0.0 ) { cout << " ( Discriminant ) Error ! " << endl; return -1; }
1236 
1237  Double_t lysi1 = ( -beta+sqrt( delta ) )/2.0;
1238 
1239  t=lysi1;
1240 
1241  Double_t xp = x1+t*dx;
1242 
1243  Double_t yp = y1+t*dy;
1244 
1245  Double_t zp = z1+t*dz;
1246 
1247  segx.push_back( xp ); segy.push_back( yp ); segz.push_back( zp );
1248 
1249  segL.push_back( 1.0*n_seg*1.0*seg_size+stag );
1250 
1251  x_seg[ n_seg ] = xp; y_seg[ n_seg ] = yp; z_seg[ n_seg ] = zp; n_seg++;
1252 
1253  x0 = xp; y0 = yp; z0 = zp;
1254 
1255  vx.push_back( x0 ); vy.push_back( y0 ); vz.push_back( z0 );
1256 
1257  ntot++;
1258 
1259  Double_t na = vx.size();
1260 
1261  Double_t sumx = 0.0;
1262 
1263  Double_t sumy = 0.0;
1264 
1265  Double_t sumz = 0.0;
1266 
1267  for ( Int_t ii=0; ii<na; ii++ )
1268  {
1269  Double_t xxw1 = vx.at( ii );
1270 
1271  Double_t yyw1 = vy.at( ii );
1272 
1273  Double_t zzw1 = vz.at( ii );
1274 
1275  sumx += xxw1; sumy += yyw1; sumz += zzw1;
1276 
1277  }
1278 
1279  sumx = sumx/na; sumy = sumy/na; sumz = sumz/na;
1280 
1281  std::vector<Double_t> mx;
1282 
1283  std::vector<Double_t> my;
1284 
1285  std::vector<Double_t> mz;
1286 
1287  TMatrixDSym m( 3 );
1288 
1289  for ( Int_t ii=0; ii<na; ii++ )
1290  {
1291  Double_t xxw1 = vx.at( ii ); Double_t yyw1 = vy.at( ii ); Double_t zzw1 = vz.at( ii );
1292 
1293  mx.push_back( xxw1-sumx ); my.push_back( yyw1-sumy ); mz.push_back( zzw1-sumz );
1294 
1295  Double_t xxw0 = mx.at( ii ); Double_t yyw0 = my.at( ii ); Double_t zzw0 = mz.at( ii );
1296 
1297  m( 0, 0 ) += xxw0*xxw0/na; m( 0, 1 ) += xxw0*yyw0/na; m( 0, 2 ) += xxw0*zzw0/na;
1298 
1299  m( 1, 0 ) += yyw0*xxw0/na; m( 1, 1 ) += yyw0*yyw0/na; m( 1, 2 ) += yyw0*zzw0/na;
1300 
1301  m( 2, 0 ) += zzw0*xxw0/na; m( 2, 1 ) += zzw0*yyw0/na; m( 2, 2 ) += zzw0*zzw0/na;
1302 
1303  }
1304 
1305  TMatrixDSymEigen me(m);
1306 
1307  TVectorD eigenval = me.GetEigenValues();
1308 
1309  TMatrixD eigenvec = me.GetEigenVectors();
1310 
1311  Double_t max1 = -666.0;
1312 
1313  Double_t ind1 = 0;
1314 
1315  for ( Int_t ii=0; ii<3; ii++ )
1316  {
1317  Double_t p1 = eigenval( ii );
1318 
1319  if ( p1>max1 ) { max1=p1; ind1=ii; }
1320 
1321  }
1322 
1323  // cout << ind1 << endl;
1324 
1325  Double_t ax = eigenvec( 0, ind1 );
1326 
1327  Double_t ay = eigenvec( 1, ind1 );
1328 
1329  Double_t az = eigenvec( 2, ind1 );
1330 
1331  if ( n_seg>1 )
1332  {
1333  if ( segx.at(n_seg-1)-segx.at(n_seg-2) > 0 ) ax = TMath::Abs( ax );
1334 
1335  else ax = -1.0*TMath::Abs( ax );
1336 
1337  if ( segy.at(n_seg-1)-segy.at(n_seg-2) > 0 ) ay = TMath::Abs( ay );
1338 
1339  else ay = -1.0*TMath::Abs( ay );
1340 
1341  if ( segz.at(n_seg-1)-segz.at(n_seg-2) > 0 ) az = TMath::Abs( az );
1342 
1343  else az = -1.0*TMath::Abs( az );
1344 
1345  segnx.push_back( ax ); segny.push_back( ay ); segnz.push_back( az );
1346 
1347  // Double_t angx = find_angle( 1.0, ax ); Double_t angy = find_angle( 1.0, ay );
1348 
1349  // cout << angx*0.001*180.0/3.14 << endl;
1350 
1351  }
1352 
1353  else return -1;
1354 
1355  ntot = 0;
1356 
1357  vx.clear(); vy.clear(); vz.clear();
1358 
1359  vx.push_back( x0 ); vy.push_back( y0 ); vz.push_back( z0 );
1360 
1361  ntot++;
1362 
1363  }
1364 
1365  else if ( dr1>seg_size )
1366  {
1367  // ..
1368 
1369  // cout << " 2 " << endl;
1370 
1371  Double_t t = 0.0;
1372 
1373  Double_t dx = x1-x0; Double_t dy = y1-y0; Double_t dz = z1-z0;
1374 
1375  Double_t dr = sqrt( dx*dx+dy*dy+dz*dz );
1376 
1377  if ( dr==0 ) { cout << " ( Zero ) Error ! " << endl; return -1; }
1378 
1379  if ( dr!=0 ) t = seg_size/dr;
1380 
1381  Double_t xp = x0+t*dx;
1382 
1383  Double_t yp = y0+t*dy;
1384 
1385  Double_t zp = z0+t*dz;
1386 
1387  segx.push_back( xp ); segy.push_back( yp ); segz.push_back( zp );
1388 
1389  segL.push_back( 1.0*n_seg*1.0*seg_size+stag );
1390 
1391  x_seg[ n_seg ] = xp; y_seg[ n_seg ] = yp; z_seg[ n_seg ] = zp; n_seg++;
1392 
1393  x0 = xp; y0 = yp; z0 = zp;
1394 
1395  i=( i-1 );
1396 
1397  // ..
1398 
1399  vx.push_back( x0 ); vy.push_back( y0 ); vz.push_back( z0 );
1400 
1401  ntot++;
1402 
1403  Double_t na = vx.size();
1404 
1405  Double_t sumx = 0.0;
1406 
1407  Double_t sumy = 0.0;
1408 
1409  Double_t sumz = 0.0;
1410 
1411  for ( Int_t ii=0; ii<na; ii++ )
1412  {
1413  Double_t xxw1 = vx.at( ii );
1414 
1415  Double_t yyw1 = vy.at( ii );
1416 
1417  Double_t zzw1 = vz.at( ii );
1418 
1419  sumx += xxw1; sumy += yyw1; sumz += zzw1;
1420 
1421  }
1422 
1423  sumx = sumx/na; sumy = sumy/na; sumz = sumz/na;
1424 
1425  std::vector<Double_t> mx;
1426 
1427  std::vector<Double_t> my;
1428 
1429  std::vector<Double_t> mz;
1430 
1431  TMatrixDSym m( 3 );
1432 
1433  for ( Int_t ii=0; ii<na; ii++ )
1434  {
1435  Double_t xxw1 = vx.at( ii ); Double_t yyw1 = vy.at( ii ); Double_t zzw1 = vz.at( ii );
1436 
1437  mx.push_back( xxw1-sumx ); my.push_back( yyw1-sumy ); mz.push_back( zzw1-sumz );
1438 
1439  Double_t xxw0 = mx.at( ii ); Double_t yyw0 = my.at( ii ); Double_t zzw0 = mz.at( ii );
1440 
1441  m( 0, 0 ) += xxw0*xxw0/na; m( 0, 1 ) += xxw0*yyw0/na; m( 0, 2 ) += xxw0*zzw0/na;
1442 
1443  m( 1, 0 ) += yyw0*xxw0/na; m( 1, 1 ) += yyw0*yyw0/na; m( 1, 2 ) += yyw0*zzw0/na;
1444 
1445  m( 2, 0 ) += zzw0*xxw0/na; m( 2, 1 ) += zzw0*yyw0/na; m( 2, 2 ) += zzw0*zzw0/na;
1446 
1447  }
1448 
1449  TMatrixDSymEigen me(m);
1450 
1451  TVectorD eigenval = me.GetEigenValues();
1452 
1453  TMatrixD eigenvec = me.GetEigenVectors();
1454 
1455  Double_t max1 = -666.0;
1456 
1457  Double_t ind1 = 0;
1458 
1459  for ( Int_t ii=0; ii<3; ii++ )
1460  {
1461  Double_t p1 = eigenval( ii );
1462 
1463  if ( p1>max1 ) { max1=p1; ind1=ii; }
1464 
1465  }
1466 
1467  // cout << ind1 << endl;
1468 
1469  Double_t ax = eigenvec( 0, ind1 );
1470 
1471  Double_t ay = eigenvec( 1, ind1 );
1472 
1473  Double_t az = eigenvec( 2, ind1 );
1474 
1475  if ( n_seg>1 )
1476  {
1477 
1478  if ( segx.at(n_seg-1)-segx.at(n_seg-2) > 0 ) ax = TMath::Abs( ax );
1479 
1480  else ax = -1.0*TMath::Abs( ax );
1481 
1482  if ( segy.at(n_seg-1)-segy.at(n_seg-2) > 0 ) ay = TMath::Abs( ay );
1483 
1484  else ay = -1.0*TMath::Abs( ay );
1485 
1486  if ( segz.at(n_seg-1)-segz.at(n_seg-2) > 0 ) az = TMath::Abs( az );
1487 
1488  else az = -1.0*TMath::Abs( az );
1489 
1490  segnx.push_back( ax ); segny.push_back( ay ); segnz.push_back( az );
1491 
1492  // Double_t angx = find_angle( 1.0, ax ); Double_t angy = find_angle( 1.0, ay );
1493 
1494  // cout << angx*0.001*180.0/3.14 << endl;
1495 
1496  }
1497 
1498  else return -1;
1499 
1500  ntot = 0;
1501 
1502  vx.clear(); vy.clear(); vz.clear();
1503 
1504  vx.push_back( x0 ); vy.push_back( y0 ); vz.push_back( z0 );
1505 
1506  ntot++;
1507 
1508  }
1509 
1510  if ( n_seg>=( stopper+1.0 ) && seg_stop!=-1 ) break;
1511 
1512  }
1513 
1514  delete gr_seg_xyz;
1515  delete gr_seg_xy;
1516  delete gr_seg_yz;
1517  delete gr_seg_xz;
1518 
1519  gr_seg_xyz = new TPolyLine3D( n_seg, z_seg, x_seg, y_seg );
1520 
1521  gr_seg_yz = new TGraph( n_seg, z_seg, y_seg ); gr_seg_xz = new TGraph( n_seg, z_seg, x_seg ); gr_seg_xy = new TGraph( n_seg, x_seg, y_seg );
1522 
1523  return 0;
1524 
1525  }
1526 
1527  void TrackMomentumCalculator::GetDeltaThetaRMS( Double_t &mean, Double_t &rms, Double_t &rmse, Double_t thick )
1528  {
1529  mean = 0.0; rms = 0.0; rmse = 0.0;
1530 
1531  Int_t a1 = segx.size(); Int_t a2 = segy.size(); Int_t a3 = segz.size();
1532 
1533  if ( ( a1!=a2 ) || ( a1!=a3 ) ) { cout << " ( Get RMS ) Error ! " << endl; return; }
1534 
1535  Int_t tot = a1-1;
1536 
1537  Double_t thick1 = thick+0.13;
1538 
1539  vector<Float_t> buf0; buf0.clear();
1540 
1541  for ( Int_t i=0; i<tot; i++ )
1542  {
1543  Double_t dx = segnx.at( i ); Double_t dy = segny.at( i ); Double_t dz = segnz.at( i );
1544 
1545  TVector3 vec_z( dx, dy, dz );
1546 
1547  TVector3 vec_x;
1548 
1549  TVector3 vec_y;
1550 
1551  Double_t switcher = basex.Dot( vec_z );
1552 
1553  if ( switcher<=0.995 )
1554  {
1555  vec_y = vec_z.Cross( basex ); vec_y = vec_y.Unit();
1556 
1557  vec_x = vec_y.Cross( vec_z );
1558 
1559  }
1560 
1561  else
1562  {
1563  // cout << " It switched ! Isn't this lovely !!! " << endl; // getchar();
1564 
1565  vec_y = basez.Cross( vec_z ); vec_y = vec_y.Unit();
1566 
1567  vec_x = vec_y.Cross( vec_z );
1568 
1569  }
1570 
1571  TVector3 Rx;
1572 
1573  Double_t ex = vec_x.Dot( basex ); Rx.SetX( ex );
1574 
1575  ex = vec_x.Dot( basey ); Rx.SetY( ex );
1576 
1577  ex = vec_x.Dot( basez ); Rx.SetZ( ex );
1578 
1579  TVector3 Ry;
1580 
1581  Double_t ey = vec_y.Dot( basex ); Ry.SetX( ey );
1582 
1583  ey = vec_y.Dot( basey ); Ry.SetY( ey );
1584 
1585  ey = vec_y.Dot( basez ); Ry.SetZ( ey );
1586 
1587  TVector3 Rz;
1588 
1589  Double_t ez = vec_z.Dot( basex ); Rz.SetX( ez );
1590 
1591  ez = vec_z.Dot( basey ); Rz.SetY( ez );
1592 
1593  ez = vec_z.Dot( basez ); Rz.SetZ( ez );
1594 
1595  Double_t refL = segL.at( i );
1596 
1597  for ( Int_t j=i; j<tot; j++ )
1598  {
1599  Double_t L1 = segL.at( j );
1600 
1601  Double_t L2 = segL.at( j+1 );
1602 
1603  Double_t dz1 = L1 - refL;
1604 
1605  Double_t dz2 = L2 - refL;
1606 
1607  if ( dz1<=thick1 && dz2>thick1 )
1608  {
1609  Double_t here_dx = segnx.at( j );
1610 
1611  Double_t here_dy = segny.at( j );
1612 
1613  Double_t here_dz = segnz.at( j );
1614 
1615  TVector3 here_vec; here_vec.SetXYZ( here_dx, here_dy, here_dz );
1616 
1617  TVector3 rot_here; rot_here.SetXYZ( Rx.Dot( here_vec ), Ry.Dot( here_vec ), Rz.Dot( here_vec ) );
1618 
1619  Double_t scx = rot_here.X();
1620 
1621  Double_t scy = rot_here.Y();
1622 
1623  Double_t scz = rot_here.Z();
1624 
1625  Double_t azy = find_angle( scz, scy ); azy*=1.0;
1626 
1627  Double_t azx = find_angle( scz, scx );
1628 
1629  Double_t ULim = 10000.0; Double_t LLim = -10000.0;
1630 
1631  // if ( azy<=ULim && azy>=LLim ) { buf0.push_back( azy ); } // hRMS.Fill( azy ); }
1632 
1633  if ( azx<=ULim && azx>=LLim ) { buf0.push_back( azx ); } // hRMS.Fill( azx ); }
1634 
1635  // if ( azy<=ULim && azy>=LLim && azx<=ULim && azx>=LLim && thick==5.0 ) { hSCA.Fill( azx, azy ); }
1636 
1637  // i=j-1;
1638 
1639  break; // of course !
1640 
1641  }
1642  }
1643 
1644  }
1645 
1646  Int_t nmeas = buf0.size(); Double_t nnn = 0.0;
1647 
1648  for ( Int_t i=0; i<nmeas; i++ ) { mean += buf0.at( i ); nnn++; } mean = mean/nnn;
1649 
1650  for ( Int_t i=0; i<nmeas; i++ ) rms += ( ( buf0.at( i ) )*( buf0.at( i ) ) );
1651 
1652  rms = rms/( nnn ); rms = sqrt( rms ); rmse = rms/sqrt( 2.0*tot );
1653 
1654  Double_t rms1 = rms;
1655 
1656  rms = 0.0;
1657 
1658  Double_t ntot1 = 0.0;
1659 
1660  Double_t lev1 = 2.50;
1661 
1662  for ( Int_t i=0; i<nmeas; i++ )
1663  {
1664  Double_t amp = buf0.at( i );
1665 
1666  if ( amp<( mean+lev1*rms1 ) && amp>( mean-lev1*rms1 ) )
1667  {
1668  ntot1++;
1669 
1670  rms += ( amp*amp );
1671 
1672  }
1673 
1674  }
1675 
1676  rms = rms/( ntot1 ); rms = sqrt( rms ); rmse = rms/sqrt( 2.0*ntot1 );
1677 
1678  return;
1679 
1680  }
1681 
1682  Double_t TrackMomentumCalculator::find_angle( Double_t vz, Double_t vy )
1683  {
1684  Double_t thetayz = -999.0;
1685 
1686  if ( vz>0 && vy>0 ) { Double_t ratio=TMath::Abs( vy/vz ); thetayz=TMath::ATan( ratio ); }
1687 
1688  else if ( vz<0 && vy>0 ) { Double_t ratio=TMath::Abs( vy/vz ); thetayz=TMath::ATan( ratio ); thetayz=3.14159-thetayz; }
1689 
1690  else if ( vz<0 && vy<0 ) { Double_t ratio=TMath::Abs( vy/vz ); thetayz=TMath::ATan( ratio ); thetayz=thetayz+3.14159; }
1691 
1692  else if ( vz>0 && vy<0 ) { Double_t ratio=TMath::Abs( vy/vz ); thetayz=TMath::ATan( ratio ); thetayz=2.0*3.14159-thetayz; }
1693 
1694  else if ( vz==0 && vy>0 ) { thetayz=3.14159/2.0; }
1695 
1696  else if ( vz==0 && vy<0 ) { thetayz=3.0*3.14159/2.0; }
1697 
1698  if ( thetayz>3.14159 ) { thetayz=thetayz-2.0*3.14159; }
1699 
1700  Double_t result = 1000.0*thetayz;
1701 
1702  return result;
1703 
1704  }
1705 
1706  Double_t TrackMomentumCalculator::my_g( Double_t xx, Double_t Q, Double_t s )
1707  {
1708  Double_t arg = 0.0;
1709 
1710  if ( s!=0 ) arg = ( xx - Q )/s;
1711 
1712  else cout << " Error : The code tries to divide by zero ! " << endl;
1713 
1714  Double_t result = 0.0;
1715 
1716  if ( s!=0 ) result = -0.5*TMath::Log( 2.0*TMath::Pi() ) - TMath::Log( s ) - 0.5*arg*arg;
1717 
1718  if ( isnan( float( result ) ) || isinf( float( result ) ) ) {
1719  cout << " Is nan ! my_g ! " << - TMath::Log( s ) << ", " << s << endl;
1720  getchar();
1721  }
1722 
1723  return result;
1724 
1725  }
1726 
1727  Double_t TrackMomentumCalculator::my_mcs_llhd( Double_t x0, Double_t x1 )
1728  {
1729  Double_t p = x0;
1730 
1731  Double_t theta0x = x1;
1732 
1733  Double_t result = 0.0;
1734 
1735  Double_t nnn1 = dEi.size();
1736 
1737  Double_t red_length = ( 10.0 )/14.0;
1738 
1739  Double_t addth = 0;
1740 
1741  for ( Int_t i=0; i<nnn1; i++ )
1742  {
1743  Double_t Ei = p-dEi.at( i );
1744 
1745  Double_t Ej = p-dEj.at( i );
1746 
1747  if ( Ei>0 && Ej<0 ) addth=3.14*1000.0;
1748 
1749  Ei = TMath::Abs( Ei );
1750 
1751  Ej = TMath::Abs( Ej );
1752 
1753  Double_t tH0 = ( 13.6/sqrt( Ei*Ej ) )*( 1.0+0.038*TMath::Log( red_length ) )*sqrt( red_length );
1754 
1755  Double_t rms = -1.0;
1756 
1757  if ( ind.at( i )==1 )
1758  {
1759  rms = sqrt( tH0*tH0+pow( theta0x, 2.0 ) );
1760 
1761  Double_t DT = dthij.at( i )+addth;
1762 
1763  Double_t prob = my_g( DT, 0.0, rms );
1764 
1765  result = result - 2.0*prob;
1766 
1767  }
1768 
1769  }
1770 
1771  if ( isnan( float( result ) ) || isinf( float( result ) ) ) { cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << endl; getchar(); }
1772 
1773  return result;
1774 
1775  }
1776 
1777 } // namespace track
x
Float_t x
Definition: compare.C:6
s
Float_t s
Definition: plot.C:23
recob::Track::LocationAtPoint
TVector3 LocationAtPoint(unsigned int p) const
Covariance matrices are either set or not.
Definition: Track.h:241
trkf::TrackMomentumCalculator::find_angle
Double_t find_angle(Double_t vz, Double_t vy)
Definition: TrackMomentumCalculator.cxx:1682
trkf::TrackMomentumCalculator::GetSegTracks2
Int_t GetSegTracks2(const std::vector< Float_t > &xxx, const std::vector< Float_t > &yyy, const std::vector< Float_t > &zzz)
Definition: TrackMomentumCalculator.cxx:1137
xx
Double_t xx
Definition: macro.C:12
E
Float_t E
Definition: plot.C:23
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2
Double_t GetMomentumMultiScatterChi2(const art::Ptr< recob::Track > &trk)
Definition: TrackMomentumCalculator.cxx:715
eymeas
Double_t eymeas[30]
Definition: TrackMomentumCalculator.h:45
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
trkf::TrackMomentumCalculator::GetTracks
Int_t GetTracks(const std::vector< Float_t > &xxx, const std::vector< Float_t > &yyy, const std::vector< Float_t > &zzz)
Definition: TrackMomentumCalculator.cxx:853
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:109
beta
Double_t beta
Definition: TG43_relative_dose.C:15
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:214
my_mcs_chi2
double my_mcs_chi2(const double *x)
Definition: TrackMomentumCalculator.cxx:7
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
trkf::TrackMomentumCalculator::GetSegTracks
Int_t GetSegTracks(const std::vector< Float_t > &xxx, const std::vector< Float_t > &yyy, const std::vector< Float_t > &zzz)
Definition: TrackMomentumCalculator.cxx:926
a2
#define a2
Definition: HoughBaseAlg.cxx:67
a3
#define a3
Definition: HoughBaseAlg.cxx:68
trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD
Double_t GetMomentumMultiScatterLLHD(const art::Ptr< recob::Track > &trk)
Definition: TrackMomentumCalculator.cxx:237
trkf::TrackMomentumCalculator::GetDeltaThetaRMS
void GetDeltaThetaRMS(Double_t &mean, Double_t &rms, Double_t &rmse, Double_t thick)
Definition: TrackMomentumCalculator.cxx:1527
nmeas
Int_t nmeas
Definition: TrackMomentumCalculator.h:45
ymeas
Double_t ymeas[30]
Definition: TrackMomentumCalculator.h:45
trkf::TrackMomentumCalculator::my_mcs_llhd
Double_t my_mcs_llhd(Double_t x0, Double_t x1)
Definition: TrackMomentumCalculator.cxx:1727
xmeas
Double_t xmeas[30]
Definition: TrackMomentumCalculator.h:45
trkf::TrackMomentumCalculator::my_g
Double_t my_g(Double_t xx, Double_t Q, Double_t s)
Definition: TrackMomentumCalculator.cxx:1706
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg)
Definition: TrackMomentumCalculator.cxx:117
trkf::TrackMomentumCalculator::GetDeltaThetaij
Int_t GetDeltaThetaij(std::vector< Float_t > &ei, std::vector< Float_t > &ej, std::vector< Float_t > &th, Double_t thick, std::vector< Float_t > &ind)
Definition: TrackMomentumCalculator.cxx:557
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:15
trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint
TVector3 GetMultiScatterStartingPoint(const art::Ptr< recob::Track > &trk)
Definition: TrackMomentumCalculator.cxx:325
dir
TDirectory * dir
Definition: macro.C:5
trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3
Double_t GetMuMultiScatterLLHD3(const art::Ptr< recob::Track > &trk, bool dir)
Definition: TrackMomentumCalculator.cxx:482
LOG_DEBUG
#define LOG_DEBUG(id)
Definition: MessageLogger.h:263
a4
#define a4
Definition: HoughBaseAlg.cxx:69
x2
Float_t x2[n_points_geant4]
Definition: compare.C:26
e
Float_t e
Definition: plot.C:34
trkf::TrackMomentumCalculator::GetRecoTracks
Int_t GetRecoTracks(const std::vector< Float_t > &xxx, const std::vector< Float_t > &yyy, const std::vector< Float_t > &zzz)
Definition: TrackMomentumCalculator.cxx:887
a1
#define a1
Definition: HoughBaseAlg.cxx:66