STEER/AliTrackResidualsChi2.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 //   Implementation of the derived class for track residuals
  18 //   based on Chi2 minimization
  19 //
  20 //-----------------------------------------------------------------
  21
  22 #include <TVirtualFitter.h>
  23 #include <TGeoMatrix.h>
  24
  25 #include "AliAlignObj.h"
  26 #include "AliTrackPointArray.h"
  27 #include "AliTrackResidualsChi2.h"
  28
  29
  30
  31 void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
  32
  33
  34 ClassImp(AliTrackResidualsChi2)
  35
  36
  37 //______________________________________________________________________________
  38 void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
  39 {
  40   // This function is called by minuit
  41   // The corresponding member method is called
  42   // using SetObjectFit/GetObjectFit methods of TMinuit
  43   AliTrackResidualsChi2* dummy = (AliTrackResidualsChi2 *)TVirtualFitter::GetFitter()->GetObjectFit();
  44   dummy->Chi2(npar, gin, f, par, iflag);
  45 }
  46
  47 //______________________________________________________________________________
  48 Bool_t AliTrackResidualsChi2::Minimize()
  49 {
  50   // Implementation of Chi2 based minimization
  51   // of track residuala sum
  52   //
  53   // CHOLM: Modified to use the TVirtualFitter interface only.  This
  54   // makes it possible other fitters than TMinuit (say TMinuit2,
  55   // TFumili) by simply changing ones configuration file, or setting
  56   // it up in a script.  Much more robust and flexible.
  57   Double_t arglist[10];
  58   Int_t ierflg = 0;
  59   TVirtualFitter *fitter = TVirtualFitter::Fitter(this,6);  //initialize TMinuit
  60   arglist[0] = -1;
  61   ierflg = fitter->ExecuteCommand("SET PRINT", arglist, 1);
  62
  63   fitter->SetFCN(TrackResidualsChi2Fcn);
  64
  65   arglist[0] = 1;
  66   ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
  67
  68   // Set starting values and step sizes for parameters
  69   Double_t pars[6] = {0,0,0,0,0,0};
  70   for(Int_t i=0;i<6;i++) if(fBFixed[i]) pars[i]=fFixed[i];
  71   Double_t step[6] = {0.0001,0.0001,0.0001,0.0001,0.0001,0.0001};
  72   ierflg = fitter->SetParameter(0, "dx", pars[0], step[0], 0,0);
  73   ierflg = fitter->SetParameter(1, "dy", pars[1], step[1], 0,0);
  74   ierflg = fitter->SetParameter(2, "dz", pars[2], step[2], 0,0);
  75   ierflg = fitter->SetParameter(3, "psi", pars[3], step[3], 0,0);
  76   ierflg = fitter->SetParameter(4, "theta", pars[4], step[4], 0,0);
  77   ierflg = fitter->SetParameter(5, "phi", pars[5], step[5], 0,0);
  78
  79   // Fix parameters
  80   if(fBFixed[0]) {printf("Fixing dx=%f\n",pars[0]); fitter->FixParameter(0);}
  81   if(fBFixed[1]) {printf("Fixing dy=%f\n",pars[1]); fitter->FixParameter(1);}
  82   if(fBFixed[2]) {printf("Fixing dz=%f\n",pars[2]); fitter->FixParameter(2);}
  83   if(fBFixed[3]) {printf("Fixing psi=%f\n",pars[3]); fitter->FixParameter(3);}
  84   if(fBFixed[4]) {printf("Fixing theta=%f\n",pars[4]); fitter->FixParameter(4);}
  85   if(fBFixed[5]) {printf("Fixing phi=%f\n",pars[5]); fitter->FixParameter(5);}
  86
  87   // Now ready for minimization step
  88   arglist[0] = 500;
  89   arglist[1] = 1.;
  90   fitter->ExecuteCommand("MIGRAD", arglist ,2);
  91
  92   // Print results
  93   Double_t amin,edm,errdef;
  94   Int_t nvpar,nparx,icstat;
  95   icstat = fitter->GetStats(amin,edm,errdef,nvpar,nparx);
  96
  97    //Construct the covariance matrix for AlignObj
  98   Double_t *cov=fitter->GetCovarianceMatrix();
  99   Int_t unfixedparam=6;
 100   Int_t fixedparamat[6]={0,0,0,0,0,0};
 101   if(fBFixed[0]==kTRUE){
 102     unfixedparam--;
 103     fixedparamat[0]=1;
 104   }
 105   for(Int_t j=1;j<6;j++){
 106     if(fBFixed[j]==kTRUE){
 107       unfixedparam--;
 108       fixedparamat[j]=fixedparamat[j-1]+1;
 109     }
 110     else fixedparamat[j]=fixedparamat[j-1];
 111   }
 112
 113   Double_t cov2[36];
 114   for(Int_t i=0;i<6;i++){
 115     for(Int_t j=0;j<6;j++){
 116       if(fBFixed[i]==kTRUE||fBFixed[j]==kTRUE){
 117         cov2[i+6*j]=0.;
 118       }
 119       else cov2[i+6*j]=cov[i-fixedparamat[i]+6*(j-fixedparamat[j])];
 120
 121     }
 122   }
 123
 124   Double_t covmatrarray[21];
 125   for(Int_t j=0;j<6;j++){
 126     for(Int_t i=j;i<6;i++){
 127       covmatrarray[i*(i+1)/2+j]=cov2[i+6*j];
 128     }
 129   }
 130
 131   //  printf("covar 2:  %.10f  \n   %.10f ; %.10f ;  \n; %.10f;  %.10f;  %.10f; \n  %.10f; %.10f;  %.10f;  %.10f;  \n  %.10f ; %.10f  ;  %.10f  ;  %.10f  ;  %.10f \n %.10f ;  %.10f ;  %.10f  ;  %.10f  ;  %.10f  ;  %.10f   ;\n",covmatrarray[0],covmatrarray[1],covmatrarray[2],covmatrarray[3],covmatrarray[4],covmatrarray[5],covmatrarray[6],covmatrarray[7],covmatrarray[8],covmatrarray[9],covmatrarray[10],covmatrarray[11],covmatrarray[12],covmatrarray[13],covmatrarray[14],covmatrarray[15],covmatrarray[16],covmatrarray[17],covmatrarray[18],covmatrarray[19],covmatrarray[20]);
 132
 133
 134   fAlignObj->SetCorrMatrix(covmatrarray);
 135   fChi2 = amin; fNdf -= nvpar;
 136
 137   return kTRUE;
 138 }
 139
 140 //______________________________________________________________________________
 141 void AliTrackResidualsChi2::Chi2(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */)
 142 {
 143   // Chi2 function to be minimized
 144   // Sums all the track residuals
 145   Double_t chi2 = 0;
 146
 147   fAlignObj->SetPars(par[0],par[1],par[2],par[3],par[4],par[5]);
 148
 149   AliTrackPoint p1,p2;
 150
 151   Bool_t count = kFALSE;
 152   if (fNdf == 0) count = kTRUE;
 153
 154   for (Int_t itrack = 0; itrack < fLast; itrack++) {
 155     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
 156     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
 157       fVolArray[itrack]->GetPoint(p1,ipoint);
 158       fAlignObj->Transform(p1);
 159       fTrackArray[itrack]->GetPoint(p2,ipoint);
 160       Float_t residual = p2.GetResidual(p1,kTRUE);
 161       chi2 += residual;
 162       if (count) fNdf += 3;
 163     }
 164   }
 165   f = chi2;
 166 }
 167
 168