ITS/AliITSRiemannFit.cxx

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