Coverity fixes.
[u/mrichter/AliRoot.git] / STEER / AliRelAlignerKalman.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 ///////////////////////////////////////////////////////////////////////////////
17 //
18 //    Kalman filter based aligner:
19 //    Finds alignement constants for  two tracking volumes (by default ITS
20 //    and TPC)
21 //    Determines the inverse transformation of the second volume (TPC)
22 //    with respect to the first (ITS) (how to realign TPC to ITS)
23 //    by measuring the residual between the 2 tracks.
24 //    Additionally calculates some callibration parameters for TPC
25 //    Fit parameters are:
26 //    - 3 shifts, x,y,z
27 //    - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
28 //    - TPC drift velocity correction,
29 //    - TPC time offset correction.
30 //
31 //    Basic usage:
32 //    When aligning two volumes, at any given time a single instance of
33 //    the class should be active. The fit of the parameters is updated
34 //    by adding new data using one of the Add.... methods:
35 //
36 //    In collision events add an ESD event to update the fit (adds all tracks):
37 //
38 //        Bool_t AddESDevent( AliESDevent* pTrack );
39 //    
40 //    or add each individual track
41 //
42 //        AddESDtrack( AliESDtrack* pTrack );
43 //
44 //    For cosmic data, the assumption is that the tracking is done twice:
45 //    once global and once only ITS and the tracklets are saved inside
46 //    one AliESDEvent. The method
47 //
48 //        Bool_t AddCosmicEvent( AliESDEvent* pEvent );
49 //
50 //    then searches the event for matching tracklets and upon succes it updates.
51 //    One cosmic ideally triggers two updates: for the upper and lower half of
52 //    the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
53 //
54 //    by default give misalignment parameters for TPC as they appear to be.
55 //    TPC calibration parameters are always given as correction to values used in reco.
56 //
57 //    _________________________________________________________________________
58 //    Expert options:
59 //    look at AddESDevent() and AddCosmicEvent() to get the idea of how the
60 //    aligner works, it's safe to repeat the needed steps outside of the class,
61 //    only public methods are used.
62 //
63 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
64 //
65 //////////////////////////////////////////////////////////////////////////////
66
67 #include <iostream>
68 #include <TObject.h>
69 #include <TMath.h>
70 #include <TMatrix.h>
71 #include <TVector.h>
72 #include <TVector3.h>
73 #include <TDecompLU.h>
74 #include <TArrayI.h>
75 #include <TObjArray.h>
76 #include <TH1D.h>
77 #include <TF1.h>
78
79 #include "AliESDtrack.h"
80 #include "AliESDEvent.h"
81 #include "AliExternalTrackParam.h"
82
83 #include "AliRelAlignerKalman.h"
84
85 ClassImp(AliRelAlignerKalman)
86
87 //______________________________________________________________________________
88 AliRelAlignerKalman::AliRelAlignerKalman():
89     TObject(),
90     fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
91     fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
92     fMagField(0.),
93     fNMeasurementParams(4),
94     fPX(new TVectorD( fgkNSystemParams )),
95     fPXcov(new TMatrixDSym( fgkNSystemParams )),
96     fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
97     fQ(1.e-15),
98     fPMeasurement(new TVectorD( fNMeasurementParams )),
99     fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
100     fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
101     fOutRejSigmas(1.),
102     fOutRejSigma2Median(5.),
103     fYZOnly(kFALSE),
104     fNumericalParanoia(kTRUE),
105     fRejectOutliers(kTRUE),
106     fRejectOutliersSigma2Median(kFALSE),
107     fRequireMatchInTPC(kFALSE),
108     fCuts(kFALSE),
109     fMinPointsVol1(3),
110     fMinPointsVol2(50),
111     fMinPt(0.),
112     fMaxPt(1.e100),
113     fMaxMatchingAngle(0.1),
114     fMaxMatchingDistance(10.),  //in cm
115     fCorrectionMode(kFALSE),
116     fNTracks(0),
117     fNUpdates(0),
118     fNOutliers(0),
119     fNOutliersSigma2Median(0),
120     fNMatchedCosmics(0),
121     fNMatchedTPCtracklets(0),
122     fNProcessedEvents(0),
123     fTrackInBuffer(0),
124     fTimeStamp(0),
125     fRunNumber(0),
126     fNMerges(0),
127     fNMergesFailed(0),
128     fTPCvd(2.64),
129     fTPCZLengthA(2.4972500e02),
130     fTPCZLengthC(2.4969799e02)
131 {
132   //Default constructor
133   for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
134   for (Int_t i=0; i<4;i++){fResArrSigma2Median[i]=NULL;}
135   Reset();
136 }
137
138 //______________________________________________________________________________
139 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
140     TObject(static_cast<TObject>(a)),
141     fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
142     fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
143     fMagField(a.fMagField),
144     fNMeasurementParams(a.fNMeasurementParams),
145     fPX(new TVectorD( *a.fPX )),
146     fPXcov(new TMatrixDSym( *a.fPXcov )),
147     fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
148     fQ(a.fQ),
149     fPMeasurement(new TVectorD( fNMeasurementParams )),
150     fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
151     fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
152     fOutRejSigmas(a.fOutRejSigmas),
153     fOutRejSigma2Median(a.fOutRejSigma2Median),
154     fYZOnly(a.fYZOnly),
155     fNumericalParanoia(a.fNumericalParanoia),
156     fRejectOutliers(a.fRejectOutliers),
157     fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
158     fRequireMatchInTPC(a.fRequireMatchInTPC),
159     fCuts(a.fCuts),
160     fMinPointsVol1(a.fMinPointsVol1),
161     fMinPointsVol2(a.fMinPointsVol2),
162     fMinPt(a.fMinPt),
163     fMaxPt(a.fMaxPt),
164     fMaxMatchingAngle(a.fMaxMatchingAngle),
165     fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
166     fCorrectionMode(a.fCorrectionMode),
167     fNTracks(a.fNTracks),
168     fNUpdates(a.fNUpdates),
169     fNOutliers(a.fNOutliers),
170     fNOutliersSigma2Median(a.fNOutliersSigma2Median),
171     fNMatchedCosmics(a.fNMatchedCosmics),
172     fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
173     fNProcessedEvents(a.fNProcessedEvents),
174     fTrackInBuffer(0),
175     fTimeStamp(a.fTimeStamp),
176     fRunNumber(a.fRunNumber),
177     fNMerges(a.fNMerges),
178     fNMergesFailed(a.fNMergesFailed),
179     fTPCvd(a.fTPCvd),
180     fTPCZLengthA(a.fTPCZLengthA),
181     fTPCZLengthC(a.fTPCZLengthC)
182 {
183   //copy constructor
184   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
185
186   //copy contents of the residuals array for sigma2median scheme
187   for (Int_t i=0;i<4;i++)
188   {
189     if ((a.fResArrSigma2Median)[i]) 
190     {
191       fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
192       memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
193              fgkNtracksSigma2Median*sizeof(Double_t));
194     }
195     else
196       fResArrSigma2Median[i] = NULL;
197   }
198 }
199
200 //______________________________________________________________________________
201 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
202 {
203   //assignment operator
204   fMagField=a.fMagField,
205   fNMeasurementParams=a.fNMeasurementParams;
206   *fPX = *a.fPX;
207   *fPXcov = *a.fPXcov;
208   fQ=a.fQ;
209   fOutRejSigmas=a.fOutRejSigmas;
210   fOutRejSigma2Median=a.fOutRejSigma2Median;
211   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
212   fYZOnly=a.fYZOnly;
213   fNumericalParanoia=a.fNumericalParanoia;
214   fRejectOutliers=a.fRejectOutliers;
215   fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
216   fRequireMatchInTPC=a.fRequireMatchInTPC;
217   fCuts=a.fCuts;
218   fMinPointsVol1=a.fMinPointsVol1;
219   fMinPointsVol2=a.fMinPointsVol2;
220   fMinPt=a.fMinPt;
221   fMaxPt=a.fMaxPt;
222   fMaxMatchingAngle=a.fMaxMatchingAngle;
223   fMaxMatchingDistance=a.fMaxMatchingDistance;  //in c;
224   fCorrectionMode=a.fCorrectionMode;
225   fNTracks=a.fNTracks;
226   fNUpdates=a.fNUpdates;
227   fNOutliers=a.fNOutliers;
228   fNOutliersSigma2Median=a.fNOutliersSigma2Median;
229   fNMatchedCosmics=a.fNMatchedCosmics;
230   fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
231   fNProcessedEvents=a.fNProcessedEvents;
232   fTrackInBuffer=0; //because the array is reset, not copied
233   fTimeStamp=a.fTimeStamp;
234   fRunNumber=a.fRunNumber;
235   fNMerges=a.fNMerges;
236   fTPCvd=a.fTPCvd;
237   fTPCZLengthA=a.fTPCZLengthA;
238   fTPCZLengthC=a.fTPCZLengthC;
239
240   //copy contents of the residuals array for sigma2median scheme
241   for (Int_t i=0;i<4;i++)
242   {
243     if ((a.fResArrSigma2Median)[i]) 
244     {
245       if (!(fResArrSigma2Median[i])) fResArrSigma2Median[i] = 
246                                      new Double_t[fgkNtracksSigma2Median];
247       memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
248              fgkNtracksSigma2Median*sizeof(Double_t));
249     }
250     else
251       fResArrSigma2Median[i] = NULL;
252   }
253   return *this;
254 }
255
256 //______________________________________________________________________________
257 AliRelAlignerKalman::~AliRelAlignerKalman()
258 {
259   //destructor
260   delete [] fPTrackParamArr1;
261   delete [] fPTrackParamArr2;
262   delete fPX;
263   delete fPXcov;
264   delete fPH;
265   delete fPMeasurement;
266   delete fPMeasurementCov;
267   for (Int_t i=0;i<4;i++) 
268   {
269     delete [] (fResArrSigma2Median[i]);
270   }
271 }
272
273 //______________________________________________________________________________
274 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
275 {
276   //Add all tracks in an ESD event
277
278   fNProcessedEvents++; //update the counter
279   
280   Bool_t success=kFALSE;
281   SetMagField( pEvent->GetMagneticField() );
282   AliESDtrack* track=NULL;
283   
284   for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
285   {
286     track = pEvent->GetTrack(i);
287     if (!track) continue;
288     if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
289          ((track->GetStatus()&AliESDtrack::kITSrefit)>0)&&
290          (track->GetNcls(0)>=fMinPointsVol1)&&
291          (track->GetNcls(1)>=fMinPointsVol2) )
292     { 
293       success = ( AddESDtrack( track ) || success );
294     }
295   }
296   if (success)
297   {
298     fTimeStamp = pEvent->GetTimeStamp();
299     fRunNumber = pEvent->GetRunNumber();
300   }
301   return success;
302 }
303
304 //______________________________________________________________________________
305 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
306 {
307   //Adds a full track, returns true if results in a new estimate
308   //  gets the inner TPC parameters from AliESDTrack::GetInnerParam()
309   //  gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
310
311   const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
312   if (!pconstparamsITS) return kFALSE;
313   const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
314   if (!pconstparamsTPC) return kFALSE;
315   
316   //TPC part
317   AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
318   paramsTPC.Rotate(pconstparamsITS->GetAlpha());
319   paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
320
321   return (AddTrackParams(pconstparamsITS, &paramsTPC));
322 }
323
324 //______________________________________________________________________________
325 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
326 {
327   //Update the estimate using new matching tracklets
328
329   if (!SetTrackParams(p1, p2)) return kFALSE;
330   return Update();
331 }
332
333 //______________________________________________________________________________
334 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
335 {
336   //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
337
338   fNProcessedEvents++; //update the counter
339
340   Bool_t success=kFALSE;
341   TArrayI trackTArrITS(1);
342   TArrayI trackTArrTPC(1);
343   if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
344   SetMagField( pEvent->GetMagneticField() );
345   AliESDtrack* ptrack=NULL;
346   const AliExternalTrackParam* pconstparams1;
347   const AliExternalTrackParam* pconstparams2;
348   AliExternalTrackParam params1;
349   AliExternalTrackParam params2;
350   
351   ////////////////////////////////
352   for (Int_t i=0;i<trackTArrITS.GetSize();i++)
353   {
354     //ITS track
355     ptrack = pEvent->GetTrack(trackTArrITS[i]);
356     pconstparams1 = ptrack->GetOuterParam();
357     if (!pconstparams1) continue;
358     params1 = *pconstparams1; //make copy to be safe
359     
360     //TPC track
361     ptrack = pEvent->GetTrack(trackTArrTPC[i]);
362     pconstparams2 = ptrack->GetInnerParam();
363     if (!pconstparams2) continue;
364     params2 = *pconstparams2; //make copy
365     params2.Rotate(params1.GetAlpha());
366     params2.PropagateTo( params1.GetX(), fMagField );
367
368     if (!SetTrackParams( &params1, &params2 )) continue;
369     
370     //do some accounting and update
371     if (Update())
372       success = kTRUE;
373     else
374       continue;
375   }
376   fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
377   fRunNumber=pEvent->GetRunNumber();
378   return success;
379 }
380
381 //______________________________________________________________________________
382 void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
383 {
384   fNMeasurementParams = (set)?2:4;
385   delete fPH;
386   fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
387   delete fPMeasurement;
388   fPMeasurement = new TVectorD( fNMeasurementParams );
389   delete fPMeasurementCov;
390   fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
391   delete fPMeasurementPrediction;
392   fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
393   fYZOnly = set;
394 }
395
396 //______________________________________________________________________________
397 void AliRelAlignerKalman::Print(Option_t*) const
398 {
399   //Print some useful info
400   Double_t rad2deg = 180./TMath::Pi();
401   printf("\nAliRelAlignerKalman\n");
402   if (fCorrectionMode) printf("(Correction mode)\n");
403   printf("  run: %i, timestamp: %i, magfield: %.3f\n", fRunNumber, fTimeStamp, fMagField);
404   printf("  %i(-%i) inputs, %i(-%i) updates, %i(-%i) merges\n", fNTracks, fNOutliersSigma2Median, fNUpdates, fNOutliers, fNMerges, fNMergesFailed );
405   printf("  psi(x):           % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(0), 1e3*TMath::Sqrt((*fPXcov)(0,0)),(*fPX)(0)*rad2deg,TMath::Sqrt((*fPXcov)(0,0))*rad2deg);
406   printf("  theta(y):         % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(1), 1e3*TMath::Sqrt((*fPXcov)(1,1)),(*fPX)(1)*rad2deg,TMath::Sqrt((*fPXcov)(1,1))*rad2deg);
407   printf("  phi(z):           % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(2), 1e3*TMath::Sqrt((*fPXcov)(2,2)),(*fPX)(2)*rad2deg,TMath::Sqrt((*fPXcov)(2,2))*rad2deg);
408   printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
409   printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
410   printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
411   if (fgkNSystemParams>6) printf("  vd corr           % .5g ± (%.2g)    [ vd should be %.4g (was %.4g in reco) ]\n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)), (*fPX)(6)*fTPCvd, fTPCvd);
412   if (fgkNSystemParams>7) printf("  t0                % .5g ± (%.2g) us  |  %.4g ± (%.2g) cm     [ t0_real = t0_rec+t0 ]\n",(*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)), fTPCvd*(*fPX)(7), fTPCvd*TMath::Sqrt((*fPXcov)(7,7)));
413   if (fgkNSystemParams>8) printf("  vd/dy             % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
414   printf("\n");
415   return;
416 }
417
418 //______________________________________________________________________________
419 void AliRelAlignerKalman::PrintSystemMatrix()
420 {
421   //Print the system matrix for this measurement
422   printf("Kalman system matrix:\n");
423   for ( Int_t i=0; i<fNMeasurementParams; i++ )
424   {
425     for ( Int_t j=0; j<fgkNSystemParams; j++ )
426     {
427       printf("% -2.2f  ", (*fPH)(i,j) );
428     }//for i
429     printf("\n");
430   }//for j
431   printf("\n");
432   return;
433 }
434
435 //______________________________________________________________________________
436 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
437 {
438   //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
439   fNTracks++; //count added input sets
440
441   //INPUT OUTLIER REJECTION
442   if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
443
444   fPTrackParamArr1[fTrackInBuffer] = *exparam1;
445   fPTrackParamArr2[fTrackInBuffer] = *exparam2;
446   
447   fTrackInBuffer++;
448
449   if (fTrackInBuffer == fgkNTracksPerMeasurement)
450   {
451     fTrackInBuffer = 0;
452     return kTRUE;
453   }
454   return kFALSE;
455 }
456
457 //______________________________________________________________________________
458 Bool_t AliRelAlignerKalman::Update()
459 {
460   //perform the update
461   
462   //if (fCalibrationMode) return UpdateCalibration();
463   //if (fFillHistograms)
464   //{
465   //  if (!UpdateEstimateKalman()) return kFALSE;
466   //  return UpdateCalibration(); //Update histograms only when update ok.
467   //}
468   //else return UpdateEstimateKalman();
469   if (!PrepareMeasurement()) return kFALSE;
470   if (!PrepareSystemMatrix()) return kFALSE;
471   if (!PreparePrediction()) return kFALSE;
472   return UpdateEstimateKalman();
473 }
474
475 //______________________________________________________________________________
476 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
477 {
478   //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
479   Double_t sinpsi = TMath::Sin(angles(0));
480   Double_t sintheta = TMath::Sin(angles(1));
481   Double_t sinphi = TMath::Sin(angles(2));
482   Double_t cospsi = TMath::Cos(angles(0));
483   Double_t costheta = TMath::Cos(angles(1));
484   Double_t cosphi = TMath::Cos(angles(2));
485
486   R(0,0) = costheta*cosphi;
487   R(0,1) = -costheta*sinphi;
488   R(0,2) = sintheta;
489   R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
490   R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
491   R(1,2) = -costheta*sinpsi;
492   R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
493   R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
494   R(2,2) = costheta*cospsi;
495 }
496
497 //______________________________________________________________________________
498 Bool_t AliRelAlignerKalman::PrepareMeasurement()
499 {
500   //Calculate the residuals and their covariance matrix
501   
502   for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
503   {
504     const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
505     const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
506
507     //Take the track parameters and calculate the input to the Kalman filter
508     Int_t x = i*4;
509     (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
510     (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
511     if (!fYZOnly)
512     {
513       (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
514       (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
515     }
516
517     //the covariance
518     const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
519     const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
520     (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
521     (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
522     (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
523     (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
524     (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
525     (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
526     (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
527     (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
528     if (!fYZOnly)
529     {
530       (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
531       (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
532       (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
533       (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
534       (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
535       (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
536       (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
537       (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
538       (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
539       (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
540       (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
541       (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
542       (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
543       (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
544       (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
545       (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
546       (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
547       (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
548       (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
549       (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
550       (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
551       (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
552       (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
553       (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
554     }
555     
556   }
557   //if (fApplyCovarianceCorrection)
558   //  *fPMeasurementCov += *fPMeasurementCovCorr;
559   return kTRUE;
560 }
561
562 //______________________________________________________________________________
563 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
564 {
565   //Calculate the system matrix for the Kalman filter
566   //approximate the system using as reference the track in the first volume
567
568   TVectorD z1( fNMeasurementParams );
569   TVectorD z2( fNMeasurementParams );
570   TVectorD x1( fgkNSystemParams );
571   TVectorD x2( fgkNSystemParams );
572   //get the derivatives
573   for ( Int_t i=0; i<fgkNSystemParams; i++ )
574   {
575     x1 = *fPX;
576     x2 = *fPX;
577     x1(i) = x1(i) - fDelta[i]/(2.0);
578     x2(i) = x2(i) + fDelta[i]/(2.0);
579     if (!PredictMeasurement( z1, x1 )) return kFALSE;
580     if (!PredictMeasurement( z2, x2 )) return kFALSE;
581     for (Int_t j=0; j<fNMeasurementParams; j++ )
582     {
583       (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
584     }
585   }
586   return kTRUE;
587 }
588
589 //______________________________________________________________________________
590 Bool_t AliRelAlignerKalman::PreparePrediction()
591 {
592   //Prepare the prediction of the measurement using state vector
593   return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
594 }
595
596 //______________________________________________________________________________
597 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
598 {
599   // Implements a system model for the Kalman fit
600   // pred is [dy,dz,dsinphi,dtgl]
601   // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
602   // note: the measurement is in a local frame, so the prediction also has to be
603   // note: state is the misalignment in global reference system
604
605   if (fCorrectionMode)
606   {
607     for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
608     {
609       AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
610       if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
611       
612       const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
613       const Double_t* newparam = track.GetParameter();
614
615       Int_t x = 4*i;
616       //calculate the predicted residual
617       pred(x+0) = oldparam[0] - newparam[0];
618       pred(x+1) = oldparam[1] - newparam[1];
619       if (!fYZOnly)
620       {
621         pred(x+2) = oldparam[2] - newparam[2];
622         pred(x+3) = oldparam[3] - newparam[3];
623       }
624       return kTRUE;
625     }
626   }
627   else
628   {
629     for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
630     {
631       AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
632       if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
633
634       const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
635       const Double_t* newparam = track.GetParameter();
636
637       Int_t x = 4*i;
638       //calculate the predicted residual
639       pred(x+0) = newparam[0] - oldparam[0];
640       pred(x+1) = newparam[1] - oldparam[1];
641       if (!fYZOnly)
642       {
643         pred(x+2) = newparam[2] - oldparam[2];
644         pred(x+3) = newparam[3] - oldparam[3];
645       }
646       return kTRUE;
647     }
648   }
649   return kFALSE;
650 }
651
652 //______________________________________________________________________________
653 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
654 {
655   //Kalman estimation of noisy constants: in the model A=1
656   //The arguments are (following the usual convention):
657   //  fPX - the state vector (parameters)
658   //  fPXcov - the state covariance matrix (parameter errors)
659   //  fPMeasurement - measurement vector
660   //  fPMeasurementCov - measurement covariance matrix
661   //  fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
662
663   TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov));            //unit matrix
664
665   //predict the state
666   //(*fPXcov) = (*fPXcov) + fQ*identity;  //add some process noise (diagonal)
667
668   // update prediction with measurement
669   // calculate Kalman gain
670   // K = PH/(HPH+fPMeasurementCov)
671   TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) );  //common factor (used twice)
672   TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
673   hpht += (*fPMeasurementCov);
674
675   //shit happens so protect yourself!
676 //  if (fNumericalParanoia)
677 //  {
678 //    TDecompLU lu(hpht);
679 //    if (lu.Condition() > 1e12) return kFALSE;
680 //    lu.Invert(hpht);
681 //  }
682 //  else
683 //  {
684     Double_t det=0.0;
685     hpht.Invert(&det); //since the matrix is small...
686     if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
687 //  }
688   //printf("KalmanUpdate: det(hpht): %.4g\n",det);
689
690   TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
691
692   // update the state and its covariance matrix
693   TVectorD xupdate(fgkNSystemParams);
694   xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
695
696   //SIMPLE OUTLIER REJECTION
697   if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
698   {
699     fNOutliers++;
700     //printf("AliRelAlignerKalman: outlier\n");
701     return kFALSE;
702   }
703
704   TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
705   TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
706   ikh = identity - kh;
707   TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
708   if (!IsPositiveDefinite(ikhp)) return kFALSE;
709
710   (*fPX) += xupdate;
711   TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
712
713   fNUpdates++;
714
715   return kTRUE;
716 }
717
718 //______________________________________________________________________________
719 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
720 {
721   //check whether an update is an outlier given the covariance matrix of the fit
722
723   Bool_t is=kFALSE;
724   for (Int_t i=0;i<fgkNSystemParams;i++)
725   {
726     if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
727     is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
728   }
729   return is;
730 }
731
732 //______________________________________________________________________________
733 Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS, 
734                                                    const AliExternalTrackParam* pTPC )
735 {
736   //check if the input residuals are not too far off their median
737   TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
738   TVectorD vecDeltaN(5);
739   Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
740   vecDeltaN[4]=0;
741   for (Int_t i=0;i<4;i++){
742     vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
743     (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
744   }
745   Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
746   for (Int_t i=0;i<fNMeasurementParams;i++){       //in point2trackmode just take the first 2 params (zy)
747     vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
748     vecRMS[i]    = TMath::RMS(entries,fResArrSigma2Median[i]);
749     vecDeltaN[i] = 0;
750     if (vecRMS[i]>0.){
751       vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
752       vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]);  //sum of abs residuals
753     }
754   }
755   Bool_t outlier=kFALSE;
756   if (fNTracks<3)  outlier=kTRUE;   //median and RMS still to be defined
757   if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
758   if (outlier) fNOutliersSigma2Median++;
759   return outlier;
760 }
761
762 //______________________________________________________________________________
763 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
764 {
765   //check for positive definiteness
766
767   for (Int_t i=0; i<mat.GetNcols(); i++)
768   {
769     if (mat(i,i)<=0.) return kFALSE;
770   }
771
772   if (!fNumericalParanoia) return kTRUE;
773
774   TDecompLU lu(mat);
775   return (lu.Decompose());
776 }
777
778 //______________________________________________________________________________
779 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
780 {
781   //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
782
783   for (Int_t i=0; i<mat.GetNcols(); i++)
784   {
785     matsym(i,i) = mat(i,i); //copy diagonal
786     for (Int_t j=i+1; j<mat.GetNcols(); j++)
787     {
788       //copy the rest
789       Double_t average = (mat(i,j)+mat(j,i))/2.;
790       matsym(i,j)=average;
791       matsym(j,i)=average;
792     }
793   }
794   matsym.MakeValid();
795   return;
796 }
797
798 //______________________________________________________________________________
799 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
800 {
801   //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
802   //b = R*a
803   angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
804   angles(1) = TMath::ASin( rotmat(0,2) );
805   angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
806   return;
807 }
808
809 //______________________________________________________________________________
810 void AliRelAlignerKalman::PrintCorrelationMatrix()
811 {
812   //Print the correlation matrix for the fitted parameters
813   printf("Correlation matrix for system parameters:\n");
814   for ( Int_t i=0; i<fgkNSystemParams; i++ )
815   {
816     for ( Int_t j=0; j<i+1; j++ )
817     {
818       if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf("   NaN  ");
819       else
820         printf("% -1.3f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
821     }//for j
822     printf("\n");
823   }//for i
824   printf("\n");
825   return;
826 }
827
828 //______________________________________________________________________________
829 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
830 {
831   //Find matching track segments in an event with tracks in TPC and ITS(standalone)
832
833   //Sanity cuts on tracks + check which tracks are ITS which are TPC
834   Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
835   fMagField = pEvent->GetMagneticField();
836   if (ntracks<2)
837   {
838     //printf("TrackFinder: less than 2 tracks!\n");
839     return kFALSE;
840   }
841   Float_t* phiArr = new Float_t[ntracks];
842   Float_t* thetaArr = new Float_t[ntracks];
843   Int_t* goodtracksArr = new Int_t[ntracks];
844   Int_t* candidateTPCtracksArr = new Int_t[ntracks];
845   Int_t* matchedITStracksArr = new Int_t[ntracks];
846   Int_t* matchedTPCtracksArr = new Int_t[ntracks];
847   Int_t* tracksArrITS = new Int_t[ntracks];
848   Int_t* tracksArrTPC = new Int_t[ntracks];
849   Int_t* nPointsArr = new Int_t[ntracks];
850   Int_t nITStracks = 0;
851   Int_t nTPCtracks = 0;
852   Int_t nGoodTracks = 0;
853   Int_t nCandidateTPCtracks = 0;
854   Int_t nMatchedITStracks = 0;
855   AliESDtrack* pTrack = NULL;
856   Bool_t foundMatchTPC = kFALSE;
857
858   //select and clasify tracks
859   for (Int_t itrack=0; itrack < ntracks; itrack++)
860   {
861     pTrack = pEvent->GetTrack(itrack);
862     if (!pTrack)
863     {
864       //std::cout<<"no track!"<<std::endl;
865       continue;
866     }
867     if (fCuts)
868     {
869       if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
870     }
871     goodtracksArr[nGoodTracks]=itrack;
872     Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
873     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
874     phiArr[nGoodTracks]=phi;
875     thetaArr[nGoodTracks]=theta;
876
877     //check if track is ITS
878     Int_t nClsITS = pTrack->GetNcls(0);
879     Int_t nClsTPC = pTrack->GetNcls(1);
880     if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
881          !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
882          !(nClsITS<fMinPointsVol1) )  //enough points
883     {
884       tracksArrITS[nITStracks] = nGoodTracks;
885       nITStracks++;
886       nGoodTracks++;
887       continue;
888     }
889
890     //check if track is TPC
891     if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
892          !(nClsTPC<fMinPointsVol2) )  //enough points
893     {
894       tracksArrTPC[nTPCtracks] = nGoodTracks;
895       nTPCtracks++;
896       nGoodTracks++;
897       //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
898       continue;
899     }
900   }//for itrack   -   selection fo tracks
901
902   //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
903
904   if ( nITStracks < 1 || nTPCtracks < 1 )
905   {
906     delete [] phiArr;
907     delete [] thetaArr;
908     delete [] goodtracksArr;
909     delete [] matchedITStracksArr;
910     delete [] candidateTPCtracksArr;
911     delete [] matchedTPCtracksArr;
912     delete [] tracksArrITS;
913     delete [] tracksArrTPC;
914     delete [] nPointsArr;
915     return kFALSE;
916   }
917
918   //find matching in TPC
919   if (nTPCtracks>1)  //if there is something to be matched, try and match it
920   {
921     Float_t min = 10000000.;
922     for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
923     {
924       for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
925       {
926         Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
927         if (deltatheta > fMaxMatchingAngle) continue;
928         Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
929         if (deltaphi > fMaxMatchingAngle) continue;
930         if (deltatheta+deltaphi<min) //only the best matching pair
931         {
932           min=deltatheta+deltaphi;
933           candidateTPCtracksArr[0] = tracksArrTPC[itr1];  //store the index of track in goodtracksArr[]
934           candidateTPCtracksArr[1] = tracksArrTPC[itr2];
935           nCandidateTPCtracks = 2;
936           foundMatchTPC = kTRUE;
937           //printf("TrackFinder: Matching TPC tracks candidates:\n");
938           //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
939           //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
940         }
941       }
942     }
943   }//if nTPCtracks>1
944   else //if nTPCtracks==1 - if nothing to match, take the only one we've got
945   {
946     candidateTPCtracksArr[0] = tracksArrTPC[0];
947     nCandidateTPCtracks = 1;
948     foundMatchTPC = kFALSE;
949   }
950   if (foundMatchTPC) fNMatchedTPCtracklets++;
951   //if no match but the requirement is set return kFALSE
952   if (fRequireMatchInTPC && !foundMatchTPC)
953   {
954     delete [] phiArr;
955     delete [] thetaArr;
956     delete [] goodtracksArr;
957     delete [] candidateTPCtracksArr;
958     delete [] matchedITStracksArr;
959     delete [] matchedTPCtracksArr;
960     delete [] tracksArrITS;
961     delete [] tracksArrTPC;
962     delete [] nPointsArr;
963     //printf("TrackFinder: no match in TPC && required\n");
964     return kFALSE;
965   }
966
967   //if no match and more than one track take all TPC tracks
968   if (!fRequireMatchInTPC && !foundMatchTPC)
969   {
970     for (Int_t i=0;i<nTPCtracks;i++)
971     {
972       candidateTPCtracksArr[i] = tracksArrTPC[i];
973     }
974     nCandidateTPCtracks = nTPCtracks;
975   }
976   //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
977
978   Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
979
980   //find ITS matches for good TPC tracks
981   Bool_t matchedITStracks=kFALSE;
982   for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
983   {
984     minDifferenceArr[nMatchedITStracks] = 10000000.;
985     matchedITStracks=kFALSE;
986     for (Int_t iits=0; iits<nITStracks; iits++)
987     {
988       AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
989       const AliExternalTrackParam* parits = itstrack->GetOuterParam();
990       AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
991       const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
992       AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
993       partpc.Rotate(parits->GetAlpha());
994       partpc.PropagateTo(parits->GetX(),fMagField);
995       Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
996       if (dtgl > fMaxMatchingAngle) continue;
997       Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
998       if (dsnp > fMaxMatchingAngle) continue;
999       Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
1000       Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
1001       if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
1002       if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
1003       {
1004         minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
1005         matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
1006         matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
1007         //printf("TrackFinder: Matching ITS to TPC:\n");
1008         //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
1009         //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
1010         //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
1011         matchedITStracks=kTRUE;;
1012       }
1013     }
1014     if (matchedITStracks) nMatchedITStracks++;
1015   }
1016
1017   if (nMatchedITStracks==0) //no match in ITS
1018   {
1019     delete [] phiArr;
1020     delete [] thetaArr;
1021     delete [] minDifferenceArr;
1022     delete [] goodtracksArr;
1023     delete [] matchedITStracksArr;
1024     delete [] candidateTPCtracksArr;
1025     delete [] matchedTPCtracksArr;
1026     delete [] tracksArrITS;
1027     delete [] tracksArrTPC;
1028     delete [] nPointsArr;
1029     //printf("TrackFinder: No match in ITS\n");
1030     return kFALSE;
1031   }
1032
1033   //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1034   //we found a cosmic
1035   fNMatchedCosmics++;
1036
1037   //Now we may have ended up with more matches than we want in the case there was
1038   //no TPC match and there were many TPC tracks
1039   //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1040   //so take only the best pair.
1041   //same applies when there are more matches than ITS tracks - means that one ITS track
1042   //matches more TPC tracks.
1043   if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1044   {
1045     Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1046     matchedITStracksArr[0] = matchedITStracksArr[imin];
1047     matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1048     nMatchedITStracks = 1;
1049     //printf("TrackFinder: too many matches - take only the best one\n");
1050     //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1051     //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1052     //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1053   }
1054
1055   ///////////////////////////////////////////////////////////////////////////
1056   outITSindexTArr.Set(nMatchedITStracks);
1057   outTPCindexTArr.Set(nMatchedITStracks);
1058   for (Int_t i=0;i<nMatchedITStracks;i++)
1059   {
1060     outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1061     outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1062     //printf("TrackFinder: Fill the output\n");
1063     //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1064     //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1065   }
1066   //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1067   ///////////////////////////////////////////////////////////////////////////
1068
1069   delete [] phiArr;
1070   delete [] thetaArr;
1071   delete [] minDifferenceArr;
1072   delete [] goodtracksArr;
1073   delete [] candidateTPCtracksArr;
1074   delete [] matchedITStracksArr;
1075   delete [] matchedTPCtracksArr;
1076   delete [] tracksArrITS;
1077   delete [] tracksArrTPC;
1078   delete [] nPointsArr;
1079   return kTRUE;
1080 }
1081
1082 //______________________________________________________________________________
1083 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1084 {
1085   //implements the system model -
1086   //applies correction for misalignment and calibration to track
1087   //track needs to be already propagated to the global reference plane
1088
1089   Double_t x = tr->GetX();
1090   Double_t alpha = tr->GetAlpha();
1091   Double_t point[3],dir[3];
1092   tr->GetXYZ(point);
1093   tr->GetDirection(dir);
1094   TVector3 Point(point);
1095   TVector3 Dir(dir);
1096   
1097   //Apply corrections to track
1098
1099   //Shift
1100   Point(0) -= misal(3); //add shift in x
1101   Point(1) -= misal(4); //add shift in y
1102   Point(2) -= misal(5); //add shift in z
1103   //Rotation
1104   TMatrixD rotmat(3,3);
1105   RotMat( rotmat, misal );
1106   Point = rotmat.T() * Point;
1107   Dir = rotmat * Dir;
1108   
1109   //TPC vdrift and T0 corrections
1110   TVector3 Point2; //second point of the track
1111   Point2 = Point + Dir;
1112   Double_t vdCorr = 1./misal(6);
1113   Double_t t0 = misal(7);
1114   Double_t vdY = 0.0;
1115   if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1116
1117   //my model
1118   if (Point(2)>0)
1119   {
1120     //A-Side
1121     Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (vdCorr-1.+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1122     Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1123   }
1124   else
1125   {
1126     //C-side
1127     Point(2) = Point(2)   - (fTPCZLengthC+Point(2))  * (1.-vdCorr-vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1128     Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1129   }
1130
1131   //Stefan's model
1132   //if (Point(2)>0)
1133   //{
1134   //  //A-Side
1135   //  Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (1.-vdCorr+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1136   //  Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1137   //}
1138   //else
1139   //{
1140   //  //C-side
1141   //  Point(2) = Point(2)   + (fTPCZLengthC+Point(2))  * (1.-vdCorr+vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1142   //  Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1143   //}
1144
1145   Dir = Point2-Point;
1146   Dir=Dir.Unit(); //keep unit length
1147
1148   //Turn back to local system
1149   Dir.GetXYZ(dir);
1150   Point.GetXYZ(point);
1151   tr->Global2LocalPosition(point,alpha);
1152   tr->Global2LocalPosition(dir,alpha);
1153
1154   //Calculate new intersection point with ref plane
1155   Double_t p[5],pcov[15];
1156   if (dir[0]==0.0) return kFALSE;
1157   Double_t s=(x-point[0])/dir[0];
1158   p[0] = point[1]+s*dir[1];
1159   p[1] = point[2]+s*dir[2];
1160   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1161   if (pt==0.0) return kFALSE;
1162   p[2] = dir[1]/pt;
1163   p[3] = dir[2]/pt;
1164   //insert everything back into track
1165   const Double_t* pcovtmp = tr->GetCovariance();
1166   p[4] = tr->GetSigned1Pt(); //copy the momentum
1167   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1168   tr->Set(x,alpha,p,pcov);
1169   return kTRUE;
1170
1171   ////put params back into track and propagate to ref
1172   //Double_t p[5],pcov[15];
1173   //p[0] = point[1];
1174   //p[1] = point[2];
1175   //Double_t xnew = point[0];
1176   //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1177   //if (pt==0.0) return kFALSE;
1178   //p[2] = dir[1]/pt;
1179   //p[3] = dir[2]/pt;
1180   //p[4] = tr->GetSigned1Pt(); //copy the momentum
1181   //const Double_t* pcovtmp = tr->GetCovariance();
1182   //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1183   //tr->Set(xnew,alpha,p,pcov);
1184   //return tr->PropagateTo(x,fMagField);
1185 }
1186
1187 //______________________________________________________________________________
1188 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1189 {
1190   //implements the system model -
1191   //applies misalignment and miscalibration to reference track
1192   //trackparams have to be at the global reference plane
1193
1194   Double_t x = tr->GetX();
1195   Double_t alpha = tr->GetAlpha();
1196   Double_t point[3],dir[3];
1197   tr->GetXYZ(point);
1198   tr->GetDirection(dir);
1199   TVector3 Point(point);
1200   TVector3 Dir(dir);
1201   
1202   //Apply misalignment to track
1203   
1204   //TPC vdrift and T0 corrections
1205   TVector3 Point2; //second point of the track
1206   Point2 = Point + Dir;
1207   Double_t vdCorr = 1./misal(6);
1208   Double_t t0 = misal(7);
1209   Double_t vdY = 0.0;
1210   if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1211
1212   if (Point(2)>0)
1213   {
1214     //A-Side
1215     Point(2) = Point(2)   + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1216                           * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1217     Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1218                           * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1219   }
1220   else
1221   {
1222     //C-side
1223     Point(2) = Point(2)   + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1224                           * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1225     Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1226                           * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1227   }
1228   Dir = Point2-Point;
1229   Dir=Dir.Unit(); //keep unit length
1230
1231   //Rotation
1232   TMatrixD rotmat(3,3);
1233   RotMat( rotmat, misal );
1234   Point = rotmat * Point;
1235   Dir = rotmat * Dir;
1236   //Shift
1237   Point(0) += misal(3); //add shift in x
1238   Point(1) += misal(4); //add shift in y
1239   Point(2) += misal(5); //add shift in z
1240
1241   //Turn back to local system
1242   Dir.GetXYZ(dir);
1243   Point.GetXYZ(point);
1244   tr->Global2LocalPosition(point,alpha);
1245   tr->Global2LocalPosition(dir,alpha);
1246
1247   //Calculate new intersection point with ref plane
1248   Double_t p[5],pcov[15];
1249   if (dir[0]==0.0) return kFALSE;
1250   Double_t s=(x-point[0])/dir[0];
1251   p[0] = point[1]+s*dir[1];
1252   p[1] = point[2]+s*dir[2];
1253   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1254   if (pt==0.0) return kFALSE;
1255   p[2] = dir[1]/pt;
1256   p[3] = dir[2]/pt;
1257   //insert everything back into track
1258   const Double_t* pcovtmp = tr->GetCovariance();
1259   p[4] = tr->GetSigned1Pt(); //copy the momentum
1260   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1261   tr->Set(x,alpha,p,pcov);
1262   return kTRUE;
1263
1264   ////put params back into track and propagate to ref
1265   //Double_t p[5];
1266   //Double_t pcov[15];
1267   //p[0] = point[1];
1268   //p[1] = point[2];
1269   //Double_t xnew = point[0];
1270   //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1271   //if (pt==0.0) return kFALSE;
1272   //p[2] = dir[1]/pt;
1273   //p[3] = dir[2]/pt;
1274   //p[4] = tr->GetSigned1Pt(); //copy the momentum
1275   //const Double_t* pcovtmp = tr->GetCovariance();
1276   //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1277   //printf("x before: %.5f, after: %.5f\n",x, xnew);
1278   //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1279   //printf("after:  %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1280   //tr->Set(xnew,alpha,p,pcov);
1281   //return tr->PropagateTo(x,fMagField);
1282 }
1283
1284 //______________________________________________________________________________
1285 void AliRelAlignerKalman::Reset()
1286 {
1287   //full reset to defaults
1288   fPX->Zero();
1289   (*fPX)(6)=1.;
1290   ResetCovariance();
1291
1292   //initialize the differentials per parameter
1293   for (Int_t i=0;i<4;i++) 
1294   {
1295     delete [] (fResArrSigma2Median[i]);
1296   }
1297   fRejectOutliersSigma2Median=kFALSE;
1298
1299   fNMatchedCosmics=0;
1300   fNMatchedTPCtracklets=0;
1301   fNUpdates=0;
1302   fNOutliers=0;
1303   fNTracks=0;
1304   fNProcessedEvents=0;
1305   fRunNumber=0;
1306   fTimeStamp=0;
1307 }
1308
1309 //______________________________________________________________________________
1310 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1311 {
1312   //Resets the covariance to the default if arg=0 or resets the off diagonals
1313   //to zero and releases the diagonals by factor arg.
1314   if (number!=0.)
1315   {
1316     for (Int_t z=0;z<6;z++)
1317     {
1318       for (Int_t zz=0;zz<6;zz++)
1319       {
1320         if (zz==z) continue; //don't touch diagonals
1321         (*fPXcov)(zz,z) = 0.;
1322         (*fPXcov)(z,zz) = 0.;
1323       }
1324       (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1325     }
1326   }
1327   else
1328   {
1329     //Resets the covariance of the fit to a default value
1330     fPXcov->Zero();
1331     (*fPXcov)(0,0) = .08*.08; //psi (rad)
1332     (*fPXcov)(1,1) = .08*.08; //theta (rad
1333     (*fPXcov)(2,2) = .08*.08; //phi (rad)
1334     (*fPXcov)(3,3) = .3*.3; //x (cm)
1335     (*fPXcov)(4,4) = .3*.3; //y (cm)
1336     (*fPXcov)(5,5) = .3*.3; //z (cm)
1337   }
1338   ResetTPCparamsCovariance(number); 
1339 }
1340
1341 //______________________________________________________________________________
1342 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1343 {
1344   //Resets the covariance to the default if arg=0 or resets the off diagonals
1345   //to zero and releases the diagonals by factor arg.
1346   
1347   //release diagonals
1348   if (number==0.)
1349   {
1350     if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1351     if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1352     if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1353   }
1354   else
1355   {
1356     if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1357     if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1358     if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1359   }
1360   
1361   //set crossterms to zero
1362   for (Int_t i=0;i<fgkNSystemParams;i++)
1363   {
1364     for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1365     {
1366       if (i==j) continue; //don't touch diagonals
1367       (*fPXcov)(i,j) = 0.;
1368       (*fPXcov)(j,i) = 0.;
1369     }
1370   }
1371 }
1372
1373 //______________________________________________________________________________
1374 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1375 {
1376   //Merge two aligners
1377   
1378   if (!al) return kFALSE;
1379   if (al==this) return kTRUE;
1380   if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1381   
1382   //store the pointers to current stuff
1383   TVectorD* pmes = fPMeasurement;
1384   TMatrixDSym* pmescov = fPMeasurementCov;
1385   TVectorD* pmespred = fPMeasurementPrediction;
1386   TMatrixD* ph = fPH;
1387
1388   //make a unity system matrix
1389   TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1390   fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1391
1392   //mesurement is the state of the new aligner
1393   fPMeasurement = al->fPX;
1394   fPMeasurementCov = al->fPXcov;
1395
1396   //the mesurement prediction is the state
1397   fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1398   
1399   //do the merging
1400   Bool_t success = UpdateEstimateKalman();
1401   
1402   //restore pointers to old stuff
1403   fPMeasurement = pmes;
1404   fPMeasurementCov = pmescov;
1405   fPMeasurementPrediction = pmespred;
1406   delete fPH;
1407   fPH = ph;
1408
1409   //merge stats
1410   if (!success)    
1411   {
1412     fNMergesFailed++;
1413     //printf("AliRelAlignerKalman::Merge failed\n");
1414     return kFALSE; //no point in merging stats if merge not succesful
1415   }
1416   fNProcessedEvents += al->fNProcessedEvents;
1417   fNUpdates += al->fNUpdates;
1418   fNOutliers += al->fNOutliers;
1419   fNOutliersSigma2Median += al->fNOutliersSigma2Median;
1420   fNTracks += al->fNTracks;
1421   fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1422   fNMatchedCosmics += al->fNMatchedCosmics;
1423   if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1424   else fNMerges += al->fNMerges;
1425   if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
1426
1427   return success;
1428 }
1429
1430 //______________________________________________________________________________
1431 Long64_t AliRelAlignerKalman::Merge( TCollection* list )
1432 {
1433   //merge all aligners in the collection
1434   Long64_t numberOfMerges=0;
1435   AliRelAlignerKalman* alignerFromList;
1436   if (!list) return 0;
1437   TIter next(list);
1438   while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
1439   {
1440     if (alignerFromList == this) continue;
1441     if (Merge(alignerFromList)) numberOfMerges++;
1442   }
1443   return numberOfMerges;
1444 }
1445
1446 //______________________________________________________________________________
1447 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1448 {
1449   if (this == obj) return 0;
1450   const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1451   if (!aobj) return 0;
1452   if (fTimeStamp < aobj->fTimeStamp) return -1;
1453   else if (fTimeStamp > aobj->fTimeStamp) return 1;
1454   else return 0;
1455 }
1456
1457 //________________________________________________________________________
1458 Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
1459 {
1460   //find matching tracks and return tobjarrays with the params
1461   
1462   Int_t ntracks = pESD->GetNumberOfTracks();
1463   Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
1464   for (Int_t i=0;i<ntracks;i++)
1465   {
1466     matchedArr[i]=0;
1467   }
1468
1469   Int_t iMatched=-1;
1470   for (Int_t i=0; i<ntracks; i++)
1471   {
1472     //get track1 and friend
1473     AliESDtrack* track1 = pESD->GetTrack(i);
1474     if (!track1) continue;
1475
1476     if (track1->GetNcls(0) < fMinPointsVol1) continue;
1477
1478     if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
1479            (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
1480
1481     const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
1482     if (!constfriendtrack1) continue;
1483     AliESDfriendTrack friendtrack1(*constfriendtrack1);
1484     
1485     if (!friendtrack1.GetITSOut()) continue;
1486     AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
1487
1488     Double_t bestd = 1000.; //best distance
1489     Bool_t newi = kTRUE; //whether we start with a new i
1490     for (Int_t j=0; j<ntracks; j++)
1491     {
1492       if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried 
1493       //get track2 and friend
1494       AliESDtrack* track2 = pESD->GetTrack(j);
1495       if (!track2) continue;
1496       if (track1==track2) continue;
1497       //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
1498       //       (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
1499
1500       if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
1501       if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
1502       if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
1503       if (track2->GetTgl() > 1.) continue; //acceptance
1504       //cut crossing tracks
1505       if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
1506       if (track2->GetInnerParam()->GetX()>90) continue;
1507       if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
1508
1509
1510       if (!track2->GetInnerParam()) continue;
1511       AliExternalTrackParam params2(*(track2->GetInnerParam()));
1512
1513       //bring to same reference plane
1514       if (!params2.Rotate(params1.GetAlpha())) continue;
1515       if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
1516
1517       //pt cut
1518       if (params2.Pt() < fMinPt) continue;
1519
1520       const Double32_t* p1 = params1.GetParameter();
1521       const Double32_t* p2 = params2.GetParameter();
1522
1523       //hard cuts
1524       Double_t dy = TMath::Abs(p2[0]-p1[0]);
1525       Double_t dz = TMath::Abs(p2[1]-p1[1]);
1526       Double_t dphi = TMath::Abs(p2[2]-p1[2]);
1527       Double_t dlam = TMath::Abs(p2[3]-p1[3]);
1528       if (dy > 2.0) continue;
1529       if (dz > 10.0) continue;
1530       if (dphi > 0.1 ) continue;
1531       if (dlam > 0.1 ) continue;
1532
1533       //best match only
1534       Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
1535       if ( d >= bestd) continue;
1536       bestd = d;
1537       matchedArr[j]=i; //j-th track matches i-th (ITS) track
1538       if (newi) iMatched++; newi=kFALSE; //increment at most once per i
1539       if (arrITS[iMatched] && arrTPC[iMatched])
1540       {
1541         *(arrITS[iMatched]) = params1;
1542         *(arrTPC[iMatched]) = params2;
1543       }
1544       else
1545       {
1546         arrITS[iMatched] = new AliExternalTrackParam(params1);
1547         arrTPC[iMatched] = new AliExternalTrackParam(params2);
1548       }//else
1549     }//for j
1550   }//for i
1551   delete [] matchedArr;
1552   return iMatched;
1553 }
1554
1555 //________________________________________________________________________
1556 void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
1557 {
1558   //Sets up or destroys the memory hungry array to hold the statistics
1559   //for data rejection with median
1560   if (set)
1561   {
1562     for (Int_t i=0;i<4;i++) 
1563     {
1564       if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] = 
1565                                    new Double_t[fgkNtracksSigma2Median];
1566     }
1567     fRejectOutliersSigma2Median = kTRUE;
1568   }//else
1569   else
1570   {
1571     // it probably doesn't make sense to delete the arrays, they are not streamed
1572     //if (fRejectOutliersSigma2Median)
1573     //for (Int_t i=0;i<4;i++) 
1574     //{
1575     //  delete [] (fResArrSigma2Median[i]);
1576     //}
1577     fRejectOutliersSigma2Median = kFALSE;
1578   }//if
1579 }