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
57 fMaxPointRadius(500.),
58 fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
61 // default constructor
63 fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
64 fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
68 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
69 AliTrackFitter(array,owner),
75 fRieman(new AliRieman(10000)), //allocate rieman
77 fMaxPointRadius(500.),
83 fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
84 fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
85 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
88 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
89 AliTrackFitter(rieman),
90 fBCorrection(rieman.fBCorrection),
91 fAlpha(rieman.fAlpha),
92 fNUsed(rieman.fNUsed),
94 fMaxDelta(rieman.fMaxDelta),
95 fRieman(new AliRieman(*(rieman.fRieman))),
96 fMinPointRadius(rieman.fMinPointRadius),
97 fMaxPointRadius(rieman.fMaxPointRadius),
103 fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
104 fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
105 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
108 //_____________________________________________________________________________
109 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
112 // Assignment operator
114 if(this==&rieman) return *this;
115 ((AliTrackFitter *)this)->operator=(rieman);
117 fBCorrection = rieman.fBCorrection;
118 fAlpha = rieman.fAlpha;
119 fNUsed = rieman.fNUsed;
120 fConv = rieman.fConv;
121 fMaxDelta = rieman.fMaxDelta;
122 fRieman = new AliRieman(*(rieman.fRieman));
123 fMinPointRadius = rieman.fMinPointRadius;
124 fMaxPointRadius = rieman.fMaxPointRadius;
126 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
131 AliTrackFitterRieman::~AliTrackFitterRieman(){
139 void AliTrackFitterRieman::Reset()
141 // Reset the track parameters and
144 AliTrackFitter::Reset();
151 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
152 AliGeomManager::ELayerID layerRangeMin,
153 AliGeomManager::ELayerID layerRangeMax)
155 // Fit the track points. The method takes as an input
156 // the set of id's (volids) of the volumes in which
157 // one wants to calculate the residuals.
158 // The following parameters are used to define the
159 // range of volumes to be used in the fitting
160 // As a result two AliTrackPointArray's obects are filled.
161 // The first one contains the space points with
162 // volume id's from volids list. The second array of points represents
163 // the track extrapolations corresponding to the space points
164 // in the first array. The two arrays can be used to find
165 // the residuals in the volids and consequently construct a
166 // chi2 function to be minimized during the alignment
167 // procedures. For the moment the track extrapolation is taken
168 // at the space-point reference plane. The reference plane is
169 // found using the covariance matrix of the point
170 // (assuming sigma(x)=0 at the reference coordinate system.
171 Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
173 // Float_t debugRatio = 1./(1.+debugLevel);
174 Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;
177 AliError("Track points array not available! Exiting...");
180 Int_t npoints = fPoints->GetNPoints();
181 if ( npoints<fMinNPoints) return kFALSE;
184 if (volIdsFit != 0x0) {
186 Int_t countPoint = 0;
187 for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
188 if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
190 if (volIdsFit != 0x0) {
191 if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
195 if (countPoint==0) return kFALSE;
196 if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
202 if (fPoints && volIdsFit && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
203 Int_t nVol = volIds->GetSize();
204 Int_t nVolFit = volIdsFit->GetSize();
205 Int_t volId = volIds->At(0);
206 (*fDebugStream)<<"PInput"<<
207 "NPoints="<<npoints<< // number of points
208 "VolId="<<volId<< // first vol ID
209 "NVol="<<nVol<< // number of volumes
210 "NvolFit="<<nVolFit<< // number of volumes to fit
211 "fPoints.="<<fPoints<< // input points
215 Bool_t isAlphaCalc = kFALSE;
216 AliTrackPoint p,plocal;
217 // fPoints->GetPoint(p,0);
218 // fAlpha = TMath::ATan2(p.GetY(),p.GetX());
222 Int_t *pindex = new Int_t[npoints];
223 for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
225 fPoints->GetPoint(p,ipoint);
226 UShort_t iVolId = p.GetVolumeID();
227 if (FindVolId(volIds,iVolId)) {
228 pindex[npVolId] = ipoint;
231 if (volIdsFit != 0x0) {
232 if (!FindVolId(volIdsFit,iVolId)) continue;
235 if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
236 iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
237 AliGeomManager::LayerSize(layerRangeMax))) continue;
240 fAlpha = p.GetAngle();
243 plocal = p.Rotate(fAlpha);
244 if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())<fMinPointRadius || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
245 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
248 printf("Problematic point\n");
249 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
251 AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
252 TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
254 // fNUsed++; AddPoint should be responsible
257 if (npVolId == 0 || fNUsed < fMinNPoints) {
269 if ((fParams[0] == 0) ||
270 ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
276 if (fNUsed < fMinNPoints) {
281 fPVolId = new AliTrackPointArray(npVolId);
282 fPTrack = new AliTrackPointArray(npVolId);
284 for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
286 Int_t index = pindex[ipoint];
287 fPoints->GetPoint(p,index);
288 if (GetPCA(p,p2) && (
289 TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
290 TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
291 TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
293 Float_t xyz[3],xyz2[3];
294 p.GetXYZ(xyz); p2.GetXYZ(xyz2);
295 // 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]);
296 fPVolId->AddPoint(ipoint,&p);
297 fPTrack->AddPoint(ipoint,&p2);
299 // what should be default bahavior -
309 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
313 AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
314 AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
315 AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
316 AliRieman * residual = fRieman->MakeResiduals();
317 for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
318 AliTrackPoint p0, p0local;
319 AliTrackPoint pFit, pFitlocal, pFitLocalE;
320 fPVolId->GetPoint(p0,ipoint);
321 Float_t lAngle = p0.GetAngle();
322 p0local= p0.MasterToLocal();
323 fPTrack->GetPoint(pFit,ipoint);
324 pFitlocal= pFit.Rotate(lAngle);
326 Float_t xyz[3], cov[6];
327 xyz[0] = pFitlocal.GetX();
328 xyz[1] = pFitlocal.GetY();
329 xyz[2] = pFitlocal.GetZ();
330 for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
331 cov[3] = GetErrY2at(xyz[0]);
332 cov[5] = GetErrZ2at(xyz[0]);
333 pFitLocalE.SetXYZ(xyz,cov);
335 lPVolId->AddPoint(ipoint,&p0local);
336 lPTrack->AddPoint(ipoint,&pFitlocal);
337 lPTrackE->AddPoint(ipoint,&pFitLocalE);
342 Int_t nVol = volIds->GetSize();
343 Int_t nVolFit = volIdsFit->GetSize();
344 Int_t volId = volIds->At(0);
346 Int_t layer = AliGeomManager::VolUIDToLayer(volId,modId);
347 Int_t volIdFit = volIdsFit->At(0);
349 Int_t layerFit = AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
351 (*fDebugStream)<<"Fit"<<
352 "VolId="<<volId<< // volume ID
353 "Layer="<<layer<< // layer ID
354 "Module="<<modId<< // module ID
355 "LayerFit="<<layerFit<< // layer ID fit
356 "ModuleFit="<<modIdFit<< // module ID fit
357 "NVol="<<nVol<< // number of volumes
358 "NvolFit="<<nVolFit<< // number of volumes to fit
359 "Points0.="<<fPVolId<< // original points
360 "Points1.="<<fPTrack<< // fitted points
361 "LPoints0.="<<lPVolId<< // original points - local frame
362 "LPoints1.="<<lPTrack<< // fitted points - local frame
363 "LPointsE.="<<lPTrackE<< // fitted points with ext error - local frame
364 "Rieman.="<<this<< // original rieman fit
365 "Res.="<<residual<< // residuals of rieman fit
377 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
380 // add point to rieman fitter
382 fRieman->AddPoint(x,y,z,sy,sz);
383 fNUsed = fRieman->GetN();
388 Bool_t AliTrackFitterRieman::Update(){
394 if (fRieman->IsValid()){
395 for (Int_t ipar=0; ipar<6; ipar++){
396 fParams[ipar] = fRieman->GetParam()[ipar];
398 fChi2 = fRieman->GetChi2();
399 fNdf = fRieman->GetN()- 2;
400 fNUsed = fRieman->GetN();
406 TLinearFitter fitY(3,"pol2");
407 TLinearFitter fitZ(3,"pol2");
408 for (Int_t ip=0; ip<fRieman->GetN();ip++){
409 Double_t x = fRieman->GetX()[ip];
410 fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
411 fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);
415 for (Int_t iparam=0; iparam<3; iparam++){
416 fCorrY[iparam]=fitY.GetParameter(iparam);
417 fCorrZ[iparam]=fitZ.GetParameter(iparam);
419 fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);
420 fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
427 //_____________________________________________________________________________
428 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
431 // Get the closest to a given spacepoint track trajectory point
432 // Look for details in the description of the Fit() method
434 if (!fConv) return kFALSE;
436 // First X and Y coordinates
437 Double_t sin = TMath::Sin(fAlpha);
438 Double_t cos = TMath::Cos(fAlpha);
440 // fParam[1] = -x0/y0
441 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
442 if (fParams[0] == 0) return kFALSE;
443 // Track parameters in the global coordinate system
444 Double_t x0 = -fParams[1]/fParams[0]*cos - 1./fParams[0]*sin;
445 Double_t y0 = 1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
446 if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
447 Double_t r = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
450 // Define space-point refence plane
451 Double_t alphap = p.GetAngle();
452 Double_t sinp = TMath::Sin(alphap);
453 Double_t cosp = TMath::Cos(alphap);
454 Double_t x = p.GetX()*cosp + p.GetY()*sinp;
455 Double_t y = p.GetY()*cosp - p.GetX()*sinp;
456 Double_t x0p= x0*cosp + y0*sinp;
457 Double_t y0p= y0*cosp - x0*sinp;
458 if ((r*r - (x-x0p)*(x-x0p))<0) {
459 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));
462 Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
463 Double_t y1 = y0p + temp;
464 Double_t y2 = y0p - temp;
465 Double_t yprime = y1;
466 if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
468 // Back to the global coordinate system
469 Double_t xsecond = x*cosp - yprime*sinp;
470 Double_t ysecond = yprime*cosp + x*sinp;
472 // Now Z coordinate and track angles
473 Double_t x2 = xsecond*cos + ysecond*sin;
474 Double_t zsecond = GetZat(x2);
475 Double_t dydx = GetDYat(x2);
476 Double_t dzdx = GetDZat(x2);
478 // Fill the cov matrix of the track extrapolation point
479 Double_t cov[6] = {0,0,0,0,0,0};
480 Double_t sigmax = 100*100.;
481 cov[0] = sigmax; cov[1] = sigmax*dydx; cov[2] = sigmax*dzdx;
482 cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
483 cov[5] = sigmax*dzdx*dzdx;
485 Double_t sigmay2 = GetErrY2at(x2);
486 Double_t sigmaz2 = GetErrZ2at(x2);
492 newcov[0] = cov[0]*cos*cos-
495 newcov[1] = cov[1]*(cos*cos-sin*sin)-
496 (cov[3]-cov[0])*sin*cos;
497 newcov[2] = cov[2]*cos-
499 newcov[3] = cov[0]*sin*sin+
502 newcov[4] = cov[4]*cos+
506 p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
507 Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
508 Float_t debugRatio = 1./(1.+debugLevel);
509 if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
510 AliTrackPoint lp0(p);
511 AliTrackPoint lp2(p2);
512 AliTrackPoint localp0(p);
513 AliTrackPoint localp2(p2);
514 Float_t lAngle = lp0.GetAngle();
515 localp0 = localp0.Rotate(lAngle);
516 localp2 = localp2.Rotate(lAngle);
518 (*fDebugStream)<<"PCA"<<
519 "P0.="<<&lp0<< //global position
521 "LP0.="<<&localp0<< //local position
528 Double_t AliTrackFitterRieman::GetYat(Double_t x) const {
530 // get y position at given point
532 Double_t correction=0;
533 if (fBCorrection){ // systematic effect correction
534 correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
536 return fRieman->GetYat(x)+correction;
539 Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
541 // get z position at given point
543 Double_t correction=0;
544 if (fBCorrection){ // systematic effect correction
545 correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
547 return fRieman->GetZat(x)+correction;
550 Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
552 // get estimate of extrapolation error
554 Double_t error = fRieman->GetErrY(x);
555 Double_t correction=0;
556 if (fBCorrection){ // everestimate error due systematic effect
558 correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
559 correction *=correction;
561 return TMath::Sqrt(error+correction);
564 Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
566 // get estimate of extrapolation error
568 Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
569 Double_t correction=0;
572 correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
573 correction*= correction;
575 return TMath::Sqrt(error+correction);
578 void AliTrackFitterRieman::SetParam(Int_t i, Double_t par) {
579 if (i<0 || i>5) return;
581 fRieman->GetParam()[i]=par;