//Author: Anders Strand Vestbo
//Last Modified: 28.6.01
-#include <TClonesArray.h>
-#include <TH2.h>
-#include <TH1.h>
-#include <TFile.h>
-#include <AliRun.h>
-#include <TParticle.h>
-#include <TTree.h>
-
-#include "AliTPCParam.h"
-#include "AliSimDigits.h"
-#include "AliL3TrackArray.h"
-#include "AliL3Transform.h"
+#include "AliL3HoughEval.h"
#include "AliL3HoughTransformer.h"
-#include "AliL3Defs.h"
+#include "AliL3DigitData.h"
#include "AliL3HoughTrack.h"
-#include "AliL3HoughEval.h"
+#include "AliL3Transform.h"
+#include "AliL3Histogram.h"
+#include "AliL3Defs.h"
ClassImp(AliL3HoughEval)
-
AliL3HoughEval::AliL3HoughEval()
{
- //Default constructor
- fTransform = new AliL3Transform();
- fHoughTransformer = 0;
- fNumOfRowsToMiss = 1;
- fNumOfPadsToLook = 1;
-}
+}
AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
{
- //Constructor
+
fHoughTransformer = transformer;
fTransform = new AliL3Transform();
- fNumOfRowsToMiss = 1;
+
+ fSlice = fHoughTransformer->GetSlice();
+ fPatch = fHoughTransformer->GetPatch();
+ fNrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
+ fNEtaSegments = fHoughTransformer->GetNEtaSegments();
+ fEtaMin = fHoughTransformer->GetEtaMin();
+ fEtaMax = fHoughTransformer->GetEtaMax();
+ fRemoveFoundTracks = kFALSE;
fNumOfPadsToLook = 1;
+ fNumOfRowsToMiss = 1;
+ GenerateLUT();
}
-
AliL3HoughEval::~AliL3HoughEval()
{
- //Destructor
if(fTransform)
delete fTransform;
- if(fMcTrackTable)
- delete [] fMcTrackTable;
+ if(fRowPointers)
+ delete [] fRowPointers;
}
-
-Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
+void AliL3HoughEval::GenerateLUT()
{
- //Look at rawdata along the road specified by the track candidates.
+ //Generate a LUT, to limit the access to raw data
+
+ fRowPointers = new AliL3DigitRowData*[fNrows];
+
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
+ if(!tempPt)
+ printf("AliL3HoughEval::GenerateLUT : Zero data pointer\n");
- if(!fHoughTransformer)
+ for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
{
- printf("\nAliL3HoughEval: No transformer object\n");
- return kFALSE;
+ Int_t prow = i - NRows[fPatch][0];
+ fRowPointers[prow] = tempPt;
+ tempPt = fHoughTransformer->UpdateDataPointer(tempPt);
}
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
+}
-
+Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
+{
+ //Look at rawdata along the road specified by the track candidates.
+ //If track is good, return true, if not return false.
+
Int_t sector,row;
- Int_t lut_index;
- Char_t **track_lut = fHoughTransformer->fTrackTable;
+
Int_t nrow=0,npixs=0;
Float_t xyz[3];
-
- for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
+
+ Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
+ for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
{
- Int_t prow = padrow - NRows[patch][0];
+ Int_t prow = padrow - NRows[fPatch][0];
+
if(!track->GetCrossingPoint(padrow,xyz))
{
printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
continue;
}
- fTransform->Slice2Sector(slice,padrow,sector,row);
+ fTransform->Slice2Sector(fSlice,padrow,sector,row);
fTransform->Local2Raw(xyz,sector,row);
npixs=0;
- //Look at both sides of the crossing point:
- for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
+
+ //Get the timebins for this pad
+ AliL3DigitRowData *tempPt = fRowPointers[prow];
+ if(!tempPt)
{
- if(p<0 || p>fTransform->GetNPads(padrow)) continue;
- lut_index = (prow<<8) + p;
- if(track_lut[eta_index][lut_index]>0) //There was a signal here
- npixs++;
+ printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
+ continue;
}
+ //Look at both sides of the pad:
+ for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
+ {
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
+ {
+ UChar_t pad = digPt[j].fPad;
+ if(pad < p) continue;
+ if(pad > p) break;
+ UShort_t time = digPt[j].fTime;
+ Double_t eta = fTransform->GetEta(row,pad,time);
+ Int_t pixel_index = (Int_t)(eta/etaslice);
+ if(pixel_index > eta_index) continue;
+ if(pixel_index != eta_index) break;
+ if(remove)
+ digPt[j].fCharge = 0; //Delete the track from image
+ npixs++;
+ }
+ }
+
if(npixs > 0)
{
nrow++;
}
}
- if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
+ if(remove)
+ return kTRUE;
+
+ if(nrow >= fNrows - fNumOfRowsToMiss)//this was a good track
{
track->SetEtaIndex(eta_index);
- if(remove)
- RemoveTrackFromImage(track,eta_index);
+ if(fRemoveFoundTracks)
+ LookInsideRoad(track,eta_index,kTRUE);
return kTRUE;
}
else
return kFALSE;
}
-void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
+void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
{
- //Evalaute the track candidates by looking along the trajectory.
- //If remove is on, the pixels along the track will be removed.
-
+ //Display the current raw data inside the slice
- if(!fHoughTransformer)
+ if(!hist)
{
- printf("AliL3HoughEval: No transformer object\n");
+ printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
return;
}
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
-
-
- Int_t sector,row;
- Int_t lut_index;
- Char_t **track_lut = fHoughTransformer->fTrackTable;
- Int_t nrow,npixs;
- Float_t xyz[3];
-
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
- if(!track) {printf("No track\n"); return;}
-
- nrow = 0;
-
- for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
- {
- Int_t prow = padrow - NRows[patch][0];
- if(!track->GetCrossingPoint(padrow,xyz))
- {
- printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
- continue;
- }
-
- fTransform->Slice2Sector(slice,padrow,sector,row);
- fTransform->Local2Raw(xyz,sector,row);
- npixs=0;
- //Look at both sides of the crossing point:
- for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
- {
- if(p<0 || p>fTransform->GetNPads(padrow)) continue;
- lut_index = (prow<<8) + p;
- if(track_lut[eta_index][lut_index]>0) //There was a signal here
- npixs++;
- }
-
- if(npixs > 0)
- {
- nrow++;
- }
- }
- if(nrow < NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)
- tracks->Remove(i); //this was not a good enough track
- else if(remove)
- RemoveTrackFromImage(track,eta_index); //this was a good track, so remove it from the image.
- }
-
- tracks->Compress();
-}
-
-void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
-{
- //Remove the pixels along the track in the image.
-
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
-
- Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
-
- Int_t lut_index;
- Char_t **track_lut = fHoughTransformer->fTrackTable;
- Int_t sector,row;
- Float_t xyz[3];
- for(Int_t padrow = NRows[patch][0]; padrow<=NRows[patch][1]; padrow++)
+ Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
+ for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
{
- Int_t prow = padrow - NRows[patch][0];
- if(!track->GetCrossingPoint(padrow,xyz))
+ Int_t prow = padrow - NRows[fPatch][0];
+
+ AliL3DigitRowData *tempPt = fRowPointers[prow];
+ if(!tempPt)
{
- printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
+ printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
continue;
}
- fTransform->Slice2Sector(slice,padrow,sector,row);
- fTransform->Local2Raw(xyz,sector,row);
-
- for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
{
- if(p<0 || p>fTransform->GetNPads(padrow)) continue;
- lut_index = (prow<<8) + p;
- track_lut[eta_index][lut_index] = -1; //remove it!!
- }
- }
-
-
-}
-
-void AliL3HoughEval::DisplaySlice(TH2F *hist)
-{
-
- AliL3Digits *pixel;
-
- for(Int_t padrow = 0; padrow < 174; padrow++)
- {
-
- for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
- {
-
+ UChar_t pad = digPt[j].fPad;
+ UChar_t charge = digPt[j].fCharge;
+ UShort_t time = digPt[j].fTime;
+ Float_t xyz[3];
Int_t sector,row;
- Float_t xyz_pix[3];
- fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
-
- if(hist)
- {
- hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
- }
-
- }
- }
-
-}
-
-void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
-{
- //define the particles that produce good enough signals to be recognized in the transform
-
- Int_t num_eta_segments = fHoughTransformer->fNumEtaSegments;
- Double_t eta_max = fHoughTransformer->fEtaMax;
- fMcTrackTable = new Int_t[num_eta_segments];
- for(Int_t i=0; i<num_eta_segments; i++)
- fMcTrackTable[i]=0;
- Double_t etaslice = eta_max/num_eta_segments;
-
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
-
- TFile *file = new TFile(rootfile);
- file->cd();
-
- AliRun *gAlice = (AliRun*)file->Get("gAlice");
- gAlice->GetEvent(0);
-
- TClonesArray *particles=gAlice->Particles();
- Int_t np = particles->GetEntriesFast();
-
- Int_t row_to_miss = fNumOfRowsToMiss;
- AliTPCParam *param = (AliTPCParam*)file->Get("75x40_100x60");
- Int_t zero = param->GetZeroSup();
- TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
- AliSimDigits da, *digits=&da;
- TD->GetBranch("Segment")->SetAddress(&digits); //Return pointer to branch segment.
- Int_t *good=new Int_t[np];
- Int_t *count = new Int_t[np]; //np number of particles.
- Int_t good_number = NRows[patch][1]-NRows[patch][0]-row_to_miss;
- Int_t i;
- for (i=0; i<np; i++) count[i]=0;
- Int_t sectors_by_rows=(Int_t)TD->GetEntries();
- for (i=0; i<sectors_by_rows; i++)
- {
- if (!TD->GetEvent(i)) continue;
- Int_t sec,row;
- param->AdjustSectorRow(digits->GetID(),sec,row);
- Int_t ss,sr;
- fTransform->Sector2Slice(ss,sr,sec,row);
- if(ss!=slice) continue;
- if(sr < NRows[patch][0]) continue;
- if(sr > NRows[patch][1]) break;
- digits->First();
- while (digits->Next()) {
- Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
- Short_t dig = digits->GetDigit(it,ip);
- Int_t idx0=digits->GetTrackID(it,ip,0);
- Int_t idx1=digits->GetTrackID(it,ip,1);
- Int_t idx2=digits->GetTrackID(it,ip,2);
- if (idx0>=0 && dig>=zero) count[idx0]+=1;
- if (idx1>=0 && dig>=zero) count[idx1]+=1;
- if (idx2>=0 && dig>=zero) count[idx2]+=1;
- }
- for (Int_t j=0; j<np; j++)
- {
- if (count[j]>1) {//at least two digits at this padrow
- good[j]++;
- }
- count[j]=0;
+ fTransform->Slice2Sector(fSlice,padrow,sector,row);
+ fTransform->Raw2Local(xyz,sector,row,pad,time);
+ Double_t eta = fTransform->GetEta(xyz);
+ Int_t pixel_index = (Int_t)(eta/etaslice);
+ if(pixel_index != eta_index) continue;
+ hist->Fill(xyz[0],xyz[1],charge);
}
}
- delete[] count;
-
- Int_t good_one=0;
- //TObjArray *part = new TObjArray(0,0);
- for(i=0; i<np; i++)
- {
- TParticle *p = (TParticle*)particles->UncheckedAt(i);
- if(p->GetFirstMother()>0) continue; //secondary particle
- if(good[i] < good_number) continue; //too few padrows
- Double_t ptg=p->Pt();
- if(ptg<pet) continue;
- Int_t eta_index = (Int_t)(p->Eta()/etaslice);
- if(eta_index < 0 || eta_index > num_eta_segments)
- continue;
- //fMcTrackTable[eta_index]++;
- //part->AddLast(p);
- good_one++;
- }
-
- printf("nparticles %d\n",good_one);
- file->Close();
- delete [] good;
- delete file;
- //return part;
-}
-
-
-void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
-{
-
- Int_t slice = fHoughTransformer->fSlice;
-
- TFile *file = new TFile(rootfile);
- file->cd();
-
- AliRun *gAlice = (AliRun*)file->Get("gAlice");
- gAlice->GetEvent(0);
-
- TClonesArray *particles=gAlice->Particles();
- Int_t n=particles->GetEntriesFast();
- Float_t torad=TMath::Pi()/180;
- Float_t phi_min = slice*20 - 10;
- Float_t phi_max = slice*20 + 10;
-
- for (Int_t j=0; j<n; j++) {
- TParticle *p=(TParticle*)particles->UncheckedAt(j);
- if (p->GetFirstMother()>=0) continue; //secondary particle
- if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
- Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
- Double_t phi_part = TMath::ATan2(pyg,pxg);
- if (phi_part < 0) phi_part += 2*TMath::Pi();
-
- if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
- if(ptg<0.100) continue;
-
- AliL3HoughTrack *sel_track=0;
-
- Double_t min_dist = 10000;
-
- for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
- {
- AliL3HoughTrack *track = (AliL3HoughTrack*)merged_tracks->GetCheckedTrack(t);
- if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
- Float_t phi0[2] = {track->GetPhi0(),0};
- fTransform->Local2GlobalAngle(phi0,slice);
- Double_t dPt = ptg - track->GetPt();
- Double_t dPhi0 = phi_part - phi0[0];
- Double_t distance = pow(dPt,2) + pow(dPhi0,2);
- if(distance < min_dist)
- {
- min_dist = distance;
- sel_track = track;
- }
-
- }
- if(sel_track)
- sel_track->SetBestMCid(j,min_dist);
-
- Double_t dpt = fabs(ptg-sel_track->GetPt());
- Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
- //printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
- }
- file->Close();
- delete file;
}
//Author: Anders Strand Vestbo
-//Last Modified: 28.6.01
+//Last Modified: 14.09.01
-#include <math.h>
-#include <TH2.h>
-#include <time.h>
-
-#include <TFile.h>
-#include <TTree.h>
-
-#include "AliTPCParam.h"
-#include "AliSimDigits.h"
-
-#include "AliL3Histogram.h"
-#include "AliL3Logging.h"
+#include "AliL3HoughTransformer.h"
#include "AliL3Defs.h"
#include "AliL3Transform.h"
-#include "AliL3HoughTransformer.h"
-#include "AliL3HoughTrack.h"
-#include "AliL3TrackArray.h"
#include "AliL3DigitData.h"
+#include "AliL3Histogram.h"
ClassImp(AliL3HoughTransformer)
}
-
-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;
- fNumOfPadRows=NRowsSlice;
-}
-
-AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,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;
+ fNEtaSegments = n_eta_segments;
+ fEtaMin = 0;
+ fEtaMax = fSlice < 18 ? 0.9 : -0.9;
+ fTransform = new AliL3Transform();
+ fThreshold = 0;
- fNRowsInPatch = NRows[fPatch][1]-NRows[fPatch][0] + 1;
- fBinTableBounds = (fNRowsInPatch+1)*(MaxNPads+1);
- fNDigitRowData = 0;
- fDigitRowData = 0;
- fHistoPt = 0;
}
-
AliL3HoughTransformer::~AliL3HoughTransformer()
{
- //Destructor
- if(fBinTable)
- {
- for(Int_t i=0; i<fBinTableBounds; i++)
- delete [] fBinTable[i];
- delete [] fBinTable;
- }
- if(fEtaIndex)
- delete [] fEtaIndex;
- if(fTrackTable)
- {
- for(Int_t i=0; i<fNumEtaSegments; i++)
- delete [] fTrackTable[i];
- delete [] fTrackTable;
- }
if(fTransform)
delete fTransform;
-
-
+ DeleteHistograms();
}
-void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
+void AliL3HoughTransformer::DeleteHistograms()
{
- fNDigitRowData = ndigits;
- fDigitRowData = ptr;
+ if(!fParamSpace)
+ return;
+ for(Int_t i=0; i<fNEtaSegments; i++)
+ delete fParamSpace[i];
+ delete [] fParamSpace;
}
-void AliL3HoughTransformer::InitTables()
+void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
+ Int_t nybin,Double_t ymin,Double_t ymax)
{
- //Create LUT for the circle transform.
- //the actual transformation is done in TransformTables.
-
- AliL3Histogram *hist = fHistoPt;
-
- Int_t nbinsy = hist->GetNbinsY();
-
- fBinTable = new Int_t*[fBinTableBounds];
- Int_t index;
-
- Double_t etaslice = fEtaMax/fNumEtaSegments;
-
- Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);
-
- fEtaIndex = new Int_t[etabounds];
- for(Int_t i=0; i<etabounds; i++)
- fEtaIndex[i] = -1;
-
- fTrackTable = new Char_t*[fNumEtaSegments];
- for(Int_t i=0; i<fNumEtaSegments; i++)
- {
- fTrackTable[i] = new Char_t[fBinTableBounds];
- for(Int_t j=0; j<fBinTableBounds; j++)
- fTrackTable[i][j] = 0;
- }
-
- for(Int_t r=NRows[fPatch][0]; r<=NRows[fPatch][1]; r++)
- {
- Int_t prow = r - NRows[fPatch][0];
- for(Int_t p=0; p<fTransform->GetNPads(r); p++)
- {
- Float_t xyz[3];
- Int_t sector,row;
- for(Int_t t=0; t<fTransform->GetNTimeBins(); t++)
- {
- fTransform->Slice2Sector(fSlice,r,sector,row);
- fTransform->Raw2Local(xyz,sector,row,p,t);
- Double_t eta = fTransform->GetEta(xyz);
- if(eta < fEtaMin || eta > fEtaMax) continue;
- Int_t ind = (prow<<17) + (p<<9) + t;
- if(fEtaIndex[ind]>=0)
- printf("AliL3HoughTransformer::InitTables : Overlapping indexes in eta!!\n");
- Int_t eta_index = (Int_t)(eta/etaslice);
- if(eta_index < 0 || eta_index > fNumEtaSegments)
- continue;
- fEtaIndex[ind] = eta_index;
-
- }
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- Double_t phi_pix = fTransform->GetPhi(xyz);
- index = (prow<<8) + p;
- fBinTable[index] = new Int_t[nbinsy+1];
-
- for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
- {
- Double_t phi0 = hist->GetBinCenterY(b);
- Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
- Int_t bin = hist->FindBin(kappa,phi0);
- if(fBinTable[index][b]!=0)
- printf("AliL3HoughTransformer::InitTables : Overlapping indexes %d %d b %d\n",fBinTable[index][b],index,b);
- fBinTable[index][b] = bin;
- }
- }
-
- }
-
-}
-
-void AliL3HoughTransformer::TransformTables(AliL3Histogram **histos,AliL3Histogram **images)
-{
- //Transform all the pixels while reading, and fill the corresponding histograms.
- //Transform is done using LUT created in InitTables.
- //fTrackTable : table telling whether a specific pixel is active (nonzero):
- //fTrackTable = 0 -> no track
- //fTrackindex > 0 -> track present
- //fTrackindex = -1 -> track has been removed (already found)
- //fEtaIndex : table storing the etaindex -> used to find correct histogram to fill
- //fBinTable : table storing all the bins to fill for each nonzero pixel
-
- Int_t eta_index;
- AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
- AliL3Histogram *hist;
+ fParamSpace = new AliL3Histogram*[fNEtaSegments];
- if(!tempPt)
+ Char_t histname[256];
+ for(Int_t i=0; i<fNEtaSegments; i++)
{
- LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformTables","data")<<
- "Zero datapointer"<<ENDLOG;
- return;
+ sprintf(histname,"paramspace_%d",i);
+ fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
}
-
- Int_t out_count=0,tot_count=0;
- Int_t index,ind;
-
- for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
- {
- Int_t prow = i - NRows[fPatch][0];
- AliL3DigitData *bins = tempPt->fDigitData;
- for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
- {
- index = (prow<<8) + bins[dig].fPad;
- ind = (prow<<17) + (bins[dig].fPad<<9) + bins[dig].fTime;
- eta_index = fEtaIndex[ind];
- if(eta_index < 0) continue; //pixel out of etarange
-
- if(fTrackTable[eta_index][index]<0) continue; //this pixel has already been removed.
- fTrackTable[eta_index][index]++; //this pixel is now marked as active (it is nonzero)
-
- tot_count++;
- hist = histos[eta_index];
-
- if(images)
- {
- //Display the transformed images.
- AliL3Histogram *image = images[eta_index];
- Float_t xyz_local[3];
- Int_t sector,row;
- fTransform->Slice2Sector(fSlice,i,sector,row);
- fTransform->Raw2Local(xyz_local,sector,row,bins[dig].fPad,bins[dig].fTime);
- image->Fill(xyz_local[0],xyz_local[1],bins[dig].fCharge);
- }
-
- if(!hist)
- {
- printf("Error getting histogram!\n");
- continue;
- }
- for(Int_t p=hist->GetFirstYbin(); p<=hist->GetLastYbin(); p++)
- hist->AddBinContent(fBinTable[index][p],bins[dig].fCharge);
-
- }
-
- Byte_t *tmp = (Byte_t*)tempPt;
- Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
- tmp += size;
- tempPt = (AliL3DigitRowData*)tmp;
- }
-
}
-void AliL3HoughTransformer::WriteTables()
-{
- //Write the tables to asci files.
-
- AliL3Histogram *hist = fHistoPt;
- Char_t name[100];
- sprintf(name,"histogram_table_%d.txt",fPatch);
- FILE *file = fopen(name,"w");
-
- Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);
- for(Int_t i=0; i<etabounds; i++)
- {
- if(fEtaIndex[i]<0) continue;
- fprintf(file,"%d %d\n",i,fEtaIndex[i]);
- }
- fclose(file);
-
- sprintf(name,"bin_table_%d.txt",fPatch);
- FILE *file2 = fopen(name,"w");
- for(Int_t i=0; i<fBinTableBounds; i++)
- {
- if(!fBinTable[i]) continue;
- fprintf(file2,"%d ",i);
- for(Int_t j=hist->GetFirstYbin(); j<=hist->GetLastYbin(); j++)
- {
- fprintf(file2,"%d ",fBinTable[i][j]);
- }
- fprintf(file2,"\n");
- }
- fclose(file2);
-}
-
-Double_t AliL3HoughTransformer::CpuTime()
+void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
{
- return (Double_t)(clock()) / CLOCKS_PER_SEC;
+ fNDigitRowData = ndigits;
+ fDigitRowData = ptr;
}
-/*
-void AliL3HoughTransformer::InitTemplates(TH2F *hist)
+AliL3DigitRowData *AliL3HoughTransformer::UpdateDataPointer(AliL3DigitRowData *tempPt)
{
+ //Update the data pointer to the next padrow
- AliL3Digits *pixel;
-
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
- Int_t nbinsy = hist->GetNbinsY();
-
- fIndex = new Int_t*[fNDigits];
- 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)
- {
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-
- Double_t phi_pix = fTransform->GetPhi(xyz);
- Int_t index = pixel->fIndex;
- if(index >= fNDigits)
- printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
- for(Int_t p=ymin; p<=ymax; p++)
- {
-
- 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");
- fIndex[index][p] = bin;
- }
- }
- }
-
+ Byte_t *tmp = (Byte_t*)tempPt;
+ Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
+ tmp += size;
+ return (AliL3DigitRowData*)tmp;
}
-
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
+void AliL3HoughTransformer::Transform()
{
- //Transformation is done with respect to local coordinates in slice.
- //Transform every pixel into whole phirange, using parametrisation:
- //kappa = 2*sin(phi-phi0)/R
- //Assumes you run InitTemplates first!!!!
-
- AliL3Digits *pix1;
-
- Int_t nbinsy = hist->GetNbinsY();
-
- Int_t totsignal=0,vol_index;
- if(fNumEtaSegments==1)
- eta_index=0; //only looking in one etaslice.
+ //Transform the input data with a circle HT.
+
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+ //Set pointer to the data
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
+ if(!tempPt || fNDigitRowData==0)
{
- 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)
- {
- Float_t xyz[3];
- fTransform->Raw2Global(xyz,2,padrow,pix1->fPad,pix1->fTime);
- Double_t eta = fTransform->GetEta(xyz);
- if(eta < fEtaMin || eta > fEtaMax)
- printf("\n Eta OUT OF RANGE\n");
-
- 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("\nAliL3HoughTransformer::TransformTables : No input data!!!\n\n");
+ return;
}
-
- printf("Total signal %d\n",totsignal);
-}
-
-void AliL3HoughTransformer::Transform2Circle(TH2F **histos,Int_t n_eta_segments,UInt_t ndigits,AliL3DigitRowData *ptr)
-{
- //Transform all the pixels while reading them, and fill the corresponding histograms.
- //Everything is done in one go here.
-
- Double_t etaslice = 0.9/n_eta_segments;
- Int_t eta_index;
- AliL3DigitRowData *tempPt = (AliL3DigitRowData*)ptr;
- TH2F *hist;
-
- Int_t out_count=0,tot_count=0;
+
+ Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
{
- printf("doing row %d\n",i);
- for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ if(i != (Int_t)tempPt->fRow)
{
-
- AliL3DigitData *bins = tempPt->fDigitData;
- Float_t xyz[3];
+ printf("AliL3HoughTransform::Transform : Mismatching padrow numbering\n");
+ continue;
+ }
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
+ {
+ UShort_t charge = digPt[j].fCharge;
+ UChar_t pad = digPt[j].fPad;
+ UShort_t time = digPt[j].fTime;
+ if(charge < fThreshold)
+ continue;
Int_t sector,row;
+ Float_t xyz[3];
fTransform->Slice2Sector(fSlice,i,sector,row);
- fTransform->Raw2Local(xyz,sector,row,bins[dig].fPad,bins[dig].fTime);
- Double_t eta = fTransform->GetEta(xyz);
- eta_index = (Int_t)(eta/etaslice);
- if(eta_index < 0 || eta_index >= n_eta_segments)
- {
- //printf("Eta index out of range %d\n",eta_index);
- out_count++;
- continue;
- }
- tot_count++;
- hist = histos[eta_index];
- if(!hist)
- {
- printf("Error getting histogramm!\n");
- return;
- }
-
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
+ fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
+ Float_t eta = fTransform->GetEta(xyz);
+ Int_t eta_index = (Int_t)(eta/etaslice);
+ if(eta_index < 0 || eta_index >= fNEtaSegments)
+ continue;
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- Double_t phi_pix = fTransform->GetPhi(xyz);
- for(Int_t p=ymin; p<=ymax; p++)
+ //Get the correct histogram:
+ AliL3Histogram *hist = fParamSpace[eta_index];
+ if(!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);
- hist->Fill(kappa,phi0,bins[dig].fCharge);
+ printf("AliL3HoughTransformer::Transform : Error getting histogram in index %d\n",eta_index);
+ continue;
}
- }
-
-
- Byte_t *tmp = (Byte_t*)tempPt;
- Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
- tmp += size;
- tempPt = (AliL3DigitRowData*)tmp;
-
- }
-
- printf("There were %d pixels out of range and %d inside\n", out_count,tot_count);
-}
-
-
-void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
-{
-
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
- if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
-
- Double_t xmin = fTransform->Row2X(track->GetFirstRow());
- Double_t xmax = fTransform->Row2X(track->GetLastRow());
-
- Double_t a = -1./tan(track->GetPsiLine());
- Double_t b = track->GetDLine()/sin(track->GetPsiLine());
-
- Double_t ymin = a*xmin + b;
- Double_t ymax = a*xmax + b;
-
- Double_t middle_x = xmin + (xmax-xmin)/2;
- Double_t middle_y = ymin + (ymax-ymin)/2;
-
- Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
- Double_t phi = atan2(middle_y,middle_x);
- Double_t phi0 = 2*phi - track->GetPsiLine();
- Double_t kappa = 2*sin(phi-phi0)/r_middle;
- hist->Fill(kappa,phi0,track->GetWeight());
-
- }
-
-
-}
-
-void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
-{
- //Transform every pixel into histogram, using parametrisation:
- //D = x*cos(psi) + y*sin(psi)
-
- // printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
-
- AliL3Digits *pix1;
-
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
-
- Double_t x0 = fTransform->Row2X(ref_row);
- Double_t y0 = 0;
-
- Int_t sector,row;
- //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
-
- Double_t phi_min = -10*ToRad;
- Double_t phi_max = 10*ToRad;
- Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
-
-
- Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
- Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
-
-
- Int_t index;
-
- for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
- {
- for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
- {
+ //Start transformation
+ Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1] + xyz[2]*xyz[2]);
+ Float_t phi = fTransform->GetPhi(xyz);
- index = phi*fNumOfPadRows + padrow;
- //printf("Looping index %d\n",index);
- if(index > fContainerBounds || index < 0)
- {
- 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)
+ //Fill the histogram along the phirange
+ for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
{
- //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- if(raw)
- raw->Fill(xyz[0],xyz[1],pix1->fCharge);
-
- xyz[0] = xyz[0]-x0;
- xyz[1] = xyz[1]-y0;
-
- //printf("calculating...");
- for(Int_t d=xmin+1; d<=xmax; d++)
- {
- Double_t psi = hist->GetXaxis()->GetBinCenter(d);
- Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
-
- Short_t signal = pix1->fCharge;
- hist->Fill(psi,D,signal);
- //printf("Filling histo, psi %f D %f\n",psi,D);
- }
- //printf("done\n");
+ Float_t phi0 = hist->GetBinCenterY(b);
+ Float_t kappa = 2*sin(phi - phi0)/R;
+ hist->Fill(kappa,phi0,charge);
}
- //printf(" done\n");
}
-
+ tempPt = UpdateDataPointer(tempPt);
}
-
}
-
-void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
-{
-
- //read data from rootfile. more or less obsolete code this.
-
- TFile *file = new TFile(rootfile);
- file->cd();
-
- AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
- TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");
-
- AliSimDigits da, *digarr=&da;
- t->GetBranch("Segment")->SetAddress(&digarr);
- Stat_t num_of_entries=t->GetEntries();
-
- Int_t digit_counter=0;
- Float_t xyz[3];
- Double_t eta;
-
- Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
- printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
-
- if(fRowContainer)
- delete [] fRowContainer;
- fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
- 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 = (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);
- Int_t sector;
- Int_t row;
- param->AdjustSectorRow(digarr->GetID(),sector,row);
- Int_t slice,padrow;
- fTransform->Sector2Slice(slice,padrow,sector,row);
- if(slice != fSlice) continue;
- if(padrow < NRows[fPatch][0]) continue;
- if(padrow > NRows[fPatch][1]) break;
- digarr->First();
- do {
- Int_t time=digarr->CurrentRow();
- Int_t pad=digarr->CurrentColumn();
- Short_t signal=digarr->CurrentDigit();
- if(time < param->GetMaxTBin()-1 && time > 0)
- if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
- && digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
- continue;
-
-
- 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);
- if(hist)
- hist->Fill(xyz[0],xyz[1],signal);
-
- AliL3Digits *dig = new AliL3Digits;
- dig->fIndex = digit_counter;
- digit_counter++;
- dig->fCharge = signal;
- dig->fPad = pad;
- dig->fTime = time;
-
- if(fRowContainer[padrow].first == NULL)
- fRowContainer[padrow].first = (void*)dig;
- else
- ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
- fRowContainer[padrow].last = (void*)dig;
-
- //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("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
- continue;
- }
-
- if(fVolume[index].first == NULL)
- fVolume[index].first = (void*)dig;
- else
- ((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());
-
- }
-
- fNDigits = digit_counter;
- printf("digitcounter %d\n",digit_counter);
- printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
- file->Close();
- delete file;
-
-}
-*/
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff