/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* Origin: Marian Ivanov Marian.Ivanov@cern.ch Not final version yet - fitting using histograms to be implemented Expected to be more robust. Model: TPC Kalman filter implementation for TPC dEdx alignment: The dEdx equalization for 3 types of the pad: Short (0.75 cm) Medium (1 cm) and Long(1.5 cm) Model of correction ci: corrected raw correction dEdxti = dEdxri*ci = dEdxri*p01*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2) // Matching - update using 2 tracklets ::UpdatedEdxPair dEdxti-dEdxtj=0 = dEdxri*ci-dEdxrj*cj // Matching - normalization of the signal ::UpdatedEdx possible only for identified particles dEdxti = dN/dx * sqrt(1+kY^2+kZ^2) */ #include "TMath.h" #include "TTreeStream.h" #include "TGraph.h" #include "TGraphErrors.h" #include "TVectorD.h" #include "AliTPCROC.h" #include "AliTPCkalmandEdx.h" AliTPCkalmandEdx::AliTPCkalmandEdx(): TNamed(), fState(0), fCovariance(0), fMat1(0), fNpad(3), // number of pad types fNpar(5), // number of parameters fNelem(3*5), // number of elements fSampleSize(0), fInit(0) { // // Default constructor // } AliTPCkalmandEdx::AliTPCkalmandEdx(const char* name, const char* title, Int_t sampleSize): TNamed(name,title), fState(0), fCovariance(0), fMat1(0), fNpad(3), // number of pad types fNpar(5), // number of parameters fNelem(3*5), // number of elements fSampleSize(sampleSize), fInit(0) { // // Default constructor // Init(); } AliTPCkalmandEdx::AliTPCkalmandEdx(const AliTPCkalmandEdx & kalman): TNamed(kalman), fState(0), fCovariance(0), fMat1(0), fNpad(kalman.fNpad), // number of pad types fNpar(kalman.fNpar), // number of parameters fNelem(kalman.fNelem), // number of elements fSampleSize(kalman.fSampleSize) { // // copy constructor // fState = new TMatrixD(*(kalman.fState)); fCovariance = new TMatrixD(*(kalman.fCovariance)); fMat1 = new TMatrixD(*(kalman.fMat1)); } void AliTPCkalmandEdx::Init(){ // // Default parameters setting // fState = new TMatrixD(fNelem,1); fCovariance = new TMatrixD(fNelem,fNelem); fMat1 = new TMatrixD(fNelem,fNelem); fInit=0; for (Int_t i=0;i<3; i++) { fSample[i].ResizeTo(fSampleSize); fSampleStat[i].ResizeTo(2); fCounter[i]=0; } TMatrixD &vecXk=*fState; TMatrixD &covXk=*fCovariance; TMatrixD &mat1=*fMat1; // // for (Int_t i=0;i4*fSampleStat[ip0][1]) return; // // Int_t nmeas = 1; TMatrixD vecXk=*fState; // X vector TMatrixD covXk=*fCovariance; // X covariance TMatrixD &mat1=*fMat1; // update covariance matrix Double_t length0 = TMath::Sqrt(1+kY0*kY0+kZ0*kZ0); // Double_t corr0 = vecXk(ip0*fNpar,0)* (1+vecXk(ip0*fNpar+1,0)*kY0+ vecXk(ip0*fNpar+2,0)*kZ0+vecXk(ip0*fNpar+3,0)*dR0+vecXk(ip0*fNpar+4,0)/length0); // TMatrixD vecZk(nmeas,nmeas); TMatrixD measR(nmeas,nmeas); TMatrixD matHk(nmeas,fNelem); // vector to mesurement TMatrixD vecYk(nmeas,nmeas); // Innovation or measurement residual TMatrixD matHkT(fNelem,nmeas); // helper matrix Hk transpose TMatrixD matSk(nmeas,nmeas); // Innovation (or residual) covariance TMatrixD matKk(fNelem,nmeas); // Optimal Kalman gain TMatrixD covXk2(fNelem,fNelem); // helper matrix TMatrixD covXk3(fNelem,fNelem); // helper matrix vecZk(0,0) =dedxRef/dedx0; measR(0,0) =s0*s0*dedxRef*dedxRef/(dedx0*dedx0); // // reset matHk for (Int_t iel=0;iel=fSampleSize) return; fSample[ip0][fCounter[ip0]]=dedxRef/dedx0; fCounter[ip0]++; if (fCounter[ip0]==fSampleSize){ fSampleStat[ip0].ResizeTo(2); fSampleStat[ip0][0] = TMath::Median(fSampleSize,fSample[ip0].GetMatrixArray()); fSampleStat[ip0][1] = TMath::RMS(fSampleSize,fSample[ip0].GetMatrixArray()); fInit++; (*fState)(ip0*fNpar,0) = fSampleStat[ip0][0]; (*fCovariance)(ip0*fNpar,ip0*fNpar) = fSampleStat[ip0][1]*fSampleStat[ip0][1]/fSampleSize; Double_t med2 = fSampleStat[ip0][0]*fSampleStat[ip0][0]; // // (*fCovariance)(ip0*fNpar+1,ip0*fNpar+1)=0.01*med2; (*fCovariance)(ip0*fNpar+2,ip0*fNpar+2)=0.01*med2; (*fCovariance)(ip0*fNpar+3,ip0*fNpar+3)=0.01*med2; (*fCovariance)(ip0*fNpar+4,ip0*fNpar+4)=0.01*med2; } }