Newer
Older
/** \class THcHodoscope
\ingroup Detectors
Class for a generic hodoscope consisting of multiple
planes with multiple paddles with phototubes on both ends.
This differs from Hall A scintillator class in that it is the whole
hodoscope array, not just one plane.
*/
#include "THcHodoHit.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
#include "THaSubDetector.h"
#include "THcHodoscope.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
Gabriel Niculescu
committed
#include <fstream>
using namespace std;
//_____________________________________________________________________________
THcHodoscope::THcHodoscope( const char* name, const char* description,
THaApparatus* apparatus ) :
THaNonTrackingDetector(name,description,apparatus)
{
// Constructor
//fTrackProj = new TClonesArray( "THaTrackProj", 5 );
// Construct the planes
Stephen A. Wood
committed
fNPlanes = 0; // No planes until we make them
//_____________________________________________________________________________
THcHodoscope::THcHodoscope( ) :
THaNonTrackingDetector()
{
// Constructor
}
//_____________________________________________________________________________
void THcHodoscope::Setup(const char* name, const char* description)
{
// static const char* const here = "Setup()";
// static const char* const message =
// "Must construct %s detector with valid name! Object construction failed.";
Stephen A. Wood
committed
cout << "In THcHodoscope::Setup()" << endl;
// Base class constructor failed?
if( IsZombie()) return;
fDebug = 1; // Keep this at one while we're working on the code
char prefix[2];
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
DBRequest listextra[]={
{"hodo_num_planes", &fNPlanes, kInt},
{"hodo_plane_names",&planenamelist, kString},
//fNPlanes = 4; // Default if not defined
gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);
cout << "Plane Name List : " << planenamelist << endl;
vector<string> plane_names = vsplit(planenamelist);
// Plane names
if(plane_names.size() != (UInt_t) fNPlanes) {
cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
// Should quit. Is there an official way to quit?
}
fPlaneNames = new char* [fNPlanes];
for(Int_t i=0;i<fNPlanes;i++) {
fPlaneNames[i] = new char[plane_names[i].length()+1];
strcpy(fPlaneNames[i], plane_names[i].c_str());
}
// Probably shouldn't assume that description is defined
char* desc = new char[strlen(description)+100];
fPlanes = new THcScintillatorPlane* [fNPlanes];
for(Int_t i=0;i < fNPlanes;i++) {
strcpy(desc, description);
strcat(desc, " Plane ");
strcat(desc, fPlaneNames[i]);
fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i+1, this); // Number planes starting from zero!!
cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
// --------------- To get energy from THcShower ----------------------
const char* shower_detector_name = "cal";
// THaApparatus* app;
THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
THaDetector* det = app->GetDetector( shower_detector_name );
if( dynamic_cast<THcShower*>(det) ) {
fShower = dynamic_cast<THcShower*>(det);
}
else if( !dynamic_cast<THcShower*>(det) ) {
cout << "Warining: calorimeter analysis module "
<< shower_detector_name << " not loaded for spectrometer "
<< prefix << endl;
fShower = NULL;
}
// --------------- To get energy from THcShower ----------------------
// --------------- To get NPEs from THcCherenkov -------------------
const char* chern_detector_name = "cher";
THaDetector* detc = app->GetDetector( chern_detector_name );
if( dynamic_cast<THcCherenkov*>(detc) ) {
fChern = dynamic_cast<THcCherenkov*>(detc);
}
else if( !dynamic_cast<THcCherenkov*>(detc) ) {
cout << "Warining: Cherenkov detector analysis module "
<< chern_detector_name << " not loaded for spectrometer "
<< prefix << endl;
fChern = NULL;
}
// --------------- To get NPEs from THcCherenkov -------------------
fScinShould = 0;
fScinDid = 0;
gHcParms->Define(Form("%shodo_did",prefix),"Total hodo tracks",fScinDid);
gHcParms->Define(Form("%shodo_should",prefix),"Total hodo triggers",fScinShould);
// Save the nominal particle mass
fPartMass = app->GetParticleMass();
fBetaNominal = app->GetBetaAtPcentral();
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
{
cout << "In THcHodoscope::Init()" << endl;
Setup(GetName(), GetTitle());
// Should probably put this in ReadDatabase as we will know the
// maximum number of hits after setting up the detector map
// But it needs to happen before the sub detectors are initialized
// so that they can get the pointer to the hitlist.
InitHitList(fDetMap, "THcRawHodoHit", 100);
EStatus status;
// This triggers call of ReadDatabase and DefineVariables
if( (status = THaNonTrackingDetector::Init( date )) )
return fStatus=status;
for(Int_t ip=0;ip<fNPlanes;ip++) {
if((status = fPlanes[ip]->Init( date ))) {
return fStatus=status;
}
}
// Replace with what we need for Hall C
// const DataDest tmp[NDEST] = {
// { &fRTNhit, &fRANhit, fRT, fRT_c, fRA, fRA_p, fRA_c, fROff, fRPed, fRGain },
// { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
// };
// memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
char EngineDID[]="xSCIN";
EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
static const char* const here = "Init()";
Error( Here(here), "Error filling detectormap for %s.",
EngineDID);
return kInitError;
}
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;
// fScinHit = new Double_t*[fNPlanes];
// for ( m = 0; m < fNPlanes; m++ ){
// fScinHit[m] = new Double_t[fNPaddle[0]];
// }
return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcHodoscope::ReadDatabase( const TDatime& date )
{
// Read this detector's parameters from the database file 'fi'.
// This function is called by THaDetectorBase::Init() once at the
// beginning of the analysis.
// 'date' contains the date/time of the run being analyzed.
// static const char* const here = "ReadDatabase()";
char prefix[2];
char parname[100];
// Read data from database
// Pull values from the THcParmList instead of reading a database
// file like Hall A does.
// Will need to determine which spectrometer in order to construct
// the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
strcpy(parname,prefix);
strcat(parname,"scin_");
// Int_t plen=strlen(parname);
cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;
fFPTime = new Double_t [fNPlanes];
fPlaneCenter = new Double_t[fNPlanes];
fPlaneSpacing = new Double_t[fNPlanes];
prefix[0]=tolower(GetApparatus()->GetName()[0]);
//
prefix[1]='\0';
for(Int_t i=0;i<fNPlanes;i++) {
DBRequest list[]={
{Form("scin_%s_nr",fPlaneNames[i]), &fNPaddle[i], kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
// GN added
// reading variables from *hodo.param
Gabriel Niculescu
committed
fMaxScinPerPlane=fNPaddle[0];
for (Int_t i=1;i<fNPlanes;i++) {
Gabriel Niculescu
committed
fMaxScinPerPlane=(fMaxScinPerPlane > fNPaddle[i])? fMaxScinPerPlane : fNPaddle[i];
}
// need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
fMaxHodoScin=fMaxScinPerPlane*fNPlanes;
if (fDebug>=1) cout <<"fMaxScinPerPlane = "<<fMaxScinPerPlane<<" fMaxHodoScin = "<<fMaxHodoScin<<endl;
fHodoVelLight=new Double_t [fMaxHodoScin];
fHodoPosSigma=new Double_t [fMaxHodoScin];
fHodoNegSigma=new Double_t [fMaxHodoScin];
fHodoPosMinPh=new Double_t [fMaxHodoScin];
fHodoNegMinPh=new Double_t [fMaxHodoScin];
fHodoPosPhcCoeff=new Double_t [fMaxHodoScin];
fHodoNegPhcCoeff=new Double_t [fMaxHodoScin];
fHodoPosTimeOffset=new Double_t [fMaxHodoScin];
fHodoNegTimeOffset=new Double_t [fMaxHodoScin];
fHodoPosPedLimit=new Int_t [fMaxHodoScin];
fHodoNegPedLimit=new Int_t [fMaxHodoScin];
fHodoPosInvAdcOffset=new Double_t [fMaxHodoScin];
fHodoNegInvAdcOffset=new Double_t [fMaxHodoScin];
fHodoPosInvAdcLinear=new Double_t [fMaxHodoScin];
fHodoNegInvAdcLinear=new Double_t [fMaxHodoScin];
fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin];
fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin];
fNHodoscopes = 2;
fxLoScin = new Int_t [fNHodoscopes];
fxHiScin = new Int_t [fNHodoscopes];
fyLoScin = new Int_t [fNHodoscopes];
fyHiScin = new Int_t [fNHodoscopes];
fHodoSlop = new Double_t [fNPlanes];
DBRequest list[]={
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
{"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},
{"xloscin", &fxLoScin[0], kInt, (UInt_t) fNHodoscopes},
{"xhiscin", &fxHiScin[0], kInt, (UInt_t) fNHodoscopes},
{"yloscin", &fyLoScin[0], kInt, (UInt_t) fNHodoscopes},
{"yhiscin", &fyHiScin[0], kInt, (UInt_t) fNHodoscopes},
{"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
{"cer_npe", &fNCerNPE, kDouble, 0, 1},
{"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
{"hodo_slop", fHodoSlop, kDouble, (UInt_t) fNPlanes},
{"debugprintscinraw", &fdebugprintscinraw, kInt, 0,1},
{0}
};
fdebugprintscinraw=0;
fTofUsingInvAdc = 0; // Default if not defined
Mark Jones
committed
fTofTolerance = 3.0; // Default if not defined
fNCerNPE = 2.0;
fNormETot = 0.7;
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
cout << " x1 lo = " << fxLoScin[0]
<< " x2 lo = " << fxLoScin[1]
<< " x1 hi = " << fxHiScin[0]
<< " x2 hi = " << fxHiScin[1]
<< endl;
cout << " y1 lo = " << fyLoScin[0]
<< " y2 lo = " << fyLoScin[1]
<< " y1 hi = " << fyHiScin[0]
<< " y2 hi = " << fyHiScin[1]
<< endl;
cout << "Hdososcope planes hits for trigger = " << fTrackEffTestNScinPlanes
<< " normalized energy min = " << fNormETot
<< " number of photo electrons = " << fNCerNPE
<< endl;
if (fTofUsingInvAdc) {
DBRequest list2[]={
{"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin},
{"hodo_neg_invadc_offset",&fHodoNegInvAdcOffset[0],kDouble,fMaxHodoScin},
{"hodo_pos_invadc_linear",&fHodoPosInvAdcLinear[0],kDouble,fMaxHodoScin},
{"hodo_neg_invadc_linear",&fHodoNegInvAdcLinear[0],kDouble,fMaxHodoScin},
{"hodo_pos_invadc_adc",&fHodoPosInvAdcAdc[0],kDouble,fMaxHodoScin},
{"hodo_neg_invadc_adc",&fHodoNegInvAdcAdc[0],kDouble,fMaxHodoScin},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list2,prefix);
};
Gabriel Niculescu
committed
if (fDebug >=1) {
cout <<"******* Testing Hodoscope Parameter Reading ***\n";
cout<<"StarTimeCenter = "<<fStartTimeCenter<<endl;
cout<<"StartTimeSlop = "<<fStartTimeSlop<<endl;
cout <<"ScintTdcToTime = "<<fScinTdcToTime<<endl;
cout <<"TdcMin = "<<fScinTdcMin<<" TdcMax = "<<fScinTdcMax<<endl;
cout <<"TofTolerance = "<<fTofTolerance<<endl;
cout <<"*** VelLight ***\n";
Gabriel Niculescu
committed
cout<<"Plane "<<i1<<endl;
Gabriel Niculescu
committed
cout<<fHodoVelLight[GetScinIndex(i1,i2)]<<" ";
Gabriel Niculescu
committed
cout <<endl;
Gabriel Niculescu
committed
cout <<endl<<endl;
// check fHodoPosPhcCoeff
/*
cout <<"fHodoPosPhcCoeff = ";
for (int i1=0;i1<fMaxHodoScin;i1++) {
cout<<this->GetHodoPosPhcCoeff(i1)<<" ";
}
cout<<endl;
*/
}
//
if ((fTofTolerance > 0.5) && (fTofTolerance < 10000.)) {
cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
}
else {
fTofTolerance= 3.0;
cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
fIsInit = true;
return kOK;
}
//_____________________________________________________________________________
Int_t THcHodoscope::DefineVariables( EMode mode )
{
// Initialize global variables and lookup table for decoder
cout << "THcHodoscope::DefineVariables called " << GetName() << endl;
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
// Register variables in global list
// Move these into THcHallCSpectrometer using track fTracks
{"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
{"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
{"fpHitsTime", "Time at focal plane from all hits", "fFPTime"},
{"starttime", "Hodoscope Start Time", "fStartTime"},
{"goodstarttime", "Hodoscope Good Start Time", "fGoodStartTime"},
{"goodscinhit", "Hit in fid area", "fGoodScinHits"},
// {"goodscinhitx", "Hit in fid x range", "fGoodScinHitsX"},
{"scinshould", "Total scin Hits in fid area", "fScinShould"},
{"scindid", "Total scin Hits in fid area with a track", "fScinDid"},
return DefineVarsFromList( vars, mode );
// return kOK;
}
//_____________________________________________________________________________
THcHodoscope::~THcHodoscope()
{
// Destructor. Remove variables from global list.
delete [] fPlaneCenter;
delete [] fPlaneSpacing;
if( fIsSetup )
RemoveVariables();
if( fIsInit )
DeleteArrays();
if (fTrackProj) {
fTrackProj->Clear();
delete fTrackProj; fTrackProj = 0;
}
}
//_____________________________________________________________________________
void THcHodoscope::DeleteArrays()
{
// Delete member arrays. Used by destructor.
// Int_t k;
// for( k = 0; k < fNPlanes; k++){
// delete [] fScinHit[k];
// }
// delete [] fScinHit;
delete [] fxLoScin; fxLoScin = NULL;
delete [] fxHiScin; fxHiScin = NULL;
delete [] fHodoSlop; fHodoSlop = NULL;
delete [] fNPaddle; fNPaddle = NULL;
delete [] fHodoVelLight; fHodoVelLight = NULL;
delete [] fHodoPosSigma; fHodoPosSigma = NULL;
delete [] fHodoNegSigma; fHodoNegSigma = NULL;
delete [] fHodoPosMinPh; fHodoPosMinPh = NULL;
delete [] fHodoNegMinPh; fHodoNegMinPh = NULL;
delete [] fHodoPosPhcCoeff; fHodoPosPhcCoeff = NULL;
delete [] fHodoNegPhcCoeff; fHodoNegPhcCoeff = NULL;
delete [] fHodoPosTimeOffset; fHodoPosTimeOffset = NULL;
delete [] fHodoNegTimeOffset; fHodoNegTimeOffset = NULL;
delete [] fHodoPosPedLimit; fHodoPosPedLimit = NULL;
delete [] fHodoNegPedLimit; fHodoNegPedLimit = NULL;
delete [] fHodoPosInvAdcOffset; fHodoPosInvAdcOffset = NULL;
delete [] fHodoNegInvAdcOffset; fHodoNegInvAdcOffset = NULL;
delete [] fHodoPosInvAdcLinear; fHodoPosInvAdcLinear = NULL;
delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL;
delete [] fHodoPosInvAdcAdc; fHodoPosInvAdcAdc = NULL;
delete [] fGoodPlaneTime; fGoodPlaneTime = NULL;
delete [] fNPlaneTime; fNPlaneTime = NULL;
delete [] fSumPlaneTime; fSumPlaneTime = NULL;
delete [] fNScinHits; fNScinHits = NULL;
}
//_____________________________________________________________________________
inline
void THcHodoscope::ClearEvent()
/*! \brief Clears variables
*
* Called by THcHodoscope::Decode
*
*/
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip]->Clear();
fPlaneCenter[ip]=0.;
fPlaneSpacing[ip]=0.;
fdEdX.clear();
fScinHitPaddle.clear();
fNScinHit.clear();
fNClust.clear();
fThreeScin.clear();
fGoodScinHitsX.clear();
}
//_____________________________________________________________________________
Int_t THcHodoscope::Decode( const THaEvData& evdata )
{
/*! \brief Decodes raw data and processes raw data into hits for each instance of THcScintillatorPlane
*
* - Calls THcHodoscope::ClearEvent
* - Reads raw data using THcHitList::DecodeToHitList
* - If one wants to subtract pedestals (assumed to be a set of data at beginning of run)
* + Must define "Pedestal_event" cut in the cuts definition file
* + For each "Pedestal_event" calls THcScintillatorPlane::AccumulatePedestals and returns
* + After First event which is not a "Pedestal_event" calls THcScintillatorPlane::CalculatePedestals
* - For each scintillator plane THcScintillatorPlane::ProcessHits
* - Calls THcHodoscope::EstimateFocalPlaneTime
*
*
*/
ClearEvent();
// Get the Hall C style hitlist (fRawHitList) for this event
Stephen A. Wood
committed
Int_t nhits = DecodeToHitList(evdata);
//
// GN: print event number so we can cross-check with engine
// if (evdata.GetEvNum()>1000)
fCheckEvent = evdata.GetEvNum();
if(gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
}
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return(0);
}
if(fAnalyzePedestals) {
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip]->CalculatePedestals();
}
fAnalyzePedestals = 0; // Don't analyze pedestals next event
}
// Let each plane get its hits
Int_t nexthit = 0;
fNfptimes=0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlaneCenter[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
fPlaneSpacing[ip] = fPlanes[ip]->GetSpacing();
// nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
// GN: select only events that have reasonable TDC values to start with
// as per the Engine h_strip_scin.f
nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
EstimateFocalPlaneTime();
if (fdebugprintscinraw == 1) {
Mark Jones
committed
cout << " Event number = " << evdata.GetEvNum()<<endl;
for(UInt_t ihit = 0; ihit < fNRawHits ; ihit++) {
THcRawHodoHit* hit = (THcRawHodoHit *) fRawHitList->At(ihit);
cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
<< hit->fADC_pos << " " << hit->fADC_neg << " " << hit->fTDC_pos
<< " " << hit->fTDC_neg << endl;
}
cout << endl;
Stephen A. Wood
committed
return nhits;
}
//_____________________________________________________________________________
void THcHodoscope::EstimateFocalPlaneTime( void )
{
/*! \brief Calculates the Drift Chamber start time and fBetaNoTrk (velocity determined without track info)
*
* - Called by THcHodoscope::Decode
* - selects good scintillator paddle hits
* + loops through hits in each scintillator plane and fills histogram array, "timehist", with corrected times for positive
* and negative ends of each paddle
* + Determines the peak of "timehist"
*
*/
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
Int_t timehist[200];
for (Int_t i=0;i<200;i++) {
timehist[i] = 0;
}
Int_t ihit=0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for(Int_t i=0;i<nphits;i++) {
Double_t postime=((THcHodoHit*) hodoHits->At(i))->GetPosTOFCorrectedTime();
Double_t negtime=((THcHodoHit*) hodoHits->At(i))->GetNegTOFCorrectedTime();
for (Int_t k=0;k<200;k++) {
Double_t tmin=0.5*(k+1);
if ((postime> tmin) && (postime < tmin+fTofTolerance)) {
timehist[k]++;
}
if ((negtime> tmin) && (negtime < tmin+fTofTolerance)) {
timehist[k]++;
}
}
}
}
// Find the bin with most hits
Int_t binmax=0;
Int_t maxhit=0;
for(Int_t i=0;i<200;i++) {
if(timehist[i]>maxhit) {
binmax = i+1;
maxhit = timehist[i];
}
}
ihit = 0;
Double_t fpTimeSum = 0.0;
fNfptimes=0;
Mark Jones
committed
Int_t Ngood_hits_plane=0;
Double_t Plane_fptime_sum=0.0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
fNoTrkPlaneInfo.push_back(NoTrkPlaneInfo());
fNoTrkPlaneInfo[ip].goodplanetime = kFALSE;
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Mark Jones
committed
Ngood_hits_plane=0;
Plane_fptime_sum=0.0;
for(Int_t i=0;i<nphits;i++) {
Mark Jones
committed
fNoTrkHitInfo[i].goodtwotimes = kFALSE;
fNoTrkHitInfo[i].goodscintime = kFALSE;
Double_t tmin = 0.5*binmax;
Double_t postime=((THcHodoHit*) hodoHits->At(i))->GetPosTOFCorrectedTime();
Double_t negtime=((THcHodoHit*) hodoHits->At(i))->GetNegTOFCorrectedTime();
if ((postime>tmin) && (postime<tmin+fTofTolerance) &&
(negtime>tmin) && (negtime<tmin+fTofTolerance)) {
Mark Jones
committed
fNoTrkHitInfo[i].goodtwotimes = kTRUE;
fNoTrkHitInfo[i].goodscintime = kTRUE;
// Both tubes fired
Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
Double_t fptime = ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime()
- (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
/ (29.979 * fBetaNominal);
if(TMath::Abs(fptime-fStartTimeCenter)<=fStartTimeSlop) {
Mark Jones
committed
Ngood_hits_plane++;
Plane_fptime_sum+=fptime;
fpTimeSum += fptime;
fNfptimes++;
}
}
Mark Jones
committed
fPlanes[ip]->SetFpTime(Plane_fptime_sum/float(Ngood_hits_plane));
fPlanes[ip]->SetNGoodHits(Ngood_hits_plane);
}
ihit++;
}
if(fNfptimes>0) {
fStartTime = fpTimeSum/fNfptimes;
fGoodStartTime=kTRUE;
} else {
fStartTime = fStartTimeCenter;
fGoodStartTime=kFALSE;
}
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
if ( ( fNoTrkPlaneInfo[0].goodplanetime || fNoTrkPlaneInfo[1].goodplanetime ) &&
( fNoTrkPlaneInfo[2].goodplanetime || fNoTrkPlaneInfo[3].goodplanetime ) ){
Double_t sumW = 0.;
Double_t sumT = 0.;
Double_t sumZ = 0.;
Double_t sumZZ = 0.;
Double_t sumTZ = 0.;
Int_t ihhit = 0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for(Int_t i=0;i<nphits;i++) {
Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
if ( fNoTrkHitInfo[ihhit].goodscintime ) {
Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
Double_t scinWeight = 1 / TMath::Power(sigma,2);
Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
// cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " << sigma << endl;
sumW += scinWeight;
sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
sumZ += scinWeight * zPosition;
sumZZ += scinWeight * ( zPosition * zPosition );
sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
} // condition of good scin time
ihhit ++;
} // loop over hits of plane
} // loop over planes
Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
fBetaNoTrk = tmp / tmpDenom;
fBetaNoTrkChiSq = 0.;
ihhit = 0;
for (Int_t ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for(Int_t i=0;i<nphits;i++) {
Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
if ( fNoTrkHitInfo[ihhit].goodscintime ) {
Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
Double_t timeDif = ( ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0 );
Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
fBetaNoTrkChiSq += ( ( zPosition / fBetaNoTrk - timeDif ) *
( zPosition / fBetaNoTrk - timeDif ) ) / ( sigma * sigma );
} // condition for good scin time
ihhit++;
} // loop over hits of a plane
} // loop over planes
Double_t pathNorm = 1.0;
fBetaNoTrk = fBetaNoTrk * pathNorm;
fBetaNoTrk = fBetaNoTrk / 29.979; // velocity / c
} // condition for fTmpDenom
else {
fBetaNoTrk = 0.;
fBetaNoTrkChiSq = -2.;
} // else condition for fTmpDenom
}
//_____________________________________________________________________________
Int_t THcHodoscope::ApplyCorrections( void )
{
return(0);
}
//_____________________________________________________________________________
Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle,
const ESide side)
{
return(0.0);
}
//_____________________________________________________________________________
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 timehist[200];
// -------------------------------------------------
fScinShould = 0; fScinDid = 0;
// **MAIN LOOP: Loop over all tracks and get corrected time, tof, beta...
Double_t* nPmtHit = new Double_t [ntracks];
Double_t* timeAtFP = new Double_t [ntracks];
for ( Int_t itrack = 0; itrack < ntracks; itrack++ ) { // Line 133
nPmtHit[itrack]=0;
timeAtFP[itrack]=0;
THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
if (!theTrack) return -1;
fNPlaneTime[ip] = 0;
fSumPlaneTime[ip] = 0.;
std::vector<Double_t> dedx_temp;
fdEdX.push_back(dedx_temp); // Create array of dedx per hit
std::vector<std::vector<GoodFlags> > goodflagstmp1;
fGoodFlags.push_back(goodflagstmp1);
// timeAtFP[itrack] = 0.;
Double_t sumFPTime = 0.; // Line 138
fNScinHit.push_back(0);
Double_t p = theTrack->GetP(); // Line 142
fBetaP = p/( TMath::Sqrt( p * p + fPartMass * fPartMass) );
//! Calculate all corrected hit times and histogram
//! This uses a copy of code below. Results are save in time_pos,neg
//! including the z-pos. correction assuming nominal value of betap
//! Code is currently hard-wired to look for a peak in the
//! range of 0 to 100 nsec, with a group of times that all
//! agree withing a time_tolerance of time_tolerance nsec. The normal
//! peak position appears to be around 35 nsec.
//! NOTE: if want to find farticles with beta different than
//! reference particle, need to make sure this is big enough
//! to accomodate difference in TOF for other particles
//! Default value in case user hasnt definedd something reasonable
// Line 162 to 171 is already done above in ReadDatabase
for (Int_t j=0; j<200; j++) { timehist[j]=0; } // Line 176
// Loop over scintillator planes.
// In ENGINE, its loop over good scintillator hits.
fTOFCalc.clear();
Int_t ihhit = 0; // Hit # overall
std::vector<GoodFlags> goodflagstmp2;
fGoodFlags[itrack].push_back(goodflagstmp2);
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// first loop over hits with in a single plane
fTOFPInfo.clear();
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
// iphit is hit # within a plane
fTOFPInfo.push_back(TOFPInfo());
// Can remove these as we will initialize in the constructor
fTOFPInfo[iphit].time_pos = -99.0;
fTOFPInfo[iphit].time_neg = -99.0;
fTOFPInfo[iphit].keep_pos = kFALSE;
fTOFPInfo[iphit].keep_neg = kFALSE;
fTOFPInfo[iphit].scin_pos_time = 0.0;
fTOFPInfo[iphit].scin_neg_time = 0.0;
Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
( paddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 183
Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
( paddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 184
if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 185
scinTrnsCoord = xHitCoord;
scinLongCoord = yHitCoord;
}
else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 188
scinTrnsCoord = yHitCoord;
scinLongCoord = xHitCoord;
}
else { return -1; } // Line 195
Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// Index to access the 2d arrays of paddle/scintillator properties
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
Double_t timep = ((THcHodoHit*)hodoHits->At(iphit))->GetPosCorrectedTime();
timep = timep - ( pathp / fHodoVelLight[fPIndex] ) - ( fPlanes[ip]->GetZpos() +
( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * fBetaP ) *
TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
fTOFPInfo[iphit].time_pos = timep;
for ( Int_t k = 0; k < 200; k++ ){ // Line 211
Double_t tmin = 0.5 * ( k + 1 ) ;
if ( ( timep > tmin ) && ( timep < ( tmin + fTofTolerance ) ) )
timehist[k] ++;
}
Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
Double_t timen = ((THcHodoHit*)hodoHits->At(iphit))->GetNegCorrectedTime();
timen = timen - ( pathn / fHodoVelLight[fPIndex] ) - ( fPlanes[ip]->GetZpos() +
( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * fBetaP ) *
TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
fTOFPInfo[iphit].time_neg = timen;
for ( Int_t k = 0; k < 200; k++ ){ // Line 230
Double_t tmin = 0.5 * ( k + 1 );
if ( ( timen > tmin ) && ( timen < ( tmin + fTofTolerance ) ) )
timehist[k] ++;
}
} // condition for cenetr on a paddle
} // First loop over hits in a plane <---------
//-----------------------------------------------------------------------------------------------
//------------- First large loop over scintillator hits in a plane ends here --------------------
//-----------------------------------------------------------------------------------------------
Int_t jmax = 0; // Line 240
Int_t maxhit = 0;
for ( Int_t k = 0; k < 200; k++ ){
if ( timehist[k] > maxhit ){
jmax = k+1;
maxhit = timehist[k];
Double_t tmin = 0.5 * jmax;
for(Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++) { // Loop over sinc. hits. in plane
if ( ( fTOFPInfo[iphit].time_pos > tmin ) && ( fTOFPInfo[iphit].time_pos < ( tmin + fTofTolerance ) ) ) {
fTOFPInfo[iphit].keep_pos=kTRUE;
}
if ( ( fTOFPInfo[iphit].time_neg > tmin ) && ( fTOFPInfo[iphit].time_neg < ( tmin + fTofTolerance ) ) ){
fTOFPInfo[iphit].keep_neg=kTRUE;
}
}
//---------------------------------------------------------------------------------------------
// ---------------------- Scond loop over scint. hits in a plane ------------------------------
//---------------------------------------------------------------------------------------------
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
GoodFlags flags;
fGoodFlags[itrack][ip].push_back(flags);
fGoodFlags[itrack][ip][iphit].onTrack = kFALSE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kFALSE;
fGoodFlags[itrack][ip][iphit].goodTdcNeg = kFALSE;
fGoodFlags[itrack][ip][iphit].goodTdcPos = kFALSE;
fTOFCalc.push_back(TOFCalc());
// Do we set back to false for each track, or just once per event?
fTOFCalc[ihhit].good_scin_time = kFALSE;
// These need a track index too to calculate efficiencies
fTOFCalc[ihhit].good_tdc_pos = kFALSE;
fTOFCalc[ihhit].good_tdc_neg = kFALSE;
fTOFCalc[ihhit].pindex = ip;
// ihhit ++;
// fRawIndex ++; // Is fRawIndex ever different from ihhit
Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
( fPlanes[ip]->GetZpos() + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 277
Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
( fPlanes[ip]->GetZpos() + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 278
if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 278
scinTrnsCoord = xHitCoord;
scinLongCoord = yHitCoord;
}
else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 281
scinTrnsCoord = yHitCoord;
scinLongCoord = xHitCoord;
}
else { return -1; } // Line 288
Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// ** Check if scin is on track
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
if ( fTOFPInfo[iphit].keep_pos ) { // 301
// ** Calculate time for each tube with a good tdc. 'pos' side first.
fTOFCalc[ihhit].good_tdc_pos = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
Double_t path = fPlanes[ip]->GetPosLeft() - scinLongCoord;
// * Convert TDC value to time, do pulse height correction, correction for
// * propogation of light thru scintillator, and offset.
Double_t time = ((THcHodoHit*)hodoHits->At(iphit))->GetPosCorrectedTime();
time = time - ( path / fHodoVelLight[fPIndex] );
fTOFPInfo[iphit].scin_pos_time = time;
} // check for good pos TDC condition
if ( fTOFPInfo[iphit].keep_neg ) { //
// ** Calculate time for each tube with a good tdc. 'pos' side first.
fTOFCalc[ihhit].good_tdc_neg = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
Double_t path = scinLongCoord - fPlanes[ip]->GetPosRight();
// * Convert TDC value to time, do pulse height correction, correction for
// * propogation of light thru scintillator, and offset.
Double_t time = ((THcHodoHit*)hodoHits->At(iphit))->GetNegCorrectedTime();
time = time - ( path / fHodoVelLight[fPIndex] );
fTOFPInfo[iphit].scin_neg_time = time;
} // check for good neg TDC condition
// ** Calculate ave time for scin and error.
if ( fTOFCalc[ihhit].good_tdc_pos ){
if ( fTOFCalc[ihhit].good_tdc_neg ){
fTOFCalc[ihhit].scin_time = ( fTOFPInfo[iphit].scin_pos_time +
fTOFPInfo[iphit].scin_neg_time ) / 2.;
fTOFCalc[ihhit].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.;
fTOFCalc[ihhit].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
fTOFCalc[ihhit].scin_time = fTOFPInfo[iphit].scin_pos_time;
fTOFCalc[ihhit].scin_sigma = fHodoPosSigma[fPIndex];
fTOFCalc[ihhit].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
if ( fTOFCalc[ihhit].good_tdc_neg ){
fTOFCalc[ihhit].scin_time = fTOFPInfo[iphit].scin_neg_time;
fTOFCalc[ihhit].scin_sigma = fHodoNegSigma[fPIndex];
fTOFCalc[ihhit].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
} // In h_tof.f this includes the following if condition for time at focal plane
// // because it is written in FORTRAN code
// c Get time at focal plane
if ( fTOFCalc[ihhit].good_scin_time ){
// scin_time_fp doesn't need to be an array
Double_t scin_time_fp = fTOFCalc[ihhit].scin_time -
( fPlanes[ip]->GetZpos() + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) /
TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
sumFPTime = sumFPTime + scin_time_fp;
nFPTime ++;
fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
fNPlaneTime[ip] ++;
fNScinHit[itrack] ++;
if ( ( fTOFCalc[ihhit].good_tdc_pos ) && ( fTOFCalc[ihhit].good_tdc_neg ) ){
fdEdX[itrack].push_back(0.0);
// --------------------------------------------------------------------------------------------
if ( fTOFCalc[ihhit].good_tdc_pos ){
if ( fTOFCalc[ihhit].good_tdc_neg ){
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Sqrt( TMath::Max( 0., ((THcHodoHit*)hodoHits->At(iphit))->GetPosADC() *
((THcHodoHit*)hodoHits->At(iphit))->GetNegADC() ) );
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Max( 0., ((THcHodoHit*)hodoHits->At(iphit))->GetPosADC() );
if ( fTOFCalc[ihhit].good_tdc_neg ){
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Max( 0., ((THcHodoHit*)hodoHits->At(iphit))->GetNegADC() );
fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
// --------------------------------------------------------------------------------------------
} // time at focal plane condition
} // on track else condition
// ** See if there are any good time measurements in the plane.
if ( fTOFCalc[ihhit].good_scin_time ){
fGoodPlaneTime[ip] = kTRUE;
fTOFCalc[ihhit].dedx = fdEdX[itrack][fNScinHit[itrack]-1];
} else {
fTOFCalc[ihhit].dedx = 0.0;
// Can this be done after looping over hits and planes?
if ( fGoodPlaneTime[2] ) theTrack->SetGoodPlane3( 1 );
if ( !fGoodPlaneTime[2] ) theTrack->SetGoodPlane3( 0 );
if ( fGoodPlaneTime[3] ) theTrack->SetGoodPlane4( 1 );
if ( !fGoodPlaneTime[3] ) theTrack->SetGoodPlane4( 0 );
ihhit ++;
} // Second loop over hits of a scintillator plane ends here
} // Loop over scintillator planes ends here
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// * * Fit beta if there are enough time measurements (one upper, one lower)
// From h_tof_fit
if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) &&
( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){
Double_t sumW = 0.;
Double_t sumT = 0.;
Double_t sumZ = 0.;
Double_t sumZZ = 0.;
Double_t sumTZ = 0.;
if (!fPlanes[ip])
return -1;
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
if ( fTOFCalc[ihhit].good_scin_time ) {
Double_t scinWeight = 1 / ( fTOFCalc[ihhit].scin_sigma * fTOFCalc[ihhit].scin_sigma );
Double_t zPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ihhit].hit_paddle % 2 ) *
fPlanes[ip]->GetDzpos() );
sumW += scinWeight;
sumT += scinWeight * fTOFCalc[ihhit].scin_time;
sumZ += scinWeight * zPosition;
sumZZ += scinWeight * ( zPosition * zPosition );
sumTZ += scinWeight * zPosition * fTOFCalc[ihhit].scin_time;
} // condition of good scin time
ihhit ++;
} // loop over hits of plane
} // loop over planes
Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
ihhit = 0;
for (Int_t ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes
if (!fPlanes[ip])
return -1;
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){ // Loop over hits of a plane
if ( fTOFCalc[ihhit].good_scin_time ){
Double_t zPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ihhit].hit_paddle % 2 ) *
fPlanes[ip]->GetDzpos() );
Double_t timeDif = ( fTOFCalc[ihhit].scin_time - t0 );
betaChiSq += ( ( zPosition / beta - timeDif ) *
( zPosition / beta - timeDif ) ) /
( fTOFCalc[ihhit].scin_sigma * fTOFCalc[ihhit].scin_sigma );
} // condition for good scin time
ihhit++;
} // loop over hits of a plane
} // loop over planes
Double_t pathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
} // condition for fTmpDenom
} // else condition for fTmpDenom
if ( nFPTime != 0 ){
timeAtFP[itrack] = ( sumFPTime / nFPTime );
//
// ---------------------------------------------------------------------------
Double_t FPTimeSum=0.0;
Int_t nFPTimeSum=0;
if ( fNPlaneTime[ip] != 0 ){
fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
FPTimeSum += fSumPlaneTime[ip];
nFPTimeSum += fNPlaneTime[ip];
fFPTime[ip] = 1000. * ( ip + 1 );
Double_t dedx=0.0;
for(UInt_t ih=0;ih<fTOFCalc.size();ih++) {
if(fTOFCalc[ih].good_scin_time) {
dedx = fTOFCalc[ih].dedx;
break;
}
}
theTrack->SetDedx(dedx);
theTrack->SetFPTime(fptime);
theTrack->SetBeta(beta);
theTrack->SetBetaChi2( betaChiSq );
theTrack->SetNPMT(nPmtHit[itrack]);
theTrack->SetFPTime( timeAtFP[itrack]);
} // Main loop over tracks ends here.
} // If condition for at least one track
//-----------------------------------------------------------------------
//
// Trnslation of h_track_tests.f file for tracking efficiency
//
//-----------------------------------------------------------------------
//************************now look at some hodoscope tests
// *second, we move the scintillators. here we use scintillator cuts to see
// *if a track should have been found.
for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
std::vector<Double_t> scin_temp;
fScinHitPaddle.push_back(scin_temp); // Create array of hits per plane
for (UInt_t ipaddle = 0; ipaddle < fNPaddle[0]; ipaddle++ ){
fScinHitPaddle[ip].push_back(0.0);
fScinHitPaddle[ip][ipaddle] = 0.0;
for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
if (!fPlanes[ip])
return -1;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// TClonesArray* scinPosTDC = fPlanes[ip]->GetPosTDC();
// TClonesArray* scinNegTDC = fPlanes[ip]->GetNegTDC();
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
fScinHitPaddle[ip][paddle] = 1;
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
}
}
// *next, look for clusters of hits in a scin plane. a cluster is a group of
// *adjacent scintillator hits separated by a non-firing scintillator.
// *Wwe count the number of three adjacent scintillators too. (A signle track
// *shouldn't fire three adjacent scintillators.
for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
// Planes ip = 0 = 1X
// Planes ip = 2 = 2X
if (!fPlanes[ip]) return -1;
fNClust.push_back(0);
fThreeScin.push_back(0);
}
// *look for clusters in x planes... (16 scins) !this assume both x planes have same
// *number of scintillators.
Int_t icount;
for (Int_t ip = 0; ip < 3; ip +=2 ) {
icount = 0;
if ( fScinHitPaddle[ip][0] == 1 )
icount ++;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[0] - 1; ipaddle++ ){
// !look for number of clusters of 1 or more hits
if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
icount ++;
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
} // Loop over paddles
fNClust[ip] = icount;
icount = 0;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[0] - 2; ipaddle++ ){
// !look for three or more adjacent hits
if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
icount ++;
} // Second loop over paddles
if ( icount > 0 )
fThreeScin[ip] = 1;
} // Loop over X plane
// *look for clusters in y planes... (10 scins) !this assume both y planes have same
// *number of scintillators.
for (Int_t ip = 1; ip < 4; ip +=2 ) {
// Planes ip = 1 = 1Y
// Planes ip = 3 = 2Y
if (!fPlanes[ip]) return -1;
icount = 0;
if ( fScinHitPaddle[ip][0] == 1 )
icount ++;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[1] - 1; ipaddle++ ){
// !look for number of clusters of 1 or more hits
if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
icount ++;
} // Loop over Y paddles
fNClust[ip] = icount;
icount = 0;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[1] - 2; ipaddle++ ){
// !look for three or more adjacent hits
if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
icount ++;
} // Second loop over Y paddles
if ( icount > 0 )
fThreeScin[ip] = 1;
}// Loop over Y planes
// *now put some "tracking" like cuts on the hslopes, based only on scins...
// *by "slope" here, I mean the difference in the position of scin hits in two
// *like-planes. For example, a track that those great straight through will
// *have a slope of zero. If it moves one scin over from s1x to s2x it has an
// *x-slope of 1... I pick the minimum slope if there are multiple scin hits.
Double_t bestXpScin = 100.0;
Double_t bestYpScin = 100.0;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[0]; ipaddle++ ){
for (Int_t ipaddle2 = 0; ipaddle2 < (Int_t) fNPaddle[0]; ipaddle2++ ){
if ( ( fScinHitPaddle[0][ipaddle] == 1 ) &&
( fScinHitPaddle[2][ipaddle2] == 1 ) ){
Double_t slope = TMath::Abs(ipaddle - ipaddle2);
if ( slope < bestXpScin ) {
bestXpScin = slope;
} // Second loop over X paddles
} // First loop over X paddles
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[1]; ipaddle++ ){
for (Int_t ipaddle2 = 0; ipaddle2 < (Int_t) fNPaddle[1]; ipaddle2++ ){
if ( ( fScinHitPaddle[1][ipaddle] == 1 ) &&
( fScinHitPaddle[3][ipaddle2] == 1 ) ){
Double_t slope = TMath::Abs(ipaddle - ipaddle2);
if ( slope < bestYpScin ) {
bestYpScin = slope;
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
} // Second loop over Y paddles
} // First loop over Y paddles
// *next we mask out the edge scintillators, and look at triggers that happened
// *at the center of the acceptance. To change which scins are in the mask
// *change the values of h*loscin and h*hiscin in htracking.param
// fGoodScinHits = 0;
for (Int_t ifidx = fxLoScin[0]; ifidx < (Int_t) fxHiScin[0]; ifidx ++ ){
fGoodScinHitsX.push_back(0);
}
// *first x plane. first see if there are hits inside the scin region
for (Int_t ifidx = fxLoScin[0]-1; ifidx < fxHiScin[0]; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){
fHitSweet1X = 1;
fSweet1XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fxLoScin[0]-1; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
}
for (Int_t ifidx = fxHiScin[0]; ifidx < (Int_t) fNPaddle[0]; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
}
// *second x plane. first see if there are hits inside the scin region
for (Int_t ifidx = fxLoScin[1]-1; ifidx < fxHiScin[1]; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){
fHitSweet2X = 1;
fSweet2XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fxLoScin[1]-1; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
}
for (Int_t ifidx = fxHiScin[1]; ifidx < (Int_t) fNPaddle[2]; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
}
// *first y plane. first see if there are hits inside the scin region
for (Int_t ifidx = fyLoScin[0]-1; ifidx < fyHiScin[0]; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){
fHitSweet1Y = 1;
fSweet1YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fyLoScin[0]-1; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
}
for (Int_t ifidx = fyHiScin[0]; ifidx < (Int_t) fNPaddle[1]; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
}
// *second y plane. first see if there are hits inside the scin region
for (Int_t ifidx = fyLoScin[1]-1; ifidx < fyHiScin[1]; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){
fHitSweet2Y = 1;
fSweet2YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fyLoScin[1]-1; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
}
for (Int_t ifidx = fyHiScin[1]; ifidx < (Int_t) fNPaddle[3]; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
}
fTestSum = fHitSweet1X + fHitSweet2X + fHitSweet1Y + fHitSweet2Y;
// * now define a 3/4 or 4/4 trigger of only good scintillators the value
// * is specified in htracking.param...
if ( fTestSum > fTrackEffTestNScinPlanes ){
fGoodScinHits = 1;
for (Int_t ifidx = fxLoScin[0]; ifidx < fxHiScin[0]; ifidx ++ ){
if ( fSweet1XScin == ifidx )
fGoodScinHitsX[ifidx] = 1;
}
}
// * require front/back hodoscopes be close to each other
if ( ( fGoodScinHits == 1 ) && ( fTrackEffTestNScinPlanes == 4 ) ){
if ( TMath::Abs( fSweet1XScin - fSweet2XScin ) > 3 )
fGoodScinHits = 0;
if ( TMath::Abs( fSweet1YScin - fSweet2YScin ) > 2 )
fGoodScinHits = 0;
}
if ( !fChern || !fShower ) {
return 0;
}
if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > fNormETot ) &&
( fChern->GetCerNPE() > fNCerNPE ) )
fScinShould = 1;
if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > fNormETot ) &&
( fChern->GetCerNPE() > fNCerNPE ) && ( tracks.GetLast() + 1 > 0 ) ) {
fScinDid = 1;
Gabriel Niculescu
committed
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
// GN: Return the index of a scintillator given the plane # and the paddle #
// This assumes that both planes and
// paddles start counting from 0!
// Result also counts from 0.
return fNPlanes*nPaddle+nPlane;
Gabriel Niculescu
committed
}
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nSide, Int_t nPlane, Int_t nPaddle ) {
return nSide*fMaxHodoScin+fNPlanes*nPaddle+nPlane-1;
}
//_____________________________________________________________________________
Double_t THcHodoscope::GetPathLengthCentral() {
return fPathLengthCentral;
}
ClassImp(THcHodoscope)
////////////////////////////////////////////////////////////////////////////////