diff --git a/Reconstruction/PFA/Arbor/CMakeLists.txt b/Reconstruction/PFA/Arbor/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..eb4468bd61180dea7e8e9dba8d3c0bf9e41ace6c --- /dev/null +++ b/Reconstruction/PFA/Arbor/CMakeLists.txt @@ -0,0 +1,34 @@ + +# Modules +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/MarlinArbor.cc + src/BushConnect.cc + LINK EventSeeder + GearSvc + DataHelperLib + DetInterface + Gaudi::GaudiKernel + k4FWCore::k4FWCore + ${LCContent_LIBRARIES} + ${CLHEP_LIBRARIES} + ${ROOT_LIBRARIES} + ${LCIO_LIBRARIES} + ${GEAR_LIBRARIES} + ${DD4hep_COMPONENT_LIBRARIES} + EDM4HEP::edm4hep EDM4HEP::edm4hepDict +) + + +install(TARGETS Arbor + EXPORT CEPCSWTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) + diff --git a/Reconstruction/PFA/Arbor/src/Arbor.cc b/Reconstruction/PFA/Arbor/src/Arbor.cc new file mode 100644 index 0000000000000000000000000000000000000000..c51a6fae8f51d261d6a4830b8866922d97fcdb25 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/Arbor.cc @@ -0,0 +1,923 @@ +#include "ArborHit.h" +#include "Arbor.h" +#include <TTree.h> +#include <algorithm> +#include <TMath.h> +#include "KDTreeLinkerAlgoT.h" +#include <unordered_map> + +//#define DEBUG + +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> IsoHitsIndex; +std::vector<int> SimpleSeedHitsIndex; +std::vector<int> StarSeedHitsIndex; +std::map<int, int> HitsFlag; // Tell Hits type; + +/* +HitLinkMap BackLinksMap; +HitLinkMap alliterBackLinksMap; +*/ +HitLinkMap IterBackLinks; + +linkcoll Links; +linkcoll InitLinks; +linkcoll IterInputLinks; +linkcoll alliterlinks; +linkcoll links_debug; +linkcoll IterLinks; +linkcoll IterLinks_1; +branchcoll LengthSortBranchCollection; +branchcoll Trees; + +int NHits = 0; + +void init() { + cleanedHits.clear(); + LeafHitsIndex.clear(); + JointHitsIndex.clear(); + StarJointHitsIndex.clear(); + IsoHitsIndex.clear(); + SimpleSeedHitsIndex.clear(); + StarSeedHitsIndex.clear(); + HitsFlag.clear(); + + Links.clear(); + InitLinks.clear(); + IterLinks_1.clear(); + IterLinks.clear(); + IterInputLinks.clear(); + IterBackLinks.clear(); + LengthSortBranchCollection.clear(); + Trees.clear(); + alliterlinks.clear(); + links_debug.clear(); +} + +void HitsCleaning( std::vector<ArborHit> inputHits ) +{ + cleanedHits = inputHits; //Cannot Really do much things here. Mapping before calling + NHits = cleanedHits.size(); +} + +void HitsClassification( linkcoll inputLinks ) +{ + int NLinks = inputLinks.size(); + + LeafHitsIndex.clear(); + JointHitsIndex.clear(); + StarJointHitsIndex.clear(); + IsoHitsIndex.clear(); + SimpleSeedHitsIndex.clear(); + StarSeedHitsIndex.clear(); + + std::pair<int, int> a_link; + + int BeginIndex[NHits]; + int EndIndex[NHits]; + + for(int i0 = 0; i0 < NHits; i0++) + { + BeginIndex[i0] = 0; + EndIndex[i0] = 0; + } + + for(int j0 = 0; j0 < NLinks; j0++) + { + BeginIndex[ (inputLinks[j0].first) ] ++; + EndIndex[ (inputLinks[j0].second) ] ++; + } + + for(int i1 = 0; i1 < NHits; i1++) + { + if(BeginIndex[i1] == 0 && EndIndex[i1] == 1) + { + LeafHitsIndex.push_back(i1); + HitsFlag[i1] = 5; + } + else if(BeginIndex[i1] == 1 && EndIndex[i1] == 1) + { + JointHitsIndex.push_back(i1); + HitsFlag[i1] = 4; + } + else if(BeginIndex[i1] > 1 && EndIndex[i1] == 1) + { + StarJointHitsIndex.push_back(i1); + HitsFlag[i1] = 3; + } + else if(BeginIndex[i1] == 1 && EndIndex[i1] == 0) + { + SimpleSeedHitsIndex.push_back(i1); + HitsFlag[i1] = 2; + } + else if(BeginIndex[i1] > 1 && EndIndex[i1] == 0) + { + StarSeedHitsIndex.push_back(i1); + HitsFlag[i1] = 1; + } + else if(BeginIndex[i1] == 0 && EndIndex[i1] == 0) + { + IsoHitsIndex.push_back(i1); + HitsFlag[i1] = 0; + } + else + { + cout<<"WARNING: UNCLASSIFIED HITS, Begin Index: "<<BeginIndex[i1]<<", End Index: "<<EndIndex[i1]<<endl; + } + } + +#ifdef DEBUG + cout<<"Verification of Hits Classification: "<<endl; + 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; +#endif +} + +linkcoll LinkClean( std::vector<ArborHit> allhits, linkcoll alllinks ) +{ + linkcoll cleanedlinks; + + int NLinks = alllinks.size(); + int Ncurrhitlinks = 0; + int MinAngleIndex = -1; + float MinAngle = 1E6; + float tmpOrder = 0; + float DirAngle = 0; + + std::pair<int, int> SelectedPair; + + TVector3 PosA, PosB, PosDiffAB; + + std::vector< std::vector<int> > LinkHits; + LinkHits.clear(); + for(int s1 = 0; s1 < NHits; s1++) + { + std::vector<int> hitlink; + for(int t1 = 0; t1 < NLinks; t1++) + { + if(alllinks[t1].second == s1) + { + hitlink.push_back(alllinks[t1].first); + } + } + LinkHits.push_back(hitlink); + } + + for(int i1 = 0; i1 < NHits; i1++) + { + PosB = cleanedHits[i1].GetPosition(); + MinAngleIndex = -10; + MinAngle = 1E6; + + std::vector<int> currhitlink = LinkHits[i1]; + + Ncurrhitlinks = currhitlink.size(); + + for(int k1 = 0; k1 < Ncurrhitlinks; k1++) + { + PosA = cleanedHits[ currhitlink[k1] ].GetPosition(); + DirAngle = (PosA + PosB).Angle(PosB - PosA); //Replace PosA + PosB with other order parameter ~ reference direction + tmpOrder = (PosB - PosA).Mag() * (DirAngle + 0.1); + if( tmpOrder < MinAngle ) // && DirAngle < 2.5 ) + { + MinAngleIndex = currhitlink[k1]; + MinAngle = tmpOrder; + } + } + + if(MinAngleIndex > -0.5) + { + SelectedPair.first = MinAngleIndex; + SelectedPair.second = i1; + cleanedlinks.push_back(SelectedPair); + } + } + +#ifdef DEBUG + cout<<"NStat "<<NHits<<" : "<<NLinks<<" InitLinks "<<InitLinks.size()<<endl; +#endif + return cleanedlinks; +} + +void BuildInitLink( std::vector<float>Thresholds ) +{ + 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; + + for(int i0 = 0; i0 < NHits; ++i0 ) + { + const auto& hit = cleanedHits[i0].GetPosition(); + nodes.emplace_back(i0,(float)hit.X(),(float)hit.Y(),(float)hit.Z()); + if( i0 == 0 ) + { + minpos[0] = hit.X(); minpos[1] = hit.Y(); minpos[2] = hit.Z(); + maxpos[0] = hit.X(); maxpos[1] = hit.Y(); maxpos[2] = hit.Z(); + } + else + { + minpos[0] = std::min((float)hit.X(),minpos[0]); + minpos[1] = std::min((float)hit.Y(),minpos[1]); + minpos[2] = std::min((float)hit.Z(),minpos[2]); + maxpos[0] = std::max((float)hit.X(),maxpos[0]); + maxpos[1] = std::max((float)hit.Y(),maxpos[1]); + maxpos[2] = std::max((float)hit.Z(),maxpos[2]); + } + } + KDTreeCube kdXYZCube(minpos[0],maxpos[0], + minpos[1],maxpos[1], + minpos[2],maxpos[2]); + kdtree.build(nodes,kdXYZCube); + nodes.clear(); + + Links.clear(); //all tmp links + + TVector3 PosA, PosB, PosDiffAB, PosDiffBA; + int NLayer_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 MagB = 0; + float ECCorr = 0; + int FlagTrkPS = 0; + int FlagPSEE = 0; + int FlagEH = 0; + int FlagHH = 0; + int FlagStaveSame = 0; + int FlagStaveDiff = 0; + + for(int i0 = 0; i0 < NHits; i0++) + { + + found.clear(); + + PosA = cleanedHits[i0].GetPosition(); + NLayer_A = cleanedHits[i0].GetLayer(); + NStave_A = cleanedHits[i0].GetStave(); + SubD_A = cleanedHits[i0].GetSubD(); //SubD_Index, 0 = track endpoint, 1 = Ecal, 2 = Hcal, 3 = EcalPS + Depth_A = cleanedHits[i0].GetDepth(); + + 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); + const float zplus(PosA.Z() + side), zminus(PosA.Z() - side); + const float xmin(std::min(xplus,xminus)), xmax(std::max(xplus,xminus)); + const float ymin(std::min(yplus,yminus)), ymax(std::max(yplus,yminus)); + const float zmin(std::min(zplus,zminus)), zmax(std::max(zplus,zminus)); + KDTreeCube searchcube( xmin, xmax, + ymin, ymax, + zmin, zmax ); + kdtree.search(searchcube,found); + + 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(); + 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; + if( (fabs(PosA.Z()) - 2600 )*(fabs(PosB.Z()) - 2600) < 0 ) + ECCorr = 40; + + // 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 * 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; + + if( FlagTrkPS || FlagPSEE || FlagEH || FlagHH ) + { + std::pair<int, int> a_Link; + MagA = PosA.Mag(); + MagB = PosB.Mag(); + + if( NStave_A != NStave_B || ( NLayer_A == 0 && NLayer_B != 0 ) || ( NLayer_B == 0 && NLayer_A != 0 ) ) + FlagStaveDiff = 1; + else if( NLayer_A != NLayer_B && NStave_A == NStave_B ) + FlagStaveSame = 1; + + if( FlagStaveDiff || FlagStaveSame || FlagTrkPS || FlagEH ) + { + if( MagA > MagB ) + { + a_Link.first = found[j0].data; + a_Link.second = i0; + } + else + { + a_Link.first = i0; + a_Link.second = found[j0].data; + } + Links.push_back(a_Link); + } + } + } + } + links_debug = Links; +} + +void LinkIteration( int time ) //Energy corrections, semi-local correction +{ + + 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; + + TVector3 RefDir[NHits]; + int Nin_hit[NHits]; + int Nout_hit[NHits]; + + for(int i0 = 0; i0 < NHits; i0++) + { + const auto& hit = cleanedHits[i0].GetPosition(); + RefDir[i0] = 1.0/hit.Mag() * hit; + Nin_hit[i0] = 0; + Nout_hit[i0] = 0; + + nodes.emplace_back(i0,(float)hit.X(),(float)hit.Y(),(float)hit.Z()); + if( i0 == 0 ) + { + minpos[0] = hit.X(); minpos[1] = hit.Y(); minpos[2] = hit.Z(); + maxpos[0] = hit.X(); maxpos[1] = hit.Y(); maxpos[2] = hit.Z(); + } + else + { + minpos[0] = std::min((float)hit.X(),minpos[0]); + minpos[1] = std::min((float)hit.Y(),minpos[1]); + minpos[2] = std::min((float)hit.Z(),minpos[2]); + maxpos[0] = std::max((float)hit.X(),maxpos[0]); + maxpos[1] = std::max((float)hit.Y(),maxpos[1]); + maxpos[2] = std::max((float)hit.Z(),maxpos[2]); + } + } + + KDTreeCube kdXYZCube(minpos[0],maxpos[0], + minpos[1],maxpos[1], + minpos[2],maxpos[2]); + kdtree.build(nodes,kdXYZCube); + nodes.clear(); + + IterBackLinks.clear(); + + if(time == 1) + { + IterLinks_1.clear(); + alliterlinks = InitLinks; + } + else + { + IterLinks.clear(); + alliterlinks = IterLinks_1; + } + int NcurrLinks = alliterlinks.size(); + + 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 SubD_B = 0; + int FlagEE = 0; + int FlagHH = 0; + int AngleAccIndex = 0; + float DisAB = 0; + float MagA = 0; + float MagB = 0; + 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; + + for(int i = 0; i < NHits; i++) + { + hitPos = cleanedHits[i].GetPosition(); + RefDir[i] = 1.0/hitPos.Mag() * hitPos; + Nin_hit[i] = 0; + Nout_hit[i] = 0; + } + + // std::vector<int> ---> all the hits linked to this hit + std::vector< std::vector<int> > hitLinksArray; + 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(); + RefDir[currlink.first] += 4*linkDir; //Weights... might be optimized... + RefDir[currlink.second] += 6*linkDir; + Nin_hit[currlink.first] ++; + Nout_hit[currlink.second] ++; + hitLinksArray[currlink.first].push_back(currlink.second); + } + + //Classification of cases: Branch to IsoHit, Branch, Seed + for(int i1 = 0; i1 < NHits; i1++) + { + found.clear(); + RefDir[i1] *= 1.0/RefDir[i1].Mag(); + PosA = cleanedHits[i1].GetPosition(); + NLayer_A = cleanedHits[i1].GetLayer(); + NStave_A = cleanedHits[i1].GetStave(); + SubD_A = cleanedHits[i1].GetSubD(); + + const float side = 200 + 100*time; + const float xplus(PosA.X() + side), xminus(PosA.X() - side); + const float yplus(PosA.Y() + side), yminus(PosA.Y() - side); + const float zplus(PosA.Z() + side), zminus(PosA.Z() - side); + const float xmin(std::min(xplus,xminus)), xmax(std::max(xplus,xminus)); + const float ymin(std::min(yplus,yminus)), ymax(std::max(yplus,yminus)); + const float zmin(std::min(zplus,zminus)), zmax(std::max(zplus,zminus)); + KDTreeCube searchcube( xmin, xmax, + ymin, ymax, + zmin, zmax ); + kdtree.search(searchcube,found); + + 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; + + bool isConnected = false; + + std::vector<int>& linksOfHitI = hitLinksArray[i1]; + std::vector<int>& linksOfHitJ = hitLinksArray[j1]; + + + for(auto it=linksOfHitI.begin(); it!=linksOfHitI.end(); ++it) { + int hitConnected = *it; + + if(hitConnected==(int)j1) { + isConnected = true; + break; + } + } + + // a piece of code copied from upper lines :( + if(!isConnected) { + for(auto it=linksOfHitJ.begin(); it!=linksOfHitJ.end(); ++it) { + int hitConnected = *it; + + if(hitConnected==(int)i1) { + isConnected = true; + break; + } + } + } + + + if( FlagNoJoint == 0 && FlagNoIso != 0 ) + { + 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(); + + 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; + else if( PosDiffBA.Angle(RefDir[i1]) < 0.6/time ) + AngleAccIndex = 3; + else if( PosDiffBA.Angle(RefDir[j1]) < 0.6/time ) + AngleAccIndex = 4; + + if(AngleAccIndex) + { + if(SubD_A == 1 && SubD_B == 1 && DisAB < 15*(time+1) ) // && DisAB > 15*time) + FlagEE = 1; + if(SubD_A > 2 && SubD_B == 2 && DisAB < 50*(time+1) ) // && DisAB > 50*time ) + FlagHH = 1; + } + + if(FlagEE || FlagHH) + { + MagA = PosA.Mag(); + 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 ) + { + a_Link.first = found[j1].data; + a_Link.second = i1; + } + else if( MagA < MagB && AngleAccIndex > 2 ) + { + a_Link.first = i1; + a_Link.second = found[j1].data; + } + alliterlinks.push_back(a_Link); + hitLinksArray[a_Link.first].push_back(a_Link.second); + } + else if( NLayer_A != NLayer_B && NStave_A == NStave_B) + { + if( MagA > MagB && AngleAccIndex < 3 ) + { + a_Link.first = found[j1].data; + a_Link.second = i1; + } + else if( MagA < MagB && AngleAccIndex > 2 ) + { + a_Link.first = i1; + a_Link.second = found[j1].data; + } + alliterlinks.push_back(a_Link); + hitLinksArray[a_Link.first].push_back(a_Link.second); + } + // } + } + } + } + } + } + + //Reusage of link iteration codes? + + int NLinks = alliterlinks.size(); + int MinAngleIndex = -10; + int Ncurrhitlinks = 0; + float MinAngle = 1E6; + float tmpOrder = 0; + float DirAngle = 0; + std::pair<int, int> SelectedPair; + + std::vector< std::vector<int> > LinkHits; + LinkHits.clear(); + for(int s1 = 0; s1 < NHits; s1++) + { + std::vector<int> hitlink; + for(int t1 = 0; t1 < NLinks; t1++) + { + if(alliterlinks[t1].second == s1) + { + hitlink.push_back(alliterlinks[t1].first); + } + } + LinkHits.push_back(hitlink); + } + + for(int i2 = 0; i2 < NHits; i2++) + { + PosB = cleanedHits[i2].GetPosition(); + MinAngleIndex = -10; + MinAngle = 1E6; + + std::vector<int> currhitlink = LinkHits[i2]; + + Ncurrhitlinks = currhitlink.size(); + + for(int j2 = 0; j2 < Ncurrhitlinks; j2++) + { + PosA = cleanedHits[ currhitlink[j2] ].GetPosition(); + DirAngle = (RefDir[i2]).Angle(PosA - PosB); + tmpOrder = (PosB - PosA).Mag() * (DirAngle + 0.6); + if(tmpOrder < MinAngle) // && DirAngle < 1.0) + { + MinAngleIndex = currhitlink[j2]; + MinAngle = tmpOrder; + } + } + + if(MinAngleIndex > -0.5) + { + SelectedPair.first = MinAngleIndex; + SelectedPair.second = i2; + if(SelectedPair.first == SelectedPair.second) + { + cout<<"WTTTTTTTTFFFFFFFFFFFFFF"<<endl; + continue; + } + if(time == 1) + { + IterLinks_1.push_back(SelectedPair); + } + else + { + IterLinks.push_back(SelectedPair); + IterBackLinks.emplace(SelectedPair.second,SelectedPair.first); + } + } + } + +#ifdef DEBUG + cout<<"Init-Iter Size "<<InitLinks.size()<<" : "<<IterLinks.size()<<endl; +#endif +} + +void BranchBuilding(float SeedThreshold) +{ + + int NLinks = IterLinks.size(); + int NBranches = 0; + std::map <int, int> HitBeginIndex; + std::map <int, int> HitEndIndex; + std::vector< std::vector<int> > InitBranchCollection; + std::vector< std::vector<int> > PrunedInitBranchCollection; + std::vector< std::vector<int> > TmpBranchCollection; + TVector3 PosA, PosB; + + for(int i1 = 0; i1 < NHits; i1++) + { + HitBeginIndex[i1] = 0; + HitEndIndex[i1] = 0; + } + + for(int j1 = 0; j1 < NLinks; j1++) + { + HitBeginIndex[ (IterLinks[j1].first) ] ++; + HitEndIndex[ (IterLinks[j1].second) ] ++; + } + + int iterhitindex = 0; + int LL = 0; //Uplimt to be set to twice the total layer thickness... + + for(int i2 = 0; i2 < NHits; i2++) + { + if(HitEndIndex[i2] > 1) + cout<<"WARNING OF INTERNAL LOOP with more than 1 link stopped at the same Hit"<<endl; + + if(HitBeginIndex[i2] == 0 && HitEndIndex[i2] == 1) //EndPoint + { + NBranches ++; + std::vector<int> currBranchhits; //array of indexes + + iterhitindex = i2; + currBranchhits.push_back(i2); + LL = 0; + + while(HitEndIndex[iterhitindex] != 0 && LL < 300) // 100 put by hand + { + /* + for(int j2 = 0; j2 < NLinks; j2++) + { + std::pair<int, int> PairIterator = IterLinks[j2]; + if( PairIterator.second == iterhitindex && std::find(currBranchhits.begin(), currBranchhits.end(), PairIterator.first) == currBranchhits.end() ) + { + currBranchhits.push_back(PairIterator.first); + iterhitindex = PairIterator.first; + break; + } + } + */ + auto iterlink_range = IterBackLinks.equal_range(iterhitindex); + assert( std::distance(iterlink_range.first,iterlink_range.second) == 1 ); + iterhitindex = iterlink_range.first->second; + currBranchhits.push_back(iterhitindex); + LL++; + } + InitBranchCollection.push_back(std::move(currBranchhits) ); + } + } + + PrunedInitBranchCollection.resize(InitBranchCollection.size()); + + std::vector<float> BranchSize; + std::vector<float> cutBranchSize; + std::vector<int> SortedBranchIndex; + std::vector<int> SortedcutBranchIndex; + std::vector<int> currBranch; + std::vector<int> iterBranch; + std::vector<int> touchedHits; + std::vector<bool> touchedHitsMap(NHits,false); + std::vector<bool> seedHitsMap(NHits,false); + std::vector<unsigned> seedHits; + std::vector<int> leadingbranch; + + 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; + bool needInitPosMinMax = true; + + HitLinkMap seedToBranchesMap; + std::unordered_map<int,int> branchToSeed; + std::map<branch, int> SortedBranchToOriginal; + SortedBranchToOriginal.clear(); + branchToSeed.clear(); + + int currBranchSize = 0; + int currHit = 0; + + for(int i3 = 0; i3 < NBranches; i3++) + { + currBranch = InitBranchCollection[i3]; + BranchSize.push_back( float(currBranch.size()) ); + } + + SortedBranchIndex = std::move( SortMeasure(BranchSize, 1) ); + + for(int i4 = 0; i4 < NBranches; i4++) + { + currBranch = InitBranchCollection[SortedBranchIndex[i4]]; + currBranchSize = currBranch.size(); + iterBranch.clear(); + + for(int j4 = 0; j4 < currBranchSize; j4++) + { + currHit = currBranch[j4]; + + if( !touchedHitsMap[currHit] ) + { + iterBranch.push_back(currHit); + touchedHitsMap[currHit] = true; + } + } + const auto theseed = currBranch[currBranchSize - 1]; + SortedBranchToOriginal[iterBranch] = theseed; //Map to seed... + branchToSeed.emplace(SortedBranchIndex[i4],theseed); + seedToBranchesMap.emplace(theseed,SortedBranchIndex[i4]); // map seed to branches + if( !seedHitsMap[theseed] ) + { + seedHitsMap[theseed] = true; + const auto& hit = cleanedHits[theseed].GetPosition(); + nodes.emplace_back(theseed,(float)hit.X(),(float)hit.Y(),(float)hit.Z()); + if( needInitPosMinMax ) { + needInitPosMinMax = false; + minpos[0] = hit.X(); minpos[1] = hit.Y(); minpos[2] = hit.Z(); + maxpos[0] = hit.X(); maxpos[1] = hit.Y(); maxpos[2] = hit.Z(); + } else { + minpos[0] = std::min((float)hit.X(),minpos[0]); + minpos[1] = std::min((float)hit.Y(),minpos[1]); + minpos[2] = std::min((float)hit.Z(),minpos[2]); + maxpos[0] = std::max((float)hit.X(),maxpos[0]); + maxpos[1] = std::max((float)hit.Y(),maxpos[1]); + maxpos[2] = std::max((float)hit.Z(),maxpos[2]); + } + } + + TmpBranchCollection.push_back(iterBranch); + cutBranchSize.push_back( float(iterBranch.size()) ); + PrunedInitBranchCollection[SortedBranchIndex[i4]] = std::move(iterBranch); + } + + SortedcutBranchIndex = std::move( SortMeasure(cutBranchSize, 1) ); + + for(int i6 = 0; i6 < NBranches; i6++) + { + currBranch.clear(); + currBranch = TmpBranchCollection[ SortedcutBranchIndex[i6]]; + LengthSortBranchCollection.push_back(currBranch);; + } + + std::vector<bool> link_helper(NBranches*NBranches,false); + TVector3 DisSeed; + KDTreeCube kdXYZCube(minpos[0],maxpos[0], + minpos[1],maxpos[1], + minpos[2],maxpos[2]); + kdtree.build(nodes,kdXYZCube); + nodes.clear(); + + QuickUnion qu(NBranches); + + for(int i7 = 0; i7 < NBranches; i7++) + { + auto SeedIndex_A = branchToSeed[i7]; + auto shared_branches = seedToBranchesMap.equal_range(SeedIndex_A); + + for( auto itr = shared_branches.first; itr != shared_branches.second; ++itr ) { + const auto foundSortedIdx = itr->second; + if( foundSortedIdx <= (unsigned)i7 ) continue; + if( link_helper[NBranches*i7 + foundSortedIdx] || + link_helper[NBranches*foundSortedIdx + i7] ) continue; + if( !qu.connected(i7,foundSortedIdx) ) { + qu.unite(i7,foundSortedIdx); + } + link_helper[NBranches*i7 + foundSortedIdx] = true; + link_helper[NBranches*foundSortedIdx + i7] = true; + } + + const auto& seedpos = cleanedHits[SeedIndex_A].GetPosition(); + const float side = 10*SeedThreshold; + //cout<<"Arbor SeedThreshold"<<SeedThreshold<<endl; + const float xplus(seedpos.X() + side), xminus(seedpos.X() - side); + const float yplus(seedpos.Y() + side), yminus(seedpos.Y() - side); + const float zplus(seedpos.Z() + side), zminus(seedpos.Z() - side); + const float xmin(std::min(xplus,xminus)), xmax(std::max(xplus,xminus)); + const float ymin(std::min(yplus,yminus)), ymax(std::max(yplus,yminus)); + const float zmin(std::min(zplus,zminus)), zmax(std::max(zplus,zminus)); + KDTreeCube searchcube( xmin, xmax, + ymin, ymax, + zmin, zmax ); + found.clear(); + kdtree.search(searchcube,found); + for(unsigned j7 = 0; j7 < found.size(); j7++) { + DisSeed = seedpos - cleanedHits[ found[j7].data ].GetPosition(); + + if( DisSeed.Mag() < SeedThreshold ) + { + auto seed_branches = seedToBranchesMap.equal_range(found[j7].data); + for( auto itr = seed_branches.first; itr != seed_branches.second; ++itr ){ + const auto foundSortedIdx = itr->second; + if( foundSortedIdx <= (unsigned)i7 ) continue; + if( link_helper[NBranches*i7 + foundSortedIdx] || + link_helper[NBranches*foundSortedIdx + i7] ) continue; + if( !qu.connected(i7,foundSortedIdx) ) { + qu.unite(i7,foundSortedIdx); + } + link_helper[NBranches*i7 + foundSortedIdx] = true; + link_helper[NBranches*foundSortedIdx + i7] = true; + } + } + } + } + + Trees.clear(); + std::unordered_map<unsigned,branch> merged_branches(qu.count()); + Trees.reserve(qu.count()); + + for( unsigned i = 0; i < (unsigned)NBranches; ++i ) { + unsigned root = qu.find(i); + const auto& branch = PrunedInitBranchCollection[i]; + auto& merged_branch = merged_branches[root]; + merged_branch.insert(merged_branch.end(),branch.begin(),branch.end()); + } + + unsigned total_hits = 0; + for( auto& final_branch : merged_branches ) { + total_hits += final_branch.second.size(); + Trees.push_back(std::move(final_branch.second)); + } +} + +void BushMerging() +{ + cout<<"Merging branch"<<endl; +} + +//HitThresholds: EEThreshold, EHThreshold, HHThreshold, EE_Seed_Threshold; +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); + } + + init(); + HitsCleaning(inputHits); + BuildInitLink(Thresholds); + InitLinks = LinkClean( cleanedHits, Links ); + + LinkIteration(1); + LinkIteration(2); + HitsClassification(IterLinks); + BranchBuilding(Thresholds[3]); + + return LengthSortBranchCollection; +} diff --git a/Reconstruction/PFA/Arbor/src/Arbor.h b/Reconstruction/PFA/Arbor/src/Arbor.h new file mode 100644 index 0000000000000000000000000000000000000000..d7678cf7cdbf9704e1b38aac03da26bbf4ec4a57 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/Arbor.h @@ -0,0 +1,37 @@ +#ifndef _Arbor_hh_ +#define _Arbor_hh_ + +#include "ArborHit.h" +#include <string> +#include <iostream> +#include <TVector3.h> +#include "ArborTool.h" + +void init(); + +void HitsCleaning( std::vector<ArborHit> inputHits ); + +void HitsClassification( linkcoll inputLinks ); + +void BuildInitLink(std::vector<float> Thresholds); + +void LinkIteration(int time); + +void BranchBuilding(float SeedThreshold); + +branchcoll Arbor( std::vector<ArborHit>, std::vector<float> Thresholds ); + +/* + int NLayer_A, NStave_A, SubD_A; + int NLayer_B, NStave_B, SubD_B; + float MagA, MagB, Depth_A, Depth_B, ECCorr, DisAB; + int FlagTrkPS, FlagEH; + int FlagPSEE, FlagHH; + int FlagStaveSame; + int FlagStaveDiff; + TVector3 PosA, PosB, PosDiffAB, PosDiffBA, linkDir; +*/ + +#endif + + diff --git a/Reconstruction/PFA/Arbor/src/ArborHit.cc b/Reconstruction/PFA/Arbor/src/ArborHit.cc new file mode 100644 index 0000000000000000000000000000000000000000..712928c84da6e3b5ecddbfeb02884968f9484f75 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborHit.cc @@ -0,0 +1,18 @@ +#include "ArborHit.h" +#include <iostream> +#include <vector> + +ArborHit::ArborHit( TVector3 hitPos, int hitLayer, float hitTime, float depth, int stave, int subD ) +{ + setHit( hitPos, hitLayer, hitTime, depth, stave, subD ); +} + +void ArborHit::setHit(TVector3 hitPos, int hitLayer, float hitTime, float depth, int stave, int subD ) +{ + m_hitTime = hitTime; + m_hitLayer = hitLayer; + m_hitPos = hitPos; + m_subD = subD; + m_stave = stave; + m_depth = depth; +} diff --git a/Reconstruction/PFA/Arbor/src/ArborHit.h b/Reconstruction/PFA/Arbor/src/ArborHit.h new file mode 100644 index 0000000000000000000000000000000000000000..aacf7d951f43f2615782c8fe092dd1d4b2b58acc --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborHit.h @@ -0,0 +1,49 @@ +#ifndef _ARBORHIT_H_ +#define _ARBORHIT_H_ + +#include "TVector3.h" +#include <iostream> + +class ArborHit +{ + float m_hitTime; + float m_depth; + int m_hitLayer; + TVector3 m_hitPos; + int m_subD; + int m_stave; + + public: + + // ArborHit(); + ArborHit( TVector3 hitPos, int hitLayer, float hitTime, float depth, int stave, int subD ); + + void setHit( TVector3 hitPos, int hitLayer, float hitTime, float depth, int stave, int subD ); + + float GetTime() + { + return m_hitTime; + } + int GetLayer() + { + return m_hitLayer; + } + TVector3 GetPosition() + { + return m_hitPos; + } + int GetSubD() + { + return m_subD; + } + int GetStave() + { + return m_stave; + } + float GetDepth() + { + return m_depth; + } +}; + +#endif diff --git a/Reconstruction/PFA/Arbor/src/ArborTool.cc b/Reconstruction/PFA/Arbor/src/ArborTool.cc new file mode 100644 index 0000000000000000000000000000000000000000..73e06ee63fceee75201bd4aacd01de4288482897 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborTool.cc @@ -0,0 +1,183 @@ +#include "ArborTool.h" +#include <TMath.h> + +using namespace std; + +std::vector<int>SortMeasure( std::vector<float> Measure, int ControlOrderFlag ) +{ + + std::vector<int> objindex; + int Nobj = Measure.size(); + + for(int k = 0; k < Nobj; k++) + { + objindex.push_back(k); + } + + int FlagSwapOrder = 1; + float SwapMeasure = 0; + int SwapIndex = 0; + + for(int i = 0; i < Nobj && FlagSwapOrder; i++) + { + FlagSwapOrder = 0; + for(int j = 0; j < Nobj - 1; j++) + { + if((Measure[j] < Measure[j+1] && ControlOrderFlag) || (Measure[j] > Measure[j+1] && !ControlOrderFlag) ) + { + FlagSwapOrder = 1; + SwapMeasure = Measure[j]; + Measure[j] = Measure[j+1]; + Measure[j+1] = SwapMeasure; + + SwapIndex = objindex[j]; + objindex[j] = objindex[j+1]; + objindex[j+1] = SwapIndex; + } + } + } + + return objindex; +} + +TMatrixF MatrixSummarize( TMatrixF inputMatrix ) +{ + + int Nrow = inputMatrix.GetNrows(); + int Ncol = inputMatrix.GetNcols(); + + TMatrixF tmpMatrix(Nrow, Ncol); + + for(int i0 = 0; i0 < Nrow; i0 ++) + { + for(int j0 = i0; j0 < Ncol; j0 ++) + { + // if( fabs(inputMatrix(i0, j0) - 2) < 0.2 || fabs(inputMatrix(i0, j0) - 10 ) < 0.2 ) //Case 2, 3: Begin-End connector + if(inputMatrix(i0, j0)) + { + tmpMatrix(i0, j0) = 1; + tmpMatrix(j0, i0) = 1; + } + else + { + tmpMatrix(i0, j0) = 0; + tmpMatrix(j0, i0) = 0; + } + } + } + + int PreviousLinks = -1; + int CurrentLinks = 0; + int symloopsize = 0; + vector <int> Indirectlinks; + int tmpI(0); + int tmpJ(0); + + // cout<<"Matrix Type: "<<Nrow<<" * "<<Ncol<<endl; + + if( Nrow == Ncol ) + { + + while( CurrentLinks > PreviousLinks ) + { + PreviousLinks = 0; + CurrentLinks = 0; + + for(int i = 0; i < Nrow; i ++) + { + for(int j = 0; j < Ncol; j ++) + { + if( tmpMatrix(i, j) > 0.1 ) + PreviousLinks ++; + } + } + + for(int k = 0; k < Nrow; k ++) + { + for(int l = 0; l < Ncol; l ++) + { + if( tmpMatrix(k, l) > 0.1) Indirectlinks.push_back(l); + } + symloopsize = Indirectlinks.size(); + + for(int l1(0); l1 < symloopsize; l1 ++) + { + tmpI = Indirectlinks[l1]; + for(int m1=l1 + 1; m1 < symloopsize; m1++) + { + tmpJ = Indirectlinks[m1]; + tmpMatrix(tmpI, tmpJ) = 1; + tmpMatrix(tmpJ, tmpI) = 1; + } + } + Indirectlinks.clear(); + } + + for(int u = 0; u < Nrow; u++) + { + for(int v = 0; v < Ncol; v++) + { + if( tmpMatrix(u, v) > 0.1) + CurrentLinks ++; + } + } + + // cout<<"Link Matrix Symmetrize Loop, PreviousFlag = "<<PreviousLinks<<", CurrentFlag = "<<CurrentLinks<<" of D = "<<Nrow<<" Matrix" <<endl; + + } + } + + return tmpMatrix; +} + +branchcoll ArborBranchMerge(branchcoll inputbranches, TMatrixF ConnectorMatrix) //ABM +{ + branchcoll outputbranches; + + int NinputBranch = inputbranches.size(); + int Nrow = ConnectorMatrix.GetNrows(); + int Ncol = ConnectorMatrix.GetNcols(); + int FlagBranchTouch[Nrow]; + int tmpCellID = 0; + + if(Ncol != NinputBranch || Nrow != Ncol || Nrow != NinputBranch) + { + cout<<"Size of Connector Matrix and inputClusterColl is not match"<<endl; + } + + for(int i0 = 0; i0 < Nrow; i0++) + { + FlagBranchTouch[i0] = 0; + } + + for(int i1 = 0; i1 < Nrow; i1++) + { + if(FlagBranchTouch[i1] == 0) + { + branch Mergebranch_A = inputbranches[i1]; + branch a_MergedBranch = Mergebranch_A; + FlagBranchTouch[i1] = 1; + + for(int j1 = i1 + 1; j1 < Nrow; j1++) + { + if(FlagBranchTouch[j1] == 0 && ConnectorMatrix(i1, j1) > 0.1 ) + { + branch Mergebranch_B = inputbranches[j1]; + FlagBranchTouch[j1] = 1; + for(unsigned int k1 = 0; k1 < Mergebranch_B.size(); k1 ++) + { + tmpCellID = Mergebranch_B[k1]; + if(find(a_MergedBranch.begin(), a_MergedBranch.end(), tmpCellID) == a_MergedBranch.end()) + { + a_MergedBranch.push_back(tmpCellID); + } + } + } + } + outputbranches.push_back(a_MergedBranch); + } + } + + return outputbranches; +} + diff --git a/Reconstruction/PFA/Arbor/src/ArborTool.h b/Reconstruction/PFA/Arbor/src/ArborTool.h new file mode 100644 index 0000000000000000000000000000000000000000..d8d4a48293a164457e489fdc6c5b15ccb212c4d8 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborTool.h @@ -0,0 +1,64 @@ +#ifndef ARBORTOOL_H_ +#define ARBORTOOL_H_ + +#include "TVector3.h" +#include "TMatrixF.h" +#include <iostream> +#include <vector> +#include <string> + + +class QuickUnion{ + std::vector<unsigned> _id; + std::vector<unsigned> _size; + int _count; + + public: + QuickUnion(const unsigned NBranches) { + _count = NBranches; + _id.resize(NBranches); + _size.resize(NBranches); + for( unsigned i = 0; i < NBranches; ++i ) { + _id[i] = i; + _size[i] = 1; + } + } + + int count() const { return _count; } + + unsigned find(unsigned p) { + while( p != _id[p] ) { + _id[p] = _id[_id[p]]; + p = _id[p]; + } + return p; + } + + bool connected(unsigned p, unsigned q) { return find(p) == find(q); } + + void unite(unsigned p, unsigned q) { + unsigned rootP = find(p); + unsigned rootQ = find(q); + _id[p] = q; + + if(_size[rootP] < _size[rootQ] ) { + _id[rootP] = rootQ; _size[rootQ] += _size[rootP]; + } else { + _id[rootQ] = rootP; _size[rootP] += _size[rootQ]; + } + --_count; + } + }; + +typedef std::vector< std::vector<int> > branchcoll; +typedef std::vector<int> branch; +typedef std::vector< std::pair<int, int> > linkcoll; + + +TMatrixF MatrixSummarize( TMatrixF inputMatrix ); //ArborCoreNeed + +std::vector<int>SortMeasure( std::vector<float> Measure, int ControlOrderFlag ); //ArborCoreNeed + +branchcoll ArborBranchMerge(branchcoll inputbranches, TMatrixF inputMatrix); //ArborCoreNeed + +#endif // diff --git a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc new file mode 100644 index 0000000000000000000000000000000000000000..b1e2c7363b093399d4c64fa7477065315149a261 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc @@ -0,0 +1,2062 @@ +#include "ArborToolLCIO.hh" +#include "ArborTool.h" +#include "Arbor.h" +#include "DetectorPos.hh" +#include "HelixClassD.hh" +#include <TMath.h> +#include <values.h> +#include <cmath> +#include <stdexcept> +#include <sstream> + +#include "TVector3.h" +#include <string> +#include <iostream> +#include <fstream> +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/Cluster.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/IDDescriptor.h" +#include "DD4hep/Plugins.h" + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +using namespace std; +/* +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); + LCFlagImpl flag; + flag.setBit(LCIO::CHBIT_LONG); + currbranchcoll->setFlag(flag.getFlag()); + + int NBranch = BranchOrder.size(); + int BranchSize = 0; + float currBranchEnergy = 0; + TVector3 SeedPos, currPos; + float MinMag = 1E9; + float currMag = 0; + float ECALTotalEn = 0; + float HCALTotalEn = 0; + + for(int i0 = 0; i0 < NBranch; i0++) + { + ClusterImpl* a_branch = new ClusterImpl(); + std::vector<int> currbranchorder = BranchOrder[i0]; + BranchSize = currbranchorder.size(); + currBranchEnergy = 0; + ECALTotalEn = 0; + HCALTotalEn = 0; + // CalorimeterHit *Seedhit = Hits[currbranchorder[BranchSize - 1]]; + MinMag = 1E9; + + for(int j = 0; j < BranchSize; j++) + { + CalorimeterHit * a_hit = Hits[currbranchorder[j]]; + currPos = a_hit->getPosition(); + // currMag = DisSeedSurface(currPos); + currMag = currPos.Mag(); + + if( currMag < MinMag ) + { + MinMag = currMag; + SeedPos = currPos; + } + + if(fabs(a_hit->getEnergy() - DHCALCalibrationConstant) < 1.0E-6 ) + { + HCALTotalEn += DHCALCalibrationConstant; + } + else + { + ECALTotalEn += a_hit->getEnergy(); + } + + currBranchEnergy += a_hit->getEnergy(); + + a_branch->addHit(a_hit, (float)1.0); + } + + a_branch->setEnergy(currBranchEnergy); + float ArraySeedPos[3] = { float(SeedPos.X()), float(SeedPos.Y()), float(SeedPos.Z()) }; + a_branch->setPosition( ArraySeedPos ); + + a_branch->subdetectorEnergies().resize(6) ; + a_branch->subdetectorEnergies()[0] = ECALTotalEn ; + a_branch->subdetectorEnergies()[1] = HCALTotalEn ; + a_branch->subdetectorEnergies()[2] = 0 ; + a_branch->subdetectorEnergies()[3] = 0 ; + a_branch->subdetectorEnergies()[4] = 0 ; + a_branch->subdetectorEnergies()[5] = 0 ; + + currbranchcoll -> addElement(a_branch); + } + + evtPP->addCollection(currbranchcoll, Name); +} +*/ + + +ArborToolLCIO::ArborToolLCIO(const std::string& name,ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) +{ + m_geosvc=service<IGeomSvc>("GeomSvc"); +} + + +ArborToolLCIO::~ArborToolLCIO() +{ +} +void ArborToolLCIO::ClusterBuilding( DataHandle<edm4hep::ClusterCollection>& _currbranchcoll, std::vector<edm4hep::ConstCalorimeterHit> Hits, branchcoll BranchOrder, int DHCALFlag ) +{ + //DataHandle<edm4hep::ClusterCollection> _currbranchcoll {"Name",Gaudi::DataHandle::Writer, this}; + //DataHandle<edm4hep::ClusterCollection> _currbranchcoll=new ClusterType(Name, Gaudi::DataHandle::Writer, this); + edm4hep::ClusterCollection* currbranchcoll = _currbranchcoll.createAndPut(); + + int NBranch = BranchOrder.size(); + int BranchSize = 0; + float currBranchEnergy = 0; + TVector3 SeedPos, currPos; + + float MinMag = 1E9; + float currMag = 0; + float ECALTotalEn = 0; + float HCALTotalEn = 0; + + for(int i0 = 0; i0 < NBranch; i0++) + { + auto a_branch=currbranchcoll->create(); + std::vector<int> currbranchorder = BranchOrder[i0]; + BranchSize = currbranchorder.size(); + currBranchEnergy = 0; + ECALTotalEn = 0; + HCALTotalEn = 0; + MinMag = 1E9; + + for(int j = 0; j < BranchSize; j++) + { + auto a_hit= Hits[currbranchorder[j]]; + currPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z);//a_hit.getPosition(); + currMag = currPos.Mag(); + + if( currMag < MinMag ) + { + MinMag = currMag; + SeedPos = currPos; + } + + if(fabs(a_hit.getEnergy() - DHCALCalibrationConstant) < 1.0E-6 ) + { + HCALTotalEn += DHCALCalibrationConstant; + } + else + { + ECALTotalEn += a_hit.getEnergy(); + } + + currBranchEnergy += a_hit.getEnergy(); + + a_branch.addToHits(a_hit); + } + + a_branch.setEnergy(currBranchEnergy); + float ArraySeedPos[3] = { float(SeedPos.X()), float(SeedPos.Y()), float(SeedPos.Z()) }; + a_branch.setPosition( ArraySeedPos ); + + //a_branch.getSubdetectorEnergies().resize(6) ; + //a_branch.getSubdetectorEnergies(0) = ECALTotalEn ; + //a_branch.getSubdetectorEnergies(1) = HCALTotalEn ; + //a_branch.getSubdetectorEnergies(2) = 0 ; + //a_branch.getSubdetectorEnergies(3) = 0 ; + //a_branch.getSubdetectorEnergies(4) = 0 ; + //a_branch.getSubdetectorEnergies(5) = 0 ; + + } + +} + + + +int ArborToolLCIO::NHScaleV2( std::vector<edm4hep::ConstCalorimeterHit> clu0, int RatioX, int RatioY, int RatioZ ) +{ + + int ReScaledNH = 0; + int NumHit = clu0.size(); + int tmpI = 0; + int tmpJ = 0; + int tmpK = 0; + float tmpEn = 0; + int NewCellID0 = 0; + + m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); + if(!m_decoder) m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + + + std::map <double, float> testIDtoEnergy; + + for(int i = 0; i < NumHit; i++) + { + auto hit = clu0[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(); + + NewCellID0 = (tmpK<<24) + (tmpJ<<12) + tmpI; + + if(testIDtoEnergy.find(NewCellID0) == testIDtoEnergy.end() ) + { + testIDtoEnergy[NewCellID0] = tmpEn; + } + else + { + testIDtoEnergy[NewCellID0] += tmpEn; + } + } + + ReScaledNH = testIDtoEnergy.size(); + + return ReScaledNH; +} + +float ArborToolLCIO::FDV2( std::vector<edm4hep::ConstCalorimeterHit> clu) +{ + float FractalDim = 0; + int NReSizeHit[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + int Scale[10] = {2, 3, 4, 5, 6, 7, 8, 9, 10, 20}; + int OriNHit = clu.size(); + + for(int j = 0; j < 10; j++) + { + NReSizeHit[j] = NHScaleV2(clu, Scale[j], Scale[j], 1); + FractalDim += 0.1 * TMath::Log(float(OriNHit)/NReSizeHit[j])/TMath::Log(float(Scale[j])); + } + + if(clu.size() == 0) + FractalDim = -1; + + return FractalDim; +} + + +int ArborToolLCIO::NHScaleV3( edm4hep::ConstCluster clu0, int RatioX, int RatioY, int RatioZ ) +{ + + int ReScaledNH = 0; + int NumHit = clu0.hits_size(); + int tmpI = 0; + int tmpJ = 0; + int tmpK = 0; + float tmpEn = 0; + int NewCellID0 = 0; + int NewCellID1 = 0; + //m_geosvc=service<IGeomSvc>("GeomSvc"); + + m_decoder = m_geosvc->getDecoder("EcalBarrelCollection"); + if(!m_decoder) m_decoder = m_geosvc->getDecoder("EcalEndcapsCollection"); + + 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(); + + NewCellID0 = (tmpK<<24) + (tmpJ<<12) + tmpI; + + testlongID = NewCellID1*1073741824 + NewCellID0; + if(testIDtoEnergy.find(testlongID) == testIDtoEnergy.end() ) + { + testIDtoEnergy[testlongID] = tmpEn; + } + else + { + testIDtoEnergy[testlongID] += tmpEn; + } + } + + ReScaledNH = testIDtoEnergy.size(); + + return ReScaledNH; + +} + +float ArborToolLCIO::FDV3( edm4hep::ConstCluster clu ){ + + float FractalDim = -1; + int NReSizeHit[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + int Scale[5] = {2, 3, 4, 5, 6}; + int OriNHit = clu.hits_size(); + if(OriNHit > 0) + { + FractalDim = 0.0; + for(int j = 0; j < 5; j++) + { + NReSizeHit[j] = NHScaleV3(clu, Scale[j], Scale[j], 1); + FractalDim += 0.2 * TMath::Log(float(OriNHit)/NReSizeHit[j])/TMath::Log(float(Scale[j])); + } + } + return FractalDim; +} + + +float ArborToolLCIO::BushDis( edm4hep::ConstCluster clu1, edm4hep::ConstCluster clu2) +{ + float DisBetweenBush = 1.0E10; + + int cluSize1 = clu1.hits_size(); + int cluSize2 = clu2.hits_size(); + + TVector3 HitPos1, HitPos2; + TVector3 PosDiff; + // TVector3 XXXPos; + + for(int i = 0; i < cluSize1; i++) + { + HitPos1 = TVector3((clu1.getHits(i)).getPosition().x,(clu1.getHits(i)).getPosition().y,(clu1.getHits(i)).getPosition().z); + for(int j = 0; j<cluSize2; j++) + { + HitPos2 = TVector3((clu2.getHits(j)).getPosition().x,(clu2.getHits(j)).getPosition().y,(clu2.getHits(j)).getPosition().z); + PosDiff = HitPos1 - HitPos2; + + if(PosDiff.Mag() < DisBetweenBush ) + { + DisBetweenBush = PosDiff.Mag(); + } + } + } + + + return DisBetweenBush; +} + + +float ArborToolLCIO::DisPointToBush(TVector3 Pos1, edm4hep::ConstCluster clu1) +{ + float Dis = 1.0E9; + float HitDis = 1.0E8; + int clusize = clu1.hits_size(); + + 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) + { + Dis = HitDis; + } + } + + return Dis; +} + + +TVector3 ArborToolLCIO::ClusterCoG(edm4hep::ConstCluster inputCluster) +{ + TVector3 CenterOfGravity; + + int inputClusterSize = inputCluster.hits_size(); + + TVector3 tmphitPos; + float tmphitEnergy; + float sumhitEnergy = 0; + + for(int i = 0; i < inputClusterSize; i++) + { + auto tmpHit = inputCluster.getHits(i); + tmphitPos = TVector3(tmpHit.getPosition().x,tmpHit.getPosition().y,tmpHit.getPosition().z); + tmphitEnergy = tmpHit.getEnergy(); + + CenterOfGravity += tmphitPos*tmphitEnergy; + sumhitEnergy += tmphitEnergy; + } + + CenterOfGravity = 1.0/sumhitEnergy * CenterOfGravity; + + return CenterOfGravity; +} + + +edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecColl( std::vector<edm4hep::ConstCluster> inputClusters, DataHandle<edm4hep::ClusterCollection>& m_clucol ) +{ + + 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; + + for(int i0 = 0; i0 < NClu; i0++) + { + auto a_clu = inputClusters[i0]; + CluEn.push_back(a_clu.getEnergy()); + } + CluIndex = SortMeasure(CluEn, 1); + + for(int i1 = 0; i1 < NClu; i1++) + { + auto branchtmp=vec_coll_Clusters->create(); + auto a_clu = inputClusters[CluIndex[i1]]; + + CurrBranchSize = a_clu.hits_size(); + + for(int j1 = 0; j1 < CurrBranchSize; j1 ++) + { + auto tmpHit = a_clu.getHits(j1); + branchtmp.addToHits(tmpHit); + } + + branchtmp.setPosition( a_clu.getPosition() ); + branchtmp.setEnergy(a_clu.getEnergy() ); + SeedPos = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + branchtmp.setITheta( SeedPos.Theta() ); //To be replaced, those worse than 1st order appro + branchtmp.setPhi( SeedPos.Phi() ); + + } + + return vec_coll_Clusters; +} + +std::vector<edm4hep::ConstCluster> ArborToolLCIO::CollClusterVec(const edm4hep::ClusterCollection * input_coll ) +{ + std::vector<edm4hep::ConstCluster> outputClusterVec; + + + outputClusterVec.clear(); + + for(int i = 0; i < input_coll->size(); i++) //We can have some sort here - according to depth/energy... + { + auto a_clu = (*input_coll)[i]; + outputClusterVec.push_back(a_clu); + } + + return outputClusterVec; +} + + +void ArborToolLCIO::NaiveCluConst(edm4hep::ConstCluster a0_clu,edm4hep::Cluster b0_clu) +{ + b0_clu.setPosition(a0_clu.getPosition()); + b0_clu.setEnergy(a0_clu.getEnergy()); + int NCaloHit = a0_clu.hits_size(); + float HitEn = 0; + float SubDEn[6] = {0, 0, 0, 0, 0, 0}; + + for(int t0 = 0; t0 < NCaloHit; t0++) + { + auto a0_hit = a0_clu.getHits(t0); + b0_clu.addToHits(a0_hit); + HitEn = a0_hit.getEnergy(); + if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) + { + SubDEn[1] += HitEn; + } + else + { + SubDEn[0] += HitEn; + } + } + + for(int i = 0; i < 6; i++) + { + b0_clu.addToSubdetectorEnergies(SubDEn[i]); + } + +} + + +edm4hep::Cluster ArborToolLCIO::NaiveCluImpl(edm4hep::ConstCluster a0_clu) +{ + edm4hep::Cluster b0_clu; + b0_clu.setPosition(a0_clu.getPosition()); + b0_clu.setEnergy(a0_clu.getEnergy()); + int NCaloHit = a0_clu.hits_size(); + float HitEn = 0; + float SubDEn[6] = {0, 0, 0, 0, 0, 0}; + + for(int t0 = 0; t0 < NCaloHit; t0++) + { + auto a0_hit = a0_clu.getHits(t0); + b0_clu.addToHits(a0_hit); + HitEn = a0_hit.getEnergy(); + if(fabs(HitEn - DHCALCalibrationConstant) < 1.0E-6) + { + SubDEn[1] += HitEn; + } + else + { + SubDEn[0] += HitEn; + } + } + + for(int i = 0; i < 6; i++) + { + b0_clu.addToSubdetectorEnergies(SubDEn[i]); + } + + return b0_clu; +} + +std::vector<edm4hep::ConstCalorimeterHit> ArborToolLCIO::CollHitVec(const edm4hep::CalorimeterHitCollection * input_coll, float EnergyThreshold) +{ + std::vector<edm4hep::ConstCalorimeterHit> outputHitVec; + + outputHitVec.clear(); + + for(int i = 0; i < input_coll->size(); i++) //We can have some sort here - according to depth/energy... + { + auto a_hit = (*input_coll)[i]; + if(a_hit.getEnergy() > EnergyThreshold) + { + outputHitVec.push_back(a_hit); + } + } + + return outputHitVec; +} + + +std::vector<edm4hep::Cluster> ArborToolLCIO::ClusterHitAbsorbtion( std::vector<edm4hep::ConstCluster> MainClusters, std::vector<edm4hep::ConstCalorimeterHit> IsoHits, float DisThreshold ) // Projective Distance + Hit Depth correlation; +{ + std::vector<edm4hep::Cluster> outputClusterVec; + + int N_Core = MainClusters.size(); + int N_Hit = IsoHits.size(); + TVector3 HitPos, MBSeedPos; + 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; + + 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); + MinHitCoreDis = 1.0E10; + + for(int j0 = 0; j0 < N_Core; j0++) + { + auto a_core = MainClusters[j0]; + currHitCoreDis = DisPointToBush(HitPos, a_core); + if(currHitCoreDis < MinHitCoreDis) + { + MinHitCoreDis = currHitCoreDis; + MinDisIndex = j0; + } + } + if(MinHitCoreDis < DisThreshold) + { + a_frag_core_link.first = i0; + a_frag_core_link.second = MinDisIndex; + Frag_Core_Links.push_back(a_frag_core_link); + } + } + + int N_frag_core_links = Frag_Core_Links.size(); + std::vector<edm4hep::ConstCalorimeterHit> tomerge_hits; + float ClusterEn = 0; + + for(int i2 = 0; i2 < N_Core; i2 ++) + { + auto a_core = MainClusters[i2]; + tomerge_hits.clear(); + + for(int j4 = 0; j4 < N_frag_core_links; j4 ++) + { + a_frag_core_link = Frag_Core_Links[j4]; + if(a_frag_core_link.second == i2) + { + auto a_frag = IsoHits[a_frag_core_link.first]; + tomerge_hits.push_back(a_frag); + } + } + edm4hep::Cluster a_mergedfrag_core; + ClusterEn = 0; + + for(unsigned int j2 = 0; j2 < a_core.hits_size(); j2++) + { + auto b_hit = a_core.getHits(j2); + a_mergedfrag_core.addToHits(b_hit); + ClusterEn += b_hit.getEnergy(); + } + + for(unsigned int j3 = 0; j3 < tomerge_hits.size(); j3++) + { + auto c_hit = tomerge_hits[j3]; + a_mergedfrag_core.addToHits(c_hit); + ClusterEn += c_hit.getEnergy(); + } + + a_mergedfrag_core.setPosition( a_core.getPosition() ); + a_mergedfrag_core.setEnergy(ClusterEn); + MBSeedPos = TVector3(a_core.getPosition().x,a_core.getPosition().y,a_core.getPosition().z); + a_mergedfrag_core.setITheta( MBSeedPos.Theta() ); //To be replaced, those worse than 1st order appro + a_mergedfrag_core.setPhi( MBSeedPos.Phi() ); + + outputClusterVec.push_back(a_mergedfrag_core); + } + + return outputClusterVec; +} + + +void ArborToolLCIO::NaiveMergeCluConst(std::vector<edm4hep::ConstCluster> inputCluVec,edm4hep::Cluster MergedClu) +{ + + int NClu = inputCluVec.size(); + 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... + + for(int i = 0; i < NClu; i++) + { + auto a_Clu = inputCluVec[i]; + CurrSeedPos = TVector3(a_Clu.getPosition().x,a_Clu.getPosition().y,a_Clu.getPosition().z); + MergedCluEnergy += a_Clu.getEnergy(); + + if(CurrSeedPos.Mag() < SeedDis) + { + SeedPos = CurrSeedPos; + SeedDis = CurrSeedPos.Mag(); + } + + int CurrCluSize = a_Clu.hits_size(); + + for(int j = 0; j < CurrCluSize; j++) + { + auto a_hit = a_Clu.getHits(j); + 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. + { + SubDEn[1] += HitEn; + } + else + { + SubDEn[0] += HitEn; + } + + if(CurrHitPos.Mag() > MaxDis) + { + MaxDis = CurrHitPos.Mag(); + CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + } + } + } + CluRefDir = (CluEndPos - SeedPos); + + float ClusterSeedPos[3] = {float(SeedPos.X()), float(SeedPos.Y()), float(SeedPos.Z())}; + MergedClu.setPosition(ClusterSeedPos); + MergedClu.setEnergy(MergedCluEnergy); + MergedClu.setITheta( CluRefDir.Theta() ); + MergedClu.setPhi( CluRefDir.Phi() ); + + //MergedClu.subdetectorEnergies().resize(6) ; + for(int i = 0; i < 6; i++) + { + // MergedClu.subdetectorEnergies()[i] = SubDEn[i]; + MergedClu.addToSubdetectorEnergies(SubDEn[i]); + } + +} +edm4hep::Cluster ArborToolLCIO::NaiveMergeClu(std::vector<edm4hep::ConstCluster> inputCluVec) +{ + edm4hep::Cluster MergedClu; + + int NClu = inputCluVec.size(); + 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... + + for(int i = 0; i < NClu; i++) + { + auto a_Clu = inputCluVec[i]; + CurrSeedPos = TVector3(a_Clu.getPosition().x,a_Clu.getPosition().y,a_Clu.getPosition().z); + MergedCluEnergy += a_Clu.getEnergy(); + + if(CurrSeedPos.Mag() < SeedDis) + { + SeedPos = CurrSeedPos; + SeedDis = CurrSeedPos.Mag(); + } + + int CurrCluSize = a_Clu.hits_size(); + + for(int j = 0; j < CurrCluSize; j++) + { + auto a_hit = a_Clu.getHits(j); + 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. + { + SubDEn[1] += HitEn; + } + else + { + SubDEn[0] += HitEn; + } + + if(CurrHitPos.Mag() > MaxDis) + { + MaxDis = CurrHitPos.Mag(); + CluEndPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + } + } + } + CluRefDir = (CluEndPos - SeedPos); + + float ClusterSeedPos[3] = {float(SeedPos.X()), float(SeedPos.Y()), float(SeedPos.Z())}; + MergedClu.setPosition(ClusterSeedPos); + MergedClu.setEnergy(MergedCluEnergy); + MergedClu.setITheta( CluRefDir.Theta() ); + MergedClu.setPhi( CluRefDir.Phi() ); + + //MergedClu.subdetectorEnergies().resize(6) ; + for(int i = 0; i < 6; i++) + { + // MergedClu.subdetectorEnergies()[i] = SubDEn[i]; + MergedClu.addToSubdetectorEnergies(SubDEn[i]); + } + + return MergedClu; +} + + + +std::vector<edm4hep::ConstCluster> ArborToolLCIO::ClusterAbsorbtion( std::vector<edm4hep::ConstCluster> MainClusters, std::vector<edm4hep::ConstCluster> FragClusters, float DisThreshold, float DepthSlope ) //ProjectiveDis +{ + std::vector<edm4hep::ConstCluster> outputClusterVec; + + int N_Core = MainClusters.size(); + int N_frag = FragClusters.size(); + + //tag minimal distance + + float MinFragCoreDis = 1.0E10; + float CurrFragCoreDis = 0; + 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; + + for(int i0 = 0; i0 < N_frag; i0 ++) + { + auto a_frag = FragClusters[i0]; + fragPos = TVector3(a_frag.getPosition().x,a_frag.getPosition().y,a_frag.getPosition().z); + MinFragCoreDis = 1.0E10; + for(int j0 = 0; j0 < N_Core; j0++) + { + auto a_core = MainClusters[j0]; + CurrFragCoreDis = BushDis(a_frag, a_core); + if(CurrFragCoreDis < MinFragCoreDis) + { + MinFragCoreDis = CurrFragCoreDis; + MinDisIndex = j0; + } + } + + if( MinFragCoreDis < DisThreshold + DepthSlope*DisSeedSurface(fragPos)) + { + a_frag_core_link.first = i0; + a_frag_core_link.second = MinDisIndex; + Frag_Core_Links.push_back(a_frag_core_link); + } + } + + int N_frag_core_links = Frag_Core_Links.size(); + std::vector<edm4hep::ConstCluster> tomerge_clu; + + for(int i4 = 0; i4 < N_Core; i4 ++) + { + auto a_core = MainClusters[i4]; + tomerge_clu.clear(); + tomerge_clu.push_back(a_core); + + for(int j4 = 0; j4 < N_frag_core_links; j4 ++) + { + a_frag_core_link = Frag_Core_Links[j4]; + if(a_frag_core_link.second == i4) + { + auto a_frag = FragClusters[a_frag_core_link.first]; + TouchedFrag[a_frag_core_link.first] = 1; + tomerge_clu.push_back(a_frag); + } + } + auto a_mergedfrag_core = NaiveMergeClu(tomerge_clu); + edm4hep::ConstCluster a_mergedfrag_coreCon=a_mergedfrag_core; + outputClusterVec.push_back(a_mergedfrag_core); + } + + for(int i1 = 0; i1 < N_frag; i1++) + { + if(TouchedFrag.find(i1) == TouchedFrag.end()) + { + auto a_frag = FragClusters[i1]; + //auto a_isofrag = NaiveCluImpl(a_frag);; + outputClusterVec.push_back(a_frag); + } + } + + return outputClusterVec; +} + + +edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecMerge( std::vector<edm4hep::ConstCluster> inputClusters, TMatrixF ConnectorMatrix, DataHandle<edm4hep::ClusterCollection>& clucol ) +{ + edm4hep::ClusterCollection* mergedbranches = clucol.createAndPut(); + + + int NinputClu = inputClusters.size(); + int Nrow = ConnectorMatrix.GetNrows(); + int Ncol = ConnectorMatrix.GetNcols(); + + if(Ncol != NinputClu || Nrow != Ncol || Nrow != NinputClu) + { + cout<<"Size of Connector Matrix and inputClusterColl is not match"<<endl; + } + + vector<edm4hep::ConstCluster> branchToMerge; + edm4hep::ConstCluster Mergebranch_A; + edm4hep::ConstCluster Mergebranch_B; + edm4hep::ConstCluster tmpMergebranch; + edm4hep::ConstCluster Mainbranch (0); + + TVector3 tmpClusterSeedPos, MBSeedPos; + + int CurrBranchSize = 0; + float SeedPosMin = 1.0E10; + float BranchEnergy = 0; + + int FlagBranchTouch[Nrow]; + + for(int i0 = 0; i0 < Nrow; i0++) + { + FlagBranchTouch[i0] = 0; + } + + for(int ibran = 0; ibran < Nrow ; ibran++) + { + if(FlagBranchTouch[ibran] == 0) + { + Mergebranch_A = inputClusters[ibran]; + branchToMerge.push_back(Mergebranch_A); + FlagBranchTouch[ibran] = 1; + BranchEnergy = 0; + auto branchtmp = mergedbranches->create(); + + for(int jbran = ibran + 1; jbran < Nrow; jbran++) + { + if(FlagBranchTouch[jbran] == 0) + { + Mergebranch_B = inputClusters[jbran]; + if( ConnectorMatrix(ibran, jbran) > 0.1 ) + { + branchToMerge.push_back(Mergebranch_B); + FlagBranchTouch[jbran] = 1; + } + } + } + + SeedPosMin = 1.0E10; + + for(unsigned int i1 = 0; i1 < branchToMerge.size(); i1++) + { + tmpMergebranch = branchToMerge[i1]; + tmpClusterSeedPos = TVector3(tmpMergebranch.getPosition().x,tmpMergebranch.getPosition().y,tmpMergebranch.getPosition().z); + if( tmpClusterSeedPos.Mag() < SeedPosMin) + { + + Mainbranch = tmpMergebranch; + SeedPosMin = tmpClusterSeedPos.Mag(); + } + + CurrBranchSize = tmpMergebranch.hits_size(); + BranchEnergy += tmpMergebranch.getEnergy(); + + for(int j1 = 0; j1 < CurrBranchSize; j1 ++) + { + auto tmpHit = tmpMergebranch.getHits(j1); + branchtmp.addToHits(tmpHit); + } + + } + branchtmp.setPosition( Mainbranch.getPosition() ); + branchtmp.setEnergy(BranchEnergy); + MBSeedPos = TVector3(Mainbranch.getPosition().x,Mainbranch.getPosition().y,Mainbranch.getPosition().z); + branchtmp.setITheta( MBSeedPos.Theta() ); //To be replaced, those worse than 1st order appro + branchtmp.setPhi( MBSeedPos.Phi() ); + + branchToMerge.clear(); + } + } + + return mergedbranches; + +} +int ArborToolLCIO::JointsBetweenBush(edm4hep::ConstCluster a_Clu, edm4hep::ConstCluster b_Clu, float CellSize) +{ + int NJoint = 0; + int a_CluSize = a_Clu.hits_size(); + int b_CluSize = b_Clu.hits_size(); + 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); + + for(int i = 0; i < a_CluSize; i++) + { + auto ahit = a_Clu.getHits(i); + aHitPos = TVector3(ahit.getPosition().x,ahit.getPosition().y,ahit.getPosition().z); + for(int j = 0; j < b_CluSize; j++) + { + auto bhit = b_Clu.getHits(j); + bHitPos = TVector3(bhit.getPosition().x,bhit.getPosition().y,bhit.getPosition().z); + 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) + NJoint++; //Change to NJoint Hit... + } + } + } + + return NJoint; +} + + +float* ArborToolLCIO::SimpleDisTrackClu( edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu) +{ + float* Distance = new float[3]; + Distance[0] = 1.0E9; + Distance[1] = 1.0E9; + Distance[2] = 1.0E9; + float minDis = 1.0E9; + float BushDist[3] = {0, 0 ,0}; + HelixClassD * a_Helix = new HelixClassD(); + //float refPoint[3] = a_Helix->getReferencePoint(); + a_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); + int NCaloHits = a_clu.hits_size(); + + for(int i1 = 0; i1 < NCaloHits; i1++) + { + auto a_hit = a_clu.getHits(i1); + + float hitpos[3]={a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z}; + float BushTime = a_Helix->getDistanceToPoint(hitpos, BushDist); + if(BushTime > 0 && BushDist[2] < minDis ) + { + minDis = BushDist[2]; + Distance[0] = BushDist[0]; + Distance[1] = BushDist[1]; + Distance[2] = BushDist[2]; + } + } + delete a_Helix; + + return Distance; +} +float ArborToolLCIO::SimpleBushTimeTrackClu(edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu) +{ + float Distance = 1.0E9; + float Time = 0; + float BushDist[3] = {0, 0 ,0}; + HelixClassD * a_Helix = new HelixClassD(); + a_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); + int NCaloHits = a_clu.hits_size(); + + for(int i1 = 0; i1 < NCaloHits; i1++) + { + auto a_hit = a_clu.getHits(i1); + + float hitpos[3]={a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z}; + float BushTime = a_Helix->getDistanceToPoint(hitpos, BushDist); + if(BushTime > 0 && BushDist[2] < Distance ) + { + Time = BushTime; + Distance = BushDist[2]; + } + } + delete a_Helix; + return Time; +} + +int ArborToolLCIO::SimpleBushNC(edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu) +{ + float Distance = 1.0E9; + //float Time = 0; + int NC = 0; + float BushDist[3] = {0, 0 ,0}; + HelixClassD * a_Helix = new HelixClassD(); + a_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); + int NCaloHits = a_clu.hits_size(); + + for(int i1 = 0; i1 < NCaloHits; i1++) + { + auto a_hit = a_clu.getHits(i1); + + float hitpos[3]={a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z}; + float BushTime = a_Helix->getDistanceToPoint(hitpos, BushDist); + int nCircles = a_Helix->getNCircle(hitpos); + if(BushTime > 0 && BushDist[2] < Distance ) + { + //Time = BushTime; + NC = nCircles; + Distance = BushDist[2]; + } + } + + delete a_Helix; + return NC; +} + +int ArborToolLCIO::ClusterFlag(edm4hep::ConstCluster a_tree, edm4hep::ConstTrack a_trk) +{ + // give each charged core cluster a flag + // Fragmentation: 999 + // MIP: matched 130 + // non-matched 131 + // EM: matched 110 + // non-matched 111 + // HAD: matched 211 + // non-matched 212 + + int CluIDFlag = 999; + int ClusterID = 211; + int EcalNHit, HcalNHit, NLEcal, NLHcal, NH[16]; + float avEnDisHtoL; + float EcalEn, HcalEn, EClu, cluDepth, maxDepth, minDepth, MaxDisHel, MinDisHel, FD_all, FD_ECAL, FD_HCAL; + float crdis, EEClu_L10, EEClu_R, EEClu_r, EEClu_p, rms_Ecal, rms_Hcal, rms_Ecal2, rms_Hcal2, av_NHE, av_NHH; + int AL_Ecal, AL_Hcal; + float FD_ECALF10; + float dEdx = 0; + + const float mass = 0.139; //Pion Mass + + 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]); + } + delete TrkInit_Helix; + + TrackEn = sqrt(TrackEn); + int nSubTrk = a_trk.getTracks().size(); + + //int NHit = a_tree->getHits().size(); + if(1 == 1) //if ( (NHit > 4 && TrackEn > 1) || TrackEn <= 1 ) + { + + TVector3 CluPos; + CluPos = TVector3(a_tree.getPosition().x,a_tree.getPosition().y,a_tree.getPosition().z); + TVector3 IntDir = ClusterCoG(a_tree)-CluPos; + EClu = a_tree.getEnergy(); + EcalNHit = 0; + HcalNHit = 0; + EcalEn = 0; + HcalEn = 0; + float currDepth = 0; + maxDepth = -100; + minDepth = 1E6; + MaxDisHel = -1; //maximal distance from Track to Helix + MinDisHel = 1E10; + + EEClu_R = 0; + EEClu_r = 0; + EEClu_p = 0; + EEClu_L10 = 0; + + + std::vector<edm4hep::ConstCalorimeterHit> Ecalhits; + std::vector<edm4hep::ConstCalorimeterHit> Hcalhits; + std::vector<edm4hep::ConstCalorimeterHit> allhits; + std::vector<edm4hep::ConstCalorimeterHit> EH_1; + std::vector<edm4hep::ConstCalorimeterHit> EH_2; + std::vector<edm4hep::ConstCalorimeterHit> EH_3; + std::vector<edm4hep::ConstCalorimeterHit> EH_4; + std::vector<edm4hep::ConstCalorimeterHit> EH_5; + std::vector<edm4hep::ConstCalorimeterHit> EH_6; + std::vector<edm4hep::ConstCalorimeterHit> HH_1; + std::vector<edm4hep::ConstCalorimeterHit> HH_2; + std::vector<edm4hep::ConstCalorimeterHit> HH_3; + std::vector<edm4hep::ConstCalorimeterHit> HH_4; + std::vector<edm4hep::ConstCalorimeterHit> HH_5; + std::vector<edm4hep::ConstCalorimeterHit> HH_6; + std::vector<edm4hep::ConstCalorimeterHit> HH_7; + std::vector<edm4hep::ConstCalorimeterHit> HH_8; + std::vector<edm4hep::ConstCalorimeterHit> HH_9; + std::vector<edm4hep::ConstCalorimeterHit> HH_0; + std::vector<edm4hep::ConstCalorimeterHit> Ecalf10hits; + + + + allhits.clear(); + Ecalhits.clear(); + Hcalhits.clear(); + EH_1.clear(); + EH_2.clear(); + EH_3.clear(); + EH_4.clear(); + EH_5.clear(); + EH_6.clear(); + HH_1.clear(); + HH_2.clear(); + HH_3.clear(); + HH_4.clear(); + HH_5.clear(); + HH_6.clear(); + HH_7.clear(); + HH_8.clear(); + HH_9.clear(); + HH_0.clear(); + Ecalf10hits.clear(); + + + HelixClassD * currHelix = new HelixClassD(); + currHelix->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 BushDist[3] = {0, 0, 0}; + float BushTime = 0; + + std::vector<float> hitTheta; + hitTheta.clear(); + + for(unsigned int j1 = 0; j1 < a_tree.hits_size(); j1++) + { + auto a_hit = a_tree.getHits(j1); + float hitpos[3]={a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z}; + BushTime = currHelix->getDistanceToPoint(hitpos, BushDist); + TVector3 tmpPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + hitTheta.push_back(tmpPos.Theta()); + if(BushTime > 0) + { + if(BushDist[2] > MaxDisHel) + { + MaxDisHel = BushDist[2]; + } + if(BushDist[2] < MinDisHel) + { + MinDisHel = BushDist[2]; + } + } + } + delete currHelix; + + float totTheta = 0; + float avTheta = 0; + float SDTheta; + + for(int t0 = 0; t0 < int(hitTheta.size()); t0++) + { + float tmpTheta = hitTheta[t0]; + totTheta += tmpTheta; + } + + avTheta = totTheta/float(hitTheta.size()); + SDTheta = 0; + + for(int t1 = 0; t1 < int(hitTheta.size()); t1++) + { + float tmpTheta = hitTheta[t1]; + SDTheta += pow(tmpTheta-avTheta,2); + } + SDTheta = sqrt(SDTheta/float(hitTheta.size())); + + TVector3 HitPos; + int currCluNHits = a_tree.hits_size(); + if(currCluNHits == 0) return CluIDFlag; + int index1 = 0, index2 = 0; + + for(int s1 = 0; s1 < currCluNHits; s1++) + { + auto a_hit = a_tree.getHits(s1); + allhits.push_back(a_hit); + auto cellid= a_hit.getCellID(); + int NLayer = m_decoder->get(cellid, "layer"); + + HitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + + currDepth = DisSeedSurface(HitPos); + crdis = (CluPos-HitPos).Mag()*sin((CluPos-HitPos).Angle(IntDir)); + + + if(currDepth > maxDepth) + { + maxDepth = currDepth; + index1 = s1; + } + if(currDepth < minDepth) + { + minDepth = currDepth; + index2 = s1; + } + + if( fabs(a_hit.getEnergy() - DHCALCalibrationConstant) < 1.0E-6 ) //or other fancy judgements...^M + { + HcalNHit++; + HcalEn += a_hit.getEnergy(); + Hcalhits.push_back(a_hit); + if(NLayer < 5) + { + HH_1.push_back(a_hit); + } + else if(NLayer < 10) + { + HH_2.push_back(a_hit); + } + else if(NLayer < 15) + { + HH_3.push_back(a_hit); + } + else if(NLayer < 20) + { + HH_4.push_back(a_hit); + } + else if(NLayer < 25) + { + HH_5.push_back(a_hit); + } + else if(NLayer < 30) + { + HH_6.push_back(a_hit); + } + else if(NLayer < 35) + { + HH_7.push_back(a_hit); + } + else if(NLayer < 40) + { + HH_8.push_back(a_hit); + } + else if(NLayer < 45) + { + HH_9.push_back(a_hit); + } + else + { + HH_0.push_back(a_hit); + } + } + else + { + EcalNHit++; + EcalEn += a_hit.getEnergy(); + Ecalhits.push_back(a_hit); + if(NLayer< 10) Ecalf10hits.push_back(a_hit); + if(crdis < 22) EEClu_R += a_hit.getEnergy(); + if(crdis < 11) EEClu_r += a_hit.getEnergy(); + if(crdis < 6) EEClu_p += a_hit.getEnergy(); + if(NLayer < 5) + { + EH_1.push_back(a_hit); + EEClu_L10 += a_hit.getEnergy(); + } + else if(NLayer < 10) + { + EH_2.push_back(a_hit); + EEClu_L10 += a_hit.getEnergy(); + } + 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); + } + else + { + EH_6.push_back(a_hit); + } + } + } + + auto maxdis_hit = a_tree.getHits(index1); + auto mindis_hit = a_tree.getHits(index2); + TVector3 maxpos = TVector3(maxdis_hit.getPosition().x,maxdis_hit.getPosition().y,maxdis_hit.getPosition().z); + TVector3 minpos = TVector3(mindis_hit.getPosition().x,mindis_hit.getPosition().y,mindis_hit.getPosition().z); + TVector3 GraPos = ClusterCoG(a_tree); + cluDepth = (maxpos-minpos).Mag(); + + float totHitEn = 0; + float totHitEnDis = 0; + float HitEn; + + for(int s2 = 0; s2 < currCluNHits; s2++) + { + auto a_hit2 = a_tree.getHits(s2); + HitPos = TVector3(a_hit2.getPosition().x,a_hit2.getPosition().y,a_hit2.getPosition().z); + HitEn = a_hit2.getEnergy(); + TVector3 par1 = GraPos-minpos; + TVector3 par2 = minpos-HitPos; + TVector3 par3 = par1.Cross(par2); + float disHtoL = par3.Mag()/par1.Mag(); + totHitEn+=HitEn; + totHitEnDis+=HitEn*disHtoL; + } + avEnDisHtoL = totHitEnDis/totHitEn; + FD_all = FDV2(allhits); + FD_ECAL = FDV2(Ecalhits); + FD_HCAL = FDV2(Hcalhits); + FD_ECALF10 = FDV2(Ecalf10hits); + + NLEcal = 0; + NLHcal = 0; + for(int p0 = 0; p0 < 8; p0++) + { + NH[p0] = 0; + } + + NH[0] = EH_1.size(); + NH[1] = EH_2.size(); + NH[2] = EH_3.size(); + NH[3] = EH_4.size(); + NH[4] = EH_5.size(); + NH[5] = EH_6.size(); + NH[6] = HH_1.size(); + NH[7] = HH_2.size(); + NH[8] = HH_3.size(); + NH[9] = HH_4.size(); + NH[10] = HH_5.size(); + NH[11] = HH_6.size(); + NH[12] = HH_7.size(); + NH[13] = HH_8.size(); + NH[14] = HH_9.size(); + NH[15] = HH_0.size(); + + NLEcal = ActiveLayers(Ecalhits); + NLHcal = ActiveLayers(Hcalhits); + cout<<"NLEcal "<<NLEcal<<" "<<Ecalhits.size()<<" "<<endl; + + + float sum_NHE = 0, sum_NHH = 0; + av_NHE = 0; + av_NHH = 0; + AL_Ecal = 0; + AL_Hcal = 0; + + for(int r1 = 0; r1 < 16; r1++) + { + if(r1 < 6 && NH[r1]>0) + { + sum_NHE += NH[r1]; + AL_Ecal++; + } + if(r1 >= 6 && NH[r1]>0) + { + sum_NHH += NH[r1]; + AL_Hcal++; + } + } + if(AL_Ecal > 0) + av_NHE = sum_NHE/AL_Ecal; + if(AL_Hcal > 0) + av_NHH = sum_NHH/AL_Hcal; + + rms_Ecal = 0; + rms_Hcal = 0; + rms_Ecal2 = 0; + rms_Hcal2 = 0; + for(int r0 = 0; r0 < 16; r0++) + { + if(r0 < 6) + { + if(NH[r0] > 0) + { + rms_Ecal+=pow(NH[r0]-av_NHE,2); + rms_Ecal2 += pow(NH[r0],2); + } + } + else + { + if(NH[r0] > 0) + { + rms_Hcal+=pow(NH[r0]-av_NHH,2); + rms_Hcal2 += pow(NH[r0],2); + } + } + } + if(AL_Ecal > 0) + { + rms_Ecal2 = sqrt(rms_Ecal2/AL_Ecal); + rms_Ecal = sqrt(rms_Ecal/AL_Ecal); + } + else + { + rms_Ecal2 = -1; + rms_Ecal = -1; + } + if(AL_Hcal > 0) + { + rms_Hcal2 = sqrt(rms_Hcal2/AL_Ecal); + rms_Hcal = sqrt(rms_Hcal/AL_Ecal); + } + else + { + rms_Hcal2 = -1; + rms_Hcal = -1; + } + + + + + bool cutmu1 = EcalNHit+2.*HcalNHit < 500; + bool cutmu2; + + if(cluDepth < 1100) cutmu2 = 0; + else if(cluDepth < 1400) cutmu2 = FD_all < 0.5/pow(400,2)*pow((cluDepth-1000),2); + else cutmu2 = 1; + + bool cutmu3; + if(TrackEn > 70) cutmu3 = FD_all < (600./avEnDisHtoL + 20)/100.; + else if (TrackEn > 10) cutmu3 = FD_all < (600./avEnDisHtoL -10 + 0.5*TrackEn)/100.; + else if (TrackEn > 7.5) cutmu3 = FD_all < (600./avEnDisHtoL -10)/100.; + else cutmu3 = 1; + bool cutmu3b; + if (TrackEn > 10) cutmu3b = avEnDisHtoL < 25; + else if (TrackEn > 4.5) cutmu3b = avEnDisHtoL < 25 + 10 * (10 - TrackEn); + else cutmu3b = 1; + bool cutmu10en = cutmu1 && cutmu2 && cutmu3 && cutmu3b; + + + bool cutmu4 = FD_HCAL >= 0; + bool cutmu5 = cluDepth > 750 - 20/TrackEn; + bool cutmu6 = cluDepth > 1200; + bool cutmu7 = FD_all < 0.3/sqrt(400.)*sqrt(cluDepth-1200.); + bool cutmu8 = rms_Hcal < 10 && FD_HCAL < -0.25/600.*MaxDisHel+0.25; + bool cutmu9 = FD_all < 0.35/sqrt(400)*sqrt(400-MaxDisHel) && MaxDisHel < 400 && EcalNHit+2.*HcalNHit > 85; + bool cutmu; + + if(TrackEn > 9.5) cutmu = cutmu10en; + else if(TrackEn < 1.5) + { + if(cutmu5) cutmu = 1; + else cutmu = 0; + } + else if(TrackEn < 3.5) + { + if(cutmu4 && cutmu6 && cutmu7) cutmu = 1; + else cutmu = 0; + } + else + { + if(cutmu3b && cutmu4 && cutmu6 && cutmu7 && cutmu8 && cutmu9) cutmu = 1; + else cutmu = 0; + } + + //bool cute1 = EcalEn/(EcalEn+HcalEn) > 0.9; + bool cute2 = (dEdx > 0.17e-6 && dEdx < 0.3e-6)||dEdx==0; + bool cute3 = FD_all > 0.9*sin(cluDepth*1.57/800.) && cluDepth < 800; + bool cute4 = FD_ECALF10 > 0.9*sin((cluDepth-200)*1.57/600.); + bool cute5 = EEClu_r/TrackEn > 0.8 * sin((avEnDisHtoL-15)*1.57/20.) && avEnDisHtoL < 35; + bool cute6 = FD_ECAL > 0.2 * log10(EcalNHit) && log10(EcalNHit) > 1.5; + bool cute7 = FD_all >= 0.6/1.5*(log10(EcalNHit+2.*HcalNHit)-1.2-0.4*TrackEn/100.); + bool cute8 = SDTheta < 0.012/1.5*(log10(EcalNHit+2.*HcalNHit)-1); + bool cute; + if(TrackEn < 1.5) cute = cute2; + else cute = cute2 && cute3 && cute4 && cute5 && cute6 && cute7 && cute8; + + + if(cutmu) ClusterID = 13; + else if (cute) ClusterID = 11; + + if(ClusterID == 13 ) + { + if(NLEcal+NLHcal < 30) CluIDFlag = 131; + else CluIDFlag = 130; + } + + if(ClusterID == 11 ) + { + if(EClu < 0.8*TrackEn) CluIDFlag = 111; + else CluIDFlag = 110; + } + + if(ClusterID == 211 ) + { + if(EClu < 0.8*TrackEn) CluIDFlag = 212; + else CluIDFlag = 211; + } + + } + + return CluIDFlag; + +} + + +int ArborToolLCIO::ActiveLayers( std::vector<edm4hep::ConstCalorimeterHit> clu ) +{ + std::vector<int> hitlayers; + hitlayers.clear(); + + int NHits = clu.size(); + int tmpK = 0; //Layer Number + int tmpS = 0; + int tmpID = 0; + + for(int i = 0; i < NHits; i++) + { + auto hit = clu[i]; + auto cellid= hit.getCellID(); + tmpK = m_decoder->get(cellid, "layer")+1 ; + tmpS = m_decoder->get(cellid, "stave")+1 ; + // cout<<"tmpK "<<tmpK<<endl; + tmpID = tmpS * 50 + tmpK; + + if( std::find(hitlayers.begin(), hitlayers.end(), tmpID) == hitlayers.end() ) + { + hitlayers.push_back(tmpID); + } + } + + return hitlayers.size(); +} + + +float ArborToolLCIO::ClusterT0(edm4hep::ConstCluster a_Clu) +{ + 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); + CluHitPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + + tmpTime = a_hit.getTime() - 1.0/300*CluHitPos.Mag(); + if(tmpTime < T0 && tmpTime > 0) + { + T0 = tmpTime; + } + } + return T0; +} + + +int ArborToolLCIO::newPhotonTag(edm4hep::ConstCluster a_clu) +{ + int flag=0; + + TVector3 PosClu = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + float Depth = 0; + Depth = DisSeedSurface(PosClu); + + int CluSize= a_clu.hits_size(); + float ClusFD=FDV3(a_clu); + float ClusT0=ClusterT0(a_clu); + + + 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; + + + +} + + +int ArborToolLCIO::ClusterFlag1st(edm4hep::ConstCluster a_tree) +{ + int ClusterID = 211; + int EcalNHit, HcalNHit, NH_ECALF10; + 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); + + EClu = a_tree.getEnergy(); + EcalNHit = 0; + HcalNHit = 0; + EcalEn = 0; + HcalEn = 0; + float currDepth = 0; + maxDepth = -100; + minDepth = 1E6; + + std::vector<edm4hep::ConstCalorimeterHit> allhits; + std::vector<edm4hep::ConstCalorimeterHit> Ecalhits; + std::vector<edm4hep::ConstCalorimeterHit> Hcalhits; + std::vector<edm4hep::ConstCalorimeterHit> Ecalf10hits; + + allhits.clear(); + Ecalhits.clear(); + Hcalhits.clear(); + Ecalf10hits.clear(); + + std::vector<float> hitTheta; + hitTheta.clear(); + + for(unsigned int j1 = 0; j1 < a_tree.hits_size(); j1++) + { + auto a_hit = a_tree.getHits(j1); + + TVector3 tmpPos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + hitTheta.push_back(tmpPos.Theta()); + + } + + float totTheta = 0; + float avTheta = 0; + float SDTheta; + + for(int t0 = 0; t0 < int(hitTheta.size()); t0++) + { + float tmpTheta = hitTheta[t0]; + totTheta += tmpTheta; + } + + avTheta = totTheta/float(hitTheta.size()); + SDTheta = 0; + + for(int t1 = 0; t1 < int(hitTheta.size()); t1++) + { + float tmpTheta = hitTheta[t1]; + SDTheta += pow(tmpTheta-avTheta,2); + } + SDTheta = sqrt(SDTheta/float(hitTheta.size())); + + TVector3 HitPos; + int currCluNHits = a_tree.hits_size(); + if(currCluNHits == 0) return 1; + int index1 = 0, index2 = 0; + int HitDepth50 = 1.0e4; + + for(int s1 = 0; s1 < currCluNHits; s1++) + { + auto a_hit = a_tree.getHits(s1); + allhits.push_back(a_hit); + auto cellid = a_hit.getCellID(); + int 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); + + currDepth = DisSeedSurface(HitPos); + if(tmpHitEn > 0.05 && currDepth < HitDepth50) + { + HitDepth50 = currDepth; + } + + + if(currDepth > maxDepth) + { + maxDepth = currDepth; + index1 = s1; + } + if(currDepth < minDepth) + { + minDepth = currDepth; + index2 = s1; + } + + if( fabs(a_hit.getEnergy() - DHCALCalibrationConstant) < 1.0E-6 ) //or other fancy judgements...^M + { + HcalNHit++; + HcalEn += a_hit.getEnergy(); + Hcalhits.push_back(a_hit); + } + else + { + EcalNHit++; + EcalEn += a_hit.getEnergy(); + Ecalhits.push_back(a_hit); + if(NLayer< 10) Ecalf10hits.push_back(a_hit); + } + } + auto maxdis_hit = a_tree.getHits(index1); + auto mindis_hit = a_tree.getHits(index2); + + TVector3 maxpos = TVector3(maxdis_hit.getPosition().x,maxdis_hit.getPosition().y,maxdis_hit.getPosition().z); + TVector3 minpos = TVector3(mindis_hit.getPosition().x,mindis_hit.getPosition().y,mindis_hit.getPosition().z); + TVector3 GraPos = ClusterCoG(a_tree); + cluDepth = (maxpos-minpos).Mag(); + + float totHitEn = 0; + float totHitEnDis = 0; + float HitEn; + + for(int s2 = 0; s2 < currCluNHits; s2++) + { + auto a_hit2 = a_tree.getHits(s2); + HitPos = TVector3(a_hit2.getPosition().x,a_hit2.getPosition().y,a_hit2.getPosition().z); + HitEn = a_hit2.getEnergy(); + TVector3 par1 = GraPos-minpos; + TVector3 par2 = minpos-HitPos; + TVector3 par3 = par1.Cross(par2); + float disHtoL = par3.Mag()/par1.Mag(); + totHitEn+=HitEn; + totHitEnDis+=HitEn*disHtoL; + } + avEnDisHtoL = totHitEnDis/totHitEn; + FD_all = FDV2(allhits); + FD_ECAL = FDV2(Ecalhits); + FD_HCAL = FDV2(Hcalhits); + FD_ECALF10 = FDV2(Ecalf10hits); + NH_ECALF10 = Ecalf10hits.size(); + + + bool cute1 = log10(FD_all/log10(EClu)) > (log10(cluDepth/log10(EClu))-2.9); + bool cute2 = FD_ECAL > 0.1 && FD_ECAL > 0.2*log10(avEnDisHtoL-3)+0.05*sqrt(avEnDisHtoL-3); + float x; + if (EcalNHit == 0) x = 0; + else x = float(NH_ECALF10)/float(EcalNHit); + 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 cute; + cute = cute1 && cute2 && cute3 && cute4; + + bool cutmu1 = (cluDepth > 1300 || EClu < 4/1000.*cluDepth+0.9); + bool cutmu2 = EClu < 15; + bool cutmu3 = avEnDisHtoL*SDTheta/EClu <0.02 && FD_all < 0.4/0.025*(0.025-avEnDisHtoL*SDTheta/EClu); + bool cutmu; + + cutmu = cutmu1 && cutmu2 && cutmu3; + + + if(currCluNHits <= 4) ClusterID = 1; + else if(cute) ClusterID = 11; + else if (cutmu) ClusterID = 13; + else + { + bool cutef1 = FD_HCAL == -1; + bool cutef2 = minDepth < 50; + bool cutef3 = FD_ECAL > 0.2; + bool cutef3b = 1; + if (FD_ECAL < 0.6) cutef3b = FD_ECAL > 6/3.5*(1.75-log10((EcalNHit+HcalNHit)/EClu)); + bool cutef4; + if (avEnDisHtoL <= 8) + cutef4 = FD_ECALF10 > 0.2; + else + cutef4 = FD_ECALF10 > 0.2+0.8*sqrt(log10(avEnDisHtoL/8)); + bool cutef5 = FD_ECAL/EClu > (log10(EcalNHit)-1.6)*(log10(EcalNHit)-1.6)*4+0.25 && FD_ECALF10 > (avEnDisHtoL-9)*(avEnDisHtoL-9)*0.006+0.2; + bool cutef; + + cutef = (cutef1 && cutef2 && cutef3 && cutef3b && cutef4) || cutef5; + if(cutef) ClusterID = 11; + } + + if(ClusterID == 211) + { + bool cutmuf1 = FD_all == 0 && log10((EcalNHit+2*HcalNHit)*EClu) > 1.2; + bool cutmuf2 = log10((EcalNHit+2*HcalNHit)*EClu) > 1.9 && log10((EcalNHit+2*HcalNHit)*EClu) < 3.1; + bool cutmuf3 = log10(FD_all/cluDepth) < -223.074+431.315*log10((EcalNHit+2*HcalNHit)*EClu)-331.72*pow(log10((EcalNHit+2*HcalNHit)*EClu),2)+124.588*pow(log10((EcalNHit+2*HcalNHit)*EClu),3)-22.8786*pow(log10((EcalNHit+2*HcalNHit)*EClu),4)+1.64577*pow(log10((EcalNHit+2*HcalNHit)*EClu),5); + bool cutmuf = cutmuf1 || (cutmuf2 && cutmuf3); + if(cutmuf) ClusterID = 13; + else if(minDepth < 0.77+0.23*EClu) ClusterID = 211; + else ClusterID = 2; + } + + return ClusterID; +} + +float ArborToolLCIO::ClusterEE(edm4hep::ConstCluster inputCluster) +{ + float ClusterEnergy = 0; + float tmpCluEn = 0; + int CluType = 211; + if(ClusterFlag1st(inputCluster) == 11) CluType = 22; + + int NCluHit = inputCluster.hits_size(); + float hitEn = 0; + float EnCorrector = 1; + + if(CluType == 22) + { + ClusterEnergy = EMClusterEE(inputCluster); + } + else + { + for(int i = 0; i < NCluHit; i++) + { + auto a_hit = inputCluster.getHits(i); + hitEn = a_hit.getEnergy(); + + if(CluType == 211) + { + if( hitEn < 1.5 ) // Or some function of the initCluEn + { + tmpCluEn += hitEn; + } + else + { + cout<<"Ultra Hot Hit"<<endl; // Use ShuZhen's Hot Hit Finding function... + tmpCluEn += 0.05; // MIP Hit Energy, Value to be determined + } + } + else // For EM & MIP: should veto accordingly + { + tmpCluEn += hitEn; + } + } + + EnCorrector = 1; + if(tmpCluEn > 1.5 && tmpCluEn < 22 && 1 == 0) + { + EnCorrector = 0.6*(1 + 1.0/log10(tmpCluEn)); + } + ClusterEnergy = tmpCluEn*EnCorrector; + //ClusterEnergy = HADClusterEE(tmpCluEn,inputCluster); + } + + return ClusterEnergy; +} + + +float ArborToolLCIO::EMClusterEE( edm4hep::ConstCluster inputCluster ) +{ + float aaa = -50.2288, ab = 219.398,ac =0.17679,ad = 0.00241144; + float ba = -56.6164,bbb = 162.647,bc = 0.679974,bd = 0.00423267,be = -0.324786; + float ff1a = -21.878,ff1b = 1146.3, ff1c = 0.0267898, ff1d = 0.000712903; + float ff2a = -85.8291,ff2b = 2466.59,ff2c = 0.55722, ff2d = 0.00159572; + float ArEhito10=0, ArEhite10=0, ArEhito20=0, ArEhite20m=0, ArEhite20p=0; + int ArNhito10=0, ArNhite10=0, ArNhito20=0, ArNhite20m=0, ArNhite20p=0; + float diff=1.0; + float x = 0, y= 0, f1 = 0, f2 = 0; + 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 CluTheta = CluPos.Theta(); + 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"); + if(NLayer > 20){ + if (NLayer%2 ==0){ + ArNhite10 ++; + ArEhite10 +=a_hit.getEnergy(); + } + else if (NLayer%2 ==1){ + ArNhito10 ++; + ArEhito10 +=a_hit.getEnergy(); + } + + } + else if(NLayer <= 20){ + if (NLayer%2 ==0){ + if(NLayer ==0){ + ArNhite20m ++; + ArEhite20m +=a_hit.getEnergy(); + } + if(NLayer >0){ + ArNhite20p ++; + ArEhite20p +=a_hit.getEnergy(); + } + } + else if (NLayer%2 ==1){ + ArNhito20 ++; + ArEhito20 +=a_hit.getEnergy(); + } + } + } + while(diff>0.01 && EMC > 0) + { + float temp_a=sqrt(EMC*EMC); + x=exp((-EMC+aaa)/ab)+ac/sqrt(EMC)+ad*EMC; + y=1/(exp((-EMC+ba)/bbb)+bc/sqrt(EMC)+bd*EMC)+be; + f1=1/(exp((-EMC+ff1a)/ff1b)+ff1c/sqrt(EMC)+ff1d*EMC); + f2=1/(exp((-EMC+ff2a)/ff2b)+ff2c/sqrt(EMC)+ff2d*EMC); + EMC=x*ArEhito20+f1*ArEhite20p+y*ArEhito10+f2*ArEhite10; + float temp_b=sqrt(EMC*EMC); + diff=fabs(temp_a-temp_b); + } + _costheta=cos(CluTheta); + + if(EMC/inputCluster.getEnergy() < 0.5) EMC = inputCluster.getEnergy(); + + Ethetacorr=(50/(0.294596*fabs(_costheta)*fabs(_costheta)-1.58336*fabs(_costheta)+49.9219)); + + float ModuleCluPhi = 0; + if(CluPhi > 17.5) + { + ModuleCluPhi = CluPhi - int((CluPhi -17.5)/45)*45; + } + else + { + ModuleCluPhi = CluPhi - int((CluPhi -17.5)/45)*45 + 45; + } + + if(ModuleCluPhi>=17.5 && ModuleCluPhi<20.5 ) Ephicorr=(50/(-0.122*ModuleCluPhi*ModuleCluPhi+2.76*ModuleCluPhi+38.04)); + if(ModuleCluPhi>=20.5 && ModuleCluPhi<22.5 ) Ephicorr=(50/(0.22*ModuleCluPhi*ModuleCluPhi-6.43*ModuleCluPhi+81.69)); + if(ModuleCluPhi>=22.5 && ModuleCluPhi<62.5 ) Ephicorr=50*(0.00839675/sqrt(ModuleCluPhi)+0.000369287*logf(ModuleCluPhi)+0.0174033); + + EMC=EMC*Ethetacorr*Ephicorr; + + return EMC; +} + +std::vector<float> ArborToolLCIO::ClusterTime(edm4hep::ConstCluster inputCluster) +{ + std::vector<float> CluTimeVector; + CluTimeVector.clear(); + + int NHit = inputCluster.hits_size(); + float currhitTime = 0; + float currhitoriTime = 0; + TVector3 hitpos; + + std::vector<float> Time; + Time.clear(); + std::map<int, float> CluHitToTime; + CluHitToTime.clear(); + std::vector<float> OriginalTime; + OriginalTime.clear(); + + int N_PropTimeHit = 0; + + for(int i = 0; i < NHit; i++) + { + auto a_hit =inputCluster.getHits(i); + hitpos = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + if(a_hit.getTime() == 0) + { + currhitTime = 1001; + currhitoriTime = 1001; + } + else + { + currhitTime = a_hit.getTime(); + currhitoriTime = a_hit.getTime(); + N_PropTimeHit++; + } + Time.push_back( currhitTime ); + OriginalTime.push_back( currhitoriTime ); + 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; + + for(int j0= 0; j0 < N_PropTimeHit; j0++) + { + auto a_hit = inputCluster.getHits(OriTimeOrder[j0]); + hitP = TVector3(a_hit.getPosition().x,a_hit.getPosition().y,a_hit.getPosition().z); + if(N_PropTimeHit > 10) + { + if( j0 < 5) + StartP += 0.2*hitP; + else if(j0 > 4 && j0 < 10) + EndP += 0.2*hitP; + } + else + { + if( j0 < N_PropTimeHit/2.) + StartP += N_PropTimeHit/2.*hitP; + else + EndP += N_PropTimeHit/2.*hitP; + } + } + + for(int j = 0; j < NHit; j++) + { + ptime = CluHitToTime[TimeOrder[j]]; + + if(j == 0) + { + CluTimeVector.push_back(ptime); + } + if(j == NHit - 1) + CluTimeVector.push_back(ptime); + + if(j > 0) + { + //cout<<"ptime "<<ptime<<" cluhittotime "<<CluHitToTime[TimeOrder[j - 1]]<<endl; + if( (ptime - CluHitToTime[TimeOrder[j - 1]]) < 3) + { + if( Break == 0 ) + { + PeakTime += ptime; + NCount ++; + } + } + else + { + Break = 1; + } + } + + AverageTime += ptime; + } + if(NCount) + PeakTime = float(PeakTime)/NCount; + + 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 + CluTimeVector.push_back((EndP + StartP).Angle(EndP - StartP)); + else + CluTimeVector.push_back(1001); + + return CluTimeVector; +} + + + + diff --git a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh new file mode 100644 index 0000000000000000000000000000000000000000..d3aaeea19c0d37a5fb7c3062c6ce72f74bbf3f72 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.hh @@ -0,0 +1,131 @@ +#ifndef ARBORTOOLLCIO_H_ +#define ARBORTOOLLCIO_H_ + + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/ReconstructedParticleCollection.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/TrackerHitCollection.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/VertexCollection.h" +#include "edm4hep/TrackCollection.h" +#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCRecoCaloAssociation.h" +#include "edm4hep/MCRecoTrackerAssociation.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCRecoParticleAssociation.h" +#include "edm4hep/MCRecoParticleAssociationCollection.h" + + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +#include "TVector3.h" +#include "GaudiKernel/INTupleSvc.h" +#include "GaudiKernel/NTuple.h" +#include "GaudiKernel/MsgStream.h" +#include "ArborTool.h" +#include "HelixClassD.hh" + +class CollectionMaps +{ +public: + CollectionMaps(); + void clear(); + std::map<std::string, std::vector<edm4hep::MCParticle> > collectionMap_MC; + std::map<std::string, std::vector<edm4hep::ConstCalorimeterHit> > collectionMap_CaloHit; + std::map<std::string, std::vector<edm4hep::Vertex> > collectionMap_Vertex; + std::map<std::string, std::vector<edm4hep::Track> > collectionMap_Track; + std::map<std::string, std::vector<edm4hep::MCRecoCaloAssociation> > collectionMap_CaloRel; + std::map<std::string, std::vector<edm4hep::MCRecoTrackerAssociation> > collectionMap_TrkRel; +}; + + + +class ArborToolLCIO : public GaudiAlgorithm +{ + public: + ArborToolLCIO(const std::string& name, ISvcLocator* svcLoc); + //ArborToolLCIO(ISvcLocator* svcLoc); + ~ArborToolLCIO(); + typedef DataHandle<edm4hep::Cluster> ClusterType; + + + void Clean(); + void ClusterBuilding(DataHandle<edm4hep::ClusterCollection>& _currbranchcoll, std::vector<edm4hep::ConstCalorimeterHit> Hits, branchcoll BranchOrder, int DHCALFlag); + + float ClusterT0(edm4hep::ConstCluster a_Clu); + + int NHScaleV2( std::vector<edm4hep::ConstCalorimeterHit> clu0, int RatioX, int RatioY, int RatioZ ); + float FDV2( std::vector<edm4hep::ConstCalorimeterHit> clu); + + int NHScaleV3( edm4hep::ConstCluster clu0, int RatioX, int RatioY, int RatioZ ); + + float FDV3( edm4hep::ConstCluster clu); + + int ActiveLayers( std::vector<edm4hep::ConstCalorimeterHit> clu ); + + float BushDis( edm4hep::ConstCluster clu1, edm4hep::ConstCluster clu2); + + + float* SimpleDisTrackClu(edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu); + + float SimpleBushTimeTrackClu(edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu); + + int SimpleBushNC(edm4hep::ConstTrack a_trk, edm4hep::ConstCluster a_clu); + + float DisPointToBush( TVector3 Pos1, edm4hep::ConstCluster clu1); + + TVector3 ClusterCoG(edm4hep::ConstCluster inputCluser); + + edm4hep::ClusterCollection* ClusterVecMerge( std::vector<edm4hep::ConstCluster> inputClusters, TMatrixF ConnectorMatrix, DataHandle<edm4hep::ClusterCollection>& clucol ); + + edm4hep::ClusterCollection* ClusterVecColl( std::vector<edm4hep::ConstCluster> inputClusters, DataHandle<edm4hep::ClusterCollection>& m_clucol); + + edm4hep::Cluster NaiveCluImpl(edm4hep::ConstCluster a0_clu); + void NaiveCluConst(edm4hep::ConstCluster a0_clu, edm4hep::Cluster); + + std::vector<edm4hep::ConstCluster> CollClusterVec(const edm4hep::ClusterCollection * input_coll ); + + std::vector<edm4hep::ConstCalorimeterHit> CollHitVec(const edm4hep::CalorimeterHitCollection * input_coll, float EnergyThreshold); + + std::vector<edm4hep::Cluster> ClusterHitAbsorbtion( std::vector<edm4hep::ConstCluster> MainClusters, std::vector<edm4hep::ConstCalorimeterHit> IsoHits, float DisThreshold ); + + edm4hep::Cluster NaiveMergeClu(std::vector<edm4hep::ConstCluster> inputCluVec); + + void NaiveMergeCluConst(std::vector<edm4hep::ConstCluster> inputCluVec,edm4hep::Cluster MergedClu); + std::vector<edm4hep::ConstCluster> ClusterAbsorbtion( std::vector<edm4hep::ConstCluster> MainClusters, std::vector<edm4hep::ConstCluster> FragClusters, float thresholds, float DepthSlope ); + + + int JointsBetweenBush(edm4hep::ConstCluster a_Clu, edm4hep::ConstCluster b_Clu, float CellSize); + + TVector3 CluEndP(edm4hep::ConstCluster a_clu); + + int ClusterFlag(edm4hep::ConstCluster a_tree, edm4hep::ConstTrack a_trk); + + int ClusterFlag1st(edm4hep::ConstCluster a_tree); + + int newPhotonTag(edm4hep::ConstCluster a_clu); + float ClusterEE(edm4hep::ConstCluster inputCluster); + + float EMClusterEE( edm4hep::ConstCluster inputCluster ); + + std::vector<float> ClusterTime(edm4hep::ConstCluster inputCluster); + + private: + SmartIF<IGeomSvc> m_geosvc;//=service<IGeomSvc>("GeomSvc"); + dd4hep::DDSegmentation::BitFieldCoder* m_decoder;// = m_geosvc->getDecoder("ECALBarrel"); + + +// m_geosvc= service<IGeomSvc>("GeomSvc"); +// m_decoder = m_geosvc->getDecoder("ECALBarrel"); + + +}; +#endif // diff --git a/Reconstruction/PFA/Arbor/src/BushConnect.cc b/Reconstruction/PFA/Arbor/src/BushConnect.cc new file mode 100644 index 0000000000000000000000000000000000000000..83559f37af14ff5e80a228c31ef25cdf2039fd17 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/BushConnect.cc @@ -0,0 +1,1328 @@ +#include "BushConnect.hh" +#include "ArborTool.h" +#include "ArborToolLCIO.hh" +#include "DetectorPos.hh" + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/Cluster.h" +#include "edm4hep/ReconstructedParticle.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/ReconstructedParticleCollection.h" + +#include "cellIDDecoder.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/IDDescriptor.h" +#include "DD4hep/Plugins.h" + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + + +#include <values.h> +#include <string> +#include <iostream> +#include <cmath> +#include <stdexcept> +#include <TFile.h> +#include <TTree.h> +#include <TMath.h> +#include <Rtypes.h> +#include <sstream> +#include <set> +#include <TVector3.h> +#include <vector> +#include <algorithm> + +#include "HelixClassD.hh" + +using namespace std; + +//const string ECALCellIDDecoder = "M:3,S-1:3,I:9,J:9,K-1:6"; + +DECLARE_COMPONENT(BushConnect) + +BushConnect::BushConnect(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) +{ + +} + + +StatusCode BushConnect::initialize() { + + ISvcLocator* svcloc = serviceLocator(); + m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc); +//printParameters(); + //Cluflag.setBit(LCIO::CHBIT_LONG); + return GaudiAlgorithm::initialize(); +} + +void BushConnect::Clean(){ + MCPTrack_Type.clear(); + Track_Energy.clear(); + Track_Type.clear(); + Track_Theta.clear(); + Track_EndPoint.clear(); + Track_P3.clear(); + TrackStartPoint.clear(); + SortedTracks.clear(); + ChCoreID.clear(); + chargedclustercore.clear(); + selfmergedcluster.clear(); + non_chargedclustercore.clear(); //all clusters besides charged cluster core + ClusterType_1stID.clear(); + CluFD.clear(); + CluT0.clear(); + CluCoG.clear(); + Clu_Depth.clear(); +} + +void BushConnect::TrackSort() //, &std::map<Track*, int>Track_Tpye, &std::map<Track*, float> Track_Energy) +{ + try{ + auto TrackColl = m_trkcol.get(); + + + int NTrk = TrackColl->size(); + cout<<NTrk<<endl; + float D0 = 0; + float Z0 = 0; + int NTrkHit = 0; + const float mass = 0.139; //Pion Mass + TVector3 EndPointPos, StartPointPos; + int TrackType = 0; + + std::vector<edm4hep::ConstTrack> tracks_HQ_Barrel; + std::vector<edm4hep::ConstTrack> tracks_HQ_Endcap; + std::vector<edm4hep::ConstTrack> tracks_HQ_Shoulder; + std::vector<edm4hep::ConstTrack> tracks_HQ_Forward; + std::vector<edm4hep::ConstTrack> tracks_MQ_Barrel; + std::vector<edm4hep::ConstTrack> tracks_MQ_Endcap; + std::vector<edm4hep::ConstTrack> tracks_MQ_Shoulder; + std::vector<edm4hep::ConstTrack> tracks_MQ_Forward; + std::vector<edm4hep::ConstTrack> tracks_Vtx; + std::vector<edm4hep::ConstTrack> tracks_LQ; + std::vector<edm4hep::ConstTrack> tracks_LE; + std::vector<edm4hep::ConstTrack> curr_tracks; + + Track_EndPoint.clear(); + + tracks_HQ_Barrel.clear(); + tracks_HQ_Endcap.clear(); + tracks_HQ_Shoulder.clear(); + tracks_HQ_Forward.clear(); + tracks_MQ_Barrel.clear(); + tracks_MQ_Endcap.clear(); + tracks_MQ_Shoulder.clear(); + tracks_MQ_Forward.clear(); + tracks_Vtx.clear(); + tracks_LQ.clear(); + tracks_LE.clear(); + + std::vector<edm4hep::ConstTrack> tracks_ILL; + tracks_ILL.clear(); + std::vector<edm4hep::ConstTrack> tracks_preInteraction; + 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::ConstTrack, int> Track_Index; + Track_Index.clear(); + Track_Energy.clear(); + Track_Type.clear(); + Track_P3.clear(); + Track_EndPoint.clear(); + TrackStartPoint.clear(); + + 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); + 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; + + delete TrkInit_Helix; + } + + TVector3 currEp, currSp; + float currMotherEn = 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 + { + currMotherEn = Track_Energy[a_Trk]; + sumDauEn = 0; + for(int i2 = 0; i2 < NTrk; i2++) + { + auto b_Trk = (*TrackColl)[i2]; + if(i2 != i1) + { + currSp = TrackStartPoint[b_Trk]; + if( (currEp - currSp).Mag() < 40 ) + sumDauEn += Track_Energy[b_Trk]; + } + } + if(currMotherEn + 0.1 > 0.9*sumDauEn && currMotherEn > 3 && sumDauEn > 0 ) //Some protection is always needed... + { + tracks_preInteraction.push_back(a_Trk); + } + } + } + + for(int t0 = 0; t0 < NTrk; t0++) + { + auto a_Trk = (*TrackColl)[t0]; + D0 = a_Trk.getTrackStates(0).D0; + Z0 = a_Trk.getTrackStates(0).Z0; + NTrkHit = a_Trk.getTrackerHits().size(); + auto last_hit = a_Trk.getTrackerHits(NTrkHit - 1); + EndPointPos = TVector3(last_hit.getPosition().x,last_hit.getPosition().y,last_hit.getPosition().z); + 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; + } + + 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; + } + else if( fabs(EndPointPos.Z()) > ECALHalfZ - 500 && EndPointPos.Perp() > TPCOuterRadius - 300 ) //Shoulder + { + TrackType = 30; + } + else if( fabs(EndPointPos.Z()) > LStar - 500 && EndPointPos.Perp() < TPCInnerRadius ) //Forward + { + TrackType = 40; + } + else if( EndPointPos.Perp() > TPCOuterRadius - 100 ) //Barrel + { + TrackType = 10; + } + else if( fabs(EndPointPos.Z()) > ECALHalfZ - 200 ) //Endcap + { + TrackType = 20; + } + + if( fabs(D0) < 1 && fabs(Z0) < 1 ) + { + TrackType += 1; + } + + Track_Type[a_Trk] = TrackType; + + if(TrackType == 11) + tracks_HQ_Barrel.push_back(a_Trk); + else if(TrackType == 21) + tracks_HQ_Endcap.push_back(a_Trk); + else if(TrackType == 31) + tracks_HQ_Shoulder.push_back(a_Trk); + else if(TrackType == 41) + tracks_HQ_Forward.push_back(a_Trk); + else if(TrackType == 10) + tracks_MQ_Barrel.push_back(a_Trk); + else if(TrackType == 20) + tracks_MQ_Endcap.push_back(a_Trk); + else if(TrackType == 30) + tracks_MQ_Shoulder.push_back(a_Trk); + else if(TrackType == 40) + tracks_MQ_Forward.push_back(a_Trk); + else if(TrackType == 1) + tracks_Vtx.push_back(a_Trk); + else if(TrackType == 101) + tracks_LE.push_back(a_Trk); + else if( (StartPointPos.Mag() > 50 && EndPointPos.Mag() < 1000 && NTrkHit < 50) || TrackType == 100 ) + { tracks_ILL.push_back(a_Trk); + } + else + tracks_LQ.push_back(a_Trk); + } + } + cout<<"LQ"<<tracks_LQ.size()<<" "<<tracks_ILL.size()<<endl; + + std::vector<float > currTrkMomentum; + std::vector<int> currTrkIndex; + + for(int t1 = 0; t1 < 11; t1++) + { + currTrkMomentum.clear(); + currTrkIndex.clear(); + curr_tracks.clear(); + if(t1 == 0) + curr_tracks = tracks_HQ_Endcap; + else if(t1 == 1) + curr_tracks = tracks_HQ_Barrel; + else if(t1 == 2) + curr_tracks = tracks_MQ_Endcap; + else if(t1 == 3) + curr_tracks = tracks_MQ_Barrel; + else if(t1 == 4) + curr_tracks = tracks_HQ_Shoulder; + else if(t1 == 5) + curr_tracks = tracks_MQ_Shoulder; + else if(t1 == 6) + curr_tracks = tracks_HQ_Forward; + else if(t1 == 7) + 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; + + + 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]]; + if(t1 < 9 || Track_Energy[b_tmpTrk] < 10) + TrackOrder.push_back(Track_Index[b_tmpTrk]); + } + } + + for(unsigned int t4 = 0; t4 < TrackOrder.size(); t4++) + { + auto b_Trk = (*TrackColl)[t4]; + SortedTracks.push_back(b_Trk); + } + }catch(GaudiException &e){} +} + +void BushConnect::BushSelfMerge() +{ + auto CaloClu = m_clucol.get(); + int NClu = CaloClu->size(); + + std::vector<edm4hep::ConstCluster > Core_1st; + std::vector<edm4hep::ConstCluster > Frag_1st; + std::vector<edm4hep::ConstCluster > UnId_1st; + Core_1st.clear(); + Frag_1st.clear(); + UnId_1st.clear(); + + float CluDepth = 0; + float CluEn = 0; + int CluSize = 0; + TVector3 PosCluSeed, PosSeedDiff, PosCoGDiff, PosSeedA, PosSeedB; + + 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; + + TMatrixF FlagMerge(NClu, NClu); + + cout<<NClu<<" clusters"<<endl; + for(int i0 = 0; i0 < NClu; i0++) + { + auto a_clu = (*CaloClu)[i0]; + CluFD[a_clu] = m_ArborToolLCIO->FDV3(a_clu); + CluT0[a_clu] = m_ArborToolLCIO->ClusterT0(a_clu); + 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); + } + // CluCoG_Top20Percent Might be used to improve Photon Split Remerge performance + + for(int s0 = 0; s0 < NClu; s0++) + { + auto a_clu = (*CaloClu)[s0]; + PosSeedA = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + Depth_A = Clu_Depth[a_clu]; + Size_A = a_clu.hits_size(); + + for(int s1 = s0 + 1; s1 < NClu; s1++) + { + auto b_clu = (*CaloClu)[s1]; + NJoints = m_ArborToolLCIO->JointsBetweenBush(a_clu, b_clu, 10); + PosSeedB = TVector3(b_clu.getPosition().x,b_clu.getPosition().y,b_clu.getPosition().z); + Depth_B = Clu_Depth[b_clu]; + Size_B = b_clu.hits_size(); + + PosSeedDiff = PosSeedA - PosSeedB; + DeeperDepth = std::max(Depth_A, Depth_B); + LaterT0 = std::max(CluT0[a_clu], CluT0[b_clu]); + PosCoGDiff = CluCoG[a_clu] - CluCoG[b_clu]; + + if(NJoints&& PosSeedDiff.Mag() < 40 + 0.05*DeeperDepth ) //And depth... + { + SmallCluSize = std::min( Size_A, Size_B ); + + //if( SmallCluSize < 10 || LaterT0 > 10 || DeeperDepth > 40 || (NJoints > 8 && PosCoGDiff.Mag()*PosCoGDiff.Angle(PosSeedB) < 15 )) + 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; + } + } + //Head Tail Connection. Could be more sophsticate && should be very strict. + if( PosSeedA.Angle(PosSeedB) < 0.1 && PosSeedDiff.Mag() < 1000 && PosSeedDiff.Mag()*PosSeedA.Angle(PosSeedB) < 60 + 0.02*DeeperDepth && ((CluFD[a_clu] < 0.2 && Size_A > 6) || (CluFD[b_clu] < 0.2 && Size_B > 6)) ) + { + if( (PosSeedA.Mag() > PosSeedB.Mag() && PosSeedA.Angle(PosSeedB - PosSeedA) < 0.2) || (PosSeedB.Mag() > PosSeedA.Mag() && PosSeedA.Angle(PosSeedA - PosSeedB) < 0.2) ) + { + FlagMerge[s0][s1] = 2.0; + FlagMerge[s1][s0] = 2.0; + } + } + } + } + + std::vector<edm4hep::ConstCluster> OriInputEHBushes = m_ArborToolLCIO->CollClusterVec(CaloClu); + TMatrixF MergeSYM = MatrixSummarize(FlagMerge); + auto CloseMergedCaloClu = m_ArborToolLCIO->ClusterVecMerge( OriInputEHBushes, MergeSYM, clucol_merged); + + std::map<edm4hep::ConstCluster,float> MinDisSeedToBush; + MinDisSeedToBush.clear(); + for(int i0 = 0; i0 < CloseMergedCaloClu->size(); i0++) + { + auto a_clu = (*CloseMergedCaloClu)[i0]; + PosCluSeed = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + float tmpmindis = 1e10; + for(int i1 = 0; i1 < CloseMergedCaloClu->size(); i1++) + { + auto b_clu = (*CloseMergedCaloClu)[i1]; + if(i1 != i0) + { + if(m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu) < tmpmindis) + tmpmindis = m_ArborToolLCIO->DisPointToBush(PosCluSeed,b_clu); + } + } + MinDisSeedToBush[a_clu] = tmpmindis; + } + + for(int i0 = 0; i0 < CloseMergedCaloClu->size(); i0++) + { + auto a_clu = (*CloseMergedCaloClu)[i0]; + PosCluSeed = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + CluDepth = DisSeedSurface(PosCluSeed); + CluEn = a_clu.getEnergy(); + CluSize = a_clu.hits_size(); + + if( CluSize > 10 && ( (CluEn > 2.0 + 0.002*CluDepth && CluDepth < 22 ) || CluEn > 5.0) ) + { + Core_1st.push_back(a_clu); + } + else if( (CluSize > 10 || CluEn > 0.2) && MinDisSeedToBush[a_clu] > 20 && m_ArborToolLCIO->ClusterT0(a_clu) < 1.2 ) // && (PhotonTag(a_clu) == 1 || ClusterFlag1st(a_clu) == 11 )) + { + Core_1st.push_back(a_clu); + } + else if( CluSize < 5 && CluEn < 0.3 && CluDepth > 40 ) + { + Frag_1st.push_back(a_clu); + } + else + { + UnId_1st.push_back(a_clu); + } + } + + std::vector<edm4hep::ConstCluster > 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); +} + +void BushConnect::TagCore() +{ + int NTrk = SortedTracks.size(); + int NClu = selfmergedcluster.size(); + int currTrackType = 0; + float currTrkEn = 0; + float DisMatrix_Track_Clu_E[NTrk][NClu]; + float TimeMatrix_Track_Clu_E[NTrk][NClu]; + float CluDepth = 0; + float CoreMergeDistanceDepthCorrector = 0; + + TVector3 TrkEndPoint(0, 0, 0); + TVector3 CluPos(0, 0, 0); + std::map<edm4hep::ConstCluster, int> BushTouchFlag; + std::map<edm4hep::ConstTrack, edm4hep::ConstCluster> FCMap_Track_CHCore; + std::map<edm4hep::ConstTrack, std::vector<edm4hep::ConstCluster>> FCMap_Track_CHCore_new; + std::map<int, int> Closest_Trk_Clu_Map; + std::vector<edm4hep::ConstCluster> TightLinkedCluster; + TightLinkedCluster.clear(); + Closest_Trk_Clu_Map.clear(); + BushTouchFlag.clear(); + FCMap_Track_CHCore.clear(); + FCMap_Track_CHCore_new.clear(); + + Clu_Depth.clear(); + + for(int s0 = 0; s0 < NTrk; s0++) + { + for(int s1 = 0; s1 < NClu; s1++) + { + DisMatrix_Track_Clu_E[s0][s1] = 1.0E10; + TimeMatrix_Track_Clu_E[s0][s1] = 1.0E10; + } + } + + 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); + } + + //~~~~~~~ find the closest cluster first... + + for(int g0 = 0; g0 < NTrk; g0++) + { + auto a_trk = SortedTracks[g0]; + float ClosestDis = 1.0E9; + int ClosestCluIndex = -1; + int ClosestNC = 1E9; + float TrackEn = Track_Energy[a_trk]; + TrkEndPoint = Track_EndPoint[a_trk]; + + currTrackType = Track_Type[a_trk]; + TVector3 TrkP3 = Track_P3[a_trk]; + + for(int g1 = 0; g1 < NClu; g1++) + { + auto fccand_bush = selfmergedcluster[g1]; + float* Dis = m_ArborToolLCIO->SimpleDisTrackClu(a_trk, fccand_bush); + float Time = m_ArborToolLCIO->SimpleBushTimeTrackClu(a_trk, fccand_bush); + int NC = m_ArborToolLCIO->SimpleBushNC(a_trk, fccand_bush); + if(Dis[2] > -0.1) + { + DisMatrix_Track_Clu_E[g0][g1] = Dis[2]; + TimeMatrix_Track_Clu_E[g0][g1] = Time; + if( Dis[2] < ClosestDis ) // && ThetaDiff < 0.05 + 0.1/Track_Energy[a_trk] ) + { + ClosestDis = Dis[2]; + ClosestCluIndex = g1; + ClosestNC = NC; + } + } + } + + //Diag for mimic + 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 ) ) + { + auto candiclu= selfmergedcluster[ClosestCluIndex]; + CluPos = TVector3(candiclu.getPosition().x,candiclu.getPosition().y,candiclu.getPosition().z); + float TrackEndPDis = (TrkEndPoint - CluPos).Mag(); + if( (TrackEndPDis < 400 + 400./TrackEn || (TrackEn < 1 && candiclu.getEnergy() < 2 )) && ( fabs(TrackEn - candiclu.getEnergy() ) < 3.0*sqrt(TrackEn) + 1.0 || candiclu.getEnergy() < 3) && (TrkEndPoint.Z()*CluPos.Z() > 0 || abs(TrkEndPoint.Z() - CluPos.Z()) < 100 ) ) // && AngDiff < 0.4 ) + { + Closest_Trk_Clu_Map[g0] = ClosestCluIndex; + BushTouchFlag[candiclu] = 0; + } + } + } //~~~~~~~ end of finding closest cluster + + cout<<"NTrk "<<NTrk<<endl; + for(int i0 = 0; i0 < NTrk; i0++) //Dropped Size can exist + { + auto a_trk = SortedTracks[i0]; + currTrackType = Track_Type[a_trk]; + currTrkEn = Track_Energy[a_trk]; + + TrkEndPoint = Track_EndPoint[a_trk]; + FCMap_Track_CHCore_new[a_trk].clear(); + if( Closest_Trk_Clu_Map.find(i0) != Closest_Trk_Clu_Map.end() ) + { + auto closeClu = selfmergedcluster[Closest_Trk_Clu_Map[i0]]; + FCMap_Track_CHCore_new[a_trk].push_back(closeClu); + } + + for(int j0 = 0; j0 < NClu; 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]; + CluPos = TVector3(fccand_bush.getPosition().x,fccand_bush.getPosition().y,fccand_bush.getPosition().z); + int currCluType = m_ArborToolLCIO->ClusterFlag1st(fccand_bush); + + CoreMergeDistanceDepthCorrector = 0; + if(CluDepth > 20) + CoreMergeDistanceDepthCorrector = 20; + else if(CluDepth > 10) + CoreMergeDistanceDepthCorrector = 10; + + float ProjectiveDis_TrackEndP_CluPos = (CluPos - TrkEndPoint).Mag()*(CluPos - TrkEndPoint).Angle(CluPos); + if(log10(BushTime) < 3.5 && BushTime > 0 && currTrackType != 101 && Dis < 10 + CoreMergeDistanceDepthCorrector && Dis > -0.1 && BushTouchFlag.find(fccand_bush) == BushTouchFlag.end() && (currTrkEn > 3 || fccand_bush.getEnergy() < 5 || currCluType == 13 ) && ProjectiveDis_TrackEndP_CluPos < 100 ) // && (TrackEndPDis < 400 || currTrackType != 11) && (CluPos.Z()*TrkEndPoint.Z() > 0 || fabs((CluPos-TrkEndPoint).Z()) < 20 ) && (fccanden - currTrkEn < something) ) + { + FCMap_Track_CHCore_new[a_trk].push_back(fccand_bush); + BushTouchFlag[fccand_bush] = currTrackType; + + } + } + //Diag 4 MIMIC + } + + //Added by Dan for multi track linked to one cluster + for(int i0 = 0; i0 < NTrk; i0++) //Dropped Size can exist + { + auto a_trk = SortedTracks[i0]; + TrkEndPoint = Track_EndPoint[a_trk]; + int naclu=FCMap_Track_CHCore_new[a_trk].size(); + + + for(int i1=0;i1<naclu;i1++) + { + auto aclu=FCMap_Track_CHCore_new[a_trk][i1]; + float dis_a=DisMatrix_Track_Clu_E[i0][i1]; + for(int j0=i0+1; j0<NTrk; j0++) + { + int breakflag=0; + auto b_trk = SortedTracks[j0]; + int nbclu=FCMap_Track_CHCore_new[b_trk].size(); + + for(int j1=0;j1<nbclu;j1++) + { + auto bclu=FCMap_Track_CHCore_new[b_trk][j1]; + float dis_b=DisMatrix_Track_Clu_E[j0][j1]; + if(aclu==bclu){ + if(dis_a>=dis_b){ + FCMap_Track_CHCore_new[a_trk].erase(FCMap_Track_CHCore_new[a_trk].begin()+i1); + breakflag=1; + break; + } + else{ + FCMap_Track_CHCore_new[b_trk].erase(FCMap_Track_CHCore_new[b_trk].begin()+j1); + } + } + } + if(breakflag)break; + } + } + if(FCMap_Track_CHCore_new[a_trk].size() > 0 ) // && EcalCoreEnergy + HcalCoreEnergy < 2.0*currTrkEn )... + { + auto chcorecluster_eh = m_ArborToolLCIO->NaiveMergeClu(FCMap_Track_CHCore_new[a_trk]); + edm4hep::ConstCluster chcorecluster_ehCon=chcorecluster_eh; + FCMap_Track_CHCore[a_trk] = chcorecluster_ehCon; + chargedclustercore.push_back(chcorecluster_ehCon); + } + }// End by Dan + + // Diagnosis in Tagcore + for(int t0 = 0; t0 < NClu; t0++) + { + auto a_clu = selfmergedcluster[t0]; + if( BushTouchFlag.find(a_clu) == BushTouchFlag.end() ) //Might be a neutral core + { + non_chargedclustercore.push_back(a_clu); + } + } + + int Track_Core_ID = -99; + + edm4hep::ReconstructedParticleCollection* chargeparticleCol = m_chargeparticleCol.createAndPut(); + edm4hep::ClusterCollection* chargedcoreclusterCol = m_chargedcoreclusterCol.createAndPut(); + for(int j5 = 0; j5 < NTrk; j5++) + { + auto a_trk = SortedTracks[j5]; + Track_Core_ID = -99; + currTrackType = Track_Type[a_trk]; + auto chargeparticle=chargeparticleCol->create(); + chargeparticle.setEnergy( Track_Energy[a_trk] ); + 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.addToTracks( a_trk ); + auto a_clu = FCMap_Track_CHCore[a_trk]; + if( FCMap_Track_CHCore[a_trk].hits_size()>0 ) // No really need to pertect, as quality will be controled in Part.Reco + { + edm4hep::Cluster chargedcorecluster = chargedcoreclusterCol->create(); + + m_ArborToolLCIO->NaiveCluConst(FCMap_Track_CHCore[a_trk],chargedcorecluster); + edm4hep::ConstCluster chargedcoreclusterCon=chargedcorecluster; + chargeparticle.addToClusters(chargedcoreclusterCon); + Track_Core_ID = m_ArborToolLCIO->ClusterFlag(a_clu, a_trk); + } + chargeparticle.setType(Track_Core_ID); + ChCoreID[chargeparticle] = Track_Core_ID; + } + + +} + +void BushConnect::ParticleReco() +{ + + auto col_IsoHit = m_col_IsoHit.get(); + std::vector<edm4hep::ConstCalorimeterHit> IsoHits = m_ArborToolLCIO->CollHitVec(col_IsoHit, 0); + + edm4hep::ReconstructedParticleCollection* arborrecoparticleCol = m_arborrecoparticleCol.createAndPut(); + + edm4hep::ClusterCollection* mergedclu_chCol = m_mergedclu_chCol.createAndPut(); + edm4hep::ClusterCollection* mergedclu_neCol = m_mergedclu_neCol.createAndPut(); + + 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; + + float CluDepth = 0; + std::map<edm4hep::ConstCluster, double> CluDepthMap; + CluDepthMap.clear(); + int currChargeCoreType = 0; + TVector3 CluPos; + + std::vector<edm4hep::ConstCluster> loosecandicluster; + std::vector<edm4hep::ConstCluster> tightcandicluster; //Muon potential candi? + std::vector<edm4hep::ConstCluster> mergedcluster; //tmp for each charged P + std::vector<edm4hep::ConstCluster> chargedclustercore_merged; //overall + chargedclustercore_merged.clear(); + + std::vector<double> reftightdis; + std::vector<double> refloosedis; + std::map<edm4hep::ConstCluster, int> NNCTouchFlag; + std::vector<edm4hep::ConstTrack> SecondIterTracks; + SecondIterTracks.clear(); + + TVector3 currTrkEnd, neighbourTrkEnd, LeadP; + + for(int i = 0; i < NChargedObj; i++) + { + auto a_recoP_ch = (*ChargedCore)[i]; + + loosecandicluster.clear(); + tightcandicluster.clear(); + mergedcluster.clear(); + reftightdis.clear(); + refloosedis.clear(); + auto a_chargedTrk = a_recoP_ch.getTracks(0); + currTrkEnd = Track_EndPoint[a_chargedTrk]; + currChargeCoreType = ChCoreID[a_recoP_ch]; + int currTrkType = Track_Type[a_chargedTrk]; + + float CurrClusterEnergy = 0; + float CurrTrackEnergy = Track_Energy[a_chargedTrk]; + edm4hep::ConstCluster a_chargedClu; + if(a_recoP_ch.clusters_size() != 0) + { + a_chargedClu = a_recoP_ch.getClusters(0); + CurrClusterEnergy = a_chargedClu.getEnergy(); + mergedcluster.push_back(a_chargedClu); //Actually can use this chance to question if previous energy are balance... + } + + float MinDisToNoClusterTrk = 1.0E10; + float MinDisToOtherTrack = 1.0E10; + + for( int is = 0; is < NChargedObj; is++ ) + { + if(is != i) + { + auto b_recoP_ch = (*ChargedCore)[is]; + auto a_neighbourTrk = b_recoP_ch.getTracks(0); + neighbourTrkEnd = Track_EndPoint[a_neighbourTrk]; + float currDD = (neighbourTrkEnd - currTrkEnd).Mag(); + if( currDD < MinDisToOtherTrack ) + { + MinDisToOtherTrack = currDD; + } + } + } + + for(int j = 0; j < NNeutralCluster; j++) + { + auto a_NeCandiClu = non_chargedclustercore[j]; + 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; + + if( ClusterType_1stID[a_NeCandiClu] == 22) continue; + + for(int k = 0; k < 2; k++) + { + DisMatrix_Core_Neutral[i][j][k] = 1.0E9; + } + + if(CurrClusterEnergy > 1E-6) //put by hand... + { + DisMatrix_Core_Neutral[i][j][0] = m_ArborToolLCIO->BushDis(a_chargedClu, a_NeCandiClu); + } + float* Dis = m_ArborToolLCIO->SimpleDisTrackClu(a_chargedTrk, a_NeCandiClu); + 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? + { + tightcandicluster.push_back(a_NeCandiClu); + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][1] ); //dependence on Cluster Flag & Clu Depth. use some more fancy sort algorithm... + } + } + else if( currChargeCoreType == 131 ) + { + 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] ); + } + 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 ) + { + 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 + { + tightcandicluster.push_back(a_NeCandiClu); + if(DisMatrix_Core_Neutral[i][j][0] < DisMatrix_Core_Neutral[i][j][1]) // not fully adequate. + { + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][0] ); + } + else + { + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][1] ); + } + } + 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. + { + refloosedis.push_back( DisMatrix_Core_Neutral[i][j][0] ); + } + else + { + refloosedis.push_back( DisMatrix_Core_Neutral[i][j][1] ); + } + } + } + else if( currChargeCoreType == 211 ) //Main Cluster distance oriented + { + if(DisMatrix_Core_Neutral[i][j][0] < 0.2*CluDepth) + { + tightcandicluster.push_back(a_NeCandiClu); + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][0] ); + } + } + else if( currChargeCoreType == 212 ) //Non_Matched + { + if(DisMatrix_Core_Neutral[i][j][0] < 10 + 0.5*CluDepth ) //Energy Dependence... + { + tightcandicluster.push_back(a_NeCandiClu); + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][0] ); + } + else if(DisMatrix_Core_Neutral[i][j][0] < 10 + 0.4*CluDepth || DisMatrix_Core_Neutral[i][j][1] < 20 + 0.5*CluDepth ) + { + loosecandicluster.push_back(a_NeCandiClu); + refloosedis.push_back( DisMatrix_Core_Neutral[i][j][0] ); + } + } + else if( currChargeCoreType == 0 ) // && a_recoP_ch->getEnergy() < 3 ) // && !FlagLowPtPz ) // + { + if(CluDepth < 20) + { + if(DisMatrix_Core_Neutral[i][j][1] < MinDisToNoClusterTrk) //Tag minimal distance cluster... and see if it can be potentially linked. + { + MinDisToNoClusterTrk = DisMatrix_Core_Neutral[i][j][1]; + } + if( MinDisToNoClusterTrk < 300 && abs(a_recoP_ch.getEnergy() - NeCandEn) < 1.5*a_recoP_ch.getEnergy() ) //some hard cut + { + tightcandicluster.push_back(a_NeCandiClu); + reftightdis.push_back( DisMatrix_Core_Neutral[i][j][1] ); + } + } + } + else + { + cout<<"Over balanced/Un matched/defined case: "<<a_recoP_ch.getEnergy()<<" ??? "<<currChargeCoreType<<endl; + } + } + } + + float totaltightcandiEn = 0; + float totalloosecandiEn = 0; + for(unsigned int s = 0; s < tightcandicluster.size(); s++) + { + totaltightcandiEn += tightcandicluster[s].getEnergy(); + } + + for(unsigned int s = 0; s < loosecandicluster.size(); s++) + { + totalloosecandiEn += loosecandicluster[s].getEnergy(); + } + + if( currChargeCoreType == 130 ) + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.0*sqrt(CurrTrackEnergy) ) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + else if( currChargeCoreType == 131 ) + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.0*sqrt(CurrTrackEnergy)) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + + //Maybe Some ID over here?... //layers & numbers... //BS Case ID + + for(unsigned int i2 = 0; i2 < loosecandicluster.size(); i2++) + { + auto a_clu = loosecandicluster[i2]; + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.0*sqrt(CurrTrackEnergy)) //Frags...Or some minmal hit cut + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + else if( currChargeCoreType == 110 ) + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + else if( currChargeCoreType == 111 ) + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy)) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + + for(unsigned int i2 = 0; i2 < loosecandicluster.size(); i2++) + { + auto a_clu = loosecandicluster[i2]; + if( fabs(CurrClusterEnergy + a_clu.getEnergy() - CurrTrackEnergy) < fabs(CurrClusterEnergy - CurrTrackEnergy) ) //Frags...Or some minmal hit cut + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + else if( currChargeCoreType == 211 ) // Matched + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + else if( currChargeCoreType == 212) + { + for(unsigned int i1 = 0; i1 < tightcandicluster.size(); i1++) + { + auto a_clu = tightcandicluster[i1]; + if( CurrClusterEnergy + a_clu.getEnergy() < CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy)) + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + + for(unsigned int i2 = 0; i2 < loosecandicluster.size(); i2++) + { + auto a_clu = loosecandicluster[i2]; + if( fabs(CurrClusterEnergy + a_clu.getEnergy() - CurrTrackEnergy) < fabs(CurrClusterEnergy - CurrTrackEnergy) ) //Frags...Or some minmal hit cut + { + mergedcluster.push_back( a_clu ); + CurrClusterEnergy += a_clu.getEnergy(); + } + } + } + + + else if( currChargeCoreType == 0 && reftightdis.size() > 0) + { + float mindis = 1.0E10; + int minindex = 0; + + for(unsigned int i1 = 0; i1 < reftightdis.size(); i1 ++) + { + if(reftightdis[i1] < mindis) + { + mindis = reftightdis[i1]; + minindex = i1; + } + } + + auto a_clu = tightcandicluster[minindex]; // Only 1? ... + + mergedcluster.push_back( a_clu ); + } + else + { + 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(); + } + + if( !( CHCluEnergy < 1 && CurrTrackEnergy > 5 ) || (MinDisToOtherTrack < 100) ) // Need to check if exist nearby larger charged cluster: maybe absorbed together and only left tiny MIP tail in the ECAL //* bool closeTonearByEnergeticChargedShower = 0 // MIP like; should also protect against energies + { + for(int i2 = 0; i2 < int(mergedcluster.size()); i2++) + { + auto a_TBM_clu = mergedcluster[i2]; + NNCTouchFlag[a_TBM_clu] = 2; // can make use of this intereting flag... + } + + float charge = a_chargedTrk.getTrackStates(0).omega/fabs(a_chargedTrk.getTrackStates(0).omega); + + auto chargeparticle = arborrecoparticleCol->create(); + chargeparticle.setCharge(charge); + chargeparticle.addToTracks( a_chargedTrk ); + TVector3 Ptrack = Track_P3[a_chargedTrk]; + double currTrkP[3] = { Ptrack.X(), Ptrack.Y(), Ptrack.Z() }; + + int flagEnergyFlow = 0; + int currChargeCoreType2 = -99; + + + auto chclustermerged = mergedclu_chCol->create(); + m_ArborToolLCIO->NaiveMergeCluConst(mergedcluster,chclustermerged); + edm4hep::ConstCluster chclustermergedCon=chclustermerged; + chargeparticle.addToClusters(chclustermergedCon); + chargedclustercore_merged.push_back(chclustermerged); + currChargeCoreType2 = m_ArborToolLCIO->ClusterFlag(chclustermerged, a_chargedTrk); + + if( currChargeCoreType2 == 130 || currChargeCoreType2 == 131 ) + { + chargeparticle.setType( int(-13*charge) ); + } + else if( currChargeCoreType2 == 110 || currChargeCoreType2 == 111 ) + { + chargeparticle.setType( int(-11*charge) ); + if(CHCluEnergy > CurrTrackEnergy + 0.5*sqrt(CurrTrackEnergy) + 1) + { + flagEnergyFlow = 1; + } + } + else + { + chargeparticle.setType( Track_Type[a_chargedTrk] ); + if(CHCluEnergy > CurrTrackEnergy + 1.5*sqrt(CurrTrackEnergy) + 1) + { + flagEnergyFlow = 2; + } + } + + 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] ); + chargeparticle.setMomentum( CorrMom ); + } + else + { + chargeparticle.setEnergy( CurrTrackEnergy ); + chargeparticle.setMomentum(edm4hep::Vector3f( currTrkP[0],currTrkP[1],currTrkP[2]) ); + } + + if( flagEnergyFlow ) + { + 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(); + 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<<"a charged particle reconstructed with en:"<<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; + } + } + + std::vector<edm4hep::ConstCluster> 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); + } + } + + float NAMom[3] = {0, 0, 0}; + + //Final Re-absorption + std::vector<edm4hep::ConstCluster> 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; + Depth = DisSeedSurface(PosClu); + float CoreEnCorr = m_ArborToolLCIO->ClusterEE(a_clu); + + if(m_ArborToolLCIO->newPhotonTag(a_clu)==1) + { + TVector3 BushSeedPos = TVector3(a_clu.getPosition().x,a_clu.getPosition().y,a_clu.getPosition().z); + auto neutralparticle = arborrecoparticleCol->create(); + neutralparticle.setType(22); + TVector3 PP = m_ArborToolLCIO->ClusterCoG(a_clu); + NAMom[0] = CoreEnCorr*1.0/PP.Mag()*PP.X(); + NAMom[1] = CoreEnCorr*1.0/PP.Mag()*PP.Y(); + NAMom[2] = CoreEnCorr*1.0/PP.Mag()*PP.Z(); + neutralparticle.setEnergy( CoreEnCorr ); + neutralparticle.setMass( 0.0 ); + neutralparticle.setCharge( 0.0 ); + neutralparticle.setMomentum( edm4hep::Vector3f(NAMom[0],NAMom[1],NAMom[2]) ); + auto a_neclu = mergedclu_neCol->create(); + m_ArborToolLCIO->NaiveCluConst(a_clu,a_neclu); + a_neclu.setEnergy( CoreEnCorr ); //Reset... + edm4hep::ConstCluster a_necluCon=a_neclu; + neutralparticle.addToClusters(a_neclu); + } + else // Distance to Charged Core > sth; + { + float MinDisToChCore = 1.0E9; + float currDis = 0; + int NChCore = mergedclu_chCol->size(); + float closestChCluEn = 0; + for(int t = 0; t < NChCore; t++) + { + auto a_chclu = (*mergedclu_chCol)[t]; + currDis = m_ArborToolLCIO->BushDis(a_chclu, a_clu); + if(currDis < MinDisToChCore) + { + MinDisToChCore = currDis; + closestChCluEn = a_chclu.getEnergy(); // Or the Trk En?? + } + } + if( MinDisToChCore > 0.4*(15 + closestChCluEn + Depth*0.01) || a_clu.getEnergy() > 2.0 ) //Joint Depth?? + { + NBBNeutral.push_back(a_clu); + } + } + } + + // Add: Neural Core Remerge & Energy Scale Recalculate, IsoHit Abso + std::vector<edm4hep::Cluster> NBBAbs = m_ArborToolLCIO->ClusterHitAbsorbtion(NBBNeutral, IsoHits, 100); //_HitAbsCut); // Huge?? + + std::vector<float> BBAbsEn; + BBAbsEn.clear(); + + for(unsigned s1 = 0; s1 < NBBAbs.size(); s1++) + { + BBAbsEn.push_back(NBBAbs[s1].getEnergy()); + } + + std::vector<int> BBAbsIndex = SortMeasure(BBAbsEn, 1); + + std::vector<edm4hep::ConstCluster > NeutronCore; + std::vector<edm4hep::ConstCluster > NeutronFlag; + NeutronCore.clear(); + NeutronFlag.clear(); + + for(unsigned int s2 = 0; s2 < NBBAbs.size(); s2++) //Sort it; the first one must be a neutral core? + { + auto p_clu = NBBAbs[BBAbsIndex[s2]]; + float currCluEn = p_clu.getEnergy(); + std::vector<float> CluTime = m_ArborToolLCIO->ClusterTime(p_clu); + if( (currCluEn > 1.0 || (currCluEn > 0.5 && s2 < 2) )&& CluTime[0] < 10) + { + NeutronCore.push_back(p_clu); + } + else + { + NeutronFlag.push_back(p_clu); + } + } + + std::vector<edm4hep::ConstCluster > Neutrons = m_ArborToolLCIO->ClusterAbsorbtion(NeutronCore, NeutronFlag, 200, 0.01); + + for(unsigned int s3 = 0; s3 < Neutrons.size(); s3++) + { + auto a_clu = Neutrons[s3]; + float CoreEnCorr = m_ArborToolLCIO->ClusterEE(a_clu); + 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; + int NChCore = mergedclu_chCol->size(); + for(int t = 0; t < NChCore; t++) + { + auto a_chclu=(*mergedclu_chCol)[t]; + currDis = m_ArborToolLCIO->BushDis(a_chclu, a_clu); + if(currDis < MinDisToChCore) + { + MinDisToChCore = currDis; + } + } + + if( !(RecoT0>0.1 && RecoT0<1E8 && MinDisToChCore <12) ) + { + if(m_ArborToolLCIO->newPhotonTag(a_clu)==1) + cout<<"WARNING... Photons after neutron merge merged"<<endl; + auto neutralparticle = arborrecoparticleCol->create(); + neutralparticle.setType(21120); + TVector3 PP = m_ArborToolLCIO->ClusterCoG(a_clu); + + if( RecoT0 > 0.16 && RecoT0 < 100 && MinDisToChCore > 50 && a_clu.hits_size() > 5 && abs(CoreEnCorr - 6) < 4 ) + { + neutralparticle.setType(2112); + float Dis = PP.Mag(); + float Beta = 1.0/(1 + 300*RecoT0/Dis); + float PPN = 0.941/sqrt(1-Beta*Beta); + float PPK = 0.53*PPN; + if(PPN/CoreEnCorr < 1) + { + CoreEnCorr = PPN; + neutralparticle.setType(2112); + } + else if(PPK/CoreEnCorr < 1 && PPN/CoreEnCorr > 1) + { + CoreEnCorr = 0.5*(PPN+PPK); + neutralparticle.setType(2112); + } + else if(PPK/CoreEnCorr > 1) + { + CoreEnCorr = PPK; + neutralparticle.setType(310); + } + } + + NAMom[0] = CoreEnCorr*1.0/PP.Mag()*PP.X(); + NAMom[1] = CoreEnCorr*1.0/PP.Mag()*PP.Y(); + NAMom[2] = CoreEnCorr*1.0/PP.Mag()*PP.Z(); + neutralparticle.setEnergy( CoreEnCorr ); + neutralparticle.setMass( 0.0 ); + neutralparticle.setCharge( 0.0 ); + neutralparticle.setMomentum( edm4hep::Vector3f(NAMom[0],NAMom[1],NAMom[2]) ); + auto a_neclu = mergedclu_neCol->create(); + m_ArborToolLCIO->NaiveCluConst(a_clu,a_neclu); + a_neclu.setEnergy( CoreEnCorr ); //Reset... + edm4hep::ConstCluster a_necluCon=a_neclu; + neutralparticle.addToClusters(a_necluCon); + } + } + + +} + +StatusCode BushConnect::execute() +{ + try{ + BushConnect::Clean(); + BushConnect::TrackSort( ); + BushConnect::BushSelfMerge( ); + BushConnect::TagCore( ); + BushConnect::ParticleReco( ); + }catch(GaudiException &e){} + return StatusCode::SUCCESS; +} + +StatusCode BushConnect::finalize() +{ + 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 new file mode 100644 index 0000000000000000000000000000000000000000..e0f23ba8ca625d4ef32d9ad62ca3c78b3e260e05 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/BushConnect.hh @@ -0,0 +1,124 @@ +#ifndef _BushConnect_hh_ +#define _BushConnect_hh_ + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/Cluster.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +#include "Arbor.h" +#include "ArborToolLCIO.hh" + +#include <TTree.h> +#include <TFile.h> +#include <TH1.h> +#include <TH2.h> +#include <TH3.h> + + +class BushConnect : public GaudiAlgorithm +{ + public: + + BushConnect(const std::string& name, ISvcLocator* svcLoc); + + + virtual StatusCode initialize() ; + + virtual StatusCode execute() ; + + virtual StatusCode finalize() ; + void Clean(); + void TrackSort(); + void BushSelfMerge(); + void TagCore(); + void ParticleReco(); + + + protected: + + std::vector<edm4hep::ConstCluster> SortedSMBushes; + std::vector<edm4hep::ConstTrack> SortedTracks; + std::map<edm4hep::ConstTrack, float> Track_Energy; + std::map<edm4hep::ConstTrack, TVector3> Track_P3; + std::map<edm4hep::ConstTrack, int> Track_Type; + std::map<edm4hep::ConstTrack, float> Track_Theta; + std::map<edm4hep::ConstTrack, float> Track_Phi; + + std::map<edm4hep::ConstCluster, int> ClusterType_1stID; + std::map<edm4hep::ReconstructedParticle, int> ChCoreID; + + std::vector<edm4hep::ConstCluster> ecalchcore_tight; //TightCores + std::vector<edm4hep::ConstCluster> ecalchcore_medium; + std::vector<edm4hep::ConstCluster> ecalchcore_loose; //LooseCores Let's also get + std::vector<edm4hep::ConstCluster> ecalchcore; //Above three + std::vector<edm4hep::ConstCluster> ecalnecore; + std::vector<edm4hep::ConstCluster> ecalnecore_EM; + std::vector<edm4hep::ConstCluster> ecalnecore_NonEM; + std::vector<edm4hep::ConstCluster> ecalfrag; + std::vector<edm4hep::ConstCluster> ecalundef; + std::vector<edm4hep::ConstCluster> ecalfrag_TBM_CH; + std::vector<edm4hep::ConstCluster> ecalfrag_TBM_NE_EM; + std::vector<edm4hep::ConstCluster> ecalfrag_TBM_NE_NonEM; + std::vector<edm4hep::ConstCluster> ecalundef_iso; + std::vector<edm4hep::ConstCluster> ecalpotentialbackscattering; + + std::vector<edm4hep::ConstCluster> chargedclustercore; + std::vector<edm4hep::ConstCluster> chargedclustercore_abs; + + std::vector<edm4hep::ConstCluster> selfmergedcluster; + std::vector<edm4hep::ConstCluster> non_chargedclustercore; + std::vector<edm4hep::ConstCluster> onlyNeutralCore; + + std::vector<edm4hep::ConstCluster> non_charged_pem_neutral_core; + std::vector<edm4hep::ConstCluster> pem_neutral_core; + + std::map<edm4hep::ConstTrack, int>MCPTrack_Type; + std::map<edm4hep::ConstTrack, TVector3> Track_EndPoint; //Last hit + std::map<edm4hep::ConstTrack, TVector3> TrackStartPoint; + std::map<edm4hep::ConstCluster, float> CluFD; + std::map<edm4hep::ConstCluster, float> CluT0; + std::map<edm4hep::ConstCluster, float> Clu_Depth; + std::map<edm4hep::ConstCluster, TVector3> CluCoG; + typedef DataHandle<edm4hep::MCParticleCollection> MCParticleColHandler; + MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this}; + + typedef DataHandle<edm4hep::TrackCollection> TrkColHandler; + TrkColHandler m_trkcol {"MarlinTrkTracks", Gaudi::DataHandle::Reader, this}; + + typedef DataHandle<edm4hep::ClusterCollection> CluColHandler; + CluColHandler m_clucol {"EHBushes", Gaudi::DataHandle::Reader, this}; + + DataHandle<edm4hep::CalorimeterHitCollection> m_col_IsoHit{"IsoHits",Gaudi::DataHandle::Reader, this}; + + DataHandle<edm4hep::ClusterCollection> m_1stclucol{"1stAbs",Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ClusterCollection> clucol_merged{"mergedCluC",Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ClusterCollection> m_mergedclu_chCol{"mergedCluCh",Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ClusterCollection> m_mergedclu_neCol{"mergedCluNe",Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ClusterCollection> m_chargedcoreclusterCol{"ChargedClu",Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ClusterCollection> NetralClu{"NetralClu",Gaudi::DataHandle::Writer, this}; + 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; + +}; + + +#endif + + diff --git a/Reconstruction/PFA/Arbor/src/ClusterAna.cc b/Reconstruction/PFA/Arbor/src/ClusterAna.cc new file mode 100644 index 0000000000000000000000000000000000000000..bce1ccd1e7cbd97a6f23b8e713c5fb10cd403aca --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ClusterAna.cc @@ -0,0 +1,185 @@ +#include "ClusterAna.hh" +#include <EVENT/LCCollection.h> +#include <IMPL/LCCollectionVec.h> +#include <EVENT/LCFloatVec.h> +#include <EVENT/MCParticle.h> +#include <EVENT/ReconstructedParticle.h> +#include <IMPL/MCParticleImpl.h> +#include <values.h> +#include <string> +#include <iostream> +#include <EVENT/LCFloatVec.h> +#include <EVENT/LCParameters.h> +#include <stdexcept> +#include <TFile.h> +#include <TTree.h> +#include <TH1F.h> +#include <TVector3.h> +#include <TRandom.h> +#include <Rtypes.h> +#include <sstream> +#include <cmath> +#include <vector> +#include <TMath.h> +#include "TLorentzVector.h" +#include <UTIL/CellIDDecoder.h> + +using namespace std; + +DECLARE_COMPONENT(ClusterAna) + + +const float sqrts = 240.0; //GeV + +const string ECALCellIDDecoder = "M:3,S-1:3,I:9,J:9,K-1:6"; +//TH1F *h_hit; + +ClusterAna::ClusterAna(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc), + _output(0) +{ +} + +StatusCode ClusterAna::initialize() { + + // printParameters(); + + TFile *tree_file=new TFile(_treeFileName.value().c_str(),(_overwrite ? "RECREATE" : "UPDATE")); + + if (!tree_file->IsOpen()) { + delete tree_file; + tree_file=new TFile(_treeFileName.value().c_str(),"NEW"); + } + + //h_hit=new TH1F("hit","hit",80,0,80); + _outputTree = new TTree(_treeName.value().c_str(),_treeName.value().c_str()); + _outputTree->SetAutoSave(32*1024*1024); // autosave every 32MB + _outputTree->Branch("EventNr", &_eventNr, "EventNr/I"); + _outputTree->Branch("Num", &_Num, "Num/I"); + _outputTree->Branch("NClu", &_NClu, "NClu/I"); + _outputTree->Branch("NPFO", &_NPFO, "NPFO/I"); + _outputTree->Branch("EcalTotalE", &_EcalTotalE, "EcalTotalE/F"); + _outputTree->Branch("HcalTotalE", &_HcalTotalE, "HcalTotalE/F"); + _outputTree->Branch("CluE", &_CluE, "CluE/F"); + + _outputClu = new TTree("Clu","Clu"); + _outputClu->Branch("EventNr", &_eventNr, "EventNr/I"); + _outputClu->Branch("Num", &_Num, "Num/I"); + _outputClu->Branch("Nhit", &_Nhit, "Nhit/I"); + _outputClu->Branch("CluEn", &_CluEn, "CluEn/F"); + _outputClu->Branch("CluPos", CluPos, "CluPos[3]/F"); + + _outputPFO = new TTree("PFO","PFO"); + _outputPFO->Branch("Num", &_Num, "Num/I"); + _outputPFO->Branch("PID", &_PID, "PID/I"); + _outputPFO->Branch("PFOEn", &_PFOEn, "PFOEn/F"); + + + _Num = 0; + + return GaudiAlgorithm::initialize(); + +} + +StatusCode ClusterAna::execute() +{ + + EVENT::LCEvent* evtP = nullptr; + + std::vector<CaloHitColHandler*> hdl_EcalHitColl{ + &m_ecalbarrelhitcol, + &m_ecalendcaphitcol + }; + std::vector<CaloHitColHandler*> hdl_HcalHitColl{ + &m_hcalbarrelhitcol, + &m_hcalendcaphitcol, + &m_hcalotherhitcol + }; + + std::vector<std::string> EcalHitColl; + std::vector<std::string> HcalHitColl; + EcalHitColl.push_back("ECALBarrel"); + EcalHitColl.push_back("ECALEndcap"); + HcalHitColl.push_back("HCALBarrel"); + HcalHitColl.push_back("HCALEndcap"); + HcalHitColl.push_back("HCALOther"); + _CluE=0; + + try{ + + for(int t = 0; t< int(hdl_EcalHitColl.size()); t++) { + const edm4hep::CalorimeterHitCollection* ecalcoll = hdl_EcalHitColl[t]->get(); + for(auto hit: *ecalcoll) { + // TODO + int NLayer = 0; + _EcalTotalE += hit.getEnergy(); + + } + } + +// for(int t2 = 0; t2< int(hdl_HcalHitColl.size()); t2++) { +// const edm4hep::CalorimeterHitCollection* hcalcoll = hdl_HcalHitColl[t2]->get(); +// for (auto hit: *hcalcoll) { +// // TODO +// int NLayer = 0; +// int HLayer = NLayer+30; +// // UTIL::CellIDDecoder<EVENT::CalorimeterHit> idDecoder(ECALCellIDDecoder); +// // int NLayer = idDecoder(a_hit)["K-1"]; +// +// //h_hit->Fill(HLayer,a_hit->getEnergy()); +// _HcalTotalE += hit.getEnergy(); +// +// } +// } +// +// + }catch(lcio::DataNotAvailableException err) { } + + try{ + auto col_Clu = m_Clu.get(); + _NClu=col_Clu->size(); + for(int i0 = 0; i0 < col_Clu->size(); i0++) { + auto a_clu = (*col_Clu)[i0]; + auto currPos0 = a_clu.getPosition(); + TVector3 currPos(currPos0.x, currPos0.y, currPos0.z); + _CluEn=a_clu.getEnergy(); + _CluE+=a_clu.getEnergy(); + + + _outputClu->Fill(); + + } + }catch (lcio::DataNotAvailableException err) { } + + try{ + auto col_PFO = m_PFO.get(); + _NPFO=col_PFO->size(); + for(int i0 = 0; i0 < col_PFO->size(); i0++){ + auto a_PFO=(*col_PFO)[i0]; + _PFOEn=a_PFO.getEnergy(); + _PID=a_PFO.getType(); + _outputPFO->Fill(); + } + }catch (lcio::DataNotAvailableException err) { } + _outputTree->Fill(); + _Num++; + // } + return StatusCode::SUCCESS; + +} + +StatusCode ClusterAna::finalize() +{ + + if (_outputTree) { + + TFile *tree_file = _outputTree->GetCurrentFile(); //just in case we switched to a new file + tree_file->Write(); + delete tree_file; + } + + return GaudiAlgorithm::finalize(); +} + + + diff --git a/Reconstruction/PFA/Arbor/src/ClusterAna.hh b/Reconstruction/PFA/Arbor/src/ClusterAna.hh new file mode 100644 index 0000000000000000000000000000000000000000..48ab40a1313d67bf09bc67d0893449c41f1d1bd3 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/ClusterAna.hh @@ -0,0 +1,79 @@ +#ifndef _ClusterAna_hh_ +#define _ClusterAna_hh_ + +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" + +#include "k4FWCore/DataHandle.h" + +#include <string> +#include <iostream> +#include <fstream> +#include <TNtuple.h> +#include <TObject.h> +#include <TTree.h> +#include <TFile.h> + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/ReconstructedParticleCollection.h" + +class ClusterAna : public GaudiAlgorithm +{ +public: + + ClusterAna(const std::string& name, ISvcLocator* svcLoc); + + ~ClusterAna() {}; + + StatusCode initialize() override; + + StatusCode execute() override; + + StatusCode finalize() override; + +protected: + typedef DataHandle<edm4hep::MCParticleCollection> MCParticleColHandler; + MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this}; + + typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloHitColHandler; + CaloHitColHandler m_ecalbarrelhitcol{"ECALBarrel", Gaudi::DataHandle::Reader, this}; + CaloHitColHandler m_ecalendcaphitcol{"ECALEndcap", Gaudi::DataHandle::Reader, this}; + + CaloHitColHandler m_hcalbarrelhitcol{"HCALBarrel", Gaudi::DataHandle::Reader, this}; + CaloHitColHandler m_hcalendcaphitcol{"HCALEndcap", Gaudi::DataHandle::Reader, this}; + CaloHitColHandler m_hcalotherhitcol {"HCALOther", Gaudi::DataHandle::Reader, this}; + + typedef DataHandle<edm4hep::ClusterCollection> CluColHandler; + CluColHandler m_Clu{"EHBushes", Gaudi::DataHandle::Reader, this}; + + typedef DataHandle<edm4hep::ReconstructedParticleCollection> PFOColHandler; + PFOColHandler m_PFO{"ArborPFO", Gaudi::DataHandle::Reader, this}; + +Gaudi::Property<std::string> _treeFileName{this, + "TreeOutputFile", "Ana.root", + "The name of the file to which the ROOT tree will be written"}; + Gaudi::Property<std::string> _treeName{this, + "TreeName", "Evt", + "The name of the ROOT tree"}; + std::string _colName; + std::string _colAdcVals; + + Gaudi::Property<int> _overwrite{this, + "OverwriteFile", 0, + "If zero an already existing file will not be overwritten."}; + TTree *_outputTree, *_outputClu, *_outputPFO; + + float _EcalTotalE,_HcalTotalE,_CluE,_CluEn,_PFOEn; + int _eventNr,_Num,_NClu,_Nhit,_NPFO,_PID; + float CluPos[3]; + + std::string _fileName; + std::ostream *_output; + std::string _histFileName; +}; + +#endif + + diff --git a/Reconstruction/PFA/Arbor/src/DetectorPos.cc b/Reconstruction/PFA/Arbor/src/DetectorPos.cc new file mode 100644 index 0000000000000000000000000000000000000000..cc70613b44bc445fa51e4e885d57fad0224a3106 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/DetectorPos.cc @@ -0,0 +1,191 @@ +#include <TMath.h> +#include "DetectorPos.hh" + +using namespace std; + +const double pi = acos(-1.0); + +//Geometric Parameter - ... need to be changed for different detector models + +int BarrelFlag( TVector3 inputPos) +{ + int isBarrel = 0; + + if(fabs(inputPos[2]) < ECALHalfZ) + { + isBarrel = 1; + } + + return isBarrel; +} + +int DepthFlag( TVector3 inputPos ) //Used to calculate depth of given position... +{ + + float ShiftDHCAL = 530; // 20 Layers Inside DHCAL + + float DHCALDeepInnerZ = DHCALEndCapInnerZ + ShiftDHCAL; + float DHCALDeepBarrelRadius = DHCALBarrelRadius + ShiftDHCAL; + + float DHCALInnerOctRadius = DHCALBarrelRadius/cos(pi/4.0); + float DHCALDeepOctRadius = DHCALBarrelRadius/cos(pi/4.0) + ShiftDHCAL; + float ECALInnerOctRadius = ECALRadius/cos(pi/4.0); + + int FlagD(-1); + + if( fabs(inputPos[2]) > DHCALDeepInnerZ || fabs(inputPos[1]) > DHCALDeepBarrelRadius || fabs(inputPos[0]) > DHCALDeepBarrelRadius || fabs(inputPos[0] + inputPos[1]) > DHCALDeepOctRadius || fabs(inputPos[0] - inputPos[1]) > DHCALDeepOctRadius ) + { + FlagD = 2; + } + else if( fabs(inputPos[2]) > DHCALEndCapInnerZ || fabs(inputPos[1]) > DHCALBarrelRadius || fabs(inputPos[0]) > DHCALBarrelRadius || fabs(inputPos[0] + inputPos[1]) > DHCALInnerOctRadius || fabs(inputPos[0] - inputPos[1]) > DHCALInnerOctRadius ) + { + FlagD = 1; // Position outsider than DHCAL Region + } + else if( fabs(inputPos[2]) > ECALEndCapInnerZ || fabs(inputPos[1]) > ECALRadius || fabs(inputPos[0]) > ECALRadius || fabs(inputPos[0] + inputPos[1]) > ECALInnerOctRadius || fabs(inputPos[0] - inputPos[1]) > ECALInnerOctRadius ) + { + FlagD = 0; + } + else + { + FlagD = 10; // Position inside Calo... Problematic for Seeds... But could be PreShower hits. + } + + return FlagD; + +} + +TVector3 CalVertex( TVector3 Pos1, TVector3 Dir1, TVector3 Pos2, TVector3 Dir2 ) +{ + TVector3 VertexPos; + + float tau1(0), tau2(0); + + TVector3 delP; + delP = Pos1 - Pos2; + + double Normal(0); + Normal = (Dir1.Dot(Dir2))*(Dir1.Dot(Dir2)) - Dir1.Mag()*Dir1.Mag()*Dir2.Mag()*Dir2.Mag(); + + if(Normal != 0) + { + tau1 = (Dir2.Mag()*Dir2.Mag() * delP.Dot(Dir1) - Dir1.Dot(Dir2)*delP.Dot(Dir2))/Normal; + tau2 = (Dir1.Dot(Dir2)*delP.Dot(Dir1) - Dir1.Mag()*Dir1.Mag() * delP.Dot(Dir2))/Normal; + } + + VertexPos = 0.5*(Pos1 + Pos2 + tau1*Dir1 + tau2*Dir2); + + return VertexPos; +} + +int TPCPosition( TVector3 inputPos ) +{ + 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; + + return flagPos; +} + +float DisSeedSurface( TVector3 SeedPos ) //ECAL, HCAL, EndCapRing... +{ + + float DisSS = 0; + + if( fabs(SeedPos[2]) > ECALHalfZ ) //EcalEndcap hit start from 2350 + 100 = 2450 + { + + if( SeedPos.Perp() > ECALRadius ) + { + if( fabs(SeedPos[2])/SeedPos.Perp() > (ECALHalfZ + 103)/(ECALRadius + 100) ) + { + DisSS = ( fabs(SeedPos[2]) - ECALHalfZ - 103 ) * SeedPos.Mag()/fabs(SeedPos[2]); + } + else + { + DisSS = (SeedPos.Perp() - ECALRadius - 100 )*SeedPos.Mag()/SeedPos.Perp(); + } + } + else + { + DisSS = fabs(SeedPos[2]) - ECALHalfZ - 103; + } + } + else if( SeedPos.Perp() > ECALRadius + 400 ) + { + DisSS = SeedPos.Perp() - ECALRadius - 100; + } + else if( (SeedPos.Phi() > 0 && int(SeedPos.Phi() * 4/pi + 0.5) % 2 == 0 ) || (SeedPos.Phi() < 0 && int(SeedPos.Phi() * 4/pi + 8.5) % 2 == 0 )) + { + DisSS = min( fabs(fabs(SeedPos[0]) - ECALRadius), fabs(fabs(SeedPos[1]) - ECALRadius ) ); + } + else + { + DisSS = min( fabs(fabs(SeedPos[0] + SeedPos[1])/1.414214 -ECALRadius), fabs(fabs(SeedPos[0] - SeedPos[1])/1.414214 - ECALRadius) ); + } + + return DisSS; +} + +float DisSeedSurfaceSimple( TVector3 SeedPos ) //ECAL, HCAL, EndCapRing... +{ + // TVector3 SeedPos = a_clu->getPosition(); + float DisSS = 0; + float DisZ = abs(SeedPos.Z()) - ECALHalfZ; + float DisR = SeedPos.Perp() - ECALRadius; + if(DisR < 0 && DisZ > 0) + { + DisSS = DisZ; + } + else if(DisZ < 0 && DisR > 0) + { + DisSS = DisR; + } + else if(DisR > 0 && DisZ > 0) + { + if( fabs(SeedPos.Z())/SeedPos.Perp() > ECALHalfZ/ECALRadius) + { + DisSS = ( fabs(SeedPos[2]) - ECALHalfZ) * SeedPos.Mag()/fabs(SeedPos[2]); + } + else + { + DisSS = (SeedPos.Perp() - ECALRadius)*SeedPos.Mag()/SeedPos.Perp(); + } + } + + return DisSS; +} + +/* +float DisSeedSurfaceClu( Cluster * a_clu ) //ECAL, HCAL, EndCapRing... +{ + TVector3 SeedPos = a_clu->getPosition(); + return DisSeedSurface(SeedPos); +} +*/ + +float DisTPCBoundary( TVector3 Pos ) +{ + float DisZ = TMath::Min( fabs(ECALHalfZ-Pos.Z()),fabs(ECALHalfZ+Pos.Z()) ); + float DisR = ECALRadius - Pos.Perp(); + float Dis = TMath::Min(DisZ, DisR); + + 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; + + Diff_Clu_CPRef = CluPosition - CPRefPos; + NormCPRefMom = 1.0/CPRefMom.Mag()*CPRefMom; + float ProDis = Diff_Clu_CPRef.Dot(NormCPRefMom); + + DisCPClu = sqrt(Diff_Clu_CPRef.Mag()*Diff_Clu_CPRef.Mag() - ProDis*ProDis); + + return DisCPClu; +} + diff --git a/Reconstruction/PFA/Arbor/src/DetectorPos.hh b/Reconstruction/PFA/Arbor/src/DetectorPos.hh new file mode 100644 index 0000000000000000000000000000000000000000..8bf3f4ce8e48e383211cefc163f667528c6ab4c4 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/DetectorPos.hh @@ -0,0 +1,41 @@ +#ifndef DETECTORPOS_H_ +#define DETECTORPOS_H_ + +#include "TVector3.h" +#include "TMatrixF.h" +#include <iostream> +#include <vector> +#include <string> + +const float BField = 3.0; +const float DHCALBarrelRadius = 2058.0; //Octo +const float DHCALEndCapInnerZ = 2650.0; +const float ECALEndCapInnerZ = 2450.0; +const float ECALHalfZ = 2350.0; // mm, Endcap Ente +const float ECALRadius = 1847.4; // mm... minimal part for the octagnle. +const float TPCRadius = 1808.0 ; + +const float TPCOuterRadius = 1808.0; +const float TPCInnerRadius = 325.0; +const float LStar = 1500.0; + +const float DHCALCalibrationConstant = 0.12; + +int BarrelFlag( TVector3 inputPos ); + +int DepthFlag( TVector3 inputPos ); + +TVector3 CalVertex( TVector3 Pos1, TVector3 Dir1, TVector3 Pos2, TVector3 Dir2 ); + +int TPCPosition( TVector3 inputPos ); //Used to tag MCParticle position, if generated inside TPC & Dead outside + +float DisSeedSurface( TVector3 SeedPos ); //for a given position, calculate the distance to Calo surface ( ECAL ) + +float DisSeedSurfaceSimple( TVector3 SeedPos ); + + +float DisTPCBoundary( TVector3 Pos ); + +float DistanceChargedParticleToCluster(TVector3 CPRefPos, TVector3 CPRefMom, TVector3 CluPosition); + +#endif // diff --git a/Reconstruction/PFA/Arbor/src/HelixClassD.cc b/Reconstruction/PFA/Arbor/src/HelixClassD.cc new file mode 100644 index 0000000000000000000000000000000000000000..20b79466237d97a8d7bf5314ec00d3867cb5e453 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/HelixClassD.cc @@ -0,0 +1,801 @@ +#include "HelixClassD.hh" +#include <math.h> +#include <stdlib.h> +#include <iostream> +//#include "ced_cli.h" + +HelixClassD::HelixClassD() { + _const_2pi = 2.0*M_PI; + _const_pi2 = 0.5*M_PI; + _FCT = 2.99792458E-4; +} + +HelixClassD::~HelixClassD() {} + +void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) { + _referencePoint[0] = pos[0]; + _referencePoint[1] = pos[1]; + _referencePoint[2] = pos[2]; + _momentum[0] = mom[0]; + _momentum[1] = mom[1]; + _momentum[2] = mom[2]; + _charge = q; + _bField = B; + _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]); + _radius = _pxy / (_FCT*B); + _omega = q/_radius; + _tanLambda = mom[2]/_pxy; + _phiMomRefPoint = atan2(mom[1],mom[0]); + _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q); + _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q); + _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*q + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + // _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + double pxy = double(_pxy); + 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) { + d0 = double(q)*radius - double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + else { + d0 = double(q)*radius + double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + + _d0 = float(d0); + +// if (fabs(_d0)>0.001 ) { +// std::cout << "New helix : " << std::endl; +// 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 << " D0 = " << _d0 << std::endl; +// } + + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles); + +} + +void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, + float omega, float tanLambda, float B) { + _omega = omega; + _d0 = d0; + _phi0 = phi0; + _z0 = z0; + _tanLambda = tanLambda; + _charge = omega/fabs(omega); + _radius = 1./fabs(omega); + _xAtPCA = -_d0*sin(_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; + _pxAtPCA = _momentum[0]; + _pyAtPCA = _momentum[1]; + _phiMomRefPoint = atan2(_momentum[1],_momentum[0]); + _xCentre = _referencePoint[0] + + _radius*cos(_phi0-_const_pi2*_charge); + _yCentre = _referencePoint[1] + + _radius*sin(_phi0-_const_pi2*_charge); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phiRefPoint = _phiAtPCA ; + _bField = B; +} + + +void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, + float bZ, float phi0, float B, float signPz, + float zBegin) { + + _radius = radius; + _pxy = _FCT*B*_radius; + _charge = -(bZ*signPz)/fabs(bZ*signPz); + _momentum[2] = -_charge*_pxy/(bZ*_radius); + _xCentre = xCentre; + _yCentre = yCentre; + _omega = _charge/radius; + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*_charge + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + _referencePoint[2] = zBegin; + _referencePoint[0] = xCentre + radius*cos(bZ*zBegin+phi0); + _referencePoint[1] = yCentre + radius*sin(bZ*zBegin+phi0); + _phiRefPoint = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + _phiMomRefPoint = -_const_pi2*_charge + _phiRefPoint; + _tanLambda = _momentum[2]/_pxy; + _momentum[0] = _pxy*cos(_phiMomRefPoint); + _momentum[1] = _pxy*sin(_phiMomRefPoint); + + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = bZ*_referencePoint[2] - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + 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; + _bField = B; + +} + +const float * HelixClassD::getMomentum() { + return _momentum; +} +const float * HelixClassD::getReferencePoint() { + return _referencePoint; +} +float HelixClassD::getPhi0() { + if (_phi0<0.0) + _phi0 += 2*M_PI; + return _phi0; +} +float HelixClassD::getD0() { + return _d0; +} +float HelixClassD::getZ0() { + return _z0; +} +float HelixClassD::getOmega() { + return _omega; +} +float HelixClassD::getTanLambda() { + return _tanLambda; +} +float HelixClassD::getPXY() { + return _pxy; +} +float HelixClassD::getXC() { + return _xCentre; +} + +float HelixClassD::getYC() { + return _yCentre; +} + +float HelixClassD::getRadius() { + return _radius; +} + +float HelixClassD::getBz() { + return _bZ; +} + +float HelixClassD::getPhiZ() { + return _phiZ; +} + +float HelixClassD::getCharge() { + return _charge; +} + +float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, + float * ref , float * point) { + + float time; + + float AA = sqrt(ax*ax+ay*ay); + + + if (AA <= 0) { + time = -1.0e+20; + return time; + } + + + float BB = ax*(x0-_xCentre) + ay*(y0-_yCentre); + BB = BB / AA; + + float CC = (x0-_xCentre)*(x0-_xCentre) + + (y0-_yCentre)*(y0-_yCentre) - _radius*_radius; + + CC = CC / AA; + + float DET = BB*BB - CC; + float tt1 = 0.; + float tt2 = 0.; + float xx1,xx2,yy1,yy2; + + + if (DET < 0 ) { + time = 1.0e+10; + point[0]=0.0; + point[1]=0.0; + point[2]=0.0; + return time; + } + + + tt1 = - BB + sqrt(DET); + tt2 = - BB - sqrt(DET); + + xx1 = x0 + tt1*ax; + 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); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + 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) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + time = tt1; + } + else { + point[0] = xx2; + point[1] = yy2; + time = tt2; + } + + point[2] = ref[2]+time*_momentum[2]; + + + + return time; + +} + + +float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) { + + float distCenterToIP = sqrt(_xCentre*_xCentre + _yCentre*_yCentre); + + point[0] = 0.0; + point[1] = 0.0; + point[2] = 0.0; + + if ((distCenterToIP+_radius)<Radius) { + float xx = -1.0e+20; + return xx; + } + + if ((_radius+Radius)<distCenterToIP) { + float xx = -1.0e+20; + return xx; + } + + float phiCentre = atan2(_yCentre,_xCentre); + float phiStar = Radius*Radius + distCenterToIP*distCenterToIP + - _radius*_radius; + + phiStar = 0.5*phiStar/fmax(1.0e-20,Radius*distCenterToIP); + + if (phiStar > 1.0) + phiStar = 0.9999999; + + if (phiStar < -1.0) + phiStar = -0.9999999; + + phiStar = acos(phiStar); + + float tt1,tt2,time; + + float xx1 = Radius*cos(phiCentre+phiStar); + float yy1 = Radius*sin(phiCentre+phiStar); + + float xx2 = Radius*cos(phiCentre-phiStar); + float yy2 = Radius*sin(phiCentre-phiStar); + + + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); + float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); + float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + 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) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + float time2; + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + point[3] = xx2; + point[4] = yy2; + time = tt1; + time2 = tt2; + } + else { + point[0] = xx2; + point[1] = yy2; + point[3] = xx1; + point[4] = yy1; + time = tt2; + time2 = tt1; + } + + point[2] = ref[2]+time*_momentum[2]; + point[5] = ref[2]+time2*_momentum[2]; + + + return time; + +} + + +float HelixClassD::getPointInZ(float zLine, float * ref, float * point) { + + float time = zLine - ref[2]; + + if (_momentum[2] == 0.) { + time = -1.0e+20; + point[0] = 0.; + point[1] = 0.; + point[2] = 0.; + return time; + } + + time = time/_momentum[2]; + + float phi0 = atan2(ref[1] - _yCentre , ref[0] - _xCentre); + float phi = phi0 - _charge*_pxy*time/_radius; + float xx = _xCentre + _radius*cos(phi); + float yy = _yCentre + _radius*sin(phi); + + point[0] = xx; + point[1] = yy; + point[2] = zLine; + + return time; + + +} + +int HelixClassD::getNCircle(float * xPoint) { + + int nCircles = 0; + + float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + + if (fabs(_tanLambda*_radius)>1.0e-20) { + float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius); + xCircles = xCircles/_const_2pi; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + + } + return nCircles; + +} + +float HelixClassD::getDistanceToPoint(float * xPoint, float * Distance) { + + float zOnHelix; + float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + 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) { + float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius); + xCircles = xCircles/_const_2pi; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + + } + + float DPhi = _const_2pi*((float)nCircles) + phi - phi0; + + zOnHelix = _referencePoint[2] - _charge*_radius*_tanLambda*DPhi; + + float DistZ = fabs(zOnHelix - xPoint[2]); + float Time; + + if (fabs(_momentum[2]) > 1.0e-20) { + Time = (zOnHelix - _referencePoint[2])/_momentum[2]; + } + else { + Time = _charge*_radius*DPhi/_pxy; + } + + Distance[0] = DistXY; + Distance[1] = DistZ; + Distance[2] = sqrt(DistXY*DistXY+DistZ*DistZ); + + return Time; + + +} + +//When we are not interested in the exact distance, we can check if we are +//already far enough away in XY, before we start calculating in Z as the +//distance will only increase +float HelixClassD::getDistanceToPoint(const std::vector<float>& xPoint, float distCut) { + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + float tempx = xPoint[0]-_xCentre; + float tempy = xPoint[1]-_yCentre; + float tempsq = sqrt(tempx*tempx + tempy*tempy); + float tempdf = tempsq - _radius; + float DistXY = fabs( tempdf ); + //If this is bigger than distCut, we dont have to know how much bigger this is + if( DistXY > distCut) { + return DistXY; + } + + int nCircles = 0; + float phi = atan2(tempy,tempx); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + float phidiff = phi0-phi; + float tempz = xPoint[2] - _referencePoint[2];//Yes referencePoint + float tanradius = _tanLambda*_radius; + if (fabs(tanradius)>1.0e-20) { + float xCircles = (phidiff -_charge*tempz/tanradius)/_const_2pi; + int n1,n2; + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + } + float DistZ = - tempz - _charge*tanradius*(_const_2pi*((float)nCircles) - phidiff); + return sqrt(DistXY*DistXY+DistZ*DistZ); +}//getDistanceToPoint(vector,float) + +float HelixClassD::getDistanceToPoint(const float* xPoint, float distCut) { + std::vector<float> xPosition(xPoint, xPoint + 3 );//We are expecting three coordinates, must be +3, last element is excluded! + return getDistanceToPoint(xPosition, distCut); +}//getDistanceToPoint(float*,float) + + + +void HelixClassD::setHelixEdges(float * xStart, float * xEnd) { + for (int i=0; i<3; ++i) { + _xStart[i] = xStart[i]; + _xEnd[i] = xEnd[i]; + } + +} + +float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * mom) { + + float x01 = getXC(); + float y01 = getYC(); + + float x02 = helix->getXC(); + float y02 = helix->getYC(); + + float rad1 = getRadius(); + float rad2 = helix->getRadius(); + + float distance = sqrt((x01-x02)*(x01-x02)+(y01-y02)*(y01-y02)); + + bool singlePoint = true; + + float phi1 = 0; + float phi2 = 0; + + if (rad1+rad2<distance) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = atan2(y01-y02,x01-x02); + } + else if (distance+rad2<rad1) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = phi1; + } + else if (distance+rad1<rad2) { + phi1 = atan2(y01-y02,x01-x02); + phi2 = phi1; + } + else { + singlePoint = false; + float cosAlpha = 0.5*(distance*distance+rad2*rad2-rad1*rad1)/(distance*rad2); + float alpha = acos(cosAlpha); + float phi0 = atan2(y01-y02,x01-x02); + phi1 = phi0 + alpha; + phi2 = phi0 - alpha; + } + + + float ref1[3]; + float ref2[3]; + for (int i=0;i<3;++i) { + ref1[i]=_referencePoint[i]; + ref2[i]=helix->getReferencePoint()[i]; + } + + float pos1[3]; + float pos2[3]; + float mom1[3]; + float mom2[3]; + + + if (singlePoint ) { + + float xSect1 = x01 + rad1*cos(phi1); + float ySect1 = y01 + rad1*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + float R1 = sqrt(xSect1*xSect1+ySect1*ySect1); + float R2 = sqrt(xSect2*xSect2+ySect2*ySect2); + + getPointOnCircle(R1,ref1,pos1); + helix->getPointOnCircle(R2,ref2,pos2); + + } + else { + + float xSect1 = x02 + rad2*cos(phi1); + float ySect1 = y02 + rad2*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + +// std::cout << "(xSect1,ySect1)=(" << xSect1 << "," << ySect1 << ")" << std::endl; +// std::cout << "(xSect2,ySect2)=(" << xSect2 << "," << ySect2 << ")" << std::endl; + + float temp21[3]; + float temp22[3]; + + 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; + float deltaPhi22 = phiF22 - phiI2; + float charge2 = helix->getCharge(); + float pxy2 = helix->getPXY(); + float pz2 = helix->getMomentum()[2]; + if (deltaPhi21 < 0 && charge2 < 0) { + deltaPhi21 += _const_2pi; + } + else if (deltaPhi21 > 0 && charge2 > 0) { + deltaPhi21 -= _const_2pi; + } + + if (deltaPhi22 < 0 && charge2 < 0) { + deltaPhi22 += _const_2pi; + } + else if (deltaPhi22 > 0 && charge2 > 0) { + deltaPhi22 -= _const_2pi; + } + + float time21 = -charge2*deltaPhi21*rad2/pxy2; + float time22 = -charge2*deltaPhi22*rad2/pxy2; + + float Z21 = ref2[2] + time21*pz2; + float Z22 = ref2[2] + time22*pz2; + + 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; + + + float temp11[3]; + float temp12[3]; + + 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; + float deltaPhi12 = phiF12 - phiI1; + float charge1 = _charge; + float pxy1 = _pxy; + float pz1 = _momentum[2]; + if (deltaPhi11 < 0 && charge1 < 0) { + deltaPhi11 += _const_2pi; + } + else if (deltaPhi11 > 0 && charge1 > 0) { + deltaPhi11 -= _const_2pi; + } + + if (deltaPhi12 < 0 && charge1 < 0) { + deltaPhi12 += _const_2pi; + } + else if (deltaPhi12 > 0 && charge1 > 0) { + deltaPhi12 -= _const_2pi; + } + + float time11 = -charge1*deltaPhi11*rad1/pxy1; + float time12 = -charge1*deltaPhi12*rad1/pxy1; + + float Z11 = ref1[2] + time11*pz1; + float Z12 = ref1[2] + time12*pz1; + + 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; + + float Dist1 = 0; + float Dist2 = 0; + + for (int j=0;j<3;++j) { + Dist1 += (temp11[j]-temp21[j])*(temp11[j]-temp21[j]); + Dist2 += (temp12[j]-temp22[j])*(temp12[j]-temp22[j]); + } + + if (Dist1<Dist2) { + for (int l=0;l<3;++l) { + pos1[l] = temp11[l]; + pos2[l] = temp21[l]; + } + } + else { + for (int l=0;l<3;++l) { + pos1[l] = temp12[l]; + pos2[l] = temp22[l]; + } + } + + } + + getExtrapolatedMomentum(pos1,mom1); + helix->getExtrapolatedMomentum(pos2,mom2); + + float helixDistance = 0; + + for (int i=0;i<3;++i) { + helixDistance += (pos1[i] - pos2[i])*(pos1[i] - pos2[i]); + pos[i] = 0.5*(pos1[i]+pos2[i]); + mom[i] = mom1[i]+mom2[i]; + } + helixDistance = sqrt(helixDistance); + + return helixDistance; + +} + +void HelixClassD::getExtrapolatedMomentum(float * pos, float * momentum) { + + float phi = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + if (phi <0.) phi += _const_2pi; + phi = phi - _phiAtPCA + _phi0; + momentum[0] = _pxy*cos(phi); + momentum[1] = _pxy*sin(phi); + momentum[2] = _momentum[2]; + + +} diff --git a/Reconstruction/PFA/Arbor/src/HelixClassD.hh b/Reconstruction/PFA/Arbor/src/HelixClassD.hh new file mode 100644 index 0000000000000000000000000000000000000000..c9e4d4c386addf9312eb56ac8a89e25b6260820f --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/HelixClassD.hh @@ -0,0 +1,304 @@ +#ifndef HELIXARD_HH +#define HELIXARD_HH 1 +#include <vector> +/** + * Utility class to manipulate with different parameterisations <br> + * of helix. Helix can be initialized in a three different <br> + * ways using the following public methods : <br> + * 1) Initialize_VP(float * pos, float * mom, float q, float B) : <br> + * initialization of helix is done using <br> + * - position of the reference point : pos[3]; <br> + * - momentum vector at the reference point : mom[3];<br> + * - particle charge : q;<br> + * - magnetic field (in Tesla) : B;<br> + * 2) Initialize_BZ(float xCentre, float yCentre, float radius, <br> + * float bZ, float phi0, float B, float signPz,<br> + * float zBegin):<br> + * initialization of helix is done according to the following<br> + * parameterization<br> + * x = xCentre + radius*cos(bZ*z + phi0)<br> + * y = yCentre + radius*sin(bZ*z + phi0)<br> + * where (x,y,z) - position of point on the helix<br> + * - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br> + * - radius is the radius of circumference<br> + * - bZ is the helix slope parameter<br> + * - phi0 is the initial phase of circumference<br> + * - B is the magnetic field (in Tesla)<br> + * - signPz is the sign of the z component of momentum vector<br> + * - zBegin is the z coordinate of the reference point;<br> + * 3) void Initialize_Canonical(float phi0, float d0, float z0, float omega,<br> + * float tanlambda, float B) :<br> + * canonical (LEP-wise) parameterisation with the following parameters<br> + * - phi0 - Phi angle of momentum vector at the point of<br> + * closest approach to IP in R-Phi plane; + * - d0 - signed distance of closest approach to IP in R-Phi plane;<br> + * - z0 - z coordinate of the point of closest approach in R-Phi plane;<br> + * - omega - signed curvature;<br> + * - tanlambda - tangent of dip angle;<br> + * - B - magnetic field (in Tesla);<br> + * A set of public methods (getters) provide access to <br> + * various parameters associated with helix. Helix Class contains<br> + * also several utility methods, allowing for calculation of helix<br> + * intersection points with planes parallel and perpendicular to <br> + * z (beam) axis and determination of the distance of closest approach<br> + * from arbitrary 3D point to the helix. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: HelixClassD.h,v 1.16 2008-06-05 13:47:18 rasp Exp $<br> + * + */ + +//#include "LineClass.h" +class HelixClassD; + +class HelixClassD { + public: + +/** + * Constructor. Initializations of constants which are used + * to calculate various parameters associated with helix. + */ + HelixClassD(); +/** + * Destructor. + */ + ~HelixClassD(); +/** + * Initialization of helix using <br> + * - position of the reference point : pos[3]; <br> + * - momentum vector at the reference point : mom[3];<br> + * - particle charge : q;<br> + * - magnetic field (in Tesla) : B<br> + */ + void Initialize_VP(float * pos, float * mom, float q, float B); + +/** + * Initialization of helix according to the following<br> + * parameterization<br> + * x = xCentre + radius*cos(bZ*z + phi0)<br> + * y = yCentre + radius*sin(bZ*z + phi0)<br> + * where (x,y,z) - position of point on the helix<br> + * - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br> + * - radius is the radius of circumference<br> + * - bZ is the helix slope parameter<br> + * - phi0 is the initial phase of circumference<br> + * - B is the magnetic field (in Tesla)<br> + * - signPz is the sign of the z component of momentum vector<br> + * - zBegin is the z coordinate of the reference point<br> + */ + void Initialize_BZ(float xCentre, float yCentre, float radius, + float bZ, float phi0, float B, float signPz, + float zBegin); +/** + * Canonical (LEP-wise) parameterisation with the following parameters<br> + * - phi0 - Phi angle of momentum vector at the point of<br> + * closest approach to IP in R-Phi plane; + * - d0 - signed distance of closest approach in R-Phi plane;<br> + * - z0 - z coordinate of the point of closest approach to IP + * in R-Phi plane;<br> + * - omega - signed curvature;<br> + * - tanlambda - tangent of dip angle;<br> + * - B - magnetic field (in Tesla)<br> + */ + void Initialize_Canonical(float phi0, float d0, float z0, float omega, + float tanlambda, float B); + /** + * Returns momentum of particle at the point of closest approach <br> + * to IP <br> + */ + const float * getMomentum(); + + /** + * Returns reference point of track <br> + */ + const float * getReferencePoint(); + + /** + * Returns Phi angle of the momentum vector <br> + * at the point of closest approach to IP <br> + */ + float getPhi0(); + + /** + * Returns signed distance of closest approach <br> + * to IP in the R-Phi plane <br> + */ + float getD0(); + + /** + * Returns z coordinate of the point of closest + * approach to IP in the R-Phi plane <br> + */ + float getZ0(); + + /** + * Returns signed curvature of the track <br> + */ + float getOmega(); + + /** + * Returns tangent of dip angle of the track <br> + */ + float getTanLambda(); + + /** + * Returns transverse momentum of the track <br> + */ + float getPXY(); + + + /** + * Returns x coordinate of circumference + */ + float getXC(); + + /** + * Returns y coordinate of circumference + */ + float getYC(); + + + /** + * Returns radius of circumference + */ + float getRadius(); + + + /** + * Returns helix intersection point with the plane <br> + * parallel to z axis. Plane is defined by two coordinates <br> + * of the point belonging to the plane (x0,y0) and normal <br> + * vector (ax,ay). ref[3] is the reference point of the helix. <br> + * point[3] - returned vector holding the coordinates of <br> + * intersection point <br> + */ + float getPointInXY(float x0, float y0, float ax, float ay, + float * ref , float * point); + + /** + * Returns helix intersection point with the plane <br> + * perpendicular to z axis. Plane is defined by z coordinate <br> + * of the plane. ref[3] is the reference point of the helix. <br> + * point[3] - returned vector holding the coordinates of <br> + * intersection point <br> + */ + float getPointInZ(float zLine, float * ref, float * point); + + /** + * Return distance of the closest approach of the helix to <br> + * arbitrary 3D point in space. xPoint[3] - coordinates of <br> + * space point. Distance[3] - vector of distances of helix to <br> + * a given point in various projections : <br> + * Distance[0] - distance in R-Phi plane <br> + * Distance[1] - distance along Z axis <br> + * Distance[2] - 3D distance <br> + */ + float getDistanceToPoint(float * xPoint, float * Distance); + + /** + * Return distance of the closest approach of the helix to <br> + * arbitrary 3D point in space. xPoint[3] - coordinates of <br> + * space point. distCut - limit on the distance between helix <br> + * and the point to reduce calculation time <br> + * If R-Phi is found to be greater than distCut, rPhi distance is returned <br> + * If the R-Phi distance is not too big, than the exact 3D distance is returned <br> + * This function can be used, if the exact distance is not always needed <br> + */ + float getDistanceToPoint(const float* xPoint, float distCut); + float getDistanceToPoint(const std::vector<float>& xPoint, float distCut); + + /** + * This method calculates coordinates of both intersection <br> + * of the helix with a cylinder. <br> + * Rotational symmetry with respect to z axis is assumed, <br> + * meaning that cylinder axis corresponds to the z axis <br> + * of reference frame. <br> + * Inputs : <br> + * Radius - radius of cylinder. <br> + * ref[3] - point of closest approach to the origin of the helix. <br> + * Output : <br> + * point[6] - coordinates of intersection point. <br> + * Method returns also generic time, defined as the <br> + * ratio of helix length from reference point to the intersection <br> + * point to the particle momentum <br> + */ + float getPointOnCircle(float Radius, float * ref, float * point); + + /** Returns distance between two helixes <br> + * Output : <br> + * pos[3] - position of the point of closest approach <br> + * mom[3] - momentum of V0 <br> + */ + float getDistanceToHelix(HelixClassD * helix, float * pos, float * mom); + + int getNCircle(float * XPoint); + + /** + * Set Edges of helix + */ + void setHelixEdges(float * xStart, float * xEnd); + + /** + * Returns starting point of helix + */ + float * getStartingPoint() {return _xStart;} + + /** + * Returns endpoint of helix + */ + float * getEndPoint() {return _xEnd;} + + /** + * Returns BZ for the second parameterization + */ + float getBz(); + + /** + * Returns Phi for the second parameterization + */ + float getPhiZ(); + + /** + * Returns extrapolated momentum + */ + void getExtrapolatedMomentum(float * pos, float * momentum); + + /** + * Returns charge + */ + float getCharge(); + + private: + float _momentum[3]; // momentum @ ref point + float _referencePoint[3]; // coordinates @ ref point + float _phi0; // phi0 in canonical parameterization + float _d0; // d0 in canonical parameterisation + float _z0; // z0 in canonical parameterisation + float _omega; // signed curvuture in canonical parameterisation + float _tanLambda; // TanLambda + float _pxy; // Transverse momentum + float _charge; // Particle Charge + float _bField; // Magnetic field (assumed to point to Z>0) + float _radius; // radius of circle in XY plane + float _xCentre; // X of circle centre + float _yCentre; // Y of circle centre + float _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point + float _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA + float _xAtPCA; // X @ PCA + float _yAtPCA; // Y @ PCA + float _pxAtPCA; // PX @ PCA + float _pyAtPCA; // PY @ PCA + float _phiMomRefPoint; // Phi of Momentum vector @ ref point + float _const_pi; // PI + float _const_2pi; // 2*PI + float _const_pi2; // PI/2 + float _FCT; // 2.99792458E-4 + float _xStart[3]; // Starting point of track segment + float _xEnd[3]; // Ending point of track segment + + float _bZ; + float _phiZ; + +}; + + +#endif diff --git a/Reconstruction/PFA/Arbor/src/KDTreeLinkerAlgoT.h b/Reconstruction/PFA/Arbor/src/KDTreeLinkerAlgoT.h new file mode 100644 index 0000000000000000000000000000000000000000..01ce7affa77a6820e24c6b21daa81cf8b4ecaa94 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/KDTreeLinkerAlgoT.h @@ -0,0 +1,336 @@ +#ifndef KDTreeLinkerAlgoTemplated_h +#define KDTreeLinkerAlgoTemplated_h + +#include "KDTreeLinkerToolsT.h" + +#include <cassert> +#include <vector> + +// Class that implements the KDTree partition of 2D space and +// a closest point search algorithme. + +template <typename DATA, unsigned DIM=2> +class KDTreeLinkerAlgo +{ + public: + KDTreeLinkerAlgo(); + + // Dtor calls clear() + ~KDTreeLinkerAlgo(); + + // Here we build the KD tree from the "eltList" in the space define by "region". + void build(std::vector<KDTreeNodeInfoT<DATA,DIM> > &eltList, + const KDTreeBoxT<DIM> ®ion); + + // Here we search in the KDTree for all points that would be + // contained in the given searchbox. The founded points are stored in resRecHitList. + void search(const KDTreeBoxT<DIM> &searchBox, + std::vector<KDTreeNodeInfoT<DATA,DIM> > &resRecHitList); + + // This reurns true if the tree is empty + bool empty() {return nodePoolPos_ == -1;} + + // This returns the number of nodes + leaves in the tree + // (nElements should be (size() +1)/2) + int size() { return nodePoolPos_ + 1;} + + // This method clears all allocated structures. + void clear(); + + private: + // The KDTree root + KDTreeNodeT<DATA,DIM>* root_; + + // The node pool allow us to do just 1 call to new for each tree building. + KDTreeNodeT<DATA,DIM>* nodePool_; + int nodePoolSize_; + int nodePoolPos_; + + + + std::vector<KDTreeNodeInfoT<DATA,DIM> > *closestNeighbour; + std::vector<KDTreeNodeInfoT<DATA,DIM> > *initialEltList; + + private: + + // Get next node from the node pool. + KDTreeNodeT<DATA,DIM>* getNextNode(); + + //Fast median search with Wirth algorithm in eltList between low and high indexes. + int medianSearch(int low, + int high, + int treeDepth); + + // Recursif kdtree builder. Is called by build() + KDTreeNodeT<DATA,DIM> *recBuild(int low, + int hight, + int depth, + const KDTreeBoxT<DIM> ®ion); + + // Recursif kdtree search. Is called by search() + void recSearch(const KDTreeNodeT<DATA,DIM> *current, + const KDTreeBoxT<DIM> &trackBox); + + // Add all elements of an subtree to the closest elements. Used during the recSearch(). + void addSubtree(const KDTreeNodeT<DATA,DIM> *current); + + // This method frees the KDTree. + void clearTree(); +}; + + +//Implementation + +template < typename DATA, unsigned DIM > +void +KDTreeLinkerAlgo<DATA,DIM>::build(std::vector<KDTreeNodeInfoT<DATA,DIM> > &eltList, + const KDTreeBoxT<DIM> ®ion) +{ + if (eltList.size()) { + initialEltList = &eltList; + + size_t size = initialEltList->size(); + nodePoolSize_ = size * 2 - 1; + nodePool_ = new KDTreeNodeT<DATA,DIM>[nodePoolSize_]; + + // Here we build the KDTree + root_ = recBuild(0, size, 0, region); + + initialEltList = 0; + } +} + +//Fast median search with Wirth algorithm in eltList between low and high indexes. +template < typename DATA, unsigned DIM > +int +KDTreeLinkerAlgo<DATA,DIM>::medianSearch(int low, + int high, + int treeDepth) +{ + //We should have at least 1 element to calculate the median... + //assert(low < high); + + const int nbrElts = high - low; + int median = nbrElts/2 - ( 1 - 1*(nbrElts&1) ); + median += low; + + int l = low; + int m = high - 1; + + while (l < m) { + KDTreeNodeInfoT<DATA,DIM> elt = (*initialEltList)[median]; + int i = l; + int j = m; + + do { + // The even depth is associated to dim1 dimension + // The odd one to dim2 dimension + const unsigned thedim = treeDepth % DIM; + while( (*initialEltList)[i].dims[thedim] < elt.dims[thedim] ) ++i; + while( (*initialEltList)[j].dims[thedim] > elt.dims[thedim] ) --j; + + if (i <= j){ + std::swap((*initialEltList)[i], (*initialEltList)[j]); + i++; + j--; + } + } while (i <= j); + if (j < median) l = i; + if (i > median) m = j; + } + + return median; +} + + + +template < typename DATA, unsigned DIM > +void +KDTreeLinkerAlgo<DATA,DIM>::search(const KDTreeBoxT<DIM> &trackBox, + std::vector<KDTreeNodeInfoT<DATA,DIM> > &recHits) +{ + if (root_) { + closestNeighbour = &recHits; + recSearch(root_, trackBox); + closestNeighbour = 0; + } +} + + +template < typename DATA, unsigned DIM > +void +KDTreeLinkerAlgo<DATA,DIM>::recSearch(const KDTreeNodeT<DATA,DIM> *current, + const KDTreeBoxT<DIM> &trackBox) +{ + /* + // By construction, current can't be null + assert(current != 0); + + // By Construction, a node can't have just 1 son. + assert (!(((current->left == 0) && (current->right != 0)) || + ((current->left != 0) && (current->right == 0)))); + */ + + if ((current->left == 0) && (current->right == 0)) {//leaf case + + // If point inside the rectangle/area + bool isInside = true; + for( unsigned i = 0; i < DIM; ++i ) { + const auto thedim = current->info.dims[i]; + isInside *= thedim >= trackBox.dimmin[i] && thedim <= trackBox.dimmax[i]; + } + if( isInside ) closestNeighbour->push_back(current->info); + + } else { + + bool isFullyContained = true; + bool hasIntersection = true; + //if region( v->left ) is fully contained in the rectangle + for( unsigned i = 0; i < DIM; ++i ) { + const auto regionmin = current->left->region.dimmin[i]; + const auto regionmax = current->left->region.dimmax[i]; + isFullyContained *= ( regionmin >= trackBox.dimmin[i] && + regionmax <= trackBox.dimmax[i] ); + hasIntersection *= ( regionmin < trackBox.dimmax[i] && regionmax > trackBox.dimmin[i]); + } + if( isFullyContained ) { + addSubtree(current->left); + } else if ( hasIntersection ) { + recSearch(current->left, trackBox); + } + // reset flags + isFullyContained = true; + hasIntersection = true; + //if region( v->right ) is fully contained in the rectangle + for( unsigned i = 0; i < DIM; ++i ) { + const auto regionmin = current->right->region.dimmin[i]; + const auto regionmax = current->right->region.dimmax[i]; + isFullyContained *= ( regionmin >= trackBox.dimmin[i] && + regionmax <= trackBox.dimmax[i] ); + hasIntersection *= ( regionmin < trackBox.dimmax[i] && regionmax > trackBox.dimmin[i]); + } + if( isFullyContained ) { + addSubtree(current->right); + } else if ( hasIntersection ) { + recSearch(current->right, trackBox); + } + } +} + +template < typename DATA, unsigned DIM > +void +KDTreeLinkerAlgo<DATA,DIM>::addSubtree(const KDTreeNodeT<DATA,DIM> *current) +{ + // By construction, current can't be null + // assert(current != 0); + + if ((current->left == 0) && (current->right == 0)) // leaf + closestNeighbour->push_back(current->info); + else { // node + addSubtree(current->left); + addSubtree(current->right); + } +} + + + + +template <typename DATA, unsigned DIM> +KDTreeLinkerAlgo<DATA,DIM>::KDTreeLinkerAlgo() + : root_ (0), + nodePool_(0), + nodePoolSize_(-1), + nodePoolPos_(-1) +{ +} + +template <typename DATA, unsigned DIM> +KDTreeLinkerAlgo<DATA,DIM>::~KDTreeLinkerAlgo() +{ + clear(); +} + + +template <typename DATA, unsigned DIM> +void +KDTreeLinkerAlgo<DATA,DIM>::clearTree() +{ + delete[] nodePool_; + nodePool_ = 0; + root_ = 0; + nodePoolSize_ = -1; + nodePoolPos_ = -1; +} + +template <typename DATA, unsigned DIM> +void +KDTreeLinkerAlgo<DATA,DIM>::clear() +{ + if (root_) + clearTree(); +} + + +template <typename DATA, unsigned DIM> +KDTreeNodeT<DATA,DIM>* +KDTreeLinkerAlgo<DATA,DIM>::getNextNode() +{ + ++nodePoolPos_; + + // The tree size is exactly 2 * nbrElts - 1 and this is the total allocated memory. + // If we have used more than that....there is a big problem. + // assert(nodePoolPos_ < nodePoolSize_); + + return &(nodePool_[nodePoolPos_]); +} + + +template <typename DATA, unsigned DIM> +KDTreeNodeT<DATA,DIM>* +KDTreeLinkerAlgo<DATA,DIM>::recBuild(int low, + int high, + int depth, + const KDTreeBoxT<DIM>& region) +{ + int portionSize = high - low; + + // By construction, portionSize > 0 can't happend. + // assert(portionSize > 0); + + if (portionSize == 1) { // Leaf case + + KDTreeNodeT<DATA,DIM> *leaf = getNextNode(); + leaf->setAttributs(region, (*initialEltList)[low]); + return leaf; + + } else { // Node case + + // The even depth is associated to dim1 dimension + // The odd one to dim2 dimension + int medianId = medianSearch(low, high, depth); + + // We create the node + KDTreeNodeT<DATA,DIM> *node = getNextNode(); + node->setAttributs(region); + + + // Here we split into 2 halfplanes the current plane + KDTreeBoxT<DIM> leftRegion = region; + KDTreeBoxT<DIM> rightRegion = region; + unsigned thedim = depth % DIM; + auto medianVal = (*initialEltList)[medianId].dims[thedim]; + leftRegion.dimmax[thedim] = medianVal; + rightRegion.dimmin[thedim] = medianVal; + + ++depth; + ++medianId; + + // We recursively build the son nodes + node->left = recBuild(low, medianId, depth, leftRegion); + node->right = recBuild(medianId, high, depth, rightRegion); + + return node; + } +} + +#endif diff --git a/Reconstruction/PFA/Arbor/src/KDTreeLinkerToolsT.h b/Reconstruction/PFA/Arbor/src/KDTreeLinkerToolsT.h new file mode 100644 index 0000000000000000000000000000000000000000..78e5696c9f3d1f8bb5ea846ee69a0e5b0665ebda --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/KDTreeLinkerToolsT.h @@ -0,0 +1,82 @@ +#ifndef KDTreeLinkerToolsTemplated_h +#define KDTreeLinkerToolsTemplated_h + +#include <array> + +// Box structure used to define 2D field. +// It's used in KDTree building step to divide the detector +// space (ECAL, HCAL...) and in searching step to create a bounding +// box around the demanded point (Track collision point, PS projection...). +template<unsigned DIM> +struct KDTreeBoxT +{ + std::array<float,DIM> dimmin, dimmax; + + template<typename... Ts> + KDTreeBoxT(Ts... dimargs) { + static_assert(sizeof...(dimargs) == 2*DIM,"Constructor requires 2*DIM args"); + std::vector<float> dims = {dimargs...}; + for( unsigned i = 0; i < DIM; ++i ) { + dimmin[i] = dims[2*i]; + dimmax[i] = dims[2*i+1]; + } + } + + KDTreeBoxT() {} +}; + +typedef KDTreeBoxT<2> KDTreeBox; +typedef KDTreeBoxT<3> KDTreeCube; + + +// Data stored in each KDTree node. +// The dim1/dim2 fields are usually the duplication of some PFRecHit values +// (eta/phi or x/y). But in some situations, phi field is shifted by +-2.Pi +template<typename DATA,unsigned DIM> +struct KDTreeNodeInfoT +{ + DATA data; + std::array<float,DIM> dims; + +public: + KDTreeNodeInfoT() + {} + + template<typename... Ts> + KDTreeNodeInfoT(const DATA& d,Ts... dimargs) + : data(d), dims{ {dimargs...} } + {} +}; + +// KDTree node. +template <typename DATA, unsigned DIM> +struct KDTreeNodeT +{ + // Data + KDTreeNodeInfoT<DATA,DIM> info; + + // Right/left sons. + KDTreeNodeT<DATA,DIM> *left, *right; + + // Region bounding box. + KDTreeBoxT<DIM> region; + + public: + KDTreeNodeT() + : left(0), right(0) + {} + + void setAttributs(const KDTreeBoxT<DIM>& regionBox, + const KDTreeNodeInfoT<DATA,DIM>& infoToStore) + { + info = infoToStore; + region = regionBox; + } + + void setAttributs(const KDTreeBoxT<DIM>& regionBox) + { + region = regionBox; + } +}; + +#endif diff --git a/Reconstruction/PFA/Arbor/src/MarlinArbor.cc b/Reconstruction/PFA/Arbor/src/MarlinArbor.cc new file mode 100755 index 0000000000000000000000000000000000000000..e8b6e7d21cf5ba1d69b655292d13dab780da146e --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/MarlinArbor.cc @@ -0,0 +1,219 @@ +#include "MarlinArbor.hh" +#include "ArborTool.h" +#include "ArborToolLCIO.hh" +#include "ArborHit.h" + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/Cluster.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" + +#include "cellIDDecoder.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/IDDescriptor.h" +#include "DD4hep/Plugins.h" + +#include <values.h> +#include <string> +#include <iostream> +#include <cmath> +#include <stdexcept> +#include <TFile.h> +#include <TTree.h> +#include <TMath.h> +#include <Rtypes.h> +#include <sstream> +#include <set> +#include <TVector3.h> +#include <vector> +#include <algorithm> + +#include "DetectorPos.hh" + +using namespace std; + +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; + +DECLARE_COMPONENT(MarlinArbor) + +MarlinArbor::MarlinArbor(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc), + _eventNr(0),_output(0) +{ +} + +StatusCode MarlinArbor::initialize() { + + + _cepc_thresholds.push_back(10); + _cepc_thresholds.push_back(90); + _cepc_thresholds.push_back(50); + _cepc_thresholds.push_back(7.5); + + m_geosvc = service<IGeomSvc>("GeomSvc"); + + ISvcLocator* svcloc = serviceLocator(); + m_ArborToolLCIO=new ArborToolLCIO("arborTools",svcloc); + 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); + } + for(unsigned int i = 0; i < m_hcalReadoutNames.value().size(); i++){ + m_col_readout_map[m_hcalColNames.value().at(i)] = m_hcalReadoutNames.value().at(i); + } + + for (auto& ecal : m_ecalColNames) { + _ecalCollections.push_back( new CaloType(ecal, Gaudi::DataHandle::Reader, this) ); + _calCollections.push_back( new CaloType(ecal, Gaudi::DataHandle::Reader, this) ); + } + for (auto& hcal : m_hcalColNames) { + _hcalCollections.push_back( new CaloType(hcal, Gaudi::DataHandle::Reader, this) ); + _calCollections.push_back( new CaloType(hcal, Gaudi::DataHandle::Reader, this) ); + } + return GaudiAlgorithm::initialize(); +} + +void MarlinArbor::HitsPreparation() +{ + cout<<"Start to prepare Hits"<<endl; +} + +void MarlinArbor::MakeIsoHits( std::vector<edm4hep::ConstCalorimeterHit> inputCaloHits, DataHandle<edm4hep::CalorimeterHitCollection>& m_hitcol) +{ + edm4hep::CalorimeterHitCollection* isohitcoll = m_hitcol.createAndPut(); + + int nhit = inputCaloHits.size(); + + for(int i = 0; i < nhit; i++) + { + auto a_hit = inputCaloHits[i]; + auto IsoHit = isohitcoll->create(); + IsoHit.setPosition(a_hit.getPosition()); + IsoHit.setCellID(a_hit.getCellID()); + IsoHit.setEnergy(a_hit.getEnergy()); + //isohitcoll->addElement(collhit); + } + +} + +StatusCode MarlinArbor::execute() +{ + //if(_eventNr % m_reportEvery == 0) cout<<"eventNr: "<<_eventNr<<endl; + _eventNr++; + + MarlinArbor::HitsPreparation(); //Absorb isolated hits; + + TVector3 currHitPos; + + std::vector< TVector3 > inputHitsPos; + std::vector< ArborHit > inputABHit; + std::vector< edm4hep::ConstCalorimeterHit > inputHits; + std::vector< edm4hep::ConstCalorimeterHit > inputECALHits; + std::vector< edm4hep::ConstCalorimeterHit > 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::ConstCalorimeterHit> IsoHits; + + 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; + 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() ) + SubDId = 2; + else + SubDId = 3; + + if(i1 > _EcalCalCollections.size() - 1) + KShift = 100; + else if( i1 == _calCollections.size() - 2) //HCAL Ring + KShift = 50; + + auto CaloHitColl = _calCollections[i1]->get(); + + //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 ){ + + inputECALHits.push_back(a_hit); + } + else{ + inputHCALHits.push_back(a_hit); + } + ArborHit a_abhit(currHitPos, LayerNum, 0, Depth, StaveNum, SubDId); + inputABHit.push_back(a_abhit); + inputHits.push_back(a_hit); + } + + + // cout<<i1<<" Stat "<<SubDId<<" ~~~ "<<inputABHit.size()<<endl; + + } + //cout<<"hit size"<<inputHits.size()<<endl; + + Sequence = Arbor(inputABHit, _cepc_thresholds); + + m_ArborToolLCIO->ClusterBuilding( branchCol, inputHits, Trees, 0 ); + + for(unsigned int i2 = 0; i2 < IsoHitsIndex.size(); i2++) + { + auto a_Isohit = inputHits[ IsoHitsIndex[i2] ]; + if(a_Isohit.getEnergy() > 0) //Veto Trk End Hits + { + IsoHits.push_back(a_Isohit); + } + } + + MakeIsoHits(IsoHits, m_isohitcol); + return StatusCode::SUCCESS; +} + + +StatusCode MarlinArbor::finalize() +{ + std::cout<<"Arbor Ends. Good luck"<<std::endl; + return GaudiAlgorithm::finalize(); +} diff --git a/Reconstruction/PFA/Arbor/src/MarlinArbor.hh b/Reconstruction/PFA/Arbor/src/MarlinArbor.hh new file mode 100644 index 0000000000000000000000000000000000000000..5f8a144494a3131c6bb4a4e7a732182e0b825833 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/MarlinArbor.hh @@ -0,0 +1,126 @@ +#ifndef _MarlinArbor_hh_ +#define _MarlinArbor_hh_ + +#include <string> +#include <iostream> +#include <fstream> +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "Gaudi/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/Cluster.h" +#include "edm4hep/ClusterCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +#include "Arbor.h" +#include "ArborToolLCIO.hh" + +#include <TTree.h> +#include <TFile.h> +#include <TH1.h> +#include <TH2.h> +#include <TH3.h> + +class TTree; + +//class CollectionMaps +//{ +//public: +// CollectionMaps(); +// void clear(); +// std::map<std::string, std::vector<edm4hep::MCParticle> > collectionMap_MC; +// std::map<std::string, std::vector<edm4hep::CalorimeterHit> > collectionMap_CaloHit; +//}; + +class MarlinArbor : public GaudiAlgorithm +{ + public: + + MarlinArbor(const std::string& name, ISvcLocator* svcLoc); + + + virtual StatusCode initialize() ; + + virtual StatusCode execute() ; + + virtual StatusCode finalize() ; + void HitsPreparation(); // LCEvent * evtP); + + void MakeIsoHits(std::vector<edm4hep::ConstCalorimeterHit> inputCaloHits, DataHandle<edm4hep::CalorimeterHitCollection>& m_isohitcol); + protected: + std::string _colName; + std::vector<std::string> _CalCollections; + std::vector<std::string> _SimCalCollections; + std::vector<std::string> _garlicCollections; + std::vector<std::string> _endHitTrackCollections; + std::vector<std::string> _EcalPreShowerCollections; + std::vector<std::string> _EcalCalCollections; + std::vector<std::string> _HcalCalCollections; + std::vector<float> _cepc_thresholds; + + + typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloType; + Gaudi::Property<std::vector<std::string>> m_ecalColNames{this, "ECALCollections", {"ECALBarrel", "ECALEndcap", "ECALOther"}, "Input ECAL Hits Collection Names"}; + Gaudi::Property<std::vector<std::string>> m_hcalColNames{this, "HCALCollections", {"HCALBarrel", "HCALEndcap", "HCALOther"}, "Input HCAL Hits Collection Names"}; + + + 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"}; + std::map<std::string, std::string> m_col_readout_map; + + std::vector<CaloType*> _ecalCollections; + std::vector<CaloType*> _hcalCollections; + std::vector<CaloType*> _calCollections; + + TTree *_outputTree; + std::string _treeFileName; + int _EH; + float _HitPos[3]; + float _BushP[3]; + float _CloseDis; + float _HitEnergy; + + int _CellSize; + int _CaloTrackLengthCut; + + int _Num, _Seg, _eventNr; + int numElements; + + float _DHCALFirstThreshold; + float _InitLinkDisThreshold; + + bool _DHCALSimuDigiMode; + bool _FlagInputSimHit; + bool _FlagMutePhoton; + bool _FlagMuteChargeParticle; + bool _FlagMuteGarlicHits; + 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; + std::ostream _output; + float _HLayerCut; + + SmartIF<IGeomSvc> m_geosvc; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + ArborToolLCIO * m_ArborToolLCIO; +}; + + +#endif + + diff --git a/Reconstruction/PFA/Arbor/src/cellIDDecoder.h b/Reconstruction/PFA/Arbor/src/cellIDDecoder.h new file mode 100644 index 0000000000000000000000000000000000000000..6503f352d754fdac3279dab087abd6a359994145 --- /dev/null +++ b/Reconstruction/PFA/Arbor/src/cellIDDecoder.h @@ -0,0 +1,120 @@ +#ifndef cellIDDecoder_h +#define cellIDDecoder_h 1 + +#include "EVENT/LCCollection.h" +#include "UTIL/BitField64.h" +#include "lcio.h" +#include <string> + +// fixes problem in gcc 4.0.3 +#include "EVENT/LCParameters.h" + +//##################### changed for EMD4HEP ######## + +namespace ID_UTIL{ + + + /** Convenient class for decoding cellIDs from collection parameter LCIO::CellIDEncoding. + * See UTIL::BitField64 for a description of the encoding string. + * + * @see BitField64 + * @version $Id: CellIDDecoder.h,v 1.9.16.1 2011-03-04 14:09:07 engels Exp $ + */ + template <class T> + class CellIDDecoder { + + public: + + CellIDDecoder() = default ; + CellIDDecoder(const CellIDDecoder& ) = delete ; + CellIDDecoder& operator=(const CellIDDecoder& ) = delete ; + + + /** Constructor takes encoding string as argument. + */ + CellIDDecoder( const std::string& encoder_str ) : _oldHit(0) { + + if( encoder_str.length() == 0 ){ + throw( lcio::Exception( "CellIDDecoder : string of length zero provided as encoder string" ) ) ; + } + _b = new UTIL::BitField64( encoder_str ) ; + + } + + /** Constructor reads encoding string from collection parameter LCIO::CellIDEncoding. + */ + CellIDDecoder( const EVENT::LCCollection* col ) : _oldHit(0) { + + std::string initString("") ; + + if( col !=0 ) + initString = col->getParameters().getStringVal( lcio::LCIO::CellIDEncoding ) ; + + if( initString.size() == 0 ) { + + initString = _defaultEncoding ; + + std::cout << " ----------------------------------------- " << std::endl + << " WARNING: CellIDDecoder - no CellIDEncoding parameter in collection ! " + << std::endl + << " -> using default : \"" << initString << "\"" + << std::endl + << " ------------------------------------------ " + << std::endl ; + } + + _b = new UTIL::BitField64( initString ) ; + } + + ~CellIDDecoder(){ + + delete _b ; + } + + + /** Provides access to the bit fields, e.g. <br> + * int layer = myCellIDEncoding( hit )[ "layer" ] ; + * + */ + inline const UTIL::BitField64 & operator()( const T* hit ){ + + if( hit != _oldHit && hit ) { + auto id = hit->getCellID(); + unsigned int id0 = id&0xFFFFFFFF; + unsigned int id1 = id>>32; + //lcio::long64 val = lcio::long64( hit->getCellID0() & 0xffffffff ) | ( lcio::long64( hit->getCellID1() ) << 32 ) ; + lcio::long64 val = lcio::long64( id0 & 0xffffffff ) | ( lcio::long64( id1 ) << 32 ) ; + + _b->setValue( val ) ; + + _oldHit = hit ; + } + + return *_b ; + } + + + /** This can be used to set the default encoding that is used if no + * CellIDEncoding parameter is set in the collection, e.g. in older lcio files. + */ + static void setDefaultEncoding(const std::string& defaultEncoding ) { + + _defaultEncoding = std::string( defaultEncoding ) ; + } + + protected: + UTIL::BitField64* _b{} ; + const T* _oldHit{NULL} ; + + static std::string _defaultEncoding; + } ; + + template <class T> + std::string CellIDDecoder<T>::_defaultEncoding + = std::string("byte0:8,byte1:8,byte2:8,byte3:8,byte4:8,byte5:8,byte6:8,byte7:8") ; + + +} // namespace +#endif + + diff --git a/Reconstruction/PFA/CMakeLists.txt b/Reconstruction/PFA/CMakeLists.txt index b3a1b42645d41829b472ce79e403cb99e0aa26d5..05a8afcc79cefcc218af482d3edf8e34cc1ce994 100644 --- a/Reconstruction/PFA/CMakeLists.txt +++ b/Reconstruction/PFA/CMakeLists.txt @@ -1,2 +1,3 @@ +add_subdirectory(Arbor) add_subdirectory(Pandora)