diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp index 730f61e34df936eab5d54bb4e84fff1ca6866a05..9c073a70b8a4723435abf9abc0f80e53399641a5 100644 --- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp +++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp @@ -98,10 +98,10 @@ StatusCode PlanarDigiAlg::initialize() error() << "Failed to find TrackSystemSvc ..." << endmsg; return StatusCode::FAILURE; } - + MarlinTrk::IMarlinTrkSystem* _trksystem = _trackSystemSvc->getTrackSystem(this); _trksystem->init(); - + _trackSystemSvc->removeTrackSystem(this); return GaudiAlgorithm::initialize(); diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp index cf0d6ec9b44a4014462ec2ecf9f74499d3bdb90d..4d79865881b5224c6d6afab12f607a669f53b082 100644 --- a/Service/GearSvc/src/GearSvc.cpp +++ b/Service/GearSvc/src/GearSvc.cpp @@ -87,7 +87,7 @@ StatusCode GearSvc::initialize() const std::map<std::string, dd4hep::DetElement>& subs = world.children(); for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ dd4hep::DetElement sub = it->second; - info() << it->first << " " << sub.path() << " " << sub.placementPath() << endmsg; + debug() << it->first << " " << sub.path() << " " << sub.placementPath() << endmsg; if(it->first=="Tube"||it->first=="BeamPipe"){ sc = convertBeamPipe(sub); } @@ -109,38 +109,43 @@ StatusCode GearSvc::initialize() else if(it->first=="SET"){ sc = convertSET(sub); } - else{ - info() << it->first << " will convert in future! now fake parameters" << endmsg; - } - if(sc==StatusCode::FAILURE) return sc; - } - - gear::CalorimeterParametersImpl* barrelParam = new gear::CalorimeterParametersImpl(1847.415655, 2350., 8, 0.); - gear::CalorimeterParametersImpl* endcapParam = new gear::CalorimeterParametersImpl(400., 2088.8, 2450., 2, 0.); - for(int i=0;i<29;i++){ - if(i<19){ - barrelParam->layerLayout().positionLayer(0, 5.25, 1.016666667e+01, 1.016666667e+01, 2.1); - endcapParam->layerLayout().positionLayer(0, 5.25, 1.016666667e+01, 1.016666667e+01, 2.1); - } - else if(i<20){ - barrelParam->layerLayout().positionLayer(0, 6.3, 1.016666667e+01, 1.016666667e+01, 2.1); - endcapParam->layerLayout().positionLayer(0, 6.3, 1.016666667e+01, 1.016666667e+01, 2.1); + else if(it->first=="EcalBarrel"||it->first=="EcalEndcap"||it->first=="EcalPlug"|| + it->first=="HcalBarrel"||it->first=="HcalEndcap"||it->first=="HcalRing"|| + it->first=="YokeBarrel"||it->first=="YokeEndcap"||it->first=="YokePlug"|| + it->first=="Coil"){ + sc = convertCal(sub); } else{ - barrelParam->layerLayout().positionLayer(0, 4.2, 1.016666667e+01, 1.016666667e+01, 4.2); - endcapParam->layerLayout().positionLayer(0, 4.2, 1.016666667e+01, 1.016666667e+01, 4.2); + info() << it->first << " will convert in future!" << endmsg; } + if (sc==StatusCode::FAILURE) return sc; } - m_gearMgr->setEcalBarrelParameters(barrelParam); - m_gearMgr->setEcalEndcapParameters(endcapParam); - gear::CalorimeterParametersImpl* barrelYokeParam = new gear::CalorimeterParametersImpl(4173.929932, 4072., 12, 0.0); - gear::CalorimeterParametersImpl* endcapYokeParam = new gear::CalorimeterParametersImpl(320., 7414.929932, 4072., 2, 0.0); - gear::CalorimeterParametersImpl* plugYokeParam = new gear::CalorimeterParametersImpl(320., 2849.254326, 3781.43, 2, 0.0); - plugYokeParam->setDoubleVal("YokePlugThickness", 290.57) ; - m_gearMgr->setYokeBarrelParameters(barrelYokeParam) ; - m_gearMgr->setYokeEndcapParameters(endcapYokeParam) ; - m_gearMgr->setYokePlugParameters(plugYokeParam) ; + try { + const auto& ecalB = m_gearMgr->getEcalBarrelParameters(); + } + catch (...) { + info() << "EcalBarrelParameters not create! create fake parameters (big size) now" << endmsg; + gear::CalorimeterParametersImpl* barrelParam = new gear::CalorimeterParametersImpl(2500, 4500, 8, 0.); + barrelParam->layerLayout().positionLayer(0, 5.25, 10, 10, 2.1); + m_gearMgr->setEcalBarrelParameters(barrelParam); + } + try { + const auto& ecalE = m_gearMgr->getEcalEndcapParameters(); + } + catch (...) { + info() << "EcalEndcapsParameters not create! create fake parameters (big size) now" << endmsg; + gear::CalorimeterParametersImpl* endcapParam = new gear::CalorimeterParametersImpl(400., 3000, 4510, 2, 0.); + endcapParam->layerLayout().positionLayer(0, 5.25, 10, 10, 2.1); + m_gearMgr->setEcalEndcapParameters(endcapParam); + } + //gear::CalorimeterParametersImpl* barrelYokeParam = new gear::CalorimeterParametersImpl(4173.929932, 4072., 12, 0.0); + //gear::CalorimeterParametersImpl* endcapYokeParam = new gear::CalorimeterParametersImpl(320., 7414.929932, 4072., 2, 0.0); + //gear::CalorimeterParametersImpl* plugYokeParam = new gear::CalorimeterParametersImpl(320., 2849.254326, 3781.43, 2, 0.0); + //plugYokeParam->setDoubleVal("YokePlugThickness", 290.57) ; + //m_gearMgr->setYokeBarrelParameters(barrelYokeParam) ; + //m_gearMgr->setYokeEndcapParameters(endcapYokeParam) ; + //m_gearMgr->setYokePlugParameters(plugYokeParam) ; gear::GearXML::createXMLFile(m_gearMgr, "test.xml"); } @@ -995,7 +1000,78 @@ StatusCode GearSvc::convertSET(dd4hep::DetElement& set){ return StatusCode::SUCCESS; } -TGeoNode* GearSvc::FindNode(TGeoNode* mother, char* name){ +StatusCode GearSvc::convertCal(dd4hep::DetElement& cal) { + std::string name = cal.name(); + + dd4hep::rec::LayeredCalorimeterData* calData = nullptr; + bool extensionDataValid = true; + try{ + calData = cal.extension<dd4hep::rec::LayeredCalorimeterData>(); + } + catch(std::runtime_error& e){ + extensionDataValid = false; + info() << e.what() << " " << calData << endmsg; + } + if(calData){ + double rmin = calData->extent[0]/dd4hep::mm; + double rmax = calData->extent[1]/dd4hep::mm; + double zmin = calData->extent[2]/dd4hep::mm; + double zmax = calData->extent[3]/dd4hep::mm; + int inner_symmetry = calData->inner_symmetry; + int outer_symmetry = calData->outer_symmetry; + double phi0 = calData->phi0; + gear::CalorimeterParametersImpl* param = nullptr; + if (calData->layoutType==dd4hep::rec::LayeredCalorimeterData::BarrelLayout) { + param = new gear::CalorimeterParametersImpl(rmin, zmax, inner_symmetry, phi0); + } + else if (calData->layoutType==dd4hep::rec::LayeredCalorimeterData::EndcapLayout) { + param = new gear::CalorimeterParametersImpl(rmin, rmax, zmin, inner_symmetry, phi0); + } + else { + error() << "Not BarrelLayout and EndcapLayout: " << calData->layoutType << endmsg; + } + auto layers = calData->layers; + for (auto layer : layers) { + double distance = layer.distance/dd4hep::mm; + double thickness = layer.sensitive_thickness/dd4hep::mm; + double absorberThickness = layer.absorberThickness/dd4hep::mm; + double cellsize0 = layer.cellSize0/dd4hep::mm; + double cellsize1 = layer.cellSize1/dd4hep::mm; + param->layerLayout().positionLayer(0, thickness, cellsize0, cellsize1, absorberThickness); + } + if (calData->layoutType==dd4hep::rec::LayeredCalorimeterData::BarrelLayout) { + if (name.find("Ecal")!=-1) m_gearMgr->setEcalBarrelParameters(param); + else if (name.find("Hcal")!=-1) m_gearMgr->setHcalBarrelParameters(param); + else if (name.find("Yoke")!=-1) m_gearMgr->setYokeBarrelParameters(param); + else m_gearMgr->setGearParameters("CoilParameters", param); + + info() << "BarrelParameters set " << name << endmsg; + } + else if (calData->layoutType==dd4hep::rec::LayeredCalorimeterData::EndcapLayout) { + if (name.find("Endcap")) { + if (name.find("Ecal")!=-1) m_gearMgr->setEcalEndcapParameters(param); + else if (name.find("Hcal")!=-1) m_gearMgr->setHcalEndcapParameters(param); + else if (name.find("Yoke")!=-1) m_gearMgr->setYokeEndcapParameters(param); + + info() << "EndcapParameters set" << name << endmsg; + } + else { + if (name.find("Ecal")!=-1) m_gearMgr->setEcalPlugParameters(param); + else if (name.find("Hcal")!=-1) m_gearMgr->setHcalRingParameters(param); + else if (name.find("Yoke")!=-1) m_gearMgr->setYokePlugParameters(param); + + info() << "Plug(Ring)Parameters set" << name << endmsg; + } + } + } + else{ + info() << name << " will convert in future!" << endmsg; + } + + return StatusCode::SUCCESS; +} + +TGeoNode* GearSvc::FindNode(TGeoNode* mother, char* name) { TGeoNode* next = 0; if(mother->GetNdaughters()!=0){ for(int i=0;i<mother->GetNdaughters();i++){ diff --git a/Service/GearSvc/src/GearSvc.h b/Service/GearSvc/src/GearSvc.h index 312491eb4c6f7e3950e50830f1faa87da08be119..3e2e250b9daf25feba27d6972a95f61f56b45ea0 100644 --- a/Service/GearSvc/src/GearSvc.h +++ b/Service/GearSvc/src/GearSvc.h @@ -25,6 +25,7 @@ class GearSvc : public extends<Service, IGearSvc> StatusCode convertDC (dd4hep::DetElement& dc); StatusCode convertSET(dd4hep::DetElement& set); StatusCode convertFTD(dd4hep::DetElement& ftd); + StatusCode convertCal(dd4hep::DetElement& cal); TGeoNode* FindNode(TGeoNode* mother, char* name); Gaudi::Property<std::string> m_gearFile{this, "GearXMLFile", ""}; diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/HelixFit.h b/Service/TrackSystemSvc/include/TrackSystemSvc/HelixFit.h index 6497531b35b1dd9cb4e2bf3a17761e40c6407d86..dc23662a4f36c580c0fab384ed30612209d7b7a2 100644 --- a/Service/TrackSystemSvc/include/TrackSystemSvc/HelixFit.h +++ b/Service/TrackSystemSvc/include/TrackSystemSvc/HelixFit.h @@ -48,7 +48,7 @@ namespace MarlinTrk { public: int fastHelixFit(int npt, double* xf, double* yf, float* rf, float* pf, double* wf, float* zf , float* wzf, int iopt, - float* vv0, float* ee0, float& ch2ph, float& ch2z); + float* vv0, float* ee0, float& ch2ph, float& ch2z, float MAX_CHI2=5000.0); }; diff --git a/Service/TrackSystemSvc/src/HelixFit.cc b/Service/TrackSystemSvc/src/HelixFit.cc index 070bbf6cdbd853dc5e932dddc626e2ef35918444..f68dd0355ba7cbce80b93d63fb998d38c9f0af40 100644 --- a/Service/TrackSystemSvc/src/HelixFit.cc +++ b/Service/TrackSystemSvc/src/HelixFit.cc @@ -17,12 +17,12 @@ namespace MarlinTrk{ int HelixFit::fastHelixFit(int npt, double* xf, double* yf, float* rf, float* pf, double* wf, float* zf , float* wzf,int iopt, - float* vv0, float* ee0, float& ch2ph, float& ch2z){ + float* vv0, float* ee0, float& ch2ph, float& ch2z, float MAX_CHI2){ if (npt < 3) { - //streamlog_out(ERROR) << "Cannot fit less than 3 points return 1" << std::endl; + //std::cerr << "Cannot fit less than 3 points return 1" << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; return 1; @@ -30,14 +30,23 @@ namespace MarlinTrk{ double eps = 1.0e-16; #define ITMAX 15 -#define MPT 600 -#define MAX_CHI2 5000.0 - - - float sp2[MPT], - del[MPT],deln[MPT],delzn[MPT],sxy[MPT],ss0[MPT],eee[MPT], - delz[MPT],grad[5],cov[15],vv1[5],dv[5]; - + //#define MPT 2500 + //#define MAX_CHI2 5000.0 + + + //float sp2[MPT], + //del[MPT],deln[MPT],delzn[MPT],sxy[MPT],ss0[MPT],eee[MPT], + //delz[MPT],grad[5],cov[15],vv1[5],dv[5]; + std::unique_ptr<float[], void(*)(float *)> sp2(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> del(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> deln(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> delzn(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> sxy(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> ss0(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> eee(new float[npt],[](float* p){delete p;}); + std::unique_ptr<float[], void(*)(float *)> delz(new float[npt],[](float* p){delete p;}); + float grad[5],cov[15],vv1[5],dv[5]; + // double xf[MPT],yf[MPT],wf[MPT],zf[MPT],wzf[MPT]; double alf,a0,a1,a2,a22,bet,cur, @@ -182,15 +191,10 @@ namespace MarlinTrk{ den2= 1.0/(x1*x1 + y1*y1 + gam*det*det); if(den2 <= 0.0) { - // streamlog_out(ERROR) << "den2 less than or equal to zero" - // << " x1 = " << x1 - // << " y1 = " << y1 - // << " gam = " << gam - // << " det = " << det - // << std::endl; + //std::cerr << "den2 less than or equal to zero x1 = " << x1 << " y1 = " << y1 << " gam = " << gam << " det = " << det << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; - return 1; + return 2; } den = sqrt(den2); @@ -216,10 +220,10 @@ namespace MarlinTrk{ rr0 = sst*cur; if( (alf*alf+bet*bet) <= 0.0 ){ - // streamlog_out(ERROR) << "(alf*alf+bet*bet) less than or equal to zero" << std::endl; + //std::cerr << "(alf*alf+bet*bet) less than or equal to zero" << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; - return 1; + return 3; } sa2b2 = 1.0/sqrt(alf*alf+bet*bet); @@ -343,10 +347,10 @@ namespace MarlinTrk{ double denom = sumw*sumss - sums*sums; if (fabs(denom) < eps){ - // streamlog_out(ERROR) << "fabs(denom) less than or equal to zero" << std::endl; + //std::cerr << "fabs(denom) less than or equal to zero" << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; - return 1; + return 4; } double dzds = (sumw*sumsz-sums*sumz) /denom; @@ -368,10 +372,10 @@ namespace MarlinTrk{ } if(chi2 > MAX_CHI2) { - // streamlog_out(ERROR) << "Chi2 greater than " << MAX_CHI2 << "return 1 " << std::endl; + //std::cerr << "Chi2 greater than " << MAX_CHI2 << "return 1 " << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; - return 1; + return 5; } for (int i = 0 ; i<npt; ++i) { @@ -464,10 +468,10 @@ namespace MarlinTrk{ cov0.invert(error); if( error != 0 ) { - //streamlog_out(ERROR) << "CLHEP Matrix inversion failed" << "return 1 " << std::endl; + //std::cerr << "CLHEP Matrix inversion failed" << "return 1 " << std::endl; ch2ph = 1.0e30; ch2z = 1.0e30; - return 1; + return 7; } icov = 0 ; @@ -558,20 +562,3 @@ namespace MarlinTrk{ } - - - - - - - - - - - - - - - - - diff --git a/Service/TrackSystemSvc/src/MarlinKalTest.cc b/Service/TrackSystemSvc/src/MarlinKalTest.cc index fa8d9fd4a3d0af7b0eb111109334d8a9c9aebd47..5931143ffd6bd50ba136097eb32f5bf95abf62e2 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTest.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTest.cc @@ -70,7 +70,7 @@ namespace MarlinTrk{ void MarlinKalTest::init() { - std::cout << "debug: MarlinKalTest - call this init " << std::endl ; + //std::cout << "debug: MarlinKalTest - call this init " << std::endl ; //ILDSITKalDetector* sitdet = new ILDSITKalDetector( *_gearMgr, _geoSvc ) ; MeasurementSurfaceStore& surfstore = _gearMgr->getMeasurementSurfaceStore(); @@ -84,7 +84,7 @@ namespace MarlinTrk{ } else { - std::cout << "debug: MarlinKalTest - MeasurementSurfaceStore is already full. Using store as filled by MeasurementSurfaceStoreFiller " << surfstore.getFillerName() << std::endl ; + //std::cout << "debug: MarlinKalTest - MeasurementSurfaceStore is already full. Using store as filled by MeasurementSurfaceStoreFiller " << surfstore.getFillerName() << std::endl ; } if (_gearMgr -> getDetectorName() == "LPTPC") { diff --git a/Utilities/KalDet/src/ild/common/ILDMeasurementSurfaceStoreFiller.cc b/Utilities/KalDet/src/ild/common/ILDMeasurementSurfaceStoreFiller.cc index 1d7e527aa65f823abb4cbbb852d549628704334f..16fe6fe18833afc685f325682dbe6461c56834cd 100644 --- a/Utilities/KalDet/src/ild/common/ILDMeasurementSurfaceStoreFiller.cc +++ b/Utilities/KalDet/src/ild/common/ILDMeasurementSurfaceStoreFiller.cc @@ -138,16 +138,17 @@ void ILDMeasurementSurfaceStoreFiller::getMeasurementSurfaces( std::vector<Measu if( _paramVXD ) this->storeZPlanar( _paramVXD , UTIL::ILDDetID::VXD, surface_list ); - + int nvxd = surface_list.size(); if( _paramSIT ) this->storeZPlanar( _paramSIT , UTIL::ILDDetID::SIT, surface_list ); - + int nsit = surface_list.size()-nvxd; if(_paramSET ) this->storeZPlanar( _paramSET , UTIL::ILDDetID::SET, surface_list ); - + int nset = surface_list.size()-nvxd-nsit; if( _paramFTD ) this->storeFTD( _paramFTD , surface_list); - + int nftd = surface_list.size()-nvxd-nsit-nset; + std::cout << "ILDMeasurementSurfaceStoreFiller::getMeasurementSurfaces " << nvxd << " " << nsit << " " << nset << " " << nftd << std::endl; } diff --git a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc index 64dd65463d1cf178c0c85e7d87bd079390591e21..1376b4ad50f7c070ff503b8dffd07f8a863a6d8e 100644 --- a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc +++ b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc @@ -87,7 +87,8 @@ TVKalDetector(10) // first check the sizes of the beam-pipe values r_in r_out and z if ( ( z.size() != rInner.size() ) || ( z.size() != rOuter.size() ) ) { - // streamlog_out(ERROR) << "ILDSupportKalDetector::ILDSupportKalDetector miss-match in number of double values for beam-pipe: Z = " << z.size() << " RInner = " << rInner.size() << " ROuter = " << rOuter.size() << " exit(1) called from "<< __FILE__ << " line " << __LINE__ << std::endl; + std::cout << "ILDSupportKalDetector::ILDSupportKalDetector miss-match in number of double values for beam-pipe: Z = " << z.size() + << " RInner = " << rInner.size() << " ROuter = " << rOuter.size() << " exit(1) called from "<< __FILE__ << " line " << __LINE__ << std::endl; exit(1); } @@ -127,7 +128,8 @@ TVKalDetector(10) if (i==0) { // check that z values start a zero if (fabs(zStart) > epsilon) { - // streamlog_out(ERROR) << "ILDSupportKalDetector::ILDSupportKalDetector first z value for beam-pipe must equal zero: Z = " << zStart << " exit(1) called from "<< __FILE__ << " line " << __LINE__ << std::endl; + std::cout << "ILDSupportKalDetector::ILDSupportKalDetector first z value for beam-pipe must equal zero: Z = " << zStart + << " exit(1) called from "<< __FILE__ << " line " << __LINE__ << std::endl; exit(1); } @@ -288,17 +290,17 @@ TVKalDetector(10) const gear::CalorimeterParameters& ecalB = gearMgr.getEcalBarrelParameters(); const gear::CalorimeterParameters& ecalE = gearMgr.getEcalEndcapParameters(); - - if (ecalB.getSymmetryOrder()!=8) { - // streamlog_out(ERROR) << "ILDSupportKalDetector::ILDSupportKalDetector ECal barrel is not eightfold symmetry: exit(1) called from " << __FILE__ << " line " << __LINE__ << std::endl; - exit(1); + int nside = ecalB.getSymmetryOrder(); + if (ecalB.getSymmetryOrder()<=0) { + std::cout << "ILDSupportKalDetector::ILDSupportKalDetector ECal barrel zero symmetry: exit(1) called from " << __FILE__ << " line " << __LINE__ << std::endl; + exit(1); } double phi0 = ecalB.getPhi0(); double r_min_ecal_bar = ecalB.getExtent()[0]; double z_max_ecal_bar = ecalB.getExtent()[3]; - + //std::cout << "fucdout ILDSupportKalDetector " << phi0 << " " << r_min_ecal_bar << " " << z_max_ecal_bar << std::endl; UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ; encoder.reset() ; // reset to 0 @@ -308,12 +310,12 @@ TVKalDetector(10) std::vector<int> module_ids; - for (int i=0; i<8; ++i) { + for (int i=0; i<nside; ++i) { encoder[lcio::ILDCellID0::module] = i; module_ids.push_back(encoder.lowWord()); - double segment_dphi = 2.0*M_PI / 8; + double segment_dphi = 2.0*M_PI / nside; double phi = i*segment_dphi+phi0; double width = 2.0*(r_min_ecal_bar*tan(segment_dphi*0.5)); double length = 2.0*z_max_ecal_bar;