diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDBasicKalman.h b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDBasicKalman.h
index 23dcdfbcd1eb18f7f556cfe6384325b17ce929ce..95ad0cb404d969773eb1930a8d77e2068ebe9e74 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDBasicKalman.h
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDBasicKalman.h
@@ -28,8 +28,8 @@
#define SIDSTATIONID 4
-#define MAXCHI2HIT 2.0
-#define MAXCHI2TRK 2.0
+#define MAXCHI2HIT 3.0
+#define MAXCHI2TRK 3.0
#define MAXSHAREDHITS 2
using namespace std;
diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDLocReco.h b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDLocReco.h
index 65b72e9d3919827f759f1b88c2970e8dd8260bed..af4386b6ad41afa11b5a3c4833c102e1f529691a 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDLocReco.h
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/AFP_LocReco/AFP_SIDLocReco.h
@@ -2,125 +2,127 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
- #ifndef AFP_SIDLOCRECO_h
-#define AFP_SIDLOCRECO_h
-
-#include <iostream>
-#include <string>
-#include <list>
-#include <map>
-#include <vector>
-#include <fstream>
-
-//#include "GaudiKernel/Algorithm.h"
-#include "AthenaBaseComps/AthAlgorithm.h"
-#include "GaudiKernel/MsgStream.h"
-#include "GaudiKernel/ObjectVector.h"
-#include "GaudiKernel/ServiceHandle.h"
-#include "GaudiKernel/IToolSvc.h"
-
-#include "AthenaKernel/getMessageSvc.h"
-#include "AthenaKernel/IAtRndmGenSvc.h"
-
-#include "StoreGate/StoreGateSvc.h"
-#include "StoreGate/DataHandle.h"
-
-#include "EventInfo/EventInfo.h"
-#include "EventInfo/EventID.h"
-
-#include "AthenaPoolUtilities/AthenaAttributeList.h"
-#include "AthenaPoolUtilities/CondAttrListCollection.h"
-
-//#include "AFP_RawEv/AFP_RawData.h"
-//#include "AFP_RawEv/AFP_RawDataContainer.h"
-//#include "AFP_RawEv/AFP_RawDataCollection.h"
-//#include "AFP_RawEv/AFP_SIDDigitCollection.h"
-//#include "AFP_RawEv/AFP_TDDigitCollection.h"
-
-#include "AFP_Geometry/AFP_constants.h"
-#include "AFP_Geometry/AFP_Geometry.h"
-#include "AFP_Geometry/AFP_ConfigParams.h"
-#include "CLHEP/Geometry/Point3D.h"
-
-#include "AFP_DigiEv/AFP_SiDigiCollection.h"
-#include "AFP_LocRecoEv/AFP_SIDLocRecoEvCollection.h"
-
-//#include "AFP_LocRecoEv/AFP_SIDLocRecoEvCollection.h"
-
-#include "AFP_LocReco/AFP_UserObjects.h"
-#include "AFP_LocReco/AFP_SIDBasicKalman.h"
-
-#include "TROOT.h"
-
-//for truth particles
-#include "GeneratorObjects/McEventCollection.h"
-#include "HepMC/GenEvent.h"
-#include "HepMC/GenVertex.h"
-#include "HepMC/GenParticle.h"
-
-
-
-using namespace std;
-
-#define SIDSTATIONID 4
-
-class StoreGateSvc;
-class ActiveStoreSvc;
-
-class AFP_SIDLocReco : public AthAlgorithm
-{
- public:
- AFP_SIDLocReco(const string& name, ISvcLocator* pSvcLocator);
- ~AFP_SIDLocReco();
-
- private:
- AFP_CONFIGURATION m_Config;
- AFP_Geometry* m_pGeometry;
-
- // a handle on Store Gate
- StoreGateSvc* m_storeGate;
- //StoreGateSvc* m_pDetStore;
-
- AFP_SIDLocRecoEvCollection* m_pSIDLocRecoEvCollection;
- AFP_SIDLocRecoEvent* m_pSIDLocRecoEvent;
-
- private:
-
- UInt_t m_eventNum; //real event number
- Int_t m_iRunNum;
- Int_t m_iDataType; //data type (simulation or real data) using in the local reconstruction
- Int_t m_iEvent; //event number from zero value
- Float_t m_AmpThresh; // TD signal amplitude threshold
-
- //slope and X,Y,Z-pos for SID plates [4][6]
- Float_t m_fsSID[SIDSTATIONID][SIDCNT];
- Float_t m_fxSID[SIDSTATIONID][SIDCNT];
- Float_t m_fySID[SIDSTATIONID][SIDCNT];
- Float_t m_fzSID[SIDSTATIONID][SIDCNT];
-
-
- string m_strKeyGeometryForReco;
- vector<string> m_vecListAlgoSID;
- string m_strAlgoSID;
-
- string m_strKeySIDLocRecoEvCollection;
- string m_strSIDCollectionName;
-
- public:
- StatusCode initialize();
- StatusCode execute();
- StatusCode finalize();
-
- private:
- bool ReadGeometryDetCS();
- bool StoreReconstructionGeometry(/*const char* szDataDestination*/);
- void SaveGeometry();
- void ClearGeometry();
-
- StatusCode AFPCollectionReading(list<SIDHIT> &ListSIDHits);
-
- StatusCode RecordSIDCollection();
- StatusCode ExecuteRecoMethod(const string strAlgo, const list<SIDHIT> &ListSIDHits);
-};
-
-#endif //AFP_TDLOCRECO_h
+#ifndef AFP_SIDLOCRECO_h
+#define AFP_SIDLOCRECO_h
+
+#include <iostream>
+#include <string>
+#include <list>
+#include <map>
+#include <vector>
+#include <fstream>
+
+//#include "GaudiKernel/Algorithm.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/ObjectVector.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IToolSvc.h"
+
+#include "AthenaKernel/getMessageSvc.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include "StoreGate/StoreGateSvc.h"
+#include "StoreGate/DataHandle.h"
+
+#include "EventInfo/EventInfo.h"
+#include "EventInfo/EventID.h"
+
+#include "AthenaPoolUtilities/AthenaAttributeList.h"
+#include "AthenaPoolUtilities/CondAttrListCollection.h"
+
+//#include "AFP_RawEv/AFP_RawData.h"
+//#include "AFP_RawEv/AFP_RawDataContainer.h"
+//#include "AFP_RawEv/AFP_RawDataCollection.h"
+//#include "AFP_RawEv/AFP_SIDDigitCollection.h"
+//#include "AFP_RawEv/AFP_TDDigitCollection.h"
+
+#include "AFP_Geometry/AFP_constants.h"
+#include "AFP_Geometry/AFP_Geometry.h"
+#include "AFP_Geometry/AFP_ConfigParams.h"
+#include "CLHEP/Geometry/Point3D.h"
+
+#include "AFP_DigiEv/AFP_SiDigiCollection.h"
+#include "AFP_LocRecoEv/AFP_SIDLocRecoEvCollection.h"
+
+//#include "AFP_LocRecoEv/AFP_SIDLocRecoEvCollection.h"
+
+#include "AFP_LocReco/AFP_UserObjects.h"
+#include "AFP_LocReco/AFP_SIDBasicKalman.h"
+
+#include "TROOT.h"
+
+//for truth particles
+#include "GeneratorObjects/McEventCollection.h"
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenVertex.h"
+#include "HepMC/GenParticle.h"
+
+// xAOD
+#include "xAODForward/AFPTrackContainer.h"
+
+
+using namespace std;
+
+#define SIDSTATIONID 4
+
+class StoreGateSvc;
+class ActiveStoreSvc;
+
+class AFP_SIDLocReco : public AthAlgorithm
+{
+ public:
+ AFP_SIDLocReco(const string& name, ISvcLocator* pSvcLocator);
+ ~AFP_SIDLocReco();
+
+ private:
+ AFP_CONFIGURATION m_Config;
+ AFP_Geometry* m_pGeometry;
+
+ // a handle on Store Gate
+ StoreGateSvc* m_storeGate;
+ //StoreGateSvc* m_pDetStore;
+
+ AFP_SIDLocRecoEvCollection* m_pSIDLocRecoEvCollection;
+ AFP_SIDLocRecoEvent* m_pSIDLocRecoEvent;
+
+ private:
+
+ UInt_t m_eventNum; //real event number
+ Int_t m_iRunNum;
+ Int_t m_iDataType; //data type (simulation or real data) using in the local reconstruction
+ Int_t m_iEvent; //event number from zero value
+ Float_t m_AmpThresh; // TD signal amplitude threshold
+
+ //slope and X,Y,Z-pos for SID plates [4][6]
+ Float_t m_fsSID[SIDSTATIONID][SIDCNT];
+ Float_t m_fxSID[SIDSTATIONID][SIDCNT];
+ Float_t m_fySID[SIDSTATIONID][SIDCNT];
+ Float_t m_fzSID[SIDSTATIONID][SIDCNT];
+
+
+ string m_strKeyGeometryForReco;
+ vector<string> m_vecListAlgoSID;
+ string m_strAlgoSID;
+
+ string m_strKeySIDLocRecoEvCollection;
+ string m_strSIDCollectionName;
+
+ public:
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+ private:
+ bool ReadGeometryDetCS();
+ bool StoreReconstructionGeometry(/*const char* szDataDestination*/);
+ void SaveGeometry();
+ void ClearGeometry();
+
+ StatusCode AFPCollectionReading(list<SIDHIT> &ListSIDHits);
+
+ StatusCode RecordSIDCollection();
+ StatusCode ExecuteRecoMethod(const string strAlgo, const list<SIDHIT> &ListSIDHits, xAOD::AFPTrackContainer* resultContainer);
+};
+
+#endif //AFP_TDLOCRECO_h
diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/CMakeLists.txt b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/CMakeLists.txt
index cf848f76af42ec7738afbfdbbb95eeb70b94dc7a..60e6fe773698f52a30d0a7e52680670a042dc40b 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/CMakeLists.txt
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/CMakeLists.txt
@@ -16,7 +16,10 @@ atlas_depends_on_subdirs( PUBLIC
ForwardDetectors/AFP/AFP_Geometry
ForwardDetectors/AFP/AFP_RecoEv/AFP_LocRecoEv
GaudiKernel
- Generators/GeneratorObjects )
+ Generators/GeneratorObjects
+ Event/xAOD/xAODForward
+ PRIVATE
+ Control/AthLinks)
# External dependencies:
find_package( CLHEP )
@@ -32,7 +35,7 @@ atlas_add_component( AFP_LocReco
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaBaseComps AthenaKernel StoreGateLib SGtests AthenaPoolUtilities EventInfo AFP_DigiEv AFP_Geometry AFP_LocRecoEv GaudiKernel GeneratorObjects )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaBaseComps AthenaKernel StoreGateLib SGtests AthenaPoolUtilities EventInfo AFP_DigiEv AFP_Geometry AFP_LocRecoEv GaudiKernel GeneratorObjects xAODForward AthLinks)
# Install files from the package:
atlas_install_headers( AFP_LocReco )
diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/cmt/requirements b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/cmt/requirements
index d9ec1ab3d73fd2c2a7296766ec7a70ce2228e697..d359b76043712a659e6d60fbfbf4c6d1e56aef19 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/cmt/requirements
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/cmt/requirements
@@ -20,6 +20,8 @@ use AFP_Geometry AFP_Geometry-* ForwardDetectors/AFP
#use AFP_RawEv AFP_RawEv-* ForwardDetectors/AFP
use AFP_LocRecoEv AFP_LocRecoEv-* ForwardDetectors/AFP/AFP_RecoEv
use AFP_DigiEv AFP_DigiEv-* ForwardDetectors/AFP
+use AFP_EventTPCnv AFP_EventTPCnv-* ForwardDetectors/AFP/AFP_EventCnv
+use xAODForward xAODForward-* Event/xAOD
library AFP_LocReco *.cxx components/*.cxx
@@ -35,9 +37,8 @@ apply_pattern declare_runtime files="-s=../share/mapping *.txt *.dat *.py"
apply_pattern declare_joboptions files="*.py"
private
-
#use PathResolver PathResolver-* Tools
-
+use AthLinks AthLinks-* Control
#macro cppdebugflags '$(cppdebugflags_s)'
#macro_remove componentshr_linkopts "-Wl,-s"
diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDBasicKalman.cxx b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDBasicKalman.cxx
index 161fcb89770e946f203695f0c25eda17cad7104f..a7727ead66c00658c3821dcdaeb106794fa21989 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDBasicKalman.cxx
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDBasicKalman.cxx
@@ -2,718 +2,812 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
-#include "AFP_LocReco/AFP_SIDBasicKalman.h"
-
-AFP_SIDBasicKalman::AFP_SIDBasicKalman()
-{
- MyFile = NULL;
- histS1_PixMap = NULL;
- histS2_PixMap = NULL;
- histS3_PixMap = NULL;
- histS4_PixMap = NULL;
-
- m_listResults.clear();
-
- m_AmpThresh = 0;
- m_iDataType = 0;
-
- pTrkSeeds.clear();
-
- Fk.clear();
- Ck.clear();
- CkP.clear();
- Rk.clear();
-
- xk.clear();
- xkP.clear();
- rk.clear();
- rkP.clear();
- mk.clear();
- chik.clear();
-
- zk.clear();
- HID.clear();
-
- xkS.clear();
- CkS.clear();
- chikS.clear();
-
- z0 = 0;
-}
-
-AFP_SIDBasicKalman::~AFP_SIDBasicKalman()
-{
- m_listResults.clear();
-}
-
-StatusCode AFP_SIDBasicKalman::Initialize(float fAmpThresh, int iDataType, const list<SIDHIT> &ListSIDHits, Float_t fsSID[SIDSTATIONID][SIDCNT], Float_t fxSID[SIDSTATIONID][SIDCNT], Float_t fySID[SIDSTATIONID][SIDCNT], Float_t fzSID[SIDSTATIONID][SIDCNT])
-{
- m_AmpThresh = (float)fAmpThresh;
- m_iDataType = iDataType;
- m_ListSIDHits = ListSIDHits;
-
- Float_t delta_pixel_x = 0.050*CLHEP::mm; // size of pixel along x, in mm
- Float_t delta_pixel_y = 0.250*CLHEP::mm; // size of pixel along y, in mm
-
-
- // x-y-z map of all pixels
- for(Int_t nStationID = 0; nStationID < SIDSTATIONID; nStationID++)
- {
- for (Int_t nPlateID = 0; nPlateID < SIDCNT; nPlateID++)
- {
- for (Int_t nPixel_x = 0; nPixel_x < 336; nPixel_x++)
- {
- for (Int_t nPixel_y = 0; nPixel_y < 80; nPixel_y++)
- {
- fxMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fxSID[nStationID][nPlateID]+(delta_pixel_x*nPixel_x)*cos(fsSID[nStationID][nPlateID]); //sign changed!
- fyMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fySID[nStationID][nPlateID]+(delta_pixel_y*nPixel_y); //sign changed!
- fzMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fzSID[nStationID][nPlateID] - ((nStationID<2)?1.:-1.)*(delta_pixel_x*nPixel_x)*sin(fsSID[nStationID][nPlateID]); //sign changed!
- }
- }
- }
- }
-
- // x-y-z map of hit pixels
- list<SIDHIT>::const_iterator iter;
- for (iter=ListSIDHits.begin(); iter!=ListSIDHits.end(); iter++)
- {
- if ((*iter).fADC > m_AmpThresh)
- {
- FillSIDHITPOS(*iter, m_MapSIDHitPos);
-
- /*SIDHITPOS SIDHitPos;
- Int_t nStID = (*iter).nStationID;
- Int_t nPlID = (*iter).nDetectorID;
- Int_t nPixCol = (*iter).nPixelCol;
- Int_t nPixRow = (*iter).nPixelRow;
-
-
- SIDHitPos.nStationID = nStID;
- SIDHitPos.nPlateID = nPlID;
- SIDHitPos.fPixX = fxSID[nStID][nPlID] - (delta_pixel_x*nPixCol+delta_pixel_x/2.)*cos(fsSID[nStID][nPlID]);
- SIDHitPos.fPixY = fySID[nStID][nPlID] - (delta_pixel_y*nPixRow+delta_pixel_y/2.);
- SIDHitPos.fPixZ = ((nStID<2)?1.:-1.)*fzSID[nStID][nPlID] + ((nStID<2)?1.:-1.)*(delta_pixel_x*nPixCol+delta_pixel_x/2.)*sin(fsSID[nStID][nPlID]);
-
- Int_t HitID = nStID*1000000+nPlID*100000+nPixCol*100+nPixRow;
-
- if(m_MapSIDHitPos.find(HitID-100)==m_MapSIDHitPos.end() && m_MapSIDHitPos.find(HitID+100)==m_MapSIDHitPos.end())
- {
- m_MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
- }*/
- }
- }
-
- // Kalman matrices initialization
- Hk = CLHEP::HepMatrix(2, 4, 0);
- Hk[0][0] = 1.;
- Hk[1][2] = 1.;
-
- Qk = CLHEP::HepMatrix(4, 4, 0);
-
- Vk = CLHEP::HepMatrix(2, 2, 0);
- Vk[0][0] = delta_pixel_x*delta_pixel_x;
- Vk[1][1] = delta_pixel_y*delta_pixel_y;
-
- C0 = CLHEP::HepMatrix(4, 4, 0);
- C0[0][0] = delta_pixel_x*delta_pixel_x/3.;
- C0[1][1] = delta_pixel_x*delta_pixel_x/NOMINALSIDSPACING/NOMINALSIDSPACING/3.;
- //C0[1][0] = 0;
- //C0[0][1] = 0;
-
- C0[2][2] = delta_pixel_y*delta_pixel_y/3.;
- C0[3][3] = delta_pixel_y*delta_pixel_y/NOMINALSIDSPACING/NOMINALSIDSPACING/3.;
- //C0[3][2] = 0;
- //C0[2][3] = 0;
-
- //HistInitialize();
-
- return StatusCode::SUCCESS;
-}
-
-StatusCode AFP_SIDBasicKalman::Execute()
-{
- MsgStream LogStream(Athena::getMessageSvc(), "AFP_SIDBasicKalman::Execute()");
-
- /*
- map<Int_t, SIDHITPOS>::const_iterator iter0;
- for (iter0=m_MapSIDHitPos.begin(); iter0!=m_MapSIDHitPos.end(); iter0++)
- {
- if( (*iter0).second.nStationID==1 )
- {
- histS1_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixX/CLHEP::mm);
- histS3_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixY/CLHEP::mm);
- }
- else if( (*iter0).second.nStationID==0 )
- {
- histS2_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixX/CLHEP::mm);
- histS4_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixY/CLHEP::mm);
- }
- }
- */
-
-
- ////////////////////////
- GetTrkSeeds();
- ////////////////////////
-
- vector<SIDHITSEED>::const_iterator iter1;
- for (iter1= pTrkSeeds.begin(); iter1!= pTrkSeeds.end(); iter1++)
- {
- ClearMatrix();
-
- Int_t LastPlate = (*iter1).nLastPlate;
-
- if ( FillTrkPropagators(*iter1, ++LastPlate) )
- {
- ++LastPlate;
- Int_t DoubleHole=0;
- for (Int_t i=LastPlate; i<SIDCNT; i++)
- {
- if ( FillTrkPropagators((*iter1).nStationID, i) )
- {
- if(DoubleHole==1) DoubleHole=0;
- }
- else
- {
- ++DoubleHole;
- }
- if(DoubleHole==2) break;
-
- }
- if(DoubleHole==2) continue;
- }
-
- else if ( FillTrkPropagators(*iter1, ++LastPlate) )
- {
- ++LastPlate;
- Int_t DoubleHole=0;
- for (Int_t i=LastPlate; i<SIDCNT; i++)
- {
- if ( FillTrkPropagators((*iter1).nStationID, i) )
- {
- if(DoubleHole==1) DoubleHole=0;
- }
- else
- {
- ++DoubleHole;
- }
- if(DoubleHole==2) break;
-
- }
- if(DoubleHole==2) continue;
- }
-
- else continue; // double hole after the seed pixels
-
- LogStream << MSG::DEBUG << "Found new track candidate with parameters:" << endreq;
- LogStream << MSG::DEBUG << "Hit ID's : ";
- vector<Int_t>::const_iterator iter2;
- for (iter2= HID.begin(); iter2!= HID.end(); iter2++) LogStream << MSG::DEBUG << (*iter2) << "\t";
- LogStream << endreq;
-
- LogStream << MSG::DEBUG << "Filtered parameters : X, DX, Y, DY, chi2" << endreq;
- vector< CLHEP::HepVector >::const_iterator iter3;
- vector< Float_t >::const_iterator iter4 = chik.begin();
- LogStream << MSG::DEBUG << (*iter1).fSeedX << "\t" << (*iter1).fSeedDX << "\t" << (*iter1).fSeedY << "\t" << (*iter1).fSeedDY << "\t" << 0 << endreq;
- for (iter3=xk.begin(); iter3!= xk.end(); iter3++)
- {
- LogStream << MSG::DEBUG << (*iter3)[0] << "\t" << (*iter3)[1] << "\t" << (*iter3)[2] << "\t" << (*iter3)[3] << "\t" << (*iter4) << endreq;
- iter4++;
- }
-
-
- ////////////
- Smooth();
- ////////////
-
- Float_t Chi2Sum = 0;
-
- LogStream << MSG::DEBUG << "Smoothed parameters : X, DX, Y, DY, chi2" << endreq;
- vector< CLHEP::HepVector >::const_reverse_iterator iter5;
- vector< Float_t >::const_reverse_iterator iter6 = chikS.rbegin();
- for (iter5=xkS.rbegin(); iter5!= xkS.rend(); ++iter5)
- {
- LogStream << MSG::DEBUG << (*iter5)[0] << "\t" << (*iter5)[1] << "\t" << (*iter5)[2] << "\t" << (*iter5)[3] << "\t" << (*iter6) << endreq;
- Chi2Sum += (*iter6);
- ++iter6;
- }
- LogStream << endreq;
-
-
- //////////////////////////////////////////
- // fill tracking collection
- SIDRESULT Results;
- Results.nStationID = (*iter1).nStationID;
- Results.x_pos = xkS.back()[0];
- Results.y_pos = xkS.back()[2];
- Results.z_pos = z0;
- Results.x_slope = xkS.back()[1];
- Results.y_slope = xkS.back()[3];
- Results.z_slope = 1.;
- Results.nHits = xkS.size();
- Results.nHoles = SIDCNT - xkS.size();
- Results.fChi2 = xkS.size() + (MAXCHI2TRK - Chi2Sum/(Float_t)chikS.size())/(MAXCHI2TRK+1.);
- Results.ListHitID = HID;
- m_listResults.push_back(Results);
- //////////////////////////////////////////
- }
-
- //////////////////////
- FilterTrkCollection();
- //////////////////////
-
-
- if (m_listResults.size()!=0)
- {
- list<SIDRESULT>::const_iterator iter7;
- LogStream << MSG::INFO << "Filtered tracks parameters : X, DX, Y, DY, Z, quality: " << endreq;
- for (iter7=m_listResults.begin(); iter7!=m_listResults.end(); iter7++)
- {
- LogStream << MSG::INFO << std::fixed << std::setprecision(6) \
- << (*iter7).x_pos << "\t" << (*iter7).x_slope << "\t" << (*iter7).y_pos << "\t" << (*iter7).y_slope <<"\t"<< (*iter7).z_pos << "\t" << (*iter7).fChi2 << endreq;
- }
- LogStream << endreq;
- }
-
-
- return StatusCode::SUCCESS;
-}
-
-StatusCode AFP_SIDBasicKalman::Finalize(list<SIDRESULT>* pListResults)
-{
- *pListResults = m_listResults;
-
- //HistFinalize();
-
- return StatusCode::SUCCESS;
-}
-
-void AFP_SIDBasicKalman::FillSIDHITPOS(const SIDHIT &SIDHit, map<Int_t, SIDHITPOS> &MapSIDHitPos)
-{
- SIDHITPOS SIDHitPos;
- Int_t nStID = SIDHit.nStationID;
- Int_t nPlID = SIDHit.nDetectorID;
- Int_t nPixCol = SIDHit.nPixelCol;
- Int_t nPixRow = SIDHit.nPixelRow;
-
- SIDHitPos.nStationID = nStID;
- SIDHitPos.nPlateID = nPlID;
- SIDHitPos.fPixX = fxMapSID[nStID][nPlID][nPixCol][nPixRow];
- SIDHitPos.fPixY = fyMapSID[nStID][nPlID][nPixCol][nPixRow];
- SIDHitPos.fPixZ = fzMapSID[nStID][nPlID][nPixCol][nPixRow];
- SIDHitPos.fAmp = SIDHit.fADC;
-
- Int_t HitID = nStID*1000000+nPlID*100000+nPixCol*100+nPixRow;
-
- // Remove close hit pixel in X wrt signal amplitude
- map<Int_t, SIDHITPOS>::iterator iter0 = MapSIDHitPos.find(HitID-100);
- map<Int_t, SIDHITPOS>::iterator iter1 = MapSIDHitPos.find(HitID+100);
-
- if( iter0==MapSIDHitPos.end() && iter1==MapSIDHitPos.end())
- {
- MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
- }
-
- else if( iter0!=MapSIDHitPos.end() )
- {
- if ( (*iter0).second.fAmp < SIDHitPos.fAmp )
- {
- MapSIDHitPos.erase(iter0);
- MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
- }
- }
-
- else
- {
- if ( (*iter1).second.fAmp < SIDHitPos.fAmp )
- {
- MapSIDHitPos.erase(iter1);
- MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
- }
- }
-
-}
-
-void AFP_SIDBasicKalman::GetTrkSeeds()
-{
- map<Int_t, SIDHITPOS> MapSIDHitPosSeed0;
- map<Int_t, SIDHITPOS> MapSIDHitPosSeed1;
-
- list<SIDHIT>::const_iterator iter;
- for (iter=m_ListSIDHits.begin(); iter!=m_ListSIDHits.end(); iter++)
- {
- if ((*iter).fADC > m_AmpThresh && (*iter).nDetectorID == 0)
- {
- FillSIDHITPOS(*iter, MapSIDHitPosSeed0);
- }
-
- else if ((*iter).fADC > m_AmpThresh && (*iter).nDetectorID == 1)
- {
- FillSIDHITPOS(*iter, MapSIDHitPosSeed1);
- }
- }
-
-
- map<Int_t, SIDHITPOS>::const_iterator iter0;
- map<Int_t, SIDHITPOS>::const_iterator iter1;
- for (iter0=MapSIDHitPosSeed0.begin(); iter0!=MapSIDHitPosSeed0.end(); iter0++)
- {
- for (iter1=MapSIDHitPosSeed1.begin(); iter1!=MapSIDHitPosSeed1.end(); iter1++)
- {
- if ( (*iter0).second.nStationID == (*iter1).second.nStationID && 1 /* place for additional discriminant condition*/)
- {
- SIDHITSEED SIDHitSeed;
- SIDHitSeed.nHitID1 = (*iter0).first;
- SIDHitSeed.nHitID2 = (*iter1).first;
- SIDHitSeed.nLastPlate = 1;
- SIDHitSeed.nStationID = (*iter0).second.nStationID;
- SIDHitSeed.fSeedX = (*iter1).second.fPixX;
- SIDHitSeed.fSeedDX = ((*iter1).second.fPixX-(*iter0).second.fPixX)/abs((*iter1).second.fPixZ-(*iter0).second.fPixZ);
- SIDHitSeed.fSeedY = (*iter1).second.fPixY;
- SIDHitSeed.fSeedDY = ((*iter1).second.fPixY-(*iter0).second.fPixY)/abs((*iter1).second.fPixZ-(*iter0).second.fPixZ);
- SIDHitSeed.fSeedZ = (*iter1).second.fPixZ;
-
- pTrkSeeds.push_back(SIDHitSeed);
- }
- }
- }
-
-
-}
-
-
-bool AFP_SIDBasicKalman::FillTrkPropagators(const SIDHITSEED &SIDHitSeed, Int_t plateF)
-{
- HID.push_back(SIDHitSeed.nHitID1);
- HID.push_back(SIDHitSeed.nHitID2);
-
- Float_t Zdist = (fzMapSID[SIDHitSeed.nStationID][plateF][1][0] - fzMapSID[SIDHitSeed.nStationID][plateF][0][0])/(fxMapSID[SIDHitSeed.nStationID][plateF][1][0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0])*(SIDHitSeed.fSeedX - fxMapSID[SIDHitSeed.nStationID][plateF][0][0]) + fzMapSID[SIDHitSeed.nStationID][plateF][0][0] - SIDHitSeed.fSeedZ;
-
- CLHEP::HepMatrix Fi(4, 4, 1);
- Fi[0][1] = Zdist;
- Fi[2][3] = Zdist;
-
- CLHEP::HepVector x0i(4);
- x0i[0] = SIDHitSeed.fSeedX;
- x0i[1] = SIDHitSeed.fSeedDX;
- x0i[2] = SIDHitSeed.fSeedY;
- x0i[3] = SIDHitSeed.fSeedDY;
- CLHEP::HepVector m0i(2);
- m0i[0] = SIDHitSeed.fSeedX;
- m0i[1] = SIDHitSeed.fSeedY;
-
-
- CLHEP::HepVector xiP = Fi*x0i;
- CLHEP::HepMatrix CiP = Fi*C0*Fi.T();//+Qk;
-
- //mk finder
- CLHEP::HepVector mi(2);
- Float_t XYdist=100.*CLHEP::mm;
- Int_t HitID=-1;
- map<Int_t, SIDHITPOS>::const_iterator iter;
- for (iter=m_MapSIDHitPos.begin(); iter!=m_MapSIDHitPos.end(); iter++)
- {
- if( (*iter).second.nStationID != SIDHitSeed.nStationID ) continue;
- if( (*iter).second.nPlateID != plateF ) continue;
- Float_t minXYdist = sqrt( pow((*iter).second.fPixX-xiP[0],2)+pow((*iter).second.fPixY-xiP[2],2) );
- if( minXYdist < XYdist )
- {
- XYdist = minXYdist;
- mi[0] = (*iter).second.fPixX;
- mi[1] = (*iter).second.fPixY;
- HitID = (*iter).first;
- }
- }
- /////////////////////////
- if (HitID == -1) return 0; // no hits in layer
- /////////////////////////
-
- CLHEP::HepVector riP = mi - Hk*xiP;
- CLHEP::HepMatrix RiP = Vk + Hk*CiP*Hk.T();
-
- CLHEP::HepMatrix Ki = CiP*Hk.T()*qr_inverse(RiP);
- CLHEP::HepVector xi = xiP + Ki*riP;
-
- CLHEP::HepMatrix Ci = CiP - Ki*Hk*CiP;
-
- CLHEP::HepVector ri = mi - Hk*xi;
- CLHEP::HepMatrix Ri = Vk - Hk*Ki*Vk;
-
- //chi2 statistics to remove fakes
- CLHEP::HepVector chii = ri.T()*qr_inverse(Ri)*ri;
-
- //////////////////////////////
- if (chii[0]>MAXCHI2HIT) return 0; // no good hit in layer
- //////////////////////////////
-
- Fk.push_back(Fi);
- CkP.push_back(CiP);
- Ck.push_back(Ci);
- mk.push_back(mi);
- xk.push_back(xi);
- xkP.push_back(xiP);
- chik.push_back(chii[0]);
-
- m0 = m0i;
- x0 = x0i;
-
- HID.push_back(HitID);
- zk.push_back( (fzMapSID[SIDHitSeed.nStationID][plateF][1][0] - fzMapSID[SIDHitSeed.nStationID][plateF][0][0])/(fxMapSID[SIDHitSeed.nStationID][plateF][1][0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0])*(xi[0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0]) + fzMapSID[SIDHitSeed.nStationID][plateF][0][0] );
-
-
- return 1;
-}
-
-
-bool AFP_SIDBasicKalman::FillTrkPropagators(Int_t stID, Int_t plateF)
-{
- Float_t Zdist = (fzMapSID[stID][plateF][1][0] - fzMapSID[stID][plateF][0][0])/(fxMapSID[stID][plateF][1][0] - fxMapSID[stID][plateF][0][0])*(xk.back()[0] - fxMapSID[stID][plateF][0][0]) + fzMapSID[stID][plateF][0][0] - zk.back();
-
- CLHEP::HepMatrix Fi(4, 4, 1);
- Fi[0][1] = Zdist;
- Fi[2][3] = Zdist;
-
- CLHEP::HepVector xiP = Fi*xk.back();
- CLHEP::HepMatrix CiP = Fi*Ck.back()*Fi.T();//+Qk;
-
- //mk finder
- CLHEP::HepVector mi(2);
- Float_t XYdist=100.*CLHEP::mm;
- Int_t HitID=-1;
- map<Int_t, SIDHITPOS>::const_iterator iter;
- for (iter=m_MapSIDHitPos.begin(); iter!=m_MapSIDHitPos.end(); iter++)
- {
- if( (*iter).second.nStationID != stID ) continue;
- if( (*iter).second.nPlateID != plateF ) continue;
- Float_t minXYdist = sqrt( pow((*iter).second.fPixX-xiP[0],2)+pow((*iter).second.fPixY-xiP[2],2) );
- if( minXYdist < XYdist )
- {
- XYdist = minXYdist;
- mi[0] = (*iter).second.fPixX;
- mi[1] = (*iter).second.fPixY;
- HitID = (*iter).first;
- }
- }
- /////////////////////////
- if (HitID == -1) return 0; // no hits in layer
- /////////////////////////
-
- CLHEP::HepVector riP = mi - Hk*xiP;
- CLHEP::HepMatrix RiP = Vk + Hk*CiP*Hk.T();
-
- CLHEP::HepMatrix Ki = CiP*Hk.T()*qr_inverse(RiP);
- CLHEP::HepVector xi = xiP + Ki*riP;
-
- CLHEP::HepMatrix Ci = CiP - Ki*Hk*CiP;
-
- CLHEP::HepVector ri = mi - Hk*xi;
- CLHEP::HepMatrix Ri = Vk - Hk*Ki*Vk;
-
- //chi2 statistics to remove fakes
- CLHEP::HepVector chii = ri.T()*qr_inverse(Ri)*ri;
-
- //////////////////////////////
- if (chii[0]>MAXCHI2HIT) return 0; // no good hit in layer
- //////////////////////////////
-
- //cout<<Ck.back()*Fi.T()*qr_inverse(CiP)<<endl;
-
- Fk.push_back(Fi);
- CkP.push_back(CiP);
- Ck.push_back(Ci);
- mk.push_back(mi);
- xk.push_back(xi);
- xkP.push_back(xiP);
- chik.push_back(chii[0]);
-
- HID.push_back(HitID);
- zk.push_back( (fzMapSID[stID][plateF][1][0] - fzMapSID[stID][plateF][0][0])/(fxMapSID[stID][plateF][1][0] - fxMapSID[stID][plateF][0][0])*(xi[0] - fxMapSID[stID][plateF][0][0]) + fzMapSID[stID][plateF][0][0] );
-
- return 1;
-}
-
-
-void AFP_SIDBasicKalman::Smooth()
-{
- CLHEP::HepVector xiS;
- CLHEP::HepMatrix CiS;
-
- for (Int_t i=(xk.size()-2); i>=0; i--)
- {
- if(i==(xk.size()-2))
- {
- xiS = xk[xk.size()-1]; // !!!!!!!!!
- CiS = Ck[xk.size()-1]; // !!!!!!!!!
-
- xkS.push_back(xiS);
- CkS.push_back(CiS);
- CLHEP::HepVector riiS = mk[xk.size()-1] - Hk*xiS;
- CLHEP::HepMatrix RiiS = Vk - Hk*CiS*Hk.T();
- CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
- chikS.push_back(chiiS[0]);
- }
- else
- {
- xiS = xkS.back();
- CiS = CkS.back();
- }
-
- //CLHEP::HepMatrix Ai = qr_inverse(Fk[i+1]); // when Qk == 0
- CLHEP::HepMatrix Ai = Ck[i]*Fk[i+1].T()*qr_inverse(CkP[i+1]);
- CLHEP::HepVector xiiS = xk[i] + Ai*( xiS - xkP[i+1] );
- CLHEP::HepMatrix CiiS = Ck[i] + Ai*( CiS - CkP[i+1] )*Ai.T();
-
- CLHEP::HepVector riiS = mk[i] - Hk*xiiS;
- CLHEP::HepMatrix RiiS = Vk - Hk*CiiS*Hk.T();
- CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
-
- xkS.push_back(xiiS);
- CkS.push_back(CiiS);
- chikS.push_back(chiiS[0]);
- }
-
-
- //CLHEP::HepMatrix Ai = qr_inverse(Fk[0]); // when Qk == 0
- CLHEP::HepMatrix Ai = C0*Fk[0].T()*qr_inverse(CkP[0]);
- CLHEP::HepVector xiiS = x0 + Ai*( xkS.back() - xkP[0] );
- CLHEP::HepMatrix CiiS = C0 + Ai*( CkS.back() - CkP[0] )*Ai.T();
-
- CLHEP::HepVector riiS = m0 - Hk*xiiS;
- CLHEP::HepMatrix RiiS = Vk - Hk*CiiS*Hk.T();
- CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
-
- CkS.push_back(CiiS);
- xkS.push_back(xiiS);
- chikS.push_back(chiiS[0]);
-
-
-
- // first hit !!!!!!!!!!
- CLHEP::HepMatrix Fi(4, 4, 1);
- Fi[0][1] = (*m_MapSIDHitPos.find(HID[1])).second.fPixZ - (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
- Fi[2][3] = (*m_MapSIDHitPos.find(HID[1])).second.fPixZ - (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
-
- Ai = qr_inverse(Fi); // when Qk == 0
- xiiS = Ai*( xkS.back() );
- CiiS = C0 + Ai*( CkS.back() - C0 )*Ai.T();
-
- riiS = Hk*Ai*(x0 - xkS.back());
- RiiS = Vk - Hk*CiiS*Hk.T();
- chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
-
- CkS.push_back(CiiS);
- xkS.push_back(xiiS);
- chikS.push_back(chiiS[0]);
-
- // z position of the seed hit0
- //Int_t StID = (*m_MapSIDHitPos.find(HID[0])).second.nStationID;
- //Int_t PlID = (*m_MapSIDHitPos.find(HID[0])).second.nPlateID;
- //z0 = (fzMapSID[StID][PlID][1][0] - fzMapSID[StID][PlID][0][0])/(fxMapSID[StID][PlID][1][0] - fxMapSID[StID][PlID][0][0])*(xiiS[0] - fxMapSID[StID][PlID][0][0]) + fzMapSID[StID][PlID][0][0];
- z0 = (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
-}
-
-void AFP_SIDBasicKalman::FilterTrkCollection()
-{
- //filtering tracking collection using shared hits + quality requirement
- list<SIDRESULT>::iterator iter1;
- list<SIDRESULT>::iterator iter2;
- for (iter1=m_listResults.begin(); iter1!=m_listResults.end(); iter1++)
- {
- for (iter2=m_listResults.begin(); iter2!=m_listResults.end(); )
- {
- if ( GetSharedHits((*iter1).ListHitID, (*iter2).ListHitID) > MAXSHAREDHITS && iter1!=iter2 )
- {
- if ( (*iter1).fChi2 < (*iter2).fChi2 )
- {
- //m_listResults.erase(iter1);
- //--iter1;
- ++iter2;
- }
- else
- {
- iter2 = m_listResults.erase(iter2);
- //--iter2;
- }
- }
- else ++iter2;
- }
-
- }
-
-}
-
-Int_t AFP_SIDBasicKalman::GetSharedHits(const vector<Int_t> &HID1, const vector<Int_t> &HID2)
-{
- Int_t SharedHits = 0;
-
- vector<Int_t>::const_iterator iter1;
- vector<Int_t>::const_iterator iter2;
-
- for (iter1=HID1.begin(); iter1!=HID1.end(); iter1++)
- {
- for (iter2=HID2.begin(); iter2!=HID2.end(); iter2++)
- {
- if ((*iter1) == (*iter2)) SharedHits++;
- }
- }
-
- return SharedHits;
-}
-
-void AFP_SIDBasicKalman::ClearMatrix()
-{
- Fk.clear();
- Ck.clear();
- CkP.clear();
- Rk.clear();
-
- xk.clear();
- xkP.clear();
- rk.clear();
- rkP.clear();
- mk.clear();
- chik.clear();
-
- zk.clear();
- HID.clear();
-
- xkS.clear();
- CkS.clear();
- chikS.clear();
-
-
-}
-
-void AFP_SIDBasicKalman::HistInitialize()
-{
- //supporting histograms with binning convention: z, x
- MyFile = new TFile("PixMap.root","RECREATE");
- histS1_PixMap = new TH2F("histS1x_PixMap", "", 600, fzMapSID[1][0][0][0]-5., fzMapSID[1][5][335][0]+5., 400, fxMapSID[1][5][335][0]-2., fxMapSID[1][0][0][0]+2.);
- histS2_PixMap = new TH2F("histS2x_PixMap", "", 600, fzMapSID[0][0][0][0]-5., fzMapSID[0][5][335][0]+5., 400, fxMapSID[0][5][335][0]-2., fxMapSID[0][0][0][0]+2.);
- histS3_PixMap = new TH2F("histS1y_PixMap", "", 600, fzMapSID[1][0][0][0]-5., fzMapSID[1][5][335][0]+5., 400, fyMapSID[1][0][0][79]-2., fyMapSID[1][0][0][0]+2.);
- histS4_PixMap = new TH2F("histS2y_PixMap", "", 600, fzMapSID[0][0][0][0]-5., fzMapSID[0][5][335][0]+5., 400, fyMapSID[0][0][0][79]-2., fyMapSID[0][0][0][0]+2.);
-}
-
-void AFP_SIDBasicKalman::HistFinalize()
-{
-
- MyFile->Write();
- //MyFile->Close();
- delete histS1_PixMap;
- delete histS2_PixMap;
- delete histS3_PixMap;
- delete histS4_PixMap;
- delete MyFile;
-}
-
-
-void AFP_SIDBasicKalman::GetData()
-{
- MsgStream LogStream(Athena::getMessageSvc(), "AFP_SIDBasicKalman::GetData()");
- LogStream << MSG::DEBUG << "begin AFP_SIDBasicKalman::GetData()" << endreq;
-
- /////
- /////
-
- LogStream << MSG::DEBUG << "end AFP_SIDBasicKalman::GetData()" << endreq;
-}
+#include "AFP_LocReco/AFP_SIDBasicKalman.h"
+
+AFP_SIDBasicKalman::AFP_SIDBasicKalman()
+{
+ MyFile = NULL;
+ histS1_PixMap = NULL;
+ histS2_PixMap = NULL;
+ histS3_PixMap = NULL;
+ histS4_PixMap = NULL;
+
+ m_listResults.clear();
+
+ m_AmpThresh = 0;
+ m_iDataType = 0;
+
+ pTrkSeeds.clear();
+
+ Fk.clear();
+ Ck.clear();
+ CkP.clear();
+ Rk.clear();
+
+ xk.clear();
+ xkP.clear();
+ rk.clear();
+ rkP.clear();
+ mk.clear();
+ chik.clear();
+
+ zk.clear();
+ HID.clear();
+
+ xkS.clear();
+ CkS.clear();
+ chikS.clear();
+
+ z0 = 0;
+}
+
+AFP_SIDBasicKalman::~AFP_SIDBasicKalman()
+{
+ m_listResults.clear();
+}
+
+StatusCode AFP_SIDBasicKalman::Initialize(float fAmpThresh, int iDataType, const list<SIDHIT> &ListSIDHits, Float_t fsSID[SIDSTATIONID][SIDCNT], Float_t fxSID[SIDSTATIONID][SIDCNT], Float_t fySID[SIDSTATIONID][SIDCNT], Float_t fzSID[SIDSTATIONID][SIDCNT])
+{
+
+ m_AmpThresh = (float)fAmpThresh;
+ m_AmpThresh = 1000;
+ m_iDataType = iDataType;
+ m_ListSIDHits = ListSIDHits;
+
+ Float_t delta_pixel_x = 0.050*CLHEP::mm; // size of pixel along x, in mm
+ Float_t delta_pixel_y = 0.250*CLHEP::mm; // size of pixel along y, in mm
+
+ Float_t interPlaneX[4][4] ={0,0,0,0,
+ 0,0,0,0,
+ 0,-0.187358,-0.127971,0,
+ 0,-0.101478,-0.0661546,-0.0675869};
+ Float_t interPlaneY[4][4] ={0,0,0,0,
+ 0,0,0,0,
+ 0,-0.134756, -0.204807,0,
+ 0, -0.155841, -0.334444,-0.341143};
+ Float_t Alpha[4][4] ={0,0,0,0,
+ 0,0,0,0,
+ 0.2443461,0.253198,0.227862,0,
+ 0.2443461,0.250661,0.252849,0.247896};
+
+
+ // x-y-z map of all pixels
+ for(Int_t nStationID = 0; nStationID < SIDSTATIONID; nStationID++)
+ {
+ for (Int_t nPlateID = 0; nPlateID < SIDCNT; nPlateID++)
+ {
+
+ if ( iDataType == 1) fsSID[nStationID][nPlateID] = Alpha[nStationID][nPlateID];
+
+ for (Int_t nPixel_x = 0; nPixel_x < 336; nPixel_x++)
+ {
+ for (Int_t nPixel_y = 0; nPixel_y < 80; nPixel_y++)
+ {
+ fxMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fxSID[nStationID][nPlateID]+(delta_pixel_x*(nPixel_x-168))*cos(fsSID[nStationID][nPlateID]); //sign changed!
+ fyMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fySID[nStationID][nPlateID]+(delta_pixel_y*nPixel_y); //sign changed!
+ fzMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] = fzSID[nStationID][nPlateID] - ((nStationID<2)?1.:-1.)*(delta_pixel_x*(nPixel_x-168))*sin(fsSID[nStationID][nPlateID]); //sign changed!
+ if( iDataType == 1) {
+ fxMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] += interPlaneX[nStationID][nPlateID] + 168*delta_pixel_x*cos(fsSID[nStationID][0]);
+ fyMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] += interPlaneY[nStationID][nPlateID];
+ fzMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] += -168*((nStationID<2)?1.:-1.)*delta_pixel_x*sin(fsSID[nStationID][0]); }
+ else {
+ fxMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] += 168*delta_pixel_x*cos(fsSID[nStationID][0]);
+ fzMapSID[nStationID][nPlateID][nPixel_x][nPixel_y] += -168*((nStationID<2)?1.:-1.)*delta_pixel_x*sin(fsSID[nStationID][0]); }
+ }
+ }
+ }
+ }
+
+ // x-y-z map of hit pixels
+ list<SIDHIT>::const_iterator iter;
+ for (iter=ListSIDHits.begin(); iter!=ListSIDHits.end(); iter++)
+ {
+ if ((*iter).fADC > m_AmpThresh)
+ {
+ FillSIDHITPOS(*iter, m_MapSIDHitPos);
+
+ /*SIDHITPOS SIDHitPos;
+ Int_t nStID = (*iter).nStationID;
+ Int_t nPlID = (*iter).nDetectorID;
+ Int_t nPixCol = (*iter).nPixelCol;
+ Int_t nPixRow = (*iter).nPixelRow;
+
+
+ SIDHitPos.nStationID = nStID;
+ SIDHitPos.nPlateID = nPlID;
+ SIDHitPos.fPixX = fxSID[nStID][nPlID] - (delta_pixel_x*nPixCol+delta_pixel_x/2.)*cos(fsSID[nStID][nPlID]);
+ SIDHitPos.fPixY = fySID[nStID][nPlID] - (delta_pixel_y*nPixRow+delta_pixel_y/2.);
+ SIDHitPos.fPixZ = ((nStID<2)?1.:-1.)*fzSID[nStID][nPlID] + ((nStID<2)?1.:-1.)*(delta_pixel_x*nPixCol+delta_pixel_x/2.)*sin(fsSID[nStID][nPlID]);
+
+ Int_t HitID = nStID*1000000+nPlID*100000+nPixCol*100+nPixRow;
+
+ if(m_MapSIDHitPos.find(HitID-100)==m_MapSIDHitPos.end() && m_MapSIDHitPos.find(HitID+100)==m_MapSIDHitPos.end())
+ {
+ m_MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
+ }*/
+ }
+ }
+
+ // Kalman matrices initialization
+ Hk = CLHEP::HepMatrix(2, 4, 0);
+ Hk[0][0] = 1.;
+ Hk[1][2] = 1.;
+
+ Qk = CLHEP::HepMatrix(4, 4, 0);
+
+ Vk = CLHEP::HepMatrix(2, 2, 0);
+ Vk[0][0] = delta_pixel_x*delta_pixel_x;
+ Vk[1][1] = delta_pixel_y*delta_pixel_y;
+
+ C0 = CLHEP::HepMatrix(4, 4, 0);
+ C0[0][0] = delta_pixel_x*delta_pixel_x/3.;
+ C0[1][1] = delta_pixel_x*delta_pixel_x/SID_NOMINALSPACING/SID_NOMINALSPACING/3.;
+ //C0[1][0] = 0;
+ //C0[0][1] = 0;
+
+ C0[2][2] = delta_pixel_y*delta_pixel_y/3.;
+ C0[3][3] = delta_pixel_y*delta_pixel_y/SID_NOMINALSPACING/SID_NOMINALSPACING/3.;
+ //C0[3][2] = 0;
+ //C0[2][3] = 0;
+
+ //HistInitialize();
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode AFP_SIDBasicKalman::Execute()
+{
+ MsgStream LogStream(Athena::getMessageSvc(), "AFP_SIDBasicKalman::Execute()");
+
+// /*
+ map<Int_t, SIDHITPOS>::const_iterator iter0;
+ for (iter0=m_MapSIDHitPos.begin(); iter0!=m_MapSIDHitPos.end(); iter0++)
+ {
+ /*
+ if( (*iter0).second.nStationID==1 )
+ {
+ histS1_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixX/CLHEP::mm);
+ histS3_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixY/CLHEP::mm);
+ }
+ else if( (*iter0).second.nStationID==0 )
+ {
+ histS2_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixX/CLHEP::mm);
+ histS4_PixMap->Fill((*iter0).second.fPixZ/CLHEP::mm, (*iter0).second.fPixY/CLHEP::mm);
+ }
+ */
+ }
+// */
+
+
+ ////////////////////////
+ GetTrkSeeds();
+ ////////////////////////
+
+ vector<SIDHITSEED>::const_iterator iter1;
+ for (iter1= pTrkSeeds.begin(); iter1!= pTrkSeeds.end(); iter1++)
+ {
+
+ ClearMatrix();
+
+ Int_t LastPlate = (*iter1).nLastPlate;
+ Int_t DoubleHole1=1;
+ Int_t DoubleHole2=0;
+
+ if ( FillTrkPropagators(*iter1, ++LastPlate) )
+ {
+ ++LastPlate;
+ DoubleHole1=0;
+ for (Int_t i=LastPlate; i<SIDCNT; i++)
+ {
+ if ( FillTrkPropagators((*iter1).nStationID, i) )
+ {
+ if(DoubleHole1==1) DoubleHole1=0;
+ }
+ else
+ {
+ ++DoubleHole1;
+ }
+ if(DoubleHole1==2) break;
+
+ }
+ if(DoubleHole1==2) continue;
+ }
+
+ else if ( FillTrkPropagators(*iter1, ++LastPlate) )
+ {
+ ++LastPlate;
+ DoubleHole2=0;
+ for (Int_t i=LastPlate; i<SIDCNT; i++)
+ {
+ if ( FillTrkPropagators((*iter1).nStationID, i) )
+ {
+ if(DoubleHole2==1) DoubleHole2=0;
+ }
+ else
+ {
+ ++DoubleHole2;
+ }
+ if(DoubleHole2==2) break;
+
+ }
+ if(DoubleHole2==2) continue;
+ }
+
+ else continue; // double hole after the seed pixels
+
+ if( DoubleHole1 > 1 || DoubleHole2 != 0) continue ;
+
+
+ LogStream << MSG::DEBUG << "Found new track candidate with parameters:" << endreq;
+ LogStream << MSG::DEBUG << "Hit ID's : ";
+ vector<Int_t>::const_iterator iter2;
+ for (iter2= HID.begin(); iter2!= HID.end(); iter2++) LogStream << MSG::DEBUG << (*iter2) << "\t";
+ LogStream << endreq;
+
+ LogStream << MSG::DEBUG << "Filtered parameters : X, DX, Y, DY, chi2" << endreq;
+ vector< CLHEP::HepVector >::const_iterator iter3;
+ vector< Float_t >::const_iterator iter4 = chik.begin();
+ LogStream << MSG::DEBUG << (*iter1).fSeedX << "\t" << (*iter1).fSeedDX << "\t" << (*iter1).fSeedY << "\t" << (*iter1).fSeedDY << "\t" << 0 << endreq;
+ for (iter3=xk.begin(); iter3!= xk.end(); iter3++)
+ {
+ LogStream << MSG::DEBUG << (*iter3)[0] << "\t" << (*iter3)[1] << "\t" << (*iter3)[2] << "\t" << (*iter3)[3] << "\t" << (*iter4) << endreq;
+ iter4++;
+ }
+
+
+ ////////////
+ Smooth();
+ ////////////
+
+ Float_t Chi2Sum = 0;
+
+ LogStream << MSG::DEBUG << "Smoothed parameters : X, DX, Y, DY, chi2" << endreq;
+ vector< CLHEP::HepVector >::const_reverse_iterator iter5;
+ vector< Float_t >::const_reverse_iterator iter6 = chikS.rbegin();
+ for (iter5=xkS.rbegin(); iter5!= xkS.rend(); ++iter5)
+ {
+ LogStream << MSG::DEBUG << (*iter5)[0] << "\t" << (*iter5)[1] << "\t" << (*iter5)[2] << "\t" << (*iter5)[3] << "\t" << (*iter6) << endreq;
+ Chi2Sum += (*iter6);
+ ++iter6;
+ }
+ LogStream << endreq;
+
+
+ //////////////////////////////////////////
+ // fill tracking collection
+ SIDRESULT Results;
+ Results.nStationID = (*iter1).nStationID;
+ Results.x_pos = xkS.back()[0];
+ Results.y_pos = xkS.back()[2];
+ Results.z_pos = z0;
+ Results.x_slope = xkS.back()[1];
+ Results.y_slope = xkS.back()[3];
+ Results.z_slope = 1.;
+ Results.nHits = xkS.size();
+ Results.nHoles = SIDCNT - xkS.size();
+ Results.fChi2 = xkS.size() + (MAXCHI2TRK - Chi2Sum/(Float_t)chikS.size())/(MAXCHI2TRK+1.);
+ Results.ListHitID = HID;
+ m_listResults.push_back(Results);
+ //////////////////////////////////////////
+ }
+
+
+ //////////////////////
+ FilterTrkCollection();
+ //////////////////////
+
+
+ if (m_listResults.size()!=0)
+ {
+ list<SIDRESULT>::const_iterator iter7;
+ LogStream << MSG::INFO << "Filtered tracks parameters : X, DX, Y, DY, Z, quality: " << endreq;
+ for (iter7=m_listResults.begin(); iter7!=m_listResults.end(); iter7++)
+ {
+ LogStream << MSG::INFO << std::fixed << std::setprecision(6) \
+ << (*iter7).x_pos << "\t" << (*iter7).x_slope << "\t" << (*iter7).y_pos << "\t" << (*iter7).y_slope <<"\t"<< (*iter7).z_pos << "\t" << (*iter7).fChi2 << endreq;
+ }
+ LogStream << endreq;
+ }
+
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode AFP_SIDBasicKalman::Finalize(list<SIDRESULT>* pListResults)
+{
+ *pListResults = m_listResults;
+
+ //HistFinalize();
+
+ return StatusCode::SUCCESS;
+}
+
+void AFP_SIDBasicKalman::FillSIDHITPOS(const SIDHIT &SIDHit, map<Int_t, SIDHITPOS> &MapSIDHitPos)
+{
+ SIDHITPOS SIDHitPos;
+
+
+ Int_t nStID = SIDHit.nStationID;
+ Int_t nPlID = SIDHit.nDetectorID;
+ Int_t nPixCol = SIDHit.nPixelCol;
+ Int_t nPixRow = SIDHit.nPixelRow;
+
+ // if ( m_iDataType ==1 ) {
+ // nPixCol = 336-SIDHit.nPixelRow;
+ // nPixRow = SIDHit.nPixelCol;
+ // }
+
+ SIDHitPos.nStationID = nStID;
+ SIDHitPos.nPlateID = nPlID;
+ SIDHitPos.fPixX = fxMapSID[nStID][nPlID][nPixCol][nPixRow];
+ SIDHitPos.fPixY = fyMapSID[nStID][nPlID][nPixCol][nPixRow];
+ SIDHitPos.fPixZ = fzMapSID[nStID][nPlID][nPixCol][nPixRow];
+ SIDHitPos.fAmp = SIDHit.fADC;
+
+ Int_t HitID = nStID*1000000+nPlID*100000+nPixCol*100+nPixRow;
+ // Remove close hit pixel in X wrt signal amplitude
+ // Actually do simple clustering
+ map<Int_t, SIDHITPOS>::iterator iter0 = MapSIDHitPos.find(HitID-100);
+ map<Int_t, SIDHITPOS>::iterator iter1 = MapSIDHitPos.find(HitID+100);
+
+ if( iter0==MapSIDHitPos.end() && iter1==MapSIDHitPos.end())
+ {
+ MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
+ }
+
+ else if( iter0!=MapSIDHitPos.end() )
+ {
+ if ( (*iter0).second.fAmp < SIDHitPos.fAmp )
+ {
+
+
+ SIDHitPos.fPixX = SIDHitPos.fPixX*SIDHitPos.fAmp + (*iter0).second.fPixX*(*iter0).second.fAmp;
+ SIDHitPos.fPixY = SIDHitPos.fPixY*SIDHitPos.fAmp + (*iter0).second.fPixY*(*iter0).second.fAmp;
+ SIDHitPos.fPixZ = SIDHitPos.fPixZ*SIDHitPos.fAmp + (*iter0).second.fPixZ*(*iter0).second.fAmp;
+
+ SIDHitPos.fAmp += (*iter0).second.fAmp;
+ SIDHitPos.fPixX /=SIDHitPos.fAmp;
+ SIDHitPos.fPixY /=SIDHitPos.fAmp;
+ SIDHitPos.fPixZ /=SIDHitPos.fAmp;
+
+ MapSIDHitPos.erase(iter0);
+ MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
+ }
+ }
+
+ else
+ {
+ if ( (*iter1).second.fAmp < SIDHitPos.fAmp )
+ {
+
+ SIDHitPos.fPixX = SIDHitPos.fPixX*SIDHitPos.fAmp + (*iter1).second.fPixX*(*iter1).second.fAmp;
+ SIDHitPos.fPixY = SIDHitPos.fPixY*SIDHitPos.fAmp + (*iter1).second.fPixY*(*iter1).second.fAmp;
+ SIDHitPos.fPixZ = SIDHitPos.fPixZ*SIDHitPos.fAmp + (*iter1).second.fPixZ*(*iter1).second.fAmp;
+
+ SIDHitPos.fAmp += (*iter1).second.fAmp;
+ SIDHitPos.fPixX /=SIDHitPos.fAmp;
+ SIDHitPos.fPixY /=SIDHitPos.fAmp;
+ SIDHitPos.fPixZ /=SIDHitPos.fAmp;
+
+ MapSIDHitPos.erase(iter1);
+ MapSIDHitPos.insert(pair<Int_t, SIDHITPOS>(HitID, SIDHitPos));
+ }
+ }
+
+}
+
+void AFP_SIDBasicKalman::GetTrkSeeds()
+{
+ map<Int_t, SIDHITPOS> MapSIDHitPosSeed0;
+ map<Int_t, SIDHITPOS> MapSIDHitPosSeed1;
+
+ list<SIDHIT>::const_iterator iter;
+ for (iter=m_ListSIDHits.begin(); iter!=m_ListSIDHits.end(); iter++)
+ {
+ if ((*iter).fADC > m_AmpThresh && (*iter).nDetectorID == 0)
+ {
+ FillSIDHITPOS(*iter, MapSIDHitPosSeed0);
+ }
+
+ else if ((*iter).fADC > m_AmpThresh && (*iter).nDetectorID == 1)
+ {
+ FillSIDHITPOS(*iter, MapSIDHitPosSeed1);
+ }
+ }
+
+
+ map<Int_t, SIDHITPOS>::const_iterator iter0;
+ map<Int_t, SIDHITPOS>::const_iterator iter1;
+ for (iter0=MapSIDHitPosSeed0.begin(); iter0!=MapSIDHitPosSeed0.end(); iter0++)
+ {
+ for (iter1=MapSIDHitPosSeed1.begin(); iter1!=MapSIDHitPosSeed1.end(); iter1++)
+ {
+ if ( (*iter0).second.nStationID == (*iter1).second.nStationID && 1 /* place for additional discriminant condition*/)
+ {
+ SIDHITSEED SIDHitSeed;
+ SIDHitSeed.nHitID1 = (*iter0).first;
+ SIDHitSeed.nHitID2 = (*iter1).first;
+ SIDHitSeed.nLastPlate = 1;
+ SIDHitSeed.nStationID = (*iter0).second.nStationID;
+ SIDHitSeed.fSeedX = (*iter1).second.fPixX;
+ SIDHitSeed.fSeedDX = ((*iter1).second.fPixX-(*iter0).second.fPixX)/abs((*iter1).second.fPixZ-(*iter0).second.fPixZ);
+ SIDHitSeed.fSeedY = (*iter1).second.fPixY;
+ SIDHitSeed.fSeedDY = ((*iter1).second.fPixY-(*iter0).second.fPixY)/abs((*iter1).second.fPixZ-(*iter0).second.fPixZ);
+ SIDHitSeed.fSeedZ = (*iter1).second.fPixZ;
+
+ pTrkSeeds.push_back(SIDHitSeed);
+ }
+ }
+ }
+
+
+}
+
+
+bool AFP_SIDBasicKalman::FillTrkPropagators(const SIDHITSEED &SIDHitSeed, Int_t plateF)
+{
+ // This sensoor is not active
+
+ if( SIDHitSeed.nStationID == 2 && plateF==3) return 0;
+
+ HID.clear();
+ HID.push_back(SIDHitSeed.nHitID1);
+ HID.push_back(SIDHitSeed.nHitID2);
+
+ Float_t Zdist = (fzMapSID[SIDHitSeed.nStationID][plateF][1][0] - fzMapSID[SIDHitSeed.nStationID][plateF][0][0])/(fxMapSID[SIDHitSeed.nStationID][plateF][1][0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0])*(SIDHitSeed.fSeedX - fxMapSID[SIDHitSeed.nStationID][plateF][0][0]) + fzMapSID[SIDHitSeed.nStationID][plateF][0][0] - SIDHitSeed.fSeedZ;
+
+ CLHEP::HepMatrix Fi(4, 4, 1);
+ Fi[0][1] = Zdist;
+ Fi[2][3] = Zdist;
+
+ CLHEP::HepVector x0i(4);
+ x0i[0] = SIDHitSeed.fSeedX;
+ x0i[1] = SIDHitSeed.fSeedDX;
+ x0i[2] = SIDHitSeed.fSeedY;
+ x0i[3] = SIDHitSeed.fSeedDY;
+ CLHEP::HepVector m0i(2);
+ m0i[0] = SIDHitSeed.fSeedX;
+ m0i[1] = SIDHitSeed.fSeedY;
+
+
+ CLHEP::HepVector xiP = Fi*x0i;
+ CLHEP::HepMatrix CiP = Fi*C0*Fi.T();//+Qk;
+
+ //mk finder
+ CLHEP::HepVector mi(2);
+ Float_t XYdist=100.*CLHEP::mm;
+ Int_t HitID=-1;
+ map<Int_t, SIDHITPOS>::const_iterator iter;
+ for (iter=m_MapSIDHitPos.begin(); iter!=m_MapSIDHitPos.end(); iter++)
+ {
+ if( (*iter).second.nStationID != SIDHitSeed.nStationID ) continue;
+ if( (*iter).second.nPlateID != plateF ) continue;
+ Float_t minXYdist = sqrt( pow((*iter).second.fPixX-xiP[0],2)+pow((*iter).second.fPixY-xiP[2],2) );
+ if( minXYdist < XYdist )
+ {
+ XYdist = minXYdist;
+ mi[0] = (*iter).second.fPixX;
+ mi[1] = (*iter).second.fPixY;
+ HitID = (*iter).first;
+ }
+ }
+ /////////////////////////
+ if (HitID == -1) return 0; // no hits in layer
+ /////////////////////////
+
+ CLHEP::HepVector riP = mi - Hk*xiP;
+ CLHEP::HepMatrix RiP = Vk + Hk*CiP*Hk.T();
+
+ CLHEP::HepMatrix Ki = CiP*Hk.T()*qr_inverse(RiP);
+ CLHEP::HepVector xi = xiP + Ki*riP;
+
+ CLHEP::HepMatrix Ci = CiP - Ki*Hk*CiP;
+
+ CLHEP::HepVector ri = mi - Hk*xi;
+ CLHEP::HepMatrix Ri = Vk - Hk*Ki*Vk;
+
+ //chi2 statistics to remove fakes
+ CLHEP::HepVector chii = ri.T()*qr_inverse(Ri)*ri;
+
+ //////////////////////////////
+ if (chii[0]>MAXCHI2HIT) return 0; // no good hit in layer
+ //////////////////////////////
+
+ Fk.push_back(Fi);
+ CkP.push_back(CiP);
+ Ck.push_back(Ci);
+ mk.push_back(mi);
+ xk.push_back(xi);
+ xkP.push_back(xiP);
+ chik.push_back(chii[0]);
+
+ m0 = m0i;
+ x0 = x0i;
+
+ HID.push_back(HitID);
+ zk.push_back( (fzMapSID[SIDHitSeed.nStationID][plateF][1][0] - fzMapSID[SIDHitSeed.nStationID][plateF][0][0])/(fxMapSID[SIDHitSeed.nStationID][plateF][1][0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0])*(xi[0] - fxMapSID[SIDHitSeed.nStationID][plateF][0][0]) + fzMapSID[SIDHitSeed.nStationID][plateF][0][0] );
+
+
+ return 1;
+}
+
+
+bool AFP_SIDBasicKalman::FillTrkPropagators(Int_t stID, Int_t plateF)
+{
+ Float_t Zdist = (fzMapSID[stID][plateF][1][0] - fzMapSID[stID][plateF][0][0])/(fxMapSID[stID][plateF][1][0] - fxMapSID[stID][plateF][0][0])*(xk.back()[0] - fxMapSID[stID][plateF][0][0]) + fzMapSID[stID][plateF][0][0] - zk.back();
+
+ // This sensor is not active
+ if( stID == 2 && plateF == 3) return 0;
+
+ CLHEP::HepMatrix Fi(4, 4, 1);
+ Fi[0][1] = Zdist;
+ Fi[2][3] = Zdist;
+
+ CLHEP::HepVector xiP = Fi*xk.back();
+ CLHEP::HepMatrix CiP = Fi*Ck.back()*Fi.T();//+Qk;
+
+ //mk finder
+ CLHEP::HepVector mi(2);
+ Float_t XYdist=100.*CLHEP::mm;
+ Int_t HitID=-1;
+ map<Int_t, SIDHITPOS>::const_iterator iter;
+ for (iter=m_MapSIDHitPos.begin(); iter!=m_MapSIDHitPos.end(); iter++)
+ {
+ if( (*iter).second.nStationID != stID ) continue;
+ if( (*iter).second.nPlateID != plateF ) continue;
+ Float_t minXYdist = sqrt( pow((*iter).second.fPixX-xiP[0],2)+pow((*iter).second.fPixY-xiP[2],2) );
+ if( minXYdist < XYdist )
+ {
+ XYdist = minXYdist;
+ mi[0] = (*iter).second.fPixX;
+ mi[1] = (*iter).second.fPixY;
+ HitID = (*iter).first;
+ }
+ }
+
+
+ /////////////////////////
+ if (HitID == -1) return 0; // no hits in layer
+ /////////////////////////
+
+ CLHEP::HepVector riP = mi - Hk*xiP;
+ CLHEP::HepMatrix RiP = Vk + Hk*CiP*Hk.T();
+
+ CLHEP::HepMatrix Ki = CiP*Hk.T()*qr_inverse(RiP);
+ CLHEP::HepVector xi = xiP + Ki*riP;
+
+ CLHEP::HepMatrix Ci = CiP - Ki*Hk*CiP;
+
+ CLHEP::HepVector ri = mi - Hk*xi;
+ CLHEP::HepMatrix Ri = Vk - Hk*Ki*Vk;
+
+ //chi2 statistics to remove fakes
+ CLHEP::HepVector chii = ri.T()*qr_inverse(Ri)*ri;
+
+
+ //////////////////////////////
+ if (chii[0]>MAXCHI2HIT) return 0; // no good hit in layer
+ //////////////////////////////
+
+
+ Fk.push_back(Fi);
+ CkP.push_back(CiP);
+ Ck.push_back(Ci);
+ mk.push_back(mi);
+ xk.push_back(xi);
+ xkP.push_back(xiP);
+ chik.push_back(chii[0]);
+
+ HID.push_back(HitID);
+ zk.push_back( (fzMapSID[stID][plateF][1][0] - fzMapSID[stID][plateF][0][0])/(fxMapSID[stID][plateF][1][0] - fxMapSID[stID][plateF][0][0])*(xi[0] - fxMapSID[stID][plateF][0][0]) + fzMapSID[stID][plateF][0][0] );
+
+ return 1;
+}
+
+
+void AFP_SIDBasicKalman::Smooth()
+{
+ CLHEP::HepVector xiS;
+ CLHEP::HepMatrix CiS;
+
+
+ if( xk.size() > 1) {
+ for (Int_t i=(xk.size()-2); i>=0; i--)
+ {
+ if(i==(xk.size()-2))
+ {
+ xiS = xk[xk.size()-1]; // !!!!!!!!!
+ CiS = Ck[xk.size()-1]; // !!!!!!!!!
+
+ xkS.push_back(xiS);
+ CkS.push_back(CiS);
+ CLHEP::HepVector riiS = mk[xk.size()-1] - Hk*xiS;
+ CLHEP::HepMatrix RiiS = Vk - Hk*CiS*Hk.T();
+ CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
+ chikS.push_back(chiiS[0]);
+ }
+ else
+ {
+ xiS = xkS.back();
+ CiS = CkS.back();
+ }
+
+ //CLHEP::HepMatrix Ai = qr_inverse(Fk[i+1]); // when Qk == 0
+ CLHEP::HepMatrix Ai = Ck[i]*Fk[i+1].T()*qr_inverse(CkP[i+1]);
+ CLHEP::HepVector xiiS = xk[i] + Ai*( xiS - xkP[i+1] );
+ CLHEP::HepMatrix CiiS = Ck[i] + Ai*( CiS - CkP[i+1] )*Ai.T();
+
+ CLHEP::HepVector riiS = mk[i] - Hk*xiiS;
+ CLHEP::HepMatrix RiiS = Vk - Hk*CiiS*Hk.T();
+ CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
+
+ xkS.push_back(xiiS);
+ CkS.push_back(CiiS);
+ chikS.push_back(chiiS[0]);
+ }
+ }
+ else if (xk.size() > 0) {
+ xiS = xk[xk.size()-1]; // !!!!!!!!!
+ CiS = Ck[xk.size()-1]; // !!!!!!!!!
+
+ xkS.push_back(xiS);
+ CkS.push_back(CiS);
+ CLHEP::HepVector riiS = mk[xk.size()-1] - Hk*xiS;
+ CLHEP::HepMatrix RiiS = Vk - Hk*CiS*Hk.T();
+ CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
+ chikS.push_back(chiiS[0]);
+
+
+ }
+
+
+
+
+ //CLHEP::HepMatrix Ai = qr_inverse(Fk[0]); // when Qk == 0
+ CLHEP::HepMatrix Ai = C0*Fk[0].T()*qr_inverse(CkP[0]);
+
+ CLHEP::HepVector xiiS = x0 + Ai*( xkS.back() - xkP[0] );
+
+ CLHEP::HepMatrix CiiS = C0 + Ai*( CkS.back() - CkP[0] )*Ai.T();
+
+ CLHEP::HepVector riiS = m0 - Hk*xiiS;
+ CLHEP::HepMatrix RiiS = Vk - Hk*CiiS*Hk.T();
+
+
+ CLHEP::HepVector chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
+
+ CkS.push_back(CiiS);
+ xkS.push_back(xiiS);
+ chikS.push_back(chiiS[0]);
+
+
+
+ // first hit !!!!!!!!!!
+ CLHEP::HepMatrix Fi(4, 4, 1);
+ Fi[0][1] = (*m_MapSIDHitPos.find(HID[1])).second.fPixZ - (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
+ Fi[2][3] = (*m_MapSIDHitPos.find(HID[1])).second.fPixZ - (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
+
+ Ai = qr_inverse(Fi); // when Qk == 0
+ xiiS = Ai*( xkS.back() );
+ CiiS = C0 + Ai*( CkS.back() - C0 )*Ai.T();
+
+ riiS = Hk*Ai*(x0 - xkS.back());
+ RiiS = Vk - Hk*CiiS*Hk.T();
+ chiiS = riiS.T()*qr_inverse(RiiS)*riiS;
+
+ CkS.push_back(CiiS);
+ xkS.push_back(xiiS);
+ chikS.push_back(chiiS[0]);
+
+ // z position of the seed hit0
+ //Int_t StID = (*m_MapSIDHitPos.find(HID[0])).second.nStationID;
+ //Int_t PlID = (*m_MapSIDHitPos.find(HID[0])).second.nPlateID;
+ //z0 = (fzMapSID[StID][PlID][1][0] - fzMapSID[StID][PlID][0][0])/(fxMapSID[StID][PlID][1][0] - fxMapSID[StID][PlID][0][0])*(xiiS[0] - fxMapSID[StID][PlID][0][0]) + fzMapSID[StID][PlID][0][0];
+ z0 = (*m_MapSIDHitPos.find(HID[0])).second.fPixZ;
+}
+
+void AFP_SIDBasicKalman::FilterTrkCollection()
+{
+ //filtering tracking collection using shared hits + quality requirement
+ list<SIDRESULT>::iterator iter1;
+ list<SIDRESULT>::iterator iter2;
+ for (iter1=m_listResults.begin(); iter1!=m_listResults.end(); iter1++)
+ {
+ for (iter2=m_listResults.begin(); iter2!=m_listResults.end(); )
+ {
+ if ( GetSharedHits((*iter1).ListHitID, (*iter2).ListHitID) > MAXSHAREDHITS && iter1!=iter2 )
+ {
+ if ( (*iter1).fChi2 < (*iter2).fChi2 )
+ {
+ //m_listResults.erase(iter1);
+ //--iter1;
+ ++iter2;
+ }
+ else
+ {
+ iter2 = m_listResults.erase(iter2);
+ //--iter2;
+ }
+ }
+ else ++iter2;
+ }
+
+ }
+
+}
+
+Int_t AFP_SIDBasicKalman::GetSharedHits(const vector<Int_t> &HID1, const vector<Int_t> &HID2)
+{
+ Int_t SharedHits = 0;
+
+ vector<Int_t>::const_iterator iter1;
+ vector<Int_t>::const_iterator iter2;
+
+ for (iter1=HID1.begin(); iter1!=HID1.end(); iter1++)
+ {
+ for (iter2=HID2.begin(); iter2!=HID2.end(); iter2++)
+ {
+ if ((*iter1) == (*iter2)) SharedHits++;
+ }
+ }
+
+ return SharedHits;
+}
+
+void AFP_SIDBasicKalman::ClearMatrix()
+{
+ Fk.clear();
+ Ck.clear();
+ CkP.clear();
+ Rk.clear();
+
+ xk.clear();
+ xkP.clear();
+ rk.clear();
+ rkP.clear();
+ mk.clear();
+ chik.clear();
+
+ zk.clear();
+ HID.clear();
+
+ xkS.clear();
+ CkS.clear();
+ chikS.clear();
+
+
+}
+
+void AFP_SIDBasicKalman::HistInitialize()
+{
+ //supporting histograms with binning convention: z, x
+ MyFile = new TFile("PixMap.root","RECREATE");
+ histS1_PixMap = new TH2F("histS1x_PixMap", "", 600, fzMapSID[1][0][0][0]-5., fzMapSID[1][5][335][0]+5., 400, fxMapSID[1][5][335][0]-2., fxMapSID[1][0][0][0]+2.);
+ histS2_PixMap = new TH2F("histS2x_PixMap", "", 600, fzMapSID[0][0][0][0]-5., fzMapSID[0][5][335][0]+5., 400, fxMapSID[0][5][335][0]-2., fxMapSID[0][0][0][0]+2.);
+ histS3_PixMap = new TH2F("histS1y_PixMap", "", 600, fzMapSID[1][0][0][0]-5., fzMapSID[1][5][335][0]+5., 400, fyMapSID[1][0][0][79]-2., fyMapSID[1][0][0][0]+2.);
+ histS4_PixMap = new TH2F("histS2y_PixMap", "", 600, fzMapSID[0][0][0][0]-5., fzMapSID[0][5][335][0]+5., 400, fyMapSID[0][0][0][79]-2., fyMapSID[0][0][0][0]+2.);
+}
+
+void AFP_SIDBasicKalman::HistFinalize()
+{
+
+ MyFile->Write();
+ //MyFile->Close();
+ delete histS1_PixMap;
+ delete histS2_PixMap;
+ delete histS3_PixMap;
+ delete histS4_PixMap;
+ delete MyFile;
+}
+
+
+void AFP_SIDBasicKalman::GetData()
+{
+ MsgStream LogStream(Athena::getMessageSvc(), "AFP_SIDBasicKalman::GetData()");
+ LogStream << MSG::DEBUG << "begin AFP_SIDBasicKalman::GetData()" << endreq;
+
+ /////
+ /////
+
+ LogStream << MSG::DEBUG << "end AFP_SIDBasicKalman::GetData()" << endreq;
+}
diff --git a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDLocReco.cxx b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDLocReco.cxx
index 7d4003433e1aeaf234ac3a04bb3d5743572bd087..95d704444ec4e43c8e753c69346b0b8de7d50cec 100644
--- a/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDLocReco.cxx
+++ b/ForwardDetectors/AFP/AFP_Reconstruction/AFP_LocReco/src/AFP_SIDLocReco.cxx
@@ -2,428 +2,503 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
-#include "AFP_LocReco/AFP_SIDLocReco.h"
-
-#include "AthenaKernel/errorcheck.h"
-#include "AthenaBaseComps/AthMsgStreamMacros.h"
-
-
-AFP_SIDLocReco::AFP_SIDLocReco(const string& name, ISvcLocator* pSvcLocator) :
-AthAlgorithm(name, pSvcLocator)
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::AFP_SIDLocReco");
-
- m_Config.clear();
- m_pGeometry = new AFP_Geometry(&m_Config);
-
-
- // data type using in the local reconstruction
- // for the simulation data the value is 0, for the real data the value is 1. Unset value is -1
- declareProperty("DataType", m_iDataType=0, "data type using in the local reconstruction");
-
- //reconstruction methods properties
- declareProperty("AmpThresh", m_AmpThresh=10.);
-
-
- //reconstruction method selection for TD
- declareProperty("AlgoSID", m_strAlgoSID="SIDBasicKalman");
-
-
- m_vecListAlgoSID.clear();
- declareProperty("ListAlgoSID", m_vecListAlgoSID);
-
-
-
- m_strKeyGeometryForReco = "AFP_GeometryForReco";
- m_strSIDCollectionName = "AFP_SiDigiCollection";
- m_strKeySIDLocRecoEvCollection = "AFP_SIDLocRecoEvCollection";
-
-
- m_pSIDLocRecoEvCollection = NULL;
- m_pSIDLocRecoEvent = NULL;
- m_storeGate = NULL;
-
- m_eventNum = 0;
- m_iEvent = 0;
- m_iRunNum = 0;
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::AFP_SIDLocReco");
-}
-
-AFP_SIDLocReco::~AFP_SIDLocReco()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::~AFP_SIDLocReco");
-
-// if(m_pGeometryReader!=NULL)
-// {
-// delete m_pGeometryReader;
-// m_pGeometryReader = NULL;
-// }
-
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::~AFP_SIDLocReco");
-}
-
-StatusCode AFP_SIDLocReco::initialize()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::initialize()");
-
- StatusCode sc;
- ClearGeometry();
-
- sc = service("StoreGateSvc",m_storeGate);
- if(sc.isFailure())
- {
- ATH_MSG_ERROR("reconstruction: unable to retrieve pointer to StoreGateSvc");
- return sc;
- }
-
- //read geometry
- if(ReadGeometryDetCS())
- {
- ATH_MSG_DEBUG("Geometry loaded successfully");
- }
- else
- {
- ATH_MSG_FATAL("Could not load geometry");
- return StatusCode::FAILURE;
- }
-
- //write AFP_GeometryReader to StoreGate
- /*if (StatusCode::SUCCESS!=service("DetectorStore",m_pDetStore))
- {
- LogStream << MSG::ERROR << "Detector store not found" << endreq;
- return StatusCode::FAILURE;
- }
- sc = m_pDetStore->record(m_pGeometryReader, m_strKeyGeometryForReco);
- if(sc.isFailure())
- {
- LogStream << MSG::ERROR << "m_pGeometryReader: unable to record to StoreGate" << endreq;
- return sc;
- }*/
-
- m_iEvent = 0;
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::initialize()");
- return StatusCode::SUCCESS;
-}
-
-StatusCode AFP_SIDLocReco::execute()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::execute()");
-
- StatusCode sc;
-
- list<SIDHIT> ListSIDHits;
- ListSIDHits.clear();
-
- m_eventNum = 0;
- m_iRunNum = 0;
- const EventInfo* eventInfo;
- sc = m_storeGate->retrieve( eventInfo );
- if (sc.isFailure())
- {
- ATH_MSG_ERROR("AFP_SIDLocReco, Cannot get event info.");
-// return sc;
- }
- else
- {
- // current event number
- m_eventNum = eventInfo->event_ID()->event_number();
- m_iRunNum = eventInfo->event_ID()->run_number();
- }
-
-/*
- // get truth information
- //////////////////////////////////////////
- const McEventCollection* mcTru = 0;
- sc = m_storeGate->retrieve(mcTru,"TruthEvent");
- if(sc.isFailure() || !mcTru){
- LogStream << MSG::DEBUG << "Container "<< "TruthEvent" <<" NOT FOUND !!!!!!!" << endreq;
-// return StatusCode::FAILURE;
- }
-*/
-
-
- sc = RecordSIDCollection();
- if (sc.isFailure())
- {
- ATH_MSG_ERROR("RecordCollection() failed");
- return StatusCode::SUCCESS;
- }
-
- sc = AFPCollectionReading(ListSIDHits);
- if (sc.isFailure())
- {
- ATH_MSG_WARNING("AFP_Collection_Reading Failed");
-// return StatusCode::SUCCESS;
- }
- else
- {
- string strAlgoSID;
- for(unsigned int i=0; i<m_vecListAlgoSID.size(); i++)
- {
- strAlgoSID = m_vecListAlgoSID[i];
-
- //execute SID algorithm
- sc = ExecuteRecoMethod(strAlgoSID, ListSIDHits);
-
- if (sc.isFailure())
- {
- ATH_MSG_FATAL("SID Algorithm " << strAlgoSID << " failure!");
- return sc;
- }
- }
- }
-
-
-
-
-
- m_iEvent++;
- ATH_MSG_DEBUG("end AFP_SIDLocReco::execute()");
- return StatusCode::SUCCESS;
-}
-
-StatusCode AFP_SIDLocReco::finalize()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::finalize()");
-
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::finalize()");
- return StatusCode::SUCCESS;
-}
-
-void AFP_SIDLocReco::ClearGeometry()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::ClearGeometry()");
-
- ///////
- memset(&m_fsSID, 0, sizeof(m_fsSID));
- memset(&m_fxSID, 0, sizeof(m_fxSID));
- memset(&m_fySID, 0, sizeof(m_fySID));
- memset(&m_fzSID, 0, sizeof(m_fzSID));
- ///////
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::ClearGeometry()");
-}
-
-
-StatusCode AFP_SIDLocReco::AFPCollectionReading(list<SIDHIT> &ListSIDHits)
-{
- StatusCode sc;
-
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::AFPCollectionReading()");
-
- SIDHIT SIDHit;
-
- ListSIDHits.clear();
-
-
- const AFP_SiDigiCollection* mcSIDGen = 0;
-
-
- sc = m_storeGate->retrieve(mcSIDGen,m_strSIDCollectionName);
- if(sc.isFailure() || !mcSIDGen)
- {
- return StatusCode::FAILURE;
- }
-
- AFP_SiDigiCollection::const_iterator mcSIDGenBeg = mcSIDGen->begin();
- AFP_SiDigiCollection::const_iterator mcSIDGenEnd = mcSIDGen->end();
-
- for(;mcSIDGenBeg!=mcSIDGenEnd;++mcSIDGenBeg)
- {
- SIDHit.iEvent = m_eventNum;
- SIDHit.fADC = mcSIDGenBeg->m_fADC;
- SIDHit.fTDC = mcSIDGenBeg->m_fTDC;
- SIDHit.nDetectorID = mcSIDGenBeg->m_nDetectorID;
- SIDHit.nStationID = mcSIDGenBeg->m_nStationID;
- SIDHit.nPixelRow = mcSIDGenBeg->m_nPixelRow;
- SIDHit.nPixelCol = mcSIDGenBeg->m_nPixelCol;
-
- ListSIDHits.push_back(SIDHit);
- }
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::AFPCollectionReading()");
-
- return StatusCode::SUCCESS;
-}
-
-
-bool AFP_SIDLocReco::ReadGeometryDetCS()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::ReadGeometryDetCS()");
-
- //////
- StatusCode sc;
-
- for(Int_t nStationID = 0; nStationID < SIDSTATIONID; nStationID++)
- {
- for (Int_t nPlateID = 0; nPlateID < SIDCNT; nPlateID++)
- {
-
- HepGeom::Point3D<double> LocPoint = HepGeom::Point3D<double>(-SID_SENSORXDIM+SID_DEATH_EDGE, -SID_SENSORYDIM+SID_LOWER_EDGE, 0.*CLHEP::mm); //changed! (death edge info from AFP_Geometry)
- HepGeom::Point3D<double> GloPoint = HepGeom::Point3D<double>();
- sc = m_pGeometry->GetPointInSIDSensorGlobalCS(nStationID, nPlateID, LocPoint, GloPoint);
-
- if (sc.isFailure())
- {
- ATH_MSG_WARNING("AFP_Geometry::GetPointInSIDSensorGlobalCS() Failed");
-
- m_fsSID[nStationID][nPlateID] = NOMINALSIDSLOPE;
- m_fxSID[nStationID][nPlateID] = LocPoint.x();
- m_fySID[nStationID][nPlateID] = LocPoint.y();
- m_fzSID[nStationID][nPlateID] = NOMINALSIDSPACING*nPlateID;
- }
-
- else
- {
- m_fsSID[nStationID][nPlateID] = NOMINALSIDSLOPE;
- m_fxSID[nStationID][nPlateID] = GloPoint.x();
- m_fySID[nStationID][nPlateID] = GloPoint.y();
- m_fzSID[nStationID][nPlateID] = GloPoint.z();
-
- ATH_MSG_DEBUG("Global edge position of SID sensor: " <<GloPoint.x()<< "\t" <<GloPoint.y()<< "\t" <<GloPoint.z()<< "\t");
- }
-
- }
- }
- //////
-
- //SaveGeometry();
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::ReadGeometryDetCS()");
-
- return 1;
-}
-
-bool AFP_SIDLocReco::StoreReconstructionGeometry(/*const char* szDataDestination*/)
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::StoreReconstructionGeometry()");
-
- //////
- // (...)
- //////
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::StoreReconstructionGeometry()");
-
- return true;
-}
-
-void AFP_SIDLocReco::SaveGeometry()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::SaveGeometry()");
-
- /////
- // (...)
- /////
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::SaveGeometry()");
-}
-
-StatusCode AFP_SIDLocReco::ExecuteRecoMethod(const string strAlgo, const list<SIDHIT> &ListSIDHits)
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::ExecuteRecoMethod()");
-
- StatusCode sc = StatusCode::SUCCESS;
- SIDRESULT SIDResults;
- list<SIDRESULT> listSIDResults;
- listSIDResults.clear();
-
-
- map<string, int> mapRecoMethods;
- mapRecoMethods.clear();
- mapRecoMethods.insert(pair<string, int>("SIDBasicKalman", 1));
-
-
- switch(mapRecoMethods[strAlgo])
- {
- case 1:
- {
- AFP_SIDBasicKalman* pSIDBasicKalman = new AFP_SIDBasicKalman();
-
- sc = pSIDBasicKalman->Initialize(m_AmpThresh, m_iDataType, ListSIDHits, m_fsSID, m_fxSID, m_fySID, m_fzSID);
- sc = pSIDBasicKalman->Execute();
- sc = pSIDBasicKalman->Finalize(&listSIDResults);
-
- list<SIDRESULT>::const_iterator iter;
- for(iter=listSIDResults.begin(); iter!=listSIDResults.end(); iter++)
- {
- SIDResults.nStationID = (*iter).nStationID;
- SIDResults.x_pos = (*iter).x_pos;
- SIDResults.y_pos = (*iter).y_pos;
- SIDResults.z_pos = (*iter).z_pos;
- SIDResults.x_slope = (*iter).x_slope;
- SIDResults.y_slope = (*iter).y_slope;
- SIDResults.z_slope = (*iter).z_slope;
-
- SIDResults.nHits = (*iter).nHits;
- SIDResults.nHoles = (*iter).nHoles;
- SIDResults.fChi2 = (*iter).fChi2;
-
- if ((SIDResults.nStationID != -1))
- {
-
- m_pSIDLocRecoEvCollection->push_back(new AFP_SIDLocRecoEvent(1, SIDResults.nStationID, SIDResults.x_pos, SIDResults.y_pos, SIDResults.z_pos, SIDResults.x_slope, SIDResults.y_slope, SIDResults.z_slope, SIDResults.nHits, SIDResults.nHoles, SIDResults.fChi2));
-
- }
- }
-
- delete pSIDBasicKalman;
- break;
- }
-
-
-
-
- default:
- {
- ATH_MSG_ERROR("Unable to recognize selected algorithm");
- return StatusCode::FAILURE;
- }
- }
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::ExecuteRecoMethod()");
-
- return sc;
-}
-
-StatusCode AFP_SIDLocReco::RecordSIDCollection()
-{
- ATH_MSG_DEBUG("begin AFP_SIDLocReco::RecordSIDCollection()");
-
- StatusCode sc = StatusCode::SUCCESS;
-
- m_pSIDLocRecoEvCollection = new AFP_SIDLocRecoEvCollection();
- sc = m_storeGate->record(m_pSIDLocRecoEvCollection, m_strKeySIDLocRecoEvCollection);
-
- if (sc.isFailure())
- {
- ATH_MSG_ERROR( m_strAlgoSID << "SID - Could not record the empty LocRecoEv collection in StoreGate");
- return sc;
- }
- else
- {
- ATH_MSG_DEBUG("SID - LocRecoEv collection was recorded in StoreGate");
- }
-
- ATH_MSG_DEBUG("end AFP_SIDLocReco::RecordSIDCollection()");
-
- return sc;
-}
-
-
-void SIDRESULT::clear()
-{
- nStationID=-1;
- x_pos=0.;
- y_pos=0.;
- z_pos=0.;
- x_slope=0.;
- y_slope=0.;
- z_slope=0.;
-
- nHits=0;
- nHoles=0;
- fChi2=0.;
-}
+#include "AFP_LocReco/AFP_SIDLocReco.h"
+
+#include "AthenaKernel/errorcheck.h"
+#include "AthenaBaseComps/AthMsgStreamMacros.h"
+
+#include "AthLinks/ElementLink.h"
+
+#include "xAODForward/AFPSiHit.h"
+#include "xAODForward/AFPSiHitContainer.h"
+#include "xAODForward/AFPTrack.h"
+#include "xAODForward/AFPTrackContainer.h"
+#include "xAODForward/AFPTrackAuxContainer.h"
+
+
+AFP_SIDLocReco::AFP_SIDLocReco(const string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name, pSvcLocator)
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::AFP_SIDLocReco");
+
+ m_Config.clear();
+ m_pGeometry = new AFP_Geometry(&m_Config);
+
+
+ // data type using in the local reconstruction
+ // for the simulation data the value is 0, for the real data the value is 1. Unset value is -1
+ declareProperty("DataType", m_iDataType=0, "data type using in the local reconstruction");
+
+ //reconstruction methods properties
+ declareProperty("AmpThresh", m_AmpThresh=10.);
+
+
+ //reconstruction method selection for TD
+ declareProperty("AlgoSID", m_strAlgoSID="SIDBasicKalman");
+
+
+ m_vecListAlgoSID.clear();
+ declareProperty("ListAlgoSID", m_vecListAlgoSID);
+
+
+ m_strKeyGeometryForReco = "AFP_GeometryForReco";
+ m_strSIDCollectionName = "AFPSiHitContainer";
+ m_strKeySIDLocRecoEvCollection = "AFPTrackContainer";
+
+ m_pSIDLocRecoEvCollection = NULL;
+ m_pSIDLocRecoEvent = NULL;
+ m_storeGate = NULL;
+
+ m_eventNum = 0;
+ m_iEvent = 0;
+ m_iRunNum = 0;
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::AFP_SIDLocReco");
+}
+
+AFP_SIDLocReco::~AFP_SIDLocReco()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::~AFP_SIDLocReco");
+
+ // if(m_pGeometryReader!=NULL)
+ // {
+ // delete m_pGeometryReader;
+ // m_pGeometryReader = NULL;
+ // }
+
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::~AFP_SIDLocReco");
+}
+
+StatusCode AFP_SIDLocReco::initialize()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::initialize()");
+
+ StatusCode sc;
+ ClearGeometry();
+
+ sc = service("StoreGateSvc",m_storeGate);
+ if(sc.isFailure())
+ {
+ ATH_MSG_WARNING("reconstruction: unable to retrieve pointer to StoreGateSvc");
+ return sc;
+ }
+
+ //read geometry
+ if(ReadGeometryDetCS())
+ {
+ ATH_MSG_DEBUG("Geometry loaded successfully");
+ }
+ else
+ {
+ ATH_MSG_WARNING("Could not load geometry");
+ return StatusCode::SUCCESS;
+ }
+
+ //write AFP_GeometryReader to StoreGate
+ /*if (StatusCode::SUCCESS!=service("DetectorStore",m_pDetStore))
+ {
+ LogStream << MSG::ERROR << "Detector store not found" << endreq;
+ return StatusCode::FAILURE;
+ }
+ sc = m_pDetStore->record(m_pGeometryReader, m_strKeyGeometryForReco);
+ if(sc.isFailure())
+ {
+ LogStream << MSG::ERROR << "m_pGeometryReader: unable to record to StoreGate" << endreq;
+ return sc;
+ }*/
+
+ m_iEvent = 0;
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::initialize()");
+ return StatusCode::SUCCESS;
+}
+
+StatusCode AFP_SIDLocReco::execute()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::execute()");
+
+ StatusCode sc;
+
+ list<SIDHIT> ListSIDHits;
+ ListSIDHits.clear();
+
+ m_eventNum = 0;
+ m_iRunNum = 0;
+ const EventInfo* eventInfo;
+ sc = m_storeGate->retrieve( eventInfo );
+ if (sc.isFailure())
+ {
+ ATH_MSG_WARNING("AFP_SIDLocReco, Cannot get event info.");
+ return StatusCode::SUCCESS;
+ }
+ else
+ {
+ // current event number
+ m_eventNum = eventInfo->event_ID()->event_number();
+ m_iRunNum = eventInfo->event_ID()->run_number();
+ }
+
+
+ // // get truth information
+ // //////////////////////////////////////////
+ // const McEventCollection* mcTru = 0;
+ // sc = m_storeGate->retrieve(mcTru,"TruthEvent");
+ // if(sc.isFailure() || !mcTru){
+ // LogStream << MSG::DEBUG << "Container "<< "TruthEvent" <<" NOT FOUND !!!!!!!" << endreq;
+ // // return StatusCode::FAILURE;
+ // }
+
+
+
+ sc = RecordSIDCollection();
+ if (sc.isFailure())
+ {
+ ATH_MSG_WARNING("RecordCollection() failed");
+ return StatusCode::SUCCESS;
+ }
+
+ // init xAOD
+ xAOD::AFPTrackContainer* trackContainer = new xAOD::AFPTrackContainer();
+ CHECK( evtStore()->record(trackContainer, m_strKeySIDLocRecoEvCollection) );
+ xAOD::AFPTrackAuxContainer* trackAuxContainer = new xAOD::AFPTrackAuxContainer();
+ CHECK( evtStore()->record(trackAuxContainer, m_strKeySIDLocRecoEvCollection + "Aux.") );
+ trackContainer->setStore(trackAuxContainer);
+
+ sc = AFPCollectionReading(ListSIDHits);
+ if (sc.isFailure()) {
+ ATH_MSG_WARNING("AFP_Collection_Reading Failed");
+ // return StatusCode::SUCCESS;
+ }
+ else {
+ string strAlgoSID;
+ for(unsigned int i=0; i<m_vecListAlgoSID.size(); i++)
+ {
+ strAlgoSID = m_vecListAlgoSID[i];
+
+ //execute SID algorithm
+ sc = ExecuteRecoMethod(strAlgoSID, ListSIDHits, trackContainer);
+
+ if (sc.isFailure())
+ {
+ ATH_MSG_WARNING("SID Algorithm " << strAlgoSID << " failure!");
+ return StatusCode::SUCCESS;
+ }
+ }
+ }
+
+
+ m_iEvent++;
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::execute()");
+ return StatusCode::SUCCESS;
+}
+
+StatusCode AFP_SIDLocReco::finalize()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::finalize()");
+
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::finalize()");
+ return StatusCode::SUCCESS;
+}
+
+void AFP_SIDLocReco::ClearGeometry()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::ClearGeometry()");
+
+ ///////
+ memset(&m_fsSID, 0, sizeof(m_fsSID));
+ memset(&m_fxSID, 0, sizeof(m_fxSID));
+ memset(&m_fySID, 0, sizeof(m_fySID));
+ memset(&m_fzSID, 0, sizeof(m_fzSID));
+ ///////
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::ClearGeometry()");
+}
+
+
+StatusCode AFP_SIDLocReco::AFPCollectionReading(list<SIDHIT> &ListSIDHits)
+{
+ StatusCode sc;
+
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::AFPCollectionReading()");
+
+ SIDHIT SIDHit;
+
+ ListSIDHits.clear();
+
+ const xAOD::AFPSiHitContainer* siHitContainer = 0;
+ sc = m_storeGate->retrieve(siHitContainer, m_strSIDCollectionName);
+ if(sc.isFailure() || !siHitContainer)
+ {
+ return StatusCode::SUCCESS;
+ }
+
+ xAOD::AFPSiHitContainer::const_iterator hitIter = siHitContainer->begin();
+ xAOD::AFPSiHitContainer::const_iterator mcSIDGenEnd = siHitContainer->end();
+
+ for(;hitIter!=mcSIDGenEnd;++hitIter) {
+ xAOD::AFPSiHit* hit = *hitIter;
+ SIDHit.iEvent = m_eventNum;
+ SIDHit.fADC = hit->depositedCharge();
+ SIDHit.fTDC = 0;
+ SIDHit.nDetectorID = hit->pixelLayerID();
+ SIDHit.nStationID = hit->stationID();
+ SIDHit.nPixelRow = hit->pixelVertID();
+ SIDHit.nPixelCol = hit->pixelHorizID();
+
+
+ ListSIDHits.push_back(SIDHit);
+ }
+
+ // const AFP_SiDigiCollection* mcSIDGen = 0;
+ // sc = m_storeGate->retrieve(mcSIDGen,m_strSIDCollectionName);
+ // if(sc.isFailure() || !mcSIDGen)
+
+ // AFP_SiDigiCollection::const_iterator mcSIDGenBeg = mcSIDGen->begin();
+ // AFP_SiDigiCollection::const_iterator mcSIDGenEnd = mcSIDGen->end();
+
+ // for(;mcSIDGenBeg!=mcSIDGenEnd;++mcSIDGenBeg) {
+ // SIDHit.iEvent = m_eventNum;
+ // SIDHit.fADC = mcSIDGenBeg->m_fADC;
+ // SIDHit.fTDC = mcSIDGenBeg->m_fTDC;
+ // SIDHit.nDetectorID = mcSIDGenBeg->m_nDetectorID;
+ // SIDHit.nStationID = mcSIDGenBeg->m_nStationID;
+ // SIDHit.nPixelRow= mcSIDGenBeg->m_nPixelRow;
+ // SIDHit.nPixelCol= mcSIDGenBeg->m_nPixelCol;
+
+ // ListSIDHits.push_back(SIDHit);
+ // }
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::AFPCollectionReading()");
+
+ return StatusCode::SUCCESS;
+}
+
+bool AFP_SIDLocReco::ReadGeometryDetCS()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::ReadGeometryDetCS()");
+
+ //////
+ StatusCode sc;
+
+ for(Int_t nStationID = 0; nStationID < SIDSTATIONID; nStationID++)
+ {
+ for (Int_t nPlateID = 0; nPlateID < SIDCNT; nPlateID++)
+ {
+
+ HepGeom::Point3D<double> LocPoint = HepGeom::Point3D<double>(-SID_SENSORXDIM+SID_DEATH_EDGE, -SID_SENSORYDIM+SID_LOWER_EDGE, 0.*CLHEP::mm); //changed! (death edge info from AFP_Geometry)
+ HepGeom::Point3D<double> GloPoint = HepGeom::Point3D<double>();
+ sc = m_pGeometry->GetPointInSIDSensorGlobalCS(nStationID, nPlateID, LocPoint, GloPoint);
+
+ if (sc.isFailure())
+ {
+ ATH_MSG_WARNING("AFP_Geometry::GetPointInSIDSensorGlobalCS() Failed");
+
+ m_fsSID[nStationID][nPlateID] = SID_NOMINALSLOPE;
+ m_fxSID[nStationID][nPlateID] = LocPoint.x();
+ m_fySID[nStationID][nPlateID] = LocPoint.y();
+ m_fzSID[nStationID][nPlateID] = SID_NOMINALSPACING*nPlateID;
+ }
+
+ else
+ {
+ m_fsSID[nStationID][nPlateID] = SID_NOMINALSLOPE;
+ m_fxSID[nStationID][nPlateID] = GloPoint.x();
+ m_fySID[nStationID][nPlateID] = GloPoint.y();
+ m_fzSID[nStationID][nPlateID] = GloPoint.z();
+
+ ATH_MSG_DEBUG("Global edge position of SID sensor: " <<GloPoint.x()<< "\t" <<GloPoint.y()<< "\t" <<GloPoint.z()<< "\t");
+ }
+
+ }
+ }
+ //////
+
+ //SaveGeometry();
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::ReadGeometryDetCS()");
+
+ return 1;
+}
+
+bool AFP_SIDLocReco::StoreReconstructionGeometry(/*const char* szDataDestination*/)
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::StoreReconstructionGeometry()");
+
+ //////
+ // (...)
+ //////
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::StoreReconstructionGeometry()");
+
+ return true;
+}
+
+void AFP_SIDLocReco::SaveGeometry()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::SaveGeometry()");
+
+ /////
+ // (...)
+ /////
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::SaveGeometry()");
+}
+
+StatusCode AFP_SIDLocReco::ExecuteRecoMethod(const string strAlgo, const list<SIDHIT> &ListSIDHits, xAOD::AFPTrackContainer* resultContainer)
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::ExecuteRecoMethod()");
+
+ if (!resultContainer) {
+ ATH_MSG_WARNING ("Null pointer given for a result tracks container.");
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode sc = StatusCode::SUCCESS;
+ SIDRESULT SIDResults;
+ list<SIDRESULT> listSIDResults;
+ listSIDResults.clear();
+
+ map<string, int> mapRecoMethods;
+ mapRecoMethods.clear();
+ mapRecoMethods.insert(pair<string, int>("SIDBasicKalman", 1));
+
+ switch(mapRecoMethods[strAlgo])
+ {
+ case 1:
+ {
+ const xAOD::AFPSiHitContainer* siHitContainer = 0;
+ sc = m_storeGate->retrieve(siHitContainer, m_strSIDCollectionName);
+ if(sc.isFailure() || !siHitContainer) {
+ return StatusCode::SUCCESS;
+ }
+
+
+ AFP_SIDBasicKalman* pSIDBasicKalman = new AFP_SIDBasicKalman();
+
+ sc = pSIDBasicKalman->Initialize(m_AmpThresh, m_iDataType, ListSIDHits, m_fsSID, m_fxSID, m_fySID, m_fzSID);
+ sc = pSIDBasicKalman->Execute();
+ sc = pSIDBasicKalman->Finalize(&listSIDResults);
+
+ list<SIDRESULT>::const_iterator iter;
+ for(iter=listSIDResults.begin(); iter!=listSIDResults.end(); iter++) {
+ if (iter->nStationID != -1) {
+ xAOD::AFPTrack* track = new xAOD::AFPTrack;
+ resultContainer->push_back(track);
+
+ track->setStationID(iter->nStationID);
+ track->setXLocal(iter->x_pos);
+ track->setYLocal(iter->y_pos);
+ track->setZLocal(iter->z_pos);
+ track->setXSlope(iter->x_slope);
+ track->setYSlope(iter->y_slope);
+ // track->setz_slope iter->z_slope;
+ track->setNHits(iter->nHits);
+ track->setNHoles(iter->nHoles);
+ track->setChi2(iter->fChi2);
+
+ ATH_MSG_DEBUG("Track reconstructed with "<<iter->ListHitID.size()<<" hits. In hits container there are in total "<<siHitContainer->size()<<" hits.");
+
+ std::vector<int>::const_iterator end = iter->ListHitID.end();
+ for (std::vector<int>::const_iterator hitIter = iter->ListHitID.begin(); hitIter != end; ++hitIter) {
+ const int pixelRow = (*hitIter)%100;
+ const int pixelCol = (( (*hitIter)/100) % 1000);
+ const int pixelLayer = ( (*hitIter)/100000) % 10;
+ const int pixelStation = ( (*hitIter)/1000000) % 10;
+ int result = 0;
+ auto endHits = siHitContainer->end();
+ for (auto origHitIter = siHitContainer->begin(); origHitIter != endHits; ++origHitIter) {
+ if (
+ (*origHitIter)->pixelVertID() == pixelRow
+ && (*origHitIter)->pixelHorizID() == pixelCol
+ && (*origHitIter)->pixelLayerID() == pixelLayer
+ && (*origHitIter)->stationID() == pixelStation
+ )
+ break;
+
+ result++;
+ }
+
+ // check if the hit was found
+ if (result < siHitContainer->size()) {
+ ATH_MSG_DEBUG("To the list of hits in a track adding hit "<<result<<"/"<<siHitContainer->size()<<".");
+
+ ElementLink< xAOD::AFPSiHitContainer >* hitLink = new ElementLink< xAOD::AFPSiHitContainer >;
+ hitLink->toIndexedElement(*siHitContainer, result);
+ track->addHit(*hitLink);
+
+ ElementLink< xAOD::AFPTrackContainer >* trackLink = new ElementLink< xAOD::AFPTrackContainer >;
+ trackLink->toIndexedElement(*resultContainer, resultContainer->size() -1);
+ xAOD::AFPSiHit * hit = const_cast<xAOD::AFPSiHit*> (siHitContainer->at(result));
+ hit->addTrackLink(*trackLink);
+ }
+ else
+ ATH_MSG_WARNING("Track hit not found in hits list. HitID: "<<(*hitIter)
+ <<" station: "<<pixelStation
+ <<" layer: "<<pixelLayer
+ <<" row: "<<pixelRow
+ <<" col: "<<pixelCol
+ <<" dataType: "<<m_iDataType
+ );
+ }
+ }
+ }
+
+ delete pSIDBasicKalman;
+ break;
+ }
+
+
+
+
+ default:
+ {
+ ATH_MSG_WARNING("Unable to recognize selected algorithm");
+ return StatusCode::SUCCESS;
+ }
+ }
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::ExecuteRecoMethod()");
+
+ return sc;
+}
+
+StatusCode AFP_SIDLocReco::RecordSIDCollection()
+{
+ ATH_MSG_DEBUG("begin AFP_SIDLocReco::RecordSIDCollection()");
+
+ StatusCode sc = StatusCode::SUCCESS;
+
+ m_pSIDLocRecoEvCollection = new AFP_SIDLocRecoEvCollection();
+ sc = m_storeGate->record(m_pSIDLocRecoEvCollection, m_strKeySIDLocRecoEvCollection);
+
+ if (sc.isFailure())
+ {
+ ATH_MSG_WARNING( m_strAlgoSID << "SID - Could not record the empty LocRecoEv collection in StoreGate");
+ return StatusCode::SUCCESS;
+ }
+ else
+ {
+ ATH_MSG_DEBUG("SID - LocRecoEv collection was recorded in StoreGate");
+ }
+
+ ATH_MSG_DEBUG("end AFP_SIDLocReco::RecordSIDCollection()");
+
+ return sc;
+}
+
+
+void SIDRESULT::clear()
+{
+ nStationID=-1;
+ x_pos=0.;
+ y_pos=0.;
+ z_pos=0.;
+ x_slope=0.;
+ y_slope=0.;
+ z_slope=0.;
+
+ nHits=0;
+ nHoles=0;
+ fChi2=0.;
+}