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 AliHLTTPCVertexArray::AliHLTTPCVertexArray()
  39   :
  40   fZSector(0),
  41   fZSectorErr(0),
  42   fMaxSeed(400),
  43   fNSeed(0)
  44 {
  45 }
  46
  47 AliHLTTPCVertexArray::AliHLTTPCVertexArray(AliHLTTPCVertexArray&)
  48   :
  49   fZSector(0),
  50   fZSectorErr(0),
  51   fMaxSeed(400),
  52   fNSeed(0)
  53 {
  54 }
  55
  56 AliHLTTPCVertexArray::AliHLTTPCVertexArray(Int_t maxseed)
  57   :
  58   fZSector(0),
  59   fZSectorErr(0),
  60   fMaxSeed(maxseed),
  61   fNSeed(0)
  62 {
  63 }
  64
  65 AliHLTTPCVertexArray::~AliHLTTPCVertexArray()
  66 {
  67 }
  68
  69 void AliHLTTPCVertexArray::AnalyzeSector(Float_t *vertex, Int_t *array, Int_t len)
  70 {
  71   //loop over all seeds and all vertex position
  72   LOG(AliHLTTPCLog::kInformational,"AliHLTTPCVertexArray::AnalyzeSector","Analyze")
  73   <<AliHLTTPCLog::kDec<<"Number of Seeds: "<<fNSeed<<ENDLOG;
  74   for(Int_t i =0;i<fNSeed;i++)
  75     for(Int_t bin = 0;bin<len;bin++)
  76       array[bin] += Trace(fZSeed[i],fRSeed[i],fSecSeed[i],vertex[bin]);
  77 }
  78
  79 void AliHLTTPCVertexArray::FindSectorVertex(Double_t pos, Double_t range, Int_t nbin)
  80 {
  81   //define position and range for search and
  82   //loop over all seeds
  83   //and find position of vertex and error
  84   const Int_t klen = nbin;
  85   const Double_t kwidth = range;
  86   const Double_t kxmin = pos - kwidth/2;
  87   const Double_t kstep = kwidth/klen;
  88   const Double_t kstart = kxmin + kstep/2.;
  89   Int_t * array = new Int_t[klen];
  90   Float_t * ver   = new Float_t[klen];
  91   for(Int_t i=0;i<klen;i++){
  92     ver[i] =  kstart + kstep * i;
  93     array[i] = 0;
  94   }
  95   AnalyzeSector(ver,array,klen);
  96   FindMean(ver,array,klen);
  97   delete[] array;
  98   delete[] ver;
  99 }
 100
 101 void AliHLTTPCVertexArray::FindMean(Float_t *vertex,Int_t *array, Int_t len)
 102 {
 103   //find mean and error of array and store it in
 104   //fZSector and fZSectorErr
 105   const Int_t knbin = len;
 106   Int_t xbin =0;
 107   Int_t max=0;
 108   for(Int_t i = 0;i<knbin;i++){
 109     if(array[i]>max){
 110       max = array[i];
 111       xbin =i;
 112     }
 113   }
 114   Int_t hmax = max/2;
 115   Int_t xmin,xmax;
 116   Int_t ops = 0;
 117   xmin = xbin;
 118   while(xmin--){
 119     if(xmin<0) {ops++;break;}
 120     if(array[xmin]<hmax) {
 121       break;
 122     }
 123   }
 124   xmax = xbin;
 125   while(xmax++){
 126     if(xmax>=knbin) {ops++;break;}
 127     if(array[xmax]<hmax){
 128       break;
 129     }
 130   }
 131   if(ops){
 132     if(xbin >= knbin/2){xmin = 2 * xbin - knbin +1;xmax = knbin-1;}
 133     else{xmin = 0;xmax = 2 * xbin;}
 134   }
 135   Double_t sumw=0;
 136   Double_t sumw2=0;
 137   Double_t sumwx=0;
 138   Double_t sumwx2=0;
 139   for(Int_t bin = xmin;bin<=xmax;bin++){
 140     sumw   += array[bin];
 141     sumw2  += array[bin] * array[bin];
 142     sumwx  += array[bin] * vertex[bin];
 143     sumwx2 += array[bin] * vertex[bin] * vertex[bin];
 144   }
 145   if(sumw){
 146     Double_t mean = sumwx/sumw;
 147     Double_t rms2 = fabs(sumwx2/sumw - mean*mean);
 148     fZSectorErr = sqrt(rms2/sumw);
 149     fZSector = mean;
 150   }
 151   else{fZSectorErr = fZSector = 0;}
 152   sumw=sumw2=sumwx=sumwx2=0;
 153   for(Int_t bin = 0;bin<knbin;bin++){
 154     sumw   += array[bin];
 155     sumw2  += array[bin] * array[bin];
 156     sumwx  += array[bin] * vertex[bin];
 157     sumwx2 += array[bin] * vertex[bin] * vertex[bin];
 158   }
 159   if(sumw){
 160     Double_t mean = fZSector;
 161     Double_t rms2 = fabs(sumwx2/sumw - mean*mean);
 162     fZSectorErr = sqrt(rms2/sumw);
 163   }
 164 }
 165