Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / STEER / STEER / AliTrackResidualsChi2.cxx
CommitLineData
98937d93 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
93e1dfe1 22#include <TVirtualFitter.h>
98937d93 23#include <TGeoMatrix.h>
24
25#include "AliAlignObj.h"
26#include "AliTrackPointArray.h"
27#include "AliTrackResidualsChi2.h"
28
4c84f475 29
30
a4ef2f82 31void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
46ae650f 32
93e1dfe1 33
98937d93 34ClassImp(AliTrackResidualsChi2)
35
93e1dfe1 36
98937d93 37//______________________________________________________________________________
a4ef2f82 38void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
98937d93 39{
40 // This function is called by minuit
41 // The corresponding member method is called
42 // using SetObjectFit/GetObjectFit methods of TMinuit
93e1dfe1 43 AliTrackResidualsChi2* dummy = (AliTrackResidualsChi2 *)TVirtualFitter::GetFitter()->GetObjectFit();
98937d93 44 dummy->Chi2(npar, gin, f, par, iflag);
45}
46
47//______________________________________________________________________________
48Bool_t AliTrackResidualsChi2::Minimize()
49{
50 // Implementation of Chi2 based minimization
51 // of track residuala sum
93e1dfe1 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.
46ae650f 57 Double_t arglist[10];
58 Int_t ierflg = 0;
93e1dfe1 59 TVirtualFitter *fitter = TVirtualFitter::Fitter(this,6); //initialize TMinuit
46ae650f 60 arglist[0] = -1;
93e1dfe1 61 ierflg = fitter->ExecuteCommand("SET PRINT", arglist, 1);
46ae650f 62
a4ef2f82 63 fitter->SetFCN(TrackResidualsChi2Fcn);
98937d93 64
98937d93 65 arglist[0] = 1;
93e1dfe1 66 ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
98937d93 67
68 // Set starting values and step sizes for parameters
46ae650f 69 Double_t pars[6] = {0,0,0,0,0,0};
29317b68 70 for(Int_t i=0;i<6;i++) if(fBFixed[i]) pars[i]=fFixed[i];
46ae650f 71 Double_t step[6] = {0.0001,0.0001,0.0001,0.0001,0.0001,0.0001};
93e1dfe1 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);
98937d93 78
29317b68 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
98937d93 87 // Now ready for minimization step
88 arglist[0] = 500;
89 arglist[1] = 1.;
93e1dfe1 90 fitter->ExecuteCommand("MIGRAD", arglist ,2);
98937d93 91
92 // Print results
93 Double_t amin,edm,errdef;
94 Int_t nvpar,nparx,icstat;
93e1dfe1 95 icstat = fitter->GetStats(amin,edm,errdef,nvpar,nparx);
98937d93 96
4c84f475 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
98937d93 137 return kTRUE;
138}
139
140//______________________________________________________________________________
141void 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]);
98937d93 148
149 AliTrackPoint p1,p2;
46ae650f 150
151 Bool_t count = kFALSE;
152 if (fNdf == 0) count = kTRUE;
98937d93 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);
46ae650f 158 fAlignObj->Transform(p1);
98937d93 159 fTrackArray[itrack]->GetPoint(p2,ipoint);
cc345ce3 160 Float_t residual = p2.GetResidual(p1,kTRUE);
46ae650f 161 chi2 += residual;
162 if (count) fNdf += 3;
98937d93 163 }
164 }
165 f = chi2;
166}
4c84f475 167
168