//$Id$
-// Author: Anders Vestbo <mailto:vestbo$fi.uib.no>
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright © ASV
#include <stream.h>
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
-#include "AliL3Defs.h"
+//_____________________________________________________________
+// AliL3Modeller
+//
+// Class for modeling TPC data.
+//
+// This performs the cluster finding, based on track parameters.
+// Basically it propagates the tracks to all padrows, and looks
+// for a corresponding cluster. For the moment only cog is calculated,
+// and no deconvolution is done.
ClassImp(AliL3Modeller)
{
fMemHandler=0;
fTracks=0;
- fTransform=0;
+ fTrackThreshold=0;
+ fPadOverlap=0;
+ fTimeOverlap=0;
+ fRowData=0;
}
delete fMemHandler;
if(fTracks)
delete fTracks;
- if(fTransform)
- delete fTransform;
-
}
-void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
+void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks)
{
fSlice = slice;
fPatch = patch;
- fPadOverlap=4;
- fTimeOverlap=4;
- fTransform = new AliL3Transform();
+ fPadOverlap=6;
+ fTimeOverlap=8;
+
+ sprintf(fPath,"%s",path);
+
+ AliL3Transform::Init(fPath);
+
fTracks = new AliL3TrackArray("AliL3ModelTrack");
+ Char_t fname[100];
AliL3MemHandler *file = new AliL3MemHandler();
- if(!file->SetBinaryInput("tracks.raw"))
+ sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
+ if(!file->SetBinaryInput(fname))
{
- cerr<<"AliL3Modeller::Init : Error opening trackfile"<<endl;
+ cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
return;
}
file->Binary2TrackArray(fTracks);
file->CloseBinaryInput();
delete file;
- fTracks->QSort();
+ if(houghtracks)
+ cout<<"AliL3Modeller is assuming local hough tracksegments!"<<endl;
+ else
+ cout<<"AliL3Modeller is assuming global tracks!"<<endl;
+
+ if(!houghtracks)
+ fTracks->QSort();
+
for(Int_t i=0; i<fTracks->GetNTracks(); i++)
{
AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
if(!track) continue;
track->Init(fSlice,fPatch);
- track->Rotate(fSlice,kTRUE);
+
+ //Only if the tracks has been merged across sector boundaries:
+ //if(!houghtracks)
+ //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
+
track->CalculateHelix();
}
CalculateCrossingPoints();
- CheckForOverlaps();
+
+ if(!houghtracks)
+ CheckForOverlaps();
fMemHandler = new AliL3MemHandler();
- Char_t fname[100];
- sprintf(fname,"%s/digits_%d_%d.raw",path,fSlice,fPatch);
+ sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
if(!fMemHandler->SetBinaryInput(fname))
{
cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
SetInputData(digits);
}
-void AliL3Modeller::Process()
+void AliL3Modeller::FindClusters()
{
if(!fTracks)
AliL3DigitRowData *rowPt = fRowData;
AliL3DigitData *digPt=0;
- Int_t ntimes = fTransform->GetNTimeBins()+1;
- Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
- Digit *row = new Digit[(ntimes)*(npads)];
+ Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
+ Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
+ Int_t bounds = ntimes*npads;
+ Digit *row = new Digit[bounds];
Int_t seq_charge;
- Int_t pad,time;
+ Int_t pad,time,index;
Short_t charge;
Cluster cluster;
- for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+ for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
{
memset((void*)row,0,ntimes*npads*sizeof(Digit));
digPt = (AliL3DigitData*)rowPt->fDigitData;
charge = digPt[j].fCharge;
row[ntimes*pad+time].fCharge = charge;
row[ntimes*pad+time].fUsed = kFALSE;
+ //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
}
for(Int_t k=0; k<fTracks->GetNTracks(); k++)
{
AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
if(!track) continue;
- if(track->GetOverlap()>=0) continue;//Track is overlapping
+ //if(track->GetOverlap(i)>=0) continue;//Track is overlapping
+
+ if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
+ {
+ track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch.
+ continue;
+ }
Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
- cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
+ //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
pad = hitpad;
time = hittime;
Int_t padsign=-1;
memset(&cluster,0,sizeof(Cluster));
+ Int_t npads=0;
while(1)//Process this padrow
{
+ if(pad < 0 || pad >= AliL3Transform::GetNPads(i))
+ {
+ //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
+ FillCluster(track,&cluster,i,npads);
+ break;
+ }
seq_charge=0;
timesign=-1;
time = hittime;
+
while(1) //Process sequence on this pad:
{
- charge = row[ntimes*pad+time].fCharge;
+ if(time < 0) break;
+ index = ntimes*pad + time;
+ if(index < 0 || index >= bounds)
+ {
+ cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
+ <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
+ break;
+ }
+
+ /*
+ if(row[index].fUsed == kTRUE)//Only use the digits once....
+ charge = 0;
+ else
+ charge = row[index].fCharge;
+ */
+
+ charge = row[index].fCharge;
if(charge==0 && timesign==-1)
- {time=hittime+1; timesign=1; continue;}
+ {
+ if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
+ {
+ time--;
+ continue;
+ }
+ else
+ {
+ time = hittime+1;
+ timesign=1;
+ continue;
+ }
+ }
else if(charge==0 && timesign==1)
- break;
+ {
+ if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
+ {
+ time++;
+ continue;
+ }
+ else
+ {
+ //cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
+ break;
+ }
+ }
+
//cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
seq_charge += charge;
row[ntimes*pad+time].fUsed = kTRUE;
time += timesign;
}
- cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
- if(seq_charge)
- pad += padsign;
-
+
+
+ if(seq_charge)//There was something on this pad.
+ {
+ pad += padsign;
+ npads++;
+ }
else //Nothing more on this pad, goto next pad
{
if(padsign==-1)
- {pad=hitpad+1; padsign=1; continue;}
+ {
+ if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad > 0)
+ {
+ pad--; //In this case, we haven't found anything yet,
+ } //so we will try to expand our search within the natural boundaries.
+ else
+ {
+ pad=hitpad+1;
+ padsign=1;
+ }
+ continue;
+ }
+
+ else if(padsign==1)
+ {
+ if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2 && pad < AliL3Transform::GetNPads(i)-2)
+ {
+ pad++; //In this case, we haven't found anything yet,
+ continue; //so we will try to expand our search within the natural boundaries.
+ }
+ else //We are out of range, or cluster if finished.
+ {
+ //cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
+ FillCluster(track,&cluster,i,npads);
+ break;
+ }
+ }
else //Nothing more in this cluster
{
- Float_t fcharge = (Float_t)cluster.fCharge;
- Float_t fpad = ((Float_t)cluster.fPad/fcharge);
- Float_t ftime = ((Float_t)cluster.fTime/fcharge);
- Float_t sigmaY2,sigmaZ2;
- CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
- //cout<<"row "<<i<<" pad "<<dpad<<endl;
- track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
+ //cout<<"Filling final cluster"<<endl;
+ FillCluster(track,&cluster,i,npads);
break;
}
}
- // pad += padsign;
}
+ //cout<<"done"<<endl;
}
- fMemHandler->UpdateRowPointer(rowPt);
+ FillZeros(rowPt,row);
+ fMemHandler->UpdateRowPointer(rowPt);
}
delete [] row;
+ cout<<"done processing"<<endl;
+
+
+ //Debug:
+ for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+ {
+ AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+ if(!track) continue;
+ if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
+ cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
+ }
}
+void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
+{
+ if(cluster->fCharge==0)
+ {
+ track->SetCluster(row,0,0,0,0,0,0);
+ return;
+ }
+ Float_t fcharge = (Float_t)cluster->fCharge;
+ Float_t fpad = ((Float_t)cluster->fPad/fcharge);
+ Float_t ftime = ((Float_t)cluster->fTime/fcharge);
+ Float_t sigmaY2,sigmaZ2;
+ CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
+ track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
+
+}
+
+void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
+{
+ //Fill zero where data has been used.
+
+ Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
+ AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
+ for(UInt_t j=0; j<rowPt->fNDigit; j++)
+ {
+ Int_t pad = digPt[j].fPad;
+ Int_t time = digPt[j].fTime;
+ if(row[ntimes*pad+time].fUsed==kTRUE)
+ digPt[j].fCharge = 0;
+ }
+}
+
+void AliL3Modeller::WriteRemaining()
+{
+ //Write remaining (nonzero) digits to file.
+
+ AliL3DigitRowData *rowPt;
+ rowPt = (AliL3DigitRowData*)fRowData;
+ Int_t digitcount=0;
+ Int_t ndigits[(AliL3Transform::GetNRows(fPatch))];
+ for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
+ {
+ AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
+ ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
+ for(UInt_t j=0; j<rowPt->fNDigit; j++)
+ {
+ if(digPt[j].fCharge==0) continue;
+ digitcount++;
+ ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
+ }
+ //cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
+ fMemHandler->UpdateRowPointer(rowPt);
+ }
+
+ Int_t size = digitcount*sizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)*sizeof(AliL3DigitRowData);
+ Byte_t *data = new Byte_t[size];
+ memset(data,0,size);
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
+ rowPt = (AliL3DigitRowData*)fRowData;
+
+ for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
+ {
+ Int_t localcount=0;
+ tempPt->fRow = i;
+ tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
+ AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
+ for(UInt_t j=0; j<rowPt->fNDigit; j++)
+ {
+ if(digPt[j].fCharge==0) continue;
+ if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
+ {
+ cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
+ return;
+ }
+ tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
+ tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
+ tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
+
+ localcount++;
+ }
+ if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
+ {
+ cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
+ return;
+ }
+ fMemHandler->UpdateRowPointer(rowPt);
+ Byte_t *tmp = (Byte_t*)tempPt;
+ Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
+ tmp += size;
+ tempPt = (AliL3DigitRowData*)tmp;
+ }
+
+ Char_t fname[100];
+ AliL3MemHandler *mem = new AliL3MemHandler();
+ sprintf(fname,"%s/remains_%d_%d.raw",fPath,fSlice,fPatch);
+ mem->SetBinaryOutput(fname);
+ mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
+ mem->CloseBinaryOutput();
+ delete mem;
+ delete [] data;
+}
+
+
void AliL3Modeller::CalculateCrossingPoints()
{
cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
return;
}
Float_t hit[3];
- for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
+
+
+ //Remove tracks which are no good:
+ for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+ {
+ AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+ if(!track) continue;
+ if(track->GetFirstPointX() > AliL3Transform::Row2X(AliL3Transform::GetLastRow(fPatch)) || track->GetPt()<0.1)
+ fTracks->Remove(i);
+ if(track->GetNHits() < fTrackThreshold)
+ {
+ fTracks->Remove(i);
+ continue;
+ }
+ }
+
+ Int_t sector,row;
+ for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
{
for(Int_t j=0; j<fTracks->GetNTracks(); j++)
{
AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
if(!track) continue;
- //if(!track->GetCrossingPoint(i,hit))
- // fTracks->Remove(j); //Track is bending out.
- track->CalculatePoint(fTransform->Row2X(i));
- if(!track->IsPoint())
+
+ if(!track->GetCrossingPoint(i,hit))
{
+ cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
+ "First point "<<track->GetFirstPointX()<<
+ " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
+ //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
+ //"--------"<<endl;
fTracks->Remove(j);
continue;
}
- hit[1]=track->GetPointY();
- hit[2]=track->GetPointZ();
- fTransform->Local2Raw(hit,fSlice,i);
+ //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
+
+ AliL3Transform::Slice2Sector(fSlice,i,sector,row);
+ AliL3Transform::Local2Raw(hit,sector,row);
+ if(hit[1]<0 || hit[1]>AliL3Transform::GetNPads(i) ||
+ hit[2]<0 || hit[2]>AliL3Transform::GetNTimeBins())
+ {//Track is leaving the patch, so flag the track hits (<0)
+ track->SetPadHit(i,-1);
+ track->SetTimeHit(i,-1);
+ continue;
+ }
+
+
track->SetPadHit(i,hit[1]);
track->SetTimeHit(i,hit[2]);
+
//if(hit[1]<0 || hit[2]>445)
- cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
+ //if(hit[2]<0 || hit[2]>445)
+ //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
//cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
}
}
{
//Flag the tracks that overlap
- cout<<"Checking for overlaps"<<endl;
+ cout<<"Checking for overlaps...";
for(Int_t i=0; i<fTracks->GetNTracks(); i++)
{
{
AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
if(!track2) continue;
- for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
+ for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
{
- if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
- fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
+ if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
+ track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
+ continue;
+ if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
+ fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
{
- track1->SetOverlap(j);
- track2->SetOverlap(i);
+ track2->SetOverlap(k,i);
+ track1->SetOverlap(k,j);
}
}
}
}
-
+ cout<<"done"<<endl;
}
{
Float_t padw,timew;
- if(fPatch < 3)
- padw = fTransform->GetPadPitchWidthLow();
- else
- padw = fTransform->GetPadPitchWidthUp();
+
+ padw = AliL3Transform::GetPadPitchWidth(fPatch);
+
Float_t charge = (Float_t)cl->fCharge;
Float_t pad = (Float_t)cl->fPad/charge;
Float_t time = (Float_t)cl->fTime/charge;
Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
- sigmaY2 = (s2 + 1./12)*padw*padw;
-
- if(s2 != 0)
+
+ //Save the sigmas in pad and time:
+
+ sigmaY2 = (s2);// + 1./12);//*padw*padw;
+
+ /*Constants added by offline
+ if(s2 != 0)
{
- sigmaY2 = sigmaY2*0.108;
- if(fPatch<3)
- sigmaY2 = sigmaY2*2.07;
+ sigmaY2 = sigmaY2*0.108;
+ if(fPatch<3)
+ sigmaY2 = sigmaY2*2.07;
}
-
+ */
+
s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
- timew = fTransform->GetZWidth();
- sigmaZ2 = (s2 +1./12)*timew*timew;
- if(s2 != 0)
+ timew = AliL3Transform::GetZWidth();
+ sigmaZ2 = (s2);// +1./12);//*timew*timew;
+
+
+
+ /*Constants added by offline
+ if(s2 != 0)
{
- sigmaZ2 = sigmaZ2*0.169;
- if(fPatch < 3)
- sigmaZ2 = sigmaZ2*1.77;
+ sigmaZ2 = sigmaZ2*0.169;
+ if(fPatch < 3)
+ sigmaZ2 = sigmaZ2*1.77;
}
-
+ */
}