diff --git a/.gitignore b/.gitignore index 7853163353f7e852e3f472d7a6675292940be6e0..ff7ac6feeaead84301a0b30a66a0445fe80e49cd 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,6 @@ build.* build spack* ./Generator/output/ -./Generator/options/ \ No newline at end of file +./Generator/options/ + +InstallArea/ \ No newline at end of file diff --git a/Analysis/TotalInvMass/src/TotalInvMass.cc b/Analysis/TotalInvMass/src/TotalInvMass.cc index 536739b25b1939bb841dbb54216c21671c78a28c..53eab8c6cbdfa6949d367ad46562917beaf561b3 100644 --- a/Analysis/TotalInvMass/src/TotalInvMass.cc +++ b/Analysis/TotalInvMass/src/TotalInvMass.cc @@ -11,13 +11,13 @@ #include <EVENT/LCFloatVec.h> #include <EVENT/LCParameters.h> #include <stdexcept> -#include <TFile.h> +#include <TFile.h> #include <TTree.h> #include <TH1F.h> #include <TVector3.h> #include <TRandom.h> -#include <Rtypes.h> -#include <sstream> +#include <Rtypes.h> +#include <sstream> #include <cmath> #include <vector> #include <TMath.h> @@ -47,14 +47,14 @@ TotalInvMass::TotalInvMass(const std::string& name, ISvcLocator* svcLoc) _colName , _colName); */ - + /* _leptonID = 13; registerProcessorParameter( "LeptonIDTag" , "Lepton ID that will be used in this analysis." , _leptonID , _leptonID); - */ + */ } @@ -74,7 +74,7 @@ StatusCode TotalInvMass::initialize() { _outputTree->SetAutoSave(32*1024*1024); // autosave every 32MB _outputTree->Branch("EventNr", &_eventNr, "EventNr/I"); _outputTree->Branch("Num", &_Num, "Num/I"); - + _outputTree->Branch("OriQuarkID", &_OriQuarkID, "OriQuarkID/I"); _outputTree->Branch("ISREn", &_ISREn, "ISREn/F"); @@ -88,6 +88,9 @@ StatusCode TotalInvMass::initialize() { _outputTree->Branch("N3En", &_N3En, "N3En/F"); _outputTree->Branch("N3Pt", &_N3Pt, "N3Pt/F"); + _outputTree->Branch("ArborPFO_E", &_ArborPFO_E); + _outputTree->Branch("ArborPFO_Charge", &_ArborPFO_Charge); + _outputTree->Branch("OriJ1CosTheta", &_OriJ1CosTheta, "OriJ1CosTheta/F"); _outputTree->Branch("OriJ2CosTheta", &_OriJ2CosTheta, "OriJ2CosTheta/F"); @@ -115,7 +118,7 @@ StatusCode TotalInvMass::initialize() { _outputTree->Branch("FrPh", FrPh, "FrPh[4]/F"); _outputTree->Branch("FrNe", FrNe, "FrNe[4]/F"); - + _outputTree->Branch("KPF", KPF, "KPF[4]/F"); _outputTree->Branch("UdP", UdP, "UdP[4]/F"); @@ -172,13 +175,16 @@ StatusCode TotalInvMass::initialize() { } -StatusCode TotalInvMass::execute() -{ +StatusCode TotalInvMass::execute() +{ info() << "TotalInvMass::executing..." << endmsg; EVENT::LCEvent* evtP = nullptr; - // if (evtP) { + MCParticleColHandler m_mcParticle{_MCPCollectionName, Gaudi::DataHandle::Reader, this}; + + + // if (evtP) { TLorentzVector ArborTotalP(0, 0, 0, 0); TLorentzVector ArborChP(0, 0, 0, 0); @@ -194,6 +200,9 @@ StatusCode TotalInvMass::execute() TLorentzVector ArborISR(0, 0, 0, 0); TLorentzVector PandoraISR(0, 0, 0, 0); + _ArborPFO_E.clear(); + _ArborPFO_Charge.clear(); + // TODO: using the event header // _eventNr = evtP->getEventNumber(); @@ -221,23 +230,23 @@ StatusCode TotalInvMass::execute() } } - nCHPFO_a = 0; + nCHPFO_a = 0; nCHPFO_p = 0; - nNEPFO_a = 0; + nNEPFO_a = 0; nNEPFO_p = 0; Type = -100; - Charge = -100; + Charge = -100; Energy = 0; - CluEn = 0; - TrkSumEn = 0; + CluEn = 0; + TrkSumEn = 0; - _EcalTotalE = 0; _HcalTotalE = 0; _EcalCluE = 0; _HcalCluE = 0; + _EcalTotalE = 0; _HcalTotalE = 0; _EcalCluE = 0; _HcalCluE = 0; _EcalCluE_p = 0; _HcalCluE_p = 0; _HcalEn1=0;_HcalEn2=0;_HcalEn3=0;_HcalEn4=0;_HcalEn5=0; _EcalEn1=0;_EcalEn2=0;_EcalEn3=0;_EcalEn4=0;_EcalEn5=0; - _HDPID = -1; - _OriQuarkID = 0; + _HDPID = -1; + _OriQuarkID = 0; _OQDir = -10; _HDir = -10; _visE=0; @@ -323,12 +332,12 @@ StatusCode TotalInvMass::execute() try{ auto MCPCol = m_mcParticle.get(); - TVector3 tmpP; + TVector3 tmpP; TVector3 ISRP(0, 0, 0); _N3En = 0; _N3Pt = 0; - - int NNeutrinoCount = 0; + + int NNeutrinoCount = 0; for (int s0 = 0; s0 < MCPCol->size(); ++s0) { auto MCP = (*MCPCol)[s0]; @@ -340,10 +349,20 @@ StatusCode TotalInvMass::execute() TVector3 VTX(VTX0.x, VTX0.y, VTX0.z); TVector3 EndP(EndP0.x, EndP0.y, EndP0.z); + int GenStatus = MCP.getGeneratorStatus(); + int SimuStatus = MCP.getSimulatorStatus(); + + if(0) debug()<<"MCP "<<s0<<": tmpPID = "<<tmpPID<<", NParent = "<<NParent<<", NDaughter = "<<NDaughter<<", GenStatus = "<<GenStatus<<", SimuStatus = "<<SimuStatus<<endmsg; + + // auto aDau = MCP.getDaughters(0); + // if(aDau){ + // debug()<<"Dau 0: PDG = "<<aDau->getPDG()<<endmsg; + + if(tmpPID == 22 && NParent == 0 && s0 < 4) { auto tmpP0 = MCP.getMomentum(); tmpP = TVector3(tmpP0.x, tmpP0.y, tmpP0.z); - ISRP += tmpP; + ISRP += tmpP; } if( (abs(tmpPID) == 12 || abs(tmpPID) == 14 || abs(tmpPID) == 16) && NParent != 0 && NDaughter == 0 && VTX.Mag() < 100 && EndP.Mag() > 100) { @@ -355,20 +374,25 @@ StatusCode TotalInvMass::execute() tmpP = TVector3(tmpP0.x, tmpP0.y, tmpP0.z); _N3En += tmpP.Mag(); _N3Pt += tmpP.Perp(); - cout<<"Found Neutrino: "<<NNeutrinoCount<<" En "<<_N3En<<" Pt "<<_N3Pt<<endl; + debug()<<"Found Neutrino: "<<NNeutrinoCount<<" En "<<_N3En<<" Pt "<<_N3Pt<<endmsg; } } if(tmpPID == 25 && NDaughter > 1 && NParent !=0 ) { //Higgs - // cout<<"tmpPID:HDPID"<<tmpPID<<" "<<NDaughter<<" "<<NParent<<endl; + // debug()<<"tmpPID:HDPID"<<tmpPID<<" "<<NDaughter<<" "<<NParent<<endmsg; _HDPID = abs(MCP.getDaughters(0).getPDG()); _HDir = MCP.getMomentum()[2]/MCP.getEnergy(); + // debug()<<"This is Higgs, has "<<NDaughter<<" daughters"<<endmsg; + debug()<<"This is Higgs, has "<<NDaughter<<" daughters"<<endmsg; + if(NDaughter == 2) { auto D1 = MCP.getDaughters(0); _J1CosTheta = D1.getMomentum()[2]/D1.getEnergy(); auto D2 = MCP.getDaughters(1); _J2CosTheta = D2.getMomentum()[2]/D2.getEnergy(); + + debug()<<"---> PDG: "<<D1.getPDG()<<" + "<<D2.getPDG()<<endmsg; } } if(abs(tmpPID)<7 && NParent == 0) { @@ -379,9 +403,9 @@ StatusCode TotalInvMass::execute() } if(abs(tmpPID)<7 && NParent == 0) { - _OriQuarkID = abs(tmpPID); + _OriQuarkID = abs(tmpPID); _OQDir = MCP.getMomentum()[2]/MCP.getEnergy(); - } + } if(abs(_J1CosTheta) > abs(_J2CosTheta)) @@ -413,17 +437,21 @@ StatusCode TotalInvMass::execute() auto a_RecoP = (*col_RecoNeP)[i0]; TLorentzVector currP( a_RecoP.getMomentum()[0], a_RecoP.getMomentum()[1], a_RecoP.getMomentum()[2], a_RecoP.getEnergy()); ArborTotalP += currP; + + _ArborPFO_E.push_back(a_RecoP.getEnergy()); + _ArborPFO_Charge.push_back(a_RecoP.getCharge()); + auto currMom0 = a_RecoP.getMomentum(); TVector3 currMom(currMom0.x, currMom0.y, currMom0.z); // if(a_RecoP->getType() == 22 && a_RecoP->getEnergy() > 2) //Compensate... // ArborTotalP += 0.98*currP; - // else + // else // ArborTotalP += currP; if(a_RecoP.getCharge() != 0) { ArborChP += currP; } else if(a_RecoP.getType() == 310) { - ArborKPF += currP; + ArborKPF += currP; } else if(a_RecoP.getType() == 22) { if(a_RecoP.getEnergy() < 3.0) ArborFrPh += currP; @@ -436,21 +464,36 @@ StatusCode TotalInvMass::execute() ArborNeP += currP; } else if(a_RecoP.getEnergy() < 3.0) { ArborFrP += currP; - cout<<"Undef "<<a_RecoP.getType()<<endl; + if(0) debug()<<"Undef "<<a_RecoP.getType()<<endmsg; } else { ArborUdP += currP; - cout<<"Undef "<<a_RecoP.getType() << "En "<<a_RecoP.getEnergy()<<endl; + if(0) debug()<<"Undef "<<a_RecoP.getType() << "En "<<a_RecoP.getEnergy()<<endmsg; } + + if(0) debug()<<"[YX debug - BMR] PFO "<<i0<<", E = "<<a_RecoP.getEnergy()<<", Charge = "<<a_RecoP.getCharge()<<endmsg; + + if(a_RecoP.clusters_size() > 0) { + + int nTrkHit = 0; + if(a_RecoP.getCharge()){ + auto a_Trk = a_RecoP.getTracks(0); + nTrkHit = a_Trk.trackerHits_size(); + } + auto a_clu = a_RecoP.getClusters(0); auto CluPos0 = a_clu.getPosition(); TVector3 CluPos(CluPos0.x, CluPos0.y, CluPos0.z); if( CluPos.Perp() < 300 && abs(CluPos.Z()) < 1300 && a_RecoP.getCharge() == 0 ) // 1150-1300 - ArborLCAL += currP; + ArborLCAL += currP; + + + if(0) debug()<<"[YX debug - BMR] ---> nTrkHit = "<<nTrkHit<<", CluE = "<<a_clu.getEnergy()<<", CluPos = ("<<CluPos.X()<<", "<<CluPos.Y()<<", "<<CluPos.Z()<<")"<<endmsg; + float MinAngleToCH = 1.0E6; - float MinAngleToNE = 1.0E6; + float MinAngleToNE = 1.0E6; if(a_RecoP.getEnergy() > 5 && a_RecoP.getCharge() == 0) { for(int i1 = 0; i1 < col_RecoNeP->size(); i1++) { @@ -464,7 +507,7 @@ StatusCode TotalInvMass::execute() if(tmpAngle < MinAngleToCH) { MinAngleToCH = tmpAngle; } - } else { + } else { if(tmpAngle < MinAngleToNE) { MinAngleToNE = tmpAngle; } @@ -474,19 +517,19 @@ StatusCode TotalInvMass::execute() } if( MinAngleToNE > 0.5 || MinAngleToCH > 0.5 ) { - ArborISR += currP; + ArborISR += currP; } } } - TrackHit = 0; + TrackHit = 0; for(int k = 0; k < 3; k++) { StartPos[k] = 0; EndPos[k] = 0; - } + } Charge = int(a_RecoP.getCharge()); - CluEn = 0; + CluEn = 0; CluEnCom[0] = 0; CluEnCom[1] = 0; @@ -510,7 +553,7 @@ StatusCode TotalInvMass::execute() if((!Charge) and (currClu.subdetectorEnergies_size() > 2)) { CluEnCom[0] = currClu.getSubdetectorEnergies(0); - CluEnCom[1] = currClu.getSubdetectorEnergies(1); + CluEnCom[1] = currClu.getSubdetectorEnergies(1); NeCaloE_a[0] += currClu.getSubdetectorEnergies(0); NeCaloE_a[1] += currClu.getSubdetectorEnergies(1); } @@ -522,7 +565,7 @@ StatusCode TotalInvMass::execute() if(Charge) { - TrkSumEn += Energy; + TrkSumEn += Energy; if (a_RecoP.tracks_size()>0) { auto a_Trk = a_RecoP.getTracks(0); @@ -535,19 +578,19 @@ StatusCode TotalInvMass::execute() EndPos[0] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[0]; EndPos[1] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[1]; - EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2]; + EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2]; } } if( Energy > CluEn + sqrt(Energy)) { - ElargeP[0] += Energy; - ElargeP[1] += CluEn; + ElargeP[0] += Energy; + ElargeP[1] += CluEn; } else if( fabs(Energy - CluEn) < sqrt(Energy) ) { EequP[0] += Energy; - EequP[1] += CluEn; + EequP[1] += CluEn; } else { EsmallP[0] += Energy; - EsmallP[1] += CluEn; + EsmallP[1] += CluEn; } } @@ -559,18 +602,18 @@ StatusCode TotalInvMass::execute() }catch (lcio::DataNotAvailableException err) { } try{ - //LCCollection* col_RecoPandora = evtP->getCollection( "PandoraPFOs" ); + //LCCollection* col_RecoPandora = evtP->getCollection( "PandoraPFOs" ); //for(int i2 = 0; i2 < col_RecoPandora->getNumberOfElements(); i2++) for(int s = 0; s < 1; s++) { auto col_PFO_iter = m_arbopfo.get(); - /* + /* if(s==0) col_PFO_iter = evtP->getCollection( "ArborChargedCore" ); else - col_PFO_iter = evtP->getCollection( "ArborNeutralCore" ); + col_PFO_iter = evtP->getCollection( "ArborNeutralCore" ); */ - for(int i2 = 0; i2 < col_PFO_iter->size(); i2++) { + for(int i2 = 0; i2 < col_PFO_iter->size(); i2++) { auto a_RecoP = (*col_PFO_iter)[i2]; TLorentzVector currP( a_RecoP.getMomentum()[0], a_RecoP.getMomentum()[1], a_RecoP.getMomentum()[2], a_RecoP.getEnergy()); PandoraTotalP += currP; @@ -580,7 +623,7 @@ StatusCode TotalInvMass::execute() } else { nNEPFO_p++; } - + auto currMom0 = a_RecoP.getMomentum(); TVector3 currMom(currMom0.x, currMom0.y, currMom0.z); @@ -592,7 +635,7 @@ StatusCode TotalInvMass::execute() auto currClu = a_RecoP.getClusters(0); CluEn = currClu.getEnergy(); - _EcalCluE_p += CluEn; + _EcalCluE_p += CluEn; _HcalCluE_p += currClu.getSubdetectorEnergies(1); if(a_RecoP.getEnergy() > 5 && a_RecoP.getCharge() == 0) { @@ -625,11 +668,11 @@ StatusCode TotalInvMass::execute() } }catch (lcio::DataNotAvailableException err) { } - _Mass_a = 0; - _Mass_p = 0; + _Mass_a = 0; + _Mass_p = 0; _Mass_a_Pisr = 0; - _Mass_a_Plcal = 0; - _Mass_p_Pisr = 0; + _Mass_a_Plcal = 0; + _Mass_p_Pisr = 0; _Mass_a = ArborTotalP.M(); _Mass_p = PandoraTotalP.M(); @@ -656,7 +699,7 @@ StatusCode TotalInvMass::execute() PhP[0] = ArborPhP.X(); PhP[1] = ArborPhP.Y(); PhP[2] = ArborPhP.Z(); - PhP[3] = ArborPhP.T(); + PhP[3] = ArborPhP.T(); NeP[0] = ArborNeP.X(); NeP[1] = ArborNeP.Y(); @@ -688,18 +731,18 @@ StatusCode TotalInvMass::execute() KPF[2] = ArborKPF.Z(); KPF[3] = ArborKPF.T(); - cout<<_Mass_a<<" : "<<_Mass_p<<endl; + info()<<_Mass_a<<" : "<<_Mass_p<<endmsg; _outputTree->Fill(); _Num++; - // } + // } info() << "TotalInvMass::execute done" << endmsg; return StatusCode::SUCCESS; -} +} StatusCode TotalInvMass::finalize() { diff --git a/Analysis/TotalInvMass/src/TotalInvMass.hh b/Analysis/TotalInvMass/src/TotalInvMass.hh index 97d48d5ad80a82aa59ab6c3769b04bbae016afd8..4084b6133a1cf21a813e843dccc1acb29c4df07a 100644 --- a/Analysis/TotalInvMass/src/TotalInvMass.hh +++ b/Analysis/TotalInvMass/src/TotalInvMass.hh @@ -34,7 +34,11 @@ public: protected: typedef DataHandle<edm4hep::MCParticleCollection> MCParticleColHandler; - MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this}; + // MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this}; + // MCParticleColHandler m_mcParticle{"MCParticleGen", Gaudi::DataHandle::Reader, this}; + Gaudi::Property<std::string> _MCPCollectionName{this, + "MCPCollectionName", "MCParticle", + "MCPCollectionName"}; typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloHitColHandler; CaloHitColHandler m_ecalbarrelhitcol{"ECALBarrel", Gaudi::DataHandle::Reader, this}; @@ -45,8 +49,9 @@ protected: CaloHitColHandler m_hcalotherhitcol {"HCALOther", Gaudi::DataHandle::Reader, this}; typedef DataHandle<edm4hep::ReconstructedParticleCollection> RecParticleColHandler; - RecParticleColHandler m_reconep{"AncientPFOs", Gaudi::DataHandle::Reader, this}; - RecParticleColHandler m_arbopfo{"ArborLICHPFOs", Gaudi::DataHandle::Reader, this}; + RecParticleColHandler m_reconep{"ArborPFO", Gaudi::DataHandle::Reader, this}; + RecParticleColHandler m_arbopfo{"ArborPFO", Gaudi::DataHandle::Reader, this}; + Gaudi::Property<std::string> _treeFileName{this, "TreeOutputFile", "MCTruth.root", @@ -61,36 +66,36 @@ protected: "OverwriteFile", 0, "If zero an already existing file will not be overwritten."}; int _leptonID; - float _cmsE; + float _cmsE; TTree *_outputTree, *_outputPFO; float _ISRP[3]; - float _ISREn, _ISRPt; - float _N3En, _N3Pt; - float _NEn, _NPt; + float _ISREn, _ISRPt; + float _N3En, _N3Pt; + float _NEn, _NPt; - int _HDPID, _OriQuarkID; + int _HDPID, _OriQuarkID; float _OQDir, _HDir; int _NMuP, _NMuM, _NChP, _NChM; float _P_MuP[4], _P_MuM[4], _P_DL[4]; - int _EventType; - float _InvMass, _RecoilMass; - float _J1CosTheta, _J2CosTheta, _MaxJetCosTheta; - float _OriJ1CosTheta, _OriJ2CosTheta, _MaxOriJetCosTheta; - float _Mass_a, _Mass_p; + int _EventType; + float _InvMass, _RecoilMass; + float _J1CosTheta, _J2CosTheta, _MaxJetCosTheta; + float _OriJ1CosTheta, _OriJ2CosTheta, _MaxOriJetCosTheta; + float _Mass_a, _Mass_p; int _PID1, _PID2; float _PL1[4], _PL2[4], _RPL1[4], _RPL2[4], _SM[4], _P_allCharged[4], _P_allNeutral[4], _P_Higgs[4], _P_allReco[4]; - float _Hmass; + float _Hmass; int _Num; - int _NHDaug; - int _HdaughterPID; - int _ZdaughterPID; + int _NHDaug; + int _HdaughterPID; + int _ZdaughterPID; float _Pz[4], _Ph[4], _PzD1[4], _PzD2[4], _PhD1[4], _PhD2[4], _RPzD1[4], _RPzD2[4], _RPhD1[4], _RPhD2[4]; float _P[4], _SumP[4], _VisP[4], _MissP[4]; - int _PID, _NFMCP, _MotherFlag, _NNeutrino; - float _ENeutrino, _DiPhMass, _DiPhMassCorr; + int _PID, _NFMCP, _MotherFlag, _NNeutrino; + float _ENeutrino, _DiPhMass, _DiPhMassCorr; // float _CosTheta, _Phi, _Charge; - float _Mz, _Mrecoil, _MzReco, _MhReco, _MrecoilReco; + float _Mz, _Mrecoil, _MzReco, _MhReco, _MrecoilReco; float KthEn[7][9]; unsigned int _eventNr; @@ -102,12 +107,12 @@ protected: float ElargeP[2], EequP[2], EsmallP[2]; float ChP[4], FrP[4], PhP[4], NeP[4], UdP[4], FrPh[4], FrNe[4], KPF[4]; - float _EcalTotalE, _HcalTotalE, _EcalCluE, _HcalCluE, _EcalCluE_p, _HcalCluE_p; + float _EcalTotalE, _HcalTotalE, _EcalCluE, _HcalCluE, _EcalCluE_p, _HcalCluE_p; float _HcalEn1, _HcalEn2,_HcalEn3,_HcalEn4,_HcalEn5; float _EcalEn1, _EcalEn2,_EcalEn3,_EcalEn4,_EcalEn5; int Type, Charge; - float Energy, TrkSumEn; + float Energy, TrkSumEn; float P[3], CluEnCom[2]; float CluEn; @@ -115,6 +120,17 @@ protected: float StartPos[3], EndPos[3]; float _visE; + std::vector<int> _ArborPFO_Charge; + std::vector<double> _ArborPFO_E; + // std::vector<double> _ArborPFO_CosTheta; + // std::vector<int> _ArborPFO_nTrk; + // std::vector<int> _ArborPFO_nClu; + // std::vector<TVector3> _ArborPFO_TrkP3; + // std::vector<double> _ArborPFO_CluE; + // std::vector<double> _ArborPFO_CluHitESum; + // std::vector<double> _ArborPFO_CluHitESum; + + std::string _fileName; std::ostream *_output; std::string _histFileName; diff --git a/Detector/DetCEPCv4/compact/CEPCV4_FullDet_GearOutput.xml b/Detector/DetCEPCv4/compact/CEPCV4_FullDet_GearOutput.xml new file mode 100644 index 0000000000000000000000000000000000000000..c434ae54c043be5f20f21c1881b339b255cd4b04 --- /dev/null +++ b/Detector/DetCEPCv4/compact/CEPCV4_FullDet_GearOutput.xml @@ -0,0 +1,334 @@ +<gear> + <global detectorName="CEPC_v4" /> + <!--Gear XML file automatically created with GearXML::createXMLFile ....--> + <BField type="ConstantBField" x="0.000000000e+00" y="0.000000000e+00" z="3.000000000e+00" /> + <detectors> + <detector geartype="TPCParameters" name="TPC"> + <maxDriftLength value="2.225000000e+03" /> + <driftVelocity value="0.000000000e+00" /> + <coordinateType value="polar" /> + <modules> + <module> + <moduleID value="0" /> + <readoutFrequency value="0.000000000e+00" /> + <PadRowLayout2D type="FixedPadSizeDiskLayout" rMin="3.840000000e+02" rMax="1.716000000e+03" padHeight="6.000000000e+00" padWidth="1.000000000e+00" maxRow="222" padGap="0.000000000e+00" phiMax="6.283185307e+00" /> + <offset x_r="0.000000000e+00" y_phi="0.000000000e+00" /> + <angle value="0.000000000e+00" /> + <enlargeActiveAreaBy value="0.000000000e+00" /> + </module> + </modules> + <parameter name="TPCGasProperties_RadLen" type="double" value="1.155205461e+05" /> + <parameter name="TPCGasProperties_dEdx" type="double" value="2.668179899e-07" /> + <parameter name="TPCInnerWallProperties_RadLen" type="double" value="2.740688665e+03" /> + <parameter name="TPCInnerWallProperties_dEdx" type="double" value="1.241647394e-05" /> + <parameter name="TPCOuterWallProperties_RadLen" type="double" value="6.495646008e+03" /> + <parameter name="TPCOuterWallProperties_dEdx" type="double" value="5.300932694e-06" /> + <parameter name="TPCWallProperties_RadLen" type="double" value="2.740688665e+03" /> + <parameter name="TPCWallProperties_dEdx" type="double" value="1.241647394e-05" /> + <parameter name="tpcInnerRadius" type="double" value="3.290000000e+02" /> + <parameter name="tpcInnerWallThickness" type="double" value="2.500000000e+01" /> + <parameter name="tpcOuterRadius" type="double" value="1.808000000e+03" /> + <parameter name="tpcOuterWallThickness" type="double" value="6.000000000e+01" /> + <parameter name="tpcZAnode" type="double" value="4.600000000e+03" /> + </detector> + <detector name="EcalBarrel" geartype="CalorimeterParameters"> + <layout type="Barrel" symmetry="8" phi0="0.000000000e+00" /> + <dimensions inner_r="1.847415655e+03" outer_z="2.350000000e+03" /> + <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + </detector> + <detector name="EcalEndcap" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="4.000000000e+02" outer_r="2.088800000e+03" inner_z="2.450000000e+03" /> + <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + </detector> + <detector name="EcalPlug" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="2.400000000e+02" outer_r="4.000000000e+02" inner_z="2.450000000e+03" /> + <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" /> + </detector> + <detector name="YokeBarrel" geartype="CalorimeterParameters"> + <layout type="Barrel" symmetry="12" phi0="0.000000000e+00" /> + <dimensions inner_r="4.173929932e+03" outer_z="4.072000000e+03" /> + <layer repeat="1" thickness="4.000000000e+01" absorberThickness="0.000000000e+00" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + <layer repeat="9" thickness="1.400000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + <layer repeat="3" thickness="6.000000000e+02" absorberThickness="5.600000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + <layer repeat="1" thickness="4.000000000e+01" absorberThickness="4.027623145e-320" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + </detector> + <detector name="YokeEndcap" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="3.200000000e+02" outer_r="7.414929932e+03" inner_z="4.072000000e+03" /> + <layer repeat="1" thickness="1.000000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + <layer repeat="9" thickness="1.400000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + <layer repeat="2" thickness="6.000000000e+02" absorberThickness="5.600000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" /> + </detector> + <detector name="YokePlug" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="3.200000000e+02" outer_r="2.849254326e+03" inner_z="3.781430000e+03" /> + <parameter name="YokePlugThickness" type="double" value="2.905700000e+02" /> + </detector> + <detector name="HcalBarrel" geartype="CalorimeterParameters"> + <layout type="Barrel" symmetry="8" phi0="1.570796327e+00" /> + <dimensions inner_r="2.058000000e+03" outer_z="2.350000000e+03" /> + <layer repeat="40" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" /> + <parameter name="Hcal_barrel_number_modules" type="int" value="5" /> + <parameter name="N_cells_z" type="int" value="91" /> + <parameter name="FrameWidth" type="double" value="1.000000000e+00" /> + <parameter name="Hcal_lateral_structure_thickness" type="double" value="1.000000000e+01" /> + <parameter name="Hcal_modules_gap" type="double" value="2.000000000e+00" /> + <parameter name="Hcal_outer_radius" type="double" value="3.144432447e+03" /> + <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" /> + <parameter name="InnerOctoSize" type="double" value="1.704903012e+03" /> + <parameter name="RPC_PadSeparation" type="double" value="0.000000000e+00" /> + <parameter name="TPC_Ecal_Hcal_barrel_halfZ" type="double" value="2.350000000e+03" /> + </detector> + <detector name="HcalEndcap" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="3.500000000e+02" outer_r="3.144432447e+03" inner_z="2.650000000e+03" /> + <layer repeat="40" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" /> + <parameter name="FrameWidth" type="double" value="0.000000000e+00" /> + <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" /> + </detector> + <detector name="HcalRing" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" /> + <dimensions inner_r="2.138800000e+03" outer_r="3.144432447e+03" inner_z="2.450000000e+03" /> + <layer repeat="6" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" /> + <parameter name="FrameWidth" type="double" value="0.000000000e+00" /> + <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" /> + </detector> + <detector name="Lcal" geartype="CalorimeterParameters"> + <layout type="Endcap" symmetry="1" phi0="0.000000000e+00" /> + <dimensions inner_r="3.225828541e+01" outer_r="9.880000000e+01" inner_z="9.519000000e+02" /> + <layer repeat="30" thickness="4.290000000e+00" absorberThickness="3.500000000e+00" cellSize0="1.039714290e+00" cellSize1="1.308996939e-01" /> + <parameter name="beam_crossing_angle" type="double" value="0.000000000e+00" /> + </detector> + <detector name="VXD" geartype="ZPlanarParameters"> + <type technology="HYBRID" /> + <shell halfLength="1.450000000e+02" gap="0.000000000e+00" innerRadius="6.500000000e+01" outerRadius="6.549392000e+01" radLength="3.527597571e+02" /> + <layers> + <layer nLadders="10" phi0="-1.570796327e+00"> + <ladder distance="1.600000000e+01" thickness="1.000000000e+00" width="1.150000000e+01" length="6.250000000e+01" offset="-1.874869853e+00" radLength="1.014262421e+03" /> + <sensitive distance="1.595000000e+01" thickness="5.000000000e-02" width="1.100000000e+01" length="6.250000000e+01" offset="-1.624869853e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="10" phi0="-1.570796327e+00"> + <ladder distance="1.700000000e+01" thickness="1.000000000e+00" width="1.150000000e+01" length="6.250000000e+01" offset="-1.874869853e+00" radLength="1.014262421e+03" /> + <sensitive distance="1.800000000e+01" thickness="5.000000000e-02" width="1.100000000e+01" length="6.250000000e+01" offset="-1.624869853e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="11" phi0="-1.570796327e+00"> + <ladder distance="3.700000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-1.837940563e+00" radLength="1.014262421e+03" /> + <sensitive distance="3.695000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-1.587940563e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="11" phi0="-1.570796327e+00"> + <ladder distance="3.800000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-1.837940563e+00" radLength="1.014262421e+03" /> + <sensitive distance="3.900000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-1.587940563e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="17" phi0="-1.570796327e+00"> + <ladder distance="5.800000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-2.636744400e+00" radLength="1.014262421e+03" /> + <sensitive distance="5.795000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-2.386744400e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="17" phi0="-1.570796327e+00"> + <ladder distance="5.900000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-2.636744400e+00" radLength="1.014262421e+03" /> + <sensitive distance="6.000000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-2.386744400e+00" radLength="9.366070445e+01" /> + </layer> + </layers> + </detector> + <detector name="FTD" geartype="FTDParameters"> + <layers> + <layer nPetals="16" nSensors="1" isDoubleSided="0" sensorType="PIXEL" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="2.205200000e+02"> + <support thickness="1.000000000e+00" width="1.224000000e+02" lengthMin="1.173582968e+01" lengthMax="6.042957721e+01" rInner="2.950000000e+01" radLength="2.807467352e+02" /> + <sensitive thickness="2.000000000e-02" width="1.224000000e+02" lengthMin="1.173582968e+01" lengthMax="6.042957721e+01" rInner="2.950000000e+01" radLength="9.366070445e+01" /> + </layer> + <layer nPetals="16" nSensors="1" isDoubleSided="0" sensorType="PIXEL" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="3.715200000e+02"> + <support thickness="1.000000000e+00" width="1.213600000e+02" lengthMin="1.214956740e+01" lengthMax="6.042957721e+01" rInner="3.054000000e+01" radLength="2.807467352e+02" /> + <sensitive thickness="2.000000000e-02" width="1.213600000e+02" lengthMin="1.214956740e+01" lengthMax="6.042957721e+01" rInner="3.054000000e+01" radLength="9.366070445e+01" /> + </layer> + <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="6.462000000e+02"> + <support thickness="2.000000000e+00" width="2.665000000e+02" lengthMin="1.292930388e+01" lengthMax="1.189495957e+02" rInner="3.250000000e+01" radLength="2.807467352e+02" /> + <sensitive thickness="2.000000000e-01" width="2.665000000e+02" lengthMin="1.292930388e+01" lengthMax="1.189495957e+02" rInner="3.250000000e+01" radLength="9.366070445e+01" /> + </layer> + <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="8.472000000e+02"> + <support thickness="2.000000000e+00" width="2.750000000e+02" lengthMin="1.352604098e+01" lengthMax="1.229278430e+02" rInner="3.400000000e+01" radLength="2.807467352e+02" /> + <sensitive thickness="2.000000000e-01" width="2.750000000e+02" lengthMin="1.352604098e+01" lengthMax="1.229278430e+02" rInner="3.400000000e+01" radLength="9.366070445e+01" /> + </layer> + <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="9.262000000e+02"> + <support thickness="2.000000000e+00" width="2.735000000e+02" lengthMin="1.412277808e+01" lengthMax="1.229278430e+02" rInner="3.550000000e+01" radLength="2.807467352e+02" /> + <sensitive thickness="2.000000000e-01" width="2.735000000e+02" lengthMin="1.412277808e+01" lengthMax="1.229278430e+02" rInner="3.550000000e+01" radLength="9.366070445e+01" /> + </layer> + </layers> + <parameter name="strip_angle_deg" type="double" value="5.000000000e+00" /> + <parameter name="strip_length_mm" type="double" value="1.600000000e+03" /> + <parameter name="strip_pitch_mm" type="double" value="1.000000000e-02" /> + <parameter name="strip_width_mm" type="double" value="1.000000000e-03" /> + </detector> + <detector name="SIT" geartype="ZPlanarParameters"> + <type technology="CCD" /> + <shell halfLength="0.000000000e+00" gap="0.000000000e+00" innerRadius="0.000000000e+00" outerRadius="0.000000000e+00" radLength="0.000000000e+00" /> + <layers> + <layer nLadders="10" phi0="0.000000000e+00"> + <ladder distance="1.531000000e+02" thickness="1.000000000e+00" width="9.916044311e+01" length="3.680000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="1.529000000e+02" thickness="2.000000000e-01" width="9.916044311e+01" length="3.680000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="10" phi0="0.000000000e+00"> + <ladder distance="1.544000000e+02" thickness="1.000000000e+00" width="1.001352022e+02" length="3.680000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="1.554000000e+02" thickness="2.000000000e-01" width="1.001352022e+02" length="3.680000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="19" phi0="0.000000000e+00"> + <ladder distance="3.001000000e+02" thickness="1.000000000e+00" width="9.988891763e+01" length="6.440000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="2.999000000e+02" thickness="2.000000000e-01" width="9.988891763e+01" length="6.440000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="19" phi0="0.000000000e+00"> + <ladder distance="3.014000000e+02" thickness="1.000000000e+00" width="1.003895291e+02" length="6.440000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="3.024000000e+02" thickness="2.000000000e-01" width="1.003895291e+02" length="6.440000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + </layers> + <parameter name="sensor_length_mm" type="double" value="9.200000000e+01" /> + <parameter name="strip_angle_deg" type="double" value="7.000000000e+00" /> + <parameter name="strip_length_mm" type="double" value="9.200000000e+01" /> + <parameter name="strip_pitch_mm" type="double" value="5.000000000e-02" /> + <parameter name="strip_width_mm" type="double" value="1.250000000e-02" /> + <parameter name="n_sensors_per_ladder" type="IntVec" value="8 8 14 14" /> + </detector> + <detector name="SET" geartype="ZPlanarParameters"> + <type technology="CCD" /> + <shell halfLength="0.000000000e+00" gap="0.000000000e+00" innerRadius="0.000000000e+00" outerRadius="0.000000000e+00" radLength="0.000000000e+00" /> + <layers> + <layer nLadders="24" phi0="0.000000000e+00"> + <ladder distance="1.811100000e+03" thickness="1.000000000e+00" width="4.766190158e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="1.810900000e+03" thickness="2.000000000e-01" width="4.766190158e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + <layer nLadders="24" phi0="0.000000000e+00"> + <ladder distance="1.812400000e+03" thickness="1.000000000e+00" width="4.770139733e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="2.134851878e+02" /> + <sensitive distance="1.813400000e+03" thickness="2.000000000e-01" width="4.770139733e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="9.366070445e+01" /> + </layer> + </layers> + <parameter name="sensor_length_mm" type="double" value="9.200000000e+01" /> + <parameter name="strip_angle_deg" type="double" value="7.000000000e+00" /> + <parameter name="strip_length_mm" type="double" value="9.200000000e+01" /> + <parameter name="strip_pitch_mm" type="double" value="5.000000000e-02" /> + <parameter name="strip_width_mm" type="double" value="1.250000000e-02" /> + <parameter name="n_sensors_per_ladder" type="IntVec" value="50 50" /> + </detector> + <detector name="BeamPipe" geartype="GearParameters"> + <parameter name="BeamPipeHalfZ" type="double" value="7.300000000e+02" /> + <parameter name="BeamPipeProperties_RadLen" type="double" value="3.527597571e+02" /> + <parameter name="BeamPipeProperties_dEdx" type="double" value="2.941795296e-04" /> + <parameter name="BeamPipeRadius" type="double" value="1.400000000e+01" /> + <parameter name="BeamPipeThickness" type="double" value="5.000000000e-01" /> + <parameter name="RInner" type="DoubleVec" value="1.400000000e+01 1.400000000e+01 2.500000000e+00 1.300000000e+01 1.300000000e+01 1.300000000e+01 1.550000000e+01 1.550000000e+01 1.900000000e+01 1.900000000e+01 2.500000000e+01 2.500000000e+01 1.300000000e+01 1.300000000e+01 2.050000000e+01 2.050000000e+01 2.300000000e+01 2.300000000e+01 2.600000000e+01 2.600000000e+01 3.200000000e+01 3.200000000e+01" /> + <parameter name="ROuter" type="DoubleVec" value="1.450000000e+01 1.450000000e+01 1.800000000e+01 1.800000000e+01 1.550000000e+01 1.550000000e+01 1.900000000e+01 1.900000000e+01 2.500000000e+01 2.500000000e+01 3.300000000e+01 3.300000000e+01 1.550000000e+01 1.550000000e+01 2.300000000e+01 2.300000000e+01 2.600000000e+01 2.600000000e+01 3.200000000e+01 3.200000000e+01 4.000000000e+01 4.000000000e+01" /> + <parameter name="Z" type="DoubleVec" value="0.000000000e+00 5.000000000e+02 7.000000000e+02 7.010000000e+02 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 3.950000000e+03 3.950000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03" /> + </detector> + <detector name="CoilParameters" geartype="GearParameters"> + <parameter name="Coil_cryostat_c_modules_half_z" type="double" value="1.224000000e+03" /> + <parameter name="Coil_cryostat_c_modules_inner_radius" type="double" value="3.348930000e+03" /> + <parameter name="Coil_cryostat_c_modules_outer_radius" type="double" value="3.599930000e+03" /> + <parameter name="Coil_cryostat_half_z" type="double" value="3.872000000e+03" /> + <parameter name="Coil_cryostat_inner_cyl_half_z" type="double" value="3.872000000e+03" /> + <parameter name="Coil_cryostat_inner_cyl_inner_radius" type="double" value="3.173930000e+03" /> + <parameter name="Coil_cryostat_inner_cyl_outer_radius" type="double" value="3.963930000e+03" /> + <parameter name="Coil_cryostat_inner_radius" type="double" value="3.173930000e+03" /> + <parameter name="Coil_cryostat_mandrel_half_z" type="double" value="3.675000000e+03" /> + <parameter name="Coil_cryostat_mandrel_inner_radius" type="double" value="3.599930000e+03" /> + <parameter name="Coil_cryostat_mandrel_outer_radius" type="double" value="3.827930000e+03" /> + <parameter name="Coil_cryostat_modules_half_z" type="double" value="7.960000000e+02" /> + <parameter name="Coil_cryostat_modules_inner_radius" type="double" value="3.348930000e+03" /> + <parameter name="Coil_cryostat_modules_outer_radius" type="double" value="3.599930000e+03" /> + <parameter name="Coil_cryostat_outer_cyl_half_z" type="double" value="3.872000000e+03" /> + <parameter name="Coil_cryostat_outer_cyl_inner_radius" type="double" value="3.893930000e+03" /> + <parameter name="Coil_cryostat_outer_cyl_outer_radius" type="double" value="3.923930000e+03" /> + <parameter name="Coil_cryostat_outer_radius" type="double" value="3.923930000e+03" /> + <parameter name="Coil_cryostat_scint1_inner_radius" type="double" value="3.263930000e+03" /> + <parameter name="Coil_cryostat_scint1_outer_radius" type="double" value="3.273930000e+03" /> + <parameter name="Coil_cryostat_scint1_zposend" type="double" value="3.972000000e+03" /> + <parameter name="Coil_cryostat_scint1_zposin" type="double" value="3.772000000e+03" /> + <parameter name="Coil_cryostat_scint2_inner_radius" type="double" value="3.278930000e+03" /> + <parameter name="Coil_cryostat_scint2_outer_radius" type="double" value="3.288930000e+03" /> + <parameter name="Coil_cryostat_scint2_zposend" type="double" value="3.972000000e+03" /> + <parameter name="Coil_cryostat_scint2_zposin" type="double" value="3.772000000e+03" /> + <parameter name="Coil_cryostat_scint3_inner_radius" type="double" value="3.833930000e+03" /> + <parameter name="Coil_cryostat_scint3_outer_radius" type="double" value="3.843930000e+03" /> + <parameter name="Coil_cryostat_scint3_zposend" type="double" value="3.972000000e+03" /> + <parameter name="Coil_cryostat_scint3_zposin" type="double" value="3.772000000e+03" /> + <parameter name="Coil_cryostat_scint4_inner_radius" type="double" value="3.818930000e+03" /> + <parameter name="Coil_cryostat_scint4_outer_radius" type="double" value="3.828930000e+03" /> + <parameter name="Coil_cryostat_scint4_zposend" type="double" value="3.972000000e+03" /> + <parameter name="Coil_cryostat_scint4_zposin" type="double" value="3.772000000e+03" /> + <parameter name="Coil_cryostat_side_l_half_z" type="double" value="2.500000000e+01" /> + <parameter name="Coil_cryostat_side_l_inner_radius" type="double" value="3.213930000e+03" /> + <parameter name="Coil_cryostat_side_l_outer_radius" type="double" value="3.893930000e+03" /> + <parameter name="Coil_cryostat_side_r_half_z" type="double" value="2.500000000e+01" /> + <parameter name="Coil_cryostat_side_r_inner_radius" type="double" value="3.213930000e+03" /> + <parameter name="Coil_cryostat_side_r_outer_radius" type="double" value="3.893930000e+03" /> + <parameter name="Coil_material_c_modules" type="string" value="aluminium" /> + <parameter name="Coil_material_inner_cyl" type="string" value="aluminium" /> + <parameter name="Coil_material_mandrel" type="string" value="aluminium" /> + <parameter name="Coil_material_modules" type="string" value="aluminium" /> + <parameter name="Coil_material_outer_cyl" type="string" value="aluminium" /> + <parameter name="Coil_material_scint1" type="string" value="polystyrene" /> + <parameter name="Coil_material_scint2" type="string" value="polystyrene" /> + <parameter name="Coil_material_scint3" type="string" value="polystyrene" /> + <parameter name="Coil_material_scint4" type="string" value="polystyrene" /> + <parameter name="Coil_material_side_l" type="string" value="aluminium" /> + <parameter name="Coil_material_side_r" type="string" value="aluminium" /> + </detector> + <detector name="MokkaParameters" geartype="GearParameters"> + <parameter name="Ecal_endcap_outer_radius" type="string" value="2088.8" /> + <parameter name="Ecal_endcap_plug_rmin" type="string" value="240" /> + <parameter name="Ecal_endcap_zmax" type="string" value="2635" /> + <parameter name="Ecal_endcap_zmin" type="string" value="2450" /> + <parameter name="Ecal_outer_radius" type="string" value="2028" /> + <parameter name="Hcal_R_max" type="string" value="3144.43" /> + <parameter name="Hcal_endcap_zmin" type="string" value="2650" /> + <parameter name="Lcal_z_begin" type="string" value="951.9" /> + <parameter name="Lcal_z_thickness" type="string" value="128.1" /> + <parameter name="MokkaModel" type="string" value="CEPC_v4" /> + <parameter name="MokkaVersion" type="string" value="void" /> + <parameter name="SIT1_Half_Length_Z" type="string" value="368" /> + <parameter name="SIT1_Radius" type="string" value="152.9" /> + <parameter name="SIT2_Half_Length_Z" type="string" value="644" /> + <parameter name="SIT2_Radius" type="string" value="299.9" /> + <parameter name="SiTrackerEndcap" type="string" value="FTD_PIXEL,29.5,151.9,220,16;FTD_PIXEL,30.54,151.9,371,16;FTD_STRIP,32.5,299,645,16;FTD_STRIP,34,309,846,16;FTD_STRIP,35.5,309,925,16" /> + <parameter name="SiTrackerLayerStructure" type="string" value="FTD_PIXEL,Si:-0.02,CarbonFiber:1;FTD_STRIP,Si:-0.2,CarbonFiber:2,Si:-0.2" /> + <parameter name="TPC_Ecal_Hcal_barrel_halfZ" type="string" value="2350" /> + <parameter name="Yoke_Z_start_endcaps" type="string" value="4072" /> + <parameter name="Yoke_barrel_inner_radius" type="string" value="4173.929931640625" /> + <parameter name="calorimeter_region_rmax" type="string" value="3144.43" /> + <parameter name="calorimeter_region_zmax" type="string" value="3736.43" /> + <parameter name="tracker_region_rmax" type="string" value="1842.9" /> + <parameter name="tracker_region_zmax" type="string" value="2350" /> + <parameter name="world_box_hx" type="string" value="" /> + <parameter name="world_box_hy" type="string" value="" /> + <parameter name="world_box_hz" type="string" value="" /> + </detector> + <detector name="VXDInfra" geartype="GearParameters"> + <parameter name="ActiveLayerProperties_dEdx" type="double" value="3.870163611e-04" /> + <parameter name="BeSupportEndplateThickness" type="double" value="2.000000000e+00" /> + <parameter name="BeSupport_dEdx" type="double" value="2.941795296e-04" /> + <parameter name="CryostatAlHalfZ" type="double" value="1.766000000e+02" /> + <parameter name="CryostatAlInnerR" type="double" value="2.420000000e+01" /> + <parameter name="CryostatAlRadius" type="double" value="1.000000000e+02" /> + <parameter name="CryostatAlThickness" type="double" value="5.000000000e-01" /> + <parameter name="CryostatAlZEndCap" type="double" value="1.768500000e+02" /> + <parameter name="CryostatFoamRadius" type="double" value="9.000000000e+01" /> + <parameter name="CryostatFoamThickness" type="double" value="1.000000000e+01" /> + <parameter name="Cryostat_RadLen" type="double" value="8.896317758e+01" /> + <parameter name="Cryostat_dEdx" type="double" value="4.350185478e-04" /> + <parameter name="ElectronicEndLength" type="double" value="1.000000000e+01" /> + <parameter name="ElectronicEndThickness" type="double" value="1.000000000e-01" /> + <parameter name="StripLineBeamPipeRadius" type="double" value="2.430000000e+01" /> + <parameter name="VXDEndPlateInnerRadius" type="double" value="3.000000000e+01" /> + <parameter name="VXDSupport_dEdx" type="double" value="5.431907412e-05" /> + <parameter name="LadderGaps" type="DoubleVec" value="0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00" /> + <parameter name="StripLineFinalZ" type="DoubleVec" value="1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02" /> + </detector> + </detectors> + <materials> + <material name="VXDFoamShellMaterial" A="1.043890843e+01" Z="5.612886646e+00" density="2.500000000e+01" radLength="1.751650267e+04" intLength="6.594366018e+01" /> + <material name="VXDSupportMaterial" A="2.075865162e+01" Z="1.039383117e+01" density="2.765900000e+02" radLength="1.014262421e+03" intLength="1.206635688e+02" /> + </materials> +</gear> diff --git a/Digitisers/G2CDArbor/CMakeLists.txt b/Digitisers/G2CDArbor/CMakeLists.txt index 85949dbbe28e38d8c68c1868f327ebe30f260579..2bd2057d6ace2c94b993eee6637556d037b2f82f 100644 --- a/Digitisers/G2CDArbor/CMakeLists.txt +++ b/Digitisers/G2CDArbor/CMakeLists.txt @@ -3,12 +3,13 @@ gaudi_add_module(G2CDArbor SOURCES src/G2CDArborAlg.cpp LINK k4FWCore::k4FWCore GearSvc + DecoderHelperLib DetInterface Gaudi::GaudiKernel - Gaudi::GaudiAlgLib + Gaudi::GaudiAlgLib ${CLHEP_LIBRARIES} - ${GEAR_LIBRARIES} - ${GSL_LIBRARIES} + ${GEAR_LIBRARIES} + ${GSL_LIBRARIES} ${LCIO_LIBRARIES} EDM4HEP::edm4hep EDM4HEP::edm4hepDict DD4hep::DDRec diff --git a/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp b/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp index 2fbda2cc75388448ff365edff01e8a38019a183a..0cb54fada5d0fb6ce74769e9e386790d23e9542e 100644 --- a/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp +++ b/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp @@ -9,6 +9,8 @@ #include "DD4hep/IDDescriptor.h" #include "DD4hep/Plugins.h" +#include "DecoderHelper/DD4hep2Lcio.h" + #include <values.h> #include <string> @@ -33,23 +35,23 @@ DECLARE_COMPONENT( G2CDArborAlg ) // const double pi = acos(-1.0); struct DigiHit { - int digihitCellID0; + int digihitCellID0; int digihitCellID1; - float digihitCharge; - float digihitEnergyDepo; - int digihitNum1mmCell; + float digihitCharge; + float digihitEnergyDepo; + int digihitNum1mmCell; float PosX; float PosY; - float PosZ; - float LeadChargeDepo; - float ChargeShare; - edm4hep::SimCalorimeterHit * LeadSimCaloHit; + float PosZ; + float LeadChargeDepo; + float ChargeShare; + edm4hep::SimCalorimeterHit * LeadSimCaloHit; } ; //std::map <int, std::pair<float, float> >WeightVector; float ChanceOfKink = 0.1; -float KinkHitChargeBoost = 2.2; +float KinkHitChargeBoost = 2.2; // G2CDArborAlg aG2CDArborAlg ; @@ -138,7 +140,7 @@ G2CDArborAlg::G2CDArborAlg(const std::string& name, ISvcLocator* svcLoc) // _caloTruthLinkCollection, // caloTruthLinkCollection); - // std::vector<float> ChargeSpatialDistri; + // std::vector<float> ChargeSpatialDistri; // ChargeSpatialDistri.push_back(0.1); // ChargeSpatialDistri.push_back(0.2); // ChargeSpatialDistri.push_back(0.4); @@ -149,7 +151,7 @@ G2CDArborAlg::G2CDArborAlg(const std::string& name, ISvcLocator* svcLoc) // _ChargeSpatialDistri, // ChargeSpatialDistri); - // std::vector<float> ShowerPositionShiftID; + // std::vector<float> ShowerPositionShiftID; // ShowerPositionShiftID.push_back(0); // ShowerPositionShiftID.push_back(0); // ShowerPositionShiftID.push_back(0); @@ -186,7 +188,7 @@ G2CDArborAlg::G2CDArborAlg(const std::string& name, ISvcLocator* svcLoc) // registerProcessorParameter( "PolyaParaC" , // "Polya: x^a*exp(-b*x) + C" , // _PolyaParaC, - // float(0.03) ); + // float(0.03) ); // registerProcessorParameter( "ChanceOfKink" , // "Chance of one boundary hit to create a multiple hit with boosted charge" , @@ -317,7 +319,7 @@ StatusCode G2CDArborAlg::initialize() { else { ChanceOfKink = _ChanceOfKink; - KinkHitChargeBoost = _KinkHitChargeBoost; + KinkHitChargeBoost = _KinkHitChargeBoost; } float PolyaDomain = 100; @@ -341,10 +343,10 @@ StatusCode G2CDArborAlg::initialize() { int tmpIndex = 0; float NormalWeightFactor = 0; float NormalWeight[SizeCSD]; - int IndexA = 0; + int IndexA = 0; int IndexB = 0; //Used to denote charge distributions... - float WeightA = 0; - float WeightB = 0; + float WeightA = 0; + float WeightB = 0; for( int i0 = 0; i0 < SizeCSD; i0++ ) { @@ -370,37 +372,37 @@ StatusCode G2CDArborAlg::initialize() { WeightA = 0; WeightB = 0; - if( i2 < CoverAreaLength ) + if( i2 < CoverAreaLength ) { IndexA = CoverAreaLength - i2; - SignX = -1; + SignX = -1; } else if( i2 > _DigiCellSize - CoverAreaLength - 1) { IndexA = CoverAreaLength + i2 - _DigiCellSize + 1; - SignX = 1; + SignX = 1; } - if( j2 < CoverAreaLength ) + if( j2 < CoverAreaLength ) { IndexB = CoverAreaLength - j2; - SignY = -1; + SignY = -1; } else if( j2 > _DigiCellSize - CoverAreaLength - 1) { IndexB = CoverAreaLength + j2 - _DigiCellSize + 1; - SignY = 1; + SignY = 1; } - for(int i3 = 0; i3 < CoverAreaLength; i3 ++) + for(int i3 = 0; i3 < CoverAreaLength; i3 ++) { if(i3 < IndexA) WeightA += NormalWeight[i3]; if(i3 < IndexB) WeightB += NormalWeight[i3]; } - WMatrix.first = SignX*WeightA; - WMatrix.second = SignY*WeightB; - WeightVector[tmpIndex] = WMatrix; + WMatrix.first = SignX*WeightA; + WMatrix.second = SignY*WeightB; + WeightVector[tmpIndex] = WMatrix; cout<<WMatrix.first<<"/"<<WMatrix.second<<", "; } @@ -419,41 +421,41 @@ StatusCode G2CDArborAlg::execute() _eventNr++; float HitEn = 0; - float DigiHitEn = 0; + float DigiHitEn = 0; int LayerNum = 0; // TVector3 HitPos; //to be solved int tmpM, tmpS, tmpI, tmpJ, tmpK; - int CurrI = 0; - int CurrJ = 0; + int CurrI = 0; + int CurrJ = 0; int DHIndexI = 0; - int DHIndexJ = 0; + int DHIndexJ = 0; float DHChargeWeight = 0; - int DHCellID0 = 0; + int DHCellID0 = 0; float RndCharge = 0; - int SingleMCPPID = 0; - //float SingleMCPPEn = 0; + int SingleMCPPID = 0; + //float SingleMCPPEn = 0; float RefPosX = 0; float RefPosY = 0; float RefPosZ = 0; float DeltaPosI = 0; float DeltaPosJ = 0; - std::vector<float> CurrWeightVec; - float WeiI = 0; + std::vector<float> CurrWeightVec; + float WeiI = 0; float WeiJ = 0; int DeltaI = 0; - int DeltaJ = 0; - int MapI = 0; - int MapJ = 0; - int MapIndex = 0; + int DeltaJ = 0; + int MapI = 0; + int MapJ = 0; + int MapIndex = 0; float FHitPos[3] = {0, 0, 0}; float currTime = 0; float HitStepEn = 0; float EmaxStep = 0; // LCCollectionVec *relcol = new LCCollectionVec(LCIO::LCRELATION); - // LCFlagImpl linkflag; + // LCFlagImpl linkflag; // linkflag.setBit(LCIO::CHBIT_LONG); - // relcol->setFlag(linkflag.getFlag()); + // relcol->setFlag(linkflag.getFlag()); edm4hep::MCRecoCaloAssociationCollection* relcol = _caloTruthLinkCollection.createAndPut(); // LCCollection * MCPCol = evtP->getCollection("MCParticle"); @@ -464,7 +466,7 @@ StatusCode G2CDArborAlg::execute() // { // SingleMCPPID = a_MCP->getPDG(); // //SingleMCPPEn = a_MCP->getEnergy(); - // break; + // break; // } // } @@ -485,99 +487,142 @@ StatusCode G2CDArborAlg::execute() // LCCollectionVec *ecalcol = new LCCollectionVec(LCIO::CALORIMETERHIT); // ecalcol->setFlag(flag.getFlag()); // string EcalinitString = Ecalcol->getParameters().getStringVal(LCIO::CellIDEncoding); - // ecalcol->parameters().setValue(LCIO::CellIDEncoding, EcalinitString); + // ecalcol->parameters().setValue(LCIO::CellIDEncoding, EcalinitString); // for(int k1 = 0; k1 < NumEcalhit; k1++) - // { - if(m_readLCIO==false){ - std::string tmp_readout = m_col_readout_map[m_ecalColNames.value().at(k0)]; - // get the DD4hep readout - m_decoder = m_geosvc->getDecoder(tmp_readout); - if (!m_decoder) { - error() << "Failed to get the decoder. Skip this collection:"<<m_ecalColNames.value().at(k0)<< endmsg; - continue; - } - } - - edm4hep::CalorimeterHitCollection* ecalcol = _outputEcalCollections[k0]->createAndPut(); - auto Ecalcol = _ecalCollections[k0]->get(); - for (auto SimEcalhit: *Ecalcol){ - auto cellid = SimEcalhit.getCellID(); - // SimCalorimeterHit * SimEcalhit = dynamic_cast<SimCalorimeterHit*>( Ecalcol->getElementAt( k1 ) ) ; - // HitEn = SimEcalhit->getEnergy(); - // LayerNum = idDecoder(SimEcalhit)["K-1"]; - - // edm4hep::SimCalorimeterHit aa(SimEcalhit.getCellID(), SimEcalhit.getEnergy(), SimEcalhit.getPosition()); - ID_UTIL::CellIDDecoder< decltype(SimEcalhit) > cellIdDecoder(m_encoder_str); - const std::string layerCodingString(m_encoder_str); - const std::string layerCoding(this->GetLayerCoding(layerCodingString)); - if(m_readLCIO==false) LayerNum = m_decoder->get(cellid, "layer");//from 0 - 29, 0 is preshower - else LayerNum = cellIdDecoder(&SimEcalhit)[layerCoding.c_str()] + 1 ;//now it is 0 - 29, 0 is preshower - //cout << "LayerNum = " << LayerNum << endl; - - HitEn = SimEcalhit.getEnergy(); - unsigned long long cellID = SimEcalhit.getCellID(); - - currTime = 0; - EmaxStep = 0; - HitStepEn = 0; - // for(int k=0; k<SimEcalhit->getNMCContributions(); k++) - // { - for(int k=0; k<SimEcalhit.contributions_size(); k++){ - edm4hep::CaloHitContribution hitContribution = SimEcalhit.getContributions(k); - // HitStepEn = SimEcalhit->getEnergyCont(k); - HitStepEn = hitContribution.getEnergy(); - if(HitStepEn > EmaxStep) - { - EmaxStep = HitStepEn; - // currTime = SimEcalhit->getTimeCont(k); - currTime = hitContribution.getTime(); - } - } - - // _calibCoeffEcal[0] = 48.16; - // _calibCoeffEcal[1] = 96.32; - //if(LayerNum < _NEcalThinLayer) - if(LayerNum <= _NEcalThinLayer) //layer from 0 - 20 should be thin, total 21 thin layers, _NEcalThinLayer should be 20 - DigiHitEn = HitEn * _calibCoeffEcal[0]; - else - DigiHitEn = HitEn * _calibCoeffEcal[1]; - if( LayerNum==0) DigiHitEn = HitEn; // 0 is preshower layer - - totalEnergy += DigiHitEn; - if(HitEn > _thresholdEcal) - { - // CalorimeterHitImpl * DigiEcalhit = new CalorimeterHitImpl(); - // FHitPos[0] = SimEcalhit->getPosition()[0] + _ShowerPositionShiftID[0]*10.0; - // FHitPos[1] = SimEcalhit->getPosition()[1] + _ShowerPositionShiftID[1]*10.0; - // FHitPos[2] = SimEcalhit->getPosition()[2] + _ShowerPositionShiftID[2]*10.0; - - edm4hep::Vector3f hitPos = SimEcalhit.getPosition(); - FHitPos[0] = hitPos.x + _ShowerPositionShiftID[0]*10.0; - FHitPos[1] = hitPos.y + _ShowerPositionShiftID[1]*10.0; - FHitPos[2] = hitPos.z + _ShowerPositionShiftID[2]*10.0; - edm4hep::Vector3f FHitPosition(FHitPos[0], FHitPos[1], FHitPos[2]); - - // DigiEcalhit->setTime(currTime); - // DigiEcalhit->setPosition( FHitPos ); - // DigiEcalhit->setCellID0(SimEcalhit->getCellID0()); - // DigiEcalhit->setCellID1(SimEcalhit->getCellID1()); - // DigiEcalhit->setEnergy(DigiHitEn); - // ecalcol->addElement(DigiEcalhit); - auto DigiEcalhit = ecalcol->create(); - DigiEcalhit.setTime(currTime); - DigiEcalhit.setPosition(FHitPosition); - DigiEcalhit.setCellID(cellID); - DigiEcalhit.setEnergy(DigiHitEn); - - // LCRelationImpl *rel = new LCRelationImpl(DigiEcalhit, SimEcalhit, 1.0); //only keep the leading contribution - // relcol->addElement(rel); - auto rel = relcol->create(); - rel.setRec(DigiEcalhit); - rel.setSim(SimEcalhit); - rel.setWeight(1.0); - } - } + // { + if(m_readLCIO==false){ + std::string tmp_readout = m_col_readout_map[m_ecalColNames.value().at(k0)]; + // get the DD4hep readout + m_decoder = m_geosvc->getDecoder(tmp_readout); + if (!m_decoder) { + error() << "Failed to get the decoder. Skip this collection:"<<m_ecalColNames.value().at(k0)<< endmsg; + continue; + } + } + + edm4hep::CalorimeterHitCollection* ecalcol = _outputEcalCollections[k0]->createAndPut(); + auto Ecalcol = _ecalCollections[k0]->get(); + + int iHit = 0; + double ESum_ECALSimu = 0; + double ESum_ECALDigi = 0; + cout<<"[YX debug - G2CD] ECAL Collection Simu "<<k0<<", nSimuHit = "<<Ecalcol->size()<<", "<<m_ecalColNames.value().at(k0)<<endl; + + for (auto SimEcalhit: *Ecalcol){ + auto cellid = SimEcalhit.getCellID(); + // SimCalorimeterHit * SimEcalhit = dynamic_cast<SimCalorimeterHit*>( Ecalcol->getElementAt( k1 ) ) ; + // HitEn = SimEcalhit->getEnergy(); + // LayerNum = idDecoder(SimEcalhit)["K-1"]; + + // edm4hep::SimCalorimeterHit aa(SimEcalhit.getCellID(), SimEcalhit.getEnergy(), SimEcalhit.getPosition()); + + ID_UTIL::CellIDDecoder< decltype(SimEcalhit) > cellIdDecoder(m_encoder_str); + const std::string layerCodingString(m_encoder_str); + const std::string layerCoding(this->GetLayerCoding(layerCodingString)); + if(m_readLCIO){ + // LayerNum = cellIdDecoder(&SimEcalhit)[layerCoding.c_str()] + 1 ;//now it is 0 - 29, 0 is preshower + // actually 1-30 + + LayerNum = cellIdDecoder(&SimEcalhit)[layerCoding.c_str()]; // 0 - 29 + + }else{ + // LayerNum = m_decoder->get(cellid, "layer"); //from 0 - 29, 0 is preshower + + int Raw_Layer = m_decoder->get(cellid, "layer"); // 0 - 29 + LayerNum = DD4hep2Lcio::CEPCv4::getEcalLayer(Raw_Layer); // -1 - 28, -1 is preshower + } + + + HitEn = SimEcalhit.getEnergy(); + unsigned long long cellID = SimEcalhit.getCellID(); + + currTime = 0; + EmaxStep = 0; + HitStepEn = 0; + // for(int k=0; k<SimEcalhit->getNMCContributions(); k++) + // { + for(int k=0; k<SimEcalhit.contributions_size(); k++){ + edm4hep::CaloHitContribution hitContribution = SimEcalhit.getContributions(k); + // HitStepEn = SimEcalhit->getEnergyCont(k); + HitStepEn = hitContribution.getEnergy(); + if(HitStepEn > EmaxStep) + { + EmaxStep = HitStepEn; + // currTime = SimEcalhit->getTimeCont(k); + currTime = hitContribution.getTime(); + } + } + + // _calibCoeffEcal[0] = 48.16; + // _calibCoeffEcal[1] = 96.32; + // if(LayerNum <= _NEcalThinLayer) //layer from 0 - 20 should be thin, total 21 thin layers, _NEcalThinLayer should be 20 + if(LayerNum < _NEcalThinLayer) // _NEcalThinLayer = 20 + DigiHitEn = HitEn * _calibCoeffEcal[0]; + else + DigiHitEn = HitEn * _calibCoeffEcal[1]; + + if(LayerNum==-1) DigiHitEn = HitEn; // -1 is preshower layer + + totalEnergy += DigiHitEn; + + // TVector3 HitPos = SimEcalhit.getPosition(); + edm4hep::Vector3f HitPos = SimEcalhit.getPosition(); + double HitDis = sqrt(HitPos.x*HitPos.x + HitPos.y*HitPos.y + HitPos.z*HitPos.z); + double HitTime = currTime - HitDis / 300; + + + if(0){ + cout<<"Hit Pos(x,y,z,perp) = ("<<HitPos.x<<", "<<HitPos.y<<", "<<HitPos.z<<", "<<sqrt(HitPos.x*HitPos.x+HitPos.y*HitPos.y)<<")"<<endl; + cout<<"---> LayerNum = "<<LayerNum<<endl; + } + + // if(HitEn > _thresholdEcal) + if(HitEn > _thresholdEcal && HitTime < _TimeThreshold) + { + // CalorimeterHitImpl * DigiEcalhit = new CalorimeterHitImpl(); + // FHitPos[0] = SimEcalhit->getPosition()[0] + _ShowerPositionShiftID[0]*10.0; + // FHitPos[1] = SimEcalhit->getPosition()[1] + _ShowerPositionShiftID[1]*10.0; + // FHitPos[2] = SimEcalhit->getPosition()[2] + _ShowerPositionShiftID[2]*10.0; + + edm4hep::Vector3f hitPos = SimEcalhit.getPosition(); + FHitPos[0] = hitPos.x + _ShowerPositionShiftID[0]*10.0; + FHitPos[1] = hitPos.y + _ShowerPositionShiftID[1]*10.0; + FHitPos[2] = hitPos.z + _ShowerPositionShiftID[2]*10.0; + edm4hep::Vector3f FHitPosition(FHitPos[0], FHitPos[1], FHitPos[2]); + + // DigiEcalhit->setTime(currTime); + // DigiEcalhit->setPosition( FHitPos ); + // DigiEcalhit->setCellID0(SimEcalhit->getCellID0()); + // DigiEcalhit->setCellID1(SimEcalhit->getCellID1()); + // DigiEcalhit->setEnergy(DigiHitEn); + // ecalcol->addElement(DigiEcalhit); + auto DigiEcalhit = ecalcol->create(); + DigiEcalhit.setTime(currTime); + DigiEcalhit.setPosition(FHitPosition); + DigiEcalhit.setCellID(cellID); + DigiEcalhit.setEnergy(DigiHitEn); + + // LCRelationImpl *rel = new LCRelationImpl(DigiEcalhit, SimEcalhit, 1.0); //only keep the leading contribution + // relcol->addElement(rel); + auto rel = relcol->create(); + rel.setRec(DigiEcalhit); + rel.setSim(SimEcalhit); + rel.setWeight(1.0); + + ESum_ECALSimu += SimEcalhit.getEnergy(); + ESum_ECALDigi += DigiEcalhit.getEnergy(); + // cout<<"[YX debug - G2CD] ECAL DigiHit "<<iHit<<", E = "<<DigiEcalhit.getEnergy()<<", T = "<<DigiEcalhit.getTime()<<", Pos ("<<FHitPos[0]<<", "<<FHitPos[1]<<", "<<FHitPos[2]<<")"<<endl; + + } + + iHit+=1; + } + + + int nDigiHit_ECAL = ecalcol->size(); + cout<<"[YX debug - G2CD] ECAL Collection Digi "<<k0<<", nDigiHit = "<<nDigiHit_ECAL<<", "<<m_ecalOutputColNames.value().at(k0)<<", DigiESum= "<<ESum_ECALSimu<<", SimuESum = "<<ESum_ECALDigi<<endl; + // evtP->addCollection(ecalcol,_outputEcalCollections[k0].c_str()); @@ -598,7 +643,7 @@ StatusCode G2CDArborAlg::execute() // //if(k=="CEPC_v1" || k=="ild_o2_v05") // 1cm cell simulation... // { - // cout<<"Detector Tpype: "<<k<<endl; + // cout<<"Detector Tpype: "<<k<<endl; for (unsigned int k3 = 0; k3 < _hcalCollections.size(); ++k3) { @@ -609,12 +654,20 @@ StatusCode G2CDArborAlg::execute() // hcalcol->setFlag(flag.getFlag()); // string HcalinitString = Hcalcol->getParameters().getStringVal(LCIO::CellIDEncoding); // hcalcol->parameters().setValue(LCIO::CellIDEncoding, HcalinitString); - + // for(int k4 = 0; k4 < NumHcalhit; k4++) // { + + edm4hep::CalorimeterHitCollection* hcalcol = _outputHcalCollections[k3]->createAndPut(); auto Hcalcol = _hcalCollections[k3]->get(); + + int iHit = 0; + double ESum_HCALSimu = 0; + double ESum_HCALDigi = 0; + cout<<"[YX debug - G2CD] HCAL Collection Simu "<<k3<<", nSimuHit = "<<Hcalcol->size()<<", "<<m_hcalColNames.value().at(k3)<<endl; + for(auto SimHcalhit: *Hcalcol){ // SimCalorimeterHit * SimHcalhit = dynamic_cast<SimCalorimeterHit*>( Hcalcol->getElementAt( k4 ) ) ; @@ -638,8 +691,12 @@ StatusCode G2CDArborAlg::execute() } } - // if(SimHcalhit->getEnergy() > 0) //some threshold can be added. - if(SimHcalhit.getEnergy() > 0) //some threshold can be added. + edm4hep::Vector3f HitPos = SimHcalhit.getPosition(); + double HitDis = sqrt(HitPos.x*HitPos.x + HitPos.y*HitPos.y + HitPos.z*HitPos.z); + double HitTime = currTime - HitDis / 300; + + // if(SimHcalhit.getEnergy() > 0) //some threshold can be added. + if(SimHcalhit.getEnergy() > 0 && HitTime < _TimeThreshold) //some threshold can be added. { // CalorimeterHitImpl * DigiHcalhit = new CalorimeterHitImpl(); @@ -672,9 +729,23 @@ StatusCode G2CDArborAlg::execute() rel.setRec(DigiHcalhit); rel.setSim(SimHcalhit); rel.setWeight(1.0); + + + ESum_HCALSimu += SimHcalhit.getEnergy(); + ESum_HCALDigi += DigiHcalhit.getEnergy(); + // cout<<"[YX debug - G2CD] HCAL DigiHit "<<iHit<<", E = "<<DigiHcalhit.getEnergy()<<", T = "<<DigiHcalhit.getTime()<<", Pos ("<<FHitPos[0]<<", "<<FHitPos[1]<<", "<<FHitPos[2]<<")"<<endl; + } + + iHit+=1; + } + int nDigiHit_HCAL = hcalcol->size(); + cout<<"[YX debug - G2CD] HCAL Collection Digi "<<k3<<", nDigiHit = "<<nDigiHit_HCAL<<", "<<m_hcalOutputColNames.value().at(k3)<<", DigiESum= "<<ESum_HCALSimu<<", SimuESum = "<<ESum_HCALDigi<<endl; + + + // evtP->addCollection(hcalcol,_outputHcalCollections[k3].c_str()); // }catch(lcio::DataNotAvailableException zero) { } @@ -685,7 +756,7 @@ StatusCode G2CDArborAlg::execute() // string EcalPSinitString = "M:3,S-1:3,I:9,J:9,K-1:6"; // ecalPScol->parameters().setValue(LCIO::CellIDEncoding, EcalPSinitString); - // // parameter CellIDEncoding [string]: M:3,S-1:3,I:9,J:9,K-1:6, + // // parameter CellIDEncoding [string]: M:3,S-1:3,I:9,J:9,K-1:6, // for(unsigned int s0 = 0; s0 < _EcalPreShowerCollections.size(); s0++) //I think should digitize and give a very small energy; even a new collection called PS Hits // { @@ -708,7 +779,7 @@ StatusCode G2CDArborAlg::execute() // EmaxStep = HitStepEn; // currTime = a_hit->getTimeCont(k); // } - // } + // } // CalorimeterHitImpl * PShit = new CalorimeterHitImpl(); @@ -732,16 +803,16 @@ StatusCode G2CDArborAlg::execute() // } // else // { - // for (unsigned int i(0); i < _hcalCollections.size(); ++i) + // for (unsigned int i(0); i < _hcalCollections.size(); ++i) // { //strictly follow barrel, endcap, ring order - // std::map <int, DigiHit> IDtoDigiHit; + // std::map <int, DigiHit> IDtoDigiHit; // IDtoDigiHit.clear(); - // std::map <int, TVector3> IDtoPos; + // std::map <int, TVector3> IDtoPos; // IDtoPos.clear(); // std::map <int, float> IDtoTime; // IDtoTime.clear(); - + // try{ // LCCollection * col = evtP->getCollection( _hcalCollections[i].c_str() ) ; @@ -782,12 +853,12 @@ StatusCode G2CDArborAlg::execute() // RefPosY = hit->getPosition()[1]; // RefPosZ = hit->getPosition()[2]; - // CurrI = tmpI/_DigiCellSize; + // CurrI = tmpI/_DigiCellSize; // CurrJ = tmpJ/_DigiCellSize; // MapI = tmpI % _DigiCellSize; // MapJ = tmpJ % _DigiCellSize; - // MapIndex = MapI * _DigiCellSize + MapJ; - // WeiI = WeightVector[MapIndex].first; + // MapIndex = MapI * _DigiCellSize + MapJ; + // WeiI = WeightVector[MapIndex].first; // WeiJ = WeightVector[MapIndex].second; // DeltaPosI = (float(_DigiCellSize) - 1.0)/2 - MapI; @@ -795,11 +866,11 @@ StatusCode G2CDArborAlg::execute() // // cout<<"DeltaI "<<DeltaPosI<<", "<<DeltaPosJ<<endl; - // DeltaI = 0; - // DeltaJ = 0; + // DeltaI = 0; + // DeltaJ = 0; - // if(fabs(WeiI) > 1E-9) - // { + // if(fabs(WeiI) > 1E-9) + // { // DeltaI = ((WeiI > 0) ? 1: -1); // } // if(fabs(WeiJ) > 1E-9) @@ -840,7 +911,7 @@ StatusCode G2CDArborAlg::execute() // if(DeltaI) // { // DHIndexI = CurrI + DeltaI; - // DHIndexJ = CurrJ; + // DHIndexJ = CurrJ; // } // if(DeltaJ) // { @@ -876,9 +947,9 @@ StatusCode G2CDArborAlg::execute() // { // IDtoDigiHit[DHCellID0].PosX = RefPosX + (DeltaPosI + int(DHIndexI - CurrI)*_DigiCellSize) * cos(tmpS*pi/4.0) + _ShowerPositionShiftID[0]*_DigiCellSize; //(mm in unit) // IDtoDigiHit[DHCellID0].PosY = RefPosY + (DeltaPosI + int(DHIndexI - CurrI)*_DigiCellSize)* sin(tmpS*pi/4.0) + _ShowerPositionShiftID[1]*_DigiCellSize; - // IDtoDigiHit[DHCellID0].PosZ = RefPosZ + DeltaPosJ + int(DHIndexJ - CurrJ)*_DigiCellSize + _ShowerPositionShiftID[2]*_DigiCellSize; //Rotation is needed, based on S; + // IDtoDigiHit[DHCellID0].PosZ = RefPosZ + DeltaPosJ + int(DHIndexJ - CurrJ)*_DigiCellSize + _ShowerPositionShiftID[2]*_DigiCellSize; //Rotation is needed, based on S; // } - // else //endcap or ring; + // else //endcap or ring; // { // IDtoDigiHit[DHCellID0].PosX = RefPosX + DeltaPosI + (DHIndexI - CurrI)*_DigiCellSize + _ShowerPositionShiftID[0]*_DigiCellSize; //(mm in unit) // IDtoDigiHit[DHCellID0].PosY = RefPosY + DeltaPosJ + (DHIndexJ - CurrJ)*_DigiCellSize + _ShowerPositionShiftID[1]*_DigiCellSize; @@ -888,7 +959,7 @@ StatusCode G2CDArborAlg::execute() // HitPos.SetXYZ(IDtoDigiHit[DHCellID0].PosX, IDtoDigiHit[DHCellID0].PosY, IDtoDigiHit[DHCellID0].PosZ); // IDtoPos[DHCellID0] = HitPos; - // IDtoTime[DHCellID0] = currTime; + // IDtoTime[DHCellID0] = currTime; // } // else // { @@ -898,7 +969,7 @@ StatusCode G2CDArborAlg::execute() // if(RndCharge * DHChargeWeight > IDtoDigiHit[DHCellID0].LeadChargeDepo) // { // IDtoDigiHit[DHCellID0].LeadChargeDepo = RndCharge * DHChargeWeight; - // IDtoDigiHit[DHCellID0].LeadSimCaloHit = hit; + // IDtoDigiHit[DHCellID0].LeadSimCaloHit = hit; // IDtoDigiHit[DHCellID0].ChargeShare = DHChargeWeight; // } @@ -917,14 +988,14 @@ StatusCode G2CDArborAlg::execute() // CalorimeterHitImpl * calhit = new CalorimeterHitImpl(); // LCRelationImpl *rel = new LCRelationImpl(calhit, ff->second.LeadSimCaloHit, ff->second.ChargeShare); //only keep the leading contribution // relcol->addElement(rel); - + // calhit->setCellID0( ff->first ); //Assume 100% efficiency // //calhit->setEnergy( DHCALCalibrationConstant ); //Charge // calhit->setEnergy( _thresholdHcal ); //Charge // calhit->setCellID1(SingleMCPPID); //Use ID1 & Energy Error to denote the MCP info... // calhit->setEnergyError(ff->second.digihitCharge); // (SingleMCPPEn); // /* - // DigiHitPos[0] = ff->second.PosX; + // DigiHitPos[0] = ff->second.PosX; // DigiHitPos[1] = ff->second.PosY; // DigiHitPos[2] = ff->second.PosZ; // */ @@ -933,7 +1004,7 @@ StatusCode G2CDArborAlg::execute() // DigiHitPos[2] = IDtoPos[ff->first].Z(); // calhit->setTime(IDtoTime[ff->first]); - // calhit->setPosition(DigiHitPos); + // calhit->setPosition(DigiHitPos); // hcalcol->addElement(calhit); // } diff --git a/Digitisers/G2CDArbor/src/G2CDArborAlg.h b/Digitisers/G2CDArbor/src/G2CDArborAlg.h index 08c04f0335266d3ecfbd035c0b30c6c798509161..45e1d1cfef34f6789f4044e2cf74144853d03a7b 100644 --- a/Digitisers/G2CDArbor/src/G2CDArborAlg.h +++ b/Digitisers/G2CDArbor/src/G2CDArborAlg.h @@ -38,7 +38,7 @@ public: /* G2CDArborAlg(); */ /* ~G2CDArborAlg() {}; */ - /** Called at the begin of the job before anything is read. + /** Called at the begin of the job before anything is read. * Use to initialize the processor, e.g. book histograms. */ virtual StatusCode initialize() ; @@ -98,6 +98,7 @@ protected: Gaudi::Property<int> _NEcalThinLayer{this, "NumThinEcalLayer", 20, "Num of thiner Ecal layers"}; Gaudi::Property<float> _thresholdEcal{this, "ECALThreshold", (float)5.0e-5, "Threshold for ECAL Hits in GeV"}; Gaudi::Property<float> _thresholdHcal{this, "HCALThreshold", (float)0.11, "Threshold for HCAL Hits in GeV"}; + Gaudi::Property<float> _TimeThreshold{this, "TimeThreshold", (float)1000, "Time Threshold for both ECAL and HCAL Hits in ns"}; Gaudi::Property<int> _DigiCellSize{this, "DigiCellSize", 10, "Size of Digitized Cell (in mm)"}; Gaudi::Property<float> _ShiftInX{this, "ShiftInX", (float)0.0, "Shift Distance in X directoin (in mm) NP only"}; Gaudi::Property<int> _UsingDefaultDetector{this, "UsingDefaultDetector", 0, "Flag Parameter Setting (0 ~ self definition, 1 ~ MircoMegas, 2 ~ GRPC_PS, 3 ~ GRPC_SPS)"}; @@ -133,18 +134,18 @@ protected: /* float _PolyaParaA, _PolyaParaB, _PolyaParaC; */ /* float _ChanceOfKink, _KinkHitChargeBoost; */ TTree *_outputTree; - TH1F *_NH1stLayer, *_NH8thLayer; - TF1 * _QPolya; + TH1F *_NH1stLayer, *_NH8thLayer; + TF1 * _QPolya; int _Num; - int _eventNr; + int _eventNr; int _M, _S, _I, _J, _K, _Seg; - float _PosX, _PosY, _PosZ; + float _PosX, _PosY, _PosZ; /* float _EDepo, _Charge, _ShiftInX; */ int _NHit1mm, _NHit1mmCenter, _NHit1mmCorner, _NHit1mmSide; - int _TotalNHit1mm, _TotalNHit, _TotalMultiHit; - int _N1, _N2, _N3; + int _TotalNHit1mm, _TotalNHit, _TotalMultiHit; + int _N1, _N2, _N3; std::string _fileName; std::ostream *_output; diff --git a/Examples/options/CEPCV4_simu_reco_Arbor.py b/Examples/options/CEPCV4_simu_reco_Arbor.py new file mode 100644 index 0000000000000000000000000000000000000000..dbd5e7212b64f3b4851cb3f99c490b7d8501a87f --- /dev/null +++ b/Examples/options/CEPCV4_simu_reco_Arbor.py @@ -0,0 +1,313 @@ +#!/usr/bin/env python + +import os + +from Gaudi.Configuration import * + +# NTupleSvc().Output = ["MyTuples DATAFILE='Examples/options/TPC_out.root' OPT='NEW' TYP='ROOT'"] + +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ +rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() +rndmengine.SetSingleton = True +rndmengine.Seeds = [1] + +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc") + +from Configurables import MarlinEvtSeeder +evtseeder = MarlinEvtSeeder("EventSeeder") + +geometry_option = "CepC_v4.xml" + +if not os.getenv("DETCEPCV4ROOT"): + print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import GeomSvc +geosvc = GeomSvc("GeomSvc") +geosvc.compact = geometry_path + +from Configurables import GearSvc +gearsvc = GearSvc("GearSvc") +#gearsvc.GearXMLFile = "Detector/DetCEPCv4/compact/FullDetGear.xml" + +############################################################################## +# Generator +from Configurables import GenAlgo +from Configurables import GtGunTool +from Configurables import StdHepRdr +from Configurables import SLCIORdr +from Configurables import HepMCRdr +from Configurables import GenPrinter + +gun = GtGunTool("GtGunTool") +gun.Particles = ["e-"] +#gun.EnergyMins = [1] +#gun.EnergyMaxs = [50] +#gun.ThetaMins = [50] +#gun.ThetaMaxs = [130] +#gun.PhiMins = [-90] +#gun.PhiMaxs = [90] +gun.EnergyMins = [10] +gun.EnergyMaxs = [10] +gun.ThetaMins = [90] +gun.ThetaMaxs = [90] +gun.PhiMins = [0] +gun.PhiMaxs = [0] + +stdheprdr = StdHepRdr("StdHepRdr") +stdheprdr.Input = "/cefs/data/stdhep/CEPC240/higgs/Higgs_10M/data/E240.Pnnh_gg.e0.p0.whizard195/nnh_gg.e0.p0.00001.stdhep" +# stdheprdr.Input = "GenFile" + +genprinter = GenPrinter("GenPrinter") + +genalg = GenAlgo("GenAlgo") +# genalg.GenTools = ["GtGunTool"] +genalg.GenTools = ["StdHepRdr"] + +############################################################################## +# Detector simulation +from Configurables import DetSimSvc +detsimsvc = DetSimSvc("DetSimSvc") + +from Configurables import DetSimAlg +detsimalg = DetSimAlg("DetSimAlg") +# detsimalg.VisMacs = ["vis.mac"] +detsimalg.RunCmds = [ +# "/physics_lists/factory/addOptical" +] +detsimalg.PhysicsList = "FTFP_BERT" +detsimalg.AnaElems = ["Edm4hepWriterAnaElemTool"] +detsimalg.RootDetElem = "WorldDetElemTool" + +from Configurables import AnExampleDetElemTool +example_dettool = AnExampleDetElemTool("AnExampleDetElemTool") + +############################################################################## +# Tracker + +from Configurables import TimeProjectionChamberSensDetTool +tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool") +tpc_sensdettool.TypeOption = 1 + +from Configurables import TrackSystemSvc +tracksystemsvc = TrackSystemSvc("TrackSystemSvc") + +vxdhitname = "VXDTrackerHits" +sithitname = "SITTrackerHits" +sitspname = "SITSpacePoints" +tpchitname = "TPCTrackerHits" +sethitname = "SETTrackerHits" +setspname = "SETSpacePoints" +ftdspname = "FTDSpacePoints" +ftdhitname = "FTDTrackerHits" +from Configurables import PlanarDigiAlg +digiVXD = PlanarDigiAlg("VXDDigi") +digiVXD.SimTrackHitCollection = "VXDCollection" +digiVXD.TrackerHitCollection = vxdhitname +digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation" +digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] +digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] + +digiSIT = PlanarDigiAlg("SITDigi") +digiSIT.IsStrip = 1 +digiSIT.SimTrackHitCollection = "SITCollection" +digiSIT.TrackerHitCollection = sithitname +digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" +digiSIT.ResolutionU = [0.007] +digiSIT.ResolutionV = [0.000] + +digiSET = PlanarDigiAlg("SETDigi") +digiSET.IsStrip = 1 +digiSET.SimTrackHitCollection = "SETCollection" +digiSET.TrackerHitCollection = sethitname +digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation" +digiSET.ResolutionU = [0.007] +digiSET.ResolutionV = [0.000] + +digiFTD = PlanarDigiAlg("FTDDigi") +digiFTD.SimTrackHitCollection = "FTDCollection" +digiFTD.TrackerHitCollection = ftdhitname +digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +digiFTD.ResolutionU = [0.003, 0.003, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007] +digiFTD.ResolutionV = [0.003, 0.003, 0, 0, 0, 0, 0, 0 ] +#digiFTD.OutputLevel = DEBUG + +from Configurables import SpacePointBuilderAlg +spSIT = SpacePointBuilderAlg("SITBuilder") +spSIT.TrackerHitCollection = sithitname +spSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" +spSIT.SpacePointCollection = sitspname +spSIT.SpacePointAssociationCollection = "SITSpacePointAssociation" +#spSIT.OutputLevel = DEBUG + +spFTD = SpacePointBuilderAlg("FTDBuilder") +spFTD.TrackerHitCollection = ftdhitname +spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +spFTD.SpacePointCollection = ftdspname +spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation" +#spFTD.OutputLevel = DEBUG + +from Configurables import TPCDigiAlg +digiTPC = TPCDigiAlg("TPCDigi") +digiTPC.TPCCollection = "TPCCollection" +digiTPC.TPCLowPtCollection = "TPCLowPtCollection" +digiTPC.TPCTrackerHitsCol = tpchitname +digiTPC.TPCTrackerHitAssCol = "TPCTrackerHitAssociation" +#digiTPC.OutputLevel = DEBUG + +from Configurables import ClupatraAlg +clupatra = ClupatraAlg("Clupatra") +clupatra.TPCHitCollection = tpchitname +#clupatra.OutputLevel = DEBUG + +from Configurables import SiliconTrackingAlg +tracking = SiliconTrackingAlg("SiliconTracking") +tracking.HeaderCol = "EventHeader" +tracking.VTXHitCollection = vxdhitname +tracking.SITHitCollection = sitspname +tracking.FTDPixelHitCollection = ftdhitname +tracking.FTDSpacePointCollection = ftdspname +tracking.SITRawHitCollection = sithitname +tracking.FTDRawHitCollection = ftdhitname +tracking.UseSIT = 1 +tracking.SmoothOn = 0 +#tracking.OutputLevel = DEBUG + +from Configurables import ForwardTrackingAlg +forward = ForwardTrackingAlg("ForwardTracking") +forward.FTDPixelHitCollection = ftdhitname +forward.FTDSpacePointCollection = ftdspname +forward.FTDRawHitCollection = ftdhitname +forward.Chi2ProbCut = 0.0 +forward.HitsPerTrackMin = 3 +forward.BestSubsetFinder = "SubsetSimple" +forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation", + "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"] +forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0] +forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05] +#forward.OutputLevel = DEBUG + +from Configurables import TrackSubsetAlg +subset = TrackSubsetAlg("TrackSubset") +subset.TrackInputCollections = ["ForwardTracks", "SiTracks"] +subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, sitspname, ftdspname] +subset.TrackSubsetCollection = "SubsetTracks" +#subset.OutputLevel = DEBUG + +from Configurables import FullLDCTrackingAlg +full = FullLDCTrackingAlg("FullTracking") +full.VTXTrackerHits = vxdhitname +full.SITTrackerHits = sitspname +full.TPCTrackerHits = tpchitname +full.SETTrackerHits = setspname +full.FTDPixelTrackerHits = ftdhitname +full.FTDSpacePoints = ftdspname +full.SITRawHits = sithitname +full.SETRawHits = sethitname +full.FTDRawHits = ftdhitname +full.TPCTracks = "ClupatraTracks" +full.SiTracks = "SubsetTracks" +full.OutputTracks = "MarlinTrkTracks" +#full.OutputLevel = DEBUG +''' +from Configurables import DumpMCParticleAlg +dumpMC = DumpMCParticleAlg("DumpMC") +dumpMC.MCParticleCollection = "MCParticle" + +from Configurables import DumpTrackAlg +dumpFu = DumpTrackAlg("DumpFu") +dumpFu.TrackCollection = "MarlinTrkTracks" +#dumpFu.OutputLevel = DEBUG + +dumpCl = DumpTrackAlg("DumpCl") +dumpCl.TrackCollection = "ClupatraTracks" +#dumpCl.OutputLevel = DEBUG + +dumpSu = DumpTrackAlg("DumpSu") +dumpSu.TrackCollection = "SubsetTracks" +#dumpSu.OutputLevel = DEBUG + +dumpSi = DumpTrackAlg("DumpSi") +dumpSi.TrackCollection = "SiTracks" +#dumpSi.OutputLevel = DEBUG + +dumpFo = DumpTrackAlg("DumpFo") +dumpFo.TrackCollection = "ForwardTracks" +#dumpFo.OutputLevel = DEBUG +''' +############################################################################## +# Calorimeter + +from Configurables import SimHitMergeAlg +simHitMerge = SimHitMergeAlg("SimHitMergeAlg") +simHitMerge.sanity_check = False +simHitMerge.InputCollections=["EcalBarrelCollection", "EcalEndcapsCollection","EcalEndcapRingCollection", "HcalBarrelCollection", "HcalEndcapsCollection", "HcalEndcapRingCollection"] +simHitMerge.OutputCollections=["EcalBarrelCollectionMerged", "EcalEndcapsCollectionMerged", "EcalEndcapRingCollectionMerged", "HcalBarrelCollectionMerged", "HcalEndcapsCollectionMerged", "HcalEndcapRingCollectionMerged"] + +# simHitMerge.InputCollections=["EcalBarrelCollection", "EcalEndcapsCollection"] +# simHitMerge.OutputCollections=["EcalBarrelCollectionMerged", "EcalEndcapsCollectionMerged"] +############################################################################## +from Configurables import G2CDArborAlg +caloDigi = G2CDArborAlg("G2CDArborAlg") +caloDigi.ReadLCIO = False +caloDigi.ECALCollections = ["EcalBarrelCollectionMerged", "EcalEndcapsCollectionMerged", "EcalEndcapRingCollectionMerged"] +caloDigi.HCALCollections = ["HcalBarrelCollectionMerged", "HcalEndcapsCollectionMerged", "HcalEndcapRingCollectionMerged"] +caloDigi.ECALReadOutNames = ["EcalBarrelCollection", "EcalEndcapsCollection", "EcalEndcapRingCollection"] +caloDigi.HCALReadOutNames = ["HcalBarrelCollection", "HcalEndcapsCollection", "HcalEndcapRingCollection"] +caloDigi.DigiECALCollection = ["ECALBarrel", "ECALEndcap", "ECALOther"] +caloDigi.DigiHCALCollection = ["HCALBarrel", "HCALEndcap", "HCALOther"] +caloDigi.EventReportEvery = 100 +# caloDigi.CalibrECAL = [46.538, 93.0769] # Yudan +caloDigi.CalibrECAL = [48.16, 96.32] +caloDigi.HCALThreshold = 0.12 + +caloDigi.PolyaParaA = 1.1 +caloDigi.PolyaParaB = 1.0 +caloDigi.PolyaParaC = 0.0 +############################################################################## +# Reconstruction: Arbor + +from Configurables import MarlinArbor +marlinArbor = MarlinArbor("MarlinArbor") +marlinArbor.ReadLCIO = False +marlinArbor.ECALCollections =["ECALBarrel","ECALEndcap","ECALOther"] +marlinArbor.HCALCollections =["HCALBarrel","HCALEndcap","HCALOther"] +marlinArbor.ECALReadOutNames = ["EcalBarrelCollection", "EcalEndcapsCollection", "EcalEndcapRingCollection"] +marlinArbor.HCALReadOutNames = ["HcalBarrelCollection", "HcalEndcapsCollection", "HcalEndcapRingCollection"] +############################################################################## +from Configurables import BushConnect +bushconnect = BushConnect("BushConnect") +bushconnect.ReadLCIO = False +############################################################################## +# BMR analysis +from Configurables import TotalInvMass +totalInvM = TotalInvMass("TotalInvMass") +totalInvM.TreeOutputFile = "Examples/options/CEPCV4_simu_reco_Arbor_nnHgg_BMRAna.root" +# totalInvM.TreeOutputFile = "BMRFile" +totalInvM.MCPCollectionName = "MCParticleGen" +# totalInvM.MCPCollectionName = "MCParticle" +############################################################################## +# write PODIO file +from Configurables import PodioOutput +write = PodioOutput("write") +write.filename = "Examples/options/CEPCV4_simu_reco_Arbor_nnHgg_Reco.root" +# write.filename = "RecoFile" +write.outputCommands = ["keep *"] +############################################################################## +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( + # TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSIT, spFTD, digiTPC, clupatra, tracking, forward, subset, full, simHitMerge, caloDigi, write], + TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSIT, spFTD, digiTPC, clupatra, tracking, forward, subset, full, simHitMerge, caloDigi, marlinArbor, bushconnect, totalInvM, write], + EvtSel = 'NONE', + EvtMax = 1, + ExtSvc = [rndmengine, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc], + HistogramPersistency='ROOT', + OutputLevel=INFO +) diff --git a/Examples/options/LCIO_read_Arbor.py b/Examples/options/LCIO_read_Arbor.py new file mode 100644 index 0000000000000000000000000000000000000000..df01a342338a7e5b7411cc2a490e958e9ed648c1 --- /dev/null +++ b/Examples/options/LCIO_read_Arbor.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python + +import os +from Gaudi.Configuration import * +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc") +######################################################################### +# read LCIO files +from Configurables import LCIOInput +read = LCIOInput("read") +read.inputs = [ +"/cefs/higgs/PFAData/Sample/Generator/FullGeo_Baseline/E240_nnH_gg/Reco/Reco_Baseline_E240_nnH_gg_00001_0.slcio" +] +read.mode = 1 +read.collections = [ + "MCParticle:MCParticle", + "SimCalorimeterHit:EcalBarrelSiliconCollection", + "SimCalorimeterHit:EcalBarrelSiliconPreShowerCollection", + "SimCalorimeterHit:EcalEndcapRingCollection", + "SimCalorimeterHit:EcalEndcapRingPreShowerCollection", + "SimCalorimeterHit:EcalEndcapSiliconCollection", + "SimCalorimeterHit:EcalEndcapSiliconPreShowerCollection", + "SimCalorimeterHit:HcalBarrelCollection", + "SimCalorimeterHit:HcalEndCapRingsCollection", + "SimCalorimeterHit:HcalEndCapsCollection", + "SimCalorimeterHit:LumiCalCollection", + "SimCalorimeterHit:MuonEndCapCollection", + "Track:ClupatraTrackSegments", + "Track:ClupatraTracks", + "Track:ForwardTracks", + "Track:MarlinTrkTracks", + "Track:SiTracks", + "Track:SubsetTracks", + "TrackerHit:FTDSpacePoints", + "TrackerHit:SETSpacePoints", + "TrackerHit:SITSpacePoints", + "TrackerHit:TPCTrackerHits", + "TrackerHitPlane:FTDPixelTrackerHits", + "TrackerHitPlane:FTDStripTrackerHits", + "TrackerHitPlane:SETTrackerHits", + "TrackerHitPlane:SITTrackerHits", + "TrackerHitPlane:VXDTrackerHits", + "SimTrackerHit:COILCollection", + "SimTrackerHit:FTD_PIXELCollection", + "SimTrackerHit:FTD_STRIPCollection", + "SimTrackerHit:SETCollection", + "SimTrackerHit:SITCollection", + "SimTrackerHit:TPCCollection", + "SimTrackerHit:TPCSpacePointCollection", + "SimTrackerHit:VXDCollection", + "CalorimeterHit:LCAL", + "CalorimeterHit:LHCAL" +] +######################################################################### +geometry_option = "CepC_v4.xml" + +if not os.getenv("DETCEPCV4ROOT"): + print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import GeomSvc +geosvc = GeomSvc("GeomSvc") +geosvc.compact = geometry_path +######################################################################## +from Configurables import GearSvc +gearSvc = GearSvc("GearSvc") +# gearSvc.GearXMLFile = "Detector/DetCEPCv4/compact/FullDetGear.xml" +# gearSvc.GearXMLFile = "/cefs/higgs/PFAData/Sample/Generator/FullGeo_Baseline/E240_nnH_gg/Simu/GearOutput.xml" +gearSvc.GearXMLFile = "Detector/DetCEPCv4/compact/CEPCV4_FullDet_GearOutput.xml" +############################################################################## +from Configurables import G2CDArborAlg +caloDigi = G2CDArborAlg("G2CDArborAlg") +caloDigi.ReadLCIO = True +caloDigi.ECALCollections = ["EcalBarrelSiliconCollection","EcalEndcapSiliconCollection","EcalEndcapRingCollection"] +caloDigi.HCALCollections = ["HcalBarrelCollection","HcalEndCapsCollection","HcalEndCapRingsCollection"] +caloDigi.DigiECALCollection = ["ECALBarrel","ECALEndcap","ECALOther"] +caloDigi.DigiHCALCollection = ["HCALBarrel","HCALEndcap","HCALOther"] + +caloDigi.CalibrECAL = [48.16, 96.32] +caloDigi.HCALThreshold = 0.12 + +caloDigi.PolyaParaA = 1.1 +caloDigi.PolyaParaB = 1.0 +caloDigi.PolyaParaC = 0.0 + +caloDigi.EventReportEvery = 1 +############################################################################## +from Configurables import MarlinArbor +marlinArbor = MarlinArbor("MarlinArbor") +marlinArbor.ReadLCIO = True +marlinArbor.ECALCollections =["ECALBarrel","ECALEndcap","ECALOther"] +marlinArbor.HCALCollections =["HCALBarrel","HCALEndcap","HCALOther"] +marlinArbor.ECALReadOutNames= ["EcalBarrelCollection","EcalEndcapsCollection","EcalEndcapRingCollection"] +marlinArbor.HCALReadOutNames= ["HcalBarrelCollection","HcalEndcapsCollection","HcalEndcapRingCollection"] +############################################################################## +from Configurables import BushConnect +bushconnect = BushConnect("BushConnect") +bushconnect.ReadLCIO = True +############################################################################## +from Configurables import TotalInvMass +totalInvM = TotalInvMass("TotalInvMass") +totalInvM.TreeOutputFile = "Examples/options/LCIO_read_Arbor_nnHgg_BMRAna.root" +############################################################################## +from Configurables import PodioOutput +write = PodioOutput("write") +write.filename = "Examples/options/LCIO_read_Arbor_nnHgg_Reco.root" +write.outputCommands = ["keep *"] +######################################################################### +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [read, caloDigi, marlinArbor, bushconnect, totalInvM, write], + EvtSel = 'NONE', + EvtMax = 2, + ExtSvc = [dsvc, geosvc, gearSvc], + OutputLevel=DEBUG +) diff --git a/Reconstruction/PFA/Arbor/CMakeLists.txt b/Reconstruction/PFA/Arbor/CMakeLists.txt index 6508c401d98b504e4f462baaaff87302ef95619d..82dc846271f4a64f1cf60e132ad36d135439213d 100644 --- a/Reconstruction/PFA/Arbor/CMakeLists.txt +++ b/Reconstruction/PFA/Arbor/CMakeLists.txt @@ -4,20 +4,21 @@ gaudi_add_module(Arbor SOURCES src/Arbor.cc src/ArborHit.cc src/ArborTool.cc - src/ArborToolLCIO.cc - src/ClusterAna.cc - src/DetectorPos.cc - src/HelixClassD.cc + src/ArborToolLCIO.cc + src/ClusterAna.cc + src/DetectorPos.cc + src/HelixClassD.cc src/MarlinArbor.cc src/BushConnect.cc LINK EventSeeder GearSvc DataHelperLib + DecoderHelperLib DetInterface Gaudi::GaudiKernel k4FWCore::k4FWCore ${LCContent_LIBRARIES} - ${CLHEP_LIBRARIES} + ${CLHEP_LIBRARIES} ${ROOT_LIBRARIES} ${LCIO_LIBRARIES} ${GEAR_LIBRARIES} diff --git a/Reconstruction/PFA/Arbor/src/Arbor.cc b/Reconstruction/PFA/Arbor/src/Arbor.cc index c51a6fae8f51d261d6a4830b8866922d97fcdb25..0e96c3c7e84f40f2e77bce10ec762a795118a094 100644 --- a/Reconstruction/PFA/Arbor/src/Arbor.cc +++ b/Reconstruction/PFA/Arbor/src/Arbor.cc @@ -8,19 +8,19 @@ //#define DEBUG -using namespace std; +using namespace std; typedef std::unordered_multimap<unsigned,unsigned> HitLinkMap; std::vector<ArborHit> cleanedHits; -std::vector<int> LeafHitsIndex; -std::vector<int> JointHitsIndex; -std::vector<int> StarJointHitsIndex; +std::vector<int> LeafHitsIndex; +std::vector<int> JointHitsIndex; +std::vector<int> StarJointHitsIndex; std::vector<int> IsoHitsIndex; std::vector<int> SimpleSeedHitsIndex; std::vector<int> StarSeedHitsIndex; -std::map<int, int> HitsFlag; // Tell Hits type; +std::map<int, int> HitsFlag; // Tell Hits type; /* HitLinkMap BackLinksMap; @@ -30,15 +30,15 @@ HitLinkMap IterBackLinks; linkcoll Links; linkcoll InitLinks; -linkcoll IterInputLinks; +linkcoll IterInputLinks; linkcoll alliterlinks; -linkcoll links_debug; -linkcoll IterLinks; -linkcoll IterLinks_1; +linkcoll links_debug; +linkcoll IterLinks; +linkcoll IterLinks_1; branchcoll LengthSortBranchCollection; -branchcoll Trees; +branchcoll Trees; -int NHits = 0; +int NHits = 0; void init() { cleanedHits.clear(); @@ -71,7 +71,7 @@ void HitsCleaning( std::vector<ArborHit> inputHits ) void HitsClassification( linkcoll inputLinks ) { int NLinks = inputLinks.size(); - + LeafHitsIndex.clear(); JointHitsIndex.clear(); StarJointHitsIndex.clear(); @@ -106,7 +106,7 @@ void HitsClassification( linkcoll inputLinks ) else if(BeginIndex[i1] == 1 && EndIndex[i1] == 1) { JointHitsIndex.push_back(i1); - HitsFlag[i1] = 4; + HitsFlag[i1] = 4; } else if(BeginIndex[i1] > 1 && EndIndex[i1] == 1) { @@ -130,7 +130,7 @@ void HitsClassification( linkcoll inputLinks ) } else { - cout<<"WARNING: UNCLASSIFIED HITS, Begin Index: "<<BeginIndex[i1]<<", End Index: "<<EndIndex[i1]<<endl; + cout<<"WARNING: UNCLASSIFIED HITS, Begin Index: "<<BeginIndex[i1]<<", End Index: "<<EndIndex[i1]<<endl; } } @@ -139,14 +139,14 @@ void HitsClassification( linkcoll inputLinks ) cout<<"Seed - Simple/Star: "<<SimpleSeedHitsIndex.size()<<" : "<<StarSeedHitsIndex.size()<<endl; cout<<"Joint - Simple/Star: "<<JointHitsIndex.size()<<" : "<<StarJointHitsIndex.size()<<endl; cout<<"Leaves: "<<LeafHitsIndex.size()<<endl; - cout<<"IsoHits: "<<IsoHitsIndex.size()<<endl; - cout<<"TotalHits: "<<NHits<<endl; + cout<<"IsoHits: "<<IsoHitsIndex.size()<<endl; + cout<<"TotalHits: "<<NHits<<endl; #endif } linkcoll LinkClean( std::vector<ArborHit> allhits, linkcoll alllinks ) { - linkcoll cleanedlinks; + linkcoll cleanedlinks; int NLinks = alllinks.size(); int Ncurrhitlinks = 0; @@ -158,7 +158,9 @@ linkcoll LinkClean( std::vector<ArborHit> allhits, linkcoll alllinks ) std::pair<int, int> SelectedPair; TVector3 PosA, PosB, PosDiffAB; - + + cout<<"[YX debug - LinkClean] Input NLinks = "<<NLinks<<endl; + std::vector< std::vector<int> > LinkHits; LinkHits.clear(); for(int s1 = 0; s1 < NHits; s1++) @@ -212,11 +214,15 @@ linkcoll LinkClean( std::vector<ArborHit> allhits, linkcoll alllinks ) void BuildInitLink( std::vector<float>Thresholds ) { + cout<<"[YX debug - BuildInitLink] Thresholds = "<<Thresholds[0]<<", "<<Thresholds[1]<<", "<<Thresholds[2]<<", "<<Thresholds[3]<<endl; + KDTreeLinkerAlgo<unsigned,3> kdtree; typedef KDTreeNodeInfoT<unsigned,3> KDTreeNodeInfo; std::array<float,3> minpos{ {0.0f,0.0f,0.0f} }, maxpos{ {0.0f,0.0f,0.0f} }; std::vector<KDTreeNodeInfo> nodes, found; + cout<<"[YX debug - BuildInitLink] Input NHits = "<<NHits<<endl; + for(int i0 = 0; i0 < NHits; ++i0 ) { const auto& hit = cleanedHits[i0].GetPosition(); @@ -243,26 +249,26 @@ void BuildInitLink( std::vector<float>Thresholds ) nodes.clear(); Links.clear(); //all tmp links - - TVector3 PosA, PosB, PosDiffAB, PosDiffBA; + + TVector3 PosA, PosB, PosDiffAB, PosDiffBA; int NLayer_A = 0; - int NLayer_B = 0; - int NStave_A = 0; + int NLayer_B = 0; + int NStave_A = 0; int NStave_B = 0; int SubD_A = 0; int SubD_B = 0; float Depth_A = 0; - float Depth_B = 0; - float DisAB = 0; - float MagA = 0; + float Depth_B = 0; + float DisAB = 0; + float MagA = 0; float MagB = 0; float ECCorr = 0; int FlagTrkPS = 0; - int FlagPSEE = 0; + int FlagPSEE = 0; int FlagEH = 0; int FlagHH = 0; - int FlagStaveSame = 0; - int FlagStaveDiff = 0; + int FlagStaveSame = 0; + int FlagStaveDiff = 0; for(int i0 = 0; i0 < NHits; i0++) { @@ -275,6 +281,11 @@ void BuildInitLink( std::vector<float>Thresholds ) SubD_A = cleanedHits[i0].GetSubD(); //SubD_Index, 0 = track endpoint, 1 = Ecal, 2 = Hcal, 3 = EcalPS Depth_A = cleanedHits[i0].GetDepth(); + // cout<<"[YX debug - BuildInitLink] Hit "<<i0<<", NLayer_A = "<<NLayer_A<<", NStave_A = "<<NStave_A<<", SubD_A = "<<SubD_A<<", Depth_A = "<<Depth_A + // <<", Pos = ("<<PosA.X()<<", "<<PosA.Y()<<", "<<PosA.Z()<<")"<<endl; + + + const float side = 200; //could also be sub detector dependent const float xplus(PosA.X() + side), xminus(PosA.X() - side); const float yplus(PosA.Y() + side), yminus(PosA.Y() - side); @@ -287,33 +298,40 @@ void BuildInitLink( std::vector<float>Thresholds ) zmin, zmax ); kdtree.search(searchcube,found); + // cout<<"[YX debug - BuildInitLink] Hit "<<i0<<", found.size() = "<<found.size()<<endl; + + int nLink_pushback = 0; + for(unsigned int j0 = 0; j0 < found.size(); j0++) { if( found[j0].data <= (unsigned)i0 ) continue; PosB = cleanedHits[found[j0].data].GetPosition(); NLayer_B = cleanedHits[found[j0].data].GetLayer(); - NStave_B = cleanedHits[found[j0].data].GetStave(); + NStave_B = cleanedHits[found[j0].data].GetStave(); SubD_B = cleanedHits[found[j0].data].GetSubD(); Depth_B = cleanedHits[found[j0].data].GetDepth(); PosDiffAB = PosA - PosB; - PosDiffBA = PosB - PosA; - DisAB = PosDiffAB.Mag(); - ECCorr = 0; + PosDiffBA = PosB - PosA; + DisAB = PosDiffAB.Mag(); + ECCorr = 0; if( (fabs(PosA.Z()) - 2600 )*(fabs(PosB.Z()) - 2600) < 0 ) ECCorr = 40; + // cout<<"[YX debug - BuildInitLink] ---> Hit "<<j0<<", NLayer_B = "<<NLayer_B<<", NStave_B = "<<NStave_B<<", SubD_B = "<<SubD_B<<", Depth_B = "<<Depth_B + // <<", Pos = ("<<PosB.X()<<", "<<PosB.Y()<<", "<<PosB.Z()<<"), DisAB = "<<DisAB<<endl; + // For the XX seed, using triangle method... FlagTrkPS = 0; FlagPSEE = 0; FlagEH = 0; FlagHH = 0; FlagStaveSame = 0; FlagStaveDiff = 0; if( SubD_A*SubD_B == 0 && SubD_A + SubD_B == 3 && DisAB < 180 && (PosDiffAB.Angle(PosA) < 0.1 || PosDiffBA.Angle(PosA) < 0.1) ) FlagTrkPS = 1; - else if( (SubD_A == 1 || SubD_A == 3) && (SubD_B == 1 || SubD_B == 3) && DisAB < Thresholds[0] + 0.05*(Depth_A + Depth_B) ) - FlagPSEE = 1; + else if( (SubD_A == 1 || SubD_A == 3) && (SubD_B == 1 || SubD_B == 3) && DisAB < Thresholds[0] + 0.05*(Depth_A + Depth_B) ) + FlagPSEE = 1; else if( (SubD_A * SubD_B == 2 && DisAB < Thresholds[1] + ECCorr) ) FlagEH = 1; else if( SubD_A == 2 && SubD_B == 2 && DisAB < Thresholds[2] + 0.03*(Depth_A + Depth_B) ) - FlagHH = 1; + FlagHH = 1; if( FlagTrkPS || FlagPSEE || FlagEH || FlagHH ) { @@ -343,7 +361,10 @@ void BuildInitLink( std::vector<float>Thresholds ) } } } - links_debug = Links; + links_debug = Links; + + cout<<"[YX debug - BuildInitLink] Output Links.size() = "<<Links.size()<<endl; + } void LinkIteration( int time ) //Energy corrections, semi-local correction @@ -402,25 +423,25 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction } int NcurrLinks = alliterlinks.size(); - TVector3 hitPos, PosA, PosB, PosDiffAB, PosDiffBA, linkDir; - int NLayer_A = 0; - int NLayer_B = 0; + TVector3 hitPos, PosA, PosB, PosDiffAB, PosDiffBA, linkDir; + int NLayer_A = 0; + int NLayer_B = 0; int NStave_A = 0; - int NStave_B = 0; - int SubD_A = 0; + int NStave_B = 0; + int SubD_A = 0; int SubD_B = 0; int FlagEE = 0; int FlagHH = 0; - int AngleAccIndex = 0; + int AngleAccIndex = 0; float DisAB = 0; - float MagA = 0; + float MagA = 0; float MagB = 0; - std::pair<int, int> currlink; - std::pair<int, int> a_Link; + std::pair<int, int> currlink; + std::pair<int, int> a_Link; std::pair<int, int> a_tmpLink, b_tmpLink; int FlagNoJoint = 0; int FlagNoIso = 0; -// int FlagNoExisting = 0; +// int FlagNoExisting = 0; for(int i = 0; i < NHits; i++) { @@ -432,16 +453,16 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction // std::vector<int> ---> all the hits linked to this hit std::vector< std::vector<int> > hitLinksArray; - hitLinksArray.resize(NHits); + hitLinksArray.resize(NHits); for(int j = 0; j < NcurrLinks; j++) { currlink = alliterlinks[j]; PosA = cleanedHits[ currlink.first ].GetPosition(); PosB = cleanedHits[ currlink.second ].GetPosition(); linkDir = (PosA - PosB); //Links are always from first point to second - verify - linkDir *= 1.0/linkDir.Mag(); + linkDir *= 1.0/linkDir.Mag(); RefDir[currlink.first] += 4*linkDir; //Weights... might be optimized... - RefDir[currlink.second] += 6*linkDir; + RefDir[currlink.second] += 6*linkDir; Nin_hit[currlink.first] ++; Nout_hit[currlink.second] ++; hitLinksArray[currlink.first].push_back(currlink.second); @@ -469,16 +490,16 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction zmin, zmax ); kdtree.search(searchcube,found); - for(unsigned int j1 = 0; j1 < found.size(); j1++) + for(unsigned int j1 = 0; j1 < found.size(); j1++) { if( found[j1].data <= (unsigned)i1 ) continue; FlagNoJoint = Nout_hit[j1] * Nin_hit[i1] * Nout_hit[i1] * Nin_hit[j1]; FlagNoIso = Nout_hit[j1] + Nin_hit[i1] + Nout_hit[i1] + Nin_hit[j1]; - a_tmpLink.first = i1; - a_tmpLink.second = j1; - b_tmpLink.first = j1; - b_tmpLink.second = i1; + a_tmpLink.first = i1; + a_tmpLink.second = j1; + b_tmpLink.first = j1; + b_tmpLink.second = i1; bool isConnected = false; @@ -510,24 +531,24 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction if( FlagNoJoint == 0 && FlagNoIso != 0 ) { - if(!isConnected) - { + if(!isConnected) + { PosB = cleanedHits[found[j1].data].GetPosition(); NLayer_B = cleanedHits[found[j1].data].GetLayer(); NStave_B = cleanedHits[found[j1].data].GetStave(); SubD_B = cleanedHits[found[j1].data].GetSubD(); - PosDiffAB = PosB - PosA; - PosDiffBA = PosA - PosB; - DisAB = PosDiffAB.Mag(); + PosDiffAB = PosB - PosA; + PosDiffBA = PosA - PosB; + DisAB = PosDiffAB.Mag(); - FlagEE = 0; - FlagHH = 0; + FlagEE = 0; + FlagHH = 0; AngleAccIndex = 0; if( PosDiffAB.Angle(RefDir[i1]) < 0.6/time ) AngleAccIndex = 1; else if( PosDiffAB.Angle(RefDir[j1]) < 0.6/time ) - AngleAccIndex = 2; + AngleAccIndex = 2; else if( PosDiffBA.Angle(RefDir[i1]) < 0.6/time ) AngleAccIndex = 3; else if( PosDiffBA.Angle(RefDir[j1]) < 0.6/time ) @@ -547,7 +568,7 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction MagB = PosB.Mag(); // if(NLayer_A != NLayer_B) - // { + // { if( NStave_A != NStave_B || ( NLayer_A == 0 && NLayer_B != 0 ) || ( NLayer_B == 0 && NLayer_A != 0 ) ) { if( MagA > MagB && AngleAccIndex < 3 ) @@ -581,7 +602,7 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction // } } } - } + } } } @@ -589,11 +610,11 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction int NLinks = alliterlinks.size(); int MinAngleIndex = -10; - int Ncurrhitlinks = 0; - float MinAngle = 1E6; + int Ncurrhitlinks = 0; + float MinAngle = 1E6; float tmpOrder = 0; - float DirAngle = 0; - std::pair<int, int> SelectedPair; + float DirAngle = 0; + std::pair<int, int> SelectedPair; std::vector< std::vector<int> > LinkHits; LinkHits.clear(); @@ -639,7 +660,7 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction if(SelectedPair.first == SelectedPair.second) { cout<<"WTTTTTTTTFFFFFFFFFFFFFF"<<endl; - continue; + continue; } if(time == 1) { @@ -651,7 +672,7 @@ void LinkIteration( int time ) //Energy corrections, semi-local correction IterBackLinks.emplace(SelectedPair.second,SelectedPair.first); } } - } + } #ifdef DEBUG cout<<"Init-Iter Size "<<InitLinks.size()<<" : "<<IterLinks.size()<<endl; @@ -693,7 +714,7 @@ void BranchBuilding(float SeedThreshold) if(HitBeginIndex[i2] == 0 && HitEndIndex[i2] == 1) //EndPoint { NBranches ++; - std::vector<int> currBranchhits; //array of indexes + std::vector<int> currBranchhits; //array of indexes iterhitindex = i2; currBranchhits.push_back(i2); @@ -711,7 +732,7 @@ void BranchBuilding(float SeedThreshold) iterhitindex = PairIterator.first; break; } - } + } */ auto iterlink_range = IterBackLinks.equal_range(iterhitindex); assert( std::distance(iterlink_range.first,iterlink_range.second) == 1 ); @@ -902,22 +923,22 @@ void BushMerging() std::vector< std::vector<int> > Arbor( std::vector<ArborHit> inputHits, std::vector<float>Thresholds ) { //cout<<"Thresholds"<<Thresholds[0]<<" "<<Thresholds[1]<<" "<<Thresholds[2]<<endl; - + if(Thresholds.size() != 4) { - cout<<"Threshold Set Wrong, Threshold Size is "<<Thresholds.size()<<" Should be 4"<<endl; - exit(2); + cout<<"Threshold Set Wrong, Threshold Size is "<<Thresholds.size()<<" Should be 4"<<endl; + exit(2); } - init(); + init(); HitsCleaning(inputHits); BuildInitLink(Thresholds); InitLinks = LinkClean( cleanedHits, Links ); LinkIteration(1); LinkIteration(2); - HitsClassification(IterLinks); - BranchBuilding(Thresholds[3]); + HitsClassification(IterLinks); + BranchBuilding(Thresholds[3]); return LengthSortBranchCollection; } diff --git a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc index 6a1e0cc20da795c05f1e7bc328060f5b232f91a8..4725b56c9552cab7a33b02663739cb6aae4731a1 100644 --- a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc +++ b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc @@ -9,6 +9,9 @@ #include <stdexcept> #include <sstream> +#include "DDSegmentation/MegatileLayerGridXY.h" +#include "DecoderHelper/DD4hep2Lcio.h" + #include "TVector3.h" #include <string> #include <iostream> @@ -31,6 +34,7 @@ #include "DD4hep/IDDescriptor.h" #include "DD4hep/Plugins.h" +#include "cellIDDecoder.h" #include <DDRec/DetectorData.h> #include <DDRec/CellIDPositionConverter.h> #include "DetInterface/IGeomSvc.h" @@ -38,7 +42,61 @@ #include "podio/podioVersion.h" using namespace std; -/* +// bool _USE_LCIO=0; + +// struct segInfo; +struct segInfo { + double megaTileSizeX = 0; + double megaTileSizeY = 0; + double megaTileOffsetX = 0; + double megaTileOffsetY = 0; + unsigned int nCellsX = 0; + unsigned int nCellsY = 0; + segInfo() = default; +}; + +struct SegParameters { +// unsigned int layer; +// unsigned int tile; + double sizex; + double sizey; + double offsetx; + double offsety; + unsigned int ncellsx; + unsigned int ncellsy; +}; + +class MyMegatileLayerGridXY: public dd4hep::DDSegmentation::MegatileLayerGridXY { +public: + int get_nCellsX(){ + return _unif_nCellsX; + } + + SegParameters getSegParameters(){ + SegParameters Pars; + // Pars.layer = 1; + // Pars.tile = 2; + Pars.sizex = _megaTileSizeX; + Pars.sizey = _megaTileSizeY; + Pars.offsetx = _megaTileOffsetX; + Pars.offsety = _megaTileOffsetY; + Pars.ncellsx = _unif_nCellsX; + Pars.ncellsy = _unif_nCellsY; + return Pars; + } + + std::map < std::pair < unsigned int, unsigned int > , segInfo > get_specialMegaTiles_layerWafer(){ + return specialMegaTiles_layerWafer; + } + // void setSpecialMegaTile( unsigned int layer, unsigned int tile, + // double sizex, double sizey, + // double offsetx, double offsety, + // unsigned int ncellsx, unsigned int ncellsy ); + + +}; + +/* void ClusterBuilding( LCEvent * evtPP, std::string Name, std::vector<CalorimeterHit*> Hits, std::vector< std::vector<int> > BranchOrder, int DHCALFlag ) { LCCollection *currbranchcoll = new LCCollectionVec(LCIO::CLUSTER); @@ -51,8 +109,8 @@ void ClusterBuilding( LCEvent * evtPP, std::string Name, std::vector<Calorimeter float currBranchEnergy = 0; TVector3 SeedPos, currPos; float MinMag = 1E9; - float currMag = 0; - float ECALTotalEn = 0; + float currMag = 0; + float ECALTotalEn = 0; float HCALTotalEn = 0; for(int i0 = 0; i0 < NBranch; i0++) @@ -113,12 +171,14 @@ void ClusterBuilding( LCEvent * evtPP, std::string Name, std::vector<Calorimeter */ -ArborToolLCIO::ArborToolLCIO(const std::string& name,ISvcLocator* svcLoc) +ArborToolLCIO::ArborToolLCIO(const std::string& name,ISvcLocator* svcLoc, bool m_readLCIO) : GaudiAlgorithm(name, svcLoc) { m_geosvc=service<IGeomSvc>("GeomSvc"); + m_encoder_str = "M:3,S-1:3,I:9,J:9,K-1:6"; + _USE_LCIO = m_readLCIO; } - + ArborToolLCIO::~ArborToolLCIO() { @@ -135,10 +195,12 @@ void ArborToolLCIO::ClusterBuilding( DataHandle<edm4hep::ClusterCollection>& _cu TVector3 SeedPos, currPos; float MinMag = 1E9; - float currMag = 0; - float ECALTotalEn = 0; + float currMag = 0; + float ECALTotalEn = 0; float HCALTotalEn = 0; + cout<<"[YX debug - ClusterBuilding] NBranch = "<<NBranch<<endl; + for(int i0 = 0; i0 < NBranch; i0++) { auto a_branch=currbranchcoll->create(); @@ -149,6 +211,10 @@ void ArborToolLCIO::ClusterBuilding( DataHandle<edm4hep::ClusterCollection>& _cu HCALTotalEn = 0; MinMag = 1E9; + + // cout<<"[YX debug - ClusterBuilding] Branch "<<i0<<", BranchSize = "<<BranchSize<<endl; + + for(int j = 0; j < BranchSize; j++) { auto a_hit= Hits[currbranchorder[j]]; @@ -189,6 +255,10 @@ void ArborToolLCIO::ClusterBuilding( DataHandle<edm4hep::ClusterCollection>& _cu } + + int nBranch_out = currbranchcoll->size(); + cout<<"[YX debug - ClusterBuilding] nBranch_out = "<<nBranch_out<<endl; + } @@ -215,10 +285,25 @@ int ArborToolLCIO::NHScaleV2( std::vector<edm4hep::CalorimeterHit> clu0, int Rat auto hit = clu0[i]; auto cellid = hit.getCellID(); + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + const std::string cellICoding(GetCellICoding(idCodingString)); + const std::string cellJCoding(GetCellJCoding(idCodingString)); + + tmpI=cellIdDecoder(&hit)[cellICoding.c_str()]/RatioX; + tmpJ=cellIdDecoder(&hit)[cellJCoding.c_str()]/RatioY; + tmpK=(cellIdDecoder(&hit)[layerCoding.c_str()]+1)/RatioZ; + } + else{ tmpI = m_decoder->get(cellid, "cellX")/RatioX; tmpJ = m_decoder->get(cellid, "cellY")/RatioY; tmpK = (m_decoder->get(cellid, "layer")+1)/RatioZ; + + + } tmpEn = hit.getEnergy(); NewCellID0 = (tmpK<<24) + (tmpJ<<12) + tmpI; @@ -251,8 +336,8 @@ float ArborToolLCIO::FDV2( std::vector<edm4hep::CalorimeterHit> clu) FractalDim += 0.1 * TMath::Log(float(OriNHit)/NReSizeHit[j])/TMath::Log(float(Scale[j])); } - if(clu.size() == 0) - FractalDim = -1; + if(clu.size() == 0) + FractalDim = -1; return FractalDim; } @@ -269,24 +354,169 @@ int ArborToolLCIO::NHScaleV3( edm4hep::Cluster clu0, int RatioX, int RatioY, int float tmpEn = 0; int NewCellID0 = 0; int NewCellID1 = 0; - //m_geosvc=service<IGeomSvc>("GeomSvc"); + //m_geosvc=service<IGeomSvc>("GeomSvc"); - m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); - if(!m_decoder) m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); // raw_system==20 + // if(!m_decoder) m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + + + // ---------------------------------------------------------------------------- + dd4hep::Detector* m_dd4hep = m_geosvc->lcdd(); + if ( !m_dd4hep ) throw "ArborToolLCIO: Failed to get dd4hep::Detector ..."; + dd4hep::Readout readout = m_dd4hep->readout("EcalBarrelCollection"); + + // auto m_segmentation = dynamic_cast<dd4hep::DDSegmentation::MegatileLayerGridXY*>(readout.segmentation().segmentation()); + auto m_segmentation = static_cast<MyMegatileLayerGridXY*>(readout.segmentation().segmentation()); + auto m_segPars = m_segmentation->getSegParameters(); + auto m_seg = m_segmentation->get_specialMegaTiles_layerWafer(); + + // std::cout<<m_segmentation<<std::endl; + if(0) std::cout<<"get_nCellsX = "<<m_segmentation->get_nCellsX()<<std::endl; + if(0) std::cout<<"m_segPars->ncellsx = "<<m_segPars.ncellsx<<std::endl; + // ---------------------------------------------------------------------------- std::map <double, float> testIDtoEnergy; double testlongID = 0; - + for(int i = 0; i < NumHit; i++) { auto hit = clu0.getHits(i); auto cellid = hit.getCellID(); - - tmpI = m_decoder->get(cellid, "cellX")/RatioX; - tmpJ = m_decoder->get(cellid, "cellY")/RatioY; - tmpK = (m_decoder->get(cellid, "layer")+1)/RatioZ; tmpEn = hit.getEnergy(); + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + const std::string cellICoding(GetCellICoding(idCodingString)); + const std::string cellJCoding(GetCellJCoding(idCodingString)); + + tmpI=cellIdDecoder(&hit)[cellICoding.c_str()]/RatioX; + tmpJ=cellIdDecoder(&hit)[cellJCoding.c_str()]/RatioY; + tmpK=(cellIdDecoder(&hit)[layerCoding.c_str()]+1)/RatioZ; + } + else{ + // initial decoder: ECAL barrel + // <id>system:5,module:3,stave:4,tower:5,layer:6,wafer:6,cellX:32:-16,cellY:-16</id> + m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); // raw_system==20 + int raw_system = m_decoder->get(cellid, "system"); + int raw_module = m_decoder->get(cellid, "module"); + int raw_stave = m_decoder->get(cellid, "stave"); + int raw_layer = m_decoder->get(cellid, "layer"); + int raw_cellX = m_decoder->get(cellid, "cellX"); + int raw_cellY = m_decoder->get(cellid, "cellY"); + int raw_tower = m_decoder->get(cellid, "tower"); + int raw_wafer = m_decoder->get(cellid, "wafer"); + + if(raw_system==29){ // ECAL endcap + // <id>system:5,module:3,stave:4,tower:5,layer:6,wafer:6,x:32:-16,y:-16</id> + m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + raw_stave = m_decoder->get(cellid, "stave"); + raw_layer = m_decoder->get(cellid, "layer"); + raw_cellX = m_decoder->get(cellid, "x"); + raw_cellY = m_decoder->get(cellid, "y"); + raw_tower = m_decoder->get(cellid, "tower"); + raw_wafer = m_decoder->get(cellid, "wafer"); + }else if(raw_system==22){ // HCAL barrel + // <id>system:5,module:3,stave:3,tower:5,layer:6,slice:4,x:32:-16,y:-16</id> + m_decoder = m_geosvc->getDecoder("HcalBarrelCollection"); + raw_stave = m_decoder->get(cellid, "stave"); + raw_layer = m_decoder->get(cellid, "layer"); + raw_cellX = m_decoder->get(cellid, "x"); + raw_cellY = m_decoder->get(cellid, "y"); + }else if(raw_system==30){ // HCAL endcap + // <id>system:5,module:3,stave:3,tower:5,layer:6,y:32:-16,x:-16</id> + m_decoder = m_geosvc->getDecoder("HcalEndcapsCollection"); + raw_stave = m_decoder->get(cellid, "stave"); + raw_layer = m_decoder->get(cellid, "layer"); + raw_cellX = m_decoder->get(cellid, "y"); + raw_cellY = m_decoder->get(cellid, "x"); + } + + + + // --------------------------------------------- + std::pair<unsigned int, unsigned int> layer_wafer = std::make_pair(raw_layer, raw_wafer); + // segInfo seginfo = m_seg.find(layer_wafer); + // --------------------------------------------- + + TVector3 currHitPos = TVector3(hit.getPosition().x, hit.getPosition().y, hit.getPosition().z); + double hitposx = currHitPos.X(); + double hitposy = currHitPos.Y(); + double hitposz = currHitPos.Z(); + double hitposp = currHitPos.Perp(); + + // --------------------------------------------- + + int new_layer = DD4hep2Lcio::CEPCv4::getEcalLayer(raw_layer); + int new_stave = DD4hep2Lcio::CEPCv4::getEcalBarrelStave(raw_stave); + if(raw_system==29){// ECAL endcap + new_stave = DD4hep2Lcio::CEPCv4::getEcalEndcapStave(raw_stave); + } + if(raw_system==22 || raw_system==30){ // HCAL, barrel 22, endcap 30 + new_layer = DD4hep2Lcio::CEPCv4::getHcalLayer(raw_layer); + new_stave = DD4hep2Lcio::CEPCv4::getHcalStave(raw_stave); + } + + int new_I = raw_cellX; + int new_J = raw_cellY; + int new_K = new_layer; + + int nwaferx = (raw_wafer-1)/2; + int nwafery = raw_wafer%2 ==0 ? 1 : 0; + int ncellsx = 9; + int ncellsy = 9; + // int cellX_constant = raw_cellX<0 ? 4 : 0; + // int cellY_constant = raw_cellY<0 ? 4 : 0; + int cellX_constant = 0; + int cellY_constant = 0; + + + // if(raw_system==20 || raw_system==22){ // barrel, ECAL 20, HCAL 22 + if(raw_system==20){ // barrel, ECAL 20 + new_I = raw_cellX + nwaferx * 9 + cellX_constant; + new_J = (8-raw_cellY) + (4-raw_tower) * 2 * 9 + (1-nwafery) * 9 + cellY_constant; + + // }else if(raw_system==29 || raw_system==30){ // endcap, ECAL 29, HCAL 30 + }else if(raw_system==29){ // endcap, ECAL 29 + // new_I = raw_cellX + raw_tower * 2 * 9 + nwaferx * 9 + 4; + // new_J = raw_cellY + nwafery * 9 + 4; + + new_I = raw_cellY + raw_tower * 2 * 9 + nwafery * 9 + 4; + new_J = raw_cellX + nwaferx * 9 + 4; + } + + if(0){ + + cout<<"\n[ArborToolLCIO::NHScaleV3] Hit Pos ("<<hitposx<<", "<<hitposy<<", "<<hitposz<<", "<<hitposp<<")mm:"<<endl; + + cout<<"[ArborToolLCIO::NHScaleV3] ---> Raw M = "<<raw_module<<", stave = "<<raw_stave<<", layer = "<<raw_layer<<", I(x) = "<<raw_cellX<<", J(y) = "<<raw_cellY<<", wafer = "<<raw_wafer<<", tower = "<<raw_tower<<", system = "<<raw_system<<endl; + + cout<<"[ArborToolLCIO::NHScaleV3] ------> system = "<<raw_system<<", layer = "<<raw_layer<<", tower = "<<raw_tower<<", wafer = "<<raw_wafer<<", nwaferx = "<<nwaferx<<", nwafery = "<<nwafery<<", ncellsx = "<<ncellsx<<", ncellsy = "<<ncellsx<<", const x = "<<cellX_constant<<", y = "<<cellY_constant<<endl; + + cout<<"[ArborToolLCIO::NHScaleV3] ---> New M = "<<raw_module<<", stave = "<<new_stave<<", layer = "<<new_layer<<", I(x) = "<<new_I<<", J(y) = "<<new_J<<", wafer = "<<raw_wafer<<", tower = "<<raw_tower<<", system = "<<raw_system<<endl; + + + // cout<<"[ArborToolLCIO::NHScaleV3] ---> megaTileSizeX = "<<seginfo.megaTileSizeX<<", megaTileSizeY = "<<seginfo.megaTileSizeY<<", megaTileOffsetX = "<<seginfo.megaTileOffsetX<<", megaTileOffsetY = "<<seginfo.megaTileOffsetY<<", nCellsX = "<<seginfo.nCellsX<<", nCellsY = "<<seginfo.nCellsY<<endl; + // cout<<"[ArborToolLCIO::NHScaleV3] ---> sizex = "<<m_segPars.sizex<<", sizey = "<<m_segPars.sizey<<", offsetx = "<<m_segPars.offsetx<<", offsety = "<<m_segPars.offsety<<", ncellsx = "<<m_segPars.ncellsx<<", ncellsy = "<<m_segPars.ncellsy<<endl; + } + + // ---------------------------------------------------------------------------- + // tmpI = m_decoder->get(cellid, "cellX")/RatioX; + // tmpJ = m_decoder->get(cellid, "cellY")/RatioY; + // tmpK = (m_decoder->get(cellid, "layer")+1)/RatioZ; + + tmpI = new_I/RatioX; + tmpJ = new_J/RatioY; + tmpK = (new_K+1)/RatioZ; + + } + + if(0) cout<<"[ArborToolLCIO::NHScaleV3] ---> RatioX = "<<RatioX<<", RatioY = "<<RatioY<<", RatioZ = "<<RatioZ<<endl; + if(0) cout<<"[ArborToolLCIO::NHScaleV3] ---> tmpI = "<<tmpI<<", tmpJ = "<<tmpJ<<", tmpK = "<<tmpK<<endl; + + + NewCellID0 = (tmpK<<24) + (tmpJ<<12) + tmpI; testlongID = NewCellID1*1073741824 + NewCellID0; @@ -327,14 +557,14 @@ float ArborToolLCIO::FDV3( edm4hep::Cluster clu ){ float ArborToolLCIO::BushDis( edm4hep::Cluster clu1, edm4hep::Cluster clu2) { - float DisBetweenBush = 1.0E10; + float DisBetweenBush = 1.0E10; int cluSize1 = clu1.hits_size(); int cluSize2 = clu2.hits_size(); - TVector3 HitPos1, HitPos2; - TVector3 PosDiff; - // TVector3 XXXPos; + TVector3 HitPos1, HitPos2; + TVector3 PosDiff; + // TVector3 XXXPos; for(int i = 0; i < cluSize1; i++) { @@ -352,41 +582,41 @@ float ArborToolLCIO::BushDis( edm4hep::Cluster clu1, edm4hep::Cluster clu2) } - return DisBetweenBush; + return DisBetweenBush; } float ArborToolLCIO::DisPointToBush(TVector3 Pos1, edm4hep::Cluster clu1) { - float Dis = 1.0E9; + float Dis = 1.0E9; float HitDis = 1.0E8; int clusize = clu1.hits_size(); - TVector3 HitPos; + TVector3 HitPos; for(int s = 0; s < clusize; s++) { HitPos = TVector3((clu1.getHits(s)).getPosition().x,(clu1.getHits(s)).getPosition().y,(clu1.getHits(s)).getPosition().z); HitDis = (HitPos - Pos1).Mag(); - if(HitDis < Dis) + if(HitDis < Dis) { - Dis = HitDis; + Dis = HitDis; } } - return Dis; + return Dis; } TVector3 ArborToolLCIO::ClusterCoG(edm4hep::Cluster inputCluster) { - TVector3 CenterOfGravity; + TVector3 CenterOfGravity; int inputClusterSize = inputCluster.hits_size(); - TVector3 tmphitPos; + TVector3 tmphitPos; float tmphitEnergy; - float sumhitEnergy = 0; + float sumhitEnergy = 0; for(int i = 0; i < inputClusterSize; i++) { @@ -395,12 +625,12 @@ TVector3 ArborToolLCIO::ClusterCoG(edm4hep::Cluster inputCluster) tmphitEnergy = tmpHit.getEnergy(); CenterOfGravity += tmphitPos*tmphitEnergy; - sumhitEnergy += tmphitEnergy; + sumhitEnergy += tmphitEnergy; } - CenterOfGravity = 1.0/sumhitEnergy * CenterOfGravity; + CenterOfGravity = 1.0/sumhitEnergy * CenterOfGravity; - return CenterOfGravity; + return CenterOfGravity; } @@ -410,10 +640,10 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecColl( std::vector<edm4hep:: edm4hep::ClusterCollection* vec_coll_Clusters = m_clucol.createAndPut(); int NClu = inputClusters.size(); - int CurrBranchSize = 0; - TVector3 SeedPos; - std::vector<float> CluEn; - std::vector<int> CluIndex; + int CurrBranchSize = 0; + TVector3 SeedPos; + std::vector<float> CluEn; + std::vector<int> CluIndex; for(int i0 = 0; i0 < NClu; i0++) { @@ -448,7 +678,7 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecColl( std::vector<edm4hep:: std::vector<edm4hep::Cluster> ArborToolLCIO::CollClusterVec(const edm4hep::ClusterCollection * input_coll ) { - std::vector<edm4hep::Cluster> outputClusterVec; + std::vector<edm4hep::Cluster> outputClusterVec; outputClusterVec.clear(); @@ -459,7 +689,7 @@ std::vector<edm4hep::Cluster> ArborToolLCIO::CollClusterVec(const edm4hep::Clust outputClusterVec.push_back(a_clu); } - return outputClusterVec; + return outputClusterVec; } @@ -468,7 +698,7 @@ void ArborToolLCIO::NaiveCluConst(edm4hep::MutableCluster a0_clu,edm4hep::Mutabl b0_clu.setPosition(a0_clu.getPosition()); b0_clu.setEnergy(a0_clu.getEnergy()); int NCaloHit = a0_clu.hits_size(); - float HitEn = 0; + float HitEn = 0; float SubDEn[6] = {0, 0, 0, 0, 0, 0}; for(int t0 = 0; t0 < NCaloHit; t0++) @@ -481,7 +711,7 @@ void ArborToolLCIO::NaiveCluConst(edm4hep::MutableCluster a0_clu,edm4hep::Mutabl SubDEn[1] += HitEn; } else - { + { SubDEn[0] += HitEn; } } @@ -490,7 +720,7 @@ void ArborToolLCIO::NaiveCluConst(edm4hep::MutableCluster a0_clu,edm4hep::Mutabl { b0_clu.addToSubdetectorEnergies(SubDEn[i]); } - + } @@ -500,7 +730,7 @@ edm4hep::Cluster ArborToolLCIO::NaiveCluImpl(edm4hep::MutableCluster a0_clu) b0_clu.setPosition(a0_clu.getPosition()); b0_clu.setEnergy(a0_clu.getEnergy()); int NCaloHit = a0_clu.hits_size(); - float HitEn = 0; + float HitEn = 0; float SubDEn[6] = {0, 0, 0, 0, 0, 0}; for(int t0 = 0; t0 < NCaloHit; t0++) @@ -513,7 +743,7 @@ edm4hep::Cluster ArborToolLCIO::NaiveCluImpl(edm4hep::MutableCluster a0_clu) SubDEn[1] += HitEn; } else - { + { SubDEn[0] += HitEn; } } @@ -522,8 +752,8 @@ edm4hep::Cluster ArborToolLCIO::NaiveCluImpl(edm4hep::MutableCluster a0_clu) { b0_clu.addToSubdetectorEnergies(SubDEn[i]); } - - return b0_clu; + + return b0_clu; } std::vector<edm4hep::CalorimeterHit> ArborToolLCIO::CollHitVec(const edm4hep::CalorimeterHitCollection * input_coll, float EnergyThreshold) @@ -545,23 +775,28 @@ std::vector<edm4hep::CalorimeterHit> ArborToolLCIO::CollHitVec(const edm4hep::Ca } -std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::vector<edm4hep::Cluster> MainClusters, std::vector<edm4hep::CalorimeterHit> IsoHits, float DisThreshold ) // Projective Distance + Hit Depth correlation; +std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::vector<edm4hep::Cluster> MainClusters, std::vector<edm4hep::CalorimeterHit> IsoHits, float DisThreshold ) // Projective Distance + Hit Depth correlation; { std::vector<edm4hep::MutableCluster> outputClusterVec; int N_Core = MainClusters.size(); int N_Hit = IsoHits.size(); TVector3 HitPos, MBSeedPos; - float currHitCoreDis = 0; - float MinHitCoreDis = 1.0E10; - int MinDisIndex = -1; + float currHitCoreDis = 0; + float MinHitCoreDis = 1.0E10; + int MinDisIndex = -1; std::vector<std::pair<int, int> > Frag_Core_Links; std::pair<int, int> a_frag_core_link; + + double ESum_IsoHit = 0; + double ESum_Core = 0; + + for(int i0 = 0; i0 < N_Hit; i0++) { auto a_hit = IsoHits[i0]; - HitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + HitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); MinHitCoreDis = 1.0E10; for(int j0 = 0; j0 < N_Core; j0++) @@ -570,7 +805,7 @@ std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::v currHitCoreDis = DisPointToBush(HitPos, a_core); if(currHitCoreDis < MinHitCoreDis) { - MinHitCoreDis = currHitCoreDis; + MinHitCoreDis = currHitCoreDis; MinDisIndex = j0; } } @@ -580,11 +815,14 @@ std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::v a_frag_core_link.second = MinDisIndex; Frag_Core_Links.push_back(a_frag_core_link); } + + ESum_IsoHit += a_hit.getEnergy(); + } int N_frag_core_links = Frag_Core_Links.size(); std::vector<edm4hep::CalorimeterHit> tomerge_hits; - float ClusterEn = 0; + float ClusterEn = 0; for(int i2 = 0; i2 < N_Core; i2 ++) { @@ -601,7 +839,7 @@ std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::v } } edm4hep::MutableCluster a_mergedfrag_core; - ClusterEn = 0; + ClusterEn = 0; for(unsigned int j2 = 0; j2 < a_core.hits_size(); j2++) { @@ -624,9 +862,17 @@ std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterHitAbsorbtion( std::v a_mergedfrag_core.setPhi( MBSeedPos.Phi() ); outputClusterVec.push_back(a_mergedfrag_core); + + + ESum_Core += ClusterEn; + // cout<<"[YX debug - ClusterHitAbsorbtion] Core "<<i2<<", merged E = "<<ClusterEn<<endl; + } - return outputClusterVec; + cout<<"[YX debug - ClusterHitAbsorbtion] N_Core = "<<N_Core<<", N_IsoHit = "<<N_Hit<<", ESum_IsoHit = "<<ESum_IsoHit<<", ESum_Core = "<<ESum_Core<<endl; + + + return outputClusterVec; } @@ -634,10 +880,10 @@ void ArborToolLCIO::NaiveMergeCluConst(std::vector<edm4hep::Cluster> inputCluVec { int NClu = inputCluVec.size(); - float SeedDis = 1E9; - float MaxDis = 0; - float MergedCluEnergy = 0; - float HitEn = 0; + float SeedDis = 1E9; + float MaxDis = 0; + float MergedCluEnergy = 0; + float HitEn = 0; float SubDEn[6] = {0, 0, 0, 0, 0, 0}; TVector3 CurrSeedPos, SeedPos, CurrHitPos, CluEndPos, CluRefDir; //Seed Depth... CoG Comp... @@ -650,7 +896,7 @@ void ArborToolLCIO::NaiveMergeCluConst(std::vector<edm4hep::Cluster> inputCluVec if(CurrSeedPos.Mag() < SeedDis) { - SeedPos = CurrSeedPos; + SeedPos = CurrSeedPos; SeedDis = CurrSeedPos.Mag(); } @@ -662,19 +908,19 @@ void ArborToolLCIO::NaiveMergeCluConst(std::vector<edm4hep::Cluster> inputCluVec MergedClu.addToHits(a_hit); CurrHitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); HitEn = a_hit.getEnergy(); - if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) // ECAL, HCAL, Should use better criteria. + if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) // ECAL, HCAL, Should use better criteria. { - SubDEn[1] += HitEn; + SubDEn[1] += HitEn; } else { - SubDEn[0] += HitEn; + SubDEn[0] += HitEn; } if(CurrHitPos.Mag() > MaxDis) { MaxDis = CurrHitPos.Mag(); - CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); } } } @@ -699,10 +945,10 @@ edm4hep::MutableCluster ArborToolLCIO::NaiveMergeClu(std::vector<edm4hep::Cluste edm4hep::MutableCluster MergedClu; int NClu = inputCluVec.size(); - float SeedDis = 1E9; - float MaxDis = 0; - float MergedCluEnergy = 0; - float HitEn = 0; + float SeedDis = 1E9; + float MaxDis = 0; + float MergedCluEnergy = 0; + float HitEn = 0; float SubDEn[6] = {0, 0, 0, 0, 0, 0}; TVector3 CurrSeedPos, SeedPos, CurrHitPos, CluEndPos, CluRefDir; //Seed Depth... CoG Comp... @@ -715,7 +961,7 @@ edm4hep::MutableCluster ArborToolLCIO::NaiveMergeClu(std::vector<edm4hep::Cluste if(CurrSeedPos.Mag() < SeedDis) { - SeedPos = CurrSeedPos; + SeedPos = CurrSeedPos; SeedDis = CurrSeedPos.Mag(); } @@ -727,19 +973,19 @@ edm4hep::MutableCluster ArborToolLCIO::NaiveMergeClu(std::vector<edm4hep::Cluste MergedClu.addToHits(a_hit); CurrHitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); HitEn = a_hit.getEnergy(); - if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) // ECAL, HCAL, Should use better criteria. + if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) // ECAL, HCAL, Should use better criteria. { - SubDEn[1] += HitEn; + SubDEn[1] += HitEn; } else { - SubDEn[0] += HitEn; + SubDEn[0] += HitEn; } if(CurrHitPos.Mag() > MaxDis) { MaxDis = CurrHitPos.Mag(); - CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); } } } @@ -772,14 +1018,14 @@ std::vector<edm4hep::MutableCluster> ArborToolLCIO::ClusterAbsorbtion( std::vect //tag minimal distance - float MinFragCoreDis = 1.0E10; + float MinFragCoreDis = 1.0E10; float CurrFragCoreDis = 0; - int MinDisIndex = -1; + int MinDisIndex = -1; std::vector<std::pair<int, int> > Frag_Core_Links; std::pair<int, int> a_frag_core_link; std::map<int, int> TouchedFrag; TouchedFrag.clear(); - TVector3 fragPos; + TVector3 fragPos; for(int i0 = 0; i0 < N_frag; i0 ++) { @@ -868,7 +1114,7 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecMerge( std::vector<edm4hep: #endif - TVector3 tmpClusterSeedPos, MBSeedPos; + TVector3 tmpClusterSeedPos, MBSeedPos; int CurrBranchSize = 0; float SeedPosMin = 1.0E10; @@ -925,7 +1171,7 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecMerge( std::vector<edm4hep: auto tmpHit = tmpMergebranch.getHits(j1); branchtmp.addToHits(tmpHit); } - + } branchtmp.setPosition( Mainbranch.getPosition() ); branchtmp.setEnergy(BranchEnergy); @@ -942,10 +1188,10 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecMerge( std::vector<edm4hep: } int ArborToolLCIO::JointsBetweenBush(edm4hep::Cluster a_Clu, edm4hep::Cluster b_Clu, float CellSize) { - int NJoint = 0; + int NJoint = 0; int a_CluSize = a_Clu.hits_size(); int b_CluSize = b_Clu.hits_size(); - TVector3 aHitPos, bHitPos, PosDiff, aCluPos, bCluPos; + TVector3 aHitPos, bHitPos, PosDiff, aCluPos, bCluPos; aCluPos = TVector3(a_Clu.getPosition().x,a_Clu.getPosition().y,a_Clu.getPosition().z); bCluPos = TVector3(b_Clu.getPosition().x,b_Clu.getPosition().y,b_Clu.getPosition().z); @@ -957,7 +1203,7 @@ int ArborToolLCIO::JointsBetweenBush(edm4hep::Cluster a_Clu, edm4hep::Cluster b_ { auto bhit = b_Clu.getHits(j); bHitPos = TVector3(bhit.getPosition().x,bhit.getPosition().y,bhit.getPosition().z); - PosDiff = aHitPos - bHitPos; + PosDiff = aHitPos - bHitPos; if(PosDiff.Mag() < 1.5*CellSize) //allow Diag connect... else use 1.2 { // if((aCluPos - bHitPos).Mag() < 60 || (bCluPos - aHitPos).Mag() < 60) @@ -966,7 +1212,7 @@ int ArborToolLCIO::JointsBetweenBush(edm4hep::Cluster a_Clu, edm4hep::Cluster b_ } } - return NJoint; + return NJoint; } @@ -1029,7 +1275,7 @@ float ArborToolLCIO::SimpleBushTimeTrackClu(edm4hep::Track a_trk, edm4hep::Clust int ArborToolLCIO::SimpleBushNC(edm4hep::Track a_trk, edm4hep::Cluster a_clu) { float Distance = 1.0E9; - //float Time = 0; + //float Time = 0; int NC = 0; float BushDist[3] = {0, 0 ,0}; HelixClassD * a_Helix = new HelixClassD(); @@ -1065,7 +1311,7 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) // non-matched 111 // HAD: matched 211 // non-matched 212 - + int CluIDFlag = 999; int ClusterID = 211; int EcalNHit, HcalNHit, NLEcal, NLHcal, NH[16]; @@ -1220,7 +1466,19 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) auto a_hit = a_tree.getHits(s1); allhits.push_back(a_hit); auto cellid= a_hit.getCellID(); - int NLayer = m_decoder->get(cellid, "layer"); + // int NLayer = m_decoder->get(cellid, "layer"); + int NLayer; + + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + NLayer=cellIdDecoder(&a_hit)[layerCoding.c_str()]+1; + } + + else{ + NLayer = m_decoder->get(cellid, "layer"); + } HitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); @@ -1307,11 +1565,11 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) else if(NLayer < 15) { EH_3.push_back(a_hit); - } + } else if(NLayer < 20) { EH_4.push_back(a_hit); - } + } else if(NLayer < 25) { EH_5.push_back(a_hit); @@ -1378,7 +1636,7 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) NLEcal = ActiveLayers(Ecalhits); NLHcal = ActiveLayers(Hcalhits); - cout<<"NLEcal "<<NLEcal<<" "<<Ecalhits.size()<<" "<<endl; + if(0) cout<<"NLEcal "<<NLEcal<<" "<<Ecalhits.size()<<" "<<endl; float sum_NHE = 0, sum_NHH = 0; @@ -1408,7 +1666,7 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) rms_Ecal = 0; rms_Hcal = 0; rms_Ecal2 = 0; - rms_Hcal2 = 0; + rms_Hcal2 = 0; for(int r0 = 0; r0 < 16; r0++) { if(r0 < 6) @@ -1539,21 +1797,36 @@ int ArborToolLCIO::ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk) int ArborToolLCIO::ActiveLayers( std::vector<edm4hep::CalorimeterHit> clu ) { - std::vector<int> hitlayers; + std::vector<int> hitlayers; hitlayers.clear(); - int NHits = clu.size(); + int NHits = clu.size(); int tmpK = 0; //Layer Number - int tmpS = 0; + int tmpS = 0; int tmpID = 0; - + for(int i = 0; i < NHits; i++) { auto hit = clu[i]; auto cellid= hit.getCellID(); + + + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + const std::string staveCoding(GetStaveCoding(idCodingString)); + + + tmpK=cellIdDecoder(&hit)[layerCoding.c_str()]+1; + tmpS=cellIdDecoder(&hit)[staveCoding.c_str()]+1; + } + + else{ tmpK = m_decoder->get(cellid, "layer")+1 ; tmpS = m_decoder->get(cellid, "stave")+1 ; - // cout<<"tmpK "<<tmpK<<endl; + } + // cout<<"tmpK "<<tmpK<<endl; tmpID = tmpS * 50 + tmpK; if( std::find(hitlayers.begin(), hitlayers.end(), tmpID) == hitlayers.end() ) @@ -1568,9 +1841,9 @@ int ArborToolLCIO::ActiveLayers( std::vector<edm4hep::CalorimeterHit> clu ) float ArborToolLCIO::ClusterT0(edm4hep::Cluster a_Clu) { - float T0 = 1E9; - float tmpTime = 0; - TVector3 CluHitPos; + float T0 = 1E9; + float tmpTime = 0; + TVector3 CluHitPos; for(unsigned int i = 0; i < a_Clu.hits_size(); i++) { auto a_hit = a_Clu.getHits(i); @@ -1579,7 +1852,7 @@ float ArborToolLCIO::ClusterT0(edm4hep::Cluster a_Clu) tmpTime = a_hit.getTime() - 1.0/300*CluHitPos.Mag(); if(tmpTime < T0 && tmpTime > 0) { - T0 = tmpTime; + T0 = tmpTime; } } return T0; @@ -1591,19 +1864,21 @@ int ArborToolLCIO::newPhotonTag(edm4hep::Cluster a_clu) int flag=0; TVector3 PosClu = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); - float Depth = 0; + float Depth = 0; Depth = DisSeedSurface(PosClu); int CluSize= a_clu.hits_size(); float ClusFD=FDV3(a_clu); float ClusT0=ClusterT0(a_clu); - + // cout<<"[YX debug - newPhotonTag] CluSize = "<<CluSize<<", ClusFD = "<<ClusFD<<", ClusT0 = "<<ClusT0<<", Depth = "<<Depth<<endl; + + if(ClusFD>0.18*TMath::Log((float)CluSize)-0.53&&ClusFD<0.16*TMath::Log((float)CluSize)+0.025&&ClusFD>-0.2*TMath::Log((float)CluSize)+0.4&&((log10(ClusT0)<-2&&log10(Depth)<2&&log10(CluSize)>2)||(log10(ClusT0)<-1.5&&log10(CluSize)<2))) { flag=1; } - + return flag; @@ -1619,7 +1894,7 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) float avEnDisHtoL; float EcalEn, HcalEn, EClu, cluDepth, maxDepth, minDepth, FD_all, FD_ECAL, FD_HCAL; float FD_ECALF10; - + TVector3 CluPos; CluPos = TVector3(a_tree.getPosition().x,a_tree.getPosition().y,a_tree.getPosition().z); @@ -1686,7 +1961,18 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) auto a_hit = a_tree.getHits(s1); allhits.push_back(a_hit); auto cellid = a_hit.getCellID(); - int NLayer = m_decoder->get(cellid, "layer"); + // int NLayer = m_decoder->get(cellid, "layer"); + int NLayer; + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + NLayer=cellIdDecoder(&a_hit)[layerCoding.c_str()]+1; + } + + else{ + NLayer = m_decoder->get(cellid, "layer"); + } float tmpHitEn = a_hit.getEnergy(); HitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); @@ -1762,7 +2048,7 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) bool cute3; if (x < 0.9) cute3 = FD_ECALF10 > 0.3/sqrt(sqrt(0.9))*sqrt(sqrt(0.9-x)); else cute3 = 1; - bool cute4 = HcalNHit/EClu < 0.3; + bool cute4 = HcalNHit/EClu < 0.3; bool cute; cute = cute1 && cute2 && cute3 && cute4; @@ -1777,7 +2063,7 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) if(currCluNHits <= 4) ClusterID = 1; else if(cute) ClusterID = 11; else if (cutmu) ClusterID = 13; - else + else { bool cutef1 = FD_HCAL == -1; bool cutef2 = minDepth < 50; @@ -1793,7 +2079,7 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) bool cutef; cutef = (cutef1 && cutef2 && cutef3 && cutef3b && cutef4) || cutef5; - if(cutef) ClusterID = 11; + if(cutef) ClusterID = 11; } if(ClusterID == 211) @@ -1804,7 +2090,7 @@ int ArborToolLCIO::ClusterFlag1st(edm4hep::Cluster a_tree) bool cutmuf = cutmuf1 || (cutmuf2 && cutmuf3); if(cutmuf) ClusterID = 13; else if(minDepth < 0.77+0.23*EClu) ClusterID = 211; - else ClusterID = 2; + else ClusterID = 2; } return ClusterID; @@ -1849,11 +2135,11 @@ float ArborToolLCIO::ClusterEE(edm4hep::Cluster inputCluster) tmpCluEn += hitEn; } } - + EnCorrector = 1; if(tmpCluEn > 1.5 && tmpCluEn < 22 && 1 == 0) { - EnCorrector = 0.6*(1 + 1.0/log10(tmpCluEn)); + EnCorrector = 0.6*(1 + 1.0/log10(tmpCluEn)); } ClusterEnergy = tmpCluEn*EnCorrector; //ClusterEnergy = HADClusterEE(tmpCluEn,inputCluster); @@ -1876,17 +2162,29 @@ float ArborToolLCIO::EMClusterEE( edm4hep::Cluster inputCluster ) float _costheta = 0; float Ethetacorr = 1; float Ephicorr = 1; - float EMC = inputCluster.getEnergy(); - TVector3 CluPos = TVector3(inputCluster.getPosition().x,inputCluster.getPosition().y,inputCluster.getPosition().z); + float EMC = inputCluster.getEnergy(); + TVector3 CluPos = TVector3(inputCluster.getPosition().x,inputCluster.getPosition().y,inputCluster.getPosition().z); float CluTheta = CluPos.Theta(); - float CluPhi = CluPos.Phi()*5.72957795130823229e+01; - + float CluPhi = CluPos.Phi()*5.72957795130823229e+01; + int NCluHits = inputCluster.hits_size(); for(int s1 = 0; s1 < NCluHits; s1++) { auto a_hit = inputCluster.getHits(s1); auto cellid=a_hit.getCellID(); - int NLayer = m_decoder->get(cellid, "layer"); + // int NLayer = m_decoder->get(cellid, "layer"); + + int NLayer; + if(_USE_LCIO){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string idCodingString(m_encoder_str); + const std::string layerCoding(GetLayerCoding(idCodingString)); + NLayer=cellIdDecoder(&a_hit)[layerCoding.c_str()]+1; + } + + else{ + NLayer = m_decoder->get(cellid, "layer"); + } if(NLayer > 20){ if (NLayer%2 ==0){ ArNhite10 ++; @@ -1953,16 +2251,16 @@ float ArborToolLCIO::EMClusterEE( edm4hep::Cluster inputCluster ) std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::Cluster inputCluster) { - std::vector<float> CluTimeVector; + std::vector<float> CluTimeVector; CluTimeVector.clear(); int NHit = inputCluster.hits_size(); float currhitTime = 0; - float currhitoriTime = 0; - TVector3 hitpos; + float currhitoriTime = 0; + TVector3 hitpos; std::vector<float> Time; - Time.clear(); + Time.clear(); std::map<int, float> CluHitToTime; CluHitToTime.clear(); std::vector<float> OriginalTime; @@ -1976,8 +2274,8 @@ std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::Cluster inputCluster) hitpos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); if(a_hit.getTime() == 0) { - currhitTime = 1001; - currhitoriTime = 1001; + currhitTime = 1001; + currhitoriTime = 1001; } else { @@ -1987,17 +2285,17 @@ std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::Cluster inputCluster) } Time.push_back( currhitTime ); OriginalTime.push_back( currhitoriTime ); - CluHitToTime[i] = currhitTime; + CluHitToTime[i] = currhitTime; } std::vector<int> TimeOrder = SortMeasure(Time, 0); std::vector<int> OriTimeOrder = SortMeasure(OriginalTime, 0); float PeakTime = 0; - float AverageTime = 0; - int NCount = 0; - float ptime = 0; - int Break = 0; - TVector3 StartP, EndP, hitP; + float AverageTime = 0; + int NCount = 0; + float ptime = 0; + int Break = 0; + TVector3 StartP, EndP, hitP; for(int j0= 0; j0 < N_PropTimeHit; j0++) { @@ -2037,22 +2335,22 @@ std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::Cluster inputCluster) { if( Break == 0 ) { - PeakTime += ptime; - NCount ++; + PeakTime += ptime; + NCount ++; } - } + } else { Break = 1; } } - AverageTime += ptime; + AverageTime += ptime; } if(NCount) - PeakTime = float(PeakTime)/NCount; + PeakTime = float(PeakTime)/NCount; - CluTimeVector.push_back(float(AverageTime)/NHit); + CluTimeVector.push_back(float(AverageTime)/NHit); CluTimeVector.push_back(PeakTime); CluTimeVector.push_back(NCount); if( N_PropTimeHit > 1 ) // Direction: Cos Theta Between Time Flow And Position @@ -2060,9 +2358,65 @@ std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::Cluster inputCluster) else CluTimeVector.push_back(1001); - return CluTimeVector; + return CluTimeVector; } +std::string ArborToolLCIO::GetStaveCoding(const std::string &encodingString) const +{ + if (encodingString.find("stave") != std::string::npos) + return std::string("stave"); + + if (encodingString.find("S-1") != std::string::npos) + return std::string("S-1"); + + if (encodingString.find("S") != std::string::npos) + return std::string("S"); + + return std::string("unknown_stave_encoding"); +} + +std::string ArborToolLCIO::GetLayerCoding(const std::string &encodingString) const +{ + if (encodingString.find("layer") != std::string::npos) + return std::string("layer"); + + if (encodingString.find("K-1") != std::string::npos) + return std::string("K-1"); + + if (encodingString.find("K") != std::string::npos) + return std::string("K"); + + return std::string("unknown_layer_encoding"); +} + + +std::string ArborToolLCIO::GetCellICoding(const std::string &encodingString) const +{ + if (encodingString.find("cell_x") != std::string::npos) + return std::string("cell_x"); + + if (encodingString.find("I-1") != std::string::npos) + return std::string("I-1"); + + if (encodingString.find("I") != std::string::npos) + return std::string("I"); + + return std::string("unknown_cellI_encoding"); +} + +std::string ArborToolLCIO::GetCellJCoding(const std::string &encodingString) const +{ + if (encodingString.find("cell_y") != std::string::npos) + return std::string("cell_y"); + + if (encodingString.find("J-1") != std::string::npos) + return std::string("J-1"); + + if (encodingString.find("J") != std::string::npos) + return std::string("J"); + + return std::string("unknown_cellJ_encoding"); +} diff --git a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh index 8a5a13d9fb3c9e9361f6f4c22613d3d298aae221..b61955b46336a3556c4c3000288f514e48bfb0d3 100644 --- a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh +++ b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh @@ -12,7 +12,7 @@ #include "edm4hep/CalorimeterHitCollection.h" #include "edm4hep/VertexCollection.h" #include "edm4hep/TrackCollection.h" -#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticle.h" #include "edm4hep/MCParticleCollection.h" #include "edm4hep/MCRecoCaloAssociation.h" #include "edm4hep/MCRecoTrackerAssociation.h" @@ -51,81 +51,89 @@ public: class ArborToolLCIO : public GaudiAlgorithm { public: - ArborToolLCIO(const std::string& name, ISvcLocator* svcLoc); + ArborToolLCIO(const std::string& name, ISvcLocator* svcLoc, bool m_readLCIO); //ArborToolLCIO(ISvcLocator* svcLoc); ~ArborToolLCIO(); typedef DataHandle<edm4hep::Cluster> ClusterType; + // bool _USE_LCIO; - void Clean(); + + void Clean(); void ClusterBuilding(DataHandle<edm4hep::ClusterCollection>& _currbranchcoll, std::vector<edm4hep::CalorimeterHit> Hits, branchcoll BranchOrder, int DHCALFlag); - + float ClusterT0(edm4hep::Cluster a_Clu); int NHScaleV2( std::vector<edm4hep::CalorimeterHit> clu0, int RatioX, int RatioY, int RatioZ ); float FDV2( std::vector<edm4hep::CalorimeterHit> clu); - + int NHScaleV3( edm4hep::Cluster clu0, int RatioX, int RatioY, int RatioZ ); - + float FDV3( edm4hep::Cluster clu); - + int ActiveLayers( std::vector<edm4hep::CalorimeterHit> clu ); - + float BushDis( edm4hep::Cluster clu1, edm4hep::Cluster clu2); - - + + float* SimpleDisTrackClu(edm4hep::Track a_trk, edm4hep::Cluster a_clu); - + float SimpleBushTimeTrackClu(edm4hep::Track a_trk, edm4hep::Cluster a_clu); - + int SimpleBushNC(edm4hep::Track a_trk, edm4hep::Cluster a_clu); - + float DisPointToBush( TVector3 Pos1, edm4hep::Cluster clu1); - + TVector3 ClusterCoG(edm4hep::Cluster inputCluser); - + edm4hep::ClusterCollection* ClusterVecMerge( std::vector<edm4hep::Cluster> inputClusters, TMatrixF ConnectorMatrix, DataHandle<edm4hep::ClusterCollection>& clucol ); - + edm4hep::ClusterCollection* ClusterVecColl( std::vector<edm4hep::MutableCluster> inputClusters, DataHandle<edm4hep::ClusterCollection>& m_clucol); edm4hep::Cluster NaiveCluImpl(edm4hep::MutableCluster a0_clu); void NaiveCluConst(edm4hep::MutableCluster a0_clu, edm4hep::MutableCluster b0_clu); - + std::vector<edm4hep::Cluster> CollClusterVec(const edm4hep::ClusterCollection * input_coll ); - + std::vector<edm4hep::CalorimeterHit> CollHitVec(const edm4hep::CalorimeterHitCollection * input_coll, float EnergyThreshold); - + std::vector<edm4hep::MutableCluster> ClusterHitAbsorbtion( std::vector<edm4hep::Cluster> MainClusters, std::vector<edm4hep::CalorimeterHit> IsoHits, float DisThreshold ); - + edm4hep::MutableCluster NaiveMergeClu(std::vector<edm4hep::Cluster> inputCluVec); void NaiveMergeCluConst(std::vector<edm4hep::Cluster> inputCluVec,edm4hep::MutableCluster MergedClu); std::vector<edm4hep::MutableCluster> ClusterAbsorbtion( std::vector<edm4hep::Cluster> MainClusters, std::vector<edm4hep::MutableCluster> FragClusters, float thresholds, float DepthSlope ); - - + + int JointsBetweenBush(edm4hep::Cluster a_Clu, edm4hep::Cluster b_Clu, float CellSize); - + TVector3 CluEndP(edm4hep::Cluster a_clu); - + int ClusterFlag(edm4hep::Cluster a_tree, edm4hep::Track a_trk); - + int ClusterFlag1st(edm4hep::Cluster a_tree); - + int newPhotonTag(edm4hep::Cluster a_clu); float ClusterEE(edm4hep::Cluster inputCluster); - + float EMClusterEE( edm4hep::Cluster inputCluster ); - + std::vector<float> ClusterTime(edm4hep::Cluster inputCluster); + std::string GetCellICoding(const std::string &encodingString) const; + std::string GetCellJCoding(const std::string &encodingString) const; + std::string GetStaveCoding(const std::string &encodingString) const; + std::string GetLayerCoding(const std::string &encodingString) const; + private: SmartIF<IGeomSvc> m_geosvc;//=service<IGeomSvc>("GeomSvc"); dd4hep::DDSegmentation::BitFieldCoder* m_decoder;// = m_geosvc->getDecoder("ECALBarrel"); + std::string m_encoder_str; // m_geosvc= service<IGeomSvc>("GeomSvc"); // m_decoder = m_geosvc->getDecoder("ECALBarrel"); - + bool _USE_LCIO; }; #endif // diff --git a/Reconstruction/PFA/Arbor/src/BushConnect.cc b/Reconstruction/PFA/Arbor/src/BushConnect.cc index 052f49db319d300059ca6da58751798102b4c557..5ba4be966f37ac3d2946ef0c6342ff2e8be5fe0b 100644 --- a/Reconstruction/PFA/Arbor/src/BushConnect.cc +++ b/Reconstruction/PFA/Arbor/src/BushConnect.cc @@ -2,6 +2,7 @@ #include "ArborTool.h" #include "ArborToolLCIO.hh" #include "DetectorPos.hh" +#include "HelixClassD.hh" #include "k4FWCore/DataHandle.h" #include "GaudiAlg/GaudiAlgorithm.h" @@ -45,7 +46,6 @@ #include <vector> #include <algorithm> -#include "HelixClassD.hh" using namespace std; @@ -63,13 +63,16 @@ BushConnect::BushConnect(const std::string& name, ISvcLocator* svcLoc) StatusCode BushConnect::initialize() { ISvcLocator* svcloc = serviceLocator(); - m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc); + m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc, m_readLCIO); + // m_ArborToolLCIO->_USE_LCIO = m_readLCIO; //printParameters(); //Cluflag.setBit(LCIO::CHBIT_LONG); return GaudiAlgorithm::initialize(); } void BushConnect::Clean(){ + cout<<"[YX debug] BushConnect::Clean Begin:"<<endl; + MCPTrack_Type.clear(); Track_Energy.clear(); Track_Type.clear(); @@ -87,38 +90,43 @@ void BushConnect::Clean(){ CluT0.clear(); CluCoG.clear(); Clu_Depth.clear(); + + cout<<"[YX debug] BushConnect::Clean End."<<endl; + } void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Track*, float> Track_Energy) { + // cout<<"[YX debug] TrackSort Begin:"<<endl; try{ auto TrackColl = m_trkcol.get(); int NTrk = TrackColl->size(); - cout<<NTrk<<endl; + // cout<<"[YX debug] NTrk = "<<NTrk<<endl; + float D0 = 0; float Z0 = 0; int NTrkHit = 0; const float mass = 0.139; //Pion Mass - TVector3 EndPointPos, StartPointPos; - int TrackType = 0; + TVector3 EndPointPos, StartPointPos; + int TrackType = 0; - std::vector<edm4hep::Track> tracks_HQ_Barrel; + std::vector<edm4hep::Track> tracks_HQ_Barrel; std::vector<edm4hep::Track> tracks_HQ_Endcap; std::vector<edm4hep::Track> tracks_HQ_Shoulder; - std::vector<edm4hep::Track> tracks_HQ_Forward; + std::vector<edm4hep::Track> tracks_HQ_Forward; std::vector<edm4hep::Track> tracks_MQ_Barrel; std::vector<edm4hep::Track> tracks_MQ_Endcap; std::vector<edm4hep::Track> tracks_MQ_Shoulder; std::vector<edm4hep::Track> tracks_MQ_Forward; - std::vector<edm4hep::Track> tracks_Vtx; - std::vector<edm4hep::Track> tracks_LQ; - std::vector<edm4hep::Track> tracks_LE; + std::vector<edm4hep::Track> tracks_Vtx; + std::vector<edm4hep::Track> tracks_LQ; + std::vector<edm4hep::Track> tracks_LE; std::vector<edm4hep::Track> curr_tracks; - + Track_EndPoint.clear(); - + tracks_HQ_Barrel.clear(); tracks_HQ_Endcap.clear(); tracks_HQ_Shoulder.clear(); @@ -137,57 +145,71 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr tracks_preInteraction.clear(); //Used to denote pion and electron interaction inside TPC/Tracker. Simply vetoed for avoid double counting... but muon may still be problematic. Better way of treating would be find the cascade photons & tracks - clusters, and veto all the daughters instead of mother. Similar can done for Kshort... // Condition, tracks_head to others tail. head position far from boundary. and, track energy >= sum of cascade - std::vector<int> TrackOrder; - TrackOrder.clear(); - std::map<edm4hep::Track, int> Track_Index; + std::vector<int> TrackOrder; + TrackOrder.clear(); + std::map<edm4hep::Track, int> Track_Index; Track_Index.clear(); Track_Energy.clear(); Track_Type.clear(); Track_P3.clear(); Track_EndPoint.clear(); TrackStartPoint.clear(); - + + // cout<<"[YX debug] i0 for loop begin."<<endl; + for(int i0 = 0; i0 < NTrk; i0++) { auto a_Trk = (*TrackColl)[i0]; - NTrkHit = a_Trk.getTrackerHits().size(); - EndPointPos = TVector3((a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().z); + NTrkHit = a_Trk.getTrackerHits().size(); + + // cout<<"[YX debug] Trk "<<i0<<", NTrkHit = "<<NTrkHit<<", a_Trk = "<<a_Trk<<endl; + + // if(NTrkHit==0) continue; + + EndPointPos = TVector3((a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(a_Trk.getTrackerHits(NTrkHit - 1)).getPosition().z); StartPointPos = TVector3((a_Trk.getTrackerHits(0)).getPosition().x,(a_Trk.getTrackerHits(0)).getPosition().y,(a_Trk.getTrackerHits(0)).getPosition().z); Track_EndPoint[a_Trk] = EndPointPos; TrackStartPoint[a_Trk] = StartPointPos; - + HelixClassD * TrkInit_Helix = new HelixClassD(); TrkInit_Helix->Initialize_Canonical(a_Trk.getTrackStates(0).phi, a_Trk.getTrackStates(0).D0, a_Trk.getTrackStates(0).Z0, a_Trk.getTrackStates(0).omega, a_Trk.getTrackStates(0).tanLambda, BField); float TrackEn = mass*mass; - + for (int q3 = 0; q3 < 3; q3 ++) { TrackEn += (TrkInit_Helix->getMomentum()[q3])*(TrkInit_Helix->getMomentum()[q3]); } TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]); - + TrackEn = sqrt(TrackEn); Track_Energy[a_Trk] = TrackEn; Track_Theta[a_Trk] = TrkMom.Theta(); // Track_Phi[a_Trk] = TrkMom.Phi(); - Track_P3[a_Trk] = TrkMom; - + Track_P3[a_Trk] = TrkMom; + + // cout<<"[YX debug - TrackSort] Input Trk "<<i0<<", PMag = "<<TrkMom.Mag()<<", NTrkHit = "<<NTrkHit<<endl; // 能对上 + + delete TrkInit_Helix; } - + // cout<<"[YX debug] i0 for loop end."<<endl; + TVector3 currEp, currSp; float currMotherEn = 0; - float sumDauEn = 0; - + float sumDauEn = 0; + for(int i1 = 0; i1 < NTrk; i1++) { auto a_Trk = (*TrackColl)[i1]; currEp = Track_EndPoint[a_Trk]; - - if( currEp.Perp() < 1600 && currEp.Perp() > 400 && abs(currEp.Z()) < 2000 ) //Only check + + // cout<<"[YX debug - TrackSort] if_preInteraction Trk "<<i1<<", PMag = "<<Track_P3[a_Trk].Mag()<<", NTrkHit = "<<a_Trk.getTrackerHits().size()<<endl; + + + if( currEp.Perp() < 1600 && currEp.Perp() > 400 && abs(currEp.Z()) < 2000 ) //Only check { currMotherEn = Track_Energy[a_Trk]; - sumDauEn = 0; + sumDauEn = 0; for(int i2 = 0; i2 < NTrk; i2++) { auto b_Trk = (*TrackColl)[i2]; @@ -201,10 +223,18 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr if(currMotherEn + 0.1 > 0.9*sumDauEn && currMotherEn > 3 && sumDauEn > 0 ) //Some protection is always needed... { tracks_preInteraction.push_back(a_Trk); + + // cout<<"[YX debug - TrackSort] ---> is preInteraction push back!!!"<<endl; + } } } - + + cout<<"[YX debug - TrackSort] tracks_preInteraction.size() = "<<tracks_preInteraction.size()<<endl; + + cout<<"[YX debug - TrackSort] Geo TPCInnerRadius = "<<TPCInnerRadius<<", TPCOuterRadius = "<<TPCOuterRadius<<", ECALHalfZ = "<<ECALHalfZ<<", LStar = "<<LStar<<endl; + // cout<<"[YX debug] i1 for loop end."<<endl; + for(int t0 = 0; t0 < NTrk; t0++) { auto a_Trk = (*TrackColl)[t0]; @@ -216,18 +246,18 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr Track_EndPoint[a_Trk] = EndPointPos; StartPointPos = TVector3((a_Trk.getTrackerHits(0)).getPosition().x,(a_Trk.getTrackerHits(0)).getPosition().y,(a_Trk.getTrackerHits(0)).getPosition().z); Track_Index[a_Trk] = t0; - + if( NTrkHit > 9 || (fabs(EndPointPos.Z()) > LStar - 500 && EndPointPos.Perp() < TPCInnerRadius ) || fabs(EndPointPos.Z()) > ECALHalfZ - 200 ) // Min requirement for track quality { // LStar - 500, suppose to be the last Disk Position - + if( find(tracks_preInteraction.begin(), tracks_preInteraction.end(), a_Trk ) != tracks_preInteraction.end() ) { - cout<<"So We Drop it! "<<Track_Energy[a_Trk]<<endl; - continue; + cout<<"So We Drop it! "<<Track_Energy[a_Trk]<<endl; + continue; } - + TrackType = 0; - + if((Track_Energy[a_Trk] < 1.0 && fabs(Track_Theta[a_Trk]-1.57)< 0.4) || (fabs(Track_Theta[a_Trk]-1.57) >= 0.4 && log10(Track_Energy[a_Trk]) < -(fabs(Track_Theta[a_Trk]-1.57)-0.4)*0.2/0.3 )) { TrackType = 100; @@ -246,16 +276,16 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr } else if( fabs(EndPointPos.Z()) > ECALHalfZ - 200 ) //Endcap { - TrackType = 20; + TrackType = 20; } - + if( fabs(D0) < 1 && fabs(Z0) < 1 ) { TrackType += 1; } - - Track_Type[a_Trk] = TrackType; - + + Track_Type[a_Trk] = TrackType; + if(TrackType == 11) tracks_HQ_Barrel.push_back(a_Trk); else if(TrackType == 21) @@ -282,12 +312,18 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr else tracks_LQ.push_back(a_Trk); } + + + // cout<<"[YX debug - TrackSort] TrackType Trk "<<t0<<", TrackType = "<<TrackType<<", NTrkHit = "<<a_Trk.getTrackerHits().size()<<", PMag = "<<Track_P3[a_Trk].Mag()<<", E = "<<Track_Energy[a_Trk]<<", Theta = "<<Track_Theta[a_Trk]<<", EndPointPos (Perp, Z) = ("<<EndPointPos.Perp()<<", "<<EndPointPos.Z()<<"), D0 = "<<D0<<", Z0 = "<<Z0<<endl; // åªæ˜¯æœ‰äº› TrackType 对ä¸ä¸Šï¼Œä½†å¥½åƒä¸æ˜¯åŽé¢ TagCore 里对ä¸ä¸Šçš„ + } - cout<<"LQ"<<tracks_LQ.size()<<" "<<tracks_ILL.size()<<endl; - + if(1) cout<<"LQ "<<tracks_LQ.size()<<" "<<tracks_ILL.size()<<endl; + + // cout<<"[YX debug] t0 for loop end."<<endl; + std::vector<float > currTrkMomentum; std::vector<int> currTrkIndex; - + for(int t1 = 0; t1 < 11; t1++) { currTrkMomentum.clear(); @@ -311,22 +347,22 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr curr_tracks = tracks_MQ_Forward; else if(t1 == 8) curr_tracks = tracks_Vtx; - else if(t1 == 9) - curr_tracks = tracks_LQ; - else if(t1 == 10) - curr_tracks = tracks_LE; - - + else if(t1 == 9) + curr_tracks = tracks_LQ; + else if(t1 == 10) + curr_tracks = tracks_LE; + + int N_currTrack = curr_tracks.size(); - + for(int t2 = 0; t2 < N_currTrack; t2++) { auto tmpTrk = curr_tracks[t2]; currTrkMomentum.push_back(Track_Energy[tmpTrk]); } - + currTrkIndex = SortMeasure(currTrkMomentum, 1); - + for(int t3 = 0; t3 < N_currTrack; t3++) { auto b_tmpTrk = curr_tracks[currTrkIndex[t3]]; @@ -334,40 +370,54 @@ void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Tr TrackOrder.push_back(Track_Index[b_tmpTrk]); } } - + // cout<<"[YX debug] 11 for loop end."<<endl; + + for(unsigned int t4 = 0; t4 < TrackOrder.size(); t4++) { - auto b_Trk = (*TrackColl)[t4]; + // auto b_Trk = (*TrackColl)[t4]; /// !!! + auto b_Trk = (*TrackColl)[TrackOrder[t4]]; SortedTracks.push_back(b_Trk); } + + // cout<<"[YX debug] TrackOrder for loop end."<<endl; + }catch(GaudiException &e){} + + // cout<<"[YX debug] TrackSort End."<<endl; + } void BushConnect::BushSelfMerge() { + // cout<<"[YX debug] BushSelfMerge Begin:"<<endl; + auto CaloClu = m_clucol.get(); int NClu = CaloClu->size(); - std::vector<edm4hep::Cluster > Core_1st; + cout<<"[YX debug - BushSelfMerge] Input EHBushes nClu = "<<NClu<<endl; + + std::vector<edm4hep::Cluster > Core_1st; std::vector<edm4hep::MutableCluster > Frag_1st; - std::vector<edm4hep::Cluster > UnId_1st; + std::vector<edm4hep::Cluster > UnId_1st; Core_1st.clear(); Frag_1st.clear(); UnId_1st.clear(); - float CluDepth = 0; + float CluDepth = 0; float CluEn = 0; - int CluSize = 0; - TVector3 PosCluSeed, PosSeedDiff, PosCoGDiff, PosSeedA, PosSeedB; + int CluSize = 0; + TVector3 PosCluSeed, PosSeedDiff, PosCoGDiff, PosSeedA, PosSeedB; - int NJoints = 0; - int SmallCluSize = 0; float DeeperDepth = 0; float LaterT0 = 0; + int NJoints = 0; + int SmallCluSize = 0; float DeeperDepth = 0; float LaterT0 = 0; float Depth_A = 0; float Depth_B = 0; - int Size_A = 0; int Size_B = 0; + int Size_A = 0; int Size_B = 0; TMatrixF FlagMerge(NClu, NClu); - cout<<NClu<<" clusters"<<endl; + double ESum_Clu_Input = 0; + // cout<<NClu<<" clusters"<<endl; for(int i0 = 0; i0 < NClu; i0++) { auto a_clu = (*CaloClu)[i0]; @@ -376,8 +426,15 @@ void BushConnect::BushSelfMerge() CluCoG[a_clu] = m_ArborToolLCIO->ClusterCoG(a_clu); PosCluSeed = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); Clu_Depth[a_clu] = DisSeedSurface(PosCluSeed); + + ESum_Clu_Input += a_clu.getEnergy(); + // cout<<"[YX debug - BushSelfMerge] Input EHBushes Clu "<<i0<<", E = "<<a_clu.getEnergy()<<", Pos = ("<<PosCluSeed.X()<<", "<<PosCluSeed.Y()<<", "<<PosCluSeed.Z()<<")"<<endl; + } // CluCoG_Top20Percent Might be used to improve Photon Split Remerge performance + cout<<"[YX debug - BushSelfMerge] Input EHBushes ESum_Clu_Input = "<<ESum_Clu_Input<<endl; + + for(int s0 = 0; s0 < NClu; s0++) { @@ -407,7 +464,7 @@ void BushConnect::BushSelfMerge() if( SmallCluSize < 10 || LaterT0 > 10 || DeeperDepth > 40 || (NJoints > 8 && PosCoGDiff.Mag()*PosCoGDiff.Angle(PosSeedB) < 15 ) || NJoints>16) //if( SmallCluSize < 10 || LaterT0 > 10 || DeeperDepth > 40 || (NJoints > 8 ) ) - { + { FlagMerge[s0][s1] = 1.0; FlagMerge[s1][s0] = 1.0; } @@ -440,8 +497,8 @@ void BushConnect::BushSelfMerge() auto b_clu = (*CloseMergedCaloClu)[i1]; if(i1 != i0) { - if(m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu) < tmpmindis) - tmpmindis = m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu); + if(m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu) < tmpmindis) + tmpmindis = m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu); } } MinDisSeedToBush[a_clu] = tmpmindis; @@ -476,10 +533,17 @@ void BushConnect::BushSelfMerge() std::vector<edm4hep::MutableCluster > UndefFrag_1stAB = m_ArborToolLCIO->ClusterAbsorbtion(UnId_1st, Frag_1st, 50, 0.02); selfmergedcluster = m_ArborToolLCIO->ClusterAbsorbtion(Core_1st, UndefFrag_1stAB, 50, 0.02); auto CluAB_1st=m_ArborToolLCIO->ClusterVecColl(selfmergedcluster,m_1stclucol); + + cout<<"[YX debug - BushSelfMerge] Output CluAB_1st nClu = "<<selfmergedcluster.size()<<endl; + + // cout<<"[YX debug] BushSelfMerge End."<<endl; + } -void BushConnect::TagCore() +void BushConnect::TagCore() { + // cout<<"[YX debug] TagCore Begin:"<<endl; + int NTrk = SortedTracks.size(); int NClu = selfmergedcluster.size(); int currTrackType = 0; @@ -489,9 +553,13 @@ void BushConnect::TagCore() float CluDepth = 0; float CoreMergeDistanceDepthCorrector = 0; - TVector3 TrkEndPoint(0, 0, 0); + + cout<<"[YX debug - TagCore] Input NTrk = "<<NTrk<<", nClu = "<<NClu<<endl; + + + TVector3 TrkEndPoint(0, 0, 0); TVector3 CluPos(0, 0, 0); - std::map<edm4hep::Cluster, int> BushTouchFlag; + std::map<edm4hep::Cluster, int> BushTouchFlag; std::map<edm4hep::Track, edm4hep::MutableCluster> FCMap_Track_CHCore; std::map<edm4hep::Track, std::vector<edm4hep::Cluster>> FCMap_Track_CHCore_new; std::map<int, int> Closest_Trk_Clu_Map; @@ -509,24 +577,31 @@ void BushConnect::TagCore() for(int s1 = 0; s1 < NClu; s1++) { DisMatrix_Track_Clu_E[s0][s1] = 1.0E10; - TimeMatrix_Track_Clu_E[s0][s1] = 1.0E10; + TimeMatrix_Track_Clu_E[s0][s1] = 1.0E10; } } + double ESum_Clu_Input =0; for(int t0 = 0; t0 < NClu; t0++) { auto a_clu = selfmergedcluster[t0]; TVector3 PosCluSeed = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); Clu_Depth[a_clu] = DisSeedSurface(PosCluSeed); + + ESum_Clu_Input += a_clu.getEnergy(); + // cout<<"[YX debug - TagCore] Input AB_1st Clu "<<t0<<", E = "<<a_clu.getEnergy()<<", Pos = ("<<a_clu.getPosition().x<<", "<<a_clu.getPosition().y<<", "<<a_clu.getPosition().z<<")"<<endl; + } + cout<<"[YX debug - TagCore] Input AB_1st ESum_Clu_Input = "<<ESum_Clu_Input<<endl; + //~~~~~~~ find the closest cluster first... for(int g0 = 0; g0 < NTrk; g0++) { auto a_trk = SortedTracks[g0]; - float ClosestDis = 1.0E9; - int ClosestCluIndex = -1; + float ClosestDis = 1.0E9; + int ClosestCluIndex = -1; int ClosestNC = 1E9; float TrackEn = Track_Energy[a_trk]; TrkEndPoint = Track_EndPoint[a_trk]; @@ -534,6 +609,8 @@ void BushConnect::TagCore() currTrackType = Track_Type[a_trk]; TVector3 TrkP3 = Track_P3[a_trk]; + // cout<<"[YX debug - TagCore] SortedTrack "<<g0<<", PMag = "<<TrkP3.Mag()<<", TrackEn = "<<TrackEn<<endl; + for(int g1 = 0; g1 < NClu; g1++) { auto fccand_bush = selfmergedcluster[g1]; @@ -546,7 +623,7 @@ void BushConnect::TagCore() TimeMatrix_Track_Clu_E[g0][g1] = Time; if( Dis[2] < ClosestDis ) // && ThetaDiff < 0.05 + 0.1/Track_Energy[a_trk] ) { - ClosestDis = Dis[2]; + ClosestDis = Dis[2]; ClosestCluIndex = g1; ClosestNC = NC; } @@ -554,10 +631,10 @@ void BushConnect::TagCore() } //Diag for mimic - cout<<" Z R "<<TrkEndPoint.Z()<<" MM "<<TrkEndPoint.Perp()<< " ClosestCluIndex "<<ClosestCluIndex<<" ClosestDis "<<ClosestDis <<endl; + if(0) cout<<" Z R "<<TrkEndPoint.Z()<<" MM "<<TrkEndPoint.Perp()<< " ClosestCluIndex "<<ClosestCluIndex<<" ClosestDis "<<ClosestDis <<endl; //End Diag - if( ClosestDis < 15 + 15./TrackEn && ClosestCluIndex > -0.1 && (ClosestNC < 3 || abs(TrkP3.Theta() - 1.57) < 0.01 ) ) + if( ClosestDis < 15 + 15./TrackEn && ClosestCluIndex > -0.1 && (ClosestNC < 3 || abs(TrkP3.Theta() - 1.57) < 0.01 ) ) { auto candiclu= selfmergedcluster[ClosestCluIndex]; CluPos = TVector3(candiclu.getPosition().x,candiclu.getPosition().y,candiclu.getPosition().z); @@ -570,7 +647,7 @@ void BushConnect::TagCore() } } //~~~~~~~ end of finding closest cluster - cout<<"NTrk "<<NTrk<<endl; + // cout<<"NTrk "<<NTrk<<endl; for(int i0 = 0; i0 < NTrk; i0++) //Dropped Size can exist { auto a_trk = SortedTracks[i0]; @@ -587,7 +664,7 @@ void BushConnect::TagCore() for(int j0 = 0; j0 < NClu; j0++) { - auto fccand_bush = selfmergedcluster[j0]; + auto fccand_bush = selfmergedcluster[j0]; float Dis = DisMatrix_Track_Clu_E[i0][j0]; //SimpleDisTrackClu(a_trk, fccand_bush); float BushTime = TimeMatrix_Track_Clu_E[i0][j0]; CluDepth = Clu_Depth[fccand_bush]; @@ -670,6 +747,11 @@ void BushConnect::TagCore() edm4hep::ReconstructedParticleCollection* chargeparticleCol = m_chargeparticleCol.createAndPut(); edm4hep::ClusterCollection* chargedcoreclusterCol = m_chargedcoreclusterCol.createAndPut(); + + cout<<"[YX debug - TagCore] NTrk = "<<NTrk<<endl; + double ESum_ChCore_All = 0; + double ESum_ChCore_in = 0; + for(int j5 = 0; j5 < NTrk; j5++) { auto a_trk = SortedTracks[j5]; @@ -680,9 +762,16 @@ void BushConnect::TagCore() chargeparticle.setCharge(a_trk.getTrackStates(0).omega/fabs(a_trk.getTrackStates(0).omega)); TVector3 Ptrack = Track_P3[a_trk]; edm4hep::Vector3f currTrkP = edm4hep::Vector3f( Ptrack.X(), Ptrack.Y(), Ptrack.Z() ); - chargeparticle.setMomentum(currTrkP); + chargeparticle.setMomentum(currTrkP); chargeparticle.addToTracks( a_trk ); auto a_clu = FCMap_Track_CHCore[a_trk]; + + double Clu_ChCore_E = a_clu.getEnergy(); + int nCluHit = a_clu.hits_size(); + if(0) cout<<"[YX debug - TagCore] Trk "<<j5<<", PMag = "<<Ptrack.Mag()<<", CluCore E = "<<Clu_ChCore_E<<", nHit = "<<nCluHit<<endl; + ESum_ChCore_All += Clu_ChCore_E; + + if( FCMap_Track_CHCore[a_trk].hits_size()>0 ) // No really need to pertect, as quality will be controled in Part.Reco { auto chargedcorecluster = chargedcoreclusterCol->create(); @@ -691,16 +780,43 @@ void BushConnect::TagCore() edm4hep::Cluster chargedcoreclusterCon=chargedcorecluster; chargeparticle.addToClusters(chargedcoreclusterCon); Track_Core_ID = m_ArborToolLCIO->ClusterFlag(a_clu, a_trk); + + if(0) cout<<"[YX debug - TagCore] ---> ChCore in"<<endl; + ESum_ChCore_in += Clu_ChCore_E; + } chargeparticle.setType(Track_Core_ID); ChCoreID[chargeparticle] = Track_Core_ID; } +// auto ChCoreCluCol = m_chargedcoreclusterCol.get(); + + + + int nPFO_ChCore = chargeparticleCol->size(); + int nClu_ChCore = chargedcoreclusterCol->size(); + // int nClu_NeCore = arborneutralcorecluster.size(); + // int nPFO_NeCore = arborrecoparticle_ne.size(); + + int nClu_NonChCore = non_chargedclustercore.size(); + + cout<<"[YX debug - TagCore] Output nPFO_ChCore = "<<nPFO_ChCore<<endl; + cout<<"[YX debug - TagCore] Output nClu_ChCore = "<<nClu_ChCore<<", ESum_ChCore_All = "<<ESum_ChCore_All<<", ESum_ChCore_in = "<<ESum_ChCore_in<<endl; + cout<<"[YX debug - TagCore] Output nClu_NonChCore = "<<nClu_NonChCore<<endl; + // cout<<"[YX debug - TagCore] Output nClu_NeCore = "<<nClu_NeCore<<endl; + // cout<<"[YX debug - TagCore] Output nPFO_NeCore = "<<nPFO_NeCore<<endl; + + + + + + // cout<<"[YX debug] TagCore End."<<endl; } void BushConnect::ParticleReco() { + // cout<<"[YX debug] ParticleReco Begin:"<<endl; auto col_IsoHit = m_col_IsoHit.get(); std::vector<edm4hep::CalorimeterHit> IsoHits = m_ArborToolLCIO->CollHitVec(col_IsoHit, 0); @@ -710,30 +826,71 @@ void BushConnect::ParticleReco() edm4hep::ClusterCollection* mergedclu_chCol = m_mergedclu_chCol.createAndPut(); edm4hep::ClusterCollection* mergedclu_neCol = m_mergedclu_neCol.createAndPut(); - auto ChargedCore = m_chargeparticleCol.get(); + auto ChargedCore = m_chargeparticleCol.get(); int NChargedObj = ChargedCore->size(); int NNeutralCluster = non_chargedclustercore.size(); - double DisMatrix_Core_Neutral[NChargedObj][NNeutralCluster][2]; //Define different types of distances; + double DisMatrix_Core_Neutral[NChargedObj][NNeutralCluster][2]; //Define different types of distances; + + + // ------------------------------------------------------------------------------------------ + double ESum_ChCore = 0; + double ESum_NeClu = 0; + for(int i = 0; i < NChargedObj; i++) + { + auto a_recoP_ch = (*ChargedCore)[i]; + auto aTrk = a_recoP_ch.getTracks()[0]; + float Trk_E = Track_Energy[aTrk]; + float Clu_E = 0; + if(a_recoP_ch.clusters_size() != 0) + { + edm4hep::Cluster aClu = a_recoP_ch.getClusters()[0]; + Clu_E = aClu.getEnergy(); + } + ESum_ChCore+=Trk_E; + + // cout<<"[YX debug - ParticleReco] Input Trk "<<i<<", E = "<<Trk_E<<", Clu_E = "<<Clu_E<<endl; + + } + + for(int j = 0; j < NNeutralCluster; j++) + { + auto aClu = non_chargedclustercore[j]; + int Clu_nHit = aClu.hits_size(); + float Clu_E = aClu.getEnergy(); + TVector3 Clu_Pos = TVector3(aClu.getPosition().x,aClu.getPosition().y,aClu.getPosition().z); + double Clu_Depth = DisSeedSurface(Clu_Pos); + + ESum_NeClu+=Clu_E; + + // cout<<"[YX debug - ParticleReco] Input NeClu "<<j<<", E = "<<Clu_E<<", nHit = "<<Clu_nHit<<", Depth = "<<Clu_Depth<<endl; + + } + + cout<<"[YX debug - ParticleReco] NChargedObj = "<<NChargedObj<<", ESum_ChCore = "<<ESum_ChCore<<endl; + cout<<"[YX debug - ParticleReco] NNeutralCluster = "<<NNeutralCluster<<", ESum_NeClu = "<<ESum_NeClu<<endl; + // ------------------------------------------------------------------------------------------ + + float CluDepth = 0; - std::map<edm4hep::Cluster, double> CluDepthMap; + std::map<edm4hep::Cluster, double> CluDepthMap; CluDepthMap.clear(); - int currChargeCoreType = 0; - TVector3 CluPos; + int currChargeCoreType = 0; + TVector3 CluPos; - std::vector<edm4hep::Cluster> loosecandicluster; + std::vector<edm4hep::Cluster> loosecandicluster; std::vector<edm4hep::Cluster> tightcandicluster; //Muon potential candi? std::vector<edm4hep::Cluster> mergedcluster; //tmp for each charged P std::vector<edm4hep::Cluster> chargedclustercore_merged; //overall chargedclustercore_merged.clear(); - std::vector<double> reftightdis; - std::vector<double> refloosedis; - std::map<edm4hep::Cluster, int> NNCTouchFlag; + std::vector<double> reftightdis; + std::vector<double> refloosedis; + std::map<edm4hep::Cluster, int> NNCTouchFlag; std::vector<edm4hep::Track> SecondIterTracks; SecondIterTracks.clear(); - TVector3 currTrkEnd, neighbourTrkEnd, LeadP; + TVector3 currTrkEnd, neighbourTrkEnd, LeadP; for(int i = 0; i < NChargedObj; i++) { @@ -759,7 +916,7 @@ void BushConnect::ParticleReco() mergedcluster.push_back(a_chargedClu); //Actually can use this chance to question if previous energy are balance... } - float MinDisToNoClusterTrk = 1.0E10; + float MinDisToNoClusterTrk = 1.0E10; float MinDisToOtherTrack = 1.0E10; for( int is = 0; is < NChargedObj; is++ ) @@ -780,12 +937,12 @@ void BushConnect::ParticleReco() for(int j = 0; j < NNeutralCluster; j++) { auto a_NeCandiClu = non_chargedclustercore[j]; - float NeCandEn = a_NeCandiClu.getEnergy(); + float NeCandEn = a_NeCandiClu.getEnergy(); CluPos = TVector3(a_NeCandiClu.getPosition().x,a_NeCandiClu.getPosition().y,a_NeCandiClu.getPosition().z); CluDepth = DisSeedSurface(CluPos); - CluDepthMap[a_NeCandiClu] = CluDepth; + CluDepthMap[a_NeCandiClu] = CluDepth; - if( ClusterType_1stID[a_NeCandiClu] == 22) continue; + if( ClusterType_1stID[a_NeCandiClu] == 22) continue; for(int k = 0; k < 2; k++) { @@ -800,7 +957,7 @@ void BushConnect::ParticleReco() DisMatrix_Core_Neutral[i][j][1] = Dis[2]; if( NNCTouchFlag.find(a_NeCandiClu) == NNCTouchFlag.end() && ( currChargeCoreType == 0 || DisMatrix_Core_Neutral[i][j][0] < 1000 ) && currTrkType != 101) - { + { if( currChargeCoreType == 130 ) //Matched Muon, should ignore { if( DisMatrix_Core_Neutral[i][j][1] < 0.2*CluDepth && CluDepth > 200 ) //&& FD? @@ -814,14 +971,14 @@ void BushConnect::ParticleReco() if( DisMatrix_Core_Neutral[i][j][1] < 0.3*CluDepth && CluDepth > 150 ) { tightcandicluster.push_back(a_NeCandiClu); - reftightdis.push_back( DisMatrix_Core_Neutral[i][j][1] ); + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][1] ); } else if( DisMatrix_Core_Neutral[i][j][1] < 0.5*CluDepth && CluDepth > 100 ) { loosecandicluster.push_back(a_NeCandiClu); refloosedis.push_back( DisMatrix_Core_Neutral[i][j][1] ); } - } + } else if( currChargeCoreType == 110 ) // Electron { if( DisMatrix_Core_Neutral[i][j][0] < 0.15*CluDepth + 15 ) @@ -829,7 +986,7 @@ void BushConnect::ParticleReco() tightcandicluster.push_back(a_NeCandiClu); reftightdis.push_back( DisMatrix_Core_Neutral[i][j][0] ); } - } + } else if( currChargeCoreType == 111 ) // look behind... might be pion... { if( DisMatrix_Core_Neutral[i][j][0] < 0.1*CluDepth + 15 && DisMatrix_Core_Neutral[i][j][1] < 0.1*CluDepth + 10 ) //Define Brems Photon region for correct @@ -845,7 +1002,7 @@ void BushConnect::ParticleReco() } } else if( DisMatrix_Core_Neutral[i][j][0] < 0.2*CluDepth + 15 || DisMatrix_Core_Neutral[i][j][1] < 0.2*CluDepth + 15 ) - { + { loosecandicluster.push_back(a_NeCandiClu); if(DisMatrix_Core_Neutral[i][j][0] < DisMatrix_Core_Neutral[i][j][1]) // not fully adequate. @@ -896,13 +1053,13 @@ void BushConnect::ParticleReco() } else { - cout<<"Over balanced/Un matched/defined case: "<<a_recoP_ch.getEnergy()<<" ??? "<<currChargeCoreType<<endl; + cout<<"Over balanced/Un matched/defined case: "<<a_recoP_ch.getEnergy()<<" ??? "<<currChargeCoreType<<endl; } } } - float totaltightcandiEn = 0; - float totalloosecandiEn = 0; + float totaltightcandiEn = 0; + float totalloosecandiEn = 0; for(unsigned int s = 0; s < tightcandicluster.size(); s++) { totaltightcandiEn += tightcandicluster[s].getEnergy(); @@ -954,7 +1111,7 @@ void BushConnect::ParticleReco() for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) { auto a_clu = tightcandicluster[i1]; - if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) { mergedcluster.push_back( a_clu ); CurrClusterEnergy += a_clu.getEnergy(); @@ -967,7 +1124,7 @@ void BushConnect::ParticleReco() { auto a_clu = tightcandicluster[i1]; - if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) { mergedcluster.push_back( a_clu ); CurrClusterEnergy += a_clu.getEnergy(); @@ -982,7 +1139,7 @@ void BushConnect::ParticleReco() mergedcluster.push_back( a_clu ); CurrClusterEnergy += a_clu.getEnergy(); } - } + } } else if( currChargeCoreType == 211 ) // Matched { @@ -990,7 +1147,7 @@ void BushConnect::ParticleReco() { auto a_clu = tightcandicluster[i1]; - if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) { mergedcluster.push_back( a_clu ); CurrClusterEnergy += a_clu.getEnergy(); @@ -1002,7 +1159,7 @@ void BushConnect::ParticleReco() for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) { auto a_clu = tightcandicluster[i1]; - if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) { mergedcluster.push_back( a_clu ); CurrClusterEnergy += a_clu.getEnergy(); @@ -1024,14 +1181,14 @@ void BushConnect::ParticleReco() else if( currChargeCoreType == 0 && reftightdis.size() > 0) { float mindis = 1.0E10; - int minindex = 0; + int minindex = 0; for(unsigned int i1 = 0; i1 < reftightdis.size(); i1 ++) { if(reftightdis[i1] < mindis) { mindis = reftightdis[i1]; - minindex = i1; + minindex = i1; } } @@ -1041,13 +1198,13 @@ void BushConnect::ParticleReco() } else { - cout<<"No_match, currChargeCoreType "<<currChargeCoreType<<endl; + cout<<"No_match, currChargeCoreType "<<currChargeCoreType<<endl; } float CHCluEnergy = 0; for(int is = 0; is < int(mergedcluster.size()); is++) - { + { auto a_TBM_clu = mergedcluster[is]; CHCluEnergy +=a_TBM_clu.getEnergy(); } @@ -1079,16 +1236,16 @@ void BushConnect::ParticleReco() chargedclustercore_merged.push_back(chclustermerged); currChargeCoreType2 = m_ArborToolLCIO->ClusterFlag(chclustermerged, a_chargedTrk); - if( currChargeCoreType2 == 130 || currChargeCoreType2 == 131 ) + if( currChargeCoreType2 == 130 || currChargeCoreType2 == 131 ) { chargeparticle.setType( int(-13*charge) ); } - else if( currChargeCoreType2 == 110 || currChargeCoreType2 == 111 ) + else if( currChargeCoreType2 == 110 || currChargeCoreType2 == 111 ) { chargeparticle.setType( int(-11*charge) ); if(CHCluEnergy > CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy) + 1) { - flagEnergyFlow = 1; + flagEnergyFlow = 1; } } else @@ -1100,10 +1257,10 @@ void BushConnect::ParticleReco() } } - if( (currChargeCoreType2 != 130 && currChargeCoreType2 != 131 && CurrTrackEnergy > 15 && CHCluEnergy > 0.5 && CurrTrackEnergy > CHCluEnergy + 2.5*sqrt(CHCluEnergy) ) || (( currChargeCoreType2 == 130 || currChargeCoreType2 == 131 ) && CurrTrackEnergy > 100 && CHCluEnergy < 3 && chclustermerged.hits_size() < 20 ) ) + if( (currChargeCoreType2 != 130 && currChargeCoreType2 != 131 && CurrTrackEnergy > 15 && CHCluEnergy > 0.5 && CurrTrackEnergy > CHCluEnergy + 2.5*sqrt(CHCluEnergy) ) || (( currChargeCoreType2 == 130 || currChargeCoreType2 == 131 ) && CurrTrackEnergy > 100 && CHCluEnergy < 3 && chclustermerged.hits_size() < 20 ) ) { chargeparticle.setEnergy( CHCluEnergy ); - edm4hep::Vector3f CorrMom = edm4hep::Vector3f(CHCluEnergy/CurrTrackEnergy*currTrkP[0], CHCluEnergy/CurrTrackEnergy*currTrkP[1], CHCluEnergy/CurrTrackEnergy*currTrkP[2] ); + edm4hep::Vector3f CorrMom = edm4hep::Vector3f(CHCluEnergy/CurrTrackEnergy*currTrkP[0], CHCluEnergy/CurrTrackEnergy*currTrkP[1], CHCluEnergy/CurrTrackEnergy*currTrkP[2] ); chargeparticle.setMomentum( CorrMom ); } else @@ -1117,50 +1274,68 @@ void BushConnect::ParticleReco() auto a_Ef_Ne_particle = arborrecoparticleCol->create(); a_Ef_Ne_particle.setEnergy( CHCluEnergy - CurrTrackEnergy ); TVector3 corePos = TVector3(chclustermerged.getPosition().x,chclustermerged.getPosition().y,chclustermerged.getPosition().z); - float WFactor = (CHCluEnergy - CurrTrackEnergy)/corePos.Mag(); + float WFactor = (CHCluEnergy - CurrTrackEnergy)/corePos.Mag(); edm4hep::Vector3f PFNEMom = edm4hep::Vector3f(WFactor*float(corePos.X()), WFactor*float(corePos.Y()), WFactor*float(corePos.Z())); a_Ef_Ne_particle.setMomentum(PFNEMom); a_Ef_Ne_particle.setMass( 0.0 ); a_Ef_Ne_particle.setCharge( 0.0 ); a_Ef_Ne_particle.setType(501); - cout<<"Energy Flow Neutral Tagged "<<CHCluEnergy - CurrTrackEnergy<<endl; + cout<<"Energy Flow Neutral Tagged "<<CHCluEnergy - CurrTrackEnergy<<endl; + + cout<<"[YX debug - ParticleReco] EfPFO E = "<<a_Ef_Ne_particle.getEnergy()<<endl; + } - cout<<"a charged particle reconstructed with en:"<<chargeparticle.getEnergy()<<endl; + if(0) cout<<"a charged particle reconstructed with en:"<<chargeparticle.getEnergy()<<endl; + + if(0) cout<<"[YX debug - ParticleReco] ChPFO E = "<<chargeparticle.getEnergy()<<endl; } else // push non valid tracks, etc to second iteration, as those for PreInteracting ones { SecondIterTracks.push_back(a_chargedTrk); - cout<<"Second Iter Track Found"<<endl; - } + cout<<"Second Iter Track Found"<<endl; + } } - std::vector<edm4hep::MutableCluster> BBCore; + std::vector<edm4hep::MutableCluster> BBCore; BBCore.clear(); for(int p6 = 0; p6 < NNeutralCluster; p6 ++) { auto c_clu = non_chargedclustercore[p6]; + + if( NNCTouchFlag.find(c_clu) == NNCTouchFlag.end() ) { BBCore.push_back(c_clu); - } + + if(0) cout<<"[YX debug - BBCore] Clu "<<p6<<" push back"<<endl; + + } + + if(0) cout<<"[YX debug - BBCore] Clu "<<p6<<", NNCTouchFlag = "<<NNCTouchFlag[c_clu]<<", CluE = "<<c_clu.getEnergy()<<endl; + } float NAMom[3] = {0, 0, 0}; //Final Re-absorption - std::vector<edm4hep::Cluster> NBBNeutral; + std::vector<edm4hep::Cluster> NBBNeutral; NBBNeutral.clear(); for(int s = 0; s < int (BBCore.size()); s++) { auto a_clu = BBCore[s]; TVector3 PosClu = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); - float Depth = 0; + float Depth = 0; Depth = DisSeedSurface(PosClu); float CoreEnCorr = m_ArborToolLCIO->ClusterEE(a_clu); + int PhotonTag = m_ArborToolLCIO->newPhotonTag(a_clu); + double CluE = a_clu.getEnergy(); + if(0) cout<<"[YX debug - BBCore] Clu "<<s<<", newPhotonTag = "<<PhotonTag<<", CluE = "<<CluE<<", CoreEnCorr = "<<CoreEnCorr<<endl; + + if(m_ArborToolLCIO->newPhotonTag(a_clu)==1) { TVector3 BushSeedPos = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); @@ -1179,13 +1354,18 @@ void BushConnect::ParticleReco() a_neclu.setEnergy( CoreEnCorr ); //Reset... edm4hep::Cluster a_necluCon=a_neclu; neutralparticle.addToClusters(a_neclu); + + cout<<"[YX debug - BBCore] ---> Is photon "<<endl; + + cout<<"[YX debug - ParticleReco] Photon E = "<<neutralparticle.getEnergy()<<endl; + } else // Distance to Charged Core > sth; { float MinDisToChCore = 1.0E9; - float currDis = 0; + float currDis = 0; int NChCore = mergedclu_chCol->size(); - float closestChCluEn = 0; + float closestChCluEn = 0; for(int t = 0; t < NChCore; t++) { auto a_chclu = (*mergedclu_chCol)[t]; @@ -1196,30 +1376,50 @@ void BushConnect::ParticleReco() closestChCluEn = a_chclu.getEnergy(); // Or the Trk En?? } } + + if(0) cout<<"[YX debug - BBCore] ---> Is others: NChCore = "<<NChCore<<", MinDisToChCore = "<<MinDisToChCore<<", closestChCluEn = "<<closestChCluEn<<endl; + if( MinDisToChCore > 0.4*(15 + closestChCluEn + Depth*0.01) || a_clu.getEnergy() > 2.0 ) //Joint Depth?? { NBBNeutral.push_back(a_clu); + + if(0) cout<<"[YX debug - BBCore] ---|---> NBBNeutral push back "<<endl; + } } + + } + + cout<<"[YX debug - ParticleReco] BBCore.size() = "<<BBCore.size()<<", NBBNeutral.size() = "<<NBBNeutral.size()<<", IsoHits.size() = "<<IsoHits.size()<<endl; + + // Add: Neural Core Remerge & Energy Scale Recalculate, IsoHit Abso std::vector<edm4hep::MutableCluster> NBBAbs = m_ArborToolLCIO->ClusterHitAbsorbtion(NBBNeutral, IsoHits, 100); //_HitAbsCut); // Huge?? - std::vector<float> BBAbsEn; + std::vector<float> BBAbsEn; BBAbsEn.clear(); + // double ESum_NBBAbs = 0; + for(unsigned s1 = 0; s1 < NBBAbs.size(); s1++) { BBAbsEn.push_back(NBBAbs[s1].getEnergy()); + + // ESum_NBBAbs+=NBBAbs[s1].getEnergy(); + } + // cout<<"[YX debug - ParticleReco] NBBAbs.size() = "<<NBBAbs.size()<<", ESum_NBBAbs = "<<ESum_NBBAbs<<endl; + + std::vector<int> BBAbsIndex = SortMeasure(BBAbsEn, 1); std::vector<edm4hep::Cluster > NeutronCore; std::vector<edm4hep::MutableCluster > NeutronFlag; NeutronCore.clear(); - NeutronFlag.clear(); + NeutronFlag.clear(); for(unsigned int s2 = 0; s2 < NBBAbs.size(); s2++) //Sort it; the first one must be a neutral core? { @@ -1238,6 +1438,10 @@ void BushConnect::ParticleReco() std::vector<edm4hep::MutableCluster > Neutrons = m_ArborToolLCIO->ClusterAbsorbtion(NeutronCore, NeutronFlag, 200, 0.01); + + cout<<"[YX debug - ParticleReco] NeutronCore.size() = "<<NeutronCore.size()<<", NeutronFlag.size() = "<<NeutronFlag.size()<<endl; + + for(unsigned int s3 = 0; s3 < Neutrons.size(); s3++) { auto a_clu = Neutrons[s3]; @@ -1245,7 +1449,7 @@ void BushConnect::ParticleReco() TVector3 PosClu = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); float MinDisToChCore = 1.0E9; float RecoT0 = m_ArborToolLCIO->ClusterT0(a_clu); - float currDis = 0; + float currDis = 0; int NChCore = mergedclu_chCol->size(); for(int t = 0; t < NChCore; t++) { @@ -1257,10 +1461,14 @@ void BushConnect::ParticleReco() } } + + if(0) cout<<"[YX debug - Neutrons] Clu "<<s3<<", MinDisToChCore = "<<MinDisToChCore<<", CluE = "<<a_clu.getEnergy()<<", CoreEnCorr = "<<CoreEnCorr<<endl; + + if( !(RecoT0>0.1 && RecoT0<1E8 && MinDisToChCore <12) ) { if(m_ArborToolLCIO->newPhotonTag(a_clu)==1) - cout<<"WARNING... Photons after neutron merge merged"<<endl; + cout<<"WARNING... Photons after neutron merge merged"<<endl; auto neutralparticle = arborrecoparticleCol->create(); neutralparticle.setType(21120); TVector3 PP = m_ArborToolLCIO->ClusterCoG(a_clu); @@ -1301,27 +1509,34 @@ void BushConnect::ParticleReco() a_neclu.setEnergy( CoreEnCorr ); //Reset... edm4hep::Cluster a_necluCon=a_neclu; neutralparticle.addToClusters(a_necluCon); + + if(0) cout<<"[YX debug - ParticleReco] NePFO E = "<<neutralparticle.getEnergy()<<endl; + } } - + +// cout<<"[YX debug] ParticleReco End."<<endl; } StatusCode BushConnect::execute() { try{ - BushConnect::Clean(); + cout<<"[YX debug - BushConnect] Begin:"<<endl; BushConnect::TrackSort( ); - BushConnect::BushSelfMerge( ); - BushConnect::TagCore( ); + BushConnect::BushSelfMerge( ); + BushConnect::TagCore( ); BushConnect::ParticleReco( ); + + BushConnect::Clean(); + }catch(GaudiException &e){} return StatusCode::SUCCESS; } StatusCode BushConnect::finalize() { - std::cout<<"Bush Connection Finished, ArborObject Formed"<<std::endl; + std::cout<<"Bush Connection Finished, ArborObject Formed"<<std::endl; return GaudiAlgorithm::finalize(); } diff --git a/Reconstruction/PFA/Arbor/src/BushConnect.hh b/Reconstruction/PFA/Arbor/src/BushConnect.hh index 7cf32b9a248d449e1dbf9277827530d237f05b6c..075773b171b2bc196b0f6cafd5b633e511a1d9e3 100644 --- a/Reconstruction/PFA/Arbor/src/BushConnect.hh +++ b/Reconstruction/PFA/Arbor/src/BushConnect.hh @@ -41,11 +41,11 @@ class BushConnect : public GaudiAlgorithm virtual StatusCode execute() ; virtual StatusCode finalize() ; - void Clean(); - void TrackSort(); - void BushSelfMerge(); - void TagCore(); - void ParticleReco(); + void Clean(); + void TrackSort(); + void BushSelfMerge(); + void TagCore(); + void ParticleReco(); protected: @@ -56,10 +56,10 @@ class BushConnect : public GaudiAlgorithm std::map<edm4hep::Track, TVector3> Track_P3; std::map<edm4hep::Track, int> Track_Type; std::map<edm4hep::Track, float> Track_Theta; - std::map<edm4hep::Track, float> Track_Phi; + std::map<edm4hep::Track, float> Track_Phi; std::map<edm4hep::Cluster, int> ClusterType_1stID; - std::map<edm4hep::ReconstructedParticle, int> ChCoreID; + std::map<edm4hep::ReconstructedParticle, int> ChCoreID; std::vector<edm4hep::Cluster> ecalchcore_tight; //TightCores std::vector<edm4hep::Cluster> ecalchcore_medium; @@ -79,7 +79,7 @@ class BushConnect : public GaudiAlgorithm std::vector<edm4hep::Cluster> chargedclustercore; std::vector<edm4hep::Cluster> chargedclustercore_abs; - std::vector<edm4hep::MutableCluster> selfmergedcluster; + std::vector<edm4hep::MutableCluster> selfmergedcluster; std::vector<edm4hep::MutableCluster> non_chargedclustercore; std::vector<edm4hep::Cluster> onlyNeutralCore; @@ -89,9 +89,9 @@ class BushConnect : public GaudiAlgorithm std::map<edm4hep::Track, int>MCPTrack_Type; std::map<edm4hep::Track, TVector3> Track_EndPoint; //Last hit std::map<edm4hep::Track, TVector3> TrackStartPoint; - std::map<edm4hep::Cluster, float> CluFD; + std::map<edm4hep::Cluster, float> CluFD; std::map<edm4hep::Cluster, float> CluT0; - std::map<edm4hep::Cluster, float> Clu_Depth; + std::map<edm4hep::Cluster, float> Clu_Depth; std::map<edm4hep::Cluster, TVector3> CluCoG; typedef DataHandle<edm4hep::MCParticleCollection> MCParticleColHandler; MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this}; @@ -113,8 +113,10 @@ class BushConnect : public GaudiAlgorithm DataHandle<edm4hep::ReconstructedParticleCollection> m_chargeparticleCol{"ArborCharged",Gaudi::DataHandle::Writer, this}; DataHandle<edm4hep::ReconstructedParticleCollection> nerecoparticleCol{"ArborNeutral",Gaudi::DataHandle::Writer, this}; DataHandle<edm4hep::ReconstructedParticleCollection> m_arborrecoparticleCol{"ArborPFO",Gaudi::DataHandle::Writer, this}; - ArborToolLCIO * m_ArborToolLCIO; - + ArborToolLCIO * m_ArborToolLCIO; + + Gaudi::Property<bool> m_readLCIO{this, "ReadLCIO", true, "Read sim file with LCIO"}; + }; diff --git a/Reconstruction/PFA/Arbor/src/DetectorPos.cc b/Reconstruction/PFA/Arbor/src/DetectorPos.cc index cc70613b44bc445fa51e4e885d57fad0224a3106..ed32f0bc02161876f7377bbfacefa510ecd82249 100644 --- a/Reconstruction/PFA/Arbor/src/DetectorPos.cc +++ b/Reconstruction/PFA/Arbor/src/DetectorPos.cc @@ -1,22 +1,23 @@ #include <TMath.h> #include "DetectorPos.hh" +// #include "Arbor/DetectorPos.hh" using namespace std; -const double pi = acos(-1.0); +const double pi = acos(-1.0); //Geometric Parameter - ... need to be changed for different detector models int BarrelFlag( TVector3 inputPos) { - int isBarrel = 0; + int isBarrel = 0; if(fabs(inputPos[2]) < ECALHalfZ) { - isBarrel = 1; + isBarrel = 1; } - return isBarrel; + return isBarrel; } int DepthFlag( TVector3 inputPos ) //Used to calculate depth of given position... @@ -79,7 +80,7 @@ TVector3 CalVertex( TVector3 Pos1, TVector3 Dir1, TVector3 Pos2, TVector3 Dir2 ) int TPCPosition( TVector3 inputPos ) { - int flagPos(-1); // == 0 means inside TPC, == 1 means outside; + int flagPos(-1); // == 0 means inside TPC, == 1 means outside; if( fabs(inputPos[2]) > ECALHalfZ || sqrt( inputPos[0]*inputPos[0] + inputPos[1]*inputPos[1] ) > TPCRadius ) flagPos = 1; else flagPos = 0; @@ -135,7 +136,7 @@ float DisSeedSurfaceSimple( TVector3 SeedPos ) //ECAL, HCAL, EndCapRing.. float DisR = SeedPos.Perp() - ECALRadius; if(DisR < 0 && DisZ > 0) { - DisSS = DisZ; + DisSS = DisZ; } else if(DisZ < 0 && DisR > 0) { @@ -167,25 +168,25 @@ float DisSeedSurfaceClu( Cluster * a_clu ) //ECAL, HCAL, EndCapRing... float DisTPCBoundary( TVector3 Pos ) { float DisZ = TMath::Min( fabs(ECALHalfZ-Pos.Z()),fabs(ECALHalfZ+Pos.Z()) ); - float DisR = ECALRadius - Pos.Perp(); + float DisR = ECALRadius - Pos.Perp(); float Dis = TMath::Min(DisZ, DisR); - return Dis; + return Dis; } float DistanceChargedParticleToCluster(TVector3 CPRefPos, TVector3 CPRefMom, TVector3 CluPosition) //Extend to Track/MCP { // Line extrapolation from RefPos with RefMom, calculate the minimal distance to Cluster - float DisCPClu = 0; - TVector3 Diff_Clu_CPRef, NormCPRefMom; + float DisCPClu = 0; + TVector3 Diff_Clu_CPRef, NormCPRefMom; - Diff_Clu_CPRef = CluPosition - CPRefPos; + Diff_Clu_CPRef = CluPosition - CPRefPos; NormCPRefMom = 1.0/CPRefMom.Mag()*CPRefMom; - float ProDis = Diff_Clu_CPRef.Dot(NormCPRefMom); + float ProDis = Diff_Clu_CPRef.Dot(NormCPRefMom); DisCPClu = sqrt(Diff_Clu_CPRef.Mag()*Diff_Clu_CPRef.Mag() - ProDis*ProDis); - return DisCPClu; + return DisCPClu; } diff --git a/Reconstruction/PFA/Arbor/src/HelixClassD.cc b/Reconstruction/PFA/Arbor/src/HelixClassD.cc index 20b79466237d97a8d7bf5314ec00d3867cb5e453..023a336eebcbd70ec18b73be056726d4aaead65b 100644 --- a/Reconstruction/PFA/Arbor/src/HelixClassD.cc +++ b/Reconstruction/PFA/Arbor/src/HelixClassD.cc @@ -40,7 +40,7 @@ void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) { double radius = pxy/double(_FCT*B); double xCentre = double(pos[0]) + radius*double(cos(_phiMomRefPoint-_const_pi2*q)); double yCentre = double(pos[1]) + radius*double(sin(_phiMomRefPoint-_const_pi2*q)); - + double d0; if (q>0) { @@ -54,17 +54,17 @@ void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) { // if (fabs(_d0)>0.001 ) { // std::cout << "New helix : " << std::endl; -// std::cout << " Position : " << pos[0] +// std::cout << " Position : " << pos[0] // << " " << pos[1] // << " " << pos[2] << std::endl; // std::cout << " Radius = " << _radius << std::endl; -// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl; +// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl; // std::cout << " D0 = " << _d0 << std::endl; // } _pxAtPCA = _pxy*cos(_phi0); _pyAtPCA = _pxy*sin(_phi0); - float deltaPhi = _phiRefPoint - _phiAtPCA; + float deltaPhi = _phiRefPoint - _phiAtPCA; float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi; xCircles = xCircles/_const_2pi; int nCircles; @@ -78,7 +78,7 @@ void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) { n1 = int(xCircles) - 1; n2 = n1 + 1; } - + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { nCircles = n1; } @@ -89,7 +89,7 @@ void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) { } -void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, +void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, float omega, float tanLambda, float B) { _omega = omega; _d0 = d0; @@ -99,20 +99,20 @@ void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, _charge = omega/fabs(omega); _radius = 1./fabs(omega); _xAtPCA = -_d0*sin(_phi0); - _yAtPCA = _d0*cos(_phi0); + _yAtPCA = _d0*cos(_phi0); _referencePoint[0] = _xAtPCA; _referencePoint[1] = _yAtPCA; _referencePoint[2] = _z0; _pxy = _FCT*B*_radius; _momentum[0] = _pxy*cos(_phi0); _momentum[1] = _pxy*sin(_phi0); - _momentum[2] = _tanLambda * _pxy; + _momentum[2] = _tanLambda * _pxy; _pxAtPCA = _momentum[0]; _pyAtPCA = _momentum[1]; _phiMomRefPoint = atan2(_momentum[1],_momentum[0]); - _xCentre = _referencePoint[0] + + _xCentre = _referencePoint[0] + _radius*cos(_phi0-_const_pi2*_charge); - _yCentre = _referencePoint[1] + + _yCentre = _referencePoint[1] + _radius*sin(_phi0-_const_pi2*_charge); _phiAtPCA = atan2(-_yCentre,-_xCentre); _phiRefPoint = _phiAtPCA ; @@ -120,7 +120,7 @@ void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, } -void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, +void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, float bZ, float phi0, float B, float signPz, float zBegin) { @@ -148,8 +148,8 @@ void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, _tanLambda = _momentum[2]/_pxy; _momentum[0] = _pxy*cos(_phiMomRefPoint); _momentum[1] = _pxy*sin(_phiMomRefPoint); - - float deltaPhi = _phiRefPoint - _phiAtPCA; + + float deltaPhi = _phiRefPoint - _phiAtPCA; float xCircles = bZ*_referencePoint[2] - deltaPhi; xCircles = xCircles/_const_2pi; int nCircles; @@ -163,14 +163,14 @@ void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, n1 = int(xCircles) - 1; n2 = n1 + 1; } - + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { nCircles = n1; } else { nCircles = n2; - } - _z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ; + } + _z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ; _bField = B; } @@ -225,7 +225,7 @@ float HelixClassD::getCharge() { return _charge; } -float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, +float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, float * ref , float * point) { float time; @@ -234,7 +234,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, if (AA <= 0) { - time = -1.0e+20; + time = -1.0e+20; return time; } @@ -242,7 +242,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, float BB = ax*(x0-_xCentre) + ay*(y0-_yCentre); BB = BB / AA; - float CC = (x0-_xCentre)*(x0-_xCentre) + float CC = (x0-_xCentre)*(x0-_xCentre) + (y0-_yCentre)*(y0-_yCentre) - _radius*_radius; CC = CC / AA; @@ -250,7 +250,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, float DET = BB*BB - CC; float tt1 = 0.; float tt2 = 0.; - float xx1,xx2,yy1,yy2; + float xx1,xx2,yy1,yy2; if (DET < 0 ) { @@ -260,8 +260,8 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, point[2]=0.0; return time; } - - + + tt1 = - BB + sqrt(DET); tt2 = - BB - sqrt(DET); @@ -269,7 +269,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, yy1 = y0 + tt1*ay; xx2 = x0 + tt2*ax; yy2 = y0 + tt2*ay; - + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); @@ -280,14 +280,14 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, if (dphi1 < 0 && _charge < 0) { dphi1 = dphi1 + _const_2pi; } - else if (dphi1 > 0 && _charge > 0) { + else if (dphi1 > 0 && _charge > 0) { dphi1 = dphi1 - _const_2pi; } if (dphi2 < 0 && _charge < 0) { dphi2 = dphi2 + _const_2pi; } - else if (dphi2 > 0 && _charge > 0) { + else if (dphi2 > 0 && _charge > 0) { dphi2 = dphi2 - _const_2pi; } @@ -299,7 +299,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, std::cout << "WARNING " << tt1 << std::endl; if (tt2 < 0. ) std::cout << "WARNING " << tt2 << std::endl; - + if (tt1 < tt2) { point[0] = xx1; @@ -314,7 +314,7 @@ float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, point[2] = ref[2]+time*_momentum[2]; - + return time; @@ -340,17 +340,17 @@ float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) { } float phiCentre = atan2(_yCentre,_xCentre); - float phiStar = Radius*Radius + distCenterToIP*distCenterToIP + float phiStar = Radius*Radius + distCenterToIP*distCenterToIP - _radius*_radius; phiStar = 0.5*phiStar/fmax(1.0e-20,Radius*distCenterToIP); - - if (phiStar > 1.0) + + if (phiStar > 1.0) phiStar = 0.9999999; - + if (phiStar < -1.0) phiStar = -0.9999999; - + phiStar = acos(phiStar); float tt1,tt2,time; @@ -372,14 +372,14 @@ float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) { if (dphi1 < 0 && _charge < 0) { dphi1 = dphi1 + _const_2pi; } - else if (dphi1 > 0 && _charge > 0) { + else if (dphi1 > 0 && _charge > 0) { dphi1 = dphi1 - _const_2pi; } if (dphi2 < 0 && _charge < 0) { dphi2 = dphi2 + _const_2pi; } - else if (dphi2 > 0 && _charge > 0) { + else if (dphi2 > 0 && _charge > 0) { dphi2 = dphi2 - _const_2pi; } @@ -391,7 +391,7 @@ float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) { std::cout << "WARNING " << tt1 << std::endl; if (tt2 < 0. ) std::cout << "WARNING " << tt2 << std::endl; - + float time2; if (tt1 < tt2) { @@ -413,7 +413,7 @@ float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) { point[2] = ref[2]+time*_momentum[2]; point[5] = ref[2]+time2*_momentum[2]; - + return time; @@ -490,7 +490,7 @@ float HelixClassD::getDistanceToPoint(float * xPoint, float * Distance) { float DistXY = (_xCentre-xPoint[0])*(_xCentre-xPoint[0]) + (_yCentre-xPoint[1])*(_yCentre-xPoint[1]); DistXY = sqrt(DistXY); DistXY = fabs(DistXY - _radius); - + int nCircles = 0; if (fabs(_tanLambda*_radius)>1.0e-20) { @@ -506,7 +506,7 @@ float HelixClassD::getDistanceToPoint(float * xPoint, float * Distance) { n1 = int(xCircles) - 1; n2 = n1 + 1; } - + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { nCircles = n1; } @@ -515,7 +515,7 @@ float HelixClassD::getDistanceToPoint(float * xPoint, float * Distance) { } } - + float DPhi = _const_2pi*((float)nCircles) + phi - phi0; zOnHelix = _referencePoint[2] - _charge*_radius*_tanLambda*DPhi; @@ -601,10 +601,10 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * float x01 = getXC(); float y01 = getYC(); - + float x02 = helix->getXC(); float y02 = helix->getYC(); - + float rad1 = getRadius(); float rad2 = helix->getRadius(); @@ -635,7 +635,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * phi1 = phi0 + alpha; phi2 = phi0 - alpha; } - + float ref1[3]; float ref2[3]; @@ -643,7 +643,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * ref1[i]=_referencePoint[i]; ref2[i]=helix->getReferencePoint()[i]; } - + float pos1[3]; float pos2[3]; float mom1[3]; @@ -661,9 +661,9 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * getPointOnCircle(R1,ref1,pos1); helix->getPointOnCircle(R2,ref2,pos2); - + } - else { + else { float xSect1 = x02 + rad2*cos(phi1); float ySect1 = y02 + rad2*sin(phi1); @@ -676,7 +676,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * float temp21[3]; float temp22[3]; - float phiI2 = atan2(ref2[1]-y02,ref2[0]-x02); + float phiI2 = atan2(ref2[1]-y02,ref2[0]-x02); float phiF21 = atan2(ySect1-y02,xSect1-x02); float phiF22 = atan2(ySect2-y02,xSect2-x02); float deltaPhi21 = phiF21 - phiI2; @@ -687,14 +687,14 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * if (deltaPhi21 < 0 && charge2 < 0) { deltaPhi21 += _const_2pi; } - else if (deltaPhi21 > 0 && charge2 > 0) { + else if (deltaPhi21 > 0 && charge2 > 0) { deltaPhi21 -= _const_2pi; } if (deltaPhi22 < 0 && charge2 < 0) { deltaPhi22 += _const_2pi; } - else if (deltaPhi22 > 0 && charge2 > 0) { + else if (deltaPhi22 > 0 && charge2 > 0) { deltaPhi22 -= _const_2pi; } @@ -706,7 +706,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * temp21[0] = xSect1; temp21[1] = ySect1; temp21[2] = Z21; temp22[0] = xSect2; temp22[1] = ySect2; temp22[2] = Z22; - + // std::cout << "temp21 = " << temp21[0] << " " << temp21[1] << " " << temp21[2] << std::endl; // std::cout << "temp22 = " << temp22[0] << " " << temp22[1] << " " << temp22[2] << std::endl; @@ -715,7 +715,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * float temp11[3]; float temp12[3]; - float phiI1 = atan2(ref1[1]-y01,ref1[0]-x01); + float phiI1 = atan2(ref1[1]-y01,ref1[0]-x01); float phiF11 = atan2(ySect1-y01,xSect1-x01); float phiF12 = atan2(ySect2-y01,xSect2-x01); float deltaPhi11 = phiF11 - phiI1; @@ -726,14 +726,14 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * if (deltaPhi11 < 0 && charge1 < 0) { deltaPhi11 += _const_2pi; } - else if (deltaPhi11 > 0 && charge1 > 0) { + else if (deltaPhi11 > 0 && charge1 > 0) { deltaPhi11 -= _const_2pi; } if (deltaPhi12 < 0 && charge1 < 0) { deltaPhi12 += _const_2pi; } - else if (deltaPhi12 > 0 && charge1 > 0) { + else if (deltaPhi12 > 0 && charge1 > 0) { deltaPhi12 -= _const_2pi; } @@ -745,7 +745,7 @@ float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * temp11[0] = xSect1; temp11[1] = ySect1; temp11[2] = Z11; temp12[0] = xSect2; temp12[1] = ySect2; temp12[2] = Z12; - + // std::cout << "temp11 = " << temp11[0] << " " << temp11[1] << " " << temp11[2] << std::endl; // std::cout << "temp12 = " << temp12[0] << " " << temp12[1] << " " << temp12[2] << std::endl; diff --git a/Reconstruction/PFA/Arbor/src/MarlinArbor.cc b/Reconstruction/PFA/Arbor/src/MarlinArbor.cc index fb930a683b1f6fb98db07028e453d2f8fc08a5dc..3ebba5dcbe7f8ea7c97f15f0b3e14c99d0073027 100644 --- a/Reconstruction/PFA/Arbor/src/MarlinArbor.cc +++ b/Reconstruction/PFA/Arbor/src/MarlinArbor.cc @@ -18,6 +18,11 @@ #include "edm4hep/MCParticleCollection.h" #include "cellIDDecoder.h" +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +#include "DecoderHelper/DD4hep2Lcio.h" #include "DD4hep/Detector.h" #include "DD4hep/IDDescriptor.h" @@ -42,11 +47,11 @@ using namespace std; -extern linkcoll InitLinks; -extern linkcoll IterLinks_1; -extern linkcoll IterLinks; -extern linkcoll links_debug; -extern branchcoll Trees; +extern linkcoll InitLinks; +extern linkcoll IterLinks_1; +extern linkcoll IterLinks; +extern linkcoll links_debug; +extern branchcoll Trees; extern std::vector<int> IsoHitsIndex; //std::vector<std::string> CaloHitCollections; @@ -54,23 +59,30 @@ extern std::vector<int> IsoHitsIndex; DECLARE_COMPONENT(MarlinArbor) MarlinArbor::MarlinArbor(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc), + : GaudiAlgorithm(name, svcLoc), _eventNr(0),_output(0) { } StatusCode MarlinArbor::initialize() { - - _cepc_thresholds.push_back(10); +// ??? + // _cepc_thresholds.push_back(10); + // _cepc_thresholds.push_back(90); + // _cepc_thresholds.push_back(50); + // _cepc_thresholds.push_back(7.5); + + m_encoder_str = "M:3,S-1:3,I:9,J:9,K-1:6"; + _cepc_thresholds.push_back(20); _cepc_thresholds.push_back(90); _cepc_thresholds.push_back(50); - _cepc_thresholds.push_back(7.5); + _cepc_thresholds.push_back(11); m_geosvc = service<IGeomSvc>("GeomSvc"); ISvcLocator* svcloc = serviceLocator(); - m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc); + // m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc); + m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc,m_readLCIO); for(unsigned int i = 0; i < m_ecalReadoutNames.value().size(); i++){ m_col_readout_map[m_ecalColNames.value().at(i)] = m_ecalReadoutNames.value().at(i); } @@ -117,66 +129,181 @@ StatusCode MarlinArbor::execute() //if(_eventNr % m_reportEvery == 0) cout<<"eventNr: "<<_eventNr<<endl; _eventNr++; - MarlinArbor::HitsPreparation(); //Absorb isolated hits; + MarlinArbor::HitsPreparation(); //Absorb isolated hits; TVector3 currHitPos; std::vector< TVector3 > inputHitsPos; - std::vector< ArborHit > inputABHit; - std::vector< edm4hep::CalorimeterHit > inputHits; - std::vector< edm4hep::CalorimeterHit > inputECALHits; - std::vector< edm4hep::CalorimeterHit > inputHCALHits; - std::vector< std::vector<int> > Sequence; - int LayerNum = 0; - int StaveNum = 0; - int SubDId = -10; - float Depth = 0; - int KShift = 0; - TVector3 TrkEndPointPos; + std::vector< ArborHit > inputABHit; + std::vector< edm4hep::CalorimeterHit > inputHits; + std::vector< edm4hep::CalorimeterHit > inputECALHits; + std::vector< edm4hep::CalorimeterHit > inputHCALHits; + std::vector< std::vector<int> > Sequence; + int LayerNum = 0; + int StaveNum = 0; + int SubDId = -10; + float Depth = 0; + int KShift = 0; + TVector3 TrkEndPointPos; std::vector<edm4hep::CalorimeterHit> IsoHits; + + cout<<"[YX debug - MarlinArbor] InputCollections.size() = "<<_calCollections.size()<<endl; + + unsigned int nECALCol = m_ecalColNames.size(); + unsigned int nHCALCol = m_hcalColNames.size(); + for(unsigned int i1 = 0; i1 < _calCollections.size(); i1++) { - std::cout<<i1<<"th collection"<<m_col_readout_map[m_ecalColNames.value().at(i1)]<<std::endl; + // cout<<"[YX debug - MarlinArbor] CaloHitCollections "<<i1<<endl; + + if(i1<nECALCol) + std::cout<<"[YX debug - MarlinArbor] CaloHitCollections "<<i1<<": "<<m_col_readout_map[m_ecalColNames.value().at(i1)]<<std::endl; + else + std::cout<<"[YX debug - MarlinArbor] CaloHitCollections "<<i1<<": "<<m_col_readout_map[m_hcalColNames.value().at(i1-nECALCol)]<<std::endl; + + std::string tmp_readout; - - if(i1<2)tmp_readout = m_col_readout_map[m_ecalColNames.value().at(i1)]; - else - tmp_readout = m_col_readout_map[m_hcalColNames.value().at(i1-2)]; - - std::cout<<tmp_readout<<std::endl; - // get the DD4hep readout - m_decoder = m_geosvc->getDecoder(tmp_readout); - KShift = 0; - SubDId = -1; - if( i1 < _EcalCalCollections.size() ) - SubDId = 1; - else if( i1 < _EcalCalCollections.size() + _HcalCalCollections.size() ) + // if(i1<2)tmp_readout = m_col_readout_map[m_ecalColNames.value().at(i1)]; + if(i1<nECALCol) + tmp_readout = m_col_readout_map[m_ecalColNames.value().at(i1)]; + else + tmp_readout = m_col_readout_map[m_hcalColNames.value().at(i1-nECALCol)]; + // tmp_readout = m_col_readout_map[m_hcalColNames.value().at(i1-2)]; + + std::cout<<"tmp_readout: "<<tmp_readout<<std::endl; + + // // get the DD4hep readout + // m_decoder = m_geosvc->getDecoder(tmp_readout); + + KShift = 0; + SubDId = -1; + + // ??? + // if( i1 < _EcalCalCollections.size() ) + // SubDId = 1; + // else if( i1 < _EcalCalCollections.size() + _HcalCalCollections.size() ) + // SubDId = 2; + // else + // SubDId = 3; + + // if(i1 > _EcalCalCollections.size() - 1) + // KShift = 100; + // else if( i1 == _calCollections.size() - 2) //HCAL Ring + // KShift = 50; + + if( i1 < _ecalCollections.size() ) + SubDId = 1; + else if( i1 < _ecalCollections.size() + _hcalCollections.size() ) SubDId = 2; else - SubDId = 3; + SubDId = 3; - if(i1 > _EcalCalCollections.size() - 1) - KShift = 100; + if(i1 > _ecalCollections.size() - 1) + KShift = 100; else if( i1 == _calCollections.size() - 2) //HCAL Ring KShift = 50; + // !!! auto CaloHitColl = _calCollections[i1]->get(); - + + int i2 = 0; + //int NHitsCurrCol = CaloHitColl->getNumberOfElements(); //CellIDDecoder<CalorimeterHit> idDecoder(CaloHitColl); - for (auto a_hit: *CaloHitColl){ - currHitPos = TVector3(a_hit.getPosition().x, a_hit.getPosition().y, a_hit.getPosition().z); - Depth = DisSeedSurface(currHitPos); - - auto cellid = a_hit.getCellID(); - LayerNum = m_decoder->get(cellid, "layer")+ KShift; - StaveNum=m_decoder->get(cellid, "stave"); - - if(SubDId!=2 ){ + for(auto a_hit: *CaloHitColl){ + ID_UTIL::CellIDDecoder<edm4hep::CalorimeterHit> cellIdDecoder(m_encoder_str); + const std::string layerCodingString(m_encoder_str); + const std::string staveCodingString(m_encoder_str); + const std::string idCodingString(m_encoder_str); + + const std::string staveCoding(m_ArborToolLCIO->GetStaveCoding(staveCodingString)); + const std::string layerCoding(m_ArborToolLCIO->GetLayerCoding(layerCodingString)); + // const std::string layerCoding(m_ArborToolLCIO->GetLayerCoding(idCodingString)); + const std::string cellICoding(m_ArborToolLCIO->GetCellICoding(idCodingString)); + const std::string cellJCoding(m_ArborToolLCIO->GetCellJCoding(idCodingString)); + + if(!m_readLCIO) + m_decoder = m_geosvc->getDecoder(tmp_readout); + + // m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); + // if(!m_decoder) m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + + currHitPos = TVector3(a_hit.getPosition().x, a_hit.getPosition().y, a_hit.getPosition().z); + Depth = DisSeedSurface(currHitPos); + + if(m_readLCIO){ + LayerNum=cellIdDecoder(&a_hit)[layerCoding.c_str()] + 1 + KShift; + StaveNum=cellIdDecoder(&a_hit)[staveCoding.c_str()] + 1 ; + + if(1){ + double hitposx = currHitPos.X(); + double hitposy = currHitPos.Y(); + double hitposz = currHitPos.Z(); + double hitposp = currHitPos.Perp(); + + // int tmp_M = idDecoder(a_hit)["M"]; + int tmp_S = cellIdDecoder(&a_hit)[staveCoding.c_str()]; + int tmp_K = cellIdDecoder(&a_hit)[layerCoding.c_str()]; + int tmp_I = cellIdDecoder(&a_hit)[cellICoding.c_str()]; + int tmp_J = cellIdDecoder(&a_hit)[cellJCoding.c_str()]; + + if(0) cout<<"[MarlinArbor] Hit Pos ("<<hitposx<<", "<<hitposy<<", "<<hitposz<<", "<<hitposp<<")mm:"<<endl; + if(0) cout<<"[MarlinArbor] ---> M = xx, stave = "<<tmp_S<<", layer = "<<tmp_K<<", I(x) = "<<tmp_I<<", J(y) = "<<tmp_J<<endl; + if(0) cout<<"[MarlinArbor] ---> StaveNum = "<<StaveNum<<", LayerNum = "<<LayerNum<<endl; + } + + }else{ + auto cellid = a_hit.getCellID(); + + + // SEcal05_siw_Barrel <id>system:5,module:3,stave:4,tower:5,layer:6,wafer:6,cellX:32:-16,cellY:-16</id> + // SEcal05_siw_Endcaps <id>system:5,module:3,stave:4,tower:5,layer:6,wafer:6,x:32:-16,y:-16</id> + // SEcal05_siw_ECRing_01 <id>system:5,module:3,stave:4,tower:3,layer:6,x:32:-16,y:-16</id> + // SHcalRpc01_Barrel_01 <id>system:5,module:3,stave:3,tower:5,layer:6,slice:4,x:32:-16,y:-16</id> + // SHcalRpc01_Endcaps_01 <id>system:5,module:3,stave:3,tower:5,layer:6,y:32:-16,x:-16</id> + // SHcalRpc01_EndcapRing_01 <id>system:5,module:3,stave:4,tower:3,layer:6,y:32:-16,x:-16</id> + + int Raw_system = m_decoder->get(cellid, "system"); + int Raw_module = m_decoder->get(cellid, "module"); + int Raw_Stave = m_decoder->get(cellid, "stave"); + int Raw_Layer = m_decoder->get(cellid, "layer"); + + int New_Layer = DD4hep2Lcio::CEPCv4::getEcalLayer(Raw_Layer); + int New_Stave = DD4hep2Lcio::CEPCv4::getEcalBarrelStave(Raw_Stave); + if(Raw_system==29){// ECAL endcap + New_Stave = DD4hep2Lcio::CEPCv4::getEcalEndcapStave(Raw_Stave); + } + if(Raw_system==22 || Raw_system==30){ + New_Layer = DD4hep2Lcio::CEPCv4::getHcalLayer(Raw_Layer); + New_Stave = DD4hep2Lcio::CEPCv4::getHcalStave(Raw_Stave); + } + LayerNum = New_Layer + KShift + 1; + StaveNum = New_Stave + 1; + + // LayerNum = m_decoder->get(cellid, "layer") + KShift; + // StaveNum = m_decoder->get(cellid, "stave"); + + + if(0){ + double hitposx = currHitPos.X(); + double hitposy = currHitPos.Y(); + double hitposz = currHitPos.Z(); + double hitposp = currHitPos.Perp(); + + cout<<"[MarlinArbor] Hit Pos ("<<hitposx<<", "<<hitposy<<", "<<hitposz<<", "<<hitposp<<")mm:"<<endl; + // cout<<"---> stave = "<<Raw_Stave<<", layer = "<<Raw_Layer<<", X = "<<Raw_cellX<<", Y = "<<Raw_cellY<<", wafer = "<<Raw_Wafer<<", tower = "<<Raw_Tower<<endl; + + cout<<"[MarlinArbor] ---> Raw M = "<<Raw_module<<", stave = "<<Raw_Stave<<", layer = "<<Raw_Layer<<", system = "<<Raw_system<<endl; + cout<<"[MarlinArbor] ---> New M = "<<Raw_module<<", stave = "<<New_Stave<<", layer = "<<New_Layer<<", system = "<<Raw_system<<endl; + cout<<"[MarlinArbor] ---> StaveNum = "<<StaveNum<<", LayerNum = "<<LayerNum<<endl; + } + + } + if(SubDId!=2 ){ inputECALHits.push_back(a_hit); } else{ @@ -185,15 +312,29 @@ StatusCode MarlinArbor::execute() ArborHit a_abhit(currHitPos, LayerNum, 0, Depth, StaveNum, SubDId); inputABHit.push_back(a_abhit); inputHits.push_back(a_hit); + + + // cout<<"[YX debug - MarlinArbor] Hit "<<i2<<", KShift = "<<KShift<<", LayerNum = "<<LayerNum<<", StaveNum = "<<StaveNum<<", SubDId = "<<SubDId<<", Depth = "<<Depth<<endl; + // // auto cellID = a_hit.getCellID(); + // // cout<<"[YX debug - MarlinArbor] ---> raw layerNum = "<<m_decoder->get(cellID, "layer")<<", staveNum = "<<m_decoder->get(cellID, "stave")<<endl; + // // cout<<"[YX debug - MarlinArbor] ---> raw cellID (DigiHit) = "<<cellID<<endl; + + + + + i2+=1; } - - // cout<<i1<<" Stat "<<SubDId<<" ~~~ "<<inputABHit.size()<<endl; + + // cout<<i1<<" Stat "<<SubDId<<" ~~~ "<<inputABHit.size()<<endl; } //cout<<"hit size"<<inputHits.size()<<endl; - Sequence = Arbor(inputABHit, _cepc_thresholds); + Sequence = Arbor(inputABHit, _cepc_thresholds); + + cout<<"[YX debug - MarlinArbor] inputABHit.size() = "<<inputABHit.size()<<endl; + cout<<"[YX debug - MarlinArbor] Sequence.size() = "<<Sequence.size()<<endl; m_ArborToolLCIO->ClusterBuilding( branchCol, inputHits, Trees, 0 ); @@ -207,6 +348,9 @@ StatusCode MarlinArbor::execute() } MakeIsoHits(IsoHits, m_isohitcol); + + cout<<"[YX debug - MarlinArbor] End."<<endl; + return StatusCode::SUCCESS; } diff --git a/Reconstruction/PFA/Arbor/src/MarlinArbor.hh b/Reconstruction/PFA/Arbor/src/MarlinArbor.hh index e768ea3ff207d0dd91bbb8e1a5ddec1766190c8a..c68128af02d1b0c16588986dafa658a77d85f4ee 100644 --- a/Reconstruction/PFA/Arbor/src/MarlinArbor.hh +++ b/Reconstruction/PFA/Arbor/src/MarlinArbor.hh @@ -66,7 +66,7 @@ class MarlinArbor : public GaudiAlgorithm std::vector<std::string> _EcalPreShowerCollections; std::vector<std::string> _EcalCalCollections; std::vector<std::string> _HcalCalCollections; - std::vector<float> _cepc_thresholds; + std::vector<float> _cepc_thresholds; typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloType; @@ -76,6 +76,7 @@ class MarlinArbor : public GaudiAlgorithm Gaudi::Property<std::vector<std::string>> m_ecalReadoutNames{this, "ECALReadOutNames", {"EcalBarrelCollection", "EcalEndcapsCollection","EcalEndcapRingCollection"}, "Name of readouts"}; Gaudi::Property<std::vector<std::string>> m_hcalReadoutNames{this, "HCALReadOutNames", {"HcalBarrelCollection", "HcalEndcapsCollection","HcalEndcapRingCollection"}, "Name of readouts"}; + Gaudi::Property<bool> m_readLCIO{this, "ReadLCIO", true, "Read sim file with LCIO"}; std::map<std::string, std::string> m_col_readout_map; std::vector<CaloType*> _ecalCollections; @@ -83,34 +84,34 @@ class MarlinArbor : public GaudiAlgorithm std::vector<CaloType*> _calCollections; TTree *_outputTree; - std::string _treeFileName; - int _EH; + std::string _treeFileName; + int _EH; float _HitPos[3]; float _BushP[3]; - float _CloseDis; + float _CloseDis; float _HitEnergy; - int _CellSize; + int _CellSize; int _CaloTrackLengthCut; int _Num, _Seg, _eventNr; int numElements; - float _DHCALFirstThreshold; + float _DHCALFirstThreshold; float _InitLinkDisThreshold; - bool _DHCALSimuDigiMode; + bool _DHCALSimuDigiMode; bool _FlagInputSimHit; bool _FlagMutePhoton; bool _FlagMuteChargeParticle; bool _FlagMuteGarlicHits; - bool _FlagUseTrackerEndHit; + bool _FlagUseTrackerEndHit; std::string m_encoder_str; DataHandle<edm4hep::ClusterCollection> branchCol{"EHBushes",Gaudi::DataHandle::Writer, this}; DataHandle<edm4hep::CalorimeterHitCollection> m_isohitcol{"IsoHits",Gaudi::DataHandle::Writer, this}; - TH2F *_h1, *_h2, *_h7; - TH1F *_h3, *_h4, *_h5, *_h6; + TH2F *_h1, *_h2, *_h7; + TH1F *_h3, *_h4, *_h5, *_h6; std::ostream _output; float _HLayerCut; diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt b/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt index dade7622c8f3144db1a18916bd5af4852739efdd..9f250bd8cce7ff83c0397e3fac9ca7775ab9d466 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt +++ b/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt @@ -17,7 +17,7 @@ gaudi_add_module(GaudiPandora k4FWCore::k4FWCore ${PandoraSDK_LIBRARIES} ${LCContent_LIBRARIES} - ${CLHEP_LIBRARIES} + ${CLHEP_LIBRARIES} ${ROOT_LIBRARIES} ${LCIO_LIBRARIES} ${GEAR_LIBRARIES} diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp index bb63363fe1b6719569888bcdfeaa91d39a1a6481..503e802312d0ba15e64d4242e647f9bc840b4c03 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp @@ -1,7 +1,7 @@ /** - * + * * @brief Implementation of the calo hit creator class. - * + * * $Log: $ */ @@ -31,10 +31,10 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora m_settings(settings), m_pPandora(pPandora) { - m_encoder_str = ""; - m_encoder_str_MUON = ""; - m_encoder_str_LCal = ""; - m_encoder_str_LHCal = ""; + m_encoder_str = ""; + m_encoder_str_MUON = ""; + m_encoder_str_LCal = ""; + m_encoder_str_LHCal = ""; if(encoder_style==0) // LCIO style { m_encoder_str = "M:3,S-1:3,I:9,J:9,K-1:6"; @@ -44,7 +44,7 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora } IGearSvc* iSvc = 0; StatusCode sc = svcloc->service("GearSvc", iSvc, false); - if ( !sc ) + if ( !sc ) { throw "Failed to find GearSvc ..."; } @@ -55,7 +55,7 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora } if(m_settings.m_use_dd4hep_geo){ const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); - if(eCalBarrelExtension){ + if(eCalBarrelExtension){ m_eCalBarrelOuterZ = eCalBarrelExtension->extent[3]/dd4hep::mm; m_eCalBarrelInnerPhi0 = eCalBarrelExtension->inner_phi0/dd4hep::rad; m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry; @@ -72,9 +72,9 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0()); m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder()); } - //Get HCal Barrel extension by type, ignore plugs and rings + //Get HCal Barrel extension by type, ignore plugs and rings const dd4hep::rec::LayeredCalorimeterData * hCalBarrelExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::HADRONIC | dd4hep::DetType::BARREL),( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); - //Get HCal Endcap extension by type, ignore plugs and rings + //Get HCal Endcap extension by type, ignore plugs and rings const dd4hep::rec::LayeredCalorimeterData * hCalEndcapExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::HADRONIC | dd4hep::DetType::ENDCAP),( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); if(hCalBarrelExtension){ m_hCalBarrelOuterZ = hCalBarrelExtension->extent[3]/dd4hep::mm; @@ -113,7 +113,7 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora "Hcal_outer_polygon_order") != _GEAR->getHcalBarrelParameters().getIntKeys().end() ? _GEAR->getHcalBarrelParameters().getIntVal("Hcal_outer_polygon_order") : 0)); - const gear::LayerLayout &hCalBarrelLayerLayout(_GEAR->getHcalBarrelParameters().getLayerLayout()); + const gear::LayerLayout &hCalBarrelLayerLayout(_GEAR->getHcalBarrelParameters().getLayerLayout()); m_hCalBarrelLayerThickness = hCalBarrelLayerLayout.getThickness(hCalBarrelLayerLayout.getNLayers() - 1); } if(hCalEndcapExtension){ @@ -137,10 +137,10 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora m_hCalEndCapLayerThickness = hCalEndCapLayerLayout.getThickness(hCalEndCapLayerLayout.getNLayers() - 1); } - //Get Muon Barrel extension by type, ignore plugs and rings + //Get Muon Barrel extension by type, ignore plugs and rings const dd4hep::rec::LayeredCalorimeterData * muonBarrelExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::MUON | dd4hep::DetType::BARREL),( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); //fg: muon endcap is not used : - //Get Muon Endcap extension by type, ignore plugs and rings + //Get Muon Endcap extension by type, ignore plugs and rings // const dd4hep::rec::LayeredCalorimeterData * muonEndcapExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::MUON | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY ) ); if(muonBarrelExtension){ m_muonBarrelOuterZ = muonBarrelExtension->extent[3]/dd4hep::mm; @@ -191,7 +191,7 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora "Hcal_outer_polygon_order") != _GEAR->getHcalBarrelParameters().getIntKeys().end() ? _GEAR->getHcalBarrelParameters().getIntVal("Hcal_outer_polygon_order") : 0)); - const gear::LayerLayout &hCalBarrelLayerLayout(_GEAR->getHcalBarrelParameters().getLayerLayout()); + const gear::LayerLayout &hCalBarrelLayerLayout(_GEAR->getHcalBarrelParameters().getLayerLayout()); m_hCalBarrelLayerThickness = hCalBarrelLayerLayout.getThickness(hCalBarrelLayerLayout.getNLayers() - 1); m_hCalEndCapOuterR = (_GEAR->getHcalEndcapParameters().getExtent()[1]); m_hCalEndCapOuterZ = (_GEAR->getHcalEndcapParameters().getExtent()[3]); @@ -219,7 +219,7 @@ CaloHitCreator::~CaloHitCreator() pandora::StatusCode CaloHitCreator::CreateCaloHits(const CollectionMaps& collectionMaps) { - + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateECalCaloHits (collectionMaps)); PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateHCalCaloHits (collectionMaps)); PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateMuonCaloHits (collectionMaps)); @@ -268,7 +268,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col endcapLayers= &(PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )->layers); } else{ - barrelLayerLayout = &(_GEAR->getEcalBarrelParameters().getLayerLayout()); + barrelLayerLayout = &(_GEAR->getEcalBarrelParameters().getLayerLayout()); endcapLayerLayout = &(_GEAR->getEcalEndcapParameters().getLayerLayout()); } @@ -316,7 +316,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col caloHitParameters.m_hitType = pandora::ECAL; caloHitParameters.m_isDigital = false; caloHitParameters.m_layer = m_settings.m_use_dd4hep_decoder == false ? cellIdDecoder(pCaloHit)[layerCoding.c_str()] + 1 : m_decoder->get(pCaloHit->getCellID(), "layer");// from 0 to 29, 0 is preshower layer - int Stave = 0 ; + int Stave = 0 ; if (m_settings.m_use_dd4hep_decoder == false){ Stave = cellIdDecoder(pCaloHit)[ staveCoding]; } @@ -327,7 +327,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col } //std::cout<<"0Stave="<<Stave<<",0layer="<<caloHitParameters.m_layer.Get()<<std::endl; if (Stave<0) throw "wrong Stave"; - if (m_settings.m_use_preshower==false && caloHitParameters.m_layer.Get()<1) continue;//don't use preshower layer + if (m_settings.m_use_preshower==false && caloHitParameters.m_layer.Get()<1) continue;//don't use preshower layer //std::cout<<"Stave="<<Stave<<",layer="<<caloHitParameters.m_layer.Get()<<std::endl; caloHitParameters.m_isInOuterSamplingLayer = false; this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters); @@ -343,7 +343,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col else { if(m_settings.m_use_dd4hep_geo) this->GetEndCapCaloHitProperties(pCaloHit, *endcapLayers, caloHitParameters, absorberCorrection); - + else this->GetEndCapCaloHitProperties(pCaloHit, *endcapLayerLayout, caloHitParameters, absorberCorrection); caloHitParameters.m_hadronicEnergy = eCalToHadGeVEndCap * pCaloHit->getEnergy(); } @@ -444,26 +444,26 @@ pandora::StatusCode CaloHitCreator::CreateHCalCaloHits(const CollectionMaps& col //std::cout<<"HCAL layer="<<caloHitParameters.m_layer.Get()<<std::endl; caloHitParameters.m_isInOuterSamplingLayer = (this->GetNLayersFromEdge(pCaloHit) <= m_settings.m_nOuterSamplingLayers); this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters); - int Stave = 0 ; + int Stave = 0 ; if (m_settings.m_use_dd4hep_decoder == false){ Stave = cellIdDecoder(pCaloHit)[ staveCoding]; } else{ Stave = m_decoder->get(pCaloHit->getCellID(), "stave"); Stave = DD4hep2Lcio::CEPCv4::getHcalStave(Stave); - //Stave = Stave ==0 ? Stave+7 : Stave-1 ;//correct, same with LCIO + // Stave = Stave ==0 ? Stave+7 : Stave-1 ;//correct, same with LCIO /* 1 0 **** **** 2 * * 0 1 * * 7 * * * * - 3* * 7 ---> 2* * 6 + 3* * 7 ---> 2* * 6 * * * * - 4 * * 6 3 * * 5 - **** **** + 4 * * 6 3 * * 5 + **** **** 5 4 - - + + */ } @@ -551,7 +551,7 @@ pandora::StatusCode CaloHitCreator::CreateMuonCaloHits(const CollectionMaps& col } else{ plugLayerLayout= &(_GEAR->getYokePlugParameters().getLayerLayout()); - barrelLayerLayout = &(_GEAR->getYokeBarrelParameters().getLayerLayout()); + barrelLayerLayout = &(_GEAR->getYokeBarrelParameters().getLayerLayout()); endcapLayerLayout = &(_GEAR->getYokeEndcapParameters().getLayerLayout()); } @@ -570,7 +570,7 @@ pandora::StatusCode CaloHitCreator::CreateMuonCaloHits(const CollectionMaps& col //std::cout<<"Muon layer="<<caloHitParameters.m_layer.Get()<<std::endl; caloHitParameters.m_isInOuterSamplingLayer = true; this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters); - int Stave = 0 ; + int Stave = 0 ; if (m_settings.m_use_dd4hep_decoder == false){ Stave = cellIdDecoder(pCaloHit)[ staveCoding]; /* @@ -670,7 +670,7 @@ pandora::StatusCode CaloHitCreator::CreateLCalCaloHits(const CollectionMaps& col auto tmpCaloHit0 = pCaloHitCollection.at(0); - const gear::LayerLayout &endcapLayerLayout(_GEAR->getLcalParameters().getLayerLayout()); + const gear::LayerLayout &endcapLayerLayout(_GEAR->getLcalParameters().getLayerLayout()); ID_UTIL::CellIDDecoder<decltype(tmpCaloHit0)> cellIdDecoder(m_encoder_str_LCal); const std::string layerCodingString(m_encoder_str_LCal); @@ -896,7 +896,7 @@ void CaloHitCreator::GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *c for (unsigned int i = 0, iMax = layers.size(); i < iMax; ++i) { float absorberThickness((layers[i].inner_thickness - layers[i].sensitive_thickness/2.0 )/dd4hep::mm); - + if (i>0) absorberThickness += (layers[i-1].outer_thickness - layers[i-1].sensitive_thickness/2.0)/dd4hep::mm; @@ -996,7 +996,7 @@ void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *c nIntLengths += layers[physicalLayer-1].outer_nInteractionLengths; layerAbsorberThickness += (layers[physicalLayer-1].outer_thickness -layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm; } - + caloHitParameters.m_cellThickness = thickness; caloHitParameters.m_nCellRadiationLengths = nRadLengths; caloHitParameters.m_nCellInteractionLengths = nIntLengths; @@ -1008,11 +1008,11 @@ void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *c throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER); } absorberCorrection = 1.; - float absorberThickness_0 = 0; + float absorberThickness_0 = 0; for (unsigned int i = 0, iMax = layers.size(); i < iMax; ++i) { float absorberThickness((layers[i].inner_thickness - layers[i].sensitive_thickness/2.0 )/dd4hep::mm); - + if (i>0) absorberThickness += (layers[i-1].outer_thickness - layers[i-1].sensitive_thickness/2.0)/dd4hep::mm; @@ -1097,7 +1097,7 @@ int CaloHitCreator::GetNLayersFromEdge(const edm4hep::CalorimeterHit *const pCal float CaloHitCreator::GetMaximumRadius(const edm4hep::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const { - + const float pCaloHitPosition[3]={pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2]}; if (symmetryOrder <= 2) return std::sqrt((pCaloHitPosition[0] * pCaloHitPosition[0]) + (pCaloHitPosition[1] * pCaloHitPosition[1])); diff --git a/Utilities/DecoderHelper/CMakeLists.txt b/Utilities/DecoderHelper/CMakeLists.txt index 3c3f4ca9ec0729ade291be0428a6cb5348a61814..4b82d87aa0107f1eab214dc1ae1a7ba9bfdd7847 100644 --- a/Utilities/DecoderHelper/CMakeLists.txt +++ b/Utilities/DecoderHelper/CMakeLists.txt @@ -1,5 +1,5 @@ -gaudi_add_library(DecoderHelperLib +gaudi_add_library(DecoderHelperLib SOURCES src/DD4hep2Lcio.cc ) diff --git a/Utilities/DecoderHelper/include/DecoderHelper/DD4hep2Lcio.h b/Utilities/DecoderHelper/include/DecoderHelper/DD4hep2Lcio.h index a5a2210df3e94839bff9f1988311a01195520844..1205bc3a4bab0882749b264caafc336208313916 100644 --- a/Utilities/DecoderHelper/include/DecoderHelper/DD4hep2Lcio.h +++ b/Utilities/DecoderHelper/include/DecoderHelper/DD4hep2Lcio.h @@ -5,10 +5,12 @@ namespace DD4hep2Lcio { namespace CEPCv4 { int getEcalLayer(int dd4hep_layer); int getEcalStave(int dd4hep_stave); + int getEcalBarrelStave(int dd4hep_stave); + int getEcalEndcapStave(int dd4hep_stave); int getHcalLayer(int dd4hep_layer); int getHcalStave(int dd4hep_stave); int getMuonLayer(int dd4hep_layer); int getMuonStave(int dd4hep_stave); } } -#endif +#endif \ No newline at end of file diff --git a/Utilities/DecoderHelper/src/DD4hep2Lcio.cc b/Utilities/DecoderHelper/src/DD4hep2Lcio.cc index 30dba288134e03f2b89f5b8452f3c4e7dee24fe3..71999a83ab3d6fb6253e7a21bf2dbbc5f9e5ce13 100644 --- a/Utilities/DecoderHelper/src/DD4hep2Lcio.cc +++ b/Utilities/DecoderHelper/src/DD4hep2Lcio.cc @@ -7,6 +7,14 @@ int DD4hep2Lcio::CEPCv4::getEcalStave(int dd4hep_stave){ int lcio_stave = dd4hep_stave <=2 ? dd4hep_stave+5 : dd4hep_stave-3 ; return lcio_stave ; } +int DD4hep2Lcio::CEPCv4::getEcalBarrelStave(int dd4hep_stave){ + int lcio_stave = dd4hep_stave <=2 ? dd4hep_stave+5 : dd4hep_stave-3 ; + return lcio_stave ; +} +int DD4hep2Lcio::CEPCv4::getEcalEndcapStave(int dd4hep_stave){ + int lcio_stave = dd4hep_stave <=2 ? dd4hep_stave+1 : dd4hep_stave-3 ; + return lcio_stave ; +} int DD4hep2Lcio::CEPCv4::getHcalLayer(int dd4hep_layer){ return dd4hep_layer - 1 ; } @@ -18,19 +26,19 @@ int DD4hep2Lcio::CEPCv4::getHcalStave(int dd4hep_stave){ **** **** 2 * * 0 1 * * 7 * * * * - 3* * 7 ---> 2* * 6 + 3* * 7 ---> 2* * 6 * * * * - 4 * * 6 3 * * 5 - **** **** + 4 * * 6 3 * * 5 + **** **** 5 4 - - + + */ - return lcio_stave ; + return lcio_stave ; } int DD4hep2Lcio::CEPCv4::getMuonLayer(int dd4hep_layer){ return dd4hep_layer - 1 ; } int DD4hep2Lcio::CEPCv4::getMuonStave(int dd4hep_stave){ return 12 - dd4hep_stave ; -} +} \ No newline at end of file