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