#include "MWDC.h"
#include "WirePlane.h"
#include "SpecialWirePlane.h"
#include "Projection.h"
#include "Hit.h"
#include "Road.h"
#include "Helper.h"
#include "THaDetMap.h"
#include "THaEvData.h"
#include "THaTrack.h"
#include "TString.h"
#include "TMath.h"
#include "THashTable.h"
#include "TVector2.h"
#include "TDecompChol.h"
#include "TSystem.h"
#include "TThread.h"
#include "TCondition.h"
#include "TMutex.h"
#include <iostream>
#include <algorithm>
#include <numeric>
#include <map>
#include <string>
#ifdef TESTCODE
#include "TStopwatch.h"
#endif
using namespace std;
typedef string::size_type ssiz_t;
namespace {
class CSpair : public TObject {
public:
CSpair( UShort_t crate, UShort_t slot ) : fCrate(crate), fSlot(slot) {}
virtual ULong_t Hash() const {
UInt_t cs = static_cast<UInt_t>(fCrate)<<16 + static_cast<UInt_t>(fSlot);
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,16,0)
return TString::Hash( &cs, sizeof(cs) );
#else
return TMath::Hash( &cs, sizeof(cs) );
#endif
}
virtual Bool_t IsEqual( const TObject* obj ) const {
const CSpair* m;
if( !obj or !(m = dynamic_cast<const CSpair*>(obj)) ) return kFALSE;
return ( fCrate == m->fCrate and fSlot == m->fSlot );
}
UShort_t fCrate;
UShort_t fSlot;
};
class DAQmodule : public CSpair {
public:
DAQmodule( UShort_t crate, UShort_t slot, UInt_t model, Double_t res,
UInt_t nchan )
: CSpair(crate, slot), fModel(model), fResolution(res), fNchan(nchan) {}
virtual ~DAQmodule() {}
virtual void Copy( TObject& obj ) const {
TObject::Copy(obj);
DAQmodule* m = dynamic_cast<DAQmodule*>(&obj);
if( !m ) return;
m->fCrate = fCrate; m->fSlot = fSlot; m->fModel = fModel;
m->fResolution = fResolution, m->fNchan = fNchan;
}
virtual void Print( Option_t* ) const {
cout << "DAQmodule: "
<< " crate = " << fCrate
<< " slot = " << fSlot
<< " model = " << fModel
<< " res = " << fResolution
<< " nchan = " << fNchan
<< endl;
}
UInt_t fModel;
Double_t fResolution;
UInt_t fNchan;
};
}
namespace TreeSearch {
typedef vector<WirePlane*>::size_type vwsiz_t;
typedef vector<WirePlane*>::iterator vwiter_t;
typedef vector<Projection*>::size_type vpsiz_t;
typedef vector<vector<Int_t> >::iterator vviter_t;
const Double_t kBig = 1e38;
#define ALL(c) (c).begin(), (c).end()
static void DoTrack( void* ptr );
struct TrackThread {
Projection* proj;
Int_t* status;
Int_t* running;
TMutex* start_m;
TMutex* done_m;
TCondition* start;
TCondition* done;
TThread* thread;
TrackThread()
: proj(0), status(0), running(0), start(0), done(0), thread(0) {}
~TrackThread() {
if( thread ) {
TThread::Delete(thread);
delete thread;
}
}
};
class ThreadCtrl {
friend THaAnalysisObject::EStatus MWDC::Init(const TDatime&);
friend Int_t MWDC::CoarseTrack(TClonesArray&);
public:
ThreadCtrl() : fTrackStatus(0), fTrackToDo(0),
fTrackStartM(new TMutex), fTrackDoneM(new TMutex),
fTrackStart(new TCondition(fTrackStartM)),
fTrackDone(new TCondition(fTrackDoneM)) {}
~ThreadCtrl() {
if( !fTrack.empty() ) {
fTrackStartM->Lock();
fTrackToDo = 0;
for( vector<TrackThread>::iterator it = fTrack.begin(); it !=
fTrack.end(); ++it ) {
TThread* th = (*it).thread;
if( th and th->GetState() == TThread::kRunningState )
SETBIT(fTrackToDo, (*it).proj->GetType());
}
SETBIT(fTrackToDo, 31);
fTrackStart->Broadcast();
fTrackDoneM->Lock();
fTrackStartM->UnLock();
while( true ) {
Int_t ret = fTrackDone->Wait();
if( ret == 0 and fTrackToDo == BIT(31) )
break;
}
fTrackDoneM->UnLock();
}
delete fTrackDone;
delete fTrackStart;
delete fTrackDoneM;
delete fTrackStartM;
}
TThread* AddTrackThread( Projection* proj ) {
assert( proj );
fTrack.push_back( TrackThread() );
TrackThread* t = &fTrack.back();
t->proj = proj;
t->status = &fTrackStatus;
t->running = &fTrackToDo;
t->start_m = fTrackStartM;
t->done_m = fTrackDoneM;
t->start = fTrackStart;
t->done = fTrackDone;
string tn;
tn = "trk_"; tn.append( proj->GetName() );
t->thread = new TThread( tn.c_str(), DoTrack, (void*)t );
return t->thread;
}
private:
vector<TrackThread> fTrack;
Int_t fTrackStatus;
Int_t fTrackToDo;
TMutex* fTrackStartM;
TMutex* fTrackDoneM;
TCondition* fTrackStart;
TCondition* fTrackDone;
};
MWDC::MWDC( const char* name, const char* desc, THaApparatus* app )
: THaTrackingDetector(name,desc,app), fCrateMap(0),
fMaxThreads(1), fThreads(0),
f3dMatchvalScalefact(1), f3dMatchCut(0), fMinNdof(1),
fFailNhits(0), fFailNpat(0)
#ifdef TESTCODE
, fNcombos(0), fN3dFits(0), fEvNum(0)
#endif
, fuseshower(0), fshower(NULL)
{
fRefMap = new THaDetMap;
fProj.resize( kTypeEnd, 0 );
Double_t p_angle[] = { -60.0, 60.0, 0.0, 90.0 };
const char* p_name[] = { "u", "v", "x", "y" };
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
Projection* proj = new Projection( type,
p_name[type],
p_angle[type]*TMath::DegToRad(),
this
);
if( !proj or proj->IsZombie() ) {
Error( Here("MWDC"), "Error creating projection %s. Call expert.",
p_name[type] );
MakeZombie();
return;
}
fProj[type] = proj;
}
}
MWDC::~MWDC()
{
if (fIsSetup)
RemoveVariables();
delete fThreads;
if( fMaxThreads > 1 )
gSystem->Unload("libThread");
DeleteContainer( fPlanes );
DeleteContainer( fProj );
delete fshower;
delete fRefMap;
}
void MWDC::Clear( Option_t* opt )
{
THaTrackingDetector::Clear(opt);
if( !opt or *opt != 'I' ) {
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane )
fPlanes[iplane]->Clear(opt);
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type )
fProj[type]->Clear(opt);
if( fRefMap->GetSize() > 0 )
fRefTime.assign( fRefMap->GetSize(), kBig );
}
fFailNpat = fFailNhits = 0;
#ifdef TESTCODE
size_t nbytes = (char*)&t_coarse - (char*)&fNcombos + sizeof(t_coarse);
memset( &fNcombos, 0, nbytes );
#endif
}
Int_t MWDC::Decode( const THaEvData& evdata )
{
static const char* const here = "Decode";
#ifdef TESTCODE
fEvNum = evdata.GetEvNum();
#endif
#ifdef VERBOSE
if( fDebug > 0 ) {
cout << "============== Event ";
#ifdef TESTCODE
cout << "# " << fEvNum << " ";
#endif
cout << "==============" << endl;
}
#endif
if (fuseshower==1){
fshower->Decode(evdata);
TClonesArray tracks;
fshower->CoarseProcess(tracks);
fshower->FineProcess(tracks);
}
for( Int_t imod = 0; imod < fRefMap->GetSize(); ++imod ) {
THaDetMap::Module* d = fRefMap->GetModule(imod);
Int_t chan = d->lo;
Int_t nhits = evdata.GetNumHits( d->crate, d->slot, chan );
if( nhits > 0 ) {
Int_t data = evdata.GetData( d->crate, d->slot, chan, nhits-1 );
if( nhits > 1 ) {
Warning( Here(here), "%d hits on reference channel %d module %d",
nhits, chan, imod );
}
fRefTime[imod] = d->resolution * data;
} else {
Warning( Here(here), "No hits on reference channel %d module %d.",
chan, imod );
fRefTime[imod] = kBig;
}
}
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
Int_t nhits = fProj[type]->Decode( evdata );
if( nhits < 0 ) {
fFailNhits = 1;
continue;
}
if( TestBit(kDoCoarse) )
fProj[type]->FillHitpattern();
}
return 0;
}
class AddFitCoord
{
public:
AddFitCoord() {}
void operator() ( Road* rd, Road::Point* p, const vector<Double_t>& coef )
{
const Projection* proj = rd->GetProjection();
Double_t slope2d = rd->GetSlope();
Double_t trkpos2d = rd->GetPos() + p->z * slope2d;
Double_t cosa = proj->GetCosAngle();
Double_t sina = proj->GetSinAngle();
Double_t slope = coef[1]*cosa + coef[3]*sina;
Double_t pos = coef[0]*cosa + coef[2]*sina + slope * p->z;
p->hit->GetWirePlane()->AddFitCoord
( FitCoord( p->hit, rd, p->x, trkpos2d, slope2d, pos, slope ));
}
};
#ifdef VERBOSE
class PrintFitPoint
{
public:
PrintFitPoint() {}
void operator() ( Road*, Road::Point* p, const vector<Double_t>& )
{
cout << p->hit->GetWirePlane()->GetName() << " "
<< "z = " << p->z << " x = " << p->x << endl;
}
};
#endif
class FillFitMatrix
{
public:
FillFitMatrix( TMatrixDSym& AtA, TVectorD& Aty )
: fAtA(AtA), fAty(Aty), fNpoints(0)
{
assert( fAtA.GetNrows() == 4 && fAty.GetNrows() == 4 );
}
void operator() ( Road* rd, Road::Point* p, const vector<Double_t>& )
{
const Projection* proj = rd->GetProjection();
Double_t cosa = proj->GetCosAngle();
Double_t sina = proj->GetSinAngle();
Double_t Ai[4] = { cosa, cosa * p->z, sina, sina * p->z };
Double_t s2 = 1.0/(p->res()*p->res());
for( int j = 0; j<4; ++j ) {
for( int k = j; k<4; ++k ) {
fAtA(j,k) += Ai[j] * s2 * Ai[k];
}
fAty(j) += Ai[j] * s2 * p->x;
}
++fNpoints;
}
Int_t GetNpoints() const { return fNpoints; }
private:
TMatrixDSym& fAtA;
TVectorD& fAty;
Int_t fNpoints;
};
class CalcChisquare
{
public:
CalcChisquare() : fChi2(0) {}
void operator() ( Road* rd, Road::Point* p, const vector<Double_t>& coef )
{
const Projection* proj = rd->GetProjection();
Double_t cosa = proj->GetCosAngle();
Double_t sina = proj->GetSinAngle();
Double_t slope = coef[1]*cosa + coef[3]*sina;
Double_t pos = coef[0]*cosa + coef[2]*sina + slope*p->z;
Double_t diff = (pos - p->x) / p->res();
fChi2 += diff*diff;
}
void Clear() { fChi2 = 0; }
Double_t GetResult() const { return fChi2; }
private:
Double_t fChi2;
};
template< typename Action >
Action MWDC::ForAllTrackPoints( const Rvec_t& roads,
const vector<Double_t>& coef, Action action )
{
for( Rvec_t::const_iterator it = roads.begin(); it != roads.end(); ++it ) {
Road* rd = *it;
const Road::Pvec_t& points = rd->GetPoints();
for( Road::Pvec_t::const_iterator it2 = points.begin();
it2 != points.end(); ++it2 ) {
Road::Point* p = *it2;
action( rd, p, coef );
}
}
return action;
}
void MWDC::FitErrPrint( Int_t err ) const
{
static const char* const here = "FitTrack";
#ifdef TESTCODE
Error( Here(here), "Event %d: Failure fitting 3D track, err = %d. "
"Should never happen. Call expert.", fEvNum, err );
#else
Error( Here(here), "Failure fitting 3D track, err = %d. "
"Should never happen. Call expert.", err );
#endif
}
Int_t MWDC::FitTrack( const Rvec_t& roads, vector<Double_t>& coef,
Double_t& chi2, TMatrixDSym* coef_covar ) const
{
TMatrixDSym AtA(4);
TVectorD Aty(4);
FillFitMatrix f = ForAllTrackPoints(roads, coef, FillFitMatrix(AtA, Aty));
Int_t npoints = f.GetNpoints();
assert( npoints > 4 );
if( npoints <=4 ) return -1;
for( int j = 0; j<4; ++j )
for( int k = j+1; k<4; ++k )
AtA(k,j) = AtA(j,k);
TDecompChol chol(AtA);
Bool_t ok = chol.Solve(Aty);
assert(ok);
if( !ok ) return -2;
assert( Aty.GetNrows() == 4 );
coef.assign( Aty.GetMatrixArray(), Aty.GetMatrixArray()+Aty.GetNrows() );
#ifdef VERBOSE
if( fDebug > 2 ) {
cout << "Points in 3D fit:" << endl;
ForAllTrackPoints( roads, coef, PrintFitPoint() );
}
#endif
chi2 = ForAllTrackPoints(roads, coef, CalcChisquare()).GetResult();
if( coef_covar ) {
if( coef_covar->GetNrows() != 4 )
coef_covar->ResizeTo(4,4);
Bool_t ok = chol.Invert(*coef_covar);
assert(ok);
if( !ok ) return -3;
}
#ifdef VERBOSE
if( fDebug > 1 )
cout << "3D fit: x/y = " << coef[0] << "/" << coef[2] << " "
<< "mx/my = " << coef[1] << "/" << coef[3] << " "
<< "ndof = " << npoints-4 << " rchi2 = " << chi2/(double)(npoints-4)
<< endl;
#endif
return npoints-4;
}
void MWDC::FindNearestHits( WirePlane* wp, const THaTrack* track,
const Rvec_t& roads ) const
{
assert( !roads.empty() );
Double_t z = wp->GetZ();
Double_t cosa = wp->GetProjection()->GetCosAngle();
Double_t sina = wp->GetProjection()->GetSinAngle();
Double_t slope = track->GetTheta()*cosa + track->GetPhi()*sina;
Double_t x = track->GetX()*cosa + track->GetY()*sina + slope*z;
Double_t dmin = kBig;
Double_t pmin = kBig;
Hit* hmin = 0;
const TSeqCollection* hits = wp->GetHits();
Int_t first = 0, last = hits->GetSize();
const TObjArray* arr = dynamic_cast<const TObjArray*>(hits);
if( arr )
last = arr->GetLast()+1;
Int_t len = last - first;
Int_t half, middle;
while( len > 0 ) {
half = len >> 1;
middle = first + half;
if( static_cast<Hit*>(hits->At(middle))->GetWirePos() < x ) {
first = middle + 1;
len -= half + 1;
} else
len = half;
}
if( last > 0 ) {
assert( first <= last );
Hit *hnext = 0, *hprev = 0;
if( first < last ) {
hnext = static_cast<Hit*>(hits->At(first));
assert( hnext->GetWirePos() >= x );
}
if( first > 0 ) {
hprev = static_cast<Hit*>(hits->At(first-1));
assert( hprev->GetWirePos() < x );
if( hnext ) {
assert( hprev->GetWireNum() < hnext->GetWireNum() );
if( hprev->GetWireNum() + 1 < hnext->GetWireNum() ) {
if( x - hprev->GetWirePos() < hnext->GetWirePos() - x )
hnext = 0;
else
hprev = 0;
}
}
}
if( hnext ) {
hmin = hnext;
pmin = hnext->GetPosL();
dmin = TMath::Abs(pmin-x);
Int_t i = first;
Hit* h;
while( ++i < last and (h = static_cast<Hit*>(hits->At(i)))->GetWireNum()
== hnext->GetWireNum() ) {
Double_t d = TMath::Abs(h->GetPosL()-x);
if( d < dmin ) {
dmin = d;
hmin = h;
pmin = h->GetPosL();
}
}
}
if( hprev ) {
Double_t d = TMath::Abs(hprev->GetPosR()-x);
if( !hmin or d < dmin ) {
dmin = d;
hmin = hprev;
pmin = hprev->GetPosR();
}
Int_t i = first-1;
Hit* h;
while( --i >= 0 and (h = static_cast<Hit*>(hits->At(i)))->GetWireNum()
== hprev->GetWireNum() ) {
d = TMath::Abs(h->GetPosR()-x);
if( d < dmin ) {
dmin = d;
hmin = h;
pmin = h->GetPosR();
}
}
}
}
Road* rd = 0;
Int_t k = min( roads.size(), (Rvec_t::size_type)wp->GetType() );
do {
if( roads[k]->GetProjection()->GetType() > wp->GetType() )
--k;
else {
if( roads[k]->GetProjection() == wp->GetProjection() )
rd = roads[k];
break;
}
} while( k>=0 );
Double_t slope2d = rd ? rd->GetSlope() : kBig;
Double_t pos2d = rd ? rd->GetPos() + z * slope2d : kBig;
wp->AddFitCoord( FitCoord(hmin, rd, pmin, pos2d, slope2d, x, slope) );
}
THaTrack* MWDC::NewTrack( TClonesArray& tracks, const FitRes_t& fit_par )
{
Double_t x = fit_par.coef[0];
Double_t xp = fit_par.coef[1];
Double_t y = fit_par.coef[2];
Double_t yp = fit_par.coef[3];
THaTrack* newTrack = AddTrack( tracks, x, y, xp, yp );
assert( newTrack );
newTrack->SetD( x, y, xp, yp );
newTrack->SetChi2( fit_par.chi2, fit_par.ndof );
ForAllTrackPoints( *fit_par.roads, fit_par.coef, AddFitCoord() );
for( Wpvec_t::const_iterator it = fCalibPlanes.begin(); it !=
fCalibPlanes.end(); ++it ) {
FindNearestHits( (*it), newTrack, *fit_par.roads );
}
return newTrack;
}
class CheckTypes : public unary_function<Road*,void>
{
public:
CheckTypes( UInt_t req ) : fReq(req), fActive(0) {}
void operator() ( const Road* rd )
{ fActive |= 1U << rd->GetProjection()->GetType(); }
operator bool() const { return (fActive & fReq) == fReq; }
bool operator!() const { return (fActive & fReq) != fReq; }
UInt_t GetTypes() const { return fActive; }
private:
UInt_t fReq;
UInt_t fActive;
};
class AnySharedHits : public unary_function<Road*,bool>
{
public:
AnySharedHits() {}
bool operator() ( const Road* rd )
{
const Hset_t& test_hits = rd->GetHits();
for( Hset_t::const_iterator it = test_hits.begin(); it != test_hits.end();
++it ) {
if( fHits.find(*it) != fHits.end() )
return true;
}
return false;
}
const AnySharedHits& use( const Rset_t& tuple )
{
fHits.clear();
for( Rset_t::const_iterator it = tuple.begin(); it != tuple.end(); ++it ) {
const Road* rd = *it;
fHits.insert( ALL(rd->GetHits()) );
}
return *this;
}
private:
Hset_t fHits;
};
class MWDC::TrackFitWeight
{
public:
TrackFitWeight( const MWDC::FitRes_t& fit_par ) :
fNdof(fit_par.ndof), fChi2(fit_par.chi2) { assert(fNdof); }
bool operator<( const TrackFitWeight& rhs ) const
{
if( fNdof > rhs.fNdof ) return true;
if( fNdof < rhs.fNdof ) return false;
return ( fChi2 < rhs.fChi2 );
}
operator double() const { return fChi2/(double)fNdof; }
private:
UInt_t fNdof;
Double_t fChi2;
};
template< typename Container, typename Weight, typename TestF,
typename QuitF > Double_t
OptimalN( const Container& choices, const multimap<Weight,Container>& weights,
vector<Container>& picks, TestF testf, QuitF quitf )
{
picks.clear();
Container choices_left(choices);
double wsum = 0;
typename multimap<Weight,Container>::const_iterator it = weights.begin();
for( ; it != weights.end(); ++it ) {
const Container& tuple = (*it).second;
if( includes(ALL(choices_left), ALL(tuple)) ) {
picks.push_back( tuple );
wsum += (*it).first;
Container choices_less_tuple;
set_difference( ALL(choices_left), ALL(tuple),
inserter(choices_less_tuple, choices_less_tuple.end()) );
Container new_choices_left;
remove_copy_if( ALL(choices_less_tuple),
inserter(new_choices_left, new_choices_left.end()),
testf.use(tuple) );
if( !for_each(ALL(new_choices_left), quitf) )
break;
choices_left.swap( new_choices_left );
}
}
return wsum;
}
inline
void MWDC::Add3dMatch( const Rvec_t& selected, Double_t matchval,
list< pair<Double_t,Rvec_t> >& combos_found,
Rset_t& unique_found ) const
{
combos_found.push_back( make_pair(matchval,selected) );
unique_found.insert( ALL(selected) );
#ifdef VERBOSE
if( fDebug > 3 )
cout << "ACCEPTED" << endl;
#endif
}
UInt_t MWDC::MatchRoads( vector<Rvec_t>& roads,
list< pair<Double_t,Rvec_t> >& combos_found,
Rset_t& unique_found )
{
vector<Rvec_t>::size_type nproj = roads.size();
assert( nproj >= 2 );
combos_found.clear();
unique_found.clear();
UInt_t ncombos = accumulate( ALL(roads), (UInt_t)1, SizeMul<Rvec_t>() );
UInt_t nfound = 0;
bool fast_3d = TestBit(k3dFastMatch);
#ifdef VERBOSE
if( fDebug > 0 ) {
cout << "Matching ";
for( vector<Rvec_t>::size_type i = 0; i < nproj; ++i ) {
cout << roads[i].size();
if( i+1 < nproj ) cout << "x";
}
cout << " track projections in 3D (";
if( fast_3d ) cout << "fast";
else cout << "generic";
cout << " algo, " << ncombos << " combinations):" << endl;
}
#endif
#ifdef TESTCODE
fNcombos = ncombos;
#endif
if( ncombos == 0 )
return 0;
Rvec_t selected;
selected.reserve(nproj);
UInt_t qxnum=17;
for (UInt_t i=0;i!=fPlanes.size();i++){
if (fPlanes[i]->IsSpecial()){
}else{
qxnum = i;
}
}
Double_t zback = fPlanes[qxnum]->GetZ();
if( fast_3d ) {
assert( nproj == 3 && fProj.size() == 3 );
Prvec_t::iterator ip = fProj.begin();
Double_t su = (*ip)->GetSinAngle();
Double_t cu = (*ip)->GetCosAngle();
++ip;
Double_t sv = (*ip)->GetSinAngle();
Double_t cv = (*ip)->GetCosAngle();
Double_t inv_denom = 1.0/(sv*cu-su*cv);
su *= inv_denom; cu *= inv_denom; sv *= inv_denom; cv *= inv_denom;
++ip;
Double_t xax_x = ((*ip)->GetAxis()).X();
Double_t xax_y = ((*ip)->GetAxis()).Y();
Road* tuple[3];
UInt_t nrd[2] = { roads[0].size(), roads[1].size() };
UInt_t ird[2];
WirePlane *front_plane = fPlanes.front(), *back_plane = fPlanes[qxnum];
sort( ALL(roads[2]), Road::PosIsLess() );
Road::PosIsNear pos_near( TMath::Sqrt(f3dMatchCut) );
Double_t matchval = 0.0;
ird[0] = 0;
while( ird[0] != nrd[0] ) {
tuple[0] = roads[0][ird[0]];
Double_t uf = tuple[0]->GetPos();
Double_t ub = tuple[0]->GetPos(zback);
Double_t usf = uf * sv;
Double_t ucf = uf * cv;
Double_t usb = ub * sv;
Double_t ucb = ub * cv;
ird[1] = 0;
while( ird[1] != nrd[1] ) {
tuple[1] = roads[1][ird[1]];
Double_t v = tuple[1]->GetPos();
Double_t xf = usf - v * su;
Double_t yf = v * cu - ucf;
if( front_plane->Contains(xf,yf) ) {
v = tuple[1]->GetPos(zback);
Double_t xb = usb - v * su;
Double_t yb = v * cu - ucb;
if( back_plane->Contains(xb,yb) ) {
Double_t pf = xf*xax_x + yf*xax_y;
Double_t pb = xb*xax_x + yb*xax_y;
pair<Rvec_t::iterator,Rvec_t::iterator> range =
equal_range( ALL(roads[2]), pf, pos_near );
for( Rvec_t::iterator it=range.first; it != range.second; ++it ) {
tuple[2] = *it;
Double_t d1 = tuple[2]->GetPos() - pf;
Double_t d2 = tuple[2]->GetPos(zback) - pb;
matchval = d1*d1 + d2*d2;
#ifdef VERBOSE
if( fDebug > 3 ) {
if( matchval < f3dMatchCut || fDebug > 4 ) {
cout << tuple[0]->GetProjection()->GetName()
<< tuple[1]->GetProjection()->GetName()
<< " front = " << "(" << xf << "," << yf << ")" << endl;
cout << "front x = (" << tuple[2]->GetPos() * xax_x << ","
<< tuple[2]->GetPos() * xax_y << ")" << endl;
cout << "front dist = " << d1 << endl;
cout << tuple[0]->GetProjection()->GetName()
<< tuple[1]->GetProjection()->GetName()
<< " back = " << "(" << xb << "," << yb << ")" << endl;
cout << "back x = ("
<< tuple[2]->GetPos(zback) * xax_x << ","
<< tuple[2]->GetPos(zback) * xax_y << ")" << endl;
cout << "back dist = " << d2 << endl;
cout << "matchval = " << matchval*f3dMatchvalScalefact
<< endl;
}
}
#endif
if( matchval < f3dMatchCut ) {
++nfound;
selected.assign( tuple, tuple+3 );
Add3dMatch( selected, matchval, combos_found, unique_found );
}
}
}
}
++ird[1];
}
++ird[0];
}
} else {
vector<TVector2> fxpts, bxpts;
fxpts.reserve( nproj*(nproj-1)/2 );
bxpts.reserve( nproj*(nproj-1)/2 );
for( UInt_t i = 0; i < ncombos; ++i ) {
Double_t matchval = 0.0;
NthCombination( i, roads, selected );
assert( selected.size() == nproj );
fxpts.clear();
bxpts.clear();
TVector2 fctr, bctr;
for( Rvec_t::iterator it1 = selected.begin();
it1 != selected.end(); ++it1 ) {
for( Rvec_t::iterator it2 = it1+1; it2 != selected.end(); ++it2 ) {
fxpts.push_back( (*it1)->Intersect(*it2, 0.0) );
bxpts.push_back( (*it1)->Intersect(*it2, zback) );
fctr += fxpts.back();
bctr += bxpts.back();
#ifdef VERBOSE
if( fDebug > 3 ) {
cout << (*it1)->GetProjection()->GetName()
<< (*it2)->GetProjection()->GetName()
<< " front(" << fxpts.size() << ") = ";
fxpts.back().Print();
cout << (*it1)->GetProjection()->GetName()
<< (*it2)->GetProjection()->GetName()
<< " back (" << bxpts.size() << ") = ";
bxpts.back().Print();
}
#endif
}
}
assert( fxpts.size() <= nproj*(nproj-1)/2 );
assert( bxpts.size() == fxpts.size() );
fctr /= static_cast<Double_t>( fxpts.size() );
if( !fPlanes.front()->Contains(fctr) )
continue;
bctr /= static_cast<Double_t>( fxpts.size() );
if( !fPlanes[qxnum]->Contains(bctr) )
continue;
for( vector<TVector2>::size_type k = 0; k < fxpts.size(); ++k ) {
matchval += (fxpts[k]-fctr).Mod2() + (bxpts[k]-bctr).Mod2();
}
#ifdef VERBOSE
if( fDebug > 3 ) {
cout << "fctr = "; fctr.Print();
cout << "bctr = "; bctr.Print();
cout << "matchval = " << matchval << endl;
}
#endif
if( matchval < f3dMatchCut ) {
++nfound;
Add3dMatch( selected, matchval, combos_found, unique_found );
}
}
}
#ifdef VERBOSE
if( fDebug > 0 ) {
cout << nfound << " match";
if( nfound != 1 )
cout << "es";
cout << " found" << endl;
}
#endif
return nfound;
}
inline
Bool_t MWDC::PassTrackCuts( const FitRes_t& fit_par ) const
{
if( fit_par.ndof < fMinNdof )
return false;
pdbl_t chi2_interval = GetChisqLimits(fit_par.ndof);
if( fit_par.chi2 < chi2_interval.first or
fit_par.chi2 > chi2_interval.second )
return false;
return true;
}
static void DoTrack( void* ptr )
{
TrackThread* arg = reinterpret_cast<TrackThread*>(ptr);
bool terminate = false;
arg->start_m->Lock();
while( !terminate ) {
while( true ) {
Int_t ret = arg->start->Wait();
if( ret == 0 and TESTBIT(*arg->running, arg->proj->GetType()) )
break;
}
terminate = TESTBIT(*arg->running, 31);
arg->start_m->UnLock();
Int_t nrd = 0;
if( !terminate )
nrd = arg->proj->Track();
arg->done_m->Lock();
if( nrd < 0 ) {
*arg->status = 1;
}
assert( TESTBIT(*arg->running, arg->proj->GetType()) );
CLRBIT( *arg->running, arg->proj->GetType() );
if( (*arg->running & ~BIT(31)) == 0 )
arg->done->Signal();
if( !terminate )
arg->start_m->Lock();
arg->done_m->UnLock();
}
}
Int_t MWDC::CoarseTrack( TClonesArray& tracks )
{
if( fFailNhits )
return -1;
if( !TestBit(kDoCoarse) )
return 0;
#ifdef TESTCODE
TStopwatch timer, timer_tot;
#endif
Int_t err = 0;
bool do_thread = (fMaxThreads > 1);
if( do_thread ) {
fThreads->fTrackDoneM->Lock();
fThreads->fTrackStatus = 0;
Int_t all_todo = BIT(kTypeEnd)-1 & ~(BIT(kTypeBegin)-1);
Int_t max_todo = BIT(fMaxThreads)-1;
Int_t start = kTypeBegin;
while( start < kTypeEnd ) {
Int_t now_todo = (max_todo << start) & all_todo;
fThreads->fTrackStartM->Lock();
fThreads->fTrackToDo = now_todo;
fThreads->fTrackStart->Broadcast();
fThreads->fTrackStartM->UnLock();
while( true ) {
Int_t ret = fThreads->fTrackDone->Wait();
if( ret == 0 and fThreads->fTrackToDo == 0 )
break;
}
start += fMaxThreads;
}
err = fThreads->fTrackStatus;
fThreads->fTrackDoneM->UnLock();
} else {
for( Prvec_t::iterator it = fProj.begin(); it != fProj.end(); ++it ) {
Int_t ret = (*it)->Track();
if( ret < 0 )
err = 1;
}
}
if( err != 0 ) {
fFailNpat = 1;
return -1;
}
vector<Rvec_t>::size_type nproj = 0;
vector<Rvec_t> roads;
roads.reserve( kTypeEnd-kTypeBegin );
UInt_t found_types = 0;
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
Projection* proj = fProj[type];
Int_t nrd = proj->GetNgoodRoads();
if( nrd > 0 ) {
++nproj;
found_types |= 1U << type;
roads.push_back( Rvec_t() );
roads.back().reserve(nrd);
for( UInt_t i = 0; i < proj->GetNroads(); ++i ) {
Road* rd = proj->GetRoad(i);
assert(rd);
if( rd->IsGood() )
roads.back().push_back(rd);
}
}
}
assert( roads.size() == nproj );
#ifdef TESTCODE
t_track = 1e6*timer.RealTime();
timer.Start();
#endif
if( nproj >= 3 ) {
list< pair<Double_t,Rvec_t> > road_combos;
Rset_t unique_found;
UInt_t nfits = MatchRoads( roads, road_combos, unique_found );
#ifdef TESTCODE
t_3dmatch = 1e6*timer.RealTime();
timer.Start();
#endif
FitRes_t fit_par;
fit_par.coef.reserve(4);
if( nfits == 1 ) {
fit_par.matchval = road_combos.front().first;
Rvec_t& these_roads = road_combos.front().second;
fit_par.roads = &these_roads;
fit_par.ndof = FitTrack( these_roads, fit_par.coef, fit_par.chi2 );
if( fit_par.ndof > 0 ) {
if( PassTrackCuts(fit_par) )
NewTrack( tracks, fit_par );
}
else
FitErrPrint( fit_par.ndof );
}
else if( nfits > 0 ) {
roads.clear();
typedef map<Rset_t, FitRes_t> FitResMap_t;
FitResMap_t fit_results;
multimap< TrackFitWeight, Rset_t > fit_chi2;
for( list< pair<Double_t,Rvec_t> >::iterator it = road_combos.begin();
it != road_combos.end(); ++it ) {
fit_par.matchval = (*it).first;
Rvec_t& these_roads = (*it).second;
fit_par.roads = &these_roads;
fit_par.ndof = FitTrack( these_roads, fit_par.coef, fit_par.chi2 );
if( fit_par.ndof > 0 ) {
if( PassTrackCuts(fit_par) ) {
Rset_t road_tuple( ALL(these_roads) );
pair<FitResMap_t::iterator,bool> ins =
fit_results.insert( make_pair(road_tuple, fit_par) );
assert( ins.second );
fit_chi2.insert( make_pair( TrackFitWeight(fit_par), road_tuple) );
}
} else
FitErrPrint( fit_par.ndof );
}
#ifdef VERBOSE
if( fDebug > 2 ) {
cout << "Track candidates:" << endl;
for( multimap<TrackFitWeight,Rset_t>::iterator it = fit_chi2.begin();
it != fit_chi2.end(); ++it ) {
FitResMap_t::iterator found = fit_results.find((*it).second);
FitRes_t& r = (*found).second;
cout << "ndof = " << r.ndof
<< " rchi2 = " << r.chi2/(double)r.ndof
<< " x/y = " << r.coef[0] << "/" << r.coef[2]
<< " mx/my = " << r.coef[1] << "/" << r.coef[3]
<< endl;
}
}
#endif
vector<Rset_t> best_roads;
OptimalN( unique_found, fit_chi2, best_roads,
AnySharedHits(), CheckTypes(found_types) );
for( vector<Rset_t>::iterator it = best_roads.begin(); it !=
best_roads.end(); ++it ) {
FitResMap_t::iterator found = fit_results.find(*it);
assert( found != fit_results.end() );
NewTrack( tracks, (*found).second );
}
}
#ifdef TESTCODE
fN3dFits = nfits;
t_3dfit = 1e6*timer.RealTime();
#endif
#ifdef VERBOSE
if( fDebug > 0 ) {
Int_t ntr = tracks.GetLast()+1;
cout << ntr << " track";
if( ntr != 1 ) cout << "s";
if( nfits > 1 ) cout << " after chi2 optimization";
if( fDebug > 1 and ntr > 0 ) {
cout << ":" << endl;
for( Int_t i = 0; i < ntr; ++i ) {
THaTrack* tr = (THaTrack*)tracks.UncheckedAt(i);
cout << "3D track: x/y = " << tr->GetX() << "/" << tr->GetY()
<< " mx/my = " << tr->GetTheta() << "/" << tr->GetPhi()
<< " ndof = " << tr->GetNDoF()
<< " rchi2 = " << tr->GetChi2()/(double)tr->GetNDoF()
<< endl<<endl<<endl;
}
} else
cout << endl << endl;
}
#endif
}
else if( nproj == 2 ) {
}
#ifdef TESTCODE
t_coarse = 1e6*timer_tot.RealTime();
#endif
return 0;
}
Int_t MWDC::FineTrack( TClonesArray& tracks )
{
if( !TestBit(kDoFine) )
return 0;
return 0;;
}
Int_t MWDC::DefineVariables( EMode mode )
{
if( mode == kDefine and fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
RVarDef vars[] = {
{ "fail.nhits", "Too many hits in wire plane(s)", "fFailNhits" },
{ "fail.npat", "Too many treesearch patterns", "fFailNpat" },
#ifdef TESTCODE
{ "ncombos", "Number of road combinations", "fNcombos" },
{ "nfits", "Number of 3D track fits done", "fN3dFits" },
{ "t_track", "Time in 1st stage tracking (us)", "t_track" },
{ "t_3dmatch", "Time in MatchRoads (us)", "t_3dmatch" },
{ "t_3dfit", "Time fitting/selecting 3D tracks (us)", "t_3dfit" },
{ "t_coarse", "Total time in CoarseTrack (us)", "t_coarse" },
#endif
{ 0 }
};
DefineVarsFromList( vars, mode );
return 0;
}
THaAnalysisObject::EStatus MWDC::Init( const TDatime& date )
{
static const char* const here = "Init";
fRefMap->Reset();
fRefTime.clear();
fCrateMap = new THashTable(100);
fCrateMap->SetOwner();
EStatus status = THaTrackingDetector::Init(date);
if( !status ) {
if (fuseshower==1){
fshower->Init(date);
}
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
status = fPlanes[iplane]->Init(date);
if( status )
break;
}
}
delete fCrateMap; fCrateMap = 0;
if( status )
return fStatus = status;
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
THaDetMap* planeMap = fPlanes[iplane]->GetDetMap();
for( Int_t imod = 0; imod < planeMap->GetSize(); ++imod ) {
THaDetMap::Module* d = planeMap->GetModule(imod);
if( d->refchan < 0 ) continue;
THaDetMap::Module* refmod = fRefMap->Find( d->crate, d->slot,
d->refchan );
if( refmod ) {
d->refindex = refmod->first;
} else {
d->refindex = fRefMap->GetSize();
fRefMap->AddModule( d->crate, d->slot, d->refchan,
d->refchan, d->refindex, d->model );
refmod = fRefMap->GetModule(d->refindex);
if( !refmod ) {
Error( Here(here), "Error when adding reference channel cr/sl/ch="
"%u/%u/%u. Call expert.", d->crate, d->slot, d->refchan );
return fStatus = kInitError;
}
refmod->SetResolution( d->resolution );
refmod->MakeTDC();
}
}
}
for( Int_t imod = 0; imod < fRefMap->GetSize(); ++imod ) {
THaDetMap::Module* r = fRefMap->GetModule(imod);
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
THaDetMap* planeMap = fPlanes[iplane]->GetDetMap();
if( planeMap->Find( r->crate, r->slot, r->lo ) ) {
Error( Here(here), "Reference channel cr/sl/ch=%u/%u/%u also defined "
"as data channel in plane \"%s\". Fix database.",
r->crate, r->slot, r->lo, fPlanes[iplane]->GetName() );
return fStatus = kInitError;
}
}
}
fRefTime.assign( fRefMap->GetSize(), kBig );
sort( ALL(fPlanes), WirePlane::ZIsLess() );
if( !TestBit(kNoPartner) ) {
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
WirePlane* thePlane = fPlanes[iplane];
TString name( thePlane->GetName() );
if( name.EndsWith("p") ) {
TString other = name.Chop();
if( other.IsNull() )
continue;
vwiter_t it = find_if( ALL(fPlanes),WirePlane::NameEquals( other ) );
if( it != fPlanes.end() ) {
WirePlane* partner = *it;
if( thePlane->GetType() != partner->GetType() ) {
Error( Here(here), "Partner planes %s and %s have different types!"
" Fix database.", thePlane->GetName(), partner->GetName() );
return fStatus = kInitError;
}
if( fDebug > 0 )
Info( Here(here), "Partnering plane %s with %s",
thePlane->GetName(), partner->GetName() );
partner->SetPartner( thePlane );
}
}
}
}
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
status = fProj[type]->Init(date);
if( status )
return fStatus = status;
}
Double_t u_angle = fProj[kUPlane]->GetAngle()*TMath::RadToDeg();
Double_t v_angle = fProj[kVPlane]->GetAngle()*TMath::RadToDeg();
Int_t qu = TMath::FloorNint( u_angle/90.0 );
Int_t qv = TMath::FloorNint( v_angle/90.0 );
if( qu&1 == qv&1 ) {
Error( Here(here), "Plane misconfiguration: uangle (%6.2lf) and vangle "
"(%6.2lf) are in equivalent quadrants. Fix database.",
u_angle, v_angle );
return fStatus = kInitError;
}
Double_t du = u_angle - 90.0*qu;
Double_t dv = v_angle - 90.0*qv;
if( TMath::Abs(TMath::Abs(du)-45.0) > 44.0 or
TMath::Abs(TMath::Abs(dv)-45.0) > 44.0 ) {
Error( Here(here), "uangle (%6.2lf) or vangle (%6.2lf) too close "
"to 0 or 90 degrees. Fix database.", u_angle, v_angle );
return fStatus = kInitError;
}
ResetBit(k3dFastMatch);
if( fProj.size() == 3 ) {
assert( kUPlane < 2 && kVPlane < 2 );
Double_t uang =
TMath::Abs( TVector2::Phi_mpi_pi(fProj[kUPlane]->GetAngle()) );
if( (TMath::Abs( TVector2::Phi_mpi_pi(fProj[kVPlane]->GetAngle()) )-uang )
< 0.5*TMath::DegToRad() ) {
SetBit(k3dFastMatch);
Double_t tan = TMath::Tan( 0.5*TMath::Pi()-uang );
f3dMatchvalScalefact = 2.0 * (1.0/3.0 + tan*tan );
f3dMatchCut /= f3dMatchvalScalefact;
}
}
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
Projection* theProj = fProj[type];
Double_t width = 0.0;
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
WirePlane* thePlane = fPlanes[iplane];
if( thePlane->GetType() == type ) {
if( !thePlane->GetProjection() ) {
WirePlane* partner = thePlane->GetPartner();
assert( !partner or partner->GetProjection() == 0 );
theProj->AddPlane( thePlane, partner );
}
Double_t s = thePlane->GetWireStart();
Double_t d = thePlane->GetWireSpacing();
Double_t n = static_cast<Double_t>( thePlane->GetNelem() );
Double_t lo = s - 0.5*d;
Double_t hi = s + (n-0.5)*d;
Double_t w = max( TMath::Abs(hi), TMath::Abs(lo) );
if( w > width )
width = w;
}
}
UInt_t n = theProj->GetNplanes();
if( n < 3 ) {
Error( Here(here), "Too few planes of type \"%s\" defined. "
"Need >= 3, have %u. Fix database.", theProj->GetName(), n );
return fStatus = kInitError;
}
width *= 2.0;
if( width > 0.01 )
theProj->SetWidth( width );
else {
Error( Here(here), "Error calculating width of projection \"%s\". "
"Wire spacing too small? Fix database.", theProj->GetName() );
return fStatus = kInitError;
}
Double_t dz = TMath::Abs(theProj->GetZsize());
if( dz > 0.01 ) {
Double_t maxslope = width/dz;
if( theProj->GetMaxSlope() < 0.01 ) {
theProj->SetMaxSlope( maxslope );
} else if( theProj->GetMaxSlope() > maxslope ) {
if( fDebug > 0 ) {
Warning( Here(here), "For projection \"%s\", maxslope from "
"database = %lf exceeds geometric maximum = %lf. "
"Using smaller value.",
theProj->GetName(), theProj->GetMaxSlope(), maxslope );
}
theProj->SetMaxSlope( maxslope );
}
} else {
Error( Here(here), "Error calculating geometric maxslope for plane "
"type \"%s\". z-range of planes too small. Fix database.",
theProj->GetName() );
return fStatus = kInitError;
}
status = fProj[type]->InitLevel2(date);
if( status )
return fStatus = status;
}
if( fMaxThreads > 1 ) {
gSystem->Load("libThread");
delete fThreads;
fThreads = new ThreadCtrl;
fThreads->fTrack.reserve( fProj.size() );
fThreads->fTrackStartM->Lock();
for( vpsiz_t k = 0; k < fProj.size(); ++k ) {
TThread* t = fThreads->AddTrackThread( fProj[k] );
t->Run();
}
fThreads->fTrackStartM->UnLock();
}
return fStatus = kOK;
}
Int_t MWDC::ReadDatabase( const TDatime& date )
{
static const char* const here = "ReadDatabase";
fIsInit = kFALSE;
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane ) {
delete fPlanes[iplane];
}
fPlanes.clear();
fCalibPlanes.clear();
FILE* file = OpenFile( date );
if( !file ) return kFileError;
Int_t err = ReadGeometry( file, date );
if( err )
return err;
vector<vector<Int_t> > *cmap = new vector<vector<Int_t> >;
string planeconfig, calibconfig,specialplaneconfig;
Int_t time_cut = 1, pairs_only = 0, mc_data = 0, nopartner = 0;
f3dMatchCut = 1e-4;
Int_t event_display = 0, disable_tracking = 0, disable_finetrack = 0;
Int_t maxmiss = -1, maxthreads = -1;
Double_t conf_level = 1e-9;
DBRequest request[] = {
{ "planeconfig", &planeconfig, kString },
{ "specialplaneconfig",&specialplaneconfig,kString },
{ "usespecialplane", &fuseshower, kInt, 0, 1 },
{ "cratemap", cmap, kIntM, 6 },
{ "timecut", &time_cut, kInt, 0, 1 },
{ "pairsonly", &pairs_only, kInt, 0, 1 },
{ "MCdata", &mc_data, kInt, 0, 1 },
{ "nopartner", &nopartner, kInt, 0, 1 },
{ "3d_matchcut", &f3dMatchCut, kDouble, 0, 1 },
{ "event_display", &event_display, kInt, 0, 1 },
{ "disable_tracking", &disable_tracking, kInt, 0, 1 },
{ "disable_finetrack", &disable_finetrack, kInt, 0, 1 },
{ "calibrate", &calibconfig, kString, 0, 1 },
{ "3d_maxmiss", &maxmiss, kInt, 0, 1 },
{ "3d_chi2_conflevel", &conf_level, kDouble, 0, 1 },
{ "maxthreads", &maxthreads, kInt, 0, 1 },
{ 0 }
};
Int_t status = kInitError;
err = LoadDB( file, date, request, fPrefix );
fclose(file);
if( !err ) {
if( cmap->empty() ) {
Error(Here(here), "No cratemap defined. Set \"cratemap\" in database.");
} else {
for( vviter_t it = cmap->begin(); it != cmap->end(); ++it ) {
vector<int>& row = *it;
for( Int_t slot = row[1]; slot <= row[2]; ++slot ) {
DAQmodule* m =
new DAQmodule( row[0], slot, row[3], row[4]*1e-12, row[5] );
DAQmodule* found = static_cast<DAQmodule*>(fCrateMap->FindObject(m));
if( found ) {
m->Copy(*found);
delete m;
}
else
fCrateMap->Add(m);
}
}
status = kOK;
}
}
delete cmap; cmap = 0;
if( status != kOK )
return status;
SetBit( kDoTimeCut, time_cut );
SetBit( kPairsOnly, pairs_only );
SetBit( kMCdata, mc_data );
SetBit( kNoPartner, nopartner );
SetBit( kEventDisplay, event_display );
SetBit( kDoCoarse, !disable_tracking );
SetBit( kDoFine, !(disable_tracking or disable_finetrack) );
vector<string> planes = vsplit(planeconfig);
if( planes.empty() ) {
Error( Here(here), "No planes defined. Set \"planeconfig\" in database." );
return kInitError;
}
vector<string> specialplanes = vsplit(specialplaneconfig);
if( planes.empty() ) {
Error( Here(here), "No special planes defined. Set \"specialplaneconfig\" in database." );
return kInitError;
}
vector<string> calibplanes = vsplit(calibconfig);
for( ssiz_t i=0; i<planes.size(); ++i ) {
assert( !planes[i].empty() );
const char* name = planes[i].c_str();
vwiter_t it = find_if( ALL(fPlanes), WirePlane::NameEquals( name ) );
if( it != fPlanes.end() ) {
Error( Here(here), "Duplicate plane name: %s. Fix database.", name );
return kInitError;
}
WirePlane* newplane = new WirePlane( name, name, this );
if( !newplane or newplane->IsZombie() ) {
Error( Here(here), "Error creating wire plane %s. Call expert.", name );
return kInitError;
}
fPlanes.push_back( newplane );
newplane->SetDebug( fDebug );
if( !calibplanes.empty() ) {
vector<string>::iterator its =
find_if( ALL(calibplanes), bind2nd(equal_to<string>(), planes[i]) );
if( its != calibplanes.end() ) {
Info( Here(here), "Plane %s in calibration mode", name );
newplane->EnableCalibration();
fCalibPlanes.push_back( newplane );
calibplanes.erase( its );
}
}
}
if (fuseshower==1){
cout << "Use special planes! (Shower detector)" << endl;
assert( !specialplanes[0].empty() );
fshower = new THaBBShower(specialplanes[0].c_str(),"BigBite Shower");
string qxtemp;
qxtemp = "x" + specialplanes[0];
const char* name = qxtemp.c_str();
vwiter_t it = find_if( ALL(fPlanes), WirePlane::NameEquals( name ) );
if( it != fPlanes.end() ) {
Error( Here(here), "Duplicate plane name: %s. Fix database.", name );
return kInitError;
}
SpecialWirePlane* newplane = new SpecialWirePlane( name, name, this );
newplane->SetSpecial(1);
if( !newplane or newplane->IsZombie() ) {
Error( Here(here), "Error creating wire plane %s. Call expert.", name );
return kInitError;
}
fPlanes.push_back( newplane );
newplane->SetDebug( fDebug );
}else{
cout << "Do not use special planes! (Shower detector)" << endl;
}
if( !calibplanes.empty() ) {
string s("plane");
if( calibplanes.size() > 1 )
s.append("s");
for( vector<string>::size_type i = 0; i < calibplanes.size(); ++i ) {
s.append(" ");
s.append(calibplanes[i]);
}
Warning( Here(here), "Requested calibration for undefined %s. "
"Typo in database?", s.c_str() );
}
UInt_t nplanes = fPlanes.size();
if( fDebug > 0 )
Info( Here(here), "Loaded %u planes", nplanes );
if( maxmiss >= 0 )
fMinNdof = nplanes - 4 - maxmiss;
else
fMinNdof = 1;
if( conf_level < 0.0 or conf_level > 1.0 ) {
Error( Here(here), "Illegal fit confidence level = %lf. "
"Must be 0-1. Fix database.", conf_level );
return kInitError;
}
fChisqLimits.clear();
fChisqLimits.resize( nplanes-3, make_pair<Double_t,Double_t>(0,0) );
for( vec_pdbl_t::size_type dof = fMinNdof; dof < fChisqLimits.size();
++dof ) {
fChisqLimits[dof].first = TMath::ChisquareQuantile( conf_level, dof );
fChisqLimits[dof].second = TMath::ChisquareQuantile( 1.0-conf_level, dof );
}
if( maxthreads > 0 ) {
fMaxThreads = max(maxthreads,1);
if( fMaxThreads > 20 ) {
fMaxThreads = 20;
Warning( Here(here), "Excessive value of maxthreads = %d, "
"limited to %u", maxthreads, fMaxThreads );
}
} else {
SysInfo_t sysifo;
gSystem->GetSysInfo( &sysifo );
if( sysifo.fCpus > 0 )
fMaxThreads = sysifo.fCpus;
else
fMaxThreads = 1;
}
fIsInit = kTRUE;
return kOK;
}
void MWDC::Print(const Option_t* opt) const
{
THaTrackingDetector::Print(opt);
Int_t verbose = 0;
if( opt ) {
TString opt_s(opt);
opt_s.ToLower();
verbose = opt_s.CountChar('v');
}
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane )
fPlanes[iplane]->Print(opt);
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
if( type == kTypeBegin or verbose > 0 )
cout << endl;
fProj[type]->Print(opt);
}
}
void MWDC::SetDebug( Int_t level )
{
THaTrackingDetector::SetDebug( level );
for( vwsiz_t iplane = 0; iplane < fPlanes.size(); ++iplane )
fPlanes[iplane]->SetDebug( level );
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type )
fProj[type]->SetDebug( level );
}
void MWDC::EnableEventDisplay( Bool_t b )
{
if( !fIsInit ) {
Error( Here("EnableEventDisplay"), "Cannot enable/disable event display "
"support after initialization." );
return;
}
SetBit( kEventDisplay, b );
}
EProjType MWDC::NameToType( const char* name )
{
if( name and *name ) {
TString s(name);
s.ToLower();
for( EProjType type = kTypeBegin; type < kTypeEnd; ++type ) {
if( !fProj[type] )
continue;
TString ps( fProj[type]->GetName() );
ps.ToLower();
if( s == ps )
return type;
}
}
return kUndefinedType;
}
inline
static DAQmodule* FindDAQmodule( UShort_t crate, UShort_t slot,
const THashTable* table )
{
if( !table ) return 0;
CSpair m( crate, slot );
return static_cast<DAQmodule*>( table->FindObject(&m) );
}
UInt_t MWDC::LoadDAQmodel( THaDetMap::Module* mod ) const
{
DAQmodule* found = FindDAQmodule( mod->crate, mod->slot, fCrateMap );
UInt_t num = found ? found->fModel : 0;
mod->SetModel( num );
return num;
}
Double_t MWDC::LoadDAQresolution( THaDetMap::Module* mod ) const
{
DAQmodule* found = FindDAQmodule( mod->crate, mod->slot, fCrateMap );
Double_t res = found ? found->fResolution : 0.0;
mod->SetResolution( res );
return res;
}
UInt_t MWDC::GetDAQnchan( THaDetMap::Module* mod ) const
{
DAQmodule* found = FindDAQmodule( mod->crate, mod->slot, fCrateMap );
return found ? found->fNchan : 0;
}
}
ClassImp(TreeSearch::MWDC)
Last update: Tue Jul 7 19:26:19 2009
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.