HLT/TPCLib/AliHLTTPCVertexArray.cxx

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