diff --git a/setup.csh b/setup.csh
index d7df18c89cf06ac4bea58b32ff02b2cd2fbf27d2..8bbc4063544bafc591e3e7ca94f5b2d0b065eaf1 100644
--- a/setup.csh
+++ b/setup.csh
@@ -2,17 +2,29 @@
set called=($_)
if ("$called" != "") then
- set scriptdir=$called[2]
- set MYDIR=`dirname $scriptdir`
- set MYDIR=`cd $MYDIR && pwd` # ensure absolute path
+set scriptdir=$called[2]
+ echo $scriptdir
+# set MYDIR=`dirname $scriptdir`
+# set MYDIR=`cd $MYDIR && pwd` # ensure absolute path
+
+ set MYDIR=`pwd` # ensure absolute path
+
+# echo $MYDIR
+# echo "Are we here"
else
- set scriptdir=$1
+set scriptdir=$1
set MYDIR=$scriptdir
endif
+
+
setenv ANALYZER $MYDIR/podd
setenv HCANALYZER $MYDIR
+
# Check if LD_LIBRARY_PATH is defined
if ( ! ($?LD_LIBRARY_PATH) ) then
- setenv LD_LIBRARY_PATH ""
+setenv LD_LIBRARY_PATH ""
endif
setenv LD_LIBRARY_PATH "${LD_LIBRARY_PATH}:${ANALYZER}:${HCANALYZER}"
+
+echo $ANALYZER
+echo $HCANALYZER
diff --git a/src/THcHallCSpectrometer.cxx b/src/THcHallCSpectrometer.cxx
index 7ce0bfcdf91041a520451b440999af641ca948c3..073c516d6966a093bbf46166d667e1b0c23ed972 100644
--- a/src/THcHallCSpectrometer.cxx
+++ b/src/THcHallCSpectrometer.cxx
@@ -40,6 +40,19 @@
// 'reference distance', corresponding to the pathlength correction
// matrix.
//
+//
+// Golden track using scin. Zafar Ahmed. August 19 2014
+// Goldent track is moved to THcHallCSpectrometer::TrackCalc()
+// if fSelUsingScin == 0 then golden track is calculated just
+// like podd. i.e. it is the first track with minimum chi2/ndf
+// with sorting ON
+//
+// if fSelUsingScin == 1 then golden track is calculetd just like
+// engine/HTRACKING/h_select_best_track_using_scin.h. This method
+// gives the best track with minimum value of chi2/ndf but with
+// additional cuts on the tracks. These cuts are on dedx, beta
+// and on energy.
+//
//////////////////////////////////////////////////////////////////////////
#include "THcHallCSpectrometer.h"
@@ -51,11 +64,21 @@
#include "THaTriggerTime.h"
#include "TMath.h"
#include "TList.h"
+
+#include "THcRawShowerHit.h"
+#include "THcSignalHit.h"
#include "THcShower.h"
+#include "THcHitList.h"
+#include "THcHodoscope.h"
+#include <vector>
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
#include <iostream>
#include <fstream>
+using std::vector;
using namespace std;
//_____________________________________________________________________________
@@ -74,6 +97,13 @@ THcHallCSpectrometer::THcHallCSpectrometer( const char* name, const char* descri
THcHallCSpectrometer::~THcHallCSpectrometer()
{
// Destructor
+
+ delete [] fX2D; fX2D = NULL; // Ahmed
+ delete [] fY2D; fY2D = NULL; // Ahmed
+
+ delete [] f2XHits; f2XHits = NULL; // Ahmed
+ delete [] f2YHits; f2YHits = NULL; // Ahmed
+
}
//_____________________________________________________________________________
@@ -123,6 +153,33 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
cout << "In THcHallCSpectrometer::ReadDatabase()" << endl;
#endif
+ MAXHODHITS = 53;
+
+ fX2D = new Double_t [MAXHODHITS];
+ fY2D = new Double_t [MAXHODHITS];
+
+ f2XHits = new Int_t [16];
+ f2YHits = new Int_t [16];
+
+ // --------------- To get energy from THcShower ----------------------
+
+ const char* detector_name = "hod";
+ //THaApparatus* app = GetDetector();
+ THaDetector* det = GetDetector("hod");
+ // THaDetector* det = app->GetDetector( shower_detector_name );
+
+ if( !dynamic_cast<THcHodoscope*>(det) ) {
+ Error("THcHallCSpectrometer", "Cannot find hodoscope detector %s",
+ detector_name );
+ return fStatus = kInitError;
+ }
+
+ fHodo = static_cast<THcHodoscope*>(det); // fHodo is a membervariable
+ // fShower = static_cast<THcShower*>(det); // fShower is a membervariable
+
+ // --------------- To get energy from THcShower ----------------------
+
+
// Get the matrix element filename from the variable store
// Read in the matrix
@@ -137,19 +194,34 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
string reconCoeffFilename;
DBRequest list[]={
- {"_recon_coeff_filename",&reconCoeffFilename,kString},
- {"theta_offset",&fThetaOffset,kDouble},
- {"phi_offset",&fPhiOffset,kDouble},
- {"delta_offset",&fDeltaOffset,kDouble},
- {"thetacentral_offset",&fThetaCentralOffset,kDouble},
- {"_oopcentral_offset",&fOopCentralOffset,kDouble},
- {"pcentral_offset",&fPCentralOffset,kDouble},
- {"pcentral",&fPCentral,kDouble},
- {"theta_lab",&fTheta_lab,kDouble},
+ {"_recon_coeff_filename", &reconCoeffFilename, kString },
+ {"theta_offset", &fThetaOffset, kDouble },
+ {"phi_offset", &fPhiOffset, kDouble },
+ {"delta_offset", &fDeltaOffset, kDouble },
+ {"thetacentral_offset", &fThetaCentralOffset, kDouble },
+ {"_oopcentral_offset", &fOopCentralOffset, kDouble },
+ {"pcentral_offset", &fPCentralOffset, kDouble },
+ {"pcentral", &fPCentral, kDouble },
+ {"theta_lab", &fTheta_lab, kDouble },
+ {"sel_using_scin", &fSelUsingScin, kInt, 0, 1},
+ {"sel_ndegreesmin", &fSelNDegreesMin, kDouble, 0, 1},
+ {"sel_dedx1min", &fSeldEdX1Min, kDouble, 0, 1},
+ {"sel_dedx1max", &fSeldEdX1Max, kDouble, 0, 1},
+ {"sel_betamin", &fSelBetaMin, kDouble, 0, 1},
+ {"sel_betamax", &fSelBetaMax, kDouble, 0, 1},
+ {"sel_etmin", &fSelEtMin, kDouble, 0, 1},
+ {"sel_etmax", &fSelEtMax, kDouble, 0, 1},
+ {"hodo_num_planes", &fNPlanes, kInt },
+ {"scin_2x_zpos", &fScin2XZpos, kDouble, 0, 1},
+ {"scin_2x_dzpos", &fScin2XdZpos, kDouble, 0, 1},
+ {"scin_2y_zpos", &fScin2YZpos, kDouble, 0, 1},
+ {"scin_2y_dzpos", &fScin2YdZpos, kDouble, 0, 1},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
//
+
+ cout << "\n\n\nhodo planes = " << fNPlanes << endl;
cout << " fPcentral = " << fPCentral << " " <<fPCentralOffset << endl;
cout << " fThate_lab = " << fTheta_lab << " " <<fThetaCentralOffset << endl;
fPCentral= fPCentral*(1.+fPCentralOffset/100.);
@@ -228,6 +300,8 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
// In Hall C, we do the target traceback here since the traceback should
// not depend on which tracking detectors are used.
+ fNtracks = tracks.GetLast()+1;
+
for (Int_t it=0;it<tracks.GetLast()+1;it++) {
THaTrack* track = static_cast<THaTrack*>( tracks[it] );
@@ -281,6 +355,10 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
track->SetMomentum(fPCentral*(1+track->GetDp()/100.0));
}
+
+ // ------------------ Moving it to TrackCalc --------------------
+
+ /*
// If enabled, sort the tracks by chi2/ndof
if( GetTrSorting() )
fTracks->Sort();
@@ -305,14 +383,215 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
} else
fGoldenTrack = NULL;
+ */
+ // ------------------ Moving it to TrackCalc --------------------
+
return 0;
}
//_____________________________________________________________________________
Int_t THcHallCSpectrometer::TrackCalc()
{
- // Additioal track calculations. At present, we only calculate beta here.
+ if ( fSelUsingScin == 0 ) {
+
+ if( GetTrSorting() )
+ fTracks->Sort();
+
+ // Find the "Golden Track".
+ // if( GetNTracks() > 0 ) {
+ if( fNtracks > 0 ) {
+ // Select first track in the array. If there is more than one track
+ // and track sorting is enabled, then this is the best fit track
+ // (smallest chi2/ndof). Otherwise, it is the track with the best
+ // geometrical match (smallest residuals) between the U/V clusters
+ // in the upper and lower VDCs (old behavior).
+ //
+ // Chi2/dof is a well-defined quantity, and the track selected in this
+ // way is immediately physically meaningful. The geometrical match
+ // criterion is mathematically less well defined and not usually used
+ // in track reconstruction. Hence, chi2 sortiing is preferable, albeit
+ // obviously slower.
+
+ fGoldenTrack = static_cast<THaTrack*>( fTracks->At(0) );
+ fTrkIfo = *fGoldenTrack;
+ fTrk = fGoldenTrack;
+ } else
+ fGoldenTrack = NULL;
+ }
+
+ if ( fSelUsingScin == 1 ){
+ if( fNtracks > 0 ) {
+
+ Double_t fY2Dmin, fX2Dmin, fZap, ft, fChi2PerDeg; //, fShowerEnergy;
+ Double_t fHitPos4, fHitPos3, fHitDist3, fHitDist4; //, fChi2Min;
+ Int_t i, j, itrack, ip, ihit; //, fGoodTimeIndex = -1;
+ Int_t fHitCnt4, fHitCnt3, fRawIndex, fGoodRawPad;
+
+ fChi2Min = 10000000000.0; fGoodTrack = -1; fY2Dmin = 100.;
+ fX2Dmin = 100.; fZap = 0.;
+
+ fRawIndex = -1;
+ for ( itrack = 0; itrack < fNtracks; itrack++ ){
+
+ THaTrack* goodTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
+ if (!goodTrack) return -1;
+
+ if ( goodTrack->GetNDoF() > fSelNDegreesMin ){
+ fChi2PerDeg = goodTrack->GetChi2() / goodTrack->GetNDoF();
+
+ if( ( goodTrack->GetDedx() > fSeldEdX1Min ) &&
+ ( goodTrack->GetDedx() < fSeldEdX1Max ) &&
+ ( goodTrack->GetBeta() > fSelBetaMin ) &&
+ ( goodTrack->GetBeta() < fSelBetaMax ) &&
+ ( goodTrack->GetEnergy() > fSelEtMin ) &&
+ ( goodTrack->GetEnergy() < fSelEtMax ) )
+ {
+
+ for ( j = 0; j < 16; j++ ){
+ f2XHits[j] = -1;
+ f2YHits[j] = -1;
+ }
+
+ for ( ip = 0; ip < fNPlanes; ip++ ){
+ for ( ihit = 0; ihit < fHodo->GetNScinHits(ip); ihit++ ){
+ fRawIndex ++;
+
+ // fGoodRawPad = fHodo->GetGoodRawPad(fRawIndex)-1;
+ fGoodRawPad = fHodo->GetGoodRawPad(fRawIndex);
+
+ if ( ip == 2 ){
+ f2XHits[fGoodRawPad] = 0;
+ }
+
+ if ( ip == 3 ){
+ f2YHits[fGoodRawPad] = 0;
+
+ }
+
+ } // loop over hits of a plane
+ } // loop over planes
+
+ fHitPos4 = goodTrack->GetY() + goodTrack->GetPhi() * ( fScin2YZpos + 0.5 * fScin2YdZpos );
+ fHitCnt4 = TMath::Nint( ( fHodo->GetHodoCenter4() - fHitPos4 ) / fHodo->GetScin2YSpacing() ) + 1;
+ fHitCnt4 = TMath::Max( TMath::Min(fHitCnt4,TMath::Nint(10) ) , 1); // scin_2y_nr = 10
+ fHitDist4 = fHitPos4 - ( fHodo->GetHodoCenter4() - fHodo->GetScin2YSpacing() * ( fHitCnt4 - 1 ) );
+
+ //----------------------------------------------------------------
+
+ if ( fNtracks > 1 ){ // Plane 4
+ fZap = 0.;
+ ft = 0;
+
+ for ( i = 0; i < 10; i++ ){
+ if ( f2YHits[i] == 0 ) {
+ fY2D[itrack] = TMath::Abs(fHitCnt4-i-1);
+ ft ++;
+
+ if ( ft == 1 ) fZap = fY2D[itrack];
+ if ( ( ft == 2 ) && ( fY2D[itrack] < fZap ) ) fZap = fY2D[itrack];
+ if ( ( ft == 3 ) && ( fY2D[itrack] < fZap ) ) fZap = fY2D[itrack];
+ if ( ( ft == 4 ) && ( fY2D[itrack] < fZap ) ) fZap = fY2D[itrack];
+ if ( ( ft == 5 ) && ( fY2D[itrack] < fZap ) ) fZap = fY2D[itrack];
+ if ( ( ft == 6 ) && ( fY2D[itrack] < fZap ) ) fZap = fY2D[itrack];
+
+ } // condition for fHodScinHit[4][i]
+ } // loop over 10
+ fY2D[itrack] = fZap;
+ } // condition for track. Plane 4
+
+ if ( fNtracks == 1 ) fY2D[itrack] = 0.;
+
+ fHitPos3 = goodTrack->GetX() + goodTrack->GetTheta() * ( fScin2XZpos + 0.5 * fScin2XdZpos );
+ fHitCnt3 = TMath::Nint( ( fHitPos3 - fHodo->GetHodoCenter3() ) / fHodo->GetScin2XSpacing() ) + 1;
+ fHitCnt3 = TMath::Max( TMath::Min(fHitCnt3,TMath::Nint(16) ) , 1); // scin_2x_nr = 16
+ fHitDist3 = fHitPos3 - ( fHodo->GetScin2XSpacing() * ( fHitCnt3 - 1 ) + fHodo->GetHodoCenter3() );
+
+ //----------------------------------------------------------------
+
+ if ( fNtracks > 1 ){ // Plane 3 (2X)
+ fZap = 0.;
+ ft = 0;
+ for ( i = 0; i < 16; i++ ){
+ if ( f2XHits[i] == 0 ) {
+ fX2D[itrack] = TMath::Abs(fHitCnt3-i-1);
+
+ ft ++;
+ if ( ft == 1 ) fZap = fX2D[itrack];
+ if ( ( ft == 2 ) && ( fX2D[itrack] < fZap ) ) fZap = fX2D[itrack];
+ if ( ( ft == 3 ) && ( fX2D[itrack] < fZap ) ) fZap = fX2D[itrack];
+ if ( ( ft == 4 ) && ( fX2D[itrack] < fZap ) ) fZap = fX2D[itrack];
+ if ( ( ft == 5 ) && ( fX2D[itrack] < fZap ) ) fZap = fX2D[itrack];
+ if ( ( ft == 6 ) && ( fX2D[itrack] < fZap ) ) fZap = fX2D[itrack];
+
+ } // condition for fHodScinHit[4][i]
+ } // loop over 16
+ fX2D[itrack] = fZap;
+ } // condition for track. Plane 3
+
+ //----------------------------------------------------------------
+
+ if ( fNtracks == 1 )
+ fX2D[itrack] = 0.;
+
+ if ( fY2D[itrack] <= fY2Dmin ) {
+ if ( fY2D[itrack] < fY2Dmin ) {
+ fX2Dmin = 100.;
+ fChi2Min = 10000000000.;
+ } // fY2D min
+
+ if ( fX2D[itrack] <= fX2Dmin ){
+ if ( fX2D[itrack] < fX2Dmin ){
+ fChi2Min = 10000000000.0;
+ } // condition fX2D
+ if ( fChi2PerDeg < fChi2Min ){
+
+ fGoodTrack = itrack; // fGoodTrack = itrack
+ fY2Dmin = fY2D[itrack];
+ fX2Dmin = fX2D[itrack];
+ fChi2Min = fChi2PerDeg;
+
+ fGoldenTrack = static_cast<THaTrack*>( fTracks->At( fGoodTrack ) );
+ fTrkIfo = *fGoldenTrack;
+ fTrk = fGoldenTrack;
+ }
+ } // condition fX2D
+ } // condition fY2D
+ } // conditions for dedx, beta and trac energy
+ } // confition for fNFreeFP greater than fSelNDegreesMin
+ } // loop over tracks
+
+ if ( fGoodTrack == -1 ){
+
+ fChi2Min = 10000000000.0;
+ for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
+ THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
+ if (!aTrack) return -1;
+
+ if ( aTrack->GetNDoF() > fSelNDegreesMin ){
+ fChi2PerDeg = aTrack->GetChi2() / aTrack->GetNDoF();
+
+ if ( fChi2PerDeg < fChi2Min ){
+
+ fGoodTrack = iitrack;
+ fChi2Min = fChi2PerDeg;
+
+ fGoldenTrack = static_cast<THaTrack*>( fTracks->At(fGoodTrack) );
+ fTrkIfo = *fGoldenTrack;
+ fTrk = fGoldenTrack;
+
+ }
+ }
+ } // loop over trakcs
+
+ } // Condition for fNtrack > 0
+
+ } else // condtion fSelUsingScin == 1
+ fGoldenTrack = NULL;
+
+ }
+
+
return TrackTimes( fTracks );
}
@@ -323,6 +602,7 @@ Int_t THcHallCSpectrometer::TrackTimes( TClonesArray* Tracks ) {
//
// To be useful, a meaningful timing resolution should be assigned
// to each Scintillator object (part of the database).
+
return 0;
}
diff --git a/src/THcHallCSpectrometer.h b/src/THcHallCSpectrometer.h
index f906546ec9c187e94ec3f607d90b2089420ae0ad..0bbb5229ec79bac0fcb47e11798e7e8f36f5e1d1 100644
--- a/src/THcHallCSpectrometer.h
+++ b/src/THcHallCSpectrometer.h
@@ -9,6 +9,29 @@
#include "THaSpectrometer.h"
+#include <vector>
+
+#include "TClonesArray.h"
+#include "THaNonTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcHodoscopeHit.h"
+#include "THcScintillatorPlane.h"
+#include "THcShower.h"
+
+//#include "THaTrackingDetector.h"
+//#include "THcHitList.h"
+#include "THcRawDCHit.h"
+#include "THcSpacePoint.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
+#include "TMath.h"
+
+#include "THaSubDetector.h"
+#include "TClonesArray.h"
+#include <iostream>
+#include <fstream>
+
+
//class THaScintillator;
class THcHallCSpectrometer : public THaSpectrometer {
@@ -30,6 +53,40 @@ public:
protected:
void InitializeReconstruction();
+ Int_t MAXHODHITS;
+
+ Int_t fGoodTrack;
+ Int_t fSelUsingScin;
+ Int_t fNPlanes;
+ Int_t fNtracks;
+
+ Int_t* f2XHits;
+ Int_t* f2YHits;
+
+ Double_t* fX2D;
+ Double_t* fY2D;
+
+ Double_t fChi2Min;
+ Double_t fSelNDegreesMin;
+ Double_t fSeldEdX1Min;
+ Double_t fSeldEdX1Max;
+ Double_t fSelBetaMin;
+ Double_t fSelBetaMax;
+ Double_t fSelEtMin;
+ Double_t fSelEtMax;
+ Double_t fScin2XZpos;
+ Double_t fScin2XdZpos;
+ Double_t fScin2YZpos;
+ Double_t fScin2YdZpos;
+
+ Double_t fHodoCenter4, fHodoCenter3;
+ Double_t fScin2YSpacing, fScin2XSpacing;
+
+ // Int_t** fHodScinHit; // [4] Array
+
+ THcShower* fShower;
+ THcHodoscope* fHodo;
+
// Should look at the ThaMatrixElement class in THaVDC.h for better way
// to store matrix element data
#define fMaxReconElements 1000
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index e9a8b1d032a92e1c64f717e576395b2b6a00dbb0..e0e3933900c16fe7a6bd57007b60b4d7e80e6b41 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -15,6 +15,13 @@
///////////////////////////////////////////////////////////////////////////////
#include "THcSignalHit.h"
+#include "THcShower.h"
+
+#include "THcHitList.h"
+#include "THcRawShowerHit.h"
+#include "TClass.h"
+#include "math.h"
+#include "THaSubDetector.h"
#include "THcHodoscope.h"
#include "THaEvData.h"
@@ -104,6 +111,9 @@ void THcHodoscope::Setup(const char* name, const char* description)
fPlaneNames[i] = new char[plane_names[i].length()];
strcpy(fPlaneNames[i], plane_names[i].c_str());
}
+
+ //myShower = new THcShower("cal", "Shower" );
+
/* fPlaneNames = new char* [fNPlanes];
for(Int_t i=0;i<fNPlanes;i++) {fPlaneNames[i] = new char[3];}
// Should get the plane names from parameters.
@@ -136,6 +146,23 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
// But it needs to happen before the sub detectors are initialized
// so that they can get the pointer to the hitlist.
+ // --------------- To get energy from THcShower ----------------------
+
+ const char* shower_detector_name = "cal";
+ THaApparatus* app = GetApparatus();
+ THaDetector* det = app->GetDetector( shower_detector_name );
+
+ if( !dynamic_cast<THcShower*>(det) ) {
+ Error("THcHodoscope", "Cannot find shower detector %s",
+ shower_detector_name );
+ return fStatus = kInitError;
+ }
+
+ fShower = static_cast<THcShower*>(det); // fShower is a membervariable
+
+ // --------------- To get energy from THcShower ----------------------
+
+
InitHitList(fDetMap, "THcHodoscopeHit", 100);
EStatus status;
@@ -165,6 +192,50 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
return kInitError;
}
+ fKeepPos = new Bool_t [MAXHODHITS];
+ fKeepNeg = new Bool_t [MAXHODHITS];
+ fGoodTDCPos = new Bool_t [MAXHODHITS];
+ fGoodTDCNeg = new Bool_t [MAXHODHITS];
+
+ fGoodRawPad = new Int_t [MAXHODHITS];
+ fHitPaddle = new Int_t [MAXHODHITS];
+ fNScinHit = new Int_t [MAXHODHITS];
+ fNPmtHit = new Int_t [MAXHODHITS];
+ fTimeHist = new Int_t [200];
+
+ fTimeAtFP = new Double_t [MAXHODHITS];
+ fScinSigma = new Double_t [MAXHODHITS];
+ fGoodScinTime = new Double_t [MAXHODHITS];
+ fScinTime = new Double_t [MAXHODHITS];
+ fTime = new Double_t [MAXHODHITS];
+ adcPh = new Double_t [MAXHODHITS];
+ fPath = new Double_t [MAXHODHITS];
+ fTimePos = new Double_t [MAXHODHITS];
+ fTimeNeg = new Double_t [MAXHODHITS];
+ fScinTimefp = new Double_t [MAXHODHITS];
+ fScinPosTime = new Double_t [MAXHODHITS];
+ fScinNegTime = new Double_t [MAXHODHITS];
+
+ fNScinHits = new Int_t [fNPlanes];
+ fGoodPlaneTime = new Bool_t [fNPlanes];
+ fNPlaneTime = new Int_t [fNPlanes];
+ fSumPlaneTime = new Double_t [fNPlanes];
+
+ // Double_t fHitCnt4 = 0., fHitCnt3 = 0.;
+
+ Int_t m = 0;
+
+ fdEdX = new Double_t*[MAXHODHITS];
+ for ( m = 0; m < MAXHODHITS; m++ ){
+ fdEdX[m] = new Double_t[MAXHODHITS];
+ }
+
+ fScinHit = new Double_t*[fNPlanes];
+ for ( m = 0; m < fNPlanes; m++ ){
+ fScinHit[m] = new Double_t[fNPaddle[0]];
+ }
+
+
return fStatus = kOK;
}
//_____________________________________________________________________________
@@ -369,32 +440,58 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin];
fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin];
+
+
prefix[1]='\0';
DBRequest list[]={
- {"start_time_center", &fStartTimeCenter, kDouble},
- {"start_time_slop", &fStartTimeSlop, kDouble},
- {"scin_tdc_to_time", &fScinTdcToTime, kDouble},
- {"scin_tdc_min", &fScinTdcMin, kDouble},
- {"scin_tdc_max", &fScinTdcMax, kDouble},
- {"tof_tolerance", &fTofTolerance, kDouble,0,1},
- {"pathlength_central", &fPathLengthCentral, kDouble},
- {"hodo_vel_light",&fHodoVelLight[0],kDouble,fMaxHodoScin},
- {"hodo_pos_sigma",&fHodoPosSigma[0],kDouble,fMaxHodoScin},
- {"hodo_neg_sigma",&fHodoNegSigma[0],kDouble,fMaxHodoScin},
- {"hodo_pos_minph",&fHodoPosMinPh[0],kDouble,fMaxHodoScin},
- {"hodo_neg_minph",&fHodoNegMinPh[0],kDouble,fMaxHodoScin},
- {"hodo_pos_phc_coeff",&fHodoPosPhcCoeff[0],kDouble,fMaxHodoScin},
- {"hodo_neg_phc_coeff",&fHodoNegPhcCoeff[0],kDouble,fMaxHodoScin},
- {"hodo_pos_time_offset",&fHodoPosTimeOffset[0],kDouble,fMaxHodoScin},
- {"hodo_neg_time_offset",&fHodoNegTimeOffset[0],kDouble,fMaxHodoScin},
- {"hodo_pos_ped_limit",&fHodoPosPedLimit[0],kInt,fMaxHodoScin},
- {"hodo_neg_ped_limit",&fHodoNegPedLimit[0],kInt,fMaxHodoScin},
- {"tofusinginvadc",&fTofUsingInvAdc,kInt,0,1},
+ // {"scin_2x_zpos", &fScin2XZpos, kDouble, 0, 1},
+ // {"scin_2x_dzpos", &fScin2XdZpos, kDouble, 0, 1},
+ // {"scin_2y_zpos", &fScin2YZpos, kDouble, 0, 1},
+ // {"scin_2y_dzpos", &fScin2YdZpos, kDouble, 0, 1},
+ // {"sel_betamin", &fSelBetaMin, kDouble, 0, 1},
+ // {"sel_dedx1min", &fSeldEdX1Min, kDouble, 0, 1},
+ // {"sel_dedx1max", &fSeldEdX1Max, kDouble, 0, 1},
+ // {"sel_using_scin", &fSelUsingScin, kInt, 0, 1},
+ // {"sel_ndegreesmin", &fSelNDegreesMin, kDouble, 0, 1},
+ {"start_time_center", &fStartTimeCenter, kDouble},
+ {"start_time_slop", &fStartTimeSlop, kDouble},
+ {"scin_tdc_to_time", &fScinTdcToTime, kDouble},
+ {"scin_tdc_min", &fScinTdcMin, kDouble},
+ {"scin_tdc_max", &fScinTdcMax, kDouble},
+ {"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
+ {"pathlength_central", &fPathLengthCentral, kDouble},
+ {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_minph", &fHodoPosMinPh[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_minph", &fHodoNegMinPh[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_phc_coeff", &fHodoPosPhcCoeff[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_phc_coeff", &fHodoNegPhcCoeff[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_time_offset", &fHodoPosTimeOffset[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_time_offset", &fHodoNegTimeOffset[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin},
+ {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin},
+ {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
{0}
};
fTofUsingInvAdc = 0; // Default if not defined
fTofTolerance = 3.0; // Default if not defined
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+
+ cout << "\n\n\n\n\n\nPaddles1x = " << fNPaddle[0]
+ << "\nscin_2y_zpos = " << fScin2YZpos
+ << "\nscin_2y_dzpos = " << fScin2YdZpos
+ << endl;
+ // << "\ndedx max = " << fSeldEdX1Max
+ // << "\nbeta min = " << fSelBetaMin
+ // << "\nbeta max = " << fSelBetaMax
+ // << "\net min = " << fSelEtMin
+ // << "\net max = " << fSelEtMax
+
+ // << endl;
+
+
if (fTofUsingInvAdc) {
DBRequest list2[]={
{"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin},
@@ -480,7 +577,7 @@ Int_t THcHodoscope::DefineVariables( EMode mode )
// { "y_adc", "y-position from amplitudes (m)", "fYa" },
// { "time", "Time of hit at plane (s)", "fTime" },
// { "dtime", "Est. uncertainty of time (s)", "fdTime" },
- // { "dedx", "dEdX-like deposited in paddle", "fAmpl" },
+ { "dedx", "dEdX-like deposited in paddle", "fdEdX" },
// { "troff", "Trigger offset for paddles", "fTrigOff"},
// { "trn", "Number of tracks for hits", "GetNTracks()" },
// { "trx", "x-position of track in det plane", "fTrackProj.THaTrackProj.fX" },
@@ -517,7 +614,17 @@ THcHodoscope::~THcHodoscope()
void THcHodoscope::DeleteArrays()
{
// Delete member arrays. Used by destructor.
-
+ Int_t k;
+ for( k = 0; k < MAXHODHITS; k++){
+ delete [] fdEdX[k];
+ }
+ delete [] fdEdX;
+
+ for( k = 0; k < fNPlanes; k++){
+ delete [] fScinHit[k];
+ }
+ delete [] fScinHit;
+
delete [] fNPaddle; fNPaddle = NULL;
delete [] fHodoVelLight; fHodoVelLight = NULL;
delete [] fHodoPosSigma; fHodoPosSigma = NULL;
@@ -536,12 +643,14 @@ void THcHodoscope::DeleteArrays()
delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL;
delete [] fHodoPosInvAdcAdc; fHodoPosInvAdcAdc = NULL;
+ delete [] fGoodRawPad; fGoodRawPad = NULL; // Ahmed
delete [] fHitPaddle; fHitPaddle = NULL; // Ahmed
delete [] fNScinHit; fNScinHit = NULL; // Ahmed
delete [] fNPmtHit; fNPmtHit = NULL; // Ahmed
delete [] fTimeHist; fTimeHist = NULL; // Ahmed
delete [] fTimeAtFP; fTimeAtFP = NULL; // Ahmed
+
delete [] fScinSigma; fScinSigma = NULL; // Ahmed
delete [] fGoodScinTime; fGoodScinTime = NULL; // Ahmed
delete [] fScinTime; fScinTime = NULL; // Ahmed
@@ -561,6 +670,8 @@ void THcHodoscope::DeleteArrays()
delete [] fNPlaneTime; fNPlaneTime = NULL; // Ahmed
delete [] fSumPlaneTime; fSumPlaneTime = NULL; // Ahmed
+ delete [] fNScinHits; fNScinHits = NULL; // Ahmed
+
// delete [] fSpacing; fSpacing = NULL;
//delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete?
@@ -616,8 +727,10 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
//
// GN: print event number so we can cross-check with engine
// if (evdata.GetEvNum()>1000)
- // cout <<"\nhcana event = " << evdata.GetEvNum()<<endl;
+ // cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;
+ fCheckEvent = evdata.GetEvNum();
+
if(gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
@@ -648,6 +761,15 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
fNfptimes=0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if ( ip == 2 ){
+ fHodoCenter3 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ fScin2XSpacing = fPlanes[ip]->GetSpacing();
+ }
+ if ( ip == 3 ){
+ fHodoCenter4 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ fScin2YSpacing = fPlanes[ip]->GetSpacing();
+ }
+
// if ( !fPlanes[ip] ) // Ahmed
// return -1; // Ahmed
@@ -703,47 +825,34 @@ Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle,
//_____________________________________________________________________________
Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
+{
+
+ ApplyCorrections();
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
{
Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
- Int_t fPaddle = 0, fIndex, k, fJMax, fMaxHit, ip, ihit;
- Int_t fNfpTime, fSumfpTime, fNScinHits;
+ Int_t fPaddle = 0, fIndex, k, fJMax, fMaxHit, ip, ihit, itrack;
+ Int_t fNfpTime, fSumfpTime, fRawIndex = -1; //, fNScinHits[fNPlanes];
Double_t fScinTrnsCoord, fScinLongCoord, fScinCenter;
Double_t fP, fBetaP, fXcoord, fYcoord, fTMin;
- // Int_t fNtof, fNtofPairs;
- // -------------------------------------------------
-
- fKeepPos = new Bool_t [MAXHODHITS];
- fKeepNeg = new Bool_t [MAXHODHITS];
- fGoodTDCPos = new Bool_t [MAXHODHITS];
- fGoodTDCNeg = new Bool_t [MAXHODHITS];
-
- fHitPaddle = new Int_t [MAXHODHITS];
- fNScinHit = new Int_t [MAXHODHITS];
- fNPmtHit = new Int_t [MAXHODHITS];
- fTimeHist = new Int_t [MAXHODHITS];
- fTimeAtFP = new Double_t [MAXHODHITS];
-
- fScinSigma = new Double_t [MAXHODHITS];
- fGoodScinTime = new Double_t [MAXHODHITS];
- fScinTime = new Double_t [MAXHODHITS];
- fTime = new Double_t [MAXHODHITS];
- adcPh = new Double_t [MAXHODHITS];
-
- fPath = new Double_t [MAXHODHITS];
- fTimePos = new Double_t [MAXHODHITS];
- fTimeNeg = new Double_t [MAXHODHITS];
- fScinTimefp = new Double_t [MAXHODHITS];
-
- fTimeHist = new Int_t [200];
-
// -------------------------------------------------
Double_t hpartmass=0.00051099; // Fix it
for ( Int_t m = 0; m < MAXHODHITS; m++ ){
+
+ for(k = 0; k < MAXHODHITS; k++){
+ fdEdX[m][k] = 0.;
+ }
+ fGoodRawPad[m] = 0;
fScinSigma[m] = 0.;
fHitPaddle[m] = 0.;
fScinTime[m] = 0.;
@@ -764,24 +873,25 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
}
-
if (tracks.GetLast()+1 > 0 ) {
// **MAIN LOOP: Loop over all tracks and get corrected time, tof, beta...
- for ( Int_t itrack = 0; itrack < Ntracks; itrack++ ) { // Line 133
+ for ( itrack = 0; itrack < Ntracks; itrack++ ) { // Line 133
THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
if (!theTrack) return -1;
- fGoodPlaneTime = new Bool_t[4];
- for ( ip = 0; ip < fNPlanes; ip++ ){ fGoodPlaneTime[ip] = kFALSE; }
-
+ for ( ip = 0; ip < fNPlanes; ip++ ){
+ fGoodPlaneTime[ip] = kFALSE;
+ fNScinHits[ip] = 0;
+ }
+
fNfpTime = 0;
fBetaChisq[itrack] = -3;
fTimeAtFP[itrack] = 0.;
fSumfpTime = 0.; // Line 138
fNScinHit[itrack] = 0; // Line 140
- fP = theTrack->GetP(); // Line 142 Fix it
+ fP = theTrack->GetP(); // Line 142
fBetaP = fP/( TMath::Sqrt( fP * fP + hpartmass * hpartmass) );
//! Calculate all corrected hit times and histogram
@@ -799,12 +909,20 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
for (Int_t j=0; j<200; j++) { fTimeHist[j]=0; } // Line 176
- fNPlaneTime = new Int_t [4];
- fSumPlaneTime = new Double_t [4];
+ // fNPlaneTime = new Int_t [4];
+ // fSumPlaneTime = new Double_t [4];
for ( ip = 0; ip < fNPlanes; ip++ ){
fNPlaneTime[ip] = 0;
fSumPlaneTime[ip] = 0.;
+ // if ( ip == 2 ){
+ // fHodoCenter3 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ // fScin2XSpacing = fPlanes[ip]->GetSpacing();
+ // }
+ // if ( ip == 3 ){
+ // fHodoCenter4 = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ // fScin2YSpacing = fPlanes[ip]->GetSpacing();
+ // }
}
// Loop over scintillator planes.
// In ENGINE, its loop over good scintillator hits.
@@ -812,10 +930,10 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
fGoodTimeIndex = -1;
for( ip = 0; ip < fNPlanes; ip++ ) {
- fNScinHits = fPlanes[ip]->GetNScinHits();
+ fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
// first loop over hits with in a single plane
- for ( ihit = 0; ihit < fNScinHits; ihit++ ){
+ for ( ihit = 0; ihit < fNScinHits[ip]; ihit++ ){
fTimePos[ihit] = -99.;
fTimeNeg[ihit] = -99.;
@@ -917,7 +1035,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
if ( fJMax >= 0 ){ // Line 248. Here I followed the code of THcSCintilaltorPlane::PulseHeightCorrection
fTMin = 0.5 * fJMax;
- for( ihit = 0; ihit < fNScinHits; ihit++) { // Loop over sinc. hits. in plane
+ for( ihit = 0; ihit < fNScinHits[ip]; ihit++) { // Loop over sinc. hits. in plane
if ( ( fTimePos[ihit] > fTMin ) && ( fTimePos[ihit] < ( fTMin + fTofTolerance ) ) ) {
fKeepPos[ihit]=kTRUE;
}
@@ -927,8 +1045,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
}
} // fJMax > 0 condition
- fScinPosTime = new Double_t [MAXHODHITS];
- fScinNegTime = new Double_t [MAXHODHITS];
+ // fScinPosTime = new Double_t [MAXHODHITS];
+ // fScinNegTime = new Double_t [MAXHODHITS];
for ( k = 0; k < MAXHODHITS; k++ ){
fScinPosTime[k] = 0;
@@ -939,12 +1057,15 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// ---------------------- Scond loop over scint. hits in a plane ------------------------------
//---------------------------------------------------------------------------------------------
- for ( ihit = 0; ihit < fNScinHits; ihit++ ){
+ for ( ihit = 0; ihit < fNScinHits[ip]; ihit++ ){
fGoodTimeIndex ++;
+ fRawIndex ++;
fPaddle = ((THcSignalHit*)scinPosTDC->At(ihit))->GetPaddleNumber()-1;
- fHitPaddle[fGoodTimeIndex] = fPaddle;
+ fHitPaddle[fGoodTimeIndex] = fPaddle;
+
+ fGoodRawPad[fRawIndex] = fPaddle;
fXcoord = theTrack->GetX() + theTrack->GetTheta() *
( fPlanes[ip]->GetZpos() + ( fPaddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 277
@@ -1069,7 +1190,6 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// scinHit[itrack][fNScinHit[itrack]] = ihit;
// scinfFPTime[itrack][fNScinHit[itrack]] = fScinTimefp[ihit];
-
// ---------------------------------------------------------------------------
// Date: July 8 2014
// This counts the pmt hits. Right now we don't need it so it is commentd off
@@ -1083,29 +1203,34 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// ---------------------------------------------------------------------------
+ // fdEdX[itrack][ihit] = 5.0;
+
+
// --------------------------------------------------------------------------------------------
// Date: July 8 201 May be we need this, not sure.
//
- // if ( fGoodTDCPos[itrack][ihit] ){
- // if ( fGoodTDCNeg[itrack][ihit] ){
- // dedX[itrack][fNScinHit[itrack]] =
- // TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() *
- // ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ) );
- // }
- // else{
- // dedX[itrack][fNScinHit[itrack]] =
- // TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() );
- // }
- // }
- // else{
- // if ( fGoodTDCNeg[itrack][ihit] ){
- // dedX[itrack][fNScinHit[itrack]] =
- // TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(ihit))->GetData() );
- // }
- // else{
- // dedX[itrack][fNScinHit[itrack]] = 0.;
- // }
- // }
+
+ if ( fGoodTDCPos[fGoodTimeIndex] ){
+ if ( fGoodTDCNeg[fGoodTimeIndex] ){
+
+ fdEdX[itrack][fNScinHit[itrack]-1] =
+ TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() *
+ ((THcSignalHit*)scinNegADC->At(ihit))->GetData() ) );
+ }
+ else{
+ fdEdX[itrack][fNScinHit[itrack]-1] =
+ TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() );
+ }
+ }
+ else{
+ if ( fGoodTDCNeg[fGoodTimeIndex] ){
+ fdEdX[itrack][fNScinHit[itrack]-1] =
+ TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(ihit))->GetData() );
+ }
+ else{
+ fdEdX[itrack][fNScinHit[itrack]-1] = 0.;
+ }
+ }
// --------------------------------------------------------------------------------------------
@@ -1144,8 +1269,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
if (!fPlanes[ip])
return -1;
- fNScinHits = fPlanes[ip]->GetNScinHits();
- for (ihit = 0; ihit < fNScinHits; ihit++ ){
+ fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
+ for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){
fGoodTimeIndex ++;
if ( fGoodScinTime[fGoodTimeIndex] ) {
@@ -1178,8 +1303,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
if (!fPlanes[ip])
return -1;
- fNScinHits = fPlanes[ip]->GetNScinHits();
- for (ihit = 0; ihit < fNScinHits; ihit++ ){ // Loop over hits of a plane
+ fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
+ for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ // Loop over hits of a plane
fGoodTimeIndex ++;
if ( fGoodScinTime[fGoodTimeIndex] ){
@@ -1246,19 +1371,20 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
}
}
- } // Main loop over tracks ends here.
- } // If condition for at least one track
- ApplyCorrections();
-
- return 0;
-}
+ // fBetaChisq[itrack]
+ // fFPTime[ind]
+
+ theTrack->SetDedx(fdEdX[itrack][0]);
+ theTrack->SetBeta(fBeta[itrack]);
+ theTrack->SetBetaChi2( fBetaChisq[itrack] );
+
+ } // Main loop over tracks ends here.
+
+ } // If condition for at least one track
-//_____________________________________________________________________________
-Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
-{
-
return 0;
+
}
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index aba78e9c642d193539fdcfe77065d91e17734103..5ed3f59230ebe7e595a696c17325044cac36910d 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -14,6 +14,21 @@
#include "THcHitList.h"
#include "THcHodoscopeHit.h"
#include "THcScintillatorPlane.h"
+#include "THcShower.h"
+
+#include "THaTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcRawDCHit.h"
+#include "THcSpacePoint.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
+#include "TMath.h"
+
+#include "THaSubDetector.h"
+#include "TClonesArray.h"
+#include <iostream>
+#include <fstream>
+
class THaScCalib;
@@ -54,13 +69,22 @@ public:
Double_t GetStartTimeSlop() const {return fStartTimeSlop;}
Double_t GetBetaNotrk() const {return fBetaNotrk;}
+ Int_t GetGoodRawPad(Int_t iii){return fGoodRawPad[iii];}
+ Double_t GetNScinHits(Int_t iii){return fNScinHits[iii];}
+ Double_t GetHodoCenter3() { return fHodoCenter3;}
+ Double_t GetHodoCenter4() { return fHodoCenter4;}
+ Double_t GetScin2XSpacing() { return fScin2XSpacing;}
+ Double_t GetScin2YSpacing() { return fScin2YSpacing;}
+
// Double_t GetBeta() const {return fBeta[];}
Double_t GetBeta(Int_t iii) const {return fBeta[iii];} // Ahmed
+ // Int_t GetEvent(){ return fCheckEvent;}
Double_t GetHodoPosSigma(Int_t iii) const {return fHodoPosSigma[iii];}
Double_t GetHodoNegSigma(Int_t iii) const {return fHodoNegSigma[iii];}
+
const TClonesArray* GetTrackHits() const { return fTrackProj; }
friend class THaScCalib;
@@ -79,7 +103,7 @@ protected:
Double_t fStartTime;
Int_t fNfptimes;
- Double_t fdEdX;
+ // Double_t fdEdX;
Double_t fBetaNotrk;
// Double_t fBeta;
@@ -121,19 +145,43 @@ protected:
//-------------------------- Ahmed -----------------------------
- Int_t MAXHODHITS;
-
- Double_t* fTestArr; // [MAXHODHITS] Array
+ THcShower* fShower;
+
+ Int_t fCheckEvent;
+
+ Int_t fGoodTrack;
+ Int_t MAXHODHITS;
+ // Int_t fSelUsingScin;
+ Int_t fSelNDegreesMin;
+ Double_t fSeldEdX1Min;
+ Double_t fSeldEdX1Max;
+ Double_t fSelBetaMin;
+ Double_t fSelBetaMax;
+ Double_t fSelEtMin;
+ Double_t fSelEtMax;
+ Double_t fScin2XZpos;
+ Double_t fScin2XdZpos;
+ Double_t fScin2YZpos;
+ Double_t fScin2YdZpos;
+
+ Double_t fChi2Min;
+ Double_t fHodoCenter4, fHodoCenter3;
+ Double_t fScin2YSpacing, fScin2XSpacing;
+
+ Double_t** fdEdX; // [MAXHODHITS] Array
+ Double_t** fScinHit; // [fNPlanes] Array
+ Int_t* fGoodRawPad;
Double_t* fBeta; // [MAXHODHITS] Array
Double_t* fBetaChisq; // [MAXHODHITS] Array
Double_t* fFPTime; // [fNPlanes] Array
+
Double_t* fScinSigma;
Double_t* fGoodScinTime;
Double_t* fScinTime;
Double_t* fTime;
- Double_t* adcPh;
+ Double_t* adcPh; // Correct it
Double_t* fTimeAtFP;
Double_t* fPath;
Double_t* fTimePos;
@@ -145,6 +193,7 @@ protected:
Int_t* fHitPaddle;
Int_t* fNScinHit;
+ Int_t* fNScinHits;
Int_t* fNPmtHit;
Int_t* fTimeHist;
Int_t* fNPlaneTime;
@@ -163,6 +212,10 @@ protected:
TClonesArray* scinPosTDC;
TClonesArray* scinNegTDC;
+ //test array
+ Double_t test_arr[53];
+ Double_t test_arr1[2];
+
//----------------------------------------------------------------
// Useful derived quantities
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index e945e5105a40c10248255df7f7163b072740ff23..e6a8073c6efd392837e2fbb4cbcf5f7add9f999c 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -45,6 +45,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fNScinHits = 0;
//
fMaxHits=53;
+
fpTimes = new Double_t [fMaxHits];
fScinTime = new Double_t [fMaxHits];
fScinSigma = new Double_t [fMaxHits];
@@ -72,6 +73,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fNScinHits = 0;
//
fMaxHits=53;
+
fpTimes = new Double_t [fMaxHits];
fScinTime = new Double_t [fMaxHits];
fScinSigma = new Double_t [fMaxHits];
@@ -94,6 +96,7 @@ THcScintillatorPlane::~THcScintillatorPlane()
delete fScinTime;
delete fScinSigma;
delete fScinZpos;
+
}
//______________________________________________________________________________
@@ -303,20 +306,29 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
mintdc=((THcHodoscope *)GetParent())->GetTdcMin();
maxtdc=((THcHodoscope *)GetParent())->GetTdcMax();
Int_t ihit = nexthit;
+
+ // cout << "THcScintillatorPlane: raw htis = " << nrawhits << endl;
+
while(ihit < nrawhits) {
THcHodoscopeHit* hit = (THcHodoscopeHit *) rawhits->At(ihit);
if(hit->fPlane > fPlaneNum) {
break;
}
Int_t padnum=hit->fCounter;
+
Int_t index=padnum-1;
- if (hit->fTDC_pos > 0) ((THcSignalHit*) frPosTDCHits->ConstructedAt(nrPosTDCHits++))->Set(padnum, hit->fTDC_pos);
- if (hit->fTDC_neg > 0) ((THcSignalHit*) frNegTDCHits->ConstructedAt(nrNegTDCHits++))->Set(padnum, hit->fTDC_neg);
- if ((hit->fADC_pos-fPosPed[index]) >= 50) ((THcSignalHit*) frPosADCHits->ConstructedAt(nrPosADCHits++))->Set(padnum, hit->fADC_pos-fPosPed[index]);
- if ((hit->fADC_neg-fNegPed[index]) >= 50) ((THcSignalHit*) frNegADCHits->ConstructedAt(nrNegADCHits++))->Set(padnum, hit->fADC_neg-fNegPed[index]);
+ if (hit->fTDC_pos > 0)
+ ((THcSignalHit*) frPosTDCHits->ConstructedAt(nrPosTDCHits++))->Set(padnum, hit->fTDC_pos);
+ if (hit->fTDC_neg > 0)
+ ((THcSignalHit*) frNegTDCHits->ConstructedAt(nrNegTDCHits++))->Set(padnum, hit->fTDC_neg);
+ if ((hit->fADC_pos-fPosPed[index]) >= 50)
+ ((THcSignalHit*) frPosADCHits->ConstructedAt(nrPosADCHits++))->Set(padnum, hit->fADC_pos-fPosPed[index]);
+ if ((hit->fADC_neg-fNegPed[index]) >= 50)
+ ((THcSignalHit*) frNegADCHits->ConstructedAt(nrNegADCHits++))->Set(padnum, hit->fADC_neg-fNegPed[index]);
// check TDC values
if (((hit->fTDC_pos >= mintdc) && (hit->fTDC_pos <= maxtdc)) ||
((hit->fTDC_neg >= mintdc) && (hit->fTDC_neg <= maxtdc))) {
+
//TDC positive hit
THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
sighit->Set(padnum, hit->fTDC_pos);
@@ -332,13 +344,12 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fNScinHits++;
}
else {
- //cout <<"pos TDC "<<hit->fTDC_pos<<" "<<mintdc<<" "<<maxtdc<<endl;
- //cout <<"neg TDC "<<hit->fTDC_neg<<" "<<mintdc<<" "<<maxtdc<<endl;
- //cout <<"skipping BAD tdc event\n";
}
ihit++;
}
+ // cout << "THcScintillatorPlane: ihit = " << ihit << endl;
+
return(ihit);
}
//________________________________________________________________________________
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index 959038f6f3c926e468ee8e8f832d13e2fdc82076..fce7cd336ea326a95fee1114a023f66640700d33 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -63,12 +63,10 @@ class THcScintillatorPlane : public THaSubDetector {
TClonesArray* fParentHitList;
- TClonesArray* GetPosADC() { return fPosADCHits;}; // Ahmed
- TClonesArray* GetNegADC() { return fNegADCHits;}; // Ahmed
- TClonesArray* GetPosTDC() { return fPosTDCHits;}; // Ahmed
- TClonesArray* GetNegTDC() { return fNegTDCHits;}; // Ahmed
-
-
+ TClonesArray* GetPosADC() { return fPosADCHits;}; // Ahmed
+ TClonesArray* GetNegADC() { return fNegADCHits;}; // Ahmed
+ TClonesArray* GetPosTDC() { return fPosTDCHits;}; // Ahmed
+ TClonesArray* GetNegTDC() { return fNegTDCHits;}; // Ahmed
protected:
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 060d120c909fa87747d47ed9f4464ebd1350ba35..b2f328be0eeaad2a049cace869ccb51ffde3fd98 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -522,6 +522,8 @@ Int_t THcShower::Decode( const THaEvData& evdata )
// Get the Hall C style hitlist (fRawHitList) for this event
Int_t nhits = DecodeToHitList(evdata);
+ fEvent = evdata.GetEvNum();
+
if(gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
@@ -1014,6 +1016,11 @@ Int_t THcShower::FineProcess( TClonesArray& tracks )
Float_t energy = GetShEnergy(theTrack);
theTrack->SetEnergy(energy);
+ // if ( fEvent == 13252 )
+ // cout << "THcShower: track = " << itrk + 1
+ // << " energy = " << energy << endl;
+
+
if (fdbg_tracks_cal) {
cout << "THcShower::FineProcess, Track " << itrk << ": "
<< " X = " << theTrack->GetX()
diff --git a/src/THcShower.h b/src/THcShower.h
index 04d351dfd28b55b8969390a3188663f55b88dcdc..01de29c182d43ecb8eccc08bdd46ca98080530b6 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -118,6 +118,8 @@ public:
protected:
+ Int_t fEvent;
+
Int_t fAnalyzePedestals; // Flag for pedestal analysis.
Int_t* fShPosPedLimit; // [fNtotBlocks] ADC limits for pedestal calc.-s.