Skip to content
Snippets Groups Projects
Commit 815a0674 authored by Zafar Ahmed's avatar Zafar Ahmed Committed by Stephen A. Wood
Browse files

Hodoscope efficiencies are added

parent d0c6b0f3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/THcHodoEff.cxx
+ 215
− 133
View file @ 815a0674
......@@ -9,6 +9,17 @@
// For now trying to emulate work done in h_scin_eff/h_scin_eff_shutdown //
///////////////////////////////////////////////////////////////////////////////
#include "THaEvData.h"
#include "THaCutList.h"
#include "VarDef.h"
#include "VarType.h"
#include "TClonesArray.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include "THcHodoEff.h"
#include "THaApparatus.h"
#include "THcHodoHit.h"
......@@ -107,6 +118,120 @@ THaAnalysisObject::EStatus THcHodoEff::Init( const TDatime& run_time )
return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcHodoEff::ReadDatabase( const TDatime& date )
{
// Read database. Gets variable needed for efficiency calculation
// Get # of planes and their z positions here.
fNPlanes = fHod->GetNPlanes();
fPlanes = new THcScintillatorPlane* [fNPlanes];
fPosZ = new Double_t[fNPlanes];
fSpacing = new Double_t[fNPlanes];
fCenterFirst = new Double_t[fNPlanes];
fNCounters = new Int_t[fNPlanes];
fHodoSlop = new Double_t[fNPlanes];
fHodoPosEffi = new Int_t[100];
fHodoNegEffi = new Int_t[100];
fHodoOrEffi = new Int_t[100];
fHodoAndEffi = new Int_t[100];
fStatTrk = new Int_t[100];
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip] = fHod->GetPlane(ip);
fPosZ[ip] = fPlanes[ip]->GetZpos() + 0.5*fPlanes[ip]->GetDzpos();
fSpacing[ip] = fPlanes[ip]->GetSpacing();
fCenterFirst[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
fNCounters[ip] = fPlanes[ip]->GetNelem();
}
char prefix[2];
prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
prefix[1] = '\0';
DBRequest list[]={
{"stat_slop", &fStatSlop, kDouble},
{"stat_maxchisq",&fMaxChisq, kDouble},
{"hodo_slop", fHodoSlop, kDouble, fNPlanes},
{0}
};
// fMaxShTrk = 0.05; // For cut on fraction of momentum seen in shower
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
cout << "\n\nTHcHodoEff::ReadDatabase nplanes=" << fHod->GetNPlanes() << endl;
// Setup statistics arrays
// Better method to put this in?
// These all need to be cleared in Begin
fHitPlane = new Int_t[fNPlanes];
fStatTrkDel.resize(fNPlanes);
fStatAndHitDel.resize(fNPlanes);
fStatPosHit.resize(fNPlanes);
fStatNegHit.resize(fNPlanes);
fStatAndHit.resize(fNPlanes);
fStatOrHit.resize(fNPlanes);
fBothGood.resize(fNPlanes);
fPosGood.resize(fNPlanes);
fNegGood.resize(fNPlanes);
for(Int_t ip=0;ip<fNPlanes;ip++) {
cout << "Plane = " << ip + 1 << " counters = " << fNCounters[ip] << endl;
fStatTrkDel[ip].resize(fNCounters[ip]);
fStatAndHitDel[ip].resize(fNCounters[ip]);
fStatPosHit[ip].resize(fNCounters[ip]);
fStatNegHit[ip].resize(fNCounters[ip]);
fStatAndHit[ip].resize(fNCounters[ip]);
fStatOrHit[ip].resize(fNCounters[ip]);
fBothGood[ip].resize(fNCounters[ip]);
fPosGood[ip].resize(fNCounters[ip]);
fNegGood[ip].resize(fNCounters[ip]);
for(Int_t ic=0;ic<fNCounters[ip];ic++) {
fStatTrkDel[ip][ic].resize(20); // Max this settable
fStatAndHitDel[ip][ic].resize(20); // Max this settable
fHodoPosEffi[fHod->GetScinIndex(ip,ic)] = 0;
fHodoNegEffi[fHod->GetScinIndex(ip,ic)] = 0;
fHodoOrEffi[fHod->GetScinIndex(ip,ic)] = 0;
fHodoAndEffi[fHod->GetScinIndex(ip,ic)] = 0;
fStatTrk[fHod->GetScinIndex(ip,ic)] = 0;
}
}
// Int_t fHodPaddles = fNCounters[0];
// gHcParms->Define(Form("%shodo_pos_hits[%d][%d]",fPrefix,fNPlanes,fHodPaddles),
// "Golden track's pos pmt hit",*&fStatPosHit);
Int_t fTotalPaddles = 100;
gHcParms->Define(Form("%shodo_pos_eff[%d]", prefix,fTotalPaddles), "Hodo positive effi",*fHodoPosEffi);
gHcParms->Define(Form("%shodo_neg_eff[%d]", prefix,fTotalPaddles), "Hodo negative effi",*fHodoNegEffi);
gHcParms->Define(Form("%shodo_or_eff[%d]", prefix,fTotalPaddles), "Hodo or effi", *fHodoOrEffi);
gHcParms->Define(Form("%shodo_and_eff[%d]", prefix,fTotalPaddles), "Hodo and effi", *fHodoAndEffi);
gHcParms->Define(Form("%shodo_gold_hits[%d]",prefix,fTotalPaddles), "Hodo golden hits", *fStatTrk);
return kOK;
}
//_____________________________________________________________________________
Int_t THcHodoEff::DefineVariables( EMode mode )
{
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
fEffiTest = 0;
gHcParms->Define(Form("hodoeffi"),"Testing effi",fEffiTest);
const RVarDef vars[] = {
// Move these into THcHallCSpectrometer using track fTracks
// {"effitestvar", "efficiency test var", "fEffiTest"},
// {"goldhodposhit", "pos pmt hit in hodo", "fStatPosHit"},
{ 0 }
};
return DefineVarsFromList( vars, mode );
// return kOK;
}
//_____________________________________________________________________________
Int_t THcHodoEff::Process( const THaEvData& evdata )
{
......@@ -116,6 +241,7 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
if( !IsOK() ) return -1;
// Project the golden track to each
// plane. Need to get track at Focal Plane, not tgt.
//
......@@ -127,7 +253,7 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
if(!theTrack) return 0;
Int_t trackIndex = theTrack->GetTrkNum()-1;
// May make these member variables
Double_t hitPos[fNPlanes];
Double_t hitDistance[fNPlanes];
......@@ -140,19 +266,25 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
// Should really have plane object self identify as X or Y
if(ip%2 == 0) { // X Plane
hitPos[ip] = theTrack->GetX() + theTrack->GetTheta()*fPosZ[ip];
} else { // Y Plane
hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip];
}
// Which counter does track pass through?
hitCounter[ip] = TMath::Min(TMath::Max(TMath::Nint((hitPos[ip] - fCenterFirst[ip])/fSpacing[ip]) + 1,1),fNCounters[ip]);
// How far from paddle center is track?
if(ip%2 == 0) { // X Plane
hitCounter[ip] = TMath::Max(
TMath::Min(
TMath::Nint((hitPos[ip]-fCenterFirst[ip])/
fSpacing[ip]+1),
TMath::Nint( fNCounters[ip] )),1);
hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) +
fCenterFirst[ip]);
} else { // Y Plane
hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) -
fCenterFirst[ip]);
hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip];
hitCounter[ip] = TMath::Max(
TMath::Min(
TMath::Nint((fCenterFirst[ip]-hitPos[ip])/
fSpacing[ip]+1),
TMath::Nint( fNCounters[ip] )),1);
hitDistance[ip] = hitPos[ip] -(fCenterFirst[ip] -
fSpacing[ip]*(hitCounter[ip]-1));
}
}
// Fill dpos histograms and set checkHit for each plane.
......@@ -186,27 +318,29 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
Int_t hitcounter = hitCounter[ip];
Double_t dist = hitDistance[ip];
if(TMath::Abs(dist) <= fStatSlop &&
theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
// && hsshtrk >= 0.05
Double_t delta = theTrack->GetDp();
Int_t idel = TMath::Floor(delta+10.0);
// Should
if(idel >=0 && idel < 20) {
fStatTrkDel[ip][hitcounter][idel]++;
theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq &&
theTrack->GetEnergy() >= 0.05 )
{
fStatTrk[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
Double_t delta = theTrack->GetDp();
Int_t idel = TMath::Floor(delta+10.0);
// Should
// if(idel >=0 && idel < 20) {
// fStatTrkDel[ip][hitcounter][idel]++;
// }
// lookat[ip] = TRUE;
}
// lookat[ip] = TRUE;
}
fHitPlane[ip] = 0;
}
}
// Is there a hit on or adjacent to paddle that track
// passes through?
// May collapse this loop into last
// record the hits as a "didhit" if track is near center of
// scintillator, the chisqared of the track is good and it is the
// first "didhit" in that plane.
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t hitcounter = hitCounter[ip];
Double_t dist = hitDistance[ip];
......@@ -219,13 +353,14 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
Bool_t onTrack, goodScinTime, goodTdcNeg, goodTdcPos;
fHod->GetFlags(trackIndex,ip,ihit,
onTrack, goodScinTime, goodTdcNeg, goodTdcPos);
if(TMath::Abs(dist) <= fStatSlop &&
TMath::Abs(hitcounter-counter) <= checkHit[ip] &&
fHitPlane[ip] == 0 &&
theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
// && hsshtrk >= 0.05
theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq &&
theTrack->GetEnergy() >= 0.05 ) {
fHitPlane[ip]++;
// Need to find out hgood_tdc_pos(igoldentrack,ihit) and neg
if(goodTdcPos) {
if(goodTdcNeg) { // Both fired
......@@ -233,131 +368,78 @@ Int_t THcHodoEff::Process( const THaEvData& evdata )
fStatNegHit[ip][hitcounter]++;
fStatAndHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
fHodoPosEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
fHodoNegEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
fHodoAndEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
Double_t delta = theTrack->GetDp();
Int_t idel = TMath::Floor(delta+10.0);
if(idel >=0 && idel < 20) {
fStatAndHitDel[ip][hitcounter][idel]++;
}
// Int_t idel = TMath::Floor(delta+10.0);
// if(idel >=0 && idel < 20) {
// fStatAndHitDel[ip][hitcounter][idel]++;
// }
} else {
fStatPosHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
fHodoPosEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
}
} else if(goodTdcNeg) {
} else if (goodTdcNeg) {
fStatNegHit[ip][hitcounter]++;
fStatOrHit[ip][hitcounter]++;
fHodoNegEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
fHodoOrEffi[fHod->GetScinIndex(ip,hitCounter[ip]-1)]++;
}
}
// Increment pos/neg/both fired. Track independent, so
// no chisquared cut, but note that only scintillators on the
// track are examined.
if(goodTdcPos) {
if(goodTdcNeg) {
fBothGood[ip][hitcounter]++;
} else {
fPosGood[ip][hitcounter]++;
// Increment pos/neg/both fired. Track independent, so
// no chisquared cut, but note that only scintillators on the
// track are examined.
if(goodTdcPos) {
if(goodTdcNeg) {
fBothGood[ip][hitcounter]++;
} else {
fPosGood[ip][hitcounter]++;
}
} else if (goodTdcNeg) {
fNegGood[ip][hitcounter]++;
}
// Determine if one or both PMTs had a good tdc
if(goodTdcPos && goodTdcNeg) {
goodTdcBothSides[ip] = kTRUE;
}
if(goodTdcPos || goodTdcNeg) {
goodTdcOneSide[ip] = kTRUE;
}
} else if (goodTdcNeg) {
fNegGood[ip][hitcounter]++;
}
// Determine if one or both PMTs had a good tdc
if(goodTdcPos && goodTdcNeg) {
goodTdcBothSides[ip] = kTRUE;
}
if(goodTdcPos || goodTdcNeg) {
goodTdcOneSide[ip] = kTRUE;
}
/*
For each plane, see of other 3 fired. This means that they were enough
to form a 3/4 trigger, and so the fraction of times this plane fired is
the plane trigger efficiency. NOTE: we only require a TDC hit, not a
TDC hit within the SCIN 3/4 trigger window, so high rates will make
this seem better than it is. Also, make sure we're not near the edge
of the hodoscope (at the last plane), using the same hhodo_slop param.
as for h_tof.f
NOTE ALSO: to make this check simpler, we are assuming that all planes
have identical active areas. y_scin = y_cent + y_offset, so shift track
position by offset for comparing to edges.
*/
// Need to add calculation and cuts on
// xatback and yatback in order to set the
// htrig_hododidflag, htrig_hodoshouldflag and otherthreehit flags
//
++fNevt;
}
}
/*
For each plane, see of other 3 fired. This means that they were enough
to form a 3/4 trigger, and so the fraction of times this plane fired is
the plane trigger efficiency. NOTE: we only require a TDC hit, not a
TDC hit within the SCIN 3/4 trigger window, so high rates will make
this seem better than it is. Also, make sure we're not near the edge
of the hodoscope (at the last plane), using the same hhodo_slop param. as for h_tof.f
NOTE ALSO: to make this check simpler, we are assuming that all planes
have identical active areas. y_scin = y_cent + y_offset, so shift track
position by offset for comparing to edges.
*/
// Need to add calculation and cuts on
// xatback and yatback in order to set the
// htrig_hododidflag, htrig_hodoshouldflag and otherthreehit flags
//
++fNevt;
return 0;
}
//_____________________________________________________________________________
Int_t THcHodoEff::ReadDatabase( const TDatime& date )
{
// Read database. Gets variable needed for efficiency calculation
// Get # of planes and their z positions here.
fNPlanes = fHod->GetNPlanes();
fPlanes = new THcScintillatorPlane* [fNPlanes];
fPosZ = new Double_t[fNPlanes];
fSpacing = new Double_t[fNPlanes];
fCenterFirst = new Double_t[fNPlanes];
fNCounters = new Int_t[fNPlanes];
fHodoSlop = new Double_t[fNPlanes];
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip] = fHod->GetPlane(ip);
fPosZ[ip] = fPlanes[ip]->GetZpos() + 0.5*fPlanes[ip]->GetDzpos();
fSpacing[ip] = fPlanes[ip]->GetSpacing();\
fCenterFirst[ip] = fPlanes[ip]->GetPosCenter(0);
fNCounters[ip] = fPlanes[ip]->GetNelem();
}
char prefix[2];
prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
prefix[1] = '\0';
DBRequest list[]={
{"stat_slop", &fStatSlop, kDouble},
{"stat_maxchisq",&fMaxChisq, kDouble},
{"hodo_slop", fHodoSlop, kDouble, fNPlanes},
{0}
};
// fMaxShTrk = 0.05; // For cut on fraction of momentum seen in shower
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
cout << "THcHodoEff::ReadDatabase nplanes=" << fHod->GetNPlanes() << endl;
// Setup statistics arrays
// Better method to put this in?
// These all need to be cleared in Begin
fHitPlane = new Int_t[fNPlanes];
fStatTrkDel.resize(fNPlanes);
fStatAndHitDel.resize(fNPlanes);
fStatPosHit.resize(fNPlanes);
fStatNegHit.resize(fNPlanes);
fStatAndHit.resize(fNPlanes);
fStatOrHit.resize(fNPlanes);
fBothGood.resize(fNPlanes);
fPosGood.resize(fNPlanes);
fNegGood.resize(fNPlanes);
for(Int_t ip=0;ip<fNPlanes;ip++) {
fStatTrkDel[ip].resize(fNCounters[ip]);
fStatAndHitDel[ip].resize(fNCounters[ip]);
fStatPosHit[ip].resize(fNCounters[ip]);
fStatNegHit[ip].resize(fNCounters[ip]);
fStatAndHit[ip].resize(fNCounters[ip]);
fStatOrHit[ip].resize(fNCounters[ip]);
fBothGood[ip].resize(fNCounters[ip]);
fPosGood[ip].resize(fNCounters[ip]);
fNegGood[ip].resize(fNCounters[ip]);
for(Int_t ic=0;ic<fNCounters[ip];ic++) {
fStatTrkDel[ip][ic].resize(20); // Max this settable
fStatAndHitDel[ip][ic].resize(20); // Max this settable
}
}
return kOK;
}
//_____________________________________________________________________________
ClassImp(THcHodoEff)
////////////////////////////////////////////////////////////////////////////////
src/THcHodoEff.h
+ 20
− 1
View file @ 815a0674
......@@ -7,6 +7,17 @@
// //
///////////////////////////////////////////////////////////////////////////////
#include "THaEvData.h"
#include "THaCutList.h"
#include "VarDef.h"
#include "VarType.h"
#include "TClonesArray.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include "THaPhysicsModule.h"
#include "THcHodoscope.h"
#include "THaSpectrometer.h"
......@@ -27,6 +38,8 @@ public:
protected:
virtual Int_t ReadDatabase( const TDatime& date);
virtual Int_t DefineVariables( EMode mode = kDefine );
/* Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle); */
// Data needed for efficiency calculation for one Hodoscope paddle
......@@ -41,13 +54,19 @@ protected:
// Information about the hodoscopes that we get from the
// THcHodoscope object
Int_t fEffiTest;
Int_t fNPlanes;
THcScintillatorPlane** fPlanes;
Double_t* fPosZ;
Double_t* fSpacing;
Double_t* fCenterFirst;
Int_t* fNCounters;
// Int_t fMaxNcounters;
// Int_t* fHodoPlnContHit;
Int_t* fHodoPosEffi;
Int_t* fHodoNegEffi;
Int_t* fHodoOrEffi;
Int_t* fHodoAndEffi;
Int_t* fStatTrk;
Double_t fStatSlop;
Double_t fMaxChisq;
Double_t* fHodoSlop;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment