#include "AliL3StandardIncludes.h"
-#include "AliL3Logging.h"
#include "AliL3ClusterFitter.h"
#include "AliL3FitUtilities.h"
#include "AliL3DigitData.h"
#include "AliL3MemHandler.h"
#include "AliL3HoughTrack.h"
#include "AliL3SpacePointData.h"
-#include "AliL3Compress.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** /class AliL3ClusterFitter
+//<pre>
//_____________________________________________________________
//
// AliL3ClusterFitter
//
-
+</pre>
+*/
ClassImp(AliL3ClusterFitter)
-Int_t AliL3ClusterFitter::fBadFitError=0;
-Int_t AliL3ClusterFitter::fFitError=0;
-Int_t AliL3ClusterFitter::fResultError=0;
-Int_t AliL3ClusterFitter::fFitRangeError=0;
+Int_t AliL3ClusterFitter::fgBadFitError=0;
+Int_t AliL3ClusterFitter::fgFitError=0;
+Int_t AliL3ClusterFitter::fgResultError=0;
+Int_t AliL3ClusterFitter::fgFitRangeError=0;
AliL3ClusterFitter::AliL3ClusterFitter()
{
plane=0;
fNmaxOverlaps = 3;
- fChiSqMax=12;
+ fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
fRowMin=-1;
fRowMax=-1;
fFitted=0;
fClusters=0;
fNMaxClusters=0;
fNClusters=0;
+ fEvent=0;
}
AliL3ClusterFitter::AliL3ClusterFitter(Char_t *path)
strcpy(fPath,path);
plane=0;
fNmaxOverlaps = 3;
- fChiSqMax=12;
+ fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
fRowMin=-1;
fRowMax=-1;
fFitted=0;
fNMaxClusters=100000;
fClusters=0;
fNClusters=0;
+ fEvent=0;
}
AliL3ClusterFitter::~AliL3ClusterFitter()
Float_t hit[3];
track->GetLineCrossingPoint(j,hit);
hit[0] += AliL3Transform::Row2X(track->GetFirstRow());
- Float_t R = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
- hit[2] = R*track->GetTgl();
+ Float_t r = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
+ hit[2] = r*track->GetTgl();
Int_t se,ro;
AliL3Transform::Slice2Sector(slice,j,se,ro);
AliL3Transform::Local2Raw(hit,se,ro);
{
fSlice=slice;
fPatch=patch;
-
+
fRowMin=AliL3Transform::GetFirstRow(patch);
fRowMax=AliL3Transform::GetLastRow(patch);
}
-
-void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline)
+void AliL3ClusterFitter::LoadLocalSegments()
{
+ Char_t filename[1024];
+ sprintf(filename,"%s/hough/tracks_ho_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
+ AliL3MemHandler mem;
+ mem.SetBinaryInput(filename);
+ mem.Binary2TrackArray(fTracks);
+ mem.CloseBinaryInput();
+ for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+ {
+ AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+ if(!track) continue;
+
+ track->CalculateHelix();
+
+ track->Init(fSlice,fPatch);
+
+ for(Int_t j=fRowMin; j<=fRowMax; j++)
+ {
+ //Calculate the crossing point between track and padrow
+
+ Float_t xyzCross[3];
+ if(!track->GetCrossingPoint(j,xyzCross))
+ continue;
+
+ Int_t sector,row;
+ AliL3Transform::Slice2Sector(fSlice,j,sector,row);
+ AliL3Transform::Local2Raw(xyzCross,sector,row);
+
+ if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j) ||
+ xyzCross[2] < 0 || xyzCross[2] >= AliL3Transform::GetNTimeBins()) //track goes out of range
+ continue;
+
+ track->SetPadHit(j,xyzCross[1]);
+ track->SetTimeHit(j,xyzCross[2]);
+
+ Float_t crossingangle = track->GetCrossingAngle(j);
+ track->SetCrossingAngleLUT(j,crossingangle);
+ track->CalculateClusterWidths(j);
+ track->GetClusterModel(j)->fSlice = fSlice;
+
+ }
+ }
+}
+
+void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline,Int_t eventnr,Float_t zvertex)
+{
+ //Function assumes _global_ tracks written to a single file.
cout<<"Loading the seeds"<<endl;
Char_t fname[1024];
+ fEvent = eventnr;
if(offline)
- sprintf(fname,"%s/offline/tracks_%d.raw",fPath,0);
+ sprintf(fname,"%s/offline/tracks_%d.raw",fPath,fEvent);
else
- sprintf(fname,"%s/hough/tracks_%d.raw",fPath,0);
+ sprintf(fname,"%s/hough/tracks_%d.raw",fPath,fEvent);
cout<<"AliL3ClusterFitter::LoadSeeds : Loading input tracks from "<<fname<<endl;
if(fSeeds)
delete fSeeds;
fSeeds = new AliL3TrackArray("AliL3ModelTrack");
+
tfile.Binary2TrackArray(fSeeds);
tfile.CloseBinaryInput();
if(!offline)
{
- /*
- if(track->GetNHits() < 10 || track->GetPt() < 0.08)
+ if(i==0) cerr<<"AliL3ClusterFitter::LoadSeeds : Cutting on pt of 4 GeV!!"<<endl;
+ if(track->GetPt() > 4.)
{
fSeeds->Remove(i);
continue;
}
- */
}
clustercount += track->GetNHits();
track->CalculateHelix();
UInt_t *hitids = track->GetHitNumbers();
Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
-
+ /*
+ if(i==0) cerr<<"Cluster fitter only in HALF TPC!!!"<<endl;
+ if(origslice > 17)
+ {
+ fSeeds->Remove(i);
+ continue;
+ }
+ */
track->Init(origslice,-1);
Int_t slice = origslice;
AliL3Transform::Local2GlobalAngle(&angle,slice);
if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
{
- cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
+ //cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
continue;
//track->Print();
//exit(5);
}
- Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
+ Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
+ xyzCross[2] += zvertex;
Int_t sector,row;
AliL3Transform::Slice2Sector(slice,j,sector,row);
- AliL3Transform::Global2Raw(xyz_cross,sector,row);
- //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
- if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
+ AliL3Transform::Global2Raw(xyzCross,sector,row);
+ //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyzCross[1]<<" time "<<xyzCross[2]<<endl;
+ if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
{
newslice:
Int_t tslice=slice;
- Float_t lastcross=xyz_cross[1];
- if(xyz_cross[1] > 0)
+ Float_t lastcross=xyzCross[1];
+ if(xyzCross[1] > 0)
{
if(slice == 17)
slice=0;
//track->Print();
//exit(5);
}
- xyz_cross[0] = track->GetPointX();
- xyz_cross[1] = track->GetPointY();
- xyz_cross[2] = track->GetPointZ();
+ xyzCross[0] = track->GetPointX();
+ xyzCross[1] = track->GetPointY();
+ xyzCross[2] = track->GetPointZ();
+ xyzCross[2] += zvertex;
Int_t sector,row;
AliL3Transform::Slice2Sector(slice,j,sector,row);
- AliL3Transform::Global2Raw(xyz_cross,sector,row);
- if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
+ AliL3Transform::Global2Raw(xyzCross,sector,row);
+ if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
{
- if(xyz_cross[1] > 0 && lastcross > 0 || xyz_cross[1] < 0 && lastcross < 0)
+ if(xyzCross[1] > 0 && lastcross > 0 || xyzCross[1] < 0 && lastcross < 0)
goto newslice;
else
{
}
}
- if(xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
+ if(xyzCross[2] < 0 || xyzCross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
continue;
- if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j))
+ if(xyzCross[1] < 0 || xyzCross[1] >= AliL3Transform::GetNPads(j))
{
- cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
+ cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyzCross[1]<<" time "<<xyzCross[2]<<endl;
track->Print();
exit(5);
}
- track->SetPadHit(j,xyz_cross[1]);
- track->SetTimeHit(j,xyz_cross[2]);
+ track->SetPadHit(j,xyzCross[1]);
+ track->SetTimeHit(j,xyzCross[2]);
angle=0;
AliL3Transform::Local2GlobalAngle(&angle,slice);
Float_t crossingangle = track->GetCrossingAngle(j,slice);
track->GetClusterModel(j)->fSlice = slice;
}
- memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));
+ memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));//Reset the hitnumbers
track->SetNHits(0);
}
fSeeds->Compress();
- AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
- c->WriteFile(fSeeds,"tracks_before.raw");
- delete c;
-
cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
}
pad = digPt[j].fPad;
time = digPt[j].fTime;
charge = digPt[j].fCharge;
- if(charge > 1024)
- charge -= 1024;
fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
//cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
fFailed++;
}
}
- //FillZeros(rowPt);
AliL3MemHandler::UpdateRowPointer(rowPt);
}
cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
cout<<"Distribution:"<<endl;
- cout<<"Bad fit "<<fBadFitError<<endl;
- cout<<"Fit error "<<fFitError<<endl;
- cout<<"Result error "<<fResultError<<endl;
- cout<<"Fit range error "<<fFitRangeError<<endl;
+ cout<<"Bad fit "<<fgBadFitError<<endl;
+ cout<<"Fit error "<<fgFitError<<endl;
+ cout<<"Result error "<<fgResultError<<endl;
+ cout<<"Fit range error "<<fgFitRangeError<<endl;
}
cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
<<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
<<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
- fFitRangeError++;
+ fgFitRangeError++;
return kFALSE;
}
//Check if any other track contributes to this cluster:
- //This is done by checking if the tracks are overlapping within
- //the range defined by the track parameters
+
for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
{
AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
if(!tr) continue;
if(fSeeding)
if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
- Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))) + 1;
- Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))) + 1;
+
+ if(tr->GetPadHit(row) > padr[0] && tr->GetPadHit(row) < padr[1] &&
+ tr->GetTimeHit(row) > timer[0] && tr->GetTimeHit(row) < timer[1])
+ {
+ if(SetFitRange(tr,padr,timer))
+ track->SetOverlap(row,t);
+ }
+ /*
+ Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))*GetYWidthFactor());
+ Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))*GetZWidthFactor());
if(
(tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1])
)
{
+
if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
- track->SetOverlap(row,t); //Set overlap
+ {
+ track->SetOverlap(row,t); //Set overlap
+ }
}
+ */
}
if(fDebug)
cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
+
FitClusters(track,padr,timer);
+ //CalculateWeightedMean(track,padr,timer);
return kTRUE;
}
Int_t nt = AliL3Transform::GetNTimeBins()+1;
Int_t nsearchbins=0;
- if(row < 63)
+ if(row < AliL3Transform::GetNRowLow())
+ nsearchbins=25;
+ else if(row < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
nsearchbins=25;
else
- nsearchbins=49;
+ nsearchbins=25;
+
+ /*
Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
,0,1,2,3,3,3,3,3,3,3
,2,1,0,-1,-2,-3
,-3,-3,-3,-3,-3,-3,-1,-1};
Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
- ,-3,-3,-3,-3,-2,-1,0,1,2,3
+ ,-3,-3,-3,-3,-2,-1,0,1,2,3
,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
+ */
+
+ Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
+ ,-3,3,-3,3,-3,3,-3,3,-3,3
+ ,0,0,-1,-1,1,1,-2,-2,2,2,-3,-3,3,3};
+ Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
+ ,0,0,-1,-1,1,1,-2,-2,2,2
+ ,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3};
Int_t padhit = (Int_t)rint(track->GetPadHit(row));
Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
}
}
- /*
- if(IsMaximum(padhit,timehit))
- {
- padmax = padhit;
- timemax = timehit;
- }
- */
+
//Define the cluster region of this hit:
//The region we look for, is centered at the local maxima
//and expanded around using the parametrized cluster width
//according to track parameters.
- Int_t xyw = (Int_t)ceil(sqrt(track->GetParSigmaY2(row)))+1;
- Int_t zw = (Int_t)ceil(sqrt(track->GetParSigmaZ2(row)))+1;
+
+ Int_t xyw = 3;
+ Int_t zw = 5;
+
if(padmax>=0 && timemax>=0)
{
- if(fDebug)
- {
- cout<<"Expanding cluster range using expected cluster widths: "<<xyw<<" "<<zw
- <<" and setting local maxima pad "<<padmax<<" time "<<timemax<<endl;
- if(xyw > 10 || zw > 10)
- track->Print();
- }
-
//Set the hit to the local maxima of the cluster.
//Store the maxima in the cluster model structure,
//-only temporary, it will be overwritten when calling SetCluster.
track->GetClusterModel(row)->fDPad = padmax;
track->GetClusterModel(row)->fDTime = timemax;
-
- for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
+
+ Int_t i=padmax,j=timemax;
+ for(Int_t pdir=-1; pdir<=1; pdir+=2)
{
- for(Int_t j=timemax-zw; j<=timemax+zw; j++)
+ i=padmax;
+ while(abs(padmax-i) < xyw)
{
- if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
- if(fRow[nt*i+j].fCharge)
+ Bool_t chargeonpad=kFALSE;
+ for(Int_t tdir=-1; tdir<=1; tdir+=2)
{
- if(i < padrange[0]) padrange[0]=i;
- if(i > padrange[1]) padrange[1]=i;
- if(j < timerange[0]) timerange[0]=j;
- if(j > timerange[1]) timerange[1]=j;
+ j=timemax;
+ while(abs(timemax-j) < zw)
+ {
+ if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) break;
+ if(fRow[nt*i+j].fCharge)
+ {
+ if(i < padrange[0]) padrange[0]=i;
+ if(i > padrange[1]) padrange[1]=i;
+ if(j < timerange[0]) timerange[0]=j;
+ if(j > timerange[1]) timerange[1]=j;
+ chargeonpad=kTRUE;
+ }
+ else
+ break;
+ j+=tdir;
+ }
}
+ if(!chargeonpad)
+ break;
+ i+=pdir;
}
}
+ /*
+ for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
+ {
+ for(Int_t j=timemax-zw; j<=timemax+zw; j++)
+ {
+ if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
+ if(fRow[nt*i+j].fCharge)
+ {
+ if(i < padrange[0]) padrange[0]=i;
+ if(i > padrange[1]) padrange[1]=i;
+ if(j < timerange[0]) timerange[0]=j;
+ if(j > timerange[1]) timerange[1]=j;
+ }
+ }
+ }
+ */
+
if(fDebug)
cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
return kTRUE;
Int_t nt = AliL3Transform::GetNTimeBins()+1;
if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
Int_t charge = fRow[nt*pad+time].fCharge;
- if(charge == 1023 || charge==0) return kFALSE;
+ if(charge == 1024 || charge==0) return kFALSE;
+ //if(charge == 0) return kFALSE;
- fRow[nt*pad+time].fUsed = kTRUE;
- return kTRUE;
-
- //if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
+ //fRow[nt*pad+time].fUsed = kTRUE;
+ //return kTRUE;
+
+ if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
- //if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
+ if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
- //if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
+ if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
- //if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
+ if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
fRow[nt*pad+time].fUsed = kTRUE;
return kTRUE;
}
+void AliL3ClusterFitter::CalculateWeightedMean(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
+{
+ Float_t sum=0;
+ Int_t npads=0;
+ Float_t pad=0,time=0;
+ Int_t nt = AliL3Transform::GetNTimeBins()+1;
+ for(Int_t i=padrange[0]; i<=padrange[1]; i++)
+ {
+ Int_t lsum=0;
+ for(Int_t j=timerange[0]; j<=timerange[1]; j++)
+ {
+ lsum += fRow[nt*(i-1)+(j-1)].fCharge;
+ time += j*fRow[nt*(i-1)+(j-1)].fCharge;
+ }
+ if(lsum)
+ npads++;
+ pad += i * lsum;
+ }
+ if(sum)
+ {
+ pad /= sum;
+ time /= sum;
+ track->SetCluster(fCurrentPadRow,pad,time,sum,0,0,npads);
+ }
+ else
+ track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
+}
+
void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
{
//Handle single and overlapping clusters
-
- //Check whether this cluster has been set before:
+ /*
+ if( (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) <= 1 ||
+ (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) >= AliL3Transform::GetNPads(fCurrentPadRow)-2 ||
+ (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) <= 1 ||
+ (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) >= AliL3Transform::GetNTimeBins()-2)
+ {
+ CalculateWeightedMean(track,padrange,timerange);
+ return;
+ }
+ */
Int_t size = FIT_PTS;
- Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
- if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
+ Int_t maxTracks = FIT_MAXPAR/NUM_PARS;
+ if(track->GetNOverlaps(fCurrentPadRow) > maxTracks)
{
cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
return;
Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
//Check if at least one cluster is not already fitted
- Bool_t all_fitted=kTRUE;
+ Bool_t allFitted=kTRUE;
Int_t k=-1;
while(k < track->GetNOverlaps(fCurrentPadRow))
if(!tr) continue;
if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
{
- all_fitted = kFALSE;
+ allFitted = kFALSE;
break;
}
}
- if(all_fitted)
+ if(allFitted)
{
if(fDebug)
cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
//Fill the fit parameters:
Double_t a[FIT_MAXPAR];
Int_t lista[FIT_MAXPAR];
- Double_t dev[FIT_MAXPAR],chisq_f;
+ Double_t dev[FIT_MAXPAR],chisqF;
- Int_t fit_pars=0;
+ Int_t fitPars=0;
- Int_t n_overlaps=0;
+ Int_t nOverlaps=0;
k=-1;
//Fill the overlapping tracks:
//Use the local maxima as the input to the fitting routine.
//The local maxima is temporary stored in the cluster model:
- Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad); //rint(tr->GetPadHit(fCurrentPadRow));
- Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime); //rint(tr->GetTimeHit(fCurrentPadRow));
+ Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
+ Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
if(fDebug)
exit(5);
}
- a[n_overlaps*NUM_PARS+2] = hitpad;
- a[n_overlaps*NUM_PARS+4] = hittime;
+ a[nOverlaps*NUM_PARS+2] = hitpad;
+ a[nOverlaps*NUM_PARS+4] = hittime;
if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
{
- a[n_overlaps*NUM_PARS+1] = charge;
- a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
- a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
- a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
- lista[n_overlaps*NUM_PARS + 1] = 1;
- lista[n_overlaps*NUM_PARS + 2] = 1;
- lista[n_overlaps*NUM_PARS + 3] = 0;
- lista[n_overlaps*NUM_PARS + 4] = 1;
- lista[n_overlaps*NUM_PARS + 5] = 0;
- lista[n_overlaps*NUM_PARS + 6] = 0;
- fit_pars += 3;
+ a[nOverlaps*NUM_PARS+1] = charge;
+ a[nOverlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
+ a[nOverlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
+ //a[nOverlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
+ lista[nOverlaps*NUM_PARS + 1] = 1;
+ lista[nOverlaps*NUM_PARS + 2] = 1;
+ lista[nOverlaps*NUM_PARS + 3] = 0;
+ lista[nOverlaps*NUM_PARS + 4] = 1;
+ lista[nOverlaps*NUM_PARS + 5] = 0;
+ //lista[nOverlaps*NUM_PARS + 6] = 0;
+ fitPars += 3; //<-------------------
}
else //Cluster was fitted before
{
xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
if(fDebug)
- cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
+ cout<<
endl<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
- a[n_overlaps*NUM_PARS+2] = pad;
- a[n_overlaps*NUM_PARS+4] = time;
- a[n_overlaps*NUM_PARS+1] = charge;
- a[n_overlaps*NUM_PARS+3] = sqrt(xywidth) * GetYWidthFactor();
- a[n_overlaps*NUM_PARS+5] = sqrt(zwidth) * GetZWidthFactor();
- a[n_overlaps*NUM_PARS+6] = sqrt(zwidth) * GetZWidthFactor();
+ a[nOverlaps*NUM_PARS+2] = pad;
+ a[nOverlaps*NUM_PARS+4] = time;
+ a[nOverlaps*NUM_PARS+1] = charge;
+ a[nOverlaps*NUM_PARS+3] = xywidth * GetYWidthFactor();
+ a[nOverlaps*NUM_PARS+5] = zwidth * GetZWidthFactor();
+ //a[nOverlaps*NUM_PARS+6] = zwidth * GetZWidthFactor();
- lista[n_overlaps*NUM_PARS + 1] = 1;
- lista[n_overlaps*NUM_PARS + 2] = 0;
- lista[n_overlaps*NUM_PARS + 3] = 0;
- lista[n_overlaps*NUM_PARS + 4] = 0;
- lista[n_overlaps*NUM_PARS + 5] = 0;
- lista[n_overlaps*NUM_PARS + 6] = 0;
- fit_pars += 1;
+ lista[nOverlaps*NUM_PARS + 1] = 1;
+ lista[nOverlaps*NUM_PARS + 2] = 0;
+ lista[nOverlaps*NUM_PARS + 3] = 0;
+ lista[nOverlaps*NUM_PARS + 4] = 0;
+ lista[nOverlaps*NUM_PARS + 5] = 0;
+ //lista[nOverlaps*NUM_PARS + 6] = 0;
+ fitPars += 1;
}
- n_overlaps++;
+ nOverlaps++;
}
- if(n_overlaps==0) //No clusters here
+ if(nOverlaps==0) //No clusters here
{
delete [] plane;
return;
}
- Int_t pad_num=0;
- Int_t time_num_max=0;
+ Int_t padNum=0;
+ Int_t timeNumMax=0;
Int_t ndata=0;
- Int_t tot_charge=0;
+ Int_t totCharge=0;
if(fDebug)
cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
for(Int_t i=padrange[0]; i<=padrange[1]; i++)
{
- Int_t max_charge = 0;
- Int_t time_num=0;
+ Int_t maxCharge = 0;
+ Int_t timeNum=0;
for(Int_t j=timerange[0]; j<=timerange[1]; j++)
{
Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
if(charge <= 0) continue;
- time_num++;
- if(charge > max_charge)
+ timeNum++;
+ if(charge > maxCharge)
{
- max_charge = charge;
- //time_num++;
+ maxCharge = charge;
+ //timeNum++;
}
if(fDebug)
cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
- tot_charge += charge;
+ totCharge += charge;
ndata++;
if(ndata >= size)
{
y[ndata]=charge;
s[ndata]= 1 + sqrt((Double_t)charge);
}
- if(max_charge) //there was charge on this pad
- pad_num++;
- if(time_num_max < time_num)
- time_num_max = time_num;
+ if(maxCharge) //there was charge on this pad
+ padNum++;
+ if(timeNumMax < timeNum)
+ timeNumMax = timeNum;
}
-
- if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps || ndata <= fit_pars) //too few to do fit
+
+ if(padNum <= 1 || timeNumMax <=1 || nOverlaps > fNmaxOverlaps || ndata <= fitPars) //too few to do fit
{
SetClusterfitFalse(track);
if(fDebug)
- cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<" ndata "<<ndata<<" fit_pars "<<fit_pars<<endl;
+ cout<<"Too few digits or too many overlaps: "<<padNum<<" "<<timeNumMax<<" "<<nOverlaps<<" ndata "<<ndata<<" fitPars "<<fitPars<<endl;
delete [] plane;
return;
}
- Int_t npars = n_overlaps * NUM_PARS;
+ Int_t npars = nOverlaps * NUM_PARS;
if(fDebug)
- cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
- Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
+ cout<<"Number of overlapping clusters "<<nOverlaps<<endl;
+ Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisqF, f2gauss5 );
if(ret<0)
{
SetClusterfitFalse(track);
fFailed++;
- fFitError++;
+ fgFitError++;
delete [] plane;
return;
//exit(5);
}
- chisq_f /= (ndata-fit_pars);
+ chisqF /= (ndata-fitPars);
if(fDebug)
- cout<<"Chisq "<<chisq_f<<endl;
+ cout<<"Chisq "<<chisqF<<endl;
+
+ Bool_t overlapping=kFALSE;
+ if(track->GetNOverlaps(fCurrentPadRow) > 0)//There is a overlap
+ overlapping=kTRUE;
k=-1;
- n_overlaps=0;
+ nOverlaps=0;
while(k < track->GetNOverlaps(fCurrentPadRow))
{
AliL3ModelTrack *tr=0;
{
if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
- if(chisq_f < fChiSqMax)//cluster fit is good enough
+ Int_t lpatch;
+ if(fCurrentPadRow < AliL3Transform::GetNRowLow())
+ lpatch=0;
+ else if(fCurrentPadRow < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
+ lpatch=1;
+ else
+ lpatch=2;
+
+ //if(chisqF < fChiSqMax[(Int_t)overlapping])//cluster fit is good enough
+ if(chisqF < fChiSqMax[lpatch])//cluster fit is good enough
{
- tot_charge = (Int_t)(a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
- Float_t fpad = a[n_overlaps*NUM_PARS+2];
- Float_t ftime = a[n_overlaps*NUM_PARS+4];
- if(tot_charge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) ||
+ totCharge = (Int_t)(2*AliL3Transform::Pi() * a[nOverlaps*NUM_PARS+1] * a[nOverlaps*NUM_PARS+3] * a[nOverlaps*NUM_PARS+5]);
+ Float_t fpad = a[nOverlaps*NUM_PARS+2];
+ Float_t ftime = a[nOverlaps*NUM_PARS+4];
+ if(totCharge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) ||
ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
{
if(fDebug)
cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
- <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<tot_charge<<" xywidth "<<a[n_overlaps*NUM_PARS+3]
- <<" zwidth "<<a[n_overlaps*NUM_PARS+5]<<" peakcharge "<<a[n_overlaps*NUM_PARS+1]<<endl;
+ <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<totCharge<<" xywidth "<<a[nOverlaps*NUM_PARS+3]
+ <<" zwidth "<<a[nOverlaps*NUM_PARS+5]<<" peakcharge "<<a[nOverlaps*NUM_PARS+1]<<endl;
tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
fFailed++;
- fResultError++;
+ fgResultError++;
continue;
}
- tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,0,0,pad_num);
+ tr->SetCluster(fCurrentPadRow,fpad,ftime,totCharge,0,0,padNum);
+ /*
+ tr->SetCluster(fCurrentPadRow,fpad,ftime,totCharge,
+ pow(a[nOverlaps*NUM_PARS+3],2),
+ pow(a[nOverlaps*NUM_PARS+5],2),padNum);
+ */
if(fDebug)
- cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]<<" charge "<<tot_charge<<endl;
+ {
+ cout<<"Setting cluster in slice "<<fSlice<<" row "<<fCurrentPadRow<<" pad "<<a[nOverlaps*NUM_PARS+2]<<" time "<<a[nOverlaps*NUM_PARS+4]
+ <<" padwidth "<<a[nOverlaps*NUM_PARS+3]<<" timewidth "<<a[nOverlaps*NUM_PARS+5]
+ <<" charge "<<totCharge<<endl;
+ }
/*
//Set the digits to used:
for(Int_t i=padrange[0]; i<=padrange[1]; i++)
if(fDebug)
cout<<"Cluster fit was too bad"<<endl;
tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
- fBadFitError++;
+ fgBadFitError++;
fFailed++;
}
}
- n_overlaps++;
+ nOverlaps++;
}
delete [] plane;
i++;
if(!tr) continue;
+ //Set the digit data to unused, so it can be fitted to another bastard:
+ Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
+ Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
+ fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fUsed = kFALSE;
+
tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
}
+
}
tr->Print();
exit(5);
}
-
+
tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
- tr->GetXYWidth(i,xywidth);
- tr->GetZWidth(i,zwidth);
+ tr->GetSigmaY2(i,xywidth);
+ tr->GetSigmaZ2(i,zwidth);
Float_t xyz[3];
Int_t sector,row;
AliL3Transform::Slice2Sector(fSlice,i,sector,row);
- AliL3Transform::Raw2Global(xyz,sector,row,pad,time);
+ AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
if(fNClusters >= fNMaxClusters)
{
cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
exit(5);
}
+
fClusters[fNClusters].fX = xyz[0];
fClusters[fNClusters].fY = xyz[1];
fClusters[fNClusters].fZ = xyz[2];
pat=0;
fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
- if(nhits >= 159)
+ if(nhits >= AliL3Transform::GetNRows())
{
cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
exit(5);
//Copy back the number of assigned clusters
tr->SetNHits(nhits);
+
}
}
-void AliL3ClusterFitter::WriteTracks()
+void AliL3ClusterFitter::WriteTracks(Int_t minHits)
{
if(!fSeeds)
return;
- AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
- c->WriteFile(fSeeds,"tracks_after.raw");
- delete c;
+ AliL3TrackArray *fakes = new AliL3TrackArray();
Int_t clustercount=0;
for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
{
AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
if(!tr) continue;
- if(tr->GetNHits()==0)
- fSeeds->Remove(i);
+ if(tr->GetNHits() < minHits)
+ {
+ fakes->AddLast(tr);
+ fSeeds->Remove(i);
+ }
clustercount += tr->GetNHits();
- /*
- if(tr->GetPt() > 1 && tr->GetNPresentClusters() < 150)
- tr->Print();
- */
}
+
cout<<"Writing "<<clustercount<<" clusters"<<endl;
fSeeds->Compress();
AliL3MemHandler mem;
Char_t filename[1024];
- sprintf(filename,"%s/fitter/tracks_0.raw",fPath);
+ sprintf(filename,"%s/fitter/tracks_%d.raw",fPath,fEvent);
mem.SetBinaryOutput(filename);
mem.TrackArray2Binary(fSeeds);
mem.CloseBinaryOutput();
+ //Write the fake tracks to its own file
+ mem.Free();
+ sprintf(filename,"%s/fitter/tracks_fakes_%d.raw",fPath,fEvent);
+ mem.SetBinaryOutput(filename);
+ mem.TrackArray2Binary(fakes);
+ mem.CloseBinaryOutput();
+ delete fakes;
}
-void AliL3ClusterFitter::WriteClusters()
+void AliL3ClusterFitter::WriteClusters(Bool_t global)
{
AliL3MemHandler mem;
if(fDebug)
cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
Char_t filename[1024];
- sprintf(filename,"%s/fitter/points_0_%d_%d.raw",fPath,fSlice,fPatch);
+ sprintf(filename,"%s/fitter/points_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
mem.SetBinaryOutput(filename);
+ if(global == kTRUE)
+ mem.Transform(fNClusters,fClusters,fSlice);
mem.Memory2Binary(fNClusters,fClusters);
mem.CloseBinaryOutput();
mem.Free();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff