From d9de5556ebe593a5c6e2764d052d3a861016cafb Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Tue, 31 Aug 2021 10:20:39 +0200 Subject: [PATCH] TPC parameters changed to mm in KinkFinder --- Tracking/KinkFinder/src/KinkFinder.cc | 170 +++++++++++++------------- 1 file changed, 84 insertions(+), 86 deletions(-) diff --git a/Tracking/KinkFinder/src/KinkFinder.cc b/Tracking/KinkFinder/src/KinkFinder.cc index 2d7186cb..74cfcb57 100644 --- a/Tracking/KinkFinder/src/KinkFinder.cc +++ b/Tracking/KinkFinder/src/KinkFinder.cc @@ -16,6 +16,7 @@ #include #include "HelixClass.h" +#include "GeometryUtil.h" using namespace lcio ; using namespace marlin ; @@ -41,7 +42,7 @@ KinkFinder aKinkFinder ; KinkFinder::KinkFinder() : Processor("KinkFinder") { _description = "Kink Finder Processor : finds kinks, prongs and splits" ; - + registerInputCollection(LCIO::TRACK, "TrackCollection", "Name of input collection of reconstructed particles", @@ -188,33 +189,31 @@ void KinkFinder::init() { dd4hep::Detector& theDet = dd4hep::Detector::getInstance(); - double bfieldV[3] ; - theDet.field().magneticField( { 0., 0., 0. } , bfieldV ) ; - _bField = bfieldV[2]/dd4hep::tesla ; - + _bField = MarlinUtil::getBzAtOrigin() ; + dd4hep::DetElement tpcDE = theDet.detector("TPC") ; const dd4hep::rec::FixedPadSizeTPCData* tpc = tpcDE.extension() ; - _tpcInnerR = tpc->rMinReadout; - _tpcOuterR = tpc->rMaxReadout; - _tpcZmax = tpc->driftLength; + _tpcInnerR = tpc->rMinReadout/dd4hep::mm; + _tpcOuterR = tpc->rMaxReadout/dd4hep::mm; + _tpcZmax = tpc->driftLength/dd4hep::mm; _tpcMaxRow = tpc->maxRow; dd4hep::DetElement vtxDE = theDet.detector("VXD"); dd4hep::rec::ZPlanarData* vtx = vtxDE.extension(); - _nLayersVTX=vtx->layers.size(); + _nLayersVTX=vtx->layers.size(); for(int iL=0;iL<_nLayersVTX;iL++){ - _rVTX.push_back( vtx->layers[iL].distanceSensitive ); + _rVTX.push_back( (vtx->layers[iL].distanceSensitive)/dd4hep::mm ); } dd4hep::DetElement sitDE = theDet.detector("SIT"); dd4hep::rec::ZPlanarData* sit = sitDE.extension(); - _nLayersSIT=sit->layers.size(); + _nLayersSIT=sit->layers.size(); for(int iL=0;iL<_nLayersSIT;iL++){ - _rSIT.push_back( sit->layers[iL].distanceSensitive ); + _rSIT.push_back( (sit->layers[iL].distanceSensitive)/dd4hep::mm ); } @@ -224,14 +223,14 @@ void KinkFinder::init() { } -void KinkFinder::processRunHeader( LCRunHeader* ) { +void KinkFinder::processRunHeader( LCRunHeader* ) { _nRun++ ; _nEvt = 0; -} +} -void KinkFinder::processEvent( LCEvent * evt ) { +void KinkFinder::processEvent( LCEvent * evt ) { TrackVec tracks; @@ -242,8 +241,8 @@ void KinkFinder::processEvent( LCEvent * evt ) { LCCollection* relationTrackCollection = NULL; LCRelationNavigator* trackNavigator=NULL; try { - relationTrackCollection = evt->getCollection("LDCTracksMCP"); - trackNavigator = new LCRelationNavigator(relationTrackCollection); + relationTrackCollection = evt->getCollection("LDCTracksMCP"); + trackNavigator = new LCRelationNavigator(relationTrackCollection); } catch(DataNotAvailableException &e){ } @@ -354,13 +353,13 @@ void KinkFinder::processEvent( LCEvent * evt ) { zBegin = zmax; zEnd = zmin; } - float signPz = zEnd - zBegin; + float signPz = zEnd - zBegin; zAtEnd[itrack] = zEnd; zAtStart[itrack] = zBegin; int _nHitsInFit = 50; int nhitsFit; - if (nhits > _nHitsInFit) { + if (nhits > _nHitsInFit) { for (int iz = 0 ; iz < nhits-1; iz++) for (int jz = 0; jz < nhits-iz-1; jz++){ TrackerHit * one = hitvec[jz]; @@ -385,9 +384,9 @@ void KinkFinder::processEvent( LCEvent * evt ) { float omega = tracks[itrack]->getOmega(); float z0 = tracks[itrack]->getZ0(); float tanlambda = tracks[itrack]->getTanLambda(); - helix.Initialize_Canonical(phi0, d0, z0, omega, tanlambda, _bField); - float totMomentum =0; - float Mom[3]; + helix.Initialize_Canonical(phi0, d0, z0, omega, tanlambda, _bField); + float totMomentum =0; + float Mom[3]; for (int j=0; j < 3; ++j) { Mom[j] = helix.getMomentum()[j]; totMomentum += Mom[j]*Mom[j]; @@ -425,13 +424,13 @@ void KinkFinder::processEvent( LCEvent * evt ) { float * xh = new float[nhitsFit]; float * yh = new float[nhitsFit]; float * zh = new float[nhitsFit]; - float * ah = new float[nhitsFit]; + float * ah = new float[nhitsFit]; for (int i=0; igetPosition()[0]); yh[i] = float(hitvec[i]->getPosition()[1]); zh[i] = float(hitvec[i]->getPosition()[2]); - } + } ClusterShapes * shapesS = new ClusterShapes(nhitsFit,ah,xh,yh,zh); float par[5]; @@ -445,9 +444,9 @@ void KinkFinder::processEvent( LCEvent * evt ) { float y0 = par[1]; float r0 = par[2]; float bz = par[3]; - float phiH = par[4]; + float phiH = par[4]; helixStart[itrack] = new HelixClass(); - helixStart[itrack]->Initialize_BZ(x0, y0, r0, + helixStart[itrack]->Initialize_BZ(x0, y0, r0, bz, phiH, _bField,signPz, zBegin); @@ -464,7 +463,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { xh[i] = float(hitvec[nhits-1-i]->getPosition()[0]); yh[i] = float(hitvec[nhits-1-i]->getPosition()[1]); zh[i] = float(hitvec[nhits-1-i]->getPosition()[2]); - } + } ClusterShapes * shapesE = new ClusterShapes(nhitsFit,ah,xh,yh,zh); shapesE->FitHelix(500, 0, 1, par, dpar, chi2, distmax); delete shapesE; @@ -472,9 +471,9 @@ void KinkFinder::processEvent( LCEvent * evt ) { y0 = par[1]; r0 = par[2]; bz = par[3]; - phiH = par[4]; + phiH = par[4]; helixEnd[itrack] = new HelixClass(); - helixEnd[itrack]->Initialize_BZ(x0, y0, r0, + helixEnd[itrack]->Initialize_BZ(x0, y0, r0, bz, phiH, _bField,signPz, zEnd); @@ -485,7 +484,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { } - for(unsigned int i=0;i< tracks.size();++i){ + for(unsigned int i=0;i< tracks.size();++i){ // Track* tracki = tracks[i]; // float d0i = tracki->getD0(); // float z0i = tracki->getZ0(); @@ -543,7 +542,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { refs[0] = helixEnd[j]->getReferencePoint()[0]; refs[1] = helixEnd[j]->getReferencePoint()[1]; refs[2] = helixEnd[j]->getReferencePoint()[2]; - helixEnd[j]->getPointInZ(z, refs, seedj); + helixEnd[j]->getPointInZ(z, refs, seedj); deltaz = zAtEnd[i] - zAtEnd[j]; ddx = (zout[i].x-zout[j].x); ddy = (zout[i].y-zout[j].y); @@ -557,7 +556,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { float dy = seedi[1]-seedj[1]; float dz = seedi[2]-seedj[2]; float dr = sqrt(dx*dx+dy*dy+dz*dz); - + float xkink = 0; float ykink = 0; float zkink = 0; @@ -586,8 +585,8 @@ void KinkFinder::processEvent( LCEvent * evt ) { } } } - - + + if(dr<100 || mcKink){ dr = 999999.; @@ -609,7 +608,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { } float zstep = (ze-zs)/10.; if(zstep<1)zstep=1; - + for(float zf = zs; zf < ze; zf+= zstep){ helixEnd[ip]->getPointInZ(zf, refsp, seedp); helixStart[id]->getPointInZ(zf, refsd, seedd); @@ -628,11 +627,11 @@ void KinkFinder::processEvent( LCEvent * evt ) { } } } - + bool ok = true; if(fabs(deltaz)>200)ok=false; if(fabs(deltaz)>100 && dr > 5.0)ok=false; - + float deltaRxyCut = -100; float drCut = -100; // bool goodRadialSep = false; @@ -672,13 +671,13 @@ void KinkFinder::processEvent( LCEvent * evt ) { if(fix>deltaRxyCut)deltaRxyCut = fix; } // add some slop to account for geometry - deltaRxyCut = deltaRxyCut*1.1; + deltaRxyCut = deltaRxyCut*1.1; } } - + // looser cuts for debug - // std::cout << i << " : " << j << " dr = " << dr << " ( " << drCut << " ) deltaRxy = " << deltarxy << " ( " << deltaRxyCut << " ) " << std::endl; + // std::cout << i << " : " << j << " dr = " << dr << " ( " << drCut << " ) deltaRxy = " << deltarxy << " ( " << deltaRxyCut << " ) " << std::endl; if( (dr0){ - float deltaPion = fabs(massMuNu-mPion)/_pionDecayMassCut; + float deltaPion = fabs(massMuNu-mPion)/_pionDecayMassCut; if(deltaPion<1 && tPion > 0.001)probPion = 3.125*deltaPion*deltaPion+tPion; } if(_kaonDecayMassCut>0){ - float deltaKaonMuNu = fabs(massMuNu-mKaon)/_kaonDecayMassCut; + float deltaKaonMuNu = fabs(massMuNu-mKaon)/_kaonDecayMassCut; float deltaKaonPiPi = fabs(massPiPi-mKaon)/_kaonDecayMassCut; float deltaKaon = deltaKaonMuNu; if(deltaKaonPiPi 0.005)probKaon = 3.125*deltaKaon*deltaKaon+tKaon; } if(_sigmaDecayMassCut>0 && momentum[i] > _minELambda){ - float deltaSigmaPiN = fabs(massPiN-mSigma)/_sigmaDecayMassCut; + float deltaSigmaPiN = fabs(massPiN-mSigma)/_sigmaDecayMassCut; float deltaSigmaPPi0 = fabs(massPPi0-mSigma)/_sigmaDecayMassCut; float deltaSigma = deltaSigmaPiN; if(deltaSigmaPPi00 && momentum[i] > _minELambda){ - float deltaHyperon = fabs(massPiL-mHyperon)/_hyperonDecayMassCut; + float deltaHyperon = fabs(massPiL-mHyperon)/_hyperonDecayMassCut; if(deltaHyperon<1 && tHyperon < _hyperonTimeCut)probHyperon = 3.125*deltaHyperon*deltaHyperon+tHyperon; } bool goodKinkMass = false; if(probPion>0.0001||probKaon>0.0001||probSigma>0.0001||probHyperon>0.0001)goodKinkMass=true; - float dpop = 2*fabs(momentum[i]-momentum[j])/(momentum[i]+momentum[j]); + float dpop = 2*fabs(momentum[i]-momentum[j])/(momentum[i]+momentum[j]); float dp = fabs(momentum[i]-momentum[j]); if(dr0){ - std::cout << " CAND SPLIT I : " << nhitsi << " ntpc " << ntpci << " nclose " << nclosei << " max " << maxdisti << " fclose : " << fclosei << std::endl; - std::cout << " CAND SPLIT J : " << nhitsj << " ntpc " << ntpcj << " nclose " << nclosej << " max " << maxdistj << " fclose : " << fclosej << std::endl; + std::cout << " CAND SPLIT I : " << nhitsi << " ntpc " << ntpci << " nclose " << nclosei << " max " << maxdisti << " fclose : " << fclosei << std::endl; + std::cout << " CAND SPLIT J : " << nhitsj << " ntpc " << ntpcj << " nclose " << nclosej << " max " << maxdistj << " fclose : " << fclosej << std::endl; } if(maxdistj<50 && maxdisti < 50 && fclosej > 0.95 && fclosei > 0.95 && ntpcj+ntpci < _tpcMaxRow+10.)split = true; splitDaughters[i].push_back(kinkij); @@ -820,28 +819,28 @@ void KinkFinder::processEvent( LCEvent * evt ) { if(charge[i]*charge[j]>0 && goodKinkMass){ if(_debugPrinting>0)std::cout << "Found kink candidate " << i << " " << j << std::endl; if(probPion>0||probKaon>0||probSigma>0||probHyperon>0)goodKinkMass=true; - if(probPion > probKaon && + if(probPion > probKaon && probPion > probSigma && probPion > probHyperon){ kinkij.pdgCode = 211; if(charge[i]<0)kinkij.pdgCode = -211; kinkij.mass = mPion; } - if(probKaon > probPion && + if(probKaon > probPion && probKaon > probSigma && probKaon > probHyperon){ kinkij.mass = mKaon; kinkij.pdgCode = 321; if(charge[i]<0)kinkij.pdgCode = -321; } - if(probSigma > probPion && + if(probSigma > probPion && probSigma > probKaon && probSigma > probHyperon){ kinkij.mass = mSigma; kinkij.pdgCode = 3222; if(charge[i]<0)kinkij.pdgCode = 3222; } - if(probHyperon > probPion && + if(probHyperon > probPion && probHyperon > probKaon && probHyperon > probSigma){ kinkij.mass = mHyperon; @@ -854,7 +853,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { // prongs if(deltarxy0)std::cout << "Found prong candidate " << i << " " << j << std::endl; prongDaughters[i].push_back(kinkij); kinkij.pdgCode = 211; @@ -863,7 +862,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { } } - + if(_debugPrinting>0){ if(mcParticle[j]!=NULL){ EVENT::MCParticleVec parents = mcParticle[j]->getParents(); @@ -873,24 +872,24 @@ void KinkFinder::processEvent( LCEvent * evt ) { } } } - + std::cout << " Candidate kink tracks : " << i << "," << j << " p[i] " << momentum[i] << " p[j] " << momentum[j] << std::endl; if(flipped)std::cout << " Flipped Track " << std::endl;; std::cout << " MC : "; if(mcParticle[i]!=NULL)std::cout << mcParticle[i]->getPDG(); - std::cout << " - "; + std::cout << " - "; if(mcParticle[j]!=NULL)std::cout << mcParticle[j]->getPDG(); - std::cout << std::endl; - std::cout << " Mass : " << massENu << " " << massMuNu << " " << massPiPi << " " << massPiN << " " << massPPi0 << " " << massPiL << std::endl; - if(flipped)std::cout << " MassF : " << FmassENu << " " << FmassMuNu << " " << FmassPiPi << " " << FmassPiN << " " << FmassPPi0 << FmassPiL << std::endl; + std::cout << std::endl; + std::cout << " Mass : " << massENu << " " << massMuNu << " " << massPiPi << " " << massPiN << " " << massPPi0 << " " << massPiL << std::endl; + if(flipped)std::cout << " MassF : " << FmassENu << " " << FmassMuNu << " " << FmassPiPi << " " << FmassPiN << " " << FmassPPi0 << FmassPiL << std::endl; std::cout << " ProbPion : " << probPion << " TimePion : " << tPion << " Dm = " << massMuNu-mPion << std::endl; - std::cout << " ProbKaon : " << probKaon << " TimeKaon : " << tKaon << - " Dm = " << massMuNu-mKaon << + std::cout << " ProbKaon : " << probKaon << " TimeKaon : " << tKaon << + " Dm = " << massMuNu-mKaon << " Dm = " << massPiPi-mKaon << std::endl; - std::cout << " ProbSigma : " << probSigma << " TimeSigma : " << tSigma << " Dm = " << massPiN-mSigma << + std::cout << " ProbSigma : " << probSigma << " TimeSigma : " << tSigma << " Dm = " << massPiN-mSigma << " Dm = " << massPPi0-mSigma << std::endl; std::cout << " ProbHyperon : " << probHyperon << " TimeHyperon : " << tHyperon << " Dm = " << massPiL-mHyperon << std::endl; std::cout << " Charge : " << i << "," << j << " q[i] " << charge[i] << " q[j] " << charge[j] << std::endl; @@ -913,7 +912,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { int countk = 0; int countp = 0; - for(unsigned int i=0;i< tracks.size();++i){ + for(unsigned int i=0;i< tracks.size();++i){ if(kinkDaughters[i].size()>0 || prongDaughters[i].size()>0){ @@ -931,7 +930,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { LCCollectionVec * colSplitRecoPart = NULL; LCCollectionVec * colSplitVertex = NULL; - for(unsigned int i=0;i< tracks.size();++i){ + for(unsigned int i=0;i< tracks.size();++i){ if(kinkDaughters[i].size()>0 || prongDaughters[i].size()>0 || splitDaughters[i].size()>0){ if(_debugPrinting>0){ if(kinkDaughters[i].size()>0){ @@ -963,16 +962,16 @@ void KinkFinder::processEvent( LCEvent * evt ) { vertex[iC] = splitDaughters[i][0].vtx[iC]; mom[iC] = splitDaughters[i][0].p[iC]; } - float distance = splitDaughters[i][0].distance; + float distance = splitDaughters[i][0].distance; vtx->setPosition( vertex ); vtx->addParameter( distance ); part->setMomentum( mom ); part->setType( code ); - + if(_debugPrinting>0){ std::cout << " Code = " << code << " Distance = " << distance << std::endl; - std::cout << " Vertex = (" - << vertex[0] << "," + std::cout << " Vertex = (" + << vertex[0] << "," << vertex[1] << "," << vertex[2] << ")" << std::endl; std::cout << " Momentum = (" @@ -981,7 +980,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { << mom[2] << ")" << std::endl; } float mass = splitDaughters[i][0].mass; - part->setMass( mass ); + part->setMass( mass ); vtx->setAssociatedParticle( part ); part->setStartVertex( vtx ); part->addTrack( tracks[splitDaughters[i][0].tracki] ); @@ -1007,15 +1006,15 @@ void KinkFinder::processEvent( LCEvent * evt ) { vertex[iC] = kinkDaughters[i][0].vtx[iC]; mom[iC] = kinkDaughters[i][0].p[iC]; } - float distance = kinkDaughters[i][0].distance; + float distance = kinkDaughters[i][0].distance; vtx->setPosition( vertex ); vtx->addParameter( distance ); part->setMomentum( mom ); part->setType( code ); if(_debugPrinting>0){ std::cout << " Code = " << code << " Distance = " << distance << std::endl; - std::cout << " Vertex = (" - << vertex[0] << "," + std::cout << " Vertex = (" + << vertex[0] << "," << vertex[1] << "," << vertex[2] << ")" << std::endl; std::cout << " Momentum = (" @@ -1024,7 +1023,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { << mom[2] << ")" << std::endl; } float mass = kinkDaughters[i][0].mass; - part->setMass( mass ); + part->setMass( mass ); vtx->setAssociatedParticle( part ); part->setStartVertex( vtx ); part->addTrack( tracks[kinkDaughters[i][0].tracki] ); @@ -1049,15 +1048,15 @@ void KinkFinder::processEvent( LCEvent * evt ) { vertex[iC] = prongDaughters[i][0].vtx[iC]; mom[iC] = prongDaughters[i][0].p[iC]; } - float distance = prongDaughters[i][0].distance; + float distance = prongDaughters[i][0].distance; vtx->setPosition( vertex ); vtx->addParameter( distance ); part->setMomentum( mom ); part->setType( code ); if(_debugPrinting>0){ std::cout << " Code = " << code << " Distance = " << distance << std::endl; - std::cout << " Vertex = (" - << vertex[0] << "," + std::cout << " Vertex = (" + << vertex[0] << "," << vertex[1] << "," << vertex[2] << ")" << std::endl; std::cout << " Momentum = (" @@ -1066,7 +1065,7 @@ void KinkFinder::processEvent( LCEvent * evt ) { << mom[2] << ")" << std::endl; } float mass = prongDaughters[i][0].mass; - part->setMass( mass ); + part->setMass( mass ); vtx->setAssociatedParticle( part ); part->setStartVertex( vtx ); part->addTrack( tracks[prongDaughters[i][0].tracki] ); @@ -1082,14 +1081,14 @@ void KinkFinder::processEvent( LCEvent * evt ) { // trackUsed[secondTrack] = 1; } } - + } if(colKinkRecoPart!=NULL){ evt->addCollection(colKinkRecoPart, _kinkRecoPartColName.c_str() ); evt->addCollection(colKinkVertex, _kinkVertexColName.c_str() ); - } + } if(colProngRecoPart!=NULL){ evt->addCollection(colProngRecoPart, _prongRecoPartColName.c_str() ); evt->addCollection(colProngVertex, _prongVertexColName.c_str() ); @@ -1111,14 +1110,14 @@ void KinkFinder::processEvent( LCEvent * evt ) { void KinkFinder::check( LCEvent * ) { } - -void KinkFinder::end(){ } + +void KinkFinder::end(){ } void KinkFinder::Sorting( TrackPairVec & trkPairVec ) { int sizeOfVector = int(trkPairVec.size()); TrackPair *one,*two,*Temp; - + for (int i = 0 ; i < sizeOfVector-1; i++) for (int j = 0; j < sizeOfVector-i-1; j++) { @@ -1130,14 +1129,14 @@ void KinkFinder::Sorting( TrackPairVec & trkPairVec ) { trkPairVec[j] = trkPairVec[j+1]; trkPairVec[j+1] = Temp; } - } + } } float KinkFinder::kinkMass(HelixClass* parent, HelixClass* daughter, float daughterMass, float neutralMass){ // Calculate the invariant mass for a decaying charged particle - + float parentMom[4]; parentMom[0] = parent->getMomentum()[0]; parentMom[1] = parent->getMomentum()[1]; @@ -1153,11 +1152,10 @@ float KinkFinder::kinkMass(HelixClass* parent, HelixClass* daughter, float daugh pnu[0] = parentMom[0]-daughterMom[0]; pnu[1] = parentMom[1]-daughterMom[1]; pnu[2] = parentMom[2]-daughterMom[2]; - + float Enu = sqrt(pnu[0]*pnu[0]+pnu[1]*pnu[1]+pnu[2]*pnu[2]+neutralMass*neutralMass); float mx = sqrt((daughterE+Enu)*(daughterE+Enu)-parentMom[0]*parentMom[0]- parentMom[1]*parentMom[1]-parentMom[2]*parentMom[2]); return mx; } -