and src, as well as little changes concering the build process.
@echo "ALIHLT_NOLOGGING = $(ALIHLT_NOLOGGING)"
@echo "ALIHLT_DOMC = $(ALIHLT_DOMC)"
@echo "ALIHLT_ALIDETECT = $(ALIHLT_ALIDETECT)"
+ @echo "ALIHLT_ROWHOUGH = $(ALIHLT_ROWHOUGH)"
help:
cat doc/README
ifeq ($(ALICE_LEVEL),ali-head)
ALIHLT_USENEWIO = true
endif
-ifeq ($(ALICE_LEVEL),ali-v4-01-00)
-ALIHLT_USENEWIO = true
endif
+
+ifeq ($(ALIHLT_ROWHOUGH),true)
+USEROWHOUGH = 1
endif
ifeq ($(ALIHLT_USENEWIO),true)
USENEWIO = 1
endif
-
#----------------------------------------------------
#Some compiler flags or defines: You can use your
#own setting by defining them outside (make -e)
#----------------------------------------------------
-GCCVERSION = $(shell $(CXX) --version | head -n 1 | cut -d" " -f 3 | cut -d. -f 1 | cut -d" " -f1)
-CXXGCC3FLAGS = -DGCCVERSION=$(GCCVERSION)
-PROFILEFLAGS = -g -pg
-
ifeq ($(ARCH),Darwin)
FINKDIR = /sw
CXX = g++
DYFLAGS = -dynamiclib -flat_namespace -undefined suppress \
-compatibility_version 1 -current_version 1.0.0
LDFLAGS = -O $(EXTRALDFLAGS) -L/sw/lib -ldl
-LDSTATIC = ar
-STATICFLAGS = rucs
else
CXX = g++
CXXFLAGS = -O2 -fPIC -Wall $(CXXGCC3FLAGS) $(EXTRACXXFLAGS)
LD = $(CXX)
LDFLAGS = -O2 $(EXTRALDFLAGS)
SOFLAGS = -shared
-LDSTATIC = ar
-STATICFLAGS = rucs
endif
+
+#static flags for profiling
+PROFILEFLAGS = -g -pg
+LDSTATIC = ar
+STATICFLAGS = rucs
DEFSTR += -Ddo_mc
endif
+ifeq ($(USEROWHOUGH),1)
+DEFSTR += -DROWHOUGH
+endif
+
ifneq ($(NOLOGGING),1)
DEFSTR += -Duse_logging
ifdef ALIHLT_MLUCDIR
$(DICT): $(HDRS)
@echo "Generating dictionary..."
- rootcint -f $(DICT) -c $(CINTCXXFLAGS) $(INCLUDES) -DGCCVERSION=$(GCCVERSION) \
+ rootcint -f $(DICT) -c $(CINTCXXFLAGS) $(INCLUDES) \
$(DEFSTR) -include AliL3StandardIncludes.h $(HDRS)
$(OBJDIR)/%.o: %.cxx
@echo "ALIHLT_NOLOGGING = $(ALIHLT_NOLOGGING)"
@echo "ALIHLT_DOMC = $(ALIHLT_DOMC)"
@echo "ALIHLT_ALIDETECT = $(ALIHLT_ALIDETECT)"
+ @echo "ALIHLT_ROWHOUGH = $(ALIHLT_ROWHOUGH)"
@echo "ROOTSTR = $(ROOTSTR)"
@echo "ALIROOTSTR = $(ALIROOTSTR)"
setenv ALIHLT_NOLOGGING false
setenv ALIHLT_DOMC true
setenv ALIHLT_ALIDETECT true
-
+setenv ALIHLT_ROWHOUGH false
setenv ALIHLT_MLUCDIR ${ALIHLT_BASEDIR}/kip/MLUC
setenv ALIHLT_DATADIR /data1/head
export ALIHLT_NOLOGGING=false
export ALIHLT_DOMC=true
export ALIHLT_ALIDETECT=true
-
+export ALIHLT_ROWHOUGH=false
export ALIHLT_MLUCDIR=/usr/local/kip/MLUC
#export ALIHLT_DATADIR=/mnt/local/alidata/head
#export ALIHLT_TRANSFORMFILE=$ALIHLT_DATADIR/l3transform.config
#export ALIHLT_GEOPATH=$ALIDATADIR
-
if test -z "$LD_LIBRARY_PATH"; then
export LD_LIBRARY_PATH=$ALIHLT_MLUCDIR/lib:$ALIHLT_LIBDIR
elif test -z "`echo $LD_LIBRARY_PATH | grep $ALIHLT_MLUCDIR/lib`";
setenv ALIHLT_NOLOGGING false
setenv ALIHLT_DOMC true
-setenv ALIHLT_ALIDETECT true
+setenv ALIHLT_ALIDETECT false
+setenv ALIHLT_ROWHOUGH true
setenv ALIHLT_MLUCDIR ${ALIHLT_BASEDIR}/kip/MLUC
-setenv ALIHLT_DATADIR /data1/head
+#setenv ALIHLT_DATADIR /data1/head
#setenv ALIHLT_TRANSFORMFILE ${ALIHLT_DATADIR}/l3transform.config
#setenv ALIHLT_GEOPATH ${ALIHLT_DATADIR}
export ALIHLT_NOLOGGING=true
export ALIHLT_DOMC=true
export ALIHLT_ALIDETECT=false
-
+export ALIHLT_ROWHOUGH=true
#export ALIHLT_MLUCDIR=/usr/local/kip/MLUC
#export ALIHLT_DATADIR=/mnt/local/alidata/head
#include "AliL3SpacePointData.h"
#include "AliL3Compress.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** /class AliL3ClusterFitter
+//<pre>
//_____________________________________________________________
//
// AliL3ClusterFitter
//
-
+</pre>
+*/
ClassImp(AliL3ClusterFitter)
}
}
-void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline,Int_t eventnr)
+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;
//exit(5);
}
Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
+ xyz_cross[2] += zvertex;
Int_t sector,row;
AliL3Transform::Slice2Sector(slice,j,sector,row);
xyz_cross[0] = track->GetPointX();
xyz_cross[1] = track->GetPointY();
xyz_cross[2] = track->GetPointZ();
+ xyz_cross[2] += zvertex;
Int_t sector,row;
AliL3Transform::Slice2Sector(slice,j,sector,row);
AliL3Transform::Global2Raw(xyz_cross,sector,row);
fClusters=0;
fNClusters=0;
}
-
void Init(Int_t slice,Int_t patch,Int_t *rowrange,AliL3TrackArray *tracks);
void Init(Int_t slice,Int_t patch);
- void LoadSeeds(Int_t *rowrange,Bool_t offline=kTRUE,Int_t eventnr=0);
+ void LoadSeeds(Int_t *rowrange,Bool_t offline=kTRUE,Int_t eventnr=0,Float_t zvertex=0.0);
void LoadLocalSegments();
void FindClusters();
void AddClusters();
Float_t GetZWidthFactor() {return fCurrentPadRow < AliL3Transform::GetLastRow(1) ? fZInnerWidthFactor : fZOuterWidthFactor;}
AliL3TrackArray *GetSeeds() {return fSeeds;}
-
ClassDef(AliL3ClusterFitter,1)
};
#include "AliL3Compress.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3DataCompressorHelper.h"
#include "AliL3DataCompressor.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3DataCompressorHelper.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3ModelTrack.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3FileHandler.h"
#endif
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3OfflineDataCompressor.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
Compiliation flags:
$ALIHLT_USENEWIO (false) --> switch to NEWIO
-$ALIHLT_NOLOGGING (false) --> use HLT logger classes (requires MLUC lib)
+$ALIHLT_ROWHOUGH (false) --> use row hough transformer
$ALIHLT_DOMC (true) --> store Monte Carlo ids (for ALIROOT package)
+$ALIHLT_NOLOGGING (false) --> use HLT logger classes (requires MLUC lib)
$ALIHLT_ALIDETECT (true) --> detect ALIROOT version using cvs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/bin/sethlt_cern.csh, /alice/cvs/hltcvs/level3code/bin/sethlt_cern.sh, /alice/cvs/hltcvs/level3code/bin/sethlt.csh, /alice/cvs/hltcvs/level3code/bin/sethlt.sh:
+ Added ALIHLT_ROWHOUGH variable.
+
2003-07-29 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/bin/sethlt_cern.csh, /alice/cvs/hltcvs/level3code/bin/sethlt_cern.sh, /alice/cvs/hltcvs/level3code/bin/sethlt.csh, /alice/cvs/hltcvs/level3code/bin/sethlt.sh, /alice/cvs/hltcvs/level3code/bin/usehlt.csh, /alice/cvs/hltcvs/level3code/bin/usehlt.sh:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/comp/AliL3ClusterFitter.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3Compress.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3DataCompressor.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3DataCompressorHelper.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3Modeller.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3ModelTrack.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3OfflineDataCompressor.cxx:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
+ * /alice/cvs/hltcvs/level3code/comp/AliL3ClusterFitter.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3ClusterFitter.h:
+ Added zvertex.
+
2004-02-25 Anders Strand Vestbo <vestbo@hansa00>
* /alice/cvs/hltcvs/level3code/comp/AliL3Compress.cxx, /alice/cvs/hltcvs/level3code/comp/AliL3Compress.h:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/doc/README: Added ALIHLT_ROWHOUGH info.
+
+ * /alice/cvs/hltcvs/level3code/doc/taginfo:
+ Added new tag before big Cvetan changes.
+
2004-01-22 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/doc/taginfo: Added new tag.
+2004-03-30 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/exa/evaltracker.C: Make it compilable.
+
+ * /alice/cvs/hltcvs/level3code/exa/runrowhough.C:
+ Little changes to have default params.
+
+ * /alice/cvs/hltcvs/level3code/exa/runtracker.C:
+ Save benchmark with bname.
+
+ * /alice/cvs/hltcvs/level3code/exa/evalrowhough.C:
+ Tried to make it compilable, but failled.
+
+2004-03-28 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/exa/runtracker.C:
+ Added name for tracker eval.
+
+ * /alice/cvs/hltcvs/level3code/exa/runrowhough.C:
+ Changes to compile macro.
+
+ * /alice/cvs/hltcvs/level3code/exa/rootlogon.C: Added runrowhough.
+
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/exa/evalrowhough.C, /alice/cvs/hltcvs/level3code/exa/runrowhough.C:
+ Added examples for new row hough transform.
+
2004-02-02 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/exa/deconvclusters.C, /alice/cvs/hltcvs/level3code/exa/rootlogon.C, /alice/cvs/hltcvs/level3code/exa/runhough.C, /alice/cvs/hltcvs/level3code/exa/runtracker.C, /alice/cvs/hltcvs/level3code/exa/SetFitParameters.C, /alice/cvs/hltcvs/level3code/exa/SetHoughParameters.C:
+2004-03-28 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughBaseTransformer.h:
+ Added old type for label.
+
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3Histogram1D.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HistogramAdaptive.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3Histogram.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughClusterTransformer.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3Hough.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughDisplay.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughEval.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughIntMerger.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughMaxFinder.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughMerger.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTest.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTrack.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformer.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerGlobal.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerLUT.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerNew.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerRow.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerVhdl.cxx:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughEval.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughEval.h:
+ Added zvertex.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3Hough.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3Hough.h:
+ Added new row transformer (using special peakfinder method)
+ Added in addition support for reading data directly from
+ DATE, so that this version of hough transform can run directly
+ on the GDCs or during some data challenges for testing.
+ (changed methods init and constructor)
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3Histogram.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3Histogram.h:
+ Added GetPreciseBinCenter functions.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughBaseTransformer.cxx:
+ Cosmetics.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughMaxFinder.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughMaxFinder.h:
+ Added FindAdaptedRowPeaks function for row transformer.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughTrack.cxx:
+ Added ROWHOUGH def for mc label.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughBaseTransformer.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughBaseTransformer.h:
+ Added zvertex information.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughLinkDef.h, /alice/cvs/hltcvs/level3code/hough/Makefile:
+ Added new row transformer and removed old gap transformer.
+
+ * /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerGap.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerGap.h, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerRow.cxx, /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerRow.h:
+ Added Cvetans new version of fast Hough Transform called "Counting
+ Gaps and Rows." The idea is to count the number of consecutive rows
+ over gaps per entry of the Hough space. (The old version called
+ AliL3HoughTransformGap is removed.)
+
2004-02-12 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/hough/AliL3HoughTransformerGap.cxx:
+2004-03-18 Thomas Vik <tvik@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/kalman/AliL3KalmanTrack.cxx:
+ Change in calling AliL3Transfrom::Global2Local command.
+
+ * /alice/cvs/hltcvs/level3code/kalman/AliL3KalmanTrack.h, /alice/cvs/hltcvs/level3code/kalman/AliL3KalmanTrack.cxx, /alice/cvs/hltcvs/level3code/kalman/AliL3Kalman.h, /alice/cvs/hltcvs/level3code/kalman/AliL3Kalman.cxx:
+ Changed global coordinates of padrows to local. Removed obsolete MakeSeed function.
+
2004-02-02 Thomas Vik <tvik@hansa00>
* /alice/cvs/hltcvs/level3code/kalman/AliL3Kalman.cxx:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/misc/AliL3DataHandler.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3DDLDataFileHandler.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3DDLRawReaderFile.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3Stopwatch.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3TPCMapping.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3TransBit.cxx, /alice/cvs/hltcvs/level3code/misc/AliL3VHDLClusterFinder.cxx:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
+ * /alice/cvs/hltcvs/level3code/misc/AliL3DDLRawReaderFile.cxx:
+ Added namespace for gcc-3.x
+
2004-02-03 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/misc/AliL3DDLRawReaderFile.cxx:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/programs/ali2raw.cxx, /alice/cvs/hltcvs/level3code/programs/convcosmicsfile.cxx, /alice/cvs/hltcvs/level3code/programs/gettransform.cxx, /alice/cvs/hltcvs/level3code/programs/Makefile, /alice/cvs/hltcvs/level3code/programs/read.cxx, /alice/cvs/hltcvs/level3code/programs/runhough.cxx, /alice/cvs/hltcvs/level3code/programs/runit.cxx, /alice/cvs/hltcvs/level3code/programs/runtracker.cxx, /alice/cvs/hltcvs/level3code/programs/runvhdlcf.cxx, /alice/cvs/hltcvs/level3code/programs/runvhdlhough.cxx:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
2004-02-25 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/programs/gettransform.cxx, /alice/cvs/hltcvs/level3code/programs/Makefile, /alice/cvs/hltcvs/level3code/programs/runhough.cxx:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/src/AliL3Transform.cxx:
+ Update version information.
+
+ * /alice/cvs/hltcvs/level3code/src/AliL3ClustFinderNew.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Display.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Evaluate.cxx, /alice/cvs/hltcvs/level3code/src/AliL3FileHandler.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Logger.h, /alice/cvs/hltcvs/level3code/src/AliL3Logging.h, /alice/cvs/hltcvs/level3code/src/AliL3RawDataFileHandler.cxx, /alice/cvs/hltcvs/level3code/src/AliL3RootTypes.h, /alice/cvs/hltcvs/level3code/src/AliL3StandardIncludes.h, /alice/cvs/hltcvs/level3code/src/AliL3Track.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Transform.cxx:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
+ * /alice/cvs/hltcvs/level3code/src/AliL3TrackArray.cxx, /alice/cvs/hltcvs/level3code/src/AliL3TrackSegmentData.h:
+ Added ROWHOUGH flag to store weight and track id
+ in the track segment data structure. This is used
+ for the new row transformer and probably only
+ a temporary solution, as the proper way would be
+ to have a new AliL3HoughTrackRow class.
+
+ By default the flag $ALIHLT_ROWHOUGH points to false,
+ so there is no change to earlier versions.
+
+2004-03-15 Anders Strand Vestbo <vestbo@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/src/AliL3ConfMapFit.cxx, /alice/cvs/hltcvs/level3code/src/AliL3ConfMapTrack.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Fitter.cxx, /alice/cvs/hltcvs/level3code/src/AliLevel3.cxx, /alice/cvs/hltcvs/level3code/src/AliLevel3.h:
+ Bugfix related to track fit parameters <-> global track merging.
+
+ It turned out that there was a slight problem with the global track
+ merging in the case a vertex constraint was imposed in the circle fit.
+ The problem is quite involved, and is due to inconsistency when storing
+ the socalled first point on the track in AliL3Track::*fFirstPoint. This
+ point is set to the first point lying on the fit, but in the case of vertex
+ constraint this point is correspondingly set to the point of closest
+ approach to the vertex. HOWEVER, the global merger always assumes that
+ this point is the first associated cluster on the track, and thus we have
+ a conflict.
+
+ The (temporary) solution to the problem is to set the *fFirstPoint to the
+ innermost cluster on the track ALWAYS, i.e. in both vertexconstraint and
+ no vertex constraint. This means also that the linear fit in (s,z) space
+ does not currently include the vertex in the fit, in order to also set
+ Z0 to the innermost cluster of the track.
+
2004-02-12 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/src/AliL3Transform.cxx, /alice/cvs/hltcvs/level3code/src/AliL3Transform.h:
+2004-03-20 Constantin Loizides <loizides@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/Makefile.conf, /alice/cvs/hltcvs/level3code/Makefile.rules:
+ Removed GCCVERSION string from code. Replaced by compiler internal
+ macro __GNUC__ which is 3 for gcc version >3.
+
+ * /alice/cvs/hltcvs/level3code/Makefile.conf, /alice/cvs/hltcvs/level3code/Makefile.rules:
+ Added ROWHOUGH flag to store weight and track id
+ in the track segment data structure. This is used
+ for the new row transformer and probably only
+ a temporary solution, as the proper way would be
+ to have a new AliL3HoughTrackRow class.
+
+ By default the flag $ALIHLT_ROWHOUGH points to false,
+ so there is no change to earlier versions.
+
2004-02-25 Constantin Loizides <loizides@hansa00>
* /alice/cvs/hltcvs/level3code/Makefile.conf, /alice/cvs/hltcvs/level3code/Makefile.rules:
+2004-03-12 Gaute Ovrebekk <ovrebekk@hansa00>
+
+ * /alice/cvs/hltcvs/level3code/trigger/RunD0Trigger.C: small chages
+
+ * /alice/cvs/hltcvs/level3code/trigger/AliD0Trigger.cxx, /alice/cvs/hltcvs/level3code/trigger/AliD0Trigger.h:
+ Added more methodes
+
+ * /alice/cvs/hltcvs/level3code/trigger/RunD0offline.C:
+ Running with offline methids
+
2004-02-24 Gaute Ovrebekk <ovrebekk@hansa00>
* /alice/cvs/hltcvs/level3code/trigger/AliD0Trigger.cxx:
v1-2: First real release to CERN (with AliRoot 3-09-Release)
v1-3: Deconvoluter and Hough changes tested (with AliRoot 3-09-Release)
v1-4: Merge with cern-hlt tree to support NEWIO
+v1-5: Fitter and Track updates and status before intro of MACOSX handling
+v1-6: Compression and tracking updates as of March, 20.
+
+
--- /dev/null
+// $Id$
+
+struct GoodTrack
+{
+ Int_t label;
+ Double_t eta;
+ Int_t code;
+ Double_t px,py,pz;
+ Double_t x,y,z;
+ Int_t nhits;
+ Int_t sector;
+};
+typedef struct GoodTrack GoodTrack;
+
+//Histograms
+TNtuple *fNtuppel;
+TH1F *fPtRes;
+TH1F *fNGoodTracksPt;
+TH1F *fNFoundTracksPt;
+TH1F *fNFakeTracksPt;
+TH1F *fTrackEffPt;
+TH1F *fFakeTrackEffPt;
+TH1F *fNGoodTracksEta;
+TH1F *fNFoundTracksEta;
+TH1F *fNFakeTracksEta;
+TH1F *fTrackEffEta;
+TH1F *fFakeTrackEffEta;
+TH1F *fNGoodTracksPhi;
+TH1F *fNFoundTracksPhi;
+TH1F *fNFakeTracksPhi;
+TH1F *fTrackEffPhi;
+TH1F *fFakeTrackEffPhi;
+TProfile *fFakesPt;
+TProfile *fFakesPhi;
+TProfile *fSlicesPt;
+TProfile *fSlicesPhi;
+TProfile2D *fFakesPtvsPhi;
+TNtuple *fNtupleRes;
+TH1F *fPtRes2;
+TH1F *fEtaRes;
+TH1F *fPsiRes;
+
+TH1F *fNGoodTracksRad;
+TH1F *fNFoundTracksRad;
+TH1F *fNGoodTracksZ;
+TH1F *fNFoundTracksZ;
+TH1F *fRadEff;
+TH1F *fZEff;
+
+void evalrowhough(Char_t *path="./fitter",Char_t *offlinepath="./",Int_t ievent=0, Float_t ptmin=0.4)
+{
+ CreateHistos(25,0.1,4.1);
+
+ Char_t fname[1024];
+ sprintf(fname,"%s/tracks_%d.raw",path,ievent);
+ AliL3FileHandler *tfile = new AliL3FileHandler();
+ if(!tfile->SetBinaryInput(fname)){
+ LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
+ <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
+ return;
+ }
+
+ AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
+ tfile->Binary2TrackArray(fTracks);
+ //fTracks->QSort();
+ tfile->CloseBinaryInput();
+ delete tfile;
+
+ for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+ {
+ AliL3HoughTrack *track = fTracks->GetCheckedTrack(i);
+ if(!track) continue;
+
+ track->SetEta(fTracks->GetCheckedTrack(i)->GetPseudoRapidity());
+ UInt_t *ids = track->GetHitNumbers();
+ Int_t slice = (ids[0]>>25)&0x7f;
+ track->SetSlice(slice);
+ cout<<"Track with label "<<ids<<" "<<track->GetNHits()<<" "<<track->GetMCid()<<" "<<ids[0]<<" "<<track->GetWeight()<<" "<<track->GetSlice()<<endl;
+ }
+
+ Char_t filename[1024];
+ sprintf(filename,"%s/good_tracks_tpc_%d",offlinepath,ievent);
+ ifstream in(filename);
+ if(!in)
+ {
+ cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
+ return;
+ }
+
+ Int_t MaxTracks=20000;
+ Int_t fGoodGen=0;
+ fGoodTracks = new GoodTrack[MaxTracks];
+
+ while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
+ fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
+ fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
+ {
+ fGoodTracks[fGoodGen].nhits=-1;
+ fGoodTracks[fGoodGen].sector=-1;
+ fGoodGen++;
+ if (fGoodGen==MaxTracks)
+ {
+ cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
+ break;
+ }
+ }
+
+ Int_t fGoodGen2 = 0;
+ Int_t fGoodFound = 0;
+ Int_t fGoodFound2 = 0;
+ Float_t fMinGoodPt = ptmin;
+ Float_t fMaxGoodPt = 99999.;
+
+ if(!fGoodTracks)
+ {
+ cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
+ return;
+ }
+ cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
+ for(Int_t i=0; i<fGoodGen; i++)
+ {
+ //cout<<"Checking particle "<<i<<endl;
+ Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
+ Float_t pg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py + fGoodTracks[i].pz*fGoodTracks[i].pz);
+ if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
+ Float_t pzg=fGoodTracks[i].pz;
+ Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
+ Float_t etag=0.5 * TMath::Log((pg+pzg)/(pg-pzg));
+ Float_t phig=TMath::ATan2(fGoodTracks[i].py,fGoodTracks[i].px);
+ Float_t rg = TMath::Sqrt(fGoodTracks[i].x*fGoodTracks[i].x + fGoodTracks[i].y*fGoodTracks[i].y);
+ Float_t zg = TMath::Abs(fGoodTracks[i].z);
+ if (phig<0) phig+=2*TMath::Pi();
+
+ fNGoodTracksPt->Fill(ptg);
+ fNGoodTracksEta->Fill(etag);
+ fNGoodTracksPhi->Fill(phig);
+ if((zg != 0) || (rg != 0)) {
+ fNGoodTracksRad->Fill(rg);
+ fNGoodTracksZ->Fill(zg);
+ }
+ fGoodGen2++;
+ Int_t found = 0;
+ Int_t nslices = 1;
+ Int_t org_slice;
+ for(Int_t k=0; k<fTracks->GetNTracks(); k++)
+ {
+ AliL3HoughTrack *track = (AliL3HoughTrack *)fTracks->GetCheckedTrack(k);
+ if(!track) continue;
+ // Int_t nHits = track->GetNumberOfPoints();
+ // if(nHits < fMinPointsOnTrack) break;
+ Int_t nHits = 1;
+ Int_t tracklabel;
+ tracklabel = track->GetMCid();
+
+ if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
+ if(found == 0) {
+ found=1;
+ Float_t pt=track->GetPt();
+ org_slice = track->GetSlice();
+ Float_t eta = track->GetEta();
+ Float_t phi = track->GetPsi();
+ //cout<<eta<<" "<<phi<<" "<<etag<<" "<<phig<<" "<<track->GetEtaIndex()<<endl;
+ if(tracklabel == fGoodTracks[i].label)
+ {
+ fGoodFound++;
+
+ fNFoundTracksPt->Fill(ptg);
+ fNFoundTracksEta->Fill(etag);
+ fNFoundTracksPhi->Fill(phig);
+ if((zg != 0) || (rg != 0)) {
+ fNFoundTracksRad->Fill(rg);
+ fNFoundTracksZ->Fill(zg);
+ }
+ fNtuppel->Fill(ptg,pt,nHits);
+ fPtRes->Fill((pt-ptg)/ptg*100.);
+ fEtaRes->Fill(eta-etag);
+ fPsiRes->Fill(phi-phig);
+ }
+ else
+ {
+ fNFakeTracksPt->Fill(ptg);
+ fNFakeTracksEta->Fill(etag);
+ fNFakeTracksPhi->Fill(phig);
+ }
+ //fPtRes->Fill((pt-ptg)/ptg*100.);
+ //fNtuppel->Fill(ptg,pt,nHits);
+ }
+ else {
+ found++;
+ if(track->GetSlice() != org_slice)
+ nslices = 2;
+ }
+ }
+ if(!found)
+ //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
+ continue;
+ else {
+ fFakesPt->Fill(ptg,(Double_t)found-1.0);
+ fFakesPhi->Fill(phig,(Double_t)found-1.0);
+ fSlicesPt->Fill(ptg,(Double_t)nslices);
+ fSlicesPhi->Fill(phig,(Double_t)nslices);
+ //Float_t local_phi = phi-(phi/)
+ //fFakesPtvsPhi->Fill(1.0/ptg,
+ }
+ fGoodFound2 += found;
+ }
+
+ cout<<fGoodFound<<"("<<fGoodFound2<<") tracks found out of "<<fGoodGen2<<" generated tracks"<<endl;
+ cout<<" The overall efficiency is "<<(Float_t)fGoodFound/(Float_t)fGoodGen2<<endl;
+
+ CalcEffHistos();
+ plotptres(ptmin);
+ Write2File();
+}
+
+void CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
+{
+ //Create the histograms
+
+ fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
+ fNtuppel->SetDirectory(0);
+ fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
+ fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
+ fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
+ fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
+ fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
+ fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
+
+ fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",nbin,-1,1);
+ fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",nbin,-1,1);
+ fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",nbin,-1,1);
+ fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",nbin,-1,1);
+ fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",nbin,-1,1);
+
+
+ fNGoodTracksPhi = new TH1F("fNGoodTracksPhi","Good tracks vs phi",nbin,0,6.28);
+ fNFoundTracksPhi = new TH1F("fNFoundTracksPhi","Found tracks vs phi",nbin,0,6.28);
+ fNFakeTracksPhi = new TH1F("fNFakeTracksPhi","Fake tracks vs phi",nbin,0,6.28);
+ fTrackEffPhi = new TH1F("fTrackEffPhi","Tracking efficienct vs phi",nbin,0,6.28);
+ fFakeTrackEffPhi = new TH1F("fFakeTrackEffPhi","Efficiency for fake tracks vs phi",nbin,0,6.28);
+ fFakesPt = new TProfile("fFakesPt","Number of ghosts vs pt",nbin,xlow,xup,0,100);
+ fFakesPhi = new TProfile("fFakesPhi","Number of ghosts vs phi",nbin,0,6.28,0,100);
+ fSlicesPt = new TProfile("fSlicesPt","Number of slices vs pt",nbin,xlow,xup,0,100);
+ fSlicesPhi = new TProfile("fSlicesPhi","Number of slices vs phi",nbin,0,6.28,0,100);
+ fPtRes2 = new TH1F("fPtRes2","Relative Pt resolution (Pt>2 GeV)",20,-100.,100.);
+ fEtaRes = new TH1F("fEtaRes","Eta resolution",30,-0.1,0.1);
+ fPsiRes = new TH1F("fPsiRes","Psi resolution",30,-0.1,0.1);
+
+ fNGoodTracksRad = new TH1F("fNGoodTracksRad","Good tracks vs distance(x,y)",1,0,0.2);
+ fNFoundTracksRad = new TH1F("fNFoundTracksRad","Found tracks vs distance(x,y)",1,0,0.2);
+ fRadEff = new TH1F("fRadEff","Tracking efficiency vs distance(x,y)",1,0,0.2);
+ fNGoodTracksZ = new TH1F("fNGoodTracksZ","Good tracks vs distance(z)",1,0,0.2);
+ fNFoundTracksZ = new TH1F("fNFoundTracksZ","Found tracks vs distance(z)",1,0,0.2);
+ fZEff = new TH1F("fZEff","Tracking efficiency vs distance(z)",1,0,0.2);
+}
+
+void Write2File()
+{
+ //Write histograms to file:
+
+ TFile *of = TFile::Open("hough_efficiencies.root","RECREATE");
+ if(!of->IsOpen())
+ {
+ cout<<"Problems opening rootfile"<<endl;
+ return;
+ }
+
+ of->cd();
+ fNtuppel->Write();
+ fPtRes->Write();
+ fNGoodTracksPt->Write();
+ fNFoundTracksPt->Write();
+ fNFakeTracksPt->Write();
+ fTrackEffPt->Write();
+ fFakeTrackEffPt->Write();
+ fNGoodTracksEta->Write();
+ fNFoundTracksEta->Write();
+ fNFakeTracksEta->Write();
+ fTrackEffEta->Write();
+ fFakeTrackEffEta->Write();
+
+ fNGoodTracksPhi->Write();
+ fNFoundTracksPhi->Write();
+ fNFakeTracksPhi->Write();
+ fTrackEffPhi->Write();
+ fFakeTrackEffPhi->Write();
+
+ fFakesPt->Write();
+ fFakesPhi->Write();
+ fSlicesPt->Write();
+ fSlicesPhi->Write();
+
+ fPtRes2->Write();
+ fEtaRes->Write();
+ fPsiRes->Write();
+
+ fNGoodTracksRad->Write();
+ fNFoundTracksRad->Write();
+ fRadEff->Write();
+ fNGoodTracksZ->Write();
+ fNFoundTracksZ->Write();
+ fZEff->Write();
+
+ of->Close();
+}
+
+void CalcEffHistos()
+{
+ fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
+ fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
+ fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
+ fTrackEffPt->SetMaximum(1.4);
+ fTrackEffPt->SetXTitle("P_{T} [GeV]");
+ fTrackEffPt->SetLineWidth(2);
+ fFakeTrackEffPt->SetFillStyle(3013);
+ fTrackEffPt->SetLineColor(4);
+ fFakeTrackEffPt->SetFillColor(2);
+
+ fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
+ fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
+ fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
+ fTrackEffEta->SetMaximum(1.4);
+ fTrackEffEta->SetXTitle("#eta");
+ fTrackEffEta->SetLineWidth(2);
+ fFakeTrackEffEta->SetFillStyle(3013);
+ fTrackEffEta->SetLineColor(4);
+ fFakeTrackEffEta->SetFillColor(2);
+
+ fNFoundTracksPhi->Sumw2(); fNGoodTracksPhi->Sumw2();
+ fTrackEffPhi->Divide(fNFoundTracksPhi,fNGoodTracksPhi,1,1.,"b");
+ fFakeTrackEffPhi->Divide(fNFakeTracksPhi,fNGoodTracksPhi,1,1.,"b");
+ fTrackEffPhi->SetMaximum(1.4);
+ fTrackEffPhi->SetXTitle("#psi [rad]");
+ fTrackEffPhi->SetLineWidth(2);
+ fFakeTrackEffPhi->SetFillStyle(3013);
+ fTrackEffPhi->SetLineColor(4);
+ fFakeTrackEffPhi->SetFillColor(2);
+
+ fEtaRes->SetXTitle("#eta");
+ fPsiRes->SetXTitle("#psi [rad]");
+
+ fRadEff->Divide(fNFoundTracksRad,fNGoodTracksRad,1,1.,"b");
+ fRadEff->SetXTitle("Distance [cm]");
+ fZEff->Divide(fNFoundTracksZ,fNGoodTracksZ,1,1.,"b");
+ fZEff->SetXTitle("Distance [cm]");
+}
+
+void plotptres(Float_t ptmin)
+{
+ //Plot the relative pt resolution vs pt.
+
+ const Int_t n=8;
+
+ Double_t hltx[9] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,2.0};
+ Double_t hltxerr[n];
+ Double_t hltavx[n];
+ Double_t hlty[n];
+ Double_t hltyerr[n];
+ Char_t string[1024];
+
+ for(Int_t i=0; i<n; i++)
+ {
+ TH1F *hist = new TH1F("hist","",100,-50.*(float)(i+1),50.*(float)(i+1));
+ sprintf(string,"pt_gen > %f && pt_gen <= %f",hltx[i],hltx[i+1]);
+ fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>hist",string,"goff");
+ TF1 *f1 = new TF1("f1","gaus");
+ hist->Fit("f1","E");
+ hlty[i] = f1->GetParameter(2);
+ hltyerr[i] = f1->GetParError(2);
+ //hlty[i] = hist->GetRMS();
+ //hltyerr[i] = 0;
+ hltavx[i] = (hltx[i]+hltx[i+1])/2.0;
+ hltxerr[i] = (hltx[i+1]-hltx[i])/2.0;
+ cout<<string<<" "<<i<<" "<<hlty[i]<<" "<<hltyerr[i]<<endl;
+ delete f1;
+ delete hist;
+ }
+
+ sprintf(string,"pt_gen > %f",2.0);
+ fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>fPtRes2",string,"goff");
+ TF1 *f1 = new TF1("f1","gaus");
+ fPtRes2->Fit("f1","E");
+ fPtRes2->SetXTitle("#Delta P_{t} / P_{t} [%]");
+ fPtRes2->SetYTitle("#Delta N");
+ delete f1;
+
+ TGraphErrors *gr1 = new TGraphErrors(n,hltavx,hlty,hltxerr,hltyerr);
+ TCanvas *c1 = new TCanvas("c1","",1);
+ fPtRes = c1->DrawFrame(0.1,0,2,15);
+ gr1->SetTitle("");
+ gr1->Draw("P");
+ gr1->SetMarkerStyle(20);
+ gr1->SetLineWidth(2);
+ gr1->SetMarkerSize(1.3);
+ f1 = new TF1("f1","pol1",ptmin,2.0);
+ gr1->Fit("f1","ER");
+ fPtRes->SetXTitle("p_{t} [GeV]");
+ fPtRes->SetYTitle("#Delta P_{t} / P_{t} [%]");
+ delete f1;
+}
and written to file.
*/
-void evaltracker(char *path,char *trackpath,char *offlinepath,int nevent=1)
+#ifndef __CINT__
+#include "AliL3Logger.h"
+#include "AliL3Evaluate.h"
+#include "AliL3FileHandler.h"
+#include "AliL3DigitData.h"
+#include "AliL3Transform.h"
+#include "AliLevel3.h"
+#include <TROOT.h>
+#include <TNtuple.h>
+#include <TRandom.h>
+#include <TSystem.h>
+#include <TStyle.h>
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TCanvas.h>
+#include <TF1.h>
+#include <stdio.h>
+#include <iostream.h>
+#include <time.h>
+#endif
+
+void evaltracker(Char_t *path="./",Char_t *trackpath="./tracker/",char *offlinepath="./",int nevent=1)
{
//Make sure you got the correct parameters:
AliL3Transform::Init(path,kTRUE);
//Define which slices to include:
- int slicerange[2]={0,35};
+ Int_t slicerange[2]={0,35};
//Define which padrows to include (should normally be all)
- int rowrange[2]={AliL3Transform::GetFirstRow(-1),AliL3Transform::GetLastRow(-1)};
+ //Int_t rowrange[2]={AliL3Transform::GetFirstRow(-1),AliL3Transform::GetLastRow(-1)};
//Define the minimum number of clusters on a simulated track to be found:
- int good = (int)(0.4*AliL3Transform::GetNRows());
+ Int_t good = (Int_t)(0.4*AliL3Transform::GetNRows());
//Define the minumum number of clusters on a found track to be found (should normally be the same as above)
- int nclusters = (int)(0.4*AliL3Transform::GetNRows());
+ Int_t nclusters = (Int_t)(0.4*AliL3Transform::GetNRows());
//Define which pt range to include
- float ptmin = 0.1;
- float ptmax = 4.;
+ Float_t ptmin = 0.1;
+ Float_t ptmax = 4.;
//Define the maximum ratio of false clusters which are allowed on a found track (default=0.1 -> 10%)
- float maxfalseratio = 0.1;
-
- a = new AliL3Evaluate(trackpath,nclusters,good,ptmin,ptmax,slicerange);
+ Float_t maxfalseratio = 0.1;
+ AliL3Evaluate *a = new AliL3Evaluate(trackpath,nclusters,good,ptmin,ptmax,slicerange);
a->CreateHistos(20,0.1,4.1);//(nbins in pt,minpt,maxpt) -> the same as used by standard offline
a->SetMaxFalseClusters(maxfalseratio);
//a->SetStandardComparison(kFALSE); //use AliTPCComparison_HLT
//Loop over the number of events:
- for(int event=0; event<nevent; event++)
+ for(Int_t event=0; event<nevent; event++)
{
cout<<"Processing event "<<event<<endl;
}
-void ploteff(char *rootfile="hlt_efficiencies.root")
+void ploteff(Char_t *rootfile="hlt_efficiencies.root")
{
//Plot the efficiency vs pt.
gStyle->SetOptStat(0);
- char filename[1024];
-
- double hltx[22];
- double hlty[22];
- double hltxerr[22];
- double hltyerr[22];
- int i;
- file = TFile::Open(rootfile);
-
- int hltcounter=0;
- h = (TH1*)fTrackEffPt;
- for(i=0; i<22; i++)
+ Double_t hltx[22];
+ Double_t hlty[22];
+ Double_t hltxerr[22];
+ Double_t hltyerr[22];
+ TFile *file = TFile::Open(rootfile);
+ Int_t hltcounter=0;
+ TH1* h = (TH1*)file->Get("fTrackEffPt");
+ for(Int_t i=0; i<22; i++)
{
if(h->GetBinContent(i)==0) continue;
hltx[hltcounter] = h->GetBinCenter(i);
}
file->Close();
- gr1 = new TGraphErrors(hltcounter,hltx,hlty,hltxerr,hltyerr);
+ TGraphErrors *gr1 = new TGraphErrors(hltcounter,hltx,hlty,hltxerr,hltyerr);
- c1 = new TCanvas("c1","",1);
- hist = c1->DrawFrame(0.1,0,3,1.4);
+ TCanvas *c1 = new TCanvas("c1","",1);
+ TH1F *hist = c1->DrawFrame(0.1,0,3,1.4);
c1->SetGridx();
c1->SetGridy();
gr1->SetTitle("");
gr1->SetMarkerStyle(20);
gr1->SetLineWidth(2);
gr1->SetMarkerSize(1.3);
-
}
-void plotptres(char *rootfile="hlt_efficiencies.root")
+void plotptres(Char_t *rootfile="hlt_efficiencies.root")
{
//Plot the relative pt resolution vs pt.
- const int n=5;
+ const Int_t n=15;
- double hltx[15] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0};
- double hlty[n];
- double hltyerr[n];
- int i;
+ Double_t hltx[15] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0};
+ Double_t hlty[n];
+ Double_t hltyerr[n];
+ Char_t string[1024];
- char string[1024];
+ TFile *file = TFile::Open(rootfile);
+ TH1F *hist = new TH1F("hist","",100,-10,10);
+ TNtuple *fNtuppel=(TNtuple*)file->Get("fNtuppel");
- file = TFile::Open(rootfile);
-
- hist = new TH1F("hist","",100,-10,10);
-
- for(i=0; i<n; i++)
+ for(Int_t i=0; i<n; i++)
{
sprintf(string,"pt_gen > %f && pt_gen <= %f && nHits>63",hltx[i]-0.1,hltx[i]+0.1);
fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>hist",string,"goff");
- f1 = new TF1("f1","gaus");
+ TF1 *f1 = new TF1("f1","gaus");
hist->Fit("f1","QN");
hlty[i] = f1->GetParameter(2);
hltyerr[i] = f1->GetParError(2);
delete hist;
file->Close();
- gr1 = new TGraphErrors(n,hltx,hlty,0,hltyerr);
-
- c1 = new TCanvas("c1","",1);
- hist = c1->DrawFrame(0.1,0,3,5);
+ TGraphErrors *gr1 = new TGraphErrors(n,hltx,hlty,0,hltyerr);
+ TCanvas *c1 = new TCanvas("c1","",1);
+ hist = c1->DrawFrame(0.1,0,3,15);
c1->SetGridx();
c1->SetGridy();
gr1->SetTitle("");
gr1->SetMarkerSize(1.3);
hist->SetXTitle("p_{t} [GeV]");
hist->SetYTitle("#Delta P_{T} / P_{T} [%]");
-
}
gSystem->Load("MakePileup_C.so");
gSystem->Load("Read_C.so");
gSystem->Load("runhough_C.so");
+ gSystem->Load("runrowhough_C.so");
gSystem->Load("deconvclusters_C.so");
gSystem->Load("runtracker_C.so");
gErrorIgnoreLevel=saveErrIgLevel;
--- /dev/null
+// $Id$
+
+#ifndef __CINT__
+#include "AliL3Logger.h"
+#include "AliL3FileHandler.h"
+#include "AliL3DigitData.h"
+#include "AliL3Transform.h"
+#include "AliL3Hough.h"
+#include "AliL3TrackArray.h"
+#include "AliL3Track.h"
+#include "AliL3HoughTrack.h"
+#include "AliL3Fitter.h"
+#include "AliL3ClusterFitter.h"
+#include "AliL3Vertex.h"
+#include "AliL3Benchmark.h"
+#include <AliRunLoader.h>
+#include <AliStack.h>
+#include <TParticle.h>
+#include <TNtuple.h>
+#include <TRandom.h>
+#include <TSystem.h>
+#include <TStopwatch.h>
+#include <TBenchmark.h>
+#include <stdio.h>
+#include <iostream.h>
+#include <time.h>
+#endif
+
+Int_t runrowhough(Char_t *path="./",Char_t *outpath="./fitter",int s1=0,int s2=35,int nevent=1,Bool_t skip=kTRUE)
+{
+ Bool_t isinit=AliL3Transform::Init(path,kTRUE);
+ if(!isinit){
+ cerr << "Could not create transform settings, please check log for error messages!" << endl;
+ return 1;
+ }
+ Float_t ptmin = 0.1*AliL3Transform::GetSolenoidField();
+ Float_t zvertex;
+
+ {
+ AliRunLoader *rl = AliRunLoader::Open("galice.root");
+ rl->LoadHeader();
+ rl->LoadKinematics();
+ AliStack* stack = rl->Stack();
+ TParticle *orig = (TParticle*)stack->Particle(0);
+ Float_t xori = orig->Vx();
+ Float_t yori = orig->Vy();
+ zvertex = orig->Vz();
+ cout<<" Primary vertex at ("<<xori<<","<<yori<<","<<zvertex<<")"<<endl;
+ if (rl->LoadgAlice()) {
+ cerr<<"Error occured while loading gAlice"<<endl;
+ return 1;
+ }
+ delete rl;
+ }
+
+ cout<<" Hough Tranform will run with ptmin="<<ptmin<<" and zvertex="<<zvertex<<endl;
+
+ AliL3Benchmark *fBenchmark = new AliL3Benchmark();
+ AliL3Hough *hough = new AliL3Hough();
+ hough->SetThreshold(4);
+ hough->SetTransformerParams(140,76,ptmin,-1);
+ hough->SetPeakThreshold(50,-1);
+ hough->Init(path, kFALSE, 100, kFALSE,4,0,0,zvertex);
+ hough->SetAddHistograms();
+
+ for(int ev=0; ev<nevent; ev++)
+ {
+ for(int slice=s1; slice<=s2; slice++)
+ {
+ cout<<"Processing slice "<<slice<<endl;
+ hough->ReadData(slice,ev);
+ hough->Transform();
+ hough->AddAllHistogramsRows();
+ hough->FindTrackCandidates();
+ hough->AddTracks();
+ }
+ hough->WriteTracks(outpath);
+ Char_t bname[100];
+ sprintf(bname,"rowhough_%d",ev);
+ hough->DoBench(bname);
+
+ if(!skip) {
+ // Run cluster fitter
+ AliL3ClusterFitter *fitter = new AliL3ClusterFitter(path);
+
+ // Set debug flag for the cluster fitter
+ // fitter->Debug();
+
+ // Setting fitter parameters
+ fitter->SetInnerWidthFactor(1,1.5);
+ fitter->SetOuterWidthFactor(1,1.5);
+ fitter->SetNmaxOverlaps(5);
+
+ // fitter->SetChiSqMax(5,kFALSE); //isolated clusters
+ fitter->SetChiSqMax(5,kTRUE); //overlapping clusters
+
+ Int_t rowrange[2] = {0,AliL3Transform::GetNRows()-1};
+
+ // Takes input from global hough tracks produced by HT
+ fitter->LoadSeeds(rowrange,kFALSE,ev,zvertex);
+
+ UInt_t ndigits;
+
+ for(int slice=s1; slice<=s2; slice++)
+ {
+ for(Int_t ipatch = 0; ipatch < AliL3Transform::GetNPatches(); ipatch++)
+ {
+ // Read digits
+ hough->GetMemHandler(ipatch)->Free();
+ hough->GetMemHandler(ipatch)->Init(slice,ipatch);
+ AliL3DigitRowData *digits = (AliL3DigitRowData *)hough->GetMemHandler(ipatch)->AliAltroDigits2Memory(ndigits,ev);
+
+ fBenchmark->Start("Fitter Init");
+ fitter->Init(slice,ipatch);
+ fBenchmark->Stop("Fitter Init");
+ fitter->SetInputData(digits);
+ fBenchmark->Start("Fitter cluster finder");
+ fitter->FindClusters();
+ fBenchmark->Stop("Fitter cluster finder");
+ fitter->WriteClusters();
+ }
+ }
+
+ // Refit of the clusters
+ AliL3Vertex vertex;
+ //The seeds are the input tracks from circle HT
+ AliL3TrackArray *tracks = fitter->GetSeeds();
+ AliL3Fitter *ft = new AliL3Fitter(&vertex,1);
+
+ ft->LoadClusters("./fitter/",ev,kFALSE);
+ fBenchmark->Start("Track fitter");
+ for(Int_t i=0; i<tracks->GetNTracks(); i++)
+ {
+ AliL3Track *track = tracks->GetCheckedTrack(i);
+ if(!track) continue;
+ if(track->GetNHits() < 20) continue;
+ ft->SortTrackClusters(track);
+ ft->FitHelix(track);
+ track->UpdateToFirstPoint();
+ }
+ fBenchmark->Stop("Track fitter");
+ delete ft;
+
+ //Write the final tracks
+ fitter->WriteTracks(20);
+
+ delete fitter;
+ }
+ }
+
+ if(!skip) {
+ fBenchmark->Analyze("fitter");
+ }
+
+ hough->DoBench("rowhough");
+ delete hough;
+ delete fBenchmark;
+ return 0;
+}
//a->DoMc(); /*do monte carlo identification*/
a->WriteFiles(opath); /*enable output*/
a->ProcessEvent(minslice,maxslice);
- a->DoBench("benchmark_0");
-
+ Char_t bname[100];
+ sprintf(bname,"benchmark_tracker_%d",ev);
+ a->DoBench(bname);
delete a;
} // event loop
}
#include "AliL3Logging.h"
#include "AliL3Histogram.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** \class AliL3Histogram
+<pre>
//_____________________________________________________________
// AliL3Histogram
//
// 2D histogram class
+//
+</pre>
+*/
+
+//uncomment if you want overflow checks
+//#define _IFON_
ClassImp(AliL3Histogram)
return fYmin + (ybin-0.5) * fBinwidthY;
}
+Double_t AliL3Histogram::GetPreciseBinCenterX(Float_t xbin) const
+{
+ if(xbin < fFirstXbin || xbin > fLastXbin)
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterX","xbin")
+ <<"Bin-value out of range "<<xbin<<ENDLOG;
+ return -1;
+ }
+ // return fXmin + (xbin-1) * fBinwidthX + 0.5*fBinwidthX;
+ return fXmin + (xbin-0.5) * fBinwidthX;
+}
+
+Double_t AliL3Histogram::GetPreciseBinCenterY(Float_t ybin) const
+{
+ if(ybin < fFirstYbin || ybin > fLastYbin)
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterY","ybin")
+ <<"Bin-value out of range "<<ybin<<ENDLOG;
+ return -1;
+ }
+ // return fYmin + (ybin-1) * fBinwidthY + 0.5*fBinwidthY;
+ return fYmin + (ybin-0.5) * fBinwidthY;
+}
+
void AliL3Histogram::Draw(Char_t *option)
{
#ifdef use_root
Double_t GetYmax() const {return fYmax;}
virtual Double_t GetBinCenterX(Int_t xbin) const;
virtual Double_t GetBinCenterY(Int_t ybin) const;
+ Double_t GetPreciseBinCenterX(Float_t xbin) const;
+ Double_t GetPreciseBinCenterY(Float_t ybin) const;
Double_t GetBinWidthX() const {return fBinwidthX;}
Double_t GetBinWidthY() const {return fBinwidthY;}
Int_t GetFirstXbin() const {return fFirstXbin;}
#include "AliL3Logging.h"
#include "AliL3Histogram1D.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3Transform.h"
#include "AliL3Track.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3HoughClusterTransformer.h"
#include "AliL3HoughTransformerLUT.h"
#include "AliL3HoughTransformerVhdl.h"
-#include "AliL3HoughTransformerGap.h"
+#include "AliL3HoughTransformerRow.h"
#include "AliL3HoughMaxFinder.h"
#include "AliL3Benchmark.h"
#ifdef use_aliroot
#include "AliL3HoughTrack.h"
#include "AliL3DDLDataFileHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
// hough->FindTrackCandidates();
//
// AliL3TrackArray *tracks = hough->GetTracks(patch);
+//
//</pre>
*/
fCurrentSlice = 0;
fEvent = 0;
- fKappaSpread=6;
- fPeakRatio=0.5;
+ fKappaSpread = 6;
+ fPeakRatio = 0.5;
+ fInputFile = 0;
+ fInputPtr = 0;
SetTransformerParams();
SetThreshold();
#endif
}
-AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile)
+AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
{
fBinary = binary;
strcpy(fPath,path);
fVersion = tv;
fKappaSpread=6;
fPeakRatio=0.5;
- if(!fBinary)
- fInputFile = infile;
- else
+ if(!fBinary) {
+ if(infile) {
+ fInputFile = infile;
+ fInputPtr = 0;
+ }
+ else {
+ fInputFile = 0;
+ fInputPtr = ptr;
+ }
+ }
+ else {
fInputFile = 0;
-
+ fInputPtr = 0;
+ }
#ifdef use_aliroot
//just be sure that index is empty for new event
AliL3FileHandler::CleanStaticIndex();
//cout << "Cleaned class mem " << endl;
}
-void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile)
+void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
{
fBinary = binary;
strcpy(fPath,path);
fWriteDigits = kFALSE;
fUse8bits = bit8;
fVersion = tv;
- if(!fBinary)
- fInputFile = infile;
- else
+ if(!fBinary) {
+ if(infile) {
+ fInputFile = infile;
+ fInputPtr = 0;
+ }
+ else {
+ fInputFile = 0;
+ fInputPtr = ptr;
+ }
+ }
+ else {
fInputFile = 0;
+ fInputPtr = 0;
+ }
+ fZVertex = zvertex;
Init(); //do the rest
}
fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
break;
case 4:
- fHoughTransformer[i] = new AliL3HoughTransformerGap(0,i,fNEtaSegments);
+ fHoughTransformer[i] = new AliL3HoughTransformerRow(0,i,fNEtaSegments,kFALSE,fZVertex);
break;
default:
fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments,kFALSE,kFALSE);
#ifdef use_aliroot
{
if(!fInputFile) {
- /* In case of reading digits file */
- fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
- if(!fBinary) {
- Char_t filename[1024];
- sprintf(filename,"%s/digitfile.root",fPath);
- fMemHandler[i]->SetAliInput(filename);
+ if(!fInputPtr) {
+ /* In case of reading digits file */
+ fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
+ if(!fBinary) {
+ Char_t filename[1024];
+ sprintf(filename,"%s/digitfile.root",fPath);
+ fMemHandler[i]->SetAliInput(filename);
+ }
+ }
+ else {
+ /* In case of reading from DATE */
+ fMemHandler[i] = new AliL3DDLDataFileHandler();
+ fMemHandler[i]->SetReaderInput(fInputPtr,-1);
}
}
else {
{
ReadData(i);
Transform();
- if(fAddHistograms)
- AddAllHistograms();
+ if(fAddHistograms) {
+ if(fVersion != 4)
+ AddAllHistograms();
+ else
+ AddAllHistogramsRows();
+ }
FindTrackCandidates();
//Evaluate();
//fGlobalMerger->FillTracks(fTracks[0],i);
initTime = GetCpuTime();
for(Int_t i=0; i<fNPatches; i++)
{
- fHoughTransformer[i]->Reset();//Reset the histograms
+ // In case of Row transformer reset the arrays only once
+ if((fVersion != 4) || (i == 0))
+ fHoughTransformer[i]->Reset();//Reset the histograms
fBenchmark->Start("Hough Transform");
if(!row_range)
fHoughTransformer[i]->TransformCircle();
<<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
}
+void AliL3Hough::AddAllHistogramsRows()
+{
+ //Add the histograms within one etaslice.
+ //Resulting histogram are in patch=0.
+
+ Double_t initTime,cpuTime;
+ initTime = GetCpuTime();
+ fBenchmark->Start("Add HistogramsRows");
+
+ for(Int_t i=0; i<fNEtaSegments; i++)
+ {
+ UChar_t *rowcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetRowCount(i);
+ UChar_t *gapcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
+
+ AliL3Histogram *hist = fHoughTransformer[0]->GetHistogram(i);
+ Int_t xmin = hist->GetFirstXbin();
+ Int_t xmax = hist->GetLastXbin();
+ Int_t ymin = hist->GetFirstYbin();
+ Int_t ymax = hist->GetLastYbin();
+
+ for(Int_t ybin=ymin; ybin<=ymax; ybin++)
+ {
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ Int_t bin = hist->GetBin(xbin,ybin);
+ if(gapcount[bin] < 3)
+ // hist->AddBinContent(bin,rowcount[bin]);
+ hist->AddBinContent(bin,rowcount[bin]/gapcount[bin]);
+ }
+ }
+ }
+
+ fBenchmark->Stop("Add HistogramsRows");
+ fAddHistograms = kTRUE;
+ cpuTime = GetCpuTime() - initTime;
+ LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistogramsRows()","Timing")
+ <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
+}
+
void AliL3Hough::AddTracks()
{
if(!fTracks[0])
fPeakFinder->SetHistogram(hist);
fPeakFinder->SetThreshold(fPeakThreshold[i]);
- fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
-
- //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
+ if(fVersion != 4) {
+ fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
+ //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
+ }
+ else /* row transformer needs different peak finder method */
+ fPeakFinder->FindAdaptedRowPeaks(1,2,0);
for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
{
track->SetEtaIndex(j);
track->SetEta(tr->GetEta(j,fCurrentSlice));
track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
+#ifdef do_mc
+ Int_t label = tr->GetTrackID(j,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
}
}
cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
private:
Char_t *fInputFile;//!
-
+ Char_t *fInputPtr;//!
Char_t fPath[1024];
Bool_t fBinary;
Bool_t fAddHistograms;
Int_t fKappaSpread;
Float_t fPeakRatio;
+ Float_t fZVertex;
+
AliL3MemHandler **fMemHandler; //!
AliL3HoughBaseTransformer **fHoughTransformer; //!
AliL3HoughEval **fEval; //!
public:
AliL3Hough();
- AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0);
+ AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0);
virtual ~AliL3Hough();
- void Init(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0);
+ void Init(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0,Float_t zvertex=0.0);
void Init(Bool_t doit=kFALSE, Bool_t addhists=kFALSE);
void Process(Int_t minslice,Int_t maxslice);
void FindTrackCandidates();
void AddAllHistograms();
+ void AddAllHistogramsRows();
Int_t Evaluate(Int_t road_width=1,Int_t nrowstomiss=1);
void EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss);
void WriteTracks(Int_t slice,Char_t *path="./");
#include "AliL3Histogram.h"
#include "AliL3HoughBaseTransformer.h"
-
+/** \class AliL3HoughBaseTransformer
+<pre>
//_____________________________________________________________
// AliL3HoughBaseTransformer
//
//
// This is an abstract class, and is only meant to provide the interface
// to the different implementations.
+//
+</pre>
+*/
ClassImp(AliL3HoughBaseTransformer)
fEtaMax = 0;
fLowerThreshold = 0;
fUpperThreshold = 1023;
+ fZVertex = 0.0;
}
-AliL3HoughBaseTransformer::AliL3HoughBaseTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments)
+AliL3HoughBaseTransformer::AliL3HoughBaseTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Float_t zvertex)
{
fDigitRowData = 0;
fEtaMax = 0;
fLowerThreshold = 3;
fUpperThreshold = 1023;
+ fZVertex = zvertex;
Init(slice,patch,n_eta_segments);
}
#include "AliL3RootTypes.h"
#ifdef do_mc
-const UInt_t MaxTrack=35;
+const UInt_t MaxTrack=100;
struct TrackIndex {
Int_t fLabel[MaxTrack];
+#ifdef ROWHOUGH
+ UChar_t fNHits[MaxTrack];
+ UChar_t fCurrentRow[MaxTrack];
+#else
Int_t fNHits[MaxTrack];
+#endif
};
typedef struct TrackIndex TrackIndex;
#endif
Int_t fUpperThreshold;
AliL3DigitRowData *fDigitRowData; //!
-
+
+ Float_t fZVertex;
+
public:
AliL3HoughBaseTransformer();
- AliL3HoughBaseTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments);
+ AliL3HoughBaseTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Float_t zvertex=0.0);
virtual ~AliL3HoughBaseTransformer();
void SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr) {fDigitRowData = ptr;}
Int_t GetUpperThreshold() {return fUpperThreshold;}
Double_t GetEtaMin() {return fEtaMin;}
Double_t GetEtaMax() {return fEtaMax;}
-
+ Float_t GetZVertex() {return fZVertex;}
+
AliL3DigitRowData *GetDataPointer() {return fDigitRowData;}
virtual Int_t GetEtaIndex(Double_t eta) = 0;
#include "AliL3Histogram.h"
#include "AliL3ClustFinderNew.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3MemHandler.h"
#include "AliL3HoughDisplay.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3Histogram.h"
#include "AliL3Histogram1D.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** /class AliL3HoughEval
+//<pre>
//_____________________________________________________________
// AliL3HoughEval
//
// Evaluation class for tracklets produced by the Hough transform.
+//
+</pre>
+*/
ClassImp(AliL3HoughEval)
fNEtaSegments = fHoughTransformer->GetNEtaSegments();
fEtaMin = fHoughTransformer->GetEtaMin();
fEtaMax = fHoughTransformer->GetEtaMax();
+ fZVertex = fHoughTransformer->GetZVertex();
GenerateLUT();
}
Int_t sector,row;
AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
+ xyz[2] -= fZVertex;
Double_t eta = AliL3Transform::GetEta(xyz);
Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
if(pixel_index != eta_index) continue;
Int_t fNumOfPadsToLook;
Int_t fNumOfRowsToMiss;
AliL3Histogram1D **fEtaHistos; //!
+ Float_t fZVertex;
//Flags
Bool_t fRemoveFoundTracks;
void SetNumOfRowsToMiss(Int_t i) {fNumOfRowsToMiss = i;}
void SetNumOfPadsToLook(Int_t i) {fNumOfPadsToLook = i;}
void SetSlice(Int_t i) {fSlice=i;}
+ void SetZVertex(Float_t zvertex) {fZVertex=zvertex;}
ClassDef(AliL3HoughEval,1) //Hough transform verfication class
#include "AliL3Transform.h"
#include "AliL3TrackArray.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#pragma link C++ class AliL3HoughTransformerLUT;
#pragma link C++ class AliL3HoughTransformerVhdl;
#pragma link C++ class AliL3HoughTransformerNew;
-#pragma link C++ class AliL3HoughTransformerGap;
+#pragma link C++ class AliL3HoughTransformerRow;
#ifndef Darwin
#pragma link C++ class AliL3HoughTrack;
#endif
#endif
#endif
-
#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** \class AliL3HoughMaxFinder
+<pre>
//_____________________________________________________________
// AliL3HoughMaxFinder
//
// Maximum finder
+//
+</pre>
+*/
ClassImp(AliL3HoughMaxFinder)
#endif
}
-
AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
{
//Constructor
fThreshold=0;
}
-
AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
{
//Destructor
value[bin_index]=hist->GetBinContent(bin[bin_index]);
bin_index++;
}
-
}
if(value[12]==0) continue;
Int_t b=0;
if(hist->GetNEntries() == 0)
return;
-
+
Int_t xmin = hist->GetFirstXbin();
Int_t xmax = hist->GetLastXbin();
Int_t ymin = hist->GetFirstYbin();
Int_t ymax = hist->GetLastYbin();
-
+
+
//Start by looking for pre-peaks:
Window **local_maxima = new Window*[hist->GetNbinsY()];
delete [] maxs;
}
+struct PreYPeak
+{
+ Int_t start_position;
+ Int_t end_position;
+ Int_t value;
+ Int_t prev_value;
+};
+
+struct Pre2DPeak
+{
+ Int_t start_x_position;
+ Int_t end_x_position;
+ Int_t start_y_position;
+ Int_t end_y_position;
+ Int_t value;
+};
+
+void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
+{
+
+ AliL3Histogram *hist = fCurrentHisto;
+
+ if(!hist)
+ {
+ cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
+ return;
+ }
+
+ if(hist->GetNEntries() == 0)
+ return;
+
+ Int_t xmin = hist->GetFirstXbin();
+ Int_t xmax = hist->GetLastXbin();
+ Int_t ymin = hist->GetFirstYbin();
+ Int_t ymax = hist->GetLastYbin();
+
+ //Start by looking for pre-peaks:
+
+ PreYPeak **local_maxima = new PreYPeak*[hist->GetNbinsY()];
+
+ Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
+ Int_t n2dmax = 0;
+ Int_t last_value,value;
+ for(Int_t ybin=ymin; ybin<=ymax; ybin++)
+ {
+ local_maxima[ybin-ymin] = new PreYPeak[hist->GetNbinsX()];
+ nmaxs[ybin-ymin] = 0;
+ last_value = 0;
+ Bool_t found = 0;
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ value = hist->GetBinContent(hist->GetBin(xbin,ybin));
+ if(value >= fThreshold)
+ {
+ if(value > last_value)
+ // if(value > (last_value + 1))
+ {
+ local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start_position = xbin;
+ local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
+ local_maxima[ybin-ymin][nmaxs[ybin-ymin]].value = value;
+ local_maxima[ybin-ymin][nmaxs[ybin-ymin]].prev_value = 0;
+ found = 1;
+ }
+ if(value == last_value)
+ // if(abs(value - last_value) <= 1)
+ if(found)
+ local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
+ }
+
+ if((value < fThreshold) || (value < last_value))
+ // if((value < fThreshold) || (value < (last_value - 1)))
+ {
+ if(found)
+ {
+ nmaxs[ybin-ymin]++;
+ found = 0;
+ }
+ }
+ last_value = value;
+
+ }
+ n2dmax += nmaxs[ybin-ymin];
+ }
+
+ // Pre2DPeak *maxima = new Pre2DPeak[n2dmax];
+ for(Int_t ybin=ymax; ybin >= ymin; ybin--)
+ {
+ for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
+ {
+ Int_t local_value = local_maxima[ybin-ymin][i].value;
+ Int_t local_prev_value = local_maxima[ybin-ymin][i].prev_value;
+ Int_t local_next_value = 0;
+
+ if(local_value<0)
+ continue; //already used
+
+ //start expanding in the psi-direction:
+
+ Int_t local_x_start = local_maxima[ybin-ymin][i].start_position;
+ Int_t local_x_end = local_maxima[ybin-ymin][i].end_position;
+ Int_t temp_x_start = local_maxima[ybin-ymin][i].start_position;
+ Int_t temp_x_end = local_maxima[ybin-ymin][i].end_position;
+
+ Int_t local_y=ybin-1,nybins=1;
+
+ while(local_y >= ymin)
+ {
+ Bool_t found=0;
+ for(Int_t j=0; j<nmaxs[local_y-ymin]; j++)
+ {
+ Int_t adapted_kappawindow;
+ Float_t adapted_x,adapted_y;
+ adapted_x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0;
+ adapted_y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0;
+ // if(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56) && adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))
+ // if(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))
+ // if((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)))
+ if(((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))) ||
+ ((adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))&& !(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56))))
+ adapted_kappawindow = 1;
+ else
+ adapted_kappawindow = 0;
+ //temprorary here
+ // adapted_kappawindow = 0;
+
+ if( (local_maxima[local_y-ymin][j].start_position <= (temp_x_end + kappawindow + adapted_kappawindow)) && (local_maxima[local_y-ymin][j].end_position >= temp_x_start))
+ {
+ if( local_maxima[local_y-ymin][j].value == local_value )
+ {
+ local_x_end = local_maxima[local_y-ymin][j].end_position;
+ temp_x_start = local_maxima[local_y-ymin][j].start_position;
+ temp_x_end = local_maxima[local_y-ymin][j].end_position;
+ local_maxima[local_y-ymin][j].value = -1;
+ found = 1;
+ nybins++;
+ break;
+ }
+ else
+ {
+ local_maxima[local_y-ymin][j].prev_value = local_value;
+ local_next_value = local_maxima[local_y-ymin][j].value;
+ }
+ }
+ }
+ if(!found || local_y == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
+ {
+ Int_t adapted_xsize;
+ Float_t adapted_x,adapted_y;
+ adapted_x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0;
+ adapted_y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0;
+ // if(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56) && adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))
+ // if(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))
+ if((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)))
+ adapted_xsize = 1;
+ else
+ adapted_xsize = 2;
+ //temprorary here
+ // adapted_xsize = 1;
+ if((nybins > ysize) && ((local_x_end-local_x_start+1) > adapted_xsize) && (local_value > local_prev_value) && (local_value > local_next_value))
+ {
+ fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(adapted_x);
+ fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(adapted_y);
+ fWeight[fNPeaks] = local_value;
+#ifdef do_mc
+ cout<<"Peak found at: "<<((Float_t)local_x_start+(Float_t)local_x_end)/2.0<<" "<<((Float_t)ybin+(Float_t)(local_y+1))/2.0<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<local_x_end-local_x_start+1<<" by "<<nybins<<endl;
+#endif
+ fNPeaks++;
+ }
+ break;
+ }
+ else
+ local_y--;//Search continues...
+ }
+ }
+ }
+
+ for(Int_t i=0; i<hist->GetNbinsY(); i++)
+ delete local_maxima[i];
+
+ delete [] local_maxima;
+ delete [] nmaxs;
+ //delete [] maxima;
+}
void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
{
delete windowPt[i];
delete [] windowPt;
delete [] anotherPt;
-
}
void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
if(a->weight < b->weight) return 1;
if(a->weight > b->weight) return -1;
return 0;
-
}
void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
delete [] m_up;
//return track;
-
-
}
void FindBigMaxima();
void FindMaxima(Int_t threshold=0);
void FindAdaptedPeaks(Int_t nkappawindow,Float_t cut_ratio);
+ //Peak finder for HoughTransformerRow
+ void FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize);
//More sophisticated peak finders:
void FindPeak(Int_t t1,Double_t t2,Int_t t3);
void FindPeak1(Int_t y_window=2,Int_t x_bin_sides=1);
#include "AliL3HoughMerger.h"
#include "AliL3HoughTransformer.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include <TH3.h>
#include <TAxis.h>
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3HoughTrack.h"
#include "AliL3Transform.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+/** \class AliL3HoughTrack
+<pre>
//_____________________________________________________________
// AliL3HoughTrack
//
// Track class for Hough tracklets
+//
+</pre>
+*/
ClassImp(AliL3HoughTrack)
fEta = 0;
}
-
AliL3HoughTrack::~AliL3HoughTrack()
{
}
void AliL3HoughTrack::Set(AliL3Track *track)
{
-
AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
SetEtaIndex(tpt->GetEtaIndex());
SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
SetSlice(tpt->GetSlice());
SetHits(tpt->GetNHits(),(UInt_t *)tpt->GetHitNumbers());
- return;
-
- fWeight = tpt->GetWeight();
- fDLine = tpt->GetDLine();
- fPsiLine = tpt->GetPsiLine();
- SetNHits(tpt->GetWeight());
- SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
- fIsHelix = false;
+#ifdef ROWHOUGH
+ SetMCid(tpt->GetMCid());
+#endif
+ /*
+ fWeight = tpt->GetWeight();
+ fDLine = tpt->GetDLine();
+ fPsiLine = tpt->GetPsiLine();
+ SetNHits(tpt->GetWeight());
+ SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
+ fIsHelix = false;
+ */
}
Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
SetTgl(tgl);
}
-
void AliL3HoughTrack::UpdateToFirstRow()
{
//Update the track parameters to the point where track cross
void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
{
-
fWeight = weight;
fMinDist = 100000;
SetKappa(kappa);
SetNHits(1);
SetRowRange(rowrange[0],rowrange[1]);
fIsHelix = false;
-
}
void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t min_dist)
{
-
if(min_dist < fMinDist)
{
fMinDist = min_dist;
SetMCid(mcid);
}
-
}
void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Float_t *xy)
{
-
if(fIsHelix)
{
printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
Float_t yhit = a*xhit + b;
xy[0] = xhit;
xy[1] = yhit;
-
}
#include "AliL3DigitData.h"
#include "AliL3HistogramAdaptive.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
+++ /dev/null
-// @(#) $Id$
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include "AliL3StandardIncludes.h"
-
-#include "AliL3Logging.h"
-#include "AliL3MemHandler.h"
-#include "AliL3Transform.h"
-#include "AliL3DigitData.h"
-#include "AliL3HistogramAdaptive.h"
-#include "AliL3HoughTransformerGap.h"
-
-#if GCCVERSION == 3
-using namespace std;
-#endif
-
-//_____________________________________________________________
-// AliL3HoughTransformerGap
-//
-// Hough transformation class
-//
-
-ClassImp(AliL3HoughTransformerGap)
-
-#ifdef COUNTINGROWS
-UChar_t **AliL3HoughTransformerGap::fRowCount = 0;
-UChar_t **AliL3HoughTransformerGap::fGapCount = 0;
-UChar_t **AliL3HoughTransformerGap::fCurrentRowCount = 0;
-#endif
-
-AliL3HoughTransformerGap::AliL3HoughTransformerGap()
-{
- //Default constructor
- fParamSpace = 0;
- fDoMC = kFALSE;;
-
-#ifdef do_mc
- fTrackID = 0;
-#endif
- fZVertex = 0.;
-#ifdef SUBSLICES
- fEtaSpace = 0;
- fEtaSpaceSize = 0;
-#endif
-#ifdef COUNTINGROWS
- fLUT2sinphi0up=0;
- fLUT2sinphi0low=0;
- fLUT2cosphi0up=0;
- fLUT2cosphi0low=0;
-#endif
-}
-
-AliL3HoughTransformerGap::AliL3HoughTransformerGap(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap,Bool_t DoMC,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
-{
- //Normal constructor
- fParamSpace = 0;
- fDoMC = kFALSE;
- fDoMC=kFALSE;
-
-#ifdef do_mc
- fTrackID = 0;
-#endif
- fZVertex = zvertex;
-#ifdef SUBSLICES
- fEtaSpace = 0;
- fEtaSpaceSize = 800/n_eta_segments;
-#endif
-#ifdef COUNTINGROWS
- fLUT2sinphi0up=0;
- fLUT2sinphi0low=0;
- fLUT2cosphi0up=0;
- fLUT2cosphi0low=0;
-#endif
-}
-
-AliL3HoughTransformerGap::~AliL3HoughTransformerGap()
-{
- DeleteHistograms();
-#ifdef do_mc
- if(fTrackID)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fTrackID[i]) continue;
- delete fTrackID[i];
- }
- delete [] fTrackID;
- }
-#endif
-#ifdef SUBSLICES
- if(fEtaSpace)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fEtaSpace[i]) continue;
- delete fEtaSpace[i];
- }
- delete [] fEtaSpace;
- }
-#endif
-#ifdef COUNTINGROWS
- if(fRowCount)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fRowCount[i]) continue;
- delete fRowCount[i];
- }
- delete [] fRowCount;
- }
- if(fGapCount)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fGapCount[i]) continue;
- delete fGapCount[i];
- }
- delete [] fGapCount;
- }
- if(fCurrentRowCount)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(fCurrentRowCount[i])
- delete fCurrentRowCount[i];
- }
- delete [] fCurrentRowCount;
- }
-#endif
-}
-
-void AliL3HoughTransformerGap::DeleteHistograms()
-{
-#ifdef COUNTINGROWS
- delete[] fLUT2sinphi0up;
- delete[] fLUT2cosphi0up;
- delete[] fLUT2sinphi0low;
- delete[] fLUT2cosphi0low;
-
- fLUT2sinphi0up=0;
- fLUT2cosphi0up=0;
- fLUT2sinphi0low=0;
- fLUT2cosphi0low=0;
-#endif
- if(!fParamSpace)
- return;
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fParamSpace[i]) continue;
- delete fParamSpace[i];
- }
- delete [] fParamSpace;
-}
-
-void AliL3HoughTransformerGap::CreateHistograms(Int_t nxbin,Float_t pt_min,
- Int_t nybin,Float_t phimin,Float_t phimax)
-{
- //Create the histograms (parameter space).
- //These are 2D histograms, span by kappa (curvature of track)
- //and phi0 (emission angle with x-axis).
- //The arguments give the range and binning;
- //nxbin = #bins in kappa
- //nybin = #bins in phi0
- //pt_min = mimium Pt of track (corresponding to maximum kappa)
- //phi_min = mimimum phi0
- //phi_max = maximum phi0
-
- Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
- //Double_t torad = AliL3Transform::Pi()/180;
- CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
-}
-
-void AliL3HoughTransformerGap::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
- Int_t nybin,Float_t ymin,Float_t ymax)
-{
-
- fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
-
- Char_t histname[256];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- sprintf(histname,"paramspace_%d",i);
- fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
- }
-
-#ifdef do_mc
- if(fDoMC)
- {
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
- fTrackID = new TrackIndex*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fTrackID[i] = new TrackIndex[ncells];
- }
-#endif
-#ifdef SUBSLICES
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*fEtaSpaceSize*sizeof(Int_t)<<" bytes to fEtaSpace"<<endl;
- fEtaSpace = new Int_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fEtaSpace[i] = new Int_t[ncells*fEtaSpaceSize];
-#endif
-#ifdef COUNTINGROWS
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- if(!fRowCount)
- {
- cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fRowCount"<<endl;
- fRowCount = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fRowCount[i] = new UChar_t[ncells];
- }
- if(!fGapCount)
- {
- cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<endl;
- fGapCount = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fGapCount[i] = new UChar_t[ncells];
- }
- if(!fCurrentRowCount)
- {
- cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<endl;
- fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fCurrentRowCount[i] = new UChar_t[ncells];
- }
-
- //create lookup table for sin and cos
- Int_t nbins = hist->GetNbinsY()+2;
- fLUT2sinphi0up=new Float_t[nbins];
- fLUT2cosphi0up=new Float_t[nbins];
- fLUT2sinphi0low=new Float_t[nbins];
- fLUT2cosphi0low=new Float_t[nbins];
- Float_t hist_bin=hist->GetBinWidthY()/2.0;
- for(Int_t i=hist->GetFirstYbin(); i<=hist->GetLastYbin(); i++){
- Float_t phi0=hist->GetBinCenterY(i);
- fLUT2sinphi0low[i]=2.*sin(phi0+hist_bin);
- fLUT2cosphi0low[i]=2.*cos(phi0+hist_bin);
- fLUT2sinphi0up[i]=2.*sin(phi0-hist_bin);
- fLUT2cosphi0up[i]=2.*cos(phi0-hist_bin);
- }
-#endif
-
-}
-
-void AliL3HoughTransformerGap::Reset()
-{
- //Reset all the histograms
-
- if(!fParamSpace)
- {
- LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
- <<"No histograms to reset"<<ENDLOG;
- return;
- }
-
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fParamSpace[i]->Reset();
-
-#ifdef do_mc
- if(fDoMC)
- {
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
- }
-#endif
-#ifdef SUBSLICES
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- memset(fEtaSpace[i],0,ncells*fEtaSpaceSize*sizeof(Int_t));
-#endif
-#ifdef COUNTINGROWS
- AliL3Histogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- memset(fRowCount[i],0,ncells*sizeof(UChar_t));
- memset(fGapCount[i],1,ncells*sizeof(UChar_t));
- memset(fCurrentRowCount[i],160,ncells*sizeof(UChar_t));
- }
-#endif
-}
-
-Int_t AliL3HoughTransformerGap::GetEtaIndex(Double_t eta)
-{
- //Return the histogram index of the corresponding eta.
-
- Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
- Double_t index = (eta-GetEtaMin())/etaslice;
- return (Int_t)index;
-}
-
-inline AliL3Histogram *AliL3HoughTransformerGap::GetHistogram(Int_t eta_index)
-{
- if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
- return 0;
- if(!fParamSpace[eta_index])
- return 0;
- return fParamSpace[eta_index];
-}
-
-Double_t AliL3HoughTransformerGap::GetEta(Int_t eta_index,Int_t slice)
-{
- Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
- Double_t eta=0;
- eta=(Double_t)((eta_index+0.5)*eta_slice);
- return eta;
-}
-
-#ifdef SUBSLICES
-Int_t AliL3HoughTransformerGap::GetEtaSubIndex(Double_t eta,Int_t eta_index)
-{
- //Return the eta space sub index of the corresponding eta.
-
- Double_t sub_eta = GetEta(eta_index,GetSlice());
- Int_t eta_subindex = (Int_t)((eta-sub_eta+0.5/(Double_t)GetNEtaSegments())*800.);
- return eta_subindex;
-}
-
-Double_t AliL3HoughTransformerGap::GetMaxSubEta(Float_t kappa,Float_t phi0,Int_t eta_index,Int_t slice)
-{
- AliL3Histogram *hist = fParamSpace[eta_index];
- Int_t bin = hist->FindBin(kappa,phi0);
- Int_t maxbin=fEtaSpaceSize/2;
- Int_t maxvalue=0;
- for(Int_t i=0;i<fEtaSpaceSize;i++) {
- if(fEtaSpace[eta_index][bin*fEtaSpaceSize+i]>maxvalue) {
- maxbin=i;
- maxvalue=fEtaSpace[eta_index][bin*fEtaSpaceSize+i];
- }
- }
- Double_t eta = GetEta(eta_index,slice)+((Double_t)maxbin-0.5*(Double_t)fEtaSpaceSize+0.5)/800.;
- cout<<"DEBUG maxetasubindex "<<eta<<" "<<eta_index<<" "<<maxbin<<endl;
- return eta;
-}
-#endif
-
-#ifdef notused
-void AliL3HoughTransformerGap::TransformCircle()
-{
-
- AliL3DigitRowData *tempPt = GetDataPointer();
- if(!tempPt)
- {
- LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
- <<"No input data "<<ENDLOG;
- return;
- }
-
- Int_t ipatch = GetPatch();
- //Loop over the padrows:
- for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
- {
- //Get the data on this padrow:
- AliL3DigitData *digPt = tempPt->fDigitData;
- if(i != (Int_t)tempPt->fRow)
- {
- cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
- continue;
- }
-
- //Loop over the data on this padrow:
- 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((Int_t)charge <= GetLowerThreshold())
- continue;
-
- if((Int_t)charge > GetUpperThreshold())
- charge = GetUpperThreshold();
-
- Int_t sector,row;
- Float_t xyz[3];
-
- //Transform data to local cartesian coordinates:
- AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
- AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
-
- //Calculate the eta:
- xyz[2] -= fZVertex;
- Double_t eta = AliL3Transform::GetEta(xyz);
-
- //Get the corresponding index, which determines which histogram to fill:
- Int_t eta_index = GetEtaIndex(eta);
-
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
- continue;
-
- //Get the correct histogrampointer:
- AliL3Histogram *hist = fParamSpace[eta_index];
- if(!hist)
- {
- cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl;
- continue;
- }
-
- //Do the transformation:
-#ifndef COUNTINGROWS
- Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- Float_t phi = AliL3Transform::GetPhi(xyz);
-#else
- xyz[1] -= AliL3Transform::GetPadPitchWidth(ipatch)/2.0;
- Float_t R1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- Float_t phi1 = AliL3Transform::GetPhi(xyz);
- xyz[1] += AliL3Transform::GetPadPitchWidth(ipatch);
- Float_t R2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- Float_t phi2 = AliL3Transform::GetPhi(xyz);
-#endif
-#ifdef SUBSLICES
- Int_t eta_subindex = GetEtaSubIndex(eta,eta_index);
-#endif
- //Fill the histogram along the phirange
- for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
- {
- Float_t phi0 = hist->GetBinCenterY(b);
-#ifndef COUNTINGROWS
- Float_t kappa = 2*sin(phi - phi0)/R;
-#else
- Float_t kappa1 = 2*sin(phi1 - phi0 - hist->GetBinWidthY()/2.)/R1;
- Float_t kappa2 = 2*sin(phi2 - phi0 + hist->GetBinWidthY()/2.)/R2;
-#endif
- //hist->Fill(kappa,phi0,(int)rint(log((Float_t)charge)));
- // hist->Fill(kappa,phi0,charge);
-#ifndef COUNTINGROWS
- hist->Fill(kappa,phi0,30);
- //hist->Fill(kappa,phi0,1);
-#else
- Int_t binx1,binx2;
- binx1 = hist->FindXbin(kappa1);
- binx2 = hist->FindXbin(kappa2);
-
- // cout<<"DEBUG counting rows "<<ipatch<<" "<<i<<" "<<" bins "<<binx1<<" "<<binx2<<" "<<" kappas "<<kappa1<<" "<<kappa2<<endl;
-
- Int_t bin1 = hist->GetBin(binx1,b);
- Int_t bin2 = hist->GetBin(binx2,b);
- if((bin1==0) || (bin2==0)) continue;
- // cout<<"DEBUG counting rows "<<ipatch<<" "<<i<<" "<<" bins "<<binx1<<" "<<binx2<<" "<<bin1<<" "<<bin2<<" kappas "<<kappa1<<" "<<kappa2<<endl;
- for(Int_t bin=bin1;bin<=bin2;bin++) {
- fRowCount[eta_index][bin]=fRowCount[eta_index][bin]|(1<<(i-AliL3Transform::GetFirstRow(ipatch)));
- }
-#endif
-#ifndef COUNTINGROWS
-#ifdef SUBSLICES
- //Fill Eta space histogram
- Int_t bin = hist->FindBin(kappa,phi0);
- // cout<<" DEBUG etasubindex "<<eta<<" "<<eta_index<<" "<<eta_subindex<<endl;
- fEtaSpace[eta_index][bin*fEtaSpaceSize+eta_subindex]++;
-#endif
-#ifdef do_mc
- if(fDoMC)
- {
- Int_t bin = hist->FindBin(kappa,phi0);
- for(Int_t t=0; t<3; t++)
- {
- Int_t label = digPt[j].fTrackID[t];
- if(label < 0) break;
- UInt_t c;
- for(c=0; c<MaxTrack; c++)
- if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
- break;
- if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
- fTrackID[eta_index][bin].fLabel[c] = label;
- fTrackID[eta_index][bin].fNHits[c]++;
- }
- }
-#endif
-#endif
- }
- }
-
- //Move the data pointer to the next padrow:
- AliL3MemHandler::UpdateRowPointer(tempPt);
- }
-}
-#endif
-
-#ifdef COUNTINGROWS
-
-struct EtaRow {
- UChar_t start_pad;
- UChar_t end_pad;
- Bool_t found;
- Float_t start_y;
-};
-
-void AliL3HoughTransformerGap::TransformCircle()
-{
- //Assumes that all the etaslice histos are the same!
- AliL3Histogram *h = fParamSpace[0];
- //Float_t hist_bin = h->GetBinWidthY()/2.;
- Int_t first_bin = h->GetFirstYbin();
- Int_t last_bin = h->GetLastYbin();
- Float_t x_min = h->GetXmin();
- Float_t x_max = h->GetXmax();
- Float_t x_bin = (x_max-x_min)/h->GetNbinsX();
- Int_t first_binx = h->GetFirstXbin()-1;
- Int_t last_binx = h->GetLastXbin()+1;
- Int_t nbinx = h->GetNbinsX()+2;
-
- UChar_t last_pad;
- Int_t last_eta_index=-1;;
- EtaRow *eta_clust = new EtaRow[GetNEtaSegments()];
-
- AliL3DigitRowData *tempPt = GetDataPointer();
- if(!tempPt)
- {
- LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
- <<"No input data "<<ENDLOG;
- return;
- }
-
- Int_t ipatch = GetPatch();
- Double_t pad_pitch = AliL3Transform::GetPadPitchWidth(ipatch);
- Int_t islice = GetSlice();
- //Loop over the padrows:
- for(Int_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
- {
- Int_t npads = AliL3Transform::GetNPads(i)-1;
-
- last_pad = 255;
- //Flush eta clusters array
- memset(eta_clust,0,GetNEtaSegments()*sizeof(EtaRow));
-
- Float_t x = AliL3Transform::Row2X(i);
- Float_t x2 = x*x;
- Float_t y=0.,r2=0.;
-
- UChar_t lrow = i-AliL3Transform::GetFirstRow(ipatch);
-
- //Get the data on this padrow:
- AliL3DigitData *digPt = tempPt->fDigitData;
- if(i != (Int_t)tempPt->fRow)
- {
- cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
- continue;
- }
- // cout<<" Starting row "<<i<<endl;
- //Loop over the data on this padrow:
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- UShort_t charge = digPt[j].fCharge;
- if((Int_t)charge <= GetLowerThreshold())
- continue;
- UChar_t pad = digPt[j].fPad;
- UShort_t time = digPt[j].fTime;
- if(pad != last_pad)
- {
- y = (pad-0.5*npads)*pad_pitch;
- Float_t y2 = y*y;
- r2 = x2 + y2;
- last_eta_index = -1;
- }
-
- //Transform data to local cartesian coordinates:
- Float_t z = AliL3Transform::GetZFast(islice,(Int_t)time,fZVertex);
- Float_t z2 = z*z;
- //Calculate the eta:
- Double_t r = sqrt(r2+z2);
- Double_t eta = 0.5 * log((r+z)/(r-z));
- //Get the corresponding index, which determines which histogram to fill:
- Int_t eta_index = GetEtaIndex(eta);
-
- if(eta_index == last_eta_index) continue;
-
- // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<eta_index<<endl;
-
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
- continue;
-
- if(!eta_clust[eta_index].found)
- {
- eta_clust[eta_index].start_pad = pad;
- eta_clust[eta_index].end_pad = pad;
- eta_clust[eta_index].start_y = y - pad_pitch/2.0;
- eta_clust[eta_index].found = 1;
- continue;
- }
- else
- {
- if(pad <= (eta_clust[eta_index].end_pad + 1))
- eta_clust[eta_index].end_pad = pad;
- else
- {
-
- // cout<<" Cluster found at etaslice "<<eta_index<<" from pad "<<(Int_t)eta_clust[eta_index].start_pad<<" to pad "<<(Int_t)eta_clust[eta_index].end_pad<<endl;
- //Get the correct histogrampointer:
- AliL3Histogram *hist = fParamSpace[eta_index];
- if(!hist)
- {
- cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl;
- continue;
- }
- UChar_t *nrows = fRowCount[eta_index];
- UChar_t *ngaps = fGapCount[eta_index];
- UChar_t *currentrow = fCurrentRowCount[eta_index];
-
- //Do the transformation:
- Float_t start_y = eta_clust[eta_index].start_y;
- Float_t R1 = x2 + start_y*start_y;
- Float_t end_y = (eta_clust[eta_index].end_pad-0.5*npads)*pad_pitch + pad_pitch/2.0;
- Float_t R2 = x2 + end_y*end_y;
-
- //Fill the histogram along the phirange
- for(Int_t b=first_bin; b<=last_bin; b++)
- {
-#ifdef WOLUT
- Float_t phi0 = hist->GetBinCenterY(b);
- // Float_t kappa1 = 2*sin(phi1 - phi0 - hist_bin)/R1;
- // Float_t kappa2 = 2*sin(phi2 - phi0 + hist_bin)/R2;
- Float_t kappa1 = 2*(start_y*cos(phi0 + hist_bin)-x*sin(phi0 + hist_bin))/R1;
- if(kappa1>x_max) continue;
- Float_t kappa2 = 2*(end_y*cos(phi0 - hist_bin)-x*sin(phi0 - hist_bin))/R2;
- if(kappa2<x_min) break;
-#else
- Float_t kappa1 = (start_y*fLUT2cosphi0low[b]-x*fLUT2sinphi0low[b])/R1;
- if(kappa1>x_max) continue;
- Float_t kappa2 = (end_y*fLUT2cosphi0up[b]-x*fLUT2sinphi0up[b])/R2;
- if(kappa2<x_min) break;
-#endif
- Int_t binx1,binx2;
- if(kappa1<x_min)
- binx1 = first_binx;
- else
- binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
- if(kappa2>x_max)
- binx2 = last_binx;
- else
- binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
-
- Int_t temp_bin = b*nbinx + binx1;
- UChar_t *nrows2 = nrows + temp_bin;
- UChar_t *ngaps2 = ngaps + temp_bin;
- UChar_t *currentrow2 = currentrow + temp_bin;
- for(Int_t bin=binx1;bin<=binx2;bin++)
- {
- //Assumes threshold about 32
- if(*ngaps2 < 5) {
- if(lrow != *currentrow2)
- {
- (*nrows2)++;
- if(lrow > (*currentrow2 + 1))
- (*ngaps2)++;
- *currentrow2=lrow;
- }
- }
- nrows2++;
- ngaps2++;
- currentrow2++;
- }
- }
- //End of the transformation
-
- eta_clust[eta_index].start_pad = pad;
- eta_clust[eta_index].end_pad = pad;
- eta_clust[eta_index].start_y = y - pad_pitch/2.0;
- }
- }
- last_pad = pad;
- last_eta_index = eta_index;
- }
- //Write remaining clusters
- for(Int_t eta_index = 0;eta_index < GetNEtaSegments();eta_index++)
- {
- //Check for empty row
- if((eta_clust[eta_index].start_pad == 0) && (eta_clust[eta_index].end_pad == 0)) continue;
-
- //Get the correct histogrampointer:
- AliL3Histogram *hist = fParamSpace[eta_index];
- if(!hist)
- {
- cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl;
- continue;
- }
- UChar_t *nrows = fRowCount[eta_index];
- UChar_t *ngaps = fGapCount[eta_index];
- UChar_t *currentrow = fCurrentRowCount[eta_index];
-
- //Do the transformation:
- Float_t start_y = eta_clust[eta_index].start_y;
- Float_t R1 = x2 + start_y*start_y;
- Float_t end_y = (eta_clust[eta_index].end_pad-0.5*npads)*pad_pitch + pad_pitch/2.0;
- Float_t R2 = x2 + end_y*end_y;
-
- //Fill the histogram along the phirange
- for(Int_t b=first_bin; b<=last_bin; b++)
- {
-#ifdef WOLUT
- Float_t phi0 = hist->GetBinCenterY(b);
- // Float_t kappa1 = 2*sin(phi1 - phi0 - hist_bin)/R1;
- // Float_t kappa2 = 2*sin(phi2 - phi0 + hist_bin)/R2;
- Float_t kappa1 = 2*(start_y*cos(phi0 + hist_bin)-x*sin(phi0 + hist_bin))/R1;
- if(kappa1>x_max) continue;
- Float_t kappa2 = 2*(end_y*cos(phi0 - hist_bin)-x*sin(phi0 - hist_bin))/R2;
- if(kappa2<x_min) break;
-#else
- Float_t kappa1 = (start_y*fLUT2cosphi0low[b]-x*fLUT2sinphi0low[b])/R1;
- if(kappa1>x_max) continue;
- Float_t kappa2 = (end_y*fLUT2cosphi0up[b]-x*fLUT2sinphi0up[b])/R2;
- if(kappa2<x_min) break;
-#endif
- Int_t binx1,binx2;
- if(kappa1<x_min)
- binx1 = first_binx;
- else
- binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
- if(kappa2>x_max)
- binx2 = last_binx;
- else
- binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
-
- Int_t temp_bin = b*nbinx + binx1;
- UChar_t *nrows2 = nrows + temp_bin;
- UChar_t *ngaps2 = ngaps + temp_bin;
- UChar_t *currentrow2 = currentrow + temp_bin;
- for(Int_t bin=binx1;bin<=binx2;bin++)
- {
- //Assumes threshold about 32
- if(*ngaps2 < 5) {
- if(lrow != *currentrow2)
- {
- (*nrows2)++;
- if(lrow > (*currentrow2 + 1))
- (*ngaps2)++;
- *currentrow2=lrow;
- }
- }
- nrows2++;
- ngaps2++;
- currentrow2++;
- }
- }
- //End of the transformation
- }
-
- //Move the data pointer to the next padrow:
- AliL3MemHandler::UpdateRowPointer(tempPt);
- }
- delete [] eta_clust;
-}
-#endif
-
-#ifdef COUNTINGROWS
-/*
-void AliL3HoughTransformer::CalculateNRows()
-{
- Int_t ipatch = GetPatch();
-
- for(Int_t i=0; i<GetNEtaSegments(); i++) {
- AliL3Histogram *hist = fParamSpace[i];
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t sum,ngaps;
- Bool_t isPrevRow;
-
- for(Int_t ybin=ymin; ybin<=ymax; ybin++)
- {
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- sum = 0;
- ngaps = 1;
- isPrevRow = kFALSE;
- Int_t bin = hist->GetBin(xbin,ybin);
- for(Int_t irow=0;irow<=(AliL3Transform::GetLastRow(ipatch)-AliL3Transform::GetFirstRow(ipatch));irow++) {
- if(((fRowCount[i][bin]>>irow)&1)==1) {
- sum++;
- isPrevRow = kTRUE;
- }
- else {
- if(isPrevRow)
- ngaps++;
- isPrevRow = kFALSE;
- }
- }
- hist->SetBinContent(bin,sum/ngaps);
- }
- }
- hist->SetNEntries(1);
- }
-}
-*/
-
-UChar_t *AliL3HoughTransformerGap::GetRowCount(Int_t eta_index)
-{
- return fRowCount[eta_index];
-}
-
-UChar_t *AliL3HoughTransformerGap::GetGapCount(Int_t eta_index)
-{
- return fGapCount[eta_index];
-}
-
-#endif
#include <TH2.h>
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3Histogram.h"
#include "AliL3HistogramAdaptive.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3Vertex.h"
#include "AliL3Fitter.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
--- /dev/null
+// @(#) $Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright © ALICE HLT Group
+
+#include "AliL3StandardIncludes.h"
+
+#include "AliL3Logging.h"
+#include "AliL3MemHandler.h"
+#include "AliL3Transform.h"
+#include "AliL3DigitData.h"
+#include "AliL3HistogramAdaptive.h"
+#include "AliL3HoughTransformerRow.h"
+
+#if __GNUC__ == 3
+using namespace std;
+#endif
+
+//_____________________________________________________________
+// AliL3HoughTransformerRow
+//
+// TPC rows Hough transformation class
+//
+
+ClassImp(AliL3HoughTransformerRow)
+
+UChar_t **AliL3HoughTransformerRow::fRowCount = 0;
+UChar_t **AliL3HoughTransformerRow::fGapCount = 0;
+UChar_t **AliL3HoughTransformerRow::fCurrentRowCount = 0;
+#ifdef do_mc
+TrackIndex **AliL3HoughTransformerRow::fTrackID = 0;
+#endif
+
+
+AliL3HoughTransformerRow::AliL3HoughTransformerRow()
+{
+ //Default constructor
+ fParamSpace = 0;
+ fDoMC = kFALSE;;
+
+ fLUT2sinphi0up=0;
+ fLUT2sinphi0low=0;
+ fLUT2cosphi0up=0;
+ fLUT2cosphi0low=0;
+
+ fLUTforwardZ=0;
+ fLUTforwardZ2=0;
+ fLUTbackwardZ=0;
+ fLUTbackwardZ2=0;
+}
+
+AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoMC,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments,zvertex)
+{
+ //Normal constructor
+ fParamSpace = 0;
+ fDoMC = kFALSE;
+ fDoMC=kFALSE;
+#ifdef do_mc
+ fDoMC = kTRUE;
+#endif
+
+ fLUT2sinphi0up=0;
+ fLUT2sinphi0low=0;
+ fLUT2cosphi0up=0;
+ fLUT2cosphi0low=0;
+
+ fLUTforwardZ=0;
+ fLUTforwardZ2=0;
+ fLUTbackwardZ=0;
+ fLUTbackwardZ2=0;
+}
+
+AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
+{
+ DeleteHistograms();
+#ifdef do_mc
+ if(fTrackID)
+ {
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fTrackID[i]) continue;
+ delete fTrackID[i];
+ }
+ delete [] fTrackID;
+ fTrackID = 0;
+ }
+#endif
+
+ if(fRowCount)
+ {
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fRowCount[i]) continue;
+ delete [] fRowCount[i];
+ }
+ delete [] fRowCount;
+ fRowCount = 0;
+ }
+ if(fGapCount)
+ {
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fGapCount[i]) continue;
+ delete [] fGapCount[i];
+ }
+ delete [] fGapCount;
+ fGapCount = 0;
+ }
+ if(fCurrentRowCount)
+ {
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(fCurrentRowCount[i])
+ delete [] fCurrentRowCount[i];
+ }
+ delete [] fCurrentRowCount;
+ fCurrentRowCount = 0;
+ }
+}
+
+void AliL3HoughTransformerRow::DeleteHistograms()
+{
+ delete[] fLUT2sinphi0up;
+ delete[] fLUT2cosphi0up;
+ delete[] fLUT2sinphi0low;
+ delete[] fLUT2cosphi0low;
+
+ fLUT2sinphi0up=0;
+ fLUT2cosphi0up=0;
+ fLUT2sinphi0low=0;
+ fLUT2cosphi0low=0;
+
+ delete[] fLUTforwardZ;
+ delete[] fLUTforwardZ2;
+ delete[] fLUTbackwardZ;
+ delete[] fLUTbackwardZ2;
+
+ fLUTforwardZ=0;
+ fLUTforwardZ2=0;
+ fLUTbackwardZ=0;
+ fLUTbackwardZ2=0;
+
+ if(!fParamSpace)
+ return;
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fParamSpace[i]) continue;
+ delete fParamSpace[i];
+ }
+ delete [] fParamSpace;
+}
+
+void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t pt_min,
+ Int_t nybin,Float_t phimin,Float_t phimax)
+{
+ //Create the histograms (parameter space).
+ //These are 2D histograms, span by kappa (curvature of track)
+ //and phi0 (emission angle with x-axis).
+ //The arguments give the range and binning;
+ //nxbin = #bins in kappa
+ //nybin = #bins in phi0
+ //pt_min = mimium Pt of track (corresponding to maximum kappa)
+ //phi_min = mimimum phi0
+ //phi_max = maximum phi0
+
+ Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
+ //Double_t torad = AliL3Transform::Pi()/180;
+ CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
+}
+
+void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
+ Int_t nybin,Float_t ymin,Float_t ymax)
+{
+
+ fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
+
+ Char_t histname[256];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ sprintf(histname,"paramspace_%d",i);
+ fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
+ }
+
+#ifdef do_mc
+ if(fDoMC)
+ {
+ AliL3Histogram *hist = fParamSpace[0];
+ Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
+ if(!fTrackID)
+ {
+ cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
+ fTrackID = new TrackIndex*[GetNEtaSegments()];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fTrackID[i] = new TrackIndex[ncells];
+ }
+ }
+#endif
+ AliL3Histogram *hist = fParamSpace[0];
+ Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
+ if(!fRowCount)
+ {
+ cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fRowCount"<<endl;
+ fRowCount = new UChar_t*[GetNEtaSegments()];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fRowCount[i] = new UChar_t[ncells];
+ }
+ if(!fGapCount)
+ {
+ cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<endl;
+ fGapCount = new UChar_t*[GetNEtaSegments()];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fGapCount[i] = new UChar_t[ncells];
+ }
+ if(!fCurrentRowCount)
+ {
+ cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<endl;
+ fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fCurrentRowCount[i] = new UChar_t[ncells];
+ }
+
+ //create lookup table for sin and cos
+ Int_t nbins = hist->GetNbinsY()+2;
+ fLUT2sinphi0up=new Float_t[nbins];
+ fLUT2cosphi0up=new Float_t[nbins];
+ fLUT2sinphi0low=new Float_t[nbins];
+ fLUT2cosphi0low=new Float_t[nbins];
+ Float_t hist_bin=hist->GetBinWidthY()/2.0;
+ for(Int_t i=hist->GetFirstYbin(); i<=hist->GetLastYbin(); i++){
+ Float_t phi0=hist->GetBinCenterY(i);
+ fLUT2sinphi0low[i]=2.*sin(phi0+hist_bin);
+ fLUT2cosphi0low[i]=2.*cos(phi0+hist_bin);
+ fLUT2sinphi0up[i]=2.*sin(phi0-hist_bin);
+ fLUT2cosphi0up[i]=2.*cos(phi0-hist_bin);
+ }
+
+ //create lookup table for z of the digits
+ Int_t ntimebins = AliL3Transform::GetNTimeBins();
+ fLUTforwardZ = new Float_t[ntimebins];
+ fLUTforwardZ2 = new Float_t[ntimebins];
+ fLUTbackwardZ = new Float_t[ntimebins];
+ fLUTbackwardZ2 = new Float_t[ntimebins];
+ for(Int_t i=0; i<ntimebins; i++){
+ Float_t z;
+ z=AliL3Transform::GetZFast(0,i,GetZVertex());
+ fLUTforwardZ[i]=z;
+ fLUTforwardZ2[i]=z*z;
+ z=AliL3Transform::GetZFast(18,i,GetZVertex());
+ fLUTbackwardZ[i]=z;
+ fLUTbackwardZ2[i]=z*z;
+ }
+
+}
+
+void AliL3HoughTransformerRow::Reset()
+{
+ //Reset all the histograms
+
+ if(!fParamSpace)
+ {
+ LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
+ <<"No histograms to reset"<<ENDLOG;
+ return;
+ }
+
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fParamSpace[i]->Reset();
+
+#ifdef do_mc
+ if(fDoMC)
+ {
+ AliL3Histogram *hist = fParamSpace[0];
+ Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
+ }
+#endif
+ AliL3Histogram *hist = fParamSpace[0];
+ Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ memset(fRowCount[i],0,ncells*sizeof(UChar_t));
+ memset(fGapCount[i],1,ncells*sizeof(UChar_t));
+ memset(fCurrentRowCount[i],160,ncells*sizeof(UChar_t));
+ }
+}
+
+Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta)
+{
+ //Return the histogram index of the corresponding eta.
+
+ Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
+ Double_t index = (eta-GetEtaMin())/etaslice;
+ return (Int_t)index;
+}
+
+inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t eta_index)
+{
+ if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
+ return 0;
+ if(!fParamSpace[eta_index])
+ return 0;
+ return fParamSpace[eta_index];
+}
+
+Double_t AliL3HoughTransformerRow::GetEta(Int_t eta_index,Int_t slice)
+{
+ Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
+ Double_t eta=0;
+ eta=(Double_t)((eta_index+0.5)*eta_slice);
+ return eta;
+}
+
+struct EtaRow {
+ UChar_t start_pad;
+ UChar_t end_pad;
+ Bool_t found;
+ Float_t start_y;
+#ifdef do_mc
+ Int_t mc_labels[MaxTrack];
+#endif
+};
+
+void AliL3HoughTransformerRow::TransformCircle()
+{
+ Int_t n_eta_segments = GetNEtaSegments();
+ Double_t eta_min = GetEtaMin();
+ Double_t eta_slice = (GetEtaMax() - eta_min)/n_eta_segments;
+
+ Int_t lower_threshold = GetLowerThreshold();
+
+ //Assumes that all the etaslice histos are the same!
+ AliL3Histogram *h = fParamSpace[0];
+ Int_t first_bin = h->GetFirstYbin();
+ Int_t last_bin = h->GetLastYbin();
+ Float_t x_min = h->GetXmin();
+ Float_t x_max = h->GetXmax();
+ Float_t x_bin = (x_max-x_min)/h->GetNbinsX();
+ Int_t first_binx = h->GetFirstXbin()-1;
+ Int_t last_binx = h->GetLastXbin()+1;
+ Int_t nbinx = h->GetNbinsX()+2;
+
+ UChar_t last_pad;
+ Int_t last_eta_index;
+ EtaRow *eta_clust = new EtaRow[n_eta_segments];
+
+ AliL3DigitRowData *tempPt = GetDataPointer();
+ if(!tempPt)
+ {
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
+ <<"No input data "<<ENDLOG;
+ return;
+ }
+
+ Int_t ipatch = GetPatch();
+ Double_t pad_pitch = AliL3Transform::GetPadPitchWidth(ipatch);
+ Int_t islice = GetSlice();
+ Float_t *fLUTz;
+ Float_t *fLUTz2;
+ if(islice < 18) {
+ fLUTz = fLUTforwardZ;
+ fLUTz2 = fLUTforwardZ2;
+ }
+ else {
+ fLUTz = fLUTbackwardZ;
+ fLUTz2 = fLUTbackwardZ2;
+ }
+
+ Int_t npads_in_patch = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(ipatch))+2;
+ Bool_t **IsPrevRow = new Bool_t*[n_eta_segments];
+ Int_t nrows_in_patch = AliL3Transform::GetLastRow(ipatch)-AliL3Transform::GetFirstRow(ipatch)+2;
+ for(Int_t i=0; i<n_eta_segments; i++) {
+ IsPrevRow[i] = new Bool_t[npads_in_patch*nrows_in_patch];
+ memset(IsPrevRow[i],1,npads_in_patch*nrows_in_patch*sizeof(Bool_t));
+ }
+
+ //Loop over the padrows:
+ for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
+ {
+ Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
+
+ last_pad = 255;
+ //Flush eta clusters array
+ memset(eta_clust,0,n_eta_segments*sizeof(EtaRow));
+
+ Float_t x = AliL3Transform::Row2X((Int_t)i);
+ Float_t x2 = x*x;
+ Float_t y,r2;
+
+ Int_t lrow = (i-AliL3Transform::GetFirstRow(ipatch))*npads_in_patch;
+
+ //Get the data on this padrow:
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ if((Int_t)i != (Int_t)tempPt->fRow)
+ {
+ cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<endl;
+ continue;
+ }
+ // cout<<" Starting row "<<i<<endl;
+ //Loop over the data on this padrow:
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
+ {
+ UShort_t charge = digPt[j].fCharge;
+ if((Int_t)charge <= lower_threshold)
+ continue;
+ UChar_t pad = digPt[j].fPad;
+ UShort_t time = digPt[j].fTime;
+ if(pad != last_pad)
+ {
+ y = (pad-0.5*npads)*pad_pitch;
+ Float_t y2 = y*y;
+ r2 = x2 + y2;
+ last_eta_index = -1;
+ }
+
+ //Transform data to local cartesian coordinates:
+ Float_t z = fLUTz[(Int_t)time];
+ Float_t z2 = fLUTz2[(Int_t)time];
+ //Calculate the eta:
+ Double_t r = sqrt(r2+z2);
+ Double_t eta = 0.5 * log((r+z)/(r-z));
+ //Get the corresponding index, which determines which histogram to fill:
+ Int_t eta_index = (Int_t)((eta-eta_min)/eta_slice);
+
+#ifndef do_mc
+ if(eta_index == last_eta_index) continue;
+#endif
+ // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<eta_index<<endl;
+
+ if(eta_index < 0 || eta_index >= n_eta_segments)
+ continue;
+
+ if(!eta_clust[eta_index].found)
+ {
+ eta_clust[eta_index].start_pad = pad;
+ eta_clust[eta_index].end_pad = pad;
+ eta_clust[eta_index].start_y = y - pad_pitch/2.0;
+ eta_clust[eta_index].found = 1;
+#ifdef do_mc
+ if(fDoMC)
+ {
+ for(Int_t t=0; t<3; t++)
+ {
+ Int_t label = digPt[j].fTrackID[t];
+ if(label < 0) break;
+ UInt_t c;
+ for(c=0; c<(MaxTrack-1); c++)
+ if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ break;
+ // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
+ eta_clust[eta_index].mc_labels[c] = label;
+ }
+ }
+#endif
+ continue;
+ }
+ else
+ {
+ if(pad <= (eta_clust[eta_index].end_pad + 1))
+ {
+ eta_clust[eta_index].end_pad = pad;
+#ifdef do_mc
+ if(fDoMC)
+ {
+ for(Int_t t=0; t<3; t++)
+ {
+ Int_t label = digPt[j].fTrackID[t];
+ if(label < 0) break;
+ UInt_t c;
+ for(c=0; c<(MaxTrack-1); c++)
+ if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ break;
+ // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
+ eta_clust[eta_index].mc_labels[c] = label;
+ }
+ }
+#endif
+ }
+ else
+ {
+
+ Bool_t fill_cluster = kFALSE;
+ Bool_t *IsPrevRow2 = &IsPrevRow[eta_index][lrow+eta_clust[eta_index].start_pad];
+ for(Int_t ipad = 0; ipad <= (eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 2); ipad++)
+ {
+ if(*IsPrevRow2)
+ {
+ fill_cluster = kTRUE;
+ break;
+ }
+ IsPrevRow2++;
+ }
+
+ // if(eta_index == 51)
+ // cout<<" Cluster to fill:"<<" row "<<(Int_t)i<<" pads from "<<(Int_t)eta_clust[eta_index].start_pad<<"("<<eta_clust[eta_index].start_y<<") to "<<(Int_t)eta_clust[eta_index].end_pad<<"("<<(eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch<<") fill "<<fill_cluster<<endl;
+
+ if(fill_cluster) {
+ // cout<<" Cluster found at etaslice "<<eta_index<<" from pad "<<(Int_t)eta_clust[eta_index].start_pad<<" to pad "<<(Int_t)eta_clust[eta_index].end_pad<<endl;
+ // Bool_t fill_next_row = kFALSE;
+
+ UChar_t *nrows = fRowCount[eta_index];
+ UChar_t *ngaps = fGapCount[eta_index];
+ UChar_t *currentrow = fCurrentRowCount[eta_index];
+
+ //Do the transformation:
+ Float_t start_y = eta_clust[eta_index].start_y;
+ if(eta_clust[eta_index].start_pad == 0)
+ start_y -= 3.0*pad_pitch;
+ Float_t R1 = x2 + start_y*start_y;
+ Float_t x_over_R1 = x/R1;
+ Float_t start_y_over_R1 = start_y/R1;
+ Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
+ if(eta_clust[eta_index].end_pad == npads)
+ end_y += 3.0*pad_pitch;
+ Float_t R2 = x2 + end_y*end_y;
+ Float_t x_over_R2 = x/R2;
+ Float_t end_y_over_R2 = end_y/R2;
+
+ //Fill the histogram along the phirange
+ for(Int_t b=first_bin; b<=last_bin; b++)
+ {
+ Float_t kappa1 = start_y_over_R1*fLUT2cosphi0low[b]-x_over_R1*fLUT2sinphi0low[b];
+ if(kappa1>x_max) continue;
+ Float_t kappa2 = end_y_over_R2*fLUT2cosphi0up[b]-x_over_R2*fLUT2sinphi0up[b];
+ if(kappa2<x_min) break;
+
+ Int_t binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
+ if(binx1<first_binx) binx1 = first_binx;
+ Int_t binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
+ if(binx2>last_binx) binx2 = last_binx;
+ // if(eta_index == 51)
+ // cout<<" phi bin "<<b<<" kappa1 "<<kappa1<<" kappa2 "<<kappa2<<" from bin "<<binx1<<" to "<<binx2<<endl;
+ Int_t temp_bin = b*nbinx + binx1;
+ UChar_t *nrows2 = nrows + temp_bin;
+ UChar_t *ngaps2 = ngaps + temp_bin;
+ UChar_t *currentrow2 = currentrow + temp_bin;
+ for(Int_t bin=binx1;bin<=binx2;bin++)
+ {
+ //Assumes threshold about 32
+ if(*ngaps2 < 3) {
+ if(i != *currentrow2)
+ {
+ (*nrows2)++;
+ if(i > (*currentrow2 + 1))
+ (*ngaps2)++;
+ *currentrow2=i;
+ }
+#ifdef do_mc
+ if(fDoMC)
+ {
+ for(UInt_t t=0;t<(MaxTrack-1); t++)
+ {
+ Int_t label = eta_clust[eta_index].mc_labels[t];
+ if(label == 0) break;
+ UInt_t c;
+ Int_t temp_bin2 = b*nbinx + bin;
+ for(c=0; c<(MaxTrack-1); c++)
+ if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
+ break;
+ // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
+ fTrackID[eta_index][temp_bin2].fLabel[c] = label;
+ if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
+ fTrackID[eta_index][temp_bin2].fNHits[c]++;
+ fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
+ }
+ }
+ }
+#endif
+ }
+ nrows2++;
+ ngaps2++;
+ currentrow2++;
+ }
+ }
+ //End of the transformation
+
+ /* if(!fill_next_row) {
+ Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch];
+ memset(&IsCurrentRow[eta_clust[eta_index].start_pad + 1],0,eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 1);
+ } */
+ }
+ else {
+ Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].start_pad+1];
+ memset(IsCurrentRow,0,eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 1);
+ }
+
+ Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].end_pad+2];
+ memset(IsCurrentRow,0,pad - eta_clust[eta_index].end_pad - 1);
+
+ eta_clust[eta_index].start_pad = pad;
+ eta_clust[eta_index].end_pad = pad;
+ eta_clust[eta_index].start_y = y - pad_pitch/2.0;
+
+#ifdef do_mc
+ if(fDoMC)
+ {
+ memset(eta_clust[eta_index].mc_labels,0,MaxTrack);
+ for(Int_t t=0; t<3; t++)
+ {
+ Int_t label = digPt[j].fTrackID[t];
+ if(label < 0) break;
+ UInt_t c;
+ for(c=0; c<(MaxTrack-1); c++)
+ if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ break;
+ // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
+ eta_clust[eta_index].mc_labels[c] = label;
+ }
+ }
+#endif
+ }
+ }
+ last_pad = pad;
+ last_eta_index = eta_index;
+ }
+ //Write remaining clusters
+ for(Int_t eta_index = 0;eta_index < n_eta_segments;eta_index++)
+ {
+ //Check for empty row
+ if((eta_clust[eta_index].start_pad == 0) && (eta_clust[eta_index].end_pad == 0)) continue;
+
+ UChar_t *nrows = fRowCount[eta_index];
+ UChar_t *ngaps = fGapCount[eta_index];
+ UChar_t *currentrow = fCurrentRowCount[eta_index];
+
+ //Do the transformation:
+ Float_t start_y = eta_clust[eta_index].start_y;
+ if(eta_clust[eta_index].start_pad == 0)
+ start_y -= 3.0*pad_pitch;
+ Float_t R1 = x2 + start_y*start_y;
+ Float_t x_over_R1 = x/R1;
+ Float_t start_y_over_R1 = start_y/R1;
+ Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
+ if(eta_clust[eta_index].end_pad == npads)
+ end_y += 3.0*pad_pitch;
+ Float_t R2 = x2 + end_y*end_y;
+ Float_t x_over_R2 = x/R2;
+ Float_t end_y_over_R2 = end_y/R2;
+
+ //Fill the histogram along the phirange
+ for(Int_t b=first_bin; b<=last_bin; b++)
+ {
+ Float_t kappa1 = start_y_over_R1*fLUT2cosphi0low[b]-x_over_R1*fLUT2sinphi0low[b];
+ if(kappa1>x_max) continue;
+ Float_t kappa2 = end_y_over_R2*fLUT2cosphi0up[b]-x_over_R2*fLUT2sinphi0up[b];
+ if(kappa2<x_min) break;
+
+ Int_t binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
+ if(binx1<first_binx) binx1 = first_binx;
+ Int_t binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
+ if(binx2>last_binx) binx2 = last_binx;
+
+ Int_t temp_bin = b*nbinx + binx1;
+ UChar_t *nrows2 = nrows + temp_bin;
+ UChar_t *ngaps2 = ngaps + temp_bin;
+ UChar_t *currentrow2 = currentrow + temp_bin;
+ for(Int_t bin=binx1;bin<=binx2;bin++)
+ {
+ //Assumes threshold about 32
+ if(*ngaps2 < 3) {
+ if(i != *currentrow2)
+ {
+ (*nrows2)++;
+ if(i > (*currentrow2 + 1))
+ (*ngaps2)++;
+ *currentrow2=i;
+ }
+#ifdef do_mc
+ if(fDoMC)
+ {
+ for(UInt_t t=0;t<(MaxTrack-1); t++)
+ {
+ Int_t label = eta_clust[eta_index].mc_labels[t];
+ if(label == 0) break;
+ UInt_t c;
+ Int_t temp_bin2 = b*nbinx + bin;
+ for(c=0; c<(MaxTrack-1); c++)
+ if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
+ break;
+ // if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
+ fTrackID[eta_index][temp_bin2].fLabel[c] = label;
+ if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
+ fTrackID[eta_index][temp_bin2].fNHits[c]++;
+ fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
+ }
+ }
+ }
+#endif
+ }
+ nrows2++;
+ ngaps2++;
+ currentrow2++;
+ }
+ }
+ //End of the transformation
+
+ if(eta_clust[eta_index].end_pad < npads)
+ {
+ Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].end_pad+2];
+ memset(IsCurrentRow,0,npads - eta_clust[eta_index].end_pad);
+ }
+
+ }
+
+ //Move the data pointer to the next padrow:
+ AliL3MemHandler::UpdateRowPointer(tempPt);
+ }
+
+
+ for(Int_t i=0; i<n_eta_segments; i++)
+ {
+ delete [] IsPrevRow[i];
+ }
+ delete [] IsPrevRow;
+
+ delete [] eta_clust;
+}
+
+Int_t AliL3HoughTransformerRow::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
+{
+ if(!fDoMC)
+ {
+ cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
+ return -1;
+ }
+
+#ifdef do_mc
+ if(eta_index < 0 || eta_index > GetNEtaSegments())
+ {
+ cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
+ return -1;
+ }
+ AliL3Histogram *hist = fParamSpace[eta_index];
+ Int_t bin = hist->FindBin(kappa,psi);
+ Int_t label=-1;
+ Int_t max=0;
+ for(UInt_t i=0; i<(MaxTrack-1); i++)
+ {
+ Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
+ if(nhits == 0) break;
+ if(nhits > max)
+ {
+ max = nhits;
+ label = fTrackID[eta_index][bin].fLabel[i];
+ }
+ }
+ Int_t label2=-1;
+ Int_t max2=0;
+ for(UInt_t i=0; i<(MaxTrack-1); i++)
+ {
+ Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
+ if(nhits == 0) break;
+ if(nhits > max2)
+ {
+ if(fTrackID[eta_index][bin].fLabel[i]!=label) {
+ max2 = nhits;
+ label2 = fTrackID[eta_index][bin].fLabel[i];
+ }
+ }
+ }
+ // cout<<" TrackID label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fTrackID[eta_index][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fTrackID[eta_index][bin].fNHits[MaxTrack-1]<<endl;
+ return label;
+#endif
+ cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
+ return -1;
+}
+
+UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t eta_index)
+{
+ return fRowCount[eta_index];
+}
+
+UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t eta_index)
+{
+ return fGapCount[eta_index];
+}
// @(#) $Id$
-#ifndef ALIL3_HOUGHTRANSFORMERGAP
-#define ALIL3_HOUGHTRANSFORMERGAP
+#ifndef ALIL3_HOUGHTRANSFORMERROW
+#define ALIL3_HOUGHTRANSFORMERROW
#include "AliL3RootTypes.h"
#include "AliL3HoughBaseTransformer.h"
class AliL3Histogram;
-#define COUNTINGROWS
-//#define SUBSLICES
-
-class AliL3HoughTransformerGap : public AliL3HoughBaseTransformer {
+class AliL3HoughTransformerRow : public AliL3HoughBaseTransformer {
private:
AliL3Histogram **fParamSpace; //!
#ifdef do_mc
- TrackIndex **fTrackID; //!
+ static TrackIndex **fTrackID; //!
#endif
Bool_t fDoMC;
- Float_t fZVertex;
-
void DeleteHistograms();
-#ifdef SUBSLICES
- Int_t **fEtaSpace; //!
- Int_t fEtaSpaceSize;
-#endif
-#ifdef COUNTINGROWS
static UChar_t **fRowCount; //!
static UChar_t **fGapCount; //!
static UChar_t **fCurrentRowCount; //!
Float_t *fLUT2cosphi0up; //!
Float_t *fLUT2sinphi0low; //!
Float_t *fLUT2cosphi0low; //!
-#endif
+
+ Float_t *fLUTforwardZ; //!
+ Float_t *fLUTforwardZ2; //!
+ Float_t *fLUTbackwardZ; //!
+ Float_t *fLUTbackwardZ2; //!
public:
- AliL3HoughTransformerGap();
- AliL3HoughTransformerGap(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap=kFALSE,Bool_t DoMC=kFALSE,Float_t zvertex=0.0);
- virtual ~AliL3HoughTransformerGap();
+ AliL3HoughTransformerRow();
+ AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoMC=kFALSE,Float_t zvertex=0.0);
+ virtual ~AliL3HoughTransformerRow();
//void CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,Int_t nybin,Float_t psi);
void CreateHistograms(Int_t nxbin,Float_t ptmin,Int_t nybin,Float_t phimin,Float_t phimax);
Int_t GetEtaIndex(Double_t eta);
AliL3Histogram *GetHistogram(Int_t eta_index);
Double_t GetEta(Int_t eta_index,Int_t slice);
-
-#ifdef SUBSLICES
- Int_t GetEtaSubIndex(Double_t eta,Int_t eta_index);
- Double_t GetMaxSubEta(Float_t kappa,Float_t phi0,Int_t eta_index,Int_t slice);
-#endif
-
- //Int_t GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi);
-
+ Int_t GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi);
UChar_t *GetRowCount(Int_t eta_index);
UChar_t *GetGapCount(Int_t eta_index);
-#ifdef COUNTINGROWS
- //void CalculateNRows();
-#endif
-
- ClassDef(AliL3HoughTransformerGap,1) //Hough transformation class
+ ClassDef(AliL3HoughTransformerRow,1) //TPC Rows Hough transformation class
};
#include "AliL3HoughTransformerVhdl.h"
#include "AliL3FFloat.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
AliL3HoughEval.cxx AliL3HoughMerger.cxx AliL3HoughBaseTransformer.cxx \
AliL3HoughIntMerger.cxx AliL3HoughGlobalMerger.cxx AliL3HoughTransformerVhdl.cxx \
AliL3Histogram.cxx AliL3Histogram1D.cxx AliL3HoughMaxFinder.cxx AliL3Hough.cxx \
- AliL3HoughTransformerLUT.cxx AliL3HoughTransformerGap.cxx AliL3HistogramAdaptive.cxx
+ AliL3HoughTransformerLUT.cxx AliL3HoughTransformerRow.cxx AliL3HistogramAdaptive.cxx
ifeq ($(ARCH),Darwin)
## AliL3HoughTrack put into src as symbolic link
fWriteOut = kTRUE;
fBenchmark = 0;
- fMakeSeed = kFALSE;
}
AliL3Kalman::~AliL3Kalman()
fKalmanTracks = new AliL3TrackArray();
- // Make a tree to store state vector, covariance matrix and chisquare
+ // Make a ntuple to store state vector, covariance matrix and chisquare
// Will eventually not need a TTree??
- TTree *kalmanTree = new TTree("kalmanTree","kalmantracks");
+ TNtuple *kalmanTree = new TNtuple("kalmanTree","kalmantracks","x0:x1:x2:x3:x4:c0:c1:c2:c3:c4:c5:c6:c7:c8:c9:c10:c11:c12:c13:c14:chisq");
Float_t meas[21];
- Int_t lab = 123456789;
- kalmanTree->Branch("x0",&meas[0],"x0/F");
- kalmanTree->Branch("x1",&meas[1],"x1/F");
- kalmanTree->Branch("x2",&meas[2],"x2/F");
- kalmanTree->Branch("x3",&meas[3],"x3/F");
- kalmanTree->Branch("x4",&meas[4],"x4/F");
- kalmanTree->Branch("c0",&meas[5],"c0/F");
- kalmanTree->Branch("c1",&meas[6],"c1/F");
- kalmanTree->Branch("c2",&meas[7],"c2/F");
- kalmanTree->Branch("c3",&meas[8],"c3/F");
- kalmanTree->Branch("c4",&meas[9],"c4/F");
- kalmanTree->Branch("c5",&meas[10],"c5/F");
- kalmanTree->Branch("c6",&meas[11],"c6/F");
- kalmanTree->Branch("c7",&meas[12],"c7/F");
- kalmanTree->Branch("c8",&meas[13],"c8/F");
- kalmanTree->Branch("c9",&meas[14],"c9/F");
- kalmanTree->Branch("c10",&meas[15],"c10/F");
- kalmanTree->Branch("c11",&meas[16],"c11/F");
- kalmanTree->Branch("c12",&meas[17],"c12/F");
- kalmanTree->Branch("c13",&meas[18],"c13/F");
- kalmanTree->Branch("c14",&meas[19],"c14/F");
- kalmanTree->Branch("chisq",&meas[20],"chisq/F");
- kalmanTree->Branch("lab",&lab,"lab/I");
// Go through the tracks from conformal or hough tracker
for (Int_t iTrack = 0; iTrack < fTracks->GetNTracks(); iTrack++)
if (!track) continue;
if (track->GetNumberOfPoints() < fMinPointsOnTrack) continue;
- UInt_t *hitnum = track->GetHitNumbers();
- UInt_t id;
-
- id = hitnum[0];
- Int_t slice = (id>>25) & 0x7f;
- Int_t patch = (id>>22) & 0x7;
- UInt_t pos = id&0x3fffff;
- AliL3SpacePointData *points = fClusters[slice][patch];
- lab=points[pos].fTrackID[0];
- //cout << lab << endl;
-
AliL3KalmanTrack *kalmantrack = new AliL3KalmanTrack();
Bool_t save = kTRUE;
- if (fMakeSeed == kTRUE)
- {
- if (MakeSeed(kalmantrack, track) == 0)
- {
- save = kFALSE;
- continue;
- }
- }
- else
+
+ if (InitKalmanTrack(kalmantrack, track) == 0)
{
- if (InitKalmanTrack(kalmantrack, track) == 0)
- {
- save = kFALSE;
- continue;
- }
+ save = kFALSE;
+ continue;
}
+
if (Propagate(kalmantrack, track) == 0)
{
save = kFALSE;
- continue;
}
if (save) {
outtrack->Set(track);
// Fill the ntuple with the state vector, covariance matrix and
- // chisquare and track label
- kalmanTree->Fill();
+ // chisquare
+ kalmanTree->Fill(meas);
}
delete track;
UInt_t pos = id&0x3fffff;
AliL3SpacePointData *points = fClusters[slice][patch];
- return kalmantrack->Init(track, points, pos);
-}
-
-Int_t AliL3Kalman::MakeSeed(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
-{
- // Makes a rough state vector and covariance matrix based on three
- // space points of the loaded track
- // Will eventually be removed?? Will use track parameters from conformal
- // tracker or HT as seeds for the kalman filter.
-
- UInt_t *hitnum = track->GetHitNumbers();
- UInt_t id1, id2, id3;
-
- id1 = hitnum[0];
- Int_t slice1 = (id1>>25) & 0x7f;
- Int_t patch1 = (id1>>22) & 0x7;
- UInt_t pos1 = id1&0x3fffff;
- AliL3SpacePointData *points1 = fClusters[slice1][patch1];
-
- id2 = hitnum[Int_t(track->GetNHits()/2)];
- Int_t slice2 = (id2>>25) & 0x7f;
- Int_t patch2 = (id2>>22) & 0x7;
- UInt_t pos2 = id2&0x3fffff;
- AliL3SpacePointData *points2 = fClusters[slice2][patch2];
-
- id3 = hitnum[track->GetNHits()-1];
- Int_t slice3 = (id3>>25) & 0x7f;
- Int_t patch3 = (id3>>22) & 0x7;
- UInt_t pos3 = id3&0x3fffff;
- AliL3SpacePointData *points3 = fClusters[slice3][patch3];
-
- return kalmantrack->MakeTrackSeed(points1,pos1,points2,pos2,points3,pos3);
-
+ return kalmantrack->Init(track, points, pos, slice);
}
Int_t AliL3Kalman::Propagate(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
AliL3SpacePointData *points = fClusters[slice][patch];
if (!points) continue;
- if (kalmantrack->Propagate(points,pos) == 0)
+ if (kalmantrack->Propagate(points,pos,slice) == 0)
{
badpoint++;
continue;
Char_t fWriteOutPath[256];
Bool_t fWriteOut;
Int_t fEvent;
- Bool_t fMakeSeed;
public:
- AliL3Kalman(Char_t *datapath, Int_t *slice=0, Int_t min_clusters=0);
+ AliL3Kalman(Char_t *datapath, Int_t *slice, Int_t min_clusters);
virtual ~AliL3Kalman();
void Init();
void LoadTracks(Int_t event, Bool_t sp);
void ProcessTracks();
Int_t InitKalmanTrack(AliL3KalmanTrack *kalmantrack, AliL3Track *track);
- Int_t MakeSeed(AliL3KalmanTrack *kalmantrack, AliL3Track *track);
Int_t Propagate(AliL3KalmanTrack *kalmantrack, AliL3Track *track);
Int_t Update(AliL3SpacePointData *points, UInt_t pos, AliL3KalmanTrack *kalmantrack);
void WriteFiles(Char_t *path="data"){fWriteOut = kTRUE; sprintf(fWriteOutPath,"%s",path);}
- void DoMakeSeed(){fMakeSeed = kTRUE;}
Double_t GetCpuTime();
AliL3TrackArray *GetTracks() {return fKalmanTracks;}
};
#include "AliL3SpacePointData.h"
#include "AliL3Logging.h"
#include "AliL3Transform.h"
+#include "AliL3StandardIncludes.h"
+
+// includes for offline comparison, will be removed
+#include "AliTPCtrack.h"
+// includes for offline comparison, will be removed
+
+#include "Riostream.h"
+
+
ClassImp(AliL3KalmanTrack)
// Class for kalman tracks
AliL3KalmanTrack::AliL3KalmanTrack()
{
+ fX = 0;
+
fMaxChi2 = 12;
// Constructor
}
// Destructor
}
-Int_t AliL3KalmanTrack::Init(AliL3Track *track, AliL3SpacePointData *points, UInt_t pos)
+Int_t AliL3KalmanTrack::Init(AliL3Track *track, AliL3SpacePointData *points, UInt_t pos, Int_t slice)
{
- fX = points[pos].fX;
- fP0 = points[pos].fY;
+ Float_t xyz[3];
+
+ xyz[0] = points[pos].fX;
+ xyz[1] = points[pos].fY;
+
+ // NB! I think boolean variable in the command under is true if single slice.
+ // Better do a test if it is signle slice, and if so set it true.??.
+ AliL3Transform::Global2Local(xyz,slice);
+
+ fX = xyz[0];
+
+ fP0 = xyz[1];
fP1 = points[pos].fZ;
fP3 = track->GetTgl(); //cout << fP3;
if (TMath::Abs(fP3) > 1.2) return 0; //From AliTPCtrackerMI
fP4 = 0.5*(-track->GetCharge()*1./track->GetPt())*0.0029980*AliL3Transform::GetBField(); // Calculation of the curvature
-
if (TMath::Abs(fP4) >= 0.0066) return 0; // From AliTPCtrackerMI
- Float_t centerX = track->GetFirstPointX() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Cos(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
- Float_t centerY = track->GetFirstPointY() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Sin(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
+
+ Float_t firstXY[2];
+ firstXY[0] = track->GetFirstPointX();
+ firstXY[1] = track->GetFirstPointY();
+
+ AliL3Transform::Global2Local(firstXY,slice);
+
+ //cout << firstXY[0] << endl;
+
+ //Float_t centerX = track->GetFirstPointX() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Cos(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
+ Float_t centerX = firstXY[0] - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Cos(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
+ //Float_t centerY = track->GetFirstPointY() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Sin(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
+
fP2 = fP4*centerX; // Curvature times center of curvature
-
+
if (TMath::Abs(fP4*fX - fP2) >= 0.9) // What's this??
{
return 0;
}
//cout << track->GetCharge() << endl;
//cout << "AliL3KalmanTrack::Init, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
- //cout << fP2 << " " << fP4 << endl;
- Float_t num = 12;
+ //cout << fP4 << endl;
+ //Float_t num = 12;
fC00 = points[pos].fSigmaY2;
fC10 = 0; fC11 = points[pos].fSigmaZ2;
return 1;
}
-Int_t AliL3KalmanTrack::MakeTrackSeed(AliL3SpacePointData *points1, UInt_t pos1, AliL3SpacePointData *points2, UInt_t pos2, AliL3SpacePointData *points3, UInt_t pos3)
-{
- // Make track seed based on three clusters
- fChisq = 0;
- fX = points1[pos1].fX;
-
- fP0 = points1[pos1].fY;
- fP1 = points1[pos1].fZ;
-
- Float_t X2 = points2[pos2].fX;
- Float_t Y2 = points2[pos2].fY;
- Float_t Z2 = points2[pos2].fZ;
-
- Float_t X3 = points3[pos3].fX;
- Float_t Y3 = points3[pos3].fY;
- Float_t Z3 = points3[pos3].fZ;
-
- Float_t ZZ = fP1 - ((fP1 - Z3)/(fX-X3))*(fX-X2);
- if (TMath::Abs(ZZ - Z2) > 10) return 0; //What's this?? (fP1 - Z3)/(fX-X3)*(fX-X2) is an angle
-
- // It may make no difference. Check on a big event??.
- if ((X2-fX)*(0-Y2)-(0-X2)*(Y2-fP0) == 0) return 0; //Straight seed
-
- // Initial approximation of the state vector
- fP2 = f2(fX,fP0,X2,Y2,X3,Y3);
- fP3 = f3(fX,fP0,X2,Y2,fP1,Z2);
- fP4 = f4(fX,fP0,X2,Y2,X3,Y3);
- if (TMath::Abs(fP3) > 1.2) return 0; //From AliTPCtrackerMI
- if (TMath::Abs(fP4) >= 0.0066) return 0; // From AliTPCtrackerMI
-
- //cout << "AliL3KalmanTrack::MakeTrackSeed, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
- //cout << fP2 << " " << fP4 << endl;
-
- // Initial approximation of the covariance matrix
- Float_t sy1=points1[pos1].fSigmaY2;
- Float_t sz1=points1[pos1].fSigmaZ2;
- Float_t sy2=points2[pos2].fSigmaY2;
- Float_t sz2=points2[pos2].fSigmaZ2;
- Float_t sy3=25000*fP4*fP4+0.1;
- Float_t sy=0.1;
- Float_t sz=0.1;
-
- Float_t f40=(f4(fX,fP0+sy,X2,Y2,X3,Y3)-fP4)/sy;
- Float_t f42=(f4(fX,fP0,X2,Y2+sy,X3,Y3)-fP4)/sy;
- Float_t f43=(f4(fX,fP0,X2,Y2,X3,Y3+sy)-fP4)/sy;
- Float_t f20=(f2(fX,fP0+sy,X2,Y2,X3,Y3)-fP2)/sy;
- Float_t f22=(f2(fX,fP0,X2,Y2+sy,X3,Y3)-fP2)/sy;
- Float_t f23=(f2(fX,fP0,X2,Y2,X3,Y3+sy)-fP2)/sy;
- Float_t f30=(f3(fX,fP0+sy,X2,Y2,fP1,Z2)-fP3)/sy;
- Float_t f31=(f3(fX,fP0,X2,Y2,fP1+sz,Z2)-fP3)/sz;
- Float_t f32=(f3(fX,fP0,X2,Y2+sy,fP1,Z2)-fP3)/sy;
- Float_t f34=(f3(fX,fP0,X2,Y2,fP1,Z2+sz)-fP3)/sz;
-
- fC00 = sy1;
- fC10 = 0;
- fC11 = sz1;
- fC20 = f20*sy1;
- fC21 = 0;
- fC22 = f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
- fC30 = f30*sy1;
- fC31 = f31*sz1;
- fC32 = f30*sy1*f20+f32*sy2*f22;
- fC33 = f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
- fC40 = f40*sy1;
- fC41 = 0;
- fC42 = f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
- fC43 = f30*sy1*f40+f32*sy2*f42;
- fC44 = f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
- //cout << "Errors: " << fC00 << " " << fC11 << " " << fC20 << " " << fC22 << " " << " " << fC30 << " " << fC31 << " " << fC32 << " " << fC33 << " " << fC40 << " " << fC42 << " " << fC43 << " " << fC44 << endl;
-
- return 1;
-}
-
-Int_t AliL3KalmanTrack::Propagate(AliL3SpacePointData *points, UInt_t pos)
+Int_t AliL3KalmanTrack::Propagate(AliL3SpacePointData *points, UInt_t pos, Int_t slice)
{
// Propagates track to the plane of the next found cluster
Float_t Xold = fX; // X position for previous space point
- Float_t Xnew = points[pos].fX; // X position of current space point
- Float_t dx = Xnew - Xold;
+ //Float_t Xnew = points[pos].fX; // X position of current space point
+ //Float_t dx = Xnew - Xold; cout << Xnew << endl;
Float_t Yold = fP0; // Y position of old point
Float_t Zold = fP1; // Z position of old point
Float_t par2old = fP2;
Float_t oldfC43 = fC43;
Float_t oldfC44 = fC44;
+ Float_t xyz[3];
+
+ xyz[0] = points[pos].fX;
+ xyz[1] = points[pos].fY;
+
+ AliL3Transform::Global2Local(xyz,slice);
+
+ Float_t Xnew = xyz[0];
+ Float_t dx = Xnew - Xold; //cout << Xnew << endl;
+
if (TMath::Abs(fP4*Xnew - fP2) >= 0.9) // What's this??
{
//cout << "Propagation failed! Stiff track. " << pos << endl;
Float_t f02 = -dx*(2*RR + CC*(Cold/Rold + Cnew/Rnew))/(RR*RR);
Float_t f04 = dx*(RR*XX + CC*(Cold*Xold/Rold + Cnew*Xnew/Rnew))/(RR*RR);
Float_t CR = Cold*Rnew + Cnew*Rold;
+
Float_t f12 = -dx*fP3*(2*CR + CC*(Cnew*(Cold/Rold)-Rold + Cold*(Cnew/Rnew)-Rnew))/(CR*CR);
Float_t f13 = dx*CC/CR;
Float_t f14 = dx*fP3*(CR*XX-CC*(Rold*Xnew-Cnew*Cold*Xold/Rold + Rnew*Xold-Cold*Cnew*Xnew/Rnew))/(CR*CR);
fP2 += Xnew*(fP4-CC);
*/
// Update the track parameters with the measured values of the new point
- Int_t value = UpdateTrack(points, pos);
- /*if (value == 0)
+ Int_t value = UpdateTrack(points, pos, slice);
+
+ if (value == 0)
{
fP0 = Yold;
fP1 = Zold;
fC44 = oldfC44;
return value;
- }*/
+ }
- //else
-return value;
+ else
+ return value;
//return UpdateTrack(points, pos);
}
-Int_t AliL3KalmanTrack::UpdateTrack(AliL3SpacePointData *points, UInt_t pos)
+Int_t AliL3KalmanTrack::UpdateTrack(AliL3SpacePointData *points, UInt_t pos, Int_t slice)
{
// Update the track parameters with the measured values
- fX = points[pos].fX;
+
+ Float_t xyz[3];
+
+ xyz[0] = points[pos].fX;
+ xyz[1] = points[pos].fY;
+
+ AliL3Transform::Global2Local(xyz,slice);
+
+ //fX = points[pos].fX;
+ fX = xyz[0];
// The errors from the measurement of the spacepoint
Float_t sigmaY2 = points[pos].fSigmaY2;
Float_t k41 = fC40*sigmaYZ + fC41*sigmaZ2;
// Deviation between the predicted and measured values of y and z
- Float_t dy = points[pos].fY-fP0; //cout << "dy = " << dy;
+ //Float_t dy = points[pos].fY-fP0; //cout << "dy = " << dy;
+ Float_t dy = xyz[1] - fP0; //cout << "dy = " << dy;
Float_t dz = points[pos].fZ-fP1; //cout << ", dz = " << dz << endl;
- //cout << "Measured " << points[pos].fY << " " << points[pos].fZ << endl;
+ //cout << "Measured " << xyz[2] << " " << points[pos].fZ << endl;
// Prediction of fP2 and fP4
Float_t cur = fP4 + k40*dy + k41*dz;
//cout << "AliL3KalmanTrack::Update, Chi2 = " << GetChisq() << endl;
//cout << "AliKalmanTrack::Update, sigmaY2 = " << sigmaY2 << " sigmaZ2 = " << sigmaZ2 << " sigmaYZ = " << sigmaYZ << " dy = " << dy << " dz = " << dz << endl;
// Calculate increase of chisquare
- fChisq = GetChisq() + GetChisqIncrement(points,pos);
+ fChisq = GetChisq() + GetChisqIncrement(xyz[1],sigmaY2,points[pos].fZ,sigmaZ2);
//cout << "fChisq = " << fChisq << endl;
//(dy*sigmaY2*dy + 2*sigmaYZ*dy*dz + dz*sigmaZ2*dz) / (sigmaY2*sigmaZ2 - sigmaYZ*sigmaYZ);
// Must at some point make an cut on chisq. Here?
- // if (fChisq > fMaxChi2) return 0;
+ //if (fChisq > fMaxChi2) return 0;
return 1;
}
-Float_t AliL3KalmanTrack::GetChisqIncrement(AliL3SpacePointData *points, UInt_t pos)
+//Float_t AliL3KalmanTrack::GetChisqIncrement(AliL3SpacePointData *points, UInt_t pos)
+Float_t AliL3KalmanTrack::GetChisqIncrement(Float_t y, Float_t error_y, Float_t z, Float_t error_z)
{
// This function calculates a predicted chi2 increment.
//-----------------------------------------------------------------
- Float_t r00=points[pos].fSigmaY2, r01=0., r11=points[pos].fSigmaZ2;
+ //Float_t r00=points[pos].fSigmaY2, r01=0., r11=points[pos].fSigmaZ2;
+ Float_t r00=error_y, r01=0., r11=error_z;
r00+=fC00; r01+=fC10; r11+=fC11;
Double_t det=r00*r11 - r01*r01;
}*/
Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
- Double_t dy=points[pos].fY - fP0, dz=points[pos].fZ - fP1;
+ Double_t dy=y - fP0, dz=z - fP1;
//cout << "AliTPCtrack::GetPredictedChi2, r00 = " << r00 << " r11 = " << r11 << " ro1 = " << r01 << " dy = " << dy << " dz = " << dz << endl;
return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
}
SetNHits(tpt->GetNHits());
}
+
+Int_t AliL3KalmanTrack::PropagateOfflineTrack(Double_t x, Double_t y, Double_t z, Double_t ey, Double_t ez)
+{
+ // Propagates track to the plane of the next found cluster
+ Float_t Xold = fX; // X position for previous space point
+ Float_t Xnew = x; // X position of current space point
+ Float_t dx = Xnew - Xold; //cout << dx << endl;
+ Float_t Yold = fP0; // Y position of old point
+ Float_t Zold = fP1; // Z position of old point
+ Float_t par2old = fP2;
+ Float_t par3old = fP3;
+ Float_t par4old = fP4;
+ Float_t oldfC00 = fC00;
+ Float_t oldfC10 = fC10;
+ Float_t oldfC11 = fC11;
+ Float_t oldfC20 = fC20;
+ Float_t oldfC21 = fC21;
+ Float_t oldfC22 = fC22;
+ Float_t oldfC30 = fC30;
+ Float_t oldfC31 = fC31;
+ Float_t oldfC32 = fC32;
+ Float_t oldfC33 = fC33;
+ Float_t oldfC40 = fC40;
+ Float_t oldfC41 = fC41;
+ Float_t oldfC42 = fC42;
+ Float_t oldfC43 = fC43;
+ Float_t oldfC44 = fC44;
+
+ if (TMath::Abs(fP4*Xnew - fP2) >= 0.9) // What's this??
+ {
+ //cout << "Propagation failed! Stiff track. " << pos << endl;
+ return 0;
+ }
+
+ if (TMath::Abs(dx) > 5) return 0;
+
+ /*if (TMath::Abs(fP4*Xold - fP2) >= 0.9) // What's this??
+ {
+ return 0;
+ }*/
+
+ // R must be approx. of the radius of the circle based on
+ // the old and new spacepoint. What then is Cod and Cnew??
+ // C has something to do with curvature at least.
+
+ Float_t Cold = fP4*Xold - fP2;
+
+ //if (TMath::Abs(Cold) >= 1.0) return 0;
+ /*{
+ cout << "Cold " << endl << fP4*Xnew - fP2 << endl;
+ return 0;
+ }*/
+ Float_t Rold = TMath::Sqrt(1 - Cold*Cold);
+ Float_t Cnew = fP4*Xnew - fP2;
+
+ //if (TMath::Abs(Cnew) >= 1.0) Cnew = 0.9;
+ /*{
+ cout << "Cnew " << endl;
+ return 0;
+ }*/
+ Float_t Rnew = TMath::Sqrt(1 - Cnew*Cnew);
+ //if (Rold < 0.9)
+ //cout << "Cold = " << Cold << " , Rold = " << Rold << " , Cnew = " << Cnew << " , Rnew = " << Rnew << endl;
+ // Prediction of the y- and z- coordinate in the next plane
+ fP0 += dx*(Cold+Cnew)/(Rold+Rnew);
+ fP1 += dx*(Cold+Cnew)/(Cold*Rnew + Cnew*Rold)*fP3;
+ //cout << "Old point " << Yold << " " << Zold << endl;
+ //cout << "Propagate " << fP0 << " " << fP1 << endl;
+ //cout << "Measured " << points[pos].fY << " " << points[pos].fZ << endl;
+
+ fX = Xnew;
+
+ // f = F - 1 //What is this??
+ // Must be the f-matrix for the prediction, as in eq 1 in ALICE Kalman paper
+ Float_t RR = Rold + Rnew;
+ Float_t CC = Cold + Cnew;
+ Float_t XX = Xold + Xnew;
+
+ Float_t f02 = -dx*(2*RR + CC*(Cold/Rold + Cnew/Rnew))/(RR*RR);
+ Float_t f04 = dx*(RR*XX + CC*(Cold*Xold/Rold + Cnew*Xnew/Rnew))/(RR*RR);
+ Float_t CR = Cold*Rnew + Cnew*Rold;
+ Float_t f12 = -dx*fP3*(2*CR + CC*(Cnew*(Cold/Rold)-Rold + Cold*(Cnew/Rnew)-Rnew))/(CR*CR);
+ Float_t f13 = dx*CC/CR;
+ Float_t f14 = dx*fP3*(CR*XX-CC*(Rold*Xnew-Cnew*Cold*Xold/Rold + Rnew*Xold-Cold*Cnew*Xnew/Rnew))/(CR*CR);
+
+ // b = C*ft // This?
+ Float_t b00=f02*fC20 + f04*fC40;
+ Float_t b01=f12*fC20 + f14*fC40 + f13*fC30;
+ Float_t b10=f02*fC21 + f04*fC41;
+ Float_t b11=f12*fC21 + f14*fC41 + f13*fC31;
+ Float_t b20=f02*fC22 + f04*fC42;
+ Float_t b21=f12*fC22 + f14*fC42 + f13*fC32;
+ Float_t b30=f02*fC32 + f04*fC43;
+ Float_t b31=f12*fC32 + f14*fC43 + f13*fC33;
+ Float_t b40=f02*fC42 + f04*fC44;
+ Float_t b41=f12*fC42 + f14*fC44 + f13*fC43;
+
+ //a = f*b = f*C*ft
+ Float_t a00 = f02*b20 + f04*b40;
+ Float_t a01 = f02*b21 + f04*b41;
+ Float_t a11 = f12*b21 + f14*b41+f13*b31;
+
+ //F*C*Ft = C + (a + b + bt) /This is the covariance matrix, the samll t
+ // means transform. Then F must be df/dx
+ fC00 += a00 + 2*b00;
+ fC10 += a01 + b01 + b10;
+ fC20 += b20;
+ fC30 += b30;
+ fC40 += b40;
+ fC11 += a11 + 2*b11;
+ fC21 += b21;
+ fC31 += b31;
+ fC41 += b41;
+
+ // Update the track parameters with the measured values of the new point
+ Int_t value = UpdateOfflineTrack(x,y,z,ey,ez);
+
+ if (value == 0)
+ {
+ fP0 = Yold;
+ fP1 = Zold;
+ fP2 = par2old;
+ fP3 = par3old;
+ fP4 = par4old;
+ fC00 = oldfC00;
+ fC10 = oldfC10;
+ fC11 = oldfC11;
+ fC20 = oldfC20;
+ fC21 = oldfC21;
+ fC22 = oldfC22;
+ fC30 = oldfC30;
+ fC31 = oldfC31;
+ fC32 = oldfC32;
+ fC33 = oldfC33;
+ fC40 = oldfC40;
+ fC41 = oldfC41;
+ fC42 = oldfC42;
+ fC43 = oldfC43;
+ fC44 = oldfC44;
+
+ return value;
+ }
+
+ else
+ return value;
+ //return 1;
+ //return UpdateTrack(points, pos);
+
+}
+
+Int_t AliL3KalmanTrack::UpdateOfflineTrack(Double_t x, Double_t y, Double_t z, Double_t ey, Double_t ez)
+{
+ // Update the track parameters with the measured values
+ //fX = x;
+
+ // The errors from the measurement of the spacepoint
+ Float_t sigmaY2 = ey;
+ Float_t sigmaZ2 = ez;
+ Float_t sigmaYZ = 0;
+
+ sigmaY2 += fC00;
+ sigmaZ2 += fC11;
+ sigmaYZ += fC10;
+
+ Float_t det = sigmaY2*sigmaZ2 - sigmaYZ*sigmaYZ;
+ Float_t tmp = sigmaY2;
+ sigmaY2 = sigmaZ2/det;
+ sigmaZ2 = tmp/det;
+ sigmaYZ = -sigmaYZ/det;
+
+ // Kalman gain matrix
+ Float_t k00 = fC00*sigmaY2 + fC10*sigmaYZ, k01 = fC00*sigmaYZ + fC10*sigmaZ2;
+ Float_t k10 = fC10*sigmaY2 + fC11*sigmaYZ, k11 = fC10*sigmaYZ + fC11*sigmaZ2;
+ Float_t k20 = fC20*sigmaY2 + fC21*sigmaYZ, k21 = fC20*sigmaYZ + fC21*sigmaZ2;
+ Float_t k30 = fC30*sigmaY2 + fC31*sigmaYZ, k31 = fC30*sigmaYZ + fC31*sigmaZ2;
+ Float_t k40 = fC40*sigmaY2 + fC41*sigmaYZ, k41 = fC40*sigmaYZ + fC41*sigmaZ2;
+ //cout << "x = " << fX << endl;
+ // Deviation between the predicted and measured values of y and z
+ Float_t dy = y-fP0; //cout << "dy = " << dy;
+ Float_t dz = z-fP1; //cout << ", dz = " << dz << endl;
+ //cout << "Measured " << points[pos].fY << " " << points[pos].fZ << endl;
+
+ // Prediction of fP2 and fP4
+ Float_t cur = fP4 + k40*dy + k41*dz;
+ Float_t eta = fP2 + k20*dy + k21*dz;
+
+ if (TMath::Abs(cur*fX - eta) >= 0.9)
+ {
+ //cout << "Update failed! Stiff track. " << pos << endl;
+ return 0;
+ }
+
+ // Filtered state vector
+ fP0 += k00*dy + k01*dz;
+ fP1 += k10*dy + k11*dz;
+ fP2 = eta;
+ fP3 += k30*dy + k31*dz; //cout << "update " << fP3 << endl;
+ fP4 = cur;
+ //cout << "AliL3KalmanTrack::Update, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
+ //cout << "Measured, " << points[pos].fY << " " << points[pos].fZ << endl;
+
+ Float_t c01 = fC10;
+ Float_t c02 = fC20;
+ Float_t c03 = fC30;
+ Float_t c04 = fC40;
+ Float_t c12 = fC21;
+ Float_t c13 = fC31;
+ Float_t c14 = fC41;
+
+ // Filtered covariance matrix
+ fC00 -= k00*fC00 + k01*fC10; fC10 -= k00*c01 + k01*fC11;
+ fC20 -= k00*c02 + k01*c12; fC30 -= k00*c03 + k01*c13;
+ fC40 -= k00*c04 + k01*c14;
+
+ fC11 -= k10*c01 + k11*fC11;
+ fC21 -= k10*c02 + k11*c12; fC31 -= k10*c03 + k11*c13;
+ fC41 -= k10*c04 + k11*c14;
+
+ fC22 -= k20*c02 + k21*c12; fC32 -= k20*c03 + k21*c13;
+ fC42 -= k20*c04 + k21*c14;
+
+ fC33 -= k30*c03 + k31*c13;
+ fC43 -= k40*c03 + k41*c13;
+
+ fC44 -= k40*c04 + k41*c14;
+
+ //cout << "AliL3KalmanTrack::Update, error " << fC00 << " " << fC11 << " " << fC22 << " " << fC33 << " " << fC44 << endl;
+
+ sigmaY2 = sigmaY2*det;
+ sigmaZ2 = sigmaZ2*det;
+ sigmaYZ = sigmaYZ*det;
+ //cout << "AliL3KalmanTrack::Update, Chi2 = " << GetChisq() << endl;
+
+ fChisq = GetChisq() + GetChisqIncrementOfflineTrack(y,z,ey,ez);
+ //cout << "fChisq = " << fChisq << endl;
+
+ // Must at some point make an cut on chisq. Here?
+ //if (fChisq > fMaxChi2) return 0;
+
+ return 1;
+
+}
+
+Float_t AliL3KalmanTrack::GetChisqIncrementOfflineTrack(Double_t y, Double_t z, Double_t ey, Double_t ez)
+{
+ // This function calculates a predicted chi2 increment.
+ //-----------------------------------------------------------------
+ Float_t r00=ey, r01=0., r11=ez;
+ r00+=fC00; r01+=fC10; r11+=fC11;
+
+ Double_t det=r00*r11 - r01*r01;
+ /*if (TMath::Abs(det) < 1.e-10) {
+ Int_t n=GetNumberOfClusters();
+ if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+ return 1e10;
+ }*/
+ Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+
+ Double_t dy=y - fP0, dz=z - fP1;
+ //cout << "AliTPCtrack::GetPredictedChi2, r00 = " << r00 << " r11 = " << r11 << " ro1 = " << r01 << " dy = " << dy << " dz = " << dz << endl;
+ return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
+}
#include "AliL3RootTypes.h"
#include "AliL3Track.h"
+
+// includes for offline comparison, will be removed
+#include "AliTPCtrack.h"
+// includes for offline comparison, will be removed
+
class AliL3SpacePointData;
class AliL3KalmanTrack : public AliL3Track {
AliL3KalmanTrack();
virtual ~AliL3KalmanTrack();
- Int_t Init(AliL3Track *track, AliL3SpacePointData *points, UInt_t pos);
- Int_t MakeTrackSeed(AliL3SpacePointData *points1, UInt_t pos1, AliL3SpacePointData *points2, UInt_t pos2, AliL3SpacePointData *points3, UInt_t pos3);
- Int_t Propagate(AliL3SpacePointData *points, UInt_t pos);
- Int_t UpdateTrack(AliL3SpacePointData *points, UInt_t pos);
+ Int_t Init(AliL3Track *track, AliL3SpacePointData *points, UInt_t pos,Int_t slice);
+ Int_t Propagate(AliL3SpacePointData *points, UInt_t pos, Int_t slice);
+ Int_t UpdateTrack(AliL3SpacePointData *points, UInt_t pos, Int_t slice);
Int_t UpdateTrackII(AliL3SpacePointData *points, UInt_t pos);
void AddTrack();
// Float_t GetStateVector(Float_t xx[5]) const {
Float_t GetC13() {return fC43;}
Float_t GetC14() {return fC44;}
+ void SetX(Float_t f) {fX = f;}
void SetX0(Float_t f) {fP0 = f;}
void SetX1(Float_t f) {fP1 = f;}
void SetX2(Float_t f) {fP2 = f;}
fC44 = f[14];}
void SetChisq(Float_t f) {fChisq = f;}
void SetMaxChi2(Float_t f) {fMaxChi2 = f;}
- Float_t GetChisqIncrement(AliL3SpacePointData *points, UInt_t pos);
+ //Float_t GetChisqIncrement(AliL3SpacePointData *points, UInt_t pos);
+ Float_t GetChisqIncrement(Float_t y, Float_t error_y, Float_t z, Float_t error_z);
Float_t f2(Float_t x1,Float_t y1, Float_t x2,Float_t y2, Float_t x3,Float_t y3);
Float_t f3(Float_t x1,Float_t y1, Float_t x2,Float_t y2, Float_t z1,Float_t z2);
Float_t f4(Float_t x1,Float_t y1, Float_t x2,Float_t y2, Float_t x3,Float_t y3);
Int_t GetNHits() const {return fNHits;}
void SetNHits(Int_t f) {fNHits = f;}
+
+ Int_t PropagateOfflineTrack(Double_t x, Double_t y, Double_t z, Double_t ey, Double_t ez);
+ Int_t UpdateOfflineTrack(Double_t x, Double_t y, Double_t z, Double_t ey, Double_t ez);
+ Float_t GetChisqIncrementOfflineTrack(Double_t y, Double_t z, Double_t ey, Double_t ez);
};
#endif
#endif
#include "AliL3DDLDataFileHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3DDLRawReaderFile.h"
+#if __GNUC__ == 3
+using namespace std;
+#endif
+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
</pre>
*/
-
ClassImp(AliL3DDLRawReaderFile)
-
AliL3DDLRawReaderFile::AliL3DDLRawReaderFile(const Char_t* name, Bool_t addnum)
{
// create an object to read digits from the given input file(s)
#include "AliL3Transform.h"
#include "AliL3DataHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3StandardIncludes.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3Transform.h"
#include "AliL3TPCMapping.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3TransBit.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
//#define VHDLDEBUG
#include "AliL3VHDLClusterFinder.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
CXX = g++
LD = $(CXX)
-GCCVERSION = $(shell $(CXX) --version | head -n 1 | cut -d" " -f 3 | cut -d. -f 1 | cut -d" " -f1)
-
-CXXGCC3FLAGS = -DGCCVERSION=$(GCCVERSION)
#CXXGCC3FLAGS += -pedantic
#CXXGCC3FLAGS += -Wno-deprecated
#CXXGCC3FLAGS += -Woverloaded-virtual
-
CXXFLAGS = -O2 -Wall -ggdb $(CXXGCC3FLAGS) $(EXTRACXXFLAGS)
LDFLAGS = -O2 $(EXTRALDFLAGS) $(LIBS)
STATICLDFLAGS = -O2 $(PROFILEFLAGS) $(EXTRALDFLAGS) $(STATICLIBS)
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
#include "AliL3StandardIncludes.h"
#include "AliL3RootTypes.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include <AliL3MemHandler.h>
#include <AliL3HoughTransformerVhdl.h>
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
}
AliL3Transform::Init(path);
- //cerr << "Transform version: " << AliL3Transform::GetVersion() << endl;
+ cerr << "Transform version: " << AliL3Transform::GetVersion() << endl;
AliL3HoughTransformerVhdl vtest(slice,patch,100,10);
vtest.CreateHistograms(64,0.1,64,-30,30);
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
#include <TGraphErrors.h>
#endif
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3MemHandler.h"
#include "AliL3SpacePointData.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
//Standalone program to run the track follower for benchmark tests.
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
#include "AliL3AltroMemHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
#include "AliL3Logging.h"
#include "AliL3Logger.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#else
#include <stream.h>
#include "AliL3SpacePointData.h"
#include "AliL3MemHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
fTrack->SetCharge(q);
-//
-// Get other track parameters
-//
- Double_t x0, y0,phi0,r0,psi,pt ;
- if ( fTrack->ComesFromMainVertex() == true )
- {
- //flag = 1 ; // primary track flag
- x0 = fVertex->GetX() ;
- y0 = fVertex->GetY() ;
- phi0 = fVertex->GetPhi() ;
- r0 = fVertex->GetR() ;
- fTrack->SetPhi0(phi0);
- fTrack->SetR0(r0);
- }
- else
- {
- //AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)hits->Last();
- AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)fTrack->lastHit;
- //flag = 0 ; // primary track flag
- x0 = lHit->GetX() ;
- y0 = lHit->GetY() ;
- phi0 = atan2(lHit->GetY(),lHit->GetX());
- if ( phi0 < 0 ) phi0 += AliL3Transform::TwoPi();
- r0 = sqrt ( lHit->GetX() * lHit->GetX() + lHit->GetY() * lHit->GetY() ) ;
- fTrack->SetPhi0(phi0);
- fTrack->SetR0(r0);
-
- }
//Set the first point on the track to the space point coordinates of the innermost track
//This will be updated to lie on the fit later on (AliL3Track::UpdateToFirstPoint).
+ Double_t x0,y0,psi,pt ;
+ AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)fTrack->lastHit;
+ x0 = lHit->GetX();
+ y0 = lHit->GetY();
fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
psi = (Double_t)atan2(bcent-y0,acent-x0) ;
//TObjArray *hits = fTrack->GetHits();
//Int_t num_of_hits = fTrack->GetNumberOfPoints();
- if ( fTrack->ComesFromMainVertex() == true )
+ if (0)// fTrack->ComesFromMainVertex() == true )
{
dx = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetX() - fVertex->GetX();
dy = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetY() - fVertex->GetY() ;
AliL3ConfMapPoint *fHit = (AliL3ConfMapPoint*)firstHit;
SetLastPoint(fHit->GetX(),fHit->GetY(),fHit->GetZ());
- if(AliLevel3::IsTracksAtFirstPoint())//Set the track fit parameter at the innermost most point
- {
- UpdateToFirstPoint();
- }
+ UpdateToFirstPoint();
delete fit;
}
#include "AliL3SpacePointData.h"
#include "AliL3MemHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include <AliComplexCluster.h>
#include <AliStack.h>
-#if GCCVERSION == 3
+#if __GNUC__ == 3
#include <fstream>
#include <iosfwd>
#else
#include "AliL3TrackArray.h"
#include "AliL3Evaluate.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#include "AliL3TrackArray.h"
#include "AliL3FileHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
fTrack->SetCharge(q);
-//
-// Get other track parameters
-//
- Double_t x0, y0,phi0,r0,psi,pt ;
- if ( fVertexConstraint == kTRUE)
- {
- //flag = 1 ; // primary track flag
- x0 = fVertex->GetX() ;
- y0 = fVertex->GetY() ;
- phi0 = fVertex->GetPhi() ;
- r0 = fVertex->GetR() ;
- fTrack->SetPhi0(phi0);
- fTrack->SetR0(r0);
- }
- else
- {
- Int_t lastid=fTrack->GetNHits()-1;
- UInt_t id = hitnum[lastid];
- Int_t slice = (id>>25) & 0x7f;
- Int_t patch = (id>>22) & 0x7;
- UInt_t pos = id&0x3fffff;
- AliL3SpacePointData *points = fClusters[slice][patch];
-
- //flag = 0 ; // primary track flag
- x0 = points[pos].fX;
- y0 = points[pos].fY;
- phi0 = atan2(points[pos].fY,points[pos].fX);
- if ( phi0 < 0 ) phi0 += 2*AliL3Transform::Pi();
- r0 = sqrt ( points[pos].fX * points[pos].fX + points[pos].fY*points[pos].fY);
- fTrack->SetPhi0(phi0);
- fTrack->SetR0(r0);
- }
-
//Set the first point on the track to the space point coordinates of the innermost track
//This will be updated to lie on the fit later on (AliL3Track::UpdateToFirstPoint).
+ Double_t x0,y0,psi,pt ;
+ Int_t lastid=fTrack->GetNHits()-1;
+ UInt_t id = hitnum[lastid];
+ Int_t slice = (id>>25) & 0x7f;
+ Int_t patch = (id>>22) & 0x7;
+ UInt_t pos = id&0x3fffff;
+ AliL3SpacePointData *points = fClusters[slice][patch];
+ x0 = points[pos].fX;
+ y0 = points[pos].fY;
fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
-
+
+ //Set the remaining fit parameters
psi = (Double_t)atan2(bcent-y0,acent-x0) ;
psi = psi + q * 0.5F * AliL3Transform::Pi() ;
if ( psi < 0 ) psi = psi + 2*AliL3Transform::Pi();
Double_t fS[(fTrack->GetNHits())];
Double_t *fZWeight = new Double_t[fTrack->GetNHits()];
UInt_t *hitnum = fTrack->GetHitNumbers();
- if (fVertexConstraint==kTRUE)
+ if (0)//fVertexConstraint==kTRUE)
{
UInt_t id = hitnum[0];
Int_t slice = (id>>25) & 0x7f;
MLUCLogServer *so; //!
MLUCLogServer *se; //!
MLUCLogServer *sm; //!
-#if GCCVERSION == 3
+#if __GNUC__ == 3
std::ofstream *of; //!
#else
ofstream *of; //!
enum TLogCmd { kEnd, kPrec, kHex, kDec };
};
-#if GCCVERSION == 3
+#if __GNUC__ == 3
#define LOG( lvl, origin, keyword ) std::cerr<<"["<<origin<<": "<<keyword<<"] "
#define ENDLOG std::endl
#else
#include "AliL3RawDataFileHandler.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
#else
-#if GCCVERSION == 3
+#if __GNUC__ == 3
#include <cstdio>
#include <cmath>
#else
#ifndef ALIL3STANDARDINCLUDESH
#define ALIL3STANDARDINCLUDESH
-#if GCCVERSION == 3
+#if __GNUC__ == 3
#include <iostream>
#include <fstream>
#define STDCERR cerr
#define STDENDL endl
-#endif //GCCVERSION
+#endif //__GNUC__
#endif
#include "AliL3Transform.h"
#include "AliL3Vertex.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
track->SetFirstPoint(trs->fX,trs->fY,trs->fZ);
track->SetLastPoint(trs->fLastX,trs->fLastY,trs->fLastZ);
track->SetHits( trs->fNPoints, trs->fPointIDs );
+#ifdef ROWHOUGH
+ if(GetTrackType()=='h') {
+ ((AliL3HoughTrack *)track)->SetWeight(trs->fWeight);
+ track->SetMCid(trs->fTrackID);
+ }
+#endif
UChar_t *tmpP = (UChar_t*)trs;
tmpP += sizeof(AliL3TrackSegmentData)+trs->fNPoints*sizeof(UInt_t);
trs = (AliL3TrackSegmentData*)tmpP;
UChar_t *tmpP = (UChar_t*)trs;
tmpP += sizeof(AliL3TrackSegmentData)+trs->fNPoints*sizeof(UInt_t);
trs = (AliL3TrackSegmentData*)tmpP;
+#ifdef ROWHOUGH
+ if(GetTrackType()=='h') {
+ ((AliL3HoughTrack *)track)->SetWeight(trs->fWeight);
+ track->SetMCid(trs->fTrackID);
+ }
+#endif
}
}
tP->fTgl = track->GetTgl();
tP->fCharge = track->GetCharge();
tP->fNPoints = track->GetNHits();
-
+#ifdef ROWHOUGH
+ if(GetTrackType()=='h') {
+ tP->fWeight = ((AliL3HoughTrack *)track)->GetWeight();
+ tP->fTrackID = track->GetMCid();
+ }
+#endif
pP = (UInt_t*)track->GetHitNumbers();
for (UInt_t j=0;j<tP->fNPoints;j++){
tP->fPointIDs[j] = pP[j];
Double_t fPsi;
Double_t fTgl;
Int_t fCharge;
+#ifdef ROWHOUGH
+ UInt_t fWeight;
+ Int_t fTrackID;
+#endif
UInt_t fNPoints;
UInt_t fPointIDs[0];
};
#include "AliL3Logging.h"
#include "AliL3Transform.h"
-#if GCCVERSION == 3
+#if __GNUC__ == 3
using namespace std;
#endif
</pre>
*/
-
ClassImp(AliL3Transform)
-
const Double_t AliL3Transform::fAnodeWireSpacing = 0.25; //Taken from the TDR
const Double_t AliL3Transform::fBFACT = 0.0029980; //Conversion Factor
const Double_t AliL3Transform::fPi = 3.141592653589793;
void AliL3Transform::PrintCompileOptions()
{
- cout << "Compiler (g++) version used: " << GCCVERSION << endl;
+ cout << "Compiler (g++) version used: " << __GNUC__ << endl;
#ifdef no_root
cout << "STANDALONE version: -Dno_root was given." << endl;
cout << "NOT using Monte Carlo Info: -Ddo_mc was not given." << endl;
#endif
+#ifdef ROWHOUGH
+ cout << "Using row transformer version: -DROWHOUGH was given." << endl;
+#else
+ cout << "NOT using row transformer version: -DROWHOUGH was not given." << endl;
+#endif
+
+#ifdef use_newio
+ cout << "Using NEWIO version: -Duse_newio was given." << endl;
+#else
+ cout << "NOT using NEWIO version: -Duse_newio was not given." << endl;
+#endif
+
#ifdef use_logging
cout << "Using logging classes (MLUC): -Duse_logging was given." << endl;
#else
cout << "NOT using logging classes (MLUC): -Duse_logging not was given." << endl;
#endif
+
}
ClassImp(AliLevel3)
-Bool_t AliLevel3::fSetTracks2FirstPoint = kTRUE;//Define track parameters at first point
Bool_t AliLevel3::fDoVertexFit = kTRUE;//Include the vertex in the final track fit
AliLevel3::AliLevel3()
AliL3Track *tr = tracks->GetCheckedTrack(i);
if(!tr) continue;
fitter->FitHelix(tr);
- if(AliLevel3::IsTracksAtFirstPoint())
- tr->UpdateToFirstPoint();
+ tr->UpdateToFirstPoint();
}
fBenchmark->Stop("Global track fitter");
delete fitter;
Bool_t fUseBinary;
Bool_t fWriteOut;
- //Define whether track parameters should be given at first point on track (default)
- //If not, the parameters will be given at the vertex.
- static Bool_t fSetTracks2FirstPoint;
static Bool_t fDoVertexFit;
Bool_t fClusterDeconv;
void DoRoi(Float_t e0=0.4,Float_t e1=0.5){fEta[0]=e0;fEta[1]=e1;fDoRoi=kTRUE;}
void WriteFiles(Char_t *path="./"){fWriteOut = kTRUE; sprintf(fWriteOutPath,"%s",path);}
- static void SetTracks2FirstPoint() {fSetTracks2FirstPoint = kTRUE;}
- static void SetTracks2Vertex() {fSetTracks2FirstPoint = kFALSE;}
static void SetVertexFit(Bool_t f) {fDoVertexFit=f;}
- static Bool_t IsTracksAtFirstPoint() {return fSetTracks2FirstPoint;}
static Bool_t DoVertexFit() {return fDoVertexFit;}
ClassDef(AliLevel3,1) //Interface class for Level3-tracking
#include "AliITStrackV2.h"
#include "AliL3Transform.h"
#include <TVector3.h>
+#include <iostream.h>
ClassImp(AliD0Trigger)
posTrack=0;
negTrack=0;
}
-AliD0Trigger::AliD0Trigger(double c[4],double b,double pv[3])
+AliD0Trigger::AliD0Trigger(double c[6],double b,double pv[3])
{
posTrack=0;
negTrack=0;
cutV0high=c[1];
cutInvMass=c[2];
cutPointAngle=c[3];
+ cutd0d0=c[4];
+ cutCosThetaStar=c[5];
+ cutpTchild=c[6];
Bfield=b;
primaryVertex[0]=pv[0];
primaryVertex[1]=pv[1];
}
bool AliD0Trigger::FindV0(){
- bool goodV0=false;
+ bool goodV0=false;
+
+ double Gxpos=0,Gypos=0,Gxneg=0,Gyneg=0;
+ double r1=0,r2=0,a1=0,a2=0,b1=0,b2=0;
+ double q=0,p=0,t=0,Fa=0,Fb=0,Fc=0;
+ double y1=0,y2=0,x1=0,x2=0,x3=0,x4=0;
+ double V01=0,V02=0,V03=0,V04=0;
+ double V0[4]={0,0,0,0};
+ bestV0[0]=0; bestV0[1]=0; bestV0[2]=0;
- double Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
- double Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
- double Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
- double Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
+ Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
+ Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
+ Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
+ Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
- double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
- double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
- //double a1=posTrack->GetX()-(r1*cos(asin(posTrack->GetSnp())+TMath::PiOver2()));
- //double a2=negTrack->GetX()-(r1*cos(asin(negTrack->GetSnp())-TMath::PiOver2()));
- //double b1=posTrack->GetY()-(r1*sin(asin(posTrack->GetSnp())+TMath::PiOver2()));
- //double b2=negTrack->GetY()-(r1*sin(asin(negTrack->GetSnp())-TMath::PiOver2()));
-
- double a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
- double a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
- double b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
- double b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
-
- double a3=a1-a2;
- double b3=b1-b2;
- double r3=sqrt(pow(a3,2)+pow(b3,2));
-
- if(r3<r1+r2){
- double q=a1-a2;
- double p=2*(b1-b2);
- double t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
-
- double Fa=pow(p,2)+(4*pow(q,2));
- double Fb=-(2*t*p+8*pow(q,2)*b1);
- double Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
-
- if(pow(Fb,2)-(4*Fa*Fc)>=0){
- double y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //noe feil her. floating point
- double y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //trolig negativ under rot
-
- double x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
- double x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
- double x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
- double x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
-
- double V01=sqrt(pow(x1,2)+pow(y1,2));
- double V02=sqrt(pow(x2,2)+pow(y2,2));
- double V03=sqrt(pow(x3,2)+pow(y1,2));
- double V04=sqrt(pow(x4,2)+pow(y2,2));
-
- double V0[4]={V01,V02,V03,V04};
- int index=0;
-
- double nearestV0=V0[0];
- for(int i=1;i<4;i++){
- if(nearestV0>V0[i]){
- nearestV0=V0[i];
- index=i;
- }
+ r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
+ r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
+
+ a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
+ a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
+ b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
+ b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
+
+ //double a3=a1-a2;
+ //double b3=b1-b2;
+ //double r3=sqrt(pow(a3,2)+pow(b3,2));
+
+ //if(r3<r1+r2){
+ q=a1-a2;
+ p=2*(b1-b2);
+ t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
+
+ Fa=pow(p,2)+(4*pow(q,2));
+ Fb=-(2*t*p+8*pow(q,2)*b1);
+ Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
+
+ if(pow(Fb,2)-(4*Fa*Fc)>=0){
+ y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //noe feil her. floating point
+ y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //trolig negativ under rot
+
+ x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
+ x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
+ x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
+ x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
+
+ V01=sqrt(pow(x1-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
+ V02=sqrt(pow(x2-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
+ V03=sqrt(pow(x3-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
+ V04=sqrt(pow(x4-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
+
+ V0[0]=V01;V0[1]=V02;V0[2]=V03;V0[3]=V04;
+
+ int index=0;
+
+ double nearestV0=V0[0];
+
+ for(int i=0;i<4;i++){
+ if(nearestV0>V0[i]){
+ nearestV0=V0[i];
+ index=i;
}
- if(nearestV0<cutV0high && nearestV0>cutV0low){
- goodV0=true;
- switch (index){
- case 0 : bestV0[0]=x1;bestV0[1]=y1;bestV0[0]=0.;break;
- case 1 : bestV0[0]=x2;bestV0[1]=y2;bestV0[0]=0.;break;
- case 2 : bestV0[0]=x3;bestV0[1]=y1;bestV0[0]=0.;break;
- case 3 : bestV0[0]=x4;bestV0[1]=y2;bestV0[0]=0.;break;
- default : printf("Can't set V0");
+ }
+ cout<<"index: "<<index<<endl;
+ if(nearestV0<cutV0high){
+ if(nearestV0>cutV0low){
+ goodV0=true;
+ if(index==0){
+ bestV0[0]=x1; bestV0[1]=y1; bestV0[2]=0.;
+ }
+ if(index==1){
+ bestV0[0]=x2; bestV0[1]=y2; bestV0[2]=0.;
+ }
+ if(index==2){
+ bestV0[0]=x3; bestV0[1]=y1; bestV0[2]=0.;
+ }
+ if(index==3){
+ bestV0[0]=x4; bestV0[1]=y2; bestV0[2]=0.;
}
}
}
+ //}
}
+
+ cout<<"My V0: x: "<<bestV0[0]<<" y: "<<bestV0[1]<<" z: "<<bestV0[2]<<endl;
+
return goodV0;
}
void AliD0Trigger::FindMomentaAtVertex(){
else
return false;
}
+void AliD0Trigger::SetMomenta(double m[6])
+{
+ momenta[0] = m[0];
+ momenta[1] = m[1];
+ momenta[2] = m[2];
+ momenta[3] = m[3];
+ momenta[4] = m[4];
+ momenta[5] = m[5];
+}
+void AliD0Trigger::FindMomentaOffline()
+{
+ double ptP,alphaP,phiP,ptN,alphaN,phiN;
+ // momenta of the tracks at the vertex
+ ptP = 1./TMath::Abs(posTrack->Get1Pt());
+ alphaP = posTrack->GetAlpha();
+ phiP = alphaP+TMath::ASin(posTrack->GetSnp());
+ momenta[0] = ptP*TMath::Cos(phiP);
+ momenta[1] = ptP*TMath::Sin(phiP);
+ momenta[2] = ptP*posTrack->GetTgl();
+
+ ptN = 1./TMath::Abs(negTrack->Get1Pt());
+ alphaN = negTrack->GetAlpha();
+ phiN = alphaN+TMath::ASin(negTrack->GetSnp());
+ momenta[3] = ptN*TMath::Cos(phiN);
+ momenta[4] = ptN*TMath::Sin(phiN);
+ momenta[5] = ptN*negTrack->GetTgl();
+}
+bool AliD0Trigger::FindV0offline(double v[3])
+{
+
+ bool goodV0=false;
+ double r=0;
+
+ bestV0[0]=v[0];
+ bestV0[1]=v[1];
+ bestV0[2]=v[2];
+
+ r=sqrt(pow(bestV0[0]-primaryVertex[0],2)+pow(bestV0[1]-primaryVertex[1],2)+pow(bestV0[2]-primaryVertex[2],2));
+
+ if(r<cutV0high){
+ if(r>cutV0low){
+ goodV0=true;
+ }
+ }
+
+ return goodV0;
+
+}
+bool AliD0Trigger::d0d0()
+{
+ bool goodd0=false;
+ double d00=0, d01=0;
+
+ d00 = 10000.*posTrack->GetD(primaryVertex[0],primaryVertex[1]);
+ d01 = -10000.*negTrack->GetD(primaryVertex[0],primaryVertex[1]);
+
+ if(d00*d01<cutd0d0){
+ goodd0=true;
+ }
+
+ return goodd0;
+}
+double AliD0Trigger::Energy()
+{
+ double kMD0 = 1.8645; // D0 mass
+ return sqrt(P()*P()+kMD0*kMD0);
+}
+bool AliD0Trigger::CosThetaStar()
+{
+ bool goodtheta=false;
+ double kMD0 = 1.8645; // D0 mass
+ double kMK = 0.49368; // K+ mass
+ double kMPi = 0.13957; // Pi+ mass
+ double qL=0,qL2=0,ctsD0=0,ctsD0bar=0;
+
+ Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
+
+ Double_t beta = P()/Energy();
+ Double_t gamma = Energy()/kMD0;
+
+ TVector3 mom(momenta[0],momenta[1],momenta[2]);
+ TVector3 mom2(momenta[3],momenta[4],momenta[5]);
+ TVector3 momD(Px(),Py(),Pz());
+
+ qL = mom.Dot(momD)/momD.Mag();
+
+ ctsD0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
+
+ qL2 = mom.Dot(momD)/momD.Mag();
+
+ ctsD0bar = (qL2/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
+
+ if(TMath::Abs(ctsD0) > cutCosThetaStar){
+ if(TMath::Abs(ctsD0bar) > cutCosThetaStar){
+ goodtheta=true;
+ }
+ }
+ return goodtheta;
+}
+bool AliD0Trigger::pTchild()
+{
+ bool goodpT=false;
+ double pT1=0,pT2=0;
+
+ pT1=sqrt(pow(momenta[0],2)+pow(momenta[1],2));
+ pT2=sqrt(pow(momenta[3],2)+pow(momenta[4],2));
+
+ if(pT1>cutpTchild){
+ if(pT2>cutpTchild){
+ goodpT=true;
+ }
+ }
+ return goodpT;
+}
double momenta[6];
double bestV0[3],primaryVertex[3];
- double cutV0low, cutV0high, cutInvMass, cutPointAngle;
+ double cutV0low, cutV0high, cutInvMass, cutPointAngle, cutd0d0,cutCosThetaStar,cutpTchild;
double Bfield;
public:
AliD0Trigger();
- AliD0Trigger(double c[4],double Bfield,double pv[3]);
+ AliD0Trigger(double c[7],double Bfield,double pv[3]);
AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT);
virtual ~AliD0Trigger();
void SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT);
bool FindInvMass();
bool FindV0();
+ bool FindV0offline(double v[3]);
void FindMomentaAtVertex();
+ void FindMomentaOffline();
bool PointingAngle();
-
+ void SetMomenta(double m[6]);
+ bool d0d0();
+ bool CosThetaStar();
+ double P(){return sqrt(Pt()*Pt()+Pz()*Pz());}
+ double Pt(){return sqrt(Px()*Px()+Py()*Py());}
+ double Px(){return (momenta[0]+momenta[3]);}
+ double Py(){return (momenta[1]+momenta[4]);}
+ double Pz(){return (momenta[2]+momenta[5]);}
+ double Energy();
+ bool pTchild();
+
ClassDef(AliD0Trigger,1)
};
Double_t primaryvertex[3] = {0.,0.,0,};
// sigle track cuts
-const Double_t kPtCut = 0.5; // GeV/c
-const Double_t kd0Cut = 50.; // micron
+const Double_t kPtCut = 0.5; // GeV/c
+const Double_t kd0Cut = 50.; // micron
const Double_t kd0CutHigh = 200.; // micron
//cuts for combined tracks
-// cuts[0] = lowest V0 cut
-// cuts[1] = highest V0 cut
-// cuts[2] = inv. mass cut (diiferense)
-// cuts[3] = max value for pointing angle
-const Double_t cuts[4] = {0.005,
- 0.02,
- 0.5,
- 0.95};
-
-// this function applies single track cuts
-Bool_t TrkCuts(const AliITStrackV2& trk);
-
-// this function creates TObjArrays with positive and negative tracks
-void SelectTracks(TTree& itsTree,
- TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
- TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
-
-//void GetPrimaryVertex(int i,Char_t* path="./");
+const Double_t cuts[5] = {0.005, // cuts[0] = lowest V0 cut (cm)
+ 200, // cuts[1] = highest V0 cut (cm) //0.02
+ 0.1, // cuts[2] = inv. mass cut (diferense) (Gev/c)
+ 0.95, // cuts[3] = max cosine value for pointing angle
+ -5000}; // cuts[4] = max cosine value for pointing angle
Int_t iTrkP,iTrkN,itsEntries;
Char_t trksName[100];
Int_t nD0=0;
int ev=0;
-void RunD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
+void RunD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
const Char_t *name="AliD0Trigger";
cerr<<'\n'<<name<<"...\n";
--- /dev/null
+//#define __COMPILE__
+//#ifdef __COMPILE__
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers-------------
+#include <iostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TVector3.h>
+#include <TTree.h>
+#include <TParticle.h>
+#include <TArray.h>
+//----- AliRoot headers ---------------
+#include "alles.h"
+#include "AliRun.h"
+#include "AliKalmanTrack.h"
+#include "AliITStrackV2.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliV0vertex.h"
+#include "AliV0vertexer.h"
+#include "AliITSVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+#include "AliD0Trigger.h"
+#endif
+//-------------------------------------
+// field (T)
+const Double_t kBz = 0.4;
+
+// primary vertex
+Double_t primaryvertex[3] = {0.,0.,0,};
+
+//sec. vertex
+double v2[3]={0,0,0};
+
+// sigle track cuts
+const Double_t kPtCut = 0.5; // GeV/c
+const Double_t kd0Cut = 50.; // micron
+const Double_t kd0CutHigh = 200.; // micron
+
+
+//cuts for combined tracks
+const Double_t cuts[7] = {0.005, // cuts[0] = lowest V0 cut (cm)
+ 0.015, // cuts[1] = highest V0 cut (cm)
+ 0.05, // cuts[2] = inv. mass cut (diferense) (Gev/c)
+ 0.95, // cuts[3] = max cosine value for pointing angle
+ -5000, // cuts[4] = d0d0
+ 0.8, // cuts[5] = costhetastar
+ 0.5}; // cuts[6] = ptchild
+//cut for distance of closest aprach
+double cutDCA=0.01;
+
+// this function applies single track cuts
+Bool_t TrkCuts(const AliITStrackV2& trk);
+
+// this function creates TObjArrays with positive and negative tracks
+void SelectTracks(TTree& itsTree,
+ TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
+ TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
+
+//void GetPrimaryVertex(int i,Char_t* path="./");
+
+Int_t iTrkP,iTrkN,itsEntries;
+Char_t trksName[100];
+Int_t nTrksP=0,nTrksN=0;
+Int_t nD0=0;
+int ev=0;
+double mom[6];
+
+void RunD0offline(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
+
+ const Char_t *name="AliD0offline";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
+
+ // Open file with ITS tracks
+ Char_t fnameTrack[1024];
+ sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
+ TFile* itstrks = TFile::Open(fnameTrack);
+
+ // tracks from ITS
+ sprintf(trksName,"TreeT_ITS_%d",ev);
+ TTree *itsTree=(TTree*)itstrks->Get(trksName);
+ if(!itsTree) continue;
+ itsEntries = (Int_t)itsTree->GetEntries();
+ printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
+
+ //Getting primary Vertex
+ GetPrimaryVertex(0,path);
+
+ // call function which applies sigle track selection and
+ // separetes positives and negatives
+ TObjArray trksP(itsEntries/2);
+ Int_t *itsEntryP = new Int_t[itsEntries];
+ TObjArray trksN(itsEntries/2);
+ Int_t *itsEntryN = new Int_t[itsEntries];
+ SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
+
+ cout<<"#pos: "<<nTrksP<<endl;
+ cout<<"#neg: "<<nTrksN<<endl;
+
+ //the offline stuff
+ // define the cuts for vertexing
+ Double_t vtxcuts[]={33., // max. allowed chi2
+ 0.0, // min. allowed negative daughter's impact param
+ 0.0, // min. allowed positive daughter's impact param
+ 1.0, // max. allowed DCA between the daughter tracks
+ -1.0, // min. allowed cosine of V0's pointing angle
+ 0.0, // min. radius of the fiducial volume
+ 2.9};// max. radius of the fiducial volume
+
+ // create the AliV0vertexer object
+ AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
+
+ AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
+
+ double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
+
+ for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
+ postrack = (AliITStrackV2*)trksP.At(iTrkP);
+ for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
+ negtrack = (AliITStrackV2*)trksN.At(iTrkN);
+ D0.SetTracks(postrack,negtrack);
+ //
+ // ----------- DCA MINIMIZATION ------------------
+ //
+ // find the DCA and propagate the tracks to the DCA
+ double dca = vertexer2->PropagateToDCA(negtrack,postrack);
+
+ if(dca<cutDCA){
+ // define the AliV0vertex object
+ AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
+ // get position of the vertex
+ vertex2->GetXYZ(v2[0],v2[1],v2[2]);
+ delete vertex2;
+ if(D0.pTchild()){
+ if(D0.d0d0()){
+ if(D0.FindV0offline(v2)){
+
+ // momenta of the tracks at the vertex
+ ptP = 1./TMath::Abs(postrack->Get1Pt());
+ alphaP = postrack->GetAlpha();
+ phiP = alphaP+TMath::ASin(postrack->GetSnp());
+ mom[0] = ptP*TMath::Cos(phiP);
+ mom[1] = ptP*TMath::Sin(phiP);
+ mom[2] = ptP*postrack->GetTgl();
+
+ ptN = 1./TMath::Abs(negtrack->Get1Pt());
+ alphaN = negtrack->GetAlpha();
+ phiN = alphaN+TMath::ASin(negtrack->GetSnp());
+ mom[3] = ptN*TMath::Cos(phiN);
+ mom[4] = ptN*TMath::Sin(phiN);
+ mom[5] = ptN*negtrack->GetTgl();
+
+ D0.SetMomenta(mom);
+
+ if(D0.FindInvMass()){
+ if(D0.CosThetaStar()){
+ if(D0.PointingAngle()){
+ nD0++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ cout<<"#D0: "<<nD0<<endl;
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+}
+//___________________________________________________________________________
+void SelectTracks(TTree& itsTree,
+ TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
+ TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
+//
+// this function creates two TObjArrays with positive and negative tracks
+//
+ nTrksP=0,nTrksN=0;
+
+
+ Int_t entr = (Int_t)itsTree.GetEntries();
+
+ // trasfer tracks from tree to arrays
+ for(Int_t i=0; i<entr; i++) {
+
+ AliITStrackV2 *itstrack = new AliITStrackV2;
+ itsTree.SetBranchAddress("tracks",&itstrack);
+
+ itsTree.GetEvent(i);
+
+ // single track selection
+ if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
+
+ if(itstrack->Get1Pt()>0.) { // negative track
+ trksN.AddLast(itstrack);
+ itsEntryN[nTrksN] = i;
+ nTrksN++;
+ } else { // positive track
+ trksP.AddLast(itstrack);
+ itsEntryP[nTrksP] = i;
+ nTrksP++;
+ }
+
+ }
+
+ return;
+}
+//____________________________________________________________________________
+Bool_t TrkCuts(const AliITStrackV2& trk) {
+//
+// this function tells if track passes some kinematical cuts
+//
+ if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
+ if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
+ if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
+
+ return kTRUE;
+}
+//____________________________________________________________________________
+void GetPrimaryVertex(int i,Char_t* path="./") {
+
+ int event=i;
+
+ Char_t falice[1024];
+ sprintf(falice,"%s/galice.root",path);
+ TFile * galice = new TFile(falice);
+
+ TDirectory * curdir;
+
+ Char_t vname[20];
+ galice->cd();
+
+ sprintf(vname,"Vertex_%d",event);
+ TArrayF o = 0;
+ o.Set(3);
+ AliHeader * header = 0;
+
+ TTree * treeE = (TTree*)gDirectory->Get("TE");
+ treeE->SetBranchAddress("Header",&header);
+ treeE->GetEntry(event);
+ AliGenEventHeader* genHeader = header->GenEventHeader();
+ if(genHeader){
+ genHeader->PrimaryVertex(o);
+ primaryvertex[0] = (Double_t)o[0];
+ primaryvertex[1] = (Double_t)o[1];
+ primaryvertex[2] = (Double_t)o[2];
+ }
+ else{
+ printf("Can't find Header");
+ }
+ delete header;
+ delete galice;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff