HLT/hough/AliL3HistogramAdaptive.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <string.h>
   7 #include "AliL3Logging.h"
   8 #include "AliL3HistogramAdaptive.h"
   9 #include "AliL3Transform.h"
  10 #include "AliL3Track.h"
  11
  12 //_____________________________________________________________
  13 // AliL3HistogramAdaptive
  14 //
  15 // 2D histogram class
  16
  17 ClassImp(AliL3HistogramAdaptive)
  18
  19 AliL3HistogramAdaptive::AliL3HistogramAdaptive() : AliL3Histogram()
  20 {
  21
  22 }
  23
  24
  25 AliL3HistogramAdaptive::AliL3HistogramAdaptive(Char_t *name,Double_t minpt,Double_t maxpt,Double_t ptres,
  26                                                Int_t nybins,Double_t ymin,Double_t ymax)
  27 {
  28   strcpy(fName,name);
  29
  30   fPtres = ptres;
  31   fXmin = -1*AliL3Transform::GetBFact()*AliL3Transform::GetBField()/minpt;
  32   fXmax = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/minpt;
  33
  34   fPtstep = 0.1;
  35   fMinPt = minpt;
  36   fMaxPt = maxpt;
  37   fNxbins = InitPtBins();
  38   //cout<<"Setting "<<fNxbins<<" bins on x"<<endl;
  39
  40   fNybins = nybins;
  41   fYmin = ymin;
  42   fYmax = ymax;
  43   fFirstXbin=1;
  44   fFirstYbin=1;
  45   fLastXbin = fNxbins;
  46   fLastYbin = fNybins;
  47   fNcells = (fNxbins+2)*(fNybins+2);
  48
  49   fThreshold=0;
  50   fContent = new Int_t[fNcells];
  51   Reset();
  52 }
  53
  54 AliL3HistogramAdaptive::~AliL3HistogramAdaptive()
  55 {
  56
  57 }
  58
  59 Int_t AliL3HistogramAdaptive::InitPtBins()
  60 {
  61
  62   Double_t pt = fMinPt,delta_pt,local_pt;
  63   Int_t bin=0;
  64
  65   while(pt < fMaxPt + fPtstep)
  66     {
  67       local_pt = pt;
  68       delta_pt = fPtres*(local_pt+0.5*fPtstep);
  69
  70       while(local_pt < pt + fPtstep)
  71         {
  72           if(local_pt > fMaxPt) break;
  73           local_pt += delta_pt;
  74           bin++;
  75         }
  76       //cout<<"Setting "<<bin<<" at step "<<pt<<endl;
  77       pt += fPtstep;
  78     }
  79
  80   return (bin+1)*2; //Both negative and positive kappa.
  81 }
  82
  83
  84 void AliL3HistogramAdaptive::Fill(Double_t x,Double_t y,Int_t weight)
  85 {
  86   Int_t bin = FindBin(x,y);
  87   if(bin < 0)
  88     return;
  89   AddBinContent(bin,weight);
  90
  91 }
  92
  93 Int_t AliL3HistogramAdaptive::FindBin(Double_t x,Double_t y)
  94 {
  95
  96   Int_t xbin = FindXbin(x);
  97   Int_t ybin = FindYbin(y);
  98
  99   if(xbin < 0)
 100     return -1;
 101   return GetBin(xbin,ybin);
 102 }
 103
 104 Int_t AliL3HistogramAdaptive::FindXbin(Double_t x)
 105 {
 106
 107   Double_t ptfind = fabs(AliL3Transform::GetBFact()*AliL3Transform::GetBField()/x);
 108   if(ptfind < fMinPt || ptfind > fMaxPt) return -1;
 109
 110   Double_t delta_pt,ptplus1,pt = fMinPt,local_pt;
 111   Bool_t found=kFALSE;
 112   Int_t bin=0;
 113   while(pt < fMaxPt + fPtstep)
 114     {
 115       local_pt = pt;
 116       delta_pt = fPtres*(local_pt+0.5*fPtstep);
 117
 118       while(local_pt < pt + fPtstep)
 119         {
 120           ptplus1 = local_pt + delta_pt;
 121           if(ptfind >= local_pt && ptfind < ptplus1)
 122             {
 123               found=kTRUE;
 124               break;
 125             }
 126           local_pt = ptplus1;
 127           bin++;
 128         }
 129       if(found) break;
 130       pt += fPtstep;
 131     }
 132   if(bin >= fNxbins/2)
 133     cerr<<"AliL3HistogramAdaptive::FindXbin : Bin out of range : "<<bin<<endl;
 134
 135   if(x < 0)
 136     return bin;
 137   else
 138     return fNxbins - 1 - bin;
 139 }
 140
 141 Int_t AliL3HistogramAdaptive::FindYbin(Double_t y)
 142 {
 143   if(y < fYmin || y > fYmax)
 144     return 0;
 145
 146   return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
 147
 148 }
 149
 150 Double_t AliL3HistogramAdaptive::GetBinCenterX(Int_t xbin)
 151 {
 152   if(xbin < 0 || xbin > fNxbins)
 153     {
 154       cerr<<"AliL3HistogramAdaptive::GetBinCenterX : Xbin out of range "<<xbin<<endl;
 155       return 0;
 156     }
 157   Double_t delta_pt,ptplus1,local_pt,pt=fMinPt;
 158   Bool_t found=kFALSE;
 159   Int_t bin=0;
 160   while(pt < fMaxPt + fPtstep)
 161     {
 162       local_pt = pt;
 163       delta_pt = fPtres*(local_pt+0.5*fPtstep);
 164       while(local_pt < pt + fPtstep)
 165         {
 166           ptplus1 = local_pt + delta_pt;
 167           if(xbin == bin || xbin == fNxbins - bin)
 168             {
 169               found=kTRUE;
 170               break;
 171             }
 172           local_pt = ptplus1;
 173           bin++;
 174         }
 175       if(found) break;
 176       pt += fPtstep;
 177     }
 178   Double_t binwidth = ptplus1 - local_pt;
 179   Double_t kappa = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/(local_pt + binwidth*0.5);
 180   if(xbin <= fNxbins/2 - 1)
 181     return -1.*kappa;
 182   else
 183     return kappa;
 184 }
 185
 186 Double_t AliL3HistogramAdaptive::GetBinCenterY(Int_t ybin)
 187 {
 188
 189   Double_t binwidth = (fYmax - fYmin) / fNybins;
 190   return fYmin + (ybin-1) * binwidth + 0.5*binwidth;
 191
 192 }
 193
 194 void AliL3HistogramAdaptive::Draw(Char_t *option)
 195 {
 196 #ifdef use_root
 197   if(!fRootHisto)
 198     CreateRootHisto();
 199
 200   Double_t kappa,psi;
 201   Int_t content,bin;
 202   for(Int_t i=0; i<fNxbins; i++)
 203     {
 204       kappa = GetBinCenterX(i);
 205       for(Int_t j=0; j<fNybins; j++)
 206         {
 207           psi = GetBinCenterY(j);
 208           bin = GetBin(i,j);
 209           content = GetBinContent(bin);
 210           fRootHisto->Fill(kappa,psi,content);
 211         }
 212     }
 213   fRootHisto->Draw(option);
 214   return;
 215 #endif
 216   cerr<<"AliL3HistogramAdaptive::Draw : You need to compile with ROOT in order to draw histogram"<<endl;
 217 }
 218
 219 void AliL3HistogramAdaptive::Print()
 220 {
 221   cout<<"Printing content of histogram "<<fName<<endl;
 222   for(Int_t i=0; i<fNcells; i++)
 223     {
 224       if(GetBinContent(i)==0) continue;
 225       cout<<"Bin "<<i<<": "<<GetBinContent(i)<<endl;
 226     }
 227
 228 }