When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / ITS / AliITSNeuralTrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // AliITSNeuralTrack
19 //
20 // The format of output data from Neural Tracker
21 // It can export data in the format of AliITSIOTrack
22 // (V1) tracking.
23 // Compatibility adaptation to V2 tracking is on the way.
24 // Author: A. Pulvirenti
25
26 #include <Riostream.h>
27 //#include <cstdlib>
28 //#include <cstring>
29
30 //#include <TObject.h>
31 //#include <TROOT.h>
32 #include <TMath.h>
33 #include <TString.h>
34 //#include <TObjArray.h>
35 //#include <TH1.h>
36 #include <TMatrixD.h>
37 #if ROOT_VERSION_CODE >= 262146
38 #include <TMatrixDEigen.h>
39 #endif
40
41 //#include "AliITSVertex.h"
42 #include "AliITSIOTrack.h"
43 #include "AliITSNeuralPoint.h"
44
45 #include "AliITSNeuralTrack.h"
46
47
48
49 ClassImp(AliITSNeuralTrack)
50 //
51 //
52 //
53 AliITSNeuralTrack::AliITSNeuralTrack() :  
54 fXC(0.0),
55 fYC(0.0),
56 fR(0.0),
57 fC(0.0),
58 fTanL(0.0),
59 fG0(0.0),
60 fDt(0.0),
61 fDz(0.0),
62 fStateR(0.0),
63 fStatePhi(0.0),
64 fStateZ(0.0),
65 fMatrix(5,5),
66 fChi2(0.0),
67 fNSteps(0.0),
68 fMass(0.1396),// default assumption: pion
69 fField(2.0),// default assumption: B = 0.4 Tesla
70 fPDG(0),
71 fLabel(0),
72 fCount(0),
73 fVertex(){
74         // Default constructor
75         
76         Int_t i;
77
78         for (i = 0; i < 6; i++) fPoint[i] = 0;
79
80         fVertex.X() = 0.0;
81         fVertex.Y() = 0.0;
82         fVertex.Z() = 0.0;
83         fVertex.ErrX() = 0.0;
84         fVertex.ErrY() = 0.0;
85         fVertex.ErrZ() = 0.0;
86 }
87 //
88 //
89 //
90 AliITSNeuralTrack::AliITSNeuralTrack(const AliITSNeuralTrack &track) 
91 : TObject((TObject&)track), 
92 fXC(0.0),
93 fYC(0.0),
94 fR(0.0),
95 fC(0.0),
96 fTanL(0.0),
97 fG0(0.0),
98 fDt(0.0),
99 fDz(0.0),
100 fStateR(0.0),
101 fStatePhi(0.0),
102 fStateZ(0.0),
103 fMatrix(5,5),
104 fChi2(0.0),
105 fNSteps(0.0),
106 fMass(0.1396),// default assumption: pion
107 fField(2.0),// default assumption: B = 0.4 Tesla
108 fPDG(0),
109 fLabel(0),
110 fCount(0),
111 fVertex(){
112 // copy constructor
113
114         Int_t i;
115
116         fMass  = 0.1396;   // default assumption: pion
117         fField = 2.0;      // default assumption: B = 0.4 Tesla
118
119         fXC = fYC = fR = fC = 0.0;
120         fTanL = fG0 = fDt = fDz = 0.0;
121         fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
122
123         fLabel = 0;
124         fCount = 0;
125         for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i];
126
127         fVertex.X() = 0.0;
128         fVertex.Y() = 0.0;
129         fVertex.Z() = 0.0;
130         fVertex.ErrX() = 0.0;
131         fVertex.ErrY() = 0.0;
132         fVertex.ErrZ() = 0.0;   
133 }
134 //
135 //
136 //
137
138 AliITSNeuralTrack& AliITSNeuralTrack::operator=(const AliITSNeuralTrack& track){
139   //assignment operator
140   this->~AliITSNeuralTrack();
141   new(this) AliITSNeuralTrack(track);
142   return *this;
143 }
144 AliITSNeuralTrack::~AliITSNeuralTrack()
145 {
146         Int_t l;
147         for (l = 0; l < 6; l++) fPoint[l] = 0;
148 }
149 //
150 //
151 //
152 void AliITSNeuralTrack::AssignLabel()
153 {       
154 // Assigns a GEANT label to the found track. 
155 // Every cluster has up to three labels (it can have less). Then each label is 
156 // recorded for each point. Then, counts are made to check if some of the labels
157 // appear more than once. Finally, the label which appears most times is assigned
158 // to the track in the field fLabel.
159 // The number of points containing that label is counted in the fCount data-member. 
160
161         Bool_t found;
162         Int_t i, j, l, lab, max = 0;
163         
164         // We have up to 6 points for 3 labels each => up to 18 possible different values
165         Int_t idx[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
166         Int_t count[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
167         
168         for (l = 0; l < 6; l++) {
169                 if (!fPoint[l]) continue;
170                 // Sometimes the same label appears two times in the same recpoint.
171                 // With these if statements, such problem is solved by turning
172                 // one of them to -1.
173                 if (fPoint[l]->GetLabel(1) >= 0 && fPoint[l]->GetLabel(1) == fPoint[l]->GetLabel(0))
174                         fPoint[l]->SetLabel(1, -1);
175                 if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(0))
176                         fPoint[l]->SetLabel(2, -1);
177                 if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(1))
178                         fPoint[l]->SetLabel(2, -1);
179                 for (i = 0; i < 3; i++) {
180                         lab = fPoint[l]->GetLabel(i);
181                         if (lab < 0) continue;
182                         found = kFALSE;
183                         for (j = 0; j < max; j++) {
184                                 if (idx[j] == lab) {
185                                         count[j]++;
186                                         found = kTRUE;
187                                 }
188                         }
189                         if(!found) {
190                                 max++;
191                                 idx[max - 1] = lab;
192                                 count[max - 1] = 1;
193                         }
194                 }
195         }
196         
197         j = 0, max = count[0];
198         for (i = 0; i < 18; i++) {
199                 if (count[i] > max) {
200                         j = i;
201                         max = count[i];
202                 }
203         }
204         fLabel = idx[j];
205         fCount = count[j];
206 }
207 //
208 //
209 //
210 void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
211 {
212 // Removes a point from the corresponding layer slot in the found track.
213 // If the argument is TRUE, the point object is also deleted from heap.
214
215         if (i >= 0 && i < 6) {
216                 if (del) delete fPoint[i];
217                 fPoint[i] = 0;
218         }
219 }
220 //
221 //
222 //
223 void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
224 {
225 // Returns the point coordinates according to the TreeR philosophy in galice.root files
226 // that consist in the module number (mod) and the position in the TClonesArray of
227 // the points reconstructed in that module for the run being examined.
228
229         if (layer < 0 || layer > 5) {
230                 Error("GetModuleData", "Layer out of range: %d", layer);
231                 return;
232         }
233         mod = fPoint[layer]->GetModule();
234         pos = fPoint[layer]->GetIndex();
235 }
236 //
237 //
238 //
239 void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
240 {
241 // A trivial method to insert a point in the tracks;
242 // the point is inserted to the slot corresponding to its ITS layer.
243
244         if (!point) return;
245         
246         Int_t layer = point->GetLayer();
247         if (layer < 0 || layer > 6) {
248                 Error("Insert", "Layer index %d out of range", layer);
249                 return;
250         } 
251         
252         fPoint[layer] = point;
253 }
254 //
255 //
256 //
257 Int_t AliITSNeuralTrack::OccupationMask() const
258 {
259 // Returns a byte which maps the occupied slots. 
260 // Each bit represents a layer going from the less significant on.
261
262         Int_t i, check, mask = 0;
263         for (i = 0; i < 6; i++) {
264                 check = 1 << i;
265                 if (fPoint[i]) mask |= check;
266         }
267         return mask;
268 }
269 //
270 //
271 //
272 void AliITSNeuralTrack::PrintLabels()
273 {
274 // Prints the results of the AssignLabel() method, together with
275 // the GEANT labels assigned to each point, in order to evaluate 
276 // how the assigned label is distributed among points.
277
278         cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
279         for (Int_t i = 0; i < 6; i++) {
280                 cout << "Point #" << i + 1 << " --> ";
281                 if (fPoint[i]) {
282                         cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
283                         cout << fPoint[i]->GetLabel(1) << ", ";
284                         cout << fPoint[i]->GetLabel(2) << endl;
285                 }
286                 else {
287                         cout << "not assigned" << endl;
288                 }
289         }
290         cout << endl;
291 }
292 //
293 //
294 //
295 Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
296 {
297 // Calculates the correction for energy loss
298
299         Double_t width = 0.0;
300         switch (layer) {
301                 case 0: width = 0.00260 + 0.00283; break;
302                 case 1: width = 0.0180; break;
303                 case 2: width = 0.0094; break;
304                 case 3: width = 0.0095; break;
305                 case 4: width = 0.0091; break;
306                 case 5: width = 0.0087; break;
307                 default:
308                         Error("AddEL", "Layer value %d out of range!", layer);
309                         return kFALSE;
310         }
311         width *= 1.7;
312
313         if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
314                 width += 0.012;
315         }
316
317         Double_t invSqCosL = 1. + fTanL * fTanL;            // = 1 / (cos(lambda)^2) = 1 + tan(lambda)^2
318         Double_t invCosL   = TMath::Sqrt(invSqCosL);        // = 1  / cos(lambda)
319         Double_t pt = GetPt();                              // = transverse momentum
320         Double_t p2 = pt *pt * invSqCosL;                   // = square modulus of momentum
321         Double_t energy = TMath::Sqrt(p2 + fMass * fMass);  // = energy
322         Double_t beta2 = p2 / (p2 + fMass * fMass);         // = (v / c) ^ 2
323         if (beta2 == 0.0) {
324                 printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
325                 return kFALSE;
326         }
327
328         Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
329         dE = sign * dE * 0.001;
330
331         energy += dE;
332         p2 = energy * energy - fMass * fMass;
333         pt = TMath::Sqrt(p2) / invCosL;
334         if (fC < 0.) pt = -pt;
335         fC = (0.299792458 * 0.2 * fField) / (pt * 100.);
336
337         return kTRUE;
338 }
339 //
340 //
341 //
342 Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
343 {
344 // Calculates the noise perturbation due to multiple scattering
345
346         Double_t width = 0.0;
347         switch (layer) {
348                 case 0: width = 0.00260 + 0.00283; break;
349                 case 1: width = 0.0180; break;
350                 case 2: width = 0.0094; break;
351                 case 3: width = 0.0095; break;
352                 case 4: width = 0.0091; break;
353                 case 5: width = 0.0087; break;
354                 default:
355                         Error("AddEL", "Layer value %d out of range!", layer);
356                         return kFALSE;
357         }
358         width *= 1.7;
359
360         if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
361                 width += 0.012;
362         }
363
364         Double_t cosL = TMath::Cos(TMath::ATan(fTanL));
365         Double_t halfC = fC / 2.;
366         Double_t q20 = 1. / (cosL * cosL);
367         Double_t q30 = fC * fTanL;
368
369         Double_t q40 = halfC * (fStateR * fStateR - fDt * fDt) / (1. + 2. * halfC * fDt);
370         Double_t dd  = fDt + halfC * fDt * fDt - halfC * fStateR * fStateR;
371         Double_t dprova = fStateR * fStateR - dd * dd;
372         Double_t q41 = 0.;
373         if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
374
375         Double_t p2 = (GetPt()*GetPt()) / (cosL * cosL);
376         Double_t beta2 = p2 / (p2 + fMass * fMass);
377         Double_t theta2 = 14.1 * 14.1 / (beta2 * p2 * 1.e6) * (width / TMath::Abs(cosL));
378
379         fMatrix(2,2) += theta2 * (q40 * q40 + q41 * q41);
380         fMatrix(3,2) += theta2 * q20 * q40;
381         fMatrix(2,3) += theta2 * q20 * q40;
382         fMatrix(3,3) += theta2 * q20 * q20;
383         fMatrix(4,2) += theta2 * q30 * q40;
384         fMatrix(2,4) += theta2 * q30 * q40;
385         fMatrix(4,3) += theta2 * q30 * q20;
386         fMatrix(3,4) += theta2 * q30 * q20;
387         fMatrix(4,4) += theta2 * q30 * q30;
388         
389         return kTRUE;
390 }
391 //
392 //
393 //
394 Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
395 {
396         // Propagation method.
397         // Changes the state vector according to a new radial position
398         // which is specified by the passed 'r' value (in cylindircal coordinates).
399         // The covariance matrix is also propagated (and enlarged) according to
400         // the FCFt technique, where F is the jacobian of the new parameters
401         // w.r.t. their old values.
402         // The option argument forces the method to add also the energy loss
403         // and the multiple scattering effects, which respectively have the effect
404         // of changing the curvature and widening the covariance matrix.
405
406         if (rk < TMath::Abs(fDt)) {
407                 Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
408                 return 0;
409         }
410         
411         Double_t duepi = 2. * TMath::Pi();
412         Double_t rkm1 = fStateR;
413         Double_t aAk = ArgPhi(rk), aAkm1 = ArgPhi(rkm1);
414         Double_t ak = ArgZ(rk), akm1 = ArgZ(rkm1);
415
416         fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
417         if(fStatePhi > duepi) fStatePhi -= duepi;
418         if(fStatePhi < 0.) fStatePhi += duepi;
419
420         Double_t halfC = 0.5 * fC;
421         fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
422         
423         Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
424         Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
425         
426         Double_t f02 = ck / TMath::Sqrt((1.-aAk)*(1.+aAk)) - ckm1 / TMath::Sqrt((1.-aAkm1)*(1.+aAkm1));
427         Double_t f04 = bk / TMath::Sqrt((1.-aAk)*(1.+aAk)) - bkm1 / TMath::Sqrt((1.-aAkm1)*(1.+aAkm1));
428         Double_t f12 = fTanL * fDt * (1. / rk - 1. / rkm1);
429         Double_t f13 = rk - rkm1;
430         
431         Double_t c00 = fMatrix(0,0);
432         Double_t c10 = fMatrix(1,0);
433         Double_t c11 = fMatrix(1,1);
434         Double_t c20 = fMatrix(2,0);
435         Double_t c21 = fMatrix(2,1);
436         Double_t c22 = fMatrix(2,2);
437         Double_t c30 = fMatrix(3,0);
438         Double_t c31 = fMatrix(3,1);
439         Double_t c32 = fMatrix(3,2);
440         Double_t c33 = fMatrix(3,3);
441         Double_t c40 = fMatrix(4,0);
442         Double_t c41 = fMatrix(4,1);
443         Double_t c42 = fMatrix(4,2);
444         Double_t c43 = fMatrix(4,3);
445         Double_t c44 = fMatrix(4,4);
446
447         Double_t r10 = c10 + c21*f02 + c41*f04;
448         Double_t r20 = c20 + c22*f02 + c42*f04;
449         Double_t r30 = c30 + c32*f02 + c43*f04;
450         Double_t r40 = c40 + c42*f02 + c44*f04;
451         Double_t r21 = c21 + c22*f12 + c32*f13;
452         Double_t r31 = c31 + c32*f12 + c33*f13;
453         Double_t r41 = c41 + c42*f12 + c43*f13;
454
455         fMatrix(0,0) = c00 + c20*f02 + c40*f04 + f02*r20 + f04*r40;
456         fMatrix(1,0) = fMatrix(0,1) = r10 + f12*r20 + f13*r30;
457         fMatrix(1,1) = c11 + c21*f12 + c31*f13 + f12*r21 + f13*r31;
458         fMatrix(2,0) = fMatrix(0,2) = r20;
459         fMatrix(2,1) = fMatrix(1,2) = r21;
460         fMatrix(3,0) = fMatrix(0,3) = r30;
461         fMatrix(3,1) = fMatrix(1,3) = r31;
462         fMatrix(4,0) = fMatrix(0,4) = r40;
463         fMatrix(4,1) = fMatrix(1,4) = r41;
464
465         fStateR = rk;
466
467         if (rkm1 < fStateR)  // going to greater R --> energy LOSS
468                 return -1;
469         else                 // going to smaller R --> energy GAIN
470                 return 1;
471 }
472 //
473 //
474 //
475 Bool_t AliITSNeuralTrack::SeedCovariance()
476 {
477         // generate a covariance matrix depending on the results obtained from
478         // the preliminary seeding fit procedure.
479         // It calculates the variances for C, D ans TanL, according to the
480         // differences of the fitted values from the requested ones necessary
481         // to make the curve exactly pass through each point.
482
483         /*
484         Int_t i, j;
485         AliITSNeuralPoint *p = 0;
486         Double_t r, argPhi, phiC, phiD, argZ, zL;
487         Double_t sumC = 0.0, sumD = 0.0, sumphi = 0., sumz = 0., sumL = 0.;
488         for (i = 0; i < fNum; i++) {
489                 p = At(i);
490                 if (!p) continue;
491                 r = p->GetR2();
492                 // weight and derivatives of phi and zeta w.r.t. various params
493                 sumphi += 1./ p->ErrorGetPhi();
494                 argPhi = ArgPhi(r);
495                 argZ = ArgZ(r);
496                 if (argPhi > 100.0 || argZ > 100.0) {
497                         Error("InitCovariance", "Argument error");
498                         return kFALSE;
499                 }
500                 phiC = DerArgPhiC(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
501                 phiD = DerArgPhiD(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
502                 if (phiC > 100.0 || phiD > 100.0) {
503                         Error("InitCovariance", "Argument error");
504                         return kFALSE;
505                 }
506                 zL = asin(argZ) / fC;
507                 sumL += zL * zL;
508                 sumC += phiC * phiC;
509                 sumD += phiD * phiD;
510                 sumz += 1.0 / (p->fError[2] * p->fError[2]);
511         }
512
513         for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) fMatrix(i,j) = 0.;
514         fMatrix(0,0) = 1. / sumphi;
515         fMatrix(1,1) = 1. / sumz;
516         fMatrix(2,2) = 1. / sumD;
517         fMatrix(3,3) = 1. / sumL;
518         fMatrix(4,4) = 1. / sumC;
519         fMatrix.Print();
520         */
521         
522         AliITSNeuralPoint *p = 0;
523         Double_t delta, cs, sn, r, argz;
524         Double_t diffC, diffD, diffL, calcC, calcD, calcL;
525
526         Int_t l;
527         for (l = 0; l < 6; l++) {
528                 p = fPoint[l];
529                 if (!p) break;
530                 sn = TMath::Sin(p->GetPhi() - fG0);
531                 cs = TMath::Cos(p->GetPhi() - fG0);
532                 r  = p->GetR2();
533                 calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
534                 argz = ArgZ(r);
535                 if (argz > 1000.0) {
536                         Error("Covariance", "Value too high");
537                         return kFALSE;
538                 }
539                 calcL = (p->Z() - fDz) * fC / asin(argz);
540                 delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
541                 if (delta < 0.E0) {
542                         if (delta >= -0.5)
543                                 delta = 0.;
544                         else {
545                                 Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
546                                 return kFALSE;
547                         }
548                 }
549                 delta = sqrt(delta);
550                 if (fC >= 0)
551                         calcD = delta - fR;
552                 else
553                         calcD = fR - delta;
554                 diffD = calcD - fDt;
555                 diffL = calcL - fTanL;
556                 diffC = fC - calcC;
557                 fMatrix(0,0) += 100000000.0 * p->GetError("phi") * p->GetError("phi");
558                 fMatrix(1,1) += 10000.0 * p->ErrZ() * p->ErrZ();
559                 fMatrix(2,2) += 100000.0 * diffD * diffD;
560                 fMatrix(3,3) += diffL * diffL;
561                 fMatrix(4,4) += 100000000.0 * diffC * diffC;
562         }
563         Double_t n = 0.;
564         for (l = 0; l < 6; l++) if (fPoint[l]) n++;
565         fMatrix *= 1./(n * (n+1));
566         return kTRUE;
567 }
568 //
569 //
570 //
571 Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
572 {
573         // Makes all calculations which apply the Kalman filter to the
574         // stored guess of the state vector, after propagation to a new layer
575
576         if (!test) {
577                 Error("Filter", "Null pointer passed");
578                 return kFALSE;
579         }
580         
581         Double_t m[2];
582         Double_t rk, phik, zk;
583         rk = test->GetR2();
584         phik = test->GetPhi();
585         zk = test->Z();
586         m[0]=phik;
587         m[1]=zk;
588         
589         //////////////////////// Evaluation of the error matrix V  /////////
590         Double_t v00 = test->GetError("phi") * rk;
591         Double_t v11 = test->ErrZ();
592         ////////////////////////////////////////////////////////////////////  
593         
594         // Get the covariance matrix
595         Double_t cin00, cin10, cin20, cin30, cin40;
596         Double_t cin11, cin21, cin31, cin41, cin22;
597         Double_t cin32, cin42, cin33, cin43, cin44;
598         cin00 = fMatrix(0,0);
599         cin10 = fMatrix(1,0);
600         cin20 = fMatrix(2,0);
601         cin30 = fMatrix(3,0);
602         cin40 = fMatrix(4,0);
603         cin11 = fMatrix(1,1);
604         cin21 = fMatrix(2,1);
605         cin31 = fMatrix(3,1);
606         cin41 = fMatrix(4,1);
607         cin22 = fMatrix(2,2);
608         cin32 = fMatrix(3,2);
609         cin42 = fMatrix(4,2);
610         cin33 = fMatrix(3,3);
611         cin43 = fMatrix(4,3);
612         cin44 = fMatrix(4,4);
613         
614         // Calculate R matrix
615         Double_t rold00 = cin00 + v00;
616         Double_t rold10 = cin10;
617         Double_t rold11 = cin11 + v11;
618         
619         ////////////////////// R matrix inversion  /////////////////////////
620         Double_t det = rold00*rold11 - rold10*rold10;
621         Double_t r00 = rold11/det;
622         Double_t r10 = -rold10/det;
623         Double_t r11 = rold00/det;
624         ////////////////////////////////////////////////////////////////////
625         
626         // Calculate Kalman matrix
627         Double_t k00 = cin00*r00 + cin10*r10;
628         Double_t k01 = cin00*r10 + cin10*r11;
629         Double_t k10 = cin10*r00 + cin11*r10;  
630         Double_t k11 = cin10*r10 + cin11*r11;
631         Double_t k20 = cin20*r00 + cin21*r10;  
632         Double_t k21 = cin20*r10 + cin21*r11;  
633         Double_t k30 = cin30*r00 + cin31*r10;  
634         Double_t k31 = cin30*r10 + cin31*r11;  
635         Double_t k40 = cin40*r00 + cin41*r10;
636         Double_t k41 = cin40*r10 + cin41*r11;
637         
638         // Get state vector (will keep the old values for phi and z)
639         Double_t x0, x1, x2, x3, x4, savex0, savex1;
640         x0 = savex0 = fStatePhi;
641         x1 = savex1 = fStateZ;
642         x2 = fDt;
643         x3 = fTanL;
644         x4 = fC;
645
646         // Update the state vector
647         x0 += k00*(m[0]-savex0) + k01*(m[1]-savex1);
648         x1 += k10*(m[0]-savex0) + k11*(m[1]-savex1);
649         x2 += k20*(m[0]-savex0) + k21*(m[1]-savex1);
650         x3 += k30*(m[0]-savex0) + k31*(m[1]-savex1);
651         x4 += k40*(m[0]-savex0) + k41*(m[1]-savex1);
652         
653         // Update the covariance matrix
654         Double_t cout00, cout10, cout20, cout30, cout40;
655         Double_t cout11, cout21, cout31, cout41, cout22;
656         Double_t cout32, cout42, cout33, cout43, cout44;
657         
658         cout00 = cin00 - k00*cin00 - k01*cin10;
659         cout10 = cin10 - k00*cin10 - k01*cin11;
660         cout20 = cin20 - k00*cin20 - k01*cin21;
661         cout30 = cin30 - k00*cin30 - k01*cin31;
662         cout40 = cin40 - k00*cin40 - k01*cin41;
663         cout11 = cin11 - k10*cin10 - k11*cin11;
664         cout21 = cin21 - k10*cin20 - k11*cin21;
665         cout31 = cin31 - k10*cin30 - k11*cin31;
666         cout41 = cin41 - k10*cin40 - k11*cin41;
667         cout22 = cin22 - k20*cin20 - k21*cin21;
668         cout32 = cin32 - k20*cin30 - k21*cin31;
669         cout42 = cin42 - k20*cin40 - k21*cin41;
670         cout33 = cin33 - k30*cin30 - k31*cin31;
671         cout43 = cin43 - k30*cin40 - k31*cin41;
672         cout44 = cin44 - k40*cin40 - k41*cin41;
673         
674         // Store the new covariance matrix
675         fMatrix(0,0) = cout00;
676         fMatrix(1,0) = fMatrix(0,1) = cout10;
677         fMatrix(2,0) = fMatrix(0,2) = cout20;
678         fMatrix(3,0) = fMatrix(0,3) = cout30;
679         fMatrix(4,0) = fMatrix(0,4) = cout40;
680         fMatrix(1,1) = cout11;
681         fMatrix(2,1) = fMatrix(1,2) = cout21;
682         fMatrix(3,1) = fMatrix(1,3) = cout31;
683         fMatrix(4,1) = fMatrix(1,4) = cout41;
684         fMatrix(2,2) = cout22;
685         fMatrix(3,2) = fMatrix(2,3) = cout32;
686         fMatrix(4,2) = fMatrix(2,4) = cout42;
687         fMatrix(3,3) = cout33;
688         fMatrix(4,3) = fMatrix(3,4) = cout43;
689         fMatrix(4,4) = cout44;
690         
691         // Calculation of the chi2 increment
692         Double_t vmcold00 = v00 - cout00;
693         Double_t vmcold10 = -cout10;
694         Double_t vmcold11 = v11 - cout11;
695         ////////////////////// Matrix vmc inversion  ///////////////////////
696         det = vmcold00*vmcold11 - vmcold10*vmcold10;
697         Double_t vmc00=vmcold11/det;
698         Double_t vmc10 = -vmcold10/det;
699         Double_t vmc11 = vmcold00/det;
700         ////////////////////////////////////////////////////////////////////
701         Double_t chi2 = (m[0] - x0)*( vmc00*(m[0] - x0) + 2.*vmc10*(m[1] - x1) ) + (m[1] - x1)*vmc11*(m[1] - x1);
702         fChi2 += chi2;
703         fNSteps++;
704         
705         return kTRUE;
706 }
707 //
708 //
709 //
710 Bool_t AliITSNeuralTrack::KalmanFit()
711 {
712 // Applies the Kalman Filter to improve the track parameters resolution.
713 // First, thre point which lies closer to the estimated helix is chosen.
714 // Then, a fit is performed towards the 6th layer
715 // Finally, the track is refitted to the 1st layer
716
717         Double_t rho;
718         Int_t l, layer, sign;
719         
720         fStateR = fPoint[0]->GetR2();
721         fStatePhi = fPoint[0]->GetPhi();
722         fStateZ = fPoint[0]->Z();
723         
724         if (!PropagateTo(3.0)) {
725                 Error("KalmanFit", "Unsuccessful initialization");
726                 return kFALSE;
727         }
728         l=0;
729
730         // Performs a Kalman filter going from the actual state position
731         // towards layer 6 position
732         // Now, the propagation + filtering operations can be performed
733         Double_t argPhi = 0.0, argZ = 0.0;
734         while (l <= 5) {
735                 if (!fPoint[l]) {
736                         Error("KalmanFit", "Not six points!");
737                         return kFALSE;
738                 }
739                 layer = fPoint[l]->GetLayer();
740                 rho = fPoint[l]->GetR2();
741                 sign = PropagateTo(rho);
742                 if (!sign) return kFALSE;
743                 AddEL(layer, -1.0);
744                 AddMS(layer);
745                 if (!Filter(fPoint[l])) return kFALSE;
746                 // these two parameters are update according to the filtered values
747                 argPhi = ArgPhi(fStateR);
748                 argZ = ArgZ(fStateR);
749                 if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
750                         Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
751                         return kFALSE;
752                 }
753                 fG0 = fStatePhi - asin(argPhi);
754                 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
755                 l++;
756         }
757
758         // Now a Kalman filter i performed going from the actual state position
759         // towards layer 1 position and then propagates to vertex
760         if (l >= 5) l = 5;
761         while (l >= 1) {
762                 layer = fPoint[l]->GetLayer();
763                 rho = fPoint[l]->GetR2();
764                 AddEL(layer, 1.0);
765                 sign = PropagateTo(rho);
766                 if (!sign) return kFALSE;
767                 AddMS(layer);
768                 if (!Filter(fPoint[l])) return kFALSE;
769                 // these two parameters are update according to the filtered values
770                 argPhi = ArgPhi(fStateR);
771                 argZ = ArgZ(fStateR);
772                 if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
773                         Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
774                         return kFALSE;
775                 }
776                 fG0 = fStatePhi - asin(argPhi);
777                 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
778                 l--;
779         }
780         return kTRUE;
781 }
782 //
783 //
784 //
785 Bool_t AliITSNeuralTrack::RiemannFit()
786 {
787         // Method which executes the circle fit via a Riemann Sphere projection
788         // with the only improvement of a weighted mean, due to different errors
789         // over different point measurements.
790         // As an output, it returns kTRUE or kFALSE respectively if the fit succeeded or not
791         // in fact, if some variables assume strange values, the fit is aborted,
792         // in order to prevent the class from raising a floating point error;
793
794         Int_t i, j;
795
796         // M1 - matrix of ones
797         TMatrixD m1(6,1);
798         for (i = 0; i < 6; i++) m1(i,0) = 1.0;
799
800         // X - matrix of Rieman projection coordinates
801         TMatrixD mX(6,3);
802         for (i = 0; i < 6; i++) {
803                 mX(i,0) = fPoint[i]->X();
804                 mX(i,1) = fPoint[i]->Y();
805                 mX(i,2) = fPoint[i]->GetR2sq();
806         }
807
808         // W - matrix of weights
809         Double_t xterm, yterm, ex, ey;
810         TMatrixD mW(6,6);
811         for (i = 0; i < 6; i++) {
812                 xterm = fPoint[i]->X() * fPoint[i]->GetPhi() - fPoint[i]->Y() / fPoint[i]->GetR2();
813                 ex = fPoint[i]->ErrX();
814                 yterm = fPoint[i]->Y() * fPoint[i]->GetPhi() + fPoint[i]->X() / fPoint[i]->GetR2();
815                 ey = fPoint[i]->ErrY();
816                 mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
817         }
818
819         // Xm - weighted sample mean
820         Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
821         for (i = 0; i < 6; i++) {
822                 meanX += mW(i,i) * mX(i,0);
823                 meanY += mW(i,i) * mX(i,1);
824                 meanW += mW(i,i) * mX(i,2);
825                 sw += mW(i,i);
826         }
827         meanX /= sw;
828         meanY /= sw;
829         meanW /= sw;
830
831         // V - sample covariance matrix
832         for (i = 0; i < 6; i++) {
833                 mX(i,0) -= meanX;
834                 mX(i,1) -= meanY;
835                 mX(i,2) -= meanW;
836         }
837         TMatrixD mXt(TMatrixD::kTransposed, mX);
838         TMatrixD mWX(mW, TMatrixD::kMult, mX);
839         TMatrixD mV(mXt, TMatrixD::kMult, mWX);
840         for (i = 0; i < 3; i++) {
841                 for (j = i + 1; j < 3; j++) {
842                         mV(i,j)  = mV(j,i)  = (mV(i,j) + mV(j,i)) * 0.5;
843                 }
844         }
845
846         // Eigenvalue problem solving for V matrix
847         Int_t ileast = 0;
848         TVectorD eval(3), n(3);
849         //      TMatrixD evec = mV.EigenVectors(eval);
850 #if ROOT_VERSION_CODE >= 262146
851         TMatrixDEigen ei(mV);
852         TMatrixD evec = ei.GetEigenVectors();
853         eval = ei.GetEigenValues();
854 #else
855         TMatrixD evec = mV.EigenVectors(eval);
856 #endif
857
858         if (eval(1) < eval(ileast)) ileast = 1;
859         if (eval(2) < eval(ileast)) ileast = 2;
860         n(0) = evec(0, ileast);
861         n(1) = evec(1, ileast);
862         n(2) = evec(2, ileast);
863
864         // c - known term in the plane intersection with Riemann axes
865         Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
866
867         fXC = -n(0) / (2. * n(2));
868         fYC = -n(1) / (2. * n(2));
869         fR  = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
870
871         if (fR <= 0.E0) {
872                 Error("RiemannFit", "Radius comed less than zero!!!");
873                 return kFALSE;
874         }
875         fR = TMath::Sqrt(fR);
876         fC = 1.0 / fR;
877
878         // evaluating signs for curvature and others
879         Double_t phi1 = 0.0, phi2, temp1, temp2, sumdphi = 0.0, ref = TMath::Pi();
880         AliITSNeuralPoint *p = fPoint[0];
881         phi1 = p->GetPhi();
882         for (i = 1; i < 6; i++) {
883                 p = (AliITSNeuralPoint*)fPoint[i];
884                 if (!p) break;
885                 phi2 = p->GetPhi();
886                 temp1 = phi1;
887                 temp2 = phi2;
888                 if (temp1 > ref && temp2 < ref)
889                         temp2 += 2.0 * ref;
890                 else if (temp1 < ref && temp2 > ref)
891                         temp1 += 2.0 * ref;
892                 sumdphi += temp2 - temp1;
893                 phi1 = phi2;
894         }
895         if (sumdphi < 0.E0) fC = -fC;
896         Double_t diff, angle = TMath::ATan2(fYC, fXC);
897         if (fC < 0.E0)
898                 fG0 = angle + 0.5 * TMath::Pi();
899         else
900                 fG0 = angle - 0.5 * TMath::Pi();
901         diff = angle - fG0;
902
903         Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
904         if (fC >= 0.E0)
905                 fDt = d;
906         else
907                 fDt = -d;
908
909         Int_t nn = 6;
910         Double_t halfC = 0.5 * fC;
911         Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn];
912         for (j = 0; j < 6; j++) {
913                 p = fPoint[j];
914                 if (!p) break;
915                 s[j] = ArgZ(p->GetR2());
916                 if (s[j] > 100.0) return kFALSE;
917                 z[j] = p->Z();
918                 s[j] = asin(s[j]) / halfC;
919                 ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
920         }
921
922         // second tep final fit
923         Double_t sums2 = 0.0, sumz = 0.0, sumsz = 0.0, sums = 0.0, sumw = 0.0;
924         for (i = 0; i < nn; i++) {
925                 sums2 += ws[i] * s[i] * s[i];
926                 sumz  += ws[i] * z[i];
927                 sums  += ws[i] * s[i];
928                 sumsz += ws[i] * s[i] * z[i];
929                 sumw += ws[i];
930         }
931         sums2 /= sumw;
932         sumz /= sumw;
933         sums /= sumw;
934         sumsz /= sumw;
935         d = sums2 - sums*sums;
936
937         fDz = (sums2*sumz - sums*sumsz) / d;
938         fTanL = (sumsz - sums*sumz) / d;
939
940         delete [] s;
941         delete [] z;
942         delete [] ws;
943
944         return kTRUE;
945 }
946 //
947 //
948 //
949 void AliITSNeuralTrack::PrintState(Bool_t matrix)
950 {
951 // Prints the state vector values.
952 // The argument switches on or off the printing of the covariance matrix.
953
954         cout << "\nState vector: " << endl;
955         cout << " Rho = " << fStateR << "\n";
956         cout << " Phi = " << fStatePhi << "\n";
957         cout << "   Z = " << fStateZ << "\n";
958         cout << "  Dt = " << fDt << "\n";
959         cout << "  Dz = " << fDz << "\n";
960         cout << "TanL = " << fTanL << "\n";
961         cout << "   C = " << fC << "\n";
962         cout << "  G0 = " << fG0 << "\n";
963         cout << "  XC = " << fXC << "\n";
964         cout << "  YC = " << fYC << "\n";
965         if (matrix) {
966                 cout << "\nCovariance Matrix: " << endl;
967                 fMatrix.Print();
968         }
969         cout << "Actual square chi = " << fChi2;
970 }
971 //
972 //
973 //
974 Double_t AliITSNeuralTrack::GetDz() const
975 {
976 //      Double_t argZ = ArgZ(fStateR);
977 //      if (argZ > 9.9) {
978 //              Error("GetDz", "Too large value: %g", argZ);
979 //              return 0.0;
980 //      }
981 //      fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
982         return fDz;
983 }
984 //
985 //
986 //
987 Double_t AliITSNeuralTrack::GetGamma() const
988 {
989 // these two parameters are update according to the filtered values
990 //      Double_t argPhi = ArgPhi(fStateR);
991 //      if (argPhi > 9.9) {
992 //              Error("Filter", "Too large value: %g", argPhi);
993 //              return kFALSE;
994 //      }
995 //      fG0 = fStatePhi - asin(argPhi);
996         return fG0;
997 }
998 //
999 //
1000 //
1001 Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
1002 {
1003 // Gives the value of azymuthal coordinate in the helix
1004 // as a function of cylindric radius
1005
1006         Double_t arg = ArgPhi(r);
1007         if (arg > 0.9) return 0.0;
1008         arg = fG0 + asin(arg);
1009         while (arg >= 2.0 * TMath::Pi()) { arg -= 2.0 * TMath::Pi(); }
1010         while (arg < 0.0) { arg += 2.0 * TMath::Pi(); }
1011         return arg;
1012 }
1013 //
1014 //
1015 //
1016 Double_t AliITSNeuralTrack::GetZ(Double_t r) const
1017 {
1018 // gives the value of Z in the helix
1019 // as a function of cylindric radius
1020
1021         Double_t arg = ArgZ(r);
1022         if (arg > 0.9) return 0.0;
1023         return fDz + fTanL * asin(arg) / fC;
1024 }
1025 //
1026 //
1027 //
1028 Double_t AliITSNeuralTrack::GetdEdX()
1029 {
1030 // total energy loss of the track
1031
1032         Double_t q[4] = {0., 0., 0., 0.}, dedx = 0.0;
1033         Int_t i = 0, swap = 0;
1034         for (i = 2; i < 6; i++) {
1035                 if (!fPoint[i]) continue;
1036                 q[i - 2] = (Double_t)fPoint[i]->GetCharge();
1037                 q[i - 2] /= (1 + fTanL*fTanL);
1038         }
1039         q[0] /= 280.;
1040         q[1] /= 280.;
1041         q[2] /= 38.;
1042         q[3] /= 38.;
1043         do {
1044                 swap = 0;
1045                 for (i = 0; i < 3; i++) {
1046                         if (q[i] <= q[i + 1]) continue;
1047                         Double_t tmp = q[i];
1048                         q[i] = q[i + 1]; 
1049                         q[i+1] = tmp;
1050                         swap++;
1051                 }
1052         } while(swap); 
1053         if(q[0] < 0.) {
1054                 q[0] = q[1];
1055                 q[1] = q[2];
1056                 q[2] = q[3];
1057                 q[3] = -1.;
1058         } 
1059         dedx = (q[0] + q[1]) / 2.;
1060         return dedx;
1061 }
1062 //
1063 //
1064 //
1065 void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
1066 {
1067         // Stores vertex data
1068
1069         if (!pos || !err) return;
1070         fVertex.ErrX() = err[0];
1071         fVertex.ErrY() = err[1];
1072         fVertex.ErrZ() = err[2];
1073         fVertex.SetLayer(0);
1074         fVertex.SetModule(0);
1075         fVertex.SetIndex(0);
1076         fVertex.SetLabel(0, -1);
1077         fVertex.SetLabel(1, -1);
1078         fVertex.SetLabel(2, -1);
1079         fVertex.SetUser(1);
1080 }
1081 //
1082 //
1083 //
1084 AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
1085 {
1086 // Exports an object in the standard format for reconstructed tracks
1087
1088         Int_t layer = 0;
1089         AliITSIOTrack *track = new AliITSIOTrack;
1090
1091         // covariance matrix
1092         track->SetCovMatrix(fMatrix(0,0), fMatrix(1,0), fMatrix(1,1),
1093                             fMatrix(2,0), fMatrix(2,1), fMatrix(2,2),
1094                             fMatrix(3,0), fMatrix(3,1), fMatrix(3,2),
1095                             fMatrix(3,3), fMatrix(4,0), fMatrix(4,1),
1096                             fMatrix(4,2), fMatrix(4,3), fMatrix(4,4));
1097         
1098         // labels
1099         track->SetLabel(IsGood(min) ? fLabel : -fLabel);
1100         track->SetTPCLabel(-1);
1101
1102         // points characteristics
1103         for (layer = 0; layer < 6; layer++) {
1104                 if (fPoint[layer]) {
1105                         track->SetIdModule(layer, fPoint[layer]->GetModule());
1106                         track->SetIdPoint(layer, fPoint[layer]->GetIndex());
1107                 }
1108         }
1109         
1110         // state vector
1111         track->SetStatePhi(fStatePhi);
1112         track->SetStateZ(fStateZ);
1113         track->SetStateD(fDt);
1114         track->SetStateTgl(fTanL);
1115         track->SetStateC(fC);
1116         track->SetRadius(fStateR);
1117         track->SetCharge((fC > 0.0) ? -1 : 1);
1118         track->SetDz(fDz);
1119
1120         // track parameters in the closest point
1121         track->SetX(fStateR * cos(fStatePhi));
1122         track->SetY(fStateR * cos(fStatePhi));
1123         track->SetZ(fStateZ);
1124         track->SetPx(GetPt() * cos(fG0));
1125         track->SetPy(GetPt() * sin(fG0));
1126         track->SetPz(GetPt() * fTanL);
1127
1128         // PID
1129         track->SetPid(fPDG);
1130         track->SetMass(fMass);
1131
1132         return track;
1133 }
1134 //
1135 //
1136 //====================================================================================
1137 //============================       PRIVATE METHODS      ============================
1138 //====================================================================================
1139 //
1140 //
1141 Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
1142 {
1143         // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
1144
1145         Double_t arg, num, den;
1146         num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
1147         den = 1. + fC * fDt;
1148         if (den == 0.) {
1149                 Error("ArgPhi", "Denominator = 0!");
1150                 return 10.0;
1151         }
1152         arg = num / den;
1153         if (TMath::Abs(arg) < 1.) return arg;
1154         if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
1155         Error("ArgPhi", "Value too large: %17.15g", arg);
1156         return 10.0;
1157 }
1158 //
1159 //
1160 //
1161 Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
1162 {
1163         // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
1164
1165         Double_t arg;
1166         arg = (r * r - fDt * fDt) / (1. + fC * fDt);
1167         if (arg < 0.) {
1168                 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
1169                 else {
1170                         Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
1171                         return 10.;
1172                 }
1173         }
1174         arg = 0.5 * fC * TMath::Sqrt(arg);
1175         if (TMath::Abs(arg) < 1.) return arg;
1176         if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
1177         Error("ArgZ", "Value too large: %17.15g", arg);
1178         return 10.0;
1179 }
1180 //
1181 //
1182 //
1183 Double_t AliITSNeuralTrack::ArgB(Double_t r) const 
1184 {
1185 // UTILITY FUNCTION 
1186
1187         Double_t arg;
1188         arg = (r*r - fDt*fDt);
1189         arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
1190         return arg;
1191 }
1192 //
1193 //
1194 //
1195 Double_t AliITSNeuralTrack::ArgC(Double_t r) const 
1196 {
1197 // UTILITY FUNCTION
1198
1199         Double_t arg;
1200         arg = (1./r - fC * ArgPhi(r));
1201         arg /= 1.+ fC*fDt;
1202         return arg;
1203 }