Version of the tracking algorithm for March 2022 White Paper

dq_install.sh

@@ -0,0 +1,45 @@
+#!/bin/bash
+
+#sudo mount -t cvmfs seaquest.opensciencegrid.org /cvmfs/seaquest.opensciencegrid.org
+#sudo yum install sqlite-devel
+#source  /cvmfs/seaquest.opensciencegrid.org/seaquest/software/current/setup.sh
+
+#git clone https://github.com/cmantill/e1039-core.git
+git clone https://github.com/wpmccormack/e1039-core.git
+#git clone git@github.com:wpmccormack/e1039-core.git
+cd e1039-core/
+git checkout patrick_new_tracking_withEmbedding
+cd script/
+
+cp /home/submit/pcharris/DQ/setup-install.sh
+cd ..
+./script/setup-install.sh auto
+source ../core-inst/this-e1039.sh
+./build.sh
+
+#git clone https://github.com/DarkQuest-FNAL/DarkQuest.git
+git clone https://github.com/wpmccormack/DarkQuest.git
+#git clone git@github.com:wpmccormack/DarkQuest.git
+cd DarkQuest
+#git checkout add_ana_flag_withEmbedding
+cd e1039-analysis/SimHits/
+sed "s@XXX@$PWD@g" /home/submit/pcharris/DQ/setup_mye1039_tmp.sh > setup_mye1039.sh
+#sed "s@XXX@$PWD@g" /data/t3home000/pharris/DQ/DarkQuest/e1039-analysis/SimHits/setup_mye1039_tmp.sh > setup_mye1039.sh
+
+source setup_mye1039.sh
+mkdir work
+mkdir install
+cd work
+cmake ../src/ -DCMAKE_INSTALL_PREFIX=../install
+make clean
+make
+make install
+cd ../../
+pwd
+
+cd HitEmbedding
+sed "s@XXX@$PWD@g" /work/submit/wmccorma/DQ_Embed/setup_my1039_tmpHitEmbed.sh > setup_my1039.sh
+source setup_my1039.sh
+cmake-this
+make-this
+cd ..

framework/phool/recoConsts.cc

@@ -52,10 +52,13 @@ void recoConsts::set_defaults()
   set_BoolFlag("MC_MODE", false);
   set_BoolFlag("COSMIC_MODE", false);
 
+  set_BoolFlag("NOT_DISPLACED", true);
   set_BoolFlag("TARGETONLY", false);
   set_BoolFlag("DUMPONLY", false);
   set_BoolFlag("TRACK_ELECTRONS", false); /**This flag is used to turn off hit requirements in proportional tubes and tracking layers after the DarkQuest ECal*/
   set_BoolFlag("TRACK_DISPLACED", false); /**When this flag is turned on, an expected-position window is not used when connecting station2+3 tracklets to station 1 tracklets, AND both charges are checked for every station1+2+3 tracklet.  The charge that yields a lower chi2 in the tracklet fit is used.  This is needed for displaced tracks, where the tracklet's charge cannot be determined from x0 position alone.*/
+  set_BoolFlag("OLD_TRACKING", false); /**This flag can be used to turn on the old version of tracking, where tracklets in station 2 and 3 are matched based on the minimization of the sum of the DCA^2 <- that's distance of closest approach*/
+
 
   //Following values are fed to GeomSvc
   set_BoolFlag("OnlineAlignment", false);

packages/calibrator/CalibDriftDist.cc

@@ -75,9 +75,14 @@ int CalibDriftDist::process_event(PHCompositeNode* topNode)
     int det = hit->get_detector_id();
     if (!geom->isChamber(det) && !geom->isPropTube(det)) continue;
 
-    CalibParamXT::Set* xt = m_cal_xt->GetParam(det);
-    if (! det) {
-      cerr << "  WARNING:  Cannot find the in-time parameter for det=" << det << " in CalibDriftDist.\n";
+    int ele = hit->get_element_id();
+    TGraphErrors* gr_t2x;
+    TGraphErrors* gr_t2dx;
+    double center, width;
+    if (! m_cal_xt ->Find(det, gr_t2x, gr_t2dx) || 
+        ! m_cal_int->Find(det, ele, center, width) ) {
+      cerr << "  WARNING:  Cannot find the in-time parameter for det=" << det << " ele=" << ele << " in CalibDriftDist.\n";
+
       continue;
       //return Fun4AllReturnCodes::ABORTEVENT;
     }
@@ -90,6 +95,7 @@ int CalibDriftDist::process_event(PHCompositeNode* topNode)
 
     int ele = hit->get_element_id();
     hit->set_pos(geom->getMeasurement(det, ele));
+    std::cout<<"still still in CalibDriftDist process_event.  index = "<<hit->get_hit_id()<<", detID = "<<hit->get_detector_id()<<", elementID = "<<hit->get_element_id()<<", tdcTime = "<<hit->get_tdc_time()<<", driftDistance = "<<fabs(hit->get_drift_distance())<<", pos = "<<hit->get_pos()<<", is_in_time = "<<hit->is_in_time()<<", get_tdc_time = "<<hit->get_tdc_time()<<", t1 = "<<t1<<", t0 = "<<t0<<std::endl; //WPM
   }
   return Fun4AllReturnCodes::EVENT_OK;
 }

packages/reco/interface/FastTracklet.h

@@ -176,6 +176,7 @@ class Tracklet : public PHObject
     double Eval(const double* par);
     double Eval4(const double* par);
     double calcChisq();
+    double calcChisq_noDrift();
 
     //Add dummy hits
     void addDummyHits();
@@ -226,7 +227,14 @@ class Tracklet : public PHObject
 
     //List of signed hits
     std::list<SignedHit> hits;
-
+  SignedHit getHit(int _i){
+    std::list<SignedHit>::iterator it = hits.begin();
+    for(int i=0; i<_i; i++){
+      ++it;
+    }
+    return *it;
+  }
+
     //Corresponding prop. tube segments
     PropSegment seg_x;
     PropSegment seg_y;
@@ -238,6 +246,71 @@ class Tracklet : public PHObject
     double y0;
     double invP;
 
+    //This is an admittedly messy way of keeping track of various bits of needed information to describe the possible particle trajectories in a single-station tracklet
+    struct linedef {
+      double slopeX;
+      double slopeY;
+      double initialX;
+      double initialY;
+      double initialZ;
+      double slopeU;
+      double initialU;
+      double slopeV;
+      double initialV;
+
+      double wireHit1Pos;
+      double wireHit2Pos;
+      double wireHit1PosX;
+      double wireHit2PosX;
+      double wireHit1PosY;
+      double wireHit2PosY;
+      double wireHit1PosZ;
+      double wireHit2PosZ;
+      double wire1Slope;
+      double wire2Slope;
+      double wireIntercept1;
+      double wireIntercept2;
+
+      void print(){
+	std::cout<<"slopeX: "<<slopeX<<" slopeY: "<<slopeY<<" initialX: "<<initialX<<" initialY: "<<initialY<<" initialZ: "<<initialZ<<" slopeU: "<<slopeU<<" initialU: "<<initialU<<" slopeV: "<<slopeV<<" initialV: "<<initialV<<std::endl;
+      }
+
+    } ;
+
+  double st2X;
+  double st3X;
+  double st2Xsl;
+  double st3Xsl;
+  double st2U;
+  double st3U;
+  double st3V;
+  double st2V;
+  double st2Z;
+  double st3Z;
+  double st2Y;
+  double st3Y;
+  double st2Usl;
+  double st2Vsl;
+  double st3Usl;
+  double st3Vsl;
+  double st2UZ;
+  double st2VZ;
+
+  linedef acceptedXLine2;
+  linedef acceptedULine2;
+  linedef acceptedVLine2;
+  linedef acceptedXLine3;
+  linedef acceptedULine3;
+  linedef acceptedVLine3;
+
+    std::vector<linedef> possibleXLines;
+    std::vector<linedef> possibleULines;
+    std::vector<linedef> possibleVLines;
+
+    void getSlopesX(Hit hit1, Hit hit2);
+    void getSlopesU(Hit hit1, Hit hit2);
+    void getSlopesV(Hit hit1, Hit hit2);
+
     double err_tx;
     double err_ty;
     double err_x0;

packages/reco/interface/SRawEvent.cxx

@@ -206,7 +206,8 @@ std::list<Int_t> SRawEvent::getHitsIndexInDetector(Short_t detectorID, Double_t
 
     for(Int_t i = 0; i < fNHits[0]; i++)
     {
-        if(fAllHits[i].detectorID != detectorID) continue;
+      if(fAllHits[i].detectorID != detectorID) continue;
+      //std::cout<<"in getHitsIndexInDetector: index is "<<fAllHits[i].index<<", pos is "<<fAllHits[i].pos<<" and x_exp "<<x_exp<<"; detID is "<<fAllHits[i].detectorID<<std::endl; //WPM
         if(fabs(fAllHits[i].pos - x_exp) > win) continue;
 
         hit_list.push_back(i);
@@ -275,6 +276,8 @@ std::list<SRawEvent::hit_pair> SRawEvent::getPartialHitPairsInSuperDetector(Shor
         for(std::list<int>::iterator jter = _hitlist2.begin(); jter != _hitlist2.end(); ++jter)
         {
             index2++;
+	    //std::cout<<"in getPartialHitPairsInSuperDetector: index 1 is "<<fAllHits[*iter].index<<" and index 2 is "<<fAllHits[*jter].index<<"; detID 1 is "<<fAllHits[*iter].detectorID<<" and detID 2 is "<<fAllHits[*jter].detectorID<<std::endl; //WPM
+	    //std::cout<<"The difference between the two is "<<fabs(fAllHits[*iter].pos - fAllHits[*jter].pos)<<".  The allowed spacing is "<<spacing[detectorID]<<std::endl; //WPM
             if(fabs(fAllHits[*iter].pos - fAllHits[*jter].pos) > spacing[detectorID]) continue;
             _hitpairs.push_back(std::make_pair(*iter, *jter));
 
@@ -324,6 +327,8 @@ std::list<SRawEvent::hit_pair> SRawEvent::getPartialHitPairsInSuperDetector(Shor
         for(std::list<int>::iterator jter = _hitlist2.begin(); jter != _hitlist2.end(); ++jter)
         {
             index2++;
+	    //std::cout<<"in getPartialHitPairsInSuperDetector: index 1 is "<<fAllHits[*iter].index<<" and index 2 is "<<fAllHits[*jter].index<<"; detID 1 is "<<fAllHits[*iter].detectorID<<" and detID 2 is "<<fAllHits[*jter].detectorID<<std::endl; //WPM
+	    //std::cout<<"The difference between the two is "<<fabs(fAllHits[*iter].pos - fAllHits[*jter].pos)<<".  The allowed spacing is "<<spacing[detectorID]<<std::endl; //WPM
             if(fabs(fAllHits[*iter].pos - fAllHits[*jter].pos) > spacing[detectorID]) continue;
             _hitpairs.push_back(std::make_pair(*iter, *jter));