From 52a2a604748eab2530f341f89e183d0eca358c1e Mon Sep 17 00:00:00 2001 From: vestbo Date: Fri, 30 Mar 2001 14:55:23 +0000 Subject: [PATCH] Updating changes before the weekend --- HLT/hough/AliL3Defs.h | 1 + HLT/hough/AliL3HoughEval.cxx | 7 +- HLT/hough/AliL3HoughEval.h | 2 +- HLT/hough/AliL3HoughLinkDef.h | 3 +- HLT/hough/AliL3HoughMaxFinder.cxx | 330 +++++++++++++++++++++++++++- HLT/hough/AliL3HoughMaxFinder.h | 5 + HLT/hough/AliL3HoughMerge.cxx | 8 - HLT/hough/AliL3HoughTransformer.cxx | 198 ++++++++--------- HLT/hough/AliL3HoughTransformer.h | 10 +- HLT/hough/Makefile | 2 +- HLT/hough/hough_merge.C | 108 ++++----- HLT/hough/testPF.C | 91 ++++---- 12 files changed, 528 insertions(+), 237 deletions(-) diff --git a/HLT/hough/AliL3Defs.h b/HLT/hough/AliL3Defs.h index c9bfdf052e5..a59c4c5525c 100644 --- a/HLT/hough/AliL3Defs.h +++ b/HLT/hough/AliL3Defs.h @@ -5,6 +5,7 @@ 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.; diff --git a/HLT/hough/AliL3HoughEval.cxx b/HLT/hough/AliL3HoughEval.cxx index b7c608f2522..5bdfcc92d34 100644 --- a/HLT/hough/AliL3HoughEval.cxx +++ b/HLT/hough/AliL3HoughEval.cxx @@ -40,7 +40,7 @@ AliL3HoughEval::~AliL3HoughEval() } -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. @@ -75,8 +75,6 @@ void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fak continue; } - if(fake) - fake->Fill(xyz[0],xyz[1],1); Int_t npixs = 0; //Get the pixels along the track candidate @@ -93,6 +91,9 @@ void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fak 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); diff --git a/HLT/hough/AliL3HoughEval.h b/HLT/hough/AliL3HoughEval.h index 1d64d91f74c..eb56a8acefc 100644 --- a/HLT/hough/AliL3HoughEval.h +++ b/HLT/hough/AliL3HoughEval.h @@ -25,7 +25,7 @@ class AliL3HoughEval : public TObject { 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); diff --git a/HLT/hough/AliL3HoughLinkDef.h b/HLT/hough/AliL3HoughLinkDef.h index 534862b547f..18a9a234a4e 100644 --- a/HLT/hough/AliL3HoughLinkDef.h +++ b/HLT/hough/AliL3HoughLinkDef.h @@ -4,9 +4,8 @@ #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; diff --git a/HLT/hough/AliL3HoughMaxFinder.cxx b/HLT/hough/AliL3HoughMaxFinder.cxx index 9a502f0a03b..ee8473a4e3b 100644 --- a/HLT/hough/AliL3HoughMaxFinder.cxx +++ b/HLT/hough/AliL3HoughMaxFinder.cxx @@ -35,6 +35,56 @@ AliL3HoughMaxFinder::~AliL3HoughMaxFinder() } +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; xbinGetBin(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) { @@ -112,17 +162,283 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_ 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; iGetNTracks(); 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; ybinGetBin(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; iGetBin(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]; @@ -270,11 +586,13 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I 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; - + } diff --git a/HLT/hough/AliL3HoughMaxFinder.h b/HLT/hough/AliL3HoughMaxFinder.h index 4eca14eb5d9..d7788136aeb 100644 --- a/HLT/hough/AliL3HoughMaxFinder.h +++ b/HLT/hough/AliL3HoughMaxFinder.h @@ -4,6 +4,7 @@ #include "AliL3RootTypes.h" class AliL3TrackArray; +class AliL3HoughTrack; class AliL3HoughMaxFinder : public TObject { @@ -19,7 +20,11 @@ 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;} diff --git a/HLT/hough/AliL3HoughMerge.cxx b/HLT/hough/AliL3HoughMerge.cxx index 534e329c8c3..81604dae3fc 100644 --- a/HLT/hough/AliL3HoughMerge.cxx +++ b/HLT/hough/AliL3HoughMerge.cxx @@ -61,11 +61,3 @@ void AliL3HoughMerge::FillHisto(TH2F *merge_hist) } } -/* -void AliL3HoughMerge::MergeLines() -{ - - - -} -*/ diff --git a/HLT/hough/AliL3HoughTransformer.cxx b/HLT/hough/AliL3HoughTransformer.cxx index 39bf6d71691..8208830515e 100644 --- a/HLT/hough/AliL3HoughTransformer.cxx +++ b/HLT/hough/AliL3HoughTransformer.cxx @@ -26,6 +26,7 @@ AliL3HoughTransformer::AliL3HoughTransformer() fPatch = 0; fRowContainer = 0; fPhiRowContainer = 0; + fVolume=0; fNumOfPadRows=0; fContainerBounds=0; fNDigits=0; @@ -33,21 +34,41 @@ AliL3HoughTransformer::AliL3HoughTransformer() } -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() { @@ -58,6 +79,8 @@ AliL3HoughTransformer::~AliL3HoughTransformer() delete fTransform; if(fPhiRowContainer) delete [] fPhiRowContainer; + if(fVolume) + delete [] fVolume; if(fIndex) { for(Int_t i=0; inextRowPixel) { @@ -100,7 +124,7 @@ void AliL3HoughTransformer::InitTemplates(TH2F *hist) 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"); @@ -210,81 +234,7 @@ void AliL3HoughTransformer::CountBins() } -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: @@ -295,15 +245,21 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist) 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); @@ -311,7 +267,7 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist) } - + printf("Total signal %d\n",totsignal); } /* @@ -431,8 +387,8 @@ void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowra 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]; @@ -479,7 +435,7 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist) 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); @@ -490,19 +446,31 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist) 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; iGetEvent(i); @@ -527,11 +495,12 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist) 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); @@ -548,20 +517,37 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist) ((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()); } diff --git a/HLT/hough/AliL3HoughTransformer.h b/HLT/hough/AliL3HoughTransformer.h index 2833e6a4023..3e6b79045ce 100644 --- a/HLT/hough/AliL3HoughTransformer.h +++ b/HLT/hough/AliL3HoughTransformer.h @@ -9,7 +9,8 @@ struct AliL3Digits UChar_t fPad; UShort_t fTime; Int_t fIndex; - AliL3Digits *nextPhiRowPixel; + AliL3Digits *fNextVolumePixel; + //AliL3Digits *nextPhiRowPixel; AliL3Digits *nextRowPixel; }; @@ -40,6 +41,7 @@ class AliL3HoughTransformer : public TObject { AliL3HoughContainer *fRowContainer; //! AliL3HoughContainer *fPhiRowContainer; //! + AliL3HoughContainer *fVolume; //! Int_t fContainerBounds; Int_t fNDigits; Int_t **fIndex; //! @@ -50,11 +52,11 @@ class AliL3HoughTransformer : public TObject { 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); diff --git a/HLT/hough/Makefile b/HLT/hough/Makefile index a4b97f7efb6..6268a76fefe 100644 --- a/HLT/hough/Makefile +++ b/HLT/hough/Makefile @@ -10,7 +10,7 @@ PACKAGE = AliL3Hough # C++ sources -SRCS = AliL3HoughTransformer.cxx AliL3HoughPixel.cxx \ +SRCS = AliL3HoughTransformer.cxx AliL3Hough.cxx \ AliL3HoughMaxFinder.cxx AliL3HoughEval.cxx AliL3HoughMerge.cxx diff --git a/HLT/hough/hough_merge.C b/HLT/hough/hough_merge.C index 74c7c79b877..db1f017f7e8 100644 --- a/HLT/hough/hough_merge.C +++ b/HLT/hough/hough_merge.C @@ -2,8 +2,8 @@ void hough_merge(char *rootfile,Bool_t MC=false) { 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 @@ -11,91 +11,61 @@ void hough_merge(char *rootfile,Bool_t MC=false) 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; trGetNTracks(); 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; diff --git a/HLT/hough/testPF.C b/HLT/hough/testPF.C index a00161c453b..cdbbd1abf45 100644 --- a/HLT/hough/testPF.C +++ b/HLT/hough/testPF.C @@ -22,26 +22,59 @@ void testPF(char *rootfile,int patch) 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; kReset(); + 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 "<GetNTracks()<<" peaks"<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 "<GetNTracks()<GetNTracks(); i++) { track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i); @@ -50,37 +83,21 @@ void testPF(char *rootfile,int patch) 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; -- 2.39.3