-//-----------------------------------------------------------------------------
-void AliITSRiemannFit::Streamer(TBuffer &lRb){
-////////////////////////////////////////////////////////////////////////
-// The default Streamer function "written by ROOT" doesn't write out
-// the arrays referenced by pointers. Therefore, a specific Streamer function
-// has to be written. This function should not be modified but instead added
-// on to so that older versions can still be read.
-////////////////////////////////////////////////////////////////////////
- // Stream an object of class AliITSRiemannFit.
- Int_t i,j,n;
- Int_t ii,jj;
- n=20;
-
- if (lRb.IsReading()) {
- Version_t lRv = lRb.ReadVersion(); if (lRv) { }
- TObject::Streamer(lRb);
- lRb >> fSizeEvent;
- lRb >> fPrimaryTracks;
- lRb >> fPoints;
- for(i=0;i<6;i++) lRb >> fPLay[i];
- if(fPointRecs!=0){
- for(i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
- delete[] fPointRecs;
- } // end if fPointRecs!=0
- fPointRecs = new Point_tl*[fSizeEvent];
- for(i=0;i<fSizeEvent;i++){
- fPointRecs[i] = new Point_tl[n];
- for(j=0;j<n;j++){
- lRb >> fPointRecs[i][j].lay;
- lRb >> fPointRecs[i][j].lad;
- lRb >> fPointRecs[i][j].det;
- lRb >> fPointRecs[i][j].track;
- lRb >> fPointRecs[i][j].fx;
- lRb >> fPointRecs[i][j].fy;
- lRb >> fPointRecs[i][j].fz;
- lRb >> fPointRecs[i][j].fr;
- lRb >> fPointRecs[i][j].fdE;
- lRb >> fPointRecs[i][j].fdx;
- lRb >> fPointRecs[i][j].fdy;
- lRb >> fPointRecs[i][j].fdz;
- lRb >> fPointRecs[i][j].fPt;
- lRb >> (Char_t*)fPointRecs[i][j].fName;
- for (ii=0;ii<4;ii++)
- lRb << fPointRecs[i][j].fOrigin[ii];
- for (jj=0;jj<4;jj++)
- lRb << fPointRecs[i][j].fMomentum[jj];
- lRb >> fPointRecs[i][j].fCode;
- lRb >> fPointRecs[i][j].phi;
- lRb >> fPointRecs[i][j].eta;
- lRb >> fPointRecs[i][j].vertexPhi;
- } //end for j
- } //end for i
-// if(fspdi!=0){
-// for(i=0;i<fSizeEvent/6;i++) delete[] fspdi[i];
-// delete[] fspdi;
-// } // end if fspdi!=0
-// fspdi = new Point_tl*[fSizeEvent/6];
-// for(i=0;i<fSizeEvent/6;i++){
-// fspdi[i] = new Point_tl[n];
-// for(j=0;j<n;j++){
-// lRb >> fspdi[i][j].lay;
-// lRb >> fspdi[i][j].lad;
-// lRb >> fspdi[i][j].det;
-// lRb >> fspdi[i][j].track;
-// lRb >> fspdi[i][j].fx;
-// lRb >> fspdi[i][j].fy;
-// lRb >> fspdi[i][j].fz;
-// lRb >> fspdi[i][j].fr;
-// lRb >> fspdi[i][j].fdE;
-// lRb >> fspdi[i][j].fdx;
-// lRb >> fspdi[i][j].fdy;
-// lRb >> fspdi[i][j].fdz;
-// lRb >> fspdi[i][j].fPt;
-// for (ii=0;ii<4;ii++)
-// lRb << fspdi[i][j].fOrigin[ii];
-// for (jj=0;jj<4;jj++)
-// lRb << fspdi[i][j].fMomentum[jj];
-// lRb >> fspdi[i][j].fCode;
-// lRb >> (Char_t*)fspdi[i][j].fName;
-// lRb >> fspdi[i][j].phi;
-// lRb >> fspdi[i][j].eta;
-// lRb >> fspdi[i][j].vertexPhi;
-// } //end for j
-// } //end for i
-// if(fspdo!=0){
-// for(i=0;i<fSizeEvent/6;i++) delete[] fspdo[i];
-// delete[] fspdo;
-// } // end if fspdo!=0
-// fspdo = new Point_tl*[fSizeEvent/6];
-// for(i=0;i<fSizeEvent/6;i++){
-// fspdo[i] = new Point_tl[n];
-// for(j=0;j<n;j++){
-// lRb >> fspdo[i][j].lay;
-// lRb >> fspdo[i][j].lad;
-// lRb >> fspdo[i][j].det;
-// lRb >> fspdo[i][j].track;
-// lRb >> fspdo[i][j].fx;
-// lRb >> fspdo[i][j].fy;
-// lRb >> fspdo[i][j].fz;
-// lRb >> fspdo[i][j].fr;
-// lRb >> fspdo[i][j].fdE;
-// lRb >> fspdo[i][j].fdx;
-// lRb >> fspdo[i][j].fdy;
-// lRb >> fspdo[i][j].fdz;
-// lRb >> fspdo[i][j].fPt;
-// for (ii=0;ii<4;ii++)
-// lRb << fspdo[i][j].fOrigin[ii];
-// for (jj=0;jj<4;jj++)
-// lRb << fspdo[i][j].fMomentum[jj];
-// lRb >> fspdo[i][j].fCode;
-// lRb >> (Char_t*)fspdo[i][j].fName;
-// lRb >> fspdo[i][j].phi;
-// lRb >> fspdo[i][j].eta;
-// lRb >> fspdo[i][j].vertexPhi;
-// } //end for j
-// } //end for i
- } else {
- lRb.WriteVersion(AliITSRiemannFit::IsA());
- TObject::Streamer(lRb);
- lRb << fSizeEvent;
- lRb << fPrimaryTracks;
- lRb << fPoints;
- for(i=0;i<6;i++) lRb >> fPLay[i];
- for(i=0;i<fSizeEvent;i++) for(j=0;j<n;j++){
- lRb << fPointRecs[i][j].lay;
- lRb << fPointRecs[i][j].lad;
- lRb << fPointRecs[i][j].det;
- lRb << fPointRecs[i][j].track;
- lRb << fPointRecs[i][j].fx;
- lRb << fPointRecs[i][j].fy;
- lRb << fPointRecs[i][j].fz;
- lRb << fPointRecs[i][j].fr;
- lRb << fPointRecs[i][j].fdE;
- lRb << fPointRecs[i][j].fdx;
- lRb << fPointRecs[i][j].fdy;
- lRb << fPointRecs[i][j].fdz;
- lRb << fPointRecs[i][j].fPt;
- for (ii=0;ii<4;ii++)
- lRb << fPointRecs[i][j].fOrigin[ii];
- for (jj=0;jj<4;jj++)
- lRb << fPointRecs[i][j].fMomentum[jj];
- lRb << fPointRecs[i][j].fCode;
- lRb << fPointRecs[i][j].fName;
- lRb << fPointRecs[i][j].phi;
- lRb << fPointRecs[i][j].eta;
- lRb << fPointRecs[i][j].vertexPhi;
+ ///////////////////////////////////////////////////////////////////////
+ // This function finds the helix parameters
+ // d0 = impact parameter
+ // rho = radius of circle
+ // phi = atan(y0/x0)
+ // for the xy plane
+ // starting from the momentum and the outcome of
+ // the fit on the Riemann sphere (i.e. u0,v0,rho)
+ //
+ // MIND !!!! Here we assume both angular velocities be 1.0e-2 (yes, 0.01 !)
+ //
+ //
+ // Also linear fit in (z,s) is performed, so it's 3-D !
+ // z0 and vpar are calculated (intercept and z-component of velocity, but
+ // in units... you guess.
+ //
+ //
+ // Values calculated in addition:
+ //
+ // - transverse impact parameter fd0
+ // - sum of residuals in (x,y) plane fFit
+ // - chisquare of linear fit chisql
+ // - correlation coefficient fCorrLin
+ //
+ //
+ //
+ //
+ //
+ ///////////////////////////////////////////////////////////////////////
+ //
+ // All this stuff relies on this hypothesis !!!
+ //
+ Int_t ierr = 0, ierrl=0;
+ const Double_t kOmega = 1.0e-2;
+
+
+
+
+ Double_t fd0 = -9999; // fake values
+ Double_t u0 = -9999.9999; // for eventual
+ Double_t v0 = -9999.9999; // debugging
+ Double_t rho = -9999.9999; //
+ Double_t fphi, fFit, chisql, z0, vpar, fCorrLin;
+
+ //
+ // This info is no more there... to be re-considered... maybe
+ //
+ // Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
+ // Double_t pt = sqrt(Px*Px+Py*Py);
+
+ TVector3 zError;
+ TVector2 zData;
+ Double_t corrLin;
+ TVector3 *ori = new TVector3[NPoints];
+ TVector3 **original = new TVector3*[NPoints];
+ TVector3 *rie = new TVector3[NPoints];
+ TVector3 **riemann = new TVector3*[NPoints];
+ TVector3 *err = new TVector3[NPoints];
+ TVector3 **errors = new TVector3*[NPoints];
+ TVector3 *linerr = new TVector3[NPoints];
+ TVector3 **linerrors = new TVector3*[NPoints];
+ Double_t * weight = new Double_t[NPoints];
+
+ for(Int_t i=0; i<NPoints; i++){
+
+ original[i] = &(ori[i]);
+ riemann[i] = &(rie[i]);
+ errors[i] = &(err[i]);
+ linerrors[i] = &(linerr[i]);
+
+ original[i]->SetXYZ(9999,9999,9999);
+ }
+
+ //
+ // Riemann fit parameters
+ //
+ Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
+ a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
+ Double_t xbar = 0;
+ Double_t ybar = 0;
+ Double_t zbar = 0;
+ //
+ Double_t a,b,c,d; // cubic parameters
+ Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
+ Double_t value = 0.0; // minimum eigenvalue
+ Double_t x1,x2,x3; // eigenvector component
+ Double_t n1,n2,n3,nr= 0; // unit eigenvector
+ TVector3 vVec,vVecNor;
+
+ for (Int_t ip=0; ip<NPoints; ip++) {
+original[ip]->SetXYZ(0.1*fPointRecs[ip]->X(),0.1*fPointRecs[ip]->Y(),0.1*fPointRecs[ip]->Z());
+
+errors[ip]->SetXYZ(0.1*fPointRecErrors[ip]->X(),0.1*fPointRecErrors[ip]->Y(),0.1*fPointRecErrors[ip]->Z());
+ weight[ip] = (1+original[ip]->Perp2())*(1+original[ip]->Perp2())/(errors[ip]->Perp2());
+ linerrors[ip]->SetXYZ(errors[ip]->X(),errors[ip]->Y(),errors[ip]->Z());
+ }
+
+
+ //
+ //
+ // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
+ //
+
+ RiemannTransf(NPoints,original,riemann);
+
+ Double_t sumWeights = 0.0; // sum of weights factor
+
+ for(Int_t j=0;j<NPoints;j++){ // mean values for x[i],y[i],z[i]
+ xbar+=weight[j]*riemann[j]->X();
+ ybar+=weight[j]*riemann[j]->Y();
+ zbar+=weight[j]*riemann[j]->Z();
+ sumWeights+=weight[j];
+ }
+
+ xbar /= sumWeights;
+ ybar /= sumWeights;
+ zbar /= sumWeights;
+
+ for(Int_t j=0;j<NPoints;j++) { // Calculate the matrix elements
+ a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
+ a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
+ a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
+ a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
+ a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
+ a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
+ }
+ //
+ // this doesn't seem to work...
+ //
+ // a11 /= sumWeights;
+ // a12 /= sumWeights;
+ // a22 /= sumWeights;
+ // a23 /= sumWeights;
+ // a13 /= sumWeights;
+ // a33 /= sumWeights;
+
+ a11 /= NPoints;
+ a12 /= NPoints;
+ a22 /= NPoints;
+ a23 /= NPoints;
+ a13 /= NPoints;
+ a33 /= NPoints;
+ a21 = a12;
+ a32 = a23;
+ a31 = a13;
+
+// ************** Determinant parameters ********************
+// n.b. simplifications done keeping in mind symmetry of A
+//
+ a = 1;
+ b = (-a11-a33-a22);
+ c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
+ d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
+
+// ************** Find the 3 eigenvalues *************************
+ Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
+
+ if(checkCubic !=1 ){
+ printf("No real solution. Check data.\n");
+ delete [] weight;
+ return 999;
+ }
+
+// **************** Find the lowest eigenvalue *****************
+ if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
+ if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
+ if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
+
+ // ************ Eigenvector relative to value **************
+ x3 = 1;
+ x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
+ x1 = (value-a33-a32*x2)/a31;
+ vVec.SetXYZ(x1,x2,x3);
+ vVecNor = vVec.Unit();
+ n1 = vVecNor.X();
+ n2 = vVecNor.Y();
+ n3 = vVecNor.Z();
+ nr = -n1*xbar-n2*ybar-n3*zbar;
+
+ u0 = -0.5*n1/(nr+n3);
+ v0 = -0.5*n2/(nr+n3);
+ rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
+
+
+ fFit = 0.0;
+ for (Int_t i=0; i<NPoints; i++) {
+ fFit += 10.*TMath::Abs(sqrt((original[i]->X()-u0)*(original[i]->X()-u0)+(original[i]->Y()-v0)*(original[i]->Y()-v0))-rho);
+ }
+ fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
+ fphi = TMath::ATan2(v0,u0);
+
+//**************************************************************************
+// LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
+// strictly linear (no approximation)
+//**************************************************************************
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
+ // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
+ // //
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ rho *= 10.0;
+ u0 *= 10.0;
+ v0 *= 10.0;
+
+ ierrl=LinearFit(NPoints,original,linerrors,kOmega,u0,v0,fphi,zData,zError,corrLin);
+ if(ierrl==33) return 0;
+ chisql=zError.Z();
+// fprintf(pout,"%f \n",chisql);
+ z0=zData.X();
+ vpar=zData.Y();
+ fCorrLin = corrLin;
+ ierr = (ierrl > ierr ? ierrl : ierr);
+// fclose(pout);
+ delete [] weight;
+
+ f1=fphi;
+ f2=vpar/(kOmega*TMath::Abs(rho));
+ f3=1/rho;
+ delete[] ori;
+ delete[] rie;
+ delete[] err;
+ delete[] linerr;
+ delete[] original;
+ delete[] riemann;
+ delete[] errors;
+ delete[] linerrors;
+
+ return 1;
+
+
+}
+
+//____________________________________________________________
+
+Int_t AliITSRiemannFit::LinearFit(Int_t npoints, TVector3 **input,
+ TVector3 **errors, Double_t omega,
+ Double_t &thu0, Double_t &thv0, Double_t &phi,TVector2 &zData, TVector3 &zError,
+ Double_t &corrLin){
+ ///////////////////////////////////////////////////////////////////////
+ // Fit the points in the (z,s) plane - helix 3rd equation
+ //
+ ///////////////////////////////////////////////////////////////////////
+ //By R.Turrisi
+
+ Int_t direction=0;
+ //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
+ //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
+ Double_t * z = new Double_t[npoints];
+ Double_t * x = new Double_t[npoints];
+ Double_t * y = new Double_t[npoints];
+ Double_t * s = new Double_t[npoints];
+ Double_t * ez = new Double_t[npoints];
+ Double_t * ex = new Double_t[npoints];
+ Double_t * ey = new Double_t[npoints];
+ Double_t * es = new Double_t[npoints];
+ Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
+
+
+ // Double_t chi=TMath::Pi()/2.0+phi;
+ Double_t chi=-TMath::Pi()-phi;
+ Double_t angold=0.0, tpang=0.0;
+ for(Int_t k = 0; k<npoints; k++) {
+ x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
+ y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
+ z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
+ if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
+ chisquare=9999.99;
+ cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
+ delete [] z;
+ delete [] x;
+ delete [] y;
+ delete [] s;
+ delete [] ez;
+ delete [] ex;
+ delete [] ey;
+ delete [] es;
+ return 12;