Removing spurious files
authorfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Jul 2000 16:26:01 +0000 (16:26 +0000)
committerfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Jul 2000 16:26:01 +0000 (16:26 +0000)
ITS/AnalysisITSAlignment.cxx [deleted file]
ITS/AnalysisITSAlignment.h [deleted file]
ITS/Makefile
ITS/vertex.cxx [deleted file]
ITS/vertex.h [deleted file]

diff --git a/ITS/AnalysisITSAlignment.cxx b/ITS/AnalysisITSAlignment.cxx
deleted file mode 100644 (file)
index 3aae967..0000000
+++ /dev/null
@@ -1,1695 +0,0 @@
-//////////////////////////////////////////////////////////////////////////
-//  Alice ITS first detector alignment program.                         //
-//                                                                      //
-// version: 0.0.0 Draft.                                                //
-// Date: April 18 1999                                                  //
-// By: Bjorn S. Nilsen                                                  //
-//                                                                      //
-//////////////////////////////////////////////////////////////////////////
-#include <stdio.h>
-#include <math.h>
-#include <fstream.h>
-#include <stdio.h>
-#include <time.h>
-#include "AliITS.h"
-#include "TTree.h"
-#include "TClonesArray.h"
-#include "TProfile.h"
-#include "TH2.h"
-#include "TArray.h"
-#include "TRandom.h"
-#include "TCanvas.h"
-#include "TFile.h"
-#include "TParticle.h"
-#include "AliRun.h"
-#include "AliITSgeom.h"
-#include "AnalysisITSAlignment.h"
-#include "AliITSstatistics2.h"
-#include "AliITSstatistics.h"
-#include "AliITSAlignmentTrack.h"
-#include "AliITSAlignmentModule.h"
-
-#define PRIMARYONLY 0  // if PIMARY_ONLY ==0 then all tracks
-#define MOMENTUMCUT 1.0 // Sets the momentum cut for the tracks.
-
-//Double_t resz[]={0.00680,0.00650,0.00265,0.00265,0.0830,0.0830};
-//Double_t resx[]={0.00105,0.00115,0.00340,0.00340,0.0020,0.0020};
-// From table 1.3 of ITS TDR. Spatial precisions and Cell size
-Double_t  resz[]={0.00700,0.00700,0.00280,0.00280,0.08300,0.08300};
-Double_t  resx[]={0.00120,0.00120,0.00380,0.00380,0.00200,0.00200};
-Double_t cellx[]={0.00500,0.00500,0.01500,0.01500,0.00950,0.00950};
-Double_t celly[]={0.00150,0.00150,0.00300,0.00300,0.00300,0.00300};
-Double_t cellz[]={0.03000,0.03000,0.03000,0.03000,4.00000,4.00000};
-
-//_________________________________________________________________________
-void HitsTo(ClustAl_tl *trk,Int_t &ntrk,Int_t nt,TTree *TH,AliITS *ITS,
-           Float_t nsigmaT1,Float_t nsigmaT2,Float_t nsigmaT3,
-           Float_t nsigmaR1,Float_t nsigmaR2,Float_t nsigmaR3){
-
-   // Local variables
-   Int_t        nClustMin = 3;
-   const Int_t  jmax = 50;
-   Int_t        t,nb,nh,h,i,ly=-1,ld=-1,dt=-1,tt=-1,j=0;
-   ClustAl_sl   ht[jmax];
-   Int_t        id[3],ie[3],nht[jmax],track;
-   Float_t      gp[3],lp[3],p[3],tof,gp2[3];
-   Float_t      ns;
-   Float_t      px[jmax],py[jmax],pz[jmax];
-   TClonesArray *ITShits = ITS->Hits();
-   AliITShit    *itsHit;
-   // initilize and creat AliITSgeom data
-   AliITSgeom   *gm  = ITS->GetITSgeom();
-   AliITSgeom   &gm2 = *gm;
-
-//   nt = 100; // lets speed things up.
-
-   Float_t tran[]={90.e-4/3.,30.e-4/3.,780.e-4/3.}; // r,rphi,z
-   Float_t  rot[]={0.4e-3/3.,3.0e-3/3.,100.e-3/3.}; // about r, rphi, z
-   tran[0] *= nsigmaT1;
-   tran[1] *= nsigmaT2;
-   tran[2] *= nsigmaT3;
-   rot[0]  *= nsigmaR1;
-   rot[1]  *= nsigmaR2;
-   rot[2]  *= nsigmaR3;
-
-   printf("Randomize by tran=%fcm %fcm %fcm\n",tran[0],tran[1],tran[2]);
-   printf("Randomize by rot=%frad %frad %frad\n",rot[0],rot[1],rot[2]);
-
-   gm2.RandomCylindericalChange(tran,rot);
-
-   for(i=0;i<jmax;i++){// zero for for first track.
-       ht[i].xl = ht[i].xg = 0.0;
-       ht[i].yl = ht[i].yg = 0.0;
-       ht[i].zl = ht[i].zg = 0.0;
-       ht[i].lay = ht[i].lad = ht[i].det = 0;
-       px[i]    = 0.0;
-       py[i]    = 0.0;
-       pz[i]    = 0.0;
-       nht[i]   = 0;
-   } // end init to zero.
-
-   j = 0;
-   for(t=0;t<nt;t++){ // Loop over tracks.
-      gAlice->ResetHits();
-      nb = TH->GetEvent(t);
-      nh = ITShits->GetEntriesFast();
-      for(h=0;h<nh;h++){ // hits
-         itsHit = (AliITShit *)ITShits->UncheckedAt(h);
-        track  = itsHit->fTrack;
-        if(track != tt){ // if new track
-            if(j>nClustMin)if(tt<nt||PRIMARYONLY==0){ 
-                 // if there are enought points save data
-                trk[ntrk].track  = tt;
-                trk[ntrk].nclust = j+1;
-                ns               = 1./(Float_t (nht[0]));
-                trk[ntrk].px     = px[0]*ns;
-                trk[ntrk].py     = py[0]*ns;
-                trk[ntrk].pz     = pz[0]*ns;
-                trk[ntrk].pt     = sqrt(px[0]*px[0]+py[0]*py[0])*ns;
-                trk[ntrk].p      = sqrt(px[0]*px[0]+py[0]*py[0]+
-                                        pz[0]*pz[0]            )*ns;
-                trk[ntrk].clust  = new ClustAl_sl[j+1];
-                for(i=0;i<=j;i++){ // loop over all detector hit.
-                    ns                     = 1./(Float_t(nht[i]));
-                    trk[ntrk].clust[i].lay = ie[0] = ht[i].lay;
-                    trk[ntrk].clust[i].lad = ie[1] = ht[i].lad;
-                    trk[ntrk].clust[i].det = ie[2] = ht[i].det;
-                    lp[0]                  = ht[i].xl*ns;
-                    lp[1]                  = ht[i].yl*ns;
-                    lp[2]                  = ht[i].zl*ns;
-                    // set in detector resolution on local coordiante system
-                    // z->zl rphi->xl yl-> thickness leave it alone
-                    lp[0] = lp[0]/resx[ie[0]-1];
-                    lp[0] = (Float_t)((Int_t)lp[0]);
-                    lp[0] = lp[0]*resx[ie[0]-1] + 0.5*resx[ie[0]-1];
-                    lp[2] = lp[2]/resz[ie[0]-1];
-                    lp[2] = (Float_t)((Int_t)lp[2]);
-                    lp[2] = lp[2]*resz[ie[0]-1] + 0.5*resz[ie[0]-1];
-                    //
-                    trk[ntrk].clust[i].xl  = lp[0];
-                    trk[ntrk].clust[i].yl  = lp[1];
-                    trk[ntrk].clust[i].zl  = lp[2];
-                    gm2.LtoG(ie,lp,gp);
-                    gm->LtoG(ie,lp,gp2);
-                    trk[ntrk].clust[i].xg  = gp[0];
-                    trk[ntrk].clust[i].yg  = gp[1];
-                    trk[ntrk].clust[i].zg  = gp[2];
-                } // end for i: detectors clust loop
-                ntrk++; // set up for next track
-            } // end if j>nClustMin
-            for(i=0;i<jmax;i++){ // zero out structure
-                ht[i].xl = ht[i].xg = 0.0;
-                ht[i].yl = ht[i].yg = 0.0;
-                ht[i].zl = ht[i].zg = 0.0;
-                ht[i].lay = ht[i].lad = ht[i].det = 0;
-                px[i]    = 0.0;
-                py[i]    = 0.0;
-                pz[i]    = 0.0;
-                nht[i]   = 0;
-            } // end for i
-            j  = -1; // since ly=ld=dt=0 forces a j++ call.
-            ly = 0; // zero old detector values
-            ld = 0;
-            dt = 0;
-            tt = track;
-        } // end if track != tt
-         itsHit->GetPositionG(gp[0],gp[1],gp[2],tof);
-         itsHit->GetDetectorID(id[0],id[1],id[2]);
-         gm->GtoL(id,gp,lp);
-        itsHit->GetMomentumG(p[0],p[1],p[2]);
-         if(!(id[0]==ly && id[1]==ld && id[2]==dt)){//if Not the same detector
-            j++; // incriment detector counter j.
-            ly = id[0];
-            ld = id[1];
-            dt = id[2];
-        } // end if id==idold
-        ht[j].lay = id[0];
-        ht[j].lad = id[1];
-        ht[j].det = id[2];
-        ht[j].xl += lp[0];
-        ht[j].yl += lp[1];
-        ht[j].zl += lp[2];
-        ht[j].xg  = gp[0];
-        ht[j].yg  = gp[1];
-        ht[j].zg  = gp[2];
-        px[j]    += p[0];
-        py[j]    += p[1];
-        pz[j]    += p[2];
-        nht[j]++;
-      } // end for h
-  } // end for t: track loop
-  return;
-}
-//______________________________________________________________________
-void HitsToClustAl(ClustAl_tl *trk,Int_t &ntrk,Int_t nt,TTree *TH,AliITS *ITS,
-                  Float_t fraction){
-
-   // Local variables
-   Int_t        nClustMin = 3,Iseed;
-//   Int_t        icount=0,icountMAX=100;
-   const Int_t  jmax = 50;
-   Int_t        t,nb,nh,h,i,ly=-1,ld=-1,dt=-1,tt=-1,j=0;
-   ClustAl_sl   ht[jmax];
-   Int_t        id[3],ie[3],nht[jmax],track;
-   Float_t      gp[3],lp[3],tof;
-   Float_t      ns;
-   TClonesArray *ITShits = ITS->Hits(),*Prt = gAlice->Particles();
-   AliITShit    *itsHit;
-   TParticle    *prt;
-   // initilize and creat AliITSgeom data
-   AliITSgeom   *gm  = ITS->GetITSgeom();
-
-   for(i=0;i<jmax;i++){// zero for for first track.
-       ht[i].xl  = 0.0;
-       ht[i].yl  = 0.0;
-       ht[i].zl  = 0.0;
-       ht[i].lay = ht[i].lad = ht[i].det = 0;
-       nht[i]    = 0;
-   } // end init to zero.
-
-   j = 0;
-   Iseed = time(0);
-   printf("HitsToclustAl: Iseed=%d ",Iseed);
-   gRandom->SetSeed(Iseed);
-   printf("gRandom->Rndm(1)=%f\n",gRandom->Rndm(1));
-   for(t=0;t<nt;t++){ // Loop over tracks.
-      if(fraction<gRandom->Rndm(1)) continue; // skip some tracks
-      gAlice->ResetHits();
-      nb = TH->GetEvent(t);
-      nh = ITShits->GetEntriesFast();
-      for(h=0;h<nh;h++){ // hits
-         itsHit = (AliITShit *)ITShits->UncheckedAt(h);
-        track  = itsHit->fTrack;
-         prt= (TParticle *)Prt->UncheckedAt(track);
-         if(prt->P()<MOMENTUMCUT) continue;
-        if(track != tt){ // if new track
-            if(j>nClustMin)if(tt<nt||PRIMARYONLY==0){ 
-//              if(icount>icountMAX) return;
-//              icount++;
-                 // if there are enought points save data
-                prt              = (TParticle *)Prt->UncheckedAt(tt);
-                trk[ntrk].track  = tt;
-                trk[ntrk].nclust = j+1;
-                ns               = 1./(Float_t (nht[0]));
-                trk[ntrk].px     = prt->Px();
-                trk[ntrk].py     = prt->Py();
-                trk[ntrk].pz     = prt->Pz();
-                trk[ntrk].pt     = prt->Pt();
-                trk[ntrk].p      = prt->P();
-                trk[ntrk].clust  = new ClustAl_sl[j+1];
-                for(i=0;i<=j;i++){ // loop over all detector hit.
-                    ns                     = 1./(Float_t(nht[i]));
-                    trk[ntrk].clust[i].lay = ie[0] = ht[i].lay;
-                    trk[ntrk].clust[i].lad = ie[1] = ht[i].lad;
-                    trk[ntrk].clust[i].det = ie[2] = ht[i].det;
-                    lp[0]                  = ht[i].xl*ns;
-                    lp[1]                  = ht[i].yl*ns;
-                    lp[2]                  = ht[i].zl*ns;
-                    // set in detector resolution on local coordiante system
-                    // z->zl rphi->xl yl-> thickness leave it alone
-                    lp[0] = lp[0]/resx[ie[0]-1];
-                    lp[0] = (Float_t)((Int_t)lp[0]);
-                    lp[0] = lp[0]*resx[ie[0]-1] + 0.5*resx[ie[0]-1];
-                    lp[2] = lp[2]/resz[ie[0]-1];
-                    lp[2] = (Float_t)((Int_t)lp[2]);
-                    lp[2] = lp[2]*resz[ie[0]-1] + 0.5*resz[ie[0]-1];
-                    //
-                    trk[ntrk].clust[i].xl  = lp[0];
-                    trk[ntrk].clust[i].yl  = lp[1];
-                    trk[ntrk].clust[i].zl  = lp[2];
-                } // end for i: detectors clust loop
-                ntrk++; // set up for next track
-            } // end if j>nClustMin
-            for(i=0;i<jmax;i++){ // zero out structure
-                ht[i].xl = ht[i].xg = 0.0;
-                ht[i].yl = ht[i].yg = 0.0;
-                ht[i].zl = ht[i].zg = 0.0;
-                ht[i].lay = ht[i].lad = ht[i].det = 0;
-                nht[i]   = 0;
-            } // end for i
-            j  = -1; // since ly=ld=dt=0 forces a j++ call.
-            ly = 0; // zero old detector values
-            ld = 0;
-            dt = 0;
-            tt = track;
-        } // end if track != tt
-         itsHit->GetPositionG(gp[0],gp[1],gp[2],tof);
-         itsHit->GetDetectorID(id[0],id[1],id[2]);
-         gm->GtoL(id,gp,lp);
-         if(!(id[0]==ly && id[1]==ld && id[2]==dt)){//if Not the same detector
-            j++; // incriment detector counter j.
-            ly = id[0];
-            ld = id[1];
-            dt = id[2];
-        } // end if id==idold
-        ht[j].lay = id[0];
-        ht[j].lad = id[1];
-        ht[j].det = id[2];
-        ht[j].xl += lp[0];
-        ht[j].yl += lp[1];
-        ht[j].zl += lp[2];
-        nht[j]++;
-      } // end for h
-  } // end for t: track loop
-  return;
-}
-//______________________________________________________________________
-void SetDetectorResolusion(Int_t l,Float_t *xl){
-    //set in detector resolution on local coordiante system
-    // z->zl rphi->xl yl-> thickness leave it alone
-    // changed from res to cell for detector simulations.
-
-    l--;
-    xl[0] = xl[0]/cellx[l];
-    xl[0] = (Float_t)((Int_t)xl[0]);
-    xl[0] = xl[0]*cellx[l] + 0.5*cellx[l];
-    xl[2] = xl[2]/cellz[l];
-    xl[2] = (Float_t)((Int_t)xl[2]);
-    xl[2] = xl[2]*cellz[l] + 0.5*cellz[l];
-    xl[1] = 0.0;
-}
-//______________________________________________________________________
-void FillAliITSAlignmentTrack(AliITSAlignmentTrack *trk,Int_t &ntrk,Int_t nt,
-                             TTree *TH,AliITS *ITS,Float_t fraction){
-    const Int_t iMAX=20,jMAX=50,nClustMin=3;
-    Int_t   ist[iMAX],imax;
-    Int_t   t,i,h,j,nh;
-    Int_t   lay,lad,det,index,track,Iseed;
-    Int_t   indexl[jMAX];
-    Float_t tof,xl[3],exl[3][3];
-    AliITSstatistics *sxl[jMAX],*szl[jMAX];
-    TParticle    *prt;
-    AliITShit    *itshit;
-    TClonesArray *ITShits = ITS->Hits();
-    AliITSgeom   *gm = ITS->GetITSgeom();
-    TClonesArray *Prt = gAlice->Particles();
-
-    printf("Entered FillAliITSAlignmentTrack\n");
-    for(i=0;i<jMAX;i++){
-       sxl[i] = new AliITSstatistics(2);
-       szl[i] = new AliITSstatistics(2);
-    }  // end for i
-    for(i=0;i<3;i++)for(j=0;j<3;j++) exl[i][j] = 0.0;
-    Iseed = time(0);
-    gRandom->SetSeed(Iseed);
-    printf("Entering track loop. Iseed=%d\n",Iseed);
-    for(t=0;t<nt;t++) if(fraction>=gRandom->Rndm(1)){
-//     gAlice->ResetHits();
-       TH->GetEvent(t);
-       nh = ITShits->GetEntriesFast();
-       i = 0;
-       ist[0] = 0;
-       printf("Entering hit loop for track=%d\n",t);
-       for(h=1;h<nh&&i<iMAX-1;h++) 
-           if(((AliITShit *)ITShits->UncheckedAt(ist[i]))->GetTrack() !=
-              ((AliITShit *)ITShits->UncheckedAt(h))->GetTrack()) ist[++i]=h;
-       ist[++i]=nh;
-       imax = i+1;
-       printf("entering loop i nh=%d imax=%d\n",nh,imax);
-       for(i=1;i<imax;i++){ // loop over tracks from primary track
-           printf("Getting hit ist[%d-1]=%d\n",i,ist[i-1]);
-           itshit = (AliITShit *)ITShits->UncheckedAt(ist[i-1]);
-           track  = itshit->GetTrack();
-           prt    = (TParticle *)Prt->UncheckedAt(track);
-           if(prt->P()<MOMENTUMCUT) continue;
-           printf("pass P cut P=%f\n",prt->P());
-           for(j=0;j<jMAX;j++){
-               sxl[j]->Reset();
-               szl[j]->Reset();
-           } // end for j
-           printf("Exiting the reset loop, ist[%d-1]=%d\n",i,ist[i-1]);
-           j = 0;
-           indexl[0]=((AliITShit *)ITShits->UncheckedAt(ist[i-1]))->
-                                                               GetDetector();
-           printf("looping over this track's hits ist[%d-1]=%d to ist[%d]=%d\n",i,ist[i-1],i,ist[i]);
-           for(h=ist[i-1];h<ist[i];h++){ // loop over hits
-               itshit = (AliITShit *)ITShits->UncheckedAt(h);
-               index = itshit->GetDetector();
-               if(indexl[j]!=index) j++;
-               indexl[j] = index;
-               itshit->GetPositionL(xl[0],xl[1],xl[2],tof);
-               sxl[j]->AddValue((Double_t) xl[0],1.0);
-               szl[j]->AddValue((Double_t) xl[2],1.0);
-           } // end for h
-           if(j<nClustMin) continue; // get next track
-           printf("Setting up tracks ntrk=%d #digits for track=%d\n",ntrk,j+1);
-           trk[ntrk].CreatePoints(j+1);
-           printf("CreatePoints with ntrk=%d and j+1=%d\n",ntrk,j+1);
-           trk[ntrk].SetTParticle(prt);
-           printf("SetTParticle\n");
-           trk[ntrk].SetTrackNumber(track);
-           printf("SetTrackNumber=%d\n",track);
-           for(h=0;h<=j;h++){
-               xl[0] = sxl[h]->GetMean();
-               xl[1] = 0.0;
-               xl[2] = szl[h]->GetMean();
-               gm->GetModuleId(indexl[h],lay,lad,det);
-               printf("setting detector resolusion for index[%d]=%d (%d,%d,%d)\n",h,indexl[h],lay,lad,det);
-               SetDetectorResolusion(lay,xl);
-               printf("Detector resolusion set\n");
-               exl[0][0] = 1./(resx[lay-1]*resx[lay-1]);
-               exl[1][1] = 12./(celly[lay-1]*celly[lay-1]);
-               exl[2][2] = 1./(resz[lay-1]*resz[lay-1]);
-               printf("Adding point to trk[%d] indexl[%d]=%d\n",ntrk,h,indexl[h]);
-               trk[ntrk].AddPointLastL(indexl[h],(Float_t *)xl,
-                                                            (Float_t **)exl);
-           } // end for h
-           ntrk++; // set up for next track
-       } // end for i
-    } // end for t // end if(fraction>=gRandom->Rndm(1))
-    for(i=0;i<jMAX;i++){
-       delete sxl[i];
-       delete szl[i];
-    }  // end for i
-}
-//______________________________________________________________________
-void FillGlobalPositions(ClustAl_tl *trk,Int_t ntrk,AliITSgeom *g){
-    Int_t   i,j,id[3];
-    Float_t lx[3],gx[3];
-
-    for(i=0;i<ntrk;i++) for(j=0;j<trk[i].nclust;j++){
-       lx[0] = trk[i].clust[j].xl;
-       lx[1] = trk[i].clust[j].yl;
-       lx[2] = trk[i].clust[j].zl;
-       id[0] = trk[i].clust[j].lay;
-       id[1] = trk[i].clust[j].lad;
-       id[2] = trk[i].clust[j].det;
-       g->LtoG(id,lx,gx);
-       trk[i].clust[j].xg = gx[0];
-       trk[i].clust[j].yg = gx[1];
-       trk[i].clust[j].zg = gx[2];
-    } // end for i,j
-    return;
-}
-//______________________________________________________________________
-void PlotGeomChanges(AliITSgeom *gt,AliITSgeom *gc,TFile *Hfile,Float_t *Rdta){
-    Int_t    ly,ld,dt;
-    Float_t  x,y,z;
-    Double_t xd,yd,zd;
-    Double_t rt,rc,phit,phic,zt,zc;
-    Float_t  A0,B0,C0,A1,B1,C1;
-    Double_t dr,drphi,dz,dphi;
-    Double_t drmi,drma,drphimi,drphima,dzmi,dzma;
-
-    AliITSstatistics *R    = new AliITSstatistics(2);
-    AliITSstatistics *RPhi = new AliITSstatistics(2);
-    AliITSstatistics *Z    = new AliITSstatistics(2);
-    AliITSstatistics *A    = new AliITSstatistics(2);
-    AliITSstatistics *B    = new AliITSstatistics(2);
-    AliITSstatistics *C    = new AliITSstatistics(2);
-    TH1F *Gr  = new TH1F("Gr","Radial Displacement (cm)",500,-1.0,1.0);
-    drmi = -1.0;
-    drma = +1.0;
-    Gr->Sumw2();
-    Gr->SetXTitle("Displacement (cm)");
-    TH1F *Grphi  = new TH1F("Grphi","RPhi Displacement (cm)",500,-1.0,1.0);
-    drphimi = -1.0;
-    drphima = +1.0;
-    Grphi->Sumw2();
-    Grphi->SetXTitle("Displacement (cm)");
-    TH1F *Gz  = new TH1F("Gz","Z Displacement (cm)",500,-1.0,1.0);
-    dzmi = -1.0;
-    dzma = +1.0;
-    Gz->Sumw2();
-    Gz->SetXTitle("Displacement (cm)");
-
-    Float_t weight=1.0;
-    for(ly=1;ly<=gt->GetNlayers();ly++)for(ld=1;ld<=gt->GetNladders(ly);ld++)
-    for(dt=1;dt<=gt->GetNdetectors(ly);dt++){
-       gt->GetTrans(ly,ld,dt,x,y,z);xd=x;yd=y;zd=z;
-       gt->GetAngles(ly,ld,dt,A0,B0,C0);
-       rt    = TMath::Hypot(yd,xd);
-       phit  = TMath::ATan2(yd,xd);
-       if(phit<0.0) phit += 2.0*TMath::Pi();
-       zt    = zd;
-       gc->GetTrans(ly,ld,dt,x,y,z);xd=x;yd=y;zd=z;
-       gc->GetAngles(ly,ld,dt,A1,B1,C1);
-       rc    = TMath::Hypot(yd,xd);
-       phic  = TMath::ATan2(yd,xd);
-       if(phic<0.0) phic += 2.0*TMath::Pi();
-       zc    = zd;
-       dr    = rt-rc;
-       dphi  = phit - phic;
-       if(dphi>2.0*TMath::Pi()) dphi -= 2.0*TMath::Pi();
-       drphi = 0.5*(rt + rc)*dphi; // change in phi as measured in rphi
-       dz    = zt - zc;
-       Gr->Fill(dr,1.0);
-       Grphi->Fill(drphi,1.0);
-       Gz->Fill(dz,1.0);
-       if(dr>=drmi&&dr<drma) R->AddValue(dr,1.0);
-       if(drphi>=drphimi&&drphi<drphima) RPhi->AddValue(drphi,1.0);
-       if(dz>=dzmi&&dz<dzma) Z->AddValue(dz,1.0);
-       A->AddValue(A0-A1,weight);
-       B->AddValue(B0-B1,weight);
-       C->AddValue(C0-C1,weight);
-    } // end for ly,ld,dt
-    Hfile->Write();
-    Rdta[0] = R->GetRMS();
-    Rdta[1] = RPhi->GetRMS();
-    Rdta[2] = Z->GetRMS();
-    Rdta[3] = A->GetRMS();
-    Rdta[4] = B->GetRMS();
-    Rdta[5] = C->GetRMS();
-    printf("PlotGeomChanges: RMS(r) = %f RMS(rphi) = %f RMS(z) = %f\n",
-          R->GetRMS(),RPhi->GetRMS(),Z->GetRMS());
-    delete Gr;
-    delete Grphi;
-    delete Gz;
-    delete R;
-    delete RPhi;
-    delete Z;
-    delete A;
-    delete B;
-    delete C;
-    return;
-}
-//______________________________________________________________________
-Int_t FitTrackToLine(ClustAl_tl &trk){
-   // Local Variables
-   Int_t   i,l,j=0;
-   Double_t x,y,z,wx,wy,dx,dy;
-   Double_t a,b,c,d;
-   Double_t xb,sxb,zxb;
-   Double_t yb,syb,zyb;
-   Double_t xzb,yzb,z2xb,z2yb;
-   Double_t sum,phi=0.0,b0,d0;
-
-   trk.qual = -1.0;
-   if(trk.nclust<3) return -1;
-
-   b  = (trk.clust[0].xg-trk.clust[trk.nclust].xg)/
-        (trk.clust[0].zg-trk.clust[trk.nclust].zg);
-   d  = (trk.clust[0].yg-trk.clust[trk.nclust].yg)/
-        (trk.clust[0].zg-trk.clust[trk.nclust].zg);
-   do{
-     b0 = b;
-     d0 = d;
-     xb  = 0.0; sxb = 0.0; zxb  = 0.0;
-     yb  = 0.0; syb = 0.0; zyb  = 0.0;
-     xzb = 0.0; yzb = 0.0; z2xb = 0.0; z2yb = 0.0;
-     for(i=0;i<trk.nclust;i++){
-      l    = trk.clust[i].lay - 1;  // zero based
-      x    = trk.clust[i].xg;
-      y    = trk.clust[i].yg;
-      z    = trk.clust[i].zg;
-      phi  = atan2(y,x);
-      wx   = 1.0/(resx[l]*resx[l]*cos(phi)*cos(phi)+b0*b0*resz[l]*resz[l]);// 1.0/rms^2
-      wy   = 1.0/(resx[l]*resx[l]*sin(phi)*sin(phi)+d0*d0*resz[l]*resz[l]);// 1.0/rms^2
-      xb   += x*wx;
-      sxb  += wx;
-      zxb  += z*wx;
-      yb   += y*wy;
-      syb  += wy;
-      zyb  += z*wy;
-      xzb  += x*z*wx;
-      yzb  += y*z*wy;
-      z2xb += z*z*wx;
-      z2yb += z*z*wy;
-     } // end for i
-     dx = zxb*zxb - z2xb*sxb;
-     if(dx!=0.0){
-        a  = zxb*xzb - xb*z2xb;
-        a  = a/dx;
-        b  = xb*zxb  - sxb*xzb;
-        b  = b/dx;
-     }else{
-         printf("FitTrackToLine: Error dx=0 track=%d zxb=%f z2xb=%f sxb=%f\n",
-                trk.track,zxb,z2xb,sxb);
-         a = 0.0;
-         b = 1E10;
-         j++;
-         if(j>100) return(0);
-     } // end if dx!=0.0
-     dy = zyb*zyb - z2yb*syb;
-     if(dy!=0.0){
-        c  = zyb*yzb - yb*z2yb;
-        c  = c/dy;
-        d  = yb*zyb  - syb*yzb;
-        d  = d/dy;
-     }else{
-         printf("FitTrackToLine: Error dy=0 track=%d zyb=%f z2yb=%f syb=%f\n",
-                trk.track,zyb,z2yb,syb);
-         c = 0.0;
-         d = 1.E10;
-         j++;
-         if(j>100) return(0);
-     } // end if dy!=0.0
-   } while(fabs(b0-b)<1.E-5 && fabs(d0-d)<1.E-5);
-   trk.a = a;
-   trk.b = b;
-   trk.c = c;
-   trk.d = d;
-   sum   = 0.0;
-   for(i=0;i<trk.nclust;i++){
-      l  = trk.clust[i].lay - 1;  // zero based
-      x  = trk.clust[i].xg;
-      y  = trk.clust[i].yg;
-      z  = trk.clust[i].zg;
-      x  = x - (a + b*z);  // x=a+bz
-      y  = y - (c + d*z);  // y=c+dz
-      wx =  resx[l]*resx[l]*cos(phi)*cos(phi) + b*b * resz[l]*resz[l];
-      wy =  resx[l]*resx[l]*sin(phi)*sin(phi) + d*d * resz[l]*resz[l];
-      if(wx==0.0||wy==0.0) {
-         printf("FitTrackToLine: error wx=%f wy=%f trk.clust[%d].lay=%d\n",
-                 wx,wy,i,trk.clust[i].lay);
-      }// end if sxz or syz == 0.
-      sum += x*x/wx + y*y/wy;
-   } // end for i #2
-   trk.qual = sum;
-   if(trk.nclust>2) trk.qual /= Double_t (trk.nclust-2);
-   // per degree of freedom 2 in xz 2 in yz.
-   return 0;
-}
-//______________________________________________________________________
-void LtoLline(const Int_t *id1,const Int_t *id2,
-             Double_t a1,Double_t b1,Double_t c1,Double_t d1,
-             Double_t &a2,Double_t &b2,Double_t &c2,Double_t &d2,
-             AliITSgeom *gm){
-    Double_t x11[3],x12[3],x21[3],x22[3],h;
-
-    x11[0] = a1;
-    x11[1] = 0.0;
-    x11[2] = c1;
-    x12[0] = a1+b1;
-    x12[1] = 1.0;
-    x12[2] = c1+d1;
-    gm->LtoL(id1,id2,x11,x21);
-    gm->LtoL(id1,id2,x12,x22);
-    h = x21[1] - x22[1];
-    if(h!=0.0){
-       b2 = (x21[0] - x22[0])/h;
-       d2 = (x21[2] - x22[2])/h;
-       a2 = x21[0] - b2*x21[1];
-       c2 = x21[2] - d2*x21[1];
-    }else{
-       printf("LtoLline: error line in plane of detector\n");
-    } // end if h!=0
-    return;
-}
-//______________________________________________________________________
-Int_t FitTrackToLineL(ClustAl_tl &trk,AliITSgeom *gm){
-   // Local Variables
-   Int_t   i,id0[3],id[3];
-   Double_t wx/*,wy*/,wz;
-   Double_t a,b,c,d;
-   Double_t xg[3],xl[3],x2g[3],x2l[3];
-   AliITSstatistics2 *Fx  = new AliITSstatistics2(2);
-   AliITSstatistics2 *Fz  = new AliITSstatistics2(2);
-
-   trk.qual = -1.0;
-   if(trk.nclust<3) return -1;
-
-   Int_t Npts = trk.nclust;
-   id0[0] = trk.clust[0].lay;
-   id0[1] = trk.clust[0].lad;
-   id0[2] = trk.clust[0].det;
-   for(i=0;i<Npts;i++){
-       id[0] = trk.clust[i].lay;
-       id[1] = trk.clust[i].lad;
-       id[2] = trk.clust[i].det;
-       xg[0] = trk.clust[i].xg;
-       xg[1] = trk.clust[i].yg;
-       xg[2] = trk.clust[i].zg;
-       gm->GtoL(id0,xg,xl);
-       wx = 1.0/(resx[id[0]-1]*resx[id[0]-1]);
-       wz = 1.0/(resz[id[0]-1]*resz[id[0]-1]);
-       Fx->AddValue(xl[0],xl[1],wx);
-       Fz->AddValue(xl[2],xl[1],wz);
-   } // end for i
-   trk.qual  = Fx->FitToLine(a,b);
-   trk.qual += Fz->FitToLine(c,d);
-   trk.a0 = a;
-   trk.b0 = b;
-   trk.c0 = c;
-   trk.d0 = d;
-   xl[0]  = a;
-   xl[1]  = 0.0;
-   xl[2]  = c;
-   x2l[0] = a+b;
-   x2l[1] = 1.0;
-   x2l[2] = c+d;
-   gm->LtoG(id0,xl,xg);
-   gm->LtoG(id0,x2l,x2g);
-   c = xg[2] - x2g[2];
-   if(c!=0.0){
-       b = (xg[0] - x2g[0])/c;
-       d = (xg[1] - x2g[1])/c;
-       a = xg[0] - b*xg[2];
-       c = xg[1] - d*xg[2];
-       trk.a = a;
-       trk.b = b;
-       trk.c = c;
-       trk.d = d;
-   }else{
-       return -1;
-   }// end if c!=0.0
-   return 0;
-}
-//______________________________________________________________________
-Int_t FindCircleCenter(Double_t &x0,Double_t &y0,Double_t x1,Double_t y1,
-                      Double_t  x2,Double_t  y2,Double_t x3,Double_t y3){
-////////////////////////////////////////////////////////////////////////
-//     This was derived as folows. Given three non-linear points find
-// the circle that is therefor defined by those three non-linar points.
-// Assume that the circle is centers at x0,y0 and has a radous R. Then
-// (1) R^2 = (x1-x0)^2 + (y1-y0)^2
-// (2) R^2 = (x2-x0)^2 + (y2-y0)^2
-// (3) R^2 = (x3-x0)^2 + (y3-y0)^2.
-// Now consider the two equations derived from the above
-// (1) - (2) = x1^2 - x2^2 -2x0(x1-x2) + y1^2 - y2Y2 -2Y0(y1-y2) = 0
-// (1) - (3) = x1^2 - x3^2 -2x0(x1-x3) + y1^2 - y3Y2 -2Y0(y1-y3) = 0
-// solving these two equations for x0 and y0 gives
-// x0 = +{(y1-y2)*(y1-y3)*(y1-y3)+x1*x1*(y1-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2)}/2d
-// y0 = -{(x1-x2)*(x1-x3)*(x1-x3)+y1*y1*(x1-x3)+y2*y2*(x3-x1)+y3*y3*(x1-x2)}/2d
-// with d = (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2)
-////////////////////////////////////////////////////////////////////////
-    Double_t d;
-
-    d = (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2);
-    if(d==0.0) return 0;  // fits to a line!
-
-    x0 = (y1-y2)*(y1-y3)*(y1-y3)+x1*x1*(y1-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2);
-    y0 = (x1-x2)*(x1-x3)*(x1-x3)+y1*y1*(x1-x3)+y2*y2*(x3-x1)+y3*y3*(x1-x2);
-    x0 = +0.5*x0/d;
-    y0 = -0.5*y0/d;
-
-    return 1;
-}
-//______________________________________________________________________
-void FitAllTracks(ClustAl_tl *trk,Int_t ntrk,Float_t *v0,AliITSgeom *gm,
-                 const char *sfile,TFile *Hfile,Float_t *Fdta,Int_t *Ndta){
-   // Local Variables
-   Int_t    i,j,k,id0[3],id[3];
-   Double_t xm,ym,zm,dt,ad/*,rh,rp,rphih,rphip*/;
-   Double_t tp,tx,rx,rz;
-   Double_t a,b,c,d,qual;
-   Double_t xl[3],xg[3];
-   Double_t trkqualMAX = 50.0;
-   Bool_t   printit = kFALSE;
-   char     filename[80],hid[10];
-   Int_t    Nqualmax = 10,NqualgtMAX=0;
-   Int_t    Nqualless[10];
-// Distrobution statisitics
-   AliITSstatistics *Sad    = new AliITSstatistics(4);
-   AliITSstatistics *SrxzL  = new AliITSstatistics(4);
-   AliITSstatistics *SrxL[6];
-   AliITSstatistics *SrzL[6];
-   for(i=0;i<6;i++) {
-       SrxL[i] = new AliITSstatistics(4);
-       SrzL[i] = new AliITSstatistics(4);
-   } // end for i
-
-   for(i=0;i<Nqualmax;i++) Nqualless[i] = 0;
-
-// Setup histograms
-
-   TH2F *Ttqp  = new TH2F("Ttqp","Track quality vs. momentum",
-                        500,0.0,50.0,500,0.0,10.0);
-   Ttqp->Sumw2();
-   Ttqp->SetXTitle("Chi^2/degree freedom");
-   Ttqp->SetYTitle("Momentum GeV/c");
-   TH2F *Tdttq = new TH2F("Tdttq",
-                         "Distance from true vertex vs. track quality",
-                         500,0.0,0.10,100,0.0,50.0);
-   Tdttq->Sumw2();
-   Tdttq->SetXTitle("Distanct to true vertex (cm)");
-   Tdttq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Tadtq = new TH2F("Tadtq",
-                     "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. track quality",
-                         500,-0.03,0.03,500,0.0,50.0);
-   Tadtq->Sumw2();
-   Tadtq->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Tadtq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Taddt = new TH2F("Taddt",
-         "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. Distance from true vertex",
-                         500,-0.03,0.03,200,0.0,0.10);
-   Taddt->Sumw2();
-   Taddt->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Taddt->SetYTitle("distance to true vertex (cm)");
-   TH2F *Tadp = new TH2F("Tadp",
-         "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. momentum",
-                         500,-0.03,0.03,200,0.0,10.0);
-   Tadp->Sumw2();
-   Tadp->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Tadp->SetYTitle("Momentum (GeV/c)");
-
-   TH2F *TrxrzL[7];
-//   TH2F *TrrprzG[7];
-   for(i=0;i<7;i++){
-       sprintf(filename,"Layer %1.1d: xi-x(i) vs. zi-z(i) local",i+1);
-       if(i==6) sprintf(filename,"Sum of all layers: xi-x(i) "
-                       "vs. zi-z(i) local");
-       sprintf(hid,"TrxrzL%1.1d",i+1);
-       TrxrzL[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
-       TrxrzL[i]->Sumw2();
-       TrxrzL[i]->SetXTitle("xi-x(i) (cm) local");
-       TrxrzL[i]->SetYTitle("zi-z(i) (cm) local");
-//
-//       sprintf(filename,"Layer %1.1d: rPhii-rPhi(i) vs. zi-z(i) local",i+1);
-//       if(i==6) sprintf(filename,"Sum of all layers: rPhii-rPhi(i) "
-//                     "vs. zi-z(i) local");
-//       sprintf(hid,"TrrprzG%1.1d",i+1);
-//       TrrprzG[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
-//       TrrprzG[i]->Sumw2();
-//       TrrprzG[i]->SetXTitle("rPhii-rPhi(i) (cm) local");
-//       TrrprzG[i]->SetYTitle("zi-z(i) (cm) local");
-   } // end for i
-
-   TH1D *Tptqp[10];
-   TH1D *Ttqall = new TH1D("Ttqall","Track quality for all tracks fit",
-                          500,0.0,1000.0);
-   Ttqall->Sumw2();
-   Ttqall->SetXTitle("Track quality: Chi squared per degree freedom");
-
-// Fit each track and fill histograms and the like.
-
-   for(i=0;i<ntrk;i++){
-       if(trk[i].p < MOMENTUMCUT) continue;
-      if(FitTrackToLineL(trk[i],gm)!=0) continue;  // fit track to a line
-      qual = trk[i].qual;
-      if(qual<0.0) continue;
-      // Fill histograms and statistics before cource chi squared cut
-      Ttqall->Fill(qual,1.0);
-      for(j=0;j<Nqualmax;j++) if(qual<(Double_t) j+1) Nqualless[j]++;
-//      if(qual>trkqualMAX) {
-//       NqualgtMAX++;
-//       continue;
-//      } // end if iqual>trkqualMAX=50.0
-      a = trk[i].a;
-      b = trk[i].b;
-      c = trk[i].c;
-      d = trk[i].d;
-      zm  = b*(v0[0]-a) + d*(v0[1]-c);
-      dt  = (b*b + d*d);
-      if(dt!=0.0) zm /= dt;
-      else{
-         printf("FitAllTracks: trk[%d].b=%f trk[%d].d=%f\n",i,b,i,d);
-         zm=0.0;
-      }// end if else dt!=0.0
-      xm  = a + b*zm;
-      ym  = c + d*zm;
-      dt   = sqrt((xm-v0[0])*(xm-v0[0]) + (ym-v0[1])*(ym-v0[1]) +
-                (zm-v0[2])*(zm-v0[2]));
-      tx  = atan(sqrt(b*b + d*d));
-      tp  = atan2(trk[i].pt,fabs(trk[i].pz));
-      ad  = tx-tp;
-      // Fill histograms and statistics
-      Tadtq->Fill(ad,qual,1.0);
-      Sad->AddValue(ad,1.0);
-      Taddt->Fill(ad,dt,1.0);
-      Tadp->Fill(ad,trk[i].p,1.0);
-      Tdttq->Fill(dt,qual,1.0);
-      Ttqp->Fill(qual,trk[i].p,1.0);
-      //
-      id0[0] = trk[i].clust[0].lay;
-      id0[1] = trk[i].clust[0].lad;
-      id0[2] = trk[i].clust[0].det;
-      xg[0] = trk[i].clust[0].xg;
-      xg[1] = trk[i].clust[0].yg;
-      xg[2] = trk[i].clust[0].zg;
-      gm->GtoL(id0,xg,xl);
-      rx = xl[0] - trk[i].a0 - trk[i].b0 * xl[1];
-      rz = xl[2] - trk[i].c0 - trk[i].d0 * xl[1];
-      SrxL[id0[0]-1]->AddValue(rx,1.0);
-      SrzL[id0[0]-1]->AddValue(rz,1.0);
-      SrxzL->AddValue(rx,1.0);
-      SrxzL->AddValue(rz,1.0);
-      TrxrzL[6]->Fill(rx,rz,1.0);
-      TrxrzL[id0[0]-1]->Fill(rx,rz,1.0);
-      for(j=1;j<trk[i].nclust;j++){
-         id[0] = trk[i].clust[j].lay;
-         id[1] = trk[i].clust[j].lad;
-         id[2] = trk[i].clust[j].det;
-         LtoLline(id0,id,trk[i].a0,trk[i].b0,trk[i].c0,trk[i].d0,a,b,c,d,gm);
-         xg[0] = trk[i].clust[j].xg;
-         xg[1] = trk[i].clust[j].yg;
-         xg[2] = trk[i].clust[j].zg;
-         gm->GtoL(id,xg,xl);
-         rx = xl[0] - a - b * xl[1];
-         rz = xl[2] - c - d * xl[1];
-         SrxL[id[0]-1]->AddValue(rx,1.0);
-         SrzL[id[0]-1]->AddValue(rz,1.0);
-         SrxzL->AddValue(rx,1.0);
-         SrxzL->AddValue(rz,1.0);
-         TrxrzL[6]->Fill(rx,rz,1.0);
-         TrxrzL[id[0]-1]->Fill(rx,rz,1.0);
-      } // end for j
-   } // end for i
-
-// Write out information
-
-   for(i=0;i<6;i++){
-       Fdta[4*i+0] = SrxL[i]->GetRMS();
-       Fdta[4*i+1] = SrxL[i]->GetErrorRMS();
-       Fdta[4*i+2] = SrzL[i]->GetRMS();
-       Fdta[4*i+3] = SrzL[i]->GetErrorRMS();
-   } // end for i
-   for(i=0;i<10;i++) Ndta[i] = Nqualless[i];
-   Ndta[10] = Sad->GetN();
-   Ndta[11] = NqualgtMAX;
-
-   printf("FitAllTracks: %d tracks cut leaving %d, with a chi squared <%f\n",
-         NqualgtMAX,Sad->GetN(),trkqualMAX);
-   printf("The number of tracks with chi squared <1");
-   for(i=1;i<Nqualmax;i++) printf(",%d",i+1);
-   printf("\nFitAllTracks: %d",Nqualless[0]);
-   for(i=1;i<Nqualmax;i++) printf(",%d",Nqualless[i]);
-   printf("\n");
-   printf("FitAllTracks: RMSs of ad rxz = %e+-%e %e+-%e\n",
-         Sad->GetRMS(),Sad->GetErrorRMS(),
-         SrxzL->GetRMS(),SrxzL->GetErrorRMS());
-
-   printf("FitAllTracks: Residuals by layer x=-rphi z ");
-   for(i=0;i<6;i++) printf(":%d:%e+-%e %e+-%e ",i+1,
-                          SrxL[i]->GetRMS(),SrxL[i]->GetErrorRMS(),
-                          SrzL[i]->GetRMS(),SrzL[i]->GetErrorRMS());
-   printf("\n");
-
-// Setup and fill projections
-
-   for(i=0;i<10;i++){
-       xm = 0.5*Double_t(i);
-       ym = xm+0.5;
-       j  = Ttqp->GetYaxis()->FindBin(xm);
-       k  = Ttqp->GetYaxis()->FindBin(ym);
-       sprintf(filename,"Track Quality for %3.1f<p<%3.1f",xm,ym);
-       Tptqp[i] = Ttqp->ProjectionX(filename,j,k,"E");
-       Tptqp[i]->SetXTitle("Track Quality (Chi squared/df)");
-   } // end for i
-
-   TH1D *Ttq  = Ttqp->ProjectionX("Ttq",0,Ttqp->GetNbinsY()+1,"E");
-   Ttq->SetXTitle("Chi^2/degree freedom");
-   TH1D *Tp   = Ttqp->ProjectionY("Tp", 0,Ttqp->GetNbinsX()+1,"E");
-   Tp->SetXTitle("Momentum GeV/c");
-   TH1D *Tdt  = Tdttq->ProjectionX("Tdt",0,Tdttq->GetNbinsY()+1,"E");
-   Tdt->SetXTitle("Distanct to true vertex (cm)");
-   TH1D *Tad  = Tadtq->ProjectionX("Tad",0,Tadtq->GetNbinsY()+1,"E");
-   Tad->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   TH1D *TrxL[7],*TrzL[7];
-//   TH1D *TrrphiG[7],*TrzG[7];
-   for(i=0;i<7;i++){
-       sprintf(hid,"TrxL%1.1d",i+1);
-       TrxL[i] = TrxrzL[i]->ProjectionX(hid,0,TrxrzL[i]->GetNbinsY()+1,"E");
-       TrxL[i]->SetXTitle("xi-x(i) (cm)");
-       sprintf(hid,"TrzL%1.1d",i+1);
-       TrzL[i] = TrxrzL[i]->ProjectionY(hid,0,TrxrzL[i]->GetNbinsX()+1,"E");
-       TrzL[i]->SetXTitle("zi-z(i) (cm)");
-//
-//       sprintf(hid,"TrrphiG%1.1d",i+1);
-//       TrrphiG[i] = TrrprzG[i]->ProjectionX(hid,0,
-//                                           TrrprzG[i]->GetNbinsY()+1,"E");
-//       TrrphiG[i]->SetXTitle("rPhii-rPhi(i) (cm)");
-//       sprintf(hid,"TrzG%1.1d",i+1);
-//       TrzG[i] = TrrprzG[i]->ProjectionY(hid,0,
-//                                        TrrprzG[i]->GetNbinsX()+1,"E");
-//       TrzG[i]->SetXTitle("zGi-zG(i) (cm)");
-   } // end for i
-
-
-   Hfile->Write();
-
-   if(printit){
-       TCanvas *c0 = new TCanvas("c0","Track quality distribution",
-                                500,100,600,700);
-       Ttqall->Draw();
-       sprintf(filename,"%s_T_tqall.ps",sfile);
-       if(printit) c0->Print(filename);
-       Ttqp->Draw("COL");
-       sprintf(filename,"%s_T_tq_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Ttq->Draw();
-       sprintf(filename,"%s_T_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tp->Draw();
-       sprintf(filename,"%s_T_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tdttq->Draw("COL");
-       sprintf(filename,"%s_T_dt_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tdt->Draw();
-       sprintf(filename,"%s_T_dt.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tadtq->Draw("COL");
-       sprintf(filename,"%s_T_ad_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tad->Draw();
-       sprintf(filename,"%s_T_ad.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tadp->Draw("COL");
-       sprintf(filename,"%s_T_ad_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Taddt->Draw("COL");
-       sprintf(filename,"%s_T_ad_dt.ps",sfile);
-       if(printit) c0->Print(filename);
-       for(i=0;i<7;i++){
-          TrxrzL[i]->Draw("COL");
-          sprintf(filename,"%s_T%1.1d_rx_rz.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-          TrxL[i]->Draw();
-          sprintf(filename,"%s_T%1.1d_rx.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-          TrzL[i]->Draw();
-          sprintf(filename,"%s_T%1.1d_rz.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-       } // end for i
-       for(i=0;i<10;i++){
-          Tptqp[i]->Draw();
-          sprintf(filename,"%s_T_tq_p%1.1d.ps",sfile,i);
-          if(printit) c0->Print(filename);
-       } // end for i
-   } // end if printit
-
-// Delet allocated stuff.
-   for(i=0;i<7;i++) {
-       delete TrxL[i];
-       delete TrzL[i];
-//       delete TrrphiG[i];
-//       delete TrzG[i];
-   } // end for i
-   for(i=0;i<10;i++) delete Tptqp[i];
-   delete Ttqp;
-   delete Tdttq;
-   delete Tadtq;
-   delete Taddt;
-   delete Tadp;
-   for(i=0;i<7;i++) delete TrxrzL[i];
-//   for(i=0;i<7;i++) delete TrrprzG[i];
-   delete Ttqall;
-   delete Ttq;
-   delete Tp;
-   delete Tdt;
-   delete Tad;
-//   printf("finished with track fitting\n");
-   delete Sad;
-   delete SrxzL;
-   for(i=0;i<6;i++) delete SrxL[i];
-   for(i=0;i<6;i++) delete SrzL[i];
-   return;
-}
-//______________________________________________________________________
-void FitAllTracksG(ClustAl_tl *trk,Int_t ntrk,Float_t *v0,AliITSgeom *gm,
-                 const char *sfile,TFile *Hfile){
-   // Local Variables
-   Int_t    i,j,k,lay,lad,det;
-   Double_t xm,ym,zm,dt,ad,rh,rp,rphih,rphip;
-   Double_t tp,tx,rx,rz/*,rr*/,rrphi;
-   Double_t xg[3],xl[3],xp,yp,xh,yh,a,b,c,d,x1[3],x2[3],qual;
-   Double_t trkqualMAX = 50.0;
-   Bool_t   printit = kFALSE;
-   char     filename[80],hid[10];
-   Int_t    Nqualmax = 10,NqualgtMAX=0;
-   Int_t    Nqualless[10];
-// Distrobution statisitics
-   AliITSstatistics *Sad    = new AliITSstatistics(4);
-   AliITSstatistics *SrxzL  = new AliITSstatistics(4);
-   AliITSstatistics *SrxL[6];
-   AliITSstatistics *SrzL[6];
-   for(i=0;i<6;i++) {
-       SrxL[i] = new AliITSstatistics(4);
-       SrzL[i] = new AliITSstatistics(4);
-   } // end for i
-
-   for(i=0;i<Nqualmax;i++) Nqualless[i] = 0;
-
-// Setup histograms
-
-   TH2F *Ttqp  = new TH2F("Ttqp","Track quality vs. momentum",
-                        500,0.0,50.0,500,0.0,10.0);
-   Ttqp->Sumw2();
-   Ttqp->SetXTitle("Chi^2/degree freedom");
-   Ttqp->SetYTitle("Momentum GeV/c");
-   TH2F *Tdttq = new TH2F("Tdttq",
-                         "Distance from true vertex vs. track quality",
-                         500,0.0,0.10,100,0.0,50.0);
-   Tdttq->Sumw2();
-   Tdttq->SetXTitle("Distanct to true vertex (cm)");
-   Tdttq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Tadtq = new TH2F("Tadtq",
-                     "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. track quality",
-                         500,-0.03,0.03,500,0.0,50.0);
-   Tadtq->Sumw2();
-   Tadtq->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Tadtq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Taddt = new TH2F("Taddt",
-         "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. Distance from true vertex",
-                         500,-0.03,0.03,200,0.0,0.10);
-   Taddt->Sumw2();
-   Taddt->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Taddt->SetYTitle("distance to true vertex (cm)");
-   TH2F *Tadp = new TH2F("Tadp",
-         "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. momentum",
-                         500,-0.03,0.03,200,0.0,10.0);
-   Tadp->Sumw2();
-   Tadp->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   Tadp->SetYTitle("Momentum (GeV/c)");
-
-   TH2F *TrxrzL[7],*TrrprzG[7];
-   for(i=0;i<7;i++){
-       sprintf(filename,"Layer %1.1d: xi-x(i) vs. zi-z(i) local",i+1);
-       if(i==6) sprintf(filename,"Sum of all layers: xi-x(i) "
-                       "vs. zi-z(i) local");
-       sprintf(hid,"TrxrzL%1.1d",i+1);
-       TrxrzL[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
-       TrxrzL[i]->Sumw2();
-       TrxrzL[i]->SetXTitle("xi-x(i) (cm) local");
-       TrxrzL[i]->SetYTitle("zi-z(i) (cm) local");
-//
-       sprintf(filename,"Layer %1.1d: rPhii-rPhi(i) vs. zi-z(i) local",i+1);
-       if(i==6) sprintf(filename,"Sum of all layers: rPhii-rPhi(i) "
-                       "vs. zi-z(i) local");
-       sprintf(hid,"TrrprzG%1.1d",i+1);
-       TrrprzG[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
-       TrrprzG[i]->Sumw2();
-       TrrprzG[i]->SetXTitle("rPhii-rPhi(i) (cm) local");
-       TrrprzG[i]->SetYTitle("zi-z(i) (cm) local");
-   } // end for i
-
-   TH1D *Tptqp[10];
-   TH1D *Ttqall = new TH1D("Ttqall","Track quality for all tracks fit",
-                          500,0.0,1000.0);
-   Ttqall->Sumw2();
-   Ttqall->SetXTitle("Track quality: Chi squared per degree freedom");
-
-// Fit each track and fill histograms and the like.
-
-   for(i=0;i<ntrk;i++){
-       if(FitTrackToLine(trk[i])!=0) continue;  // fit track to a line
-      qual = trk[i].qual;
-      if(qual<0.0) continue;
-      // Fill histograms and statistics before cource chi squared cut
-      Ttqall->Fill(qual,1.0);
-      for(j=0;j<Nqualmax;j++) if(qual<(Double_t) j+1) Nqualless[j]++;
-      if(qual>trkqualMAX) {
-         NqualgtMAX++;
-         continue;
-      } // end if iqual>trkqualMAX=50.0
-      a = trk[i].a;
-      b = trk[i].b;
-      c = trk[i].c;
-      d = trk[i].d;
-      zm  = b*(v0[0]-a) + d*(v0[1]-c);
-      dt  = (b*b + d*d);
-      if(dt!=0.0) zm /= dt;
-      else{
-         printf("FitAllTracks: trk[%d].b=%f trk[%d].d=%f\n",i,b,i,d);
-         zm=0.0;
-      }// end if else dt!=0.0
-      xm  = a + b*zm;
-      ym  = c + d*zm;
-      dt   = sqrt((xm-v0[0])*(xm-v0[0]) + (ym-v0[1])*(ym-v0[1]) +
-                (zm-v0[2])*(zm-v0[2]));
-      tx  = atan(sqrt(b*b + d*d));
-      tp  = atan2(trk[i].pt,fabs(trk[i].pz));
-      ad  = tx-tp;
-      // Fill histograms and statistics
-      Tadtq->Fill(ad,qual,1.0);
-      Sad->AddValue(ad,1.0);
-      Taddt->Fill(ad,dt,1.0);
-      Tadp->Fill(ad,trk[i].p,1.0);
-      Tdttq->Fill(dt,qual,1.0);
-      Ttqp->Fill(qual,trk[i].p,1.0);
-      for(j=0;j<trk[i].nclust;j++){
-         lay = trk[i].clust[j].lay;
-         lad = trk[i].clust[j].lad;
-         det = trk[i].clust[j].det;
-         xh = trk[i].clust[j].xg;
-         yh = trk[i].clust[j].yg;
-         rh    = hypot(xh,yh);
-         rphih = atan2(yh,xh);
-         xp = a+b*trk[i].clust[j].zg;
-         yp = c+d*trk[i].clust[j].zg;
-         rp    = hypot(xp,yp);
-         rphip = atan2(yp,xp);
-         //
-         if(fabs(rphih-rphip)>TMath::Pi()){
-             if(rphih<rphip) rphih += 2.0*TMath::Pi();
-             else rphip = 2.0*TMath::Pi();
-         } //
-         rrphi = rh*rphih-rp*rphip;
-         //
-         x1[0] = a; x1[1] = b; x1[2] = 0.0;
-         gm->GtoL(lay,lad,det,x1,xl);
-         x1[0] = xl[0]; x1[1] = xl[1]; x1[2] = xl[2];
-         if(trk[i].clust[j].zg!=0.0){
-             x2[0] = xp; x2[1] = yp; x2[2] = trk[i].clust[j].zg;
-         }else{
-             x2[0] = a+b; x2[1] = c+d; x2[2] = 1.0;
-         } // end if trk[i].clst[j].zg!=0
-         gm->GtoL(lay,lad,det,x2,xl);
-         x2[0] = xl[0]; x2[1] = xl[1]; x2[2] = xl[2];
-         xg[0] = trk[i].clust[j].xg;
-         xg[1] = trk[i].clust[j].yg;
-         xg[2] = trk[i].clust[j].zg;
-         xl[0] = trk[i].clust[j].xl;
-         xl[1] = trk[i].clust[j].yl;
-         xl[2] = trk[i].clust[j].zl;
-         if(x1[1]!=x2[1]){
-             rx = x1[0] - x1[1]*(x1[0]-x2[0])/(x1[1]-x2[1]);
-             rz = x1[2] - x1[1]*(x1[2]-x2[2])/(x1[1]-x2[1]);
-         }else{
-             continue;
-         } // end if x1[1]!=x2[1]
-         rx = rx - xl[0];
-         rz = rz - xl[2];
-         // Fill histograms and statistics
-         TrxrzL[lay-1]->Fill(rx,rz,1.0);
-         TrxrzL[6]->Fill(rx,rz,1.0);
-         TrrprzG[lay-1]->Fill(rrphi,rz,1.0);
-         TrrprzG[6]->Fill(rrphi,rz,1.0);
-         SrxL[lay-1]->AddValue(rx,1.0);
-         SrzL[lay-1]->AddValue(rz,1.0);
-         SrxzL->AddValue(rx,1.0);
-         SrxzL->AddValue(rz,1.0);
-      } // end for j
-   } // end for i
-
-// Write out information
-
-   printf("FitAllTracks: %d tracks cut leaving %d, with a chi squared >%f\n",
-         NqualgtMAX,Sad->GetN(),trkqualMAX);
-   printf("The number of tracks with chi squared <1");
-   for(i=1;i<Nqualmax;i++) printf(",%d",i+1);
-   printf("\nFitAllTracks: %d",Nqualless[0]);
-   for(i=1;i<Nqualmax;i++) printf(",%d",Nqualless[i]);
-   printf("\n");
-   printf("FitAllTracks: RMSs of ad, rxz = %f+-%f %f+-%f\n",
-         Sad->GetRMS(),Sad->GetErrorRMS(),
-         SrxzL->GetRMS(),SrxzL->GetErrorRMS());
-
-   printf("FitAllTracks: Residuals by layer x=-rphi z ");
-   for(i=0;i<6;i++) printf("%f+-%f %f+-%f ",
-                          SrxL[i]->GetRMS(),SrxL[i]->GetErrorRMS(),
-                          SrzL[i]->GetRMS(),SrzL[i]->GetErrorRMS());
-   printf("\n");
-
-// Setup and fill projections
-
-   for(i=0;i<10;i++){
-       xm = 0.5*Double_t(i);
-       ym = xm+0.5;
-       j  = Ttqp->GetYaxis()->FindBin(xm);
-       k  = Ttqp->GetYaxis()->FindBin(ym);
-       sprintf(filename,"Track Quality for %3.1f<p<%3.1f",xm,ym);
-       Tptqp[i] = Ttqp->ProjectionX(filename,j,k,"E");
-       Tptqp[i]->SetXTitle("Track Quality (Chi squared/df)");
-   } // end for i
-
-   TH1D *Ttq  = Ttqp->ProjectionX("Ttq",0,Ttqp->GetNbinsY()+1,"E");
-   Ttq->SetXTitle("Chi^2/degree freedom");
-   TH1D *Tp   = Ttqp->ProjectionY("Tp", 0,Ttqp->GetNbinsX()+1,"E");
-   Tp->SetXTitle("Momentum GeV/c");
-   TH1D *Tdt  = Tdttq->ProjectionX("Tdt",0,Tdttq->GetNbinsY()+1,"E");
-   Tdt->SetXTitle("Distanct to true vertex (cm)");
-   TH1D *Tad  = Tadtq->ProjectionX("Tad",0,Tadtq->GetNbinsY()+1,"E");
-   Tad->SetXTitle("theta fitted line - atan(pt/p) (rad)");
-   TH1D *TrxL[7],*TrzL[7],*TrrphiG[7],*TrzG[7];
-   for(i=0;i<7;i++){
-       sprintf(hid,"TrxL%1.1d",i+1);
-       TrxL[i] = TrxrzL[i]->ProjectionX(hid,0,TrxrzL[i]->GetNbinsY()+1,"E");
-       TrxL[i]->SetXTitle("xi-x(i) (cm)");
-       sprintf(hid,"TrzL%1.1d",i+1);
-       TrzL[i] = TrxrzL[i]->ProjectionY(hid,0,TrxrzL[i]->GetNbinsX()+1,"E");
-       TrzL[i]->SetXTitle("zi-z(i) (cm)");
-//
-       sprintf(hid,"TrrphiG%1.1d",i+1);
-       TrrphiG[i] = TrrprzG[i]->ProjectionX(hid,0,
-                                             TrrprzG[i]->GetNbinsY()+1,"E");
-       TrrphiG[i]->SetXTitle("rPhii-rPhi(i) (cm)");
-       sprintf(hid,"TrzG%1.1d",i+1);
-       TrzG[i] = TrrprzG[i]->ProjectionY(hid,0,
-                                          TrrprzG[i]->GetNbinsX()+1,"E");
-       TrzG[i]->SetXTitle("zGi-zG(i) (cm)");
-   } // end for i
-
-
-   Hfile->Write();
-
-   if(printit){
-       TCanvas *c0 = new TCanvas("c0","Track quality distribution",
-                                500,100,600,700);
-       Ttqall->Draw();
-       sprintf(filename,"%s_T_tqall.ps",sfile);
-       if(printit) c0->Print(filename);
-       Ttqp->Draw("COL");
-       sprintf(filename,"%s_T_tq_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Ttq->Draw();
-       sprintf(filename,"%s_T_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tp->Draw();
-       sprintf(filename,"%s_T_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tdttq->Draw("COL");
-       sprintf(filename,"%s_T_dt_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tdt->Draw();
-       sprintf(filename,"%s_T_dt.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tadtq->Draw("COL");
-       sprintf(filename,"%s_T_ad_tq.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tad->Draw();
-       sprintf(filename,"%s_T_ad.ps",sfile);
-       if(printit) c0->Print(filename);
-       Tadp->Draw("COL");
-       sprintf(filename,"%s_T_ad_p.ps",sfile);
-       if(printit) c0->Print(filename);
-       Taddt->Draw("COL");
-       sprintf(filename,"%s_T_ad_dt.ps",sfile);
-       if(printit) c0->Print(filename);
-       for(i=0;i<7;i++){
-          TrxrzL[i]->Draw("COL");
-          sprintf(filename,"%s_T%1.1d_rx_rz.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-          TrxL[i]->Draw();
-          sprintf(filename,"%s_T%1.1d_rx.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-          TrzL[i]->Draw();
-          sprintf(filename,"%s_T%1.1d_rz.ps",sfile,i+1);
-          if(printit) c0->Print(filename);
-       } // end for i
-       for(i=0;i<10;i++){
-          Tptqp[i]->Draw();
-          sprintf(filename,"%s_T_tq_p%1.1d.ps",sfile,i);
-          if(printit) c0->Print(filename);
-       } // end for i
-   } // end if printit
-
-// Delet allocated stuff.
-   for(i=0;i<7;i++) {
-       delete TrxL[i];
-       delete TrzL[i];
-       delete TrrphiG[i];
-       delete TrzG[i];
-   } // end for i
-   for(i=0;i<10;i++) delete Tptqp[i];
-   delete Ttqp;
-   delete Tdttq;
-   delete Tadtq;
-   delete Taddt;
-   delete Tadp;
-   for(i=0;i<7;i++) delete TrxrzL[i];
-   for(i=0;i<7;i++) delete TrrprzG[i];
-   delete Ttqall;
-   delete Ttq;
-   delete Tp;
-   delete Tdt;
-   delete Tad;
-//   printf("finished with track fitting\n");
-   delete Sad;
-   delete SrxzL;
-   for(i=0;i<6;i++) delete SrxL[i];
-   for(i=0;i<6;i++) delete SrzL[i];
-   return;
-}
-//______________________________________________________________________
-void FindVertex2(ClustAl_tl &trk1,ClustAl_tl &trk2,Double_t *vt,Double_t &d){
-   // Local Variables
-   Double_t x1,y1,x2,y2;
-   Double_t z1,z2;
-   Double_t a1,b1,c1,d1;
-   Double_t a2,b2,c2,d2;
-   Double_t da,db,dc,dd,dbd;
-   Double_t den,num1,num2;
-
-   a1 = trk1.a; b1 = trk1.b; c1 = trk1.c; d1 = trk1.c;
-   a2 = trk2.a; b2 = trk2.b; c2 = trk2.c; d2 = trk2.c;
-
-   // Find z1 and z2 of points of closest approch.
-   da = a1-a2;
-   db = b1-b2;
-   dc = c1-c2;
-   dd = d1-d2;
-   dbd = b1*d2-d1*b2;
-   den = -db*db - dd*dd - dbd*dbd;
-   num1 = da*(db+d1*dbd) + dc*(dd-b1*dbd);
-   num2 = da*(db+d2*dbd) + dc*(dd-b2*dbd);
-   if(den!=0.0){
-       z1 = num1/den;
-       z2 = num2/den;
-   }else{ // parallel lines
-       z1 = 0.0;
-       z2 = 0.0;
-   } // end if den!-0.0
-
-   // find coordinate of closest approch and distance between.
-   x1    = a1+b1*z1; y1 = c1+d1*z1;
-   x2    = a2+b2*z2; y2 = c2*d2*z2;
-   vt[0] = 0.5*(x1+x2);
-   vt[1] = 0.5*(y1+y2);
-   vt[2] = 0.5*(z1+z2);
-   d     = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) );
-
-   return;
-}
-//______________________________________________________________________
-void FitVertexAll(ClustAl_tl *trk,Int_t ntrk,const char *sfile,TFile *Hfile){
-   // Local Variables
-   Int_t    i,j,k;
-   Double_t d,vz,r,q;
-   Double_t vt[3];
-   Bool_t   printit = kFALSE;
-   char     filename[80],hid[10];
-   Double_t trkqualMAX = 50.0;
-   const Int_t    IntQualMax=10;
-   Int_t    QualInt[IntQualMax];
-   AliITSstatistics *Svz   = new AliITSstatistics(4);
-   AliITSstatistics *Svx   = new AliITSstatistics(4);
-   AliITSstatistics *Svy   = new AliITSstatistics(4);
-   //   Double_t vzmi=-0.20,vzma=0.20,vxmi=-1.0,vxma=1.0,vymi=-1.0,vyma=1.0;
-
-   for(i=0;i<IntQualMax;i++) QualInt[i] = 0;
-
-   TH1F *Vvzbtq[IntQualMax];
-   for(i=0;i<IntQualMax;i++){
-       sprintf(hid,"Vvzbta%2.2d",i+1);
-       sprintf(filename,"Z of vertex of pairs for tracks with chi squared <%d",
-              i+1);
-       Vvzbtq[i] = new TH1F(hid,filename,100,-1.0,1.0);
-       Vvzbtq[i]->Sumw2();
-       Vvzbtq[i]->SetXTitle("Z of vertex (cm)");
-   } // end for i
-
-   TH2F *Vvztq  = new TH2F("Vvztq","Z of vertex of pairs vs. track quality cut",
-                          500,-0.20,0.20,200,0.0,20.0);
-   Vvztq->Sumw2();
-   Vvztq->SetXTitle("Z of vertex (cm)");
-   Vvztq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Vvzvr  = new TH2F("Vvzvr","Z vs. R of vertex of pairs",
-                                              200,-0.20,0.20,400,0.0,5.0);
-   Vvzvr->Sumw2();
-   Vvzvr->SetXTitle("Z of vertex (cm)");
-   Vvzvr->SetYTitle("R of vertex (cm)");
-   TH2F *Vdtq   = new TH2F("Vdtq","Distance between lines vs. track quality cut",
-                          500,0.0,0.2,200,0.0,50.0);
-   Vdtq->Sumw2();
-   Vdtq->SetXTitle("minimum distance between lines (cm)");
-   Vdtq->SetYTitle("Chi^2/degree freedom");
-   TH2F *Vvzd   = new TH2F("Vvzd","Z vertex vs. Dist between lines",
-                                               500,-5.0,5.0,200,0.0,5.0);
-   Vvzd->Sumw2();
-   Vvzd->SetXTitle("Z of vertex (cm)");
-   Vvzd->SetYTitle("minimmum distance between lined (cm)");
-   TH2F *Vvxvy  = new TH2F("Vvxvy","X vertex vs. Y vertex",
-                                               200,-10.0,10.0,200,-10.0,10.0);
-   Vvxvy->Sumw2();
-   Vvxvy->SetXTitle("X of vertex (cm)");
-   Vvxvy->SetYTitle("Y of vertex (cm)");
-   TH2F *Vvxvz  = new TH2F("Vvxvz","X vertex vs. Z vertex",
-                                               200,-10.0,10.0,200,-10.0,10.0);
-   Vvxvz->Sumw2();
-   Vvxvz->SetXTitle("X of vertex (cm)");
-   Vvxvz->SetYTitle("Z of vertex (cm)");
-   TH2F *Vvyvz  = new TH2F("Vvyvz","Y vertex vs. Z vertex",
-                                               200,-10.0,10.0,200,-10.0,10.0);
-   Vvyvz->Sumw2();
-   Vvyvz->SetXTitle("Y of vertex (cm)");
-   Vvyvz->SetYTitle("Z of vertex (cm)");
-
-   for(i=0;i<ntrk-1;i++){
-       if(trk[i].qual>trkqualMAX||trk[i].qual<0.0) continue;
-       for(j=i+1;j<ntrk;j++){
-          if(trk[j].qual>trkqualMAX||trk[i].qual<0.0) continue;
-          FindVertex2(trk[i],trk[j],vt,d);
-          q  = TMath::Max(trk[i].qual,trk[j].qual);
-          r  = vt[0]*vt[0] + vt[1]*vt[1];
-          r  = sqrt(r);
-          vz = vt[2];
-           Vvxvy->Fill(vt[0],vt[1],1.0);
-           Vvxvz->Fill(vt[0],vt[2],1.0);
-           Vvyvz->Fill(vt[1],vt[2],1.0);
-          Svx->AddValue(vt[0],1.0);
-          Svy->AddValue(vt[1],1.0);
-          Svz->AddValue(vt[2],1.0);
-          Vvztq->Fill(vz,q,1.0);
-          Vvzvr->Fill(vz,r,1.0);
-          Vdtq->Fill(d,q,1.0);
-          Vvzd->Fill(vz,d,1.0);
-          for(k=0;k<IntQualMax;k++)if(q<(Double_t)k+1.)QualInt[k]++;
-          for(k=0;k<IntQualMax;k++)if(q<(Double_t)k+1.)Vvzbtq[k]->Fill(vz,1.);
-       } // end for j
-   } // end for i
-
-   TH1D *Vd  = Vvzd->ProjectionY ("Vd", 0,Vvzd->GetNbinsX()+1,"E");
-   Vd->SetXTitle("minimum distance between lines (cm)");
-   TH1D *Vvr = Vvzvr->ProjectionY("Vvr",0,Vvzvr->GetNbinsX()+1,"E");
-   Vvr->SetXTitle("R of vertex (cm)");
-   TH1D *Vtq = Vdtq->ProjectionY ("Vtq",0,Vdtq->GetNbinsX()+1,"E");
-   Vtq->SetXTitle("Chi^2/degree freedom");
-//
-   TH1D *Vvx = Vvxvy->ProjectionX("Vvx",0,Vvxvy->GetNbinsY()+1,"E");
-   Vvx->SetXTitle("X of vertex (cm)");
-   TH1D *Vvy = Vvxvy->ProjectionY("Vvy",0,Vvxvy->GetNbinsX()+1,"E");
-   Vvy->SetXTitle("Y of vertex (cm)");
-   TH1D *Vvz = Vvxvz->ProjectionX ("Vvz",0,Vvxvz->GetNbinsY()+1,"E");
-   Vvz->SetXTitle("Z of vertex (cm)");
-
-   printf("FitVertexAll: N(track pairs Qual <=1");
-   for(k=1;k<IntQualMax;k++) printf(",%d",k+1);
-   printf("\nFitVertexAll: %d",QualInt[0]);
-   for(k=1;k<IntQualMax;k++) printf(",%d",QualInt[k]);
-   printf("FitVertexAll: RMS of vx vy vz=%e+-%e %e+-%e %e+-%e\n",
-         Svx->GetRMS(),Svx->GetErrorRMS(),
-         Svy->GetRMS(),Svy->GetErrorRMS(),
-         Svz->GetRMS(),Svy->GetErrorRMS());
-   Hfile->Write();
-
-   if(printit){
-       TCanvas *c1 = new TCanvas("c1","Vertex info",500,100,600,700);
-       for(k=0;k<IntQualMax;k++){
-          Vvzbtq[k]->Draw();
-          sprintf(filename,"%s_V_vz_tq%2.2d.ps",sfile,k+1);
-          if(printit) c1->Print(filename);
-       } // end for k
-       Vvzd->Draw("COL");
-       sprintf(filename,"%s_V_vz_d.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvzvr->Draw("COL");
-       sprintf(filename,"%s_V_vz_vr.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vdtq->Draw("COL");
-       sprintf(filename,"%s_V_d_tq.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvztq->Draw("COL");
-       sprintf(filename,"%s_V_vz_tq.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvz->Draw();
-       sprintf(filename,"%s_V_vz.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vd->Draw();
-       sprintf(filename,"%s_V_d.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvr->Draw();
-       sprintf(filename,"%s_V_vr.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vtq->Draw();
-       sprintf(filename,"%s_V_tq.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvxvy->Draw("COL");
-       sprintf(filename,"%s_V_vx_vy.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvxvz->Draw("COL");
-       sprintf(filename,"%s_V_vx_vz.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvyvz->Draw("COL");
-       sprintf(filename,"%s_V_vy_vz.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvx->Draw();
-       sprintf(filename,"%s_V_vx.ps",sfile);
-       if(printit) c1->Print(filename);
-       Vvy->Draw();
-       sprintf(filename,"%s_V_vy.ps",sfile);
-       if(printit) c1->Print(filename);
-   } // end if printit
-   for(k=0;k<IntQualMax;k++) delete Vvzbtq[k];
-   delete Vvzd;
-   delete Vvzvr;
-   delete Vdtq;
-   delete Vvztq;
-   delete Vvz;
-   delete Vd;
-   delete Vvr;
-   delete Vtq;
-   delete Vvxvy;
-   delete Vvxvz;
-   delete Vvyvz;
-   delete Vvx;
-   delete Vvy;
-   delete Svz;
-   delete Svx;
-   delete Svy;
-
-   return;
-}
-//______________________________________________________________________
-void OnlyOneGeometry(char *filename,AliITSgeom *gm,AliITSgeom &gm2,
-                    Float_t trans[],Float_t rot[]){
-
-    ifstream f_in(filename);
-    if(f_in){
-       printf("reading from geometry file %s\n",filename);
-       gm2.ReadGeom(f_in);
-       f_in.close();
-    }else{
-       f_in.close();
-       gm2 = *gm;
-       gm2.RandomCylindericalChange(trans,rot);
-       printf("writting to geometry file %s\n",filename);
-       ofstream f_out(filename);
-       gm2.PrintGeom(f_out);
-       f_out.close();
-    } // end if
-    return;
-}
-//______________________________________________________________________
-void deleteClustAl(ClustAl_tl *trk,Int_t ntrk){
-
-    Int_t i;
-
-    for(i=0;i<ntrk;i++) delete[] trk[i].clust;
-    return;
-}
-//______________________________________________________________________
-void RunAlignment(Int_t evnt,Float_t fraction=1.0){
-// define some variables for later use.
-//
-    printf("gAlice=%p and evnt=%d\n",gAlice,evnt);
-//
-    Int_t nparticles = gAlice->GetEvent(evnt);
-    printf("nparticles %d\n",nparticles);
-    if (nparticles <= 0) return; /* get next event */
-
-// Get pointers to Alice detectors and Clusts containers
-    AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
-    if(!ITS) return;          /* error no ITS data exit */
-    TTree *TH   = gAlice->TreeH();
-    Int_t Ntrkp =  (Int_t) TH->GetEntries(),ntrk;
-    AliITSgeom gm2,gm3,*gm = ITS->GetITSgeom();
-    Int_t Nmods = gm->GetIndexMax();
-    Float_t   trans[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
-                         2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
-                         2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // cm
-    Float_t   rots[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
-                        2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
-                        2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // rad
-   Float_t     tran[3] = {0.0,0.0,0.0},rot[3] = {0.0,0.0,0.0};
-   TFile *Hfile;
-   Float_t     Rdta[6];
-   Int_t Itimes=0,i,j,badmod;
-   Double_t Chi2b;
-
-// Array (stucture) of clusts for the first and second layer
-// this should be replaced with either clusters or digits
-// when they are proporly defined.
-    AliITSAlignmentTrack *trk = new AliITSAlignmentTrack[Ntrkp];
-    TObjArray *mods = new TObjArray(Nmods);
-    AliITSAlignmentModule *mod;
-
-    printf("Ntrkp=%d\n",Ntrkp);
-
-    FillAliITSAlignmentTrack(trk,ntrk,Ntrkp,TH,ITS,fraction);
-
-    for(Int_t Isigmas=0;Isigmas<1;Isigmas++){
-//
-//      tran[0] = sigma1;
-//     tran[1] = sigma2;
-//     tran[2] = sigma3;
-       if(Itimes==0){ tran[0] = trans[Isigmas];
-       }else tran[0] = 0.0;
-       if(Itimes==1){ tran[1] = trans[Isigmas];
-       }else tran[1] = 0.0;
-       if(Itimes==2){ tran[2] = trans[Isigmas];
-       }else tran[2] = 0.0;
-       if(Itimes==3){ rot[0] = rots[Isigmas];
-       }else rot[0] = 0.0;
-       if(Itimes==4){ rot[1] = rots[Isigmas];
-       }else rot[1] = 0.0;
-       if(Itimes==5){ rot[2] = rots[Isigmas];
-       }else rot[2] = 0.0;
-       printf("tran= %e %e %e (cm), rot=%e %e %e (rad)\n",
-              tran[0],tran[1],tran[2],rot[0],rot[1],rot[2]);
-//
-       gm2 = *gm;
-       gm2.RandomCylindericalChange(tran,rot);
-       gm3 = gm2;
-       Hfile = new TFile("Alignment_geom_0_2.root","RECREATE",
-                         "comparison of geometry before refitting");
-        PlotGeomChanges(gm,&gm2,Hfile,Rdta);
-       Hfile -> Close();
-
-       for(i=0;i<Nmods;i++) mods->AddAt(new AliITSAlignmentModule(i,&gm3,
-                                                                  ntrk,trk),i);
-//
-       j = 0;
-       do{
-           for(i=0;i<ntrk;i++) trk[i].SetGlobalPosition(&gm2);
-           //
-           for(i=0;i<ntrk;i++){
-               trk[i].FitToFunction(2,&gm3);
-           } // end for i
-           // find detector with the worst Chi2.
-           j++;
-           Chi2b  = 0.0;
-           badmod = 0;
-           for(i=0;i<Nmods;i++){
-               if(((AliITSAlignmentModule *)(mods->At(i)))->GetChi2()>Chi2b){
-                   Chi2b =((AliITSAlignmentModule *)(mods->At(i)))->GetChi2();
-                   badmod = i;
-               } // end if
-           } // end for i
-           mod = (AliITSAlignmentModule *)(mods->At(badmod));
-           mod->AlignModule();
-           gm3.SetTrans(badmod,mod->GetTranslationVector());
-           gm3.SetByAngles(badmod,mod->GetRotationAngles());
-       }while(j<Nmods);
-    } // end for Isigmas
-//
-    Hfile = new TFile("Alignment_geom_0_3.root","RECREATE",
-                     "comparison of geometry after refitting");
-    PlotGeomChanges(gm,&gm3,Hfile,Rdta);
-    Hfile->Close();
-    printf("Event %d done\n",evnt);
-//
-    delete[] trk;            // now delet memory allocated above.
-}
-//______________________________________________________________________
diff --git a/ITS/AnalysisITSAlignment.h b/ITS/AnalysisITSAlignment.h
deleted file mode 100644 (file)
index aee7640..0000000
+++ /dev/null
@@ -1,48 +0,0 @@
-#ifndef ALIGNITS_H
-#define ALIGNITS_H
-//////////////////////////////////////////////////////////////////////////
-//  Alice ITS first detector alignment program.                         //
-//                                                                      //
-// version: 0.0.0 Draft.                                                //
-// Date: April 18 1999                                                  //
-// By: Bjorn S. Nilsen                                                  //
-//                                                                      //
-//////////////////////////////////////////////////////////////////////////
-#include "AliITS.h"
-#include "TTree.h"
-#include "TClonesArray.h"
-#include "AliRun.h"
-#include "AliITSAlignmentTrack.h"
-
-// Data structure to hold averaged clusts.
-struct ClustAl_sl{
-    Int_t lay,lad,det;
-    Float_t xg,yg,zg,xl,yl,zl;
-};
-struct ClustAl_tl{
-    Int_t    track,nclust;  // track number and number of data points.
-    ClustAl_sl *clust;        // data points to fit.
-    Float_t  a,b,c,d,a0,b0,c0,d0,qual;  // fit parameters and fit quality.
-    Float_t  px,py,pz,p,pt;
-    // x=a+b*z and y=c+d*z;
-    // x=a0+b0*y and z=c0+d0*y in coordinate system of clust[0].lay,lad,det
-};
-
-void HitsTo(ClustAl_tl *trk,Int_t &ntrk,Int_t,TTree *TH,AliITS *ITS,
-           Float_t nsigmaT1,Float_t nsigmaT2,Float_t nsigmaT3,
-           Float_t nsigmaR1,Float_t nsigmaR2,Float_t nsigmaR3);
-void HitsToClustAl(ClustAl_tl *trk,Int_t &ntrk,Int_t nt,TTree *TH,
-                  AliITS *ITS,Float_t fraction);
-void PlotGeomChanges(AliITSgeom *gt,AliITSgeom *gc,TFile *Hfile,Float_t *Rdta);
-void FillGlobalPositions(ClustAl_tl *trk,Int_t ntrk,AliITSgeom *g);
-void FitAllTracks(ClustAl_tl *trk,Int_t ntrk,Float_t *v0,AliITSgeom *gm,
-                 const char *sfile,TFile *Hfile,Float_t *Fdta,Int_t *Ndta);
-void FitVertexAll(ClustAl_tl *trk,Int_t ntrk,const char *sfile,TFile *Hfile);
-void OnlyOneGeometry(char *filename,AliITSgeom *gm,AliITSgeom &gm2,
-                    Float_t trans[],Float_t rot[]);
-void deleteClustAl(ClustAl_tl *trk,Int_t ntrk);
-void FillAliITSAlignmentTrack(AliITSAlignmentTrack *trk,Int_t &ntrk,Int_t nt,
-                             TTree *TH,AliITS *ITS,Float_t fraction);
-void RunAlignment(Int_t evnt,Float_t fraction);
-
-#endif
index f4e2c16..a34b2fc 100644 (file)
@@ -28,9 +28,8 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                AliITSHuffman.cxx AliITSClusterFinderSSD.cxx \
                AliITSclusterSSD.cxx AliITSpackageSSD.cxx \
                AliITSdictSSD.cxx AliITSgeomSPD300.cxx AliITSgeomSPD425.cxx \
-               vertex.cxx AliITSstatistics.cxx AliITSstatistics2.cxx \
-               AnalysisITSAlignment.cxx  AliITSAlignmentTrack.cxx \
-               AliITSAlignmentModule.cxx 
+               AliITSstatistics.cxx AliITSstatistics2.cxx \
+               AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx 
 
 # Fortran sources
 
diff --git a/ITS/vertex.cxx b/ITS/vertex.cxx
deleted file mode 100755 (executable)
index dcbd3c6..0000000
+++ /dev/null
@@ -1,718 +0,0 @@
-#include <stdio.h>
-#include <math.h>
-#include <time.h>
-#include "TMinuit.h"
-#include "TRandom.h"
-#include "TTree.h"
-#include "TClonesArray.h"
-#include "TProfile.h"
-#include "TH2.h"
-#include "TArray.h"
-#include "TCanvas.h"
-#include "TString.h"
-#include "AliRun.h"
-#include "AliITS.h"
-#include "AliITSgeom.h"
-#include "TParticle.h"
-#include "vertex.h"
-
-// Global in this file only. Needed for Minuit
-Int_t gstrt,gend;  // starting and ending bin values
-TH1D  *gH1;        // histogram to be fit.
-
-void FindVertexs(Int_t evnt,Float_t frac=1.0,Float_t len=16.82){
-
-      // Get pointers to Alice detectors and Hits containers
-      AliDetector  *ITS       = gAlice->GetDetector("ITS");
-      if(!ITS) return;          // error no ITS data exit
-//      TClonesArray *Particles = gAlice->Particles();
-      TClonesArray *ITShits   = ITS->Hits();
-      TTree        *TH        = gAlice->TreeH();
-      Int_t        ntracks    =  (Int_t)TH->GetEntries();
-//      ntracks = 10;
-
-      printf(" %d primary tracks in the file\n",ntracks);
-      // Array (stucture) of hits for the first and second layer
-      // this should be replaced with either clusters or digits
-      // when they are proporly defined.
-      Hit_tl *spd1 = new Hit_tl[ntracks];
-      Hit_tl *spd2 = new Hit_tl[ntracks];
-      Hit_tl **spdi = new Hit_tl*[ntracks];
-      Hit_tl **spdo = new Hit_tl*[ntracks];
-      for(Int_t i=0;i<ntracks;i++){
-         spdi[i] = &(spd1[i]);
-//       printf("spdi[%d]=%p spd1[%d]=%p  ",i,spdi[i],i,&(spd1[i]));
-         spdo[i] = &(spd2[i]);
-//       printf("spdo[%d]=%p spd2[%d]=%p\n",i,spdo[i],i,&(spd2[i]));
-      } // end for i
-      Int_t   i1max,i2max;
-      Float_t vz;
-
-      HitsToV(spdi,i1max,spdo,i2max,ntracks,TH,ITShits,frac,len);
-      printf("back in Macro i1max=%d i2max=%d\n",i1max,i2max);
-
-//      Float_t r1=0.0,r2=0.0;
-//      for(i=0;i<TMath::Max(i1max,i2max);i++){
-//       if(i<i1max) r1 += spdi[i]->r;
-//       if(i<i2max) r2 += spdo[i]->r;
-//      } // end for i
-//      printf("<r1>=%f i1max=%d\t<r2>=%f i2max=%d\n",
-//          r1/Float_t(i1max),i1max,r2/Float_t(i2max),i2max);
-
-      vz = vertexSlow(spdi,i1max,spdo,i2max);
-
-      printf("Slow Sorted event=%d Zvertex=%f\n",evnt,vz);
-
-//      vz = vertex(spdi,i1max,spdo,i2max);
-
-//      printf("Phi sorted event=%d Zvertex=%f\n",evnt,vz);
-
-//      vz = vertexEta(spdi,i1max,spdo,i2max);
-
-//      printf("Eta sorted event=%d Zvertex=%f\n",evnt,vz);
-
-      return;
-}
-
-
-Bool_t L_SortPhi(const Hit_tl *s1,const Hit_tl *s2){
-  // Phi sorting function for qsort.
-   Float_t a;
-
-   a = s1->phir - s2->phir;
-   if(a<0.0) return kTRUE;
-   if(a>0.0) return kFALSE;
-   a = s1->etar - s2->etar;
-   if(a<0.0) return kTRUE;
-   if(a>0.0) return kFALSE;
-   return kFALSE;
-}
-
-Bool_t L_SortEta(const Hit_tl *s1,const Hit_tl *s2){
-  // Eta sorting function for qsort.
-   Float_t a;
-
-   a = s1->etar - s2->etar;
-   if(a<0.0) return kTRUE;
-   if(a>0.0) return kFALSE;
-   a = s1->phir - s2->phir;
-   if(a<0.0) return kTRUE;
-   if(a>0.0) return kFALSE;
-   return kFALSE;
-}
-
-void hpsortPhi(Hit_tl **ra,Int_t n){
-   Int_t i,ir,j,l;
-   Hit_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 && L_SortPhi(ra[j],ra[j+1]) ) j++;
-        if( L_SortPhi(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 hpsortEta(Hit_tl **ra,Int_t n){
-   Int_t i,ir,j,l;
-   Hit_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 && L_SortEta(ra[j],ra[j+1]) ) j++;
-        if( L_SortEta(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 fillStructure(Hit_tl **spd,Int_t &nspd,
-                  Float_t *xv,Int_t *id,Int_t track,Int_t nhits){
-    Float_t PI2 = 2.0*TMath::Pi();
-//    Float_t PI  = TMath::Pi();
-    Int_t i;
-    Float_t x,y,z,zr,r,phi,eta,rphi;
-    Float_t dz    = 0.0300; // 300 microns
-    Float_t drphi = 0.0050; //  50 microns
-
-    i = nspd;
-//   if(nspd<5) printf("fill: spd=%p spd[%d]=%p i=%d spd[i]=%p id=%d %d %d\n",
-//                   spd,nspd,spd[nspd],i,spd[i],id[0],id[1],id[2]);
-    x = xv[0];
-    y = xv[1];
-    z = xv[2];
-    r = sqrt(x*x+y*y);
-    phi = atan2(y,x);
-    if(phi<0.0) phi += PI2;
-    rphi = r*phi;
-    eta  = -0.5*tan(0.5*atan2(r,z));
-
-    spd[i]->track = track;
-    spd[i]->n     = nhits;
-    spd[i]->lad   = id[1];
-    spd[i]->det   = id[2];
-    spd[i]->x     = x;
-    spd[i]->y     = y;
-    spd[i]->z     = z;
-    spd[i]->r     = r;
-    spd[i]->phi   = phi;
-    spd[i]->eta   = eta;
-
-    zr   = dz    * Float_t(Int_t(z   /   dz)) + 0.5*dz;
-    rphi = drphi * Float_t(Int_t(rphi/drphi)) + 0.5*drphi;
-
-    spd[i]->xr   =  cos(phi)*rphi/r;
-    spd[i]->yr   = -sin(phi)*rphi/r;
-    spd[i]->zr   = zr;
-    spd[i]->phir = rphi/r;
-    spd[i]->etar = -0.5*tan(0.5*atan2(r,zr));
-
-    nspd++;
-    return;
-}
-
-void HitsToV(Hit_tl **spdi,Int_t &nspdi,Hit_tl **spdo,Int_t &nspdo,
-           Int_t ntracks,TTree *TH,TClonesArray *ITShits,
-            Float_t fraction=1.0,Float_t len=16.82){
-    Int_t     i,t,h,n,nb,nhits,id[4],idold[4],Iseed,ieta=0,itrk=0,ieta2=0;
-    Float_t   x[3],xb[3],xl[3];
-    AliITShit *itsHit;
-    TClonesArray *Part = gAlice->Particles();
-    TParticle *part;
-
-    nspdi = nspdo = 0;
-
-    idold[0] = idold[1] = idold[2] = idold[3] = 0;
-    xb[0]    = xb[1]    = xb[2]    = 0.0;
-    n = 0;
-    Iseed = time(0);
-    printf("HitsToV: Iseed=%d fraction=%f Pixel length=%fcm\n",
-          Iseed,fraction,len);
-    gRandom->SetSeed(Iseed);
-//    printf("HitsToV: gRandom->Rndm(1)=%f\n",gRandom->Rndm(1));
-    for(t=0;t<ntracks;t++){
-       if(fraction<gRandom->Rndm(1)) continue; // skip some tracks
-       itrk++;
-       gAlice->ResetHits();
-       nb       = TH->GetEvent(t);
-       nhits    = ITShits->GetEntriesFast();
-       for(h=0;h<nhits;h++){
-           itsHit = (AliITShit *) ITShits->UncheckedAt(h);
-           itsHit->GetDetectorID(id[0],id[1],id[2]); id[3]=t;
-           itsHit->GetPositionG(x[0],x[1],x[2]);
-           itsHit->GetPositionL(xl[0],xl[1],xl[2]);
-           if(h==0){
-               part = (TParticle *) (Part->UncheckedAt(itsHit->GetTrack()));
-               if(TMath::Abs(part->Eta())<=1.0) ieta++;
-               if(TMath::Abs(part->Eta())<=0.5) ieta2++;
-           } // end if h==0
-           if(TMath::Abs(x[2]/len) >= 1.0) continue;
-           if(x[0]==0.0&&x[1]==0.0){
-               printf("Hitsto: t=%d h=%d/%d id=%d,%d,%d x=%f,%f,%f\n",
-                      t,h,nhits,id[0],id[1],id[2],x[0],x[1],x[2]); 
-               continue;
-           } // end if x[0]==x[1]==0.0
-           if(!(id[0]==idold[0]&&id[1]==idold[1]&&
-                id[2]==idold[2]&&id[3]==idold[3])){
-               if(!(n<=0 || (xb[0]==0.0&&xb[1]==0.0))){
-                   for(i=0;i<3;i++) xb[i]   /= Float_t(n);
-                   if(idold[0]==1){
-                       fillStructure(spdi,nspdi,xb,idold,t,n);
-                   }
-                   if(idold[0]==2){
-                       fillStructure(spdo,nspdo,xb,idold,t,n);
-                   }
-                   if(nspdi>ntracks || nspdo>ntracks){
-                       printf("Hitsto: fill error,"
-                              " nspdi=%d nspdo=%d ntracks=%d\n",
-                              nspdi,nspdo,ntracks);
-                   } // end if fill error
-               } // end if n>0
-               for(i=0;i<4;i++) idold[i] = id[i];
-               for(i=0;i<3;i++) xb[i]    = 0.0;
-               n = 0;
-           } // end if id != idold
-           for(i=0;i<3;i++) xb[i] += x[i];
-           n++;
-       } // end for h
-    } // end for t
-    printf("exiting HitsToV: %d primary tracks in eta=+-1\n",ieta);
-    printf("exiting HitsToV: %d primary tracks #2 in eta=+-0.5\n",ieta2);
-    printf("exiting HitsToV: %d primary tracks in file used\n",itrk);
-    return;
-}
-
-Float_t vertex(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
-   Float_t r1    = 3.910078; // radius at which hit is from beam axis for spd1
-   Float_t r2    = 6.955933; // radius at which hit is from beam axis for spd2
-   Float_t DPhi  = 0.005;    // maximum allowed difference in phi angle
-   Float_t DZ    = 12.5;      // maximum allowed difference in z position
-   Float_t avt   = 0.0;
-   Int_t   nbinx = 2000;
-   Int_t   start = 0;
-   Bool_t  mod   = kFALSE;
-   Int_t   i,j;
-   Float_t zv,av,su,i0,i1,i2,x0,x1,dphi,dz;
-   Float_t PI2 = 2.0*TMath::Pi();
-   Float_t PI  = TMath::Pi();
-
-   // sort according to phi.
-   hpsortPhi(spdi,i1max);
-   hpsortPhi(spdo,i2max);
-
-   // find best vertex allong z.
-   TH2S *Pvzphi   = new TH2S("Pvzphi","Phi: Posible Z vertecies vs. Phi",
-                              nbinx,-1.0,1.0,32,0.0,PI2);
-   Pvzphi->Sumw2(); // collect good statitics
-   TH1F *Pvzfl   = new TH1F("Pvzfl","Phi: Posible Z vertecies flattened",
-                        nbinx,-1.0,1.0);
-   Pvzfl->Sumw2();
-   TH1F *Pvztr  = new TH1F("Pvztr","Phi: Z Vertex found for True Tracks",
-                        nbinx,-1.0,1.0);
-   Pvztr->Sumw2(); // collect good statitics
-
-   for(i=0;i<i1max;i++){
-       if(spdi[i]->r==0.0) {printf("spdi[%d]->r=0.0\n",i);continue;}
-       for(;spdo[start]->phir<spdi[i]->phir-1.5*DPhi;start++);
-       for(j=start;(spdo[j]->phir < spdi[i]->phir+DPhi) && (j<i2max);j++){
-          dphi = fabs(spdo[j]->phir - spdi[i]->phir); if(dphi>PI) dphi -= PI;
-          dz   = fabs(spdo[j]->zr   - spdi[i]->zr);
-          if(dphi>DPhi) continue;
-           if(dz>DZ)     continue; // If outside dz range skip it
-          r1   = spdi[i]->r;
-          r2   = spdo[j]->r;
-          if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
-          else{
-              printf("vertex: spdi[%d]->r=%f = spdo[%d]->r=%f\n",i,r1,j,r2);
-              continue;
-          } // end if else r2-r1!=0.0
-          Pvzphi->Fill(zv,spdi[i]->phir);
-          Pvzfl->Fill(zv);
-          if(spdi[i]->track == spdo[j]->track) Pvztr->Fill(zv);
-       } // end for j
-   } // end for i
-
-   TH1D *Pvzdef  = Pvzphi->ProjectionX("Phi: Posible Z vertecies",1,nbinx,"E");
-
-   i0 = Pvzfl->GetBinContent(1);
-   x0 = Float_t(1);
-   i1 = Pvzfl->GetBinContent(nbinx);
-   x1 = Float_t(nbinx);
-   su = x0-x1; if(su==0.0) return -200.0;
-   su = (i0-i1)/su;
-   x0 = x0-su*i0;
-   x1 = su;
-   for(i=1;i<=nbinx;i++){
-      su = x1*Float_t(i) + x0;
-      Pvzfl->AddBinContent(i,-su);
-   } // end for i
-
-   printf("Phi:            mean=%f RMS=%f w2=%f\n",
-               Pvzdef->GetMean(),Pvzdef->GetRMS(),Pvzdef->GetSumOfWeights());
-   printf("Phi: Flattened mean=%f RMS=%f w2=%f\n",
-                 Pvzfl->GetMean(),Pvzfl->GetRMS(),Pvzfl->GetSumOfWeights());
-
-   av = 3.0 * Pvzfl->Integral(1,nbinx)/Float_t(nbinx);
-   printf("Phi: Flattened av=%f Pvzfl->Max=%f nbinx=%d\n",
-         av,Pvzfl->GetMaximum(),nbinx);
-
-   su  = 0.0;
-   avt = 0.0;
-   i1  = Pvzfl->GetBinContent(1);
-   i2  = Pvzfl->GetBinContent(2);
-   for(i=2;i<nbinx;i++){
-       i0 = i1; // old i-1 value
-       i1 = i2; // old i value
-       i2 = Pvzfl->GetBinContent(i+1); // new i+1 value
-       if(i1 > av && (i0 > av || i2 > av ) ){
-          if(!mod) mod = kTRUE;
-       } else if(mod) break;
-       if(mod){ // inside peak
-          su  += i1;
-          avt += i1 * Pvzfl->GetXaxis()->GetBinCenter(i);
-       } // end if mod
-   } // end for i
-
-   if(su!=0.0) zv = avt/su;
-   else        zv = -100.0; // an unphysical value
-
-   TCanvas *c0 = new TCanvas("c0","Alice ITS vertex finder", 400,10,600,700);
-   Pvzphi->Draw("col");
-   c0->Print("vertex5_P_vz_phi.ps");
-   Pvzdef->Draw();
-   c0->Print("vertex5_P_vz_def.ps");
-   Pvzfl->Draw();
-   c0->Print("vertex5_P_vz_flat.ps");
-   Pvztr->Draw();
-   c0->Print("vertex5_P_vz_true.ps");
-
-   return zv;
-}
-
-void Chi2Gauss(Int_t &npar,Double_t *gin,Double_t &f,
-              Double_t *par,Int_t iflag){
-    Double_t chi2 = 0.0;
-    Double_t delta,h,x,eh;
-
-    for(Int_t i=gstrt;i<gend;i++){
-       h = gH1->GetBinContent(i);
-       eh = TMath::Sqrt(h);
-       if(eh <= 0.0) eh = 1.0;
-       x = gH1->GetXaxis()->GetBinCenter(i);
-       delta = (h - par[0]*TMath::Gaus(x,par[1],par[2]) - par[3])/eh;
-       chi2 += delta*delta;
-    } // end for i
-    f = chi2;
-    return;
-}
-
-Float_t vertexSlow(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
-   Float_t r1    = 3.910078; // radius at which hit is from beam axis for spd1
-   Float_t r2    = 6.955933; // radius at which hit is from beam axis for spd2
-   Float_t DPhi  = 0.005;    // maximum allowed difference in phi angle
-   Float_t BDphi;
-   Float_t DZ    = 12.5;      // maximum allowed difference in z position
-   Float_t avt   = 0.0;
-   Int_t   nbinx = 2000;
-   Int_t   nbst;
-   Bool_t  mod   = kFALSE;
-   Int_t   i,j;
-   Float_t zv,av,su,i0,i1,i2,x0,x1,dphi,dz;
-   Float_t PI2 = 2.0*TMath::Pi();
-   Float_t PI  = TMath::Pi();
-
-   if(i1max<=0||i2max<=0) return -1.0E10;
-
-   // find best vertex allong z.
-   TH2S *Svzphi  = new TH2S("Svzphi","Slow: Posible Z vertecies vs. Phi",
-                           nbinx,+0.0,10.0,32,0.0,PI2);
-   Svzphi->Sumw2(); // collect good statitics
-   TH2S *Svzdphi  = new TH2S("Svzdpii","Slow: Posible Z vertecies vs. DPhi",
-                            200,+4.0,6.0,20,0.0,10*DPhi);
-   Svzdphi->Sumw2(); // collect good statitics
-   TH1F *Svzfl  = new TH1F("Svzfl","Slow: Posible Z vertecies Flattened",
-                           nbinx,+4.0,6.0);
-   Svzfl->Sumw2();
-   TH1F *Svztr = new TH1F("Svztr","Slow: Z Vertex found by True Tracks",
-                         nbinx,+4.0,6.0);
-   Svztr->Sumw2(); // collect good statitics
-
-   printf("Svertex: i1max=%d i2max=%d\n",i1max,i2max);
-
-   for(i=0;i<i1max;i++) for(j=0;j<i2max;j++) {
-       dphi = fabs(spdo[j]->phir - spdi[i]->phir); if(dphi>PI) dphi -= PI;
-       dz   = fabs(spdo[j]->zr   - spdi[i]->zr);
-       r1   = spdi[i]->r;
-       r2   = spdo[j]->r;
-       if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
-       else{
-          printf("vertex_slow: spdi[%d]->r=%f = spdo[%d]->r=%f\n",i,r1,j,r2);
-          continue;
-       } // end if else r1-r2!=0.0
-//       if(j<10&&i<10) printf("zv=%e dphi=%e,r1=%e,r2=%e,dz=%e\n",
-//                          zv,dphi,r1,r2,dz);
-       Svzdphi->Fill(zv,dphi);
-       if(dphi>DPhi) continue; // If outside DPhi (momentum) range, skip it.
-       if(dz>DZ)     continue; // If outside dz range, skip it.
-       Svzphi ->Fill(zv,spdi[i]->phir);
-       Svzfl->Fill(zv);
-       if(spdi[i]->track == spdo[j]->track) Svztr->Fill(zv);
-   } // end for i,j
-
-   TH1D *Svzdef = Svzphi->ProjectionX("Slow: Posible Z vertecies",0,nbinx,"E");
-
-   i0 = Svzfl->GetBinContent(1);
-   x0 = Float_t(1);
-   i1 = Svzfl->GetBinContent(nbinx);
-   x1 = Float_t(nbinx);
-   su = x0-x1; if(su==0.0) return -200.0;
-   su = (i0-i1)/su;
-   x0 = x0-su*i0;
-   x1 = su;
-   for(i=1;i<=nbinx;i++){
-      su = x1*Float_t(i) + x0;
-      Svzfl->AddBinContent(i,-su);
-   } // end for i
-
-   printf("Slow:           mean=%f RMS=%f w2=%f\n",
-         Svzdef->GetMean(),Svzdef->GetRMS(),Svzdef->GetSumOfWeights());
-   printf("Slow: Flattened mean=%f RMS=%f w2=%f\n",
-         Svzfl->GetMean(),Svzfl->GetRMS(),Svzfl->GetSumOfWeights());
-
-   av = 3.0 * Svzfl->Integral(1,nbinx)/Float_t(nbinx);
-   printf("Slow: Flattened av=%f Tvxps->Max=%f nbinx=%d\n",
-         av,Svzfl->GetMaximum(),nbinx);
-
-   {// find best dPhi that masimizes the signal/noise ration.
-       Float_t sn=0.0,sig=0.0,nois=0.0,ns;
-       Int_t   iend = Svzdphi->GetYaxis()->GetNbins();
-       Int_t   jend = Svzdphi->GetXaxis()->GetNbins();
-       Int_t   ipeak = Svzdphi->GetXaxis()->FindBin(Svzdef->GetMean());
-       nbst = 0;
-       for(i=0;i<=iend;i++){
-          sig += Svzdphi->GetCellContent(ipeak,i);
-          ns   = 0.0;
-          for(j=1;j<6;j++) ns += Svzdphi->GetCellContent(j,i);
-          for(j=jend-6;j<jend;j++) ns += Svzdphi->GetCellContent(j,i);
-          nois += 0.1*ns;
-          if(nois<=0.0) continue;
-          if((sig-nois)/nois>sn){
-              sn = (sig-nois)/nois;
-              nbst = i;
-          } // end if
-//        printf("Svertex: bin=%d signal/noise=%e\n",i,(sig-nois)/nois);
-       } // end for i
-   } // end find best dPhi
-
-   if(nbst<=0) nbst = 1; // must have some data.
-   BDphi = Svzdphi->GetYaxis()->GetBinUpEdge(nbst);
-   TH1D *Svzbst = Svzdphi->ProjectionX("Slow: Best, Z vertecies",0,nbst,"E");
-
-   { // Start Muinuit fitting
-       // initilize Minuit package
-       gH1   = Svzbst; // histogram to be fit
-       gstrt = Svzbst->GetXaxis()->FindBin(
-        Svzdef->GetMean()-2.0*(Svzdef->GetRMS()));//histogram start bin value
-       gend  = Svzbst->GetXaxis()->FindBin(
-        Svzdef->GetMean()+2.0*(Svzdef->GetRMS()));//histogram end   bin value
-       TMinuit *gMinuit = new TMinuit(4); // init Minuit
-       gMinuit->SetFCN(Chi2Gauss);  // chi^2 function with Gaussian
-       Double_t arglist[10];  // Munuit parameter array
-       Int_t ierflg = 0; // Muniut error flag
-       //arglist[0] = 1.;
-       //gMinuit->mnexcm("SET ERR",arglist,1,ierglg);
-       // Set starting values and step size for parameters
-       Double_t vstart[4],step[4];
-       { // find background
-          Float_t ns   = 0.0;
-          Int_t jend = Svzbst->GetXaxis()->GetNbins();
-          for(j=1;j<6;j++) ns += Svzbst->GetBinContent(j);
-          for(j=jend-6;j<jend;j++) ns += Svzbst->GetBinContent(j);
-          vstart[3] = 0.1*ns;
-       } // end find backgrount
-       vstart[1] = Svzbst->GetMean();
-       vstart[2] = Svzbst->GetRMS();
-       if(vstart[2] == 0.0) vstart[2] = 0.04;
-       vstart[0] = (Svzbst->GetEntries() - (Svzbst->GetNbinsX())*vstart[3])*
-          vstart[2]/TMath::Sqrt(TMath::Pi());
-       if(vstart[0]<=0.0) vstart[0] = 1.0;
-       for(i=0;i<4;i++) step[i] = 0.05*vstart[i];
-       if(vstart[3] <= 0.0) step[3] = 0.1;
-       step[1] = 0.01; // mean expected at about zero set step my hand
-       gMinuit->mnparm(0,"Const",vstart[0],step[0],0,0,ierflg);
-       gMinuit->mnparm(1,"Mean" ,vstart[1],step[1],0,0,ierflg);
-       gMinuit->mnparm(2,"Sigma",vstart[2],step[2],0,0,ierflg);
-       gMinuit->mnparm(3,"Offst",vstart[3],step[3],0,0,ierflg);
-       // Now ready for minimization step
-       //arglist[0] = 500.; // Maximum number of calls
-       //arglist[1] = 1.;   // Tolorance
-       gMinuit->mnexcm("MIGRAD",arglist,0,ierflg); // do minimization
-       gMinuit->mnexcm("MINO",arglist,0,ierflg);   // find best errors
-       // Get results
-       Double_t parmin,edm,errdef;
-       Int_t nvpar,nparx,icstat;
-       gMinuit->mnstat(parmin,edm,errdef,nvpar,nparx,icstat);
-       printf("Svertex: chi2gauss=%e edist=%e istat=%d dPhi=%e\n",
-             parmin,edm,icstat,BDphi);
-       TString chnam = TString(10);
-       Double_t par[4],epar[4],empar[4],eppar[4];
-       for(i=0;i<4;i++){ 
-          gMinuit->mnpout(i,chnam,par[i],epar[i],empar[i],eppar[i],ierflg);
-          gMinuit->mnerrs(i,eppar[i],empar[i],epar[i],arglist[i]);
-          printf("Svertex: %s = %e +- %e (%e,%e)\n",
-                 chnam.Data(),par[i],epar[i],empar[i],eppar[i]);
-       } // end for i
-   } // End Muinuit fitting
-   su  = 0.0;
-   avt = 0.0;
-   i1  = Svzfl->GetBinContent(1);
-   i2  = Svzfl->GetBinContent(2);
-   for(i=2;i<nbinx;i++){
-       i0 = i1; // old i-1 value
-       i1 = i2; // old i value
-       i2 = Svzfl->GetBinContent(i+1); // new i+1 value
-       if(i1 > av && (i0 > av || i2 > av ) ){
-          if(!mod) mod = kTRUE;
-       } else if(mod) break;
-       if(mod){ // inside peak
-          su  += i1;
-          avt += i1 * Svzfl->GetXaxis()->GetBinCenter(i);
-       } // end if mod
-   } // end for i
-
-   if(su!=0.0) zv = avt/su;
-   else        zv = -100.0; // an unphysical value
-
-   TCanvas *c1 = new TCanvas("c1","Slow Alice ITS vertex finder",
-                            450,10,600,700);
-   Svzphi->Draw("col");
-   c1->Print("vertex5_S_vz_phi.eps");
-   Svzdphi->Draw();
-   c1->Print("vertex5_S_vz_dphi.eps");
-   Svzdef->Draw();
-   c1->Print("vertex5_S_vz_def.eps");
-   Svzbst->Draw();
-   c1->Print("vertex5_S_vz_bst.eps");
-   Svztr->Draw();
-   c1->Print("vertex5_S_vz_true.eps");
-
-   return zv;
-}
-
-Float_t vertexEta(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
-   Float_t r1    = 3.910078;// radius at which hit is from beam axis for spd1
-   Float_t r2    = 6.955933;// radius at which hit is from beam axis for spd2
-   Float_t DPhi  = 0.005;   // maximum allowed difference in phi angle
-   Float_t DEta  = 0.100;   // maximum allowed difference in eta/pseudorapidity
-   Float_t DZ    = 12.5;     // maximum allowed difference in z position
-   Float_t avt   = 0.0;
-   Int_t   nbinx = 2000;
-   Int_t   start = 0;
-   Bool_t  mod   = kFALSE;
-   Int_t   i,j;
-   Float_t zv=0.0,av,su,i0,i1,i2,x0,x1,dphi,dz;
-   Float_t PI2 = 2.0*TMath::Pi();
-   Float_t PI  = TMath::Pi();
-
-   // sort according to phi.
-   hpsortEta(spdi,i1max);
-   hpsortEta(spdo,i2max);
-
-   // find best vertex allong z.
-   TH2S *Evzphi   = new TH2S("Evzphi","Eta: Posible Z vertecies vs. Phi",
-                            nbinx,-5.0,5.0,32,0.0,PI2);
-   Evzphi->Sumw2(); // collect good statitics
-   TH2S *Evzdphi  = new TH2S("Evzdphi","Eta: Posible Z vertecies vs. DPhi",
-                            200,-1.0,1.0,20,0.0,10*DPhi);
-   Evzdphi->Sumw2(); // collect good statitics
-   TH1F *Evzfl   = new TH1F("Evzfl","Eta: Posible Z vertecies Flattened",
-                         nbinx,-5.0,5.0);
-   Evzfl->Sumw2();
-   TH1F *Evztr  = new TH1F("Evztr","Eta: Z Vertex found by True Tracks",
-                          nbinx,-5.0,5.0);
-   Evztr->Sumw2(); // collect good statitics
-
-   for(i=0;i<i1max;i++){
-       for(;spdo[start]->etar < spdi[i]->etar-1.5*DEta;start++);
-       for(j=start;(spdo[j]->etar < spdi[i]->etar+DEta) && (j < i2max);j++){
-          dphi = fabs(spdi[i]->phir - spdo[j]->phir); if(dphi>PI) dphi -= PI;
-          dz   = fabs(spdi[i]->zr   - spdo[j]->zr);
-          r1   = spdi[i]->r;
-          r2   = spdo[j]->r;
-          if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
-          else printf("vertex_Eta: spdi[%d]->r=%f = spdo[%d]->r=%f\n",
-                      i,r1,j,r2);
-          Evzdphi->Fill(zv,dphi);
-          if(dphi>DPhi) continue;
-          if(dz>DZ)     continue;
-          Evzphi->Fill(zv,spdi[i]->phir);
-          Evzfl->Fill(zv);
-          if(spdi[i]->track == spdo[j]->track) Evztr ->Fill(zv);
-       } // end for j
-   } // end for i
-
-   TH1D *Evzdef  = Evzphi->ProjectionX("Eta: Z vertecies Eta",1,nbinx,"E");
-
-   i0 = Evzfl->GetBinContent(1);
-   x0 = Float_t(1);
-   i1 = Evzfl->GetBinContent(nbinx);
-   x1 = Float_t(nbinx);
-   su = x0-x1; if(su==0.0) return -200.0;
-   su = (i0-i1)/su;
-   x0 = x0-su*i0;
-   x1 = su;
-   for(i=1;i<=nbinx;i++){
-      su = x1*Float_t(i) + x0;
-      Evzfl->AddBinContent(i,-su);
-   } // end for i
-
-   printf("Eta:           mean=%f RMS=%f w2=%f\n",
-         Evzdef->GetMean(),Evzdef->GetRMS(),Evzdef->GetSumOfWeights());
-   printf("Eta: Flattened mean=%f RMS=%f w2=%f\n",
-         Evzfl->GetMean(),Evzfl->GetRMS(),Evzfl->GetSumOfWeights());
-
-   av = 3.0 * Evzfl->Integral(1,nbinx)/Float_t(nbinx);
-   printf("Eta: Flattened av=%f TvxpE->Max=%f nbinx=%d\n",
-         av,Evzfl->GetMaximum(),nbinx);
-
-   su  = 0.0;
-   avt = 0.0;
-   i1  = Evzfl->GetBinContent(1);
-   i2  = Evzfl->GetBinContent(2);
-   for(i=2;i<nbinx;i++){
-       i0 = i1; // old i-1 value
-       i1 = i2; // old i value
-       i2 = Evzfl->GetBinContent(i+1); // new i+1 value
-       if(i1 > av && (i0 > av || i2 > av ) ){
-          if(!mod) mod = kTRUE;
-       } else if(mod) break;
-       if(mod){ // inside peak
-          su  += i1;
-          avt += i1 * Evzfl->GetXaxis()->GetBinCenter(i);
-       } // end if mod
-   } // end for i
-
-   if(su!=0.0) zv = avt/su;
-   else        zv = -100.0; // an unphysical value
-
-   TCanvas *c2 = new TCanvas("c2","Alice ITS vertex finder Eta", 
-                            500,10,600,700);
-   Evzphi->Draw("col");
-   c2->Print("vertex5_E_vz_phi.ps");
-   Evzdphi->Draw();
-   c2->Print("vertex5_E_vz_dphi.ps");
-   Evzdef->Draw();
-   c2->Print("vertex5_E_vz_def.ps");
-   Evzfl->Draw();
-   c2->Print("vertex5_E_vz_flat.ps");
-   Evztr->Draw();
-   c2->Print("vertex5_E_vz_true.ps");
-
-   return zv;
-}
diff --git a/ITS/vertex.h b/ITS/vertex.h
deleted file mode 100755 (executable)
index a562f6d..0000000
+++ /dev/null
@@ -1,19 +0,0 @@
-#ifndef VERTEX_H
-#define VERTEX_H
-struct Hit_tl {
-   Int_t   track,n,lad,det;
-   Float_t x,y,z,xl,yl,zl,xr,yr,zr;
-   Float_t r,phi,eta,phir,etar;
-};
-
-void FindVertexs(Int_t evnt,Float_t frac,Float_t len);
-
-void HitsToV(Hit_tl **spdi,Int_t &nspdi,Hit_tl **spdo,Int_t &nspdo,
-            Int_t ntracks,TTree *TH,TClonesArray *ITShits,
-            Float_t fraction,Float_t len);
-
-Float_t vertex(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max);
-Float_t vertexSlow(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max);
-Float_t vertexEta(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max);
-
-#endif