TPC/AliTPCkalmandEdx.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17    Origin: Marian Ivanov   Marian.Ivanov@cern.ch
  18    Not final version yet - fitting using histograms to be implemented
  19    Expected to be more robust.
  20
  21
  22    Model:
  23    TPC Kalman filter implementation for TPC dEdx alignment:
  24    The dEdx equalization for 3 types of the pad:
  25    Short (0.75 cm) Medium (1 cm) and Long(1.5 cm)
  26
  27    Model of correction ci:
  28    corrected   raw    correction
  29    dEdxti    = dEdxri*ci = dEdxri*p01*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
  30    //
  31    Matching - update using 2 tracklets  ::UpdatedEdxPair
  32    dEdxti-dEdxtj=0 = dEdxri*ci-dEdxrj*cj
  33    //
  34    Matching - normalization of the signal ::UpdatedEdx
  35    possible only for identified particles
  36    dEdxti = dN/dx * sqrt(1+kY^2+kZ^2)
  37 */
  38
  39 #include "TMath.h"
  40 #include "TTreeStream.h"
  41 #include "TGraph.h"
  42 #include "TGraphErrors.h"
  43 #include "TVectorD.h"
  44 #include "AliTPCROC.h"
  45 #include "AliTPCkalmandEdx.h"
  46
  47
  48 AliTPCkalmandEdx::AliTPCkalmandEdx():
  49   TNamed(),
  50   fState(0),
  51   fCovariance(0),
  52   fMat1(0),
  53   fNpad(3),                  // number of pad types
  54   fNpar(5),                // number of parameters
  55   fNelem(3*5),                // number of elements
  56   fSampleSize(0),
  57   fInit(0)
  58 {
  59   //
  60   // Default constructor
  61   //
  62 }
  63
  64 AliTPCkalmandEdx::AliTPCkalmandEdx(const char* name, const char* title, Int_t sampleSize):
  65   TNamed(name,title),
  66   fState(0),
  67   fCovariance(0),
  68   fMat1(0),
  69   fNpad(3),                  // number of pad types
  70   fNpar(5),                // number of parameters
  71   fNelem(3*5),             // number of elements
  72   fSampleSize(sampleSize),
  73   fInit(0)
  74 {
  75   //
  76   // Default constructor
  77   //
  78   Init();
  79 }
  80
  81 AliTPCkalmandEdx::AliTPCkalmandEdx(const AliTPCkalmandEdx & kalman):
  82   TNamed(kalman),
  83   fState(0),
  84   fCovariance(0),
  85   fMat1(0),
  86   fNpad(kalman.fNpad),                  // number of pad types
  87   fNpar(kalman.fNpar),                // number of parameters
  88   fNelem(kalman.fNelem),                // number of elements
  89   fSampleSize(kalman.fSampleSize),
  90   fInit(kalman.fInit)
  91 {
  92   //
  93   // copy constructor
  94   //
  95   fState      = new TMatrixD(*(kalman.fState));
  96   fCovariance = new TMatrixD(*(kalman.fCovariance));
  97   fMat1       = new TMatrixD(*(kalman.fMat1));
  98 }
  99
 100 void AliTPCkalmandEdx::Init(){
 101   //
 102   // Default parameters setting
 103   //
 104   fState      = new TMatrixD(fNelem,1);
 105   fCovariance = new TMatrixD(fNelem,fNelem);
 106   fMat1       = new TMatrixD(fNelem,fNelem);
 107
 108   fInit=0;
 109   for (Int_t i=0;i<3; i++) {
 110     fSample[i].ResizeTo(fSampleSize);
 111     fSampleStat[i].ResizeTo(2);
 112     fCounter[i]=0;
 113   }
 114
 115   TMatrixD &vecXk=*fState;
 116   TMatrixD &covXk=*fCovariance;
 117   TMatrixD &mat1=*fMat1;
 118   //
 119   //
 120   for (Int_t i=0;i<fNelem;i++)
 121     for(Int_t j=0;j<fNelem;j++) {covXk(i,j)=0;  mat1(i,j)=0;}
 122   //
 123   for (Int_t i=0;i<fNelem;i++)  {covXk(i,i)=1.; mat1(i,i)=1; vecXk(i,0)=0;}
 124   for (Int_t ipad=0;ipad<3;ipad++){
 125     vecXk(ipad*fNpar,0)=1;
 126   }
 127 //   //
 128 //   // set reference ipad=0
 129 //   vecXk(1*fNpar+0,0)=1;
 130 //   vecXk(1*fNpar+1,0)=0;
 131 //   vecXk(1*fNpar+2,0)=0;
 132 //   vecXk(1*fNpar+3,0)=0;
 133 //   vecXk(1*fNpar+4,0)=0;
 134
 135
 136 }
 137
 138 void AliTPCkalmandEdx::UpdatedEdxPair(Int_t ip0, Int_t ip1,
 139                                       Double_t dedx0, Double_t dedx1,
 140                                       Double_t s0, Double_t s1,
 141                                       Double_t kY0, Double_t kY1,
 142                                       Double_t kZ0, Double_t kZ1,
 143                                       Double_t dR0, Double_t dR1,
 144                                       TTreeSRedirector *debug){
 145   //
 146   // update relative measurement
 147   // 0 = dEdxti-dEdxtj  = dEdxri*ci-dEdxrj*cj
 148   //
 149   // Model of correction ci:
 150   // dEdxti = dEdxri*ci = dEdxri*p01*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
 151   if (fInit<3) return;  // not initialized parameters
 152   const Double_t kchi2Cut = 3.*3.;
 153
 154   Int_t nmeas = 1;
 155   TMatrixD vecXk=*fState;           // X vector
 156   TMatrixD covXk=*fCovariance;      // X covariance
 157   TMatrixD &mat1=*fMat1;            // update covariance matrix
 158
 159   Double_t length0 = TMath::Sqrt(1+kY0*kY0+kZ0*kZ0);
 160   Double_t length1 = TMath::Sqrt(1+kY1*kY1+kZ1*kZ1);
 161
 162   Double_t corr0 = vecXk(ip0*fNpar,0)*
 163     (1+vecXk(ip0*fNpar+1,0)*kY0+
 164      vecXk(ip0*fNpar+2,0)*kZ0+vecXk(ip0*fNpar+3,0)*dR0+vecXk(ip0*fNpar+4,0)/length0);
 165   Double_t corr1 = vecXk(ip1*fNpar,0)*
 166     (1+vecXk(ip1*fNpar+1,0)*kY1+
 167      vecXk(ip1*fNpar+2,0)*kZ1+vecXk(ip1*fNpar+3,0)*dR1+vecXk(ip1*fNpar+4,0)/length1);
 168   //
 169   TMatrixD vecZk(nmeas,nmeas);
 170   TMatrixD measR(nmeas,nmeas);
 171   TMatrixD matHk(nmeas,fNelem);     // vector to mesurement
 172   TMatrixD vecYk(nmeas,nmeas);          // Innovation or measurement residual
 173   TMatrixD matHkT(fNelem,nmeas);    // helper matrix Hk transpose
 174   TMatrixD matSk(nmeas,nmeas);      // Innovation (or residual) covariance
 175   TMatrixD matKk(fNelem,nmeas);     // Optimal Kalman gain
 176   TMatrixD covXk2(fNelem,fNelem);   // helper matrix
 177   TMatrixD covXk3(fNelem,fNelem);   // helper matrix
 178
 179   //
 180   vecZk(0,0) =dedx1/dedx0;          // dedx ratio
 181   measR(0,0) =(s0*s0+s1*s1);
 182   //
 183   // reset matHk
 184   for (Int_t iel=0;iel<fNelem;iel++)
 185     for (Int_t ip=0;ip<nmeas;ip++) matHk(ip,iel)=0;
 186
 187   //
 188   matHk(0, ip0*fNpar+0) = corr0/(corr1*vecXk(ip0*fNpar+0,0));
 189   matHk(0, ip0*fNpar+1) = (vecXk(ip0*fNpar+0,0)*kY0)/corr1;
 190   matHk(0, ip0*fNpar+2) = (vecXk(ip0*fNpar+0,0)*kZ0)/corr1;
 191   matHk(0, ip0*fNpar+3) = (vecXk(ip0*fNpar+0,0)*dR0)/corr1;
 192   matHk(0, ip0*fNpar+4) = (vecXk(ip0*fNpar+0,0)/length0)/corr1;
 193   //
 194   //  matHk(0, ip1*fNpar+0) = ;
 195   //matHk(0, ip1*fNpar+1) = ;
 196   //matHk(0, ip1*fNpar+2) = ;
 197   //matHk(0, ip1*fNpar+3) = ;
 198   //matHk(0, ip1*fNpar+4) = ;
 199   //
 200   //
 201   vecYk(0,0) = vecZk(0,0)- corr0/corr1;               // Innovation or measurement residual
 202   matHkT=matHk.T(); matHk.T();
 203   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 204   matSk.Invert();
 205   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 206   vecXk += matKk*vecYk;                    //  updated vector
 207   covXk2= (mat1-(matKk*matHk));
 208   covXk3 =  covXk2*covXk;
 209   covXk = covXk3;
 210   Double_t chi2 = vecYk(0,0)* vecYk(0,0)*matSk(0,0);
 211   //
 212   //
 213   //
 214   if (debug){ // dump debug info
 215     char chname[1000];
 216     sprintf(chname,"updatePair_%s",GetName());
 217     (*debug)<<chname<<
 218       // input parameters
 219       "ip0="<<ip0<<
 220       "ip1="<<ip1<<
 221       "dedx0="<<dedx0<<
 222       "dedx1="<<dedx1<<
 223       "s0="<<s0<<
 224       "s1="<<s1<<
 225       "kY0="<<kY0<<
 226       "kY1="<<kY1<<
 227       "kZ0="<<kZ0<<
 228       "kZ1="<<kZ1<<
 229       "dR0="<<dR0<<
 230       "dR1="<<dR1<<
 231       "chi2="<<chi2<<
 232       "corr0="<<corr0<<
 233       "corr1="<<corr1<<
 234       "length0="<<length0<<
 235       "length1="<<length1<<
 236       //
 237       "vecYk.="<<&vecYk<<
 238       "vecZk.="<<&vecZk<<
 239       "matHk.="<<&matHk<<
 240       "matHkT.="<<&matHkT<<
 241       "matSk.="<<&matSk<<
 242       "matKk.="<<&matKk<<
 243       "covXk2.="<<&covXk2<<
 244       "covXk.="<<&covXk<<
 245       "cov.="<<fCovariance<<
 246       "vecXk.="<<&vecXk<<
 247       "vec.="<<fState<<
 248       "\n";
 249   }
 250   if (chi2<kchi2Cut){
 251     //    (*fCovariance)=covXk;
 252     //(*fState)=vecXk;
 253   }
 254   else{
 255     printf("chi2=%f\n",chi2);
 256   }
 257 }
 258
 259 void AliTPCkalmandEdx::UpdatedEdx(Int_t ip0,
 260                                   Double_t dedx0,
 261                                   Double_t dedxRef,
 262                                   Double_t s0,
 263                                   Double_t kY0,
 264                                   Double_t kZ0,
 265                                   Double_t dR0,
 266                                   TTreeSRedirector *debug){
 267   //
 268   // update relative measurement
 269   // dEdx  = dEdxti = dEdxri*ci
 270   //
 271   // Model of correction ci:
 272   // dEdxti = dEdxri*ci = dEdxri*p0i*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
 273   const Double_t kchi2Cut = 3.*3.;
 274   // removing of "big outliers
 275   if (fInit<3) return AdddEdx(ip0,dedx0,dedxRef);
 276   if (TMath::Abs(dedxRef/dedx0-fSampleStat[ip0][0])>4*fSampleStat[ip0][1]) return;
 277   //
 278   //
 279   Int_t nmeas = 1;
 280   TMatrixD vecXk=*fState;           // X vector
 281   TMatrixD covXk=*fCovariance;      // X covariance
 282   TMatrixD &mat1=*fMat1;            // update covariance matrix
 283
 284   Double_t length0 = TMath::Sqrt(1+kY0*kY0+kZ0*kZ0);
 285   //
 286   Double_t corr0 = vecXk(ip0*fNpar,0)*
 287     (1+vecXk(ip0*fNpar+1,0)*kY0+
 288      vecXk(ip0*fNpar+2,0)*kZ0+vecXk(ip0*fNpar+3,0)*dR0+vecXk(ip0*fNpar+4,0)/length0);
 289
 290
 291   //
 292   TMatrixD vecZk(nmeas,nmeas);
 293   TMatrixD measR(nmeas,nmeas);
 294   TMatrixD matHk(nmeas,fNelem);     // vector to mesurement
 295   TMatrixD vecYk(nmeas,nmeas);          // Innovation or measurement residual
 296   TMatrixD matHkT(fNelem,nmeas);    // helper matrix Hk transpose
 297   TMatrixD matSk(nmeas,nmeas);      // Innovation (or residual) covariance
 298   TMatrixD matKk(fNelem,nmeas);     // Optimal Kalman gain
 299   TMatrixD covXk2(fNelem,fNelem);   // helper matrix
 300   TMatrixD covXk3(fNelem,fNelem);   // helper matrix
 301
 302   vecZk(0,0) =dedxRef/dedx0;
 303   measR(0,0) =s0*s0*dedxRef*dedxRef/(dedx0*dedx0);
 304   //
 305   // reset matHk
 306   for (Int_t iel=0;iel<fNelem;iel++)
 307     for (Int_t ip=0;ip<nmeas;ip++) matHk(ip,iel)=0;
 308
 309   //
 310   //
 311   matHk(0, ip0*fNpar+0) = corr0/vecXk(ip0*fNpar,0);
 312   matHk(0, ip0*fNpar+1) = kY0*vecXk(ip0*fNpar,0);
 313   matHk(0, ip0*fNpar+2) = kZ0*vecXk(ip0*fNpar,0);
 314   matHk(0, ip0*fNpar+3) = dR0*vecXk(ip0*fNpar,0);
 315   matHk(0, ip0*fNpar+4) = vecXk(ip0*fNpar,0)/length0;
 316   //
 317   //
 318   vecYk(0,0) = vecZk(0,0)-corr0;               // Innovation or measurement residual
 319   matHkT=matHk.T(); matHk.T();
 320   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 321   matSk.Invert();
 322   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 323   vecXk += matKk*vecYk;                    //  updated vector
 324   covXk2= (mat1-(matKk*matHk));
 325   covXk3 =  covXk2*covXk;
 326   covXk = covXk3;
 327   Double_t chi2 = vecYk(0,0)*vecYk(0,0)*matSk(0,0);
 328   //
 329   //
 330   //
 331   if (debug){ // dump debug info
 332     char chname[1000];
 333     sprintf(chname,"update_%s",GetName());
 334     (*debug)<<chname<<
 335       //
 336       // input parameters
 337       "ip0="<<ip0<<
 338       "dedxRef="<<dedxRef<<
 339       "dedx0="<<dedx0<<
 340       "s0="<<s0<<
 341       "kY0="<<kY0<<
 342       "kZ0="<<kZ0<<
 343       "dR0="<<dR0<<
 344       "chi2="<<chi2<<
 345       "corr0="<<corr0<<
 346       "length0="<<length0<<
 347       //
 348       "vecYk.="<<&vecYk<<
 349       "vecZk.="<<&vecZk<<
 350       "matHk.="<<&matHk<<
 351       "matHkT.="<<&matHkT<<
 352       "matSk.="<<&matSk<<
 353       "matKk.="<<&matKk<<
 354       "covXk2.="<<&covXk2<<
 355       "covXk.="<<&covXk<<
 356       "cov.="<<fCovariance<<
 357       "vecXk.="<<&vecXk<<
 358       "vec.="<<fState<<
 359       "\n";
 360   }
 361   if (chi2<kchi2Cut){
 362     (*fCovariance)=covXk;
 363     (*fState)=vecXk;
 364   }else{
 365     printf("chi2=%f\n",chi2);
 366   }
 367 }
 368
 369 void AliTPCkalmandEdx::AdddEdx(Int_t ip0,Double_t dedx0, Double_t dedxRef){
 370   //
 371   //
 372   //
 373   if (fCounter[ip0]>=fSampleSize) return;
 374   fSample[ip0][fCounter[ip0]]=dedxRef/dedx0;
 375   fCounter[ip0]++;
 376   if (fCounter[ip0]==fSampleSize){
 377     fSampleStat[ip0].ResizeTo(2);
 378     fSampleStat[ip0][0] = TMath::Median(fSampleSize,fSample[ip0].GetMatrixArray());
 379     fSampleStat[ip0][1] = TMath::RMS(fSampleSize,fSample[ip0].GetMatrixArray());
 380     fInit++;
 381     (*fState)(ip0*fNpar,0)    = fSampleStat[ip0][0];
 382     (*fCovariance)(ip0*fNpar,ip0*fNpar) = fSampleStat[ip0][1]*fSampleStat[ip0][1]/fSampleSize;
 383     Double_t med2  = fSampleStat[ip0][0]*fSampleStat[ip0][0];
 384     //
 385     //
 386     (*fCovariance)(ip0*fNpar+1,ip0*fNpar+1)=0.01*med2;
 387     (*fCovariance)(ip0*fNpar+2,ip0*fNpar+2)=0.01*med2;
 388     (*fCovariance)(ip0*fNpar+3,ip0*fNpar+3)=0.01*med2;
 389     (*fCovariance)(ip0*fNpar+4,ip0*fNpar+4)=0.01*med2;
 390   }
 391 }
 392
 393
 394