From 4a7af64cd3ab3a7ab0e56630cd404c23b1ae02e4 Mon Sep 17 00:00:00 2001
From: Carlos Yero <cyero002@fiu.edu>
Date: Thu, 23 Aug 2018 14:09:12 -0400
Subject: [PATCH] Updated DC Classes to have DC Sigma per wire
1) THcDCWire.h
a) added sigma to THcDCWire constructor which fills fSigmaWire
b) added GetSimga method
c) added SetSigma method
2) THcDCHit.h
a) Added GetWireSigma method
3) THcDriftChamberPlane.h
a) Added flag fUsingSigmaPerWire
b) Added array fSigmaWire
4) THcDriftChamberPlane.cxx
a) In ReadDatabase, set fUsingSigmaPerWire=0 by default
with reading of parameter h_using_sigma_per_wire optional
for setting fUsingSigmaPerWire.
b) If fUsingSigmaPerWire=0 then sets all wires in array fSigmaWire to
fSigma ( which is the simga per plane).
c) If fUsingSigmaPerWire=1 then reads in the sigma per wire from
a parameter file.
d) Used new THcDCWire constructor to set the sigma for wire to fSigmaWire[nwire]
e) Also change code so that the tzero offset is set in the THcDCWire constructor
f) In ProcessHits , use wire->GetTOffSet instead of the fTzeroWire array when setting
the time.
5) THcDC.cxx
a) In TrackFit, replace fSigma with hit->GetWireSigma()
---
src/THcDC.cxx | 66 +++++++++++++++++++++++++++---------
src/THcDC.h | 1 -
src/THcDCHit.h | 1 +
src/THcDCWire.h | 7 ++--
src/THcDriftChamberPlane.cxx | 43 +++++++++++++++--------
src/THcDriftChamberPlane.h | 2 ++
6 files changed, 86 insertions(+), 34 deletions(-)
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 350453b..9ec08cf 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -357,6 +357,7 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
fReadoutTB[ip] = 0.0;
}
+
gHcParms->LoadParmValues((DBRequest*)&list,fPrefix);
//Set the default plane x,y positions to those of the chamber
@@ -958,7 +959,9 @@ void THcDC::TrackFit()
coords[ihit] = hit->GetCoord();
}
}
- }
+
+
+ } //end loop over hits
theDCTrack->SetNFree(theDCTrack->GetNHits() - NUM_FPRAY);
Double_t chi2 = dummychi2;
@@ -967,11 +970,19 @@ void THcDC::TrackFit()
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++) {
- TT[irayp] += (coords[ihit]*
- fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])
- /pow(fSigma[planes[ihit]],2);
- }
+ 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++) {
@@ -979,11 +990,16 @@ void THcDC::TrackFit()
if(jrayp<irayp) { // Symmetric
AA[irayp][jrayp] = AA[jrayp][irayp];
} else {
- 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);
- }
+ 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
}
}
}
@@ -1012,12 +1028,17 @@ void THcDC::TrackFit()
// 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);
+
+ 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/fSigma[planes[ihit]],2);
+ chi2 += pow(residual/hit->GetWireSigma(),2);
+
}
theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
@@ -1026,16 +1047,23 @@ void THcDC::TrackFit()
// 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(fSigma[planes[ihit]],2);
+ /pow(hit->GetWireSigma(),2);
+
}
}
}
@@ -1045,11 +1073,17 @@ void THcDC::TrackFit()
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(fSigma[planes[ihit]],2);
+ pow(hit->GetWireSigma(),2);
+
}
}
}
diff --git a/src/THcDC.h b/src/THcDC.h
index dd37637..c4fe61a 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -92,7 +92,6 @@ protected:
Int_t fdebugtrackprint;
Int_t fdebugprintdecodeddc;
Int_t fHMSStyleChambers;
-
Int_t fTDC_RefTimeCut;
UInt_t fNDCTracks;
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index 69eeb86..25536be 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -32,6 +32,7 @@ public:
// Get and Set Functions
THcDCWire* GetWire() const { return fWire; }
+ Double_t GetWireSigma() const { return fWire->GetSigma(); }
Int_t GetWireNum() const { return fWire->GetNum(); }
Int_t GetRawTime() const { return fRawTime; }
Int_t GetRawNoRefCorrTime() const { return fRawNoRefCorrTime; }
diff --git a/src/THcDCWire.h b/src/THcDCWire.h
index 5c89cd6..efbcc74 100644
--- a/src/THcDCWire.h
+++ b/src/THcDCWire.h
@@ -14,9 +14,9 @@ class THcDCWire : public TObject {
public:
- THcDCWire( Int_t num=0, Double_t pos=0.0, Double_t offset=0.0,
+ THcDCWire( Int_t num=0, Double_t pos=0.0, Double_t offset=0.0, Double_t sigma=0.0,
THcDCTimeToDistConv* ttd=NULL ) :
- fNum(num), fFlag(0), fPos(pos), fTOffset(offset), fTTDConv(ttd) {}
+ fNum(num), fFlag(0), fPos(pos), fTOffset(offset), fSigmaWire(sigma), fTTDConv(ttd) {}
virtual ~THcDCWire() {}
// Get and Set Functions
@@ -24,12 +24,14 @@ public:
Int_t GetFlag() const { return fFlag; }
Double_t GetPos() const { return fPos; }
Double_t GetTOffset() const { return fTOffset; }
+ Double_t GetSigma() const { return fSigmaWire; }
THcDCTimeToDistConv * GetTTDConv() { return fTTDConv; }
void SetNum (Int_t num) {fNum = num;}
void SetFlag (Int_t flag) {fFlag = flag;}
void SetPos (Double_t pos) { fPos = pos; }
void SetTOffset (Double_t tOffset){ fTOffset = tOffset; }
+ void SetSigma(Double_t tSigma){ fSigmaWire = tSigma; }
void SetTTDConv (THcDCTimeToDistConv * ttdConv){ fTTDConv = ttdConv;}
protected:
@@ -37,6 +39,7 @@ protected:
Int_t fFlag; //Flag for errors (e.g. Bad wire)
Double_t fPos; //Position within the plane
Double_t fTOffset; //Timing Offset
+ Double_t fSigmaWire; //Added SIgma per Wire --Carlos
THcDCTimeToDistConv* fTTDConv; //!Time to Distance Converter
private:
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index fa9bef1..bc09d4f 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -38,6 +38,7 @@ THcDriftChamberPlane::THcDriftChamberPlane( const char* name,
// Normal constructor with name and description
fHits = new TClonesArray("THcDCHit",100);
fWires = new TClonesArray("THcDCWire", 100);
+
fTTDConv = NULL;
fPlaneNum = planenum;
@@ -97,6 +98,7 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
Double_t DriftMapFirstBin;
Double_t DriftMapBinSize;
fUsingTzeroPerWire=0;
+ fUsingSigmaPerWire=0;
prefix[0]=tolower(GetParent()->GetPrefix()[0]);
prefix[1]='\0';
DBRequest list[]={
@@ -104,6 +106,7 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
{"drift1stbin", &DriftMapFirstBin, kDouble},
{"driftbinsz", &DriftMapBinSize, kDouble},
{"_using_tzero_per_wire", &fUsingTzeroPerWire, kInt,0,1},
+ {"_using_sigma_per_wire", &fUsingSigmaPerWire, kInt,0,1},
{0}
};
@@ -142,25 +145,35 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fNSperChan = fParent->GetNSperChan();
- if (fUsingTzeroPerWire==1) {
+
fTzeroWire = new Double_t [fNWires];
- DBRequest list3[]={
+ fSigmaWire = new Double_t [fNWires];
+
+
+ if (fUsingTzeroPerWire==1) {
+ DBRequest list3[]={
{Form("tzero%s",GetName()),fTzeroWire,kDouble,(UInt_t) fNWires},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list3,prefix);
- // printf(" using tzero per wire plane = %s nwires = %d \n",GetName(),fNWires);
- // for (Int_t iw=0;iw < fNWires;iw++) {
- // // printf("%d %f ",iw+1,fTzeroWire[iw]) ;
- // if ( iw!=0 && iw%8 == 0) printf("\n") ;
- // }
+
} else {
- fTzeroWire = new Double_t [fNWires];
for (Int_t iw=0;iw < fNWires;iw++) {
fTzeroWire[iw]=0.0;
}
}
+ if (fUsingSigmaPerWire==1) {
+ DBRequest list4[]={
+ {Form("wire_sigma%s",GetName()),fSigmaWire,kDouble,(UInt_t) fNWires},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list4,prefix);
+ } else {
+ for (Int_t iw=0;iw < fNWires;iw++) {
+ fSigmaWire[iw]=fSigma;
+ }
+ }
// Calculate Geometry Constants
// Do we want to move all this to the Chamber of DC Package leve
// as that is where these things will be needed?
@@ -207,7 +220,7 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
+ pow(hxpsi*fPsi0+hxchi*hchi0,2)
+ pow(hypsi*fPsi0+hychi*hchi0,2) );
if(z0 < 0.0) hphi0 = -hphi0;
-
+
Double_t denom2 = stubxpsi*stubychi - stubxchi*stubypsi;
// Why are there 4, but only 3 used?
@@ -243,10 +256,9 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
for (int i=0; i<nWires; i++) {
Double_t pos = fPitch*( (fWireOrder==0?(i+1):fNWires-i)
- fCentralWire) - fCenter;
- new((*fWires)[i]) THcDCWire( i+1, pos , 0.0, fTTDConv);
- //if( something < 0 ) wire->SetFlag(1);
+ new((*fWires)[i]) THcDCWire( i+1, pos , fTzeroWire[i], fSigmaWire[i], fTTDConv); //added fTzeroWire/fSigmaWire to be read in as fTOffset --Carlos
}
-
+
THaApparatus* app = GetApparatus();
const char* nm = "hod";
if( !app ||
@@ -255,9 +267,9 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
Warning(Here(here),"Hodoscope \"%s\" not found. "
"Event-by-event time offsets will NOT be used!!",nm);
}
-
+
return kOK;
-}
+}
//_____________________________________________________________________________
Int_t THcDriftChamberPlane::DefineVariables( EMode mode )
{
@@ -354,7 +366,8 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else if (rawtdc > fTdcWinMax) {
// Increment late count
} else {
- Double_t time = - rawtdc*fNSperChan + fPlaneTimeZero - fTzeroWire[wireNum-1]; // fNSperChan > 0 for 1877
+ Double_t time = - rawtdc*fNSperChan + fPlaneTimeZero - wire->GetTOffset(); // fNSperChan > 0 for 1877
+
new( (*fHits)[nextHit++] ) THcDCHit(wire, rawnorefcorrtdc,rawtdc, time, this);
break; // Take just the first hit in the time window
}
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index db23cba..b074cb4 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -89,6 +89,7 @@ protected:
Int_t fPlaneIndex; /* Index of this plane within it's chamber */
Int_t fChamberNum;
Int_t fUsingTzeroPerWire;
+ Int_t fUsingSigmaPerWire;
Int_t fNRawhits;
Int_t fNWires;
Int_t fTdcWinMin;
@@ -116,6 +117,7 @@ protected:
Double_t fNSperChan; /* TDC bin size */
Double_t* fTzeroWire;
+ Double_t* fSigmaWire;
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
--
GitLab