#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"
#include "AliL3DDLDataFileHandler.h"
+#include "AliL3HoughKalmanTrack.h"
#include "TThread.h"
fPeakRatio = 0.5;
fInputFile = 0;
fInputPtr = 0;
+ fRawEvent = 0;
SetTransformerParams();
SetThreshold();
//cout << "Cleaned class mem " << endl;
}
+void AliL3Hough::Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,Float_t zvertex)
+{
+ //Normal constructor
+ fNEtaSegments = netasegments;
+ fVersion = tv;
+ fRawEvent = rawevent;
+ fZVertex = zvertex;
+
+ Init();
+}
+
void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
{
//Normal init of the AliL3Hough
else
#ifdef use_aliroot
{
- if(!fInputFile) {
- if(!fInputPtr) {
- /* In case of reading digits file */
- fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
- if(!fBinary) {
+ if(!fRawEvent) {
+ if(!fInputFile) {
+ if(!fInputPtr) {
+ /* In case of reading digits file */
+ fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
+ if(!fBinary) {
#if use_newio
- if(!fRunLoader) {
+ if(!fRunLoader) {
#endif
- Char_t filename[1024];
- sprintf(filename,"%s/digitfile.root",fPath);
- fMemHandler[i]->SetAliInput(filename);
+ Char_t filename[1024];
+ sprintf(filename,"%s/digitfile.root",fPath);
+ fMemHandler[i]->SetAliInput(filename);
#if use_newio
- }
- else {
- fMemHandler[i]->SetAliInput(fRunLoader);
- }
+ }
+ else {
+ fMemHandler[i]->SetAliInput(fRunLoader);
+ }
#endif
+ }
+ }
+ else {
+ /* In case of reading from DATE */
+ fMemHandler[i] = new AliL3DDLDataFileHandler();
+ fMemHandler[i]->SetReaderInput(fInputPtr,-1);
}
}
else {
- /* In case of reading from DATE */
+ /* In case of reading rawdata from ROOT file */
fMemHandler[i] = new AliL3DDLDataFileHandler();
- fMemHandler[i]->SetReaderInput(fInputPtr,-1);
+ fMemHandler[i]->SetReaderInput(fInputFile);
}
}
else {
- /* In case of reading rawdata from ROOT file */
+ /* In case of reading rawdata using AliRawEvent */
fMemHandler[i] = new AliL3DDLDataFileHandler();
- fMemHandler[i]->SetReaderInput(fInputFile);
+ fMemHandler[i]->SetReaderInput(fRawEvent);
}
}
#else
}
fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",50000);
- fMerger = new AliL3HoughMerger(fNPatches);
- fInterMerger = new AliL3HoughIntMerger();
+ if(fVersion!=4) {
+ fMerger = new AliL3HoughMerger(fNPatches);
+ fInterMerger = new AliL3HoughIntMerger();
+ }
+ else {
+ fMerger = 0;
+ fInterMerger = 0;
+ }
fGlobalMerger = 0;
fBenchmark = new AliL3Benchmark();
}
void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
{
// Setup the parameters for the Hough Transformer
+ // This includes the bin size and limits for
+ // the parameter space histograms
+
Int_t mrow;
Float_t psi=0;
if(patch==-1)
{
// Setup the parameters for the Hough Transformer
-
- Int_t mrow=79;
- Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
- Double_t alpha1 = AliL3Transform::GetMaxY(mrow)/pow(lineradius,2);
- Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
+ Double_t lineradius = 1.0/(AliL3HoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliL3Transform::Pi()*10/180)*tan(AliL3Transform::Pi()*10/180)));
+ Double_t alpha1 = AliL3HoughTransformerRow::GetBeta1()*tan(AliL3Transform::Pi()*10/180);
+ Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/(ptmin*0.9);
Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
// cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
- AliL3HoughTrack track;
- track.SetTrackParameters(kappa,psi,1);
- Float_t hit[3];
- Int_t mrow2 = 158;
- track.GetCrossingPoint(mrow2,hit);
- Double_t lineradius2 = sqrt(pow(AliL3Transform::Row2X(mrow2),2) + pow(AliL3Transform::GetMaxY(mrow2),2));
- Double_t alpha2 = hit[1]/pow(lineradius2,2);
+ Double_t alpha2 = alpha1 - (AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())*tan(psi);
// cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
Int_t i=0;
while(i < 6)
{
- fLowPt[i] = 1.15*alpha1;
+ fLowPt[i] = 1.1*alpha1;
fNBinY[i] = ny;
fNBinX[i] = nx;
- fPhi[i] = 1.15*alpha2;
+ fPhi[i] = alpha2;
+ i++;
+ }
+}
+
+void AliL3Hough::CalcTransformerParams(Float_t ptmin)
+{
+ // Setup the parameters for the Row Hough Transformer
+ // Automatically adjusts the number of bins in X and Y in a way
+ // that the size of the hough bin is 2x (in X) and 2.5 (in Y) the
+ // size of the tpc pads
+
+ Double_t lineradius = 1.0/(AliL3HoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliL3Transform::Pi()*10/180)*tan(AliL3Transform::Pi()*10/180)));
+ Double_t alpha1 = AliL3HoughTransformerRow::GetBeta1()*tan(AliL3Transform::Pi()*10/180);
+ Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/(ptmin*0.9);
+ Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
+ // cout<<"Calculated psi range "<<psi<<endl;
+ Double_t alpha2 = alpha1 - (AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())*tan(psi);
+ alpha1 *= 1.1;
+ // cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<endl;
+
+ Double_t sizex = 2.0*AliL3Transform::GetPadPitchWidthLow()*AliL3HoughTransformerRow::GetBeta1()*AliL3HoughTransformerRow::GetBeta1();
+ Double_t sizey = 2.5*AliL3Transform::GetPadPitchWidthUp()*AliL3HoughTransformerRow::GetBeta2()*AliL3HoughTransformerRow::GetBeta2();
+
+ Int_t nx = 2*(Int_t)(alpha1/sizex)+1;
+ Int_t ny = 2*(Int_t)(alpha2/sizey)+1;
+ // cout<<"Calculated number of bins "<<nx<<" "<<ny<<endl;
+
+ Int_t i=0;
+ while(i < 6)
+ {
+ fLowPt[i] = alpha1;
+ fNBinY[i] = ny;
+ fNBinX[i] = nx;
+ fPhi[i] = alpha2;
i++;
}
}
void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
{
+ // Prepare the parameter space for the processing of
+ // the next read patch. According to the already
+ // accumulated number of gaps in parameter space
+ // bins, the routine updates the dynamic
+ // pointers used in order to jump rapidly during the
+ // filling of the parameter space.
+
char buf[256];
sprintf(buf,"Prepare For Patch %d",nextpatch);
fBenchmark->Start(buf);
Int_t ymax = hist->GetLastYbin();
Int_t nxbins = hist->GetNbinsX()+2;
- UChar_t lastyvalue = 0;
- Int_t endybin = ymin - 1;
- for(Int_t ybin=ymin; ybin<=ymax; ybin++)
- {
- UChar_t lastxvalue = 0;
- UChar_t maxvalue = 0;
- Int_t endxbin = xmin - 1;
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- Int_t bin = xbin + ybin*nxbins;
- UChar_t value = 0;
- if(gapcount[bin] < MAX_N_GAPS) {
- value = 1;
- maxvalue = 1;
- if(tracklastrow[bin] > lastpatchlastrow) {
- if(lastpatchlastrow > currentrowcount[bin])
- gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
- }
- else {
- if(tracklastrow[bin] > currentrowcount[bin])
- gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
+ if(fLastPatch != -1) {
+ UChar_t lastyvalue = 0;
+ Int_t endybin = ymin - 1;
+ for(Int_t ybin=nextrow[ymin]; ybin<=ymax; ybin = nextrow[++ybin])
+ {
+ UChar_t lastxvalue = 0;
+ UChar_t maxvalue = 0;
+ Int_t endxbin = xmin - 1;
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ Int_t bin = xbin + ybin*nxbins;
+ UChar_t value = 0;
+ if(gapcount[bin] < MAX_N_GAPS) {
+ if(tracklastrow[bin] > lastpatchlastrow) {
+ if(lastpatchlastrow > currentrowcount[bin])
+ gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
+ }
+ else {
+ if(tracklastrow[bin] > currentrowcount[bin])
+ gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
+ }
+ if(gapcount[bin] < MAX_N_GAPS) {
+ value = 1;
+ maxvalue = 1;
+ if(trackfirstrow[bin] < nextpatchfirstrow)
+ currentrowcount[bin] = nextpatchfirstrow;
+ else
+ currentrowcount[bin] = trackfirstrow[bin];
+ }
}
- if(trackfirstrow[bin] < nextpatchfirstrow)
- currentrowcount[bin] = nextpatchfirstrow;
- else
- currentrowcount[bin] = trackfirstrow[bin];
- }
- if(fLastPatch != -1) {
if(value > 0)
{
nextbin[xbin + ybin*nxbins] = (UChar_t)xbin;
}
lastxvalue = value;
}
- }
- if(fLastPatch != -1) {
UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
memset(tempnextbin,(UChar_t)(xmax+1),xmax-endxbin);
if(maxvalue > 0)
}
lastyvalue = maxvalue;
}
- }
- if(fLastPatch != -1) {
UChar_t *tempnextrow = nextrow + endybin + 1;
- memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin);
+ memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
+ }
+ else {
+ UChar_t lastyvalue = 0;
+ Int_t endybin = ymin - 1;
+ for(Int_t ybin=ymin; ybin<=ymax; ybin++)
+ {
+ UChar_t maxvalue = 0;
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ Int_t bin = xbin + ybin*nxbins;
+ if(gapcount[bin] < MAX_N_GAPS) {
+ maxvalue = 1;
+ if(trackfirstrow[bin] < nextpatchfirstrow)
+ currentrowcount[bin] = nextpatchfirstrow;
+ else
+ currentrowcount[bin] = trackfirstrow[bin];
+ }
+ }
+ if(maxvalue > 0)
+ {
+ nextrow[ybin] = (UChar_t)ybin;
+ if(maxvalue > lastyvalue)
+ {
+ UChar_t *tempnextrow = nextrow + endybin + 1;
+ memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
+ }
+ endybin = ybin;
+ }
+ lastyvalue = maxvalue;
+ }
+ UChar_t *tempnextrow = nextrow + endybin + 1;
+ memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
}
}
for(Int_t i=0; i<npatches; i++)
{
AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
+ AliL3Histogram *h = tr->GetHistogram(0);
+ Float_t deltax = h->GetBinWidthX()*AliL3HoughTransformerRow::GetDAlpha();
+ Float_t deltay = h->GetBinWidthY()*AliL3HoughTransformerRow::GetDAlpha();
+ Float_t deltaeta = (tr->GetEtaMax()-tr->GetEtaMin())/tr->GetNEtaSegments()*AliL3HoughTransformerRow::GetDEta();
+ Float_t zvertex = tr->GetZVertex();
fTracks[i]->Reset();
fPeakFinder->Reset();
if(hist->GetNEntries()==0) continue;
fPeakFinder->SetHistogram(hist);
fPeakFinder->SetEtaSlice(j);
- fPeakFinder->SetTrackLUTs(((AliL3HoughTransformerRow *)tr)->fTrackNRows,((AliL3HoughTransformerRow *)tr)->fTrackFirstRow,((AliL3HoughTransformerRow *)tr)->fTrackLastRow);
+ fPeakFinder->SetTrackLUTs(((AliL3HoughTransformerRow *)tr)->GetTrackNRows(),((AliL3HoughTransformerRow *)tr)->GetTrackFirstRow(),((AliL3HoughTransformerRow *)tr)->GetTrackLastRow(),((AliL3HoughTransformerRow *)tr)->GetNextRow(j));
#ifdef do_mc
LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
<<"Starting "<<j<<" etaslice"<<ENDLOG;
#endif
fPeakFinder->SetThreshold(fPeakThreshold[i]);
fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
-
- //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
}
for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
{
- if(fPeakFinder->GetWeight(k) < 0) continue;
+ // if(fPeakFinder->GetWeight(k) < 0) continue;
AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
- Float_t psi = atan((fPeakFinder->GetXPeak(k)-fPeakFinder->GetYPeak(k))/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
- Float_t kappa = 2.0*(fPeakFinder->GetXPeak(k)*cos(psi)-AliL3HoughTransformerRow::GetBeta1()*sin(psi));
- // track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
- track->SetTrackParameters(kappa,psi,fPeakFinder->GetWeight(k));
+ Double_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
+ Double_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
+ Double_t eta = (starteta+endeta)/2.0;
+ track->SetTrackParametersRow(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),eta,fPeakFinder->GetWeight(k));
+ track->SetPterr(deltax); track->SetPsierr(deltay); track->SetTglerr(deltaeta);
track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
+ track->SetZ0(zvertex);
Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
track->SetEtaIndex(etaindex);
- Float_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
- Float_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
- track->SetEta((starteta+endeta)/2.0);
- track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
+ Int_t rows[2];
+ ((AliL3HoughTransformerRow *)tr)->GetTrackLength(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),rows);
+ track->SetRowRange(rows[0],rows[1]);
track->SetSector(fCurrentSlice);
track->SetSlice(fCurrentSlice);
#ifdef do_mc
Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
track->SetMCid(label);
- // cout<<"Track found with label "<<label<<" at "<<fPeakFinder->GetXPeak(k)<<" "<<fPeakFinder->GetYPeak(k)<<" with weight "<<fPeakFinder->GetWeight(k)<<endl;
#endif
}
LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
}
}
+#ifdef use_aliroot
+Int_t AliL3Hough::FillESD(AliESD *esd)
+{
+ // Fill the found hough transform tracks
+ // into the ESD. The tracks are stored as
+ // AliESDHLTtrack objects.
+
+ // No TPC PID so far,assuming pions
+ Double_t prob[AliPID::kSPECIES];
+ for(Int_t i=0;i<AliPID::kSPECIES;i++) {
+ if(i==AliPID::kPion) prob[i]=1.0;
+ else prob[i]=0.1;
+ }
+
+ if(!fGlobalTracks) return 0;
+ Int_t nglobaltracks = 0;
+ for(Int_t i=0; i<fGlobalTracks->GetNTracks(); i++)
+ {
+ AliL3HoughTrack *tpt = (AliL3HoughTrack *)fGlobalTracks->GetCheckedTrack(i);
+ if(!tpt) continue;
+
+ if(tpt->GetWeight()<0) continue;
+ AliL3HoughKalmanTrack *tpctrack = new AliL3HoughKalmanTrack(*tpt);
+ if(!tpctrack) continue;
+ AliESDtrack *esdtrack2 = new AliESDtrack() ;
+ esdtrack2->UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
+ esdtrack2->SetESDpid(prob);
+ esd->AddTrack(esdtrack2);
+ nglobaltracks++;
+ delete esdtrack2;
+ delete tpctrack;
+ }
+ return nglobaltracks;
+}
+#endif
+
void AliL3Hough::WriteDigits(Char_t *outfile)
{
//Write the current data to a new rootfile.
void *AliL3Hough::ProcessInThread(void *args)
{
+ // Called in case Hough transform tracking
+ // is executed in a thread
+
AliL3Hough *instance = (AliL3Hough *)args;
Int_t minslice = instance->GetMinSlice();
Int_t maxslice = instance->GetMaxSlice();
void AliL3Hough::StartProcessInThread(Int_t minslice,Int_t maxslice)
{
+ // Starts the Hough transform tracking as a
+ // separate thread. Takes as parameters the
+ // range of TPC slices (sectors) to be reconstructed
+
if(!fThread) {
char buf[255];
sprintf(buf,"houghtrans_%d_%d",minslice,maxslice);
Int_t AliL3Hough::WaitForThreadFinish()
{
+ // Routine is used in case we run the
+ // Hough transform tracking in several
+ // threads and want to sync them before
+ // writing the results to the ESD
+
+#if ROOT_VERSION_CODE < 262403
return TThread::Join(fThread->GetId());
+#else
+ return fThread->Join(fThread->GetId());
+#endif
}