HLT/src/AliL3VertexArray.cxx

   1 // Author: Uli Frankenfeld <mailto:franken@fi.uib.no>
   2 //*-- Copyright &copy Uli
   3
   4 #include "AliL3Logging.h"
   5 //#include <TMath.h>
   6 #include <iostream.h>
   7 #include "AliL3VertexArray.h"
   8
   9
  10 //_____________________________________________________________
  11 // AliL3VertexArray
  12 //
  13 // The L3 Fast Vertex Finder Base Class
  14 //
  15 // usage:
  16 //
  17 //for(Int_t sec=0;sec<NSEC;sec++){
  18 //  ResetSector();
  19 //  FillSectorSeed3D(x,y,z);
  20 //  FillSector3D(x,y,z);
  21 //  FindSectorVertex();
  22 //  Double_t z = GetZSector();
  23 //  Double_t zerr = GetZSectorErr();
  24 //// do somethink with z, zerr
  25 //}
  26
  27 ClassImp(AliL3VertexArray)
  28 void AliL3VertexArray::AnalyzeSector(Float_t *vertex, Int_t *array, Int_t len){
  29 //loop over all seeds and all vertex position
  30   LOG(AliL3Log::kInformational,"AliL3VertexArray::AnalyzeSector","Analyze")
  31   <<AliL3Log::kDec<<"Number of Seeds: "<<fNSeed<<ENDLOG;
  32   for(Int_t i =0;i<fNSeed;i++)
  33     for(Int_t bin = 0;bin<len;bin++)
  34       array[bin] += Trace(fZSeed[i],fRSeed[i],fSecSeed[i],vertex[bin]);
  35 }
  36
  37 void AliL3VertexArray::FindSectorVertex(Double_t pos, Double_t range, Int_t nbin){
  38 //define position and range for search and
  39 //loop over all seeds
  40 //and find position of vertex and error
  41   const Int_t len = nbin;
  42   const Double_t width = range;
  43   const Double_t  xmin = pos - width/2;
  44   const Double_t step = width/len;
  45   const Double_t start = xmin + step/2.;
  46   Int_t * array = new Int_t[len];
  47   Float_t * ver   = new Float_t[len];
  48   for(Int_t i=0;i<len;i++){
  49     ver[i] =  start + step * i;
  50     array[i] = 0;
  51   }
  52   AnalyzeSector(ver,array,len);
  53   FindMean(ver,array,len);
  54   delete[] array;
  55   delete[] ver;
  56 }
  57
  58 void AliL3VertexArray::FindMean(Float_t *vertex,Int_t *array, Int_t len){
  59 //find mean and error of array and store it in
  60 //fZSector and fZSectorErr
  61   const Int_t nbin = len;
  62   Int_t xbin =0;
  63   Int_t max=0;
  64   for(Int_t i = 0;i<nbin;i++){
  65     if(array[i]>max){
  66       max = array[i];
  67       xbin =i;
  68     }
  69   }
  70   Int_t hmax = max/2;
  71   Int_t xmin,xmax;
  72   Int_t ops = 0;
  73   xmin = xbin;
  74   while(xmin--){
  75     if(xmin<0) {ops++;break;}
  76     if(array[xmin]<hmax) {
  77       break;
  78     }
  79   }
  80   xmax = xbin;
  81   while(xmax++){
  82     if(xmax>=nbin) {ops++;break;}
  83     if(array[xmax]<hmax){
  84       break;
  85     }
  86   }
  87   if(ops){
  88     if(xbin >= nbin/2){xmin = 2 * xbin - nbin +1;xmax = nbin-1;}
  89     else{xmin = 0;xmax = 2 * xbin;}
  90   }
  91   Double_t sumw=0;
  92   Double_t sumw2=0;
  93   Double_t sumwx=0;
  94   Double_t sumwx2=0;
  95   for(Int_t bin = xmin;bin<=xmax;bin++){
  96     sumw   += array[bin];
  97     sumw2  += array[bin] * array[bin];
  98     sumwx  += array[bin] * vertex[bin];
  99     sumwx2 += array[bin] * vertex[bin] * vertex[bin];
 100   }
 101   if(sumw){
 102     Double_t mean = sumwx/sumw;
 103     Double_t rms2 = fabs(sumwx2/sumw - mean*mean);
 104     fZSectorErr = sqrt(rms2/sumw);
 105     fZSector = mean;
 106   }
 107   else{fZSectorErr = fZSector = 0;}
 108   sumw=sumw2=sumwx=sumwx2=0;
 109   for(Int_t bin = 0;bin<nbin;bin++){
 110     sumw   += array[bin];
 111     sumw2  += array[bin] * array[bin];
 112     sumwx  += array[bin] * vertex[bin];
 113     sumwx2 += array[bin] * vertex[bin] * vertex[bin];
 114   }
 115   if(sumw){
 116     Double_t mean = fZSector;
 117     Double_t rms2 = fabs(sumwx2/sumw - mean*mean);
 118     fZSectorErr = sqrt(rms2/sumw);
 119   }
 120 }
 121