introducing come internal EDM

Tracker/TrackerRecAlgs/TrackerSeedFinder/src/TrackerSeedFinder.cxx

@@ -277,60 +277,50 @@ namespace Tracker
     ATH_MSG_VERBOSE( "  2nd station " << N_2_0+N_2_1+N_2_2 << " (" << vsp2.at(0).size() << ";" << vsp2.at(1).size() << ";" << vsp2.at(2).size() << ")");
     ATH_MSG_VERBOSE( "  3rd station " << N_3_0+N_3_1+N_3_2 << " (" << vsp3.at(0).size() << ";" << vsp3.at(1).size() << ";" << vsp3.at(2).size() << ")");
 
-    double a, b;
-    double xs=0, ys=0, zs=0, x2s=0, y2s=0, z2s=0, xys=0, xzs=0; 
-    int n=0, j=0;
-    vector<int> poss;
+    double xs=0, ys=0, zs=0, z2s=0, xzs=0, yzs=0; 
+    vector<seed> vseed;
 
     for (unsigned int i=0; i<vsp1[0].size(); i++) {
       for (unsigned int j=0; j<vsp1[1].size(); j++) {
     	for (unsigned int k=0; k<vsp1[2].size(); k++) {
 
-	  const Trk::SpacePoint* sp1 = vsp1[0].at(i);
-	  const Trk::SpacePoint* sp2 = vsp1[1].at(j);
-	  const Trk::SpacePoint* sp3 = vsp1[2].at(k);
+    	  ATH_MSG_VERBOSE( " station 1 / list of space points for triplets " 
+			   << vsp1[0].at(i)->clusterList().first->identify() << " " 
+			   << vsp1[1].at(j)->clusterList().first->identify() << " " 
+			   << vsp1[2].at(k)->clusterList().first->identify());
+
+	  vector<double> x; vector<double> y; vector<double> z;
+	  x.push_back(vsp1[0].at(i)->globalPosition().x()); x.push_back(vsp1[1].at(j)->globalPosition().x()); x.push_back(vsp1[2].at(k)->globalPosition().x());
+	  y.push_back(vsp1[0].at(i)->globalPosition().y()); y.push_back(vsp1[1].at(j)->globalPosition().y()); y.push_back(vsp1[2].at(k)->globalPosition().y());
+	  z.push_back(vsp1[0].at(i)->globalPosition().z()); z.push_back(vsp1[1].at(j)->globalPosition().z()); z.push_back(vsp1[2].at(k)->globalPosition().z());
+
+	  xs=x.at(0)+x.at(1)+x.at(2); ys=y.at(0)+y.at(1)+y.at(2); zs=z.at(0)+z.at(1)+z.at(2);
+	  z2s=pow(z.at(0),2)+pow(z.at(1),2)+pow(z.at(2),2);
+	  xzs=x.at(0)*z.at(0)+x.at(1)*z.at(1)+x.at(2)*z.at(2);
+	  yzs=y.at(0)*z.at(0)+y.at(1)*z.at(1)+y.at(2)*z.at(2);
+
+	  seed m_seed;
+	  m_seed.axz=(3*xzs-xs*zs)/(3*z2s-zs*zs);
+	  m_seed.bxz=(z2s*xs-xzs*zs)/(3*z2s-zs*zs);
+	  m_seed.ayz=(3*yzs-ys*zs)/(3*z2s-zs*zs);
+	  m_seed.byz=(z2s*ys-yzs*zs)/(3*z2s-zs*zs);
+	  m_seed.add_sp(vsp1[0].at(i)); m_seed.add_sp(vsp1[1].at(j)); m_seed.add_sp(vsp1[2].at(k));
+
+	  for(int pos=0; pos<3; pos++) {
+	    m_seed.chi2_xz+=pow(x.at(pos)-(m_seed.axz*z.at(pos)+m_seed.bxz),2);
+	    m_seed.chi2_yz+=pow(x.at(pos)-(m_seed.ayz*z.at(pos)+m_seed.byz),2);
+	  }
 
-	  double x1 = vsp1[0].at(i)->globalPosition().x();
-	  double x2 = vsp1[1].at(j)->globalPosition().x();
-	  double x3 = vsp1[2].at(k)->globalPosition().x();
-	  double y1 = vsp1[0].at(i)->globalPosition().y();
-	  double y2 = vsp1[1].at(j)->globalPosition().y();
-	  double y3 = vsp1[2].at(k)->globalPosition().y();
-	  double z1 = vsp1[0].at(i)->globalPosition().z();
-	  double z2 = vsp1[1].at(j)->globalPosition().z();
-	  double z3 = vsp1[2].at(k)->globalPosition().z();
+	  ATH_MSG_VERBOSE( " linear fit on plane xz / slope " << m_seed.axz << "; intercept " << m_seed.bxz);
+	  ATH_MSG_VERBOSE( " linear fit on plane yz / slope " << m_seed.ayz << "; intercept " << m_seed.byz);
+	  ATH_MSG_VERBOSE( " chi2 / on plane xz " << m_seed.chi2_xz << "; on plane yz " << m_seed.chi2_yz);
+
+	  vseed.push_back(m_seed);
 
-    	  ATH_MSG_VERBOSE( " station 1 / list of space points for triplets " 
-			   << sp1->clusterList().first->identify() << " " 
-			   << sp2->clusterList().first->identify() << " " 
-			   << sp3->clusterList().first->identify());
-
-	  xs=x1+x2+x3; ys=y1+y2+y3; zs=z1+z2+z3;
-	  x2s=pow(x1,2)+pow(x2,2)+pow(x3,2); y2s=pow(y1,2)+pow(y2,2)+pow(y3,2); z2s=pow(z1,2)+pow(z2,2)+pow(z3,2);
-	  xys=x1*y1+x2*y2+x3*y3; xzs=x1*z1+x2*z2+x3*z3;
-
-	  a=(n*xys-xs*ys)/(n*x2s-xs*xs);
-	  b=(x2s*ys-xs*xys)/(x2s*n-xs*xs);
-	  
-	  ATH_MSG_VERBOSE( " linear fit / slope " << a << "; intercept " << b);
     	}
       }
     }
 
-    vector<pair<int, const Trk::SpacePoint*>> mseed;
-    for(unsigned int i=0; i<poss.size(); i++) {
-      ATH_MSG_VERBOSE( "  with spacepoints " << poss.at(i));
-      mseed.push_back(msp.at(poss.at(i)));
-    }
-
-    double chi2=0;
-    for (vector<pair<int, const Trk::SpacePoint*>>::iterator it = mseed.begin(); it != mseed.end(); ++it) {
-      chi2+=fabs((*it).second->globalPosition().y() - a*(*it).second->globalPosition().x()+b);
-      ATH_MSG_VERBOSE( "  chi2 " << chi2);
-    }
-
-    mseeds.push_back(mseed);
-
     SG::WriteHandle<TrackerSeedCollection> seedContainer(m_trackerSeedContainerKey, ctx);
     ATH_CHECK(seedContainer.record( std::make_unique<TrackerSeedCollection>() ) );
     ATH_MSG_INFO("Created track seed container " << m_trackerSeedContainerKey.key());
@@ -338,13 +328,10 @@ namespace Tracker
     Tracker::TrackerSeed* trackerSeed = new Tracker::TrackerSeed();
     trackerSeed->set_id(TrackerSeed::TRIPLET_SP_FIRSTSTATION);
 
-    for (vector<vector<pair<int, const Trk::SpacePoint*>>>::iterator it = mseeds.begin(); it != mseeds.end(); ++it) {
-      vector<const Trk::SpacePoint*> mseed;
-      for (vector<pair<int, const Trk::SpacePoint*>>::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) {
-	mseed.push_back((*it2).second);
-      }
+    for (vector<seed>::iterator it = vseed.begin(); it != vseed.end(); ++it) {
+      vector<const Trk::SpacePoint*> mseed = (*it).vsp;
       trackerSeed->add(mseed);
-      seedContainer->push_back(trackerSeed);
+      seedContainer->push_back(trackerSeed);      
     }
 
     return StatusCode::SUCCESS;

Tracker/TrackerRecAlgs/TrackerSeedFinder/src/TrackerSeedFinder.h

@@ -61,7 +61,6 @@ namespace Tracker {
     
   private:
 
-    StatusCode printAllPair(vector<pair<int, vector<string> > >&, int, string) const;
     StatusCode makeTriplets(vector< vector<const Trk::SpacePoint*> >&, unsigned int, string) const;
 
     TrackerSeedFinder() = delete;
@@ -106,6 +105,16 @@ namespace Tracker {
     ServiceHandle<ITHistSvc>  m_thistSvc;
     
     mutable int N_1_0=0, N_1_1=0, N_1_2=0, N_2_0=0, N_2_1=0, N_2_2=0, N_3_0=0, N_3_1=0, N_3_2=0;
+
+    struct seed {
+      vector<const Trk::SpacePoint*> vsp;
+      double axz, bxz, ayz, byz;
+      double chi2_xz, chi2_yz;
+
+      void add_sp(const Trk::SpacePoint* sp) {
+	vsp.push_back(sp);
+      }
+    };
     
   };