diff --git a/src/THcDC.cxx b/src/THcDC.cxx index 350453bb861d0c7e5ac7cbe9ff721f8d6acc2530..9ec08cf5317a1bffe394515bea224a8bed0f36e0 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 dd37637164cc9c7bd3fa99777dbb1fff3df444df..c4fe61a863ae8d79db9008d9abc99236672e7bd5 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 69eeb8659e1c67b91e4e2d7adb83a125094a98e0..25536be1676b6c53c5bde0608ab9e848aff458c8 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 5c89cd649427ab574e2388e1bac51f10b2d234dd..efbcc745f63e3768ed8a3fa3d96f76003633f9b4 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 fa9bef12879d5102ed1c84107bf1587c1da6d3f7..bc09d4fb699644834f4798f60093fecd5a743fd0 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 db23cbaeb2db77c537de3a89682f7a3b83807900..b074cb46134ee39e2d29fee49369eaac5719616b 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 );