THcDriftChamber.cxx

				  break;
			  default:
				  wireDistance = 0.0;
			  }

			  Double_t timeCorrection = wireDistance/fWireVelocity;

			  if(fFixPropagationCorrection==0) { // ENGINE behavior
				  Double_t time=hit->GetTime();
				  hit->SetTime(time - timeCorrection);
		          hit->ConvertTimeToDist();
			  }
			  else{
				  Double_t time=hit->GetTime();
				  Double_t dist=hit->GetDist();

				  hit->SetTime(time - timeCorrection);
				  //double usingOldTime = time-time_corr;
				  //hit->SetTime(time- time_corr);

				  hit->ConvertTimeToDist();
				  sp->SetHitDist(ihit, hit->GetDist());

				  hit->SetTime(time);	// Restore time
				  hit->SetDist(dist);	// Restore distance
			  }

		  }
		  /////////////////////////////////////////////////////////////

		  //     if (fhdebugflagpr) cout << "Correcting hit " << hit << " " << plane->GetPlaneNum() << " " << isp << "/" << ihit << "  " << x << "," << y << endl;
		  // Fortran ENGINE does not do this check, so hits can get "corrected"
		  // multiple times if they belong to multiple space points.
		  // To reproduce the precise ENGINE behavior, remove this if condition.
		  else{
			  if(fFixPropagationCorrection==0) { // ENGINE behavior
				  hit->SetTime(hit->GetTime() - plane->GetCentralTime()
					  + plane->GetDriftTimeSign()*time_corr);
				  hit->ConvertTimeToDist();
			  //      hit->SetCorrectedStatus(1);
			  } else {
			  // New behavior: Save corrected distance with the hit in the space point
			  // so that the same hit can have a different correction depending on
			  // which space point it is in.
			  //
			  // This is a hack now because the converttimetodist method is connected to the hit
			  // so I compute the corrected time and distance, and then restore the original
			  // time and distance.  Can probably add a method to hit that does a conversion on a time
			  // but does not modify the hit data.
			  Double_t time=hit->GetTime();
			  Double_t dist=hit->GetDist();
			  hit->SetTime(time - plane->GetCentralTime()
					  + plane->GetDriftTimeSign()*time_corr);
			  hit->ConvertTimeToDist();
			  sp->SetHitDist(ihit, hit->GetDist());
			  hit->SetTime(time);	// Restore time
			  hit->SetDist(dist);	// Restore distance
			  }
		  }
	  }
  }
}
UInt_t THcDriftChamber::Count1Bits(UInt_t x)
// From http://graphics.stanford.edu/~seander/bithacks.html
{
  x = x - ((x >> 1) & 0x55555555);
  x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
  return (((x + (x >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
}

//_____________________________________________________________________________
void THcDriftChamber::LeftRight()
{
  // For each space point,
  // Fit stubs to all possible left-right combinations of drift distances
  // and choose the set with the minimum chi**2.

  for(Int_t isp=0; isp<fNSpacePoints; isp++) {
    // Build a bit pattern of which planes are hit
    THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
    Int_t nhits = sp->GetNHits();
    UInt_t bitpat  = 0;		// Bit pattern of which planes are hit
    Double_t maxchi2= 1.0e10;
    Double_t minchi2 = maxchi2;
    Double_t tmp_minchi2=maxchi2;
    Double_t minxp = 0.25;
    Int_t hasy1 = -1;
    Int_t hasy2 = -1;
    Int_t plusminusknown[nhits];
    Int_t plusminusbest[nhits];
    Int_t plusminus[nhits];	// ENGINE makes this array float.  Why?
    Int_t tmp_plusminus[nhits];
    Int_t plane_list[nhits];
    Double_t stub[4];
    Double_t tmp_stub[4];
    Int_t nplusminus;

    if(nhits < 0) {
      if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() nhits < 0" << endl;
    } else if (nhits==0) {
      if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() nhits = 0" << endl;
    }
    for(Int_t ihit=0;ihit < nhits;ihit++) {
      THcDCHit* hit = sp->GetHit(ihit);
      Int_t pindex = hit->GetPlaneIndex();
      plane_list[ihit] = pindex;

      bitpat |= 1<<pindex;

      plusminusknown[ihit] = 0;

      if(pindex == YPlaneInd) hasy1 = ihit;
      if(pindex == YPlanePInd) hasy2 = ihit;
    }
    nplusminus = 1<<nhits;
    if(fHMSStyleChambers) {
      Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
      if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes

	if(sp->GetHit(hasy2)->GetPos() <=
	   sp->GetHit(hasy1)->GetPos()) {
	  plusminusknown[hasy1] = -1;
	  plusminusknown[hasy2] = 1;
	} else {
	  plusminusknown[hasy1] = 1;
	  plusminusknown[hasy2] = -1;
	}
	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 && TMath::Abs(hit2->GetPos()-hit1->GetPos())<0.51) { // 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);
      }//end if fSmallAngleApprox!=0
    }//end else SOS style
    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;
	 chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub);
	if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl;
	if(chi2 < minchi2) {
	  if (fStubMaxXPDiff<100. ) {
	    // 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]]);
	    Double_t xp_expect = sp->GetX()*fRatio_xpfp_to_xfp;
	    if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
	      minchi2 = chi2;
	      for(Int_t ihit=0;ihit<nhits;ihit++) {
		plusminusbest[ihit] = plusminus[ihit];
	      }
              sp->SetStub(stub);
	    } else {		// Record best stub failing angle cut
              if (chi2 < tmp_minchi2) {
	      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];
	    }
            sp->SetStub(stub);
	  }
	}
	///////////////
      } 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];
	  }
          sp->SetStub(stub);
	}
      }//////////
      else {
	if (fhdebugflagpr) cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
      }
    } // End loop of pm combinations

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

    // 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 spaleftce 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];
    sp->SetStub(stub);
    //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->GetHitDist(ihit)
      - 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.
  delete fCosBeta; fCosBeta = NULL;
  delete fSinBeta; fSinBeta = NULL;
  delete fTanBeta; fTanBeta = NULL;
  delete fSigma; fSigma = NULL;
  delete fPsi0; fPsi0 = NULL;
  delete fStubCoefs; fStubCoefs = NULL;

  // Need to delete each element of the fAA3Inv map

}

//_____________________________________________________________________________
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)
////////////////////////////////////////////////////////////////////////////////