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