Add Fwd tracks to StMuDST (#538)

- Adds StMuFwdTrack to MuDST, - minor updates to StFwdTrackMaker (control material effects + debug tools/ttree) - Add StRoot/StFwdUtils : contains example fwd analysis maker --------- Co-authored-by: Dmitri Smirnov <[email protected]>
star-bnl · Jul 27, 2023 · c677c84 · c677c84
1 parent 7fd7f97
commit c677c84
Show file tree

Hide file tree

Showing 27 changed files with 1,434 additions and 252 deletions.
diff --git a/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx b/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx
diff --git a/StRoot/StFwdTrackMaker/StFwdTrackMaker.h b/StRoot/StFwdTrackMaker/StFwdTrackMaker.h
@@ -62,6 +62,10 @@ struct FwdTreeData {
     vector<float> fstX, fstY, fstZ;
     vector<int> fstTrackId;
 
+    int fcsN;
+    vector<float> fcsX, fcsY, fcsZ;
+    vector<int> fcsDet;
+
     // RC tracks
     int rcN;
     vector<float> rcPt, rcEta, rcPhi, rcQuality;
@@ -89,6 +93,7 @@ struct FwdTreeData {
     int thdN;
     vector<float> thdX, thdY, thaX, thaY, thaZ;
 
+    bool saveCrit = false;
     std::map<string, std::vector<float>> Crits;
     std::map<string, std::vector<int>> CritTrackIds;
 

diff --git a/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h b/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h
@@ -73,7 +73,7 @@ class GenfitTrackResult {
         this->trackSeed = seedTrack;
 
         try {
-            this->track = new genfit::Track(*track);
+            this->track = track;
             this->status = *(this->track->getFitStatus());
             this->trackRep = this->track->getCardinalRep();
 
@@ -86,9 +86,10 @@ class GenfitTrackResult {
             this->nPV = this->track->getNumPoints() - (nFTT + nFST);
 
             this->momentum = this->trackRep->getMom( this->track->getFittedState(0, this->trackRep) );
+            LOG_DEBUG << "GenfitTrackResult::set Track successful" << endm;
 
         } catch ( genfit::Exception &e ) {
-            LOG_ERROR << "CANNOT GET TRACK" << endm;
+            LOG_ERROR << "GenfitTrackResult cannot get track" << endm;
             this->track = nullptr;
             this->trackRep = nullptr;
 
@@ -109,11 +110,14 @@ class GenfitTrackResult {
     }
 
     void setFst( Seed_t &seedFst, genfit::Track *track ){
+        LOG_DEBUG << "GenfitTrackResult::setFSt" << endm;
         nFST = seedFst.size();
         fstSeed = seedFst;
 
         try {
             this->fstTrack = new genfit::Track(*track);
+            // make sure the McTrackId is set correctly 
+            this->fstTrack->setMcTrackId( this->track->getMcTrackId() );
             this->fstStatus = *(this->fstTrack->getFitStatus());
             this->fstTrackRep = this->fstTrack->getCardinalRep();
 
@@ -125,7 +129,7 @@ class GenfitTrackResult {
             this->fstMomentum = this->fstTrackRep->getMom( this->fstTrack->getFittedState(0, this->fstTrackRep) );
 
         } catch ( genfit::Exception &e ) {
-            LOG_ERROR << "CANNOT GET TRACK" << endm;
+            LOG_ERROR << "CANNOT GET FST TRACK" << endm;
             this->fstTrack = nullptr;
             this->fstTrackRep = nullptr;
 
@@ -528,6 +532,7 @@ class ForwardTrackMaker {
             if (mConfig.get<bool>("TrackFitter:refitSi", true)) {
                 addSiHitsMc();
             } else {
+                LOG_DEBUG << "Skipping FST Hits" << endm;
                 // skip Si refit
             }
             /***********************************************/
@@ -656,13 +661,13 @@ class ForwardTrackMaker {
             }
 
 
-            GenfitTrackResult gtr( track.size(), 0, track, mTrackFitter->getTrack() );
+            genTrack = new genfit::Track(*mTrackFitter->getTrack());
+            genTrack->setMcTrackId(idt);
+            GenfitTrackResult gtr( track.size(), 0, track, genTrack );
 
             // assign the fit results to be saved
             fitStatus = mTrackFitter->getStatus();
-            genTrack = new genfit::Track(*mTrackFitter->getTrack());
             trackRep = mTrackFitter->getTrackRep()->clone(); // Clone the track rep
-            genTrack->setMcTrackId(idt);
 
             if ( mGenHistograms && genTrack->getFitStatus(genTrack->getCardinalRep())->isFitConverged() && p.Perp() > 1e-3) {
                 mHist["FitStatus"]->Fill("GoodCardinal", 1);
@@ -679,6 +684,7 @@ class ForwardTrackMaker {
 
             mTrackResults.push_back( gtr );
 
+            LOG_DEBUG << "FwdTracker::fitTrack complete" << endm;
         } // if (mDoTrackFitting && !bailout)
     }
 
@@ -774,6 +780,9 @@ class ForwardTrackMaker {
      */
     vector<Seed_t> doTrackingOnHitmapSubset( size_t iIteration, FwdDataSource::HitMap_t &hitmap  ) {
         long long itStart = FwdTrackerUtils::nowNanoSecond();
+
+        std::vector<Seed_t> acceptedTracks;
+        std::vector<Seed_t> rejectedTracks;
         /*************************************************************/
         // Step 2
         // build 2-hit segments (setup parent child relationships)
@@ -812,6 +821,11 @@ class ForwardTrackMaker {
         LOG_DEBUG << TString::Format( "nSegments=%lu", automaton.getSegments().size() ).Data() << endm;
         LOG_DEBUG << TString::Format( "nConnections=%u", automaton.getNumberOfConnections() ).Data() << endm;
 
+        if (automaton.getNumberOfConnections() > 900 ){
+            LOG_ERROR << "Got too many connections, bailing out of tracking" << endm;
+            return acceptedTracks;
+        }
+
         // at any point we can get a list of tracks out like this:
         // std::vector < std::vector< KiTrack::IHit* > > tracks = automaton.getTracks();
         // we can apply an optional parameter <nHits> to only get tracks with >=nHits in them
@@ -853,7 +867,7 @@ class ForwardTrackMaker {
         duration = (FwdTrackerUtils::nowNanoSecond() - itStart) * 1e-6; // milliseconds
         if (mGenHistograms)
             mHist["Step3Duration"]->Fill( duration );
-        if (duration > 200){
+        if (duration > 200 || automaton.getNumberOfConnections() > 900){
             LOG_WARN << "The Three Hit Criteria took more than 200ms to process, duration: " << duration << " ms" << endm;
             LOG_WARN << "bailing out (skipping subset HNN)" << endm;
             std::vector<Seed_t> acceptedTracks;
@@ -878,8 +892,6 @@ class ForwardTrackMaker {
 
         //  only for debug really
         bool findSubsets = mConfig.get<bool>(subsetPath + ":active", true);
-        std::vector<Seed_t> acceptedTracks;
-        std::vector<Seed_t> rejectedTracks;
 
         if (findSubsets) {
             size_t minHitsOnTrack = mConfig.get<size_t>(subsetPath + ":min-hits-on-track", 7);
@@ -1023,13 +1035,14 @@ class ForwardTrackMaker {
     void addSiHitsMc() {
         FwdDataSource::HitMap_t hitmap = mDataSource->getFstHits();
 
-        for (size_t i = 0; i < mGlobalTracks.size(); i++) {
-
-            if (mGlobalTracks[i]->getFitStatus(mGlobalTracks[i]->getCardinalRep())->isFitConverged() == false || mFitMoms[i].Perp() < 1e-3) {
+        for (size_t i = 0; i < mTrackResults.size(); i++) {
+            GenfitTrackResult &gtr = mTrackResults[i];
+
+            if ( gtr.status.isFitConverged() == false || gtr.momentum.Perp() < 1e-3) {
+                LOG_DEBUG << "Skipping addSiHitsMc, fit failed" << endm;
                 return;
             }
 
-
             if ( mGenHistograms){
                 mHist["FitStatus"]->Fill("PossibleReFit", 1);
             }
@@ -1038,7 +1051,7 @@ class ForwardTrackMaker {
 
             for (size_t j = 0; j < 3; j++) {
                 for (auto h0 : hitmap[j]) {
-                    if (dynamic_cast<FwdHit *>(h0)->_tid == mGlobalTracks[i]->getMcTrackId()) {
+                    if (dynamic_cast<FwdHit *>(h0)->_tid == gtr.track->getMcTrackId()) {
                         si_hits_for_this_track[j] = h0;
                         break;
                     }
@@ -1049,6 +1062,7 @@ class ForwardTrackMaker {
             if ( si_hits_for_this_track[0] != nullptr ) nSiHitsFound++;
             if ( si_hits_for_this_track[1] != nullptr ) nSiHitsFound++;
             if ( si_hits_for_this_track[2] != nullptr ) nSiHitsFound++;
+            LOG_DEBUG << "Found " << nSiHitsFound << " FST Hits on this track (MC lookup)" << endm;
 
             if ( mGenHistograms ){
                 this->mHist[ "nSiHitsFound" ]->Fill( 1, ( si_hits_for_this_track[0] != nullptr ? 1 : 0 ) );
@@ -1060,14 +1074,15 @@ class ForwardTrackMaker {
                 if ( mGenHistograms ){
                     mHist["FitStatus"]->Fill("AttemptReFit", 1);
                 }
-                TVector3 p = mTrackFitter->refitTrackWithSiHits(mGlobalTracks[i], si_hits_for_this_track);
+                TVector3 p = mTrackFitter->refitTrackWithSiHits(gtr.track, si_hits_for_this_track);
 
                 if ( mGenHistograms ){
                     if (p.Perp() == mFitMoms[i].Perp()) {
                         mHist["FitStatus"]->Fill("BadReFit", 1);
-
+                        LOG_DEBUG << "refitTrackWithSiHits failed refit" << endm;
                     } else {
                         mHist["FitStatus"]->Fill("GoodReFit", 1);
+                        gtr.setFst( si_hits_for_this_track, mTrackFitter->getTrack() );
                     }
                 }
 

diff --git a/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h b/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h
@@ -74,8 +74,23 @@ class TrackFitter {
         gMan = gGeoManager;
         // Set up the material interface and set material effects on/off from the config
         genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
-        genfit::MaterialEffects::getInstance()->setNoEffects(mConfig.get<bool>("TrackFitter::noMaterialEffects", false)); // false means defaul ON
 
+        // Set Material Stepper debug level
+        genfit::MaterialEffects::getInstance()->setDebugLvl( mConfig.get<int>("TrackFitter.MaterialEffects:DebugLvl", 0) );
+
+        genfit::MaterialEffects::getInstance()->setEnergyLossBetheBloch( mConfig.get<int>("TrackFitter.MaterialEffects.EnergyLossBetheBloch", true) );
+        genfit::MaterialEffects::getInstance()->setNoiseBetheBloch( mConfig.get<int>("TrackFitter.MaterialEffects.NoiseBetheBloch", true) );
+        genfit::MaterialEffects::getInstance()->setNoiseCoulomb( mConfig.get<int>("TrackFitter.MaterialEffects.NoiseCoulomb", true) );
+        genfit::MaterialEffects::getInstance()->setEnergyLossBrems( mConfig.get<int>("TrackFitter.MaterialEffects.EnergyLossBrems", true) );
+        genfit::MaterialEffects::getInstance()->setNoiseBrems( mConfig.get<int>("TrackFitter.MaterialEffects.NoiseBrems", true) );
+        genfit::MaterialEffects::getInstance()->ignoreBoundariesBetweenEqualMaterials( mConfig.get<int>("TrackFitter.MaterialEffects.ignoreBoundariesBetweenEqualMaterials", true) );
+
+        // do this last to override
+        genfit::MaterialEffects::getInstance()->setNoEffects( !mConfig.get<bool>("TrackFitter:MaterialEffects", false)); // negated, true means defaul is all effects on (noEffects off)
+        if (!mConfig.get<bool>("TrackFitter:MaterialEffects", false)){
+            LOG_INFO << "Turning OFF GenFit Material Effects in stepper" << endm;
+        }
+
         // Determine which Magnetic field to use
         // Either constant field or real field from StarFieldAdaptor
         if (mConfig.get<bool>("TrackFitter:constB", false)) {
@@ -171,7 +186,7 @@ class TrackFitter {
         // Now load FTT
         // mConfig.getVector<>(...) requires a default, hence the 
         mFTTZLocations = fwdGeoUtils.fttZ(
-            mConfig.getVector<double>("TrackFitter.Geometry:ftt", {0.0f, 0.0f, 0.0f, 0.0f})
+            mConfig.getVector<double>("TrackFitter.Geometry:ftt", {281.082,304.062,325.058,348.068})
             );
 
         if ( fwdGeoUtils.fttZ( 0 ) < 1.0 ) { // returns 0.0 on failure
@@ -772,8 +787,8 @@ class TrackFitter {
         }
 
         // create the track representations
-        auto trackRepPos = new genfit::RKTrackRep(mPdgPositron);
-        auto trackRepNeg = new genfit::RKTrackRep(mPdgElectron);
+        auto trackRepPos = new genfit::RKTrackRep(mPdgMuon);
+        auto trackRepNeg = new genfit::RKTrackRep(mPdgAntiMuon);
 
         // Create the track
         mFitTrack = new genfit::Track(trackRepPos, seedPos, seedMom);
@@ -950,6 +965,8 @@ class TrackFitter {
     const int mPdgPiMinus = -211;
     const int mPdgPositron = 11;
     const int mPdgElectron = -11;
+    const int mPdgMuon = 13;
+    const int mPdgAntiMuon = -13;
 
 
     // det z locations loaded from geom or config

diff --git a/StRoot/StFwdTrackMaker/macro/sim/fast_track.xml b/StRoot/StFwdTrackMaker/macro/sim/fast_track.xml
@@ -1,11 +1,19 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <config>
     <Output url="fast_track.root" />
-	<Geometry>fGeom.root</Geometry>
-    <Source ftt="GEANT" fttFilter="true" />
+    <Source ftt="GEANT" fttFilter="false" />
+
     <!-- Missing <TrackFinder> means use MC Track finding -->
 
-	<TrackFitter refitSi="true" mcSeed="true" refitGBL="true" >
-        <Vertex sigmaXY="0.002" sigmaZ="1.0" includeInFit="true" smearMcVertex="true" />
+	<TrackFitter refitSi="true" mcSeed="true" refitGBL="false" MaterialEffects="false" >
+        <MaterialEffects DebugLvl="0" >
+            <EnergyLossBetheBloch>false</EnergyLossBetheBloch>
+            <NoiseBetheBloch>false</NoiseBetheBloch>
+            <NoiseCoulomb>false</NoiseCoulomb>
+            <EnergyLossBrems>false</EnergyLossBrems>
+            <NoiseBrems>false</NoiseBrems>
+            <ignoreBoundariesBetweenEqualMaterials>true</ignoreBoundariesBetweenEqualMaterials>
+        </MaterialEffects>
+        <Vertex sigmaXY="0.001" sigmaZ="0.1" includeInFit="true" smearMcVertex="false" />
     </TrackFitter>
 </config>
diff --git a/StRoot/StFwdTrackMaker/macro/sim/full_track.xml b/StRoot/StFwdTrackMaker/macro/sim/full_track.xml
@@ -1,7 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <config>
     <Output url="full_track.root" />
-	<Geometry>fGeom.root</Geometry>
+	<!-- <Geometry>fGeom.root</Geometry> -->
     <Source ftt="GEANT" fttFilter="true" />
 
     <TrackFinder nIterations="1">
@@ -33,7 +33,7 @@
         <HitRemover active="true" />
     </TrackFinder>
 
-	<TrackFitter refitSi="true" mcSeed="true" >
-        <Vertex sigmaXY="0.002" sigmaZ="1.0" includeInFit="false" smearMcVertex="true" />
+	<TrackFitter refitSi="true" mcSeed="false" MaterialEffects="false" >
+        <Vertex sigmaXY="1" sigmaZ="70.0" includeInFit="true" smearMcVertex="true" />
     </TrackFitter>
 </config>