THcDC.cxx

        if (fdebuglinkstubs)
          cout << "EPIC FAIL 3:  Too many tracks found in THcDC::LinkStubs" << endl;
        fNDCTracks = 0;
        // Do something here to fail this event
        return; // Max # of allowed tracks
      }
    }
  }
  ///
  if (fdebuglinkstubs) {
    cout << " Number of tracks from link stubs = " << fNDCTracks << endl;
    printf("%s %s \n", "Track", "Plane Wire ");
    for (UInt_t itrack = 0; itrack < fNDCTracks; itrack++) {
      THcDCTrack* tempTrack = (THcDCTrack*)(fDCTracks->At(itrack));
      printf("%-5d  ", itrack + 1);
      for (Int_t ihit = 0; ihit < tempTrack->GetNHits(); ihit++) {
        THcDCHit* hit = (THcDCHit*)(tempTrack->GetHit(ihit));
        printf("%3d %3d", hit->GetPlaneNum(), hit->GetWireNum());
      }
      printf("\n");
    }
  }
}

//_____________________________________________________________________________
void THcDC::TrackFit() {
  /**
     Primary track fitting routine
  */

  // Number of ray parameters in focal plane.
  const Int_t raycoeffmap[] = {4, 5, 2, 3};

  Double_t dummychi2 = 1.0E4;
  for (UInt_t itrack = 0; itrack < fNDCTracks; itrack++) {
    //    Double_t chi2 = dummychi2;
    //    Int_t htrack_fit_num = itrack;
    THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(itrack));

    Double_t coords[theDCTrack->GetNHits()];
    Int_t    planes[theDCTrack->GetNHits()];
    for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {
      THcDCHit* hit = theDCTrack->GetHit(ihit);
      planes[ihit]  = hit->GetPlaneNum() - 1;
      if (fFixLR) {
        if (fFixPropagationCorrection) {
          coords[ihit] = hit->GetPos() + theDCTrack->GetHitLR(ihit) * theDCTrack->GetHitDist(ihit);
        } else {
          coords[ihit] = hit->GetPos() + theDCTrack->GetHitLR(ihit) * hit->GetDist();
        }

      } else {
        if (fFixPropagationCorrection) {
          coords[ihit] = hit->GetPos() + hit->GetLR() * theDCTrack->GetHitDist(ihit);
        } else {
          coords[ihit] = hit->GetCoord();
        }
      }

    } // end loop over hits

    theDCTrack->SetNFree(theDCTrack->GetNHits() - NUM_FPRAY);
    Double_t chi2 = dummychi2;
    if (theDCTrack->GetNFree() > 0) {
      TVectorD TT(NUM_FPRAY);
      TMatrixD AA(NUM_FPRAY, NUM_FPRAY);
      for (Int_t irayp = 0; irayp < NUM_FPRAY; irayp++) {
        TT[irayp] = 0.0;
        for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {

          THcDCHit* hit = theDCTrack->GetHit(ihit);

          TT[irayp] += (coords[ihit] * fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]) /
                       pow(hit->GetWireSigma(), 2);
          //	  if (hit->GetPlaneNum()==5)
          //	    {
          //	      //	cout << "Plane: " << hit->GetPlaneNum() << endl;
          //	      //cout << "Wire: " <<hit->GetWireNum() << endl;
          //	      //cout << "Sigma: " << hit->GetWireSigma() << endl;
          //	    }

        } // end hit loop
      }
      for (Int_t irayp = 0; irayp < NUM_FPRAY; irayp++) {
        for (Int_t jrayp = 0; jrayp < NUM_FPRAY; jrayp++) {
          AA[irayp][jrayp] = 0.0;
          if (jrayp < irayp) { // Symmetric
            AA[irayp][jrayp] = AA[jrayp][irayp];
          } else {
            for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {

              THcDCHit* hit = theDCTrack->GetHit(ihit);

              AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]] *
                                  fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]] /
                                  pow(hit->GetWireSigma(), 2);

            } // end ihit loop
          }
        }
      }

      // Solve 4x4 equations
      TVectorD dray(NUM_FPRAY);
      // Should check that it is invertable
      AA.Invert();
      dray = AA * TT;
      //      cout << "DRAY: " << dray[0] << " "<< dray[1] << " "<< dray[2] << " "<< dray[3] << " "
      //      << endl; if(bad_determinant) {
      //	dray[0] = dray[1] = 10000.; dray[2] = dray[3] = 2.0;
      //      }
      // Calculate hit coordinate for each plane for chi2 and efficiency
      // calculations

      // Make sure fCoords, fResiduals, and fDoubleResiduals are clear
      for (Int_t iplane = 0; iplane < fNPlanes; iplane++) {
        Double_t coord = 0.0;
        for (Int_t ir = 0; ir < NUM_FPRAY; ir++) {
          coord += fPlaneCoeffs[iplane][raycoeffmap[ir]] * dray[ir];
          // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
        }
        theDCTrack->SetCoord(iplane, coord);
      }
      // Compute Chi2 and residuals
      chi2 = 0.0;
      for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {

        THcDCHit* hit = theDCTrack->GetHit(ihit);

        Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
        theDCTrack->SetResidual(planes[ihit], residual);

        //  	  double track_coord = theDCTrack->GetCoord(planes[ihit]);
        // cout<<planes[ihit]<<"\t"<<track_coord<<"\t"<<coords[ihit]<<"\t"<<residual<<endl;
        chi2 += pow(residual / hit->GetWireSigma(), 2);
      }

      theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
    }
    theDCTrack->SetChisq(chi2);

    // calculate ray without a plane in track
    for (Int_t ipl_hit = 0; ipl_hit < theDCTrack->GetNHits(); ipl_hit++) {

      if (theDCTrack->GetNFree() > 0) {
        TVectorD TT(NUM_FPRAY);
        TMatrixD AA(NUM_FPRAY, NUM_FPRAY);
        for (Int_t irayp = 0; irayp < NUM_FPRAY; irayp++) {
          TT[irayp] = 0.0;
          for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {

            THcDCHit* hit = theDCTrack->GetHit(ihit);

            if (ihit != ipl_hit) {
              TT[irayp] += (coords[ihit] * fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]) /
                           pow(hit->GetWireSigma(), 2);
            }
          }
        }
        for (Int_t irayp = 0; irayp < NUM_FPRAY; irayp++) {
          for (Int_t jrayp = 0; jrayp < NUM_FPRAY; jrayp++) {
            AA[irayp][jrayp] = 0.0;
            if (jrayp < irayp) { // Symmetric
              AA[irayp][jrayp] = AA[jrayp][irayp];
            } else {

              for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {

                THcDCHit* hit = theDCTrack->GetHit(ihit);

                if (ihit != ipl_hit) {
                  AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]] *
                                      fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]] /
                                      pow(hit->GetWireSigma(), 2);
                }
              }
            }
          }
        }
        //
        // Solve 4x4 equations
        // Should check that it is invertable
        TVectorD dray(NUM_FPRAY);
        AA.Invert();
        dray           = AA * TT;
        Double_t coord = 0.0;
        for (Int_t ir = 0; ir < NUM_FPRAY; ir++) {
          coord += fPlaneCoeffs[planes[ipl_hit]][raycoeffmap[ir]] * dray[ir];
          // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
        }
        Double_t residual = coords[ipl_hit] - coord;
        theDCTrack->SetResidualExclPlane(planes[ipl_hit], residual);
      }
    }
  }
  // Calculate residual without plane

  // Calculate residuals for each chamber if in single stub mode
  // and there was a track found in each chamber
  // Specific for two chambers.  Can/should it be generalized?

  if (fSingleStub != 0) {
    if (fNDCTracks == 2) {
      THcDCTrack* theDCTrack1 = static_cast<THcDCTrack*>(fDCTracks->At(0));
      THcDCTrack* theDCTrack2 = static_cast<THcDCTrack*>(fDCTracks->At(1));
      //      Int_t itrack=0;
      Int_t     ihit    = 0;
      THcDCHit* hit     = theDCTrack1->GetHit(ihit);
      Int_t     plane   = hit->GetPlaneNum() - 1;
      Int_t     chamber = fNChamber[plane];
      if (chamber == 1) {
        //	itrack=1;
        hit     = theDCTrack2->GetHit(ihit);
        plane   = hit->GetPlaneNum() - 1;
        chamber = fNChamber[plane];
        if (chamber == 2) {
          Double_t ray1[4];
          Double_t ray2[4];
          theDCTrack1->GetRay(ray1);
          theDCTrack2->GetRay(ray2);
          //	  itrack = 1;
          // Loop over hits in second chamber
          for (Int_t ihit = 0; ihit < theDCTrack2->GetNHits(); ihit++) {
            // Calculate residual in second chamber from first chamber track
            THcDCHit* hit   = theDCTrack2->GetHit(ihit);
            Int_t     plane = hit->GetPlaneNum() - 1;
            Double_t  pos   = DpsiFun(ray1, plane);
            Double_t  coord;
            if (fFixLR) {
              if (fFixPropagationCorrection) {
                coord = hit->GetPos() + theDCTrack2->GetHitLR(ihit) * theDCTrack2->GetHitDist(ihit);
              } else {
                coord = hit->GetPos() + theDCTrack2->GetHitLR(ihit) * hit->GetDist();
              }
            } else {
              if (fFixPropagationCorrection) {
                coord = hit->GetPos() + hit->GetLR() * theDCTrack2->GetHitDist(ihit);
              } else {
                coord = hit->GetCoord();
              }
            }
            theDCTrack1->SetDoubleResidual(plane, coord - pos);
            //  hdc_dbl_res(pln) = hdc_double_residual(1,pln)  for hists
          }
          //	  itrack=0;
          // Loop over hits in first chamber
          for (Int_t ihit = 0; ihit < theDCTrack1->GetNHits(); ihit++) {
            // Calculate residual in first chamber from second chamber track
            THcDCHit* hit   = theDCTrack1->GetHit(ihit);
            Int_t     plane = hit->GetPlaneNum() - 1;
            Double_t  pos   = DpsiFun(ray1, plane);
            Double_t  coord;
            if (fFixLR) {
              if (fFixPropagationCorrection) {
                coord = hit->GetPos() + theDCTrack1->GetHitLR(ihit) * theDCTrack1->GetHitDist(ihit);
              } else {
                coord = hit->GetPos() + theDCTrack1->GetHitLR(ihit) * hit->GetDist();
              }
            } else {
              if (fFixPropagationCorrection) {
                coord = hit->GetPos() + hit->GetLR() * theDCTrack1->GetHitDist(ihit);
              } else {
                coord = hit->GetCoord();
              }
            }
            theDCTrack2->SetDoubleResidual(plane, coord - pos);
            //  hdc_dbl_res(pln) = hdc_double_residual(1,pln)  for hists
          }
        }
      }
    }
  }
  //
  if (fdebugtrackprint) {
    printf("%5s %-14s %-14s %-14s %-14s  %-10s %-10s \n", "Track", "x_t", "y_t", "xp_t", "yp_t",
           "chi2", "DOF");
    printf("%5s %-14s %-14s %-14s %-14s  %-10s %-10s \n", "     ", "[cm]", "[cm]", "[rad]", "[rad]",
           " ", " ");
    for (UInt_t itr = 0; itr < fNDCTracks; itr++) {
      THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(itr));
      printf("%-5d %14.6e %14.6e %14.6e %14.6e %10.3e %3d \n", itr + 1, theDCTrack->GetX(),
             theDCTrack->GetY(), theDCTrack->GetXP(), theDCTrack->GetYP(), theDCTrack->GetChisq(),
             theDCTrack->GetNFree());
    }
    for (UInt_t itr = 0; itr < fNDCTracks; itr++) {
      printf("%s %5d \n", "Hit info for track number = ", itr + 1);
      printf("%5s %-15s %-15s %-15s \n", "Plane", "WIRE_COORD", "Fit postiion", "Residual");
      THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(itr));
      for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {
        THcDCHit* hit   = theDCTrack->GetHit(ihit);
        Int_t     plane = hit->GetPlaneNum() - 1;
        Double_t  coords_temp;
        if (fFixLR) {
          if (fFixPropagationCorrection) {
            coords_temp = hit->GetPos() + theDCTrack->GetHitLR(ihit) * theDCTrack->GetHitDist(ihit);
          } else {
            coords_temp = hit->GetPos() + theDCTrack->GetHitLR(ihit) * hit->GetDist();
          }
        } else {
          if (fFixPropagationCorrection) {
            coords_temp = hit->GetPos() + hit->GetLR() * theDCTrack->GetHitDist(ihit);
          } else {
            coords_temp = hit->GetCoord();
          }
        }
        printf("%-5d %15.7e %15.7e %15.7e \n", plane + 1, coords_temp, theDCTrack->GetCoord(plane),
               theDCTrack->GetResidual(plane));
      }
    }
  }

  //
}
//
//
Double_t THcDC::DpsiFun(Double_t ray[4], Int_t plane) {
  /**
    this function calculates the psi coordinate of the intersection
    of a ray (defined by ray) with a hms wire chamber plane. the geometry
    of the plane is contained in the coeff array calculated in the
    array hplane_coeff
    Note it is call by MINUIT via H_FCNCHISQ and so uses double precision
    variables

    the ray is defined by
    x = (z-zt)*tan(xp) + xt
    y = (z-zt)*tan(yp) + yt
     at some fixed value of zt*
    ray(1) = xt
    ray(2) = yt
    ray(3) = tan(xp)
    ray(4) = tan(yp)
  */

  Double_t infinity  = 1.0E+20;
  Double_t cinfinity = 1 / infinity;
  Double_t DpsiFun   = ray[2] * ray[1] * fPlaneCoeffs[plane][0] +
                     ray[3] * ray[0] * fPlaneCoeffs[plane][1] + ray[2] * fPlaneCoeffs[plane][2] +
                     ray[3] * fPlaneCoeffs[plane][3] + ray[0] * fPlaneCoeffs[plane][4] +
                     ray[1] * fPlaneCoeffs[plane][5];
  Double_t denom =
      ray[2] * fPlaneCoeffs[plane][6] + ray[3] * fPlaneCoeffs[plane][7] + fPlaneCoeffs[plane][8];
  if (TMath::Abs(denom) < cinfinity) {
    DpsiFun = infinity;
  } else {
    DpsiFun = DpsiFun / denom;
  }
  return (DpsiFun);
}

//_____________________________________________________________________________
Int_t THcDC::End(THaRunBase* run) {
  //  EffCalc();
  MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
  return 0;
}

//_____________________________________________________________________________
void THcDC::EffInit() {
  /**
     Create, and initialize counters used to calculate
     efficiencies.  Register the counters in gHcParms so that the
     variables can be used in end of run reports.
  */

  delete[] fNChamHits;
  fNChamHits = new Int_t[fNChambers];
  delete[] fPlaneEvents;
  fPlaneEvents = new Int_t[fNPlanes];

  fTotEvents = 0;
  for (UInt_t i = 0; i < fNChambers; i++) {
    fNChamHits[i] = 0;
  }
  for (Int_t i = 0; i < fNPlanes; i++) {
    fPlaneEvents[i] = 0;
  }
  gHcParms->Define(Form("%sdc_tot_events", fPrefix), "Total DC Events", fTotEvents);
  gHcParms->Define(Form("%sdc_cham_hits[%d]", fPrefix, fNChambers),
                   "N events with hits per chamber", *fNChamHits);
  gHcParms->Define(Form("%sdc_events[%d]", fPrefix, fNPlanes), "N events with hits per plane",
                   *fPlaneEvents);
}

//_____________________________________________________________________________
void THcDC::Eff() {
  /**
     Accumulate statistics for efficiency calculations
  */

  fTotEvents++;
  for (UInt_t i = 0; i < fNChambers; i++) {
    if (fChambers[i]->GetNHits() > 0)
      fNChamHits[i]++;
  }
  for (Int_t i = 0; i < fNPlanes; i++) {
    if (fPlanes[i]->GetNHits() > 0)
      fPlaneEvents[i]++;
  }
  return;
}

ClassImp(THcDC)
    ////////////////////////////////////////////////////////////////////////////////