HLT/TPCLib/AliHLTTPCVertexArray.cxx

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