HLT/src/AliL3Fitter.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <math.h>
   7
   8 #include "AliL3Defs.h"
   9 #include "AliL3Logging.h"
  10 #include "AliL3Fitter.h"
  11 #include "AliL3Vertex.h"
  12 #include "AliL3Track.h"
  13 #include "AliL3SpacePointData.h"
  14 #include "AliL3MemHandler.h"
  15
  16 //_____________________________________________________________
  17 // AliL3Fitter
  18 //
  19 // Fit class HLT
  20
  21 ClassImp(AliL3Fitter)
  22
  23 AliL3Fitter::AliL3Fitter(AliL3Vertex *vertex)
  24 {
  25   //constructor
  26   fTrack=0;
  27   fVertex = vertex;
  28   BFACT = 0.0029980;
  29   fVertexConstraint=kTRUE;
  30 }
  31
  32 void AliL3Fitter::LoadClusters(Char_t *path)
  33 {
  34   Char_t fname[256];
  35   AliL3MemHandler *clusterfile[36][6];
  36   for(Int_t s=0; s<=35; s++)
  37     {
  38       for(Int_t p=0; p<6; p++)
  39         {
  40           clusterfile[s][p] = new AliL3MemHandler();
  41           sprintf(fname,"%spoints_%d_%d.raw",path,s,p);
  42           if(!clusterfile[s][p]->SetBinaryInput(fname))
  43             {
  44               delete clusterfile[s][p];
  45               clusterfile[s][p] = 0;
  46               continue;
  47             }
  48           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
  49           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
  50           clusterfile[s][p]->CloseBinaryInput();
  51         }
  52     }
  53 }
  54
  55 Int_t AliL3Fitter::FitHelix(AliL3Track *track)
  56 {
  57   fTrack = track;
  58   if(FitCircle())
  59     {
  60       LOG(AliL3Log::kError,"AliL3Fitter::FitHelix","TrackFit")<<AliL3Log::kDec<<
  61         "Problems during circle fit"<<ENDLOG;
  62       return 1;
  63     }
  64   if(FitLine())
  65     {
  66       LOG(AliL3Log::kError,"AliL3Fitter::FitHelix","TrackFit")<<AliL3Log::kDec<<
  67         "Problems during line fit"<<ENDLOG;
  68       return 1;
  69     }
  70   return 0;
  71 }
  72
  73 Int_t AliL3Fitter::FitCircle()
  74 {
  75   //-----------------------------------------------------------------
  76   //Fits circle parameters using algorithm
  77   //described by ChErnov and Oskov in Computer Physics
  78   //Communications.
  79   //
  80   //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
  81   //Moved to C by Pablo Yepes
  82   //Moved to AliROOT by ASV.
  83   //------------------------------------------------------------------
  84
  85   Double_t wsum  = 0.0 ;
  86   Double_t xav   = 0.0 ;
  87   Double_t yav   = 0.0 ;
  88
  89   //
  90   //     Loop over hits calculating average
  91
  92   Double_t fXYWeight[(fTrack->GetNHits())];
  93   UInt_t *hitnum = fTrack->GetHitNumbers();
  94   for(Int_t i=0; i<fTrack->GetNHits(); i++)
  95     {
  96       UInt_t id = hitnum[i];
  97       Int_t slice = (id>>25) & 0x7f;
  98       Int_t patch = (id>>22) & 0x7;
  99       UInt_t pos = id&0x3fffff;
 100
 101       AliL3SpacePointData *points = fClusters[slice][patch];
 102       fXYWeight[i] = 1./ (Double_t)(points[pos].fXYErr*points[pos].fXYErr + points[pos].fXYErr*points[pos].fXYErr);
 103       wsum += fXYWeight[i];
 104       xav += fXYWeight[i]*points[pos].fX;
 105       yav += fXYWeight[i]*points[pos].fY;
 106
 107     }
 108   if (fVertexConstraint == kTRUE)
 109     {
 110       wsum += fVertex->GetXYWeight() ;
 111       xav  += fVertex->GetX() ;
 112       yav  += fVertex->GetY() ;
 113     }
 114
 115   xav = xav / wsum ;
 116   yav = yav / wsum ;
 117 //
 118 //  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
 119 //
 120   Double_t xxav  = 0.0 ;
 121   Double_t xyav  = 0.0 ;
 122   Double_t yyav  = 0.0 ;
 123   Double_t xi, yi ;
 124
 125   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 126     {
 127       UInt_t id = hitnum[i];
 128       Int_t slice = (id>>25) & 0x7f;
 129       Int_t patch = (id>>22) & 0x7;
 130       UInt_t pos = id&0x3fffff;
 131       AliL3SpacePointData *points = fClusters[slice][patch];
 132       xi = points[pos].fX -xav;
 133       yi        = points[pos].fY - yav ;
 134       xxav     += xi * xi * fXYWeight[i];
 135       xyav     += xi * yi * fXYWeight[i];
 136       yyav     += yi * yi * fXYWeight[i];
 137     }
 138
 139   if (fVertexConstraint == kTRUE)
 140     {
 141       xi        = fVertex->GetX() - xav ;
 142       yi        = fVertex->GetY() - yav ;
 143       xxav     += xi * xi * fVertex->GetXYWeight() ;
 144       xyav     += xi * yi * fVertex->GetXYWeight() ;
 145       yyav     += yi * yi * fVertex->GetXYWeight() ;
 146     }
 147   xxav = xxav / wsum ;
 148   xyav = xyav / wsum ;
 149   yyav = yyav / wsum ;
 150 //
 151 //-->  ROTATE COORDINATES SO THAT <XY> = 0
 152 //
 153 //-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
 154 //-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
 155 //-->  SIGN(S) = SIGN(XYAV)                  >
 156
 157   Double_t a = fabs( xxav - yyav ) ;
 158   Double_t b = 4.0 * xyav * xyav ;
 159
 160   Double_t asqpb  = a * a + b  ;
 161   Double_t rasqpb = sqrt ( asqpb) ;
 162
 163   Double_t splus  = 1.0 + a / rasqpb ;
 164   Double_t sminus = b / (asqpb * splus) ;
 165
 166   splus  = sqrt (0.5 * splus ) ;
 167   sminus = sqrt (0.5 * sminus) ;
 168 //
 169 //->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
 170 //
 171   Double_t sinrot, cosrot ;
 172   if ( xxav <= yyav ) {
 173          cosrot = sminus ;
 174          sinrot = splus  ;
 175   }
 176   else {
 177           cosrot = splus ;
 178           sinrot = sminus ;
 179   }
 180 //
 181 //->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
 182 //
 183   if ( xyav < 0.0 ) sinrot = - sinrot ;
 184 //
 185 //-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
 186 //-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
 187 //-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
 188 //-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
 189 //
 190 //-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
 191 //
 192   if ( cosrot*xav+sinrot*yav < 0.0 ) {
 193           cosrot = -cosrot ;
 194           sinrot = -sinrot ;
 195   }
 196 //
 197 //->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
 198 //
 199   Double_t rrav   = xxav + yyav ;
 200   Double_t rscale = sqrt(rrav) ;
 201
 202   xxav   = 0.0 ;
 203   yyav   = 0.0 ;
 204   xyav   = 0.0 ;
 205   Double_t xrrav = 0.0 ;
 206   Double_t yrrav = 0.0 ;
 207   Double_t rrrrav  = 0.0 ;
 208
 209   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
 210
 211   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 212     {
 213       UInt_t id = hitnum[i];
 214       Int_t slice = (id>>25) & 0x7f;
 215       Int_t patch = (id>>22) & 0x7;
 216       UInt_t pos = id&0x3fffff;
 217       AliL3SpacePointData *points = fClusters[slice][patch];
 218
 219       xold = points[pos].fX - xav ;
 220       yold = points[pos].fY - yav ;
 221       //
 222       //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 223       //
 224       xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 225       yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 226
 227       xixi   = xi * xi ;
 228       yiyi   = yi * yi ;
 229       riri   = xixi + yiyi ;
 230       wiriri = fXYWeight[i] * riri ;
 231
 232       xyav   += fXYWeight[i] * xi * yi ;
 233       xxav   += fXYWeight[i] * xixi ;
 234       yyav   += fXYWeight[i] * yiyi ;
 235
 236       xrrav  += wiriri * xi ;
 237       yrrav  += wiriri * yi ;
 238       rrrrav += wiriri * riri ;
 239     }
 240   //
 241 //   Include vertex if required
 242 //
 243   if (fVertexConstraint == kTRUE)
 244     {
 245       xold = fVertex->GetX() - xav ;
 246         yold = fVertex->GetY() - yav ;
 247         //
 248         //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 249         //
 250         xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 251         yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 252
 253         xixi   = xi * xi ;
 254         yiyi   = yi * yi ;
 255         riri   = xixi + yiyi ;
 256         wiriri = fVertex->GetXYWeight() * riri ;
 257
 258         xyav   += fVertex->GetXYWeight() * xi * yi ;
 259         xxav   += fVertex->GetXYWeight() * xixi ;
 260         yyav   += fVertex->GetXYWeight() * yiyi ;
 261
 262         xrrav  += wiriri * xi ;
 263         yrrav  += wiriri * yi ;
 264         rrrrav += wiriri * riri ;
 265   }
 266   //
 267   //
 268   //
 269   //-->  DIVIDE BY WSUM TO MAKE AVERAGES
 270   //
 271   xxav    = xxav   / wsum ;
 272   yyav    = yyav   / wsum ;
 273   xrrav   = xrrav  / wsum ;
 274   yrrav   = yrrav  / wsum ;
 275   rrrrav  = rrrrav / wsum ;
 276   xyav    = xyav   / wsum ;
 277
 278   Int_t const ntry = 5 ;
 279 //
 280 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
 281 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
 282 //
 283   Double_t xrrxrr = xrrav * xrrav ;
 284   Double_t yrryrr = yrrav * yrrav ;
 285   Double_t rrrrm1 = rrrrav - 1.0  ;
 286   Double_t xxyy   = xxav  * yyav  ;
 287
 288   Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
 289   Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;
 290   Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;
 291   Double_t c4  = - 4.0  ;
 292 //
 293 //-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
 294 //
 295   Double_t c2d =   2.0 * c2 ;
 296   Double_t c4d =   4.0 * c4 ;
 297 //
 298 //-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
 299 //
 300 //   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
 301   Double_t lamda  = 0.0 ;
 302   Double_t dlamda = 0.0 ;
 303 //
 304   Double_t chiscl = wsum * rscale * rscale ;
 305   Double_t dlamax = 0.001 / chiscl ;
 306
 307   Double_t p, pd ;
 308   for ( int itry = 1 ; itry <= ntry ; itry++ ) {
 309      p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
 310      pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
 311      dlamda = -p / pd ;
 312      lamda  = lamda + dlamda ;
 313      if (fabs(dlamda)<   dlamax) break ;
 314   }
 315
 316   Double_t chi2 = (Double_t)(chiscl * lamda) ;
 317
 318   //fTrack->SetChiSq1(chi2);
 319   // Double_t dchisq = chiscl * dlamda ;
 320 //
 321 //-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
 322 //
 323   Double_t h11   = xxav  -     lamda ;
 324   Double_t h14   = xrrav ;
 325   Double_t h22   = yyav  -     lamda ;
 326   Double_t h24   = yrrav ;
 327   Double_t h34   = 1.0   + 2.0*lamda ;
 328   if ( h11 == 0.0 || h22 == 0.0 ){
 329     LOG(AliL3Log::kError,"AliL3Fitter::FitCircle","TrackFit")<<AliL3Log::kDec<<
 330       "Problems fitting circle"<<ENDLOG;
 331     return 1 ;
 332   }
 333   Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
 334
 335   Double_t ratio, kappa, beta ;
 336   if ( fabs(h22) > fabs(h24) ) {
 337      ratio  = h24 / h22 ;
 338      rootsq = ratio * ratio + rootsq ;
 339      kappa = 1.0 / sqrt(rootsq) ;
 340      beta  = - ratio * kappa ;
 341   }
 342   else {
 343      ratio  = h22 / h24 ;
 344      rootsq = 1.0 + ratio * ratio * rootsq ;
 345      beta  = 1.0 / sqrt(rootsq) ;
 346      if ( h24 > 0 ) beta = - beta ;
 347      kappa = -ratio * beta ;
 348   }
 349   Double_t alpha = - (h14/h11) * kappa ;
 350 //
 351 //-->  transform these into the lab coordinate system
 352 //-->  first get kappa and back to real dimensions
 353 //
 354   Double_t kappa1 = kappa / rscale ;
 355   Double_t dbro   = 0.5   / kappa1 ;
 356 //
 357 //-->  next rotate alpha and beta and scale
 358 //
 359   Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
 360   Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
 361 //
 362 //-->  then translate by (xav,yav)
 363 //
 364   Double_t acent  = (double)(xav - alphar) ;
 365   Double_t bcent  = (double)(yav - betar ) ;
 366   Double_t radius = (double)dbro ;
 367 //
 368 //   Get charge
 369 //
 370   Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
 371
 372   fTrack->SetCharge(q);
 373
 374 //
 375 //    Get other track parameters
 376 //
 377   Double_t x0, y0,phi0,r0,psi,pt ;
 378   if ( fVertexConstraint == kTRUE)
 379     {
 380       //flag = 1 ; // primary track flag
 381       x0   = fVertex->GetX() ;
 382       y0   = fVertex->GetY() ;
 383       phi0 = fVertex->GetPhi() ;
 384       r0   = fVertex->GetR() ;
 385       fTrack->SetPhi0(phi0);
 386       fTrack->SetR0(r0);
 387     }
 388   else
 389     {
 390       Int_t lastid=fTrack->GetNHits()-1;
 391       UInt_t id = hitnum[lastid];
 392       Int_t slice = (id>>25) & 0x7f;
 393       Int_t patch = (id>>22) & 0x7;
 394       UInt_t pos = id&0x3fffff;
 395       AliL3SpacePointData *points = fClusters[slice][patch];
 396
 397       //flag =  0 ; // primary track flag
 398       x0   =  points[pos].fX;
 399       y0   =  points[pos].fY;
 400       phi0 =  atan2(points[pos].fY,points[pos].fX);
 401       if ( phi0 < 0 ) phi0 += 2*Pi;
 402       r0   =  sqrt ( points[pos].fX * points[pos].fX + points[pos].fY*points[pos].fY);
 403       fTrack->SetPhi0(phi0);
 404       fTrack->SetR0(r0);
 405     }
 406   //
 407   psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
 408   psi  = psi + q * 0.5F * Pi ;
 409   if ( psi < 0 ) psi = psi + 2*Pi;
 410
 411   pt   = (Double_t)(BFACT * BField * radius ) ;
 412   fTrack->SetPsi(psi);
 413   fTrack->SetPt(pt);
 414   fTrack->SetFirstPoint(x0,y0,0);
 415   //
 416 //    Get errors from fast fit
 417 //
 418   //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
 419 //
 420   return 0 ;
 421
 422 }
 423
 424 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 425 //    Fit Line in s-z plane
 426 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 427 Int_t AliL3Fitter::FitLine ( )
 428 {
 429   //
 430   //Initialization
 431   //
 432   Double_t sum = 0.F ;
 433   Double_t ss  = 0.F ;
 434   Double_t sz  = 0.F ;
 435   Double_t sss = 0.F ;
 436   Double_t ssz = 0.F ;
 437   //
 438   //find sum , sums ,sumz, sumss
 439   //
 440   Double_t dx, dy ;
 441   Double_t radius = (Double_t)(fTrack->GetPt() / ( BFACT * BField ) ) ;
 442
 443   //TObjArray *hits = fTrack->GetHits();
 444   //Int_t num_of_hits = fTrack->GetNumberOfPoints();
 445
 446   Double_t fS[(fTrack->GetNHits())];
 447   UInt_t *hitnum = fTrack->GetHitNumbers();
 448   if (fVertexConstraint==kTRUE)
 449     {
 450       UInt_t id = hitnum[0];
 451       Int_t slice = (id>>25) & 0x7f;
 452       Int_t patch = (id>>22) & 0x7;
 453       UInt_t pos = id&0x3fffff;
 454       AliL3SpacePointData *points = fClusters[slice][patch];
 455
 456       dx = points[pos].fX - fVertex->GetX();
 457       dy = points[pos].fY - fVertex->GetY() ;
 458     }
 459   else
 460     {
 461       UInt_t id = hitnum[0];
 462       Int_t slice = (id>>25) & 0x7f;
 463       Int_t patch = (id>>22) & 0x7;
 464       UInt_t posf = id&0x3fffff;
 465       AliL3SpacePointData *pointsf = fClusters[slice][patch];
 466       id = hitnum[(fTrack->GetNHits()-1)];
 467       slice = (id>>25) & 0x7f;
 468       patch = (id>>22) & 0x7;
 469       UInt_t posl = id&0x3fffff;
 470       AliL3SpacePointData *pointsl = fClusters[slice][patch];
 471       dx = pointsf[posf].fX - pointsl[posl].fX;
 472       dy = pointsf[posf].fY - pointsl[posl].fY;
 473
 474     }
 475
 476   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
 477   Double_t total_s ;
 478
 479   if ( fabs(localPsi) < 1. )
 480     {
 481       total_s = 2.0 * radius * asin ( localPsi ) ;
 482     }
 483   else
 484     {
 485       total_s = 2.0 * radius * Pi ;
 486     }
 487
 488   Double_t dpsi,s;
 489
 490   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 491     {
 492       UInt_t id = hitnum[i];
 493       Int_t slice = (id>>25) & 0x7f;
 494       Int_t patch = (id>>22) & 0x7;
 495       UInt_t pos = id&0x3fffff;
 496       AliL3SpacePointData *points = fClusters[slice][patch];
 497
 498       if(i>0)
 499         {
 500           id = hitnum[i-1];
 501           slice = (id>>25) & 0x7f;
 502           patch = (id>>22) & 0x7;
 503           UInt_t lastpos = id&0x3fffff;
 504           AliL3SpacePointData *lastpoints = fClusters[slice][patch];
 505           dx = points[pos].fX -lastpoints[lastpos].fX;
 506           dy = points[pos].fY -lastpoints[lastpos].fY;
 507           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
 508           fTrack->SetPsierr(dpsi);
 509           s = fS[i-1] - 2.0 * radius * (Double_t)asin ( dpsi ) ;
 510           fS[i]=s;
 511         }
 512       else
 513         fS[i]=total_s;
 514
 515       sum += 1/(points[pos].fZErr*points[pos].fZErr);
 516       ss  += 1/(points[pos].fZErr*points[pos].fZErr) * fS[i];
 517       sz  += 1/(points[pos].fZErr*points[pos].fZErr)*points[pos].fZ;
 518       sss += 1/(points[pos].fZErr*points[pos].fZErr)* fS[i] * fS[i];
 519       ssz += 1/(points[pos].fZErr*points[pos].fZErr) * fS[i] * points[pos].fZ;
 520
 521     }
 522
 523
 524   Double_t chi2,det = sum * sss - ss * ss;
 525   if ( fabs(det) < 1e-20)
 526     {
 527       chi2 = 99999.F ;
 528       //fTrack->SetChiSq2(chi2);
 529       return 0 ;
 530     }
 531
 532   //Compute the best fitted parameters A,B
 533   Double_t tanl,z0,dtanl,dz0;
 534
 535   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
 536   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
 537
 538   fTrack->SetTgl(tanl);
 539   fTrack->SetZ0(z0);
 540
 541   //     calculate chi-square
 542
 543   chi2 = 0.;
 544   Double_t r1 ;
 545
 546   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 547     {
 548       UInt_t id = hitnum[i];
 549       Int_t slice = (id>>25) & 0x7f;
 550       Int_t patch = (id>>22) & 0x7;
 551       UInt_t pos = id&0x3fffff;
 552       AliL3SpacePointData *points = fClusters[slice][patch];
 553       r1   = points[pos].fZ - tanl * fS[i] - z0 ;
 554       chi2 += (Double_t) ( (Double_t)(1/(points[pos].fZErr*points[pos].fZErr)) * (r1 * r1) );
 555     }
 556
 557   //fTrack->SetChiSq2(chi2);
 558   //
 559   //     calculate estimated variance
 560   //      varsq=chi/(double(n)-2.)
 561   //     calculate covariance matrix
 562   //      siga=sqrt(varsq*sxx/det)
 563   //      sigb=sqrt(varsq*sum/det)
 564   //
 565   dtanl = (Double_t) ( sum / det );
 566   dz0   = (Double_t) ( sss / det );
 567
 568   fTrack->SetTglerr(dtanl);
 569   fTrack->SetZ0err(dz0);
 570
 571   return 0 ;
 572 }