THcDC.cxx

    			  for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
    				  AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
    						  fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
							  pow(fSigma[planes[ihit]],2);
    			  }
    		  }
    	  }
      }

      // 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++) {
    	  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/fSigma[planes[ihit]],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++) {
          if (ihit != ipl_hit) {
	  TT[irayp] += (coords[ihit]*
			fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])
	    /pow(fSigma[planes[ihit]],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++) {
              if (ihit != ipl_hit) {
	      AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
		fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
		pow(fSigma[planes[ihit]],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)
////////////////////////////////////////////////////////////////////////////////