New Class added to do, Tracking and vertexing for ITS Stand alone using
authornilsen <nilsen@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Mar 2002 11:32:33 +0000 (11:32 +0000)
committernilsen <nilsen@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Mar 2002 11:32:33 +0000 (11:32 +0000)
Riemann Sphere mapping, with ITSRF.C macro and updated files for compiling.

ITS/AliITSRiemannFit.cxx [new file with mode: 0644]
ITS/AliITSRiemannFit.h [new file with mode: 0644]
ITS/ITSLinkDef.h
ITS/ITSRF.C [new file with mode: 0644]
ITS/Makefile
ITS/libITS.pkg

diff --git a/ITS/AliITSRiemannFit.cxx b/ITS/AliITSRiemannFit.cxx
new file mode 100644 (file)
index 0000000..7a0d32c
--- /dev/null
@@ -0,0 +1,1064 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ *                                                                        *
+ *                                                                        *
+ * /////////////////////////////////////////////////////////////////////  *
+ *                                                                        *
+ * This class performs a fast fit of helices going through the <=6        *
+ * points of the ITS, with the goal of studying tracking and              *
+ * vertexing performances.                                                *
+ * Generated kinematics is used to take into account different weights    *
+ * associated to points in different layers (with different multiple      *
+ * scattering-originated errors).                                         *
+ *                                                                        *
+ *   Based on the work by A. Strandlie, R. Fruhwirth                      *
+ *                                                                        *
+ *   First implementation by N. Bustreo, R. Turrisi - July 2000           *
+ *                                                                        *
+ *   Further modifications by A. Dainese, R. Turrisi                      *
+ *                                                                        *
+ *   Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it                 *
+ *                                                                        *
+ * **************************************************************************/
+//
+//
+//       Modified November, 7th 2001 by Rosario Turrisi 
+//       (rosario.turrisi@pd.infn.it)
+//
+//       FitHelix returns different values. 0=ok, >0 =problem
+//       void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
+//
+//
+//       Modified July, 30th 2001 by Rosario Turrisi 
+//       (rosario.turrisi@pd.infn.it)
+//       
+//        Fit for z now in (z,s) plane.
+//        Returns parameters in order to write the helix equation
+//        and find the right phase/initial point.
+//
+//     "PROPER WEIGHTS":  (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
+//
+#include "AliITSRiemannFit.h"
+#include "AliRun.h"
+#include "TClonesArray.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "TH2.h"
+#include "TMath.h"
+#include "TF1.h"
+#include "TGraphErrors.h"
+#include "TMinuit.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TRandom.h"
+#include "TParticle.h"
+#include "TFile.h"
+#include "AliITSRecPoint.h"
+#include "AliITSgeom.h"
+#include "AliITSmodule.h"
+
+ClassImp(AliITSRiemannFit)
+
+
+AliITSRiemannFit::AliITSRiemannFit() {
+  ///////////////////////////////////////////////////////////
+  // Default constructor.
+  // Set everything to zero.
+  ////////////////////////////////////////////////////////////
+
+  fSizeEvent     = 0;
+  fPoints        = 0;
+  fPrimaryTracks = 0;
+  fPointRecs     = 0;
+  //
+  //  test erase
+//    fspdi          = 0;
+//    fspdo          = 0;
+  for(Int_t i=0;i<6;i++)fPLay[i] = 0;
+  
+}
+//----------------------------------------------------------------------
+
+AliITSRiemannFit::~AliITSRiemannFit() {
+  ///////////////////////////////////////////////////////////
+  // Default destructor.
+  // if arrays exist delete them. Then set everything to zero.
+  ////////////////////////////////////////////////////////////
+   if(fPointRecs!=0){
+      for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
+      delete[] fPointRecs;
+   } // end if fPointRecs!=0
+  fSizeEvent     = 0;  
+  fPointRecs     = 0;
+  fPoints        = 0;
+  fPrimaryTracks = 0;
+  //
+  // test erase
+//    fspdi          = 0;
+//    fspdo          = 0;
+  for(Int_t i=0;i<6;i++)fPLay[i] = 0;
+  return;
+}
+//----------------------------------------------------------------------
+
+AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) {
+  ///////////////////////////////////////////////////////////
+  // Constructor.
+  // Set fSizeEvent to size and fPrimaryTracks to ntracks.
+  // Others to zero.
+  ////////////////////////////////////////////////////////////
+
+  fSizeEvent     = size;
+  fPoints        = 0;
+  fPrimaryTracks = ntracks;
+  //
+  // test erase
+//    fspdi          = 0;
+//    fspdo          = 0;
+  Point_tl *first = new Point_tl[fSizeEvent];
+  Point_tl **PointRecs = new Point_tl*[fSizeEvent];
+  for(Int_t i=0;i<6;i++)fPLay[i] = 0;
+  for(Int_t j=0;j<fSizeEvent;j++)  // create an array of struct
+    PointRecs[j] = &(first[j]);   
+}
+// ---------------------------------------------------------------------
+
+void FillPoints(Point_tl **Points,Int_t &index,Float_t *xpoint,
+               Float_t *error,
+               TLorentzVector PE,TLorentzVector OT,Int_t *id,
+               Int_t track,const Char_t *name,Int_t code,
+               Float_t phiorigin){
+  ///////////////////////////////////////////////////////////////////////
+  // Fill the structure Point_tl with the proper data
+  //
+  //////////////////////////////////////////////////////////////////////
+  Float_t PI2 = 2.0*TMath::Pi();
+  Float_t phi,r,x,y,z;
+  Int_t i;
+  i = index;
+  x = xpoint[0];
+  y = xpoint[1];
+  z = xpoint[2];
+  r = sqrt(x*x+y*y);
+  phi = TMath::ATan2(y,x);
+  if(phi<0.0) phi += PI2;
+  Points[i]->phi = phi;
+  Points[i]->eta  = -0.5*tan(0.5*TMath::ATan2(r,z));
+  Points[i]->fx = x;
+  Points[i]->fy = y;
+  Points[i]->fz = z;
+  Points[i]->fdx = error[0];
+  Points[i]->fdy = error[1];
+  Points[i]->fdz = error[2];
+  Points[i]->fr = r;
+  Points[i]->track  = track;
+  Points[i]->lay    = id[0],
+  Points[i]->lad    = id[1];
+  Points[i]->det    = id[2];
+  Points[i]->fMomentum = PE;
+  Points[i]->fOrigin   = OT;
+  Points[i]->fPt   = sqrt(PE.X()*PE.X()+PE.Y()*PE.Y());
+  Points[i]->fCode = code;
+  Points[i]->fName = name;
+  Points[i]->vertexPhi = phiorigin;
+  index++;
+  return;
+  
+}
+// -----------------------------------------------------------------------
+
+void AliITSRiemannFit::InitPoints(Int_t evnt,Int_t ntracks,AliITS *ITS,
+                                   TTree *TR,Int_t nparticles){
+  //////////////////////////////////////////////////////////////////////
+  // Fill the class member fPointRecs with the reconstructed points
+  // Set All other members to the real values
+  //
+  /////////////////////////////////////////////////////////////////////
+  printf("\n ************* Starting Init Points *************\n");
+  TParticle *part;
+  AliITSgeom *gm = (AliITSgeom*)ITS->GetITSgeom();
+  //get pointer to modules array
+  TObjArray *ITSmodules = ITS->GetModules();
+  Int_t nmodules=ITSmodules->GetEntriesFast();
+  printf("nmodules = %d \n",nmodules);
+  // Get the points from points file
+  AliITSmodule *itsModule;
+  Int_t mod,irec;
+  Stat_t nent;
+  AliITSRecPoint  *recp;
+  nent=TR->GetEntries();
+  TClonesArray *ITSrec  = ITS->RecPoints();
+
+  Int_t TotRP=0;
+  for (mod=0; mod<nmodules; mod++) {
+    itsModule=(AliITSmodule*)ITSmodules->At(mod);
+    ITS->ResetRecPoints();
+    TR->GetEvent(mod);
+    Int_t nrecp = ITSrec->GetEntries();
+    if(!nrecp) continue;
+    TotRP += nrecp;
+  }
+
+  Int_t iMAX = TotRP;
+  fPrimaryTracks = ntracks;
+  fParticles     = nparticles;
+  Point_tl *global = new Point_tl[iMAX];
+  fPointRecs = new Point_tl*[iMAX];
+  //
+  // test erase
+//    Point_tl *first = new Point_tl[iMAX];
+//    Point_tl *second = new Point_tl[iMAX];
+//    fspdi = new Point_tl*[iMAX];
+//    fspdo = new Point_tl*[iMAX];
+  for(Int_t j=0;j<iMAX;j++) {
+    fPointRecs[j] = &(global[j]);
+    //
+    // test erase
+//      fspdi[j]      = &(first[j]);
+//      fspdo[j]      = &(second[j]);
+  }
+  
+  Int_t ieta=0,ieta2=0;
+  Int_t i,id[4],idold[4];
+  Int_t track=0;//         // track  of hit
+  Float_t xpoint[3],error_plus[3],error_minus[3],global_error[3];       // position and error of the point
+  TLorentzVector OT,PE;
+  Float_t locals[3],locals_error[3],locals_plus[3],locals_minus[3]; // local position and local errors
+  Float_t Phi;
+  Int_t code;
+  const char *name;
+  Int_t layer,ladder,detector;
+  Float_t xcluster,zcluster;
+  Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
+  for (mod=0; mod<nmodules; mod++) {
+    itsModule=(AliITSmodule*)ITSmodules->At(mod);
+    ITS->ResetRecPoints();
+    TR->GetEvent(mod);
+    Int_t nrecp = ITSrec->GetEntries();
+    if (!nrecp) continue;
+    itsModule->GetID(layer,ladder,detector);
+
+    for (irec=0;irec<nrecp;irec++) {
+      recp   = (AliITSRecPoint*)ITSrec->UncheckedAt(irec);
+      track=recp->fTracks[0];
+      if(track <0 ) continue;
+      xcluster=recp->GetX();     // x on cluster
+      zcluster=recp->GetZ();     // z on cluster
+      part   = (TParticle*) gAlice->Particle(track);    
+      part->ProductionVertex(OT);  // set the vertex 
+      part->Momentum(PE);          // set the vertex momentum
+      name      = part->GetName();
+      code      = part->GetPdgCode();
+      Phi       = part->Phi();
+      id[0]=layer;
+      id[1]=ladder;
+      id[2]=detector;
+      id[3]=irec;
+      locals[0]=xcluster;     // x on cluster
+      locals[1]=0.0;          // y on cluster
+      locals[2]=zcluster;     // z on cluster
+      locals_error[0]=sqrt(recp->GetSigmaX2());
+      locals_error[1]=0.0;
+      locals_error[2]=sqrt(recp->GetSigmaZ2());
+      locals_plus[0]=xcluster+sqrt(recp->GetSigmaX2());       // x on cluster
+      if(layer==1||layer==2) locals_plus[1]=0.0150/2;         // y on cluster
+      else if(layer==3||layer==4) locals_plus[1]=0.0280/2;    // y on cluster
+      else if(layer==5||layer==6) locals_plus[1]=0.0300/2;    // y on cluster
+      locals_plus[2]=zcluster+sqrt(recp->GetSigmaZ2());       // z on cluster
+      locals_minus[0]=xcluster-sqrt(recp->GetSigmaX2());      // x on cluster
+      locals_minus[1]=0.0;                                    // y on cluster
+      locals_minus[2]=zcluster-sqrt(recp->GetSigmaZ2());      // z on cluster
+
+      gm->LtoG(layer,ladder,detector,locals,xpoint);
+      gm->LtoG(layer,ladder,detector,locals_plus,error_plus);
+      gm->LtoG(layer,ladder,detector,locals_minus,error_minus);
+      global_error[0]=0.5*TMath::Abs(error_plus[0]-error_minus[0]);
+      global_error[1]=0.5*TMath::Abs(error_plus[1]-error_minus[1]);
+      global_error[2]=0.5*TMath::Abs(error_plus[2]-error_minus[2]);
+      if(track<ntracks) {
+       if(TMath::Abs(part->Eta())<=1.0) ieta++;
+       if(TMath::Abs(part->Eta())<=0.5) ieta2++;
+      }
+      if(!(id[0]==idold[0]&&id[1]==idold[1]&&
+          id[2]==idold[2]&&id[3]==idold[3])) {
+       FillPoints(fPointRecs,num,xpoint,global_error,PE,OT,id,track,name,code,Phi);
+       //
+       // test erase   
+       switch (idold[0]) {
+       case 1:
+         nspdi++;
+         break;
+       case 2:
+         nspdo++;
+         break;
+       case 3:
+         nsddi++;
+         break;
+       case 4:
+         nsddo++;
+         break;
+       case 5:
+         nssdi++;
+         break;
+       case 6:
+         nssdo++;
+         break;
+       }
+//     if(idold[0]==1){
+//       FillPoints(fspdi,nspdi,xpoint,global_error,PE,OT,id,track,name,code,Phi);  
+//     }
+//     if(idold[0]==2){
+        
+//       FillPoints(fspdo,nspdo,xpoint,global_error,PE,OT,id,track,name,code,Phi);
+//     }
+//     if(idold[0]==3){
+//       nsddi++;
+//     }
+//     if(idold[0]==4){
+//       nsddo++;
+//     }
+//     if(idold[0]==5){
+//       nssdi++;
+//     }
+//     if(idold[0]==6){
+//       nssdo++;
+//     }
+       for(i=0;i<4;i++) idold[i] = id[i];
+       for(i=0;i<3;i++) xpoint[i]    = 0.0;
+      } // end if id != idold
+    } // end for irec
+  }// end for mod
+
+  fPoints = num;
+  fSizeEvent = num;
+  fPLay[0] = nspdi ;
+  fPLay[1] = nspdo ;
+  fPLay[2] = nsddi ;
+  fPLay[3] = nsddo ;
+  fPLay[4] = nssdi ;
+  fPLay[5] = nssdo ;
+  printf("%d primary tracks in eta=+-1\n",ieta);
+  printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
+  printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
+  printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
+  printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
+  printf("Points on all Layers: %d\n",num);
+  printf("\n ************* Init Points Finished *************\n");
+  return;
+}
+// ------------------------------------------------------------------------
+///////////////////////////////////////////////////////////
+// Functions for sorting the fPointRecs array
+///////////////////////////////////////////////////////////
+Bool_t SortZ(const Point_tl *s1,const Point_tl *s2){
+  // Z sorting function for qsort.
+   Float_t a;
+
+   a = s1->fz - s2->fz;
+   if(a<0.0) return kTRUE;
+   if(a>0.0) return kFALSE;
+   return kFALSE;
+}
+Bool_t SortTrack(const Point_tl *s1,const Point_tl *s2){
+  // track sorting function for qsort.
+   Float_t a;
+
+   a = s1->track - s2->track;
+   if(a<0.0) return kTRUE;
+   if(a>0.0) return kFALSE;
+   return kFALSE;
+}
+void hpsortTrack(Point_tl **ra,Int_t n){
+   Int_t i,ir,j,l;
+   Point_tl *rra;
+
+   if(n<2) return;
+
+   l  = ((n-1) >> 1) +1; // divide 2 + 1
+   ir = n-1;
+   for(;;){
+     if(l>0){
+        rra = ra[--l];  // decrement first
+     }else{
+        rra    = ra[ir];
+        ra[ir] = ra[0];
+        if(--ir == 0){  // decrement first
+           ra[0] = rra;
+           break;
+        } // if --ra == 0 
+     } // end l>0 
+     i = l;
+     j = l+1;
+     while(j<=ir){
+        if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
+        if( SortTrack(rra,ra[j]) ){
+           ra[i] = ra[j];
+           i = j;
+           j <<= 1; // time 2.
+        }else{
+           break;
+        } // end if func() 
+     } // end while 
+     ra[i] = rra;
+   } // end for ever 
+}
+void hpsortZ(Point_tl **ra,Int_t n){
+   Int_t i,ir,j,l;
+   Point_tl *rra;
+
+   if(n<2) return;
+
+   l  = ((n-1) >> 1) +1; // devide 2 + 1
+   ir = n-1;
+   for(;;){
+     if(l>0){
+        rra = ra[--l];  // decrament first
+     }else{
+        rra    = ra[ir];
+        ra[ir] = ra[0];
+        if(--ir == 0){  // decrament first
+           ra[0] = rra;
+           break;
+        } // if --ra == 0 
+     } // end l>0 
+     i = l;
+     j = l+1;
+     while(j<=ir){
+        if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
+        if( SortZ(rra,ra[j]) ){
+           ra[i] = ra[j];
+           i = j;
+           j <<= 1; // time 2.
+        }else{
+           break;
+        } // end if func() 
+     } // end while 
+     ra[i] = rra;
+   } // end for ever 
+}
+//-----------------------------------------------------------------------
+////////////////////////////////////////////////////////////////////
+//      Sorting functions
+///////////////////////////////////////////////////////////////////
+Int_t Partition(Int_t array[],Int_t left,Int_t right){
+  Int_t val = array[left];
+  Int_t lm = left - 1;
+  Int_t rm = right + 1;
+  for(;;) {
+    do 
+      rm--;
+    while
+      (array[rm]>val);
+    do 
+      lm++;
+    while
+      (array[lm]<val);
+    if(lm<rm){
+      Int_t tempr = array[rm];
+      array[rm]=array[lm];
+      array[lm]=tempr;
+    }
+    else
+      return rm;
+  }
+
+  return 1;
+}
+
+///////////////////////////////////////////////////////////////////////
+
+void AliITSRiemannFit::WritePoints(void) {
+  /////////////////////////////////////////////////////////////////////
+  // write the data in a file (temporary ascii)
+  /////////////////////////////////////////////////////////////////////
+  FILE *ascii= fopen("AsciiPoints.dat","w");
+  for(Int_t i=0;i<fPoints;i++) {
+    fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->lay,
+           fPointRecs[i]->track,fPointRecs[i]->fx,fPointRecs[i]->fy,
+           fPointRecs[i]->fz);
+  }
+  fclose(ascii);
+  return;
+}
+//-----------------------------------------------------------------------
+
+void AliITSRiemannFit::ReadPoints(void) {
+  //////////////////////////////////////////////////////////////////////
+  // read the filled array
+  /////////////////////////////////////////////////////////////////////
+  hpsortTrack(fPointRecs,fPoints);
+  for(Int_t i=0;i<fPoints;i++) 
+    printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
+          i,fPointRecs[i]->lay,fPointRecs[i]->track,fPointRecs[i]->fx,
+          fPointRecs[i]->fy,fPointRecs[i]->fz,fPointRecs[i]->fOrigin.X(),
+          fPointRecs[i]->fOrigin.Y(),fPointRecs[i]->fOrigin.Z(),
+          fPointRecs[i]->fPt,fPointRecs[i]->fName);
+  return; 
+}
+//-----------------------------------------------------------------------
+
+Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
+                                   Double_t &x1,Double_t &x2,Double_t &x3){
+   //////////////////////////////////////////////
+   ///  Solve cubic equation:
+   ///  x^3 + a*x^2 +b*x + c
+   ///  
+   ///  returns  x1 , x2 , x3
+   ////////////////////////////////////////
+   
+   Double_t Q = ((a*a - 3*b)/9);
+   Double_t R = ((2*a*a*a - 9*a*b +27*c)/54);
+   Double_t theta;
+   Double_t F = -2*sqrt(Q);
+   Double_t g = a/3;
+   Double_t PI2 = TMath::Pi()*2;
+
+  if( R*R>Q*Q*Q ) {
+    cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
+    x1 = 9999999;
+    x2 = 9999999;
+    x3 = 9999999;
+    return 0;
+  }
+
+  theta = TMath::ACos( R/sqrt(Q*Q*Q));
+
+  x1 = (F*TMath::Cos(theta/3))-g;
+  x2 = (F*TMath::Cos((theta+PI2)/3))-g;
+  x3 = (F*TMath::Cos((theta-PI2)/3))-g;
+  
+  return 1;
+}
+//-----------------------------------------------------------------
+
+void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
+  ///////////////////////////////////////////////////////////////////////
+  //   This function apllies the transformation in the Riemann sphere
+  //   for xy plane
+  ///////////////////////////////////////////////////////////////////////
+  Float_t *R = new Float_t[npoints];
+  Float_t *Theta = new Float_t[npoints];
+  Float_t PI2 = 2*TMath::Pi();
+  Float_t x=0,y=0,z=0;
+  
+  for(Int_t i=0;i<npoints;i++) {
+    R[i]     = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
+    Theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
+    if(Theta[i]<0) Theta[i]+=PI2;
+    x     = R[i]*cos(Theta[i])/(1+R[i]*R[i]);
+    y     = R[i]*sin(Theta[i])/(1+R[i]*R[i]);
+    z     = R[i]*R[i]/(1+R[i]*R[i]);
+    To[i]->SetXYZ(x,y,z);
+  }
+  delete[] R;
+  delete[] Theta;
+  return;
+}
+
+
+//---------------------------------------------------------------------
+
+Int_t FitLinear(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
+  //
+  /////////////////////////////////////////////////////////////////////// 
+  Int_t direction=0;
+  Double_t z[npoints],x[npoints],y[npoints],s[npoints];
+  Double_t ez[npoints],ex[npoints],ey[npoints],es[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; 
+      return 12;
+    }
+    Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
+    if( (x[k]-thu0)<0 ) {
+      if (ang1*angold<0) {
+       tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
+       ang1=tpang;
+      }
+    }
+    angold=ang1;
+    if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
+    s[k] = (ang1+chi)/omega;
+    es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
+  }
+  if ( TMath::Abs(direction) != (npoints-1) ) {return 11;} 
+
+  TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
+  fitHist->Fit("pol1","Q");
+  z0  = fitHist->GetFunction("pol1")->GetParameter(0);
+  vpar  = fitHist->GetFunction("pol1")->GetParameter(1);
+  ez0 = fitHist->GetFunction("pol1")->GetParError(0);
+  evpar = fitHist->GetFunction("pol1")->GetParError(1);
+  chisquare = fitHist->GetFunction("pol1")->GetChisquare();
+  zData.Set(z0,vpar);
+  zError.SetXYZ(ez0,evpar,chisquare);
+
+  Double_t Sigmas=0.; 
+  Double_t Sigmaz=0.;
+  Double_t Avs=0.;
+  Double_t Avz=0.;
+  Double_t Avsz=0.;
+
+  for(Int_t j = 0; j < npoints; j++) {
+    Avs  += s[j];
+    Avz  += z[j]; 
+    Avsz += s[j]*z[j];
+  }
+    Avs  /= (Double_t)npoints;
+    Avz  /= (Double_t)npoints; 
+    Avsz /= (Double_t)npoints;
+
+  for(Int_t l = 0; l < npoints; l++) {
+    Sigmas += (s[l]-Avs)*(s[l]-Avs);
+    Sigmaz += (z[l]-Avz)*(z[l]-Avz);
+  }
+  Sigmas /=(Double_t)npoints;
+  Sigmaz /=(Double_t)npoints;
+
+  Sigmas = sqrt(Sigmas);
+  Sigmaz = sqrt(Sigmaz);
+
+  CorrLin = (Avsz-Avs*Avz)/(Sigmas*Sigmaz);
+
+  return 0;
+}
+
+//-------------------------------------------------------------------
+Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Int_t charge,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
+                                  Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
+                                  Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
+                                  Int_t first=1,Int_t second=1,Int_t third=1,Int_t fourth=1,Int_t fifth=1,Int_t sixth=1) {
+  ///////////////////////////////////////////////////////////////////////
+  //  This function finds the helix paramenters 
+  //  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.0  (yes, one-dot-zero !)
+  //
+  //
+  ///////////////////////////////////////////////////////////////////////
+  //
+  //  All this stuff relies on this hypothesis !!!
+  //
+//   FILE *pout=fopen("chisql.dat","a");
+  Int_t ierr = 0, ierrl=0;
+  omega = 1.0e-2;
+
+  Int_t bitlay[6]={1,1,1,1,1,1}; 
+  bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
+  fd0 = -9999;   // No phisycs value
+  u0 = -9999.9999; // parameters of helix - strange value...
+  v0 = -9999.9999; // parameters of helix - strange value...
+  rho = -9999.9999; // parameters of helix -unphysical strange value...
+  Int_t Layer = 0;
+  const Char_t* name = 0;
+  Int_t i=0,k=0;
+  Int_t iMAX = 50;
+  Int_t N = 0;
+  Int_t npl[6]={0,0,0,0,0,0};
+  Double_t P = 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[iMAX];
+  TVector3 **original  = new TVector3*[iMAX];
+  TVector3 *rie       = new TVector3[iMAX];
+  TVector3 **Riemann  = new TVector3*[iMAX];
+  TVector3 *err        = new TVector3[iMAX];
+  TVector3 **errors    = new TVector3*[iMAX];
+  TVector3 *linerr        = new TVector3[iMAX];
+  TVector3 **linerrors    = new TVector3*[iMAX];
+  Double_t Weight[iMAX];
+
+  for(i=0;i<iMAX;i++){
+    original[i]   = &(ori[i]);
+    Riemann[i]   = &(rie[i]);
+    errors[i]     = &(err[i]);
+    linerrors[i]  = &(linerr[i]);
+  }
+  for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
+  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 
+  Double_t Radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
+  Double_t sigma_MS = 0;
+  TVector3 Vec,VecNor;
+
+// Select RecPoints belonging to the track
+  for(k =0;k<fPoints;k++){
+    if(fPointRecs[k]->track==tracknumber) {
+      name = fPointRecs[k]->fName;
+      Pt = fPointRecs[k]->fPt;
+      Layer = fPointRecs[k]->lay;
+      Int_t ilay = Layer-1;
+      if(npl[ilay]!=0) continue;
+      if(bitlay[ilay] == 1) {
+       original[N]->SetXYZ(0.1*fPointRecs[k]->fx,0.1*fPointRecs[k]->fy,0.1*fPointRecs[k]->fz);
+       errors[N]->SetXYZ(0.1*fPointRecs[k]->fdx,0.1*fPointRecs[k]->fdy,0.1*fPointRecs[k]->fdz);
+        sigma_MS = (Radiusdm[Layer]-Radiusdm[0])*0.000724/P;// beam pipe contribution
+       for(Int_t j=1;j<Layer;j++) {
+         sigma_MS += (Radiusdm[Layer]-Radiusdm[j])*0.00136/P;
+       }
+       Weight[N] = ( 1 + original[N]->Perp2() )*( 1+ original[N]->Perp2() )/
+         ( errors[N]->Perp2() + sigma_MS*sigma_MS );
+       linerrors[N]->SetXYZ(errors[N]->X(),errors[N]->Y(),sqrt(errors[N]->Z()*errors[N]->Z()+sigma_MS*sigma_MS));
+       N++;
+       npl[ilay]++;
+      }  // end if on layer 
+    }    //end if track==tracknumber
+  }      //end for k
+  //
+  //    6 points, no more, no less
+  //
+  if(original[5]->X() == 9999 || original[6]->X() != 9999)  
+    {
+      return 1;   // not enough points
+    }
+  
+  //
+  //
+  //
+  //  FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
+  //  
+  
+  RiemannTransf(N,original,Riemann);  
+
+  Double_t Sum_Weights = 0.0; // sum of weights factor
+
+  for(Int_t j=0;j<N;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();
+    Sum_Weights+=Weight[j];
+  }
+  
+  xbar /= Sum_Weights;
+  ybar /= Sum_Weights;
+  zbar /= Sum_Weights;
+  
+  for(Int_t j=0;j<N;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);
+  }
+  
+  A11 /= N;
+  A12 /= N;
+  A22 /= N;
+  A23 /= N;
+  A13 /= N;
+  A33 /= N;
+  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 Check_Cubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
+  if(Check_Cubic !=1 ){
+    printf("Track %d Has no real solution continuing ...\n",tracknumber);
+    return 2;
+  }
+  
+// **************** 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;
+  Vec.SetXYZ(x1,x2,x3);
+  VecNor = Vec.Unit();
+  n1 = VecNor.X();
+  n2 = VecNor.Y();
+  n3 = VecNor.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;
+  fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
+  fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
+  fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
+  fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
+  fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
+  fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->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=FitLinear(N,original,linerrors,omega,u0,v0,fphi,zData,zError,CorrLin);
+  chisql=zError.Z();
+//   fprintf(pout,"%f \n",chisql); 
+  z0=zData.X();
+  vpar=zData.Y();
+  fCorrLin = CorrLin;
+  ierr = (ierrl > ierr ? ierrl : ierr);
+//   fclose(pout);
+  return ierr;
+}
+
+//-----------------------------------------------------------------------------
+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;
+    }
+//      for(i=0;i<fSizeEvent/6;i++) 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 << fspdi[i][j].fName;
+//        lRb << fspdi[i][j].phi;
+//        lRb << fspdi[i][j].eta;
+//        lRb << fspdi[i][j].vertexPhi;
+//      }
+//      for(i=0;i<fSizeEvent/6;i++) 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 << fspdo[i][j].fName;
+//        lRb << fspdo[i][j].phi;
+//        lRb << fspdo[i][j].eta;
+//        lRb << fspdo[i][j].vertexPhi;
+//    }
+  } // end if reading
+}
diff --git a/ITS/AliITSRiemannFit.h b/ITS/AliITSRiemannFit.h
new file mode 100644 (file)
index 0000000..beb880e
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef ALIITSRIEMANNFIT_H
+#define ALIITSRIEMANNFIT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "TLorentzVector.h"
+#include "TTree.h"
+#include "AliITS.h"
+#include "TVector3.h"
+
+struct Point_tl{
+  Int_t lay,lad,det,track;
+  Float_t fx,fy,fz,fr;               // global position of point 
+  Float_t fdE,fdx,fdy,fdz;               // Errors
+  TLorentzVector fOrigin,fMomentum;  // position and momentum of 
+                                     //  particle at its origin
+  Int_t fCode;                       // Geant code of particle
+  const Char_t *fName;
+  Float_t fPt;                       // Pt at the origin
+  Float_t phi,eta,vertexPhi;         // phi eta on layer and phi on vertex
+};
+
+
+class AliITSRiemannFit : public TObject{
+ public:
+  AliITSRiemannFit();
+  AliITSRiemannFit(Int_t size,Int_t ntracks);
+  ~AliITSRiemannFit();
+  
+  Int_t GetSize() const {return this->fSizeEvent;}
+  Int_t GetPrimaryTracks() const {return this->fPrimaryTracks;}
+  Int_t GetPoints() const {return this->fPoints;}
+  Int_t GetParticles() const {return this->fParticles;}
+  Int_t GetLayPoints(Int_t layer) const {return this->fPLay[layer-1];}
+  Point_tl **GetPointRecs() const {return this->fPointRecs;}
+  Float_t GetX(Int_t i) const {return this->fPointRecs[i]->fx;}
+  Float_t GetY(Int_t i) const {return this->fPointRecs[i]->fy;}
+  Float_t GetZ(Int_t i) const {return this->fPointRecs[i]->fz;}
+  Float_t GetdX(Int_t i) const {return this->fPointRecs[i]->fdx;}
+  Float_t GetdY(Int_t i) const {return this->fPointRecs[i]->fdy;}
+  Float_t GetdZ(Int_t i) const {return this->fPointRecs[i]->fdz;}
+  
+  void     InitPoints(Int_t evnt,Int_t ntracks,AliITS *ITS,TTree *TR,Int_t nparticles);
+  void     WritePoints(void);
+  void     ReadPoints(void);
+  static Int_t SolveCubic(Double_t a,Double_t b,Double_t c,Double_t&,Double_t&,Double_t&);
+  Int_t FitHelix(Int_t tracknumber,Int_t charge,Double_t Px,Double_t Py,Double_t Pz,
+                Double_t& fd0,Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,
+                Double_t& omega, Double_t& z0,
+                Double_t& vpar,Double_t& chisql,Double_t& fCorrLin,Double_t& fFit,
+                Int_t first=1,Int_t second=1,Int_t third=1,Int_t fourth=1,Int_t fifth=1,Int_t sixth=1);  
+ private:
+  Int_t fSizeEvent;      // size of array 
+  Int_t fPrimaryTracks;  // number of primary tracks in the event
+  Int_t fPoints;         // number of Reconstructed Points in the event
+  Int_t fParticles;      // number of particles in the event
+  Int_t fPLay[6];           // number of points in each layer
+  Point_tl **fPointRecs;
+  //
+  // test erase
+/*    Point_tl **fspdi,**fspdo; // This are for the first two layers and vertex analysis */
+  
+  ClassDef(AliITSRiemannFit,1)  // Fast fit of helices on ITS RecPoints
+    };
+#endif
+
index 1698c43..7f4b748 100644 (file)
 #pragma link C++ class  AliITSIOTrack+;
 #pragma link C++ class  AliITSTrackerV1+;
 #pragma link C++ class  AliITSgeoinfo+;
+#pragma link C++ class  AliITSRiemannFit-;
 // New used for Alignment studdies
 //#pragma link C++ class  AliITSAlignmentTrack-;
 //#pragma link C++ class  AliITSAlignmentModule-;
diff --git a/ITS/ITSRF.C b/ITS/ITSRF.C
new file mode 100644 (file)
index 0000000..f4b8b22
--- /dev/null
@@ -0,0 +1,141 @@
+#define __COMPILED__
+#ifdef __COMPILED__
+#include <iostream.h>
+#include <AliITS.h>
+#include <AliITSRiemannFit.h>
+#include <AliRun.h>
+#include <TROOT.h>
+#include <TClonesArray.h>
+#include <TFitter.h>
+#include <TMinuit.h>
+#include <TF2.h>
+#include <TFile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TObjArray.h>
+#include <TParticle.h>
+#include <TTree.h>
+#include <TVector3.h>
+#include <TVirtualFitter.h>
+#endif
+
+//    Writes the data of helices and the position of the true II-ary, to perform fit 
+//    of II-ary and calculation of resolution (in another macro)
+//
+void ITSRF(const char *filename="galice.root", 
+          const char *outfile="RiemannFit.dat",
+          Int_t evnt=0) {
+  
+  // Connect the Root Galice file containing Geometry, Kine and Hits
+  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+  if (!file){ file = new TFile(filename);printf("Opening new file\n");}
+  printf("Root input file is %s\n",filename);
+  
+  // Get AliRun object from file or create it if not on file
+  if (!gAlice) {
+    gAlice = (AliRun*)file->Get("gAlice");
+    if (gAlice) printf("AliRun object found on file\n");
+    if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+  }
+  
+  // Get pointers to Alice detectors and Hits containers
+  AliITS  *ITS       = (AliITS*)gAlice->GetModule("ITS");
+  if(!ITS) {cout<<"not ITS"<<endl; return;}
+  Int_t nparticles = gAlice->GetEvent(evnt);   // number of particles generated (including daughters)
+
+  //
+  //  Treat ITS information per module
+  //
+  
+  Int_t nmodules;
+  ITS->InitModules(-1,nmodules); 
+  cout<<" Filling modules..."<<endl;
+  ITS->FillModules(evnt,-1,nmodules," "," ");
+
+  //
+  //  Get the recontruction tree
+  //
+  TTree *TR = gAlice->TreeR();
+  Int_t nent = (Int_t)TR->GetEntries();
+  
+  //
+  // Initialize Riemann-fit class
+  //
+  AliITSRiemannFit *rfit = new AliITSRiemannFit();
+  rfit->InitPoints(evnt,nent,ITS,TR,nparticles);
+  
+  const TVector3 Bfield(0.0,0.0,0.2); // ALICE magnetic field
+  //
+  //   Setting the variables
+  //
+  
+  TParticle *MPart;
+  
+  Double_t r0[3]={0.0,0.0,0.0}, w, chi;
+  Double_t Pt, Px, Py, Pz, charge, Ptot;
+  Double_t x0, y0, rho, fd0, fphi, z0, vpar, chisql;
+  Double_t CorrLin, ResSum; 
+  Int_t status[4]={0,0,0,0}, gfit_status, PdgPart;
+  
+  FILE *outdat=fopen(outfile,"w");
+
+  cout<<endl<<" Looping on "<< nparticles <<" particles..."<<endl;
+
+  Int_t nAcc=0;
+  for(Int_t part = 0; part < nparticles ; part++) {
+    if (!(part%100)) cout<<part<<"\r";
+    MPart = (TParticle*)gAlice->Particle(part);
+    PdgPart = MPart->GetPdgCode();
+
+    //
+    //  get track fit
+    // 
+    Pt = MPart->Pt();
+    Px = MPart->Px();
+    Py = MPart->Py();
+    Pz = MPart->Pz();
+    Ptot = MPart->P();
+    charge = PdgPart/TMath::Abs(PdgPart);  // +/- 1
+    gfit_status = (Int_t)rfit->FitHelix(part, charge, Px, Py, Pz, fd0, fphi, x0, y0, 
+                                       rho, w, z0, vpar, chisql, CorrLin, ResSum);
+    // With this offset chi = 0  at the 
+    // distance of closest approach in bending plane
+    //
+    chi=-TMath::Pi()-fphi;  
+    if(gfit_status > 0) {
+      switch(gfit_status) {
+      case 1:  status[0]++;break;
+      case 2:  status[1]++;break;
+      case 11:  status[2]++;break;
+      case 12:  status[3]++; break;
+      default: cout<<"CASO NON PREVISTO IN K:   "<<gfit_status<<endl;break;
+      }
+      continue; 
+    }
+
+    //
+    //
+    //
+    Int_t written=0;
+    written=fprintf(outdat,"%d ",   part); 
+    written=fprintf(outdat,"%f ",   r0[0]); 
+    written=fprintf(outdat,"%f ",   r0[1]);
+    written=fprintf(outdat,"%f \n", r0[2]);
+    written=fprintf(outdat,"%f \n", ResSum); 
+    written=fprintf(outdat," %f ",  rho);
+    written=fprintf(outdat,"%f ",   w);
+    written=fprintf(outdat,"%f ",   fphi); 
+    written=fprintf(outdat,"%f \n", vpar); 
+    nAcc++;
+
+  } // loop on particles
+  fclose(outdat);
+
+  printf(" missed fits because of (resp.):\n no six points, no solution to cubic equation, \n z-direction inversion, tangent limit \n");
+  printf(" %d %d %d %d %d \n",status[0],status[1],status[2],status[3]);
+
+
+  return;
+} 
+
+
index fb123d2..45aa9f0 100644 (file)
@@ -46,7 +46,7 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                AliV0vertex.cxx AliV0vertexer.cxx \
                AliCascadeVertex.cxx AliCascadeVertexer.cxx \
                AliITSglobalRecPoint.cxx AliITSneuralTrack.cxx \
-               AliITSneuralTracker.cxx 
+               AliITSneuralTracker.cxx AliITSRiemannFit.cxx 
 #              AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
 
 # Fortran sources
index 53bb5ab..7deb851 100644 (file)
@@ -35,7 +35,8 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                 AliITSPid.cxx AliITStrackV2Pid.cxx AliCascadeVertex.cxx \
                 AliCascadeVertexer.cxx  AliITSneuralTrack.cxx \
                 AliITSneuralTracker.cxx AliITSglobalRecPoint.cxx \
-                AliITSDigitizer.cxx AliITSreconstruction.cxx 
+                AliITSDigitizer.cxx AliITSreconstruction.cxx \
+               AliITSRiemannFit.cxx 
 #              AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
 
 HDRS:=  $(SRCS:.cxx=.h)