Update of the TPC-ITS alignment code (Mikolaj)

diff --git a/STEER/AliRelAlignerKalman.cxx b/STEER/AliRelAlignerKalman.cxx
index b9f538ec03b57f276c0c0490ada242c1fbaf0f60..20c626a1cbfbf08feda52297cadfc9eadab13df3 100644 (file)
--- a/STEER/AliRelAlignerKalman.cxx
+++ b/STEER/AliRelAlignerKalman.cxx
@@ -103,6 +103,7 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fQ(1.e-15),
     fPMeasurement(new TVectorD( fgkNMeasurementParams )),
     fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
+    fPMeasurementPrediction(new TVectorD( fgkNMeasurementParams )),
     fOutRejSigmas(1.),
     fDelta(new Double_t[fgkNSystemParams]),
     fNumericalParanoia(kTRUE),
@@ -123,6 +124,7 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fNMatchedTPCtracklets(0),
     fNProcessedEvents(0),
     fTrackInBuffer(0),
+    fTimeStamp(0),
     fTPCvd(2.64),
     fTPCZLengthA(2.4972500e02),
     fTPCZLengthC(2.4969799e02)
@@ -147,6 +149,7 @@ AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
     fQ(a.fQ),
     fPMeasurement(new TVectorD( *a.fPMeasurement )),
     fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
+    fPMeasurementPrediction(new TVectorD( *a.fPMeasurement )),
     fOutRejSigmas(a.fOutRejSigmas),
     fDelta(new Double_t[fgkNSystemParams]),
     fNumericalParanoia(a.fNumericalParanoia),
@@ -167,6 +170,7 @@ AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
     fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
     fNProcessedEvents(a.fNProcessedEvents),
     fTrackInBuffer(a.fTrackInBuffer),
+    fTimeStamp(a.fTimeStamp),
     fTPCvd(a.fTPCvd),
     fTPCZLengthA(a.fTPCZLengthA),
     fTPCZLengthC(a.fTPCZLengthC)
@@ -223,6 +227,7 @@ AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a
   fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
   fNProcessedEvents=a.fNProcessedEvents;
   fTrackInBuffer=a.fTrackInBuffer;
+  fTimeStamp=a.fTimeStamp;
   //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
   //fCalibrationMode=a.fCalibrationMode;
   //fFillHistograms=a.fFillHistograms;
@@ -288,33 +293,28 @@ Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
       success = ( AddESDtrack( track ) || success );
     }
   }
+  if (success) fTimeStamp = pEvent->GetTimeStamp();
   return success;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
 {
-  ////Adds a full track, returns true if results in a new estimate
-
-  //const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
-  //if (!pconstparamsITS) return kFALSE;
-  //const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
-  //if (!pconstparamsTPC) return kFALSE;
-  //
-  ////TPC part
-  //AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
-  //paramsTPC.Rotate(pconstparamsITS->GetAlpha());
-  //paramsTPC.PropagateTo(pconstparamsITS->GetX(), pconstparamsITS->GetAlpha());
-
-  //if (!SetTrackParams( pconstparamsITS, &paramsTPC )) return kFALSE;
-  //
-  ////do some accounting and update
-  //return (Update());
-
-  const AliESDtrack* p;
-  p=pTrack; //shuts up the compiler
+  //Adds a full track, returns true if results in a new estimate
+  //  gets the inner TPC parameters from AliESDTrack::GetInnerParam()
+  //  gets the outer ITS parameters from AliESDTrack::GetOuterParam()
+
+  const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
+  if (!pconstparamsITS) return kFALSE;
+  const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
+  if (!pconstparamsTPC) return kFALSE;
   
-  return kFALSE;
+  //TPC part
+  AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
+  paramsTPC.Rotate(pconstparamsITS->GetAlpha());
+  paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
+
+  return (AddTrackParams(pconstparamsITS, &paramsTPC));
 }
 
 //______________________________________________________________________________
@@ -359,7 +359,7 @@ Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
     if (!pconstparams2) continue;
     params2 = *pconstparams2; //make copy
     params2.Rotate(params1.GetAlpha());
-    params2.PropagateTo( params1.GetX(), params1.GetAlpha() );
+    params2.PropagateTo( params1.GetX(), fMagField );
 
     if (!SetTrackParams( &params1, &params2 )) continue;
     
@@ -369,6 +369,7 @@ Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
     else
       continue;
   }
+  if (success) fTimeStamp=pEvent->GetTimeStamp();
   return success;
 }
 
@@ -453,6 +454,7 @@ Bool_t AliRelAlignerKalman::Update()
   //else return UpdateEstimateKalman();
   if (!PrepareMeasurement()) return kFALSE;
   if (!PrepareSystemMatrix()) return kFALSE;
+  if (!PreparePrediction()) return kFALSE;
   return UpdateEstimateKalman();
 }
 
@@ -565,6 +567,13 @@ Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
   return kTRUE;
 }
 
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PreparePrediction()
+{
+  //Prepare the prediction of the measurement using state vector
+  return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
+}
+
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
 {
@@ -657,9 +666,7 @@ Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 
   // update the state and its covariance matrix
   TVectorD xupdate(fgkNSystemParams);
-  TVectorD hx(fgkNMeasurementParams);
-  PredictMeasurement( hx, (*fPX) );
-  xupdate = k*((*fPMeasurement)-hx);
+  xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
 
   //SIMPLE OUTLIER REJECTION
   if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
@@ -1280,7 +1287,7 @@ void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
   {
     if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
     if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
-    if (fgkNSystemParams>8) (*fPXcov)(8,8) = 1.*1.;
+    if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
   }
   else
   {
@@ -1301,6 +1308,52 @@ void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
   }
 }
 
+Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
+{
+  //Merge two aligners
+  
+  if (!al) return kFALSE;
+  if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE;
+  
+  //store the pointers to current stuff
+  TVectorD* pmes = fPMeasurement;
+  TMatrixDSym* pmescov = fPMeasurementCov;
+  TVectorD* pmespred = fPMeasurementPrediction;
+  TMatrixD* ph = fPH;
+
+  //make a unity system matrix
+  TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
+  fPH = new TMatrixD(TMatrixD::kUnit, tmp);
+
+  //mesurement is the state of the new aligner
+  fPMeasurement = al->fPX;
+  fPMeasurementCov = al->fPXcov;
+
+  //the mesurement prediction is the state
+  fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
+  
+  //do the merging
+  Bool_t success = UpdateEstimateKalman();
+  
+  //restore pointers to old stuff
+  fPMeasurement = pmes;
+  fPMeasurementCov = pmescov;
+  fPMeasurementPrediction = pmespred;
+  delete fPH;
+  fPH = ph;
+
+  //merge stats
+  fNProcessedEvents += al->fNProcessedEvents;
+  fNUpdates += al->fNUpdates;
+  fNOutliers += al->fNOutliers;
+  fNTracks += al->fNTracks;
+  fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
+  fNMatchedCosmics += al->fNMatchedCosmics;
+  if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //a bit arbitrary: take the older one
+
+  return success;
+}
+
 //______________________________________________________________________________
 //void AliRelAlignerKalman::PrintCovarianceCorrection()
 //{

diff --git a/STEER/AliRelAlignerKalman.h b/STEER/AliRelAlignerKalman.h
index cc16cf78cf68098cf3b13ef3ca8be8363b7a5801..95c5a5d2ba1025935eeb08b06bd3c584bccc0e94 100644 (file)
--- a/STEER/AliRelAlignerKalman.h
+++ b/STEER/AliRelAlignerKalman.h
@@ -68,6 +68,7 @@ public:
     TMatrixDSym* GetStateCov() const { return fPXcov; }
     void GetSeed( TVectorD& seed, TMatrixDSym& seedCov ) const { seed = *fPX; seedCov = *fPXcov; }
     void SetSeed( const TVectorD& seed, const TMatrixDSym& seedCov ) {*fPX = seed; *fPXcov = seedCov; }
+    Bool_t Merge( const AliRelAlignerKalman* al );
 
     //Expert methods:
     Bool_t AddESDevent( const AliESDEvent* pEvent );
@@ -111,11 +112,14 @@ public:
     static void RotMat( TMatrixD& R, const TVectorD& angles );
     static void TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat );
     Bool_t IsPositiveDefinite( const TMatrixD& mat ) const;
+    void SetTimeStamp( const UInt_t ts ) { fTimeStamp = ts; }
+    UInt_t GetTimeStamp() const {return fTimeStamp;}
     
 protected:
     Bool_t UpdateEstimateKalman();
     Bool_t PrepareMeasurement();
     Bool_t PrepareSystemMatrix();
+    Bool_t PreparePrediction();
     Bool_t PredictMeasurement( TVectorD& z, const TVectorD& x );
     Bool_t IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix );
 
@@ -129,7 +133,7 @@ private:
     Double_t fLocalX;      //!local x coordinate of reference plane = r(global)
     AliExternalTrackParam* fPTrackParamArr1;   //!local track parameters
     AliExternalTrackParam* fPTrackParamArr2;   //!local track parameters
-    Double_t fMagField; //magnetic field
+    Double_t fMagField; //!magnetic field
 
     //Kalman filter related stuff
     TVectorD* fPX; //System (fit) parameters (phi, theta, psi, x, y, z, driftcorr, driftoffset )
@@ -138,20 +142,21 @@ private:
     Double_t fQ;        //!measure for system noise
     TVectorD* fPMeasurement; //!the measurement vec for Kalman filter (theta,phi,x,z)
     TMatrixDSym* fPMeasurementCov; //!measurement vec cvariance
+    TVectorD* fPMeasurementPrediction; //!prediction of the measurement
     Double_t fOutRejSigmas; //number of sigmas for outlier rejection
     Double_t* fDelta; //array with differentials for calculating derivatives for every parameter(see PrepareSystemMatrix())
 
     //Control
-    Bool_t fNumericalParanoia; //whether to perform additional checks for numerical stability
-    Bool_t fRejectOutliers; //whether to do outlier rejection in the Kalman filter
-    Bool_t fRequireMatchInTPC;  //when looking for a cosmic in event, require that TPC has 2 matching segments
-    Bool_t fCuts;    //track cuts?
-    Int_t fMinPointsVol1;   //mininum number of points in volume 1
-    Int_t fMinPointsVol2;   //mininum number of points in volume 2
-    Double_t fMinMom;       //min momentum of track for track cuts
-    Double_t fMaxMom;       //max momentum of track for track cuts
-    Double_t fMaxMatchingAngle; //cuts
-    Double_t fMaxMatchingDistance; //cuts
+    Bool_t fNumericalParanoia; //!whether to perform additional checks for numerical stability
+    Bool_t fRejectOutliers; //!whether to do outlier rejection in the Kalman filter
+    Bool_t fRequireMatchInTPC;  //!when looking for a cosmic in event, require that TPC has 2 matching segments
+    Bool_t fCuts;    //!track cuts?
+    Int_t fMinPointsVol1;   //!mininum number of points in volume 1
+    Int_t fMinPointsVol2;   //!mininum number of points in volume 2
+    Double_t fMinMom;       //!min momentum of track for track cuts
+    Double_t fMaxMom;       //!max momentum of track for track cuts
+    Double_t fMaxMatchingAngle; //!cuts
+    Double_t fMaxMatchingDistance; //!cuts
     Bool_t fCorrectionMode; //calculate corrective transform for TPC (or monitor actual TPC misal params)
     
     //Counters
@@ -161,7 +166,8 @@ private:
     Int_t fNMatchedCosmics; //number of cosmic events with matching tracklets (good cosmics)
     Int_t fNMatchedTPCtracklets;//number of cosmic events with 2 matching TPC tracklets
     Int_t fNProcessedEvents; //number of processed events
-    Int_t fTrackInBuffer; //number of tracks in buffer
+    Int_t fTrackInBuffer; //!number of tracks in buffer
+    UInt_t fTimeStamp;
 
     //TPC stuff
     Double_t fTPCvd; //TPC drift velocity