const Int_t NRowsSlice = 174;
const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+//const Int_t NRows[5][2] = {{ 0, 31},{46,173},{78,109},{110,141},{142,173}};
const Double_t Pi = 3.14159265358979323846;
const Double_t ToRad = Pi/180.;
}
-void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fake)
+void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove,TH2F *hist)
{
//Look at rawdata along the road specified by the track candidates.
continue;
}
- if(fake)
- fake->Fill(xyz[0],xyz[1],1);
Int_t npixs = 0;
//Get the pixels along the track candidate
if(fabs(xyz_pix[1] - xyz[1]) > num_of_pads_to_look*fTransform->GetPadPitchWidthLow()) continue;
npixs++;
+ if(remove) //Remove this pixel from image
+ pixel->fIndex = -1;
+
if(hist)
{
//fTransform->Local2Global(xyz_pix,slice);
virtual ~AliL3HoughEval();
TClonesArray *GetParticles(Char_t *rootfile);
- void LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist=0,TH2F *fake=0);
+ void LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove=(Bool_t)false,TH2F *hist=0);
void DisplaySlice(TH2F *hist);
void CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta);
#pragma link off all classes;
#pragma link off all functions;
+#pragma link C++ class AliL3Hough;
#pragma link C++ class AliL3HoughTransformer;
-#pragma link C++ class AliL3HoughPixel;
-//#pragma link C++ class AliL3HoughTrack;
#pragma link C++ class AliL3HoughMaxFinder;
#pragma link C++ class AliL3HoughEval;
#pragma link C++ class AliL3HoughMerge;
}
+AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
+{
+ Int_t xmin = hist->GetXaxis()->GetFirst();
+ Int_t xmax = hist->GetXaxis()->GetLast();
+ Int_t ymin = hist->GetYaxis()->GetFirst();
+ Int_t ymax = hist->GetYaxis()->GetLast();
+ Int_t bin[25],bin_index;
+ Stat_t value[25];
+
+ AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+ AliL3HoughTrack *track;
+
+ for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
+ {
+ for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
+ {
+ bin_index=0;
+ for(Int_t xb=xbin-2; xb<xbin+3; xb++)
+ {
+ for(Int_t yb=ybin-2; yb<ybin+3; yb++)
+ {
+ bin[bin_index]=hist->GetBin(xb,yb);
+ value[bin_index]=hist->GetBinContent(bin[bin_index]);
+ bin_index++;
+ }
+
+ }
+ if(value[12]==0) continue;
+ Int_t b=0;
+ while(1)
+ {
+ if(value[b] > value[12] || b==bin_index) break;
+ b++;
+ printf("b %d\n",b);
+ }
+ if(b == bin_index)
+ {
+ //Found maxima
+ Double_t max_x = hist->GetXaxis()->GetBinCenter(xbin);
+ Double_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
+ track = (AliL3HoughTrack*)tracks->NextTrack();
+ track->SetTrackParameters(max_x,max_y,value[12]);
+ }
+ }
+ }
+
+ tracks->QSort();
+ return tracks;
+}
+
AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row)
{
return tracks;
}
+
+AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *hist,Int_t nbins)
+{
+
+ AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+ AliL3HoughTrack *track;
+ Int_t xmin = hist->GetXaxis()->GetFirst();
+ Int_t xmax = hist->GetXaxis()->GetLast();
+ Int_t ymin = hist->GetYaxis()->GetFirst();
+ Int_t ymax = hist->GetYaxis()->GetLast();
+
+ Int_t weight_loc;
+ Int_t candidate=0;
+ for(Int_t xbin=xmin+nbins; xbin <= xmax - nbins; xbin++)
+ {
+ for(Int_t ybin=ymin+nbins; ybin <= ymax - nbins; ybin++)
+ {
+ weight_loc=0;
+ for(Int_t xbin_loc = xbin-nbins; xbin_loc <= xbin+nbins; xbin_loc++)
+ {
+ for(Int_t ybin_loc = ybin-nbins; ybin_loc <= ybin+nbins; ybin_loc++)
+ {
+ Int_t bin_loc = hist->GetBin(xbin_loc,ybin_loc);
+ weight_loc += (Int_t)hist->GetBinContent(bin_loc);
+ }
+ }
+
+ if(weight_loc > 0)
+ {
+ track = (AliL3HoughTrack*)tracks->NextTrack();
+ track->SetTrackParameters(hist->GetXaxis()->GetBinCenter(xbin),hist->GetYaxis()->GetBinCenter(ybin),weight_loc);
+ }
+
+ }
+ }
+ tracks->QSort();
+
+ AliL3HoughTrack *track1,*track2;
+
+ for(Int_t i=1; i<tracks->GetNTracks(); i++)
+ {
+ track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+ if(!track1) continue;
+ Int_t xbin1 = hist->GetXaxis()->FindBin(track1->GetKappa());
+ Int_t ybin1 = hist->GetYaxis()->FindBin(track1->GetPhi0());
+ for(Int_t j=0; j < i; j++)
+ {
+ track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
+ if(!track2) continue;
+ Int_t xbin2 = hist->GetXaxis()->FindBin(track2->GetKappa());
+ Int_t ybin2 = hist->GetYaxis()->FindBin(track2->GetPhi0());
+ if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
+ {
+ tracks->Remove(i);
+ break;
+ }
+ }
+
+ }
+ tracks->Compress();
+ return tracks;
+}
+
+AliL3TrackArray *AliL3HoughMaxFinder::LookInWindows(TH2F *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3)
+{
+
+ AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+ AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
+
+ Int_t xmin = hist->GetXaxis()->GetFirst();
+ Int_t xmax = hist->GetXaxis()->GetLast();
+ Int_t ymin = hist->GetYaxis()->GetFirst();
+ Int_t ymax = hist->GetYaxis()->GetLast();
+
+ for(Int_t xbin=xmin+nbins; xbin < xmax - nbins; xbin += nbins)
+ {
+ for(Int_t ybin=ymin+nbins; ybin<ymax-nbins; ybin += nbins)
+ {
+ //Int_t bin = hist->GetBin(xbin,ybin);
+ //if((Int_t)hist->GetBinContent(bin)==0) continue;
+
+ Int_t xrange[2] = {xbin-nbins,xbin+nbins};
+ Int_t yrange[2] = {ybin-nbins,ybin+nbins};
+ if(LocatePeak(hist,track,xrange,yrange,t1,t2,t3))
+ track = (AliL3HoughTrack*)tracks->NextTrack();
+ else
+ continue;
+ }
+ }
+ tracks->QSort();
+ return tracks;
+
+}
+
+Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3)
+{
+
+ /* Int_t xmin = hist->GetXaxis()->FindBin(track->GetKappa()) - nbins;
+ Int_t xmax = hist->GetXaxis()->FindBin(track->GetKappa()) + nbins;
+ Int_t ymin = hist->GetYaxis()->FindBin(track->GetPhi0()) - nbins;
+ Int_t ymax = hist->GetYaxis()->FindBin(track->GetPhi0()) + nbins;
+ Int_t nbinsx = nbins*2 + 1;
+
+ if(xmin < 0 || xmax > hist->GetXaxis()->GetNbins() || ymin < 0 || ymax > hist->GetYaxis()->GetNbins())
+ return false;
+ */
+ Int_t xmin = xrange[0];
+ Int_t xmax = xrange[1];
+ Int_t ymin = yrange[0];
+ Int_t ymax = yrange[1];
+
+ if(xmin < 0)
+ xmin=0;
+ if(xmax > hist->GetXaxis()->GetNbins())
+ xmax = hist->GetXaxis()->GetNbins();
+ if(ymin < 0)
+ ymin=0;
+ if(ymax > hist->GetYaxis()->GetNbins())
+ ymax = hist->GetYaxis()->GetNbins();
+
+ printf("Defined xrange %d %d yrange %d %d\n",xmin,xmax,ymin,ymax);
+
+
+ Int_t totbins = hist->GetXaxis()->GetNbins();
+ Int_t *m = new Int_t[totbins];
+ Int_t *m_low = new Int_t[totbins];
+ Int_t *m_up = new Int_t[totbins];
+
+
+ for(Int_t i=0; i<totbins; i++)
+ {
+ m[i]=0;
+ m_low[i]=0;
+ m_up[i]=0;
+ }
+
+ Int_t max_x=0,sum=0,max_xbin=0,bin;
+
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ for(Int_t ybin=ymin; ybin < ymax - t1; ybin++)
+ {
+ sum = 0;
+ for(Int_t y=ybin; y < ybin+t1; y++)
+ {
+ //Inside window
+ bin = hist->GetBin(xbin,y);
+ sum += (Int_t)hist->GetBinContent(bin);
+
+ }
+ if(sum > m[xbin]) //Max value locally in this xbin
+ {
+ m[xbin]=sum;
+ m_low[xbin]=ybin;
+ m_up[xbin]=ybin + t1 - 1;
+ }
+
+ }
+
+ if(m[xbin] > max_x) //Max value globally in x-direction
+ {
+ max_xbin = xbin;
+ max_x = m[xbin];//sum;
+ }
+ }
+
+ if(max_x == 0) return false;
+
+ //printf("max_xbin %d max_x %d m_low %d m_up %d\n",max_xbin,max_x,m_low[max_xbin],m_up[max_xbin]);
+ //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[max_xbin]),hist->GetYaxis()->GetBinCenter(m_up[max_xbin]));
+
+ //Determine a width in the x-direction
+ Int_t x_low=0,x_up=0;
+ for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
+ {
+ if(m[xbin]*t2 < max_x)
+ {
+ x_low = xbin+1;
+ break;
+ }
+ }
+ for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
+ {
+ if(m[xbin]*t2 < max_x)
+ {
+ x_up = xbin-1;
+ break;
+ }
+ }
+
+ Double_t top=0,butt=0,value,x_peak;
+
+ // printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(x_up));
+ //printf("Spread in x %d\n",x_up-x_low +1);
+
+ //Now, calculate the center of mass in x-direction
+ for(Int_t xbin=x_low; xbin <= x_up; xbin++)
+ {
+ value = hist->GetXaxis()->GetBinCenter(xbin);
+ top += value*m[xbin];
+ butt += m[xbin];
+ }
+ x_peak = top/butt;
+
+ //printf("FOund x_peak %f\n",x_peak);
+
+ //Find the peak in y direction:
+ Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
+ if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak)
+ x_l--;
+
+ Int_t x_u = x_l + 1;
+
+ if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak || hist->GetXaxis()->GetBinCenter(x_u) <= x_peak)
+ printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetXaxis()->GetBinCenter(x_l),x_peak,hist->GetXaxis()->GetBinCenter(x_u));
+
+ //printf("\nxlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_l),hist->GetXaxis()->GetBinCenter(x_u));
+
+ value=top=butt=0;
+
+ //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_l]),hist->GetYaxis()->GetBinCenter(m_up[x_l]));
+ //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_u]),hist->GetYaxis()->GetBinCenter(m_up[x_u]));
+
+ for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
+ {
+ value = hist->GetYaxis()->GetBinCenter(ybin);
+ bin = hist->GetBin(x_l,ybin);
+ top += value*hist->GetBinContent(bin);
+ butt += hist->GetBinContent(bin);
+ }
+ Double_t y_peak_low = top/butt;
+
+ //printf("y_peak_low %f\n",y_peak_low);
+
+ value=top=butt=0;
+ for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
+ {
+ value = hist->GetYaxis()->GetBinCenter(ybin);
+ bin = hist->GetBin(x_u,ybin);
+ top += value*hist->GetBinContent(bin);
+ butt += hist->GetBinContent(bin);
+ }
+ Double_t y_peak_up = top/butt;
+
+ //printf("y_peak_up %f\n",y_peak_up);
+
+ Double_t x_value_up = hist->GetXaxis()->GetBinCenter(x_u);
+ Double_t x_value_low = hist->GetXaxis()->GetBinCenter(x_l);
+
+ Double_t y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
+
+ //Find the weight:
+ bin = hist->FindBin(x_peak,y_peak);
+ Int_t weight = (Int_t)hist->GetBinContent(bin);
+
+ if(weight==0) return false;
+
+ track->SetTrackParameters(x_peak,y_peak,weight);
+
+ delete [] m;
+ delete [] m_up;
+ delete [] m_low;
+ return true;
+}
+
+
AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3)
{
//Attempt of a more sophisticated peak finder.
- AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+ AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack",1);
Int_t xmin = hist->GetXaxis()->GetFirst();
Int_t xmax = hist->GetXaxis()->GetLast();
Int_t ymin = hist->GetYaxis()->GetFirst();
Int_t ymax = hist->GetYaxis()->GetLast();
- Int_t nbinsx = hist->GetXaxis()->GetNbins();
+ Int_t nbinsx = hist->GetXaxis()->GetNbins()+1;
Int_t *m = new Int_t[nbinsx];
Int_t *m_low = new Int_t[nbinsx];
AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
track->SetTrackParameters(x_peak,y_peak,weight);
+
+ delete [] m;
+ delete [] m_low;
+ delete [] m_up;
+
return tracks;
- //delete [] m;
- //delete [] m_low;
- //delete [] m_up;
-
+
}
#include "AliL3RootTypes.h"
class AliL3TrackArray;
+class AliL3HoughTrack;
class AliL3HoughMaxFinder : public TObject {
AliL3HoughMaxFinder(Char_t *histotype);
virtual ~AliL3HoughMaxFinder();
+ AliL3TrackArray *FindBigMaxima(TH2F *hist);
AliL3TrackArray *FindMaxima(TH2F *hist,Int_t *rowrange=0,Int_t ref_row=0);
+ AliL3TrackArray *LookForPeaks(TH2F *hist,Int_t nbins);
+ AliL3TrackArray *LookInWindows(TH2F *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3);
+ Bool_t LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3);
AliL3TrackArray *FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3);
void SetThreshold(Int_t f) {fThreshold = f;}
}
}
-/*
-void AliL3HoughMerge::MergeLines()
-{
-
-
-
-}
-*/
fPatch = 0;
fRowContainer = 0;
fPhiRowContainer = 0;
+ fVolume=0;
fNumOfPadRows=0;
fContainerBounds=0;
fNDigits=0;
}
-AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
{
//Constructor
fTransform = new AliL3Transform();
-
fEtaMin = etarange[0];
fEtaMax = etarange[1];
fSlice = slice;
fPatch = patch;
- fNPhiSegments = phi_segments;
fNumOfPadRows=NRowsSlice;
}
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
+{
+
+ fTransform = new AliL3Transform();
+ if(etarange)
+ {
+ //assumes you want to look at a given etaslice only
+ fEtaMin = etarange[0];
+ fEtaMax = etarange[1];
+ }
+ else
+ {
+ //transform in all etaslices
+ fEtaMin = 0;
+ fEtaMax = slice < 18 ? 0.9 : -0.9;
+ }
+ fNumEtaSegments = n_eta_segments;
+ fSlice = slice;
+ fPatch = patch;
+ fNumOfPadRows = NRowsSlice;
+}
+
AliL3HoughTransformer::~AliL3HoughTransformer()
{
delete fTransform;
if(fPhiRowContainer)
delete [] fPhiRowContainer;
+ if(fVolume)
+ delete [] fVolume;
if(fIndex)
{
for(Int_t i=0; i<fNDigits; i++)
fIndex[i] = new Int_t[nbinsy+1];
Int_t sector,row;
+
for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- {
+ {
for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
{
Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
-
+ //printf("kappa %f phi0 %f\n",kappa,phi0);
Int_t bin = hist->FindBin(kappa,phi0);
if(fIndex[index][p]!=0)
printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
}
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
-{
- //Transformation is done with respect to local coordinates in slice.
-
-
- AliL3Digits *pix1;
- Int_t sector,row;
-
- //Define a common point
-
- /*
- Double_t rowdist1 = fTransform->Row2X(middle_row-1);
- Double_t rowdist2 = fTransform->Row2X(middle_row);
- Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
- */
-
- Double_t r_in_bundle = fTransform->Row2X(middle_row);
- //Make overlap between slices
- Double_t phi_min = -15*ToRad;
- Double_t phi_max = 15*ToRad;
-
- Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
-
- for(Int_t p=0; p <= fNPhiSegments; p++)
- {
- Double_t phi_in_bundle = phi_min + p*phi_slice;
- //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
-
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
- //for(Int_t padrow = 0; padrow <= middle_row; padrow++)
- {
-
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
- {
-
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
- //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-
- Double_t phi_pix = fTransform->GetPhi(xyz);
-
- //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
- Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
-
- Double_t phi0 = atan(tanPhi0);
- //if(padrow > middle_row)
- //phi0 = -phi0;
- //if(phi0 < 0.55 || phi0 > 0.88) continue;
-
- Double_t angle = phi_pix - phi0;
- Double_t kappa = 2*sin(angle)/r_pix;
-
- //Double_t angle = phi_in_bundle - phi0;
- //Double_t kappa = 2*sin(angle)/r_in_bundle;
-
- //if(kappa < -0.006 || kappa > 0.006) continue;
-
- Short_t signal = pix1->fCharge;
-
- hist->Fill(kappa,phi0,signal);
-
-
- }
-
- }
- }
-
-}
-
-
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
+void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
{
//Transformation is done with respect to local coordinates in slice.
//Transform every pixel into whole phirange, using parametrisation:
Int_t nbinsy = hist->GetNbinsY();
+ Int_t totsignal=0,vol_index;
+ if(fNumEtaSegments==1)
+ eta_index=0; //only looking in one etaslice.
+
for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
{
-
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+ vol_index = eta_index*fNumOfPadRows + padrow;
+ //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+ for(pix1=(AliL3Digits*)fVolume[vol_index].first; pix1!=0; pix1=(AliL3Digits*)pix1->fNextVolumePixel)
{
Short_t signal = pix1->fCharge;
Int_t index = pix1->fIndex;
-
+ if(index < 0) continue; //This pixel has been removed.
+ totsignal += signal;
for(Int_t p=0; p <= nbinsy; p++)
hist->AddBinContent(fIndex[index][p],signal);
}
-
+ printf("Total signal %d\n",totsignal);
}
/*
printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
return;
}
- //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
- for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
+ for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+ //for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
{
//printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
Float_t xyz[3];
Int_t digit_counter=0;
Float_t xyz[3];
- Double_t eta,phi;
+ Double_t eta;
Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
- fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
- printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
- if(fPhiRowContainer)
- delete [] fPhiRowContainer;
- fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
- memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
+ //fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
+ //printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
- Double_t phi_min = -10*ToRad;
- Double_t phi_max = 10*ToRad;
- Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+ /*
+ if(fPhiRowContainer)
+ delete [] fPhiRowContainer;
+ fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
+ memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
+ */
+ fContainerBounds = (fNumEtaSegments+1)*(fNumOfPadRows+1);
+ if(fVolume)
+ delete [] fVolume;
+ fVolume = new AliL3HoughContainer[fContainerBounds];
+ memset(fVolume,0,fContainerBounds*sizeof(AliL3HoughContainer));
+ /*
+ Double_t phi_min = -10*ToRad;
+ Double_t phi_max = 10*ToRad;
+ Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+ */
+ Double_t eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
Int_t index;
digit_counter=0;
+ printf("\nLoading ALL pixels in slice\n\n");
+
for (Int_t i=0; i<num_of_entries; i++)
{
t->GetEvent(i);
fTransform->Raw2Global(xyz,sector,row,pad,time);
eta = fTransform->GetEta(xyz);
+
if(eta < fEtaMin || eta > fEtaMax) continue;
fTransform->Global2Local(xyz,sector);
- phi = fTransform->GetPhi(xyz);
+ //phi = fTransform->GetPhi(xyz);
if(hist)
hist->Fill(xyz[0],xyz[1],signal);
((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
fRowContainer[padrow].last = (void*)dig;
- Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
- index = phi_index*fNumOfPadRows + padrow;
- if(phi_index > fContainerBounds || phi_index < 0)
+ //thisHit->etaIndex=(Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1);
+ Int_t eta_index = (Int_t)((eta-fEtaMin)/eta_slice);
+ index = eta_index*fNumOfPadRows + padrow;
+ if(index > fContainerBounds || index < 0)
{
- printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
+
+ //printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
continue;
}
-
- if(fPhiRowContainer[index].first == NULL)
- fPhiRowContainer[index].first = (void*)dig;
+
+ if(fVolume[index].first == NULL)
+ fVolume[index].first = (void*)dig;
else
- ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
- fPhiRowContainer[index].last=(void*)dig;
+ ((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
+ fVolume[index].last = (void*)dig;
+ /*
+ Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
+ index = phi_index*fNumOfPadRows + padrow;
+ if(phi_index > fContainerBounds || phi_index < 0)
+ {
+ printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
+ continue;
+ }
+
+ if(fPhiRowContainer[index].first == NULL)
+ fPhiRowContainer[index].first = (void*)dig;
+ else
+ ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
+ fPhiRowContainer[index].last=(void*)dig;
+ */
}while (digarr->Next());
}
UChar_t fPad;
UShort_t fTime;
Int_t fIndex;
- AliL3Digits *nextPhiRowPixel;
+ AliL3Digits *fNextVolumePixel;
+ //AliL3Digits *nextPhiRowPixel;
AliL3Digits *nextRowPixel;
};
AliL3HoughContainer *fRowContainer; //!
AliL3HoughContainer *fPhiRowContainer; //!
+ AliL3HoughContainer *fVolume; //!
Int_t fContainerBounds;
Int_t fNDigits;
Int_t **fIndex; //!
public:
AliL3HoughTransformer();
- AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments);
+ AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange);
+ AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange=0,Int_t n_eta_segments=1);
virtual ~AliL3HoughTransformer();
- void Transform2Circle(TH2F *hist,Int_t middle_row);
- void Transform2Circle(TH2F *hist);
+ void Transform2Circle(TH2F *hist,Int_t eta_index);
void Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw=0);
void TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks);
void GetPixels(Char_t *rootfile,TH2F *hist=0);
# C++ sources
-SRCS = AliL3HoughTransformer.cxx AliL3HoughPixel.cxx \
+SRCS = AliL3HoughTransformer.cxx AliL3Hough.cxx \
AliL3HoughMaxFinder.cxx AliL3HoughEval.cxx AliL3HoughMerge.cxx
{
gStyle->SetOptStat(0);
- Int_t xbin = 60;
- Int_t ybin = 60;
+ Int_t xbin = 70;
+ Int_t ybin = 70;
Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006
//Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
Int_t xr[2] = {0,250};
Int_t yr[2] = {-125,125};
- TH1F *ntracks = new TH1F("ntracks","",100,0,200);
- TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
- TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
- TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125);
- TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+ ntracks = new TH1F("ntracks","",100,0,200);
+ raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+ road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+ fake = new TH2F("fake","",300,0,300,250,-125,125);
+ peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
peaks->SetMarkerStyle(3);
peaks->SetMarkerColor(2);
+ road->SetMarkerStyle(6);
+
+ int slice = 2;
+ //float eta[2] = {0.3,0.4};
+ float eta[2] = {0.04,0.05};
- int slice = 2,patch=0,n_phi_segments=40;
- float eta[2] = {0.3,0.4};
- // float eta[2] = {0.04,0.05};
-
- AliL3HoughTransformer *a;
- AliL3HoughMaxFinder *b;
- AliL3HoughEval *c = new AliL3HoughEval(slice);
- AliL3HoughMerge *d = new AliL3HoughMerge(slice);
- TH2F *hist;
+ b = new AliL3HoughMaxFinder("KappaPhi");
+ //hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+ TH2F **histos = new TH2F*[5];
for(int pat=0; pat<5; pat++)
{
- a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
- hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+ a = new AliL3HoughTransformer(slice,pat,eta);
+ histos[pat] = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
a->GetPixels(rootfile,raw);
- a->Transform2Circle(hist,87);
-
- b = new AliL3HoughMaxFinder("KappaPhi");
- AliL3TrackArray *tracks = b->FindMaxima(hist);
- c->SetTransformer(a);
- c->LookInsideRoad(tracks);
- d->FillTracks(tracks,pat);
+ a->InitTemplates(hist);
+ a->Transform2Circle(hist);
delete a;
- delete b;
}
- return;
-
-
- TH2F *merge_hist = new TH2F("merge_hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
-
- d->FillHisto(merge_hist);
-
- b = new AliL3HoughMaxFinder(slice,pat);
- b->FindMaxima(merge_hist);
-
- AliL3TrackArray *mtrs = b->GetTracks();
+ hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+ for(Int_t i=0; i<5; i++)
+ {
+ hist->Add(histos[i],1);
+ }
+ tracks = (AliL3TrackArray*)b->FindPeak(hist,4,0.95,5);
+ track = (AliL3HoughTrack*)tracks->GetCheckedTrack(0);
+ peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
- c->CompareMC(rootfile,mtrs,eta);
- AliL3Transform *transform = new AliL3Transform();
- for(int tr=0; tr<mtrs->GetNTracks(); tr++)
+ for(Int_t padrow=0; padrow<174; padrow++)
{
- AliL3HoughTrack *t = (AliL3HoughTrack*)mtrs->GetCheckedTrack(tr);
- if(!t) {printf("NO TRACK11\n"); continue;}
- if(t->GetMCid() < 0) {printf("Track %d was fake, weigth %d\n",tr,t->GetNHits()); continue;}
- printf("Found pt %f weigth %d\n",t->GetPt(),t->GetNHits());
- //printf("charge %f\n",t->GetCharge());
- for(Int_t j=0; j<174; j++)
- {
- Float_t xyz[3];
- if(!t->GetCrossingPoint(j,xyz))
- {
- printf("Track does not cross line\n");
- continue;
- }
- road->Fill(xyz[0],xyz[1],1);
-
- }
-
-
+ Float_t xyz[3];
+ track->GetCrossingPoint(padrow,xyz);
+ road->Fill(xyz[0],xyz[1],1);
}
-
- TCanvas *c1 = new TCanvas("c1","",800,800);
- SetCanvasOptions(c1);
+ c1 = new TCanvas("c1","",1000,1000);
c1->Divide(2,2);
c1->cd(1);
- merge_hist->Draw("box");
-
- //TCanvas *c2 = new TCanvas("c2","",2);
- c1->cd(3);
+ hist->Draw("box");
+ peaks->Draw("same");
+ c1->cd(2);
raw->Draw();
-
- //TCanvas *c3 = new TCanvas("c3","",2);
- c1->cd(4);
- road->SetMarkerStyle(7);
+ c1->cd(3);
road->Draw();
+ printf("Pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
+ return;
+
delete b;
delete d;
TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
peaks->SetMarkerStyle(3);
peaks->SetMarkerColor(2);
+ road->SetMarkerStyle(6);
+ road->SetMarkerColor(2);
real_peaks = new TH2F("real_peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
real_peaks->SetMarkerStyle(3);
real_peaks->SetMarkerColor(4);
int slice = 2,n_phi_segments=30;
- float eta[2] = {0.3,0.4};
- //float eta[2] = {0.03,0.04};
+ //double eta[2] = {0.3,0.4};
+ double eta[2] = {0.04,0.05};
- a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
+ a = new AliL3HoughTransformer(slice,0,eta,1);
a->GetPixels(rootfile,raw);
a->InitTemplates(hist);
- a->Transform2Circle(hist);
-
+ a->Transform2Circle(hist,3);
+
b = new AliL3HoughMaxFinder("KappaPhi");
- tracks = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
-
c = new AliL3HoughEval(a);
- c->LookInsideRoad(tracks);
-
+ // tracks = (AliL3TrackArray*)b->FindMaxima(hist);
+ /*
+ int n_iter=1;
+ AliL3TrackArray **track_cand = new AliL3TrackArray*[n_iter];
+ for(int k=0; k<n_iter; k++)
+ {
+ hist->Reset();
+ a->Transform2Circle(hist);
+ track_cand[k] = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
+ c->LookInsideRoad(track_cand[k],true);
+ }
+
+ //tracks = (AliL3TrackArray*)b->LookInWindows(hist,4,4,0.95,5);
+
+ //tracks = (AliL3TrackArray*)b->FindMaxima(hist);
+ //cout << "Found "<<tracks->GetNTracks()<<" peaks"<<endl;
+
+ for(int j=0; j<n_iter; j++)
+ {
+ track = (AliL3HoughTrack*)track_cand[j]->GetCheckedTrack(0);
+ if(!track) continue;
+ printf("Pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
+ for(int padrow=0; padrow<174; padrow++)
+ {
+ Float_t xyz[3];
+ track->GetCrossingPoint(padrow,xyz);
+ road->Fill(xyz[0],xyz[1],1);
+ }
+ peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
+ }
+ */
+
+ printf("Looking for big maxima\n");
+ tracks = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
+ cout<<"NTracks after checking "<<tracks->GetNTracks()<<endl;
for(int i=0; i<tracks->GetNTracks(); i++)
{
track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
printf("Pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
}
- c2 = new TCanvas("c2","",800,800);
- // c2->Divide(2,2);
- //c2->cd(1);
- hist->Draw("box");
- peaks->Draw("same");
- real_peaks->Draw("same");
-
- /*
- AliL3TrackArray *tracks = b->FindMaxima(hist);
-
- track = (AliL3HoughTrack*)tracks->GetCheckedTrack(0);
- peaks->Fill(track->GetKappa(),track->GetPhi0(),track->GetWeight());
- int xbin = hist->GetXaxis()->FindBin(track->GetKappa());
- int ybin = hist->GetYaxis()->FindBin(track->GetPhi0());
-
- TH1D *proj_kappa = hist->ProjectionX("proj_kappa",ybin,ybin);
- TH1D *proj_phi = hist->ProjectionY("proj_phi",xbin,xbin);
-
- TCanvas *c1 = new TCanvas("c1","",800,800);
- c1->Divide(2,2);
- c1->cd(1);
+ c2 = new TCanvas("c2","",900,900);
+ c2->Divide(2,2);
+ c2->cd(1);
+ raw->Draw("");
+ road->Draw("same");
+ c2->cd(2);
+ raw->DrawCopy();
+ c2->cd(3);
hist->Draw("box");
peaks->Draw("same");
-
- c1->cd(3);
- proj_kappa->Draw();
- c1->cd(4);
- proj_phi->Draw();
- return;
- */
+ //real_peaks->Draw("same");
+ c2->cd(4);
+ hist->DrawCopy("lego1");
+
delete a;
delete b;
delete c;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff