Setting of aliases to rawReader done only once. Base decision on cosmic or calib...
[u/mrichter/AliRoot.git] / STEER / 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(),
a6585a16 90 fPTrackParam1(new AliExternalTrackParam()),
91 fPTrackParam2(new AliExternalTrackParam()),
3ebfbf52 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),
a80268df 123 fTimeStamp(0),
329b5744 124 fRunNumber(0),
52d1b7b0 125 fNMerges(0),
126 fNMergesFailed(0),
55d5da9e 127 fTPCvd(2.64),
3ebfbf52 128 fTPCZLengthA(2.4972500e02),
129 fTPCZLengthC(2.4969799e02)
043badeb 130{
55d5da9e 131 //Default constructor
fffc83e6 132 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
133 for (Int_t i=0; i<4;i++){fResArrSigma2Median[i]=NULL;}
55d5da9e 134 Reset();
55d5da9e 135}
136
137//______________________________________________________________________________
138AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
3ebfbf52 139 TObject(static_cast<TObject>(a)),
a6585a16 140 fPTrackParam1(new AliExternalTrackParam()),
141 fPTrackParam2(new AliExternalTrackParam()),
3ebfbf52 142 fMagField(a.fMagField),
16cbbc2c 143 fNMeasurementParams(a.fNMeasurementParams),
55d5da9e 144 fPX(new TVectorD( *a.fPX )),
145 fPXcov(new TMatrixDSym( *a.fPXcov )),
52d1b7b0 146 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
55d5da9e 147 fQ(a.fQ),
52d1b7b0 148 fPMeasurement(new TVectorD( fNMeasurementParams )),
149 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
150 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
55d5da9e 151 fOutRejSigmas(a.fOutRejSigmas),
fffc83e6 152 fOutRejSigma2Median(a.fOutRejSigma2Median),
16cbbc2c 153 fYZOnly(a.fYZOnly),
3ebfbf52 154 fNumericalParanoia(a.fNumericalParanoia),
55d5da9e 155 fRejectOutliers(a.fRejectOutliers),
fffc83e6 156 fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
55d5da9e 157 fRequireMatchInTPC(a.fRequireMatchInTPC),
55d5da9e 158 fCuts(a.fCuts),
159 fMinPointsVol1(a.fMinPointsVol1),
160 fMinPointsVol2(a.fMinPointsVol2),
fffc83e6 161 fMinPt(a.fMinPt),
162 fMaxPt(a.fMaxPt),
55d5da9e 163 fMaxMatchingAngle(a.fMaxMatchingAngle),
164 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
165 fCorrectionMode(a.fCorrectionMode),
166 fNTracks(a.fNTracks),
167 fNUpdates(a.fNUpdates),
168 fNOutliers(a.fNOutliers),
fffc83e6 169 fNOutliersSigma2Median(a.fNOutliersSigma2Median),
55d5da9e 170 fNMatchedCosmics(a.fNMatchedCosmics),
171 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
172 fNProcessedEvents(a.fNProcessedEvents),
a80268df 173 fTimeStamp(a.fTimeStamp),
329b5744 174 fRunNumber(a.fRunNumber),
52d1b7b0 175 fNMerges(a.fNMerges),
176 fNMergesFailed(a.fNMergesFailed),
55d5da9e 177 fTPCvd(a.fTPCvd),
178 fTPCZLengthA(a.fTPCZLengthA),
179 fTPCZLengthC(a.fTPCZLengthC)
180{
181 //copy constructor
182 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
fffc83e6 183
184 //copy contents of the residuals array for sigma2median scheme
185 for (Int_t i=0;i<4;i++)
186 {
187 if ((a.fResArrSigma2Median)[i])
188 {
189 fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
190 memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
191 fgkNtracksSigma2Median*sizeof(Double_t));
192 }
193 else
194 fResArrSigma2Median[i] = NULL;
195 }
55d5da9e 196}
197
198//______________________________________________________________________________
199AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
200{
201 //assignment operator
e99fb5c9 202 if(&a == this) return *this;
203 TObject::operator=(a);
204
205 fMagField=a.fMagField;
52d1b7b0 206 fNMeasurementParams=a.fNMeasurementParams;
55d5da9e 207 *fPX = *a.fPX;
208 *fPXcov = *a.fPXcov;
55d5da9e 209 fQ=a.fQ;
55d5da9e 210 fOutRejSigmas=a.fOutRejSigmas;
fffc83e6 211 fOutRejSigma2Median=a.fOutRejSigma2Median;
55d5da9e 212 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
16cbbc2c 213 fYZOnly=a.fYZOnly;
3ebfbf52 214 fNumericalParanoia=a.fNumericalParanoia;
55d5da9e 215 fRejectOutliers=a.fRejectOutliers;
fffc83e6 216 fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
55d5da9e 217 fRequireMatchInTPC=a.fRequireMatchInTPC;
55d5da9e 218 fCuts=a.fCuts;
219 fMinPointsVol1=a.fMinPointsVol1;
220 fMinPointsVol2=a.fMinPointsVol2;
fffc83e6 221 fMinPt=a.fMinPt;
222 fMaxPt=a.fMaxPt;
55d5da9e 223 fMaxMatchingAngle=a.fMaxMatchingAngle;
224 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
225 fCorrectionMode=a.fCorrectionMode;
226 fNTracks=a.fNTracks;
227 fNUpdates=a.fNUpdates;
228 fNOutliers=a.fNOutliers;
fffc83e6 229 fNOutliersSigma2Median=a.fNOutliersSigma2Median;
55d5da9e 230 fNMatchedCosmics=a.fNMatchedCosmics;
231 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
232 fNProcessedEvents=a.fNProcessedEvents;
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
a6585a16 260 delete fPTrackParam1;
261 delete fPTrackParam2;
fffc83e6 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
044eb03e 324//______________________________________________________________________________
18f65700 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
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
044eb03e 381//______________________________________________________________________________
16cbbc2c 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
396//______________________________________________________________________________
044eb03e 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
a6585a16 444 *fPTrackParam1 = *exparam1;
445 *fPTrackParam2 = *exparam2;
3ebfbf52 446
a6585a16 447 return kTRUE;
043badeb 448}
449
044eb03e 450//______________________________________________________________________________
3ebfbf52 451Bool_t AliRelAlignerKalman::Update()
043badeb 452{
3ebfbf52 453 //perform the update
454
455 //if (fCalibrationMode) return UpdateCalibration();
456 //if (fFillHistograms)
457 //{
458 // if (!UpdateEstimateKalman()) return kFALSE;
459 // return UpdateCalibration(); //Update histograms only when update ok.
460 //}
461 //else return UpdateEstimateKalman();
55d5da9e 462 if (!PrepareMeasurement()) return kFALSE;
463 if (!PrepareSystemMatrix()) return kFALSE;
a80268df 464 if (!PreparePrediction()) return kFALSE;
3ebfbf52 465 return UpdateEstimateKalman();
043badeb 466}
467
044eb03e 468//______________________________________________________________________________
043badeb 469void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
470{
55d5da9e 471 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
472 Double_t sinpsi = TMath::Sin(angles(0));
473 Double_t sintheta = TMath::Sin(angles(1));
474 Double_t sinphi = TMath::Sin(angles(2));
475 Double_t cospsi = TMath::Cos(angles(0));
476 Double_t costheta = TMath::Cos(angles(1));
477 Double_t cosphi = TMath::Cos(angles(2));
478
479 R(0,0) = costheta*cosphi;
480 R(0,1) = -costheta*sinphi;
481 R(0,2) = sintheta;
482 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
483 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
484 R(1,2) = -costheta*sinpsi;
485 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
486 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
487 R(2,2) = costheta*cospsi;
043badeb 488}
489
044eb03e 490//______________________________________________________________________________
491Bool_t AliRelAlignerKalman::PrepareMeasurement()
043badeb 492{
55d5da9e 493 //Calculate the residuals and their covariance matrix
3ebfbf52 494
a6585a16 495 const Double_t* pararr1 = fPTrackParam1->GetParameter();
496 const Double_t* pararr2 = fPTrackParam2->GetParameter();
3ebfbf52 497
498 //Take the track parameters and calculate the input to the Kalman filter
a6585a16 499 (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
500 (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
16cbbc2c 501 if (!fYZOnly)
502 {
a6585a16 503 (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
504 (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
16cbbc2c 505 }
3ebfbf52 506
507 //the covariance
a6585a16 508 const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
509 const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
510 (*fPMeasurementCov)(0,0)=parcovarr1[0];
511 (*fPMeasurementCov)(0,1)=parcovarr1[1];
512 (*fPMeasurementCov)(1,0)=parcovarr1[1];
513 (*fPMeasurementCov)(1,1)=parcovarr1[2];
514 (*fPMeasurementCov)(0,0)+=parcovarr2[0];
515 (*fPMeasurementCov)(0,1)+=parcovarr2[1];
516 (*fPMeasurementCov)(1,0)+=parcovarr2[1];
517 (*fPMeasurementCov)(1,1)+=parcovarr2[2];
16cbbc2c 518 if (!fYZOnly)
519 {
a6585a16 520 (*fPMeasurementCov)(0,2)=parcovarr1[3];
521 (*fPMeasurementCov)(0,3)=parcovarr1[6];
522 (*fPMeasurementCov)(1,2)=parcovarr1[4];
523 (*fPMeasurementCov)(1,3)=parcovarr1[7];
524 (*fPMeasurementCov)(2,0)=parcovarr1[3];
525 (*fPMeasurementCov)(2,1)=parcovarr1[4];
526 (*fPMeasurementCov)(2,2)=parcovarr1[5];
527 (*fPMeasurementCov)(2,3)=parcovarr1[8];
528 (*fPMeasurementCov)(3,0)=parcovarr1[6];
529 (*fPMeasurementCov)(3,1)=parcovarr1[7];
530 (*fPMeasurementCov)(3,2)=parcovarr1[8];
531 (*fPMeasurementCov)(3,3)=parcovarr1[9];
532 (*fPMeasurementCov)(0,2)+=parcovarr2[3];
533 (*fPMeasurementCov)(0,3)+=parcovarr2[6];
534 (*fPMeasurementCov)(1,2)+=parcovarr2[4];
535 (*fPMeasurementCov)(1,3)+=parcovarr2[7];
536 (*fPMeasurementCov)(2,0)+=parcovarr2[3];
537 (*fPMeasurementCov)(2,1)+=parcovarr2[4];
538 (*fPMeasurementCov)(2,2)+=parcovarr2[5];
539 (*fPMeasurementCov)(2,3)+=parcovarr2[8];
540 (*fPMeasurementCov)(3,0)+=parcovarr2[6];
541 (*fPMeasurementCov)(3,1)+=parcovarr2[7];
542 (*fPMeasurementCov)(3,2)+=parcovarr2[8];
543 (*fPMeasurementCov)(3,3)+=parcovarr2[9];
16cbbc2c 544 }
3ebfbf52 545
3ebfbf52 546 //if (fApplyCovarianceCorrection)
547 // *fPMeasurementCov += *fPMeasurementCovCorr;
55d5da9e 548 return kTRUE;
043badeb 549}
550
044eb03e 551//______________________________________________________________________________
552Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
043badeb 553{
55d5da9e 554 //Calculate the system matrix for the Kalman filter
555 //approximate the system using as reference the track in the first volume
556
16cbbc2c 557 TVectorD z1( fNMeasurementParams );
558 TVectorD z2( fNMeasurementParams );
55d5da9e 559 TVectorD x1( fgkNSystemParams );
560 TVectorD x2( fgkNSystemParams );
55d5da9e 561 //get the derivatives
562 for ( Int_t i=0; i<fgkNSystemParams; i++ )
563 {
564 x1 = *fPX;
565 x2 = *fPX;
3ebfbf52 566 x1(i) = x1(i) - fDelta[i]/(2.0);
567 x2(i) = x2(i) + fDelta[i]/(2.0);
55d5da9e 568 if (!PredictMeasurement( z1, x1 )) return kFALSE;
569 if (!PredictMeasurement( z2, x2 )) return kFALSE;
16cbbc2c 570 for (Int_t j=0; j<fNMeasurementParams; j++ )
3ebfbf52 571 {
572 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
573 }
55d5da9e 574 }
575 return kTRUE;
043badeb 576}
577
044eb03e 578//______________________________________________________________________________
a80268df 579Bool_t AliRelAlignerKalman::PreparePrediction()
580{
581 //Prepare the prediction of the measurement using state vector
582 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
583}
584
585//______________________________________________________________________________
043badeb 586Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
587{
55d5da9e 588 // Implements a system model for the Kalman fit
589 // pred is [dy,dz,dsinphi,dtgl]
590 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
591 // note: the measurement is in a local frame, so the prediction also has to be
592 // note: state is the misalignment in global reference system
593
55d5da9e 594 if (fCorrectionMode)
595 {
a6585a16 596 AliExternalTrackParam track(*fPTrackParam2); //make a copy track
3ebfbf52 597 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
598
a6585a16 599 const Double_t* oldparam = fPTrackParam2->GetParameter();
3ebfbf52 600 const Double_t* newparam = track.GetParameter();
601
3ebfbf52 602 //calculate the predicted residual
a6585a16 603 pred(0) = oldparam[0] - newparam[0];
604 pred(1) = oldparam[1] - newparam[1];
16cbbc2c 605 if (!fYZOnly)
606 {
a6585a16 607 pred(2) = oldparam[2] - newparam[2];
608 pred(3) = oldparam[3] - newparam[3];
16cbbc2c 609 }
3ebfbf52 610 return kTRUE;
55d5da9e 611 }
612 else
613 {
b2d08a09 614 AliExternalTrackParam track(*fPTrackParam1); //make a copy track
3ebfbf52 615 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
616
a6585a16 617 const Double_t* oldparam = fPTrackParam1->GetParameter();
3ebfbf52 618 const Double_t* newparam = track.GetParameter();
619
3ebfbf52 620 //calculate the predicted residual
a6585a16 621 pred(0) = newparam[0] - oldparam[0];
622 pred(1) = newparam[1] - oldparam[1];
16cbbc2c 623 if (!fYZOnly)
624 {
a6585a16 625 pred(2) = newparam[2] - oldparam[2];
626 pred(3) = newparam[3] - oldparam[3];
16cbbc2c 627 }
3ebfbf52 628 return kTRUE;
629 }
55d5da9e 630 return kFALSE;
043badeb 631}
632
044eb03e 633//______________________________________________________________________________
043badeb 634Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
635{
55d5da9e 636 //Kalman estimation of noisy constants: in the model A=1
637 //The arguments are (following the usual convention):
638 // fPX - the state vector (parameters)
639 // fPXcov - the state covariance matrix (parameter errors)
640 // fPMeasurement - measurement vector
641 // fPMeasurementCov - measurement covariance matrix
642 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
643
644 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
645
646 //predict the state
647 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
648
649 // update prediction with measurement
650 // calculate Kalman gain
651 // K = PH/(HPH+fPMeasurementCov)
652 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
653 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
654 hpht += (*fPMeasurementCov);
3ebfbf52 655
55d5da9e 656 //shit happens so protect yourself!
3ebfbf52 657// if (fNumericalParanoia)
658// {
659// TDecompLU lu(hpht);
660// if (lu.Condition() > 1e12) return kFALSE;
661// lu.Invert(hpht);
662// }
663// else
664// {
da6f4665 665 Double_t det=0.0;
a8f68d07 666 hpht.Invert(&det); //since the matrix is small...
667 if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
3ebfbf52 668// }
669 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
670
671 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
55d5da9e 672
673 // update the state and its covariance matrix
674 TVectorD xupdate(fgkNSystemParams);
a80268df 675 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
55d5da9e 676
677 //SIMPLE OUTLIER REJECTION
678 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
679 {
680 fNOutliers++;
fffc83e6 681 //printf("AliRelAlignerKalman: outlier\n");
55d5da9e 682 return kFALSE;
683 }
684
55d5da9e 685 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
686 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
687 ikh = identity - kh;
688 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
3ebfbf52 689 if (!IsPositiveDefinite(ikhp)) return kFALSE;
690
691 (*fPX) += xupdate;
692 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
693
55d5da9e 694 fNUpdates++;
3ebfbf52 695
55d5da9e 696 return kTRUE;
043badeb 697}
698
044eb03e 699//______________________________________________________________________________
700Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
701{
55d5da9e 702 //check whether an update is an outlier given the covariance matrix of the fit
703
704 Bool_t is=kFALSE;
705 for (Int_t i=0;i<fgkNSystemParams;i++)
706 {
3ebfbf52 707 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
55d5da9e 708 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
709 }
710 return is;
044eb03e 711}
712
713//______________________________________________________________________________
fffc83e6 714Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS,
715 const AliExternalTrackParam* pTPC )
716{
717 //check if the input residuals are not too far off their median
718 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
719 TVectorD vecDeltaN(5);
720 Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
721 vecDeltaN[4]=0;
722 for (Int_t i=0;i<4;i++){
723 vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
724 (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
725 }
726 Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
727 for (Int_t i=0;i<fNMeasurementParams;i++){ //in point2trackmode just take the first 2 params (zy)
728 vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
729 vecRMS[i] = TMath::RMS(entries,fResArrSigma2Median[i]);
730 vecDeltaN[i] = 0;
731 if (vecRMS[i]>0.){
732 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
733 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
734 }
735 }
736 Bool_t outlier=kFALSE;
737 if (fNTracks<3) outlier=kTRUE; //median and RMS still to be defined
738 if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
739 if (outlier) fNOutliersSigma2Median++;
740 return outlier;
741}
742
743//______________________________________________________________________________
3ebfbf52 744Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
745{
18f65700 746 //check for positive definiteness
747
3ebfbf52 748 for (Int_t i=0; i<mat.GetNcols(); i++)
749 {
750 if (mat(i,i)<=0.) return kFALSE;
751 }
752
753 if (!fNumericalParanoia) return kTRUE;
754
755 TDecompLU lu(mat);
756 return (lu.Decompose());
757}
758
759//______________________________________________________________________________
55d5da9e 760void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
043badeb 761{
55d5da9e 762 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
043badeb 763
55d5da9e 764 for (Int_t i=0; i<mat.GetNcols(); i++)
765 {
766 matsym(i,i) = mat(i,i); //copy diagonal
767 for (Int_t j=i+1; j<mat.GetNcols(); j++)
043badeb 768 {
55d5da9e 769 //copy the rest
770 Double_t average = (mat(i,j)+mat(j,i))/2.;
771 matsym(i,j)=average;
772 matsym(j,i)=average;
043badeb 773 }
55d5da9e 774 }
3ebfbf52 775 matsym.MakeValid();
55d5da9e 776 return;
043badeb 777}
778
044eb03e 779//______________________________________________________________________________
043badeb 780void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
781{
55d5da9e 782 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
783 //b = R*a
784 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
785 angles(1) = TMath::ASin( rotmat(0,2) );
786 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
787 return;
043badeb 788}
789
044eb03e 790//______________________________________________________________________________
043badeb 791void AliRelAlignerKalman::PrintCorrelationMatrix()
792{
55d5da9e 793 //Print the correlation matrix for the fitted parameters
794 printf("Correlation matrix for system parameters:\n");
795 for ( Int_t i=0; i<fgkNSystemParams; i++ )
796 {
797 for ( Int_t j=0; j<i+1; j++ )
043badeb 798 {
3ebfbf52 799 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
800 else
801 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
55d5da9e 802 }//for j
043badeb 803 printf("\n");
55d5da9e 804 }//for i
805 printf("\n");
806 return;
043badeb 807}
808
044eb03e 809//______________________________________________________________________________
c3b5bfc1 810Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
043badeb 811{
55d5da9e 812 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
043badeb 813
55d5da9e 814 //Sanity cuts on tracks + check which tracks are ITS which are TPC
815 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
3ebfbf52 816 fMagField = pEvent->GetMagneticField();
55d5da9e 817 if (ntracks<2)
818 {
819 //printf("TrackFinder: less than 2 tracks!\n");
820 return kFALSE;
821 }
822 Float_t* phiArr = new Float_t[ntracks];
823 Float_t* thetaArr = new Float_t[ntracks];
824 Int_t* goodtracksArr = new Int_t[ntracks];
825 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
826 Int_t* matchedITStracksArr = new Int_t[ntracks];
827 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
828 Int_t* tracksArrITS = new Int_t[ntracks];
829 Int_t* tracksArrTPC = new Int_t[ntracks];
830 Int_t* nPointsArr = new Int_t[ntracks];
831 Int_t nITStracks = 0;
832 Int_t nTPCtracks = 0;
833 Int_t nGoodTracks = 0;
834 Int_t nCandidateTPCtracks = 0;
835 Int_t nMatchedITStracks = 0;
836 AliESDtrack* pTrack = NULL;
837 Bool_t foundMatchTPC = kFALSE;
838
839 //select and clasify tracks
840 for (Int_t itrack=0; itrack < ntracks; itrack++)
841 {
842 pTrack = pEvent->GetTrack(itrack);
843 if (!pTrack)
043badeb 844 {
3ebfbf52 845 //std::cout<<"no track!"<<std::endl;
55d5da9e 846 continue;
043badeb 847 }
55d5da9e 848 if (fCuts)
c3b5bfc1 849 {
fffc83e6 850 if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
043badeb 851 }
55d5da9e 852 goodtracksArr[nGoodTracks]=itrack;
853 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
854 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
855 phiArr[nGoodTracks]=phi;
856 thetaArr[nGoodTracks]=theta;
857
858 //check if track is ITS
859 Int_t nClsITS = pTrack->GetNcls(0);
860 Int_t nClsTPC = pTrack->GetNcls(1);
861 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
862 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
863 !(nClsITS<fMinPointsVol1) ) //enough points
043badeb 864 {
55d5da9e 865 tracksArrITS[nITStracks] = nGoodTracks;
866 nITStracks++;
867 nGoodTracks++;
868 continue;
043badeb 869 }
870
55d5da9e 871 //check if track is TPC
872 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
873 !(nClsTPC<fMinPointsVol2) ) //enough points
043badeb 874 {
55d5da9e 875 tracksArrTPC[nTPCtracks] = nGoodTracks;
876 nTPCtracks++;
877 nGoodTracks++;
878 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
879 continue;
044eb03e 880 }
55d5da9e 881 }//for itrack - selection fo tracks
044eb03e 882
55d5da9e 883 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
884
885 if ( nITStracks < 1 || nTPCtracks < 1 )
886 {
887 delete [] phiArr;
888 delete [] thetaArr;
889 delete [] goodtracksArr;
890 delete [] matchedITStracksArr;
891 delete [] candidateTPCtracksArr;
892 delete [] matchedTPCtracksArr;
893 delete [] tracksArrITS;
894 delete [] tracksArrTPC;
895 delete [] nPointsArr;
896 return kFALSE;
897 }
898
899 //find matching in TPC
900 if (nTPCtracks>1) //if there is something to be matched, try and match it
901 {
902 Float_t min = 10000000.;
903 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
044eb03e 904 {
55d5da9e 905 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
906 {
907 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
908 if (deltatheta > fMaxMatchingAngle) continue;
909 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
910 if (deltaphi > fMaxMatchingAngle) continue;
911 if (deltatheta+deltaphi<min) //only the best matching pair
043badeb 912 {
55d5da9e 913 min=deltatheta+deltaphi;
914 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
915 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
916 nCandidateTPCtracks = 2;
917 foundMatchTPC = kTRUE;
918 //printf("TrackFinder: Matching TPC tracks candidates:\n");
919 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
920 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
043badeb 921 }
55d5da9e 922 }
043badeb 923 }
55d5da9e 924 }//if nTPCtracks>1
925 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
926 {
927 candidateTPCtracksArr[0] = tracksArrTPC[0];
928 nCandidateTPCtracks = 1;
929 foundMatchTPC = kFALSE;
930 }
931 if (foundMatchTPC) fNMatchedTPCtracklets++;
932 //if no match but the requirement is set return kFALSE
933 if (fRequireMatchInTPC && !foundMatchTPC)
934 {
935 delete [] phiArr;
936 delete [] thetaArr;
937 delete [] goodtracksArr;
938 delete [] candidateTPCtracksArr;
939 delete [] matchedITStracksArr;
940 delete [] matchedTPCtracksArr;
941 delete [] tracksArrITS;
942 delete [] tracksArrTPC;
943 delete [] nPointsArr;
944 //printf("TrackFinder: no match in TPC && required\n");
945 return kFALSE;
946 }
043badeb 947
55d5da9e 948 //if no match and more than one track take all TPC tracks
949 if (!fRequireMatchInTPC && !foundMatchTPC)
950 {
951 for (Int_t i=0;i<nTPCtracks;i++)
c3b5bfc1 952 {
55d5da9e 953 candidateTPCtracksArr[i] = tracksArrTPC[i];
c3b5bfc1 954 }
55d5da9e 955 nCandidateTPCtracks = nTPCtracks;
956 }
957 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
958
959 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
960
961 //find ITS matches for good TPC tracks
962 Bool_t matchedITStracks=kFALSE;
963 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
964 {
965 minDifferenceArr[nMatchedITStracks] = 10000000.;
966 matchedITStracks=kFALSE;
967 for (Int_t iits=0; iits<nITStracks; iits++)
c3b5bfc1 968 {
55d5da9e 969 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
970 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
971 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
972 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
973 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
974 partpc.Rotate(parits->GetAlpha());
3ebfbf52 975 partpc.PropagateTo(parits->GetX(),fMagField);
55d5da9e 976 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
977 if (dtgl > fMaxMatchingAngle) continue;
978 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
979 if (dsnp > fMaxMatchingAngle) continue;
980 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
981 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
982 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
983 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
984 {
985 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
986 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
987 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
988 //printf("TrackFinder: Matching ITS to TPC:\n");
989 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
990 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
991 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
992 matchedITStracks=kTRUE;;
993 }
c3b5bfc1 994 }
55d5da9e 995 if (matchedITStracks) nMatchedITStracks++;
996 }
044eb03e 997
55d5da9e 998 if (nMatchedITStracks==0) //no match in ITS
999 {
c3b5bfc1 1000 delete [] phiArr;
1001 delete [] thetaArr;
1002 delete [] minDifferenceArr;
044eb03e 1003 delete [] goodtracksArr;
c3b5bfc1 1004 delete [] matchedITStracksArr;
55d5da9e 1005 delete [] candidateTPCtracksArr;
c3b5bfc1 1006 delete [] matchedTPCtracksArr;
55d5da9e 1007 delete [] tracksArrITS;
1008 delete [] tracksArrTPC;
044eb03e 1009 delete [] nPointsArr;
55d5da9e 1010 //printf("TrackFinder: No match in ITS\n");
1011 return kFALSE;
1012 }
1013
1014 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1015 //we found a cosmic
1016 fNMatchedCosmics++;
1017
1018 //Now we may have ended up with more matches than we want in the case there was
1019 //no TPC match and there were many TPC tracks
1020 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1021 //so take only the best pair.
1022 //same applies when there are more matches than ITS tracks - means that one ITS track
1023 //matches more TPC tracks.
1024 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1025 {
1026 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1027 matchedITStracksArr[0] = matchedITStracksArr[imin];
1028 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1029 nMatchedITStracks = 1;
1030 //printf("TrackFinder: too many matches - take only the best one\n");
1031 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1032 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1033 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1034 }
1035
1036 ///////////////////////////////////////////////////////////////////////////
1037 outITSindexTArr.Set(nMatchedITStracks);
1038 outTPCindexTArr.Set(nMatchedITStracks);
1039 for (Int_t i=0;i<nMatchedITStracks;i++)
1040 {
1041 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1042 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1043 //printf("TrackFinder: Fill the output\n");
1044 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1045 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1046 }
1047 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1048 ///////////////////////////////////////////////////////////////////////////
1049
1050 delete [] phiArr;
1051 delete [] thetaArr;
1052 delete [] minDifferenceArr;
1053 delete [] goodtracksArr;
1054 delete [] candidateTPCtracksArr;
1055 delete [] matchedITStracksArr;
1056 delete [] matchedTPCtracksArr;
1057 delete [] tracksArrITS;
1058 delete [] tracksArrTPC;
1059 delete [] nPointsArr;
1060 return kTRUE;
043badeb 1061}
55d5da9e 1062
044eb03e 1063//______________________________________________________________________________
3ebfbf52 1064Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
043badeb 1065{
55d5da9e 1066 //implements the system model -
1067 //applies correction for misalignment and calibration to track
3ebfbf52 1068 //track needs to be already propagated to the global reference plane
55d5da9e 1069
1070 Double_t x = tr->GetX();
1071 Double_t alpha = tr->GetAlpha();
1072 Double_t point[3],dir[3];
1073 tr->GetXYZ(point);
1074 tr->GetDirection(dir);
1075 TVector3 Point(point);
1076 TVector3 Dir(dir);
1077
1078 //Apply corrections to track
1079
1080 //Shift
3ebfbf52 1081 Point(0) -= misal(3); //add shift in x
1082 Point(1) -= misal(4); //add shift in y
1083 Point(2) -= misal(5); //add shift in z
55d5da9e 1084 //Rotation
1085 TMatrixD rotmat(3,3);
1086 RotMat( rotmat, misal );
3ebfbf52 1087 Point = rotmat.T() * Point;
55d5da9e 1088 Dir = rotmat * Dir;
1089
1090 //TPC vdrift and T0 corrections
3ebfbf52 1091 TVector3 Point2; //second point of the track
1092 Point2 = Point + Dir;
1093 Double_t vdCorr = 1./misal(6);
1094 Double_t t0 = misal(7);
1095 Double_t vdY = 0.0;
1096 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1097
1098 //my model
55d5da9e 1099 if (Point(2)>0)
1100 {
1101 //A-Side
1102 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1103 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1104 }
1105 else
1106 {
1107 //C-side
1108 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1109 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1110 }
3ebfbf52 1111
1112 //Stefan's model
1113 //if (Point(2)>0)
1114 //{
1115 // //A-Side
1116 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1117 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1118 //}
1119 //else
1120 //{
1121 // //C-side
1122 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1123 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1124 //}
1125
55d5da9e 1126 Dir = Point2-Point;
1127 Dir=Dir.Unit(); //keep unit length
1128
1129 //Turn back to local system
1130 Dir.GetXYZ(dir);
1131 Point.GetXYZ(point);
1132 tr->Global2LocalPosition(point,alpha);
1133 tr->Global2LocalPosition(dir,alpha);
1134
1135 //Calculate new intersection point with ref plane
1136 Double_t p[5],pcov[15];
1137 if (dir[0]==0.0) return kFALSE;
1138 Double_t s=(x-point[0])/dir[0];
1139 p[0] = point[1]+s*dir[1];
1140 p[1] = point[2]+s*dir[2];
1141 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1142 if (pt==0.0) return kFALSE;
1143 p[2] = dir[1]/pt;
1144 p[3] = dir[2]/pt;
55d5da9e 1145 //insert everything back into track
1146 const Double_t* pcovtmp = tr->GetCovariance();
3ebfbf52 1147 p[4] = tr->GetSigned1Pt(); //copy the momentum
55d5da9e 1148 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1149 tr->Set(x,alpha,p,pcov);
55d5da9e 1150 return kTRUE;
3ebfbf52 1151
1152 ////put params back into track and propagate to ref
1153 //Double_t p[5],pcov[15];
1154 //p[0] = point[1];
1155 //p[1] = point[2];
1156 //Double_t xnew = point[0];
1157 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1158 //if (pt==0.0) return kFALSE;
1159 //p[2] = dir[1]/pt;
1160 //p[3] = dir[2]/pt;
1161 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1162 //const Double_t* pcovtmp = tr->GetCovariance();
1163 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1164 //tr->Set(xnew,alpha,p,pcov);
1165 //return tr->PropagateTo(x,fMagField);
043badeb 1166}
1167
044eb03e 1168//______________________________________________________________________________
3ebfbf52 1169Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
043badeb 1170{
55d5da9e 1171 //implements the system model -
1172 //applies misalignment and miscalibration to reference track
3ebfbf52 1173 //trackparams have to be at the global reference plane
55d5da9e 1174
1175 Double_t x = tr->GetX();
1176 Double_t alpha = tr->GetAlpha();
1177 Double_t point[3],dir[3];
1178 tr->GetXYZ(point);
1179 tr->GetDirection(dir);
1180 TVector3 Point(point);
1181 TVector3 Dir(dir);
1182
1183 //Apply misalignment to track
1184
1185 //TPC vdrift and T0 corrections
3ebfbf52 1186 TVector3 Point2; //second point of the track
1187 Point2 = Point + Dir;
1188 Double_t vdCorr = 1./misal(6);
1189 Double_t t0 = misal(7);
1190 Double_t vdY = 0.0;
1191 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1192
55d5da9e 1193 if (Point(2)>0)
1194 {
1195 //A-Side
1196 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1197 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1198 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1199 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1200 }
1201 else
1202 {
1203 //C-side
1204 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1205 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1206 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1207 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1208 }
1209 Dir = Point2-Point;
1210 Dir=Dir.Unit(); //keep unit length
1211
1212 //Rotation
1213 TMatrixD rotmat(3,3);
1214 RotMat( rotmat, misal );
1215 Point = rotmat * Point;
1216 Dir = rotmat * Dir;
1217 //Shift
1218 Point(0) += misal(3); //add shift in x
1219 Point(1) += misal(4); //add shift in y
1220 Point(2) += misal(5); //add shift in z
1221
1222 //Turn back to local system
1223 Dir.GetXYZ(dir);
1224 Point.GetXYZ(point);
1225 tr->Global2LocalPosition(point,alpha);
1226 tr->Global2LocalPosition(dir,alpha);
1227
1228 //Calculate new intersection point with ref plane
1229 Double_t p[5],pcov[15];
3ebfbf52 1230 if (dir[0]==0.0) return kFALSE;
55d5da9e 1231 Double_t s=(x-point[0])/dir[0];
1232 p[0] = point[1]+s*dir[1];
1233 p[1] = point[2]+s*dir[2];
1234 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
3ebfbf52 1235 if (pt==0.0) return kFALSE;
55d5da9e 1236 p[2] = dir[1]/pt;
1237 p[3] = dir[2]/pt;
55d5da9e 1238 //insert everything back into track
1239 const Double_t* pcovtmp = tr->GetCovariance();
3ebfbf52 1240 p[4] = tr->GetSigned1Pt(); //copy the momentum
55d5da9e 1241 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1242 tr->Set(x,alpha,p,pcov);
55d5da9e 1243 return kTRUE;
3ebfbf52 1244
1245 ////put params back into track and propagate to ref
1246 //Double_t p[5];
1247 //Double_t pcov[15];
1248 //p[0] = point[1];
1249 //p[1] = point[2];
1250 //Double_t xnew = point[0];
1251 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1252 //if (pt==0.0) return kFALSE;
1253 //p[2] = dir[1]/pt;
1254 //p[3] = dir[2]/pt;
1255 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1256 //const Double_t* pcovtmp = tr->GetCovariance();
1257 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1258 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1259 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1260 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1261 //tr->Set(xnew,alpha,p,pcov);
1262 //return tr->PropagateTo(x,fMagField);
043badeb 1263}
1264
044eb03e 1265//______________________________________________________________________________
55d5da9e 1266void AliRelAlignerKalman::Reset()
043badeb 1267{
3ebfbf52 1268 //full reset to defaults
55d5da9e 1269 fPX->Zero();
1270 (*fPX)(6)=1.;
1271 ResetCovariance();
3ebfbf52 1272
1273 //initialize the differentials per parameter
fffc83e6 1274 for (Int_t i=0;i<4;i++)
1275 {
1276 delete [] (fResArrSigma2Median[i]);
1277 }
1278 fRejectOutliersSigma2Median=kFALSE;
3ebfbf52 1279
1280 fNMatchedCosmics=0;
1281 fNMatchedTPCtracklets=0;
1282 fNUpdates=0;
1283 fNOutliers=0;
1284 fNTracks=0;
1285 fNProcessedEvents=0;
329b5744 1286 fRunNumber=0;
1287 fTimeStamp=0;
043badeb 1288}
1289
044eb03e 1290//______________________________________________________________________________
c3b5bfc1 1291void AliRelAlignerKalman::ResetCovariance( const Double_t number )
043badeb 1292{
55d5da9e 1293 //Resets the covariance to the default if arg=0 or resets the off diagonals
1294 //to zero and releases the diagonals by factor arg.
1295 if (number!=0.)
1296 {
3ebfbf52 1297 for (Int_t z=0;z<6;z++)
c3b5bfc1 1298 {
3ebfbf52 1299 for (Int_t zz=0;zz<6;zz++)
55d5da9e 1300 {
1301 if (zz==z) continue; //don't touch diagonals
1302 (*fPXcov)(zz,z) = 0.;
1303 (*fPXcov)(z,zz) = 0.;
1304 }
1305 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
c3b5bfc1 1306 }
55d5da9e 1307 }
1308 else
1309 {
1310 //Resets the covariance of the fit to a default value
1311 fPXcov->Zero();
3ebfbf52 1312 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1313 (*fPXcov)(1,1) = .08*.08; //theta (rad
1314 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1315 (*fPXcov)(3,3) = .3*.3; //x (cm)
1316 (*fPXcov)(4,4) = .3*.3; //y (cm)
1317 (*fPXcov)(5,5) = .3*.3; //z (cm)
55d5da9e 1318 }
3ebfbf52 1319 ResetTPCparamsCovariance(number);
55d5da9e 1320}
1321
1322//______________________________________________________________________________
1323void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1324{
1325 //Resets the covariance to the default if arg=0 or resets the off diagonals
1326 //to zero and releases the diagonals by factor arg.
1327
1328 //release diagonals
1329 if (number==0.)
1330 {
3ebfbf52 1331 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1332 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
a80268df 1333 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
55d5da9e 1334 }
1335 else
1336 {
3ebfbf52 1337 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1338 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1339 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
55d5da9e 1340 }
1341
1342 //set crossterms to zero
1343 for (Int_t i=0;i<fgkNSystemParams;i++)
1344 {
3ebfbf52 1345 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
c3b5bfc1 1346 {
55d5da9e 1347 if (i==j) continue; //don't touch diagonals
1348 (*fPXcov)(i,j) = 0.;
1349 (*fPXcov)(j,i) = 0.;
c3b5bfc1 1350 }
55d5da9e 1351 }
043badeb 1352}
55d5da9e 1353
329b5744 1354//______________________________________________________________________________
a80268df 1355Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1356{
1357 //Merge two aligners
1358
1359 if (!al) return kFALSE;
329b5744 1360 if (al==this) return kTRUE;
329b5744 1361 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
a80268df 1362
1363 //store the pointers to current stuff
1364 TVectorD* pmes = fPMeasurement;
1365 TMatrixDSym* pmescov = fPMeasurementCov;
1366 TVectorD* pmespred = fPMeasurementPrediction;
1367 TMatrixD* ph = fPH;
1368
1369 //make a unity system matrix
1370 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1371 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1372
1373 //mesurement is the state of the new aligner
1374 fPMeasurement = al->fPX;
1375 fPMeasurementCov = al->fPXcov;
1376
1377 //the mesurement prediction is the state
1378 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1379
1380 //do the merging
1381 Bool_t success = UpdateEstimateKalman();
1382
1383 //restore pointers to old stuff
1384 fPMeasurement = pmes;
1385 fPMeasurementCov = pmescov;
1386 fPMeasurementPrediction = pmespred;
1387 delete fPH;
1388 fPH = ph;
1389
1390 //merge stats
52d1b7b0 1391 if (!success)
1392 {
1393 fNMergesFailed++;
fffc83e6 1394 //printf("AliRelAlignerKalman::Merge failed\n");
52d1b7b0 1395 return kFALSE; //no point in merging stats if merge not succesful
1396 }
a80268df 1397 fNProcessedEvents += al->fNProcessedEvents;
1398 fNUpdates += al->fNUpdates;
1399 fNOutliers += al->fNOutliers;
fffc83e6 1400 fNOutliersSigma2Median += al->fNOutliersSigma2Median;
a80268df 1401 fNTracks += al->fNTracks;
1402 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1403 fNMatchedCosmics += al->fNMatchedCosmics;
fffc83e6 1404 if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1405 else fNMerges += al->fNMerges;
52d1b7b0 1406 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
a80268df 1407
1408 return success;
1409}
1410
3ebfbf52 1411//______________________________________________________________________________
52d1b7b0 1412Long64_t AliRelAlignerKalman::Merge( TCollection* list )
1413{
1414 //merge all aligners in the collection
1415 Long64_t numberOfMerges=0;
1416 AliRelAlignerKalman* alignerFromList;
1417 if (!list) return 0;
1418 TIter next(list);
1419 while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
1420 {
1421 if (alignerFromList == this) continue;
1422 if (Merge(alignerFromList)) numberOfMerges++;
1423 }
1424 return numberOfMerges;
1425}
1426
1427//______________________________________________________________________________
329b5744 1428Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1429{
1430 if (this == obj) return 0;
1431 const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1432 if (!aobj) return 0;
1433 if (fTimeStamp < aobj->fTimeStamp) return -1;
1434 else if (fTimeStamp > aobj->fTimeStamp) return 1;
1435 else return 0;
1436}
1437
fffc83e6 1438//________________________________________________________________________
1439Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
1440{
1441 //find matching tracks and return tobjarrays with the params
1442
1443 Int_t ntracks = pESD->GetNumberOfTracks();
1444 Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
1445 for (Int_t i=0;i<ntracks;i++)
1446 {
1447 matchedArr[i]=0;
1448 }
1449
1450 Int_t iMatched=-1;
1451 for (Int_t i=0; i<ntracks; i++)
1452 {
1453 //get track1 and friend
1454 AliESDtrack* track1 = pESD->GetTrack(i);
1455 if (!track1) continue;
1456
1457 if (track1->GetNcls(0) < fMinPointsVol1) continue;
1458
1459 if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
1460 (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
1461
1462 const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
1463 if (!constfriendtrack1) continue;
1464 AliESDfriendTrack friendtrack1(*constfriendtrack1);
1465
1466 if (!friendtrack1.GetITSOut()) continue;
1467 AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
1468
1469 Double_t bestd = 1000.; //best distance
1470 Bool_t newi = kTRUE; //whether we start with a new i
1471 for (Int_t j=0; j<ntracks; j++)
1472 {
1473 if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
1474 //get track2 and friend
1475 AliESDtrack* track2 = pESD->GetTrack(j);
1476 if (!track2) continue;
1477 if (track1==track2) continue;
1478 //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
1479 // (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
1480
1481 if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
1482 if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
1483 if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
1484 if (track2->GetTgl() > 1.) continue; //acceptance
1485 //cut crossing tracks
1486 if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
1487 if (track2->GetInnerParam()->GetX()>90) continue;
1488 if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
1489
1490
1491 if (!track2->GetInnerParam()) continue;
1492 AliExternalTrackParam params2(*(track2->GetInnerParam()));
1493
1494 //bring to same reference plane
1495 if (!params2.Rotate(params1.GetAlpha())) continue;
1496 if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
1497
1498 //pt cut
1499 if (params2.Pt() < fMinPt) continue;
1500
1501 const Double32_t* p1 = params1.GetParameter();
1502 const Double32_t* p2 = params2.GetParameter();
1503
1504 //hard cuts
1505 Double_t dy = TMath::Abs(p2[0]-p1[0]);
1506 Double_t dz = TMath::Abs(p2[1]-p1[1]);
1507 Double_t dphi = TMath::Abs(p2[2]-p1[2]);
1508 Double_t dlam = TMath::Abs(p2[3]-p1[3]);
1509 if (dy > 2.0) continue;
1510 if (dz > 10.0) continue;
1511 if (dphi > 0.1 ) continue;
1512 if (dlam > 0.1 ) continue;
1513
1514 //best match only
1515 Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
1516 if ( d >= bestd) continue;
1517 bestd = d;
1518 matchedArr[j]=i; //j-th track matches i-th (ITS) track
1519 if (newi) iMatched++; newi=kFALSE; //increment at most once per i
1520 if (arrITS[iMatched] && arrTPC[iMatched])
1521 {
1522 *(arrITS[iMatched]) = params1;
1523 *(arrTPC[iMatched]) = params2;
1524 }
1525 else
1526 {
1527 arrITS[iMatched] = new AliExternalTrackParam(params1);
1528 arrTPC[iMatched] = new AliExternalTrackParam(params2);
1529 }//else
1530 }//for j
1531 }//for i
2c16d467 1532 delete [] matchedArr;
fffc83e6 1533 return iMatched;
1534}
1535
1536//________________________________________________________________________
1537void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
1538{
1539 //Sets up or destroys the memory hungry array to hold the statistics
1540 //for data rejection with median
1541 if (set)
1542 {
1543 for (Int_t i=0;i<4;i++)
1544 {
1545 if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] =
1546 new Double_t[fgkNtracksSigma2Median];
1547 }
1548 fRejectOutliersSigma2Median = kTRUE;
1549 }//else
1550 else
1551 {
1552 // it probably doesn't make sense to delete the arrays, they are not streamed
1553 //if (fRejectOutliersSigma2Median)
1554 //for (Int_t i=0;i<4;i++)
1555 //{
1556 // delete [] (fResArrSigma2Median[i]);
1557 //}
1558 fRejectOutliersSigma2Median = kFALSE;
1559 }//if
1560}
d01710c6 1561
1562
1563Int_t AliRelAlignerKalman::CheckCovariance(){
1564 //
1565 // check covariance matrix
1566 // Return values:
1567 // -1 - everything OK
1568 // >=0 - index of "corrupted" element in the covariance matrix
1569 TMatrixDSym& mat =(*fPXcov);
1570 for (Int_t irow=0; irow<9; irow++){
1571 if (mat(irow,irow)<=0.) return irow*9;
1572 for (Int_t icol=0; icol<irow; icol++){
1573 Double_t corel = mat(irow,irow)*mat(icol,icol);
1574 if (corel<=0.) return irow*9+icol;
1575 corel=mat(irow,icol)/TMath::Sqrt(corel);
1576 if (TMath::Abs(corel)>=1.) return irow*9+icol;
1577 }
1578 }
1579 return -1;
1580}