Skip to content
Snippets Groups Projects
Commit 58c45eee authored by Stephen A. Wood's avatar Stephen A. Wood
Browse files

Add computation of drift distance from drift time.

  It has not been check that we are getting a good start time yet.
  Drift time and distance added to tree
  Change DC plane names from 1, 2, 3, ... to 1x1, 1y1, ...
    This is so that the parameters holding the time to distance maps
       can be found.  (The parameters are e.g. hwc1x1fract)
    Need to find a way not to have to hard code these plane names.
       Either use wire angles (alpha) or some kind of parameter name mapping
  Changed output.def to match new plane names
parent 3129a4bd
No related branches found
No related tags found
No related merge requests found
......@@ -25,4 +25,4 @@ raddeg=3.14159265/180
#include "PARAM/hhodo.param"
#include "PARAM/haero.param"
#include "PARAM/hdc.param"
#include "PARAM/hdriftmap.param"
......@@ -53,15 +53,15 @@ TH1F calposadc1 'HMS Cal ADC1' H.Cal.1pr.posadc1 150 -50 400
# Can we use variables for the constants. In CTP we used hdc_nwire(i)
#
TH1F hdc1x1_wm 'HDC 1X1 Wiremap' H.dc.1.tdchits 113 0.5 113.5
TH1F hdc1y1_wm 'HDC 1Y1 Wiremap' H.dc.2.tdchits 52 0.5 52.5
TH1F hdc1u1_wm 'HDC 1U1 Wiremap' H.dc.3.tdchits 107 0.5 107.5
TH1F hdc1v1_wm 'HDC 1V1 Wiremap' H.dc.4.tdchits 107 0.5 107.5
TH1F hdc1y2_wm 'HDC 1Y2 Wiremap' H.dc.5.tdchits 52 0.5 52.5
TH1F hdc1x2_wm 'HDC 1X2 Wiremap' H.dc.6.tdchits 113 0.5 113.5
TH1F hdc2x1_wm 'HDC 2X1 Wiremap' H.dc.7.tdchits 113 0.5 113.5
TH1F hdc2y1_wm 'HDC 2Y1 Wiremap' H.dc.8.tdchits 52 0.5 52.5
TH1F hdc2u1_wm 'HDC 2U1 Wiremap' H.dc.9.tdchits 107 0.5 107.5
TH1F hdc2v1_wm 'HDC 2V1 Wiremap' H.dc.10.tdchits 107 0.5 107.5
TH1F hdc2y2_wm 'HDC 2Y2 Wiremap' H.dc.11.tdchits 52 0.5 52.5
TH1F hdc2x2_wm 'HDC 2X2 Wiremap' H.dc.12.tdchits 113 0.5 113.5
TH1F hdc1x1_wm 'HDC 1X1 Wiremap' H.dc.1x1.tdchits 113 0.5 113.5
TH1F hdc1y1_wm 'HDC 1Y1 Wiremap' H.dc.1y1.tdchits 52 0.5 52.5
TH1F hdc1u1_wm 'HDC 1U1 Wiremap' H.dc.1u1.tdchits 107 0.5 107.5
TH1F hdc1v1_wm 'HDC 1V1 Wiremap' H.dc.1v1.tdchits 107 0.5 107.5
TH1F hdc1y2_wm 'HDC 1Y2 Wiremap' H.dc.1y2.tdchits 52 0.5 52.5
TH1F hdc1x2_wm 'HDC 1X2 Wiremap' H.dc.1x2.tdchits 113 0.5 113.5
TH1F hdc2x1_wm 'HDC 2X1 Wiremap' H.dc.2x1.tdchits 113 0.5 113.5
TH1F hdc2y1_wm 'HDC 2Y1 Wiremap' H.dc.2y1.tdchits 52 0.5 52.5
TH1F hdc2u1_wm 'HDC 2U1 Wiremap' H.dc.2u1.tdchits 107 0.5 107.5
TH1F hdc2v1_wm 'HDC 2V1 Wiremap' H.dc.2v1.tdchits 107 0.5 107.5
TH1F hdc2y2_wm 'HDC 2Y2 Wiremap' H.dc.2y2.tdchits 52 0.5 52.5
TH1F hdc2x2_wm 'HDC 2X2 Wiremap' H.dc.2x2.tdchits 113 0.5 113.5
......@@ -14,7 +14,7 @@ ClassImp(THcDCHit)
const Double_t THcDCHit::kBig = 1.e38; // Arbitrary large value
//_____________________________________________________________________________
Double_t THcDCHit::ConvertTimeToDist(Double_t slope)
Double_t THcDCHit::ConvertTimeToDist()
{
// Converts TDC time to drift distance
// Takes the (estimated) slope of the track as an argument
......@@ -24,7 +24,7 @@ Double_t THcDCHit::ConvertTimeToDist(Double_t slope)
if (ttdConv) {
// If a time to distance algorithm exists, use it to convert the TDC time
// to the drift distance
fDist = ttdConv->ConvertTimeToDist(fTime, slope, &fdDist);
fDist = ttdConv->ConvertTimeToDist(fTime);
return fDist;
}
......
......@@ -15,10 +15,12 @@ class THcDCHit : public TObject {
public:
THcDCHit( THcDCWire* wire=NULL, Int_t rawtime=0, Double_t time=0.0 ) :
fWire(wire), fRawTime(rawtime), fTime(time), fDist(0.0), ftrDist(kBig) {}
fWire(wire), fRawTime(rawtime), fTime(time), fDist(0.0), ftrDist(kBig) {
ConvertTimeToDist();
}
virtual ~THcDCHit() {}
virtual Double_t ConvertTimeToDist(Double_t slope);
virtual Double_t ConvertTimeToDist();
Int_t Compare ( const TObject* obj ) const;
Bool_t IsSortable () const { return kTRUE; }
......
///////////////////////////////////////////////////////////////////////////////
// //
// THcDCLookupTTDConv //
// THcDCLookupTTDConv //
// //
// Upon initialization needs to be given the lookup table for time //
// to distance conversion. //
// //
///////////////////////////////////////////////////////////////////////////////
......@@ -10,32 +13,21 @@ ClassImp(THcDCLookupTTDConv)
//______________________________________________________________________________
THcDCLookupTTDConv::THcDCLookupTTDConv()
THcDCLookupTTDConv::THcDCLookupTTDConv(Double_t T0, Double_t MaxDriftDistance,
Double_t BinSize, Int_t NumBins,
Double_t* Table) :
fT0(T0), fMaxDriftDistance(MaxDriftDistance), fBinSize(BinSize),
fNumBins(NumBins)
{
//Normal constructor
}
//______________________________________________________________________________
THcDCLookupTTDConv::THcDCLookupTTDConv( Double_t vel)
{
// Normal constructor
fDriftVel = vel;
fTable = new Double_t[fNumBins];
for(Int_t i=0;i<fNumBins;i++) {
fTable[i] = Table[i];
}
// TODO: This should be read from database!!
fA1tdcCor[0] = 2.12e-3;
fA1tdcCor[1] = 0.0;
fA1tdcCor[2] = 0.0;
fA1tdcCor[3] = 0.0;
fA2tdcCor[0] = -4.20e-4;
fA2tdcCor[1] = 1.3e-3;
fA2tdcCor[2] = 1.06e-4;
fA2tdcCor[3] = 0.0;
fdtime = 4.e-9; // 4ns -> 200 microns
}
//______________________________________________________________________________
THcDCLookupTTDConv::~THcDCLookupTTDConv()
{
......@@ -44,54 +36,25 @@ THcDCLookupTTDConv::~THcDCLookupTTDConv()
}
//______________________________________________________________________________
Double_t THcDCLookupTTDConv::ConvertTimeToDist(Double_t time,
Double_t tanTheta,
Double_t *ddist)
Double_t THcDCLookupTTDConv::ConvertTimeToDist(Double_t time)
{
// Drift Velocity in m/s
// time in s
// Return m
// printf("Converting Drift Time to Drift Distance!\n");
Double_t a1 = 0.0, a2 = 0.0;
// Find the values of a1 and a2 by evaluating the proper polynomials
// a = A_3 * x^3 + A_2 * x^2 + A_1 * x + A_0
tanTheta = 1.0 / tanTheta; // I assume this has to do w/ making the
// polynomial have the proper variable...
for (Int_t i = 3; i >= 1; i--) {
a1 = tanTheta * (a1 + fA1tdcCor[i]);
a2 = tanTheta * (a2 + fA2tdcCor[i]);
}
a1 += fA1tdcCor[0];
a2 += fA2tdcCor[0];
// printf("a1(%e) = %e\n", tanTheta, a1);
// printf("a2(%e) = %e\n", tanTheta, a2);
// ESPACE software includes corrections to the time for
// 1. Cluster t0 (offset applied to entire cluster)
// 2. Time of flight to scintillators
Double_t dist = fDriftVel * time;
Double_t unc = fDriftVel * fdtime; // watch uncertainty in the timing
if (dist < 0) {
// something screwy is going on
} else if (dist < a1 ) {
// dist = fDriftVel * time * (1 + 1 / (a1/a2 + 1));
dist *= ( 1 + a2 / a1);
unc *= ( 1 + a2 / a1);
} else {
dist += a2;
// Lookup in table
Int_t ib = (time-fT0)/fBinSize;
Double_t frac = 0;
if(ib >= 0 && ib+1 < fNumBins) {
Double_t tfrac = (time - (ib*fBinSize + fT0)) / fBinSize;
frac = fTable[ib]*(1-tfrac) + fTable[ib+1]*tfrac;
} else if (ib+1 >= fNumBins) {
frac = 1.0;
}
if (ddist) *ddist = unc;
// printf("D(%e) = %e\nUncorrected D = %e\n", time, dist, fDriftVel * time);
return dist;
}
Double_t drift_distance = fMaxDriftDistance * frac;
// Engine subtracts a hdc_card_delay from this. Seems
// to be zero in the PARAM files, bit is it always? Delay implies
// time? Whis is a time subtracted from a distance?
return(drift_distance);
}
////////////////////////////////////////////////////////////////////////////////
......@@ -13,31 +13,21 @@
class THcDCLookupTTDConv : public THcDCTimeToDistConv{
public:
THcDCLookupTTDConv( );
THcDCLookupTTDConv(Double_t vel);
THcDCLookupTTDConv(Double_t T0, Double_t MaxDriftDistance, Double_t BinSize,
Int_t NumBins, Double_t* Table);
virtual ~THcDCLookupTTDConv();
virtual Double_t ConvertTimeToDist(Double_t time, Double_t tanTheta,
Double_t *ddist=0);
virtual Double_t ConvertTimeToDist(Double_t time);
// Get and Set Functions
Double_t GetDriftVel() { return fDriftVel; }
void SetDriftVel(Double_t v) {fDriftVel = v; }
protected:
Double_t fDriftVel; // Drift velocity (m / s)
// Coefficients for a polynomial yielding correction parameters
// For now, hard code these values from db_eh845
// Eventually, this need to be read directly from the database
Double_t fA1tdcCor[4];
Double_t fA2tdcCor[4];
Double_t fdtime; // uncertainty in the measured time
Double_t fT0;
Double_t fBinSize;
Int_t fNumBins;
Double_t fMaxDriftDistance;
Double_t* fTable;
ClassDef(THcDCLookupTTDConv,0) // VDC Analytic TTD Conv class
};
......
......@@ -17,8 +17,8 @@ public:
THcDCTimeToDistConv() {}
virtual ~THcDCTimeToDistConv();
virtual Double_t ConvertTimeToDist(Double_t time, Double_t tanTheta,
Double_t *ddist=0) = 0;
virtual Double_t ConvertTimeToDist(Double_t time) = 0;
private:
THcDCTimeToDistConv( const THcDCTimeToDistConv& );
......
......@@ -73,13 +73,48 @@ void THcDriftChamber::Setup(const char* name, const char* description)
cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
fPlaneNames = new char* [fNPlanes];
for(Int_t i=0;i<fNPlanes;i++) {fPlaneNames[i] = new char[4];}
// Big hack needed because the drift map parameters are in parameters with
// names like hwc1x1fract. Need to do some kind of parameter name mapping
// or generate the names from the wire angle (alpha)
if(prefix[0] == 'h') {
strcpy(fPlaneNames[0],"1x1");
strcpy(fPlaneNames[1],"1y1");
strcpy(fPlaneNames[2],"1u1");
strcpy(fPlaneNames[3],"1v1");
strcpy(fPlaneNames[4],"1y2");
strcpy(fPlaneNames[5],"1x2");
strcpy(fPlaneNames[6],"2x1");
strcpy(fPlaneNames[7],"2y1");
strcpy(fPlaneNames[8],"2u1");
strcpy(fPlaneNames[9],"2v1");
strcpy(fPlaneNames[10],"2y2");
strcpy(fPlaneNames[11],"2x2");
} else if (prefix[0] == 's') {
strcpy(fPlaneNames[0],"1u1");
strcpy(fPlaneNames[1],"1u2");
strcpy(fPlaneNames[2],"1x1");
strcpy(fPlaneNames[3],"1x2");
strcpy(fPlaneNames[4],"1v1");
strcpy(fPlaneNames[5],"1v2");
strcpy(fPlaneNames[6],"2u1");
strcpy(fPlaneNames[7],"2u2");
strcpy(fPlaneNames[8],"2x1");
strcpy(fPlaneNames[9],"2x2");
strcpy(fPlaneNames[10],"2v1");
strcpy(fPlaneNames[11],"2v2");
} else {
cout << "Unknown Spectrometer Prefix '" << prefix << "' Guessing names" << endl;
for(Int_t i=0;i<fNPlanes;i++) {
sprintf(fPlaneNames[i],"%d",i+1);
}
}
char *desc = new char[strlen(description)+100];
fPlanes = new THcDriftChamberPlane* [fNPlanes];
for(Int_t i=0;i<fNPlanes;i++) {
fPlaneNames[i] = new char[5];
sprintf(fPlaneNames[i],"%d",i+1);
strcpy(desc, description);
strcat(desc, " Plane ");
strcat(desc, fPlaneNames[i]);
......
......@@ -80,11 +80,23 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
static const char* const here = "ReadDatabase()";
char prefix[2];
char parname[100];
Int_t NumDriftMapBins;
Double_t DriftMapFirstBin;
Double_t DriftMapBinSize;
Double_t DriftMap[1000];
prefix[0]=tolower(GetParent()->GetPrefix()[0]);
prefix[1]='\0';
DBRequest list[]={
{"driftbins", &NumDriftMapBins, kInt},
{"drift1stbin", &DriftMapFirstBin, kDouble},
{"driftbinsz", &DriftMapBinSize, kDouble},
{Form("wc%sfract",GetName()),DriftMap,kDouble,1000},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
// Retrieve parameters we need
// Retrieve parameters we need from parent class
THcDriftChamber* fParent;
fParent = (THcDriftChamber*) GetParent();
......@@ -103,7 +115,8 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
cout << fPlaneNum << " " << fNWires << endl;
fTTDConv = new THcDCLookupTTDConv();// Need to pass the lookup table
fTTDConv = new THcDCLookupTTDConv(DriftMapFirstBin,fPitch/2,DriftMapBinSize,
NumDriftMapBins,DriftMap);
Int_t nWires = fParent->GetNWires(fPlaneNum);
// For HMS, wire numbers start with one, but arrays start with zero.
......@@ -142,6 +155,10 @@ Int_t THcDriftChamberPlane::DefineVariables( EMode mode )
"fHits.THcDCHit.GetWireNum()"},
{"rawtdc", "Raw TDC Values",
"fHits.THcDCHit.GetRawTime()"},
{"time","Drift times",
"fHits.THcDCHit.GetTime()"},
{"dist","Drift distancess",
"fHits.THcDCHit.GetDist()"},
{ 0 }
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment