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