ITS/AnalysisITSAlignment.cxx

   1 //////////////////////////////////////////////////////////////////////////
   2 //  Alice ITS first detector alignment program.                         //
   3 //                                                                      //
   4 // version: 0.0.0 Draft.                                                //
   5 // Date: April 18 1999                                                  //
   6 // By: Bjorn S. Nilsen                                                  //
   7 //                                                                      //
   8 //////////////////////////////////////////////////////////////////////////
   9 #include <stdio.h>
  10 #include <math.h>
  11 #include <fstream.h>
  12 #include <stdio.h>
  13 #include <time.h>
  14 #include "AliITS.h"
  15 #include "TTree.h"
  16 #include "TClonesArray.h"
  17 #include "TProfile.h"
  18 #include "TH2.h"
  19 #include "TArray.h"
  20 #include "TRandom.h"
  21 #include "TCanvas.h"
  22 #include "TFile.h"
  23 #include "TParticle.h"
  24 #include "AliRun.h"
  25 #include "AliITSgeom.h"
  26 #include "AnalysisITSAlignment.h"
  27 #include "AliITSstatistics2.h"
  28 #include "AliITSstatistics.h"
  29 #include "AliITSAlignmentTrack.h"
  30 #include "AliITSAlignmentModule.h"
  31
  32 #define PRIMARYONLY 0  // if PIMARY_ONLY ==0 then all tracks
  33 #define MOMENTUMCUT 1.0 // Sets the momentum cut for the tracks.
  34
  35 //Double_t resz[]={0.00680,0.00650,0.00265,0.00265,0.0830,0.0830};
  36 //Double_t resx[]={0.00105,0.00115,0.00340,0.00340,0.0020,0.0020};
  37 // From table 1.3 of ITS TDR. Spatial precisions and Cell size
  38 Double_t  resz[]={0.00700,0.00700,0.00280,0.00280,0.08300,0.08300};
  39 Double_t  resx[]={0.00120,0.00120,0.00380,0.00380,0.00200,0.00200};
  40 Double_t cellx[]={0.00500,0.00500,0.01500,0.01500,0.00950,0.00950};
  41 Double_t celly[]={0.00150,0.00150,0.00300,0.00300,0.00300,0.00300};
  42 Double_t cellz[]={0.03000,0.03000,0.03000,0.03000,4.00000,4.00000};
  43
  44 //_________________________________________________________________________
  45 void HitsTo(ClustAl_tl *trk,Int_t &ntrk,Int_t nt,TTree *TH,AliITS *ITS,
  46             Float_t nsigmaT1,Float_t nsigmaT2,Float_t nsigmaT3,
  47             Float_t nsigmaR1,Float_t nsigmaR2,Float_t nsigmaR3){
  48
  49    // Local variables
  50    Int_t        nClustMin = 3;
  51    const Int_t  jmax = 50;
  52    Int_t        t,nb,nh,h,i,ly=-1,ld=-1,dt=-1,tt=-1,j=0;
  53    ClustAl_sl   ht[jmax];
  54    Int_t        id[3],ie[3],nht[jmax],track;
  55    Float_t      gp[3],lp[3],p[3],tof,gp2[3];
  56    Float_t      ns;
  57    Float_t      px[jmax],py[jmax],pz[jmax];
  58    TClonesArray *ITShits = ITS->Hits();
  59    AliITShit    *itsHit;
  60    // initilize and creat AliITSgeom data
  61    AliITSgeom   *gm  = ITS->GetITSgeom();
  62    AliITSgeom   &gm2 = *gm;
  63
  64 //   nt = 100; // lets speed things up.
  65
  66    Float_t tran[]={90.e-4/3.,30.e-4/3.,780.e-4/3.}; // r,rphi,z
  67    Float_t  rot[]={0.4e-3/3.,3.0e-3/3.,100.e-3/3.}; // about r, rphi, z
  68    tran[0] *= nsigmaT1;
  69    tran[1] *= nsigmaT2;
  70    tran[2] *= nsigmaT3;
  71    rot[0]  *= nsigmaR1;
  72    rot[1]  *= nsigmaR2;
  73    rot[2]  *= nsigmaR3;
  74
  75    printf("Randomize by tran=%fcm %fcm %fcm\n",tran[0],tran[1],tran[2]);
  76    printf("Randomize by rot=%frad %frad %frad\n",rot[0],rot[1],rot[2]);
  77
  78    gm2.RandomCylindericalChange(tran,rot);
  79
  80    for(i=0;i<jmax;i++){// zero for for first track.
  81        ht[i].xl = ht[i].xg = 0.0;
  82        ht[i].yl = ht[i].yg = 0.0;
  83        ht[i].zl = ht[i].zg = 0.0;
  84        ht[i].lay = ht[i].lad = ht[i].det = 0;
  85        px[i]    = 0.0;
  86        py[i]    = 0.0;
  87        pz[i]    = 0.0;
  88        nht[i]   = 0;
  89    } // end init to zero.
  90
  91    j = 0;
  92    for(t=0;t<nt;t++){ // Loop over tracks.
  93       gAlice->ResetHits();
  94       nb = TH->GetEvent(t);
  95       nh = ITShits->GetEntriesFast();
  96       for(h=0;h<nh;h++){ // hits
  97          itsHit = (AliITShit *)ITShits->UncheckedAt(h);
  98          track  = itsHit->fTrack;
  99          if(track != tt){ // if new track
 100              if(j>nClustMin)if(tt<nt||PRIMARYONLY==0){
 101                  // if there are enought points save data
 102                  trk[ntrk].track  = tt;
 103                  trk[ntrk].nclust = j+1;
 104                  ns               = 1./(Float_t (nht[0]));
 105                  trk[ntrk].px     = px[0]*ns;
 106                  trk[ntrk].py     = py[0]*ns;
 107                  trk[ntrk].pz     = pz[0]*ns;
 108                  trk[ntrk].pt     = sqrt(px[0]*px[0]+py[0]*py[0])*ns;
 109                  trk[ntrk].p      = sqrt(px[0]*px[0]+py[0]*py[0]+
 110                                          pz[0]*pz[0]            )*ns;
 111                  trk[ntrk].clust  = new ClustAl_sl[j+1];
 112                  for(i=0;i<=j;i++){ // loop over all detector hit.
 113                      ns                     = 1./(Float_t(nht[i]));
 114                      trk[ntrk].clust[i].lay = ie[0] = ht[i].lay;
 115                      trk[ntrk].clust[i].lad = ie[1] = ht[i].lad;
 116                      trk[ntrk].clust[i].det = ie[2] = ht[i].det;
 117                      lp[0]                  = ht[i].xl*ns;
 118                      lp[1]                  = ht[i].yl*ns;
 119                      lp[2]                  = ht[i].zl*ns;
 120                      // set in detector resolution on local coordiante system
 121                      // z->zl rphi->xl yl-> thickness leave it alone
 122                      lp[0] = lp[0]/resx[ie[0]-1];
 123                      lp[0] = (Float_t)((Int_t)lp[0]);
 124                      lp[0] = lp[0]*resx[ie[0]-1] + 0.5*resx[ie[0]-1];
 125                      lp[2] = lp[2]/resz[ie[0]-1];
 126                      lp[2] = (Float_t)((Int_t)lp[2]);
 127                      lp[2] = lp[2]*resz[ie[0]-1] + 0.5*resz[ie[0]-1];
 128                      //
 129                      trk[ntrk].clust[i].xl  = lp[0];
 130                      trk[ntrk].clust[i].yl  = lp[1];
 131                      trk[ntrk].clust[i].zl  = lp[2];
 132                      gm2.LtoG(ie,lp,gp);
 133                      gm->LtoG(ie,lp,gp2);
 134                      trk[ntrk].clust[i].xg  = gp[0];
 135                      trk[ntrk].clust[i].yg  = gp[1];
 136                      trk[ntrk].clust[i].zg  = gp[2];
 137                  } // end for i: detectors clust loop
 138                  ntrk++; // set up for next track
 139              } // end if j>nClustMin
 140              for(i=0;i<jmax;i++){ // zero out structure
 141                  ht[i].xl = ht[i].xg = 0.0;
 142                  ht[i].yl = ht[i].yg = 0.0;
 143                  ht[i].zl = ht[i].zg = 0.0;
 144                  ht[i].lay = ht[i].lad = ht[i].det = 0;
 145                  px[i]    = 0.0;
 146                  py[i]    = 0.0;
 147                  pz[i]    = 0.0;
 148                  nht[i]   = 0;
 149              } // end for i
 150              j  = -1; // since ly=ld=dt=0 forces a j++ call.
 151              ly = 0; // zero old detector values
 152              ld = 0;
 153              dt = 0;
 154              tt = track;
 155          } // end if track != tt
 156          itsHit->GetPositionG(gp[0],gp[1],gp[2],tof);
 157          itsHit->GetDetectorID(id[0],id[1],id[2]);
 158          gm->GtoL(id,gp,lp);
 159          itsHit->GetMomentumG(p[0],p[1],p[2]);
 160          if(!(id[0]==ly && id[1]==ld && id[2]==dt)){//if Not the same detector
 161              j++; // incriment detector counter j.
 162              ly = id[0];
 163              ld = id[1];
 164              dt = id[2];
 165          } // end if id==idold
 166          ht[j].lay = id[0];
 167          ht[j].lad = id[1];
 168          ht[j].det = id[2];
 169          ht[j].xl += lp[0];
 170          ht[j].yl += lp[1];
 171          ht[j].zl += lp[2];
 172          ht[j].xg  = gp[0];
 173          ht[j].yg  = gp[1];
 174          ht[j].zg  = gp[2];
 175          px[j]    += p[0];
 176          py[j]    += p[1];
 177          pz[j]    += p[2];
 178          nht[j]++;
 179       } // end for h
 180   } // end for t: track loop
 181   return;
 182 }
 183 //______________________________________________________________________
 184 void HitsToClustAl(ClustAl_tl *trk,Int_t &ntrk,Int_t nt,TTree *TH,AliITS *ITS,
 185                    Float_t fraction){
 186
 187    // Local variables
 188    Int_t        nClustMin = 3,Iseed;
 189 //   Int_t        icount=0,icountMAX=100;
 190    const Int_t  jmax = 50;
 191    Int_t        t,nb,nh,h,i,ly=-1,ld=-1,dt=-1,tt=-1,j=0;
 192    ClustAl_sl   ht[jmax];
 193    Int_t        id[3],ie[3],nht[jmax],track;
 194    Float_t      gp[3],lp[3],tof;
 195    Float_t      ns;
 196    TClonesArray *ITShits = ITS->Hits(),*Prt = gAlice->Particles();
 197    AliITShit    *itsHit;
 198    TParticle    *prt;
 199    // initilize and creat AliITSgeom data
 200    AliITSgeom   *gm  = ITS->GetITSgeom();
 201
 202    for(i=0;i<jmax;i++){// zero for for first track.
 203        ht[i].xl  = 0.0;
 204        ht[i].yl  = 0.0;
 205        ht[i].zl  = 0.0;
 206        ht[i].lay = ht[i].lad = ht[i].det = 0;
 207        nht[i]    = 0;
 208    } // end init to zero.
 209
 210    j = 0;
 211    Iseed = time(0);
 212    printf("HitsToclustAl: Iseed=%d ",Iseed);
 213    gRandom->SetSeed(Iseed);
 214    printf("gRandom->Rndm(1)=%f\n",gRandom->Rndm(1));
 215    for(t=0;t<nt;t++){ // Loop over tracks.
 216       if(fraction<gRandom->Rndm(1)) continue; // skip some tracks
 217       gAlice->ResetHits();
 218       nb = TH->GetEvent(t);
 219       nh = ITShits->GetEntriesFast();
 220       for(h=0;h<nh;h++){ // hits
 221          itsHit = (AliITShit *)ITShits->UncheckedAt(h);
 222          track  = itsHit->fTrack;
 223          prt= (TParticle *)Prt->UncheckedAt(track);
 224          if(prt->P()<MOMENTUMCUT) continue;
 225          if(track != tt){ // if new track
 226              if(j>nClustMin)if(tt<nt||PRIMARYONLY==0){
 227 //               if(icount>icountMAX) return;
 228 //               icount++;
 229                  // if there are enought points save data
 230                  prt              = (TParticle *)Prt->UncheckedAt(tt);
 231                  trk[ntrk].track  = tt;
 232                  trk[ntrk].nclust = j+1;
 233                  ns               = 1./(Float_t (nht[0]));
 234                  trk[ntrk].px     = prt->Px();
 235                  trk[ntrk].py     = prt->Py();
 236                  trk[ntrk].pz     = prt->Pz();
 237                  trk[ntrk].pt     = prt->Pt();
 238                  trk[ntrk].p      = prt->P();
 239                  trk[ntrk].clust  = new ClustAl_sl[j+1];
 240                  for(i=0;i<=j;i++){ // loop over all detector hit.
 241                      ns                     = 1./(Float_t(nht[i]));
 242                      trk[ntrk].clust[i].lay = ie[0] = ht[i].lay;
 243                      trk[ntrk].clust[i].lad = ie[1] = ht[i].lad;
 244                      trk[ntrk].clust[i].det = ie[2] = ht[i].det;
 245                      lp[0]                  = ht[i].xl*ns;
 246                      lp[1]                  = ht[i].yl*ns;
 247                      lp[2]                  = ht[i].zl*ns;
 248                      // set in detector resolution on local coordiante system
 249                      // z->zl rphi->xl yl-> thickness leave it alone
 250                      lp[0] = lp[0]/resx[ie[0]-1];
 251                      lp[0] = (Float_t)((Int_t)lp[0]);
 252                      lp[0] = lp[0]*resx[ie[0]-1] + 0.5*resx[ie[0]-1];
 253                      lp[2] = lp[2]/resz[ie[0]-1];
 254                      lp[2] = (Float_t)((Int_t)lp[2]);
 255                      lp[2] = lp[2]*resz[ie[0]-1] + 0.5*resz[ie[0]-1];
 256                      //
 257                      trk[ntrk].clust[i].xl  = lp[0];
 258                      trk[ntrk].clust[i].yl  = lp[1];
 259                      trk[ntrk].clust[i].zl  = lp[2];
 260                  } // end for i: detectors clust loop
 261                  ntrk++; // set up for next track
 262              } // end if j>nClustMin
 263              for(i=0;i<jmax;i++){ // zero out structure
 264                  ht[i].xl = ht[i].xg = 0.0;
 265                  ht[i].yl = ht[i].yg = 0.0;
 266                  ht[i].zl = ht[i].zg = 0.0;
 267                  ht[i].lay = ht[i].lad = ht[i].det = 0;
 268                  nht[i]   = 0;
 269              } // end for i
 270              j  = -1; // since ly=ld=dt=0 forces a j++ call.
 271              ly = 0; // zero old detector values
 272              ld = 0;
 273              dt = 0;
 274              tt = track;
 275          } // end if track != tt
 276          itsHit->GetPositionG(gp[0],gp[1],gp[2],tof);
 277          itsHit->GetDetectorID(id[0],id[1],id[2]);
 278          gm->GtoL(id,gp,lp);
 279          if(!(id[0]==ly && id[1]==ld && id[2]==dt)){//if Not the same detector
 280              j++; // incriment detector counter j.
 281              ly = id[0];
 282              ld = id[1];
 283              dt = id[2];
 284          } // end if id==idold
 285          ht[j].lay = id[0];
 286          ht[j].lad = id[1];
 287          ht[j].det = id[2];
 288          ht[j].xl += lp[0];
 289          ht[j].yl += lp[1];
 290          ht[j].zl += lp[2];
 291          nht[j]++;
 292       } // end for h
 293   } // end for t: track loop
 294   return;
 295 }
 296 //______________________________________________________________________
 297 void SetDetectorResolusion(Int_t l,Float_t *xl){
 298     //set in detector resolution on local coordiante system
 299     // z->zl rphi->xl yl-> thickness leave it alone
 300     // changed from res to cell for detector simulations.
 301
 302     l--;
 303     xl[0] = xl[0]/cellx[l];
 304     xl[0] = (Float_t)((Int_t)xl[0]);
 305     xl[0] = xl[0]*cellx[l] + 0.5*cellx[l];
 306     xl[2] = xl[2]/cellz[l];
 307     xl[2] = (Float_t)((Int_t)xl[2]);
 308     xl[2] = xl[2]*cellz[l] + 0.5*cellz[l];
 309     xl[1] = 0.0;
 310 }
 311 //______________________________________________________________________
 312 void FillAliITSAlignmentTrack(AliITSAlignmentTrack *trk,Int_t &ntrk,Int_t nt,
 313                               TTree *TH,AliITS *ITS,Float_t fraction){
 314     const Int_t iMAX=20,jMAX=50,nClustMin=3;
 315     Int_t   ist[iMAX],imax;
 316     Int_t   t,i,h,j,nh;
 317     Int_t   lay,lad,det,index,track,Iseed;
 318     Int_t   indexl[jMAX];
 319     Float_t tof,xl[3],exl[3][3];
 320     AliITSstatistics *sxl[jMAX],*szl[jMAX];
 321     TParticle    *prt;
 322     AliITShit    *itshit;
 323     TClonesArray *ITShits = ITS->Hits();
 324     AliITSgeom   *gm = ITS->GetITSgeom();
 325     TClonesArray *Prt = gAlice->Particles();
 326
 327     printf("Entered FillAliITSAlignmentTrack\n");
 328     for(i=0;i<jMAX;i++){
 329         sxl[i] = new AliITSstatistics(2);
 330         szl[i] = new AliITSstatistics(2);
 331     }  // end for i
 332     for(i=0;i<3;i++)for(j=0;j<3;j++) exl[i][j] = 0.0;
 333     Iseed = time(0);
 334     gRandom->SetSeed(Iseed);
 335     printf("Entering track loop. Iseed=%d\n",Iseed);
 336     for(t=0;t<nt;t++) if(fraction>=gRandom->Rndm(1)){
 337 //      gAlice->ResetHits();
 338         TH->GetEvent(t);
 339         nh = ITShits->GetEntriesFast();
 340         i = 0;
 341         ist[0] = 0;
 342         printf("Entering hit loop for track=%d\n",t);
 343         for(h=1;h<nh&&i<iMAX-1;h++)
 344             if(((AliITShit *)ITShits->UncheckedAt(ist[i]))->GetTrack() !=
 345                ((AliITShit *)ITShits->UncheckedAt(h))->GetTrack()) ist[++i]=h;
 346         ist[++i]=nh;
 347         imax = i+1;
 348         printf("entering loop i nh=%d imax=%d\n",nh,imax);
 349         for(i=1;i<imax;i++){ // loop over tracks from primary track
 350             printf("Getting hit ist[%d-1]=%d\n",i,ist[i-1]);
 351             itshit = (AliITShit *)ITShits->UncheckedAt(ist[i-1]);
 352             track  = itshit->GetTrack();
 353             prt    = (TParticle *)Prt->UncheckedAt(track);
 354             if(prt->P()<MOMENTUMCUT) continue;
 355             printf("pass P cut P=%f\n",prt->P());
 356             for(j=0;j<jMAX;j++){
 357                 sxl[j]->Reset();
 358                 szl[j]->Reset();
 359             } // end for j
 360             printf("Exiting the reset loop, ist[%d-1]=%d\n",i,ist[i-1]);
 361             j = 0;
 362             indexl[0]=((AliITShit *)ITShits->UncheckedAt(ist[i-1]))->
 363                                                                 GetDetector();
 364             printf("looping over this track's hits ist[%d-1]=%d to ist[%d]=%d\n",i,ist[i-1],i,ist[i]);
 365             for(h=ist[i-1];h<ist[i];h++){ // loop over hits
 366                 itshit = (AliITShit *)ITShits->UncheckedAt(h);
 367                 index = itshit->GetDetector();
 368                 if(indexl[j]!=index) j++;
 369                 indexl[j] = index;
 370                 itshit->GetPositionL(xl[0],xl[1],xl[2],tof);
 371                 sxl[j]->AddValue((Double_t) xl[0],1.0);
 372                 szl[j]->AddValue((Double_t) xl[2],1.0);
 373             } // end for h
 374             if(j<nClustMin) continue; // get next track
 375             printf("Setting up tracks ntrk=%d #digits for track=%d\n",ntrk,j+1);
 376             trk[ntrk].CreatePoints(j+1);
 377             printf("CreatePoints with ntrk=%d and j+1=%d\n",ntrk,j+1);
 378             trk[ntrk].SetTParticle(prt);
 379             printf("SetTParticle\n");
 380             trk[ntrk].SetTrackNumber(track);
 381             printf("SetTrackNumber=%d\n",track);
 382             for(h=0;h<=j;h++){
 383                 xl[0] = sxl[h]->GetMean();
 384                 xl[1] = 0.0;
 385                 xl[2] = szl[h]->GetMean();
 386                 gm->GetModuleId(indexl[h],lay,lad,det);
 387                 printf("setting detector resolusion for index[%d]=%d (%d,%d,%d)\n",h,indexl[h],lay,lad,det);
 388                 SetDetectorResolusion(lay,xl);
 389                 printf("Detector resolusion set\n");
 390                 exl[0][0] = 1./(resx[lay-1]*resx[lay-1]);
 391                 exl[1][1] = 12./(celly[lay-1]*celly[lay-1]);
 392                 exl[2][2] = 1./(resz[lay-1]*resz[lay-1]);
 393                 printf("Adding point to trk[%d] indexl[%d]=%d\n",ntrk,h,indexl[h]);
 394                 trk[ntrk].AddPointLastL(indexl[h],(Float_t *)xl,
 395                                                              (Float_t **)exl);
 396             } // end for h
 397             ntrk++; // set up for next track
 398         } // end for i
 399     } // end for t // end if(fraction>=gRandom->Rndm(1))
 400     for(i=0;i<jMAX;i++){
 401         delete sxl[i];
 402         delete szl[i];
 403     }  // end for i
 404 }
 405 //______________________________________________________________________
 406 void FillGlobalPositions(ClustAl_tl *trk,Int_t ntrk,AliITSgeom *g){
 407     Int_t   i,j,id[3];
 408     Float_t lx[3],gx[3];
 409
 410     for(i=0;i<ntrk;i++) for(j=0;j<trk[i].nclust;j++){
 411         lx[0] = trk[i].clust[j].xl;
 412         lx[1] = trk[i].clust[j].yl;
 413         lx[2] = trk[i].clust[j].zl;
 414         id[0] = trk[i].clust[j].lay;
 415         id[1] = trk[i].clust[j].lad;
 416         id[2] = trk[i].clust[j].det;
 417         g->LtoG(id,lx,gx);
 418         trk[i].clust[j].xg = gx[0];
 419         trk[i].clust[j].yg = gx[1];
 420         trk[i].clust[j].zg = gx[2];
 421     } // end for i,j
 422     return;
 423 }
 424 //______________________________________________________________________
 425 void PlotGeomChanges(AliITSgeom *gt,AliITSgeom *gc,TFile *Hfile,Float_t *Rdta){
 426     Int_t    ly,ld,dt;
 427     Float_t  x,y,z;
 428     Double_t xd,yd,zd;
 429     Double_t rt,rc,phit,phic,zt,zc;
 430     Float_t  A0,B0,C0,A1,B1,C1;
 431     Double_t dr,drphi,dz,dphi;
 432     Double_t drmi,drma,drphimi,drphima,dzmi,dzma;
 433
 434     AliITSstatistics *R    = new AliITSstatistics(2);
 435     AliITSstatistics *RPhi = new AliITSstatistics(2);
 436     AliITSstatistics *Z    = new AliITSstatistics(2);
 437     AliITSstatistics *A    = new AliITSstatistics(2);
 438     AliITSstatistics *B    = new AliITSstatistics(2);
 439     AliITSstatistics *C    = new AliITSstatistics(2);
 440     TH1F *Gr  = new TH1F("Gr","Radial Displacement (cm)",500,-1.0,1.0);
 441     drmi = -1.0;
 442     drma = +1.0;
 443     Gr->Sumw2();
 444     Gr->SetXTitle("Displacement (cm)");
 445     TH1F *Grphi  = new TH1F("Grphi","RPhi Displacement (cm)",500,-1.0,1.0);
 446     drphimi = -1.0;
 447     drphima = +1.0;
 448     Grphi->Sumw2();
 449     Grphi->SetXTitle("Displacement (cm)");
 450     TH1F *Gz  = new TH1F("Gz","Z Displacement (cm)",500,-1.0,1.0);
 451     dzmi = -1.0;
 452     dzma = +1.0;
 453     Gz->Sumw2();
 454     Gz->SetXTitle("Displacement (cm)");
 455
 456     Float_t weight=1.0;
 457     for(ly=1;ly<=gt->GetNlayers();ly++)for(ld=1;ld<=gt->GetNladders(ly);ld++)
 458     for(dt=1;dt<=gt->GetNdetectors(ly);dt++){
 459         gt->GetTrans(ly,ld,dt,x,y,z);xd=x;yd=y;zd=z;
 460         gt->GetAngles(ly,ld,dt,A0,B0,C0);
 461         rt    = TMath::Hypot(yd,xd);
 462         phit  = TMath::ATan2(yd,xd);
 463         if(phit<0.0) phit += 2.0*TMath::Pi();
 464         zt    = zd;
 465         gc->GetTrans(ly,ld,dt,x,y,z);xd=x;yd=y;zd=z;
 466         gc->GetAngles(ly,ld,dt,A1,B1,C1);
 467         rc    = TMath::Hypot(yd,xd);
 468         phic  = TMath::ATan2(yd,xd);
 469         if(phic<0.0) phic += 2.0*TMath::Pi();
 470         zc    = zd;
 471         dr    = rt-rc;
 472         dphi  = phit - phic;
 473         if(dphi>2.0*TMath::Pi()) dphi -= 2.0*TMath::Pi();
 474         drphi = 0.5*(rt + rc)*dphi; // change in phi as measured in rphi
 475         dz    = zt - zc;
 476         Gr->Fill(dr,1.0);
 477         Grphi->Fill(drphi,1.0);
 478         Gz->Fill(dz,1.0);
 479         if(dr>=drmi&&dr<drma) R->AddValue(dr,1.0);
 480         if(drphi>=drphimi&&drphi<drphima) RPhi->AddValue(drphi,1.0);
 481         if(dz>=dzmi&&dz<dzma) Z->AddValue(dz,1.0);
 482         A->AddValue(A0-A1,weight);
 483         B->AddValue(B0-B1,weight);
 484         C->AddValue(C0-C1,weight);
 485     } // end for ly,ld,dt
 486     Hfile->Write();
 487     Rdta[0] = R->GetRMS();
 488     Rdta[1] = RPhi->GetRMS();
 489     Rdta[2] = Z->GetRMS();
 490     Rdta[3] = A->GetRMS();
 491     Rdta[4] = B->GetRMS();
 492     Rdta[5] = C->GetRMS();
 493     printf("PlotGeomChanges: RMS(r) = %f RMS(rphi) = %f RMS(z) = %f\n",
 494            R->GetRMS(),RPhi->GetRMS(),Z->GetRMS());
 495     delete Gr;
 496     delete Grphi;
 497     delete Gz;
 498     delete R;
 499     delete RPhi;
 500     delete Z;
 501     delete A;
 502     delete B;
 503     delete C;
 504     return;
 505 }
 506 //______________________________________________________________________
 507 Int_t FitTrackToLine(ClustAl_tl &trk){
 508    // Local Variables
 509    Int_t   i,l,j=0;
 510    Double_t x,y,z,wx,wy,dx,dy;
 511    Double_t a,b,c,d;
 512    Double_t xb,sxb,zxb;
 513    Double_t yb,syb,zyb;
 514    Double_t xzb,yzb,z2xb,z2yb;
 515    Double_t sum,phi=0.0,b0,d0;
 516
 517    trk.qual = -1.0;
 518    if(trk.nclust<3) return -1;
 519
 520    b  = (trk.clust[0].xg-trk.clust[trk.nclust].xg)/
 521         (trk.clust[0].zg-trk.clust[trk.nclust].zg);
 522    d  = (trk.clust[0].yg-trk.clust[trk.nclust].yg)/
 523         (trk.clust[0].zg-trk.clust[trk.nclust].zg);
 524    do{
 525      b0 = b;
 526      d0 = d;
 527      xb  = 0.0; sxb = 0.0; zxb  = 0.0;
 528      yb  = 0.0; syb = 0.0; zyb  = 0.0;
 529      xzb = 0.0; yzb = 0.0; z2xb = 0.0; z2yb = 0.0;
 530      for(i=0;i<trk.nclust;i++){
 531       l    = trk.clust[i].lay - 1;  // zero based
 532       x    = trk.clust[i].xg;
 533       y    = trk.clust[i].yg;
 534       z    = trk.clust[i].zg;
 535       phi  = atan2(y,x);
 536       wx   = 1.0/(resx[l]*resx[l]*cos(phi)*cos(phi)+b0*b0*resz[l]*resz[l]);// 1.0/rms^2
 537       wy   = 1.0/(resx[l]*resx[l]*sin(phi)*sin(phi)+d0*d0*resz[l]*resz[l]);// 1.0/rms^2
 538       xb   += x*wx;
 539       sxb  += wx;
 540       zxb  += z*wx;
 541       yb   += y*wy;
 542       syb  += wy;
 543       zyb  += z*wy;
 544       xzb  += x*z*wx;
 545       yzb  += y*z*wy;
 546       z2xb += z*z*wx;
 547       z2yb += z*z*wy;
 548      } // end for i
 549      dx = zxb*zxb - z2xb*sxb;
 550      if(dx!=0.0){
 551          a  = zxb*xzb - xb*z2xb;
 552          a  = a/dx;
 553          b  = xb*zxb  - sxb*xzb;
 554          b  = b/dx;
 555      }else{
 556           printf("FitTrackToLine: Error dx=0 track=%d zxb=%f z2xb=%f sxb=%f\n",
 557                  trk.track,zxb,z2xb,sxb);
 558           a = 0.0;
 559           b = 1E10;
 560           j++;
 561           if(j>100) return(0);
 562      } // end if dx!=0.0
 563      dy = zyb*zyb - z2yb*syb;
 564      if(dy!=0.0){
 565          c  = zyb*yzb - yb*z2yb;
 566          c  = c/dy;
 567          d  = yb*zyb  - syb*yzb;
 568          d  = d/dy;
 569      }else{
 570           printf("FitTrackToLine: Error dy=0 track=%d zyb=%f z2yb=%f syb=%f\n",
 571                  trk.track,zyb,z2yb,syb);
 572           c = 0.0;
 573           d = 1.E10;
 574           j++;
 575           if(j>100) return(0);
 576      } // end if dy!=0.0
 577    } while(fabs(b0-b)<1.E-5 && fabs(d0-d)<1.E-5);
 578    trk.a = a;
 579    trk.b = b;
 580    trk.c = c;
 581    trk.d = d;
 582    sum   = 0.0;
 583    for(i=0;i<trk.nclust;i++){
 584       l  = trk.clust[i].lay - 1;  // zero based
 585       x  = trk.clust[i].xg;
 586       y  = trk.clust[i].yg;
 587       z  = trk.clust[i].zg;
 588       x  = x - (a + b*z);  // x=a+bz
 589       y  = y - (c + d*z);  // y=c+dz
 590       wx =  resx[l]*resx[l]*cos(phi)*cos(phi) + b*b * resz[l]*resz[l];
 591       wy =  resx[l]*resx[l]*sin(phi)*sin(phi) + d*d * resz[l]*resz[l];
 592       if(wx==0.0||wy==0.0) {
 593           printf("FitTrackToLine: error wx=%f wy=%f trk.clust[%d].lay=%d\n",
 594                  wx,wy,i,trk.clust[i].lay);
 595       }// end if sxz or syz == 0.
 596       sum += x*x/wx + y*y/wy;
 597    } // end for i #2
 598    trk.qual = sum;
 599    if(trk.nclust>2) trk.qual /= Double_t (trk.nclust-2);
 600    // per degree of freedom 2 in xz 2 in yz.
 601    return 0;
 602 }
 603 //______________________________________________________________________
 604 void LtoLline(const Int_t *id1,const Int_t *id2,
 605               Double_t a1,Double_t b1,Double_t c1,Double_t d1,
 606               Double_t &a2,Double_t &b2,Double_t &c2,Double_t &d2,
 607               AliITSgeom *gm){
 608     Double_t x11[3],x12[3],x21[3],x22[3],h;
 609
 610     x11[0] = a1;
 611     x11[1] = 0.0;
 612     x11[2] = c1;
 613     x12[0] = a1+b1;
 614     x12[1] = 1.0;
 615     x12[2] = c1+d1;
 616     gm->LtoL(id1,id2,x11,x21);
 617     gm->LtoL(id1,id2,x12,x22);
 618     h = x21[1] - x22[1];
 619     if(h!=0.0){
 620         b2 = (x21[0] - x22[0])/h;
 621         d2 = (x21[2] - x22[2])/h;
 622         a2 = x21[0] - b2*x21[1];
 623         c2 = x21[2] - d2*x21[1];
 624     }else{
 625         printf("LtoLline: error line in plane of detector\n");
 626     } // end if h!=0
 627     return;
 628 }
 629 //______________________________________________________________________
 630 Int_t FitTrackToLineL(ClustAl_tl &trk,AliITSgeom *gm){
 631    // Local Variables
 632    Int_t   i,id0[3],id[3];
 633    Double_t wx/*,wy*/,wz;
 634    Double_t a,b,c,d;
 635    Double_t xg[3],xl[3],x2g[3],x2l[3];
 636    AliITSstatistics2 *Fx  = new AliITSstatistics2(2);
 637    AliITSstatistics2 *Fz  = new AliITSstatistics2(2);
 638
 639    trk.qual = -1.0;
 640    if(trk.nclust<3) return -1;
 641
 642    Int_t Npts = trk.nclust;
 643    id0[0] = trk.clust[0].lay;
 644    id0[1] = trk.clust[0].lad;
 645    id0[2] = trk.clust[0].det;
 646    for(i=0;i<Npts;i++){
 647        id[0] = trk.clust[i].lay;
 648        id[1] = trk.clust[i].lad;
 649        id[2] = trk.clust[i].det;
 650        xg[0] = trk.clust[i].xg;
 651        xg[1] = trk.clust[i].yg;
 652        xg[2] = trk.clust[i].zg;
 653        gm->GtoL(id0,xg,xl);
 654        wx = 1.0/(resx[id[0]-1]*resx[id[0]-1]);
 655        wz = 1.0/(resz[id[0]-1]*resz[id[0]-1]);
 656        Fx->AddValue(xl[0],xl[1],wx);
 657        Fz->AddValue(xl[2],xl[1],wz);
 658    } // end for i
 659    trk.qual  = Fx->FitToLine(a,b);
 660    trk.qual += Fz->FitToLine(c,d);
 661    trk.a0 = a;
 662    trk.b0 = b;
 663    trk.c0 = c;
 664    trk.d0 = d;
 665    xl[0]  = a;
 666    xl[1]  = 0.0;
 667    xl[2]  = c;
 668    x2l[0] = a+b;
 669    x2l[1] = 1.0;
 670    x2l[2] = c+d;
 671    gm->LtoG(id0,xl,xg);
 672    gm->LtoG(id0,x2l,x2g);
 673    c = xg[2] - x2g[2];
 674    if(c!=0.0){
 675        b = (xg[0] - x2g[0])/c;
 676        d = (xg[1] - x2g[1])/c;
 677        a = xg[0] - b*xg[2];
 678        c = xg[1] - d*xg[2];
 679        trk.a = a;
 680        trk.b = b;
 681        trk.c = c;
 682        trk.d = d;
 683    }else{
 684        return -1;
 685    }// end if c!=0.0
 686    return 0;
 687 }
 688 //______________________________________________________________________
 689 Int_t FindCircleCenter(Double_t &x0,Double_t &y0,Double_t x1,Double_t y1,
 690                        Double_t  x2,Double_t  y2,Double_t x3,Double_t y3){
 691 ////////////////////////////////////////////////////////////////////////
 692 //     This was derived as folows. Given three non-linear points find
 693 // the circle that is therefor defined by those three non-linar points.
 694 // Assume that the circle is centers at x0,y0 and has a radous R. Then
 695 // (1) R^2 = (x1-x0)^2 + (y1-y0)^2
 696 // (2) R^2 = (x2-x0)^2 + (y2-y0)^2
 697 // (3) R^2 = (x3-x0)^2 + (y3-y0)^2.
 698 // Now consider the two equations derived from the above
 699 // (1) - (2) = x1^2 - x2^2 -2x0(x1-x2) + y1^2 - y2Y2 -2Y0(y1-y2) = 0
 700 // (1) - (3) = x1^2 - x3^2 -2x0(x1-x3) + y1^2 - y3Y2 -2Y0(y1-y3) = 0
 701 // solving these two equations for x0 and y0 gives
 702 // x0 = +{(y1-y2)*(y1-y3)*(y1-y3)+x1*x1*(y1-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2)}/2d
 703 // y0 = -{(x1-x2)*(x1-x3)*(x1-x3)+y1*y1*(x1-x3)+y2*y2*(x3-x1)+y3*y3*(x1-x2)}/2d
 704 // with d = (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2)
 705 ////////////////////////////////////////////////////////////////////////
 706     Double_t d;
 707
 708     d = (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2);
 709     if(d==0.0) return 0;  // fits to a line!
 710
 711     x0 = (y1-y2)*(y1-y3)*(y1-y3)+x1*x1*(y1-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2);
 712     y0 = (x1-x2)*(x1-x3)*(x1-x3)+y1*y1*(x1-x3)+y2*y2*(x3-x1)+y3*y3*(x1-x2);
 713     x0 = +0.5*x0/d;
 714     y0 = -0.5*y0/d;
 715
 716     return 1;
 717 }
 718 //______________________________________________________________________
 719 void FitAllTracks(ClustAl_tl *trk,Int_t ntrk,Float_t *v0,AliITSgeom *gm,
 720                   const char *sfile,TFile *Hfile,Float_t *Fdta,Int_t *Ndta){
 721    // Local Variables
 722    Int_t    i,j,k,id0[3],id[3];
 723    Double_t xm,ym,zm,dt,ad/*,rh,rp,rphih,rphip*/;
 724    Double_t tp,tx,rx,rz;
 725    Double_t a,b,c,d,qual;
 726    Double_t xl[3],xg[3];
 727    Double_t trkqualMAX = 50.0;
 728    Bool_t   printit = kFALSE;
 729    char     filename[80],hid[10];
 730    Int_t    Nqualmax = 10,NqualgtMAX=0;
 731    Int_t    Nqualless[10];
 732 // Distrobution statisitics
 733    AliITSstatistics *Sad    = new AliITSstatistics(4);
 734    AliITSstatistics *SrxzL  = new AliITSstatistics(4);
 735    AliITSstatistics *SrxL[6];
 736    AliITSstatistics *SrzL[6];
 737    for(i=0;i<6;i++) {
 738        SrxL[i] = new AliITSstatistics(4);
 739        SrzL[i] = new AliITSstatistics(4);
 740    } // end for i
 741
 742    for(i=0;i<Nqualmax;i++) Nqualless[i] = 0;
 743
 744 // Setup histograms
 745
 746    TH2F *Ttqp  = new TH2F("Ttqp","Track quality vs. momentum",
 747                          500,0.0,50.0,500,0.0,10.0);
 748    Ttqp->Sumw2();
 749    Ttqp->SetXTitle("Chi^2/degree freedom");
 750    Ttqp->SetYTitle("Momentum GeV/c");
 751    TH2F *Tdttq = new TH2F("Tdttq",
 752                           "Distance from true vertex vs. track quality",
 753                           500,0.0,0.10,100,0.0,50.0);
 754    Tdttq->Sumw2();
 755    Tdttq->SetXTitle("Distanct to true vertex (cm)");
 756    Tdttq->SetYTitle("Chi^2/degree freedom");
 757    TH2F *Tadtq = new TH2F("Tadtq",
 758                      "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. track quality",
 759                           500,-0.03,0.03,500,0.0,50.0);
 760    Tadtq->Sumw2();
 761    Tadtq->SetXTitle("theta fitted line - atan(pt/p) (rad)");
 762    Tadtq->SetYTitle("Chi^2/degree freedom");
 763    TH2F *Taddt = new TH2F("Taddt",
 764          "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. Distance from true vertex",
 765                           500,-0.03,0.03,200,0.0,0.10);
 766    Taddt->Sumw2();
 767    Taddt->SetXTitle("theta fitted line - atan(pt/p) (rad)");
 768    Taddt->SetYTitle("distance to true vertex (cm)");
 769    TH2F *Tadp = new TH2F("Tadp",
 770          "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. momentum",
 771                           500,-0.03,0.03,200,0.0,10.0);
 772    Tadp->Sumw2();
 773    Tadp->SetXTitle("theta fitted line - atan(pt/p) (rad)");
 774    Tadp->SetYTitle("Momentum (GeV/c)");
 775
 776    TH2F *TrxrzL[7];
 777 //   TH2F *TrrprzG[7];
 778    for(i=0;i<7;i++){
 779        sprintf(filename,"Layer %1.1d: xi-x(i) vs. zi-z(i) local",i+1);
 780        if(i==6) sprintf(filename,"Sum of all layers: xi-x(i) "
 781                         "vs. zi-z(i) local");
 782        sprintf(hid,"TrxrzL%1.1d",i+1);
 783        TrxrzL[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
 784        TrxrzL[i]->Sumw2();
 785        TrxrzL[i]->SetXTitle("xi-x(i) (cm) local");
 786        TrxrzL[i]->SetYTitle("zi-z(i) (cm) local");
 787 //
 788 //       sprintf(filename,"Layer %1.1d: rPhii-rPhi(i) vs. zi-z(i) local",i+1);
 789 //       if(i==6) sprintf(filename,"Sum of all layers: rPhii-rPhi(i) "
 790 //                      "vs. zi-z(i) local");
 791 //       sprintf(hid,"TrrprzG%1.1d",i+1);
 792 //       TrrprzG[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
 793 //       TrrprzG[i]->Sumw2();
 794 //       TrrprzG[i]->SetXTitle("rPhii-rPhi(i) (cm) local");
 795 //       TrrprzG[i]->SetYTitle("zi-z(i) (cm) local");
 796    } // end for i
 797
 798    TH1D *Tptqp[10];
 799    TH1D *Ttqall = new TH1D("Ttqall","Track quality for all tracks fit",
 800                            500,0.0,1000.0);
 801    Ttqall->Sumw2();
 802    Ttqall->SetXTitle("Track quality: Chi squared per degree freedom");
 803
 804 // Fit each track and fill histograms and the like.
 805
 806    for(i=0;i<ntrk;i++){
 807        if(trk[i].p < MOMENTUMCUT) continue;
 808       if(FitTrackToLineL(trk[i],gm)!=0) continue;  // fit track to a line
 809       qual = trk[i].qual;
 810       if(qual<0.0) continue;
 811       // Fill histograms and statistics before cource chi squared cut
 812       Ttqall->Fill(qual,1.0);
 813       for(j=0;j<Nqualmax;j++) if(qual<(Double_t) j+1) Nqualless[j]++;
 814 //      if(qual>trkqualMAX) {
 815 //        NqualgtMAX++;
 816 //        continue;
 817 //      } // end if iqual>trkqualMAX=50.0
 818       a = trk[i].a;
 819       b = trk[i].b;
 820       c = trk[i].c;
 821       d = trk[i].d;
 822       zm  = b*(v0[0]-a) + d*(v0[1]-c);
 823       dt  = (b*b + d*d);
 824       if(dt!=0.0) zm /= dt;
 825       else{
 826           printf("FitAllTracks: trk[%d].b=%f trk[%d].d=%f\n",i,b,i,d);
 827           zm=0.0;
 828       }// end if else dt!=0.0
 829       xm  = a + b*zm;
 830       ym  = c + d*zm;
 831       dt   = sqrt((xm-v0[0])*(xm-v0[0]) + (ym-v0[1])*(ym-v0[1]) +
 832                  (zm-v0[2])*(zm-v0[2]));
 833       tx  = atan(sqrt(b*b + d*d));
 834       tp  = atan2(trk[i].pt,fabs(trk[i].pz));
 835       ad  = tx-tp;
 836       // Fill histograms and statistics
 837       Tadtq->Fill(ad,qual,1.0);
 838       Sad->AddValue(ad,1.0);
 839       Taddt->Fill(ad,dt,1.0);
 840       Tadp->Fill(ad,trk[i].p,1.0);
 841       Tdttq->Fill(dt,qual,1.0);
 842       Ttqp->Fill(qual,trk[i].p,1.0);
 843       //
 844       id0[0] = trk[i].clust[0].lay;
 845       id0[1] = trk[i].clust[0].lad;
 846       id0[2] = trk[i].clust[0].det;
 847       xg[0] = trk[i].clust[0].xg;
 848       xg[1] = trk[i].clust[0].yg;
 849       xg[2] = trk[i].clust[0].zg;
 850       gm->GtoL(id0,xg,xl);
 851       rx = xl[0] - trk[i].a0 - trk[i].b0 * xl[1];
 852       rz = xl[2] - trk[i].c0 - trk[i].d0 * xl[1];
 853       SrxL[id0[0]-1]->AddValue(rx,1.0);
 854       SrzL[id0[0]-1]->AddValue(rz,1.0);
 855       SrxzL->AddValue(rx,1.0);
 856       SrxzL->AddValue(rz,1.0);
 857       TrxrzL[6]->Fill(rx,rz,1.0);
 858       TrxrzL[id0[0]-1]->Fill(rx,rz,1.0);
 859       for(j=1;j<trk[i].nclust;j++){
 860           id[0] = trk[i].clust[j].lay;
 861           id[1] = trk[i].clust[j].lad;
 862           id[2] = trk[i].clust[j].det;
 863           LtoLline(id0,id,trk[i].a0,trk[i].b0,trk[i].c0,trk[i].d0,a,b,c,d,gm);
 864           xg[0] = trk[i].clust[j].xg;
 865           xg[1] = trk[i].clust[j].yg;
 866           xg[2] = trk[i].clust[j].zg;
 867           gm->GtoL(id,xg,xl);
 868           rx = xl[0] - a - b * xl[1];
 869           rz = xl[2] - c - d * xl[1];
 870           SrxL[id[0]-1]->AddValue(rx,1.0);
 871           SrzL[id[0]-1]->AddValue(rz,1.0);
 872           SrxzL->AddValue(rx,1.0);
 873           SrxzL->AddValue(rz,1.0);
 874           TrxrzL[6]->Fill(rx,rz,1.0);
 875           TrxrzL[id[0]-1]->Fill(rx,rz,1.0);
 876       } // end for j
 877    } // end for i
 878
 879 // Write out information
 880
 881    for(i=0;i<6;i++){
 882        Fdta[4*i+0] = SrxL[i]->GetRMS();
 883        Fdta[4*i+1] = SrxL[i]->GetErrorRMS();
 884        Fdta[4*i+2] = SrzL[i]->GetRMS();
 885        Fdta[4*i+3] = SrzL[i]->GetErrorRMS();
 886    } // end for i
 887    for(i=0;i<10;i++) Ndta[i] = Nqualless[i];
 888    Ndta[10] = Sad->GetN();
 889    Ndta[11] = NqualgtMAX;
 890
 891    printf("FitAllTracks: %d tracks cut leaving %d, with a chi squared <%f\n",
 892           NqualgtMAX,Sad->GetN(),trkqualMAX);
 893    printf("The number of tracks with chi squared <1");
 894    for(i=1;i<Nqualmax;i++) printf(",%d",i+1);
 895    printf("\nFitAllTracks: %d",Nqualless[0]);
 896    for(i=1;i<Nqualmax;i++) printf(",%d",Nqualless[i]);
 897    printf("\n");
 898    printf("FitAllTracks: RMSs of ad rxz = %e+-%e %e+-%e\n",
 899           Sad->GetRMS(),Sad->GetErrorRMS(),
 900           SrxzL->GetRMS(),SrxzL->GetErrorRMS());
 901
 902    printf("FitAllTracks: Residuals by layer x=-rphi z ");
 903    for(i=0;i<6;i++) printf(":%d:%e+-%e %e+-%e ",i+1,
 904                            SrxL[i]->GetRMS(),SrxL[i]->GetErrorRMS(),
 905                            SrzL[i]->GetRMS(),SrzL[i]->GetErrorRMS());
 906    printf("\n");
 907
 908 // Setup and fill projections
 909
 910    for(i=0;i<10;i++){
 911        xm = 0.5*Double_t(i);
 912        ym = xm+0.5;
 913        j  = Ttqp->GetYaxis()->FindBin(xm);
 914        k  = Ttqp->GetYaxis()->FindBin(ym);
 915        sprintf(filename,"Track Quality for %3.1f<p<%3.1f",xm,ym);
 916        Tptqp[i] = Ttqp->ProjectionX(filename,j,k,"E");
 917        Tptqp[i]->SetXTitle("Track Quality (Chi squared/df)");
 918    } // end for i
 919
 920    TH1D *Ttq  = Ttqp->ProjectionX("Ttq",0,Ttqp->GetNbinsY()+1,"E");
 921    Ttq->SetXTitle("Chi^2/degree freedom");
 922    TH1D *Tp   = Ttqp->ProjectionY("Tp", 0,Ttqp->GetNbinsX()+1,"E");
 923    Tp->SetXTitle("Momentum GeV/c");
 924    TH1D *Tdt  = Tdttq->ProjectionX("Tdt",0,Tdttq->GetNbinsY()+1,"E");
 925    Tdt->SetXTitle("Distanct to true vertex (cm)");
 926    TH1D *Tad  = Tadtq->ProjectionX("Tad",0,Tadtq->GetNbinsY()+1,"E");
 927    Tad->SetXTitle("theta fitted line - atan(pt/p) (rad)");
 928    TH1D *TrxL[7],*TrzL[7];
 929 //   TH1D *TrrphiG[7],*TrzG[7];
 930    for(i=0;i<7;i++){
 931        sprintf(hid,"TrxL%1.1d",i+1);
 932        TrxL[i] = TrxrzL[i]->ProjectionX(hid,0,TrxrzL[i]->GetNbinsY()+1,"E");
 933        TrxL[i]->SetXTitle("xi-x(i) (cm)");
 934        sprintf(hid,"TrzL%1.1d",i+1);
 935        TrzL[i] = TrxrzL[i]->ProjectionY(hid,0,TrxrzL[i]->GetNbinsX()+1,"E");
 936        TrzL[i]->SetXTitle("zi-z(i) (cm)");
 937 //
 938 //       sprintf(hid,"TrrphiG%1.1d",i+1);
 939 //       TrrphiG[i] = TrrprzG[i]->ProjectionX(hid,0,
 940 //                                            TrrprzG[i]->GetNbinsY()+1,"E");
 941 //       TrrphiG[i]->SetXTitle("rPhii-rPhi(i) (cm)");
 942 //       sprintf(hid,"TrzG%1.1d",i+1);
 943 //       TrzG[i] = TrrprzG[i]->ProjectionY(hid,0,
 944 //                                         TrrprzG[i]->GetNbinsX()+1,"E");
 945 //       TrzG[i]->SetXTitle("zGi-zG(i) (cm)");
 946    } // end for i
 947
 948
 949    Hfile->Write();
 950
 951    if(printit){
 952        TCanvas *c0 = new TCanvas("c0","Track quality distribution",
 953                                  500,100,600,700);
 954        Ttqall->Draw();
 955        sprintf(filename,"%s_T_tqall.ps",sfile);
 956        if(printit) c0->Print(filename);
 957        Ttqp->Draw("COL");
 958        sprintf(filename,"%s_T_tq_p.ps",sfile);
 959        if(printit) c0->Print(filename);
 960        Ttq->Draw();
 961        sprintf(filename,"%s_T_tq.ps",sfile);
 962        if(printit) c0->Print(filename);
 963        Tp->Draw();
 964        sprintf(filename,"%s_T_p.ps",sfile);
 965        if(printit) c0->Print(filename);
 966        Tdttq->Draw("COL");
 967        sprintf(filename,"%s_T_dt_tq.ps",sfile);
 968        if(printit) c0->Print(filename);
 969        Tdt->Draw();
 970        sprintf(filename,"%s_T_dt.ps",sfile);
 971        if(printit) c0->Print(filename);
 972        Tadtq->Draw("COL");
 973        sprintf(filename,"%s_T_ad_tq.ps",sfile);
 974        if(printit) c0->Print(filename);
 975        Tad->Draw();
 976        sprintf(filename,"%s_T_ad.ps",sfile);
 977        if(printit) c0->Print(filename);
 978        Tadp->Draw("COL");
 979        sprintf(filename,"%s_T_ad_p.ps",sfile);
 980        if(printit) c0->Print(filename);
 981        Taddt->Draw("COL");
 982        sprintf(filename,"%s_T_ad_dt.ps",sfile);
 983        if(printit) c0->Print(filename);
 984        for(i=0;i<7;i++){
 985            TrxrzL[i]->Draw("COL");
 986            sprintf(filename,"%s_T%1.1d_rx_rz.ps",sfile,i+1);
 987            if(printit) c0->Print(filename);
 988            TrxL[i]->Draw();
 989            sprintf(filename,"%s_T%1.1d_rx.ps",sfile,i+1);
 990            if(printit) c0->Print(filename);
 991            TrzL[i]->Draw();
 992            sprintf(filename,"%s_T%1.1d_rz.ps",sfile,i+1);
 993            if(printit) c0->Print(filename);
 994        } // end for i
 995        for(i=0;i<10;i++){
 996            Tptqp[i]->Draw();
 997            sprintf(filename,"%s_T_tq_p%1.1d.ps",sfile,i);
 998            if(printit) c0->Print(filename);
 999        } // end for i
1000    } // end if printit
1001
1002 // Delet allocated stuff.
1003    for(i=0;i<7;i++) {
1004        delete TrxL[i];
1005        delete TrzL[i];
1006 //       delete TrrphiG[i];
1007 //       delete TrzG[i];
1008    } // end for i
1009    for(i=0;i<10;i++) delete Tptqp[i];
1010    delete Ttqp;
1011    delete Tdttq;
1012    delete Tadtq;
1013    delete Taddt;
1014    delete Tadp;
1015    for(i=0;i<7;i++) delete TrxrzL[i];
1016 //   for(i=0;i<7;i++) delete TrrprzG[i];
1017    delete Ttqall;
1018    delete Ttq;
1019    delete Tp;
1020    delete Tdt;
1021    delete Tad;
1022 //   printf("finished with track fitting\n");
1023    delete Sad;
1024    delete SrxzL;
1025    for(i=0;i<6;i++) delete SrxL[i];
1026    for(i=0;i<6;i++) delete SrzL[i];
1027    return;
1028 }
1029 //______________________________________________________________________
1030 void FitAllTracksG(ClustAl_tl *trk,Int_t ntrk,Float_t *v0,AliITSgeom *gm,
1031                   const char *sfile,TFile *Hfile){
1032    // Local Variables
1033    Int_t    i,j,k,lay,lad,det;
1034    Double_t xm,ym,zm,dt,ad,rh,rp,rphih,rphip;
1035    Double_t tp,tx,rx,rz/*,rr*/,rrphi;
1036    Double_t xg[3],xl[3],xp,yp,xh,yh,a,b,c,d,x1[3],x2[3],qual;
1037    Double_t trkqualMAX = 50.0;
1038    Bool_t   printit = kFALSE;
1039    char     filename[80],hid[10];
1040    Int_t    Nqualmax = 10,NqualgtMAX=0;
1041    Int_t    Nqualless[10];
1042 // Distrobution statisitics
1043    AliITSstatistics *Sad    = new AliITSstatistics(4);
1044    AliITSstatistics *SrxzL  = new AliITSstatistics(4);
1045    AliITSstatistics *SrxL[6];
1046    AliITSstatistics *SrzL[6];
1047    for(i=0;i<6;i++) {
1048        SrxL[i] = new AliITSstatistics(4);
1049        SrzL[i] = new AliITSstatistics(4);
1050    } // end for i
1051
1052    for(i=0;i<Nqualmax;i++) Nqualless[i] = 0;
1053
1054 // Setup histograms
1055
1056    TH2F *Ttqp  = new TH2F("Ttqp","Track quality vs. momentum",
1057                          500,0.0,50.0,500,0.0,10.0);
1058    Ttqp->Sumw2();
1059    Ttqp->SetXTitle("Chi^2/degree freedom");
1060    Ttqp->SetYTitle("Momentum GeV/c");
1061    TH2F *Tdttq = new TH2F("Tdttq",
1062                           "Distance from true vertex vs. track quality",
1063                           500,0.0,0.10,100,0.0,50.0);
1064    Tdttq->Sumw2();
1065    Tdttq->SetXTitle("Distanct to true vertex (cm)");
1066    Tdttq->SetYTitle("Chi^2/degree freedom");
1067    TH2F *Tadtq = new TH2F("Tadtq",
1068                      "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. track quality",
1069                           500,-0.03,0.03,500,0.0,50.0);
1070    Tadtq->Sumw2();
1071    Tadtq->SetXTitle("theta fitted line - atan(pt/p) (rad)");
1072    Tadtq->SetYTitle("Chi^2/degree freedom");
1073    TH2F *Taddt = new TH2F("Taddt",
1074          "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. Distance from true vertex",
1075                           500,-0.03,0.03,200,0.0,0.10);
1076    Taddt->Sumw2();
1077    Taddt->SetXTitle("theta fitted line - atan(pt/p) (rad)");
1078    Taddt->SetYTitle("distance to true vertex (cm)");
1079    TH2F *Tadp = new TH2F("Tadp",
1080          "atan(sqrt(b*b+d*d))-atan(pt/|pz|) vs. momentum",
1081                           500,-0.03,0.03,200,0.0,10.0);
1082    Tadp->Sumw2();
1083    Tadp->SetXTitle("theta fitted line - atan(pt/p) (rad)");
1084    Tadp->SetYTitle("Momentum (GeV/c)");
1085
1086    TH2F *TrxrzL[7],*TrrprzG[7];
1087    for(i=0;i<7;i++){
1088        sprintf(filename,"Layer %1.1d: xi-x(i) vs. zi-z(i) local",i+1);
1089        if(i==6) sprintf(filename,"Sum of all layers: xi-x(i) "
1090                         "vs. zi-z(i) local");
1091        sprintf(hid,"TrxrzL%1.1d",i+1);
1092        TrxrzL[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
1093        TrxrzL[i]->Sumw2();
1094        TrxrzL[i]->SetXTitle("xi-x(i) (cm) local");
1095        TrxrzL[i]->SetYTitle("zi-z(i) (cm) local");
1096 //
1097        sprintf(filename,"Layer %1.1d: rPhii-rPhi(i) vs. zi-z(i) local",i+1);
1098        if(i==6) sprintf(filename,"Sum of all layers: rPhii-rPhi(i) "
1099                         "vs. zi-z(i) local");
1100        sprintf(hid,"TrrprzG%1.1d",i+1);
1101        TrrprzG[i] = new TH2F(hid,filename,500,-1.0,1.0,500,-1.0,1.0);
1102        TrrprzG[i]->Sumw2();
1103        TrrprzG[i]->SetXTitle("rPhii-rPhi(i) (cm) local");
1104        TrrprzG[i]->SetYTitle("zi-z(i) (cm) local");
1105    } // end for i
1106
1107    TH1D *Tptqp[10];
1108    TH1D *Ttqall = new TH1D("Ttqall","Track quality for all tracks fit",
1109                            500,0.0,1000.0);
1110    Ttqall->Sumw2();
1111    Ttqall->SetXTitle("Track quality: Chi squared per degree freedom");
1112
1113 // Fit each track and fill histograms and the like.
1114
1115    for(i=0;i<ntrk;i++){
1116        if(FitTrackToLine(trk[i])!=0) continue;  // fit track to a line
1117       qual = trk[i].qual;
1118       if(qual<0.0) continue;
1119       // Fill histograms and statistics before cource chi squared cut
1120       Ttqall->Fill(qual,1.0);
1121       for(j=0;j<Nqualmax;j++) if(qual<(Double_t) j+1) Nqualless[j]++;
1122       if(qual>trkqualMAX) {
1123           NqualgtMAX++;
1124           continue;
1125       } // end if iqual>trkqualMAX=50.0
1126       a = trk[i].a;
1127       b = trk[i].b;
1128       c = trk[i].c;
1129       d = trk[i].d;
1130       zm  = b*(v0[0]-a) + d*(v0[1]-c);
1131       dt  = (b*b + d*d);
1132       if(dt!=0.0) zm /= dt;
1133       else{
1134           printf("FitAllTracks: trk[%d].b=%f trk[%d].d=%f\n",i,b,i,d);
1135           zm=0.0;
1136       }// end if else dt!=0.0
1137       xm  = a + b*zm;
1138       ym  = c + d*zm;
1139       dt   = sqrt((xm-v0[0])*(xm-v0[0]) + (ym-v0[1])*(ym-v0[1]) +
1140                  (zm-v0[2])*(zm-v0[2]));
1141       tx  = atan(sqrt(b*b + d*d));
1142       tp  = atan2(trk[i].pt,fabs(trk[i].pz));
1143       ad  = tx-tp;
1144       // Fill histograms and statistics
1145       Tadtq->Fill(ad,qual,1.0);
1146       Sad->AddValue(ad,1.0);
1147       Taddt->Fill(ad,dt,1.0);
1148       Tadp->Fill(ad,trk[i].p,1.0);
1149       Tdttq->Fill(dt,qual,1.0);
1150       Ttqp->Fill(qual,trk[i].p,1.0);
1151       for(j=0;j<trk[i].nclust;j++){
1152           lay = trk[i].clust[j].lay;
1153           lad = trk[i].clust[j].lad;
1154           det = trk[i].clust[j].det;
1155           xh = trk[i].clust[j].xg;
1156           yh = trk[i].clust[j].yg;
1157           rh    = hypot(xh,yh);
1158           rphih = atan2(yh,xh);
1159           xp = a+b*trk[i].clust[j].zg;
1160           yp = c+d*trk[i].clust[j].zg;
1161           rp    = hypot(xp,yp);
1162           rphip = atan2(yp,xp);
1163           //
1164           if(fabs(rphih-rphip)>TMath::Pi()){
1165               if(rphih<rphip) rphih += 2.0*TMath::Pi();
1166               else rphip = 2.0*TMath::Pi();
1167           } //
1168           rrphi = rh*rphih-rp*rphip;
1169           //
1170           x1[0] = a; x1[1] = b; x1[2] = 0.0;
1171           gm->GtoL(lay,lad,det,x1,xl);
1172           x1[0] = xl[0]; x1[1] = xl[1]; x1[2] = xl[2];
1173           if(trk[i].clust[j].zg!=0.0){
1174               x2[0] = xp; x2[1] = yp; x2[2] = trk[i].clust[j].zg;
1175           }else{
1176               x2[0] = a+b; x2[1] = c+d; x2[2] = 1.0;
1177           } // end if trk[i].clst[j].zg!=0
1178           gm->GtoL(lay,lad,det,x2,xl);
1179           x2[0] = xl[0]; x2[1] = xl[1]; x2[2] = xl[2];
1180           xg[0] = trk[i].clust[j].xg;
1181           xg[1] = trk[i].clust[j].yg;
1182           xg[2] = trk[i].clust[j].zg;
1183           xl[0] = trk[i].clust[j].xl;
1184           xl[1] = trk[i].clust[j].yl;
1185           xl[2] = trk[i].clust[j].zl;
1186           if(x1[1]!=x2[1]){
1187               rx = x1[0] - x1[1]*(x1[0]-x2[0])/(x1[1]-x2[1]);
1188               rz = x1[2] - x1[1]*(x1[2]-x2[2])/(x1[1]-x2[1]);
1189           }else{
1190               continue;
1191           } // end if x1[1]!=x2[1]
1192           rx = rx - xl[0];
1193           rz = rz - xl[2];
1194           // Fill histograms and statistics
1195           TrxrzL[lay-1]->Fill(rx,rz,1.0);
1196           TrxrzL[6]->Fill(rx,rz,1.0);
1197           TrrprzG[lay-1]->Fill(rrphi,rz,1.0);
1198           TrrprzG[6]->Fill(rrphi,rz,1.0);
1199           SrxL[lay-1]->AddValue(rx,1.0);
1200           SrzL[lay-1]->AddValue(rz,1.0);
1201           SrxzL->AddValue(rx,1.0);
1202           SrxzL->AddValue(rz,1.0);
1203       } // end for j
1204    } // end for i
1205
1206 // Write out information
1207
1208    printf("FitAllTracks: %d tracks cut leaving %d, with a chi squared >%f\n",
1209           NqualgtMAX,Sad->GetN(),trkqualMAX);
1210    printf("The number of tracks with chi squared <1");
1211    for(i=1;i<Nqualmax;i++) printf(",%d",i+1);
1212    printf("\nFitAllTracks: %d",Nqualless[0]);
1213    for(i=1;i<Nqualmax;i++) printf(",%d",Nqualless[i]);
1214    printf("\n");
1215    printf("FitAllTracks: RMSs of ad, rxz = %f+-%f %f+-%f\n",
1216           Sad->GetRMS(),Sad->GetErrorRMS(),
1217           SrxzL->GetRMS(),SrxzL->GetErrorRMS());
1218
1219    printf("FitAllTracks: Residuals by layer x=-rphi z ");
1220    for(i=0;i<6;i++) printf("%f+-%f %f+-%f ",
1221                            SrxL[i]->GetRMS(),SrxL[i]->GetErrorRMS(),
1222                            SrzL[i]->GetRMS(),SrzL[i]->GetErrorRMS());
1223    printf("\n");
1224
1225 // Setup and fill projections
1226
1227    for(i=0;i<10;i++){
1228        xm = 0.5*Double_t(i);
1229        ym = xm+0.5;
1230        j  = Ttqp->GetYaxis()->FindBin(xm);
1231        k  = Ttqp->GetYaxis()->FindBin(ym);
1232        sprintf(filename,"Track Quality for %3.1f<p<%3.1f",xm,ym);
1233        Tptqp[i] = Ttqp->ProjectionX(filename,j,k,"E");
1234        Tptqp[i]->SetXTitle("Track Quality (Chi squared/df)");
1235    } // end for i
1236
1237    TH1D *Ttq  = Ttqp->ProjectionX("Ttq",0,Ttqp->GetNbinsY()+1,"E");
1238    Ttq->SetXTitle("Chi^2/degree freedom");
1239    TH1D *Tp   = Ttqp->ProjectionY("Tp", 0,Ttqp->GetNbinsX()+1,"E");
1240    Tp->SetXTitle("Momentum GeV/c");
1241    TH1D *Tdt  = Tdttq->ProjectionX("Tdt",0,Tdttq->GetNbinsY()+1,"E");
1242    Tdt->SetXTitle("Distanct to true vertex (cm)");
1243    TH1D *Tad  = Tadtq->ProjectionX("Tad",0,Tadtq->GetNbinsY()+1,"E");
1244    Tad->SetXTitle("theta fitted line - atan(pt/p) (rad)");
1245    TH1D *TrxL[7],*TrzL[7],*TrrphiG[7],*TrzG[7];
1246    for(i=0;i<7;i++){
1247        sprintf(hid,"TrxL%1.1d",i+1);
1248        TrxL[i] = TrxrzL[i]->ProjectionX(hid,0,TrxrzL[i]->GetNbinsY()+1,"E");
1249        TrxL[i]->SetXTitle("xi-x(i) (cm)");
1250        sprintf(hid,"TrzL%1.1d",i+1);
1251        TrzL[i] = TrxrzL[i]->ProjectionY(hid,0,TrxrzL[i]->GetNbinsX()+1,"E");
1252        TrzL[i]->SetXTitle("zi-z(i) (cm)");
1253 //
1254        sprintf(hid,"TrrphiG%1.1d",i+1);
1255        TrrphiG[i] = TrrprzG[i]->ProjectionX(hid,0,
1256                                               TrrprzG[i]->GetNbinsY()+1,"E");
1257        TrrphiG[i]->SetXTitle("rPhii-rPhi(i) (cm)");
1258        sprintf(hid,"TrzG%1.1d",i+1);
1259        TrzG[i] = TrrprzG[i]->ProjectionY(hid,0,
1260                                            TrrprzG[i]->GetNbinsX()+1,"E");
1261        TrzG[i]->SetXTitle("zGi-zG(i) (cm)");
1262    } // end for i
1263
1264
1265    Hfile->Write();
1266
1267    if(printit){
1268        TCanvas *c0 = new TCanvas("c0","Track quality distribution",
1269                                  500,100,600,700);
1270        Ttqall->Draw();
1271        sprintf(filename,"%s_T_tqall.ps",sfile);
1272        if(printit) c0->Print(filename);
1273        Ttqp->Draw("COL");
1274        sprintf(filename,"%s_T_tq_p.ps",sfile);
1275        if(printit) c0->Print(filename);
1276        Ttq->Draw();
1277        sprintf(filename,"%s_T_tq.ps",sfile);
1278        if(printit) c0->Print(filename);
1279        Tp->Draw();
1280        sprintf(filename,"%s_T_p.ps",sfile);
1281        if(printit) c0->Print(filename);
1282        Tdttq->Draw("COL");
1283        sprintf(filename,"%s_T_dt_tq.ps",sfile);
1284        if(printit) c0->Print(filename);
1285        Tdt->Draw();
1286        sprintf(filename,"%s_T_dt.ps",sfile);
1287        if(printit) c0->Print(filename);
1288        Tadtq->Draw("COL");
1289        sprintf(filename,"%s_T_ad_tq.ps",sfile);
1290        if(printit) c0->Print(filename);
1291        Tad->Draw();
1292        sprintf(filename,"%s_T_ad.ps",sfile);
1293        if(printit) c0->Print(filename);
1294        Tadp->Draw("COL");
1295        sprintf(filename,"%s_T_ad_p.ps",sfile);
1296        if(printit) c0->Print(filename);
1297        Taddt->Draw("COL");
1298        sprintf(filename,"%s_T_ad_dt.ps",sfile);
1299        if(printit) c0->Print(filename);
1300        for(i=0;i<7;i++){
1301            TrxrzL[i]->Draw("COL");
1302            sprintf(filename,"%s_T%1.1d_rx_rz.ps",sfile,i+1);
1303            if(printit) c0->Print(filename);
1304            TrxL[i]->Draw();
1305            sprintf(filename,"%s_T%1.1d_rx.ps",sfile,i+1);
1306            if(printit) c0->Print(filename);
1307            TrzL[i]->Draw();
1308            sprintf(filename,"%s_T%1.1d_rz.ps",sfile,i+1);
1309            if(printit) c0->Print(filename);
1310        } // end for i
1311        for(i=0;i<10;i++){
1312            Tptqp[i]->Draw();
1313            sprintf(filename,"%s_T_tq_p%1.1d.ps",sfile,i);
1314            if(printit) c0->Print(filename);
1315        } // end for i
1316    } // end if printit
1317
1318 // Delet allocated stuff.
1319    for(i=0;i<7;i++) {
1320        delete TrxL[i];
1321        delete TrzL[i];
1322        delete TrrphiG[i];
1323        delete TrzG[i];
1324    } // end for i
1325    for(i=0;i<10;i++) delete Tptqp[i];
1326    delete Ttqp;
1327    delete Tdttq;
1328    delete Tadtq;
1329    delete Taddt;
1330    delete Tadp;
1331    for(i=0;i<7;i++) delete TrxrzL[i];
1332    for(i=0;i<7;i++) delete TrrprzG[i];
1333    delete Ttqall;
1334    delete Ttq;
1335    delete Tp;
1336    delete Tdt;
1337    delete Tad;
1338 //   printf("finished with track fitting\n");
1339    delete Sad;
1340    delete SrxzL;
1341    for(i=0;i<6;i++) delete SrxL[i];
1342    for(i=0;i<6;i++) delete SrzL[i];
1343    return;
1344 }
1345 //______________________________________________________________________
1346 void FindVertex2(ClustAl_tl &trk1,ClustAl_tl &trk2,Double_t *vt,Double_t &d){
1347    // Local Variables
1348    Double_t x1,y1,x2,y2;
1349    Double_t z1,z2;
1350    Double_t a1,b1,c1,d1;
1351    Double_t a2,b2,c2,d2;
1352    Double_t da,db,dc,dd,dbd;
1353    Double_t den,num1,num2;
1354
1355    a1 = trk1.a; b1 = trk1.b; c1 = trk1.c; d1 = trk1.c;
1356    a2 = trk2.a; b2 = trk2.b; c2 = trk2.c; d2 = trk2.c;
1357
1358    // Find z1 and z2 of points of closest approch.
1359    da = a1-a2;
1360    db = b1-b2;
1361    dc = c1-c2;
1362    dd = d1-d2;
1363    dbd = b1*d2-d1*b2;
1364    den = -db*db - dd*dd - dbd*dbd;
1365    num1 = da*(db+d1*dbd) + dc*(dd-b1*dbd);
1366    num2 = da*(db+d2*dbd) + dc*(dd-b2*dbd);
1367    if(den!=0.0){
1368        z1 = num1/den;
1369        z2 = num2/den;
1370    }else{ // parallel lines
1371        z1 = 0.0;
1372        z2 = 0.0;
1373    } // end if den!-0.0
1374
1375    // find coordinate of closest approch and distance between.
1376    x1    = a1+b1*z1; y1 = c1+d1*z1;
1377    x2    = a2+b2*z2; y2 = c2*d2*z2;
1378    vt[0] = 0.5*(x1+x2);
1379    vt[1] = 0.5*(y1+y2);
1380    vt[2] = 0.5*(z1+z2);
1381    d     = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) );
1382
1383    return;
1384 }
1385 //______________________________________________________________________
1386 void FitVertexAll(ClustAl_tl *trk,Int_t ntrk,const char *sfile,TFile *Hfile){
1387    // Local Variables
1388    Int_t    i,j,k;
1389    Double_t d,vz,r,q;
1390    Double_t vt[3];
1391    Bool_t   printit = kFALSE;
1392    char     filename[80],hid[10];
1393    Double_t trkqualMAX = 50.0;
1394    const Int_t    IntQualMax=10;
1395    Int_t    QualInt[IntQualMax];
1396    AliITSstatistics *Svz   = new AliITSstatistics(4);
1397    AliITSstatistics *Svx   = new AliITSstatistics(4);
1398    AliITSstatistics *Svy   = new AliITSstatistics(4);
1399    //   Double_t vzmi=-0.20,vzma=0.20,vxmi=-1.0,vxma=1.0,vymi=-1.0,vyma=1.0;
1400
1401    for(i=0;i<IntQualMax;i++) QualInt[i] = 0;
1402
1403    TH1F *Vvzbtq[IntQualMax];
1404    for(i=0;i<IntQualMax;i++){
1405        sprintf(hid,"Vvzbta%2.2d",i+1);
1406        sprintf(filename,"Z of vertex of pairs for tracks with chi squared <%d",
1407                i+1);
1408        Vvzbtq[i] = new TH1F(hid,filename,100,-1.0,1.0);
1409        Vvzbtq[i]->Sumw2();
1410        Vvzbtq[i]->SetXTitle("Z of vertex (cm)");
1411    } // end for i
1412
1413    TH2F *Vvztq  = new TH2F("Vvztq","Z of vertex of pairs vs. track quality cut",
1414                            500,-0.20,0.20,200,0.0,20.0);
1415    Vvztq->Sumw2();
1416    Vvztq->SetXTitle("Z of vertex (cm)");
1417    Vvztq->SetYTitle("Chi^2/degree freedom");
1418    TH2F *Vvzvr  = new TH2F("Vvzvr","Z vs. R of vertex of pairs",
1419                                                200,-0.20,0.20,400,0.0,5.0);
1420    Vvzvr->Sumw2();
1421    Vvzvr->SetXTitle("Z of vertex (cm)");
1422    Vvzvr->SetYTitle("R of vertex (cm)");
1423    TH2F *Vdtq   = new TH2F("Vdtq","Distance between lines vs. track quality cut",
1424                            500,0.0,0.2,200,0.0,50.0);
1425    Vdtq->Sumw2();
1426    Vdtq->SetXTitle("minimum distance between lines (cm)");
1427    Vdtq->SetYTitle("Chi^2/degree freedom");
1428    TH2F *Vvzd   = new TH2F("Vvzd","Z vertex vs. Dist between lines",
1429                                                500,-5.0,5.0,200,0.0,5.0);
1430    Vvzd->Sumw2();
1431    Vvzd->SetXTitle("Z of vertex (cm)");
1432    Vvzd->SetYTitle("minimmum distance between lined (cm)");
1433    TH2F *Vvxvy  = new TH2F("Vvxvy","X vertex vs. Y vertex",
1434                                                200,-10.0,10.0,200,-10.0,10.0);
1435    Vvxvy->Sumw2();
1436    Vvxvy->SetXTitle("X of vertex (cm)");
1437    Vvxvy->SetYTitle("Y of vertex (cm)");
1438    TH2F *Vvxvz  = new TH2F("Vvxvz","X vertex vs. Z vertex",
1439                                                200,-10.0,10.0,200,-10.0,10.0);
1440    Vvxvz->Sumw2();
1441    Vvxvz->SetXTitle("X of vertex (cm)");
1442    Vvxvz->SetYTitle("Z of vertex (cm)");
1443    TH2F *Vvyvz  = new TH2F("Vvyvz","Y vertex vs. Z vertex",
1444                                                200,-10.0,10.0,200,-10.0,10.0);
1445    Vvyvz->Sumw2();
1446    Vvyvz->SetXTitle("Y of vertex (cm)");
1447    Vvyvz->SetYTitle("Z of vertex (cm)");
1448
1449    for(i=0;i<ntrk-1;i++){
1450        if(trk[i].qual>trkqualMAX||trk[i].qual<0.0) continue;
1451        for(j=i+1;j<ntrk;j++){
1452            if(trk[j].qual>trkqualMAX||trk[i].qual<0.0) continue;
1453            FindVertex2(trk[i],trk[j],vt,d);
1454            q  = TMath::Max(trk[i].qual,trk[j].qual);
1455            r  = vt[0]*vt[0] + vt[1]*vt[1];
1456            r  = sqrt(r);
1457            vz = vt[2];
1458            Vvxvy->Fill(vt[0],vt[1],1.0);
1459            Vvxvz->Fill(vt[0],vt[2],1.0);
1460            Vvyvz->Fill(vt[1],vt[2],1.0);
1461            Svx->AddValue(vt[0],1.0);
1462            Svy->AddValue(vt[1],1.0);
1463            Svz->AddValue(vt[2],1.0);
1464            Vvztq->Fill(vz,q,1.0);
1465            Vvzvr->Fill(vz,r,1.0);
1466            Vdtq->Fill(d,q,1.0);
1467            Vvzd->Fill(vz,d,1.0);
1468            for(k=0;k<IntQualMax;k++)if(q<(Double_t)k+1.)QualInt[k]++;
1469            for(k=0;k<IntQualMax;k++)if(q<(Double_t)k+1.)Vvzbtq[k]->Fill(vz,1.);
1470        } // end for j
1471    } // end for i
1472
1473    TH1D *Vd  = Vvzd->ProjectionY ("Vd", 0,Vvzd->GetNbinsX()+1,"E");
1474    Vd->SetXTitle("minimum distance between lines (cm)");
1475    TH1D *Vvr = Vvzvr->ProjectionY("Vvr",0,Vvzvr->GetNbinsX()+1,"E");
1476    Vvr->SetXTitle("R of vertex (cm)");
1477    TH1D *Vtq = Vdtq->ProjectionY ("Vtq",0,Vdtq->GetNbinsX()+1,"E");
1478    Vtq->SetXTitle("Chi^2/degree freedom");
1479 //
1480    TH1D *Vvx = Vvxvy->ProjectionX("Vvx",0,Vvxvy->GetNbinsY()+1,"E");
1481    Vvx->SetXTitle("X of vertex (cm)");
1482    TH1D *Vvy = Vvxvy->ProjectionY("Vvy",0,Vvxvy->GetNbinsX()+1,"E");
1483    Vvy->SetXTitle("Y of vertex (cm)");
1484    TH1D *Vvz = Vvxvz->ProjectionX ("Vvz",0,Vvxvz->GetNbinsY()+1,"E");
1485    Vvz->SetXTitle("Z of vertex (cm)");
1486
1487    printf("FitVertexAll: N(track pairs Qual <=1");
1488    for(k=1;k<IntQualMax;k++) printf(",%d",k+1);
1489    printf("\nFitVertexAll: %d",QualInt[0]);
1490    for(k=1;k<IntQualMax;k++) printf(",%d",QualInt[k]);
1491    printf("FitVertexAll: RMS of vx vy vz=%e+-%e %e+-%e %e+-%e\n",
1492           Svx->GetRMS(),Svx->GetErrorRMS(),
1493           Svy->GetRMS(),Svy->GetErrorRMS(),
1494           Svz->GetRMS(),Svy->GetErrorRMS());
1495    Hfile->Write();
1496
1497    if(printit){
1498        TCanvas *c1 = new TCanvas("c1","Vertex info",500,100,600,700);
1499        for(k=0;k<IntQualMax;k++){
1500            Vvzbtq[k]->Draw();
1501            sprintf(filename,"%s_V_vz_tq%2.2d.ps",sfile,k+1);
1502            if(printit) c1->Print(filename);
1503        } // end for k
1504        Vvzd->Draw("COL");
1505        sprintf(filename,"%s_V_vz_d.ps",sfile);
1506        if(printit) c1->Print(filename);
1507        Vvzvr->Draw("COL");
1508        sprintf(filename,"%s_V_vz_vr.ps",sfile);
1509        if(printit) c1->Print(filename);
1510        Vdtq->Draw("COL");
1511        sprintf(filename,"%s_V_d_tq.ps",sfile);
1512        if(printit) c1->Print(filename);
1513        Vvztq->Draw("COL");
1514        sprintf(filename,"%s_V_vz_tq.ps",sfile);
1515        if(printit) c1->Print(filename);
1516        Vvz->Draw();
1517        sprintf(filename,"%s_V_vz.ps",sfile);
1518        if(printit) c1->Print(filename);
1519        Vd->Draw();
1520        sprintf(filename,"%s_V_d.ps",sfile);
1521        if(printit) c1->Print(filename);
1522        Vvr->Draw();
1523        sprintf(filename,"%s_V_vr.ps",sfile);
1524        if(printit) c1->Print(filename);
1525        Vtq->Draw();
1526        sprintf(filename,"%s_V_tq.ps",sfile);
1527        if(printit) c1->Print(filename);
1528        Vvxvy->Draw("COL");
1529        sprintf(filename,"%s_V_vx_vy.ps",sfile);
1530        if(printit) c1->Print(filename);
1531        Vvxvz->Draw("COL");
1532        sprintf(filename,"%s_V_vx_vz.ps",sfile);
1533        if(printit) c1->Print(filename);
1534        Vvyvz->Draw("COL");
1535        sprintf(filename,"%s_V_vy_vz.ps",sfile);
1536        if(printit) c1->Print(filename);
1537        Vvx->Draw();
1538        sprintf(filename,"%s_V_vx.ps",sfile);
1539        if(printit) c1->Print(filename);
1540        Vvy->Draw();
1541        sprintf(filename,"%s_V_vy.ps",sfile);
1542        if(printit) c1->Print(filename);
1543    } // end if printit
1544    for(k=0;k<IntQualMax;k++) delete Vvzbtq[k];
1545    delete Vvzd;
1546    delete Vvzvr;
1547    delete Vdtq;
1548    delete Vvztq;
1549    delete Vvz;
1550    delete Vd;
1551    delete Vvr;
1552    delete Vtq;
1553    delete Vvxvy;
1554    delete Vvxvz;
1555    delete Vvyvz;
1556    delete Vvx;
1557    delete Vvy;
1558    delete Svz;
1559    delete Svx;
1560    delete Svy;
1561
1562    return;
1563 }
1564 //______________________________________________________________________
1565 void OnlyOneGeometry(char *filename,AliITSgeom *gm,AliITSgeom &gm2,
1566                      Float_t trans[],Float_t rot[]){
1567
1568     ifstream f_in(filename);
1569     if(f_in){
1570         printf("reading from geometry file %s\n",filename);
1571         gm2.ReadGeom(f_in);
1572         f_in.close();
1573     }else{
1574         f_in.close();
1575         gm2 = *gm;
1576         gm2.RandomCylindericalChange(trans,rot);
1577         printf("writting to geometry file %s\n",filename);
1578         ofstream f_out(filename);
1579         gm2.PrintGeom(f_out);
1580         f_out.close();
1581     } // end if
1582     return;
1583 }
1584 //______________________________________________________________________
1585 void deleteClustAl(ClustAl_tl *trk,Int_t ntrk){
1586
1587     Int_t i;
1588
1589     for(i=0;i<ntrk;i++) delete[] trk[i].clust;
1590     return;
1591 }
1592 //______________________________________________________________________
1593 void RunAlignment(Int_t evnt,Float_t fraction=1.0){
1594 // define some variables for later use.
1595 //
1596     printf("gAlice=%p and evnt=%d\n",gAlice,evnt);
1597 //
1598     Int_t nparticles = gAlice->GetEvent(evnt);
1599     printf("nparticles %d\n",nparticles);
1600     if (nparticles <= 0) return; /* get next event */
1601
1602 // Get pointers to Alice detectors and Clusts containers
1603     AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
1604     if(!ITS) return;          /* error no ITS data exit */
1605     TTree *TH   = gAlice->TreeH();
1606     Int_t Ntrkp =  (Int_t) TH->GetEntries(),ntrk;
1607     AliITSgeom gm2,gm3,*gm = ITS->GetITSgeom();
1608     Int_t Nmods = gm->GetIndexMax();
1609     Float_t   trans[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
1610                           2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
1611                           2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // cm
1612     Float_t   rots[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
1613                          2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
1614                          2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // rad
1615    Float_t     tran[3] = {0.0,0.0,0.0},rot[3] = {0.0,0.0,0.0};
1616    TFile *Hfile;
1617    Float_t     Rdta[6];
1618    Int_t Itimes=0,i,j,badmod;
1619    Double_t Chi2b;
1620
1621 // Array (stucture) of clusts for the first and second layer
1622 // this should be replaced with either clusters or digits
1623 // when they are proporly defined.
1624     AliITSAlignmentTrack *trk = new AliITSAlignmentTrack[Ntrkp];
1625     TObjArray *mods = new TObjArray(Nmods);
1626     AliITSAlignmentModule *mod;
1627
1628     printf("Ntrkp=%d\n",Ntrkp);
1629
1630     FillAliITSAlignmentTrack(trk,ntrk,Ntrkp,TH,ITS,fraction);
1631
1632     for(Int_t Isigmas=0;Isigmas<1;Isigmas++){
1633 //
1634 //      tran[0] = sigma1;
1635 //      tran[1] = sigma2;
1636 //      tran[2] = sigma3;
1637         if(Itimes==0){ tran[0] = trans[Isigmas];
1638         }else tran[0] = 0.0;
1639         if(Itimes==1){ tran[1] = trans[Isigmas];
1640         }else tran[1] = 0.0;
1641         if(Itimes==2){ tran[2] = trans[Isigmas];
1642         }else tran[2] = 0.0;
1643         if(Itimes==3){ rot[0] = rots[Isigmas];
1644         }else rot[0] = 0.0;
1645         if(Itimes==4){ rot[1] = rots[Isigmas];
1646         }else rot[1] = 0.0;
1647         if(Itimes==5){ rot[2] = rots[Isigmas];
1648         }else rot[2] = 0.0;
1649         printf("tran= %e %e %e (cm), rot=%e %e %e (rad)\n",
1650                tran[0],tran[1],tran[2],rot[0],rot[1],rot[2]);
1651 //
1652         gm2 = *gm;
1653         gm2.RandomCylindericalChange(tran,rot);
1654         gm3 = gm2;
1655         Hfile = new TFile("Alignment_geom_0_2.root","RECREATE",
1656                           "comparison of geometry before refitting");
1657         PlotGeomChanges(gm,&gm2,Hfile,Rdta);
1658         Hfile -> Close();
1659
1660         for(i=0;i<Nmods;i++) mods->AddAt(new AliITSAlignmentModule(i,&gm3,
1661                                                                   ntrk,trk),i);
1662 //
1663         j = 0;
1664         do{
1665             for(i=0;i<ntrk;i++) trk[i].SetGlobalPosition(&gm2);
1666             //
1667             for(i=0;i<ntrk;i++){
1668                 trk[i].FitToFunction(2,&gm3);
1669             } // end for i
1670             // find detector with the worst Chi2.
1671             j++;
1672             Chi2b  = 0.0;
1673             badmod = 0;
1674             for(i=0;i<Nmods;i++){
1675                 if(((AliITSAlignmentModule *)(mods->At(i)))->GetChi2()>Chi2b){
1676                     Chi2b =((AliITSAlignmentModule *)(mods->At(i)))->GetChi2();
1677                     badmod = i;
1678                 } // end if
1679             } // end for i
1680             mod = (AliITSAlignmentModule *)(mods->At(badmod));
1681             mod->AlignModule();
1682             gm3.SetTrans(badmod,mod->GetTranslationVector());
1683             gm3.SetByAngles(badmod,mod->GetRotationAngles());
1684         }while(j<Nmods);
1685     } // end for Isigmas
1686 //
1687     Hfile = new TFile("Alignment_geom_0_3.root","RECREATE",
1688                       "comparison of geometry after refitting");
1689     PlotGeomChanges(gm,&gm3,Hfile,Rdta);
1690     Hfile->Close();
1691     printf("Event %d done\n",evnt);
1692 //
1693     delete[] trk;            // now delet memory allocated above.
1694 }
1695 //______________________________________________________________________