1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // Class to the track points on the Riemann sphere. Inputs are
21 // the set of id's (volids) of the volumes in which residuals are
22 // calculated to construct a chi2 function to be minimized during
23 // the alignment procedures. For the moment the track extrapolation is
24 // taken at the space-point reference plane. The reference plane is
25 // found using the covariance matrix of the point
26 // (assuming sigma(x)=0 at the reference coordinate system.
28 // Internal usage of AliRieman class for minimization
30 //////////////////////////////////////////////////////////////////////////////
33 #include <TLinearFitter.h>
36 #include <TMatrixDSym.h>
38 #include <TTreeStream.h>
42 #include "AliRieman.h"
43 #include "AliTrackFitterRieman.h"
45 ClassImp(AliTrackFitterRieman)
48 AliTrackFitterRieman::AliTrackFitterRieman():
55 fRieman(new AliRieman(10000)), // allocate rieman
56 fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
59 // default constructor
64 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
65 AliTrackFitter(array,owner),
71 fRieman(new AliRieman(10000)), //allocate rieman
77 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
80 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
81 AliTrackFitter(rieman),
82 fBCorrection(rieman.fBCorrection),
83 fAlpha(rieman.fAlpha),
84 fNUsed(rieman.fNUsed),
86 fMaxDelta(rieman.fMaxDelta),
87 fRieman(new AliRieman(*(rieman.fRieman))),
93 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
96 //_____________________________________________________________________________
97 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
100 // Assignment operator
102 if(this==&rieman) return *this;
103 ((AliTrackFitter *)this)->operator=(rieman);
105 fBCorrection = rieman.fBCorrection;
106 fAlpha = rieman.fAlpha;
107 fNUsed = rieman.fNUsed;
108 fConv = rieman.fConv;
109 fMaxDelta = rieman.fMaxDelta;
110 fRieman = new AliRieman(*(rieman.fRieman));
112 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
117 AliTrackFitterRieman::~AliTrackFitterRieman(){
125 void AliTrackFitterRieman::Reset()
127 // Reset the track parameters and
130 AliTrackFitter::Reset();
137 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
138 AliGeomManager::ELayerID layerRangeMin,
139 AliGeomManager::ELayerID layerRangeMax)
141 // Fit the track points. The method takes as an input
142 // the set of id's (volids) of the volumes in which
143 // one wants to calculate the residuals.
144 // The following parameters are used to define the
145 // range of volumes to be used in the fitting
146 // As a result two AliTrackPointArray's obects are filled.
147 // The first one contains the space points with
148 // volume id's from volids list. The second array of points represents
149 // the track extrapolations corresponding to the space points
150 // in the first array. The two arrays can be used to find
151 // the residuals in the volids and consequently construct a
152 // chi2 function to be minimized during the alignment
153 // procedures. For the moment the track extrapolation is taken
154 // at the space-point reference plane. The reference plane is
155 // found using the covariance matrix of the point
156 // (assuming sigma(x)=0 at the reference coordinate system.
157 Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
159 // Float_t debugRatio = 1./(1.+debugLevel);
160 Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;
162 Int_t npoints = fPoints->GetNPoints();
163 if ( npoints<fMinNPoints) return kFALSE;
166 if (volIdsFit != 0x0) {
168 Int_t countPoint = 0;
169 for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
170 if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
172 if (volIdsFit != 0x0) {
173 if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
177 if (countPoint==0) return kFALSE;
178 if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
184 if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
185 Int_t nVol = volIds->GetSize();
186 Int_t nVolFit = volIdsFit->GetSize();
187 Int_t volId = volIds->At(0);
188 (*fDebugStream)<<"PInput"<<
189 "NPoints="<<npoints<< // number of points
190 "VolId="<<volId<< // first vol ID
191 "NVol="<<nVol<< // number of volumes
192 "NvolFit="<<nVolFit<< // number of volumes to fit
193 "fPoints.="<<fPoints<< // input points
197 Bool_t isAlphaCalc = kFALSE;
198 AliTrackPoint p,plocal;
199 // fPoints->GetPoint(p,0);
200 // fAlpha = TMath::ATan2(p.GetY(),p.GetX());
204 Int_t *pindex = new Int_t[npoints];
205 for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
207 fPoints->GetPoint(p,ipoint);
208 UShort_t iVolId = p.GetVolumeID();
209 if (FindVolId(volIds,iVolId)) {
210 pindex[npVolId] = ipoint;
213 if (volIdsFit != 0x0) {
214 if (!FindVolId(volIdsFit,iVolId)) continue;
217 if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
218 iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
219 AliGeomManager::LayerSize(layerRangeMax))) continue;
222 fAlpha = p.GetAngle();
225 plocal = p.Rotate(fAlpha);
226 if (TMath::Abs(plocal.GetX())>500 || TMath::Abs(plocal.GetX())<2 || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
227 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
230 printf("Problematic point\n");
231 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
233 AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
234 TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
236 // fNUsed++; AddPoint should be responsible
239 if (npVolId == 0 || fNUsed < fMinNPoints) {
251 if ((fParams[0] == 0) ||
252 ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
258 if (fNUsed < fMinNPoints) {
263 fPVolId = new AliTrackPointArray(npVolId);
264 fPTrack = new AliTrackPointArray(npVolId);
266 for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
268 Int_t index = pindex[ipoint];
269 fPoints->GetPoint(p,index);
270 if (GetPCA(p,p2) && (
271 TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
272 TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
273 TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
275 Float_t xyz[3],xyz2[3];
276 p.GetXYZ(xyz); p2.GetXYZ(xyz2);
277 // printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
278 fPVolId->AddPoint(ipoint,&p);
279 fPTrack->AddPoint(ipoint,&p2);
281 // what should be default bahavior -
291 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
295 AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
296 AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
297 AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
298 AliRieman * residual = fRieman->MakeResiduals();
299 for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
300 AliTrackPoint p0, p0local;
301 AliTrackPoint pFit, pFitlocal, pFitLocalE;
302 fPVolId->GetPoint(p0,ipoint);
303 Float_t lAngle = p0.GetAngle();
304 p0local= p0.MasterToLocal();
305 fPTrack->GetPoint(pFit,ipoint);
306 pFitlocal= pFit.Rotate(lAngle);
308 Float_t xyz[3], cov[6];
309 xyz[0] = pFitlocal.GetX();
310 xyz[1] = pFitlocal.GetY();
311 xyz[2] = pFitlocal.GetZ();
312 for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
313 cov[3] = GetErrY2at(xyz[0]);
314 cov[5] = GetErrZ2at(xyz[0]);
315 pFitLocalE.SetXYZ(xyz,cov);
317 lPVolId->AddPoint(ipoint,&p0local);
318 lPTrack->AddPoint(ipoint,&pFitlocal);
319 lPTrackE->AddPoint(ipoint,&pFitLocalE);
324 Int_t nVol = volIds->GetSize();
325 Int_t nVolFit = volIdsFit->GetSize();
326 Int_t volId = volIds->At(0);
328 Int_t layer = AliGeomManager::VolUIDToLayer(volId,modId);
329 Int_t volIdFit = volIdsFit->At(0);
331 Int_t layerFit = AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
333 (*fDebugStream)<<"Fit"<<
334 "VolId="<<volId<< // volume ID
335 "Layer="<<layer<< // layer ID
336 "Module="<<modId<< // module ID
337 "LayerFit="<<layerFit<< // layer ID fit
338 "ModuleFit="<<modIdFit<< // module ID fit
339 "NVol="<<nVol<< // number of volumes
340 "NvolFit="<<nVolFit<< // number of volumes to fit
341 "Points0.="<<fPVolId<< // original points
342 "Points1.="<<fPTrack<< // fitted points
343 "LPoints0.="<<lPVolId<< // original points - local frame
344 "LPoints1.="<<lPTrack<< // fitted points - local frame
345 "LPointsE.="<<lPTrackE<< // fitted points with ext error - local frame
346 "Rieman.="<<this<< // original rieman fit
347 "Res.="<<residual<< // residuals of rieman fit
359 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
362 // add point to rieman fitter
364 fRieman->AddPoint(x,y,z,sy,sz);
365 fNUsed = fRieman->GetN();
370 Bool_t AliTrackFitterRieman::Update(){
376 if (fRieman->IsValid()){
377 for (Int_t ipar=0; ipar<6; ipar++){
378 fParams[ipar] = fRieman->GetParam()[ipar];
380 fChi2 = fRieman->GetChi2();
381 fNdf = fRieman->GetN()- 2;
382 fNUsed = fRieman->GetN();
388 TLinearFitter fitY(3,"pol2");
389 TLinearFitter fitZ(3,"pol2");
390 for (Int_t ip=0; ip<fRieman->GetN();ip++){
391 Double_t x = fRieman->GetX()[ip];
392 fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
393 fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);
397 for (Int_t iparam=0; iparam<3; iparam++){
398 fCorrY[iparam]=fitY.GetParameter(iparam);
399 fCorrZ[iparam]=fitZ.GetParameter(iparam);
401 fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);
402 fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
409 //_____________________________________________________________________________
410 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
413 // Get the closest to a given spacepoint track trajectory point
414 // Look for details in the description of the Fit() method
416 if (!fConv) return kFALSE;
418 // First X and Y coordinates
419 Double_t sin = TMath::Sin(fAlpha);
420 Double_t cos = TMath::Cos(fAlpha);
422 // fParam[1] = -x0/y0
423 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
424 if (fParams[0] == 0) return kFALSE;
425 // Track parameters in the global coordinate system
426 Double_t x0 = -fParams[1]/fParams[0]*cos - 1./fParams[0]*sin;
427 Double_t y0 = 1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
428 if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
429 Double_t r = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
432 // Define space-point refence plane
433 Double_t alphap = p.GetAngle();
434 Double_t sinp = TMath::Sin(alphap);
435 Double_t cosp = TMath::Cos(alphap);
436 Double_t x = p.GetX()*cosp + p.GetY()*sinp;
437 Double_t y = p.GetY()*cosp - p.GetX()*sinp;
438 Double_t x0p= x0*cosp + y0*sinp;
439 Double_t y0p= y0*cosp - x0*sinp;
440 if ((r*r - (x-x0p)*(x-x0p))<0) {
441 AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
444 Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
445 Double_t y1 = y0p + temp;
446 Double_t y2 = y0p - temp;
447 Double_t yprime = y1;
448 if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
450 // Back to the global coordinate system
451 Double_t xsecond = x*cosp - yprime*sinp;
452 Double_t ysecond = yprime*cosp + x*sinp;
454 // Now Z coordinate and track angles
455 Double_t x2 = xsecond*cos + ysecond*sin;
456 Double_t zsecond = GetZat(x2);
457 Double_t dydx = GetDYat(x2);
458 Double_t dzdx = GetDZat(x2);
460 // Fill the cov matrix of the track extrapolation point
461 Double_t cov[6] = {0,0,0,0,0,0};
462 Double_t sigmax = 100*100.;
463 cov[0] = sigmax; cov[1] = sigmax*dydx; cov[2] = sigmax*dzdx;
464 cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
465 cov[5] = sigmax*dzdx*dzdx;
467 Double_t sigmay2 = GetErrY2at(x2);
468 Double_t sigmaz2 = GetErrZ2at(x2);
474 newcov[0] = cov[0]*cos*cos-
477 newcov[1] = cov[1]*(cos*cos-sin*sin)-
478 (cov[3]-cov[0])*sin*cos;
479 newcov[2] = cov[2]*cos-
481 newcov[3] = cov[0]*sin*sin+
484 newcov[4] = cov[4]*cos+
488 p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
489 Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
490 Float_t debugRatio = 1./(1.+debugLevel);
491 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
492 AliTrackPoint lp0(p);
493 AliTrackPoint lp2(p2);
494 AliTrackPoint localp0(p);
495 AliTrackPoint localp2(p2);
496 Float_t lAngle = lp0.GetAngle();
497 localp0 = localp0.Rotate(lAngle);
498 localp2 = localp2.Rotate(lAngle);
500 (*fDebugStream)<<"PCA"<<
501 "P0.="<<&lp0<< //global position
503 "LP0.="<<&localp0<< //local position
510 Double_t AliTrackFitterRieman::GetYat(Double_t x) const {
512 // get y position at given point
514 Double_t correction=0;
515 if (fBCorrection){ // systematic effect correction
516 correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
518 return fRieman->GetYat(x)+correction;
521 Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
523 // get z position at given point
525 Double_t correction=0;
526 if (fBCorrection){ // systematic effect correction
527 correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
529 return fRieman->GetZat(x)+correction;
532 Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
534 // get estimate of extrapolation error
536 Double_t error = fRieman->GetErrY(x);
537 Double_t correction=0;
538 if (fBCorrection){ // everestimate error due systematic effect
540 correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
541 correction *=correction;
543 return TMath::Sqrt(error+correction);
546 Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
548 // get estimate of extrapolation error
550 Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
551 Double_t correction=0;
554 correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
555 correction*= correction;
557 return TMath::Sqrt(error+correction);