From eb30eb492a9f182e5976504caf3f12c0d83f595f Mon Sep 17 00:00:00 2001 From: jthaeder Date: Thu, 11 Sep 2008 22:17:43 +0000 Subject: [PATCH] Commit from Sergey: - Fixed bug in ESD converter - Some improvements and bugfixes for the CA Tracker to get better efficiency --- HLT/TPCLib/AliHLTTPCEsdWriterComponent.cxx | 47 +- HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.cxx | 20 +- HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.h | 2 +- HLT/TPCLib/tracking-ca/AliHLTTPCCAGBHit.h | 4 +- HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTrack.h | 9 +- .../tracking-ca/AliHLTTPCCAGBTracker.cxx | 718 +++++++++++------ HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.h | 22 +- HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.cxx | 26 + HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.h | 51 ++ HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.cxx | 33 +- HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.h | 12 +- HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx | 150 +++- HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h | 111 +-- .../tracking-ca/AliHLTTPCCAPerformance.cxx | 737 +++++++++++++++--- .../tracking-ca/AliHLTTPCCAPerformance.h | 63 +- .../tracking-ca/AliHLTTPCCATrackConvertor.cxx | 83 ++ .../tracking-ca/AliHLTTPCCATrackConvertor.h | 34 + .../tracking-ca/AliHLTTPCCATrackParam.cxx | 483 ++++++++++-- .../tracking-ca/AliHLTTPCCATrackParam.h | 24 +- HLT/TPCLib/tracking-ca/AliHLTTPCCATracker.cxx | 208 +++-- .../AliHLTTPCCATrackerComponent.cxx | 268 +++++-- .../tracking-ca/AliHLTTPCCATrackerComponent.h | 13 +- HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx | 248 ++++-- HLT/TPCLib/tracking-ca/AliTPCtrackerCA.h | 3 + HLT/libAliHLTTPC.pkg | 2 + 25 files changed, 2623 insertions(+), 748 deletions(-) create mode 100644 HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.cxx create mode 100644 HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.h create mode 100644 HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.cxx create mode 100644 HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.h diff --git a/HLT/TPCLib/AliHLTTPCEsdWriterComponent.cxx b/HLT/TPCLib/AliHLTTPCEsdWriterComponent.cxx index 026152ff565..e0667052c63 100644 --- a/HLT/TPCLib/AliHLTTPCEsdWriterComponent.cxx +++ b/HLT/TPCLib/AliHLTTPCEsdWriterComponent.cxx @@ -171,11 +171,11 @@ int AliHLTTPCEsdWriterComponent::ProcessBlocks(TTree* pTree, AliESDEvent* pESD, int* pMaxSlice) { // see header file for class documentation + int iResult=0; int iAddedDataBlocks=0; if (pESD && blocks) { - pESD->Reset(); - + pESD->Reset(); pESD->SetMagneticField(fSolenoidBz); const AliHLTComponentBlockData* iter = NULL; AliHLTTPCTrackletData* inPtr=NULL; @@ -197,7 +197,7 @@ int AliHLTTPCEsdWriterComponent::ProcessBlocks(TTree* pTree, AliESDEvent* pESD, if (pMinSlice && (*pMinSlice==-1 || *pMinSlice>minslice)) *pMinSlice=minslice; if (pMaxSlice && (*pMaxSlice==-1 || *pMaxSlicefSpecification, minslice); + HLTDebug("AliHLTTPCEsdWriterComponent::ProcessBlocks: dataspec %#x minslice %d", iter->fSpecification, minslice); if (minslice >=-1 && minslice0 && pTree) { + if ( iAddedDataBlocks>0 && pTree) { pTree->Fill(); } @@ -232,8 +232,7 @@ int AliHLTTPCEsdWriterComponent::Tracks2ESD(AliHLTTPCTrackArray* pTracks, AliESD { // see header file for class documentation int iResult=0; - if (pTracks && pESD) { - HLTDebug("converting %d tracks from track array", pTracks->GetNTracks()); + if (pTracks && pESD) { for (int i=0; iGetNTracks() && iResult>=0; i++) { AliHLTTPCTrack* pTrack=(*pTracks)[i]; if (pTrack) { @@ -241,9 +240,22 @@ int AliHLTTPCEsdWriterComponent::Tracks2ESD(AliHLTTPCTrackArray* pTracks, AliESD //pTrack->Print(); int iLocal=pTrack->Convert2AliKalmanTrack(); if (iLocal>=0) { + AliESDtrack iotrack; iotrack.UpdateTrackParams(pTrack,AliESDtrack::kTPCin); - iotrack.SetTPCPoints(pTrack->GetPoints()); + + Float_t points[4] = {pTrack->GetFirstPointX(), pTrack->GetFirstPointY(), pTrack->GetLastPointX(), pTrack->GetLastPointY() }; + + if(pTrack->GetSector() == -1){ // Set first and last points for global tracks + Double_t s = TMath::Sin( pTrack->GetAlpha() ); + Double_t c = TMath::Cos( pTrack->GetAlpha() ); + points[0] = pTrack->GetFirstPointX()*c + pTrack->GetFirstPointY()*s; + points[1] = -pTrack->GetFirstPointX()*s + pTrack->GetFirstPointY()*c; + points[2] = pTrack->GetLastPointX() *c + pTrack->GetLastPointY() *s; + points[3] = -pTrack->GetLastPointX() *s + pTrack->GetLastPointY() *c; + } + iotrack.SetTPCPoints(points); + //iotrack.SetTPCPoints(pTrack->GetPoints()); pESD->AddTrack(&iotrack); } else { HLTError("conversion to AliKalmanTrack failed for track %d of %d", i, pTracks->GetNTracks()); @@ -336,8 +348,8 @@ AliHLTTPCEsdWriterComponent::AliConverter::AliConverter() // or // refer to README to build package // or - // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt -} + // visit http://web.ift.uib.no/~kjeks + } AliHLTTPCEsdWriterComponent::AliConverter::~AliConverter() { @@ -387,7 +399,12 @@ int AliHLTTPCEsdWriterComponent::AliConverter::DoInit(int argc, const char** arg } else if (argument.CompareTo("-tree")==0) { fWriteTree=1; - } else { + }else if (argument.CompareTo("-solenoidBz")==0) { + TString tmp="-solenoidBz "; + tmp+=argv[++i]; + fBase->Configure(tmp.Data()); + } + else { HLTError("unknown argument %s", argument.Data()); break; } @@ -399,7 +416,7 @@ int AliHLTTPCEsdWriterComponent::AliConverter::DoInit(int argc, const char** arg if (iResult>=0) { iResult=fBase->Reconfigure(NULL, NULL); - } + } return iResult; } @@ -416,7 +433,7 @@ int AliHLTTPCEsdWriterComponent::AliConverter::DoEvent(const AliHLTComponentEven AliHLTUInt8_t* /*outputPtr*/, AliHLTUInt32_t& size, AliHLTComponentBlockDataList& /*outputBlocks*/ ) -{ +{ // see header file for class documentation int iResult=0; // no direct writing to the output buffer @@ -433,20 +450,24 @@ int AliHLTTPCEsdWriterComponent::AliConverter::DoEvent(const AliHLTComponentEven } AliESDEvent* pESD = fESD; + if (pESD && fBase) { + TTree* pTree = NULL; // TODO: Matthias 06.12.2007 // Tried to write the ESD directly instead to a tree, but this did not work // out. Information in the ESD is different, needs investigation. + if (fWriteTree) pTree = new TTree("esdTree", "Tree with HLT ESD objects"); + if (pTree) { pTree->SetDirectory(0); pESD->WriteToTree(pTree); } if ((iResult=fBase->ProcessBlocks(pTree, pESD, blocks, (int)evtData.fBlockCnt))>0) { - // TODO: set the specification correctly + // TODO: set the specification correctly if (pTree) { // the esd structure is written to the user info and is // needed in te ReadFromTree method to read all objects correctly diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.cxx index 5bd8ce4ba07..76b1f2110e1 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.cxx @@ -142,10 +142,13 @@ void AliHLTTPCCADisplay::SetCurrentSlice( AliHLTTPCCATracker *slice ) Double_t dr = .5*(slice->Param().RMax()-slice->Param().RMin()); Double_t cx = 0; Double_t cy = r0; - fYX->Range(cx-dr, cy-dr*1.05, cx+dr, cy+dr); Double_t cz = .5*(slice->Param().ZMax()+slice->Param().ZMin()); Double_t dz = .5*(slice->Param().ZMax()-slice->Param().ZMin())*1.2; - fZX->Range(cz-dz, cy-dr*1.05, cz+dz, cy+dr);//+dr); + fYX->Range(cx-dr, cy-dr*1.05, cx+dr, cy+dr); + fZX->Range(cz-dz, cy-dr*1.05, cz+dz, cy+dr); + + //fYX->Range(cx-dr/6, cy-dr, cx+dr/8, cy-dr + dr/8); + //fZX->Range(cz+dz/2+dz/6, cy-dr , cz+dz-dz/8, cy-dr + dr/8); //fYX->Range(cx-dr/3, cy-dr/3, cx+dr, cy+dr); //fZX->Range(cz-dz/3, cy-dr/3, cz+dz, cy+dr);//+dr); @@ -414,14 +417,14 @@ void AliHLTTPCCADisplay::ConnectCells( Int_t iRow1, AliHLTTPCCACell &cell1, Double_t x12 = row1.X(); Double_t y11 = h11.Y() - h11.ErrY()*3; Double_t y12 = h12.Y() + h12.ErrY()*3; - Double_t z11 = h11.Z(); - Double_t z12 = h12.Z(); + Double_t z11 = h11.Z() - h11.ErrZ()*3; + Double_t z12 = h12.Z() + h12.ErrZ()*3; Double_t x21 = row2.X(); Double_t x22 = row2.X(); Double_t y21 = h21.Y() - h21.ErrY()*3; Double_t y22 = h22.Y() + h22.ErrY()*3; - Double_t z21 = h21.Z(); - Double_t z22 = h22.Z(); + Double_t z21 = h21.Z() - h21.ErrZ()*3; + Double_t z22 = h22.Z() + h22.ErrZ()*3; Double_t vx11, vx12, vy11, vy12, vx21, vx22, vy21, vy22; @@ -451,12 +454,12 @@ void AliHLTTPCCADisplay::ConnectCells( Int_t iRow1, AliHLTTPCCACell &cell1, -void AliHLTTPCCADisplay::DrawTrack1( AliHLTTPCCATrack &track, Int_t color, Bool_t DrawCells ) +void AliHLTTPCCADisplay::DrawTrack( AliHLTTPCCATrack &track, Int_t color, Bool_t DrawCells ) { // draw track if( track.NCells()<2 ) return; - int width = 3; + int width = 1; AliHLTTPCCADisplayTmpCell *vCells = new AliHLTTPCCADisplayTmpCell[track.NCells()]; AliHLTTPCCATrackParam &t = fSlice->ID2Point(track.PointID()[0]).Param(); @@ -586,6 +589,7 @@ void AliHLTTPCCADisplay::DrawTrack1( AliHLTTPCCATrack &track, Int_t color, Bool_ delete[] vCells; } + void AliHLTTPCCADisplay::DrawTrackletPoint( AliHLTTPCCATrackParam &t, Int_t color ) { // draw tracklet point diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.h index e2c9b589a7b..79eb0d41344 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCADisplay.h @@ -68,7 +68,7 @@ class AliHLTTPCCADisplay:public TObject void ConnectCells( Int_t iRow1, AliHLTTPCCACell &cell1, Int_t iRow2, AliHLTTPCCACell &cell2, Int_t color=-1 ); - void DrawTrack1( AliHLTTPCCATrack &track, Int_t color=-1, Bool_t DrawCells=1 ); + void DrawTrack( AliHLTTPCCATrack &track, Int_t color=-1, Bool_t DrawCells=1 ); void DrawTrackletPoint( AliHLTTPCCATrackParam &t, Int_t color=-1 ); void SetSliceTransform( Double_t alpha ); diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBHit.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBHit.h index 43f6586fa3a..aee3f0e0b92 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBHit.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBHit.h @@ -22,7 +22,7 @@ class AliHLTTPCCAGBHit { public: AliHLTTPCCAGBHit() - :fX(0),fY(0),fZ(0),fErrX(0),fErrY(0),fErrZ(0), + :fX(0),fY(0),fZ(0),fErrX(0),fErrY(0),fErrZ(0),fAmp(0), fISlice(0), fIRow(0), fID(0), fIsUsed(0),fSliceHit(){} virtual ~AliHLTTPCCAGBHit(){} @@ -34,6 +34,7 @@ class AliHLTTPCCAGBHit Float_t &ErrX(){ return fErrX; } Float_t &ErrY(){ return fErrY; } Float_t &ErrZ(){ return fErrZ; } + Float_t &Amp() { return fAmp; } Int_t &ISlice(){ return fISlice; } Int_t &IRow(){ return fIRow; } @@ -62,6 +63,7 @@ class AliHLTTPCCAGBHit Float_t fErrY; //* Y position error Float_t fErrZ; //* Z position error + Float_t fAmp; //* Maximal amplitude Int_t fISlice; //* slice number Int_t fIRow; //* row number Int_t fID; //* external ID (id of AliTPCcluster) diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTrack.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTrack.h index aba0a01a9fb..42eddac2561 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTrack.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTrack.h @@ -21,14 +21,14 @@ class AliHLTTPCCAGBTrack { public: - AliHLTTPCCAGBTrack():fFirstHitRef(0),fNHits(0),fParam(),fAlpha(0){ ; } + AliHLTTPCCAGBTrack():fFirstHitRef(0),fNHits(0),fParam(),fAlpha(0),fDeDx(0){ ; } virtual ~AliHLTTPCCAGBTrack(){ ; } Int_t &NHits() { return fNHits; } Int_t &FirstHitRef() { return fFirstHitRef; } AliHLTTPCCATrackParam &Param() { return fParam; } - Double_t &Alpha() { return fAlpha; } - + Float_t &Alpha() { return fAlpha; } + Float_t &DeDx() { return fDeDx; } static Bool_t ComparePNClusters( const AliHLTTPCCAGBTrack *a, const AliHLTTPCCAGBTrack *b){ return (a->fNHits > b->fNHits); } @@ -38,7 +38,8 @@ class AliHLTTPCCAGBTrack Int_t fFirstHitRef; // index of the first hit reference in track->hit reference array Int_t fNHits; // number of track hits AliHLTTPCCATrackParam fParam;// fitted track parameters - Double_t fAlpha; //* Alpha angle of the parametrerisation + Float_t fAlpha; //* Alpha angle of the parametrerisation + Float_t fDeDx; //* DE/DX private: diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.cxx index a371623b8f6..b59d4de5f05 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.cxx @@ -25,7 +25,6 @@ #include "TMath.h" #include "TStopwatch.h" #include "Riostream.h" -#include "TROOT.h" //#define DRAW @@ -46,19 +45,11 @@ AliHLTTPCCAGBTracker::AliHLTTPCCAGBTracker() fTracks(0), fNTracks(0), fSliceTrackInfos(0), + fTime(0), fStatNEvents(0) { //* constructor - fStatTime[0] = 0; - fStatTime[1] = 0; - fStatTime[2] = 0; - fStatTime[3] = 0; - fStatTime[4] = 0; - fStatTime[5] = 0; - fStatTime[6] = 0; - fStatTime[7] = 0; - fStatTime[8] = 0; - fStatTime[9] = 0; + for( Int_t i=0; i<20; i++ ) fStatTime[i] = 0; } AliHLTTPCCAGBTracker::AliHLTTPCCAGBTracker(const AliHLTTPCCAGBTracker&) @@ -71,6 +62,7 @@ AliHLTTPCCAGBTracker::AliHLTTPCCAGBTracker(const AliHLTTPCCAGBTracker&) fTracks(0), fNTracks(0), fSliceTrackInfos(0), + fTime(0), fStatNEvents(0) { //* dummy @@ -139,8 +131,8 @@ void AliHLTTPCCAGBTracker::SetNHits( Int_t nHits ) fNHits = 0; } -void AliHLTTPCCAGBTracker::ReadHit( Double_t x, Double_t y, Double_t z, - Double_t errY, Double_t errZ, +void AliHLTTPCCAGBTracker::ReadHit( Float_t x, Float_t y, Float_t z, + Float_t errY, Float_t errZ, Float_t amp, Int_t ID, Int_t iSlice, Int_t iRow ) { //* read the hit to local array @@ -151,6 +143,7 @@ void AliHLTTPCCAGBTracker::ReadHit( Double_t x, Double_t y, Double_t z, hit.ErrX() = 1.e-4;//fSlices[iSlice].Param().ErrX(); hit.ErrY() = errY; hit.ErrZ() = errZ; + hit.Amp() = amp; hit.ID() = ID; hit.ISlice()=iSlice; hit.IRow() = iRow; @@ -167,7 +160,7 @@ void AliHLTTPCCAGBTracker::ReadHit( Double_t x, Double_t y, Double_t z, void AliHLTTPCCAGBTracker::FindTracks() { //* main tracking routine - + fTime = 0; fStatNEvents++; #ifdef DRAW @@ -222,6 +215,7 @@ void AliHLTTPCCAGBTracker::FindTracks() AliHLTTPCCATracker &slice = fSlices[iSlice]; slice.Reconstruct(); timer.Stop(); + fTime = timer.CpuTime(); fStatTime[0] += timer.CpuTime(); fStatTime[1]+=slice.Timers()[0]; fStatTime[2]+=slice.Timers()[1]; @@ -230,10 +224,7 @@ void AliHLTTPCCAGBTracker::FindTracks() fStatTime[5]+=slice.Timers()[4]; fStatTime[6]+=slice.Timers()[5]; fStatTime[7]+=slice.Timers()[6]; -#ifdef DRAW - //SlicePerformance( iSlice ); - //if( slice.NTracks()>0 ) AliHLTTPCCADisplay::Instance().Ask(); -#endif + fStatTime[8]+=slice.Timers()[7]; } //cout<<"End CA reconstruction"<=0 ){ - jtr = fSliceTrackInfos[jSlice][jtr].fNextNeighbour; - if( jSlice==fNSlices/2-1 ) jSlice = 0; - else if( jSlice==fNSlices-1 ) jSlice = fNSlices/2; - else jSlice = jSlice + 1; + do{ if( fSliceTrackInfos[jSlice][jtr].fUsed ) break; fSliceTrackInfos[jSlice][jtr].fUsed = 1; - AliHLTTPCCAOutTrack &jTr = fSlices[jSlice].OutTracks()[jtr]; + AliHLTTPCCATracker &jslice = fSlices[jSlice]; + AliHLTTPCCAOutTrack &jTr = jslice.OutTracks()[jtr]; for( int jhit=0; jhit=0 ) continue; + p.fISlice = h.ISlice(); + p.fHitID = id; + p.fX = jslice.Rows()[h.IRow()].X(); + p.fY = h.Y(); + p.fZ = h.Z(); + p.fErr2Y = h.ErrY()*h.ErrY(); + p.fErr2Z = h.ErrZ()*h.ErrZ(); + p.fAmp = h.Amp(); + nHits++; } + jtr = fSliceTrackInfos[jSlice][jtr].fNextNeighbour; + jSlice = nextSlice[jSlice]; + } while( jtr >=0 ); + + if( nHits < 30 ) continue; //SG!!! + + + Int_t firstRow = 0, lastRow = maxNRows-1; + for( firstRow=0; firstRow=0 ) break; + } + for( lastRow=maxNRows-1; lastRow>=0; lastRow-- ){ + if( fitPoints[lastRow].fISlice>=0 ) break; } - t.NHits() = ihit; - if( t.NHits()<50 ) continue; - Int_t nHitsOld = t.NHits(); + Int_t mmidRow = (firstRow + lastRow )/2; + Int_t midRow = firstRow; + for( Int_t i=firstRow+1; i<=lastRow; i++ ){ + if( fitPoints[i].fISlice<0 ) continue; + if( TMath::Abs(i-mmidRow)>=TMath::Abs(midRow-mmidRow) ) continue; + midRow = i; + } + if( midRow==firstRow || midRow==lastRow ) continue; + + Int_t searchRows[300]; + Int_t nSearchRows = 0; + for( Int_t i=firstRow; i<=lastRow; i++ ) searchRows[nSearchRows++] = i; + for( Int_t i=lastRow+1; i=0; i-- ) searchRows[nSearchRows++] = i; + // refit - { - if( t.NHits()<10 ) continue; + AliHLTTPCCATrackParam t0; - AliHLTTPCCAGBHit* *vhits = new AliHLTTPCCAGBHit*[nHitsOld]; - - for( Int_t ih=0; ihISlice(); - Int_t currrow = fSlices[currslice].Param().NRows()-1; - Double_t currd2 = 1.e10; - AliHLTTPCCAGBHit *currhit = 0; - - for( Int_t ih=0; ihY(), currhit->Z(), currhit->ErrY(), currhit->ErrZ()); - currhit = 0; - } - if( h.ISlice() != currslice ){ - //* Rotate to the new slice - currhit = 0; - if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[h.ISlice()].Param().Alpha() ) ) break; - currslice = h.ISlice(); - //currrow = 300; - currd2 = 1.e10; - } - //* search for missed hits in rows - { - Double_t factor2 = 3.5*3.5; - AliHLTTPCCATracker &cslice = fSlices[currslice]; - for( Int_t srow=currrow-1; srow>h.IRow(); srow--){ - AliHLTTPCCARow &row = cslice.Rows()[srow]; - if( !t0.TransportToX( row.X() ) ) continue; - Int_t bestsh = -1; - Double_t ds = 1.e10; - for( Int_t ish=0; ishfactor2*s2z ) continue; - Double_t s2y = /*t0.GetErr2Y() + */ sh.ErrY()*sh.ErrY(); - if( dy*dy>factor2*s2y ) continue; - //* update track - t0.Filter(sh.Y(), sh.Z(), sh.ErrY()/.33, sh.ErrZ()/.33); - fTrackHits[nTrackHits+t.NHits()] = sh.ID(); - t.NHits()++; - } - } - //* transport to the new row - currrow = h.IRow(); - Bool_t ok = t0.TransportToX( h.X() ); - if( ok ){ - currrow = h.IRow(); - Double_t dy = h.Y() - t0.GetY(); - Double_t dz = h.Z() - t0.GetZ(); - currd2 = dy*dy + dz*dz; - currhit = &h; - } + FitPoint &p0 = fitPoints[firstRow]; + FitPoint &p1 = fitPoints[midRow]; + FitPoint &p2 = fitPoints[lastRow]; + Float_t x0=p0.fX, y0=p0.fY, z0=p0.fZ; + Float_t x1=p1.fX, y1=p1.fY, z1=p1.fZ; + Float_t x2=p2.fX, y2=p2.fY, z2=p2.fZ; + if( p1.fISlice!=p0.fISlice ){ + Float_t dAlpha = fSlices[p0.fISlice].Param().Alpha() - fSlices[p1.fISlice].Param().Alpha(); + Float_t c = TMath::Cos(dAlpha); + Float_t s = TMath::Sin(dAlpha); + x1 = p1.fX*c + p1.fY*s; + y1 = p1.fY*c - p1.fX*s; } - } - if( currhit != 0 ){ // update track - - t0.Filter(currhit->Y(), currhit->Z(), currhit->ErrY(), currhit->ErrZ()); + if( p2.fISlice!=p0.fISlice ){ + Float_t dAlpha = fSlices[p0.fISlice].Param().Alpha() - fSlices[p2.fISlice].Param().Alpha(); + Float_t c = TMath::Cos(dAlpha); + Float_t s = TMath::Sin(dAlpha); + x2 = p2.fX*c + p2.fY*s; + y2 = p2.fY*c - p2.fX*s; + } + + Float_t sp0[5] = {x0, y0, z0, .5, .5 }; + Float_t sp1[5] = {x1, y1, z1, .5, .5 }; + Float_t sp2[5] = {x2, y2, z2, .5, .5 }; + t0.ConstructXYZ3(sp0,sp1,sp2,1., 0); } - //* search for missed hits in rows - { - Double_t factor2 = 3.5*3.5; - for( Int_t srow=currrow-1; srow>=0; srow--){ + Int_t currslice = fitPoints[firstRow].fISlice; + + for( Int_t rowID=0; rowID=0 ){ + + //* Existing hit + + //* Rotate to the new slice + + if( p.fISlice!=currslice ){ + if( !t0.Rotate( fSlices[p.fISlice].Param().Alpha() - fSlices[currslice].Param().Alpha() ) ) continue; + currslice = p.fISlice; + } + //* Transport to the new row + + if( !t0.TransportToX( p.fX ) ) continue; + + //* Calculate hit errors + + GetErrors2( p.fISlice, iRow, t0, p.fErr2Y, p.fErr2Z ); + + } else { + //continue; //SG!! + //* Search for the missed hit + + Float_t factor2 = 3.5*3.5; + AliHLTTPCCATracker *cslice = &(fSlices[currslice]); - AliHLTTPCCARow *row = &(cslice->Rows()[srow]); + AliHLTTPCCARow *row = &(cslice->Rows()[iRow]); if( !t0.TransportToX( row->X() ) ) continue; + if( t0.GetY() > row->MaxY() ){ //next slice - Int_t j = currslice+1; - if( currslice==fNSlices/2-1 ) j = 0; - else if( currslice==fNSlices-1 ) j = fNSlices/2; + + Int_t j = nextSlice[currslice]; + //* Rotate to the new slice - if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[j].Param().Alpha() ) ) break; + + if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[j].Param().Alpha() ) ) continue; currslice = j; cslice = &(fSlices[currslice]); - row = &(cslice->Rows()[srow]); - if( !t0.TransportToX( row->X() ) ) continue; + row = &(cslice->Rows()[iRow]); + if( !t0.TransportToX( row->X() ) ) continue; + if( TMath::Abs(t0.GetY()) > row->MaxY() ) continue; + }else if( t0.GetY() < -row->MaxY() ){ //prev slice - Int_t j = currslice-1; - if( currslice==0 ) j = fNSlices/2-1; - else if( currslice==fNSlices/2 ) j = fNSlices-1; + Int_t j = prevSlice[currslice]; //* Rotate to the new slice if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[j].Param().Alpha() ) ) break; currslice = j; cslice = &(fSlices[currslice]); - row = &(cslice->Rows()[srow]); + row = &(cslice->Rows()[iRow]); if( !t0.TransportToX( row->X() ) ) continue; + if( TMath::Abs(t0.GetY()) > row->MaxY() ) continue; } + Int_t bestsh = -1; - Double_t ds = 1.e10; + Float_t ds = 1.e10; for( Int_t ish=0; ishNHits(); ish++ ){ AliHLTTPCCAHit &sh = row->Hits()[ish]; - Double_t dy = sh.Y() - t0.GetY(); - Double_t dz = sh.Z() - t0.GetZ(); - Double_t dds = dy*dy+dz*dz; + Float_t dy = sh.Y() - t0.GetY(); + Float_t dz = sh.Z() - t0.GetZ(); + Float_t dds = dy*dy+dz*dz; if( ddsHits()[bestsh]; - Double_t dy = sh.Y() - t0.GetY(); - Double_t dz = sh.Z() - t0.GetZ(); - Double_t s2z = /*t0.GetErr2Z() + */ sh.ErrZ()*sh.ErrZ(); + Float_t dy = sh.Y() - t0.GetY(); + Float_t dz = sh.Z() - t0.GetZ(); + Float_t s2z = /*t0.GetErr2Z() + */ p.fErr2Z; if( dz*dz>factor2*s2z ) continue; - Double_t s2y = /*t0.GetErr2Y() + */ sh.ErrY()*sh.ErrY(); + Float_t s2y = /*t0.GetErr2Y() + */ p.fErr2Y; if( dy*dy>factor2*s2y ) continue; - //* update track - t0.Filter(sh.Y(), sh.Z(), sh.ErrY()/.33, sh.ErrZ()/.33); - fTrackHits[nTrackHits+t.NHits()] = sh.ID(); - t.NHits()++; - } + + p.fISlice = currslice; + p.fHitID = sh.ID(); + p.fX = row->X(); + p.fY = sh.Y(); + p.fZ = sh.Z(); + p.fAmp = fHits[p.fHitID].Amp(); + } + + //* Update the track + + t0.Filter2( p.fY, p.fZ, p.fErr2Y, p.fErr2Z ); } - if( vhits ) delete[] vhits; - //* search for missed hits in rows + //* final refit, dE/dx calculation + //cout<<"\n\nstart refit..\n"<Rows()[srow]); - if( !t0.TransportToX( row->X() ) ) continue; - if( t0.GetY() > row->MaxY() ){ //next slice - Int_t j = currslice+1; - if( currslice==fNSlices/2-1 ) j = 0; - else if( currslice==fNSlices-1 ) j = fNSlices/2; - //* Rotate to the new slice - if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[j].Param().Alpha() ) ) break; - currslice = j; - cslice = &(fSlices[currslice]); - row = &(cslice->Rows()[srow]); - if( !t0.TransportToX( row->X() ) ) continue; - }else if( t0.GetY() < -row->MaxY() ){ //prev slice - Int_t j = currslice-1; - if( currslice==0 ) j = fNSlices/2-1; - else if( currslice==fNSlices/2 ) j = fNSlices-1; - //* Rotate to the new slice - if( !t0.Rotate( -fSlices[currslice].Param().Alpha() +fSlices[j].Param().Alpha() ) ) break; - currslice = j; - cslice = &(fSlices[currslice]); - row = &(cslice->Rows()[srow]); - if( !t0.TransportToX( row->X() ) ) continue; + Double_t sumDeDx = 0; + Int_t nDeDx = 0; + t.NHits() = 0; + + AliHLTTPCCATrackFitParam fitPar; + t0.CalculateFitParameters( fitPar, fSlices[0].Param().Bz() ); + + t0.Cov()[ 0] = .1; + t0.Cov()[ 1] = 0; + t0.Cov()[ 2] = .1; + t0.Cov()[ 3] = 0; + t0.Cov()[ 4] = 0; + t0.Cov()[ 5] = .1; + t0.Cov()[ 6] = 0; + t0.Cov()[ 7] = 0; + t0.Cov()[ 8] = 0; + t0.Cov()[ 9] = .1; + t0.Cov()[10] = 0; + t0.Cov()[11] = 0; + t0.Cov()[12] = 0; + t0.Cov()[13] = 0; + t0.Cov()[14] = .1; + t0.Chi2() = 0; + t0.NDF() = -5; + bool first = 1; + for( Int_t iRow = maxNRows-1; iRow>=0; iRow-- ){ + FitPoint &p = fitPoints[iRow]; + if( p.fISlice<0 ) continue; + fTrackHits[nTrackHits+t.NHits()] = p.fHitID; + t.NHits()++; + + //* Rotate to the new slice + + if( p.fISlice!=currslice ){ + //cout<<"rotate..."<NHits(); ish++ ){ - AliHLTTPCCAHit &sh = row->Hits()[ish]; - Double_t dy = sh.Y() - t0.GetY(); - Double_t dz = sh.Z() - t0.GetZ(); - Double_t dds = dy*dy+dz*dz; - if( ddsHits()[bestsh]; - Double_t dy = sh.Y() - t0.GetY(); - Double_t dz = sh.Z() - t0.GetZ(); - Double_t s2z = /*t0.GetErr2Z() + */ sh.ErrZ()*sh.ErrZ(); - if( dz*dz>factor2*s2z ) continue; - Double_t s2y = /*t0.GetErr2Y() + */ sh.ErrY()*sh.ErrY(); - if( dy*dy>factor2*s2y ) continue; - //* update track - t0.Filter(sh.Y(), sh.Z(), sh.ErrY()/.33, sh.ErrZ()/.33); - fTrackHits[nTrackHits+t.NHits()] = sh.ID(); - t.NHits()++; - } + + //cout<<" before filtration:"<1.e-4 ){ + Float_t dLdX = TMath::Sqrt(1.+t0.DzDs()*t0.DzDs())/TMath::Abs(t0.CosPhi()); + sumDeDx+=p.fAmp/dLdX; + nDeDx++; + } + } + t.DeDx() = 0; + if( nDeDx >0 ) t.DeDx() = sumDeDx/nDeDx; + if( t0.GetErr2Y()<=0 ){ + cout<<"nhits = "<5. || c[2]>5. || c[5]>2. || c[9]>2 || c[14]>2 ) ok = 0; + if(!ok){ + nRejected++; + //cout<<"\n\nRejected: "<126) ? 1:2; + Float_t cosPhiInv = TMath::Abs(t.GetCosPhi())>1.e-2 ?1./t.GetCosPhi() :0; + Float_t angleY = t.GetSinPhi()*cosPhiInv ; + Float_t angleZ = t.GetDzDs()*cosPhiInv ; + + AliHLTTPCCATracker &slice = fSlices[iSlice]; + + Err2Y = slice.Param().GetClusterError2(0,type, z,angleY); + Err2Z = slice.Param().GetClusterError2(1,type, z,angleZ); +} + +void AliHLTTPCCAGBTracker::GetErrors2( AliHLTTPCCAGBHit &h, AliHLTTPCCATrackParam &t, Float_t &Err2Y, Float_t &Err2Z ) +{ + // + // Use calibrated cluster error from OCDB + // + + GetErrors2( h.ISlice(), h.IRow(), t, Err2Y, Err2Z ); +} + +void AliHLTTPCCAGBTracker::WriteSettings( ostream &out ) +{ + out<< NSlices()<> nSlices; + SetNSlices( nSlices ); + for( int iSlice=0; iSlice> nHits; + SetNHits(nHits); + for (Int_t i=0; i>x>>y>>z>>errY>>errZ>>amp>>id>>iSlice>>iRow; + ReadHit( x, y, z, errY, errZ, amp, id, iSlice, iRow ); + } +} + +void AliHLTTPCCAGBTracker::WriteTracks( ostream &out ) +{ + out<>fTime; + fStatTime[0]+=fTime; + fStatNEvents++; + if( fTrackHits ) delete[] fTrackHits; + fTrackHits = 0; + Int_t nTrackHits = 0; + in >> nTrackHits; + fTrackHits = new Int_t [nTrackHits]; + for( Int_t ih=0; ih> TrackHits()[ih]; + } + if( fTracks ) delete[] fTracks; + fTracks = 0; + in >> fNTracks; + fTracks = new AliHLTTPCCAGBTrack[fNTracks]; + for( Int_t itr=0; itr> t.NHits(); + in>> t.FirstHitRef(); + in>> t.Alpha(); + in>> t.DeDx(); + in>> p.X(); + in>> p.CosPhi(); + in>> p.Chi2(); + in>> p.NDF(); + for( Int_t i=0; i<5; i++ ) in>>p.Par()[i]; + for( Int_t i=0; i<15; i++ ) in>>p.Cov()[i]; + } } diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.h index a653fb75951..989f4ce03e6 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAGBTracker.h @@ -15,7 +15,7 @@ class AliHLTTPCCAGBTrack; class AliHLTTPCCAGBHit; class TParticle; class TProfile; - +class AliHLTTPCCATrackParam; /** * @class AliHLTTPCCAGBTracker @@ -46,9 +46,9 @@ public: void SetNSlices( Int_t N ); void SetNHits( Int_t nHits ); - void ReadHit( Double_t x, Double_t y, Double_t z, - Double_t ErrY, Double_t ErrZ, Int_t ID, - Int_t iSlice, Int_t iRow ); + void ReadHit( Float_t x, Float_t y, Float_t z, + Float_t ErrY, Float_t ErrZ, Float_t amp, + Int_t ID, Int_t iSlice, Int_t iRow ); void FindTracks(); void Merging(); @@ -56,11 +56,21 @@ public: AliHLTTPCCAGBHit *Hits(){ return fHits; } Int_t NHits() const { return fNHits; } Int_t NSlices() const { return fNSlices; } + Double_t Time() const { return fTime; } Double_t StatTime( Int_t iTimer ) const { return fStatTime[iTimer]; } Int_t StatNEvents() const { return fStatNEvents; } Int_t NTracks() const { return fNTracks; } AliHLTTPCCAGBTrack *Tracks(){ return fTracks; } Int_t *TrackHits() {return fTrackHits; } + void GetErrors2( AliHLTTPCCAGBHit &h, AliHLTTPCCATrackParam &t, Float_t &Err2Y, Float_t &Err2Z ); + void GetErrors2( Int_t iSlice, Int_t iRow, AliHLTTPCCATrackParam &t, Float_t &Err2Y, Float_t &Err2Z ); + + void WriteSettings( ostream &out ); + void ReadSettings( istream &in ); + void WriteEvent( ostream &out ); + void ReadEvent( istream &in ); + void WriteTracks( ostream &out ); + void ReadTracks( istream &in ); protected: @@ -79,8 +89,8 @@ protected: }; AliHLTTPCCAGBSliceTrackInfo **fSliceTrackInfos; //* additional information for slice tracks - - Double_t fStatTime[10]; //* time measured + Double_t fTime; + Double_t fStatTime[20]; //* timers Int_t fStatNEvents; //* n events proceed ClassDef(AliHLTTPCCAGBTracker,1) //Base class for conformal mapping tracking diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.cxx new file mode 100644 index 00000000000..3f05f1f08dc --- /dev/null +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.cxx @@ -0,0 +1,26 @@ +// $Id: AliHLTTPCCAMCPoint.cxx 27042 2008-07-02 12:06:02Z richterm $ +//*************************************************************************** +// This file is property of and copyright by the ALICE HLT Project * +// ALICE Experiment at CERN, All rights reserved. * +// * +// Primary Authors: Sergey Gorbunov * +// Ivan Kisel * +// for The ALICE HLT Project. * +// * +// Permission to use, copy, modify and distribute this software and its * +// documentation strictly for non-commercial purposes is hereby granted * +// without fee, provided that the above copyright notice appears in all * +// copies and that both the copyright notice and this permission notice * +// appear in the supporting documentation. The authors make no claims * +// about the suitability of this software for any purpose. It is * +// provided "as is" without express or implied warranty. * +//*************************************************************************** + +#include "AliHLTTPCCAMCPoint.h" + + +AliHLTTPCCAMCPoint::AliHLTTPCCAMCPoint() + : fX(0), fY(0), fZ(0), fSx(0), fSy(0), fSz(0), fTime(0), fISlice(0), fTrackID(0) +{ + //* Default constructor +} diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.h new file mode 100644 index 00000000000..eb9dff2cd7e --- /dev/null +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCPoint.h @@ -0,0 +1,51 @@ +//-*- Mode: C++ -*- + +//* This file is property of and copyright by the ALICE HLT Project * +//* ALICE Experiment at CERN, All rights reserved. * +//* See cxx source for full Copyright notice * + +#ifndef ALIHLTTPCCAMCPOINT_H +#define ALIHLTTPCCAMCPOINT_H + +#include "Rtypes.h" + + +/** + * @class AliHLTTPCCAMCPoint + * store MC point information for AliHLTTPCCAPerformance + */ +class AliHLTTPCCAMCPoint +{ + public: + + AliHLTTPCCAMCPoint(); + + Float_t &X() { return fX; } + Float_t &Y() { return fY; } + Float_t &Z() { return fZ; } + Float_t &Sx() { return fSx; } + Float_t &Sy() { return fSy; } + Float_t &Sz() { return fSz; } + Float_t &Time() { return fTime; } + Int_t &ISlice() { return fISlice; } + Int_t &TrackID() { return fTrackID; } + + static Bool_t Compare( const AliHLTTPCCAMCPoint &p1, const AliHLTTPCCAMCPoint &p2 ) + { + return (p1.fTrackID < p2.fTrackID); + } + + protected: + + Float_t fX; //* global X position + Float_t fY; //* global Y position + Float_t fZ; //* global Z position + Float_t fSx; //* slice X position + Float_t fSy; //* slice Y position + Float_t fSz; //* slice Z position + Float_t fTime; //* time + Int_t fISlice; //* slice number + Int_t fTrackID; //* mc track number +}; + +#endif diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.cxx index 311330aacd8..45e099a824f 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.cxx @@ -23,18 +23,19 @@ AliHLTTPCCAMCTrack::AliHLTTPCCAMCTrack() - : fP(0), fPt(0), fNHits(0), fNReconstructed(0), fSet(0), fNTurns(0) + : fPDG(0), fP(0), fPt(0), fNHits(0), fNMCPoints(0), fFirstMCPointID(0), fNReconstructed(0), fSet(0), fNTurns(0) { //* Default constructor } AliHLTTPCCAMCTrack::AliHLTTPCCAMCTrack( const TParticle *part ) - : fP(0), fPt(0), fNHits(0), fNReconstructed(0), fSet(0), fNTurns(0) + : fPDG(0), fP(0), fPt(0), fNHits(0), fNMCPoints(0), fFirstMCPointID(0), fNReconstructed(0), fSet(0), fNTurns(0) { //* Constructor from TParticle for( Int_t i=0; i<7; i++ ) fPar[i] = 0; + for( Int_t i=0; i<7; i++ ) fTPCPar[i] = 0; fP = 0; fPt = 0; @@ -52,9 +53,31 @@ AliHLTTPCCAMCTrack::AliHLTTPCCAMCTrack( const TParticle *part ) fPar[4] = part->Py()*pi; fPar[5] = part->Pz()*pi; fPar[6] = 0; - Int_t pdg = part->GetPdgCode(); - if ( TMath::Abs(pdg) < 100000 ){ - TParticlePDG *pPDG = TDatabasePDG::Instance()->GetParticle(pdg); + fPDG = part->GetPdgCode(); + if ( TMath::Abs(fPDG) < 100000 ){ + TParticlePDG *pPDG = TDatabasePDG::Instance()->GetParticle(fPDG); if( pPDG ) fPar[6] = pPDG->Charge()/3.0*pi; } } + +void AliHLTTPCCAMCTrack::SetTPCPar( Float_t X, Float_t Y, Float_t Z, + Float_t Px, Float_t Py, Float_t Pz ) +{ + //* Set parameters at TPC entrance + + for( Int_t i=0; i<7; i++ ) fTPCPar[i] = 0; + + fTPCPar[0] = X; + fTPCPar[1] = Y; + fTPCPar[2] = Z; + Double_t p = TMath::Sqrt(Px*Px + Py*Py + Pz*Pz); + Double_t pi = ( p >1.e-4 ) ?1./p :0; + fTPCPar[3] = Px*pi; + fTPCPar[4] = Py*pi; + fTPCPar[5] = Pz*pi; + fTPCPar[6] = 0; + if ( TMath::Abs(fPDG) < 100000 ){ + TParticlePDG *pPDG = TDatabasePDG::Instance()->GetParticle(fPDG); + if( pPDG ) fTPCPar[6] = pPDG->Charge()/3.0*pi; + } +} diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.h index 20398f36490..cd3d793bd4e 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAMCTrack.h @@ -23,24 +23,34 @@ class AliHLTTPCCAMCTrack AliHLTTPCCAMCTrack(); AliHLTTPCCAMCTrack( const TParticle *part ); + void SetTPCPar( Float_t X, Float_t Y, Float_t Z, Float_t Px, Float_t Py, Float_t Pz ); + Int_t &PDG() { return fPDG;} Double_t *Par() { return fPar; } + Double_t *TPCPar() { return fTPCPar; } Double_t &P() { return fP; } Double_t &Pt() { return fPt; } + Int_t &NHits() { return fNHits;} + Int_t &NMCPoints() { return fNMCPoints;} + Int_t &FirstMCPointID(){ return fFirstMCPointID;} Int_t &NReconstructed(){ return fNReconstructed; } Int_t &Set() { return fSet; } Int_t &NTurns() { return fNTurns; } protected: + Int_t fPDG; //* particle pdg code Double_t fPar[7]; //* x,y,z,ex,ey,ez,q/p + Double_t fTPCPar[7]; //* x,y,z,ex,ey,ez,q/p at TPC entrance (x=y=0 means no information) Double_t fP, fPt; //* momentum and transverse momentum Int_t fNHits; //* N TPC clusters + Int_t fNMCPoints; //* N MC points + Int_t fFirstMCPointID; //* id of the first MC point in the points array Int_t fNReconstructed; //* how many times is reconstructed Int_t fSet; //* set of tracks 0-OutSet, 1-ExtraSet, 2-RefSet Int_t fNTurns; //* N of turns in the current sector - + }; #endif diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx index e7c1c0218e3..39c59ad317c 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx @@ -27,20 +27,63 @@ AliHLTTPCCAParam::AliHLTTPCCAParam() fCosAlpha(0), fSinAlpha(0), fAngleMin(0), fAngleMax(0), fRMin(83.65), fRMax(133.3), fZMin(0.0529937), fZMax(249.778), fErrX(0), fErrY(0), fErrZ(0.228808),fPadPitch(0.4),fBz(-5.), fYErrorCorrection(0.33), fZErrorCorrection(0.45), - fCellConnectionAngleXY(35./180.*TMath::Pi()), - fCellConnectionAngleXZ(35./180.*TMath::Pi()), + fCellConnectionAngleXY(45./180.*TMath::Pi()), + fCellConnectionAngleXZ(45./180.*TMath::Pi()), fMaxTrackMatchDRow(4), fTrackConnectionFactor(3.5), fTrackChiCut(3.5), fTrackChi2Cut(10) { + fParamS0Par[0][0][0] = 0.00047013; + fParamS0Par[0][0][1] = 2.00135e-05; + fParamS0Par[0][0][2] = 0.0106533; + fParamS0Par[0][0][3] = 5.27104e-08; + fParamS0Par[0][0][4] = 0.012829; + fParamS0Par[0][0][5] = 0.000147125; + fParamS0Par[0][0][6] = 4.99432; + fParamS0Par[0][1][0] = 0.000883342; + fParamS0Par[0][1][1] = 1.07011e-05; + fParamS0Par[0][1][2] = 0.0103187; + fParamS0Par[0][1][3] = 4.25141e-08; + fParamS0Par[0][1][4] = 0.0224292; + fParamS0Par[0][1][5] = 8.27274e-05; + fParamS0Par[0][1][6] = 4.17233; + fParamS0Par[0][2][0] = 0.000745399; + fParamS0Par[0][2][1] = 5.62408e-06; + fParamS0Par[0][2][2] = 0.0151562; + fParamS0Par[0][2][3] = 5.08757e-08; + fParamS0Par[0][2][4] = 0.0601004; + fParamS0Par[0][2][5] = 7.97129e-05; + fParamS0Par[0][2][6] = 4.84913; + fParamS0Par[1][0][0] = 0.00215126; + fParamS0Par[1][0][1] = 6.82233e-05; + fParamS0Par[1][0][2] = 0.0221867; + fParamS0Par[1][0][3] = -6.27825e-09; + fParamS0Par[1][0][4] = -0.00745378; + fParamS0Par[1][0][5] = 0.000172629; + fParamS0Par[1][0][6] = 6.24987; + fParamS0Par[1][1][0] = 0.00181667; + fParamS0Par[1][1][1] = -4.17772e-06; + fParamS0Par[1][1][2] = 0.0253429; + fParamS0Par[1][1][3] = 1.3011e-07; + fParamS0Par[1][1][4] = -0.00362827; + fParamS0Par[1][1][5] = 0.00030406; + fParamS0Par[1][1][6] = 17.7775; + fParamS0Par[1][2][0] = 0.00158251; + fParamS0Par[1][2][1] = -3.55911e-06; + fParamS0Par[1][2][2] = 0.0247899; + fParamS0Par[1][2][3] = 7.20604e-08; + fParamS0Par[1][2][4] = 0.0179946; + fParamS0Par[1][2][5] = 0.000425504; + fParamS0Par[1][2][6] = 20.9294; + Update(); } void AliHLTTPCCAParam::Initialize( Int_t iSlice, - Int_t nRows, Double_t rowX[], - Double_t alpha, Double_t dAlpha, - Double_t rMin, Double_t rMax, - Double_t zMin, Double_t zMax, - Double_t padPitch, Double_t zSigma, - Double_t bz + Int_t nRows, Float_t rowX[], + Float_t alpha, Float_t dAlpha, + Float_t rMin, Float_t rMax, + Float_t zMin, Float_t zMax, + Float_t padPitch, Float_t zSigma, + Float_t bz ) { // initialization @@ -74,8 +117,8 @@ void AliHLTTPCCAParam::Update() fTrackChi2Cut = fTrackChiCut * fTrackChiCut; } -void AliHLTTPCCAParam::Slice2Global( Double_t x, Double_t y, Double_t z, - Double_t *X, Double_t *Y, Double_t *Z ) const +void AliHLTTPCCAParam::Slice2Global( Float_t x, Float_t y, Float_t z, + Float_t *X, Float_t *Y, Float_t *Z ) const { // conversion of coorinates sector->global *X = x*fCosAlpha - y*fSinAlpha; @@ -83,11 +126,94 @@ void AliHLTTPCCAParam::Slice2Global( Double_t x, Double_t y, Double_t z, *Z = z; } -void AliHLTTPCCAParam::Global2Slice( Double_t X, Double_t Y, Double_t Z, - Double_t *x, Double_t *y, Double_t *z ) const +void AliHLTTPCCAParam::Global2Slice( Float_t X, Float_t Y, Float_t Z, + Float_t *x, Float_t *y, Float_t *z ) const { // conversion of coorinates global->sector *x = X*fCosAlpha + Y*fSinAlpha; *y = Y*fCosAlpha - X*fSinAlpha; *z = Z; } + +Float_t AliHLTTPCCAParam::GetClusterError2( Int_t yz, Int_t type, Float_t z, Float_t angle ) +{ + //* recalculate the cluster error wih respect to the track slope + Float_t angle2 = angle*angle; + Float_t *c = fParamS0Par[yz][type]; + Float_t v = c[0] + z*(c[1] + c[3]*z) + angle2*(c[2] + angle2*c[4] + c[5]*z ); + return TMath::Abs(v); +} + +void AliHLTTPCCAParam::WriteSettings( ostream &out ) +{ + out << fISlice<> fISlice; + in >> fNRows; + in >> fAlpha; + in >> fDAlpha; + in >> fCosAlpha; + in >> fSinAlpha; + in >> fAngleMin; + in >> fAngleMax; + in >> fRMin; + in >> fRMax; + in >> fZMin; + in >> fZMax; + in >> fErrX; + in >> fErrY; + in >> fErrZ; + in >> fPadPitch; + in >> fBz; + in >> fYErrorCorrection; + in >> fZErrorCorrection; + in >> fCellConnectionAngleXY; + in >> fCellConnectionAngleXZ; + in >> fMaxTrackMatchDRow; + in >> fTrackConnectionFactor; + in >> fTrackChiCut; + in >> fTrackChi2Cut; + for( Int_t iRow = 0; iRow> fRowX[iRow]; + } + for( Int_t i=0; i<2; i++ ) + for( Int_t j=0; j<3; j++ ) + for( Int_t k=0; k<7; k++ ) + in >> fParamS0Par[i][j][k]; +} diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h index c2e92a0c9cf..a175a5eec98 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h @@ -9,6 +9,7 @@ #define ALIHLTTPCCAPARAM_H #include "Rtypes.h" +#include "Riostream.h" /** * @class ALIHLTTPCCAParam @@ -25,70 +26,76 @@ class AliHLTTPCCAParam AliHLTTPCCAParam(); virtual ~AliHLTTPCCAParam(){;} - void Initialize( Int_t iSlice, Int_t nRows, Double_t rowX[], - Double_t alpha, Double_t dAlpha, - Double_t rMin, Double_t rMax, Double_t zMin, Double_t zMax, - Double_t padPitch, Double_t zSigma, Double_t bz ); + void Initialize( Int_t iSlice, Int_t nRows, Float_t rowX[], + Float_t alpha, Float_t dAlpha, + Float_t rMin, Float_t rMax, Float_t zMin, Float_t zMax, + Float_t padPitch, Float_t zSigma, Float_t bz ); void Update(); - void Slice2Global( Double_t x, Double_t y, Double_t z, - Double_t *X, Double_t *Y, Double_t *Z ) const; - void Global2Slice( Double_t x, Double_t y, Double_t z, - Double_t *X, Double_t *Y, Double_t *Z ) const; + void Slice2Global( Float_t x, Float_t y, Float_t z, + Float_t *X, Float_t *Y, Float_t *Z ) const; + void Global2Slice( Float_t x, Float_t y, Float_t z, + Float_t *X, Float_t *Y, Float_t *Z ) const; Int_t &ISlice(){ return fISlice;} Int_t &NRows(){ return fNRows;} - Double_t &RowX( Int_t iRow ){ return fRowX[iRow]; } + Float_t &RowX( Int_t iRow ){ return fRowX[iRow]; } - Double_t &Alpha(){ return fAlpha;} - Double_t &DAlpha(){ return fDAlpha;} - Double_t &CosAlpha(){ return fCosAlpha;} - Double_t &SinAlpha(){ return fSinAlpha;} - Double_t &AngleMin(){ return fAngleMin;} - Double_t &AngleMax(){ return fAngleMax;} - Double_t &RMin(){ return fRMin;} - Double_t &RMax(){ return fRMax;} - Double_t &ZMin(){ return fZMin;} - Double_t &ZMax(){ return fZMax;} - Double_t &ErrZ(){ return fErrZ;} - Double_t &ErrX(){ return fErrX;} - Double_t &ErrY(){ return fErrY;} - Double_t &Bz(){ return fBz;} - - Double_t &TrackConnectionFactor(){ return fTrackConnectionFactor; } - Double_t &TrackChiCut() { return fTrackChiCut; } - Double_t &TrackChi2Cut(){ return fTrackChi2Cut; } - Int_t &MaxTrackMatchDRow(){ return fMaxTrackMatchDRow; } - Double_t &YErrorCorrection(){ return fYErrorCorrection; } - Double_t &ZErrorCorrection(){ return fZErrorCorrection; } - Double_t &CellConnectionAngleXY(){ return fCellConnectionAngleXY; } - Double_t &CellConnectionAngleXZ(){ return fCellConnectionAngleXZ; } + Float_t &Alpha(){ return fAlpha;} + Float_t &DAlpha(){ return fDAlpha;} + Float_t &CosAlpha(){ return fCosAlpha;} + Float_t &SinAlpha(){ return fSinAlpha;} + Float_t &AngleMin(){ return fAngleMin;} + Float_t &AngleMax(){ return fAngleMax;} + Float_t &RMin(){ return fRMin;} + Float_t &RMax(){ return fRMax;} + Float_t &ZMin(){ return fZMin;} + Float_t &ZMax(){ return fZMax;} + Float_t &ErrZ(){ return fErrZ;} + Float_t &ErrX(){ return fErrX;} + Float_t &ErrY(){ return fErrY;} + Float_t &Bz(){ return fBz;} + + Float_t &TrackConnectionFactor(){ return fTrackConnectionFactor; } + Float_t &TrackChiCut() { return fTrackChiCut; } + Float_t &TrackChi2Cut(){ return fTrackChi2Cut; } + Int_t &MaxTrackMatchDRow(){ return fMaxTrackMatchDRow; } + Float_t &YErrorCorrection(){ return fYErrorCorrection; } + Float_t &ZErrorCorrection(){ return fZErrorCorrection; } + Float_t &CellConnectionAngleXY(){ return fCellConnectionAngleXY; } + Float_t &CellConnectionAngleXZ(){ return fCellConnectionAngleXZ; } + + Float_t GetClusterError2(Int_t yz, Int_t type, Float_t z, Float_t angle ); + + void WriteSettings( ostream &out ); + void ReadSettings( istream &in ); protected: Int_t fISlice; // slice number Int_t fNRows; // number of rows - Double_t fAlpha, fDAlpha; // slice angle and angular size - Double_t fCosAlpha, fSinAlpha;// sign and cosine of the slice angle - Double_t fAngleMin, fAngleMax; // minimal and maximal angle - Double_t fRMin, fRMax;// slice R range - Double_t fZMin, fZMax;// slice Z range - Double_t fErrX, fErrY, fErrZ;// default cluster errors - Double_t fPadPitch; // pad pitch - Double_t fBz; // magnetic field value (only constant field can be used) - - Double_t fYErrorCorrection;// correction factor for Y error of input clusters - Double_t fZErrorCorrection;// correction factor for Z error of input clusters - - Double_t fCellConnectionAngleXY; // max phi angle between connected cells - Double_t fCellConnectionAngleXZ; // max psi angle between connected cells - Int_t fMaxTrackMatchDRow;// maximal jump in TPC row for connecting track segments - Double_t fTrackConnectionFactor; // allowed distance in Chi^2/3.5 for neighbouring tracks - Double_t fTrackChiCut; // cut for track Sqrt(Chi2/NDF); - Double_t fTrackChi2Cut;// cut for track Chi^2/NDF - - Double_t fRowX[200];// X-coordinate of rows + Float_t fAlpha, fDAlpha; // slice angle and angular size + Float_t fCosAlpha, fSinAlpha;// sign and cosine of the slice angle + Float_t fAngleMin, fAngleMax; // minimal and maximal angle + Float_t fRMin, fRMax;// slice R range + Float_t fZMin, fZMax;// slice Z range + Float_t fErrX, fErrY, fErrZ;// default cluster errors + Float_t fPadPitch; // pad pitch + Float_t fBz; // magnetic field value (only constant field can be used) + + Float_t fYErrorCorrection;// correction factor for Y error of input clusters + Float_t fZErrorCorrection;// correction factor for Z error of input clusters + + Float_t fCellConnectionAngleXY; // max phi angle between connected cells + Float_t fCellConnectionAngleXZ; // max psi angle between connected cells + Int_t fMaxTrackMatchDRow;// maximal jump in TPC row for connecting track segments + Float_t fTrackConnectionFactor; // allowed distance in Chi^2/3.5 for neighbouring tracks + Float_t fTrackChiCut; // cut for track Sqrt(Chi2/NDF); + Float_t fTrackChi2Cut;// cut for track Chi^2/NDF + + Float_t fRowX[200];// X-coordinate of rows + Float_t fParamS0Par[2][3][7]; // cluster error parameterization coeficients ClassDef(AliHLTTPCCAParam,1); }; diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx index b2d66ef9d4b..1b79e3928fc 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx @@ -19,15 +19,19 @@ #include "AliHLTTPCCAPerformance.h" #include "AliHLTTPCCAGBHit.h" #include "AliHLTTPCCAMCTrack.h" +#include "AliHLTTPCCAMCPoint.h" #include "AliHLTTPCCAOutTrack.h" +#include "AliHLTTPCCAGBTrack.h" #include "AliHLTTPCCAGBTracker.h" #include "AliHLTTPCCATracker.h" + #include "TMath.h" #include "TROOT.h" #include "Riostream.h" #include "TFile.h" #include "TH1.h" +#include "TH2.h" #include "TProfile.h" @@ -42,6 +46,9 @@ AliHLTTPCCAPerformance::AliHLTTPCCAPerformance() fNHits(0), fMCTracks(0), fNMCTracks(0), + fMCPoints(0), + fNMCPoints(0), + fDoClusterPulls(1), fStatNEvents(0), fStatNRecTot(0), fStatNRecOut(0), @@ -52,20 +59,42 @@ AliHLTTPCCAPerformance::AliHLTTPCCAPerformance() fStatNMCRef(0), fStatNRecRef(0), fStatNClonesRef(0), + fStatGBNRecTot(0), + fStatGBNRecOut(0), + fStatGBNGhost(0), + fStatGBNMCAll(0), + fStatGBNRecAll(0), + fStatGBNClonesAll(0), + fStatGBNMCRef(0), + fStatGBNRecRef(0), + fStatGBNClonesRef(0), fHistoDir(0), + fhResY(0), + fhResZ(0), + fhResSinPhi(0), + fhResDzDs(0), + fhResPt(0), + fhPullY(0), + fhPullZ(0), + fhPullSinPhi(0), + fhPullDzDs(0), + fhPullQPt(0), fhHitErrY(0), fhHitErrZ(0), - fhHitResX(0), fhHitResY(0), fhHitResZ(0), - fhHitPullX(0), fhHitPullY(0), fhHitPullZ(0), + fhHitResY1(0), + fhHitResZ1(0), + fhHitPullY1(0), + fhHitPullZ1(0), fhCellPurity(0), fhCellNHits(0), fhCellPurityVsN(0), fhCellPurityVsPt(0), - fhEffVsP(0) + fhEffVsP(0), + fhGBEffVsP(0) { //* constructor } @@ -78,6 +107,9 @@ AliHLTTPCCAPerformance::AliHLTTPCCAPerformance(const AliHLTTPCCAPerformance&) fNHits(0), fMCTracks(0), fNMCTracks(0), + fMCPoints(0), + fNMCPoints(0), + fDoClusterPulls(1), fStatNEvents(0), fStatNRecTot(0), fStatNRecOut(0), @@ -88,20 +120,42 @@ AliHLTTPCCAPerformance::AliHLTTPCCAPerformance(const AliHLTTPCCAPerformance&) fStatNMCRef(0), fStatNRecRef(0), fStatNClonesRef(0), + fStatGBNRecTot(0), + fStatGBNRecOut(0), + fStatGBNGhost(0), + fStatGBNMCAll(0), + fStatGBNRecAll(0), + fStatGBNClonesAll(0), + fStatGBNMCRef(0), + fStatGBNRecRef(0), + fStatGBNClonesRef(0), fHistoDir(0), + fhResY(0), + fhResZ(0), + fhResSinPhi(0), + fhResDzDs(0), + fhResPt(0), + fhPullY(0), + fhPullZ(0), + fhPullSinPhi(0), + fhPullDzDs(0), + fhPullQPt(0), fhHitErrY(0), - fhHitErrZ(0), - fhHitResX(0), + fhHitErrZ(0), fhHitResY(0), fhHitResZ(0), - fhHitPullX(0), fhHitPullY(0), fhHitPullZ(0), + fhHitResY1(0), + fhHitResZ1(0), + fhHitPullY1(0), + fhHitPullZ1(0), fhCellPurity(0), fhCellNHits(0), fhCellPurityVsN(0), fhCellPurityVsPt(0), - fhEffVsP(0) + fhEffVsP(0), + fhGBEffVsP(0) { //* dummy } @@ -134,6 +188,9 @@ void AliHLTTPCCAPerformance::StartEvent() if( fMCTracks ) delete[] fMCTracks; fMCTracks = 0; fNMCTracks = 0; + if( fMCPoints ) delete[] fMCPoints; + fMCPoints = 0; + fNMCPoints = 0; } void AliHLTTPCCAPerformance::SetNHits( Int_t NHits ) @@ -154,6 +211,15 @@ void AliHLTTPCCAPerformance::SetNMCTracks( Int_t NMCTracks ) fNMCTracks = NMCTracks; } +void AliHLTTPCCAPerformance::SetNMCPoints( Int_t NMCPoints ) +{ + //* set number of MC points + if( fMCPoints ) delete[] fMCPoints; + fMCPoints = 0; + fMCPoints = new AliHLTTPCCAMCPoint[ NMCPoints ]; + fNMCPoints = 0; +} + void AliHLTTPCCAPerformance::ReadHitLabel( Int_t HitID, Int_t lab0, Int_t lab1, Int_t lab2 ) { @@ -167,49 +233,81 @@ void AliHLTTPCCAPerformance::ReadHitLabel( Int_t HitID, void AliHLTTPCCAPerformance::ReadMCTrack( Int_t index, const TParticle *part ) { - //* read mc track to local array + //* read mc track to the local array fMCTracks[index] = AliHLTTPCCAMCTrack(part); } +void AliHLTTPCCAPerformance::ReadMCTPCTrack( Int_t index, Float_t X, Float_t Y, Float_t Z, + Float_t Px, Float_t Py, Float_t Pz ) +{ + //* read mc track parameters at TPC + fMCTracks[index].SetTPCPar(X,Y,Z,Px,Py,Pz); +} + +void AliHLTTPCCAPerformance::ReadMCPoint( Int_t TrackID, Float_t X, Float_t Y, Float_t Z, Float_t Time, Int_t iSlice ) +{ + //* read mc point to the local array + AliHLTTPCCAMCPoint &p = fMCPoints[fNMCPoints]; + p.TrackID() = TrackID; + p.X() = X; + p.Y() = Y; + p.Z() = Z; + p.Time() = Time; + p.ISlice() = iSlice; + fTracker->Slices()[iSlice].Param().Global2Slice( X, Y, Z, + &p.Sx(), &p.Sy(), &p.Sz() ); + if( X*X + Y*Y>10.) fNMCPoints++; +} + void AliHLTTPCCAPerformance::CreateHistos() { //* create performance histogramms TDirectory *curdir = gDirectory; fHistoDir = gROOT->mkdir("HLTTPCCATrackerPerformance"); fHistoDir->cd(); - gDirectory->mkdir("Fit"); - gDirectory->cd("Fit"); - /* - fhResY = new TH1D("resY", "Y resoltion [cm]", 100, -5., 5.); - fhResZ = new TH1D("resZ", "Z resoltion [cm]", 100, -5., 5.); - fhResTy = new TH1D("resTy", "Ty resoltion ", 100, -1., 1.); - fhResTz = new TH1D("resTz", "Tz resoltion ", 100, -1., 1.); - fhResP = new TH1D("resP", "P resoltion ", 100, -10., 10.); - fhPullY = new TH1D("pullY", "Y pull", 100, -10., 10.); - fhPullZ = new TH1D("pullZ", "Z pull", 100, -10., 10.); - fhPullTy = new TH1D("pullTy", "Ty pull", 100, -10., 10.); - fhPullTz = new TH1D("pullTz", "Tz pull", 100, -10., 10.); - fhPullQp = new TH1D("pullQp", "Qp pull", 100, -10., 10.); - */ + gDirectory->mkdir("TrackFit"); + gDirectory->cd("TrackFit"); + + fhResY = new TH1D("resY", "track Y resoltion [cm]", 30, -.5, .5); + fhResZ = new TH1D("resZ", "track Z resoltion [cm]", 30, -.5, .5); + fhResSinPhi = new TH1D("resSinPhi", "track SinPhi resoltion ", 30, -.03, .03); + fhResDzDs = new TH1D("resDzDs", "track DzDs resoltion ", 30, -.01, .01); + fhResPt = new TH1D("resPt", "track telative Pt resoltion", 30, -.2, .2); + fhPullY = new TH1D("pullY", "track Y pull", 30, -10., 10.); + fhPullZ = new TH1D("pullZ", "track Z pull", 30, -10., 10.); + fhPullSinPhi = new TH1D("pullSinPhi", "track SinPhi pull", 30, -10., 10.); + fhPullDzDs = new TH1D("pullDzDs", "track DzDs pull", 30, -10., 10.); + fhPullQPt = new TH1D("pullQPt", "track Q/Pt pull", 30, -10., 10.); + gDirectory->cd(".."); fhEffVsP = new TProfile("EffVsP", "Eff vs P", 100, 0., 5.); + fhGBEffVsP = new TProfile("GBEffVsP", "Global tracker: Eff vs P", 100, 0., 5.); - fhHitResX = new TH1D("resHitX", "X cluster resoltion [cm]", 100, -2., 2.); + gDirectory->mkdir("Clusters"); + gDirectory->cd("Clusters"); fhHitResY = new TH1D("resHitY", "Y cluster resoltion [cm]", 100, -2., 2.); fhHitResZ = new TH1D("resHitZ", "Z cluster resoltion [cm]", 100, -2., 2.); - fhHitPullX = new TH1D("pullHitX", "X cluster pull", 100, -10., 10.); - fhHitPullY = new TH1D("pullHitY", "Y cluster pull", 100, -10., 10.); - fhHitPullZ = new TH1D("pullHitZ", "Z cluster pull", 100, -10., 10.); + fhHitPullY = new TH1D("pullHitY", "Y cluster pull", 50, -10., 10.); + fhHitPullZ = new TH1D("pullHitZ", "Z cluster pull", 50, -10., 10.); + fhHitResY1 = new TH1D("resHitY1", "Y cluster resoltion [cm]", 100, -2., 2.); + fhHitResZ1 = new TH1D("resHitZ1", "Z cluster resoltion [cm]", 100, -2., 2.); + fhHitPullY1 = new TH1D("pullHitY1", "Y cluster pull", 50, -10., 10.); + fhHitPullZ1 = new TH1D("pullHitZ1", "Z cluster pull", 50, -10., 10.); fhHitErrY = new TH1D("HitErrY", "Y cluster error [cm]", 100, 0., 1.); fhHitErrZ = new TH1D("HitErrZ", "Z cluster error [cm]", 100, 0., 1.); + gDirectory->cd(".."); + + gDirectory->mkdir("Cells"); + gDirectory->cd("Cells"); fhCellPurity = new TH1D("CellPurity", "Cell Purity", 100, -0.1, 1.1); fhCellNHits = new TH1D("CellNHits", "Cell NHits", 40, 0., 40.); fhCellPurityVsN = new TProfile("CellPurityVsN", "Cell purity Vs N hits", 40, 2., 42.); fhCellPurityVsPt = new TProfile("CellPurityVsPt", "Cell purity Vs Pt", 100, 0., 5.); + gDirectory->cd(".."); curdir->cd(); } @@ -246,17 +344,17 @@ void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag ) { //* calculate slice tracker performance if( !fTracker ) return; - - int nRecTot = 0, nGhost=0, nRecOut=0; - int nMCAll = 0, nRecAll=0, nClonesAll=0; - int nMCRef = 0, nRecRef=0, nClonesRef=0; + + Int_t nRecTot = 0, nGhost=0, nRecOut=0; + Int_t nMCAll = 0, nRecAll=0, nClonesAll=0; + Int_t nMCRef = 0, nRecRef=0, nClonesRef=0; AliHLTTPCCATracker &slice = fTracker->Slices()[iSlice]; - int firstSliceHit = 0; + Int_t firstSliceHit = 0; for( ; firstSliceHitNHits(); firstSliceHit++){ if( fTracker->Hits()[firstSliceHit].ISlice()==iSlice ) break; } - int endSliceHit = firstSliceHit; + Int_t endSliceHit = firstSliceHit; for( ; endSliceHitNHits(); endSliceHit++){ if( fTracker->Hits()[endSliceHit].ISlice()!=iSlice ) break; @@ -269,8 +367,7 @@ void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag ) for (Int_t ic = 0; icNHits()<=0 ) lb[nla++]= l.fLab[1]; if( l.fLab[2]>=0 ) lb[nla++]= l.fLab[2]; } - //cout<<12<=0 && lmaxHits()[ih].ID()]; if( l.fLab[0]>=0 ) fMCTracks[l.fLab[0]].NHits()++; if( l.fLab[1]>=0 ) fMCTracks[l.fLab[1]].NHits()++; @@ -348,27 +444,31 @@ void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag ) } } - int traN = slice.NOutTracks(); - Int_t *traLabels = new Int_t[traN]; - Double_t *traPurity = new Double_t[traN]; + Int_t traN = slice.NOutTracks(); + Int_t *traLabels = 0; + Double_t *traPurity = 0; + traLabels = new Int_t[traN]; + traPurity = new Double_t[traN]; { for (Int_t itr=0; itrHits()[index].ID()]; + //cout<=0 ) lb[nla++]= l.fLab[0]; if(l.fLab[1]>=0 ) lb[nla++]= l.fLab[1]; if(l.fLab[2]>=0 ) lb[nla++]= l.fLab[2]; } sort( lb, lb+nla ); - int labmax = -1, labcur=-1, lmax = 0, lcurr=0; - for( int i=0; i=0 && lmaxHits()[index].ID()]; if( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax ) lmax++; } traLabels[itr] = labmax; traPurity[itr] = ( (nhits>0) ?double(lmax)/double(nhits) :0 ); + //cout<<"perf track "<=fNMCTracks){ + + for(Int_t itr=0; itr=fNMCTracks){ nGhost++; continue; } @@ -421,7 +523,6 @@ void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag ) if( mc.Set()>0 ) fhEffVsP->Fill(mc.P(), ( mc.NReconstructed()>0 ?1 :0)); } - if( traLabels ) delete[] traLabels; if( traPurity ) delete[] traPurity; @@ -489,10 +590,200 @@ void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag ) void AliHLTTPCCAPerformance::Performance() { // main routine for performance calculation + //SG!!! + /* + fStatNEvents=0; + fStatNRecTot=0; + fStatNRecOut=0; + fStatNGhost=0; + fStatNMCAll=0; + fStatNRecAll=0; + fStatNClonesAll=0; + fStatNMCRef=0; + fStatNRecRef=0; + fStatNClonesRef=0; + */ fStatNEvents++; - for( int islice=0; isliceNSlices(); islice++){ + for( Int_t islice=0; isliceNSlices(); islice++){ SlicePerformance(islice,0); } + + // global tracker performance + { + if( !fTracker ) return; + + Int_t nRecTot = 0, nGhost=0, nRecOut=0; + Int_t nMCAll = 0, nRecAll=0, nClonesAll=0; + Int_t nMCRef = 0, nRecRef=0, nClonesRef=0; + + // Select reconstructable MC tracks + + { + for (Int_t imc=0; imc=0 ) fMCTracks[l.fLab[0]].NHits()++; + if( l.fLab[1]>=0 ) fMCTracks[l.fLab[1]].NHits()++; + if( l.fLab[2]>=0 ) fMCTracks[l.fLab[2]].NHits()++; + } + + for (Int_t imc=0; imc= 50 && mc.P()>=.05 ){ + mc.Set() = 1; + nMCAll++; + if( mc.P()>=1. ){ + mc.Set() = 2; + nMCRef++; + } + } + } + } + + Int_t traN = fTracker->NTracks(); + Int_t *traLabels = 0; + Double_t *traPurity = 0; + traLabels = new Int_t[traN]; + traPurity = new Double_t[traN]; + { + for (Int_t itr=0; itrTracks()[itr]; + Int_t nhits = tCA.NHits(); + Int_t *lb = new Int_t[nhits*3]; + Int_t nla=0; + for( Int_t ihit=0; ihitTrackHits()[tCA.FirstHitRef()+ihit]; + AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()]; + if(l.fLab[0]>=0 ) lb[nla++]= l.fLab[0]; + if(l.fLab[1]>=0 ) lb[nla++]= l.fLab[1]; + if(l.fLab[2]>=0 ) lb[nla++]= l.fLab[2]; + } + sort( lb, lb+nla ); + Int_t labmax = -1, labcur=-1, lmax = 0, lcurr=0; + for( Int_t i=0; i=0 && lmax=0 && lmaxTrackHits()[tCA.FirstHitRef()+ihit]; + AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()]; + if( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax + ) lmax++; + } + traLabels[itr] = labmax; + traPurity[itr] = ( (nhits>0) ?double(lmax)/double(nhits) :0 ); + if( lb ) delete[] lb; + } + } + + nRecTot+= traN; + for(Int_t itr=0; itr=fNMCTracks){ + nGhost++; + continue; + } + + AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]]; + + mc.NReconstructed()++; + if( mc.Set()== 0 ) nRecOut++; + else{ + if( mc.NReconstructed()==1 ) nRecAll++; + else if(mc.NReconstructed() > mc.NTurns() ) nClonesAll++; + if( mc.Set()==2 ){ + if( mc.NReconstructed()==1 ) nRecRef++; + else if(mc.NReconstructed() > mc.NTurns() ) nClonesRef++; + } + } + + // track resolutions + while( TMath::Abs(mc.TPCPar()[0]) + TMath::Abs(mc.TPCPar()[1])>1 ){ + if( traPurity[itr]<.90 ) break; + AliHLTTPCCAGBTrack &t = fTracker->Tracks()[itr]; + AliHLTTPCCATrackParam p = t.Param(); + Double_t cosA = TMath::Cos( t.Alpha() ); + Double_t sinA = TMath::Sin( t.Alpha() ); + Double_t mcX = mc.TPCPar()[0]*cosA + mc.TPCPar()[1]*sinA; + Double_t mcY = -mc.TPCPar()[0]*sinA + mc.TPCPar()[1]*cosA; + Double_t mcZ = mc.TPCPar()[2]; + Double_t mcEx = mc.TPCPar()[3]*cosA + mc.TPCPar()[4]*sinA; + Double_t mcEy = -mc.TPCPar()[3]*sinA + mc.TPCPar()[4]*cosA; + Double_t mcEz = mc.TPCPar()[5]; + Double_t mcEt = TMath::Sqrt(mcEx*mcEx + mcEy*mcEy); + if( TMath::Abs(mcEt)<1.e-4 ) break; + Double_t mcSinPhi = mcEy / mcEt; + Double_t mcDzDs = mcEz / mcEt; + Double_t mcQPt = mc.TPCPar()[6]/ mcEt; + if( TMath::Abs(mcQPt)<1.e-4 ) break; + Double_t mcPt = 1./TMath::Abs(mcQPt); + if( mcPt<1. ) break; + if( t.NHits() < 50 ) break; + Double_t bz = fTracker->Slices()[0].Param().Bz(); + if( !p.TransportToXWithMaterial( mcX, bz ) ) break; + if( p.GetCosPhi()*mcEx < 0 ){ // change direction + mcSinPhi = -mcSinPhi; + mcDzDs = -mcDzDs; + mcQPt = -mcQPt; + } + const Double_t kCLight = 0.000299792458; + Double_t k2QPt = 100; + if( TMath::Abs(bz)>1.e-4 ) k2QPt= 1./(bz*kCLight); + Double_t qPt = p.GetKappa()*k2QPt; + Double_t pt = 100; + if( TMath::Abs(qPt) >1.e-4 ) pt = 1./TMath::Abs(qPt); + + fhResY->Fill( p.GetY() - mcY ); + fhResZ->Fill( p.GetZ() - mcZ ); + fhResSinPhi->Fill( p.GetSinPhi() - mcSinPhi ); + fhResDzDs->Fill( p.GetDzDs() - mcDzDs ); + fhResPt->Fill( ( pt - mcPt )/mcPt ); + + if( p.GetErr2Y()>0 ) fhPullY->Fill( (p.GetY() - mcY)/TMath::Sqrt(p.GetErr2Y()) ); + if( p.GetErr2Z()>0 ) fhPullZ->Fill( (p.GetZ() - mcZ)/TMath::Sqrt(p.GetErr2Z()) ); + if( p.GetErr2SinPhi()>0 ) fhPullSinPhi->Fill( (p.GetSinPhi() - mcSinPhi)/TMath::Sqrt(p.GetErr2SinPhi()) ); + if( p.GetErr2DzDs()>0 ) fhPullDzDs->Fill( (p.DzDs() - mcDzDs)/TMath::Sqrt(p.GetErr2DzDs()) ); + if( p.GetErr2Kappa()>0 ) fhPullQPt->Fill( (qPt - mcQPt)/TMath::Sqrt(p.GetErr2Kappa()*k2QPt*k2QPt) ); + break; + } + } + + for (Int_t ipart=0; ipart0 ) fhGBEffVsP->Fill(mc.P(), ( mc.NReconstructed()>0 ?1 :0)); + } + + if( traLabels ) delete[] traLabels; + if( traPurity ) delete[] traPurity; + + fStatGBNRecTot += nRecTot; + fStatGBNRecOut += nRecOut; + fStatGBNGhost += nGhost; + fStatGBNMCAll += nMCAll; + fStatGBNRecAll += nRecAll; + fStatGBNClonesAll += nClonesAll; + fStatGBNMCRef += nMCRef; + fStatGBNRecRef += nRecRef; + fStatGBNClonesRef += nClonesRef; + } + // distribution of cluster errors @@ -504,48 +795,302 @@ void AliHLTTPCCAPerformance::Performance() fhHitErrZ->Fill(hit.ErrZ()); } } + + // cluster pulls + + if( fDoClusterPulls && fNMCPoints>0 ) { + + { + for (Int_t ipart=0; ipartHits()[ih]; + AliHLTTPCCAHitLabel &l = fHitLabels[ih]; + + if( l.fLab[0]<0 || l.fLab[0]>=fNMCTracks + || l.fLab[1]>=0 || l.fLab[2]>=0 ) continue; + + Int_t lab = l.fLab[0]; + + AliHLTTPCCAMCTrack &track = fMCTracks[lab]; + //if( track.Pt()<1. ) continue; + Int_t ip1=-1, ip2=-1; + Double_t d1 = 1.e20, d2=1.e20; + for( Int_t ip=0; ip1.e-8 && TMath::Abs(p1.Sx()-hit.X())<2. && TMath::Abs(p2.Sx()-hit.X())<2. ){ + Double_t sx = hit.X(); + Double_t sy = p1.Sy() + dy/dx*(sx-p1.Sx()); + Double_t sz = p1.Sz() + dz/dx*(sx-p1.Sx()); + + Float_t errY, errZ; + { + AliHLTTPCCATrackParam t; + t.Z() = sz; + t.SinPhi() = dy/TMath::Sqrt(dx*dx+dy*dy); + t.CosPhi() = dx/TMath::Sqrt(dx*dx+dy*dy); + t.DzDs() = dz/TMath::Sqrt(dx*dx+dy*dy); + fTracker->GetErrors2(hit,t,errY, errZ ); + errY = TMath::Sqrt(errY); + errZ = TMath::Sqrt(errZ); + } + + fhHitResY->Fill((hit.Y()-sy)); + fhHitResZ->Fill((hit.Z()-sz)); + fhHitPullY->Fill((hit.Y()-sy)/errY); + fhHitPullZ->Fill((hit.Z()-sz)/errZ); + if( track.Pt()>=1. ){ + fhHitResY1->Fill((hit.Y()-sy)); + fhHitResZ1->Fill((hit.Z()-sz)); + fhHitPullY1->Fill((hit.Y()-sy)/errY); + fhHitPullZ1->Fill((hit.Z()-sz)/errZ); + } + } + } + } + + { + cout<<"\nSlice tracker performance: \n"<0 ) ? fStatNRecTot :1; + Double_t dMCAll = (fStatNMCAll>0 ) ? fStatNMCAll :1; + Double_t dMCRef = (fStatNMCRef>0 ) ? fStatNMCRef :1; + Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1; + Double_t dRecAll = (fStatNRecAll+fStatNClonesAll>0 ) ? fStatNRecAll+fStatNClonesAll :1; + Double_t dRecRef = (fStatNRecRef+fStatNClonesRef>0 ) ? fStatNRecRef+fStatNClonesRef :1; + Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1; + + cout<<" EffRef = "<< fStatNRecRef/dMCRef + <<", CloneRef = " << fStatNClonesRef/dRecRef <StatTime(2)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(3)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(4)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(5)/fTracker->StatNEvents()*1.e3<<"[" + <StatTime(6)/fTracker->StatNEvents()*1.e3<<"/" + <StatTime(7)/fTracker->StatNEvents()*1.e3<<"] " + <StatTime(8)/fTracker->StatNEvents()*1.e3<<" " + <<" msec/event "<0 ) ? fStatNRecTot :1; - Double_t dMCAll = (fStatNMCAll>0 ) ? fStatNMCAll :1; - Double_t dMCRef = (fStatNMCRef>0 ) ? fStatNMCRef :1; - Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1; - Double_t dRecAll = (fStatNRecAll+fStatNClonesAll>0 ) ? fStatNRecAll+fStatNClonesAll :1; - Double_t dRecRef = (fStatNRecRef+fStatNClonesRef>0 ) ? fStatNRecRef+fStatNClonesRef :1; - Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1; - - cout<<" EffRef = "<< fStatNRecRef/dMCRef - <<", CloneRef = " << fStatNClonesRef/dRecRef <StatTime(2)/fTracker->StatNEvents()*1.e3<<" " - <StatTime(3)/fTracker->StatNEvents()*1.e3<<" " - <StatTime(4)/fTracker->StatNEvents()*1.e3<<" " - <StatTime(5)/fTracker->StatNEvents()*1.e3<<"[" - <StatTime(6)/fTracker->StatNEvents()*1.e3<<"/" - <StatTime(7)/fTracker->StatNEvents()*1.e3<<"], merge:" - <StatTime(8)/fTracker->StatNEvents()*1.e3 - <<" msec/event "<0 ) ? fStatGBNRecTot :1; + Double_t dMCAll = (fStatGBNMCAll>0 ) ? fStatGBNMCAll :1; + Double_t dMCRef = (fStatGBNMCRef>0 ) ? fStatGBNMCRef :1; + Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1; + Double_t dRecAll = (fStatGBNRecAll+fStatGBNClonesAll>0 ) ? fStatGBNRecAll+fStatGBNClonesAll :1; + Double_t dRecRef = (fStatGBNRecRef+fStatGBNClonesRef>0 ) ? fStatGBNRecRef+fStatGBNClonesRef :1; + Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1; + + cout<<" EffRef = "<< fStatGBNRecRef/dMCRef + <<", CloneRef = " << fStatGBNClonesRef/dRecRef <StatTime(0)/fTracker->StatNEvents()*1.e3<<": " + <StatTime(1)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(2)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(3)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(4)/fTracker->StatNEvents()*1.e3<<" " + <StatTime(5)/fTracker->StatNEvents()*1.e3<<"[" + <StatTime(6)/fTracker->StatNEvents()*1.e3<<"/" + <StatTime(7)/fTracker->StatNEvents()*1.e3<<"] " + <StatTime(8)/fTracker->StatNEvents()*1.e3 + <<" msec/event "<StatTime(9)/fTracker->StatNEvents()*1.e3<<": " + <StatTime(10)/fTracker->StatNEvents()*1.e3<<", " + <StatTime(11)/fTracker->StatNEvents()*1.e3<<", " + <StatTime(12)/fTracker->StatNEvents()*1.e3<<" " + <<" msec/event "<>fNMCTracks; + fMCTracks = new AliHLTTPCCAMCTrack[fNMCTracks]; + for( Int_t it=0; it>j; + in>> t.PDG(); + for( Int_t i=0; i<7; i++ ) in>>t.Par()[i]; + for( Int_t i=0; i<7; i++ ) in>>t.TPCPar()[i]; + in>> t.P(); + in>> t.Pt(); + in>> t.NHits(); + in>> t.NMCPoints(); + in>> t.FirstMCPointID(); + in>> t.NReconstructed(); + in>> t.Set(); + in>> t.NTurns(); + } + + in>>fNHits; + fHitLabels = new AliHLTTPCCAHitLabel[fNHits]; + for( Int_t ih=0; ih>l.fLab[0]>>l.fLab[1]>>l.fLab[2]; + } +} + +void AliHLTTPCCAPerformance::ReadMCPoints( istream &in ) +{ + if( fMCPoints ) delete[] fMCPoints; + fMCPoints = 0; + fNMCPoints = 0; + + in>>fNMCPoints; + fMCPoints = new AliHLTTPCCAMCPoint[fNMCPoints]; + for( Int_t ip=0; ip> p.X(); + in>> p.Y(); + in>> p.Z(); + in>> p.Sx(); + in>> p.Sy(); + in>> p.Sz(); + in>> p.Time(); + in>> p.ISlice(); + in>> p.TrackID(); + } +} diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.h index 38982bd7c4d..97f53a5fd61 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.h @@ -13,9 +13,11 @@ class TParticle; class AliHLTTPCCAMCTrack; +class AliHLTTPCCAMCPoint; class AliHLTTPCCAGBTracker; class TDirectory; class TH1D; +class TH2D; class TProfile; /** @@ -41,15 +43,27 @@ class AliHLTTPCCAPerformance:public TObject void StartEvent(); void SetNHits( Int_t NHits ); void SetNMCTracks( Int_t NMCTracks ); + void SetNMCPoints( Int_t NMCPoints ); + void ReadHitLabel( Int_t HitID, Int_t lab0, Int_t lab1, Int_t lab2 ); void ReadMCTrack( Int_t index, const TParticle *part ); + void ReadMCTPCTrack( Int_t index, Float_t X, Float_t Y, Float_t Z, + Float_t Px, Float_t Py, Float_t Pz ); + + void ReadMCPoint( Int_t TrackID, Float_t X, Float_t Y, Float_t Z, Float_t Time, Int_t iSlice ); void CreateHistos(); void WriteHistos(); void SlicePerformance( Int_t iSlice, Bool_t PrintFlag ); void Performance(); + void WriteMCEvent( ostream &out ); + void ReadMCEvent( istream &in ); + void WriteMCPoints( ostream &out ); + void ReadMCPoints( istream &in ); + Bool_t& DoClusterPulls(){ return fDoClusterPulls; } + protected: AliHLTTPCCAGBTracker *fTracker; //* pointer to the tracker @@ -57,10 +71,14 @@ protected: struct AliHLTTPCCAHitLabel{ Int_t fLab[3]; //* array of 3 MC labels }; + AliHLTTPCCAHitLabel *fHitLabels; //* array of hit MC labels - Int_t fNHits; //* number of hits - AliHLTTPCCAMCTrack *fMCTracks; //* array of Mc tracks - Int_t fNMCTracks; //* number of MC tracks + Int_t fNHits; //* number of hits + AliHLTTPCCAMCTrack *fMCTracks; //* array of MC tracks + Int_t fNMCTracks; //* number of MC tracks + AliHLTTPCCAMCPoint *fMCPoints; //* array of MC points + Int_t fNMCPoints; //* number of MC points + Bool_t fDoClusterPulls; //* do cluster pulls (very slow) Int_t fStatNEvents; //* n of events proceed Int_t fStatNRecTot; //* total n of reconstructed tracks Int_t fStatNRecOut; //* n of reconstructed tracks in Out set @@ -72,17 +90,45 @@ protected: Int_t fStatNRecRef; //* n of reconstructed tracks in Ref set Int_t fStatNClonesRef; //* n of reconstructed clones in Ref set + Int_t fStatGBNRecTot; //* global tracker: total n of reconstructed tracks + Int_t fStatGBNRecOut; //* global tracker: n of reconstructed tracks in Out set + Int_t fStatGBNGhost;//* global tracker: n of reconstructed tracks in Ghost set + Int_t fStatGBNMCAll;//* global tracker: n of MC tracks + Int_t fStatGBNRecAll; //* global tracker: n of reconstructed tracks in All set + Int_t fStatGBNClonesAll;//* global tracker: total n of reconstructed tracks in Clone set + Int_t fStatGBNMCRef; //* global tracker: n of MC reference tracks + Int_t fStatGBNRecRef; //* global tracker: n of reconstructed tracks in Ref set + Int_t fStatGBNClonesRef; //* global tracker: n of reconstructed clones in Ref set + TDirectory *fHistoDir; //* ROOT directory with histogramms + TH1D + *fhResY, //* track Y resolution at the TPC entrance + *fhResZ, //* track Z resolution at the TPC entrance + *fhResSinPhi, //* track SinPhi resolution at the TPC entrance + *fhResDzDs, //* track DzDs resolution at the TPC entrance + *fhResPt, //* track Pt relative resolution at the TPC entrance + *fhPullY, //* track Y pull at the TPC entrance + *fhPullZ, //* track Z pull at the TPC entrance + *fhPullSinPhi, //* track SinPhi pull at the TPC entrance + *fhPullDzDs, //* track DzDs pull at the TPC entrance + *fhPullQPt; //* track Q/Pt pull at the TPC entrance + TH1D *fhHitErrY, //* hit error in Y - *fhHitErrZ,//* hit error in Z - *fhHitResX,//* hit resolution X + *fhHitErrZ,//* hit error in Z *fhHitResY,//* hit resolution Y *fhHitResZ,//* hit resolution Z - *fhHitPullX,//* hit pull X *fhHitPullY,//* hit pull Y - *fhHitPullZ,//* hit pull Z + *fhHitPullZ;//* hit pull Z + + TH1D + *fhHitResY1,//* hit resolution Y, pt>1GeV + *fhHitResZ1,//* hit resolution Z, pt>1GeV + *fhHitPullY1,//* hit pull Y, pt>1GeV + *fhHitPullZ1;//* hit pull Z, pt>1GeV + + TH1D *fhCellPurity,//* cell purity *fhCellNHits//* cell n hits ; @@ -90,7 +136,8 @@ protected: TProfile *fhCellPurityVsN, //* cell purity vs N hits *fhCellPurityVsPt,//* cell purity vs MC Pt - *fhEffVsP; //* reconstruction efficiency vs P plot + *fhEffVsP, //* reconstruction efficiency vs P plot + *fhGBEffVsP; //* global reconstruction efficiency vs P plot static void WriteDir2Current( TObject *obj ); diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.cxx new file mode 100644 index 00000000000..79608342256 --- /dev/null +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.cxx @@ -0,0 +1,83 @@ +// $Id: AliHLTTPCCATrackConvertor.cxx 27042 2008-07-02 12:06:02Z richterm $ +//*************************************************************************** +// This file is property of and copyright by the ALICE HLT Project * +// ALICE Experiment at CERN, All rights reserved. * +// * +// Primary Authors: Sergey Gorbunov * +// Ivan Kisel * +// for The ALICE HLT Project. * +// * +// Permission to use, copy, modify and distribute this software and its * +// documentation strictly for non-commercial purposes is hereby granted * +// without fee, provided that the above copyright notice appears in all * +// copies and that both the copyright notice and this permission notice * +// appear in the supporting documentation. The authors make no claims * +// about the suitability of this software for any purpose. It is * +// provided "as is" without express or implied warranty. * +//*************************************************************************** + +#include "AliHLTTPCCATrackConvertor.h" +#include "AliExternalTrackParam.h" +#include "AliHLTTPCCATrackParam.h" +#include "TMath.h" + + +void AliHLTTPCCATrackConvertor::GetExtParam( const AliHLTTPCCATrackParam &T1, AliExternalTrackParam &T2, Double_t alpha, Double_t Bz ) +{ + //* Convert from AliHLTTPCCATrackParam to AliExternalTrackParam parameterisation, + //* the angle alpha is the global angle of the local X axis + + Double_t par[5], cov[15]; + for( Int_t i=0; i<5; i++ ) par[i] = T1.GetPar()[i]; + for( Int_t i=0; i<15; i++ ) cov[i] = T1.GetCov()[i]; + + if(par[2]>.99 ) par[2]=.99; + if(par[2]<-.99 ) par[2]=-.99; + + { // kappa => 1/pt + const Double_t kCLight = 0.000299792458; + Double_t c = 1.e4; + if( TMath::Abs(Bz)>1.e-4 ) c = 1./(Bz*kCLight); + par[4] *= c; + cov[10]*= c; + cov[11]*= c; + cov[12]*= c; + cov[13]*= c; + cov[14]*= c*c; + } + if( T1.GetCosPhi()<0 ){ // change direction + par[2] = -par[2]; // sin phi + par[3] = -par[3]; // DzDs + par[4] = -par[4]; // kappa + cov[3] = -cov[3]; + cov[4] = -cov[4]; + cov[6] = -cov[6]; + cov[7] = -cov[7]; + cov[10] = -cov[10]; + cov[11] = -cov[11]; + } + T2.Set(T1.GetX(),alpha,par,cov); +} + +void AliHLTTPCCATrackConvertor::SetExtParam( AliHLTTPCCATrackParam &T1, const AliExternalTrackParam &T2, Double_t Bz ) +{ + //* Convert from AliExternalTrackParam parameterisation + + for( Int_t i=0; i<5; i++ ) T1.Par()[i] = T2.GetParameter()[i]; + for( Int_t i=0; i<15; i++ ) T1.Cov()[i] = T2.GetCovariance()[i]; + T1.X() = T2.GetX(); + if(T1.SinPhi()>.99 ) T1.SinPhi()=.99; + if(T1.SinPhi()<-.99 ) T1.SinPhi()=-.99; + T1.CosPhi() = TMath::Sqrt(1.-T1.SinPhi()*T1.SinPhi()); + const Double_t kCLight = 0.000299792458; + Double_t c = Bz*kCLight; + { // 1/pt -> kappa + T1.Par()[4] *= c; + T1.Cov()[10]*= c; + T1.Cov()[11]*= c; + T1.Cov()[12]*= c; + T1.Cov()[13]*= c; + T1.Cov()[14]*= c*c; + } +} + diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.h new file mode 100644 index 00000000000..0d3d8e27b70 --- /dev/null +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackConvertor.h @@ -0,0 +1,34 @@ +//-*- Mode: C++ -*- + +//* This file is property of and copyright by the ALICE HLT Project * +//* ALICE Experiment at CERN, All rights reserved. * +//* See cxx source for full Copyright notice * + + +#ifndef ALIHLTTPCCATRACKCONVERTOR_H +#define ALIHLTTPCCATRACKCONVERTOR_H +#include "Rtypes.h" + +class AliExternalTrackParam; +class AliHLTTPCCATrackParam; + +/** + * @class AliHLTTPCCATrackConvertor + * + * AliHLTTPCCATrackConvertor class converts tracks to different parameterisations + * it is used by the AliHLTTPCCATrackerComponent + * + */ +class AliHLTTPCCATrackConvertor +{ + public: + + AliHLTTPCCATrackConvertor(){} + + static void GetExtParam( const AliHLTTPCCATrackParam &T1, AliExternalTrackParam &T2, Double_t alpha, Double_t Bz ); + static void SetExtParam( AliHLTTPCCATrackParam &T1, const AliExternalTrackParam &T2, Double_t Bz ); + +}; + + +#endif diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx index 2b81e097e04..6a6d70fc5e8 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx @@ -18,7 +18,7 @@ #include "AliHLTTPCCATrackParam.h" #include "TMath.h" -#include "AliExternalTrackParam.h" +#include "Riostream.h" //ClassImp(AliHLTTPCCATrackParam) @@ -107,10 +107,10 @@ void AliHLTTPCCATrackParam::ConstructXY3( const Float_t x[3], const Float_t y[3] Float_t s0 = sigmaY2[0]; Float_t s1 = sigmaY2[1]; Float_t s2 = sigmaY2[2]; - + fC[0] = s0; fC[1] = 0; - fC[2] = 0; + fC[2] = 100.; fC[3] = d0[0]*s0; fC[4] = 0; @@ -119,7 +119,7 @@ void AliHLTTPCCATrackParam::ConstructXY3( const Float_t x[3], const Float_t y[3] fC[6] = 0; fC[7] = 0; fC[8] = 0; - fC[9] = 0; + fC[9] = 100.; fC[10] = d0[1]*s0; fC[11] = 0; @@ -269,7 +269,7 @@ Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x ) Float_t dSin = dl*k/2; if( dSin > 1 ) dSin = 1; if( dSin <-1 ) dSin = -1; - Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl; + Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl; Float_t dz = dS*DzDs(); Float_t cci = 0, exi = 0, ex1i = 0; @@ -281,6 +281,7 @@ Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x ) else ret = 0; if( !ret ) return ret; + X() += dx; fP[0]+= dy; fP[1]+= dz; @@ -338,34 +339,261 @@ Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x ) return ret; } +Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t X, Float_t Bz ) +{ + AliHLTTPCCATrackFitParam par; + CalculateFitParameters( par, Bz ); + return TransportToXWithMaterial(X, par ); +} -Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha ) +Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, AliHLTTPCCATrackFitParam &par ) { - //* Rotate the coordinate system in XY on the angle alpha + //* Transport the track parameters to X=x Bool_t ret = 1; + Float_t oldX=GetX(); + + Float_t x0 = X(); + //Float_t y0 = Y(); + Float_t k = Kappa(); + Float_t ex = CosPhi(); + Float_t ey = SinPhi(); + Float_t dx = x - x0; + + Float_t ey1 = k*dx + ey; + Float_t ex1; + if( TMath::Abs(ey1)>.99 ){ // no intersection -> check the border + ey1 = ( ey1>0 ) ?1 :-1; + ex1 = 0; + dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0; + + Float_t ddx = TMath::Abs(x0+dx - x)*k*k; + Float_t hx[] = {0, -k, 1+ey }; + Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5]; + if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error + ret = 0; // any case + return ret; + }else{ + ex1 = TMath::Sqrt(1 - ey1*ey1); + if( ex<0 ) ex1 = -ex1; + } + + Float_t dx2 = dx*dx; + CosPhi() = ex1; + Float_t ss = ey+ey1; + Float_t cc = ex+ex1; + Float_t tg = 0; + if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2) + else ret = 0; + Float_t dy = dx*tg; + Float_t dl = dx*TMath::Sqrt(1+tg*tg); + + if( cc<0 ) dl = -dl; + Float_t dSin = dl*k/2; + if( dSin > 1 ) dSin = 1; + if( dSin <-1 ) dSin = -1; + Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl; + Float_t dz = dS*DzDs(); + + Float_t cci = 0, exi = 0, ex1i = 0; + if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc; + else ret = 0; + if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex; + else ret = 0; + if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1; + else ret = 0; + + if( !ret ) return ret; + + X() += dx; + fP[0]+= dy; + fP[1]+= dz; + fP[2] = ey1; + fP[3] = fP[3]; + fP[4] = fP[4]; + + Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i; + Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci; + + Float_t c00 = fC[0]; + Float_t c10 = fC[1]; + Float_t c11 = fC[2]; + Float_t c20 = fC[3]; + Float_t c21 = fC[4]; + Float_t c22 = fC[5]; + Float_t c30 = fC[6]; + Float_t c31 = fC[7]; + Float_t c32 = fC[8]; + Float_t c33 = fC[9]; + Float_t c40 = fC[10]; + Float_t c41 = fC[11]; + Float_t c42 = fC[12]; + Float_t c43 = fC[13]; + Float_t c44 = fC[14]; + + //Float_t H0[5] = { 1,0, h2, 0, h4 }; + //Float_t H1[5] = { 0, 1, 0, dS, 0 }; + //Float_t H2[5] = { 0, 0, 1, 0, dx }; + //Float_t H3[5] = { 0, 0, 0, 1, 0 }; + //Float_t H4[5] = { 0, 0, 0, 0, 1 }; + + + fC[0]=( c00 + h2*h2*c22 + h4*h4*c44 + + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) ); + + fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43); + fC[2]= c11 + 2*dS*c31 + dS*dS*c33; + + fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44); + fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43); + fC[5]= c22 +2*dx*c42 + dx2*c44; + + fC[6]= c30 + h2*c32 + h4*c43; + fC[7]= c31 + dS*c33; + fC[8]= c32 + dx*c43; + fC[9]= c33; + + fC[10]= c40 + h2*c42 + h4*c44; + fC[11]= c41 + dS*c43; + fC[12]= c42 + dx*c44; + fC[13]= c43; + fC[14]= c44; + + Float_t d = TMath::Sqrt(dS*dS + dz*dz ); + + if (oldX > GetX() ) d = -d; + { + Float_t rho=0.9e-3; + Float_t radLen=28.94; + CorrectForMeanMaterial(d*rho/radLen,d*rho,par); + } + + return ret; +} + + + +Float_t AliHLTTPCCATrackParam::ApproximateBetheBloch( Float_t beta2 ) +{ + //------------------------------------------------------------------ + // This is an approximation of the Bethe-Bloch formula with + // the density effect taken into account at beta*gamma > 3.5 + // (the approximation is reasonable only for solid materials) + //------------------------------------------------------------------ + if (beta2 >= 1) return 0; + + if (beta2/(1-beta2)>3.5*3.5) + return 0.153e-3/beta2*( log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2); + return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2); +} + + +void AliHLTTPCCATrackParam::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, Float_t Bz, Float_t mass ) +{ + //*! + + const Float_t kCLight = 0.000299792458; + Float_t c = Bz*kCLight; + Float_t p2 = (1.+ fP[3]*fP[3])*c*c; + Float_t k2 = fP[4]*fP[4]; + Float_t beta2= p2 / (p2 + mass*mass*k2); + Float_t Bethe = ApproximateBetheBloch(beta2); + + Float_t P2 = (k2>1.e-8) ?p2/k2 :10000; // impuls 2 + par.fBethe = Bethe; + par.fE = TMath::Sqrt( P2 + mass*mass); + par.fTheta2 = 14.1*14.1/(beta2*p2*1e6)*k2; + par.fEP2 = par.fE/p2*k2; + + // Approximate energy loss fluctuation (M.Ivanov) + + const Float_t knst=0.07; // To be tuned. + par.fSigmadE2 = knst*par.fEP2*fP[4]; + par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2; + + par.fK22 = (1. + fP[3]*fP[3]); + par.fK33 = par.fK22*par.fK22; + par.fK43 = fP[3]*fP[4]*par.fK22; + par.fK44 = fP[3]*fP[3]*fP[4]*fP[4]; +} + + +Bool_t AliHLTTPCCATrackParam::CorrectForMeanMaterial( Float_t xOverX0, Float_t xTimesRho, AliHLTTPCCATrackFitParam &par ) +{ + //------------------------------------------------------------------ + // This function corrects the track parameters for the crossed material. + // "xOverX0" - X/X0, the thickness in units of the radiation length. + // "xTimesRho" - is the product length*density (g/cm^2). + //------------------------------------------------------------------ + + Float_t &fC22=fC[5]; + Float_t &fC33=fC[9]; + Float_t &fC40=fC[10]; + Float_t &fC41=fC[11]; + Float_t &fC42=fC[12]; + Float_t &fC43=fC[13]; + Float_t &fC44=fC[14]; + + //Energy losses************************ + + Float_t dE = par.fBethe*xTimesRho; + if ( TMath::Abs(dE) > 0.3*par.fE ) return 0; //30% energy loss is too much! + Float_t corr = (1.- par.fEP2*dE); + if( corr<0.3 ) return 0; + fP[4]*= corr; + fC40*= corr; + fC41*= corr; + fC42*= corr; + fC43*= corr; + fC44*= corr*corr; + fC44+= par.fSigmadE2*TMath::Abs(dE); + + + //Multiple scattering****************** + + Float_t theta2 = par.fTheta2*TMath::Abs(xOverX0); + fC22 += theta2*par.fK22*(1.- fP[2]*fP[2]); + fC33 += theta2*par.fK33; + fC43 += theta2*par.fK43; + fC44 += theta2*par.fK44; + + return 1; +} + + + +#include "Riostream.h" + +Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha ) +{ + //* Rotate the coordinate system in XY on the angle alpha + Float_t cA = TMath::Cos( alpha ); Float_t sA = TMath::Sin( alpha ); Float_t x = X(), y= Y(), sP= SinPhi(), cP= CosPhi(); - + Float_t cosPhi = cP*cA + sP*sA; + Float_t sinPhi =-cP*sA + sP*cA; + + if( TMath::Abs(sinPhi)>.99 || TMath::Abs(cosPhi)<1.e-2 || TMath::Abs(cP)<1.e-2 ) return 0; + + Float_t j0 = cP/cosPhi; + Float_t j2 = cosPhi/cP; + X() = x*cA + y*sA; Y() = -x*sA + y*cA; - CosPhi() = cP*cA + sP*sA; - SinPhi() = -cP*sA + sP*cA; + CosPhi() = cosPhi; + SinPhi() = sinPhi; - Float_t j0 = 0, j2 = 0; - - if( TMath::Abs(CosPhi())>1.e-4 ) j0 = cP/CosPhi(); else ret = 0; - if( TMath::Abs(cP)>1.e-4 ) j2 = CosPhi()/cP; else ret = 0; //Float_t J[5][5] = { { j0, 0, 0, 0, 0 }, // Y // { 0, 1, 0, 0, 0 }, // Z // { 0, 0, j2, 0, 0 }, // SinPhi // { 0, 0, 0, 1, 0 }, // DzDs // { 0, 0, 0, 0, 1 } }; // Kappa - + //cout<<"alpha="<.999 ) par[2]=.999; - if(par[2]<-.999 ) par[2]=-.999; - - const Double_t kCLight = 0.000299792458; - Double_t c = (TMath::Abs(Bz)>1.e-4) ?1./(Bz*kCLight) :1./(5.*kCLight); - { // kappa => 1/pt - par[4] *= c; - cov[10]*= c; - cov[11]*= c; - cov[12]*= c; - cov[13]*= c; - cov[14]*= c*c; - } - if( GetCosPhi()<0 ){ // change direction - par[2] = -par[2]; // sin phi - par[3] = -par[3]; // DzDs - par[4] = -par[4]; // kappa - cov[3] = -cov[3]; - cov[4] = -cov[4]; - cov[6] = -cov[6]; - cov[7] = -cov[7]; - cov[10] = -cov[10]; - cov[11] = -cov[11]; - } - T.Set(GetX(),alpha,par,cov); -} - -void AliHLTTPCCATrackParam::SetExtParam( const AliExternalTrackParam &T, Double_t Bz ) -{ - //* Convert from AliExternalTrackParam parameterisation - - for( Int_t i=0; i<5; i++ ) fP[i] = T.GetParameter()[i]; - for( Int_t i=0; i<15; i++ ) fC[i] = T.GetCovariance()[i]; - X() = T.GetX(); - if(SinPhi()>.999 ) SinPhi()=.999; - if(SinPhi()<-.999 ) SinPhi()=-.999; - CosPhi() = TMath::Sqrt(1.-SinPhi()*SinPhi()); - const Double_t kCLight = 0.000299792458; - Double_t c = Bz*kCLight; - { // 1/pt -> kappa - fP[4] *= c; - fC[10]*= c; - fC[11]*= c; - fC[12]*= c; - fC[13]*= c; - fC[14]*= c*c; - } -} - - -void AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t erry, Float_t errz ) +Bool_t AliHLTTPCCATrackParam::Filter2( Float_t y, Float_t z, Float_t err2Y, Float_t err2Z ) { //* Add the y,z measurement with the Kalman filter @@ -455,23 +625,19 @@ void AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t erry, Float_t z0 = y-fP[0], z1 = z-fP[1]; - Float_t v[3] = {erry*erry, 0, errz*errz}; + Float_t v[3] = {err2Y, 0, err2Z}; Float_t mS[3] = { c00+v[0], c10+v[1], c11+v[2] }; Float_t mSi[3]; Float_t det = (mS[0]*mS[2] - mS[1]*mS[1]); - if( TMath::Abs(det)<1.e-8 ) return; + if( det < 1.e-8 ) return 0; det = 1./det; mSi[0] = mS[2]*det; mSi[1] = -mS[1]*det; mSi[2] = mS[0]*det; - fNDF += 2; - fChi2 += ( +(mSi[0]*z0 + mSi[1]*z1 )*z0 - +(mSi[1]*z0 + mSi[2]*z1 )*z1 ); - // K = CHtS Float_t k00, k01 , k10, k11, k20, k21, k30, k31, k40, k41; @@ -483,11 +649,15 @@ void AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t erry, Float_t k40 = c40*mSi[0] + c41*mSi[1]; k41 = c40*mSi[1] + c41*mSi[2] ; Float_t sinPhi = fP[2] + k20*z0 + k21*z1 ; - if( TMath::Abs(sinPhi)>=0.99 ) return; + if( TMath::Abs(sinPhi)>=0.99 ) return 0; + + fNDF += 2; + fChi2 += ( +(mSi[0]*z0 + mSi[1]*z1 )*z0 + +(mSi[1]*z0 + mSi[2]*z1 )*z1 ); fP[ 0]+= k00*z0 + k01*z1 ; fP[ 1]+= k10*z0 + k11*z1 ; - fP[ 2]+= k20*z0 + k21*z1 ; + fP[ 2] = sinPhi; fP[ 3]+= k30*z0 + k31*z1 ; fP[ 4]+= k40*z0 + k41*z1 ; @@ -517,5 +687,150 @@ void AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t erry, Float_t }else{ CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); } + return 1; +} + +void AliHLTTPCCATrackParam::FilterY( Float_t y, Float_t erry ) +{ + //* Add the y measurement with the Kalman filter + + Float_t + c00 = fC[ 0], + c10 = fC[ 1], + c20 = fC[ 3], + c30 = fC[ 6], + c40 = fC[10]; + + Float_t + z0 = y-fP[0]; + + Float_t s = { c00+erry*erry }; + if( TMath::Abs(s)<1.e-4 ) return; + + Float_t si = 1/s; + + fNDF += 1; + fChi2 += si*z0*z0; + + // K = CHtS + + Float_t k0, k1 , k2, k3, k4; + k0 = c00*si; + k1 = c10*si; + k2 = c20*si; + k3 = c30*si; + k4 = c40*si; + + Float_t sinPhi = fP[2] + k2*z0 ; + if( TMath::Abs(sinPhi)>=0.99 ) return; + + fP[ 0]+= k0*z0 ; + fP[ 1]+= k1*z0 ; + fP[ 2] = sinPhi; + fP[ 3]+= k3*z0 ; + fP[ 4]+= k4*z0 ; + + fC[ 0]-= k0*c00; + + fC[ 1]-= k1*c00; + fC[ 2]-= k1*c10; + + fC[ 3]-= k2*c00; + fC[ 4]-= k2*c10; + fC[ 5]-= k2*c20; + + fC[ 6]-= k3*c00; + fC[ 7]-= k3*c10; + fC[ 8]-= k3*c20; + fC[ 9]-= k3*c30; + + fC[10]-= k4*c00; + fC[11]-= k4*c10; + fC[12]-= k4*c20; + fC[13]-= k4*c30; + fC[14]-= k4*c40; + + if( CosPhi()>=0 ){ + CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi()); + }else{ + CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); + } + +} + +void AliHLTTPCCATrackParam::FilterZ( Float_t z, Float_t errz ) +{ + //* Add the z measurement with the Kalman filter + + Float_t + c01 = fC[ 1], + c11 = fC[ 2], + c21 = fC[ 4], + c31 = fC[ 7], + c41 = fC[11]; + + Float_t + z1 = z-fP[1]; + + Float_t s = c11 + errz*errz; + if( TMath::Abs(s)<1.e-4 ) return; + + Float_t si = 1./s; + + fNDF += 1; + fChi2 += si*z1*z1; + + // K = CHtS + + Float_t k0, k1 , k2, k3, k4; + + k0 = 0;//c01*si; + k1 = c11*si; + k2 = 0;//c21*si; + k3 = c31*si; + k4 = 0;//c41*si; + + Float_t sinPhi = fP[2] + k2*z1 ; + if( TMath::Abs(sinPhi)>=0.99 ) return; + + fP[ 0]+= k0*z1 ; + fP[ 1]+= k1*z1 ; + fP[ 2] = sinPhi; + fP[ 3]+= k3*z1 ; + fP[ 4]+= k4*z1 ; + + + fC[ 0]-= k0*c01 ; + + fC[ 1]-= k1*c01 ; + fC[ 2]-= k1*c11 ; + + fC[ 3]-= k2*c01 ; + fC[ 4]-= k2*c11 ; + fC[ 5]-= k2*c21 ; + + fC[ 6]-= k3*c01 ; + fC[ 7]-= k3*c11 ; + fC[ 8]-= k3*c21 ; + fC[ 9]-= k3*c31 ; + + fC[10]-= k4*c01 ; + fC[11]-= k4*c11 ; + fC[12]-= k4*c21 ; + fC[13]-= k4*c31 ; + fC[14]-= k4*c41 ; + + if( CosPhi()>=0 ){ + CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi()); + }else{ + CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi()); + } + +} + +void AliHLTTPCCATrackParam::Print() const +{ + cout<<"track: "<=fParam.NRows() ) lastRow2 = fParam.NRows()-1; for (Int_t i1 = 0; i10 ) AliHLTTPCCADisplay::Instance().Ask(); @@ -623,7 +625,7 @@ void AliHLTTPCCATracker::FindTracks() TStopwatch timer4; - Bool_t doMerging=1; + Bool_t doMerging=1;//SG!!! Float_t factor2 = fParam.TrackConnectionFactor()*fParam.TrackConnectionFactor(); Int_t *refEndPoints = new Int_t[fNHitsTotal]; @@ -720,7 +722,11 @@ void AliHLTTPCCATracker::FindTracks() Float_t dy2, dz2; AliHLTTPCCATrackParam &jPar = jp.Param(); - + + // check direction + { + if( jPar.GetCosPhi()*iPar.GetCosPhi()>=0 ) continue; + } // check for neighbouring { float d = jPar.GetY() - iPar.GetY(); @@ -728,22 +734,27 @@ void AliHLTTPCCATracker::FindTracks() if( d*d>factor2*s2 ){ continue; } + //cout<<"\ndy="<0 ) AliHLTTPCCADisplay::Instance().Ask(); + //if( fNTracks>0 ) AliHLTTPCCADisplay::Instance().Ask(); #endif -#ifdef DRAW - AliHLTTPCCADisplay::Instance().Clear(); - AliHLTTPCCADisplay::Instance().DrawSector( this ); +#ifdef DRAWXX + AliHLTTPCCADisplay::Instance().ClearView(); + AliHLTTPCCADisplay::Instance().DrawSlice( this ); for( Int_t iRow=0; iRowLink()>=0; iID = ic->Link(), ic = &ID2Cell(iID) ){ //cout<<"itr="<NHits(); iHit++ ){ AliHLTTPCCAHit &h = row.GetCellHit(*ic,iHit); - // check for wrong hits - + // check for wrong hits { - if( !t.TransportToX( row.X() ) ) continue; - Float_t dy = t.GetY() - h.Y(); - Float_t dz = t.GetZ() - h.Z(); - if( dy*dy > 3.5*3.5*(t.GetErr2Y() + h.ErrY()*h.ErrY() ) ) continue; - if( dz*dz > 3.5*3.5*(t.GetErr2Z() + h.ErrZ()*h.ErrZ() ) ) continue; - /* - Float_t m[3] = {row.X(), h.Y(), h.Z() }; - Float_t mV[6] = {fParam.ErrX()*fParam.ErrX(), 0, h.ErrY()*h.ErrY(), 0, 0, h.ErrZ()*h.ErrZ() }; - Float_t mG[6]; - t.TransportBz(fParam.Bz(),m); - t.GetConnectionMatrix(fParam.Bz(),m, mG); - Float_t chi2 = t.GetChi2( m, mV, mG )/2.; - if( chi2>10. ) continue; - */ - } - - fOutTrackHits[fNOutTrackHits] = h.ID(); - fNOutTrackHits++; - if( fNOutTrackHits>fNHitsTotal ){ - //cout<<"fNOutTrackHits>fNHitsTotal"< 3.5*3.5*(/*t0.GetErr2Y() + */h.ErrY()*h.ErrY() ) ) continue;//SG!!! + //if( dz*dz > 3.5*3.5*(/*t0.GetErr2Z() + */h.ErrZ()*h.ErrZ() ) ) continue; + //if( !t0.Filter2( h.Y(), h.Z(), h.ErrY()*h.ErrY(), h.ErrZ()*h.ErrZ() ) ) continue; + + + if( !t.TransportToX( row.X() ) ) continue; + + //* Update the track + + if( first ){ + t.Cov()[ 0] = .5*.5; + t.Cov()[ 1] = 0; + t.Cov()[ 2] = .5*.5; + t.Cov()[ 3] = 0; + t.Cov()[ 4] = 0; + t.Cov()[ 5] = .2*.2; + t.Cov()[ 6] = 0; + t.Cov()[ 7] = 0; + t.Cov()[ 8] = 0; + t.Cov()[ 9] = .2*.2; + t.Cov()[10] = 0; + t.Cov()[11] = 0; + t.Cov()[12] = 0; + t.Cov()[13] = 0; + t.Cov()[14] = .2*.2; + t.Chi2() = 0; + t.NDF() = -5; + } + + if( t.Filter2( h.Y(), h.Z(), h.ErrY()*h.ErrY(), h.ErrZ()*h.ErrZ() ) ) first = 0; + else continue; + + fOutTrackHits[fNOutTrackHits] = h.ID(); + fNOutTrackHits++; + if( fNOutTrackHits>fNHitsTotal ){ + cout<<"fNOutTrackHits>fNHitsTotal"<= argc ) { - Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing B-field", "Missing B-field specifier." ); + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing solenoidBz", "Missing solenoidBz specifier." ); return ENOTSUP; } - fBField = strtod( argv[i+1], &cpErr ); + fSolenoidBz = strtod( argv[i+1], &cpErr ); if ( *cpErr ) { - Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert B-field specifier '%s'.", argv[i+1] ); + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert solenoidBz specifier '%s'.", argv[i+1] ); return EINVAL; } Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line", - "Magnetic field value is set to %f kG", fBField ); + "Magnetic field value is set to %f kG", fSolenoidBz ); i += 2; continue; } - if ( !strcmp( argv[i], "MinNTrackClusters" ) ){ + if ( !strcmp( argv[i], "minNTrackClusters" ) || !strcmp( argv[i], "-minNTrackClusters" ) ){ if ( i+1 >= argc ) { - Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing MinNTrackClusters", "Missing MinNTrackClusters specifier." ); + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing minNTrackClusters", "Missing minNTrackClusters specifier." ); return ENOTSUP; } fMinNTrackClusters = (Int_t ) strtod( argv[i+1], &cpErr ); if ( *cpErr ) { - Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert MinNTrackClusters '%s'.", argv[i+1] ); + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert minNTrackClusters '%s'.", argv[i+1] ); return EINVAL; } Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line", - "MinNTrackClusters is set to %i ", fMinNTrackClusters ); + "minNTrackClusters is set to %i ", fMinNTrackClusters ); i += 2; continue; } - + + if ( !strcmp( argv[i], "cellConnectionAngleXY" ) || !strcmp( argv[i], "-cellConnectionAngleXY" ) ){ + if ( i+1 >= argc ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing cellConnectionAngleXY", "Missing cellConnectionAngleXY specifier." ); + return ENOTSUP; + } + fCellConnectionAngleXY = strtod( argv[i+1], &cpErr ); + if ( *cpErr ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert cellConnectionAngleXY '%s'.", argv[i+1] ); + return EINVAL; + } + + Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line", + "cellConnectionAngleXY is set to %f ", fCellConnectionAngleXY ); + + i += 2; + continue; + } + if ( !strcmp( argv[i], "cellConnectionAngleXZ" ) || !strcmp( argv[i], "-cellConnectionAngleXZ" ) ){ + if ( i+1 >= argc ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing cellConnectionAngleXZ", "Missing cellConnectionAngleXZ specifier." ); + return ENOTSUP; + } + fCellConnectionAngleXZ = strtod( argv[i+1], &cpErr ); + if ( *cpErr ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert cellConnectionAngleXZ '%s'.", argv[i+1] ); + return EINVAL; + } + + Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line", + "cellConnectionAngleXZ is set to %f ", fCellConnectionAngleXZ ); + + i += 2; + continue; + } + if ( !strcmp( argv[i], "clusterZCut" ) || !strcmp( argv[i], "-clusterZCut" ) ){ + if ( i+1 >= argc ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing clusterZCut", "Missing clusterZCut specifier." ); + return ENOTSUP; + } + fClusterZCut = TMath::Abs(strtod( argv[i+1], &cpErr )); + if ( *cpErr ) + { + Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert clusterZCut '%s'.", argv[i+1] ); + return EINVAL; + } + + Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line", + "clusterZCut is set to %f ", fClusterZCut ); + + i += 2; + continue; + } + Logging(kHLTLogError, "HLT::TPCCATracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] ); return EINVAL; } @@ -204,6 +270,14 @@ int AliHLTTPCCATrackerComponent::DoDeinit() return 0; } +Bool_t AliHLTTPCCATrackerComponent::CompareClusters(AliHLTTPCSpacePointData *a, AliHLTTPCSpacePointData *b) +{ + //* Comparison function for sorting clusters + if( a->fPadRowfPadRow ) return 1; + if( a->fPadRow>b->fPadRow ) return 0; + return (a->fZ < b->fZ); +} + int AliHLTTPCCATrackerComponent::DoEvent ( const AliHLTComponentEventData& evtData, @@ -217,10 +291,14 @@ int AliHLTTPCCATrackerComponent::DoEvent AliHLTUInt32_t MaxBufferSize = size; size = 0; // output size + if(GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR )){ + return 0; + } + TStopwatch timer; // Event reconstruction in one TPC slice with CA Tracker - + //Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA::DoEvent()" ); if ( evtData.fBlockCnt<=0 ) { @@ -297,34 +375,34 @@ int AliHLTTPCCATrackerComponent::DoEvent // Initialize the tracker - Double_t Bz = fBField; + Float_t Bz = fSolenoidBz; { if( !fTracker ) fTracker = new AliHLTTPCCATracker; Int_t iSec = slice; - Double_t inRmin = 83.65; - // Double_t inRmax = 133.3; - // Double_t outRmin = 133.5; - Double_t outRmax = 247.7; - Double_t plusZmin = 0.0529937; - Double_t plusZmax = 249.778; - Double_t minusZmin = -249.645; - Double_t minusZmax = -0.0799937; - Double_t dalpha = 0.349066; - Double_t alpha = 0.174533 + dalpha*iSec; + Float_t inRmin = 83.65; + // Float_t inRmax = 133.3; + // Float_t outRmin = 133.5; + Float_t outRmax = 247.7; + Float_t plusZmin = 0.0529937; + Float_t plusZmax = 249.778; + Float_t minusZmin = -249.645; + Float_t minusZmax = -0.0799937; + Float_t dalpha = 0.349066; + Float_t alpha = 0.174533 + dalpha*iSec; Bool_t zPlus = (iSec<18 ); - Double_t zMin = zPlus ?plusZmin :minusZmin; - Double_t zMax = zPlus ?plusZmax :minusZmax; + Float_t zMin = zPlus ?plusZmin :minusZmin; + Float_t zMax = zPlus ?plusZmax :minusZmax; //TPCZmin = -249.645, ZMax = 249.778 - // Double_t rMin = inRmin; - // Double_t rMax = outRmax; + // Float_t rMin = inRmin; + // Float_t rMax = outRmax; Int_t NRows = AliHLTTPCTransform::GetNRows(); - Double_t padPitch = 0.4; - Double_t sigmaZ = 0.228808; + Float_t padPitch = 0.4; + Float_t sigmaZ = 0.228808; - Double_t *rowX = new Double_t [NRows]; + Float_t *rowX = new Float_t [NRows]; for( Int_t irow=0; irowInitialize( param ); delete[] rowX; } @@ -380,16 +460,12 @@ int AliHLTTPCCATrackerComponent::DoEvent fTracker->StartEvent(); - AliHLTTPCCAHit *vHits = new AliHLTTPCCAHit [nHitsTotal]; // CA hit array - Double_t *vHitStoreX = new Double_t [nHitsTotal]; // hit X coordinates - Int_t *vHitStoreID = new Int_t [nHitsTotal]; // hit ID's - Int_t *vHitRowID = new Int_t [nHitsTotal]; // hit ID's - - Int_t nHits = 0; - Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits", "Total %d hits to read for slice %d", nHitsTotal, slice ); + + AliHLTTPCSpacePointData** vOrigClusters = new AliHLTTPCSpacePointData* [ nHitsTotal]; + Int_t nClusters=0; for( std::vector::iterator pIter = patchIndices.begin(); pIter!=patchIndices.end(); pIter++ ){ @@ -402,26 +478,47 @@ int AliHLTTPCCATrackerComponent::DoEvent Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits", "Reading %d hits for slice %d - patch %d", inPtrSP->fSpacePointCnt, slice, patch ); - // Read patch hits, row by row + for (UInt_t i=0; ifSpacePointCnt; i++ ){ + vOrigClusters[nClusters++] = &(inPtrSP->fSpacePoints[i]); + } + } + // sort clusters since they are not sorted fore some reason + + sort( vOrigClusters, vOrigClusters+nClusters, CompareClusters ); + + AliHLTTPCCAHit *vHits = new AliHLTTPCCAHit [nHitsTotal]; // CA hit array + Float_t *vHitStoreX = new Float_t [nHitsTotal]; // hit X coordinates + Float_t *vHitStoreY = new Float_t [nHitsTotal]; // hit Y coordinates + Float_t *vHitStoreZ = new Float_t [nHitsTotal]; // hit Z coordinates + Int_t *vHitStoreID = new Int_t [nHitsTotal]; // hit ID's + Int_t *vHitRowID = new Int_t [nHitsTotal]; // hit ID's + + Int_t nHits = 0; + + { Int_t oldRow = -1; Int_t nRowHits = 0; Int_t firstRowHit = 0; - for (UInt_t i=0; ifSpacePointCnt; i++ ){ - AliHLTTPCSpacePointData* pSP = &(inPtrSP->fSpacePoints[i]); + for (Int_t i=0; ifPadRow != oldRow ){ - if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits ); + if( oldRow>=0 ){ + if( fTracker->Rows()[oldRow].NHits()!=0 ) HLTError("CA: clusters from row %d are readed twice",oldRow); + fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits ); + } oldRow = pSP->fPadRow; firstRowHit = nHits; nRowHits = 0; } AliHLTTPCCAHit &h = vHits[nHits]; - if( TMath::Abs(pSP->fX- fTracker->Rows()[pSP->fPadRow].X() )>1.e-4 ) cout<<"row "<<(Int_t)pSP->fPadRow<<" "<Rows()[pSP->fPadRow].X()-pSP->fX <fX- fTracker->Rows()[pSP->fPadRow].X() )>1.e-4 ) cout<<"row "<<(Int_t)pSP->fPadRow<<" "<Rows()[pSP->fPadRow].X()-pSP->fX <fX- fTracker->Rows()[pSP->fPadRow].X() )>1.e-4 ) HLTError( "row %d, %f",(Int_t)pSP->fPadRow, fTracker->Rows()[pSP->fPadRow].X()-pSP->fX ); h.Y() = pSP->fY; h.Z() = pSP->fZ; - if( TMath::Abs(h.Z())>230.) continue; + if( TMath::Abs(h.Z())>fClusterZCut) continue; h.ErrY() = TMath::Sqrt(TMath::Abs(pSP->fSigmaY2)); h.ErrZ() = TMath::Sqrt(TMath::Abs(pSP->fSigmaZ2)); if( h.ErrY()<.1 ) h.ErrY() = .1; @@ -429,15 +526,20 @@ int AliHLTTPCCATrackerComponent::DoEvent if( h.ErrY()>1. ) h.ErrY() = 1.; if( h.ErrZ()>1. ) h.ErrZ() = 1.; h.ID() = nHits; - vHitStoreX[nHits] = pSP->fX; + vHitStoreX[nHits] = pSP->fX; + vHitStoreY[nHits] = h.Y(); + vHitStoreZ[nHits] = h.Z(); vHitStoreID[nHits] = pSP->fID; vHitRowID[nHits] = pSP->fPadRow; nHits++; nRowHits++; - nClusters++; } - if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits ); + if( oldRow>=0 ){ + if( fTracker->Rows()[oldRow].NHits()!=0 ) HLTError("CA: clusters from row %d are readed twice",oldRow); + fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits ); + } } + if( vOrigClusters ) delete[] vOrigClusters; // reconstruct the event @@ -483,38 +585,57 @@ int AliHLTTPCCATrackerComponent::DoEvent } // convert CA track parameters to HLT Track Segment - + + Int_t iFirstRow = 1000; + Int_t iLastRow = -1; Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()]; - Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1]; - + Int_t iLastHit = iFirstHit; + for( Int_t ih=0; ihOutTrackHits()[t.FirstHitRef() + ih ]; + Int_t iRow = vHitRowID[hitID]; + if( iRowiLastRow ){ iLastRow = iRow; iLastHit = hitID; } + } + AliHLTTPCCATrackParam par = t.StartPoint(); par.TransportToX( vHitStoreX[iFirstHit] ); AliExternalTrackParam tp; - par.GetExtParam( tp, 0, fBField ); + AliHLTTPCCATrackConvertor::GetExtParam( par, tp, 0, fSolenoidBz ); currOutTracklet->fX = tp.GetX(); currOutTracklet->fY = tp.GetY(); currOutTracklet->fZ = tp.GetZ(); currOutTracklet->fCharge = (Int_t ) tp.GetSign(); currOutTracklet->fPt = TMath::Abs(tp.GetSignedPt()); - Double_t snp = tp.GetSnp() ; + Float_t snp = tp.GetSnp() ; if( snp>.999 ) snp=.999; if( snp<-.999 ) snp=-.999; currOutTracklet->fPsi = TMath::ASin( snp ); currOutTracklet->fTgl = tp.GetTgl(); - Double_t h = -currOutTracklet->fPt*currOutTracklet->fPt; + + currOutTracklet->fY0err = tp.GetSigmaY2(); + currOutTracklet->fZ0err = tp.GetSigmaZ2(); + Float_t h = -currOutTracklet->fPt*currOutTracklet->fPt; currOutTracklet->fPterr = h*h*tp.GetSigma1Pt2(); h = 1./TMath::Sqrt(1-snp*snp); currOutTracklet->fPsierr = h*h*tp.GetSigmaSnp2(); currOutTracklet->fTglerr = tp.GetSigmaTgl2(); - - par.TransportToX( vHitStoreX[iLastHit] ); - currOutTracklet->fLastX = par.GetX(); - currOutTracklet->fLastY = par.GetY(); - currOutTracklet->fLastZ = par.GetZ(); - + + if( par.TransportToX( vHitStoreX[iLastHit] ) ){ + currOutTracklet->fLastX = par.GetX(); + currOutTracklet->fLastY = par.GetY(); + currOutTracklet->fLastZ = par.GetZ(); + } else { + currOutTracklet->fLastX = vHitStoreX[iLastHit]; + currOutTracklet->fLastY = vHitStoreY[iLastHit]; + currOutTracklet->fLastZ = vHitStoreZ[iLastHit]; + } + //if( currOutTracklet->fLastX<10. ) { + //HLTError("CA last point: hitxyz=%f,%f,%f, track=%f,%f,%f, tracklet=%f,%f,%f, nhits=%d",vHitStoreX[iLastHit],vHitStoreY[iLastHit],vHitStoreZ[iLastHit], + //par.GetX(), par.GetY(),par.GetZ(),currOutTracklet->fLastX,currOutTracklet->fLastY ,currOutTracklet->fLastZ, t.NHits()); + //} #ifdef INCLUDE_TPC_HOUGH #ifdef ROWHOUGHPARAMS currOutTracklet->fTrackID = 0; @@ -537,10 +658,12 @@ int AliHLTTPCCATrackerComponent::DoEvent outPtr->fTrackletCnt++; } - delete[] vHits; - delete[] vHitStoreX; - delete[] vHitStoreID; - delete[] vHitRowID; + if( vHits ) delete[] vHits; + if( vHitStoreX ) delete[] vHitStoreX; + if( vHitStoreY ) delete[] vHitStoreY; + if( vHitStoreZ ) delete[] vHitStoreZ; + if( vHitStoreID ) delete[] vHitStoreID; + if( vHitRowID ) delete[] vHitRowID; AliHLTComponentBlockData bd; FillBlockData( bd ); @@ -553,15 +676,16 @@ int AliHLTTPCCATrackerComponent::DoEvent timer.Stop(); - fFullTime+= timer.CpuTime(); - fRecoTime+= timerReco.CpuTime(); + fFullTime+= timer.RealTime(); + fRecoTime+= timerReco.RealTime(); fNEvents++; // Set log level to "Warning" for on-line system monitoring - + Int_t hz = (Int_t) (fFullTime>1.e-10 ?fNEvents/fFullTime :100000); + Int_t hz1 = (Int_t) (fRecoTime>1.e-10 ?fNEvents/fRecoTime :100000); Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Tracks", - "CATracker slice %d: output %d tracks; input %d clusters, patches %d..%d, rows %d..%d; reco time %d/%d us", - slice, ntracks, nClusters, minPatch, maxPatch, row[0], row[1], (Int_t) (fFullTime/fNEvents*1.e6), (Int_t) (fRecoTime/fNEvents*1.e6) ); + "CATracker slice %d: output %d tracks; input %d clusters, patches %d..%d, rows %d..%d; reco time %d/%d Hz", + slice, ntracks, nClusters, minPatch, maxPatch, row[0], row[1], hz, hz1 ); return ret; diff --git a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerComponent.h b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerComponent.h index 83883edde81..eeb1d729405 100644 --- a/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerComponent.h +++ b/HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerComponent.h @@ -11,6 +11,7 @@ #include "AliHLTProcessor.h" class AliHLTTPCCATracker; +class AliHLTTPCSpacePointData; /** * @class AliHLTTPCCATrackerComponent @@ -72,12 +73,16 @@ private: AliHLTTPCCATracker* fTracker; //! transient /** magnetic field */ - Double_t fBField; // see above + Double_t fSolenoidBz; // see above Int_t fMinNTrackClusters; //* required min number of clusters on the track + Double_t fCellConnectionAngleXY; //* max phi angle between connected cells (deg) + Double_t fCellConnectionAngleXZ; //* max psi angle between connected cells (deg) + Double_t fClusterZCut; //* cut on cluster Z position (for noise rejection at the age of TPC) + Double_t fFullTime; //* total time for DoEvent() [s] + Double_t fRecoTime; //* total reconstruction time [s] + Long_t fNEvents; //* number of reconstructed events - Double_t fFullTime; //! total time for DoEvent() [s] - Double_t fRecoTime; //! total reconstruction time [s] - Long_t fNEvents; //! number of reconstructed events + static Bool_t CompareClusters(AliHLTTPCSpacePointData *a, AliHLTTPCSpacePointData *b); ClassDef(AliHLTTPCCATrackerComponent, 0); diff --git a/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx b/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx index b31be58473c..d92ff85c3b4 100644 --- a/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx +++ b/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx @@ -35,6 +35,7 @@ #include "AliHLTTPCCAOutTrack.h" #include "AliHLTTPCCAPerformance.h" #include "AliHLTTPCCAParam.h" +#include "AliHLTTPCCATrackConvertor.h" #include "TMath.h" #include "AliTPCLoader.h" @@ -46,18 +47,20 @@ #include "AliTPCtrack.h" #include "AliESDtrack.h" #include "AliESDEvent.h" +#include "AliTrackReference.h" +#include ClassImp(AliTPCtrackerCA) AliTPCtrackerCA::AliTPCtrackerCA() - :AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0) + :AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0) { //* default constructor } AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCtrackerCA &): - AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0) + AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0) { //* dummy } @@ -78,11 +81,12 @@ AliTPCtrackerCA::~AliTPCtrackerCA() } AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par): - AliTracker(),fParam(par), fClusters(0), fNClusters(0), fHLTTracker(0), fHLTPerformance(0),fDoHLTPerformance(0) + AliTracker(),fParam(par), fClusters(0), fNClusters(0), fHLTTracker(0), fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0) { //* constructor - DoHLTPerformance() = 0; + DoHLTPerformance() = 1; + DoHLTPerformanceClusters() = 1; fHLTTracker = new AliHLTTPCCAGBTracker; fHLTTracker->SetNSlices( fParam->GetNSector()/2 ); @@ -94,32 +98,32 @@ AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par): for( int iSlice=0; iSliceNSlices(); iSlice++ ){ - Double_t bz = AliTracker::GetBz(); - - Double_t inRmin = fParam->GetInnerRadiusLow(); - //Double_t inRmax = fParam->GetInnerRadiusUp(); - //Double_t outRmin = fParam->GetOuterRadiusLow(); - Double_t outRmax = fParam->GetOuterRadiusUp(); - Double_t plusZmin = 0.0529937; - Double_t plusZmax = 249.778; - Double_t minusZmin = -249.645; - Double_t minusZmax = -0.0799937; - Double_t dalpha = 0.349066; - Double_t alpha = 0.174533 + dalpha*iSlice; + Float_t bz = AliTracker::GetBz(); + + Float_t inRmin = fParam->GetInnerRadiusLow(); + //Float_t inRmax = fParam->GetInnerRadiusUp(); + //Float_t outRmin = fParam->GetOuterRadiusLow(); + Float_t outRmax = fParam->GetOuterRadiusUp(); + Float_t plusZmin = 0.0529937; + Float_t plusZmax = 249.778; + Float_t minusZmin = -249.645; + Float_t minusZmax = -0.0799937; + Float_t dalpha = 0.349066; + Float_t alpha = 0.174533 + dalpha*iSlice; Bool_t zPlus = (iSlice<18 ); - Double_t zMin = zPlus ?plusZmin :minusZmin; - Double_t zMax = zPlus ?plusZmax :minusZmax; + Float_t zMin = zPlus ?plusZmin :minusZmin; + Float_t zMax = zPlus ?plusZmax :minusZmax; //TPCZmin = -249.645, ZMax = 249.778 - //Double_t rMin = inRmin; - //Double_t rMax = outRmax; + //Float_t rMin = inRmin; + //Float_t rMax = outRmax; - Double_t padPitch = 0.4; - Double_t sigmaZ = 0.228808; + Float_t padPitch = 0.4; + Float_t sigmaZ = 0.228808; Int_t NRows = fParam->GetNRowLow()+fParam->GetNRowUp(); - Double_t rowX[200]; + Float_t rowX[200]; for( Int_t irow=0; irowGetNRowLow(); irow++){ rowX[irow] = fParam->GetPadRowRadiiLow(irow); } @@ -129,10 +133,10 @@ AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par): AliHLTTPCCAParam param; param.Initialize( iSlice, NRows, rowX, alpha, dalpha, inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, bz ); - param.YErrorCorrection() = .33;//1; - param.ZErrorCorrection() = .33;//2; + param.YErrorCorrection() = 1.;//.33; + param.ZErrorCorrection() = 1.;//.33; param.MaxTrackMatchDRow() = 5; - param.TrackConnectionFactor() = 5.; + param.TrackConnectionFactor() = 3.5; fHLTTracker->Slices()[iSlice].Initialize( param ); } } @@ -150,13 +154,13 @@ Int_t AliTPCtrackerCA::LoadClusters (TTree * tree) if( !fParam ) return 1; // load mc tracks - if( fDoHLTPerformance ){ - if( !gAlice ) return 0; + while( fDoHLTPerformance ){ + if( !gAlice ) break; AliRunLoader *rl = gAlice->GetRunLoader(); - if( !rl ) return 0; + if( !rl ) break; rl->LoadKinematics(); AliStack *stack = rl->Stack(); - if( !stack ) return 0 ; + if( !stack ) break; fHLTPerformance->SetNMCTracks( stack->GetNtrack() ); @@ -164,7 +168,73 @@ Int_t AliTPCtrackerCA::LoadClusters (TTree * tree) TParticle *part = stack->Particle(itr); fHLTPerformance->ReadMCTrack( itr, part ); } - } + + { // check for MC tracks at the TPC entrance + Bool_t *isTPC = 0; + isTPC = new Bool_t [stack->GetNtrack()]; + for( Int_t i=0; iGetNtrack(); i++ ) isTPC[i] = 0; + rl->LoadTrackRefs(); + TTree *TR = rl->TreeTR(); + if( !TR ) break; + TBranch *branch=TR->GetBranch("TrackReferences"); + if (!branch ) break; + TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs=&tpcdummy; + branch->SetAddress(&tpcRefs); + Int_t nr=(Int_t)TR->GetEntries(); + for (Int_t r=0; rGetEvent(r); + Int_t nref = tpcRefs->GetEntriesFast(); + if (!nref) continue; + AliTrackReference *tpcRef= 0x0; + for (Int_t iref=0; irefUncheckedAt(iref); + if (tpcRef->DetectorId() == AliTrackReference::kTPC) break; + tpcRef = 0x0; + } + if (!tpcRef) continue; + + if( isTPC[tpcRef->Label()] ) continue; + + fHLTPerformance->ReadMCTPCTrack(tpcRef->Label(), + tpcRef->X(),tpcRef->Y(),tpcRef->Z(), + tpcRef->Px(),tpcRef->Py(),tpcRef->Pz() ); + isTPC[tpcRef->Label()] = 1; + tpcRefs->Clear(); + } + if( isTPC ) delete[] isTPC; + } + + while( fDoHLTPerformanceClusters ){ + AliTPCLoader *tpcl = (AliTPCLoader*) rl->GetDetectorLoader("TPC"); + if( !tpcl ) break; + if( tpcl->TreeH() == 0x0 ){ + if( tpcl->LoadHits() ) break; + } + if( tpcl->TreeH() == 0x0 ) break; + + AliTPC *tpc = (AliTPC*) gAlice->GetDetector("TPC"); + Int_t nEnt=(Int_t)tpcl->TreeH()->GetEntries(); + Int_t nPoints = 0; + for (Int_t iEnt=0; iEntResetHits(); + tpcl->TreeH()->GetEvent(iEnt); + AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1); + for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ) nPoints++; + } + fHLTPerformance->SetNMCPoints( nPoints ); + + for (Int_t iEnt=0; iEntResetHits(); + tpcl->TreeH()->GetEvent(iEnt); + AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1); + for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ){ + fHLTPerformance->ReadMCPoint( phit->GetTrack(),phit->X(), phit->Y(),phit->Z(),phit->Time(), phit->fSector%36); + } + } + break; + } + break; + } TBranch * br = tree->GetBranch("Segment"); if( !br ) return 1; @@ -196,7 +266,7 @@ Int_t AliTPCtrackerCA::LoadClusters (TTree * tree) Int_t sec,row; fParam->AdjustSectorRow(clrow->GetID(),sec,row); int NClu = clrow->GetArray()->GetEntriesFast(); - Double_t x = fParam->GetPadRowRadii(sec,row); + Float_t x = fParam->GetPadRowRadii(sec,row); for (Int_t icl=0; iclSetY(posC[1]); cluster->SetZ(posC[2]); - Double_t y = cluster->GetY(); - Double_t z = cluster->GetZ(); + Float_t y = cluster->GetY(); + Float_t z = cluster->GetZ(); if( sec>=36 ){ sec = sec - 36; @@ -246,8 +316,8 @@ Int_t AliTPCtrackerCA::LoadClusters (TTree * tree) Int_t index = ind++; fClusters[index] = *cluster; fHLTTracker->ReadHit( x, y, z, - TMath::Sqrt(cluster->GetSigmaY2()), TMath::Sqrt(cluster->GetSigmaZ2()), - index, sec, row ); + TMath::Sqrt(cluster->GetSigmaY2()), TMath::Sqrt(cluster->GetSigmaZ2()), + cluster->GetMax(), index, sec, row ); if( fDoHLTPerformance ) fHLTPerformance->ReadHitLabel(index, lab0, lab1, lab2 ); } } @@ -267,34 +337,80 @@ Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event ) fHLTTracker->FindTracks(); if( fDoHLTPerformance ) fHLTPerformance->Performance(); + if( 0 ) {// Write Event + if( fStatNEvents == 0 ){ + fstream geo; + geo.open("CAEvents/settings.dat", ios::out); + if( geo.is_open() ){ + fHLTTracker->WriteSettings(geo); + } + geo.close(); + } + + fstream event; + char name[255]; + sprintf( name,"CAEvents/%i.event.dat",fStatNEvents ); + event.open(name, ios::out); + if( event.is_open() ){ + fHLTTracker->WriteEvent(event); + fstream tracks; + sprintf( name,"CAEvents/%i.tracks.dat",fStatNEvents ); + tracks.open(name, ios::out); + fHLTTracker->WriteTracks(tracks); + } + event.close(); + if( fDoHLTPerformance ){ + fstream mcevent, mcpoints; + char mcname[255]; + sprintf( mcname,"CAEvents/%i.mcevent.dat",fStatNEvents ); + mcevent.open(mcname, ios::out); + if( mcevent.is_open() ){ + fHLTPerformance->WriteMCEvent(mcevent); + } + if(fDoHLTPerformanceClusters ){ + sprintf( mcname,"CAEvents/%i.mcpoints.dat",fStatNEvents ); + mcpoints.open(mcname, ios::out); + if( mcpoints.is_open() ){ + fHLTPerformance->WriteMCPoints(mcpoints); + } + mcpoints.close(); + } + mcevent.close(); + } + } + fStatNEvents++; + if( event ){ + Float_t bz = fHLTTracker->Slices()[0].Param().Bz(); + for( Int_t itr=0; itrNTracks(); itr++ ){ AliTPCtrack tTPC; AliHLTTPCCAGBTrack &tCA = fHLTTracker->Tracks()[itr]; AliHLTTPCCATrackParam &par = tCA.Param(); - - par.GetExtParam( tTPC, tCA.Alpha(), fHLTTracker->Slices()[0].Param().Bz() ); - + AliHLTTPCCATrackConvertor::GetExtParam( par, tTPC, tCA.Alpha(), bz ); tTPC.SetMass(0.13957); + tTPC.SetdEdx( tCA.DeDx() ); + if( TMath::Abs(tTPC.GetSigned1Pt())>1./0.1 ) continue; int nhits = tCA.NHits(); if( nhits>kMaxRow ) nhits = kMaxRow; - tTPC.SetNumberOfClusters(nhits); + tTPC.SetNumberOfClusters(nhits); for( Int_t ih=0; ihTrackHits()[tCA.FirstHitRef()+ih]; Int_t ext_index = fHLTTracker->Hits()[index].ID(); tTPC.SetClusterIndex(ih, ext_index); } - CookLabel(&tTPC,0.1); + CookLabel(&tTPC,0.1); { - Double_t xTPC=83.65; - if (tTPC.AliExternalTrackParam::PropagateTo(xTPC,5)) { + Double_t xTPC=fParam->GetInnerRadiusLow(); + Double_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi(); + if (tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz)) { Double_t y=tTPC.GetY(); - Double_t ymax=xTPC*TMath::Tan(1.74532920122146606e-01); + Double_t ymax=xTPC*TMath::Tan(dAlpha/2.); if (y > ymax) { - if (tTPC.Rotate(2*1.74532920122146606e-01)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,5); + if (tTPC.Rotate(dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz); } else if (y <-ymax) { - if (tTPC.Rotate(-2*1.74532920122146606e-01)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,5); + if (tTPC.Rotate(-dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz); } } } @@ -303,7 +419,7 @@ Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event ) tESD.UpdateTrackParams( &(tTPC),AliESDtrack::kTPCin); //tESD.SetStatus( AliESDtrack::kTPCrefit ); //tESD.SetTPCPoints(tTPC.GetPoints()); - //tESD.myTPC = tTPC; + //tESD.myTPC = tTPC; event->AddTrack(&tESD); } } @@ -313,10 +429,42 @@ Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event ) } -Int_t AliTPCtrackerCA::RefitInward (AliESDEvent *) +Int_t AliTPCtrackerCA::RefitInward (AliESDEvent *event) { - //* not implemented yet - return 0; + //* back propagation of ESD tracks (not fully functional) + + Float_t bz = fHLTTracker->Slices()[0].Param().Bz(); + Float_t xTPC = fParam->GetInnerRadiusLow(); + Float_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi(); + Float_t yMax = xTPC*TMath::Tan(dAlpha/2.); + + Int_t nentr=event->GetNumberOfTracks(); + + for (Int_t i=0; iGetTrack(i); + ULong_t status=esd->GetStatus(); + if (!(status&AliESDtrack::kTPCin)) continue; + AliHLTTPCCATrackParam t0; + AliHLTTPCCATrackConvertor::SetExtParam(t0,*esd, bz ); + Float_t alpha = esd->GetAlpha(); + if( t0.TransportToXWithMaterial( xTPC, bz) ){ + if (t0.GetY() > yMax) { + if (t0.Rotate(dAlpha)){ + alpha+=dAlpha; + t0.TransportToXWithMaterial( xTPC, bz); + } + } else if (t0.GetY() <-yMax) { + if (t0.Rotate(-dAlpha)){ + alpha+=-dAlpha; + t0.TransportToXWithMaterial( xTPC, bz); + } + } + } + AliTPCtrack tt(*esd); + AliHLTTPCCATrackConvertor::GetExtParam(t0,tt,alpha,bz); + esd->UpdateTrackParams( &tt,AliESDtrack::kTPCrefit); + } + return 0; } Int_t AliTPCtrackerCA::PropagateBack(AliESDEvent *) diff --git a/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.h b/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.h index 05ad6dfb6b9..79886cb4079 100644 --- a/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.h +++ b/HLT/TPCLib/tracking-ca/AliTPCtrackerCA.h @@ -43,6 +43,7 @@ public: void UnloadClusters(){ return ; } AliCluster * GetCluster(Int_t index) const; Bool_t &DoHLTPerformance(){ return fDoHLTPerformance; } + Bool_t &DoHLTPerformanceClusters(){ return fDoHLTPerformanceClusters; } // protected: @@ -52,6 +53,8 @@ public: AliHLTTPCCAGBTracker *fHLTTracker; //* pointer to the HLT tracker AliHLTTPCCAPerformance *fHLTPerformance; //* performance calculations Bool_t fDoHLTPerformance; //* flag for call AliHLTTPCCAPerformance + Bool_t fDoHLTPerformanceClusters; //* flag for call AliHLTTPCCAPerformance with cluster pulls (takes some time to load TPC MC points) + Int_t fStatNEvents; //* N of reconstructed events ClassDef(AliTPCtrackerCA,1) }; diff --git a/HLT/libAliHLTTPC.pkg b/HLT/libAliHLTTPC.pkg index 8f3d5f7455a..ed6be7f2326 100644 --- a/HLT/libAliHLTTPC.pkg +++ b/HLT/libAliHLTTPC.pkg @@ -57,11 +57,13 @@ CLASS_HDRS:= AliHLTTPCTransform.h \ tracking-ca/AliHLTTPCCAHit.h \ tracking-ca/AliHLTTPCCAOutTrack.h \ tracking-ca/AliHLTTPCCATrackParam.h \ + tracking-ca/AliHLTTPCCATrackConvertor.h \ tracking-ca/AliHLTTPCCAParam.h \ tracking-ca/AliHLTTPCCARow.h \ tracking-ca/AliHLTTPCCAGrid.h \ tracking-ca/AliHLTTPCCAGBHit.h \ tracking-ca/AliHLTTPCCAMCTrack.h \ + tracking-ca/AliHLTTPCCAMCPoint.h \ tracking-ca/AliHLTTPCCAGBTrack.h \ tracking-ca/AliHLTTPCCAGBTracker.h \ tracking-ca/AliHLTTPCCAPerformance.h \ -- 2.43.0