THcDriftChamber.cxx

	nplusminus = 1<<(nhits-2);
	//	if (fhdebugflagpr) cout << " Small angle approx = " << smallAngOK << " " << plusminusknown[hasy1] << endl;
	//if (fhdebugflagpr) cout << "pm =  " << plusminusknown[hasy2] << " " << hasy1 << " " << hasy2 << endl;
	//if (fhdebugflagpr) cout << " Plane index " << YPlaneInd << " " << YPlanePInd << endl;
      }
    } else {			// SOS Style
      if(fSmallAngleApprox !=0) {
	// Brookhaven style chamber L/R code
	Int_t npaired=0;
	for(Int_t ihit1=0;ihit1 < nhits;ihit1++) {
	  THcDCHit* hit1 = sp->GetHit(ihit1);
	  Int_t pindex1=hit1->GetPlaneIndex();
	  if((pindex1%2)==0) { // Odd plane (or even index)
	    for(Int_t ihit2=0;ihit2<nhits;ihit2++) {
	      THcDCHit* hit2 = sp->GetHit(ihit2);
	      if(hit2->GetPlaneIndex()-pindex1 == 1) { // Adjacent plane
		if(hit2->GetPos() <= hit1->GetPos()) {
		  plusminusknown[ihit1] = -1;
		  plusminusknown[ihit2] = 1;
		} else {
		  plusminusknown[ihit1] = 1;
		  plusminusknown[ihit2] = -1;
		}
		npaired+=2;
	      }
	    }
	  }
	}
	nplusminus = 1 << (nhits-npaired);
      }
    }
    if(nhits < 2) {
      if (fdebugstubchisq) cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
    } else if (nhits == 2) {
      if (fdebugstubchisq) cout << "THcDriftChamber::LeftRight: numhits-2 = 0" << endl;
    }
    Int_t nplaneshit = Count1Bits(bitpat);
    //if (fhdebugflagpr) cout << " num of pm = " << nplusminus << " num of hits =" << nhits << endl;
    // Use bit value of integer word to set + or -
    // Loop over all combinations of left right.
    for(Int_t pmloop=0;pmloop<nplusminus;pmloop++) {
      Int_t iswhit = 1;
      for(Int_t ihit=0;ihit<nhits;ihit++) {
	if(plusminusknown[ihit]!=0) {
	  plusminus[ihit] = plusminusknown[ihit];
	} else {
	  // Max hits per point has to be less than 32.  
	  if(pmloop & iswhit) {
	    plusminus[ihit] = 1;
	  } else {
	    plusminus[ihit] = -1;
	  }
	  iswhit <<= 1;
	}
      }
      if (nplaneshit >= fNPlanes-1) {
	Double_t chi2 = FindStub(nhits, sp,
				     plane_list, bitpat, plusminus, stub);
	if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl;
	if(chi2 < minchi2) {
	  if(fHMSStyleChambers) { // Perhaps a different flag here
	    // Take best chi2 IF x' of the stub agrees with x' as expected from x.
	    // Sometimes an incorrect x' gives a good chi2 for the stub, even though it is
	    // not the correct left/right combination for the real track.
	    // Rotate x'(=stub(3)) to hut coordinates and compare to x' expected from x.
	    // THIS ASSUMES STANDARD HMS TUNE!!!!, for which x' is approx. x/875.
	    if(stub[2]*fTanBeta[plane_list[0]]==-1.0) {
	      if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() Error 3" << endl;
	    }
	    Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
	      /(1+stub[2]*fTanBeta[plane_list[0]]);
	    // Tune depdendent.  Definitely HMS specific
	    Double_t xp_expect = sp->GetX()/875.0;
	    if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
	      minchi2 = chi2;
	      for(Int_t ihit=0;ihit<nhits;ihit++) {
		plusminusbest[ihit] = plusminus[ihit];
	      }
	      Double_t *spstub = sp->GetStubP();
	      for(Int_t i=0;i<4;i++) {
		spstub[i] = stub[i];
	      }
	    } else {		// Record best stub failing angle cut
	      tmp_minchi2 = chi2;
	      for(Int_t ihit=0;ihit<nhits;ihit++) {
		tmp_plusminus[ihit] = plusminus[ihit];
	      }
	      for(Int_t i=0;i<4;i++) {
		tmp_stub[i] = stub[i];
	      }
	    }
	  } else { // Not HMS specific
	    minchi2 = chi2;
	    for(Int_t ihit=0;ihit<nhits;ihit++) {
	      plusminusbest[ihit] = plusminus[ihit];
	    }
	    Double_t *spstub = sp->GetStubP();
	    for(Int_t i=0;i<4;i++) {
	      spstub[i] = stub[i];
	    }
	  }
	}
      } else if (nplaneshit >= fNPlanes-2 && fHMSStyleChambers) { // Two planes missing
	Double_t chi2 = FindStub(nhits, sp,
				     plane_list, bitpat, plusminus, stub); 
	//if(debugging)
	//if (fhdebugflagpr) cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
	// Isn't this a bad idea, doing == with reals
	if(stub[2]*fTanBeta[plane_list[0]] == -1.0) {
	  if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() Error 3" << endl;
	}
	Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
	  /(1+stub[2]*fTanBeta[plane_list[0]]);
	if(TMath::Abs(xp_fit) <= minxp) {
	  minxp = TMath::Abs(xp_fit);
	  minchi2 = chi2;
	  for(Int_t ihit=0;ihit<nhits;ihit++) {
	    plusminusbest[ihit] = plusminus[ihit];
	  }
	  Double_t *spstub = sp->GetStubP();
	  for(Int_t i=0;i<4;i++) {
	    spstub[i] = stub[i];
	  }
	}
      } else {
	if (fhdebugflagpr) cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
      }
    } // End loop of pm combinations

    Double_t *spstub = sp->GetStubP();
    if(minchi2 > 9.9e9) {	// No track passed angle cut
      minchi2 = tmp_minchi2;
      for(Int_t ihit=0;ihit<nhits;ihit++) {
	plusminusbest[ihit] = tmp_plusminus[ihit];
      }
      for(Int_t i=0;i<4;i++) {
	spstub[i] = tmp_stub[i];
      }
      
    }

    // Calculate final coordinate based on plusminusbest
    // Update the hit positions in the space points
    for(Int_t ihit=0; ihit<nhits; ihit++) {
      // Save left/right status with the hit and in the space point
      // In THcDC will decide which to used based on fix_lr flag
      sp->GetHit(ihit)->SetLeftRight(plusminusbest[ihit]);
      sp->SetHitLR(ihit, plusminusbest[ihit]);
    }

    // Stubs are calculated in rotated coordinate system
    // (I think this rotates in case chambers not perpendicular to central ray)
    Int_t pindex=plane_list[0];
    if(spstub[2] - fTanBeta[pindex] == -1.0) {
      if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight(): stub3 error" << endl;
    }
    stub[2] = (spstub[2] - fTanBeta[pindex])
      / (1.0 + spstub[2]*fTanBeta[pindex]);
    if(spstub[2]*fSinBeta[pindex] ==  -fCosBeta[pindex]) {
      if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight(): stub4 error" << endl;
    }
    stub[3] = spstub[3]
      / (spstub[2]*fSinBeta[pindex]+fCosBeta[pindex]);
    stub[0] = spstub[0]*fCosBeta[pindex]
      - spstub[0]*stub[2]*fSinBeta[pindex];
    stub[1] = spstub[1]
      - spstub[1]*stub[3]*fSinBeta[pindex];
    for(Int_t i=0;i<4;i++) {
      spstub[i] = stub[i];
    }
    //if (fhdebugflagpr) cout << " Left/Right space pt " << isp+1 << " " << stub[0]<< " " << stub[1] << " " << stub[2]<< " " << stub[3] << endl;
      }
  // Option to print stubs
}
//_____________________________________________________________________________
Double_t THcDriftChamber::FindStub(Int_t nhits, THcSpacePoint *sp,
				       Int_t* plane_list, UInt_t bitpat,
				       Int_t* plusminus, Double_t* stub)
{
  // For a given combination of L/R, fit a stub to the space point
  // This method does a linear least squares fit of a line to the
  // hits in an individual chamber.  It assumes that the y slope is 0 
  // The wire coordinate is calculated by
  //          wire center + plusminus*(drift distance).
  // Method is called in a loop over all combinations of plusminus
  Double_t zeros[] = {0.0,0.0,0.0};
  TVectorD TT; TT.Use(3, zeros); // X, X', Y
  Double_t dpos[nhits];

  for(Int_t ihit=0;ihit<nhits; ihit++) {
    dpos[ihit] = sp->GetHit(ihit)->GetPos()
      + plusminus[ihit]*sp->GetHit(ihit)->GetDist()
      - fPsi0[plane_list[ihit]];
    for(Int_t index=0;index<3;index++) {
      TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
	/fSigma[plane_list[ihit]];
    }
  }
  //  fAA3Inv[bitpat]->Print();
  //  if (fhdebugflagpr) cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
  //  TT->Print();

  TT *= *fAA3Inv[bitpat];
  // if (fhdebugflagpr) cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;

  // Calculate Chi2.  Remember one power of sigma is in fStubCoefs
  stub[0] = TT[0];
  stub[1] = TT[1];
  stub[2] = TT[2];
  stub[3] = 0.0;
  Double_t chi2=0.0;
  for(Int_t ihit=0;ihit<nhits; ihit++) {
    chi2 += pow( dpos[ihit]/fSigma[plane_list[ihit]]
		 - fStubCoefs[plane_list[ihit]][0]*stub[0]
		 - fStubCoefs[plane_list[ihit]][1]*stub[1]
		 - fStubCoefs[plane_list[ihit]][2]*stub[2]
		 , 2);
  }
  return(chi2);
}

//_____________________________________________________________________________
THcDriftChamber::~THcDriftChamber()
{
  // Destructor. Remove variables from global list.

  if (fhdebugflagpr) cout << fChamberNum << ": delete fSpacePoints: " << hex << fSpacePoints << endl;
  delete fSpacePoints;
  // Should delete the matrices

  if( fIsSetup )
    RemoveVariables();
  if( fIsInit )
    DeleteArrays();
  if (fTrackProj) {
    fTrackProj->Clear();
    delete fTrackProj; fTrackProj = 0;
  }
}

//_____________________________________________________________________________
void THcDriftChamber::DeleteArrays()
{
  // Delete member arrays. Used by destructor.

}

//_____________________________________________________________________________
inline 
void THcDriftChamber::Clear( const Option_t* )
{
  // Reset per-event data.
  //  fNhits = 0;

  //  fTrackProj->Clear();
  fNhits = 0;

}

//_____________________________________________________________________________
Int_t THcDriftChamber::ApplyCorrections( void )
{
  return(0);
}

ClassImp(THcDriftChamber)
////////////////////////////////////////////////////////////////////////////////