Coverity fixes.
[u/mrichter/AliRoot.git] / STEER / AliRelAlignerKalman.cxx
CommitLineData
043badeb 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:
044eb03e 19// Finds alignement constants for two tracking volumes (by default ITS
20// and TPC)
55d5da9e 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.
044eb03e 24// Additionally calculates some callibration parameters for TPC
043badeb 25// Fit parameters are:
26// - 3 shifts, x,y,z
27// - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
c3b5bfc1 28// - TPC drift velocity correction,
55d5da9e 29// - TPC time offset correction.
043badeb 30//
31// Basic usage:
3ebfbf52 32// When aligning two volumes, at any given time a single instance of
043badeb 33// the class should be active. The fit of the parameters is updated
55d5da9e 34// by adding new data using one of the Add.... methods:
043badeb 35//
3ebfbf52 36// In collision events add an ESD event to update the fit (adds all tracks):
55d5da9e 37//
3ebfbf52 38// Bool_t AddESDevent( AliESDevent* pTrack );
39//
40// or add each individual track
41//
42// AddESDtrack( AliESDtrack* pTrack );
55d5da9e 43//
044eb03e 44// For cosmic data, the assumption is that the tracking is done twice:
55d5da9e 45// once global and once only ITS and the tracklets are saved inside
044eb03e 46// one AliESDEvent. The method
55d5da9e 47//
48// Bool_t AddCosmicEvent( AliESDEvent* pEvent );
49//
044eb03e 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)
55d5da9e 53//
3ebfbf52 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//
044eb03e 57// _________________________________________________________________________
58// Expert options:
3ebfbf52 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.
043badeb 62//
63// Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
64//
65//////////////////////////////////////////////////////////////////////////////
66
55d5da9e 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>
fffc83e6 75#include <TObjArray.h>
55d5da9e 76#include <TH1D.h>
77#include <TF1.h>
78
79#include "AliESDtrack.h"
55d5da9e 80#include "AliESDEvent.h"
55d5da9e 81#include "AliExternalTrackParam.h"
82
043badeb 83#include "AliRelAlignerKalman.h"
84
85ClassImp(AliRelAlignerKalman)
86
044eb03e 87//______________________________________________________________________________
043badeb 88AliRelAlignerKalman::AliRelAlignerKalman():
3ebfbf52 89 TObject(),
3ebfbf52 90 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
91 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
92 fMagField(0.),
16cbbc2c 93 fNMeasurementParams(4),
782e5230 94 fPX(new TVectorD( fgkNSystemParams )),
95 fPXcov(new TMatrixDSym( fgkNSystemParams )),
16cbbc2c 96 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
55d5da9e 97 fQ(1.e-15),
16cbbc2c 98 fPMeasurement(new TVectorD( fNMeasurementParams )),
99 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
100 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
044eb03e 101 fOutRejSigmas(1.),
fffc83e6 102 fOutRejSigma2Median(5.),
16cbbc2c 103 fYZOnly(kFALSE),
3ebfbf52 104 fNumericalParanoia(kTRUE),
044eb03e 105 fRejectOutliers(kTRUE),
fffc83e6 106 fRejectOutliersSigma2Median(kFALSE),
c3b5bfc1 107 fRequireMatchInTPC(kFALSE),
043badeb 108 fCuts(kFALSE),
c3b5bfc1 109 fMinPointsVol1(3),
044eb03e 110 fMinPointsVol2(50),
fffc83e6 111 fMinPt(0.),
112 fMaxPt(1.e100),
c3b5bfc1 113 fMaxMatchingAngle(0.1),
55d5da9e 114 fMaxMatchingDistance(10.), //in cm
115 fCorrectionMode(kFALSE),
044eb03e 116 fNTracks(0),
117 fNUpdates(0),
118 fNOutliers(0),
fffc83e6 119 fNOutliersSigma2Median(0),
044eb03e 120 fNMatchedCosmics(0),
c3b5bfc1 121 fNMatchedTPCtracklets(0),
044eb03e 122 fNProcessedEvents(0),
3ebfbf52 123 fTrackInBuffer(0),
a80268df 124 fTimeStamp(0),
329b5744 125 fRunNumber(0),
52d1b7b0 126 fNMerges(0),
127 fNMergesFailed(0),
55d5da9e 128 fTPCvd(2.64),
3ebfbf52 129 fTPCZLengthA(2.4972500e02),
130 fTPCZLengthC(2.4969799e02)
043badeb 131{
55d5da9e 132 //Default constructor
fffc83e6 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;}
55d5da9e 135 Reset();
55d5da9e 136}
137
138//______________________________________________________________________________
139AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
3ebfbf52 140 TObject(static_cast<TObject>(a)),
3ebfbf52 141 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
142 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
143 fMagField(a.fMagField),
16cbbc2c 144 fNMeasurementParams(a.fNMeasurementParams),
55d5da9e 145 fPX(new TVectorD( *a.fPX )),
146 fPXcov(new TMatrixDSym( *a.fPXcov )),
52d1b7b0 147 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
55d5da9e 148 fQ(a.fQ),
52d1b7b0 149 fPMeasurement(new TVectorD( fNMeasurementParams )),
150 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
151 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
55d5da9e 152 fOutRejSigmas(a.fOutRejSigmas),
fffc83e6 153 fOutRejSigma2Median(a.fOutRejSigma2Median),
16cbbc2c 154 fYZOnly(a.fYZOnly),
3ebfbf52 155 fNumericalParanoia(a.fNumericalParanoia),
55d5da9e 156 fRejectOutliers(a.fRejectOutliers),
fffc83e6 157 fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
55d5da9e 158 fRequireMatchInTPC(a.fRequireMatchInTPC),
55d5da9e 159 fCuts(a.fCuts),
160 fMinPointsVol1(a.fMinPointsVol1),
161 fMinPointsVol2(a.fMinPointsVol2),
fffc83e6 162 fMinPt(a.fMinPt),
163 fMaxPt(a.fMaxPt),
55d5da9e 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),
fffc83e6 170 fNOutliersSigma2Median(a.fNOutliersSigma2Median),
55d5da9e 171 fNMatchedCosmics(a.fNMatchedCosmics),
172 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
173 fNProcessedEvents(a.fNProcessedEvents),
52d1b7b0 174 fTrackInBuffer(0),
a80268df 175 fTimeStamp(a.fTimeStamp),
329b5744 176 fRunNumber(a.fRunNumber),
52d1b7b0 177 fNMerges(a.fNMerges),
178 fNMergesFailed(a.fNMergesFailed),
55d5da9e 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));
fffc83e6 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 }
55d5da9e 198}
199
200//______________________________________________________________________________
201AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
202{
203 //assignment operator
3ebfbf52 204 fMagField=a.fMagField,
52d1b7b0 205 fNMeasurementParams=a.fNMeasurementParams;
55d5da9e 206 *fPX = *a.fPX;
207 *fPXcov = *a.fPXcov;
55d5da9e 208 fQ=a.fQ;
55d5da9e 209 fOutRejSigmas=a.fOutRejSigmas;
fffc83e6 210 fOutRejSigma2Median=a.fOutRejSigma2Median;
55d5da9e 211 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
16cbbc2c 212 fYZOnly=a.fYZOnly;
3ebfbf52 213 fNumericalParanoia=a.fNumericalParanoia;
55d5da9e 214 fRejectOutliers=a.fRejectOutliers;
fffc83e6 215 fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
55d5da9e 216 fRequireMatchInTPC=a.fRequireMatchInTPC;
55d5da9e 217 fCuts=a.fCuts;
218 fMinPointsVol1=a.fMinPointsVol1;
219 fMinPointsVol2=a.fMinPointsVol2;
fffc83e6 220 fMinPt=a.fMinPt;
221 fMaxPt=a.fMaxPt;
55d5da9e 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;
fffc83e6 228 fNOutliersSigma2Median=a.fNOutliersSigma2Median;
55d5da9e 229 fNMatchedCosmics=a.fNMatchedCosmics;
230 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
231 fNProcessedEvents=a.fNProcessedEvents;
52d1b7b0 232 fTrackInBuffer=0; //because the array is reset, not copied
a80268df 233 fTimeStamp=a.fTimeStamp;
329b5744 234 fRunNumber=a.fRunNumber;
52d1b7b0 235 fNMerges=a.fNMerges;
55d5da9e 236 fTPCvd=a.fTPCvd;
237 fTPCZLengthA=a.fTPCZLengthA;
238 fTPCZLengthC=a.fTPCZLengthC;
fffc83e6 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 }
55d5da9e 253 return *this;
043badeb 254}
255
044eb03e 256//______________________________________________________________________________
257AliRelAlignerKalman::~AliRelAlignerKalman()
043badeb 258{
55d5da9e 259 //destructor
fffc83e6 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 }
3ebfbf52 271}
272
273//______________________________________________________________________________
274Bool_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() );
da6f4665 282 AliESDtrack* track=NULL;
3ebfbf52 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)&&
fffc83e6 289 ((track->GetStatus()&AliESDtrack::kITSrefit)>0)&&
3ebfbf52 290 (track->GetNcls(0)>=fMinPointsVol1)&&
291 (track->GetNcls(1)>=fMinPointsVol2) )
292 {
293 success = ( AddESDtrack( track ) || success );
294 }
295 }
329b5744 296 if (success)
297 {
298 fTimeStamp = pEvent->GetTimeStamp();
299 fRunNumber = pEvent->GetRunNumber();
300 }
3ebfbf52 301 return success;
043badeb 302}
303
044eb03e 304//______________________________________________________________________________
3ebfbf52 305Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
043badeb 306{
a80268df 307 //Adds a full track, returns true if results in a new estimate
308 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
fffc83e6 309 // gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
a80268df 310
311 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
312 if (!pconstparamsITS) return kFALSE;
313 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
314 if (!pconstparamsTPC) return kFALSE;
3ebfbf52 315
a80268df 316 //TPC part
317 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
318 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
319 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
320
321 return (AddTrackParams(pconstparamsITS, &paramsTPC));
043badeb 322}
323
18f65700 324//______________________________________________________________________________
325Bool_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
044eb03e 333//______________________________________________________________________________
55d5da9e 334Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
043badeb 335{
55d5da9e 336 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
337
3ebfbf52 338 fNProcessedEvents++; //update the counter
339
55d5da9e 340 Bool_t success=kFALSE;
341 TArrayI trackTArrITS(1);
342 TArrayI trackTArrTPC(1);
343 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
3ebfbf52 344 SetMagField( pEvent->GetMagneticField() );
da6f4665 345 AliESDtrack* ptrack=NULL;
3ebfbf52 346 const AliExternalTrackParam* pconstparams1;
347 const AliExternalTrackParam* pconstparams2;
348 AliExternalTrackParam params1;
349 AliExternalTrackParam params2;
350
55d5da9e 351 ////////////////////////////////
352 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
353 {
354 //ITS track
355 ptrack = pEvent->GetTrack(trackTArrITS[i]);
3ebfbf52 356 pconstparams1 = ptrack->GetOuterParam();
357 if (!pconstparams1) continue;
358 params1 = *pconstparams1; //make copy to be safe
044eb03e 359
55d5da9e 360 //TPC track
361 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
3ebfbf52 362 pconstparams2 = ptrack->GetInnerParam();
363 if (!pconstparams2) continue;
364 params2 = *pconstparams2; //make copy
365 params2.Rotate(params1.GetAlpha());
a80268df 366 params2.PropagateTo( params1.GetX(), fMagField );
3ebfbf52 367
368 if (!SetTrackParams( &params1, &params2 )) continue;
55d5da9e 369
370 //do some accounting and update
55d5da9e 371 if (Update())
372 success = kTRUE;
373 else
374 continue;
375 }
329b5744 376 fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
377 fRunNumber=pEvent->GetRunNumber();
55d5da9e 378 return success;
043badeb 379}
380
16cbbc2c 381//______________________________________________________________________________
382void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
383{
384 fNMeasurementParams = (set)?2:4;
fffc83e6 385 delete fPH;
16cbbc2c 386 fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
fffc83e6 387 delete fPMeasurement;
16cbbc2c 388 fPMeasurement = new TVectorD( fNMeasurementParams );
fffc83e6 389 delete fPMeasurementCov;
16cbbc2c 390 fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
fffc83e6 391 delete fPMeasurementPrediction;
16cbbc2c 392 fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
393 fYZOnly = set;
394}
395
044eb03e 396//______________________________________________________________________________
397void AliRelAlignerKalman::Print(Option_t*) const
043badeb 398{
55d5da9e 399 //Print some useful info
400 Double_t rad2deg = 180./TMath::Pi();
3ebfbf52 401 printf("\nAliRelAlignerKalman\n");
402 if (fCorrectionMode) printf("(Correction mode)\n");
fffc83e6 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 );
55d5da9e 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)));
3ebfbf52 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)));
55d5da9e 414 printf("\n");
415 return;
043badeb 416}
417
044eb03e 418//______________________________________________________________________________
419void AliRelAlignerKalman::PrintSystemMatrix()
043badeb 420{
55d5da9e 421 //Print the system matrix for this measurement
422 printf("Kalman system matrix:\n");
16cbbc2c 423 for ( Int_t i=0; i<fNMeasurementParams; i++ )
55d5da9e 424 {
425 for ( Int_t j=0; j<fgkNSystemParams; j++ )
044eb03e 426 {
55d5da9e 427 printf("% -2.2f ", (*fPH)(i,j) );
428 }//for i
044eb03e 429 printf("\n");
55d5da9e 430 }//for j
431 printf("\n");
432 return;
043badeb 433}
434
044eb03e 435//______________________________________________________________________________
3ebfbf52 436Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
043badeb 437{
3ebfbf52 438 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
fffc83e6 439 fNTracks++; //count added input sets
440
441 //INPUT OUTLIER REJECTION
442 if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
443
3ebfbf52 444 fPTrackParamArr1[fTrackInBuffer] = *exparam1;
445 fPTrackParamArr2[fTrackInBuffer] = *exparam2;
446
447 fTrackInBuffer++;
043badeb 448
3ebfbf52 449 if (fTrackInBuffer == fgkNTracksPerMeasurement)
450 {
451 fTrackInBuffer = 0;
452 return kTRUE;
453 }
454 return kFALSE;
043badeb 455}
456
044eb03e 457//______________________________________________________________________________
3ebfbf52 458Bool_t AliRelAlignerKalman::Update()
043badeb 459{
3ebfbf52 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();
55d5da9e 469 if (!PrepareMeasurement()) return kFALSE;
470 if (!PrepareSystemMatrix()) return kFALSE;
a80268df 471 if (!PreparePrediction()) return kFALSE;
3ebfbf52 472 return UpdateEstimateKalman();
043badeb 473}
474
044eb03e 475//______________________________________________________________________________
043badeb 476void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
477{
55d5da9e 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;
043badeb 495}
496
044eb03e 497//______________________________________________________________________________
498Bool_t AliRelAlignerKalman::PrepareMeasurement()
043badeb 499{
55d5da9e 500 //Calculate the residuals and their covariance matrix
3ebfbf52 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];
16cbbc2c 511 if (!fYZOnly)
512 {
513 (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
514 (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
515 }
3ebfbf52 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];
3ebfbf52 522 (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
523 (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
3ebfbf52 524 (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
525 (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
3ebfbf52 526 (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
527 (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
16cbbc2c 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 }
3ebfbf52 555
3ebfbf52 556 }
557 //if (fApplyCovarianceCorrection)
558 // *fPMeasurementCov += *fPMeasurementCovCorr;
55d5da9e 559 return kTRUE;
043badeb 560}
561
044eb03e 562//______________________________________________________________________________
563Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
043badeb 564{
55d5da9e 565 //Calculate the system matrix for the Kalman filter
566 //approximate the system using as reference the track in the first volume
567
16cbbc2c 568 TVectorD z1( fNMeasurementParams );
569 TVectorD z2( fNMeasurementParams );
55d5da9e 570 TVectorD x1( fgkNSystemParams );
571 TVectorD x2( fgkNSystemParams );
55d5da9e 572 //get the derivatives
573 for ( Int_t i=0; i<fgkNSystemParams; i++ )
574 {
575 x1 = *fPX;
576 x2 = *fPX;
3ebfbf52 577 x1(i) = x1(i) - fDelta[i]/(2.0);
578 x2(i) = x2(i) + fDelta[i]/(2.0);
55d5da9e 579 if (!PredictMeasurement( z1, x1 )) return kFALSE;
580 if (!PredictMeasurement( z2, x2 )) return kFALSE;
16cbbc2c 581 for (Int_t j=0; j<fNMeasurementParams; j++ )
3ebfbf52 582 {
583 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
584 }
55d5da9e 585 }
586 return kTRUE;
043badeb 587}
588
a80268df 589//______________________________________________________________________________
590Bool_t AliRelAlignerKalman::PreparePrediction()
591{
592 //Prepare the prediction of the measurement using state vector
593 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
594}
595
044eb03e 596//______________________________________________________________________________
043badeb 597Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
598{
55d5da9e 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
55d5da9e 605 if (fCorrectionMode)
606 {
3ebfbf52 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];
16cbbc2c 619 if (!fYZOnly)
620 {
621 pred(x+2) = oldparam[2] - newparam[2];
622 pred(x+3) = oldparam[3] - newparam[3];
623 }
3ebfbf52 624 return kTRUE;
625 }
55d5da9e 626 }
627 else
628 {
3ebfbf52 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];
16cbbc2c 641 if (!fYZOnly)
642 {
643 pred(x+2) = newparam[2] - oldparam[2];
644 pred(x+3) = newparam[3] - oldparam[3];
645 }
3ebfbf52 646 return kTRUE;
647 }
55d5da9e 648 }
649 return kFALSE;
043badeb 650}
651
044eb03e 652//______________________________________________________________________________
043badeb 653Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
654{
55d5da9e 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);
3ebfbf52 674
55d5da9e 675 //shit happens so protect yourself!
3ebfbf52 676// if (fNumericalParanoia)
677// {
678// TDecompLU lu(hpht);
679// if (lu.Condition() > 1e12) return kFALSE;
680// lu.Invert(hpht);
681// }
682// else
683// {
da6f4665 684 Double_t det=0.0;
a8f68d07 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....
3ebfbf52 687// }
688 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
689
690 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
55d5da9e 691
692 // update the state and its covariance matrix
693 TVectorD xupdate(fgkNSystemParams);
a80268df 694 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
55d5da9e 695
696 //SIMPLE OUTLIER REJECTION
697 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
698 {
699 fNOutliers++;
fffc83e6 700 //printf("AliRelAlignerKalman: outlier\n");
55d5da9e 701 return kFALSE;
702 }
703
55d5da9e 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
3ebfbf52 708 if (!IsPositiveDefinite(ikhp)) return kFALSE;
709
710 (*fPX) += xupdate;
711 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
712
55d5da9e 713 fNUpdates++;
3ebfbf52 714
55d5da9e 715 return kTRUE;
043badeb 716}
717
044eb03e 718//______________________________________________________________________________
719Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
720{
55d5da9e 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 {
3ebfbf52 726 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
55d5da9e 727 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
728 }
729 return is;
044eb03e 730}
731
fffc83e6 732//______________________________________________________________________________
733Bool_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
3ebfbf52 762//______________________________________________________________________________
763Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
764{
18f65700 765 //check for positive definiteness
766
3ebfbf52 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
044eb03e 778//______________________________________________________________________________
55d5da9e 779void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
043badeb 780{
55d5da9e 781 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
043badeb 782
55d5da9e 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++)
043badeb 787 {
55d5da9e 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;
043badeb 792 }
55d5da9e 793 }
3ebfbf52 794 matsym.MakeValid();
55d5da9e 795 return;
043badeb 796}
797
044eb03e 798//______________________________________________________________________________
043badeb 799void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
800{
55d5da9e 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;
043badeb 807}
808
044eb03e 809//______________________________________________________________________________
043badeb 810void AliRelAlignerKalman::PrintCorrelationMatrix()
811{
55d5da9e 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++ )
043badeb 817 {
3ebfbf52 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) ) );
55d5da9e 821 }//for j
043badeb 822 printf("\n");
55d5da9e 823 }//for i
824 printf("\n");
825 return;
043badeb 826}
827
044eb03e 828//______________________________________________________________________________
c3b5bfc1 829Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
043badeb 830{
55d5da9e 831 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
043badeb 832
55d5da9e 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);
3ebfbf52 835 fMagField = pEvent->GetMagneticField();
55d5da9e 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)
043badeb 863 {
3ebfbf52 864 //std::cout<<"no track!"<<std::endl;
55d5da9e 865 continue;
043badeb 866 }
55d5da9e 867 if (fCuts)
c3b5bfc1 868 {
fffc83e6 869 if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
043badeb 870 }
55d5da9e 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
043badeb 883 {
55d5da9e 884 tracksArrITS[nITStracks] = nGoodTracks;
885 nITStracks++;
886 nGoodTracks++;
887 continue;
043badeb 888 }
889
55d5da9e 890 //check if track is TPC
891 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
892 !(nClsTPC<fMinPointsVol2) ) //enough points
043badeb 893 {
55d5da9e 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;
044eb03e 899 }
55d5da9e 900 }//for itrack - selection fo tracks
044eb03e 901
55d5da9e 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++)
044eb03e 923 {
55d5da9e 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
043badeb 931 {
55d5da9e 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]);
043badeb 940 }
55d5da9e 941 }
043badeb 942 }
55d5da9e 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 }
043badeb 966
55d5da9e 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++)
c3b5bfc1 971 {
55d5da9e 972 candidateTPCtracksArr[i] = tracksArrTPC[i];
c3b5bfc1 973 }
55d5da9e 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++)
c3b5bfc1 987 {
55d5da9e 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());
3ebfbf52 994 partpc.PropagateTo(parits->GetX(),fMagField);
55d5da9e 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 }
c3b5bfc1 1013 }
55d5da9e 1014 if (matchedITStracks) nMatchedITStracks++;
1015 }
044eb03e 1016
55d5da9e 1017 if (nMatchedITStracks==0) //no match in ITS
1018 {
c3b5bfc1 1019 delete [] phiArr;
1020 delete [] thetaArr;
1021 delete [] minDifferenceArr;
044eb03e 1022 delete [] goodtracksArr;
c3b5bfc1 1023 delete [] matchedITStracksArr;
55d5da9e 1024 delete [] candidateTPCtracksArr;
c3b5bfc1 1025 delete [] matchedTPCtracksArr;
55d5da9e 1026 delete [] tracksArrITS;
1027 delete [] tracksArrTPC;
044eb03e 1028 delete [] nPointsArr;
55d5da9e 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;
043badeb 1080}
55d5da9e 1081
044eb03e 1082//______________________________________________________________________________
3ebfbf52 1083Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
043badeb 1084{
55d5da9e 1085 //implements the system model -
1086 //applies correction for misalignment and calibration to track
3ebfbf52 1087 //track needs to be already propagated to the global reference plane
55d5da9e 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
3ebfbf52 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
55d5da9e 1103 //Rotation
1104 TMatrixD rotmat(3,3);
1105 RotMat( rotmat, misal );
3ebfbf52 1106 Point = rotmat.T() * Point;
55d5da9e 1107 Dir = rotmat * Dir;
1108
1109 //TPC vdrift and T0 corrections
3ebfbf52 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
55d5da9e 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 }
3ebfbf52 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
55d5da9e 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;
55d5da9e 1164 //insert everything back into track
1165 const Double_t* pcovtmp = tr->GetCovariance();
3ebfbf52 1166 p[4] = tr->GetSigned1Pt(); //copy the momentum
55d5da9e 1167 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1168 tr->Set(x,alpha,p,pcov);
55d5da9e 1169 return kTRUE;
3ebfbf52 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);
043badeb 1185}
1186
044eb03e 1187//______________________________________________________________________________
3ebfbf52 1188Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
043badeb 1189{
55d5da9e 1190 //implements the system model -
1191 //applies misalignment and miscalibration to reference track
3ebfbf52 1192 //trackparams have to be at the global reference plane
55d5da9e 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
3ebfbf52 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
55d5da9e 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];
3ebfbf52 1249 if (dir[0]==0.0) return kFALSE;
55d5da9e 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]);
3ebfbf52 1254 if (pt==0.0) return kFALSE;
55d5da9e 1255 p[2] = dir[1]/pt;
1256 p[3] = dir[2]/pt;
55d5da9e 1257 //insert everything back into track
1258 const Double_t* pcovtmp = tr->GetCovariance();
3ebfbf52 1259 p[4] = tr->GetSigned1Pt(); //copy the momentum
55d5da9e 1260 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1261 tr->Set(x,alpha,p,pcov);
55d5da9e 1262 return kTRUE;
3ebfbf52 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);
043badeb 1282}
1283
044eb03e 1284//______________________________________________________________________________
55d5da9e 1285void AliRelAlignerKalman::Reset()
043badeb 1286{
3ebfbf52 1287 //full reset to defaults
55d5da9e 1288 fPX->Zero();
1289 (*fPX)(6)=1.;
1290 ResetCovariance();
3ebfbf52 1291
1292 //initialize the differentials per parameter
fffc83e6 1293 for (Int_t i=0;i<4;i++)
1294 {
1295 delete [] (fResArrSigma2Median[i]);
1296 }
1297 fRejectOutliersSigma2Median=kFALSE;
3ebfbf52 1298
1299 fNMatchedCosmics=0;
1300 fNMatchedTPCtracklets=0;
1301 fNUpdates=0;
1302 fNOutliers=0;
1303 fNTracks=0;
1304 fNProcessedEvents=0;
329b5744 1305 fRunNumber=0;
1306 fTimeStamp=0;
043badeb 1307}
1308
044eb03e 1309//______________________________________________________________________________
c3b5bfc1 1310void AliRelAlignerKalman::ResetCovariance( const Double_t number )
043badeb 1311{
55d5da9e 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 {
3ebfbf52 1316 for (Int_t z=0;z<6;z++)
c3b5bfc1 1317 {
3ebfbf52 1318 for (Int_t zz=0;zz<6;zz++)
55d5da9e 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;
c3b5bfc1 1325 }
55d5da9e 1326 }
1327 else
1328 {
1329 //Resets the covariance of the fit to a default value
1330 fPXcov->Zero();
3ebfbf52 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)
55d5da9e 1337 }
3ebfbf52 1338 ResetTPCparamsCovariance(number);
55d5da9e 1339}
1340
1341//______________________________________________________________________________
1342void 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 {
3ebfbf52 1350 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1351 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
a80268df 1352 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
55d5da9e 1353 }
1354 else
1355 {
3ebfbf52 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);
55d5da9e 1359 }
1360
1361 //set crossterms to zero
1362 for (Int_t i=0;i<fgkNSystemParams;i++)
1363 {
3ebfbf52 1364 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
c3b5bfc1 1365 {
55d5da9e 1366 if (i==j) continue; //don't touch diagonals
1367 (*fPXcov)(i,j) = 0.;
1368 (*fPXcov)(j,i) = 0.;
c3b5bfc1 1369 }
55d5da9e 1370 }
043badeb 1371}
55d5da9e 1372
329b5744 1373//______________________________________________________________________________
a80268df 1374Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1375{
1376 //Merge two aligners
1377
1378 if (!al) return kFALSE;
329b5744 1379 if (al==this) return kTRUE;
329b5744 1380 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
a80268df 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
52d1b7b0 1410 if (!success)
1411 {
1412 fNMergesFailed++;
fffc83e6 1413 //printf("AliRelAlignerKalman::Merge failed\n");
52d1b7b0 1414 return kFALSE; //no point in merging stats if merge not succesful
1415 }
a80268df 1416 fNProcessedEvents += al->fNProcessedEvents;
1417 fNUpdates += al->fNUpdates;
1418 fNOutliers += al->fNOutliers;
fffc83e6 1419 fNOutliersSigma2Median += al->fNOutliersSigma2Median;
a80268df 1420 fNTracks += al->fNTracks;
1421 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1422 fNMatchedCosmics += al->fNMatchedCosmics;
fffc83e6 1423 if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1424 else fNMerges += al->fNMerges;
52d1b7b0 1425 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
a80268df 1426
1427 return success;
1428}
1429
52d1b7b0 1430//______________________________________________________________________________
1431Long64_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
329b5744 1446//______________________________________________________________________________
1447Int_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
fffc83e6 1457//________________________________________________________________________
1458Int_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
2c16d467 1551 delete [] matchedArr;
fffc83e6 1552 return iMatched;
1553}
1554
1555//________________________________________________________________________
1556void 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}