diff --git a/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKDataMatrixFromMCAlgo.h b/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKDataMatrixFromMCAlgo.h
new file mode 100644
index 0000000000000000000000000000000000000000..7da2c085c50be0e41a400529064a4257d8f6b0be
--- /dev/null
+++ b/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKDataMatrixFromMCAlgo.h
@@ -0,0 +1,43 @@
+#ifndef FTKDataMatrixFromMCAlgo_h
+#define FTKDataMatrixFromMCAlgo_h
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "TrigFTKSim/FTKLogging.h"
+#include "TrigFTKBankGen/FTKGhostHitCalculator.h"
+
+#include <TFile.h>
+#include <string>
+
+class FTKDataMatrixFromMCAlgo : public AthAlgorithm, public FTKLogger {
+public:
+ FTKDataMatrixFromMCAlgo(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~FTKDataMatrixFromMCAlgo ();
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+ protected:
+ virtual void PostMessage(void); // divert FTKLogger messages to Athena
+ std::string m_MCmatrixFileName;
+ std::string m_DATAmatrixFileName;
+ std::string m_debugFileName;
+ std::string m_MCmodulePositionFileName;
+ std::string m_DATAmodulePositionFileName;
+ std::string m_pmap_path;
+ std::string m_rmap_path;
+ int m_HWMODEID;
+ int m_tower;
+ int m_invertIBLphiMatrix;
+ int m_invertIBLphiData;
+ std::string m_modulelut_path;
+ double m_matrixVtxX,m_matrixVtxY,m_matrixVtxZ;
+ double m_dataVtxX,m_dataVtxY,m_dataVtxZ;
+
+ FTKRegionMap *m_rmap;
+ FTKGhostHitCalculator *m_MCmodulePositionCalculator;
+ FTKGhostHitCalculator *m_DATAmodulePositionCalculator;
+ TFile *m_MCmatrixFile; // input: FTK matrix file for MC geometry
+ TFile *m_DATAmatrixFile; // output: FTK matrix file for data geometry
+ TFile *m_debugFile;
+};
+
+#endif
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKGhostHitCalculator.h b/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKGhostHitCalculator.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0b66a06d461124085a94c381cd8fc30631178ef
--- /dev/null
+++ b/Trigger/TrigFTK/TrigFTKBankGen/TrigFTKBankGen/FTKGhostHitCalculator.h
@@ -0,0 +1,143 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef FTKGHOSTHITCALCULATOR_H
+#define FTKGHOSTHITCALCULATOR_H
+
+/* class FTKGhostHitCalculator
+ *
+ * extrapolate truth tracks to dead modules and estimate ghost hit positions
+ *
+ * the mapping of modules to FTK planes and towers is defined by
+ * the method setFTKGeometry()
+ *
+ * the module positions are read from a TTree loadModuleGeometry()
+ * in that three, dead modules are defined
+ * overwrite the dead module positions using
+ * clearBadModuleList() and addBadModules()
+ *
+ * for a given module ID, check whether the module is on the "bad" list
+ * isBad()
+ *
+ * for a given truth track, module positions are added
+ * using the method addGhostHit()
+ *
+ * Author: Stefan Schmitt, DESY
+ * Date: 2017/06/28
+ *
+ */
+
+#include "TrigFTKSim/FTKLogging.h"
+#include <TVector3.h>
+#include <TVector2.h>
+#include <vector>
+#include <map>
+#include <set>
+
+class FTKRegionMap;
+class FTKTruthTrack;
+class TDirectory;
+class FTKHit;
+
+struct FTK_ModuleGeometry {
+ // geometry data of a single module
+ int m_isPixel;
+ int m_id;
+ int m_plane;
+ int m_hitSector;
+ TVector3 m_center;
+ TVector3 m_phiAxis;
+ TVector3 m_etaAxis;
+ TVector2 m_pitch;
+ TVector2 m_width;
+ // extrapolation method
+ bool extrapolate(double rho,double cotTh,double z0,double d0,
+ double sinPhi0,double cosPhi0,TVector3 const *trackVertex,
+ double *xLocPtr,TVector2 *xPtr,double *sPtr) const;
+};
+
+class FTK_ModuleLocator {
+ public:
+ // locate a module
+ //
+ // for barrel geometry,
+ // determine R[min,max] order all modules by Z for fast lookup
+ // for endcap geometry,
+ // determine Z[min,max] order all modules by R for fast lookup
+ //
+ void insert(FTK_ModuleGeometry const *module);
+ FTKHit *locate(FTKTruthTrack const &track,FTK_ModuleGeometry *module=0) const;
+ void print(void) const;
+ protected:
+ struct FTK_ModuleSet {
+ // barrel or disk geometry
+ int m_isBarrel;
+ // m_x: Rmin,Rmax (barrel) Zmin,Zmax (endcap)
+ double m_x[2];
+ // m_dy: maximum module width in Z (Barrel), in R (endcap)
+ double m_dy;
+ // m_modules: modules orderd by Zavg (Barrel), Ravg (Endcap)
+ std::multimap<double,FTK_ModuleGeometry const *> m_modules;
+ // method to intersect track with modules
+ // if a module matches, the intersection and the module
+ // are added to the output
+ // the candidates are ordered by the arc length
+ void getCandidates(double rho,double cotTh,double z0,double d0,
+ double sinPhi0,double cosPhi0,std::multimap
+ <double,std::pair<TVector2,FTK_ModuleGeometry
+ const *> > &candidate) const;
+ };
+ // barrel,endcap have different module sets
+ // also, if R(barrel) or z(endcap) differs too much, start a new set
+ std::vector<FTK_ModuleSet> m_moduleSetVector;
+};
+
+class FTKGhostHitCalculator : public FTKLogging {
+ public:
+ FTKGhostHitCalculator(const FTKRegionMap *rmap,
+ std::string const &name="FTKGhostHitCalculator");
+ int loadModuleGeometry(std::string const &name,bool markAllBad,bool invertIBLphi);
+ virtual int loadModuleGeometry(TDirectory *source,bool markAllBad,bool invertIBLphi);
+
+ virtual void clearBadModuleList(void);
+ int addBadModules(std::string const &name);
+ int addBadModules(std::istream &in);
+
+ virtual FTKHit *addGhostHit( FTKTruthTrack const &track,int plane,
+ FTK_ModuleGeometry *module=0 ) const;
+ virtual bool isBad(int plane,int moduleID) const;
+
+ // returns:
+ // 0: all ok
+ // 1: bad plane number
+ // 2: bad module id
+ // 3: extrapolation failed
+ virtual int extrapolateTo
+ (size_t plane,int moduleID,double rho,double cotTh,double z0,double d0,
+ double sinPhi0,double cosPhi0,TVector3 const *trackVertex,
+ double *x) const;
+ virtual FTK_ModuleGeometry const *findModule(size_t plane,int moduleID) const;
+ protected:
+ int addBadPixelModules(std::set<int> const &bad);
+ int addBadSCTModules(std::set<int> const &bad);
+
+ void updateCalculator(void);
+
+ FTKRegionMap const *m_rmap;
+
+ // for each module ID give the plane number
+ std::map<int,int> m_planeFromPIXid,m_plane_From_SCTid;
+
+ // bad modules
+ std::set<int> m_badPixel,m_badSCT;
+
+ // give geometry data by plane and module ID
+ // the vector index maps the modules to FTK planes
+ // the map index is the module id
+ std::vector<std::map<int,FTK_ModuleGeometry> > m_geometryData;
+ //
+ // each plane has its own ModuleLocator
+ std::vector<FTK_ModuleLocator> m_moduleLocatorByPlane;
+};
+
+#endif
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKConstGenAlgo.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKConstGenAlgo.cxx
index 39664e44eebd36834d7be2a12e9cafde78c72676..4ac5a195897f6471f9dfd35f76599e0292bd0049 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKConstGenAlgo.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKConstGenAlgo.cxx
@@ -1504,8 +1504,8 @@ void FTKConstGenAlgo::extract_1stStage()
double *tmpxZ;
double *tmpcovx;
- int ndim8 = 11;
- int ndim2_8 = 11*11;
+ int ndim8 = m_pmap_8L->getTotalDim(); //11;
+ int ndim2_8 = ndim8*ndim8; //11*11;
int nth_dim8 = 0;
int nth_dim8_2 = 0;
int idx_cov8 = 0;
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKDataMatrixFromMCAlgo.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKDataMatrixFromMCAlgo.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5674d3d20e46bc459a9bd034e07b83da19d2f8e3
--- /dev/null
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKDataMatrixFromMCAlgo.cxx
@@ -0,0 +1,696 @@
+#include "TrigFTKBankGen/FTKDataMatrixFromMCAlgo.h"
+#include "TrigFTKSim/FTKRegionMap.h"
+#include "TrigFTKSim/FTKSetup.h"
+#include <TTree.h>
+#include <TH2F.h>
+#include <TVectorD.h>
+#include <TMatrixD.h>
+#include <TMatrixDSymEigen.h>
+#include <TError.h>
+#include <vector>
+
+using namespace std;
+
+#define INCLUDE_SHORT_VECTOR
+
+FTKDataMatrixFromMCAlgo::FTKDataMatrixFromMCAlgo
+(const std::string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name,pSvcLocator),
+ m_HWMODEID(2),
+ m_tower(-1),
+ m_invertIBLphiMatrix(1),
+ m_invertIBLphiData(1),
+ m_matrixVtxX(0),
+ m_matrixVtxY(0),
+ m_matrixVtxZ(0),
+ m_rmap(0),
+ m_MCmodulePositionCalculator(0),
+ m_DATAmodulePositionCalculator(0),
+ m_MCmatrixFile(0),
+ m_DATAmatrixFile(0),
+ m_debugFile(0)
+{
+ SetLogger(this);
+ declareProperty("MCmatrixFileName",m_MCmatrixFileName);
+ declareProperty("DATAmatrixFileName",m_DATAmatrixFileName);
+ declareProperty("debugFileName",m_debugFileName);
+ declareProperty("MCmodulePositionFileName",m_MCmodulePositionFileName);
+ declareProperty("DATAmodulePositionFileName",m_DATAmodulePositionFileName);
+ declareProperty("pmap_path", m_pmap_path);
+ declareProperty("rmap_path",m_rmap_path);
+ declareProperty("hwmodeid", m_HWMODEID);
+ declareProperty("tower", m_tower);
+ /* declareProperty("IBLMode",m_IBLMode);
+ declareProperty("FixEndcapL0",m_fixEndcapL0);
+ declareProperty("ITkMode",m_ITkMode); */
+ declareProperty("ModuleLUTPath", m_modulelut_path);
+ declareProperty("matrixVtxX",m_matrixVtxX);
+ declareProperty("matrixVtxY",m_matrixVtxY);
+ declareProperty("matrixVtxZ",m_matrixVtxZ);
+ declareProperty("dataVtxX",m_dataVtxX);
+ declareProperty("dataVtxY",m_dataVtxY);
+ declareProperty("dataVtxZ",m_dataVtxZ);
+ declareProperty("invertIBLphiMatrix",m_invertIBLphiMatrix);
+ declareProperty("invertIBLphiData",m_invertIBLphiData);
+}
+
+FTKDataMatrixFromMCAlgo::~FTKDataMatrixFromMCAlgo() {
+ if(m_DATAmatrixFile) delete m_DATAmatrixFile;
+ if(m_MCmatrixFile) delete m_MCmatrixFile;
+ if(m_debugFile) delete m_debugFile;
+ if(m_MCmodulePositionCalculator) delete m_MCmodulePositionCalculator;
+ if(m_DATAmodulePositionCalculator) delete m_DATAmodulePositionCalculator;
+ if(m_rmap) delete m_rmap;
+}
+
+void FTKDataMatrixFromMCAlgo::PostMessage(void) {
+ if (FTKLogger::m_type==0) ATH_MSG_FATAL(m_buffer->str());
+ else if(FTKLogger::m_type==1) ATH_MSG_ERROR(m_buffer->str());
+ else if(FTKLogger::m_type==2) ATH_MSG_WARNING(m_buffer->str());
+ else if(FTKLogger::m_type==3) ATH_MSG_INFO(m_buffer->str());
+ else if(FTKLogger::m_type==4) ATH_MSG_DEBUG(m_buffer->str());
+}
+
+StatusCode FTKDataMatrixFromMCAlgo::initialize() {
+ // open root files
+ int error=0;
+ FTKSetup &ftkset = FTKSetup::getFTKSetup();
+ ftkset.setHWModeSS(m_HWMODEID);
+ /*ftkset.setIBLMode(m_IBLMode);
+ ftkset.setfixEndcapL0(m_fixEndcapL0);
+ ftkset.setITkMode(m_ITkMode); */
+ if(m_pmap_path.empty()) {
+ ATH_MSG_FATAL("no plane map specified");
+ }
+ if(m_rmap_path.empty()) {
+ ATH_MSG_FATAL("no region map specified");
+ }
+
+ ATH_MSG_INFO("Loading plane map: "+m_pmap_path);
+ FTKPlaneMap* pmap = new FTKPlaneMap(m_pmap_path.c_str());
+ if (!(*pmap)) {
+ ATH_MSG_FATAL("Error using plane map: "+ m_pmap_path);
+ }
+
+ ATH_MSG_INFO("Loading region map: "+m_rmap_path);
+ m_rmap = new FTKRegionMap(pmap, m_rmap_path.c_str());
+
+ if(m_HWMODEID==2) {
+ if(m_modulelut_path.empty()) {
+ ATH_MSG_FATAL("A module LUT is required with HWMODEID==2");
+ return StatusCode::FAILURE;
+ } else {
+ m_rmap->loadModuleIDLUT(m_modulelut_path.c_str());
+ }
+ }
+
+ m_MCmatrixFile=new TFile(m_MCmatrixFileName.c_str());
+ m_DATAmatrixFile=new TFile(m_DATAmatrixFileName.c_str(),"recreate");
+ m_MCmodulePositionCalculator=
+ new FTKGhostHitCalculator(m_rmap,"FTKGhostHitCalculator_MC");
+ m_DATAmodulePositionCalculator=
+ new FTKGhostHitCalculator(m_rmap,"FTKGhostHitCalculator_data");
+ if(!m_MCmatrixFile->IsOpen()) {
+ ATH_MSG_FATAL("Can not open input matrix file \""+
+ m_MCmatrixFileName+"\" for reading");
+ error++;
+ } else {
+ ATH_MSG_INFO("opening input matrix file \""+m_MCmatrixFileName+"\"");
+ }
+ if(!m_DATAmatrixFile->IsOpen()) {
+ ATH_MSG_FATAL("Can not open DATA matrix file \""+
+ m_DATAmatrixFileName+"\" for writing");
+ error++;
+ } else {
+ ATH_MSG_INFO("opening output matrix file \""+m_DATAmatrixFileName+"\"");
+ }
+ if(!m_MCmodulePositionCalculator->loadModuleGeometry
+ (m_MCmodulePositionFileName,true,m_invertIBLphiMatrix)) {
+ ATH_MSG_FATAL("Can not load MC module position file \""+
+ m_MCmodulePositionFileName+"\" for reading");
+ error++;
+ }
+ if(!m_DATAmodulePositionCalculator->loadModuleGeometry
+ (m_DATAmodulePositionFileName,true,m_invertIBLphiData)) {
+ ATH_MSG_FATAL("Can not open DATA module position file \""+
+ m_DATAmodulePositionFileName+"\" for reading");
+ error++;
+ }
+ if(!m_debugFileName.empty()) {
+ m_debugFile=new TFile(m_debugFileName.c_str(),"recreate");
+ if(m_debugFile->IsOpen()) {
+ ATH_MSG_INFO("Debugging information is stored in "+
+ m_debugFileName);
+ } else {
+ ATH_MSG_WARNING("File \""+m_debugFileName+
+ "\"for storing debug information can not be opened");
+ delete m_debugFile;
+ m_debugFile=0;
+ }
+ }
+ if(error) {
+ return StatusCode::FAILURE;
+ } else {
+ return StatusCode::SUCCESS;
+ }
+}
+
+StatusCode FTKDataMatrixFromMCAlgo::execute() {
+ if(!m_DATAmatrixFile) return StatusCode::FAILURE;
+ // TTRee with debugging information
+ TVector3 matrixVtx(m_matrixVtxX,m_matrixVtxY,m_matrixVtxZ);
+ TVector3 dataVtx(m_dataVtxX,m_dataVtxY,m_dataVtxZ);
+ ATH_MSG_INFO("matrix file assumed to be relative to "+
+ TString::Format("(%lf,%lf,%lf)",
+ matrixVtx(0),matrixVtx(1),matrixVtx(2)));
+ ATH_MSG_INFO("output matrix file will be relative to "+
+ TString::Format("(%lf,%lf,%lf)",
+ dataVtx(0),dataVtx(1),dataVtx(2)));
+#ifdef WRITE_DEBUG_TREE
+ TTree *debugPix=0,*debugSCT=0;
+ int debug_tower,debug_hash,debug_plane;
+ float debug_track[5],debug_hit[2],debug_ghost[2],debug_data[2];
+#endif
+ vector<TH2 *> debug_coordHM(16),debug_coordMD(16),debug_resCurv(16);
+ if(m_debugFile) {
+ m_debugFile->cd();
+#ifdef WRITE_DEBUG_TREE
+ debugPix=new TTree("debugPix","debugPix");
+ debugSCT=new TTree("debugSCT","debugSCT");
+ debugPix->Branch("tower",&debug_tower,"tower/I");
+ debugPix->Branch("plane",&debug_plane,"plane/I");
+ debugPix->Branch("hash",&debug_hash,"hash/I");
+ debugPix->Branch("track",debug_track,"track[5]/F");
+ debugPix->Branch("hit",debug_hit,"hit[2]/F");
+ debugPix->Branch("ghost",debug_ghost,"ghost[2]/F");
+ debugPix->Branch("data",debug_data,"data[2]/F");
+ debugSCT->Branch("tower",&debug_tower,"tower/I");
+ debugSCT->Branch("plane",&debug_plane,"plane/I");
+ debugSCT->Branch("hash",&debug_hash,"hash/I");
+ debugSCT->Branch("track",debug_track,"track[5]/F");
+ debugSCT->Branch("hit",debug_hit,"hit/F");
+ debugSCT->Branch("ghost",debug_ghost,"ghost/F");
+ debugSCT->Branch("data",debug_data,"data/F");
+#endif
+ for(size_t i=0;i<debug_coordHM.size();i++) {
+ double w=800;
+ if(i<8) {
+ if(i&1) w=160;
+ else w=360.;
+ }
+ double x0=-0.1*w;
+ double x1=1.1*w;
+ debug_coordHM[i]=new TH2F(TString::Format("coordHM%d",(int)i),
+ ";hit;matrix",50,x0,x1,50,x0,x1);
+ debug_coordMD[i]=new TH2F(TString::Format("coordMD%d",(int)i),
+ ";matrix;output",50,x0,x1,50,x0,x1);
+ debug_resCurv[i]=new TH2F(TString::Format("resCurv%d",(int)i),
+ ";matrix;output",50,-0.001,0.001,50,-0.3*w,0.3*w);
+ }
+ }
+
+ // output directory
+ m_DATAmatrixFile->cd();
+ // copy tree "slice"
+ TTree *slice_tree;
+ m_MCmatrixFile->GetObject("slice",slice_tree);
+ if(slice_tree) {
+ double tmp_cmax;
+ double tmp_cmin;
+ int tmp_cslices;
+
+ double tmp_phimax;
+ double tmp_phimin;
+ int tmp_phislices;
+
+ double tmp_d0max;
+ double tmp_d0min;
+ int tmp_d0slices;
+
+ double tmp_z0max;
+ double tmp_z0min;
+ int tmp_z0slices;
+
+ double tmp_etamax;
+ double tmp_etamin;
+ int tmp_etaslices;
+
+ slice_tree->SetBranchAddress("c_max",&tmp_cmax);
+ slice_tree->SetBranchAddress("c_min", &tmp_cmin);
+ slice_tree->SetBranchAddress("c_slices", &tmp_cslices);
+ slice_tree->SetBranchAddress("phi_max", &tmp_phimax);
+ slice_tree->SetBranchAddress("phi_min", &tmp_phimin);
+ slice_tree->SetBranchAddress("phi_slices", &tmp_phislices);
+ slice_tree->SetBranchAddress("d0_max", &tmp_d0max);
+ slice_tree->SetBranchAddress("d0_min", &tmp_d0min);
+ slice_tree->SetBranchAddress("d0_slices", &tmp_d0slices);
+ slice_tree->SetBranchAddress("z0_max", &tmp_z0max);
+ slice_tree->SetBranchAddress("z0_min", &tmp_z0min);
+ slice_tree->SetBranchAddress("z0_slices", &tmp_z0slices);
+ slice_tree->SetBranchAddress("eta_max", &tmp_etamax);
+ slice_tree->SetBranchAddress("eta_min", &tmp_etamin);
+ slice_tree->SetBranchAddress("eta_slices", &tmp_etaslices);
+ slice_tree->GetEntry(0);
+ delete slice_tree;
+ slice_tree=0;
+
+ m_DATAmatrixFile->cd();
+ slice_tree=new TTree("slice","slice");
+ slice_tree->Branch("c_max",&tmp_cmax);
+ slice_tree->Branch("c_min", &tmp_cmin);
+ slice_tree->Branch("c_slices", &tmp_cslices);
+ slice_tree->Branch("phi_max", &tmp_phimax);
+ slice_tree->Branch("phi_min", &tmp_phimin);
+ slice_tree->Branch("phi_slices", &tmp_phislices);
+ slice_tree->Branch("d0_max", &tmp_d0max);
+ slice_tree->Branch("d0_min", &tmp_d0min);
+ slice_tree->Branch("d0_slices", &tmp_d0slices);
+ slice_tree->Branch("z0_max", &tmp_z0max);
+ slice_tree->Branch("z0_min", &tmp_z0min);
+ slice_tree->Branch("z0_slices", &tmp_z0slices);
+ slice_tree->Branch("eta_max", &tmp_etamax);
+ slice_tree->Branch("eta_min", &tmp_etamin);
+ slice_tree->Branch("eta_slices", &tmp_etaslices);
+ slice_tree->Fill();
+ slice_tree->Write(0,TObject::kOverwrite);
+
+ } else {
+ ATH_MSG_ERROR("no slice tree found in \""+m_MCmatrixFileName+"\"");
+ }
+ // copy tree "montree"
+ /*TTree *montree;
+ m_MCmatrixFile->GetObject("montree",montree);
+ if(montree) {
+ montree->Write();
+ } else {
+ ATH_MSG_WARNING("no montree tree found in \""+m_MCmatrixFileName+"\"");
+ } */
+ gErrorAbortLevel=kError;
+ // loop over all am trees
+ for(int tower=0;tower<64;tower++) {
+ if((m_tower>=0)&&(tower!=m_tower)) continue;
+ TTree *mcTree;
+ TString amName=TString::Format("am%d",tower);
+ m_MCmatrixFile->GetObject(amName,mcTree);
+#ifdef WRITE_DEBUG_TREE
+ debug_tower=tower;
+#endif
+ if(!mcTree) {
+ ATH_MSG_ERROR("no tree \""+amName+"\" found in \""+
+ m_MCmatrixFileName+"\"");
+ } else {
+ ATH_MSG_INFO(TString::Format("Processing tower %d",tower) );
+ vector<int> errorCountPlane(m_rmap->getPlaneMap()->getNPlanes());
+ map<int,int> badModulePix,badModuleSCT;
+ int errorCountE=0;
+ int errorCountH=0;
+ // rename input tree
+ mcTree->SetName("tmptree");
+ int ndim,ndim2,nplanes;
+ mcTree->SetBranchAddress("ndim",&ndim);
+ mcTree->SetBranchAddress("ndim2",&ndim2);
+ mcTree->SetBranchAddress("nplanes",&nplanes);
+ mcTree->GetEntry(0);
+ ATH_MSG_INFO(TString::Format("ndim=%d ndim2=%d nplanes=%d",
+ ndim,ndim2,nplanes) );
+ double *Vec=new double[ndim];
+ mcTree->SetBranchAddress("Vec",Vec);
+ int *sectorID=new int[nplanes];
+ mcTree->SetBranchAddress("sectorID",sectorID);
+ int *hashID=new int[nplanes];
+ mcTree->SetBranchAddress("hashID",hashID);
+ float nhit;
+ mcTree->SetBranchAddress("nhit",&nhit);
+ double tmp_C_D_Phi_Coto_Z[5];
+ mcTree->SetBranchAddress("tmpC",tmp_C_D_Phi_Coto_Z+0);
+ mcTree->SetBranchAddress("tmpD",tmp_C_D_Phi_Coto_Z+1);
+ mcTree->SetBranchAddress("tmpPhi",tmp_C_D_Phi_Coto_Z+2);
+ mcTree->SetBranchAddress("tmpCoto",tmp_C_D_Phi_Coto_Z+3);
+ mcTree->SetBranchAddress("tmpZ",tmp_C_D_Phi_Coto_Z+4);
+ double *tmpx_C_D_Phi_Coto_Z=new double[5*ndim];
+ mcTree->SetBranchAddress("tmpxC",tmpx_C_D_Phi_Coto_Z+0*ndim);
+ mcTree->SetBranchAddress("tmpxD",tmpx_C_D_Phi_Coto_Z+1*ndim);
+ mcTree->SetBranchAddress("tmpxPhi",tmpx_C_D_Phi_Coto_Z+2*ndim);
+ mcTree->SetBranchAddress("tmpxCoto",tmpx_C_D_Phi_Coto_Z+3*ndim);
+ mcTree->SetBranchAddress("tmpxZ",tmpx_C_D_Phi_Coto_Z+4*ndim);
+ double *tmpcovx=new double[ndim2];
+ mcTree->SetBranchAddress("tmpcovx",tmpcovx);
+#ifdef INCLUDE_SHORT_VECTOR
+ std::vector<short>
+ *tmpintc=new std::vector<short>,
+ *tmpintphi=new std::vector<short>,
+ *tmpintd0=new std::vector<short>,
+ *tmpintz0=new std::vector<short>,
+ *tmpinteta=new std::vector<short>;
+ mcTree->SetBranchAddress("tmpintc",&tmpintc);
+ mcTree->SetBranchAddress("tmpintphi",&tmpintphi);
+ mcTree->SetBranchAddress("tmpintd0",&tmpintd0);
+ mcTree->SetBranchAddress("tmpintz0",&tmpintz0);
+ mcTree->SetBranchAddress("tmpinteta",&tmpinteta);
+#endif
+ // create output tree
+ m_DATAmatrixFile->cd();
+ TTree *dataTree=new TTree
+ (amName,TString::Format("Ambank %d para",tower));
+
+ dataTree->Branch("ndim",&ndim,"ndim/I");
+ dataTree->Branch("ndim2",&ndim2,"ndim2/I");
+ dataTree->Branch("nplanes",&nplanes,"nplanes/I");
+ dataTree->Branch("Vec", Vec,"Vec[ndim]/D");
+ dataTree->Branch("sectorID",sectorID,"sectorID[nplanes]/I");
+ dataTree->Branch("hashID",hashID,"hashID[nplanes]/I");
+ dataTree->Branch("nhit", &nhit,"nhit/F");\
+ dataTree->Branch("tmpC",tmp_C_D_Phi_Coto_Z+0,"tmpC/D");
+ dataTree->Branch("tmpD",tmp_C_D_Phi_Coto_Z+1,"tmpD/D");
+ dataTree->Branch("tmpPhi",tmp_C_D_Phi_Coto_Z+2,"tmpPhi/D");
+ dataTree->Branch("tmpCoto",tmp_C_D_Phi_Coto_Z+3,"tmpCoto/D");
+ dataTree->Branch("tmpZ",tmp_C_D_Phi_Coto_Z+4,"tmpZ/D");
+ dataTree->Branch("tmpxC",tmpx_C_D_Phi_Coto_Z+0*ndim,"tmpxC[ndim]/D");
+ dataTree->Branch("tmpxD",tmpx_C_D_Phi_Coto_Z+1*ndim,"tmpxD[ndim]/D");
+ dataTree->Branch("tmpxPhi",tmpx_C_D_Phi_Coto_Z+2*ndim,"tmpxPhi[ndim]/D");
+ dataTree->Branch("tmpxCoto",tmpx_C_D_Phi_Coto_Z+3*ndim,"tmpxCoto[ndim]/D");
+ dataTree->Branch("tmpxZ",tmpx_C_D_Phi_Coto_Z+4*ndim,"tmpxZ[ndim]/D");
+ dataTree->Branch("tmpcovx",tmpcovx,"tmpcovx[ndim2]/D");
+ // drop these for tests
+#ifdef INCLUDE_SHORT_VECTOR
+ dataTree->Branch("tmpintc", &tmpintc);
+ dataTree->Branch("tmpintphi", &tmpintphi);
+ dataTree->Branch("tmpintd0", &tmpintd0);
+ dataTree->Branch("tmpintz0", &tmpintz0);
+ dataTree->Branch("tmpinteta", &tmpinteta);
+#endif
+ // loop over all entries (c.f. sectors)
+ for(int iEnt=0;iEnt<mcTree->GetEntries(); iEnt++) {
+ mcTree->GetEntry(iEnt);
+ if(nhit<ndim) {
+ errorCountH++;
+ continue;
+ }
+ // transform here the coordinates to the new module geometry
+ TVectorD avgT(5);
+ TVectorD avgX(ndim);
+ TMatrixD avgTX(5,ndim);
+ TMatrixDSym avgXX(ndim);
+ // <t>
+ for(int i=0;i<5;i++) {
+ avgT(i)=tmp_C_D_Phi_Coto_Z[i]/nhit;
+ }
+ // <x>
+ for(int i=0;i<ndim;i++) {
+ avgX(i)=Vec[i]/nhit;
+ }
+ // <tx>
+ for(int i=0;i<5;i++) {
+ for(int j=0;j<ndim;j++) {
+ avgTX(i,j)=tmpx_C_D_Phi_Coto_Z[j+ndim*i]/nhit-avgX(j)*avgT(i);
+ }
+ }
+ // <xx>-<x>*<x>
+ for(int j=0;j<ndim;j++) {
+ for(int i=0;i<=j;i++) {
+ avgXX(i,j)=tmpcovx[i*ndim+j]/nhit-avgX(i)*avgX(j);
+ avgXX(j,i)=avgXX(i,j);
+ }
+ }
+ TMatrixDSymEigen EXX(avgXX);
+ TVectorD eval=EXX.GetEigenValues();
+ TMatrixD evec=EXX.GetEigenVectors();
+ TMatrixDSym H(ndim);
+ int error=0;
+ for(int i=0;i<ndim;i++) {
+ for(int j=0;j<ndim;j++) {
+ for(int k=0;k<ndim;k++) {
+ if(eval(k)<=0.0) {
+ error++;
+ break;
+ }
+ H(i,j) += evec(i,k)*evec(j,k)/eval(k);
+ }
+ if(error) break;
+ }
+ if(error) break;
+ }
+ // if there is an error, skip further corrections
+ if(error) {
+ errorCountH++;
+ /* cout<<"nhit="<<nhit<<"\n";
+ cout<<"avgT\n";
+ avgT.Print();
+ cout<<"avgX\n";
+ avgX.Print();
+ cout<<"avgTX\n";
+ avgTX.Print();
+ cout<<"avgXX\n";
+ avgXX.Print();
+ exit(0); */
+ } else {
+ // matrix A: t = A*(x-t0) + t0
+ // TMatrixD A(avgTX,TMatrixD::kMult,H);
+ //
+ // <tt>-<t>*<t>
+ TMatrixDSym avgTT(H.Similarity(avgTX));
+ // eigenvalue decomposition
+ TMatrixDSymEigen ETT(avgTT);
+ TMatrixD T=ETT.GetEigenVectors();
+ TMatrixD Tinv(T);
+ TVectorD evT=ETT.GetEigenValues();
+ for(int i=0;i<5;i++) {
+ for(int j=0;j<5;j++) {
+ double sqrtEV=sqrt(evT(j));
+ T(i,j) *= sqrtEV;
+ Tinv(i,j) /= sqrtEV;
+ }
+ }
+ TMatrixD BTinv(avgTX,TMatrixD::kTransposeMult,Tinv);
+ //
+ // calculate alignment differences for DATA-MC
+ // for central track positions and shifted
+ // by +/- eigenvectors
+ TVectorD d_mcX(ndim);
+ TMatrixD d_mcXT(ndim,5);
+ TMatrixD d_mcXX(ndim,ndim);
+ int planeError=0;
+ for(int n=0;n<11;n++) {
+ // n=0: <t>
+ // n=1,2 <t> +/- <T0>
+ // n=3,4 <t> +/- <T1>
+ // ...
+ // n=9,10 <t> +/- <T4>
+ TVectorD track(avgT);
+ double sign= (n&1) ? 1. : -1.;
+ int iEV=(n-1)/2;
+ if(n) {
+ for(int j=0;j<5;j++) {
+ track(j) +=sign*T(j,iEV);
+ }
+ }
+ //cout<<"track\n";
+ //track.Print();
+ static double const B_FACTOR=-0.3*2.083;
+ double Q_over_pt=2.0*track(0);
+ double rho=B_FACTOR* Q_over_pt;
+ double d0=track(1);
+ double phi=track(2);
+ double cotTh=track(3);
+ double z0=track(4);
+ double sinPhi0=sin(phi);
+ double cosPhi0=cos(phi);
+ int icoord=0;
+ TVectorD d_mc(ndim);
+ for(int plane=0;plane<nplanes;plane++) {
+ int hash=hashID[plane];
+ double xData[2],xMC[2];
+ int errorMC=m_MCmodulePositionCalculator->
+ extrapolateTo(plane,hash,rho,cotTh,z0,d0,
+ sinPhi0,cosPhi0,&matrixVtx,xMC);
+ int errorData=m_DATAmodulePositionCalculator->
+ extrapolateTo(plane,hash,rho,cotTh,z0,d0,
+ sinPhi0,cosPhi0,&dataVtx,xData);
+ int ncoord=m_rmap->getPlaneMap()->isSCT(plane) ? 1 : 2;
+ if(errorMC||errorData) {
+ if(ncoord==2) {
+ badModulePix[hash]|=(1<<(errorMC+8))|(1<<errorData);
+ } else {
+ badModuleSCT[hash]|=(1<<(errorMC+8))|(1<<errorData);
+ }
+ planeError |= (1<<plane);
+ error++;
+ } else {
+ for(int j=0;j<ncoord;j++) {
+ d_mc(icoord+j) = xData[j]-xMC[j];
+ }
+#ifdef WRITE_DEBUG_TREE
+ if((n==0) && (m_debugFile)) {
+ debug_hash=hash;
+ debug_plane=plane;
+ for(int j=0;j<ncoord;j++) {
+ debug_hit[j]=avgX(icoord+j);
+ debug_ghost[j]=xMC[j];
+ debug_data[j]=xData[j];
+ }
+ for(int i=0;i<5;i++) {
+ debug_track[i]=track(i);
+ }
+ /*FTK_ModuleGeometry const *module=
+ m_MCmodulePositionCalculator->
+ findModule(plane,hash); */
+ if(ncoord==2) {
+ debugPix->Fill();
+ } else {
+ debugSCT->Fill();
+ }
+ }
+#endif
+ if(m_debugFile) {
+ // fill monitoring histograms
+ for(int j=0;j<ncoord;j++) {
+ if(debug_coordHM[icoord+j])
+ debug_coordHM[icoord+j]->Fill
+ (avgX(icoord+j),xMC[j]);
+ if(debug_coordMD[icoord+j])
+ debug_coordMD[icoord+j]->Fill
+ (xMC[j],xData[j]);
+ if(debug_resCurv[icoord+j])
+ debug_resCurv[icoord+j]->Fill
+ (xMC[j]-avgX(icoord+j),track(0));
+ }
+ }
+ }
+ icoord += ncoord;
+ }
+ if(n==0) {
+ // average difference data-mc
+ d_mcX=d_mc;
+ } else {
+ TVectorD deltaX(d_mc-d_mcX);
+ for(int i=0;i<ndim;i++) {
+ for(int k=0;k<5;k++) {
+ // sum of deltaX*deltaT
+ d_mcXT(i,k) += 0.5*deltaX(i)*sign*T(k,iEV);
+ }
+ for(int j=0;j<ndim;j++) {
+ // sums of deltaX*deltaX and correlations
+ // involving B=avgTX
+ d_mcXX(i,j) += 0.5*
+ (deltaX(i)*sign*BTinv(j,iEV)+
+ deltaX(j)*sign*BTinv(i,iEV)+
+ deltaX(i)*deltaX(j));
+ }
+ }
+ }
+ }
+ if(!error) {
+ // apply alignment changes to the coordinate sums
+
+ // coordinate shifts (vector d)
+ float nhit0=nhit;
+ //nhit=0;
+ for(int i=0;i<ndim;i++) {
+ Vec[i]+= d_mcX(i)*nhit;
+ }
+ // rotations (matrices M and E)
+ for(int i=0;i<5;i++) {
+ for(int j=0;j<ndim;j++) {
+ tmpx_C_D_Phi_Coto_Z[j+ndim*i]+=nhit*
+ (d_mcX(j)*avgT(i)+d_mcXT(j,i));
+ }
+ }
+ for(int j=0;j<ndim;j++) {
+ for(int i=0;i<=j;i++) {
+ tmpcovx[i*ndim+j]+= nhit*
+ (avgX(i)*d_mcX(j)+avgX(j)*d_mcX(i)+
+ d_mcX(i)*d_mcX(j)+d_mcXX(i,j));
+ tmpcovx[j*ndim+i]=tmpcovx[i*ndim+j];
+ }
+ }
+ nhit=nhit0;
+ } else {
+ for(size_t k=0;k<errorCountPlane.size();k++) {
+ if(planeError & (1<<k)) errorCountPlane[k]++;
+ //cout<<"plane "<<k<<" error\n";
+ }
+ errorCountE++;
+ //exit(0);
+ }
+ }
+
+ if((errorCountH |errorCountE)==0) {
+ dataTree->Fill();
+ }
+
+ /*if(//(iEnt<10)||
+ //((iEnt<100)&&(iEnt%10 ==0))||
+ //((iEnt<1000)&&(iEnt%100 ==0))||
+ ((iEnt<10000)&&(iEnt%1000 ==0))||
+ ((iEnt%10000 ==0))) {
+ ATH_MSG_INFO(TString::Format
+ ("sector %d nErrorE=%d nErrorH=%d",
+ iEnt,errorCountE,errorCountH));
+ } */
+ }
+
+ // save to output file
+ m_DATAmatrixFile->cd();
+ dataTree->Write(0,TObject::kOverwrite);
+
+ if(errorCountE|errorCountH) {
+ TString errorMsg=TString::Format
+ ("tower=%d errors (%d+%d)/%lld sectors with error",
+ tower,errorCountE,errorCountH,
+ mcTree->GetEntries());
+ if(errorCountE) {
+ errorMsg+=" by plane:";
+ for(size_t k=0;k<errorCountPlane.size();k++) {
+ errorMsg += TString::Format(" %d",errorCountPlane[k]);
+ }
+ }
+ ATH_MSG_WARNING(errorMsg);
+ if(badModulePix.size()) {
+ TString pixErrors("pixel module errors");
+ for(auto m=badModulePix.begin();m!=badModulePix.end();m++) {
+ pixErrors+=TString::Format(" %d:%x",(*m).first,(*m).second);
+ }
+ ATH_MSG_WARNING(pixErrors);
+ }
+ if(badModuleSCT.size()) {
+ TString sctErrors("SCT module errors");
+ for(auto m=badModuleSCT.begin();m!=badModuleSCT.end();m++) {
+ sctErrors+=TString::Format(" %d:%x",(*m).first,(*m).second);
+ }
+ ATH_MSG_WARNING(sctErrors);
+ }
+ }
+
+ delete [] tmpcovx;
+ delete [] tmpx_C_D_Phi_Coto_Z;
+ delete [] hashID;
+ delete [] sectorID;
+ delete [] Vec;
+ }
+ }
+ if(m_DATAmatrixFile) delete m_DATAmatrixFile;
+ m_DATAmatrixFile=0;
+#ifdef WRITE_DEBUG_TREE
+ if(debugPix) {
+ m_debugFile->cd();
+ debugPix->Write(0,TObject::kOverwrite);
+ }
+ if(debugSCT) {
+ m_debugFile->cd();
+ debugSCT->Write(0,TObject::kOverwrite);
+ }
+#endif
+ for(size_t i=0;i<debug_coordHM.size();i++) {
+ if(m_debugFile) m_debugFile->cd();
+ if(debug_coordHM[i]) debug_coordHM[i]->Write(0,TObject::kOverwrite);
+ if(debug_coordMD[i]) debug_coordMD[i]->Write(0,TObject::kOverwrite);
+ if(debug_resCurv[i]) debug_resCurv[i]->Write(0,TObject::kOverwrite);
+ }
+ if(m_debugFile) delete m_debugFile;
+ m_debugFile=0;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FTKDataMatrixFromMCAlgo::finalize() {
+ return StatusCode::SUCCESS;
+}
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKGhostHitCalculator.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKGhostHitCalculator.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c4039c692efd708c16b0d25824c64b63c788804b
--- /dev/null
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKGhostHitCalculator.cxx
@@ -0,0 +1,647 @@
+#include "TrigFTKBankGen/FTKGhostHitCalculator.h"
+#include "TrigFTKSim/FTKRootFile.h"
+#include "TrigFTKSim/FTKRegionMap.h"
+#include "TrigFTKSim/ftk_dcap.h"
+#include "TrigFTKSim/FTKTruthTrack.h"
+#include "TrigFTKSim/FTKSetup.h"
+#include <TTree.h>
+
+
+bool FTK_ModuleGeometry::extrapolate
+(double rho,double cotTh,double z0,double d0,double sinPhi0,double cosPhi0,
+ TVector3 const *trackVertex,double *xLocPtr,TVector2 *xPtr,double *sPtr)
+ const {
+ // rough track extrapolation in phi
+ TVector3 d(-d0*sinPhi0,d0*cosPhi0,z0);
+ if(trackVertex) d += *trackVertex;
+ TVector3 cd(m_center-d);
+ double r=TMath::Hypot(cd.X(),cd.Y());
+ double product=(cd.X()*cosPhi0+cd.Y()*sinPhi0)/r;
+ // skip modules if track is not within +/-25 degree in phi
+ if(product<0.5) {
+ /*std::cout<<"id="<<m_id<<"/"<<m_plane<<" rho="<<rho<<" cotTh="<<cotTh<<" z0="<<z0<<" d0="<<d0
+ <<" sinPhi="<<sinPhi0<<" cosPhi="<<cosPhi0
+ <<" center="<<m_center.X()<<","<<m_center.Y()<<","<<m_center.Z()
+ <<"\n";
+ //exit(0); */
+ return false;
+ }
+ TVector2 x;
+ TVector3 const &dir0(m_phiAxis);
+ TVector3 const &dir1(m_etaAxis);
+ // iterate to get accurate extrapolation to module plane
+ double s;
+ bool haveS=false;
+ for(int iter=0;iter<3;iter++) {
+ // 1/rho may be large (or infinite)
+ // but sin(s*rho)/rho and (1-cos(s*rho))/rho are finite
+ // further complication:
+ // 1-cos(s*rho) may be small -> calculation using cos() is inaccurate
+ // -> calculate using sin(s*rho/2)
+ double sRho;
+ if(haveS) {
+ sRho=s*rho;
+ } else {
+ double rRho=r*rho;
+ if(fabs(rRho)>1.E-10) {
+ sRho=2.*asin(0.5*rRho);
+ s=sRho/rho;
+ } else {
+ // 2*arcsin(x/2) ~= x*(1+x^2/24)
+ double corr=(1.+rRho*rRho/24.);
+ s=r*corr;
+ sRho=rRho*corr;
+ }
+ }
+ double sinSrho_byRho;
+ double one_minus_cosSrho_byRho;
+ if(fabs(sRho)>1.E-10) {
+ sinSrho_byRho=sin(sRho)/rho;
+ double h=sin(0.5*sRho);
+ // h*h = sin^2(0.5*sRho)
+ // = 1/2* (sin^2(0.5*sRho)+1-cos^2(0.5*sRho))
+ // = 1/2* (1 - cos(sRho))
+ one_minus_cosSrho_byRho = 2.*h*h/rho; //(1-cos(sRho))/rho;
+ } else {
+ // sin(x) ~= x*(1-x^2/6)
+ sinSrho_byRho=s-sRho*sRho*s/6.;
+ // 1-cos(x) ~= x^2*(1-x^2/24)
+ one_minus_cosSrho_byRho=s*sRho*(0.5-sRho*sRho/24.);
+ }
+ // h is the difference vector
+ // of the extrapolated helix and the module centre
+ TVector3 h(cosPhi0*sinSrho_byRho-sinPhi0*one_minus_cosSrho_byRho,
+ sinPhi0*sinSrho_byRho+cosPhi0*one_minus_cosSrho_byRho,
+ s*cotTh);
+ h-=cd;
+ // projection on strip phi axis
+ // projection on strip eta axis
+ x=TVector2(dir0.Dot(h),dir1.Dot(h));
+ // updated estimate of impact point on module plane minus
+ TVector3 impact=cd+dir0*x.X()+dir1*x.Y();
+ if(fabs(dir1.Z())>0.5) {
+ // barrel module
+ r=TMath::Hypot(impact.X(),impact.Y());
+ haveS=false;
+ } else {
+ // endcap module
+ s=impact.Z()/cotTh;
+ haveS=true;
+ }
+ }
+ if(xLocPtr) {
+ xLocPtr[0]=(x.X()+m_width.X()*0.5)/m_pitch.X();
+ xLocPtr[1]=(x.Y()+m_width.Y()*0.5)/m_pitch.Y();
+ }
+ if(xPtr) {
+ *xPtr = x;
+ }
+ if(sPtr) {
+ *sPtr =s;
+ }
+ return true;
+}
+
+FTKGhostHitCalculator::FTKGhostHitCalculator
+(const FTKRegionMap *rmap,std::string const &name)
+ : FTKLogging(name),m_rmap(rmap) {
+ // extract module to plane mapping
+ // m_idByPlane gives for each plane the set of modules
+ m_planeFromPIXid.clear();
+ m_plane_From_SCTid.clear();
+ for(int tower=0;tower<m_rmap->getNumRegions();tower++) {
+ for(int plane=0;plane<m_rmap->getNPlanes();plane++) {
+ std::map<int,int> &planeByID=m_rmap->getPlaneMap()->isPixel(plane) ?
+ m_planeFromPIXid : m_plane_From_SCTid;
+ std::map<unsigned int,unsigned int> tpMap=
+ m_rmap->getGlobalToLocalMapping(tower,plane);
+ for(std::map<unsigned int,unsigned int>::const_iterator
+ iID=tpMap.begin();iID!=tpMap.end();iID++) {
+ planeByID[(*iID).first]=plane;
+ }
+ }
+ }
+}
+
+int FTKGhostHitCalculator::loadModuleGeometry
+(std::string const &name,bool markAllBad,bool invertIBLphi) {
+ TDirectory *file=FTKRootFile::Instance()->OpenRootFileReadonly(name.c_str());
+ if(file) {
+ Info("LoadModueGeoemetry")
+ <<"reading \""<<name<<"\""<<" markAllBad="<<markAllBad
+ <<" invertIBLphi="<<invertIBLphi<<"\n";
+ int r=loadModuleGeometry(file,markAllBad,invertIBLphi);
+ delete file;
+ if(!r) {
+ Error("loadModuleGeometry")
+ <<"failed to load module geometry from "<<name<<"\n";
+ }
+ return r;
+ } else {
+ Error("loadModuleGeometry")<<"failed to open "<<name<<" as root file\n";
+ return 0;
+ }
+}
+
+int FTKGhostHitCalculator::addBadModules(std::string const &name) {
+ ftk_dcap::istream moduleList;
+ const bool ok_read = ftk_dcap::open_for_read(name,moduleList);
+ if (!moduleList || !ok_read) {
+ Error("addBadModules")<<"failed to open file "<<name<<" for reading\n";
+ return 0;
+ } else {
+ int r=addBadModules(moduleList);
+ if(!r) {
+ Error("addBadModules")<<"failed to read bad modules from "<<name<<"\n";
+ }
+ return r;
+ }
+}
+
+int FTKGhostHitCalculator::addBadModules(std::istream &in) {
+ int r=1;
+ std::set<int> SCTmodules,PIXmodules;
+ std::string line;
+ int nline=0;
+ while(getline(in,line)) {
+ nline++;
+ std::istringstream lineRead(line);
+ int tmpisPixel,tmpBEC,tmpSector,tmpPlane,tmpEtaModule,tmpPhiModule,
+ tmpSection,tmpidhash;
+ lineRead >>tmpisPixel >>tmpBEC >>tmpSector >>tmpPlane >>tmpEtaModule
+ >>tmpPhiModule >>tmpSection>>tmpidhash;
+ if(lineRead.fail()) {
+ Error("addBadModules")
+ <<"problem to read bad modules, line="<<nline
+ <<" : "<<line<<"\n";
+ r=0;
+ continue;
+ }
+ std::set<int> &moduleSet=(tmpisPixel ? PIXmodules : SCTmodules);
+ if(!moduleSet.insert(tmpidhash).second) {
+ Warning("addBadModules")
+ <<"skipping duplicate module isPixel="<<tmpisPixel<<" id="
+ <<tmpidhash<<"\n";
+ }
+ }
+ int nPixRead=PIXmodules.size();
+ int nPix=addBadPixelModules(PIXmodules);
+ int nSCTRead=SCTmodules.size();
+ int nSCT=addBadSCTModules(SCTmodules);
+ if(nPixRead==nPix) {
+ Info("addBadModules")<<"Using "<<nPix<<" dead pixel modules\n";
+ } else {
+ Info("addBadModules")
+ <<"Using "<<nPix<<" dead pixel modules ("<<nPixRead<<" requested)\n";
+ }
+ if(nSCTRead==nSCT) {
+ Info("addBadModules")<<"Using "<<nSCT<<" dead SCT modules\n";
+ } else {
+ Info("addBadModules")
+ <<"Using "<<nSCT<<" dead SCT modules ("<<nSCTRead<<" requested)\n";
+ }
+ return r;
+}
+
+FTK_ModuleGeometry const *FTKGhostHitCalculator::findModule(size_t plane,int moduleID) const {
+ FTK_ModuleGeometry const *r=0;
+ if(plane<m_geometryData.size()) {
+ auto modulePtr=m_geometryData[plane].find(moduleID);
+ if(modulePtr!=m_geometryData[plane].end()) {
+ r= &(*modulePtr).second;
+ }
+ }
+ return r;
+}
+
+int FTKGhostHitCalculator::extrapolateTo
+(size_t plane,int moduleID,double rho,double cotTh,double z0,double d0,
+ double sinPhi0,double cosPhi0,TVector3 const *trackVertex,double *x) const {
+ int r=0;
+ if(plane<m_geometryData.size()) {
+ auto modulePtr=m_geometryData[plane].find(moduleID);
+ if(modulePtr!=m_geometryData[plane].end()) {
+ if(!(*modulePtr).second.extrapolate
+ (rho,cotTh,z0,d0,sinPhi0,cosPhi0,trackVertex,x,0,0)) {
+ r=3;
+ }
+ } else {
+ r=2;
+ //Warning("extrapolateTo")
+ // <<"module "<<moduleID<<" not found in plane "<<plane<<"\n";
+ }
+ } else {
+ r=1;
+ Error("extrapolateTo")
+ <<"bad plane number "<<plane<<" max="<<m_geometryData.size()<<"\n";
+ }
+ return r;
+}
+
+int FTKGhostHitCalculator::loadModuleGeometry
+(TDirectory *source,bool markAllBad,bool invertIBLphi) {
+ // initialize empty list of modules
+ m_geometryData.resize(0);
+ m_geometryData.resize(m_rmap->getPlaneMap()->getNPlanes());
+
+ // read modules from TTree
+ TTree *tree;
+ source->GetObject("modulePositions",tree);
+ int r=1;
+
+ std::set<int> SCTmodules,PIXmodules;
+ if(tree) {
+ Int_t idhash;
+ Int_t isbad,isPixel,hitSector,isBLayer;
+ Int_t swapPhi,swapEta;
+ Float_t center[3],phiAxis[3],etaAxis[3],width,length,phiPitch,etaPitch;
+ Float_t sinTilt;
+ tree->SetBranchAddress("id",&idhash);
+ tree->SetBranchAddress("isPixel",&isPixel);
+ tree->SetBranchAddress("isBLayer",&isBLayer);
+ tree->SetBranchAddress("hitSector",&hitSector);
+ tree->SetBranchAddress("width",&width);
+ tree->SetBranchAddress("length",&length);
+ tree->SetBranchAddress("phiPitch",&phiPitch);
+ tree->SetBranchAddress("etaPitch",&etaPitch);
+ tree->SetBranchAddress("phiAxis",phiAxis);
+ tree->SetBranchAddress("etaAxis",etaAxis);
+ tree->SetBranchAddress("center",center);
+ tree->SetBranchAddress("isbad",&isbad);
+ tree->SetBranchAddress("swapPhi",&swapPhi);
+ tree->SetBranchAddress("swapEta",&swapEta);
+ tree->SetBranchAddress("sinTilt",&sinTilt);
+ m_geometryData.resize(m_rmap->getNPlanes());
+ int nLoad[2],nFTK[2],nDuplicate[2];
+ std::vector<int> rejected(4);
+ for(int i=0;i<2;i++) {
+ nLoad[i]=0;
+ nFTK[i]=0;
+ nDuplicate[i]=0;
+ }
+ for(int i=0;i<tree->GetEntries();i++) {
+ tree->GetEntry(i);
+ nLoad[isPixel ? 1 : 0] ++;
+ int plane=-1;
+ std::map<int,int> const &id2plane=
+ (isPixel ? m_planeFromPIXid : m_plane_From_SCTid);
+ std::map<int,int>::const_iterator iID=id2plane.find(idhash);
+ if(iID != id2plane.end()) {
+ plane=(*iID).second;
+ } else {
+ plane=isPixel ? -2 : -1;
+ }
+ if((plane>=0)&&(plane<(int)m_geometryData.size())) {
+ nFTK[isPixel ? 1 : 0] ++;
+ std::map<int,FTK_ModuleGeometry>::iterator iModule=
+ m_geometryData[plane].find(idhash);
+ if(iModule!=m_geometryData[plane].end()) {
+ nDuplicate[isPixel ? 1 : 0] ++;
+ Warning("loadModuleGeometry")
+ <<"skipping duplicate module isPixel="
+ <<isPixel<<" id="<<idhash<<"\n";
+ } else {
+ FTK_ModuleGeometry &m=m_geometryData[plane][idhash];
+ m.m_center.SetXYZ(center[0],center[1],center[2]);
+ m.m_phiAxis.SetXYZ(phiAxis[0],phiAxis[1],phiAxis[2]);
+ m.m_etaAxis.SetXYZ(etaAxis[0],etaAxis[1],etaAxis[2]);
+
+ // fix for IBL inverted phi axis
+ if(isBLayer && invertIBLphi) {
+ m.m_phiAxis *= -1.0;
+ }
+ //m.m_phiAxis.SetMag(1.0);
+ //m.m_etaAxis.SetMag(1.0);
+ double test=m.m_phiAxis.Dot(m.m_etaAxis);
+ if(TMath::Abs(test)>1.E-5) {
+ Error("loadModuleGeometry")
+ <<"phi axis ["<<m.m_phiAxis.X()<<","<<m.m_phiAxis.Y()
+ <<","<<m.m_phiAxis.Z()<<"] and eta axis ["
+ <<m.m_etaAxis.X()<<","<<m.m_etaAxis.Y()
+ <<","<<m.m_etaAxis.Z()<<"] not orthogonal "<<test<<"\n";
+ }
+ m.m_pitch.Set(phiPitch,etaPitch);
+ m.m_width.Set(width,length);
+ m.m_isPixel=isPixel;
+ m.m_id=idhash;
+ m.m_plane=plane;
+ m.m_hitSector=hitSector;
+ if(isbad || markAllBad) {
+ if(isPixel) {
+ PIXmodules.insert(idhash);
+ } else {
+ SCTmodules.insert(idhash);
+ }
+ }
+ }
+ } else {
+ int iReject=isPixel ? 0 : 1;
+ if(center[2]>0.) iReject+=2;
+ rejected[iReject]++;
+ }
+ }
+ delete tree;
+ Info("loadModuleGeometry")
+ <<"Modules found: SCT "<<nFTK[0]<<"/"<<nLoad[0]
+ <<" rejected: "<<rejected[0]<<"/"<<rejected[1]
+ <<"/"<<rejected[2]<<"/"<<rejected[3]
+ <<" Pixel "<<nFTK[1]<<"/"<<nLoad[1]
+ <<" Duplicates: SCT="<<nDuplicate[0]<<" Pixel="<<nDuplicate[1]<<"\n";
+ } else {
+ Error("loadModuleGeometry")<<"unable to load TTree \"modulePositions\"\n";
+ r=0;
+ }
+ addBadPixelModules(PIXmodules);
+ addBadSCTModules(SCTmodules);
+ /* for(size_t plane=0;plane<m_geometryData.size();plane++) {
+ std::cout<<"=========== modulelocator for plane="<<plane<<"\n";
+ FTK_ModuleLocator &moduleLocator=m_moduleLocatorByPlane[plane];
+ moduleLocator.print();
+ } */
+ return r;
+}
+
+void FTKGhostHitCalculator::clearBadModuleList(void) {
+ Info("clearBadModuleList")<<"reset bad modules, pixel "<<m_badPixel.size()
+ <<" SCT "<<m_badSCT.size()<<"\n";
+ m_badPixel.clear();
+ m_badSCT.clear();
+ updateCalculator();
+}
+
+void FTKGhostHitCalculator::updateCalculator(void) {
+ //
+ m_moduleLocatorByPlane.resize(0);
+ m_moduleLocatorByPlane.resize(m_geometryData.size());
+ //
+ for(size_t plane=0;plane<m_geometryData.size();plane++) {
+ std::map<int,FTK_ModuleGeometry> const &planeData(m_geometryData[plane]);
+ FTK_ModuleLocator &moduleLocator=m_moduleLocatorByPlane[plane];
+ for(std::map<int,FTK_ModuleGeometry>::const_iterator im=planeData.begin();
+ im!=planeData.end();im++) {
+ int id=(*im).first;
+ FTK_ModuleGeometry const *module=&((*im).second);
+ std::set<int> const &badList=module->m_isPixel ? m_badPixel : m_badSCT;
+ if(badList.find(id)!=badList.end()) {
+ // insert this module into the module locator
+ moduleLocator.insert(module);
+ }
+ }
+ }
+}
+
+FTKHit *FTKGhostHitCalculator::addGhostHit
+(FTKTruthTrack const &track,int plane, FTK_ModuleGeometry *module) const {
+ return m_moduleLocatorByPlane[plane].locate(track,module);
+}
+
+bool FTKGhostHitCalculator::isBad(int plane,int moduleID) const {
+ std::map<int,FTK_ModuleGeometry> const &planeData=m_geometryData[plane];
+ std::map<int,FTK_ModuleGeometry>::const_iterator iID=
+ planeData.find(moduleID);
+ if(iID!=planeData.end()) {
+ std::set<int> const &badSet=
+ (*iID).second.m_isPixel ? m_badPixel : m_badSCT;
+ return badSet.find(moduleID)!=badSet.end();
+ }
+ return false;
+}
+
+int FTKGhostHitCalculator::addBadPixelModules(std::set<int> const &bad) {
+ int n0=m_badPixel.size();
+ m_badPixel.insert(bad.begin(),bad.end());
+ Info("addBadPixelModules")
+ <<"inserting "<<bad.size()<<" bad modules, nBad increases from "
+ <<n0<<" to "<< m_badPixel.size()<<"\n";
+ updateCalculator();
+ return 0;
+}
+
+int FTKGhostHitCalculator::addBadSCTModules(std::set<int> const &bad) {
+ int n0=m_badSCT.size();
+ m_badSCT.insert(bad.begin(),bad.end());
+ Info("addBadSCTModules")
+ <<"inserting "<<bad.size()<<" bad modules, nBad increases from "
+ <<n0<<" to "<< m_badSCT.size()<<"\n";
+ updateCalculator();
+ return 0;
+}
+
+void FTK_ModuleLocator::insert(FTK_ModuleGeometry const *module) {
+ std::vector<double> ri(2),zi(2);
+ for(int ix=0;ix<2;ix++) {
+ for(int iy=0;iy<2;iy++) {
+ TVector3 pos=module->m_center
+ +module->m_width.X()*(ix-0.5)*module->m_phiAxis
+ +module->m_width.Y()*(iy-0.5)*module->m_etaAxis;
+ double r=pos.Perp();
+ double z=pos.Z();
+ if((ix==0)&&(iy==0)) {
+ ri[0]=r;
+ ri[1]=r;
+ zi[0]=z;
+ zi[1]=z;
+ } else {
+ ri[0]=TMath::Min(r,ri[0]);
+ ri[1]=TMath::Max(r,ri[1]);
+ zi[0]=TMath::Min(z,zi[0]);
+ zi[1]=TMath::Max(z,zi[1]);
+ }
+ }
+ }
+ int isBarrel;
+ double x[2],y[2],dxMax;
+ if(ri[1]-ri[0]<zi[1]-zi[0]) {
+ // barrel
+ isBarrel=1;
+ x[0]=ri[0];
+ x[1]=ri[1];
+ y[0]=zi[0];
+ y[1]=zi[1];
+ dxMax=20.;
+ } else {
+ // endcap
+ isBarrel=0;
+ x[0]=zi[0];
+ x[1]=zi[1];
+ y[0]=ri[0];
+ y[1]=ri[1];
+ dxMax=10.;
+ }
+ size_t iMS;
+ for(iMS=0;iMS<m_moduleSetVector.size();iMS++) {
+ FTK_ModuleSet &moduleSet(m_moduleSetVector[iMS]);
+ if(moduleSet.m_isBarrel != isBarrel) continue;
+ double xmin=TMath::Min(x[0],moduleSet.m_x[0]);
+ double xmax=TMath::Max(x[1],moduleSet.m_x[1]);
+ if(xmax-xmin<=dxMax) {
+ moduleSet.m_x[0]=xmin;
+ moduleSet.m_x[1]=xmax;
+ moduleSet.m_dy=
+ TMath::Max(y[1]-y[0],moduleSet.m_dy);
+ moduleSet.m_modules.insert
+ (std::make_pair(0.5*(y[1]+y[0]),module));
+ break;
+ }
+ }
+ if(iMS==m_moduleSetVector.size()) {
+ // start a new FTK_ModuleSet
+ m_moduleSetVector.resize(iMS+1);
+ FTK_ModuleSet &moduleSet(m_moduleSetVector[iMS]);
+ moduleSet.m_isBarrel=isBarrel;
+ moduleSet.m_x[0]=x[0];
+ moduleSet.m_x[1]=x[1];
+ moduleSet.m_dy=y[1]-y[0];
+ moduleSet.m_modules
+ .insert(std::make_pair(0.5*(y[1]+y[0]),module));
+ }
+}
+
+void FTK_ModuleLocator::print(void) const {
+ std::cout<<"FTK_ModuleLocator::print "<<m_moduleSetVector.size()
+ <<" module sets\n";
+ for(size_t i=0;i<m_moduleSetVector.size();i++) {
+ FTK_ModuleSet const &moduleSet=m_moduleSetVector[i];
+ std::cout<<" moduleSet["<<i<<"] ";
+ if(moduleSet.m_isBarrel) {
+ std::cout<<"barrel R=[";
+ } else {
+ std::cout<<"endcap Z=[";
+ }
+ std::cout<<moduleSet.m_x[0]<<","<<moduleSet.m_x[1]
+ <<"] n(modules)="<<moduleSet.m_modules.size()
+ <<"\n";
+ for(std::multimap<double,FTK_ModuleGeometry const *>
+ ::const_iterator ix=moduleSet.m_modules.begin();
+ ix!=moduleSet.m_modules.end();) {
+ double z=(*ix).first;
+ if(moduleSet.m_isBarrel) {
+ std::cout<<" Z=["<<z;
+ } else {
+ std::cout<<" R=["<<z;
+ }
+ std::pair<
+ std::multimap<double,FTK_ModuleGeometry const *>::const_iterator,
+ std::multimap<double,FTK_ModuleGeometry const *>::const_iterator>
+ er=moduleSet.m_modules.equal_range(z);
+ for(ix=er.first;ix!=er.second;ix++) {
+ std::cout<<" "<<(*ix).second->m_id;
+ }
+ std::cout<<" ]\n";
+ }
+ }
+}
+
+FTKHit *FTK_ModuleLocator::locate(FTKTruthTrack const &track,
+ FTK_ModuleGeometry *ret_module) const {
+ static double const B_FACTOR=-0.3*2.083;
+ double rho=B_FACTOR/track.getPT()*track.getQ();
+ double cotTh=track.getPZ()/track.getPT();
+ double z0=track.getZ();
+ double d0=track.getD0();
+ double sinPhi=track.getPY()/track.getPT();
+ double cosPhi=track.getPX()/track.getPT();
+
+ std::multimap<double,std::pair<TVector2,FTK_ModuleGeometry const *> >
+ candidates;
+
+ for(size_t im=0;im<m_moduleSetVector.size();im++) {
+ /* std::cout<<" getCandidate moduleset["<<im<<"] rho="<<rho
+ <<" cotth="<<cotTh<<" z0="<<z0<<" d0="<<d0
+ <<" sinPhi="<<sinPhi<<" cosPhi="<<cosPhi<<"\n"; */
+ FTK_ModuleSet const &moduleSet=m_moduleSetVector[im];
+ moduleSet.getCandidates(rho,cotTh,z0,d0,sinPhi,cosPhi,
+ candidates);
+ }
+ // select candidate with smallest s
+ FTKHit *r=0;
+ for(std::multimap<double,std::pair<TVector2,FTK_ModuleGeometry const *> >
+ ::const_iterator ii=candidates.begin();ii!=candidates.end();ii++) {
+ /*FTK_ModuleGeometry const *module=(*ii).second.second;
+ double s=module->m_center.Pt();
+ std::cout<<" CAND hash="<<module->m_id
+ <<" cosphi(cent)="<<module->m_center[0]/s
+ <<" sinphi(cent)="<<module->m_center[1]/s
+ <<" track: "<<cosPhi<<" "<<sinPhi<<"\n"; */
+ }
+ if(candidates.begin()!=candidates.end()) {
+ TVector2 const &x=(*candidates.begin()).second.first;
+ FTK_ModuleGeometry const *module=(*candidates.begin()).second.second;
+ r=new FTKHit(module->m_isPixel ? 2 : 1);
+ r->setSector(module->m_hitSector);
+ r->setITkMode(FTKSetup::getFTKSetup().getITkMode());
+ r->setIdentifierHash(module->m_id);
+ r->setPhiWidth(1);
+ r->setEtaWidth(1);
+ r->setPlane(module->m_plane);
+ r->setCoord(0,(x.X()+module->m_width.X()*0.5)/module->m_pitch.X());
+ if(module->m_isPixel) {
+ r->setCoord(1,(x.Y()+module->m_width.Y()*0.5)/module->m_pitch.Y());
+ }
+ if(ret_module) *ret_module=*module;
+ }
+ return r;
+}
+
+void FTK_ModuleLocator::FTK_ModuleSet::getCandidates
+(double rho,double cotTh,double z0,double d0,double sinPhi0,double cosPhi0,
+ std::multimap <double,std::pair<TVector2,FTK_ModuleGeometry const *> >
+ &candidate) const {
+ // calculate coordinates along module extent
+ // y = Z in barrel
+ // y = R im endcap
+ double y[2];
+ for(int k=0;k<2;k++) {
+ if(m_isBarrel) {
+ // barrel m_x are Rmin and Rmax
+ double r=m_x[k];
+ double rrho=r*rho;
+ double s=r*(1.+rrho*rrho/24.);
+ y[k]=s*cotTh+z0;
+ } else {
+ double z=m_x[k];
+ double s=(z-z0)/cotTh;
+ double srho=s*rho;
+ y[k]=s*(1.-srho*srho/24.);
+ }
+ }
+ if(y[0]>y[1]) {
+ double yy=y[0];
+ y[0]=y[1];
+ y[1]=yy;
+ }
+ /*std::cout<<" isBarrel="<<m_isBarrel<<" x[]="<<m_x[0]<<" "<<m_x[1]<<" y[]="<<y[0]<<" "<<y[1]<<"\n"; */
+
+ // loop over a range of modules
+ // start at module with smallest y where y>= y[0]-0.5*dy
+ for(std::multimap<double,FTK_ModuleGeometry const *>
+ ::const_iterator iModule=m_modules.lower_bound(y[0]-m_dy*0.5);
+ iModule!=m_modules.end();iModule++) {
+ // leave loop if
+ // y> y[1]+0.5*dy
+ // because all following modules have even larger y
+ /*std::cout<<" check: "<<(*iModule).first<<" >= "<<(y[0]-m_dy*0.5)
+ <<" is: "<<(*iModule).first<<" ?<= "<<(y[1]+m_dy*0.5)<<"?"; */
+ if((y[1]+m_dy*0.5)<(*iModule).first) {
+ //std::cout<<" no\n";
+ break;
+ }
+ //std::cout<<" yes\n";
+ // here we have a module with
+ // y[0]-d/2<=y and y<=y[1]+d/2
+ // -> y[0]<= y+d/2 and y[1]>=y-d/2
+ // the estimated extrapolated track yt, which is between y[0] and y[1]
+ // possibly is on the module [distance d/2 around the centre y].
+ FTK_ModuleGeometry const *module=(*iModule).second;
+
+ TVector2 x;
+ double s;
+ module->extrapolate(rho,cotTh,z0,d0,sinPhi0,cosPhi0,0,0,&x,&s);
+
+ // check whether coordinates are on the module
+ if(TMath::Abs(x.X())>0.5*module->m_width.X()) continue;
+ if(TMath::Abs(x.Y())>0.5*module->m_width.Y()) continue;
+ // add to the list of results
+ candidate.insert(std::make_pair(s,std::make_pair(x,module)));
+ }
+}
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRoot.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRoot.cxx
index 0941b695cf6e32014bb67e9e688cda7cd488095b..6cab5f923a5236e5d03c462934142d38a6fc2b59 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRoot.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRoot.cxx
@@ -1069,9 +1069,10 @@ void FTKPattGenRoot::SetModuleGeometryCheck(const std::string &fileName,
<<"select="<<" (undefined selection algorithm)\n";
}
} else {
+ Fatal("SetModuleGeometryCheck")
+ <<"module geometry "<<fileName<<" not found but select="
+ <<m_select<<"\n";
m_select= RndSector;
- Error("SetModuleGeometryCheck")
- <<"module geometry "<<fileName<<" not found\n";
Warning("SetModuleGeometryCheck")
<<"revert to select="<<m_select
<<": no boudary check, random sector selection\n";
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRootAlgo.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRootAlgo.cxx
index 5deb6d0a1485712eacab3d09fa576c5d6341c932..409544c0e3eed4bb90882f8c708c03dcf74ff365 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRootAlgo.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKPattGenRootAlgo.cxx
@@ -209,9 +209,13 @@ StatusCode FTKPattGenRootAlgo::initialize() {
ftkset.setIBLMode(m_IBLMode);
ftkset.setfixEndcapL0(m_fixEndcapL0);
ftkset.setITkMode(m_ITkMode);
+ log << MSG::INFO <<"HWMODEID="<<m_HWMODEID
+ <<" IBLMode="<<m_IBLMode<<" fixEndcapL0="<<m_fixEndcapL0
+ <<" ITKMode="<<m_ITkMode<<"\n";
+
// --- Create the pmap file object
- log << MSG::INFO << "RunPattGen() Make FTKPlaneMap." << endmsg;
+ log << MSG::INFO << "make FTKPlaneMap "<<m_pmap_path << endmsg;
FTKPlaneMap* pmap = new FTKPlaneMap(m_pmap_path.c_str());
if (!(*pmap)) {
log << MSG::FATAL << "Error using plane map: " << m_pmap_path << endmsg;
@@ -219,34 +223,38 @@ StatusCode FTKPattGenRootAlgo::initialize() {
}
// --- Create region map object
- log << MSG::INFO << "RunPattGen() Make FTKRegionMap." << endmsg;
+ log << MSG::INFO << "Make FTKRegionMap "<<m_rmap_path << endmsg;
FTKRegionMap* rmap = new FTKRegionMap(pmap, m_rmap_path.c_str());
if(m_HWMODEID==2) {
if(m_modulelut_path.empty()) {
- log << MSG::FATAL <<"RunPattGen() A module LUT is required when HW SS calculation is required"<<endmsg;
+ log << MSG::FATAL <<"A module LUT is required when HW SS calculation is required"<<endmsg;
return StatusCode::FAILURE;
} else {
+ log << MSG::INFO << "Load LUT "<<m_modulelut_path<<endmsg;
rmap->loadModuleIDLUT(m_modulelut_path.c_str());
}
}
// --- Create ssmap
const bool force_am_hashmap = false;
- log << MSG::INFO << "RunPattGen() Make FTKSSMap." << endmsg;
+ log << MSG::INFO << "Make FTKSSMap "<<m_ssmap_path << endmsg;
FTKSSMap* ssmap = new FTKSSMap(rmap, m_ssmap_path.c_str(), force_am_hashmap);
// --- Create the slices file object
- log << MSG::INFO << "RunPattGen() Make FTKSectorSlice." << endmsg;
+ log << MSG::INFO << "Make FTKSectorSlice "
+ <<m_slices_path << endmsg;
FTKSectorSlice* sectorslice = new FTKSectorSlice();
sectorslice->loadSlices(m_slices_path);
// --- create FTKConstantBank
- log << MSG::INFO << "RunPattGen() Make FTKConstantBank." << endmsg;
+ log << MSG::INFO << "Make FTKConstantBank "
+ <<m_fitconstants_path << endmsg;
FTKConstantBank* constbank = new FTKConstantBank(pmap->getTotalDim(), m_fitconstants_path.c_str());
// --- create pattgen object
- log << MSG::INFO << "RunPattGen() Make FTKPattGenRoot."
+ log << MSG::INFO << "Make FTKPattGenRoot "
+ <<" keep7of8="<<m_keep7of8<<" tolerance7of8="<<m_tolerance7of8
<< endmsg;
m_pattgen=new FTKPattGenRoot(m_curreg,ssmap,sectorslice,constbank,m_keep7of8,
@@ -255,20 +263,33 @@ StatusCode FTKPattGenRootAlgo::initialize() {
,m_propagator
#endif
);
- log << MSG::INFO << "RunPattGen() beam spot at "
+ log << MSG::INFO << "beam spot at "
<<m_beamspotX<<" "<<m_beamspotY
<< endmsg;
m_pattgen->SetRandomNumberGenerator(m_rndmSvc->GetEngine(m_rndStreamName));
m_pattgen->setBeamspot(m_beamspotX,m_beamspotY);
+ log << MSG::INFO << "Read sector path "<<m_sectors_path<<endmsg;
m_pattgen->ReadSectorFile(m_sectors_path); // set sectors path
+ log<< MSG::INFO << "Slice parameters:"
+ <<" phi: "<<m_phi_min<<" "<<m_phi_max
+ <<" c: "<<m_c_min<<" "<<m_c_max
+ <<" d0: "<<m_d0_min<<" "<<m_d0_max
+ <<" z0: "<<m_z0_min<<" "<<m_z0_max
+ <<" eta: "<<m_eta_min<<" "<<m_eta_max
+ <<endmsg;
m_pattgen->SetSliceParameters(m_phi_min,m_phi_max, m_c_min, m_c_max,
m_d0_min, m_d0_max, m_z0_min, m_z0_max, m_eta_min, m_eta_max);
+ log<< MSG::INFO << "D0 exponent: "<<m_d0_alpha<<endmsg;
m_pattgen->SetD0Exponent(m_d0_alpha);
+ log<< MSG::INFO <<"Overlap removal: "<<m_overlap<<endmsg;
m_pattgen->SetOverlapRemoval(m_overlap);
+ log<<MSG::INFO <<"Module boundary check: "<<m_sectorSelection
+ <<" module position file: \""<<m_ModuleGeometryFile<<"\""<<endmsg;
m_pattgen->SetModuleGeometryCheck
(m_ModuleGeometryFile,(FTKPattGenRoot::SectorSelection)m_sectorSelection);
// open output file
+ log<<MSG::INFO <<"Output file: "<<m_OutputFile<<endmsg;
m_pattgen->SetRootOutput(m_OutputFile);
return StatusCode::SUCCESS;
@@ -305,11 +326,12 @@ StatusCode FTKPattGenRootAlgo::finalize() {
void FTKPattGenRootAlgo::PostMessage(void) {
- if (FTKLogger::fType==0) ATH_MSG_FATAL(fBuffer->str());
- else if(FTKLogger::fType==1) ATH_MSG_ERROR(fBuffer->str());
- else if(FTKLogger::fType==2) ATH_MSG_WARNING(fBuffer->str());
- else if(FTKLogger::fType==3) ATH_MSG_INFO(fBuffer->str());
- else if(FTKLogger::fType==4) ATH_MSG_DEBUG(fBuffer->str());
+ int fType=getLoggerMsgType();
+ if (fType==0) ATH_MSG_FATAL(getLoggerMsg());
+ else if(fType==1) ATH_MSG_ERROR(getLoggerMsg());
+ else if(fType==2) ATH_MSG_WARNING(getLoggerMsg());
+ else if(fType==3) ATH_MSG_INFO(getLoggerMsg());
+ else if(fType==4) ATH_MSG_DEBUG(getLoggerMsg());
}
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKSectorSlice.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKSectorSlice.cxx
index c97883361efa3ef63749ec2d0f0d5706eddaa9bc..9b05204e2ae56953e922de5417940f1c741edb7b 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/FTKSectorSlice.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/FTKSectorSlice.cxx
@@ -237,25 +237,41 @@ vector<int>* FTKSectorSlice::searchSectors(FTKTrack& track) {
int error=0;
+ static int print=200;
if ((track.getPhi()<m_phi_min) || (track.getPhi()>m_phi_max)) {
- cout << "Error: wrong phi parameter " << track.getPhi() << '\n';
+ if(print) {
+ cout << "Error: wrong phi parameter " << track.getPhi() << '\n';
+ print--;
+ }
error++;
}
if ((2.*track.getHalfInvPt()<m_c_min) ||
(2.*track.getHalfInvPt()>m_c_max)) {
- cout << "Error: wrong c parameter " << 2.*track.getHalfInvPt() << '\n';
+ if(print) {
+ cout << "Error: wrong c parameter " << 2.*track.getHalfInvPt() << '\n';
+ print--;
+ }
error++;
}
if ((track.getIP()<m_d0_min) || (track.getIP()>m_d0_max)) {
- cout << "Error: wrong d0 parameter " << track.getIP() << '\n';
+ if(print) {
+ cout << "Error: wrong d0 parameter " << track.getIP() << '\n';
+ print--;
+ }
error++;
}
if ((track.getZ0()<m_z0_min) || (track.getZ0()>m_z0_max)) {
- cout << "Error: wrong z0 parameter " << track.getZ0() << '\n';
+ if(print) {
+ cout << "Error: wrong z0 parameter " << track.getZ0() << '\n';
+ print--;
+ }
error++;
}
if ((track.getCotTheta()<m_ctheta_min) || (track.getCotTheta()>m_ctheta_max)) {
- cout << "Error: wrong ctheta parameter " << track.getCotTheta() << '\n';
+ if(print) {
+ cout << "Error: wrong ctheta parameter " << track.getCotTheta() << '\n';
+ print--;
+ }
error++;
}
//if(error) exit(0);
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/PattMergeRootAlgo.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/PattMergeRootAlgo.cxx
index dfaab3abbcc2a7a5bec1364c6ac5b84128f7ea64..cd4dbc7b04431cecf047550ce38ea85786969a89 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/PattMergeRootAlgo.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/PattMergeRootAlgo.cxx
@@ -158,11 +158,12 @@ StatusCode PattMergeRootAlgo::RunMerging() {
void PattMergeRootAlgo::PostMessage(void) {
- if (FTKLogger::fType==0) ATH_MSG_FATAL(fBuffer->str());
- else if(FTKLogger::fType==1) ATH_MSG_ERROR(fBuffer->str());
- else if(FTKLogger::fType==2) ATH_MSG_WARNING(fBuffer->str());
- else if(FTKLogger::fType==3) ATH_MSG_INFO(fBuffer->str());
- else if(FTKLogger::fType==4) ATH_MSG_DEBUG(fBuffer->str());
+ int fType=getLoggerMsgType();
+ if (fType==0) ATH_MSG_FATAL(getLoggerMsg());
+ else if(fType==1) ATH_MSG_ERROR(getLoggerMsg());
+ else if(fType==2) ATH_MSG_WARNING(getLoggerMsg());
+ else if(fType==3) ATH_MSG_INFO(getLoggerMsg());
+ else if(fType==4) ATH_MSG_DEBUG(getLoggerMsg());
}
diff --git a/Trigger/TrigFTK/TrigFTKBankGen/src/components/TrigFTKBankGen_entries.cxx b/Trigger/TrigFTK/TrigFTKBankGen/src/components/TrigFTKBankGen_entries.cxx
index da304f527c2d2a38b1f75545f7e4cb2f7a8e0d35..d10dc998d472db342a44b628b8a4b58884aa3930 100644
--- a/Trigger/TrigFTK/TrigFTKBankGen/src/components/TrigFTKBankGen_entries.cxx
+++ b/Trigger/TrigFTK/TrigFTKBankGen/src/components/TrigFTKBankGen_entries.cxx
@@ -12,8 +12,11 @@
#include "TrigFTKBankGen/FTKCachedBankGenAlgo.h"
+#include "TrigFTKBankGen/FTKDataMatrixFromMCAlgo.h"
+
DECLARE_ALGORITHM_FACTORY(FTKBankGenAlgo)
DECLARE_ALGORITHM_FACTORY(FTKConstGenAlgo)
+DECLARE_ALGORITHM_FACTORY(FTKDataMatrixFromMCAlgo)
DECLARE_ALGORITHM_FACTORY(PattMergeRootAlgo)
@@ -24,6 +27,7 @@ DECLARE_ALGORITHM_FACTORY(FTKCachedBankGenAlgo)
DECLARE_FACTORY_ENTRIES(TrigFTKBankGen) {
DECLARE_ALGORITHM(FTKBankGenAlgo);
DECLARE_ALGORITHM(FTKConstGenAlgo);
+ DECLARE_ALGORITHM(FTKDataMatrixFromMCAlgo);
DECLARE_ALGORITHM(PattMergeRootAlgo);
diff --git a/Trigger/TrigFTK/TrigFTKSim/CMakeLists.txt b/Trigger/TrigFTK/TrigFTKSim/CMakeLists.txt
index 603f52a10e68f6fecbd7a90452eefc0960a57d57..d8c81d5b2ab21c49ddcffe55b4af421f893417c4 100644
--- a/Trigger/TrigFTK/TrigFTKSim/CMakeLists.txt
+++ b/Trigger/TrigFTK/TrigFTKSim/CMakeLists.txt
@@ -144,6 +144,11 @@ atlas_add_executable( trigftk_dataflow
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${HEPPDT_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} ${TBB_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${HEPPDT_LIBRARIES} ${EIGEN_LIBRARIES} ${HEPMC_LIBRARIES} ${TBB_LIBRARIES} AthenaBaseComps PileUpToolsLib StoreGateLib SGtests GaudiKernel InDetIdentifier InDetReadoutGeometry InDetPrepRawData Particle TrkParameters TrkParticleBase TrkTrack TrkTruthData TrkExInterfaces TrkToolInterfaces TrigInDetEvent TrigFTKPool IdDictDetDescr EventInfo NavFourMom GeneratorObjects InDetRawData InDetSimData TrkMeasurementBase TrkRIO_OnTrack TrigFTKSimLib )
+atlas_add_executable( trigftk_makeconn
+ standalone/make_conn.cc
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${HEPPDT_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} ${TBB_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${HEPPDT_LIBRARIES} ${EIGEN_LIBRARIES} ${HEPMC_LIBRARIES} ${TBB_LIBRARIES} AthenaBaseComps PileUpToolsLib StoreGateLib SGtests GaudiKernel InDetIdentifier InDetReadoutGeometry InDetPrepRawData Particle TrkParameters TrkParticleBase TrkTrack TrkTruthData TrkExInterfaces TrkToolInterfaces TrigInDetEvent TrigFTKPool IdDictDetDescr EventInfo NavFourMom GeneratorObjects InDetRawData InDetSimData TrkMeasurementBase TrkRIO_OnTrack TrigFTKSimLib )
+
# Install files from the package:
atlas_install_python_modules( python/__init__.py python/TrigFTKTruthAlgsConfig.py python/QueryFTKdb.py python/findInputs.py python/FTKSimOptions.py )
atlas_install_joboptions( share/*.py )
diff --git a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKLogging.h b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKLogging.h
index 2da782225b4ed429ad8d0f0ad2bf40ebbfd1e432..f6e9ccd6f0a4a02a094e1f02bf7d86410165762a 100644
--- a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKLogging.h
+++ b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKLogging.h
@@ -20,7 +20,7 @@
class FTKLogging {
public:
- FTKLogging(const std::string &name) : fName(name),fOut(new std::ostream(0)) { }
+ FTKLogging(const std::string &name) : m_name(name),m_out(new std::ostream(0)) { }
virtual ~FTKLogging();
std::ostream &Fatal(const char *where) const;
std::ostream &Error(const char *where) const;
@@ -28,15 +28,15 @@ public:
std::ostream &Info(const char *where) const;
std::ostream &Debug(const char *where) const;
void ShowProgress(const char *text);
- const std::string &GetName(void) const { return fName; }
+ const std::string &GetName(void) const { return m_name; }
static void SetPrintLevel(int);
static void SetAbortLevel(int);
private:
- static int sPrintLevel,sAbortLevel;
+ static int s_printLevel,s_abortLevel;
std::string GetHeader(const char *where) const;
std::ostream &GetStream(const char *name,int level) const;
- std::string fName;
- std::ostream *fOut;
+ std::string m_name;
+ std::ostream *m_out;
};
class FTKLogger : public std::streambuf {
@@ -60,18 +60,21 @@ public:
// }
static FTKLogger *GetLogger(int type,int printLevel,int abortLevel);
protected:
- static FTKLogger *logger;
+ static FTKLogger *s_logger;
FTKLogger(void);
virtual ~FTKLogger();
static void SetLogger(FTKLogger *ftkLogger);
void SetType(int type,int printLevel,int abortLevel);
+ int getLoggerMsgType(void) const { return m_type; }
+ std::string getLoggerMsg(void) const { return m_buffer->str(); }
//FTKLogger() : fLevel(0) { }
virtual std::streamsize xsputn ( const char * s, std::streamsize n );
virtual int overflow (int c);
virtual void PostMessage(void);
- std::ostringstream *fBuffer;
- int fType,fAbortLevel,fPrintLevel;
+ std::ostringstream *m_buffer;
+
+ int m_type,m_abortLevel,m_printLevel;
};
#endif
diff --git a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKPatternBySector.h b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKPatternBySector.h
index eb839a59daf8358ca307a6e17ddffc1db5a4c20d..28fcaf5e8d9cf4d3e5c71f546db6d8b854b5022e 100644
--- a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKPatternBySector.h
+++ b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTKPatternBySector.h
@@ -79,6 +79,7 @@ class FTKPatternBySectorReader : public virtual FTKPatternBySectorBase {
// read patterns (optimized IO), given a minimum number of patterns
// repeated calls to this function will result in all sectors to be read
// for sectors appended to the list, all patterns have been read
+ // i.e. the next call with return different sectors or an empty list
// (use Rewind() to start all over)
// returns: next sector to read (or -1)
virtual int ReadRaw(int firstSector,
@@ -141,7 +142,7 @@ class FTKPatternBySectorWriter : public virtual FTKPatternBySectorBase {
(TDirectory &dir,WRITERTYPE_t type=kWriterDefault);
// append patterns from file (down to minCoverage)
// the default algorithm merges one sector a time, calling
- // AppendMergedPatterns() fro each sector.
+ // AppendMergedPatterns() for each sector.
// This procedure is not optimal for all file formats
// For this reason the method is virtual, to allow for
// a more efficient implementation
diff --git a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTK_CompressedAMBank.h b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTK_CompressedAMBank.h
index f10e77d58f3a70aaf54043b7ef9e55ec970ff19e..442259f727ac1bdbf69e14e6c73d4e627a6275f6 100644
--- a/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTK_CompressedAMBank.h
+++ b/Trigger/TrigFTK/TrigFTKSim/TrigFTKSim/FTK_CompressedAMBank.h
@@ -330,6 +330,7 @@ protected:
// holds all pattern data
PatternBank m_bank;
MAP<int,std::pair<int,int> > m_SectorFirstLastPattern;
+ MAP<int,bool> m_tooFew;
//
// hold wildcards (layer mask) per sector
typedef uint8_t HitPattern_t;
diff --git a/Trigger/TrigFTK/TrigFTKSim/scripts/TrigFTKSim_tf.py b/Trigger/TrigFTK/TrigFTKSim/scripts/TrigFTKSim_tf.py
index 2642ed93ef00a4c762fcf87f5a4188b92fdc7a89..1bbbea3e3b3628464bd787fc9005e2a2f1b55860 100755
--- a/Trigger/TrigFTK/TrigFTKSim/scripts/TrigFTKSim_tf.py
+++ b/Trigger/TrigFTK/TrigFTKSim/scripts/TrigFTKSim_tf.py
@@ -80,7 +80,7 @@ def addFTKSimulationArgs(parser):
parser.add_argument('--ssmapunused_path', type=trfArgClasses.argFactory(trfArgClasses.argString, runarg=True),
help='Location of ssmapunused file', group='TrigFTKSim')
parser.add_argument('--IBLMode',type=trfArgClasses.argFactory(trfArgClasses.argInt,runarg=True),
- help='Enalbe the IBL geometry',group='TrigFTKSim')
+ help='Enable the IBL geometry',group='TrigFTKSim')
parser.add_argument('--badmap_path', type=trfArgClasses.argFactory(trfArgClasses.argString, runarg=True),
help='Location of badmap file', group='TrigFTKSim')
parser.add_argument('--badmap_path_for_hit', type=trfArgClasses.argFactory(trfArgClasses.argString, runarg=True),
diff --git a/Trigger/TrigFTK/TrigFTKSim/src/FTKLogging.cxx b/Trigger/TrigFTK/TrigFTKSim/src/FTKLogging.cxx
index 35c791a187f5a761d382627f0e9c8770d947b398..4b9ae5e372a8a83d9d6cdf642ef7af2ec3e13569 100644
--- a/Trigger/TrigFTK/TrigFTKSim/src/FTKLogging.cxx
+++ b/Trigger/TrigFTK/TrigFTKSim/src/FTKLogging.cxx
@@ -18,26 +18,26 @@
*
*/
-int FTKLogging::sPrintLevel=4;
+int FTKLogging::s_printLevel=4;
-int FTKLogging::sAbortLevel=2;
+int FTKLogging::s_abortLevel=2;
void FTKLogging::SetPrintLevel(int level) {
- sPrintLevel=level;
+ s_printLevel=level;
}
void FTKLogging::SetAbortLevel(int level) {
- sAbortLevel=level;
+ s_abortLevel=level;
}
std::string FTKLogging::GetHeader(const char *where) const {
- return fName+"::"+where+" ";
+ return m_name+"::"+where+" ";
}
std::ostream &FTKLogging::GetStream(const char *name,int level) const {
- fOut->rdbuf(FTKLogger::GetLogger(level,sPrintLevel,sAbortLevel));
- if(name) (*fOut)<<GetHeader(name);
- return *fOut;
+ m_out->rdbuf(FTKLogger::GetLogger(level,s_printLevel,s_abortLevel));
+ if(name) (*m_out)<<GetHeader(name);
+ return *m_out;
}
@@ -66,7 +66,7 @@ void FTKLogging::ShowProgress(const char *text) {
}
FTKLogging::~FTKLogging() {
- if (fOut) delete fOut;
+ if (m_out) delete m_out;
}
/* class FTKLogger
@@ -79,33 +79,33 @@ FTKLogging::~FTKLogging() {
*
*/
-FTKLogger *FTKLogger::logger=0;
+FTKLogger *FTKLogger::s_logger=0;
FTKLogger::FTKLogger(void) {
- fBuffer=0;
- fType=0;
- fAbortLevel = 0;
- fPrintLevel = 0;
+ m_buffer=0;
+ m_type=0;
+ m_abortLevel = 0;
+ m_printLevel = 0;
}
FTKLogger::~FTKLogger() {
- if(logger==this) SetLogger(0);
+ if(s_logger==this) SetLogger(0);
}
FTKLogger *FTKLogger::GetLogger(int type,int printLevel,int abortLevel) {
- if(!logger) SetLogger(new FTKLogger());
- logger->SetType(type,printLevel,abortLevel);
- return logger;
+ if(!s_logger) SetLogger(new FTKLogger());
+ s_logger->SetType(type,printLevel,abortLevel);
+ return s_logger;
}
void FTKLogger::SetLogger(FTKLogger *ftkLogger) {
- logger=ftkLogger;
+ s_logger=ftkLogger;
}
void FTKLogger::SetType(int type,int printLevel,int abortLevel) {
- fType=type;
- fPrintLevel=printLevel;
- fAbortLevel=abortLevel;
+ m_type=type;
+ m_printLevel=printLevel;
+ m_abortLevel=abortLevel;
}
int FTKLogger::overflow (int c) {
@@ -116,13 +116,13 @@ int FTKLogger::overflow (int c) {
std::streamsize FTKLogger::xsputn(const char * s, std::streamsize n) {
for(int i=0;i<n;i++) {
- if(!fBuffer) fBuffer=new std::ostringstream();
+ if(!m_buffer) m_buffer=new std::ostringstream();
if(s[i]=='\n') {
PostMessage();
- delete fBuffer;
- fBuffer=0;
+ delete m_buffer;
+ m_buffer=0;
} else {
- (*fBuffer)<<s[i];
+ (*m_buffer)<<s[i];
}
}
return n;
@@ -130,10 +130,10 @@ std::streamsize FTKLogger::xsputn(const char * s, std::streamsize n) {
void FTKLogger::PostMessage(void) {
static char const *text[]={"FATAL","ERROR","WARNING","INFO","DEBUG"};
- if(fType<fPrintLevel) {
- std::cout<<text[fType]<<"\t"<<fBuffer->str()<<"\n";
+ if(m_type<m_printLevel) {
+ std::cout<<text[m_type]<<"\t"<<m_buffer->str()<<"\n";
}
- if(fType<fAbortLevel) {
- exit(1+fType);
+ if(m_type<m_abortLevel) {
+ exit(1+m_type);
}
}
diff --git a/Trigger/TrigFTK/TrigFTKSim/src/FTKRoadFinderAlgo.cxx b/Trigger/TrigFTK/TrigFTKSim/src/FTKRoadFinderAlgo.cxx
index 987162da4273150c94cc54727bd4a619166b78ef..3d85f2623549569d8c0e56db2f1c3ac8c6e19576 100644
--- a/Trigger/TrigFTK/TrigFTKSim/src/FTKRoadFinderAlgo.cxx
+++ b/Trigger/TrigFTK/TrigFTKSim/src/FTKRoadFinderAlgo.cxx
@@ -626,10 +626,11 @@ StatusCode FTKRoadFinderAlgo::finalize() {
}
void FTKRoadFinderAlgo::PostMessage() {
- if(fType==0) ATH_MSG_FATAL(fBuffer->str());
- else if(fType==1) ATH_MSG_ERROR(fBuffer->str());
- else if(fType==2) ATH_MSG_WARNING(fBuffer->str());
- else if(fType==3) ATH_MSG_INFO(fBuffer->str());
- else if(fType==4) ATH_MSG_DEBUG(fBuffer->str());
+ int fType=getLoggerMsgType();
+ if(fType==0) ATH_MSG_FATAL(getLoggerMsg());
+ else if(fType==1) ATH_MSG_ERROR(getLoggerMsg());
+ else if(fType==2) ATH_MSG_WARNING(getLoggerMsg());
+ else if(fType==3) ATH_MSG_INFO(getLoggerMsg());
+ else if(fType==4) ATH_MSG_DEBUG(getLoggerMsg());
}
diff --git a/Trigger/TrigFTK/TrigFTKSim/src/FTK_AMsimulation_base.cxx b/Trigger/TrigFTK/TrigFTKSim/src/FTK_AMsimulation_base.cxx
index d20ec8e22fe15bf00eca16fb17908f77756ba11e..f51416c93d8bdd99502fa97de95d972c9c80633e 100644
--- a/Trigger/TrigFTK/TrigFTKSim/src/FTK_AMsimulation_base.cxx
+++ b/Trigger/TrigFTK/TrigFTKSim/src/FTK_AMsimulation_base.cxx
@@ -9,6 +9,8 @@
//#define VERBOSE_DEBUG
+#define PRINT_IBL_HITS 0
+
using namespace std;
/* class FTK_AMsimulation_base
@@ -201,7 +203,7 @@ int FTK_AMsimulation_base::passHitsUnused(const std::vector<FTKHit> &hitlist) {
std::map<int,FTKSS>::iterator issitem = ssmap.find(ssid);
if (issitem==ssmap.end()) {
- // ad the item related to this SS number
+ // add the item related to this SS number
ssmap[ssid] = FTKSS();
issitem = ssmap.find(ssid);
}
@@ -212,6 +214,43 @@ int FTK_AMsimulation_base::passHitsUnused(const std::vector<FTKHit> &hitlist) {
if(error) {
FTKSetup::PrintMessageFmt(ftk::warn,"FTK_AMsimulation_base::passHitsUnused number or errors=%d\n",error );
}
+ static int print=PRINT_IBL_HITS;
+ if(print) {
+ cout<<"PassHitsUnused\n";
+ for(std::map< int, std::map< int, FTKSS > >::const_iterator
+ ipl=m_usedssmap_ignored.begin();
+ ipl!=m_usedssmap_ignored.end();ipl++) {
+ cout<<"unused plane="<<(*ipl).first<<" nSS="<<(*ipl).second.size();
+ int nhit=0;
+ int lastHash=-1;
+ for(std::map< int, FTKSS >::const_iterator ihit=(*ipl).second.begin();
+ ihit!=(*ipl).second.end();ihit++) {
+ nhit+= (*ihit).second.getNHits();
+ if((*ipl).first==0) {
+ int localID,localX,localY;
+ getSSMapUnused()->decodeSSTowerXY((*ihit).first,
+ getBankID(),(*ipl).first,0,
+ localID,localX,localY,true);
+ int hash=getSSMapUnused()->getRegionMap()
+ ->getGlobalId(getBankID(),(*ipl).first,localID);
+ if(hash!=lastHash) {
+ cout<<"\nModule="<<hash;
+ lastHash=hash;
+ }
+ cout<<" "<<(*ihit).first;
+ //<<":"<<localX<<","<<localY;
+ for(int ih=0;ih<(*ihit).second.getNHits();ih++) {
+ const FTKHit &h=(*ihit).second.getHit(ih);
+ cout<<"["<<h[0]<<","<<h[1]<<"]";
+ }
+ }
+ }
+ cout<<" nhit="<<nhit;
+ cout<<"\n";
+ }
+ print--;
+ }
+
return res;
}
diff --git a/Trigger/TrigFTK/TrigFTKSim/src/FTK_CompressedAMBank.cxx b/Trigger/TrigFTK/TrigFTKSim/src/FTK_CompressedAMBank.cxx
index 2563a3eb40e0fce4b2df08b9ff4ef221e24e5c86..b89af49d02fe6ef7fe0acf94014e4d27e151e430 100644
--- a/Trigger/TrigFTK/TrigFTKSim/src/FTK_CompressedAMBank.cxx
+++ b/Trigger/TrigFTK/TrigFTKSim/src/FTK_CompressedAMBank.cxx
@@ -604,7 +604,7 @@ int FTK_CompressedAMBank::writePCachedBankFile
unsigned int ndc[8]; // number of DC bits per plane
unsigned int wildcard;
//unsigned lowestSSID;
- int sector,firstPattern,lastPattern,isub,lamb;
+ int sector,firstPattern,lastPattern,isub,lamb,allPatternsUsed;
auxtree->Branch("nplanes",&nplanes,"nplanes/I");
auxtree->Branch("ndc",ndc,"ndc[nplanes]/I");
auxtree->Branch("sector",§or,"sector/I");
@@ -613,6 +613,7 @@ int FTK_CompressedAMBank::writePCachedBankFile
auxtree->Branch("lamb",&lamb,"lamb/I");
auxtree->Branch("isub",&isub,"isub/I");
auxtree->Branch("wildcard",&wildcard,"wildcard/I");
+ auxtree->Branch("allPatternsUsed",&allPatternsUsed,"allPatternsUsed/I");
//auxtree->Branch("lowestSSID",&lowestSSID,"lowestSSID/I");
// set number of DC bits per plane (constant)
for(int i=0;i<getNPlanes();i++) {
@@ -627,6 +628,11 @@ int FTK_CompressedAMBank::writePCachedBankFile
isub= nsub ? (sector % nsub) : 0;
lamb= nlamb ? (sector % nlamb) : 0;
wildcard=m_SectorWC[sector];
+ MAP<int,bool>::const_iterator tooFew=m_tooFew.find(sector);
+ allPatternsUsed=-1;
+ if(tooFew!=m_tooFew.end()) {
+ allPatternsUsed=(*tooFew).second ? 1 : 0;
+ }
//lowestSSID=;
auxtree->Fill();
}
@@ -1442,7 +1448,7 @@ void FTK_CompressedAMBank::importDCpatterns
// here: top_sectorByLamb has all bits set to one and is larger or equal the largest sector/nlamb
unsigned leadingBitMask=(top_sectorByNlamb>>1)+1;
unsigned sectorByLamb=0;
- while(true) {
+ for(unsigned count=0;count<=top_sectorByNlamb;count++) {
//for(unsigned sector=ilamb;sector<sectorPointer.size();sector+=nlamb) {
unsigned sector=sectorByLamb*nlamb+ilamb;
@@ -1459,7 +1465,6 @@ void FTK_CompressedAMBank::importDCpatterns
if(sectorByLamb & mask) break;
// if the bit is zero, have to count the next bit
}
- if(!sectorByLamb) break; // all bits were counted, stop
// when counting in reverse bit order, sectors are not produced in sequence
// high (invalid) sector numbers are produced in random order and have
// to be skipped
@@ -1694,6 +1699,7 @@ void FTK_CompressedAMBank::importDCpatterns
Debug("importDCpatterns")
<<"sector "<<sector<<" nPatt="<<ipatt<<"\n";
}
+ if(!sectorByLamb) break; // all bits were counted, stop
} // end loop over sectors
} // end loop over subregions
}
@@ -3549,7 +3555,6 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
for(std::list<FTK_CompressedAMBank::Partition>::const_iterator
iPartition=partitionList.begin();iPartition!=partitionList.end();
iPartition++) {
- maxpatts+=(*iPartition).fNumPatternMax;
// the FTKPatternBySectorBlockReader is reading the patterns
// sector-by-sector, in a pre-defined coverage range
//
@@ -3569,7 +3574,8 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
Info("readSectorOrderedBank")
<<"reading from "<<name<<" nLayer="<<TSPreader->GetNLayers()
<<" partition "<<partitionCount<<" maxpattern="
- <<(*iPartition).fNumPatternMax<<" maxSector="<<(*iPartition).fNumSectorMax
+ <<(*iPartition).fNumPatternMax
+ <<" maxSector="<<(*iPartition).fNumSectorMax
<<"\n";
// the coverage map holds for each coverage the number of patterns
// it is used in order to estimate down to which coverage the patterns
@@ -3582,7 +3588,6 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
#ifdef SEARCH_MEMORY_LEAK
printVmemUsage("after GetNPatternsByCoverage");
#endif
- std::map<int,int>::const_reverse_iterator i=coverageMap.rbegin();
TSPreader->Rewind();
#ifdef SEARCH_MEMORY_LEAK
printVmemUsage("after Rewind");
@@ -3592,6 +3597,7 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
uint32_t totalPatterns=0;
uint64_t totalPatternsCoverage=0;
int nDC0=nDC;
+ maxpatts=nDC0+(*iPartition).fNumPatternMax;
for(std::map<int,int>::const_iterator j=coverageMap.begin();
j!=coverageMap.end();j++) {
totalPatterns += (*j).second;
@@ -3603,27 +3609,26 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
<<" average coverage: "
<<(totalPatterns ? (totalPatternsCoverage/(double)totalPatterns):0)
<<"\n";
- int covEnd;
std::set<int> loadedSectorList;
int maxNumSector=(*iPartition).fNumSectorMax;
if(maxNumSector<0) maxNumSector=(*iPartition).fSectorSet.size();
- do {
- std::map<int,int>::const_reverse_iterator j=i;
- int npattEstimate=nDC;
- covEnd=(*i).first;
- while((j!=coverageMap.rend())&&(npattEstimate<maxpatts)) {
- npattEstimate += (*j).second;
- covEnd=(*j).first;
- j++;
- }
- Info("readSectorOrderedBank")
- <<"adding coverage range "<<(*i).first<<"-"<<covEnd
- <<" extra patterns="<<npattEstimate-nDC<<"\n";
+ // loop over coverage ranges
+ // starting with coverage indicated by (*i)
+ bool tooFew=true;
+ int covEnd=0;
+ for(std::map<int,int>::const_reverse_iterator i=coverageMap.rbegin();
+ i!=coverageMap.rend();) {
+ // read a minimum of one coverages
+ int nAvail=0;
+ int covBegin=(*i).first;
+ do {
+ nAvail +=(*i).second;
+ covEnd=(*i).first;
+ i++;
+ } while((nDC+nAvail<=maxpatts)&&(i!=coverageMap.rend()));
// read all patterns in the given coverage range
- // (not reading the smallest coverage)
- // these patterns for sure will fit into the DC bank
- //for(int sector=TSPreader->GetFirstSector();sector>=0;
- // sector=TSPreader->GetNextSector(sector)) {
+ int nTSPold=nTSP;
+ int nDCold=nDC;
for(std::set<int>::const_iterator sectorPtr=
(*iPartition).fSectorSet.begin();
sectorPtr!=(*iPartition).fSectorSet.end();sectorPtr++) {
@@ -3632,42 +3637,46 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
if((int)loadedSectorList.size()>=maxNumSector) continue;
loadedSectorList.insert(sector);
}
- FTKPatternOneSector *patterns=TSPreader->Read(sector,covEnd+1);
+ FTKPatternOneSector *patterns=
+ TSPreader->Read(sector,covEnd);
if(patterns) {
insertPatterns(sector,patterns,maxpatts,dcPatterns,nDC,nTSP);
delete patterns;
}
+ if(nDC>=maxpatts) break;
}
- // read patterns for the smallest coverage
- // these patterns will partially fit into the bank
- for(std::set<int>::const_iterator sectorPtr=
- (*iPartition).fSectorSet.begin();
- sectorPtr!=(*iPartition).fSectorSet.end();sectorPtr++) {
- int sector=*sectorPtr;
- if(loadedSectorList.find(sector)==loadedSectorList.end()) {
- if((int)loadedSectorList.size()>=maxNumSector) continue;
- loadedSectorList.insert(sector);
- }
- FTKPatternOneSector *patterns=TSPreader->Read(sector,covEnd);
- if(patterns) {
- insertPatterns(sector,patterns,maxpatts,dcPatterns,nDC,nTSP);
- delete patterns;
- if(nDC>=maxpatts) break;
+ Info("readSectorOrderedBank")
+ <<"added coverage range "<<covBegin<<"-"<<covEnd
+ <<" DC added="<<nDC-nDCold
+ <<" TSP used="<<nTSP-nTSPold<<"/"<<nAvail
+ <<"\n";
+ if(nDC>=maxpatts){
+ if((nTSP-nTSPold<nAvail)||(i!=coverageMap.rend())) {
+ tooFew=false;
}
+ break;
}
- i=j;
- // stop if all patterns have been read
+ // stop if all patterns have been read
// or if the maximum number of patterns is reached
- // DEBUG - speed up parttern reading
+ // for DEBUGGING - speed up pattern reading
// break;
- } while((i!=coverageMap.rend())&&(nDC<maxpatts));
+ }
+ for(std::set<int>::const_iterator sectorPtr=
+ (*iPartition).fSectorSet.begin();
+ sectorPtr!=(*iPartition).fSectorSet.end();sectorPtr++) {
+ int sector=*sectorPtr;
+ m_tooFew[sector]=tooFew;
+ }
Info("readSectorOrderedBank")
<<"partition "<<partitionCount<<" number of DC patterns: "<<nDC-nDC0
+ <<"/"<<maxpatts-nDC0
<<" smallest coverage="<<covEnd
<<" #sectors="<<loadedSectorList.size()
<<" maxnsector="<<maxNumSector
+ <<" tooFew="<<tooFew
+ <<" nDC="<<nDC
<<"\n";
partitionCount++;
sectorCount += loadedSectorList.size();
@@ -3695,6 +3704,7 @@ int FTK_CompressedAMBank::readPartitionedSectorOrderedBank
pos += dcPatterns[sector].size()*offsetSSID;
numPattern+=dcPatterns[sector].size();
}
+ Info("readSectorOrderedBank")<<"numPattern="<<numPattern<<"\n";
int32_t *ssidData=new int32_t[pos];
for(unsigned sector=0;sector<dcPatterns.size();sector++) {
HitPatternMap_t const &hitMap=dcPatterns[sector];
diff --git a/Trigger/TrigFTK/TrigFTKSim/src/TrackFitter711.cxx b/Trigger/TrigFTK/TrigFTKSim/src/TrackFitter711.cxx
index 852620854d901943e36167bfdf99993e17954e3c..1d56055239b404fde45b306340d222658a58bfbc 100644
--- a/Trigger/TrigFTK/TrigFTKSim/src/TrackFitter711.cxx
+++ b/Trigger/TrigFTK/TrigFTKSim/src/TrackFitter711.cxx
@@ -17,6 +17,8 @@
#include <cassert>
using namespace std;
+static int print_guessed=0; //1000;
+
TrackFitter711::TrackFitter711() :
TrackFitter(),
m_use_guessing(true),
@@ -619,6 +621,7 @@ void TrackFitter711::processor_end(int ibank)
// m_trackoutput_pre_hw->addNFitsHWRejectedMajorityI(ibank,m_nfits_rejmajI);
// m_trackoutput_pre_hw->addNConnections(ibank,m_nconn);
// m_trackoutput_pre_hw->addNExtrapolatedTracks(ibank,m_nextrapolatedTracks);
+
}
@@ -3366,7 +3369,10 @@ void TrackFitter711::guessedHitsToSSID() {
m_ssid[ip] = m_ssmap_complete->getSSGlobal(tmphit, 1);
}
-
+ if(print_guessed && (m_idplanes_eff[ip]==0)) {
+ cout<<"guessed_hit: "<< m_ssid[ip]<<":"<<moduleID
+ <<","<<tmphit[0]<<","<<tmphit[1]<<"\n";
+ }
// std::cout << "SSID current " << m_ssid[ip] << std::endl;
diff --git a/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.cc b/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.cc
index 7b2b73706475a0abdc9e5224a9e4ea67cecd9239..a641ee76f03cd61e217ebe98bd6e52cea6c1d44f 100644
--- a/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.cc
+++ b/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.cc
@@ -72,6 +72,14 @@ void Init() {
// create instances for the FTK generic distributions
histocoordmasketa_ftk = new TH2F("histocoordmasketa_ftk","Track probability of using a coordinate;Coordinate bit;#eta;Prob.",16,0,16,nbins,etamin,etamax);
histocoordmaskphi_ftk = new TH2F("histocoordmaskphi_ftk","Track probability of using a coordinate;Coordinate bit;#phi;Prob.",16,0,16,nbins,phimin,phimax);
+ for(int i=0;i<=16;i++) {
+ histocoordmasketaphi_ftk[i] = new TH2F
+ (TString::Format("histocoordmasketaphiL%d_ftk",i),"Track probability of using a coordinate;#eta;#phi;Prob.",32,etamin,etamax,32,phimin,phimax);
+ histoeff_etaphiplane[i]= new TH2F
+ (TString::Format("histoeff_etaphiL%d",i),
+ TString::Format("efficiency plane %d;#eta;#phi;Prob.",i),
+ 32,etamin,etamax,32,phimin,phimax);
+ }
histocoordmaskz0_ftk = new TH2F("histocoordmaskz0_ftk","Track probability of using a coordinate;Coordinate bit;#z0;Prob.",16,0,16,nbins,z0min, z0max);
histonmisseta_ftk = new TH2F("histonmisseta_ftk","Track probability of missing N layers;N Miss;#eta;Prob.",6,0,6,nbins,etamin,etamax);
histochisqndfeta_ftk = new TH2F("histochisqndfeta_ftk","Track probability of #chi^{2};#chi^{2}/ndf;#eta;Prob.",20,0,20,nbins,etamin,etamax);
@@ -177,6 +185,12 @@ void Init() {
histod0res_veta = new TProfile("histod0res_veta","#Delta d_{0} (mm);#eta", nbins, etamin, etamax);
histoz0res_veta = new TProfile("histoz0res_veta","#Delta z_{0} (mm);eta", nbins, etamin, etamax);
histod0res_vphi = new TProfile("histod0res_vphi","#Delta d_{0} (mm);#phi", nbins, phimin, phimax);
+ histo2d0res_vphi = new TH2F("histo2d0res_vphi","#Delta d_{0} (mm);#phi", nbins, phimin, phimax,50,-5.*Dd0,5.*Dd0);
+ histo2d0res_veta = new TH2F("histo2d0res_veta","#Delta d_{0} (mm);#eta", nbins, etamin, etamax,50,-5.*Dd0,5.*Dd0);
+ histo2d0res_coordbit = new TH2F("histo2d0res_coordbit","#Delta d_{0} (mm);coordinate bit", 16, -0.5, 15.5,50,-5.*Dd0,5.*Dd0);
+ histo2d0truth_vphi = new TH2F("histo2d0truth_vphi",";#phi;truth d_{0} (mm)", nbins, phimin, phimax,50,-4.,4.);
+ histo2d0ftk_vphi = new TH2F("histo2d0ftk_vphi",";#phi;FTK d_{0} (mm)", nbins, phimin, phimax,50,-4.,4.);
+ histo2curvCurv = new TH2F("histo2curvCurv",";curv;curv",50,-abscurvmax ,abscurvmax ,50,-abscurvmax,abscurvmax);
histoz0res_vphi = new TProfile("histoz0res_vphi","#Delta z_{0} (mm);#phi", nbins, phimin, phimax);
histod0res_vz0 = new TProfile("histod0res_vz0","#Delta d_{0} (mm);z_{0} (mm)", nbins, z0min, z0max);
histoz0res_vz0 = new TProfile("histoz0res_vz0","#Delta z_{0} (mm);z_{0} (mm)", nbins, z0min, z0max);
@@ -263,21 +277,15 @@ void Process(Long64_t ientry) {
curtrack = (Use1stStage ? tracks->getTrackI(itrk) : tracks->getTrack(itrk));
//if ( TMath::Abs(curtrack->getPt())*1e-3 < ptmincut ) continue;
+
+ //double temp_phi = curtrack->getPhi();
+ //double temp_d0 = curtrack->getIP();
+ //double temp_z0 = curtrack->getZ0();
- double temp_phi = curtrack->getPhi();
- double temp_d0 = curtrack->getIP();
- double temp_z0 = curtrack->getZ0();
-
- // double dx = -0.42;
- // double dy = -0.53;
-
-
- // double dx = -0.5491;
- // double dy = -0.6557;
-
- double thisd0 = temp_d0 + dx*sin(temp_phi)-dy*cos(temp_phi);
- double thisz0 = temp_z0 + ((cos(temp_phi) *dx - sin(temp_phi)*dy))*curtrack->getCotTheta();
- double thisphi = temp_phi;
+ double thisd0 = curtrack->getIP();
+ //temp_d0 + dx*sin(temp_phi)-dy*cos(temp_phi);
+ double thisz0 = curtrack->getZ0();
+ //temp_z0 + ((cos(temp_phi) *dx - sin(temp_phi)*dy))*curtrack->getCotTheta();
histod0_ftk->Fill(thisd0);
@@ -356,10 +364,14 @@ void Process(Long64_t ientry) {
if (curtrack->getBitmask()&(1<<6)) histophiz0_ftk_PixL2->Fill(curtrack->getPhi(),thisz0);
+ histocoordmasketaphi_ftk[16]
+ ->Fill(curtrack->getEta(),curtrack->getPhi());
for (Int_t icoord=0;icoord!=curtrack->getNCoords();++icoord) {
if (curtrack->getBitmask()&(1<<icoord)){
histocoordmasketa_ftk->Fill(icoord,curtrack->getEta());
+ histocoordmasketaphi_ftk[icoord]
+ ->Fill(curtrack->getEta(),curtrack->getPhi());
histocoordmaskz0_ftk->Fill(icoord,thisz0);
histocoordmaskphi_ftk->Fill(icoord,curtrack->getPhi());
// std::cout << "icoord, eta, z0, phi: " << icoord << ", " << curtrack->getEta() << ", " << thisz0 << ", " << curtrack->getPhi() << std::endl;
@@ -483,9 +495,10 @@ void Process(Long64_t ientry) {
vector<FTKTruthTrack>::const_iterator itr = truthTracks->begin();
vector<FTKTruthTrack>::const_iterator itrE = truthTracks->end();
for (;itr!=itrE;++itr) { // loop over the truth tracks
- const FTKTruthTrack &curtruth = (*itr);
+ const FTKTruthTrack &curtruth_fromFile = (*itr);
+ // apply vertex shift
- int barcode = curtruth.getBarcode();
+ int barcode = curtruth_fromFile.getBarcode();
if (barcode>100000 || barcode==0) continue;
// select only good truth tracks
@@ -493,88 +506,92 @@ void Process(Long64_t ientry) {
// double dy = curtruth.getY();
// double d0 = -1e-3*(dx*py-dy*px)*invpt;
- double px = curtruth.getPX();
- double py = curtruth.getPY();
- double pt = 1e-3*TMath::Sqrt(px*px+py*py);
- double invpt = 1./pt;
-
- if ( pt < ptmincut ) continue;
-
- double d0 = curtruth.getD0();
- if (d0<d0min || d0>d0max) continue;
-
- double z0 = curtruth.getZ();
- if (z0<z0min || z0>z0max) continue;
+ double px_truth = curtruth_fromFile.getPX();
+ double py_truth = curtruth_fromFile.getPY();
+ double pt_truth = 1e-3*TMath::Hypot(px_truth,py_truth);
+ double invpt_truth = 1./pt_truth;
+ double d0_truth = curtruth_fromFile.getD0();
+ double z0_truth = curtruth_fromFile.getZ();
+ double curv_truth = .5*curtruth_fromFile.getQ()*invpt_truth;
+ double phi_truth = curtruth_fromFile.getPhi();
+ double eta_truth = curtruth_fromFile.getEta();
+ double cotTheta_truth= pt_truth/curtruth_fromFile.getPZ();
+ int pdgcode = curtruth_fromFile.getPDGCode();
+ int eventIndex_truth=curtruth_fromFile.getEventIndex();
+
+ if (d0_truth<d0min || d0_truth>d0max) continue;
+ if (z0_truth<z0min || z0_truth>z0max) continue;
+
+ if ( pt_truth < ptmincut ) continue;
+ if (curv_truth<-abscurvmax || curv_truth>abscurvmax) continue;
+ if (phi_truth<phimin || phi_truth>phimax) continue;
+ if (eta_truth<etamin || eta_truth>etamax) continue;
+ if (eventIndex_truth!=0 && curtruth_fromFile.getQ()==0) continue;
+
+ // transform variables to shifted vertex
+ // only d0,z0 is corrected here (lowest order corrections)
+ d0_truth -= vtxTruth[0]*sin(phi_truth)-vtxTruth[1]*cos(phi_truth);
+ z0_truth -= vtxTruth[2]+((cos(phi_truth) *vtxTruth[0] - sin(phi_truth)*vtxTruth[1]))*cotTheta_truth;
- double curv = .5*curtruth.getQ()*invpt;
- if (curv<-abscurvmax || curv>abscurvmax) continue;
-
- double phi = curtruth.getPhi();
- if (phi<phimin || phi>phimax) continue;
-
- double eta = curtruth.getEta();
- if (eta<etamin || eta>etamax) continue;
-
- if (curtruth.getEventIndex()!=0 && curtruth.getQ()==0) continue;
ntruth_good += 1;
// Fill the histogram for the generic truth distribution
- histod0_truth->Fill(d0);
- histoz0_truth->Fill(z0);
- histocurv_truth->Fill(curv);
- histophi_truth->Fill(curtruth.getPhi());
- histoetaphi_truth->Fill(eta, curtruth.getPhi());
- histoetaz0_truth->Fill(eta, z0);
- histoeta_truth->Fill(eta);
- histoetaabs_truth->Fill(fabs(eta));
+ histod0_truth->Fill(d0_truth);
+ histoz0_truth->Fill(z0_truth);
+ histocurv_truth->Fill(curv_truth);
+ histophi_truth->Fill(phi_truth);
+ histoetaphi_truth->Fill(eta_truth, phi_truth);
+ histoetaz0_truth->Fill(eta_truth, z0_truth);
+ histoeta_truth->Fill(eta_truth);
+ histoetaabs_truth->Fill(fabs(eta_truth));
histoeff_truth->Fill(0);
- histopt_truth->Fill(pt);
- histopt_truth_lg->Fill(pt);
- histopt_truthlo_lg->Fill(pt);
+ histopt_truth->Fill(pt_truth);
+ histopt_truth_lg->Fill(pt_truth);
+ histopt_truthlo_lg->Fill(pt_truth);
- int pdgcode = curtruth.getPDGCode();
if (pdgcode==13 || pdgcode==-13) {
// muon block
ntruth_good_muon += 1;
- histod0_truth_muon->Fill(d0);
- histoz0_truth_muon->Fill(z0);
- histocurv_truth_muon->Fill(curv);
- histophi_truth_muon->Fill(curtruth.getPhi());
- histoeta_truth_muon->Fill(eta);
- histopt_truth_muon->Fill(pt);
- histopt_truth_muon_lg->Fill(pt);
- histopt_truth_muonlo_lg->Fill(pt);
+ histod0_truth_muon->Fill(d0_truth);
+ histoz0_truth_muon->Fill(z0_truth);
+ histocurv_truth_muon->Fill(curv_truth);
+ histophi_truth_muon->Fill(phi_truth);
+ histoeta_truth_muon->Fill(eta_truth);
+ histopt_truth_muon->Fill(pt_truth);
+ histopt_truth_muon_lg->Fill(pt_truth);
+ histopt_truth_muonlo_lg->Fill(pt_truth);
}
// match the barcode and event index values
- MatchInfo reftruth(barcode,curtruth.getEventIndex());
+ MatchInfo reftruth(barcode,eventIndex_truth);
pair<FTKBarcodeMM::const_iterator,FTKBarcodeMM::const_iterator> mrange = ftkmatchinfo.equal_range(reftruth);
+ int bitmask=0;
if (mrange.first != mrange.second) {
- histod0_truthM->Fill(d0);
- histoz0_truthM->Fill(z0);
- histocurv_truthM->Fill(curv);
- histophi_truthM->Fill(curtruth.getPhi());
- histoetaphi_truthM->Fill(eta, curtruth.getPhi());
- histoetaz0_truthM->Fill(eta, z0);
- histoeta_truthM->Fill(eta);
- histoetaabs_truthM->Fill(fabs(eta));
+ histod0_truthM->Fill(d0_truth);
+ histoz0_truthM->Fill(z0_truth);
+ histocurv_truthM->Fill(curv_truth);
+ histophi_truthM->Fill(phi_truth);
+ histoetaphi_truthM->Fill(eta_truth, phi_truth);
+ histoetaz0_truthM->Fill(eta_truth, z0_truth);
+ histoeta_truthM->Fill(eta_truth);
+ histoetaabs_truthM->Fill(fabs(eta_truth));
histoeff_truthM->Fill(0);
- histopt_truthM->Fill(pt);
- histopt_truthMlo_lg->Fill(pt);
- histopt_truthM_lg->Fill(pt);
+ histopt_truthM->Fill(pt_truth);
+ histopt_truthMlo_lg->Fill(pt_truth);
+ histopt_truthM_lg->Fill(pt_truth);
if (pdgcode==13 || pdgcode==-13) {
// matched muon block
- histod0_truthM_muon->Fill(d0);
- histoz0_truthM_muon->Fill(z0);
- histocurv_truthM_muon->Fill(curv);
- histophi_truthM_muon->Fill(curtruth.getPhi());
- histoeta_truthM_muon->Fill(eta);
- histopt_truthM_muon->Fill(pt);
- histopt_truthM_muon_lg->Fill(pt);
- histopt_truthM_muonlo_lg->Fill(pt);
+ histod0_truthM_muon->Fill(d0_truth);
+ histoz0_truthM_muon->Fill(z0_truth);
+ histocurv_truthM_muon->Fill(curv_truth);
+ histophi_truthM_muon->Fill(phi_truth);
+ histoeta_truthM_muon->Fill(eta_truth);
+ histopt_truthM_muon->Fill(pt_truth);
+ histopt_truthM_muon_lg->Fill(pt_truth);
+ histopt_truthM_muonlo_lg->Fill(pt_truth);
}
const FTKTrack *bestftk(0x0);
@@ -587,22 +604,43 @@ void Process(Long64_t ientry) {
}
if (bestftk) {
- histod0res->Fill(d0-bestftk->getIP());
- histoz0res->Fill(z0-bestftk->getZ0());
- histod0res_veta->Fill(eta,d0-bestftk->getIP());
- histoz0res_veta->Fill(eta,z0-bestftk->getZ0());
- histod0res_vphi->Fill(curtruth.getPhi(),d0-bestftk->getIP());
- histoz0res_vphi->Fill(curtruth.getPhi(),z0-bestftk->getZ0());
- histod0res_vz0->Fill(z0,d0-bestftk->getIP());
- histoz0res_vz0->Fill(z0,z0-bestftk->getZ0());
-
- histocurvres->Fill(curv-bestftk->getHalfInvPt()*1e3);
- histophires->Fill(curtruth.getPhi()-bestftk->getPhi());
- histoetares->Fill(eta-bestftk->getEta());
- historelptres->Fill((pt-TMath::Abs(bestftk->getPt())*1e-3)*invpt);
- historelptrespt->Fill(pt,(pt-TMath::Abs(bestftk->getPt())*1e-3)*invpt);
+ histod0res->Fill(d0_truth-bestftk->getIP());
+ histoz0res->Fill(z0_truth-bestftk->getZ0());
+ histod0res_veta->Fill(eta_truth,d0_truth-bestftk->getIP());
+ histoz0res_veta->Fill(eta_truth,z0_truth-bestftk->getZ0());
+ histod0res_vphi->Fill(phi_truth,d0_truth-bestftk->getIP());
+ histo2d0res_vphi->Fill(phi_truth,d0_truth-bestftk->getIP());
+ histo2d0res_veta->Fill(eta_truth,d0_truth-bestftk->getIP());
+ histo2d0truth_vphi->Fill(phi_truth,d0_truth);
+ histo2d0ftk_vphi->Fill(phi_truth,bestftk->getIP());
+ histo2curvCurv->Fill(curv_truth,bestftk->getHalfInvPt()*1.E3);
+ histoz0res_vphi->Fill(phi_truth,z0_truth-bestftk->getZ0());
+ histod0res_vz0->Fill(z0_truth,d0_truth-bestftk->getIP());
+ histoz0res_vz0->Fill(z0_truth,z0_truth-bestftk->getZ0());
+
+ histocurvres->Fill(curv_truth-bestftk->getHalfInvPt()*1e3);
+ histophires->Fill(phi_truth-bestftk->getPhi());
+ histoetares->Fill(eta_truth-bestftk->getEta());
+ historelptres->Fill((pt_truth-TMath::Abs(bestftk->getPt())*1e-3)*invpt_truth);
+ historelptrespt->Fill(pt_truth,(pt_truth-TMath::Abs(bestftk->getPt())*1e-3)*invpt_truth);
+ bitmask=bestftk->getBitmask();
+ for(int icoord=0;icoord<16;icoord++) {
+ if (bitmask&(1<<icoord)){
+ histo2d0res_coordbit->Fill(icoord,d0_truth-bestftk->getIP());
+ }
+ }
}
}
+ // fill efficiency per plane in eta,phi
+ if(mrange.first!=mrange.second) {
+ bitmask |= 1<<16;
+ for(Int_t icoord=0;icoord<=16;++icoord) {
+ if (bitmask&(1<<icoord)){
+ histoeff_etaphiplane[icoord]->Fill(eta_truth,phi_truth);
+ }
+ }
+ }
+
} // end loop over truth tracks
// fill counter histograms
@@ -720,6 +758,15 @@ void Terminate(std::string& outputname) {
ofile->Add(histocoordmasketa_ftk);
ofile->Add(histocoordmaskphi_ftk);
+ for(int i=0;i<16;i++) {
+ histocoordmasketaphi_ftk[i]->Divide(histocoordmasketaphi_ftk[16]);
+ histoeff_etaphiplane[i]->Divide(histoeff_etaphiplane[16]);
+ }
+
+ for(int i=0;i<16;i++) {
+ ofile->Add(histocoordmasketaphi_ftk[i]);
+ ofile->Add(histoeff_etaphiplane[i]);
+ }
ofile->Add(histocoordmaskz0_ftk);
ofile->Add(histonmisseta_ftk);
ofile->Add(histochisqndfeta_ftk);
@@ -814,6 +861,12 @@ void Terminate(std::string& outputname) {
ofile->Add(histod0res_veta);
ofile->Add(histoz0res_veta);
ofile->Add(histod0res_vphi);
+ ofile->Add(histo2d0res_vphi);
+ ofile->Add(histo2d0res_veta);
+ ofile->Add(histo2d0res_coordbit);
+ ofile->Add(histo2d0truth_vphi);
+ ofile->Add(histo2d0ftk_vphi);
+ ofile->Add(histo2curvCurv);
ofile->Add(histoz0res_vphi);
ofile->Add(histod0res_vz0);
ofile->Add(histoz0res_vz0);
@@ -898,8 +951,11 @@ int main(int argc, char **argv) {
std::string output;
std::vector<std::string> files;
ptmincut = 0;
- dx = -0.5;
- dy = -0.5;
+ //dx = -0.5;
+ //dy = -0.5;
+ vtxTruth[0]=0.;
+ vtxTruth[1]=0.;
+ vtxTruth[2]=0.;
try {
std::string appName = boost::filesystem::basename(argv[0]);
@@ -915,8 +971,11 @@ int main(int argc, char **argv) {
("use-first-stage", po::value<int>(&Use1stStage)->default_value(0), "-1: Use roads, 1: Use 1st stage tracks, 0(default): Use 2nd stage tracks")
("files", po::value<std::vector<std::string> >(&files)->multitoken(), "FTK NTUP files")
("ptmincut", po::value<double>(&ptmincut), "min pt cut on truth tracks")
- ("dx", po::value<double>(&dx)->default_value(-0.5), "dx")
- ("dy", po::value<double>(&dy)->default_value(-0.5), "dx")
+ //("dx", po::value<double>(&dx)->default_value(-0.5), "vertex-x for d0,z0 shift")
+ //("dy", po::value<double>(&dy)->default_value(-0.5), "vertex-y for d0,z0 shift")
+ ("vxTruth", po::value<double>(vtxTruth+0)->default_value(0.), "vertex-x shift for truth tracks")
+ ("vyTruth", po::value<double>(vtxTruth+1)->default_value(0.), "vertex-y shift for truth tracks")
+ ("vzTruth", po::value<double>(vtxTruth+2)->default_value(0.), "vertex-z shift for truth tracks")
("psfile", "Produce postscript file with efficieny plots");
po::variables_map vm;
diff --git a/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.h b/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.h
index 8e5bea844dc421993628bc60c58e9f8a009da098..149bbecaf5481438ecdb772f80d90db36945e53b 100644
--- a/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.h
+++ b/Trigger/TrigFTK/TrigFTKSim/standalone/efficiency.h
@@ -44,10 +44,12 @@ double Dcurv;
double Deta;
double Dz0;
double ptmincut;
-double dx, dy;
+//double dx, dy;
+double vtxTruth[3];
// block of generic control histograms for the FTK tracks
TH2F *histocoordmasketa_ftk;
+TH2F *histocoordmasketaphi_ftk[17];
TH2F *histocoordmaskz0_ftk;
TH2F *histocoordmaskphi_ftk;
TH2F *histonmisseta_ftk;
@@ -69,6 +71,8 @@ TH1F *histocurv_ftk;
TH1F *histoeta_ftk;
TH2F *histoetaphi_ftk;
+TH2F *histoeff_etaphiplane[17];
+
TH2F *histoetaphi_ftk_IBL;
TH2F *histoetaphi_ftk_PixL0;
@@ -161,6 +165,12 @@ TH1F *histoz0res;
TProfile *histod0res_veta;
TProfile *histoz0res_veta;
TProfile *histod0res_vphi;
+TH2F *histo2d0res_vphi;
+TH2F *histo2d0res_veta;
+TH2F *histo2d0res_coordbit;
+TH2F *histo2d0truth_vphi;
+TH2F *histo2d0ftk_vphi;
+TH2F *histo2curvCurv;
TProfile *histoz0res_vphi;
TProfile *histod0res_vz0;
TProfile *histoz0res_vz0;
diff --git a/Trigger/TrigFTK/TrigFTKSim/standalone/make_conn.cc b/Trigger/TrigFTK/TrigFTKSim/standalone/make_conn.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c4d1ab49b3acf1ba724951a41602f4bb54cbf2cb
--- /dev/null
+++ b/Trigger/TrigFTK/TrigFTKSim/standalone/make_conn.cc
@@ -0,0 +1,232 @@
+#include <iostream>
+#include <string>
+#include <vector>
+#include <set>
+#include <map>
+
+#include <TString.h>
+
+#include <glob.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+
+#include <boost/iostreams/device/file.hpp>
+//#include <boost/iostreams/concepts.hpp>
+//#include <boost/iostreams/categories.hpp>
+//#include <boost/algorithm/string/classification.hpp>
+//#include <boost/filesystem.hpp>
+#include <boost/iostreams/filtering_stream.hpp>
+#include <boost/iostreams/filter/bzip2.hpp>
+
+using namespace std;
+
+class SSIDvector : public vector<int> {
+public:
+ SSIDvector(int n=0) : vector<int>(n) { }
+ int operator<(const SSIDvector &b) const {
+ if(size()<b.size()) return true;
+ if(size()>b.size()) return false;
+ for(size_t n=0;n<size();n++) {
+ if((*this)[n]<b[n]) return true;
+ if((*this)[n]>b[n]) return false;
+ }
+ return false;
+ }
+};
+
+typedef map<SSIDvector,pair<int,vector<vector<int> > > > ConnectionMap_t;
+
+int main(int argc, char const *argv[]) {
+ if(argc!=2) {
+ cout<<"usage: "<<argv[0]<<" path\n";
+ exit(0);
+ }
+ TString path(argv[1]);
+ //
+ // locate 8L sector files
+ glob_t globList;
+ struct stat statBuf;
+ memset(&globList,0,sizeof(glob_t));
+ TString sector8Lpath=path+"/sectors_raw_8L_reg*";
+ int err=glob(sector8Lpath,GLOB_TILDE,0,&globList);
+ if(err) {
+ cout<<"Error: can not glob "<<sector8Lpath<<"\n";
+ exit(1);
+ }
+ int nTowerMax=globList.gl_pathc;
+ for(int tower=0;tower<nTowerMax;tower++) {
+ sector8Lpath=path+TString::Format("/sectors_raw_8L_reg%d.patt*",tower);
+ memset(&globList,0,sizeof(glob_t));
+ err=glob(sector8Lpath,GLOB_TILDE,0,&globList);
+ if(err) {
+ cout<<"Can not glob "<<sector8Lpath<<", skipping\n";
+ continue;
+ }
+ TString sectorFile[2];
+ sectorFile[0]=globList.gl_pathv[0];
+ TString sector12Lpath=path+
+ TString::Format("/sectors_raw_12L_reg%d.patt*",tower);
+ memset(&globList,0,sizeof(glob_t));
+ err=glob(sector12Lpath,GLOB_TILDE,0,&globList);
+ if(err) {
+ cout<<"Can not glob "<<sector12Lpath<<", skipping\n";
+ continue;
+ }
+ sectorFile[1]=globList.gl_pathv[0];
+ int error=0;
+ boost::iostreams::filtering_istream *sectorStream[2];
+ for(int i=0;i<2;i++) {
+ if(stat(sectorFile[i],&statBuf)) {
+ cout<<"Can not stat "<<sectorFile[i]<<", skipping\n";
+ error++;
+ }
+ if(!S_ISREG(statBuf.st_mode)) {
+ cout<<"Not a regular file, skipping: "
+ <<sectorFile[i]<<"\n";
+ error++;
+ }
+ sectorStream[i]=new boost::iostreams::filtering_istream;
+ if(sectorFile[i].EndsWith(".bz2")) {
+ sectorStream[i]->push(boost::iostreams::bzip2_decompressor());
+ }
+ sectorStream[i]->push(boost::iostreams::file_source(sectorFile[i].Data()));
+ cout<<"open "<<sectorFile[i]<<"\n";
+ }
+ if(error) continue;
+ // read sector files
+
+ vector<set<int> > moduleId;
+ ConnectionMap_t connections;
+ int notFound=0;
+ int found=0;
+ for(int iType=0;iType<2;iType++) {
+ int nSector,nPlane;
+ (*sectorStream[iType])>>nSector>>nPlane;
+ if(iType==0) moduleId.resize(nPlane);
+ SSIDvector ssid(nPlane);
+ for(int iSector=0;iSector<nSector;iSector++) {
+ // read sector data
+ int sectorID;
+ (*sectorStream[iType])>>sectorID;
+ //cout<<sectorID;
+ if(sectorID!=iSector) {
+ cout<<"error reading sector "<<iSector
+ <<" "<<sectorID<<"\n";
+ exit(0);
+ }
+ // read module IDs
+ for(int plane=0;plane<nPlane;plane++) {
+ (*sectorStream[iType])>>ssid[plane];
+ //cout<<" "<<ssid[plane];
+ }
+ // dummy
+ int dummy,nTrack;
+ (*sectorStream[iType])>>dummy>>nTrack;
+ //cout<<" "<<dummy<<" "<<nTrack<<"\n";
+ //if(sectorID>20) exit(0);
+ SSIDvector ssid8L(moduleId.size());
+ if(iType==0) {
+ // if 8L: store known module IDs and store 8L SSID
+ for(int plane=0;plane<nPlane;plane++) {
+ moduleId[plane].insert(ssid[plane]);
+ ssid8L[plane]=ssid[plane];
+ }
+ connections[ssid8L].first=iSector;
+ } else {
+ // if 12L: locate known planes, skip others
+ vector<int> connId(nPlane-moduleId.size()+2);
+ connId[1]=iSector;
+ size_t nExtra=2;
+ size_t plane8=0;
+ for(int plane=0;plane<nPlane;plane++) {
+ if((plane8<moduleId.size())&&
+ (moduleId[plane8].find(ssid[plane])!=
+ moduleId[plane8].end())) {
+ //cout<<" "<<ssid[plane]<<" "<<plane<<" -> "<<plane8<<"\n";
+ ssid8L[plane8++]=ssid[plane];
+ } else {
+ if(nExtra>=connId.size()) {
+ cout<<"12L sector="<<iSector<<" too few planes match\n";
+ exit(0);
+ }
+ //cout<<" "<<ssid[plane]<<" extra="<<nExtra<<"\n";
+ connId[nExtra++]=ssid[plane];
+ }
+ }
+ //exit(0);
+ ConnectionMap_t::iterator ic=connections.find(ssid8L);
+
+ /* cout<<" 8L:";
+ for(int k=0;k<ssid8L.size();k++) {
+ cout<<" "<<ssid8L[k];
+ }
+ cout<<" extra ";
+ for(int k=0;k<connId.size();k++) {
+ cout<<" "<<connId[k];
+ }
+ cout<<"\n";
+ */
+ if(ic==connections.end()) {
+ notFound++;
+ //cout<<"Number of sectors: "<<connections.size()<<"\n";
+ //cout<<"12L sector="<<iSector<<" not found, skipping\n";
+ /* for(ConnectionMap_t::const_iterator ic=connections.begin();
+ ic!=connections.end();ic++) {
+ if((*ic).second.first<20) {
+ cout<<(*ic).second.first;
+ for(int k=0;k<(*ic).first.size();k++) {
+ cout<<" "<<(*ic).first[k];
+ }
+ cout<<"\n";
+ }
+ }
+ exit(0); */
+ } else {
+ found++;
+ vector<vector<int> > &conn12=(*ic).second.second;
+ conn12.push_back(connId);
+ }
+ }
+ }
+ if(iType==0) {
+ cout<<"nSector="<<nSector<<" nPlane="<<nPlane
+ <<" Number of sectors: "<<connections.size()<<"\n";
+ for(size_t plane=0;plane< moduleId.size();plane++) {
+ cout<<" 8L["<<plane<<"] n(mod)="<<moduleId[plane].size();
+ }
+ cout<<"\n";
+ } else {
+ cout<<"nSector="<<nSector<<" nPlane="<<nPlane
+ <<" notFound="<<notFound<<" found="<<found<<"\n";
+ // order connections by sector id, then print
+ map<int,ConnectionMap_t::const_iterator> orderedConnections;
+ for(ConnectionMap_t::const_iterator ic=connections.begin();
+ ic!=connections.end();ic++) {
+ orderedConnections[(*ic).second.first]=ic;
+ }
+ TString outname(TString::Format("sectors_raw_8L_reg%d.conn",tower));
+ ofstream output(outname);
+ int sectorID12=0;
+ for(map<int,ConnectionMap_t::const_iterator>::const_iterator
+ oic=orderedConnections.begin();
+ oic!=orderedConnections.end();oic++) {
+ vector<vector<int> > const &conn12=(*(*oic).second).second.second;
+ output<<(*oic).first<<" "<<conn12.size();
+ for(size_t k=0;k<conn12.size();k++) {
+ output<<" 0 "<<sectorID12;
+ sectorID12++;
+ for(size_t l=2;l<conn12[k].size();l++) {
+ output<<" "<<conn12[k][l];
+ }
+ }
+ output<<"\n";
+ }
+ }
+ }
+
+ for(int i=0;i<2;i++) {
+ delete sectorStream[i];
+ }
+ }
+}
diff --git a/Trigger/TrigFTK/TrigFTKSim/standalone/partitionbalancing.cc b/Trigger/TrigFTK/TrigFTKSim/standalone/partitionbalancing.cc
index 857bea9514a9bee8e21f30cd20bd02659bbde183..9297c2623b073bfea2b598b1ed9c0951f14a5462 100644
--- a/Trigger/TrigFTK/TrigFTKSim/standalone/partitionbalancing.cc
+++ b/Trigger/TrigFTK/TrigFTKSim/standalone/partitionbalancing.cc
@@ -52,6 +52,7 @@ public:
fSteering->AddIntPar("NPATTERN",1,0);
fSteering->AddIntPar("maxSectorId",1,0);
+ fSteering->AddIntPar("rebinSlices",1,0);
fSteering->AddDoublePar("absetarange",2,0.0);
@@ -66,6 +67,9 @@ public:
fSteering->AddStringPar("pattern_bank_path");
fSteering->AddStringPar("slices_file_path");
+ fSteering->AddStringPar("maximum_input_file");
+ fSteering->AddStringPar("maximum_output_file");
+
}
return fSteering;
}
@@ -81,6 +85,12 @@ public:
const char *GetSlicesFilePath(void) const {
return *(*this)["slices_file_path"];
}
+ const char *GetMaximumInputFile(void) const {
+ return *(*this)["maximum_input_file"];
+ }
+ const char *GetMaximumOutputFile(void) const {
+ return *(*this)["maximum_output_file"];
+ }
double GetAbsEtaRangePar(int i) const {
return (*(*this)["absetarange"])[i];
}
@@ -96,6 +106,9 @@ public:
int GetMaxSectorId(void) const {
return (*this)["maxSectorId"][0];
}
+ int GetRebinSlices(void) const {
+ return (*this)["rebinSlices"][0];
+ }
int GetNPattern(void) const {
return (*this)["NPATTERN"][0];
}
@@ -173,6 +186,9 @@ int main(int argc, char const *argv[]) {
// this is a limit on the maximum allowed sector ID
// -> sectorIDs >= max will not be written to the partition file
int maxSectorId=PartitionSteering::Instance()->GetMaxSectorId();
+ //
+ // this combines several half-slices into one bin
+ int rebinSlices=PartitionSteering::Instance()->GetRebinSlices();
cout<<"abs(eta) range: "<<absetamin<<"..."<<absetamax<<"\n";
cout<<"etabinOffset: "<<etabinOffset<<"\n";
@@ -189,8 +205,46 @@ int main(int argc, char const *argv[]) {
TH2D *originalPerSlice2(0),*partitionedPerSlice2(0);
TH1D *etaPerSlice2(0),*cosThetaPerSlice2(0);
+ vector<map<int,int> > knownMax(nreg);
+ char const *maximumInputFile=PartitionSteering::Instance()->GetMaximumInputFile();
+ if(maximumInputFile) {
+ ifstream in(maximumInputFile);
+ int line=1;
+ while(!in.eof()) {
+ int ireg,sector,maximum;
+ in>>ireg>>sector>>maximum;
+ if(in.fail()) break;
+ if((ireg<0)||(ireg>=(int)knownMax.size())) {
+ logging.Error("maximumFromFile")<<"bad region="<<ireg<<" at line="<<line<<"\n";
+ break;
+ }
+ if(knownMax[ireg].find(sector)!=knownMax[ireg].end()) {
+ logging.Error("maximumFromFile")<<"duplicate sector="<<sector<<" in region="<<ireg<<" at line="<<line<<"\n";
+ break;
+ }
+ if(maximum<=0) {
+ logging.Error("maximumFromFile")<<"max maximum="<<maximum<<" in sector="<<sector<<" region="<<ireg<<" at line="<<line<<"\n";
+ break;
+ }
+ knownMax[ireg][sector]=maximum;
+ line++;
+ }
+ }
+
+ ofstream *dumpMaximum=0;
+ char const *maximumOutputFile= PartitionSteering::Instance()->GetMaximumOutputFile();
+ if(maximumOutputFile) {
+ dumpMaximum=new ofstream(maximumOutputFile);
+ }
+
for(int reg=0;reg<nreg;reg++) {
- logging.Info("loop")<<"Processing region "<<reg<<"\n";
+
+ struct SectorPatternData {
+ int fNPattern;
+ int fNMax;
+ };
+
+ logging.Info("loop")<<"Processing region "<<reg<<" known max for "<<knownMax[reg].size()<<" sector(s)\n";
// read slices
TString sliceFileName=
PartitionSteering::Instance()->GetSlicesFilePath()+
@@ -201,7 +255,7 @@ int main(int argc, char const *argv[]) {
// ctheta contains for each slice in cos(theta) a bitvector
// the bits correspond to sector IDs
// if the bit is on, the sector may have patterns
- // in teh given cos(theta) bin
+ // in the given cos(theta) bin
TClonesArray *ctheta;
TString name("c_bits_ctheta");
slices->GetObject(name,ctheta);
@@ -218,12 +272,16 @@ int main(int argc, char const *argv[]) {
// either read one subregion or read nSub subregions
// this variable will contain the number of patterns per sector
- map<int,int> patternsPerSector;
+ // first : patterns per sector
+ // second : 1 if at limit, -1 or 0 otherwise
+ map<int,SectorPatternData> patternsPerSector;
+
// this variable will contain the total number of patterns
int nPatternRead=0;
for(int isubRead=0;isubRead<nsub;isubRead++) {
// read pattern bank(s) and determine number of patterns
// per sector
+ int nPatternSub=0;
TString bankName=
PartitionSteering::Instance()->GetPatternBankPath()+
TString::Format("/*_reg%d_sub%d.p*.root",reg,isubRead);
@@ -266,46 +324,104 @@ int main(int argc, char const *argv[]) {
//
// count sector multiplicities
// for the given bank, read the sectorID of all patterns
- // and cout the (original) number of patterns per sector ID
+ // and count the (original) number of patterns per sector ID
TTree *tree;
- patternBank->GetObject("Bank",tree);
- if(!tree) {
- logging.Fatal("subregion")<<"can not find TTree \"Bank\"\n";
- }
- int nPattern=tree->GetEntriesFast();
- FTKPattern *pattern=new FTKPattern();
- if(tree->SetBranchAddress("Pattern",&pattern)!=0) {
- // no match, assume this is a flat file
- logging.Info("subregion")<<"try to read flat format\n";
- delete pattern;
- pattern=0;
- }
- TBranch *branch_sectorID=0;
- int sectorID;
- if(!pattern) {
- // flat format
- // access branches and set addresses
- branch_sectorID=tree->GetBranch("sectorID");
- branch_sectorID->SetAddress(§orID);
+ patternBank->GetObject("Sector",tree);
+ if(tree) {
+ int sector,first,last,allPatternsUsed;
+ tree->SetBranchAddress("sector",§or);
+ tree->SetBranchAddress("first",&first);
+ tree->SetBranchAddress("last",&last);
+ tree->SetBranchAddress("allPatternsUsed",&allPatternsUsed);
+ for(int i=0;i<tree->GetEntries();i++) {
+ allPatternsUsed=0;
+ tree->GetEntry(i);
+ int n=last+1-first;
+ for(int special=((first-1) & 0xfffe0000)+0x20000;special<=last;special += 0x20000) {
+ n--; // special patterns are not counted
+ }
+ SectorPatternData &pps=patternsPerSector[sector];
+ pps.fNPattern = n;
+ pps.fNMax = (allPatternsUsed==1) ? n : PartitionSteering::Instance()->GetNPattern();
+ nPatternSub+=n;
+ }
} else {
- branch_sectorID=tree->GetBranch("m_sectorID");
- }
- for(int iPattern=0;iPattern<nPattern;++iPattern) {
- branch_sectorID->GetEntry(iPattern);
- if(pattern) {
- sectorID=pattern->getSectorID();
+ patternBank->GetObject("Bank",tree);
+ if(!tree) {
+ logging.Fatal("subregion")<<"can not find TTree \"Bank\"\n";
+ }
+ int nPattern=tree->GetEntriesFast();
+ FTKPattern *pattern=new FTKPattern();
+ if(tree->SetBranchAddress("Pattern",&pattern)!=0) {
+ // no match, assume this is a flat file
+ logging.Info("subregion")<<"try to read flat format\n";
+ delete pattern;
+ pattern=0;
+ }
+ TBranch *branch_sectorID=0;
+ int sectorID;
+ if(!pattern) {
+ // flat format
+ // access branches and set addresses
+ branch_sectorID=tree->GetBranch("sectorID");
+ branch_sectorID->SetAddress(§orID);
+ } else {
+ branch_sectorID=tree->GetBranch("m_sectorID");
+ }
+ for(int iPattern=0;iPattern<nPattern;++iPattern) {
+ branch_sectorID->GetEntry(iPattern);
+ if(pattern) {
+ sectorID=pattern->getSectorID();
+ }
+ patternsPerSector[sectorID].fNPattern++;
+ patternsPerSector[sectorID].fNMax=PartitionSteering::Instance()->GetNPattern();
+ nPatternSub++;
}
- patternsPerSector[sectorID]++;
- nPatternRead++;
}
globfree(&globStruct);
logging.Info("subregion")
- <<"read subregion="<<isubRead<<" nPattern="<<nPattern<<"\n";
+ <<"read subregion="<<isubRead<<" nPattern="<<nPatternSub<<"\n";
+ nPatternRead += nPatternSub;
}
logging.Info("readpattern")
<<"total sectors (bank): "<<patternsPerSector.size()
<<" total patterns: "<<nPatternRead<<"\n";
+ // add information abount known maxima
+ for(map<int,int>::const_iterator i=knownMax[reg].begin();
+ i!=knownMax[reg].end();i++) {
+ map<int,SectorPatternData>::iterator
+ i2=patternsPerSector.find((*i).first);
+ if(i2==patternsPerSector.end()) {
+ SectorPatternData &pps=patternsPerSector[(*i).first];
+ pps.fNPattern=0;
+ pps.fNMax=(*i).second;
+ } else if((*i2).second.fNMax==
+ PartitionSteering::Instance()->GetNPattern()) {
+ logging.Info("setMaximum")
+ <<"reg="<<reg<<" sector="<<(*i2).first
+ <<" fMax "<<(*i2).second.fNMax
+ <<" -> "<<(*i).second<<"\n";
+ (*i2).second.fNMax=(*i).second;
+ } else if((*i2).second.fNMax!=(*i).second) {
+ logging.Error("setMaximum")
+ <<"reg="<<reg<<" sector="<<(*i2).first
+ <<" fMax="<<(*i2).second.fNMax
+ <<" fMax(iteration-1)="<<(*i).second<<"\n";
+ }
+ }
+
+ // dump information about maximum per sector
+ if(dumpMaximum) {
+ for(map<int,SectorPatternData>::const_iterator i=patternsPerSector.begin();
+ i!=patternsPerSector.end();i++) {
+ if((*i).second.fNMax<
+ PartitionSteering::Instance()->GetNPattern()) {
+ (*dumpMaximum)<<reg<<" "<<(*i).first
+ <<" "<<(*i).second.fNMax<<"\n";
+ }
+ }
+ }
// here we have:
// ctheta[] : sectors per cos(theta) slice
@@ -353,49 +469,73 @@ int main(int argc, char const *argv[]) {
// This sets the number of partitions along eta
// there is one partition for each possible half-integer number
// empty partitions will be skipped.
- int nSlice2=2*nSlice1;
+
+ int nSliceTimes2=2*nSlice1;
+ int nSliceRebin=nSliceTimes2/(rebinSlices>0 ? rebinSlices : 1);
+ if((rebinSlices>0)&&(nSliceTimes2 % rebinSlices)) {
+ nSliceRebin++;
+ }
if(!originalPerSlice2) {
outputRoot->cd();
originalPerSlice2=new TH2D("originalPerSlice2",
";slice;tower",
- nSlice2,-0.5,nSlice2-0.5,
+ nSliceRebin,-0.5,nSliceRebin-0.5,
nreg,-0.5,nreg-0.5);
partitionedPerSlice2=new TH2D("partitionedPerSlice2",
";slice;tower",
- nSlice2,-0.5,nSlice2-0.5,
+ nSliceRebin,-0.5,nSliceRebin-0.5,
nreg,-0.5,nreg-0.5);
etaPerSlice2=new TH1D("etaPerSlice2",";slice;eta",
- nSlice2,-0.5,nSlice2-0.5);
+ nSliceRebin,-0.5,nSliceRebin-0.5);
cosThetaPerSlice2=new TH1D("cosThetaPerSlice2",
";slice;cos(theta)",
- nSlice2,-0.5,nSlice2-0.5);
+ nSliceRebin,-0.5,nSliceRebin-0.5);
}
// next, each sector is assigned to one partition "nslice2"
- vector<int> patternPerSlice(nSlice2);
- vector<vector<int> > sectorList(nSlice2);
+ // first : number of patterns
+ // second : if zero, all sectors are at limit
+ vector<SectorPatternData > patternPerSlice(nSliceRebin);
+ vector<vector<int> > sectorList(nSliceRebin);
for(map<int,double>::const_iterator is=sum0.begin();
is!=sum0.end();is++) {
int sector=(*is).first;
- int slice2=(int)(2.*sum1[sector]/(*is).second);
+ int sliceRebin=
+ ((int)(2.*sum1[sector]/(*is).second))/
+ (rebinSlices>0 ? rebinSlices : 0);
+ if(sliceRebin>=nSliceRebin) sliceRebin=nSliceRebin-1;
//sliceNumber[sector]=slice;
- sectorList[slice2].push_back(sector);
+ sectorList[sliceRebin].push_back(sector);
if(patternsPerSector.find(sector)!=patternsPerSector.end()) {
- patternPerSlice[slice2]+=patternsPerSector[sector];
+ SectorPatternData &ppSlice=patternPerSlice[sliceRebin];
+ SectorPatternData &ppSector=patternsPerSector[sector];
+ ppSlice.fNPattern += ppSector.fNPattern;
+ ppSlice.fNMax += ppSector.fNMax;
+ if(ppSlice.fNMax>PartitionSteering::Instance()->GetNPattern()) {
+ ppSlice.fNMax=PartitionSteering::Instance()->GetNPattern();
+ }
+ }
+ }
+ for(int i=0;i<nSliceRebin;i++) {
+ if((patternPerSlice[i].fNMax>0)&&(patternPerSlice[i].fNMax<PartitionSteering::Instance()->GetNPattern())) {
+ logging.Info("Truncation")<<"isub="<<isub<<" islice="<<i<<" nmax="<<patternPerSlice[i].fNMax
+ <<" npattern="<<patternPerSlice[i].fNPattern<<"\n";
}
}
// here:
// sectorList[partition][] : the sectors in this partition
- // patternPerSlice[partition] : number of patterns in partition
+ // patternPerSlice[partition]
+ // : number of patterns and maximum in partition
// (2) for each slice determine the efficiency and add up
// the target number of patterns
- vector<double> targetPatterns(nSlice2);
+ vector<double> targetPatterns(nSliceRebin);
double ctheta0=TMath::SinH(-3.),ctheta1=TMath::SinH(3.);
double total=0.;
- for(int i=0;i<nSlice2;i++) {
+ for(int i2=0;i2<nSliceTimes2;i2++) {
// get cos(theta) from slice number
- double ctheta=(i+etabinOffset)*(ctheta1-ctheta0)/nSlice2+ctheta0;
+ double ctheta=(i2+etabinOffset)*(ctheta1-ctheta0)/
+ nSliceTimes2+ctheta0;
// get eta coordinate (for efficiency plot)
double eta=TMath::ASinH(ctheta);
double abseta=TMath::Abs(eta);
@@ -403,11 +543,13 @@ int main(int argc, char const *argv[]) {
double epsilon=efficiencyPlot->Eval(eta);
double epsilon0=
PartitionSteering::Instance()->GetEpsilonLimit();
+ int i=i2/(rebinSlices>0 ? rebinSlices : 1);
+ if(i>=nSliceRebin) i=nSliceRebin-1;
cosThetaPerSlice2->SetBinContent(i+1,ctheta);
etaPerSlice2->SetBinContent(i+1,eta);
if(((absetamin>=absetamax)||
((abseta>=absetamin)&&(abseta<=absetamax)))&&
- (patternPerSlice[i]>0)) {
+ (patternPerSlice[i].fNPattern>0)) {
double weight=PartitionSteering::Instance()
->GetWeightLimit();
if(epsilon<epsilon0) {
@@ -417,26 +559,23 @@ int main(int argc, char const *argv[]) {
weight=PartitionSteering::Instance()->GetWeightLimit();
}
}
- // weight: desired weighting factor fro which the predicted
+ // weight: desired weighting factor for which the predicted
// efficiency reaches the target efficiency
- targetPatterns[i]+= patternPerSlice[i]*weight;
+ targetPatterns[i]+= patternPerSlice[i].fNPattern*weight;
total+= targetPatterns[i];
}
- originalPerSlice2->SetBinContent(i+1,reg+1,patternPerSlice[i]);
+ originalPerSlice2->SetBinContent
+ (i+1,reg+1,patternPerSlice[i].fNPattern);
}
// here:
// targetPatterns[partition] : patterns needed to reach target
// efficiency
// total : total number of patterns needed to reach target
- // (3) normalize and round
+ // (3) normalize
// targetPatterns[] should be scaled down by the factor
// nPattern/total
- // Problem: finally, all number have to be integer
- // fractional patterns have to be rounded up or down
- //
- int tot=0;
// target number of patterns in this subregion
int nPattern=nPatternRead/nsub;
if(PartitionSteering::Instance()->GetNPattern())
@@ -444,53 +583,105 @@ int main(int argc, char const *argv[]) {
if(total==0) {
logging.Warning("subregion")
<<"no sectors accepted, store one pattern per partition\n";
- nPattern=nSlice2;
+ nPattern=nSliceTimes2;
+ }
+ for(int i=0;i<nSliceRebin;i++) {
+ targetPatterns[i]*=((double)nPattern)/total;
}
- for(int i=0;i<nSlice2;i++) {
+ // (4) take into account truncation
+ double truncate=0.;
+ do {
+ total=0.;
+ double free=0.;
+ truncate=0.;
+ for(int i=0;i<nSliceRebin;i++) {
+ if(targetPatterns[i]>=patternPerSlice[i].fNMax) {
+ truncate += targetPatterns[i]-patternPerSlice[i].fNMax;
+ if(targetPatterns[i]>0.0) {
+ logging.Info("target")<<" slice="<<i<<" change from "<<targetPatterns[i]<<" to "<<patternPerSlice[i].fNMax<<"\n";
+ }
+ targetPatterns[i]=patternPerSlice[i].fNMax;
+ } else {
+ free+=targetPatterns[i];
+ }
+ total+=targetPatterns[i];
+ }
+ logging.Info("target")<<" total="<<total<<" nPattern="<<nPattern<<" free="<<free<<" truncate="<<truncate<<"\n";
+ if(free<=0.) break;
+ // truncate: excess patterns to be distributed
+ double scale=truncate/free;
+ for(int i=0;i<nSliceRebin;i++) {
+ if(targetPatterns[i]<patternPerSlice[i].fNMax-0.5) {
+ double n=targetPatterns[i]*scale;
+ targetPatterns[i] += n;
+ truncate -= n;
+ }
+ }
+ } while(truncate>1.0);
+
+ // (5) finally, all numbers have to be integer
+ // fractional patterns have to be rounded up or down
+ // also, take into account truncation
+
+ // round to integer
+ int tot=0;
+ for(int i=0;i<nSliceRebin;i++) {
// skip empty partitions
- if(patternPerSlice[i]==0) continue;
+ if(patternPerSlice[i].fNPattern==0) continue;
if(total>0) {
- // normalized targen number of patterns
- patternPerSlice[i]=
- (int)((targetPatterns[i]*nPattern)/total+0.5);
+ int target=(int)(targetPatterns[i]+0.5);
+ if(patternPerSlice[i].fNMax>target) {
+ // set new normalized target number of patterns
+ // (only for those partitions which are not at limit)
+ patternPerSlice[i].fNPattern=target;
+ } else {
+ patternPerSlice[i].fNPattern=patternPerSlice[i].fNMax;
+ }
+ tot+= patternPerSlice[i].fNPattern;
} else {
- patternPerSlice[i]=1;
+ patternPerSlice[i].fNPattern=1;
}
- tot+= patternPerSlice[i];
}
- // here:
- // patternPerSlice[i] : renormalized target number of patterns
- // integer
- // !!! the original content is overwritten !!!
- // tot : sum of renormalizerd target number of patterns
-
- // execute tis loop to adjust the number of patterns
+ logging.Info("adjust")<<"isub="<<isub<<" tot="<<tot<<" nPattern="<<nPattern<<"\n";
+ // execute this loop to adjust the number of patterns
// until the integer number: tot
// equals the target number: nPattern
+ int nAdjustTotal=0;
+ int nLimitTotal=0;
while(tot-nPattern) {
// try all partitions in sequence
- for(int i=0;i<nSlice2;i++) {
+ int nAdjust=0;
+ int nLimit=0;
+ for(int i=0;i<nSliceRebin;i++) {
// do not adjust empty partitions
- if(patternPerSlice[i]<=0) continue;
+ if(patternPerSlice[i].fNPattern<=0) continue;
+ // do not adjust partitions which are at limit
int delta=tot-nPattern;
if(!delta) break; // done
else if(delta>0) {
// too many patterns, subtract one in this partition
tot--;
- patternPerSlice[i]--;
+ patternPerSlice[i].fNPattern--;
+ nAdjust++;
} else if(delta<0) {
// too few patterns, add one in this partition
- // (the if is redundant???)
- if(patternPerSlice[i]>0) {
+ if(patternPerSlice[i].fNPattern<patternPerSlice[i].fNMax) {
tot++;
- patternPerSlice[i]++;
+ patternPerSlice[i].fNPattern++;
+ nAdjust++;
+ } else {
+ nLimit++;
}
}
}
+ nAdjustTotal+=nAdjust;
+ nLimitTotal+=nLimit;
+ if(nAdjust==0) break;
}
- for(int i=0;i<nSlice2;i++) {
- partitionedPerSlice2->SetBinContent(i+1,reg+1,
- patternPerSlice[i]);
+ logging.Info("adjust")<<"total="<<tot<<" nPattern="<<nPattern<<" adjustTotal="<<nAdjustTotal<<" limit="<<nLimitTotal<<"\n";
+ for(int i=0;i<nSliceRebin;i++) {
+ partitionedPerSlice2->SetBinContent
+ (i+1,reg+1,patternPerSlice[i].fNPattern);
}
// done
// write out partition file
@@ -498,10 +689,11 @@ int main(int argc, char const *argv[]) {
int nSectorUsed=0;
int nSectorUsedBank=0;
int nSliceUsed=0;
- for(int i=0;i<nSlice2;i++) {
- if(sectorList[i].size() && patternPerSlice[i]) {
- sum += patternPerSlice[i];
- (*outputFile)<<"PARTITION "<<iPart<<" "<<patternPerSlice[i]
+ for(int i=0;i<nSliceRebin;i++) {
+ if(sectorList[i].size() && patternPerSlice[i].fNPattern) {
+ sum += patternPerSlice[i].fNPattern;
+ (*outputFile)<<"PARTITION "<<iPart<<" "
+ <<patternPerSlice[i].fNPattern
<<" "<<sectorList[i].size();
for(unsigned j=0;j<sectorList[i].size();j++) {
if(!j) (*outputFile)<<"\n";
@@ -516,6 +708,11 @@ int main(int argc, char const *argv[]) {
}
nSectorUsed+=sectorList[i].size();
(*outputFile)<<"\n";
+ logging.Info("result")
+ <<"ipart="<<iPart<<" isub="<<isub <<" islice="<<i
+ <<" npattern="<<patternPerSlice[i].fNPattern
+ <<" nmax="<<patternPerSlice[i].fNMax
+ <<"\n";
iPart++;
nSliceUsed++;
}
@@ -557,6 +754,9 @@ int main(int argc, char const *argv[]) {
etaPerSlice2->Write();
cosThetaPerSlice2->Write();
delete outputRoot;
+
+ if(dumpMaximum) delete dumpMaximum;
+
logging.Info("output")
<<"maxSector="<<totalMaxSector<<" maxSectorUsed="<<totalMaxSectorUsed
<<" nSectorTotal="<<totalNsectorTotal<<" nSectorUsed="<<totalNsectorUsed