ITS/vertex.cxx

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <time.h>
   4 #include "TMinuit.h"
   5 #include "TRandom.h"
   6 #include "TTree.h"
   7 #include "TClonesArray.h"
   8 #include "TProfile.h"
   9 #include "TH2.h"
  10 #include "TArray.h"
  11 #include "TCanvas.h"
  12 #include "TString.h"
  13 #include "AliRun.h"
  14 #include "AliITS.h"
  15 #include "AliITSgeom.h"
  16 #include "TParticle.h"
  17 #include "vertex.h"
  18
  19 // Global in this file only. Needed for Minuit
  20 Int_t gstrt,gend;  // starting and ending bin values
  21 TH1D  *gH1;        // histogram to be fit.
  22
  23 void FindVertexs(Int_t evnt,Float_t frac=1.0,Float_t len=16.82){
  24
  25       // Get pointers to Alice detectors and Hits containers
  26       AliDetector  *ITS       = gAlice->GetDetector("ITS");
  27       if(!ITS) return;          // error no ITS data exit
  28 //      TClonesArray *Particles = gAlice->Particles();
  29       TClonesArray *ITShits   = ITS->Hits();
  30       TTree        *TH        = gAlice->TreeH();
  31       Int_t        ntracks    =  (Int_t)TH->GetEntries();
  32 //      ntracks = 10;
  33
  34       printf(" %d primary tracks in the file\n",ntracks);
  35       // Array (stucture) of hits for the first and second layer
  36       // this should be replaced with either clusters or digits
  37       // when they are proporly defined.
  38       Hit_tl *spd1 = new Hit_tl[ntracks];
  39       Hit_tl *spd2 = new Hit_tl[ntracks];
  40       Hit_tl **spdi = new Hit_tl*[ntracks];
  41       Hit_tl **spdo = new Hit_tl*[ntracks];
  42       for(Int_t i=0;i<ntracks;i++){
  43           spdi[i] = &(spd1[i]);
  44 //        printf("spdi[%d]=%p spd1[%d]=%p  ",i,spdi[i],i,&(spd1[i]));
  45           spdo[i] = &(spd2[i]);
  46 //        printf("spdo[%d]=%p spd2[%d]=%p\n",i,spdo[i],i,&(spd2[i]));
  47       } // end for i
  48       Int_t   i1max,i2max;
  49       Float_t vz;
  50
  51       HitsToV(spdi,i1max,spdo,i2max,ntracks,TH,ITShits,frac,len);
  52       printf("back in Macro i1max=%d i2max=%d\n",i1max,i2max);
  53
  54 //      Float_t r1=0.0,r2=0.0;
  55 //      for(i=0;i<TMath::Max(i1max,i2max);i++){
  56 //        if(i<i1max) r1 += spdi[i]->r;
  57 //        if(i<i2max) r2 += spdo[i]->r;
  58 //      } // end for i
  59 //      printf("<r1>=%f i1max=%d\t<r2>=%f i2max=%d\n",
  60 //           r1/Float_t(i1max),i1max,r2/Float_t(i2max),i2max);
  61
  62       vz = vertexSlow(spdi,i1max,spdo,i2max);
  63
  64       printf("Slow Sorted event=%d Zvertex=%f\n",evnt,vz);
  65
  66 //      vz = vertex(spdi,i1max,spdo,i2max);
  67
  68 //      printf("Phi sorted event=%d Zvertex=%f\n",evnt,vz);
  69
  70 //      vz = vertexEta(spdi,i1max,spdo,i2max);
  71
  72 //      printf("Eta sorted event=%d Zvertex=%f\n",evnt,vz);
  73
  74       return;
  75 }
  76
  77
  78 Bool_t L_SortPhi(const Hit_tl *s1,const Hit_tl *s2){
  79   // Phi sorting function for qsort.
  80    Float_t a;
  81
  82    a = s1->phir - s2->phir;
  83    if(a<0.0) return kTRUE;
  84    if(a>0.0) return kFALSE;
  85    a = s1->etar - s2->etar;
  86    if(a<0.0) return kTRUE;
  87    if(a>0.0) return kFALSE;
  88    return kFALSE;
  89 }
  90
  91 Bool_t L_SortEta(const Hit_tl *s1,const Hit_tl *s2){
  92   // Eta sorting function for qsort.
  93    Float_t a;
  94
  95    a = s1->etar - s2->etar;
  96    if(a<0.0) return kTRUE;
  97    if(a>0.0) return kFALSE;
  98    a = s1->phir - s2->phir;
  99    if(a<0.0) return kTRUE;
 100    if(a>0.0) return kFALSE;
 101    return kFALSE;
 102 }
 103
 104 void hpsortPhi(Hit_tl **ra,Int_t n){
 105    Int_t i,ir,j,l;
 106    Hit_tl *rra;
 107
 108    if(n<2) return;
 109
 110    l  = ((n-1) >> 1) +1; // devide 2 + 1
 111    ir = n-1;
 112    for(;;){
 113      if(l>0){
 114         rra = ra[--l];  // decrament first
 115      }else{
 116         rra    = ra[ir];
 117         ra[ir] = ra[0];
 118         if(--ir == 0){  // decrament first
 119            ra[0] = rra;
 120            break;
 121         } // if --ra == 0
 122      } // end l>0
 123      i = l;
 124      j = l+1;
 125      while(j<=ir){
 126         if( j<ir && L_SortPhi(ra[j],ra[j+1]) ) j++;
 127         if( L_SortPhi(rra,ra[j]) ){
 128            ra[i] = ra[j];
 129            i = j;
 130            j <<= 1; // time 2.
 131         }else{
 132            break;
 133         } // end if func()
 134      } // end while
 135      ra[i] = rra;
 136    } // end for ever
 137 }
 138
 139
 140 void hpsortEta(Hit_tl **ra,Int_t n){
 141    Int_t i,ir,j,l;
 142    Hit_tl *rra;
 143
 144    if(n<2) return;
 145
 146    l  = ((n-1) >> 1) +1; // devide 2 + 1
 147    ir = n-1;
 148    for(;;){
 149      if(l>0){
 150         rra = ra[--l];  // decrament first
 151      }else{
 152         rra    = ra[ir];
 153         ra[ir] = ra[0];
 154         if(--ir == 0){  // decrament first
 155            ra[0] = rra;
 156            break;
 157         } // if --ra == 0
 158      } // end l>0
 159      i = l;
 160      j = l+1;
 161      while(j<=ir){
 162         if( j<ir && L_SortEta(ra[j],ra[j+1]) ) j++;
 163         if( L_SortEta(rra,ra[j]) ){
 164            ra[i] = ra[j];
 165            i = j;
 166            j <<= 1; // time 2.
 167         }else{
 168            break;
 169         } // end if func()
 170      } // end while
 171      ra[i] = rra;
 172    } // end for ever
 173 }
 174
 175 void fillStructure(Hit_tl **spd,Int_t &nspd,
 176                    Float_t *xv,Int_t *id,Int_t track,Int_t nhits){
 177     Float_t PI2 = 2.0*TMath::Pi();
 178 //    Float_t PI  = TMath::Pi();
 179     Int_t i;
 180     Float_t x,y,z,zr,r,phi,eta,rphi;
 181     Float_t dz    = 0.0300; // 300 microns
 182     Float_t drphi = 0.0050; //  50 microns
 183
 184     i = nspd;
 185 //   if(nspd<5) printf("fill: spd=%p spd[%d]=%p i=%d spd[i]=%p id=%d %d %d\n",
 186 //                    spd,nspd,spd[nspd],i,spd[i],id[0],id[1],id[2]);
 187     x = xv[0];
 188     y = xv[1];
 189     z = xv[2];
 190     r = sqrt(x*x+y*y);
 191     phi = atan2(y,x);
 192     if(phi<0.0) phi += PI2;
 193     rphi = r*phi;
 194     eta  = -0.5*tan(0.5*atan2(r,z));
 195
 196     spd[i]->track = track;
 197     spd[i]->n     = nhits;
 198     spd[i]->lad   = id[1];
 199     spd[i]->det   = id[2];
 200     spd[i]->x     = x;
 201     spd[i]->y     = y;
 202     spd[i]->z     = z;
 203     spd[i]->r     = r;
 204     spd[i]->phi   = phi;
 205     spd[i]->eta   = eta;
 206
 207     zr   = dz    * Float_t(Int_t(z   /   dz)) + 0.5*dz;
 208     rphi = drphi * Float_t(Int_t(rphi/drphi)) + 0.5*drphi;
 209
 210     spd[i]->xr   =  cos(phi)*rphi/r;
 211     spd[i]->yr   = -sin(phi)*rphi/r;
 212     spd[i]->zr   = zr;
 213     spd[i]->phir = rphi/r;
 214     spd[i]->etar = -0.5*tan(0.5*atan2(r,zr));
 215
 216     nspd++;
 217     return;
 218 }
 219
 220 void HitsToV(Hit_tl **spdi,Int_t &nspdi,Hit_tl **spdo,Int_t &nspdo,
 221             Int_t ntracks,TTree *TH,TClonesArray *ITShits,
 222              Float_t fraction=1.0,Float_t len=16.82){
 223     Int_t     i,t,h,n,nb,nhits,id[4],idold[4],Iseed,ieta=0,itrk=0,ieta2=0;
 224     Float_t   x[3],xb[3],xl[3];
 225     AliITShit *itsHit;
 226     TClonesArray *Part = gAlice->Particles();
 227     TParticle *part;
 228
 229     nspdi = nspdo = 0;
 230
 231     idold[0] = idold[1] = idold[2] = idold[3] = 0;
 232     xb[0]    = xb[1]    = xb[2]    = 0.0;
 233     n = 0;
 234     Iseed = time(0);
 235     printf("HitsToV: Iseed=%d fraction=%f Pixel length=%fcm\n",
 236            Iseed,fraction,len);
 237     gRandom->SetSeed(Iseed);
 238 //    printf("HitsToV: gRandom->Rndm(1)=%f\n",gRandom->Rndm(1));
 239     for(t=0;t<ntracks;t++){
 240         if(fraction<gRandom->Rndm(1)) continue; // skip some tracks
 241         itrk++;
 242         gAlice->ResetHits();
 243         nb       = TH->GetEvent(t);
 244         nhits    = ITShits->GetEntriesFast();
 245         for(h=0;h<nhits;h++){
 246             itsHit = (AliITShit *) ITShits->UncheckedAt(h);
 247             itsHit->GetDetectorID(id[0],id[1],id[2]); id[3]=t;
 248             itsHit->GetPositionG(x[0],x[1],x[2]);
 249             itsHit->GetPositionL(xl[0],xl[1],xl[2]);
 250             if(h==0){
 251                 part = (TParticle *) (Part->UncheckedAt(itsHit->GetTrack()));
 252                 if(TMath::Abs(part->Eta())<=1.0) ieta++;
 253                 if(TMath::Abs(part->Eta())<=0.5) ieta2++;
 254             } // end if h==0
 255             if(TMath::Abs(x[2]/len) >= 1.0) continue;
 256             if(x[0]==0.0&&x[1]==0.0){
 257                 printf("Hitsto: t=%d h=%d/%d id=%d,%d,%d x=%f,%f,%f\n",
 258                        t,h,nhits,id[0],id[1],id[2],x[0],x[1],x[2]);
 259                 continue;
 260             } // end if x[0]==x[1]==0.0
 261             if(!(id[0]==idold[0]&&id[1]==idold[1]&&
 262                  id[2]==idold[2]&&id[3]==idold[3])){
 263                 if(!(n<=0 || (xb[0]==0.0&&xb[1]==0.0))){
 264                     for(i=0;i<3;i++) xb[i]   /= Float_t(n);
 265                     if(idold[0]==1){
 266                         fillStructure(spdi,nspdi,xb,idold,t,n);
 267                     }
 268                     if(idold[0]==2){
 269                         fillStructure(spdo,nspdo,xb,idold,t,n);
 270                     }
 271                     if(nspdi>ntracks || nspdo>ntracks){
 272                         printf("Hitsto: fill error,"
 273                                " nspdi=%d nspdo=%d ntracks=%d\n",
 274                                nspdi,nspdo,ntracks);
 275                     } // end if fill error
 276                 } // end if n>0
 277                 for(i=0;i<4;i++) idold[i] = id[i];
 278                 for(i=0;i<3;i++) xb[i]    = 0.0;
 279                 n = 0;
 280             } // end if id != idold
 281             for(i=0;i<3;i++) xb[i] += x[i];
 282             n++;
 283         } // end for h
 284     } // end for t
 285     printf("exiting HitsToV: %d primary tracks in eta=+-1\n",ieta);
 286     printf("exiting HitsToV: %d primary tracks #2 in eta=+-0.5\n",ieta2);
 287     printf("exiting HitsToV: %d primary tracks in file used\n",itrk);
 288     return;
 289 }
 290
 291 Float_t vertex(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
 292    Float_t r1    = 3.910078; // radius at which hit is from beam axis for spd1
 293    Float_t r2    = 6.955933; // radius at which hit is from beam axis for spd2
 294    Float_t DPhi  = 0.005;    // maximum allowed difference in phi angle
 295    Float_t DZ    = 12.5;      // maximum allowed difference in z position
 296    Float_t avt   = 0.0;
 297    Int_t   nbinx = 2000;
 298    Int_t   start = 0;
 299    Bool_t  mod   = kFALSE;
 300    Int_t   i,j;
 301    Float_t zv,av,su,i0,i1,i2,x0,x1,dphi,dz;
 302    Float_t PI2 = 2.0*TMath::Pi();
 303    Float_t PI  = TMath::Pi();
 304
 305    // sort according to phi.
 306    hpsortPhi(spdi,i1max);
 307    hpsortPhi(spdo,i2max);
 308
 309    // find best vertex allong z.
 310    TH2S *Pvzphi   = new TH2S("Pvzphi","Phi: Posible Z vertecies vs. Phi",
 311                               nbinx,-1.0,1.0,32,0.0,PI2);
 312    Pvzphi->Sumw2(); // collect good statitics
 313    TH1F *Pvzfl   = new TH1F("Pvzfl","Phi: Posible Z vertecies flattened",
 314                          nbinx,-1.0,1.0);
 315    Pvzfl->Sumw2();
 316    TH1F *Pvztr  = new TH1F("Pvztr","Phi: Z Vertex found for True Tracks",
 317                          nbinx,-1.0,1.0);
 318    Pvztr->Sumw2(); // collect good statitics
 319
 320    for(i=0;i<i1max;i++){
 321        if(spdi[i]->r==0.0) {printf("spdi[%d]->r=0.0\n",i);continue;}
 322        for(;spdo[start]->phir<spdi[i]->phir-1.5*DPhi;start++);
 323        for(j=start;(spdo[j]->phir < spdi[i]->phir+DPhi) && (j<i2max);j++){
 324            dphi = fabs(spdo[j]->phir - spdi[i]->phir); if(dphi>PI) dphi -= PI;
 325            dz   = fabs(spdo[j]->zr   - spdi[i]->zr);
 326            if(dphi>DPhi) continue;
 327            if(dz>DZ)     continue; // If outside dz range skip it
 328            r1   = spdi[i]->r;
 329            r2   = spdo[j]->r;
 330            if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
 331            else{
 332                printf("vertex: spdi[%d]->r=%f = spdo[%d]->r=%f\n",i,r1,j,r2);
 333                continue;
 334            } // end if else r2-r1!=0.0
 335            Pvzphi->Fill(zv,spdi[i]->phir);
 336            Pvzfl->Fill(zv);
 337            if(spdi[i]->track == spdo[j]->track) Pvztr->Fill(zv);
 338        } // end for j
 339    } // end for i
 340
 341    TH1D *Pvzdef  = Pvzphi->ProjectionX("Phi: Posible Z vertecies",1,nbinx,"E");
 342
 343    i0 = Pvzfl->GetBinContent(1);
 344    x0 = Float_t(1);
 345    i1 = Pvzfl->GetBinContent(nbinx);
 346    x1 = Float_t(nbinx);
 347    su = x0-x1; if(su==0.0) return -200.0;
 348    su = (i0-i1)/su;
 349    x0 = x0-su*i0;
 350    x1 = su;
 351    for(i=1;i<=nbinx;i++){
 352       su = x1*Float_t(i) + x0;
 353       Pvzfl->AddBinContent(i,-su);
 354    } // end for i
 355
 356    printf("Phi:            mean=%f RMS=%f w2=%f\n",
 357                Pvzdef->GetMean(),Pvzdef->GetRMS(),Pvzdef->GetSumOfWeights());
 358    printf("Phi: Flattened mean=%f RMS=%f w2=%f\n",
 359                   Pvzfl->GetMean(),Pvzfl->GetRMS(),Pvzfl->GetSumOfWeights());
 360
 361    av = 3.0 * Pvzfl->Integral(1,nbinx)/Float_t(nbinx);
 362    printf("Phi: Flattened av=%f Pvzfl->Max=%f nbinx=%d\n",
 363           av,Pvzfl->GetMaximum(),nbinx);
 364
 365    su  = 0.0;
 366    avt = 0.0;
 367    i1  = Pvzfl->GetBinContent(1);
 368    i2  = Pvzfl->GetBinContent(2);
 369    for(i=2;i<nbinx;i++){
 370        i0 = i1; // old i-1 value
 371        i1 = i2; // old i value
 372        i2 = Pvzfl->GetBinContent(i+1); // new i+1 value
 373        if(i1 > av && (i0 > av || i2 > av ) ){
 374            if(!mod) mod = kTRUE;
 375        } else if(mod) break;
 376        if(mod){ // inside peak
 377            su  += i1;
 378            avt += i1 * Pvzfl->GetXaxis()->GetBinCenter(i);
 379        } // end if mod
 380    } // end for i
 381
 382    if(su!=0.0) zv = avt/su;
 383    else        zv = -100.0; // an unphysical value
 384
 385    TCanvas *c0 = new TCanvas("c0","Alice ITS vertex finder", 400,10,600,700);
 386    Pvzphi->Draw("col");
 387    c0->Print("vertex5_P_vz_phi.ps");
 388    Pvzdef->Draw();
 389    c0->Print("vertex5_P_vz_def.ps");
 390    Pvzfl->Draw();
 391    c0->Print("vertex5_P_vz_flat.ps");
 392    Pvztr->Draw();
 393    c0->Print("vertex5_P_vz_true.ps");
 394
 395    return zv;
 396 }
 397
 398 void Chi2Gauss(Int_t &npar,Double_t *gin,Double_t &f,
 399                Double_t *par,Int_t iflag){
 400     Double_t chi2 = 0.0;
 401     Double_t delta,h,x,eh;
 402
 403     for(Int_t i=gstrt;i<gend;i++){
 404         h = gH1->GetBinContent(i);
 405         eh = TMath::Sqrt(h);
 406         if(eh <= 0.0) eh = 1.0;
 407         x = gH1->GetXaxis()->GetBinCenter(i);
 408         delta = (h - par[0]*TMath::Gaus(x,par[1],par[2]) - par[3])/eh;
 409         chi2 += delta*delta;
 410     } // end for i
 411     f = chi2;
 412     return;
 413 }
 414
 415 Float_t vertexSlow(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
 416    Float_t r1    = 3.910078; // radius at which hit is from beam axis for spd1
 417    Float_t r2    = 6.955933; // radius at which hit is from beam axis for spd2
 418    Float_t DPhi  = 0.005;    // maximum allowed difference in phi angle
 419    Float_t BDphi;
 420    Float_t DZ    = 12.5;      // maximum allowed difference in z position
 421    Float_t avt   = 0.0;
 422    Int_t   nbinx = 2000;
 423    Int_t   nbst;
 424    Bool_t  mod   = kFALSE;
 425    Int_t   i,j;
 426    Float_t zv,av,su,i0,i1,i2,x0,x1,dphi,dz;
 427    Float_t PI2 = 2.0*TMath::Pi();
 428    Float_t PI  = TMath::Pi();
 429
 430    if(i1max<=0||i2max<=0) return -1.0E10;
 431
 432    // find best vertex allong z.
 433    TH2S *Svzphi  = new TH2S("Svzphi","Slow: Posible Z vertecies vs. Phi",
 434                             nbinx,+0.0,10.0,32,0.0,PI2);
 435    Svzphi->Sumw2(); // collect good statitics
 436    TH2S *Svzdphi  = new TH2S("Svzdpii","Slow: Posible Z vertecies vs. DPhi",
 437                              200,+4.0,6.0,20,0.0,10*DPhi);
 438    Svzdphi->Sumw2(); // collect good statitics
 439    TH1F *Svzfl  = new TH1F("Svzfl","Slow: Posible Z vertecies Flattened",
 440                             nbinx,+4.0,6.0);
 441    Svzfl->Sumw2();
 442    TH1F *Svztr = new TH1F("Svztr","Slow: Z Vertex found by True Tracks",
 443                           nbinx,+4.0,6.0);
 444    Svztr->Sumw2(); // collect good statitics
 445
 446    printf("Svertex: i1max=%d i2max=%d\n",i1max,i2max);
 447
 448    for(i=0;i<i1max;i++) for(j=0;j<i2max;j++) {
 449        dphi = fabs(spdo[j]->phir - spdi[i]->phir); if(dphi>PI) dphi -= PI;
 450        dz   = fabs(spdo[j]->zr   - spdi[i]->zr);
 451        r1   = spdi[i]->r;
 452        r2   = spdo[j]->r;
 453        if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
 454        else{
 455            printf("vertex_slow: spdi[%d]->r=%f = spdo[%d]->r=%f\n",i,r1,j,r2);
 456            continue;
 457        } // end if else r1-r2!=0.0
 458 //       if(j<10&&i<10) printf("zv=%e dphi=%e,r1=%e,r2=%e,dz=%e\n",
 459 //                           zv,dphi,r1,r2,dz);
 460        Svzdphi->Fill(zv,dphi);
 461        if(dphi>DPhi) continue; // If outside DPhi (momentum) range, skip it.
 462        if(dz>DZ)     continue; // If outside dz range, skip it.
 463        Svzphi ->Fill(zv,spdi[i]->phir);
 464        Svzfl->Fill(zv);
 465        if(spdi[i]->track == spdo[j]->track) Svztr->Fill(zv);
 466    } // end for i,j
 467
 468    TH1D *Svzdef = Svzphi->ProjectionX("Slow: Posible Z vertecies",0,nbinx,"E");
 469
 470    i0 = Svzfl->GetBinContent(1);
 471    x0 = Float_t(1);
 472    i1 = Svzfl->GetBinContent(nbinx);
 473    x1 = Float_t(nbinx);
 474    su = x0-x1; if(su==0.0) return -200.0;
 475    su = (i0-i1)/su;
 476    x0 = x0-su*i0;
 477    x1 = su;
 478    for(i=1;i<=nbinx;i++){
 479       su = x1*Float_t(i) + x0;
 480       Svzfl->AddBinContent(i,-su);
 481    } // end for i
 482
 483    printf("Slow:           mean=%f RMS=%f w2=%f\n",
 484           Svzdef->GetMean(),Svzdef->GetRMS(),Svzdef->GetSumOfWeights());
 485    printf("Slow: Flattened mean=%f RMS=%f w2=%f\n",
 486           Svzfl->GetMean(),Svzfl->GetRMS(),Svzfl->GetSumOfWeights());
 487
 488    av = 3.0 * Svzfl->Integral(1,nbinx)/Float_t(nbinx);
 489    printf("Slow: Flattened av=%f Tvxps->Max=%f nbinx=%d\n",
 490           av,Svzfl->GetMaximum(),nbinx);
 491
 492    {// find best dPhi that masimizes the signal/noise ration.
 493        Float_t sn=0.0,sig=0.0,nois=0.0,ns;
 494        Int_t   iend = Svzdphi->GetYaxis()->GetNbins();
 495        Int_t   jend = Svzdphi->GetXaxis()->GetNbins();
 496        Int_t   ipeak = Svzdphi->GetXaxis()->FindBin(Svzdef->GetMean());
 497        nbst = 0;
 498        for(i=0;i<=iend;i++){
 499            sig += Svzdphi->GetCellContent(ipeak,i);
 500            ns   = 0.0;
 501            for(j=1;j<6;j++) ns += Svzdphi->GetCellContent(j,i);
 502            for(j=jend-6;j<jend;j++) ns += Svzdphi->GetCellContent(j,i);
 503            nois += 0.1*ns;
 504            if(nois<=0.0) continue;
 505            if((sig-nois)/nois>sn){
 506                sn = (sig-nois)/nois;
 507                nbst = i;
 508            } // end if
 509 //         printf("Svertex: bin=%d signal/noise=%e\n",i,(sig-nois)/nois);
 510        } // end for i
 511    } // end find best dPhi
 512
 513    if(nbst<=0) nbst = 1; // must have some data.
 514    BDphi = Svzdphi->GetYaxis()->GetBinUpEdge(nbst);
 515    TH1D *Svzbst = Svzdphi->ProjectionX("Slow: Best, Z vertecies",0,nbst,"E");
 516
 517    { // Start Muinuit fitting
 518        // initilize Minuit package
 519        gH1   = Svzbst; // histogram to be fit
 520        gstrt = Svzbst->GetXaxis()->FindBin(
 521         Svzdef->GetMean()-2.0*(Svzdef->GetRMS()));//histogram start bin value
 522        gend  = Svzbst->GetXaxis()->FindBin(
 523         Svzdef->GetMean()+2.0*(Svzdef->GetRMS()));//histogram end   bin value
 524        TMinuit *gMinuit = new TMinuit(4); // init Minuit
 525        gMinuit->SetFCN(Chi2Gauss);  // chi^2 function with Gaussian
 526        Double_t arglist[10];  // Munuit parameter array
 527        Int_t ierflg = 0; // Muniut error flag
 528        //arglist[0] = 1.;
 529        //gMinuit->mnexcm("SET ERR",arglist,1,ierglg);
 530        // Set starting values and step size for parameters
 531        Double_t vstart[4],step[4];
 532        { // find background
 533            Float_t ns   = 0.0;
 534            Int_t jend = Svzbst->GetXaxis()->GetNbins();
 535            for(j=1;j<6;j++) ns += Svzbst->GetBinContent(j);
 536            for(j=jend-6;j<jend;j++) ns += Svzbst->GetBinContent(j);
 537            vstart[3] = 0.1*ns;
 538        } // end find backgrount
 539        vstart[1] = Svzbst->GetMean();
 540        vstart[2] = Svzbst->GetRMS();
 541        if(vstart[2] == 0.0) vstart[2] = 0.04;
 542        vstart[0] = (Svzbst->GetEntries() - (Svzbst->GetNbinsX())*vstart[3])*
 543            vstart[2]/TMath::Sqrt(TMath::Pi());
 544        if(vstart[0]<=0.0) vstart[0] = 1.0;
 545        for(i=0;i<4;i++) step[i] = 0.05*vstart[i];
 546        if(vstart[3] <= 0.0) step[3] = 0.1;
 547        step[1] = 0.01; // mean expected at about zero set step my hand
 548        gMinuit->mnparm(0,"Const",vstart[0],step[0],0,0,ierflg);
 549        gMinuit->mnparm(1,"Mean" ,vstart[1],step[1],0,0,ierflg);
 550        gMinuit->mnparm(2,"Sigma",vstart[2],step[2],0,0,ierflg);
 551        gMinuit->mnparm(3,"Offst",vstart[3],step[3],0,0,ierflg);
 552        // Now ready for minimization step
 553        //arglist[0] = 500.; // Maximum number of calls
 554        //arglist[1] = 1.;   // Tolorance
 555        gMinuit->mnexcm("MIGRAD",arglist,0,ierflg); // do minimization
 556        gMinuit->mnexcm("MINO",arglist,0,ierflg);   // find best errors
 557        // Get results
 558        Double_t parmin,edm,errdef;
 559        Int_t nvpar,nparx,icstat;
 560        gMinuit->mnstat(parmin,edm,errdef,nvpar,nparx,icstat);
 561        printf("Svertex: chi2gauss=%e edist=%e istat=%d dPhi=%e\n",
 562               parmin,edm,icstat,BDphi);
 563        TString chnam = TString(10);
 564        Double_t par[4],epar[4],empar[4],eppar[4];
 565        for(i=0;i<4;i++){
 566            gMinuit->mnpout(i,chnam,par[i],epar[i],empar[i],eppar[i],ierflg);
 567            gMinuit->mnerrs(i,eppar[i],empar[i],epar[i],arglist[i]);
 568            printf("Svertex: %s = %e +- %e (%e,%e)\n",
 569                   chnam.Data(),par[i],epar[i],empar[i],eppar[i]);
 570        } // end for i
 571    } // End Muinuit fitting
 572    su  = 0.0;
 573    avt = 0.0;
 574    i1  = Svzfl->GetBinContent(1);
 575    i2  = Svzfl->GetBinContent(2);
 576    for(i=2;i<nbinx;i++){
 577        i0 = i1; // old i-1 value
 578        i1 = i2; // old i value
 579        i2 = Svzfl->GetBinContent(i+1); // new i+1 value
 580        if(i1 > av && (i0 > av || i2 > av ) ){
 581            if(!mod) mod = kTRUE;
 582        } else if(mod) break;
 583        if(mod){ // inside peak
 584            su  += i1;
 585            avt += i1 * Svzfl->GetXaxis()->GetBinCenter(i);
 586        } // end if mod
 587    } // end for i
 588
 589    if(su!=0.0) zv = avt/su;
 590    else        zv = -100.0; // an unphysical value
 591
 592    TCanvas *c1 = new TCanvas("c1","Slow Alice ITS vertex finder",
 593                              450,10,600,700);
 594    Svzphi->Draw("col");
 595    c1->Print("vertex5_S_vz_phi.eps");
 596    Svzdphi->Draw();
 597    c1->Print("vertex5_S_vz_dphi.eps");
 598    Svzdef->Draw();
 599    c1->Print("vertex5_S_vz_def.eps");
 600    Svzbst->Draw();
 601    c1->Print("vertex5_S_vz_bst.eps");
 602    Svztr->Draw();
 603    c1->Print("vertex5_S_vz_true.eps");
 604
 605    return zv;
 606 }
 607
 608 Float_t vertexEta(Hit_tl **spdi,Int_t i1max,Hit_tl **spdo,Int_t i2max){
 609    Float_t r1    = 3.910078;// radius at which hit is from beam axis for spd1
 610    Float_t r2    = 6.955933;// radius at which hit is from beam axis for spd2
 611    Float_t DPhi  = 0.005;   // maximum allowed difference in phi angle
 612    Float_t DEta  = 0.100;   // maximum allowed difference in eta/pseudorapidity
 613    Float_t DZ    = 12.5;     // maximum allowed difference in z position
 614    Float_t avt   = 0.0;
 615    Int_t   nbinx = 2000;
 616    Int_t   start = 0;
 617    Bool_t  mod   = kFALSE;
 618    Int_t   i,j;
 619    Float_t zv=0.0,av,su,i0,i1,i2,x0,x1,dphi,dz;
 620    Float_t PI2 = 2.0*TMath::Pi();
 621    Float_t PI  = TMath::Pi();
 622
 623    // sort according to phi.
 624    hpsortEta(spdi,i1max);
 625    hpsortEta(spdo,i2max);
 626
 627    // find best vertex allong z.
 628    TH2S *Evzphi   = new TH2S("Evzphi","Eta: Posible Z vertecies vs. Phi",
 629                              nbinx,-5.0,5.0,32,0.0,PI2);
 630    Evzphi->Sumw2(); // collect good statitics
 631    TH2S *Evzdphi  = new TH2S("Evzdphi","Eta: Posible Z vertecies vs. DPhi",
 632                              200,-1.0,1.0,20,0.0,10*DPhi);
 633    Evzdphi->Sumw2(); // collect good statitics
 634    TH1F *Evzfl   = new TH1F("Evzfl","Eta: Posible Z vertecies Flattened",
 635                           nbinx,-5.0,5.0);
 636    Evzfl->Sumw2();
 637    TH1F *Evztr  = new TH1F("Evztr","Eta: Z Vertex found by True Tracks",
 638                            nbinx,-5.0,5.0);
 639    Evztr->Sumw2(); // collect good statitics
 640
 641    for(i=0;i<i1max;i++){
 642        for(;spdo[start]->etar < spdi[i]->etar-1.5*DEta;start++);
 643        for(j=start;(spdo[j]->etar < spdi[i]->etar+DEta) && (j < i2max);j++){
 644            dphi = fabs(spdi[i]->phir - spdo[j]->phir); if(dphi>PI) dphi -= PI;
 645            dz   = fabs(spdi[i]->zr   - spdo[j]->zr);
 646            r1   = spdi[i]->r;
 647            r2   = spdo[j]->r;
 648            if(r2-r1!=0.0) zv   = (r2*spdi[i]->zr - r1*spdo[j]->zr)/(r2-r1);
 649            else printf("vertex_Eta: spdi[%d]->r=%f = spdo[%d]->r=%f\n",
 650                        i,r1,j,r2);
 651            Evzdphi->Fill(zv,dphi);
 652            if(dphi>DPhi) continue;
 653            if(dz>DZ)     continue;
 654            Evzphi->Fill(zv,spdi[i]->phir);
 655            Evzfl->Fill(zv);
 656            if(spdi[i]->track == spdo[j]->track) Evztr ->Fill(zv);
 657        } // end for j
 658    } // end for i
 659
 660    TH1D *Evzdef  = Evzphi->ProjectionX("Eta: Z vertecies Eta",1,nbinx,"E");
 661
 662    i0 = Evzfl->GetBinContent(1);
 663    x0 = Float_t(1);
 664    i1 = Evzfl->GetBinContent(nbinx);
 665    x1 = Float_t(nbinx);
 666    su = x0-x1; if(su==0.0) return -200.0;
 667    su = (i0-i1)/su;
 668    x0 = x0-su*i0;
 669    x1 = su;
 670    for(i=1;i<=nbinx;i++){
 671       su = x1*Float_t(i) + x0;
 672       Evzfl->AddBinContent(i,-su);
 673    } // end for i
 674
 675    printf("Eta:           mean=%f RMS=%f w2=%f\n",
 676           Evzdef->GetMean(),Evzdef->GetRMS(),Evzdef->GetSumOfWeights());
 677    printf("Eta: Flattened mean=%f RMS=%f w2=%f\n",
 678           Evzfl->GetMean(),Evzfl->GetRMS(),Evzfl->GetSumOfWeights());
 679
 680    av = 3.0 * Evzfl->Integral(1,nbinx)/Float_t(nbinx);
 681    printf("Eta: Flattened av=%f TvxpE->Max=%f nbinx=%d\n",
 682           av,Evzfl->GetMaximum(),nbinx);
 683
 684    su  = 0.0;
 685    avt = 0.0;
 686    i1  = Evzfl->GetBinContent(1);
 687    i2  = Evzfl->GetBinContent(2);
 688    for(i=2;i<nbinx;i++){
 689        i0 = i1; // old i-1 value
 690        i1 = i2; // old i value
 691        i2 = Evzfl->GetBinContent(i+1); // new i+1 value
 692        if(i1 > av && (i0 > av || i2 > av ) ){
 693            if(!mod) mod = kTRUE;
 694        } else if(mod) break;
 695        if(mod){ // inside peak
 696            su  += i1;
 697            avt += i1 * Evzfl->GetXaxis()->GetBinCenter(i);
 698        } // end if mod
 699    } // end for i
 700
 701    if(su!=0.0) zv = avt/su;
 702    else        zv = -100.0; // an unphysical value
 703
 704    TCanvas *c2 = new TCanvas("c2","Alice ITS vertex finder Eta",
 705                              500,10,600,700);
 706    Evzphi->Draw("col");
 707    c2->Print("vertex5_E_vz_phi.ps");
 708    Evzdphi->Draw();
 709    c2->Print("vertex5_E_vz_dphi.ps");
 710    Evzdef->Draw();
 711    c2->Print("vertex5_E_vz_def.ps");
 712    Evzfl->Draw();
 713    c2->Print("vertex5_E_vz_flat.ps");
 714    Evztr->Draw();
 715    c2->Print("vertex5_E_vz_true.ps");
 716
 717    return zv;
 718 }