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