[u/mrichter/AliRoot.git] / HLT / hough / AliL3HistogramAdaptive.cxx

//$Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ASV

#include <string.h>
#include "AliL3Logging.h"
#include "AliL3HistogramAdaptive.h"
#include "AliL3Transform.h"
#include "AliL3Track.h"

//_____________________________________________________________
// AliL3HistogramAdaptive
//
// 2D histogram class

ClassImp(AliL3HistogramAdaptive)

AliL3HistogramAdaptive::AliL3HistogramAdaptive() : AliL3Histogram()
{
  fKvalue=0;
  fPtstep=0;
}

  
AliL3HistogramAdaptive::AliL3HistogramAdaptive(Char_t *name,Int_t firstpatch,Int_t lastpatch,Double_t minpt,Double_t maxpt,
					       Int_t nybins,Double_t ymin,Double_t ymax)
{
  strcpy(fName,name);
  
  fXmin = -1*AliL3Transform::GetBFact()*AliL3Transform::GetBField()/minpt;
  fXmax = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/minpt;

  fPtstep = 0.1;
  fMinPt = minpt;
  fMaxPt = maxpt;
  fNxbins = InitPtBins(firstpatch,lastpatch,0.1);
  
  fNybins = nybins;
  fYmin = ymin;
  fYmax = ymax;
  fFirstXbin=1;
  fFirstYbin=1;
  fLastXbin = fNxbins;
  fLastYbin = fNybins;
  fNcells = (fNxbins+2)*(fNybins+2);
  
  fThreshold=0;
  fContent = new Int_t[fNcells];
  Reset();
}

AliL3HistogramAdaptive::~AliL3HistogramAdaptive()
{
  if(fKvalue)
    delete [] fKvalue;
}

Int_t AliL3HistogramAdaptive::InitPtBins(Int_t firstpatch,Int_t lastpatch,Double_t xyresolution)
{
  
  Double_t pt = fMinPt;
  AliL3Track *track=0;
  Float_t xyz1[3],xyz2[3];
  Double_t angle1,angle2,diff_angle,s_tot,local_pt,pt_nplus1;
  Int_t nrows = AliL3Transform::GetLastRow(lastpatch)-AliL3Transform::GetFirstRow(firstpatch)+1;
  Int_t nptbins=0,index=-1;
  Int_t bounds = (Int_t)rint((fMaxPt-fMinPt)/fPtstep)+1;
  fKvalue = new Double_t[bounds];
  
  while(pt < fMaxPt + fPtstep)
    {
      track = new AliL3Track();
      track->SetPt(pt);
      track->SetPsi(0);
      track->SetCharge(1);
      track->SetFirstPoint(0,0,0);
      track->CalculateHelix();
      
      track->GetCrossingPoint(AliL3Transform::GetFirstRow(firstpatch),xyz1);
      track->GetCrossingPoint(AliL3Transform::GetLastRow(lastpatch),xyz2);
      
      angle1 = atan2((xyz2[1] - track->GetCenterY()),(xyz2[0] - track->GetCenterX()));
      if(angle1 < 0) angle1 += 2.*AliL3Transform::Pi();
      angle2 = atan2((xyz1[1] - track->GetCenterY()),(xyz1[0] - track->GetCenterX()));
      if(angle2 < 0) angle2 += 2.*AliL3Transform::Pi();
      diff_angle = angle1 - angle2;
      diff_angle = fmod(diff_angle,2*AliL3Transform::Pi());
      if((track->GetCharge()*diff_angle) > 0) diff_angle = diff_angle - track->GetCharge()*2.*AliL3Transform::Pi();
      
      s_tot = fabs(diff_angle)*track->GetRadius();
      
      index = (Int_t)rint((pt-fMinPt)/fPtstep);
      if(index < 0 || index > bounds)
	{
	  cerr<<"AliL3HistogramAdaptive::InitPtBins : Index out of range "<<index<<endl;
	  return 0;
	}
      //cout<<"Filling index "<<index<<" pt "<<pt<<" fMinPt "<<fMinPt<<" fPtstep "<<fPtstep<<endl;
      fKvalue[index] = xyresolution*sqrt(720/(nrows+4))/(AliL3Transform::GetBFact()*AliL3Transform::GetBField()*s_tot*s_tot);

      local_pt = pt;
      while(local_pt < pt + fPtstep)
	{
	  if(local_pt > fMaxPt) break;
	  pt_nplus1 = fKvalue[index]*local_pt*local_pt + local_pt;
	  //  cout<<"Pt "<<local_pt<<" binsize "<<(pt_nplus1-local_pt)<<endl;
	  local_pt = pt_nplus1;
	  nptbins++;
	}
      delete track;
      pt += fPtstep;
    }
  cout<<"Number of ptbins "<<nptbins*2<<endl;
  
  return nptbins*2; //Both negative and positive kappa.
}


void AliL3HistogramAdaptive::Fill(Double_t x,Double_t y,Int_t weight)
{
  Int_t bin = FindBin(x,y);
  if(bin < 0)
    return;
  AddBinContent(bin,weight);

}

Int_t AliL3HistogramAdaptive::FindBin(Double_t x,Double_t y)
{
  
  Int_t xbin = FindXbin(x);
  Int_t ybin = FindYbin(y);
  
  if(xbin < 0) 
    return -1;
  return GetBin(xbin,ybin);
}

Int_t AliL3HistogramAdaptive::FindXbin(Double_t x)
{
  if(x < fXmin || x > fXmax)
    return 0;
  
  Double_t pt = fabs(AliL3Transform::GetBFact()*AliL3Transform::GetBField()/x);
  
  if(pt < fMinPt || pt > fMaxPt) 
    return -1;
  Double_t local_pt = fMinPt;
  Double_t pt_nplus1;
  Int_t index,bin=0;
  while(local_pt < fMaxPt)
    {
      index = (Int_t)rint((local_pt-fMinPt)/fPtstep);
      pt_nplus1 = fKvalue[index]*local_pt*local_pt + local_pt;
      //cout<<"Looking for pt "<<local_pt<<" pt_nplus1 "<<pt_nplus1<<" index "<<index<<endl;
      if(pt >= local_pt && pt < pt_nplus1)
	break;
      local_pt = pt_nplus1;
      bin++;
    }
  if(x < 0)
    return bin;
  else
    return fNxbins - bin;
}

Int_t AliL3HistogramAdaptive::FindYbin(Double_t y)
{
  if(y < fYmin || y > fYmax)
    return 0;
  
  return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));

}

Double_t AliL3HistogramAdaptive::GetBinCenterX(Int_t xbin)
{
  Double_t local_pt = fMinPt;
  Double_t pt_nplus1=0;
  Int_t index,bin=0;
  
  while(local_pt < fMaxPt)
    {
      index = (Int_t)rint((local_pt-fMinPt)/fPtstep);
      pt_nplus1 = fKvalue[index]*local_pt*local_pt + local_pt;
      if(xbin == bin || xbin == fNxbins - bin)
	break;
      local_pt = pt_nplus1;
      bin++;
    }
  
  Double_t binwidth = (pt_nplus1 - local_pt);
  Double_t kappa = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/(local_pt + binwidth*0.5);
  
  cout<<"Returning xbin "<<xbin<<" out of total "<<fNxbins<<endl;
  if(xbin < fNxbins/2)
    return -1.*kappa;
  else
    return kappa;
      
  //Double_t binwidth = (fXmax - fXmin) / fNxbins;
  //return fXmin + (xbin-1) * binwidth + 0.5*binwidth;
}

Double_t AliL3HistogramAdaptive::GetBinCenterY(Int_t ybin)
{
  
  Double_t binwidth = (fYmax - fYmin) / fNybins;
  return fYmin + (ybin-1) * binwidth + 0.5*binwidth;
  
}

void AliL3HistogramAdaptive::Draw(Char_t *option)
{
#ifdef use_root
  if(!fRootHisto)
    CreateRootHisto();
  
  Double_t kappa,psi;
  Int_t content,bin;
  for(Int_t i=0; i<fNxbins; i++)
    {
      kappa = GetBinCenterX(i);
      for(Int_t j=0; j<fNybins; j++)
	{
	  psi = GetBinCenterY(j);
	  bin = GetBin(i,j);
	  content = GetBinContent(bin);
	  fRootHisto->Fill(kappa,psi,content);
	}
    }
  fRootHisto->Draw(option);
  return;
#endif
  cerr<<"AliL3HistogramAdaptive::Draw : You need to compile with ROOT in order to draw histogram"<<endl;
}

void AliL3HistogramAdaptive::Print()
{
  cout<<"Printing content of histogram "<<fName<<endl;
  for(Int_t i=0; i<fNcells; i++)
    {
      if(GetBinContent(i)==0) continue;
      cout<<"Bin "<<i<<": "<<GetBinContent(i)<<endl;
    }

}
Commit	Line	Data
6b88e8ee	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	fKvalue=0;
	22	fPtstep=0;
	23	}
	24
	25
	26	AliL3HistogramAdaptive::AliL3HistogramAdaptive(Char_t *name,Int_t firstpatch,Int_t lastpatch,Double_t minpt,Double_t maxpt,
	27	Int_t nybins,Double_t ymin,Double_t ymax)
	28	{
	29	strcpy(fName,name);
	30
	31	fXmin = -1AliL3Transform::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(firstpatch,lastpatch,0.1);
	38
	39	fNybins = nybins;
	40	fYmin = ymin;
	41	fYmax = ymax;
	42	fFirstXbin=1;
	43	fFirstYbin=1;
	44	fLastXbin = fNxbins;
	45	fLastYbin = fNybins;
	46	fNcells = (fNxbins+2)*(fNybins+2);
	47
	48	fThreshold=0;
	49	fContent = new Int_t[fNcells];
	50	Reset();
	51	}
	52
	53	AliL3HistogramAdaptive::~AliL3HistogramAdaptive()
	54	{
	55	if(fKvalue)
	56	delete [] fKvalue;
	57	}
	58
	59	Int_t AliL3HistogramAdaptive::InitPtBins(Int_t firstpatch,Int_t lastpatch,Double_t xyresolution)
	60	{
	61
	62	Double_t pt = fMinPt;
	63	AliL3Track *track=0;
	64	Float_t xyz1[3],xyz2[3];
65	Double_t angle1,angle2,diff_angle,s_tot,local_pt,pt_nplus1;
66	Int_t nrows = AliL3Transform::GetLastRow(lastpatch)-AliL3Transform::GetFirstRow(firstpatch)+1;
67	Int_t nptbins=0,index=-1;
68	Int_t bounds = (Int_t)rint((fMaxPt-fMinPt)/fPtstep)+1;
69	fKvalue = new Double_t[bounds];
70
71	while(pt < fMaxPt + fPtstep)
72	{
73	track = new AliL3Track();
74	track->SetPt(pt);
75	track->SetPsi(0);
76	track->SetCharge(1);
77	track->SetFirstPoint(0,0,0);
78	track->CalculateHelix();
79
80	track->GetCrossingPoint(AliL3Transform::GetFirstRow(firstpatch),xyz1);
81	track->GetCrossingPoint(AliL3Transform::GetLastRow(lastpatch),xyz2);
82
83	angle1 = atan2((xyz2[1] - track->GetCenterY()),(xyz2[0] - track->GetCenterX()));
84	if(angle1 < 0) angle1 += 2.*AliL3Transform::Pi();
85	angle2 = atan2((xyz1[1] - track->GetCenterY()),(xyz1[0] - track->GetCenterX()));
86	if(angle2 < 0) angle2 += 2.*AliL3Transform::Pi();
87	diff_angle = angle1 - angle2;
88	diff_angle = fmod(diff_angle,2*AliL3Transform::Pi());
89	if((track->GetCharge()diff_angle) > 0) diff_angle = diff_angle - track->GetCharge()2.*AliL3Transform::Pi();
90
91	s_tot = fabs(diff_angle)*track->GetRadius();
92
93	index = (Int_t)rint((pt-fMinPt)/fPtstep);
94	if(index < 0 \|\| index > bounds)
95	{
96	cerr<<"AliL3HistogramAdaptive::InitPtBins : Index out of range "<<index<<endl;
97	return 0;
98	}
99	//cout<<"Filling index "<<index<<" pt "<<pt<<" fMinPt "<<fMinPt<<" fPtstep "<<fPtstep<<endl;
100	fKvalue[index] = xyresolutionsqrt(720/(nrows+4))/(AliL3Transform::GetBFact()AliL3Transform::GetBField()s_tots_tot);
101
102	local_pt = pt;
103	while(local_pt < pt + fPtstep)
104	{
105	if(local_pt > fMaxPt) break;
106	pt_nplus1 = fKvalue[index]local_ptlocal_pt + local_pt;
107	// cout<<"Pt "<<local_pt<<" binsize "<<(pt_nplus1-local_pt)<<endl;
108	local_pt = pt_nplus1;
109	nptbins++;
110	}
111	delete track;
112	pt += fPtstep;
113	}
114	cout<<"Number of ptbins "<<nptbins*2<<endl;
115
116	return nptbins*2; //Both negative and positive kappa.
117	}
118
119
120	void AliL3HistogramAdaptive::Fill(Double_t x,Double_t y,Int_t weight)
121	{
122	Int_t bin = FindBin(x,y);
123	if(bin < 0)
124	return;
125	AddBinContent(bin,weight);
126
127	}
128
129	Int_t AliL3HistogramAdaptive::FindBin(Double_t x,Double_t y)
130	{
131
132	Int_t xbin = FindXbin(x);
133	Int_t ybin = FindYbin(y);
134
135	if(xbin < 0)
136	return -1;
137	return GetBin(xbin,ybin);
138	}
139
140	Int_t AliL3HistogramAdaptive::FindXbin(Double_t x)
141	{
142	if(x < fXmin \|\| x > fXmax)
143	return 0;
144
145	Double_t pt = fabs(AliL3Transform::GetBFact()*AliL3Transform::GetBField()/x);
146
147	if(pt < fMinPt \|\| pt > fMaxPt)
148	return -1;
149	Double_t local_pt = fMinPt;
150	Double_t pt_nplus1;
151	Int_t index,bin=0;
152	while(local_pt < fMaxPt)
153	{
154	index = (Int_t)rint((local_pt-fMinPt)/fPtstep);
155	pt_nplus1 = fKvalue[index]local_ptlocal_pt + local_pt;
156	//cout<<"Looking for pt "<<local_pt<<" pt_nplus1 "<<pt_nplus1<<" index "<<index<<endl;
157	if(pt >= local_pt && pt < pt_nplus1)
158	break;
159	local_pt = pt_nplus1;
160	bin++;
161	}
162	if(x < 0)
163	return bin;
164	else
165	return fNxbins - bin;
166	}
167
168	Int_t AliL3HistogramAdaptive::FindYbin(Double_t y)
169	{
170	if(y < fYmin \|\| y > fYmax)
171	return 0;
172
173	return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
174
175	}
176
177	Double_t AliL3HistogramAdaptive::GetBinCenterX(Int_t xbin)
178	{
179	Double_t local_pt = fMinPt;
180	Double_t pt_nplus1=0;
181	Int_t index,bin=0;
182
183	while(local_pt < fMaxPt)
184	{
185	index = (Int_t)rint((local_pt-fMinPt)/fPtstep);
186	pt_nplus1 = fKvalue[index]local_ptlocal_pt + local_pt;
187	if(xbin == bin \|\| xbin == fNxbins - bin)
188	break;
189	local_pt = pt_nplus1;
190	bin++;
191	}
192
193	Double_t binwidth = (pt_nplus1 - local_pt);
194	Double_t kappa = AliL3Transform::GetBFact()AliL3Transform::GetBField()/(local_pt + binwidth0.5);
195
196	cout<<"Returning xbin "<<xbin<<" out of total "<<fNxbins<<endl;
197	if(xbin < fNxbins/2)
198	return -1.*kappa;
199	else
200	return kappa;
201
202	//Double_t binwidth = (fXmax - fXmin) / fNxbins;
203	//return fXmin + (xbin-1) * binwidth + 0.5*binwidth;
204	}
205
206	Double_t AliL3HistogramAdaptive::GetBinCenterY(Int_t ybin)
207	{
208
209	Double_t binwidth = (fYmax - fYmin) / fNybins;
210	return fYmin + (ybin-1) * binwidth + 0.5*binwidth;
211
212	}
213
214	void AliL3HistogramAdaptive::Draw(Char_t *option)
215	{
216	#ifdef use_root
217	if(!fRootHisto)
218	CreateRootHisto();
219
220	Double_t kappa,psi;
221	Int_t content,bin;
222	for(Int_t i=0; i<fNxbins; i++)
223	{
224	kappa = GetBinCenterX(i);
225	for(Int_t j=0; j<fNybins; j++)
226	{
227	psi = GetBinCenterY(j);
228	bin = GetBin(i,j);
229	content = GetBinContent(bin);
230	fRootHisto->Fill(kappa,psi,content);
231	}
232	}
233	fRootHisto->Draw(option);
234	return;
235	#endif
236	cerr<<"AliL3HistogramAdaptive::Draw : You need to compile with ROOT in order to draw histogram"<<endl;
237	}
238
239	void AliL3HistogramAdaptive::Print()
240	{
241	cout<<"Printing content of histogram "<<fName<<endl;
242	for(Int_t i=0; i<fNcells; i++)
243	{
244	if(GetBinContent(i)==0) continue;
245	cout<<"Bin "<<i<<": "<<GetBinContent(i)<<endl;
246	}
247
248	}