ITS/AliITSRiemannFit.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
  18 //
  19 //                                                                        *
  20 // This class performs a fast fit of helices going through the <=6        *
  21 // points of the ITS, with the goal of studying tracking and              *
  22 // vertexing performances.                                                *
  23 // Generated kinematics is used to take into account different weights    *
  24 // associated to points in different layers (with different multiple      *
  25 // scattering-originated errors).                                         *
  26 //                                                                        *
  27 //   Based on the work by A. Strandlie, R. Fruhwirth                      *
  28 //                                                                        *
  29 //   First implementation by N. Bustreo, R. Turrisi - July 2000           *
  30 //                                                                        *
  31 //   Further modifications by A. Dainese, R. Turrisi                      *
  32 //                                                                        *
  33 //   Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it                 *
  34 //                                                                        *
  35 // ************************************************************************
  36 //
  37 //
  38 //       Modified November, 7th 2001 by Rosario Turrisi
  39 //       (rosario.turrisi@pd.infn.it)
  40 //
  41 //       FitHelix returns different values. 0=ok, >0 =problem
  42 //       void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
  43 //
  44 //
  45 //       Modified July, 30th 2001 by Rosario Turrisi
  46 //       (rosario.turrisi@pd.infn.it)
  47 //
  48 //        Fit for z now in (z,s) plane.
  49 //        Returns parameters in order to write the helix equation
  50 //        and find the right phase/initial point.
  51 //
  52 //     "PROPER WEIGHTS":  (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
  53 //
  54
  55
  56
  57 #include "AliITSRiemannFit.h"
  58 #include "AliRun.h"
  59 #include "TClonesArray.h"
  60 #include "stdio.h"
  61 #include "stdlib.h"
  62 #include "Riostream.h"
  63 #include "TF1.h"
  64 #include "TGraphErrors.h"
  65 #include "TParticle.h"
  66 #include "TTree.h"
  67 #include "TVector3.h"
  68 #include "AliITSRecPoint.h"
  69 #include "AliITSgeom.h"
  70 #include "AliMC.h"
  71 #include "AliITSDetTypeRec.h"
  72
  73 ClassImp(AliITSRiemannFit)
  74
  75
  76 AliITSRiemannFit::AliITSRiemannFit() {
  77   ///////////////////////////////////////////////////////////
  78   // Default constructor.
  79   // Set everything to zero.
  80   ////////////////////////////////////////////////////////////
  81
  82   fSizeEvent     = 0;
  83   fPoints        = 0;
  84   fPrimaryTracks = 0;
  85   fPointRecs     = 0;
  86   //
  87   //  test erase
  88 //    fspdi          = 0;
  89 //    fspdo          = 0;
  90   for(Int_t i=0;i<6;i++)fPLay[i] = 0;
  91
  92 }
  93
  94 //______________________________________________________________________
  95 AliITSRiemannFit::AliITSRiemannFit(const AliITSRiemannFit &rf) : TObject(rf) {
  96   // Copy constructor
  97   // Copies are not allowed. The method is protected to avoid misuse.
  98   Error("AliITSRiemannFit","Copy constructor not allowed\n");
  99 }
 100
 101 //______________________________________________________________________
 102 AliITSRiemannFit& AliITSRiemannFit::operator=(const AliITSRiemannFit& /* rf */){
 103   // Assignment operator
 104   // Assignment is not allowed. The method is protected to avoid misuse.
 105   Error("= operator","Assignment operator not allowed\n");
 106   return *this;
 107 }
 108
 109 //______________________________________________________________________
 110 AliITSRiemannFit::~AliITSRiemannFit() {
 111   ///////////////////////////////////////////////////////////
 112   // Default destructor.
 113   // if arrays exist delete them. Then set everything to zero.
 114   ////////////////////////////////////////////////////////////
 115    if(fPointRecs!=0){
 116       for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
 117       delete[] fPointRecs;
 118    } // end if fPointRecs!=0
 119   fSizeEvent     = 0;
 120   fPointRecs     = 0;
 121   fPoints        = 0;
 122   fPrimaryTracks = 0;
 123   //
 124   // test erase
 125 //    fspdi          = 0;
 126 //    fspdo          = 0;
 127   for(Int_t i=0;i<6;i++)fPLay[i] = 0;
 128   return;
 129 }
 130 //----------------------------------------------------------------------
 131
 132 AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) {
 133   ///////////////////////////////////////////////////////////
 134   // Constructor.
 135   // Set fSizeEvent to size and fPrimaryTracks to ntracks.
 136   // Others to zero.
 137   ////////////////////////////////////////////////////////////
 138
 139   fSizeEvent     = size;
 140   fPoints        = 0;
 141   fPrimaryTracks = ntracks;
 142   //
 143   // test erase
 144 //    fspdi          = 0;
 145 //    fspdo          = 0;
 146   AliPointtl *first = new AliPointtl[fSizeEvent];
 147   AliPointtl **pointRecs = new AliPointtl*[fSizeEvent];
 148   for(Int_t i=0;i<6;i++)fPLay[i] = 0;
 149   for(Int_t j=0;j<fSizeEvent;j++)  // create an array of struct
 150     pointRecs[j] = &(first[j]);
 151 }
 152
 153 // ---------------------------------------------------------------------
 154 AliITSRiemannFit::AliPointtl::AliPointtl(){
 155   // default constructor
 156   SetLay();
 157   SetLad();
 158   SetDet();
 159   SetTrack();
 160   SetX();
 161   SetY();
 162   SetZ();
 163   SetR();
 164   SetdE();
 165   SetdX();
 166   SetdY();
 167   SetdZ();
 168   SetOrigin();
 169   SetMomentum();
 170   SetCode();
 171   SetName();
 172   SetPt();
 173   SetPhi();
 174   SetEta();
 175   SetVertexPhi();
 176 }
 177
 178 // ---------------------------------------------------------------------
 179
 180 void FillPoints(AliITSRiemannFit::AliPointtl **Points,Int_t &index,Float_t *xpoint,
 181                 Float_t *error,
 182                 TLorentzVector pE,TLorentzVector oT,Int_t *id,
 183                 Int_t track, Char_t *name,Int_t code,
 184                 Float_t phiorigin){
 185   ///////////////////////////////////////////////////////////////////////
 186   // Fill the structure AliPointtl with the proper data
 187   //
 188   //////////////////////////////////////////////////////////////////////
 189   Float_t pPI2 = 2.0*TMath::Pi();
 190   Float_t phi,r,x,y,z;
 191   Int_t i;
 192   i = index;
 193   x = xpoint[0];
 194   y = xpoint[1];
 195   z = xpoint[2];
 196   r = sqrt(x*x+y*y);
 197   phi = TMath::ATan2(y,x);
 198   if(phi<0.0) phi += pPI2;
 199   Points[i]->SetPhi(phi);
 200   Points[i]->SetEta(-0.5*tan(0.5*TMath::ATan2(r,z)));
 201   Points[i]->SetX(x);
 202   Points[i]->SetY(y);
 203   Points[i]->SetZ(z);
 204   Points[i]->SetdX(error[0]);
 205   Points[i]->SetdY(error[1]);
 206   Points[i]->SetdZ(error[2]);
 207   Points[i]->SetR(r);
 208   Points[i]->SetTrack(track);
 209   Points[i]->SetLay(id[0]);
 210   Points[i]->SetLad(id[1]);
 211   Points[i]->SetDet(id[2]);
 212   Points[i]->SetMomentum(&pE);
 213   Points[i]->SetOrigin(&oT);
 214   Points[i]->SetPt(sqrt(pE.X()*pE.X()+pE.Y()*pE.Y()));
 215   Points[i]->SetCode(code);
 216   Points[i]->SetName(name);
 217   Points[i]->SetVertexPhi(phiorigin);
 218   index++;
 219   return;
 220
 221 }
 222 // -----------------------------------------------------------------------
 223
 224 void AliITSRiemannFit::InitPoints(Int_t ntracks,TTree *TR,Int_t nparticles){
 225   //////////////////////////////////////////////////////////////////////
 226   // Fill the class member fPointRecs with the reconstructed points
 227   // Set All other members to the real values
 228   //
 229   /////////////////////////////////////////////////////////////////////
 230   printf("\n ************* Starting Init Points *************\n");
 231   TParticle *part;
 232
 233   AliRunLoader* rl = AliRunLoader::Open("galice.root");
 234   rl->CdGAFile();
 235   AliITSgeom* gm = (AliITSgeom*)gDirectory->Get("AliITSgeom");
 236
 237   //get pointer to modules array
 238   Int_t nmodules = gm->GetIndexMax();
 239   // Get the points from points file
 240   Int_t mod,irec;
 241   Stat_t nent;
 242   AliITSRecPoint  *recp;
 243   nent=TR->GetEntries();
 244   AliITSDetTypeRec detTypeRec;
 245   TClonesArray *iTSrec  = detTypeRec.RecPoints();
 246   Int_t totRP=0;
 247   for (mod=0; mod<nmodules; mod++) {
 248     detTypeRec.ResetRecPoints();
 249     TR->GetEvent(mod);
 250     Int_t nrecp = iTSrec->GetEntries();
 251     if(!nrecp) continue;
 252     totRP += nrecp;
 253   }
 254
 255   Int_t iMAX = totRP;
 256   fPrimaryTracks = ntracks;
 257   fParticles     = nparticles;
 258   AliITSRiemannFit::AliPointtl *global = new AliPointtl[iMAX];
 259   fPointRecs = new AliITSRiemannFit::AliPointtl*[iMAX];
 260   //
 261   for(Int_t j=0;j<iMAX;j++) {
 262     fPointRecs[j] = &(global[j]);
 263   }
 264
 265   Int_t ieta=0,ieta2=0;
 266   Int_t i,id[4],idold[4];
 267   Int_t track=0;//         // track  of hit
 268   Float_t xpoint[3],errorPlus[3],errorMinus[3],globalError[3];       // position and error of the point
 269   TLorentzVector oT,pE;
 270   Float_t locals[3],localserror[3],localsplus[3],localsminus[3]; // local position and local errors
 271   Float_t pPhi;
 272   Int_t code;
 273   const char *name;
 274   Int_t layer,ladder,detector;
 275   Float_t xcluster,zcluster;
 276   Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
 277
 278   for (mod=0; mod<nmodules; mod++) {
 279     //itsModule=(AliITSmodule*)iTSmodules->At(mod);
 280     //ITS->ResetRecPoints();
 281     detTypeRec.ResetRecPoints();
 282     TR->GetEvent(mod);
 283     Int_t nrecp = iTSrec->GetEntries();
 284     if (!nrecp) continue;
 285     //itsModule->GetID(layer,ladder,detector);
 286     gm->GetModuleId(mod,layer,ladder,detector);
 287
 288     for (irec=0;irec<nrecp;irec++) {
 289       recp   = (AliITSRecPoint*)iTSrec->UncheckedAt(irec);
 290       track=recp->fTracks[0];
 291       if(track <0 ) continue;
 292       xcluster=recp->GetX();     // x on cluster
 293       zcluster=recp->GetZ();     // z on cluster
 294       part   = (TParticle*) gAlice->GetMCApp()->Particle(track);
 295       part->ProductionVertex(oT);  // set the vertex
 296       part->Momentum(pE);          // set the vertex momentum
 297       name      = part->GetName();
 298       Char_t nam2[50];
 299       sprintf(nam2,"%s",name);
 300       code      = part->GetPdgCode();
 301       pPhi       = part->Phi();
 302       id[0]=layer;
 303       id[1]=ladder;
 304       id[2]=detector;
 305       id[3]=irec;
 306       locals[0]=xcluster;     // x on cluster
 307       locals[1]=0.0;          // y on cluster
 308       locals[2]=zcluster;     // z on cluster
 309       localserror[0]=sqrt(recp->GetSigmaX2());
 310       localserror[1]=0.0;
 311       localserror[2]=sqrt(recp->GetSigmaZ2());
 312       localsplus[0]=xcluster+sqrt(recp->GetSigmaX2());       // x on cluster
 313       if(layer==1||layer==2) localsplus[1]=0.0150/2;         // y on cluster
 314       else if(layer==3||layer==4) localsplus[1]=0.0280/2;    // y on cluster
 315       else if(layer==5||layer==6) localsplus[1]=0.0300/2;    // y on cluster
 316       localsplus[2]=zcluster+sqrt(recp->GetSigmaZ2());       // z on cluster
 317       localsminus[0]=xcluster-sqrt(recp->GetSigmaX2());      // x on cluster
 318       localsminus[1]=0.0;                                    // y on cluster
 319       localsminus[2]=zcluster-sqrt(recp->GetSigmaZ2());      // z on cluster
 320
 321       gm->LtoG(layer,ladder,detector,locals,xpoint);
 322       gm->LtoG(layer,ladder,detector,localsplus,errorPlus);
 323       gm->LtoG(layer,ladder,detector,localsminus,errorMinus);
 324       globalError[0]=0.5*TMath::Abs(errorPlus[0]-errorMinus[0]);
 325       globalError[1]=0.5*TMath::Abs(errorPlus[1]-errorMinus[1]);
 326       globalError[2]=0.5*TMath::Abs(errorPlus[2]-errorMinus[2]);
 327       if(track<ntracks) {
 328         if(TMath::Abs(part->Eta())<=1.0) ieta++;
 329         if(TMath::Abs(part->Eta())<=0.5) ieta2++;
 330       }
 331       if(!(id[0]==idold[0]&&id[1]==idold[1]&&
 332            id[2]==idold[2]&&id[3]==idold[3])) {
 333         FillPoints(fPointRecs,num,xpoint,globalError,pE,oT,id,track,nam2,code,pPhi);
 334         //
 335         // test erase
 336         switch (idold[0]) {
 337         case 1:
 338           nspdi++;
 339           break;
 340         case 2:
 341           nspdo++;
 342           break;
 343         case 3:
 344           nsddi++;
 345           break;
 346         case 4:
 347           nsddo++;
 348           break;
 349         case 5:
 350           nssdi++;
 351           break;
 352         case 6:
 353           nssdo++;
 354           break;
 355         }
 356 //      if(idold[0]==1){
 357 //        FillPoints(fspdi,nspdi,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
 358 //      }
 359 //      if(idold[0]==2){
 360
 361 //        FillPoints(fspdo,nspdo,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
 362 //      }
 363 //      if(idold[0]==3){
 364 //        nsddi++;
 365 //      }
 366 //      if(idold[0]==4){
 367 //        nsddo++;
 368 //      }
 369 //      if(idold[0]==5){
 370 //        nssdi++;
 371 //      }
 372 //      if(idold[0]==6){
 373 //        nssdo++;
 374 //      }
 375         for(i=0;i<4;i++) idold[i] = id[i];
 376         for(i=0;i<3;i++) xpoint[i]    = 0.0;
 377       } // end if id != idold
 378     } // end for irec
 379   }// end for mod
 380
 381   fPoints = num;
 382   fSizeEvent = num;
 383   fPLay[0] = nspdi ;
 384   fPLay[1] = nspdo ;
 385   fPLay[2] = nsddi ;
 386   fPLay[3] = nsddo ;
 387   fPLay[4] = nssdi ;
 388   fPLay[5] = nssdo ;
 389
 390   delete rl;
 391   printf("%d primary tracks in eta=+-1\n",ieta);
 392   printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
 393   printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
 394   printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
 395   printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
 396   printf("Points on all Layers: %d\n",num);
 397   printf("\n ************* Init Points Finished *************\n");
 398   return;
 399 }
 400 // ------------------------------------------------------------------------
 401 ///////////////////////////////////////////////////////////
 402 // Functions for sorting the fPointRecs array
 403 ///////////////////////////////////////////////////////////
 404 Bool_t SortZ(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
 405   // Z sorting function for qsort.
 406    Float_t a;
 407
 408    a = s1->GetZ() - s2->GetZ();
 409    if(a<0.0) return kTRUE;
 410    if(a>0.0) return kFALSE;
 411    return kFALSE;
 412 }
 413 Bool_t SortTrack(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
 414   // track sorting function for qsort.
 415    Float_t a;
 416
 417    a = s1->GetTrack() - s2->GetTrack();
 418    if(a<0.0) return kTRUE;
 419    if(a>0.0) return kFALSE;
 420    return kFALSE;
 421 }
 422 void hpsortTrack(AliITSRiemannFit::AliPointtl **ra,Int_t n){
 423    Int_t i,ir,j,l;
 424    AliITSRiemannFit::AliPointtl *rra;
 425
 426    if(n<2) return;
 427
 428    l  = ((n-1) >> 1) +1; // divide 2 + 1
 429    ir = n-1;
 430    for(;;){
 431      if(l>0){
 432         rra = ra[--l];  // decrement first
 433      }else{
 434         rra    = ra[ir];
 435         ra[ir] = ra[0];
 436         if(--ir == 0){  // decrement first
 437            ra[0] = rra;
 438            break;
 439         } // if --ra == 0
 440      } // end l>0
 441      i = l;
 442      j = l+1;
 443      while(j<=ir){
 444         if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
 445         if( SortTrack(rra,ra[j]) ){
 446            ra[i] = ra[j];
 447            i = j;
 448            j <<= 1; // time 2.
 449         }else{
 450            break;
 451         } // end if func()
 452      } // end while
 453      ra[i] = rra;
 454    } // end for ever
 455 }
 456 void hpsortZ(AliITSRiemannFit::AliPointtl **ra,Int_t n){
 457    Int_t i,ir,j,l;
 458    AliITSRiemannFit::AliPointtl *rra;
 459
 460    if(n<2) return;
 461
 462    l  = ((n-1) >> 1) +1; // devide 2 + 1
 463    ir = n-1;
 464    for(;;){
 465      if(l>0){
 466         rra = ra[--l];  // decrament first
 467      }else{
 468         rra    = ra[ir];
 469         ra[ir] = ra[0];
 470         if(--ir == 0){  // decrament first
 471            ra[0] = rra;
 472            break;
 473         } // if --ra == 0
 474      } // end l>0
 475      i = l;
 476      j = l+1;
 477      while(j<=ir){
 478         if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
 479         if( SortZ(rra,ra[j]) ){
 480            ra[i] = ra[j];
 481            i = j;
 482            j <<= 1; // time 2.
 483         }else{
 484            break;
 485         } // end if func()
 486      } // end while
 487      ra[i] = rra;
 488    } // end for ever
 489 }
 490 //-----------------------------------------------------------------------
 491 ////////////////////////////////////////////////////////////////////
 492 //      Sorting functions
 493 ///////////////////////////////////////////////////////////////////
 494 Int_t Partition(Int_t array[],Int_t left,Int_t right){
 495   Int_t val = array[left];
 496   Int_t lm = left - 1;
 497   Int_t rm = right + 1;
 498   for(;;) {
 499     do
 500       rm--;
 501     while
 502       (array[rm]>val);
 503     do
 504       lm++;
 505     while
 506       (array[lm]<val);
 507     if(lm<rm){
 508       Int_t tempr = array[rm];
 509       array[rm]=array[lm];
 510       array[lm]=tempr;
 511     }
 512     else
 513       return rm;
 514   }
 515
 516   return 1;
 517 }
 518
 519 ///////////////////////////////////////////////////////////////////////
 520
 521 void AliITSRiemannFit::WritePoints(void) {
 522   /////////////////////////////////////////////////////////////////////
 523   // write the data in a file (temporary ascii)
 524   /////////////////////////////////////////////////////////////////////
 525   FILE *ascii= fopen("AsciiPoints.dat","w");
 526   for(Int_t i=0;i<fPoints;i++) {
 527     fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->GetLay(),
 528             fPointRecs[i]->GetTrack(),fPointRecs[i]->GetX(),
 529             fPointRecs[i]->GetY(),fPointRecs[i]->GetZ());
 530   }
 531   fclose(ascii);
 532   return;
 533 }
 534 //-----------------------------------------------------------------------
 535
 536 void AliITSRiemannFit::ReadPoints(void) {
 537   //////////////////////////////////////////////////////////////////////
 538   // read the filled array
 539   /////////////////////////////////////////////////////////////////////
 540   hpsortTrack(fPointRecs,fPoints);
 541   for(Int_t i=0;i<fPoints;i++)
 542     printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
 543            i,fPointRecs[i]->GetLay(),fPointRecs[i]->GetTrack(),
 544            fPointRecs[i]->GetX(),fPointRecs[i]->GetY(),
 545            fPointRecs[i]->GetZ(),fPointRecs[i]->GetOrigin()->X(),
 546            fPointRecs[i]->GetOrigin()->Y(),fPointRecs[i]->GetOrigin()->Z(),
 547            fPointRecs[i]->GetPt(),fPointRecs[i]->GetName());
 548   return;
 549 }
 550 //-----------------------------------------------------------------------
 551
 552 Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
 553                                     Double_t &x1,Double_t &x2,Double_t &x3){
 554    //////////////////////////////////////////////
 555    ///  Solve cubic equation:
 556    ///  x^3 + a*x^2 +b*x + c
 557    ///
 558    ///  returns  x1 , x2 , x3
 559    ////////////////////////////////////////
 560
 561    Double_t qQ = ((a*a - 3*b)/9);
 562    Double_t rR = ((2*a*a*a - 9*a*b +27*c)/54);
 563    Double_t theta;
 564    Double_t aF = -2*sqrt(qQ);
 565    Double_t g = a/3;
 566    Double_t pPI2 = TMath::Pi()*2;
 567
 568   if( rR*rR>qQ*qQ*qQ ) {
 569     cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
 570     x1 = 9999999;
 571     x2 = 9999999;
 572     x3 = 9999999;
 573     return 0;
 574   }
 575
 576   theta = TMath::ACos(rR/sqrt(qQ*qQ*qQ));
 577
 578   x1 = (aF*TMath::Cos(theta/3))-g;
 579   x2 = (aF*TMath::Cos((theta+pPI2)/3))-g;
 580   x3 = (aF*TMath::Cos((theta-pPI2)/3))-g;
 581
 582   return 1;
 583 }
 584 //-----------------------------------------------------------------
 585
 586 void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
 587   ///////////////////////////////////////////////////////////////////////
 588   //   This function apllies the transformation in the Riemann sphere
 589   //   for xy plane
 590   ///////////////////////////////////////////////////////////////////////
 591   Float_t *rR = new Float_t[npoints];
 592   Float_t *theta = new Float_t[npoints];
 593   Float_t pPI2 = 2*TMath::Pi();
 594   Float_t x=0,y=0,z=0;
 595
 596   for(Int_t i=0;i<npoints;i++) {
 597     rR[i]     = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
 598     theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
 599     if(theta[i]<0) theta[i]+=pPI2;
 600     x     = rR[i]*cos(theta[i])/(1+rR[i]*rR[i]);
 601     y     = rR[i]*sin(theta[i])/(1+rR[i]*rR[i]);
 602     z     = rR[i]*rR[i]/(1+rR[i]*rR[i]);
 603     To[i]->SetXYZ(x,y,z);
 604   }
 605   delete[] rR;
 606   delete[] theta;
 607   return;
 608 }
 609
 610
 611 //---------------------------------------------------------------------
 612
 613 Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega,
 614                 Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError,
 615                 Double_t &corrLin){
 616   ///////////////////////////////////////////////////////////////////////
 617   //   Fit the points in the (z,s) plane - helix 3rd equation
 618   //
 619   ///////////////////////////////////////////////////////////////////////
 620   Int_t direction=0;
 621   //PH  Double_t z[npoints],x[npoints],y[npoints],s[npoints];
 622   //PH  Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
 623   Double_t * z = new Double_t[npoints];
 624   Double_t * x = new Double_t[npoints];
 625   Double_t * y = new Double_t[npoints];
 626   Double_t * s = new Double_t[npoints];
 627   Double_t * ez = new Double_t[npoints];
 628   Double_t * ex = new Double_t[npoints];
 629   Double_t * ey = new Double_t[npoints];
 630   Double_t * es = new Double_t[npoints];
 631   Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
 632
 633   //  Double_t chi=TMath::Pi()/2.0+phi;
 634   Double_t chi=-TMath::Pi()-phi;
 635   Double_t angold=0.0, tpang=0.0;
 636   for(Int_t k = 0; k<npoints; k++) {
 637     x[k] = 10.0*input[k]->X();   ex[k] = 10.0*errors[k]->X();
 638     y[k] = 10.0*input[k]->Y();   ey[k] = 10.0*errors[k]->Y();
 639     z[k] = 10.0*input[k]->Z();   ez[k] = 10.0*errors[k]->Z();
 640     if(TMath::Abs(x[k]-thu0)<1.0e-5) {  // should never happen, nor give troubles...
 641       chisquare=9999.99;
 642       cerr<<"limit for  x-x_0 "<<x[k]<<" "<<thu0<<endl;
 643       delete [] z;
 644       delete [] x;
 645       delete [] y;
 646       delete [] s;
 647       delete [] ez;
 648       delete [] ex;
 649       delete [] ey;
 650       delete [] es;
 651       return 12;
 652     }
 653     Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
 654     if( (x[k]-thu0)<0 ) {
 655       if (ang1*angold<0) {
 656         tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
 657         ang1=tpang;
 658       }
 659     }
 660     angold=ang1;
 661     if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
 662     s[k] = (ang1+chi)/omega;
 663     es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
 664   }
 665   if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
 666
 667   TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
 668   fitHist->Fit("pol1","qQ");
 669   z0  = fitHist->GetFunction("pol1")->GetParameter(0);
 670   vpar  = fitHist->GetFunction("pol1")->GetParameter(1);
 671   ez0 = fitHist->GetFunction("pol1")->GetParError(0);
 672   evpar = fitHist->GetFunction("pol1")->GetParError(1);
 673   chisquare = fitHist->GetFunction("pol1")->GetChisquare();
 674   zData.Set(z0,vpar);
 675   zError.SetXYZ(ez0,evpar,chisquare);
 676
 677   Double_t sigmas=0.;
 678   Double_t sigmaz=0.;
 679   Double_t avs=0.;
 680   Double_t avz=0.;
 681   Double_t avsz=0.;
 682
 683   for(Int_t j = 0; j < npoints; j++) {
 684     avs  += s[j];
 685     avz  += z[j];
 686     avsz += s[j]*z[j];
 687   }
 688     avs  /= (Double_t)npoints;
 689     avz  /= (Double_t)npoints;
 690     avsz /= (Double_t)npoints;
 691
 692   for(Int_t l = 0; l < npoints; l++) {
 693     sigmas += (s[l]-avs)*(s[l]-avs);
 694     sigmaz += (z[l]-avz)*(z[l]-avz);
 695   }
 696   sigmas /=(Double_t)npoints;
 697   sigmaz /=(Double_t)npoints;
 698
 699   sigmas = sqrt(sigmas);
 700   sigmaz = sqrt(sigmaz);
 701
 702   corrLin = (avsz-avs*avz)/(sigmas*sigmaz);
 703
 704   delete [] z;
 705   delete [] x;
 706   delete [] y;
 707   delete [] s;
 708   delete [] ez;
 709   delete [] ex;
 710   delete [] ey;
 711   delete [] es;
 712
 713   return 0;
 714 }
 715
 716 //-------------------------------------------------------------------
 717 Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
 718                                    Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
 719                                    Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
 720                                    Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
 721   ///////////////////////////////////////////////////////////////////////
 722   //  This function finds the helix paramenters
 723   //  d0  = impact parameter
 724   //  rho = radius of circle
 725   //  phi = atan(y0/x0)
 726   //  for the xy plane
 727   //  starting from the momentum and the outcome of
 728   //  the fit on the Riemann sphere (i.e. u0,v0,rho)
 729   //
 730   //   MIND !!!!   Here we assume both angular velocities be 1.0  (yes, one-dot-zero !)
 731   //
 732   //
 733   ///////////////////////////////////////////////////////////////////////
 734   //
 735   //  All this stuff relies on this hypothesis !!!
 736   //
 737 //   FILE *pout=fopen("chisql.dat","a");
 738   Int_t ierr = 0, ierrl=0;
 739   omega = 1.0e-2;
 740
 741   Int_t bitlay[6]={1,1,1,1,1,1};
 742   bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
 743   fd0 = -9999;   // No phisycs value
 744   u0 = -9999.9999; // parameters of helix - strange value...
 745   v0 = -9999.9999; // parameters of helix - strange value...
 746   rho = -9999.9999; // parameters of helix -unphysical strange value...
 747   Int_t pLayer = 0;
 748   const Char_t* name = 0;
 749   Int_t i=0,k=0;
 750   Int_t iMAX = 50;
 751   Int_t nN = 0;
 752   Int_t npl[6]={0,0,0,0,0,0};
 753   Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
 754   Double_t pt = sqrt(Px*Px+Py*Py);
 755   TVector3 zError;
 756   TVector2 zData;
 757   Double_t corrLin;
 758   TVector3 *ori        = new TVector3[iMAX];
 759   TVector3 **original  = new TVector3*[iMAX];
 760   TVector3 *rie       = new TVector3[iMAX];
 761   TVector3 **riemann  = new TVector3*[iMAX];
 762   TVector3 *err        = new TVector3[iMAX];
 763   TVector3 **errors    = new TVector3*[iMAX];
 764   TVector3 *linerr        = new TVector3[iMAX];
 765   TVector3 **linerrors    = new TVector3*[iMAX];
 766   //PH  Double_t weight[iMAX];
 767   Double_t * weight = new Double_t[iMAX];
 768
 769   for(i=0;i<iMAX;i++){
 770     original[i]   = &(ori[i]);
 771     riemann[i]   = &(rie[i]);
 772     errors[i]     = &(err[i]);
 773     linerrors[i]  = &(linerr[i]);
 774   }
 775   for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
 776   Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
 777   a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
 778   Double_t xbar = 0;
 779   Double_t ybar = 0;
 780   Double_t zbar = 0;
 781   Double_t a,b,c,d;       // cubic parameters
 782   Double_t roots[3]= {0.0,0.0,0.0};  // cubic solutions
 783   Double_t value = 0.0;   // minimum eigenvalue
 784   Double_t x1,x2,x3;      // eigenvector component
 785   Double_t n1,n2,n3,nr= 0;// unit eigenvector
 786   Double_t radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
 787   Double_t sigmaMS = 0;
 788   TVector3 vVec,vVecNor;
 789
 790 // Select RecPoints belonging to the track
 791   for(k =0;k<fPoints;k++){
 792     if(fPointRecs[k]->GetTrack()==tracknumber) {
 793       name = fPointRecs[k]->GetName();
 794       pt = fPointRecs[k]->GetPt();
 795       pLayer = fPointRecs[k]->GetLay();
 796       Int_t ilay = pLayer-1;
 797       if(npl[ilay]!=0) continue;
 798       if(bitlay[ilay] == 1) {
 799         original[nN]->SetXYZ(0.1*fPointRecs[k]->GetX(),0.1*fPointRecs[k]->GetY(),0.1*fPointRecs[k]->GetZ());
 800         errors[nN]->SetXYZ(0.1*fPointRecs[k]->GetdX(),0.1*fPointRecs[k]->GetdY(),0.1*fPointRecs[k]->GetdZ());
 801         sigmaMS = (radiusdm[pLayer]-radiusdm[0])*0.000724/pP;// beam pipe contribution
 802         for(Int_t j=1;j<pLayer;j++) {
 803           sigmaMS += (radiusdm[pLayer]-radiusdm[j])*0.00136/pP;
 804         }
 805         weight[nN] = ( 1 + original[nN]->Perp2() )*( 1+ original[nN]->Perp2() )/
 806           ( errors[nN]->Perp2() + sigmaMS*sigmaMS );
 807         linerrors[nN]->SetXYZ(errors[nN]->X(),errors[nN]->Y(),sqrt(errors[nN]->Z()*errors[nN]->Z()+sigmaMS*sigmaMS));
 808         nN++;
 809         npl[ilay]++;
 810       }  // end if on layer
 811     }    //end if track==tracknumber
 812   }      //end for k
 813   //
 814   //    6 points, no more, no less
 815   //
 816   if(original[5]->X() == 9999 || original[6]->X() != 9999)
 817     {
 818       delete [] weight;
 819       return 1;   // not enough points
 820     }
 821
 822   //
 823   //
 824   //
 825   //  FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
 826   //
 827
 828   RiemannTransf(nN,original,riemann);
 829
 830   Double_t sumWeights = 0.0; // sum of weights factor
 831
 832   for(Int_t j=0;j<nN;j++){ // mean values for x[i],y[i],z[i]
 833     xbar+=weight[j]*riemann[j]->X();
 834     ybar+=weight[j]*riemann[j]->Y();
 835     zbar+=weight[j]*riemann[j]->Z();
 836     sumWeights+=weight[j];
 837   }
 838
 839   xbar /= sumWeights;
 840   ybar /= sumWeights;
 841   zbar /= sumWeights;
 842
 843   for(Int_t j=0;j<nN;j++) {  // Calculate the matrix elements
 844     a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
 845     a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
 846     a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
 847     a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
 848     a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
 849     a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
 850   }
 851
 852   a11 /= nN;
 853   a12 /= nN;
 854   a22 /= nN;
 855   a23 /= nN;
 856   a13 /= nN;
 857   a33 /= nN;
 858   a21 = a12;
 859   a32 = a23;
 860   a31 = a13;
 861
 862 // **************  Determinant  parameters ********************
 863 //   n.b. simplifications done keeping in mind symmetry of A
 864 //
 865   a = 1;
 866   b = (-a11-a33-a22);
 867   c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
 868   d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
 869
 870 // **************  Find the 3 eigenvalues *************************
 871   Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
 872
 873   if(checkCubic !=1 ){
 874     printf("Track %d Has no real solution continuing ...\n",tracknumber);
 875     delete [] weight;
 876     return 2;
 877   }
 878
 879 // **************** Find the lowest eigenvalue *****************
 880   if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
 881   if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
 882   if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
 883
 884   // ************ Eigenvector relative to value **************
 885   x3 = 1;
 886   x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
 887   x1 = (value-a33-a32*x2)/a31;
 888   vVec.SetXYZ(x1,x2,x3);
 889   vVecNor = vVec.Unit();
 890   n1 = vVecNor.X();
 891   n2 = vVecNor.Y();
 892   n3 = vVecNor.Z();
 893   nr = -n1*xbar-n2*ybar-n3*zbar;
 894
 895   u0  = -0.5*n1/(nr+n3);
 896   v0  = -0.5*n2/(nr+n3);
 897   rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
 898
 899   fFit = 0.0;
 900   fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
 901   fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
 902   fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
 903   fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
 904   fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
 905   fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->Y()-v0))-rho);
 906
 907   fd0      =   100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
 908   fphi     =   TMath::ATan2(v0,u0);
 909
 910 //**************************************************************************
 911 //  LINEAR FIT IN (z,s) PLANE:    z = zData.X() + zData.Y()*s
 912 //  strictly linear (no approximation)
 913 //**************************************************************************
 914
 915        ////////////////////////////////////////////////////////////////////////////////////////////////////////////
 916       //                                                                                                        //
 917      //      REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S    //
 918     //       rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes...  //
 919    //                                                                                                        //
 920   ////////////////////////////////////////////////////////////////////////////////////////////////////////////
 921
 922   rho    *=  10.0;
 923   u0     *=  10.0;
 924   v0     *=  10.0;
 925   ierrl=FitLinear(nN,original,linerrors,omega,u0,v0,fphi,zData,zError,corrLin);
 926   chisql=zError.Z();
 927 //   fprintf(pout,"%f \n",chisql);
 928   z0=zData.X();
 929   vpar=zData.Y();
 930   fCorrLin = corrLin;
 931   ierr = (ierrl > ierr ? ierrl : ierr);
 932 //   fclose(pout);
 933   delete [] weight;
 934   return ierr;
 935 }
 936 Int_t AliITSRiemannFit::FitHelix(Int_t NPoints, TVector3** fPointRecs,TVector3** fPointRecErrors,Float_t& f1, Float_t& f2, Float_t& f3) {
 937
 938   ///////////////////////////////////////////////////////////////////////
 939   //  This function finds the helix parameters
 940   //  d0  = impact parameter
 941   //  rho = radius of circle
 942   //  phi = atan(y0/x0)
 943   //  for the xy plane
 944   //  starting from the momentum and the outcome of
 945   //  the fit on the Riemann sphere (i.e. u0,v0,rho)
 946   //
 947   //   MIND !!!!   Here we assume both angular velocities be 1.0e-2  (yes, 0.01 !)
 948   //
 949   //
 950   //   Also linear fit in (z,s) is performed, so it's 3-D !
 951   //   z0 and vpar are calculated (intercept and z-component of velocity, but
 952   //   in units... you guess.
 953   //
 954   //
 955   //   Values calculated in addition:
 956   //
 957   //              - transverse impact parameter        fd0
 958   //              - sum of residuals in (x,y) plane    fFit
 959   //              - chisquare of linear fit            chisql
 960   //              - correlation coefficient            fCorrLin
 961   //
 962   //
 963   //
 964   //
 965   //
 966   ///////////////////////////////////////////////////////////////////////
 967   //
 968   //  All this stuff relies on this hypothesis !!!
 969   //
 970   Int_t ierr = 0, ierrl=0;
 971   const Double_t kOmega = 1.0e-2;
 972
 973
 974
 975
 976   Double_t  fd0 = -9999;      //  fake values
 977   Double_t  u0 = -9999.9999;  //  for eventual
 978   Double_t  v0 = -9999.9999;  //  debugging
 979   Double_t  rho = -9999.9999; //
 980   Double_t fphi, fFit, chisql, z0, vpar, fCorrLin;
 981
 982   //
 983   //  This info is no more there... to be re-considered... maybe
 984   //
 985   //   Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
 986   //   Double_t pt = sqrt(Px*Px+Py*Py);
 987
 988   TVector3 zError;
 989   TVector2 zData;
 990   Double_t corrLin;
 991   TVector3 *ori        = new TVector3[NPoints];
 992   TVector3 **original  = new TVector3*[NPoints];
 993   TVector3 *rie        = new TVector3[NPoints];
 994   TVector3 **riemann   = new TVector3*[NPoints];
 995   TVector3 *err        = new TVector3[NPoints];
 996   TVector3 **errors    = new TVector3*[NPoints];
 997   TVector3 *linerr     = new TVector3[NPoints];
 998   TVector3 **linerrors = new TVector3*[NPoints];
 999   Double_t * weight = new Double_t[NPoints];
1000
1001   for(Int_t i=0; i<NPoints; i++){
1002
1003     original[i]   = &(ori[i]);
1004     riemann[i]   = &(rie[i]);
1005     errors[i]     = &(err[i]);
1006     linerrors[i]  = &(linerr[i]);
1007
1008     original[i]->SetXYZ(9999,9999,9999);
1009   }
1010
1011   //
1012   //  Riemann fit parameters
1013   //
1014   Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
1015   a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
1016   Double_t xbar = 0;
1017   Double_t ybar = 0;
1018   Double_t zbar = 0;
1019   //
1020   Double_t a,b,c,d;                  // cubic parameters
1021   Double_t roots[3]= {0.0,0.0,0.0};  // cubic solutions
1022   Double_t value = 0.0;              // minimum eigenvalue
1023   Double_t x1,x2,x3;                 // eigenvector component
1024   Double_t n1,n2,n3,nr= 0;           // unit eigenvector
1025   TVector3 vVec,vVecNor;
1026
1027   for (Int_t ip=0; ip<NPoints; ip++) {
1028 original[ip]->SetXYZ(0.1*fPointRecs[ip]->X(),0.1*fPointRecs[ip]->Y(),0.1*fPointRecs[ip]->Z());
1029
1030 errors[ip]->SetXYZ(0.1*fPointRecErrors[ip]->X(),0.1*fPointRecErrors[ip]->Y(),0.1*fPointRecErrors[ip]->Z());
1031     weight[ip] = (1+original[ip]->Perp2())*(1+original[ip]->Perp2())/(errors[ip]->Perp2());
1032     linerrors[ip]->SetXYZ(errors[ip]->X(),errors[ip]->Y(),errors[ip]->Z());
1033   }
1034
1035
1036   //
1037   //
1038   //  FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
1039   //
1040
1041   RiemannTransf(NPoints,original,riemann);
1042
1043   Double_t sumWeights = 0.0; // sum of weights factor
1044
1045   for(Int_t j=0;j<NPoints;j++){ // mean values for x[i],y[i],z[i]
1046     xbar+=weight[j]*riemann[j]->X();
1047     ybar+=weight[j]*riemann[j]->Y();
1048     zbar+=weight[j]*riemann[j]->Z();
1049     sumWeights+=weight[j];
1050   }
1051
1052   xbar /= sumWeights;
1053   ybar /= sumWeights;
1054   zbar /= sumWeights;
1055
1056   for(Int_t j=0;j<NPoints;j++) {  // Calculate the matrix elements
1057     a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
1058     a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
1059     a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
1060     a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
1061     a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
1062     a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
1063   }
1064   //
1065   // this doesn't seem to work...
1066   //
1067   //   a11 /= sumWeights;
1068   //   a12 /= sumWeights;
1069   //   a22 /= sumWeights;
1070   //   a23 /= sumWeights;
1071   //   a13 /= sumWeights;
1072   //   a33 /= sumWeights;
1073
1074   a11 /= NPoints;
1075   a12 /= NPoints;
1076   a22 /= NPoints;
1077   a23 /= NPoints;
1078   a13 /= NPoints;
1079   a33 /= NPoints;
1080   a21 = a12;
1081   a32 = a23;
1082   a31 = a13;
1083
1084 // **************  Determinant  parameters ********************
1085 //   n.b. simplifications done keeping in mind symmetry of A
1086 //
1087   a = 1;
1088   b = (-a11-a33-a22);
1089   c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
1090   d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
1091
1092 // **************  Find the 3 eigenvalues *************************
1093   Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
1094
1095   if(checkCubic !=1 ){
1096     printf("No real solution. Check data.\n");
1097     delete [] weight;
1098     return 999;
1099   }
1100
1101 // **************** Find the lowest eigenvalue *****************
1102   if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
1103   if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
1104   if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
1105
1106   // ************ Eigenvector relative to value **************
1107   x3 = 1;
1108   x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
1109   x1 = (value-a33-a32*x2)/a31;
1110   vVec.SetXYZ(x1,x2,x3);
1111   vVecNor = vVec.Unit();
1112   n1 = vVecNor.X();
1113   n2 = vVecNor.Y();
1114   n3 = vVecNor.Z();
1115   nr = -n1*xbar-n2*ybar-n3*zbar;
1116
1117   u0  = -0.5*n1/(nr+n3);
1118   v0  = -0.5*n2/(nr+n3);
1119   rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
1120
1121
1122   fFit = 0.0;
1123   for (Int_t i=0; i<NPoints; i++) {
1124   fFit += 10.*TMath::Abs(sqrt((original[i]->X()-u0)*(original[i]->X()-u0)+(original[i]->Y()-v0)*(original[i]->Y()-v0))-rho);
1125   }
1126   fd0      =   100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
1127   fphi     =   TMath::ATan2(v0,u0);
1128
1129 //**************************************************************************
1130 //  LINEAR FIT IN (z,s) PLANE:    z = zData.X() + zData.Y()*s
1131 //  strictly linear (no approximation)
1132 //**************************************************************************
1133
1134        ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1135       //                                                                                                        //
1136      //      REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S    //
1137     //       rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes...  //
1138    //                                                                                                        //
1139   ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1140
1141   rho    *=  10.0;
1142   u0     *=  10.0;
1143   v0     *=  10.0;
1144
1145   ierrl=LinearFit(NPoints,original,linerrors,kOmega,u0,v0,fphi,zData,zError,corrLin);
1146   if(ierrl==33) return 0;
1147   chisql=zError.Z();
1148 //   fprintf(pout,"%f \n",chisql);
1149   z0=zData.X();
1150   vpar=zData.Y();
1151   fCorrLin = corrLin;
1152   ierr = (ierrl > ierr ? ierrl : ierr);
1153 //   fclose(pout);
1154   delete [] weight;
1155
1156   f1=fphi;
1157   f2=vpar/(kOmega*TMath::Abs(rho));
1158   f3=1/rho;
1159   delete[] ori;
1160   delete[] rie;
1161   delete[] err;
1162   delete[] linerr;
1163   delete[] original;
1164   delete[] riemann;
1165   delete[] errors;
1166   delete[] linerrors;
1167
1168   return 1;
1169
1170
1171 }
1172
1173 //____________________________________________________________
1174
1175 Int_t AliITSRiemannFit::LinearFit(Int_t npoints, TVector3 **input,
1176                                  TVector3 **errors, Double_t omega,
1177                                  Double_t &thu0, Double_t &thv0, Double_t &phi,TVector2 &zData, TVector3 &zError,
1178                                  Double_t &corrLin){
1179   ///////////////////////////////////////////////////////////////////////
1180   //   Fit the points in the (z,s) plane - helix 3rd equation
1181   //
1182   ///////////////////////////////////////////////////////////////////////
1183   //By R.Turrisi
1184
1185   Int_t direction=0;
1186   //PH  Double_t z[npoints],x[npoints],y[npoints],s[npoints];
1187   //PH  Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
1188   Double_t * z = new Double_t[npoints];
1189   Double_t * x = new Double_t[npoints];
1190   Double_t * y = new Double_t[npoints];
1191   Double_t * s = new Double_t[npoints];
1192   Double_t * ez = new Double_t[npoints];
1193   Double_t * ex = new Double_t[npoints];
1194   Double_t * ey = new Double_t[npoints];
1195   Double_t * es = new Double_t[npoints];
1196   Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
1197
1198
1199   //  Double_t chi=TMath::Pi()/2.0+phi;
1200   Double_t chi=-TMath::Pi()-phi;
1201   Double_t angold=0.0, tpang=0.0;
1202   for(Int_t k = 0; k<npoints; k++) {
1203     x[k] = 10.0*input[k]->X();   ex[k] = 10.0*errors[k]->X();
1204     y[k] = 10.0*input[k]->Y();   ey[k] = 10.0*errors[k]->Y();
1205     z[k] = 10.0*input[k]->Z();   ez[k] = 10.0*errors[k]->Z();
1206     if(TMath::Abs(x[k]-thu0)<1.0e-5) {  // should never happen, nor give troubles...
1207       chisquare=9999.99;
1208       cerr<<"limit for  x-x_0 "<<x[k]<<" "<<thu0<<endl;
1209       delete [] z;
1210       delete [] x;
1211       delete [] y;
1212       delete [] s;
1213       delete [] ez;
1214       delete [] ex;
1215       delete [] ey;
1216       delete [] es;
1217       return 12;
1218     }
1219     Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
1220     if( (x[k]-thu0)<0 ) {
1221       if (ang1*angold<0) {
1222         tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
1223         ang1=tpang;
1224       }
1225     }
1226     angold=ang1;
1227     if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
1228     s[k] = (ang1+chi)/omega;
1229     es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
1230   }
1231   if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
1232
1233   //  if(s[0]>-636 && s[0]<-625) return 33;
1234
1235   TGraph* fitHist = new TGraph(npoints,s,z);
1236   TF1* f1 = new TF1("f1",Fitfunction,-100,100,2);
1237
1238   f1->SetParameter(0,1);
1239   f1->SetParameter(1,1);
1240   f1->SetLineColor(2);
1241   fitHist->Fit(f1,"qQ");
1242
1243   z0   = f1->GetParameter(0);
1244   vpar = f1->GetParameter(1);
1245   ez0  = f1->GetParError(0);
1246   evpar= f1->GetParError(1);
1247   chisquare=f1->GetChisquare();
1248   zData.Set(z0,vpar);
1249   zError.SetXYZ(ez0,evpar,chisquare);
1250
1251   Double_t sigmas=0.;
1252   Double_t sigmaz=0.;
1253   Double_t avs=0.;
1254   Double_t avz=0.;
1255   Double_t avsz=0.;
1256
1257   for(Int_t j = 0; j < npoints; j++) {
1258     avs  += s[j];
1259     avz  += z[j];
1260     avsz += s[j]*z[j];
1261   }
1262     avs  /= (Double_t)npoints;
1263     avz  /= (Double_t)npoints;
1264     avsz /= (Double_t)npoints;
1265
1266   for(Int_t l = 0; l < npoints; l++) {
1267     sigmas += (s[l]-avs)*(s[l]-avs);
1268     sigmaz += (z[l]-avz)*(z[l]-avz);
1269   }
1270   sigmas /=(Double_t)npoints;
1271   sigmaz /=(Double_t)npoints;
1272
1273   sigmas = sqrt(sigmas);
1274   sigmaz = sqrt(sigmaz);
1275
1276   corrLin = (avsz-avs*avz)/(sigmas*sigmaz);
1277
1278
1279
1280   delete [] z;
1281   delete [] x;
1282   delete [] y;
1283   delete [] s;
1284   delete [] ez;
1285   delete [] ex;
1286   delete [] ey;
1287   delete [] es;
1288   delete f1; delete fitHist;
1289   return 0;
1290 }
1291
1292
1293 //_______________________________________________________
1294
1295 Double_t AliITSRiemannFit::Fitfunction(Double_t *x, Double_t* par){
1296   // function used for fit
1297   return par[0]+(*x)*par[1];
1298
1299 }
1300