EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Nystrom.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Nystrom.cxx
1 /* Copyright 2008-2009, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 //-----------------------------------------------------------
20 // File and Version Information:
21 // $Id$
22 //
23 // Description:
24 // Implementation of class Nystrom
25 // see Nystrom.hh for details
26 //
27 // Environment:
28 // Software developed for the PANDA Detector at FAIR.
29 //
30 // Author List:
31 // Sebastian Neubert TUM (original author)
32 //
33 //
34 //-----------------------------------------------------------
35 
36 // Panda Headers ----------------------
37 
38 // This Class' Header ------------------
39 #include "Nystrom.h"
40 #include <cmath>
41 // C/C++ Headers ----------------------
42 #include <iostream>
43 #include "assert.h"
44 #include "TMath.h"
45 
46 // Collaborating Class Headers --------
47 #include "AbsNystromEQM.h"
48 
49 // Class Member definitions -----------
50 
51 
52 Nystrom::Nystrom(AbsNystromEQM* eqm)
53  : _eqm(eqm), _acc(1E-2), _adaptive(true)
54 {}
55 
56 void
57 Nystrom::step(const TVectorT<double>& u, const TVectorT<double>& uprime,
58  const TVectorT<double>& par,
59  TVectorT<double>& newu, TVectorT<double>& newuprime,
60  double h)
61 {
62  assert(h!=0);
63  double sign= h > 0 ? 1. : -1.;
64  double habs=sign*h;
65  double hh=habs*0.5;
66  double h2=h*h;
67  TVectorT<double> k1=_eqm->eval(u,uprime,par);
68  TVectorT<double> k2=_eqm->eval(u+hh*uprime+h2/8.*k1,uprime+hh*k1,par);
69  TVectorT<double> k3=_eqm->eval(u+hh*uprime+h2/8.*k2,uprime+hh*k2,par);
70  TVectorT<double> k4=_eqm->eval(u+habs*uprime+h2/2.*k3,uprime+habs*k3,par);
71 
72  TVectorT<double> ksum=k1+k2+k3;
73  newu=u+h*uprime+h2/6.*ksum;
74  newuprime=uprime+(h/6.)*(ksum+k4+k2+k3);
75  assert(TMath::Abs(newuprime[2]-1.)<1E-4);
76 }
77 
78 
79 // this will return the new h!
80 double
81 Nystrom::adaptiveStep(const TVectorT<double>& u,
82  const TVectorT<double>& uprime,
83  const TVectorT<double>& par,
84  TVectorT<double>& newu,
85  TVectorT<double>& newuprime,
86  double& stepdone,
87  double h, double sign, double reststep)
88 {
89 
90  //std::cout<<"Starting adaptive step:"<<std::endl;
91  double d0=_acc*h; // desired accuracy
92  double d=1.1*d0;
93  double newh=h;
94  double _maxstep=0.1;
95 
96 
97  assert(reststep>=0);
98 
99  TVectorT<double> newu1(u);
100  TVectorT<double> newuprime1(uprime);
101  TVectorT<double> newu2(u);
102  TVectorT<double> newuprime2(uprime);
103  TVectorT<double> utmp(u);
104  TVectorT<double> uprimetmp(uprime);
105  TVectorT<double> delta(u);
106  //newh=0.01;
107  int n=delta.GetNrows();
108  while(d>d0){
109 
110  if(newh>=reststep)newh=reststep;
111 
112  //std::cout<<"newh="<<newh<<" d="<<d<<std::endl;
113  // make one step
114  step(u,uprime,par,
115  newu1,newuprime1,newh*sign);
116 
117 
118  // make two half-steps
119  step(u,uprime,par,
120  utmp,uprimetmp,newh*0.5*sign);
121  step(utmp,uprimetmp,par,
122  newu2,newuprime2,newh*0.5*sign);
123 
124  stepdone=newh*sign;
125 
126  // calculate delta:
127  delta=newu1-newu2;
128  //delta.Print();
129 
130  // get biggest contribution:
131  int index=0;
132  d=0;
133  for(int i=0; i<n; ++i){
134  if(TMath::Abs(delta[i])>d){
135  d=TMath::Abs(delta[i]);
136  index=i;
137  }
138  }
139  d0=TMath::Abs(_acc*h*uprime[index]);
140 
141  //std::cout<<"d="<<d<<" d0="<<d0
142  // <<" zdif="<<newu1[2]-newu2[2]
143  // <<" newh="<<newh<<std::endl;
144 
145  assert(d>=0);
146  if(d>0){
147  double frac=TMath::Abs(d0/d);
148  newh=0.9*TMath::Abs(newh*(TMath::Power(frac,0.25)));
149  if(newh<0.00001){
150  newh=0.00001;
151  break;
152  }
153  if(newh>_maxstep)newh=_maxstep;
154  }
155  else {
156  //std::cout<<"small dev! newh="<<newh<<std::endl;
157  newh=2*newh;
158  if(newh>_maxstep)newh=_maxstep;
159  break;
160  }
161  //break;
162  }
163 
164  //std::cout<<"completed step d="<<d
165  // <<" d0="<<d0
166  // <<" newh="<<newh
167  // <<" stepdone="<<stepdone<<std::endl;
168  newu=newu1;//+0.06666666667*delta;
169  newuprime=newuprime1;
170  assert(newh>0);
171  return newh;
172 }
173 
174 
175 // returns length
176 double
177 Nystrom::propagate(double start, double end,
178  const TVectorT<double>& u,
179  const TVectorT<double>& uprime,
180  const TVectorT<double>& par,
181  TVectorT<double>& newu, TVectorT<double>& newuprime)
182 {
183  if(TMath::Abs(start-end)<1E-8){
184  newu=u;
185  newuprime=uprime;
186  return 0;
187  }
188 
189  double length=0;
190 
191  double sign= start < end ? 1. : -1.;
192  double h=_acc; // TODO: This should be more dynamic!
193  double s=start;
194  TVectorT<double> utmp(u);
195  TVectorT<double> uprimetmp(uprime);
196  unsigned int steps=0;
197  unsigned int maxsteps=1000000;
198  double z0=u[2];
199  double z=z0;
200  assert(fabs(z0-start)<1E-8);
201 
202  while(TMath::Abs(s-end)>h && steps<maxsteps){
203 
204  ++steps;
205 
206  double reststep=TMath::Abs(s-end);
207  double newh=h;
208  double stepdone=sign*h;
209  if(_adaptive){
210  newh=adaptiveStep(utmp,uprimetmp,par,
211  newu,newuprime,stepdone,h,sign,reststep);
212  }
213  else {
214  step(utmp,uprimetmp,par,
215  newu,newuprime,sign*h);
216  }
217  //if(fabs(newh-h)>1E-4)std::cout<<"h="<<h<<" newh="<<newh<<std::endl;
218  double dx=h*uprimetmp[0];
219  double dy=h*uprimetmp[1];
220 
221  length+=sqrt(dx*dx+dy*dy+h*h);
222 
223 
224  s+=stepdone;
225  z+=stepdone;
226 
227  h=newh;
228  utmp=newu;
229  uprimetmp=newuprime;
230 
231  // todo: remove hard limits
232  if(fabs(utmp[0])>1000 || fabs(utmp[1])>1000){
233  std::cerr<<"Nystrom:: outside of range (-1000,1000)!"<<std::endl;
234  return length;
235  }
236 
237  }
238 
239 
240 
241 
242  if(steps>=maxsteps)std::cerr<<"Nystrom:: too many steps (>"<<maxsteps<<") aborting."<<std::endl;
243  // do the rest step
244  h=end-s;
245  if(h>1E-8){
246  step(utmp,uprimetmp,par,
247  newu,newuprime,h);
248  }
249  double dx=h*uprimetmp[0];
250  double dy=h*uprimetmp[1];
251 
252  length+=sqrt(dx*dx+dy*dy+h*h);
253  return length;
254  //double avstep=fabs(z0-z)/(double)steps;
255  //std::cout<<"z0="<<z0
256  // <<" z="<<z
257  // <<" newu[2]="<<newu[2]
258  // <<" avstep="<<avstep
259  // <<" signe="<<sign<<std::endl;
260 
261 }