// $Id: BrLocalTrackDC.cxx,v 1.14 1998/07/20 17:45:06 videbaek Exp $
// $Log: BrLocalTrackDC.cxx,v $
// Revision 1.14 1998/07/20 17:45:06 videbaek
// changed 0 return to NULL
//
//
#include "BrLocalTrackDC.h"
#include "BrEventNode.h"
#include "BrDetectorParamsDC.h"
#include "BrDetectorVolume.h"
#include "BrDigHit.h"
#include "BrDCHit.h"
#include "BrDCHitPosition.h"
#include "TSonataMath.h"
#include "BrTrackHit2D.h"
#include "btwn.h"
ClassImp(BrLocalTrackDC)
////////////////////////////////////////////////////////////////
//
// The TPhDetector class provides an interface to the BRAHMS
//
////////////////////////////////////////////////////////////////
BrLocalTrackDC::BrLocalTrackDC()
:BrLocalTrackModule()
{
}
BrLocalTrackDC::BrLocalTrackDC(Char_t *Name,Char_t *Title)
:BrLocalTrackModule(Name,Title)
{
fParams_p = NULL;
for(Int_t i=0;i<4;i++) fVolumeParams_p[i] = NULL;
fNumVol = 0;
CreateDCClonesArrays();
//defaults
fMainView = 0;
fSubView = -90;
}
void BrLocalTrackDC::CreateDCClonesArrays()
{
DCHits = new BrClonesArray("BrDCHit",1000);
TrackHits2D = new BrClonesArray("BrTrackHit2D",1000);
}
void BrLocalTrackDC::SetDetectorParamsDC(BrDetectorParamsDC* par)
{
// Set the TPC parameters with the values of the external parameters
// referenced by pointer
// Note that no owner ship is transferred. A copy is made.
//
//
if( fParams_p != 0)
delete fParams_p;
fParams_p = new BrDetectorParamsDC("Dummy","DCINT");
*fParams_p = *par;
// The conversion cast is needed because fParams_p is of
// type BrDetectorParams.
}
void BrLocalTrackDC::SetDetectorParamsDC(const BrDetectorParamsDC& par)
{
// Set the TPC parameters with the values of the external parameters.
// Note that no owner ship is transferred. A copy is made.
// There is no real need for the copy code here. Both Name and Title
// will be copied (it seems)
if( fParams_p != 0)
delete fParams_p;
fParams_p = new BrDetectorParamsDC("Dummy","ParamsTPC");
*fParams_p = par;
// The conversion cast is needed because fParams_p is of
// type BrDetectorParams.
}
void BrLocalTrackDC::SetDetectorVolume(BrDetectorVolume *vol) {
Int_t i;
for(i=0;i<fNumVol;i++) {
if(!strcmp(fVolumeParams_p[i]->GetName(),vol->GetName())) {
fVolumeParams_p[i] = vol;
return;
}
}
if(fNumVol<4) {
fVolumeParams_p[fNumVol] = vol;
fNumVol++;
}
else cout<<"Too many volume parameters, this one not added"<<endl;
}
BrDetectorVolume *BrLocalTrackDC::GetDetectorVolume(const Char_t *name) {
Int_t i;
for(i=0;i<fNumVol;i++) {
if(!strcmp(fVolumeParams_p[i]->GetName(),name)) {
return fVolumeParams_p[i];
}
}
cout<<"Detector Volume "<<name<<" does not exist"<<endl;
return NULL;
}
void BrLocalTrackDC::ListDetectorParameters(){
// List the current value of the digitization parameters on
// standard out.
if( fParams_p != NULL)
fParams_p->ListParameters();
else
cout << "No parameters set for this Digitisation module" << endl;
}
void BrLocalTrackDC::ListDetectorVolume() {
cout<<"There are "<<fNumVol<<" volumes defined"<<endl;
for(Int_t i=0;i<fNumVol;i++) {
cout<<"Volume "<<i<<" has name "<<fVolumeParams_p[i]->GetName()<<endl;
}
}
BrLocalTrackDC::~BrLocalTrackDC(){
//
// Is there a problem with the fact the the fParams_p can point to different kinds
// of datastructures i.e. potential memory leak. probably yes since fParams_p is
// defined as BrDetectorParams nor BrDetectorParamsTPC.
// By second thought it is probably ok. As long as all parameters are defined/set
// Using the SetDetectorParamsTPC / DC /BB etc..
//
if( fParams_p != 0) delete fParams_p;
for(Int_t i=0;i<4;i++) {
if( fVolumeParams_p[i] != 0) delete fVolumeParams_p[i];
}
}
void BrLocalTrackDC::Init(){
// The init member function/method is used to define histograms
// The histograms is added to the list of hist objects in this
// module.
// d_Hhits = new TH1F("Hhits","Geant Hits per event",40, (float)0., (float)40.);
// HistogramList()->Add(d_Hhits);
}
void BrLocalTrackDC::Event(BrEventNode* InputTable,BrEventNode* OutputTable)
{
Int_t status;
const Float_t clustdx = (float)0.1;
const Float_t clustdy = (float)0.1;
BrDataTable* DigitizedHits;
if(DebugLevel()>0) cout<<"Entering local tracking for "<<GetName()<<endl;
//Get stuff from the input table
Char_t TableName[32];
sprintf(TableName,"%s DigitizedDCx0",GetName());
DigitizedHits = InputTable->GetDataTable(TableName);
if(DigitizedHits == 0) {
if(DebugLevel() > 0 ) cout<<"No DigitizedDC in this event for "<<GetName()<<endl;
return;
}
CombinedHits->Clear();
LocalTracks->Clear();
TrackHits->Clear();
VirtualTracks->Clear();
DCHits->Clear();
TrackHits2D->Clear();
DCCluster(DigitizedHits,clustdx,clustdy);
status = TrackLocalElement();
if(!status) FillLocalTracks(OutputTable);
}
void BrLocalTrackDC::DCCluster(BrDataTable* DigitizedHits,Float_t dx,Float_t dy)
{
Int_t i;
Int_t numplanes;
#define MaxNumDCPLs 20
Float_t wirang[MaxNumDCPLs];
Int_t npair,ntrio,nsingle,iview,ipl,inview;
Int_t idet;
#define maxpl 20
#define maxpair 10
Int_t slist[maxpl];
Int_t plist[maxpair];
Int_t tlist[maxpair];
Int_t isview[maxpl];
Int_t ipview[maxpair];
Int_t itview[maxpair];
//STUFF ADDED FOR NEW THING
BrDetectorParamsDC *par;
BrDetectorVolume *vol[4];
DCHits->Clear();
par = GetDetectorParamsDC();
for(i=0;i<3;i++) {
Char_t volname[32];
sprintf(volname,"%sA%d",GetName(),i+1);
vol[i] = GetDetectorVolume(volname); //sub-modules
}
vol[3] = GetDetectorVolume(GetName()); //main module
for(idet=0;idet<3;idet++) {
numplanes = par->GetNplane();
for(i=0;i<numplanes;i++) wirang[i] = par->GetPlaneAt(i)->GetWirang();
npair = 0;
ntrio = 0;
nsingle = 0;
iview = 1;
ipl = 0;
while(ipl < numplanes ) {
inview = 1;
while((wirang[ipl+inview] == wirang[ipl]) && (ipl+inview<numplanes)) {
inview++;
}
if(inview == 1) {
slist[nsingle] = ipl;
isview[nsingle] = iview;
nsingle++;
}
else if(inview == 2) {
plist[npair] = ipl;
ipview[npair] = iview;
npair++;
}
else if(inview == 3) {
tlist[ntrio] = ipl;
itview[ntrio] = iview;
ntrio++;
}
else {
cout<<"DC_clust> *****************************"<<endl;
cout<<"DC_clust> ERROR ERROR ERROR ERROR ERROR"<<endl;
cout<<"DC_clust> I do not know how to handle "<<endl;
cout<<"DC_clust> this plane configuration !!! "<<endl;
}
ipl += inview;
iview++;
}
if(DebugLevel() > 3) {
cout<<"Single Planes"<<endl;
for(i=0;i<nsingle;i++) cout<<slist[i]<<" "<<isview[i]<<endl;
cout<<"Double planes"<<endl;
for(i=0;i<npair;i++) cout<<plist[i]<<" "<<ipview[i]<<endl;
}
for(i=0;i<npair;i++) DCPair(DigitizedHits,idet,plist[i],ipview[i],dx,dy);
for(i=0;i<nsingle;i++) DCSingle(DigitizedHits,idet,slist[i],isview[i]);
for(i=0;i<ntrio;i++) DCTrio(DigitizedHits,idet,tlist[i],itview[i],dx,dy);
} //loop over idet
}
void BrLocalTrackDC::DCSingle(BrDataTable *DigitizedHits,Int_t idet,Int_t ipl,Int_t iview)
{
Int_t i;
BrDCPlane *dcpl;
Float_t detZ,z1;
Int_t ih1,ip1;
Float_t hitcmb_pos[3];
Float_t hitcmb_dpos[2];
Int_t NumDigHitMod;
BrDigHit *dighit_p;
BrDCHitPosition *hitpos1;
BrCombinedHit *hitcmb_p;
BrDCHit *dchits_p;
//STUFF ADDED FOR NEW THING
BrDetectorParamsDC *par;
BrDetectorVolume *vol[4];
BrDataTable DigHitMod;
par = GetDetectorParamsDC();
for(i=0;i<3;i++) {
Char_t volname[32];
sprintf(volname,"%sA%d",GetName(),i+1);
vol[i] = GetDetectorVolume(volname); //sub-modules
}
vol[3] = GetDetectorVolume(GetName()); //main module
dcpl = par->GetPlaneAt(ipl);
detZ = vol[idet]->GetPos()[2]; //Z of center of detector
z1 = detZ + par->GetPldis()*(ipl+1 - (Float_t)(par->GetNplane()+1)/(Float_t)2.);
DigHitMod.Clear();
GetDigHitMod(idet,ipl,IANY,IANY,DigitizedHits,&DigHitMod);
NumDigHitMod = DigHitMod.GetEntries();
for(ih1=0;ih1<NumDigHitMod;ih1++) {
dighit_p = (BrDigHit*)DigHitMod.At(ih1);
if(dighit_p->GetUsed() == 0 ) {
for(ip1=0;ip1<dighit_p->GetHitposEntries();ip1++) {
BrClonesArray &hitcmb = *CombinedHits;
Int_t nhitcmb = CombinedHits->GetEntries();
hitcmb_p = new(hitcmb[nhitcmb]) BrCombinedHit();
hitpos1 = dighit_p->GetHitposAt(ip1);
hitcmb_pos[0] = hitpos1->GetUhit();
hitcmb_pos[1] = (Float_t)0.0;
hitcmb_pos[2] = z1;
hitcmb_p->SetPos(hitcmb_pos);
hitcmb_dpos[0] = par->GetPosres();
hitcmb_dpos[1] = (Float_t)0.0;
hitcmb_p->SetDpos(hitcmb_dpos);
hitcmb_p->SetImod((idet + 1)*100 + iview);
hitcmb_p->SetNhit(1);
hitcmb_p->SetStat(1);
hitcmb_p->SetType(1);
hitcmb_p->SetPlace(1);
BrClonesArray &dchits = *DCHits;
Int_t ndchits = DCHits->GetEntries();
dchits_p = new(dchits[ndchits]) BrDCHit();
dchits_p->SetImod(ipl); //should we add 1???
dchits_p->SetHitpos(hitpos1);
dchits_p->SetHitcmb(hitcmb_p);
dighit_p->IncUsed();
}
}
else {
// cout<<"Used .ne. 0"<<endl;
}
}
//Clear so that entries will not be deleted when we go out of scope
DigHitMod.Clear();
}
void BrLocalTrackDC::DCPair(BrDataTable* DigitizedHits,Int_t idet,Int_t ipl,Int_t iview,Float_t ddx,Float_t ddy)
{
Int_t i;
BrDCPlane *dcpl;
BrDCPlane *dcpl2;
BrDetectorParamsDC *par;
BrDetectorVolume *vol[4];
static const Double_t rad = .01745329;
Float_t detZ,z1,z2,zc,ddu;
Int_t ih1,ih2,ip1,ip2,staggered,wir_sgn,wir2_sgn;
Float_t hitcmb_pos[3];
Float_t hitcmb_dpos[2];
Int_t NumDigHitMod;
BrDataTable DigHitMod;
Int_t NumDigHitMod2;
BrDataTable DigHitMod2;
BrDCHitPosition *hitpos1;
BrDCHitPosition *hitpos2;
BrCombinedHit *hitcmb_p;
BrDCHit *dchits_p;
BrDigHit *dighit_p1;
BrDigHit *dighit_p2;
par = GetDetectorParamsDC();
dcpl = par->GetPlaneAt(ipl);
dcpl2 = par->GetPlaneAt(ipl+1);
for(i=0;i<3;i++) {
Char_t volname[32];
sprintf(volname,"%sA%d",GetName(),i+1);
vol[i] = GetDetectorVolume(volname); //sub-modules
}
vol[3] = GetDetectorVolume(GetName()); //main module
detZ = vol[idet]->GetPos()[2]; //center of detector
staggered = dcpl->GetNwire() != dcpl2->GetNwire();
ddu = ddx*(Float_t)TMath::Cos(dcpl->GetWirang()*rad) - ddy*(Float_t)TMath::Sin(dcpl->GetWirang()*rad);
z1 = detZ + par->GetPldis()*(ipl+1-(Float_t)(par->GetNplane()+1)/(Float_t)2.0);
z2 = z1 + par->GetPldis();
zc = (Float_t)0.5 * (z1 + z2);
wir_sgn = dcpl->GetNwire() * 2 - 1;
wir2_sgn = dcpl2->GetNwire() * 2 - 1;
DigHitMod.Clear();
GetDigHitMod(idet,ipl,IANY,IANY,DigitizedHits,&DigHitMod);
NumDigHitMod = DigHitMod.GetEntries();
for(ih1=0;ih1<NumDigHitMod;ih1++) {
dighit_p1 = (BrDigHit*)DigHitMod.At(ih1);
DigHitMod2.Clear();
GetDigHitMod(idet,ipl+1,dighit_p1->GetWire(),IANY,DigitizedHits,&DigHitMod2);
NumDigHitMod2 = DigHitMod2.GetEntries();
for(ih2=0;ih2<NumDigHitMod2;ih2++) {
dighit_p2 = (BrDigHit*)DigHitMod2.At(ih2);
for(ip1=0;ip1<dighit_p1->GetHitposEntries();ip1++) {
hitpos1 = dighit_p1->GetHitposAt(ip1);
for(ip2=0;ip2<dighit_p2->GetHitposEntries();ip2++) {
hitpos2 = dighit_p2->GetHitposAt(ip2);
if((staggered &&
((hitpos1->GetSign() == wir_sgn||hitpos2->GetSign() == wir2_sgn) &&
(hitpos1->GetSign()*hitpos2->GetSign() == -1))) ||
(!staggered &&
(hitpos1->GetSign() == hitpos2->GetSign())) ){
if(TMath::Abs(hitpos1->GetUhit()-hitpos2->GetUhit())<=ddu) {
BrClonesArray &hitcmb = *CombinedHits;
Int_t nhitcmb = CombinedHits->GetEntries();
hitcmb_p = new(hitcmb[nhitcmb]) BrCombinedHit();
hitcmb_pos[0] = (Float_t)0.5*(hitpos1->GetUhit()+hitpos2->GetUhit());
hitcmb_pos[1] = (Float_t)0.0;
hitcmb_pos[2] = zc;
hitcmb_p->SetPos(hitcmb_pos);
hitcmb_dpos[0] = par->GetPosres() / (Float_t)1.41;
hitcmb_dpos[1] = (Float_t)0.0;
hitcmb_p->SetDpos(hitcmb_dpos);
hitcmb_p->SetImod((idet + 1)*100 + iview);
hitcmb_p->SetNhit(2);
hitcmb_p->SetStat(1);
hitcmb_p->SetType(1); //DC
hitcmb_p->SetPlace(2);
BrClonesArray &dchits = *DCHits;
Int_t ndchits = DCHits->GetEntries();
dchits_p = new(dchits[ndchits]) BrDCHit();
dchits_p->SetImod(ipl); //should we add 1???
dchits_p->SetHitpos(hitpos1);
dchits_p->SetHitcmb(hitcmb_p);
BrClonesArray &dchits1 = *DCHits;
ndchits = DCHits->GetEntries();
dchits_p = new(dchits1[ndchits]) BrDCHit();
dchits_p->SetImod(ipl+1); //should we add 1???
dchits_p->SetHitpos(hitpos2);
dchits_p->SetHitcmb(hitcmb_p);
dighit_p1->IncUsed();
dighit_p2->IncUsed();
}
}
}
}
}
if(staggered) {
DigHitMod2.Clear();
GetDigHitMod(idet,ipl+1,dighit_p1->GetWire()-dcpl->GetNwire()+dcpl2->GetNwire(),IANY,DigitizedHits,&DigHitMod2);
NumDigHitMod2 = DigHitMod2.GetEntries();
for(ih2=0;ih2<NumDigHitMod2;ih2++) {
dighit_p2 = (BrDigHit*)DigHitMod2.At(ih2);
for(ip1=0;ip1<dighit_p1->GetHitposEntries();ip1++) {
hitpos1 = dighit_p1->GetHitposAt(ip1);
for(ip2=0;ip2<dighit_p2->GetHitposEntries();ip2++) {
hitpos2 = dighit_p2->GetHitposAt(ip2);
if(((hitpos1->GetSign() != wir_sgn ) &&
(hitpos2->GetSign() != wir2_sgn)) &&
(hitpos1->GetSign()*hitpos2->GetSign() == -1) ) {
//Make overlap hit if accepted
if(TMath::Abs(hitpos1->GetUhit()-hitpos2->GetUhit())<=ddu) {
BrClonesArray &hitcmb = *CombinedHits;
Int_t nhitcmb = CombinedHits->GetEntries();
hitcmb_p = new(hitcmb[nhitcmb]) BrCombinedHit();
hitcmb_pos[0] = (Float_t)0.5*(hitpos1->GetUhit()+hitpos2->GetUhit());
hitcmb_pos[1] = (Float_t)0.0;
hitcmb_pos[2] = zc;
hitcmb_p->SetPos(hitcmb_pos);
hitcmb_dpos[0] = par->GetPosres() / (Float_t)1.41;
hitcmb_dpos[1] = (Float_t)0.0;
hitcmb_p->SetDpos(hitcmb_dpos);
hitcmb_p->SetImod((idet + 1)*100 + iview);
hitcmb_p->SetNhit(2);
hitcmb_p->SetStat(1);
hitcmb_p->SetType(1); //DC
hitcmb_p->SetPlace(3);
BrClonesArray &dchits = *DCHits;
Int_t ndchits = DCHits->GetEntries();
dchits_p = new(dchits[ndchits]) BrDCHit();
dchits_p->SetImod(ipl); //should we add 1???
dchits_p->SetHitpos(hitpos1);
dchits_p->SetHitcmb(hitcmb_p);
BrClonesArray &dchits1 = *DCHits;
ndchits = DCHits->GetEntries();
dchits_p = new(dchits1[ndchits]) BrDCHit();
dchits_p->SetImod(ipl+1); //should we add 1???
dchits_p->SetHitpos(hitpos2);
dchits_p->SetHitcmb(hitcmb_p);
dighit_p1->IncUsed();
dighit_p2->IncUsed();
}
}
}
}
}
}
}
//Call the DC_single for the plane two of the pair set.
//DC_single checks on used hits.
//It must also call the first plane for unused hits - which were
//not paired up.
DCSingle(DigitizedHits,idet,ipl+1,iview);
DCSingle(DigitizedHits,idet,ipl ,iview);
//Clear so that we don't delete tables when we go out of scope
DigHitMod.Clear();
DigHitMod2.Clear();
}
void BrLocalTrackDC::DCTrio(BrDataTable* DigitizedHits,Int_t idet,Int_t ipl,Int_t iview,Float_t ddx,Float_t ddy)
{
}
void BrLocalTrackDC::GetDigHitMod(Int_t idet,Int_t imod,Int_t wire,Int_t time,BrDataTable *DigitizedHits,BrDataTable *DigHitMod)
{
Int_t i;
BrDigHit* dighit_p;
Int_t NumDigHits = DigitizedHits->GetEntries();
for(i=0;i<NumDigHits;i++) {
dighit_p = (BrDigHit*)DigitizedHits->At(i);
if(( dighit_p->GetIdet() == idet || idet == IANY ) &&
( dighit_p->GetImod() == imod || imod == IANY ) &&
( dighit_p->GetWire() == wire || wire == IANY ) &&
( dighit_p->GetTime() == time || time == IANY )) {
DigHitMod->Add(dighit_p);
}
}
DigHitMod->Sort();
}
Int_t BrLocalTrackDC::TrackLocalElement()
{
static const Float_t dccut_matchdx = (Float_t)0.018;
static const Float_t dccut_matchdy = (Float_t)0.018;
static const Int_t dccut_miss[3] = {2,1,1};
static const Double_t rad = .01745329;
Float_t dddd[2];
Int_t maxmiss[3];
Int_t nmv,nsv,nov,nvv,nview,Imod;
Int_t i,j,k,i1,i2,it;
Int_t NoSaved,NoSavedPl,nmmiss,mview,nok;
Float_t u3,z3,du3,uc,duc;
static Float_t const asigma = (Float_t)3.0;
BrDCPlane *dcpl;
Int_t numplanes,idet;
Int_t vlistID[5][maxview];//????????
//Algorithm parameters
//Assume main view is x, submain view is y.
//Note: the x-view should have the most planes
//and there should be at least 3 sets of y-view planes
//static const Int_t main_view = 0;
//static const Int_t sub_view = -90;
Int_t main_view; //set to fMainView
Int_t sub_view; //goes to fSubView;
Int_t ncutoff,iv,numloctra,oview;
Float_t u1,z1,du1,u2,z2,du2,au,au1,au2,dau,th1,th2,th3;
Float_t uc1,uc2,x,y,dx2,dy2,arg1,arg2,arg3;
Float_t tra2d_pos[2];
Float_t loctra_pos[3];
Float_t loctra_vec[3];
Float_t dx = (Float_t)0.0;
Float_t dy = (Float_t)0.0;
Int_t vlist[maxview]; //list of angles of views
Int_t mvlist[maxview]; //list of main views
Int_t svlist[maxview]; //list of sub views
Int_t ovlist[maxview]; //list of other views
//static const Int_t *xvlist[2] = {mvlist,svlist}; //pointers to above lists
//static const Int_t *nxv[2] = {&nmv,&nsv}; //pointers to above indices
const Int_t *xvlist[2] = {mvlist,svlist}; //pointers to above lists
const Int_t *nxv[2] = {&nmv,&nsv}; //pointers to above indices
#define max_saved 1024
BrDataTable SavedId;
Int_t IDSave,IDCurrent,IDCurrentPl,IDCur;
BrCombinedHit* hitcmbb;
BrCombinedHit* hitcmb1;
BrCombinedHit* hitcmb2;
BrTrackHit2D* tra2d_p;
BrTrackHit2D* tra2d_p1;
BrLocalTrack* loctra_p;
BrLocalTrack* loctra_p1;
Int_t NumHitCMBMod[maxview];
BrDataTable HitCMBMod[maxview];
BrDataTable HitCMBSel;
Int_t NumTra2DMod1;
BrDataTable Tra2DMod1;
Int_t NumTra2DMod2;
BrDataTable Tra2DMod2;
Float_t sol[4];
Float_t covar[4][4];
Float_t chisq;
Int_t ifail;
Int_t nloctra,nloctraa,flg;
BrDetectorParamsDC *par;
main_view = GetMainView();
sub_view = GetSubView();
if( DebugLevel() > 4 ) {
// See what kind of hitcmb's we have.
Int_t NumHitCMB = CombinedHits->GetEntries();
for(i=0;i<NumHitCMB;i++) {
BrCombinedHit *hitcmb_p = (BrCombinedHit*)CombinedHits->At(i);
cout<<hitcmb_p;
}
}
par = GetDetectorParamsDC();
for(i=0;i<3;i++) maxmiss[i] = dccut_miss[i];
dddd[0] = (Float_t)dccut_matchdx;
dddd[1] = (Float_t)dccut_matchdy;
//Figure out the views
nvv = 0;
nmv = 0;
nsv = 0;
nov = 0;
nview = 0;
for(idet=0;idet<3;idet++) {
numplanes = par->GetNplane();
// role of nview needs to investigated carefully
for(j=0;j<numplanes;j++) {
dcpl = par->GetPlaneAt(j);
if( nview == 0 ) {
nvv = 0;
// TSonataMath::Int appears to behave as FORTRAN NINT
vlist[nview] = TSonataMath::Int(dcpl->GetWirang());
vlistID[idet][nview] = nview;
Imod = (idet + 1)*100 + nview + 1;
HitCMBMod[0].Clear();
GetHitCMBMod(Imod,&HitCMBMod[0]);
NumHitCMBMod[0] = HitCMBMod[0].GetEntries();
nview = 1;
}
else if((TSonataMath::Int(dcpl->GetWirang()) != vlist[nview-1]) || j == 0) {
// Come in here only if this wire has a different angle than the one
// in front of it. Also come in here if this is the first plane
// in a sub-module
if(j==0) {
nvv = 0;
}
else {
nvv++;
}
vlist[nview] = TSonataMath::Int(dcpl->GetWirang());
vlistID[idet][nvv] = nview;
if(nview>maxview) cout<<"nview too large, it is "<<nview<<endl;
Imod = (idet + 1)*100 + nvv + 1;
HitCMBMod[nview].Clear();
GetHitCMBMod(Imod,&HitCMBMod[nview]);
NumHitCMBMod[nview] = HitCMBMod[nview].GetEntries();
nview++;
}
}
}
//Set up list of main and sub views. Also keep track of how many there are.
for(i=0;i<nview;i++) {
if(vlist[i] == main_view) {
mvlist[nmv] = i;
nmv++;
}
else if(vlist[i] == sub_view) {
svlist[nsv] = i;
nsv++;
}
else {
ovlist[nov] = i;
nov++;
}
}
//Setup cutoff parameters based on actual plane configuration
ncutoff = nmv + nsv + nov;
ncutoff -= 5;
//Main and submain view 2d tracking
//=================================
//Define local 2d tracks by taking hitcmb for the first and last
//planes in main and subview
if(DebugLevel()>0) {
printf("Total number hitcmb is %dn",CombinedHits->GetEntries());
if(DebugLevel() > 2) {
Int_t numhitcmbtest = 0;
for(i=0;i<nview;i++) {
printf("NumHitCMB for i = %d is %dn",i,HitCMBMod[i].GetEntries());
numhitcmbtest += HitCMBMod[i].GetEntries();
}
printf("Total hitcmb calculated from parts is %dn",numhitcmbtest);
}
}
for(iv=0;iv<2;iv++) {
i1 = xvlist[iv][0]; //First plane
i2 = xvlist[iv][*nxv[iv]-1]; //Last plane
NumHitCMBMod[i1] = HitCMBMod[i1].GetEntries(); //Entries in first plane
NumHitCMBMod[i2] = HitCMBMod[i2].GetEntries();//Entries in last plane
if(DebugLevel() > 2) {
printf("i1,i2 = %d,%d NumHitCMBMod[i1,i2] = %d,%dn",i1,i2,NumHitCMBMod[i1],NumHitCMBMod[i2]);
}
for(i=0;i<NumHitCMBMod[i1];i++) {
hitcmb1 = (BrCombinedHit*)HitCMBMod[i1].At(i);
u1 = hitcmb1->GetPos()[0];
z1 = hitcmb1->GetPos()[2];
du1 = hitcmb1->GetDpos()[0];
i2 = xvlist[iv][*nxv[iv]-1]; //Last plane
NumHitCMBMod[i2] = HitCMBMod[i2].GetEntries();//Entries in last plane
for(j=0;j<NumHitCMBMod[i2];j++) {
hitcmb2 = (BrCombinedHit*)HitCMBMod[i2].At(j);
u2 = hitcmb2->GetPos()[0];
z2 = hitcmb2->GetPos()[2];
du2 = hitcmb2->GetDpos()[0];
au = (u2-u1)/(z2-z1);
dau = (Float_t)TMath::Sqrt(du1*du1 + du2*du2)/(z2-z1);
BrClonesArray &tra2d = *TrackHits2D;
Int_t ntra2d = TrackHits2D->GetEntries();
tra2d_p = new(tra2d[ntra2d]) BrTrackHit2D();
tra2d_pos[0] = hitcmb1->GetPos()[2];
tra2d_pos[1] = hitcmb1->GetPos()[0];
tra2d_p->SetPos(tra2d_pos);
tra2d_p->SetVec(au);
tra2d_p->SetDpos(du1);
tra2d_p->SetDvec(dau);
tra2d_p->SetNhit(2);
tra2d_p->SetFlg(iv);
tra2d_p->SetStat(0);
InsertTrahit(hitcmb1,tra2d_p);
InsertTrahit(hitcmb2,tra2d_p);
} //loop over j
} //loop over i
} //loop over iv
//Initial possible 2d tracks have now been identified
//Loop over intermediate positions in the local tracks and compare
//to the initially defined local tracks.
//This might create additional ghost 2d tracks
//First put into same order as ADAMO for easy comparison
//ChangeTra2DOrder();
Int_t NumTra2DD = TrackHits2D->GetEntries();
if(NumTra2DD > 300) {
if(DebugLevel()>0) {
cout<<"Number Tra2D too large, it is "<<NumTra2DD<<endl;
cout<<"Aborting tracking procedure for this part of event"<<endl;
}
goto Return;
}
for(i=0;i<NumTra2DD;i++) {
NoSaved = 1;
tra2d_p = (BrTrackHit2D*)TrackHits2D->At(i);
SavedId.Clear();
SavedId.Add(tra2d_p);
if(tra2d_p->GetStat() == 0) {
iv = tra2d_p->GetFlg(); //x view or y view
nmmiss = 0;
for(mview=1;mview<*nxv[iv]-1;mview++) {
nok = 0;
NoSavedPl = NoSaved;
i1 = xvlist[iv][mview];
for(k=0;k<NumHitCMBMod[i1];k++) {
hitcmbb = (BrCombinedHit*)HitCMBMod[i1].At(k);
u1 = tra2d_p->GetPos()[1];
z1 = tra2d_p->GetPos()[0];
u3 = hitcmbb->GetPos()[0];
z3 = hitcmbb->GetPos()[2];
au = tra2d_p->GetVec();
dau = tra2d_p->GetDvec();
du1 = tra2d_p->GetDpos();
du3 = hitcmbb->GetDpos()[0];
uc = u1 + au*(z3-z1);
duc = dddd[iv];
if(BTWN(uc-asigma*duc,u3,uc+asigma*duc)) {
nok++;
if(nok == 1) {
if(NoSaved>max_saved) printf("NoSaved = %dn",NoSaved);
for(it=0;it<NoSaved;it++) {
BrTrackHit2D* tra2d_pp = (BrTrackHit2D*)SavedId.At(it);
tra2d_pp->IncNhit();
InsertTrahit(hitcmbb,tra2d_pp);
}
}
else if(nok > 1) {
// More than one hit -- create a ghost local track
if(2*NoSaved > max_saved) {
cout<<"**** Error : cannot save GHOSTS **"<<endl;
}
else {
for(j=0;j<NoSavedPl;j++) {
BrClonesArray &tra2d = *TrackHits2D;
Int_t ntra2d = TrackHits2D->GetEntries();
tra2d_p1 = new(tra2d[ntra2d]) BrTrackHit2D(tra2d_p);
if(NoSaved+j>max_saved) printf("NoSaved1 = %dn",NoSaved+j);
SavedId.AddAt(tra2d_p1,NoSaved+j);
CopyTrahit1((BrTrackHit2D*)SavedId.At(j),tra2d_p1);
}
// Restore the current hitcmb picked at this plane
// before inserting into the trahit structure for
// this new ghost track(s)
for(j=0;j<NoSavedPl;j++) {
BrTrackHit2D* tra2d_pp = (BrTrackHit2D*)SavedId.At(NoSaved+j);
InsertTrahit(hitcmbb,tra2d_pp);
}
// Update total number of tracks which belong to
// this group of ghosts.
NoSaved += NoSavedPl;
}
}
}
} //loop over k
// No hits found for this view.
// Update the missed hit pointer
if(nok == 0) {
nmmiss++;
if(nmmiss>maxmiss[iv]) {
for(j=0;j<NoSaved;j++) RemoveTrack((BrTrackHit2D*)SavedId.At(j));
break;
}
}
}
if(DebugLevel() > 2) cout<<"nmiss for tra2d = "<<nmmiss<<endl;
}
}
//Clean up tra2d table
if(DebugLevel()>5) cout<<"Calling DeleteBadTrack()"<<endl;
DeleteBadTrack();
if(DebugLevel()>5) cout<<"Calling GetTra2dMod(...)"<<endl;
//Combine views
flg = 0;
GetTra2DMod(0,flg,&Tra2DMod1);
flg = 1;
GetTra2DMod(0,flg,&Tra2DMod2);
if(DebugLevel()>5) cout<<"Getting number of trahitmods"<<endl;
//printf("TrackLocalElement checkpoint 5 Mod1 = %d Mod2 = %dn",NumTra2DMod1,NumTra2DMod2);
NumTra2DMod1 = Tra2DMod1.GetEntries();
NumTra2DMod2 = Tra2DMod2.GetEntries();
if(DebugLevel()>1) {
cout<<"NumTra2DMod1,2 = "<<NumTra2DMod1<<","<<NumTra2DMod2<<endl;
}
for(it=0;it<NumTra2DMod1;it++) {
BrTrackHit2D* tra2dd = (BrTrackHit2D*)Tra2DMod1.At(it);
u1 = tra2dd->GetPos()[1];
du1 = tra2dd->GetDpos();
z1 = tra2dd->GetPos()[0];
au1 = tra2dd->GetVec();
th1 = (Float_t)main_view;
for(j=0;j<NumTra2DMod2;j++) {
BrTrackHit2D* tra2dd2 = (BrTrackHit2D*)Tra2DMod2.At(j);
u2 = tra2dd2->GetPos()[1];
du2 = tra2dd2->GetDpos();
z2 = tra2dd2->GetPos()[0];
au2 = tra2dd2->GetVec();
th2 = (Float_t)sub_view;
BrClonesArray &loctra = *LocalTracks;
nloctra = LocalTracks->GetEntries();
loctra_p = new(loctra[nloctra]) BrLocalTrack();
loctra_pos[0] = u1;
loctra_pos[1] = u2 + au2*(z1-z2);
loctra_pos[2] = z1;
loctra_p->SetPos(loctra_pos);
loctra_vec[0] = au1;
loctra_vec[1] = au2;
loctra_vec[2] = (Float_t)1.0;
loctra_p->SetVec(loctra_vec);
loctra_p->SetNhit(tra2dd->GetNhit() + tra2dd2->GetNhit());
loctra_p->SetFlg(0);
loctra_p->SetTr1(tra2dd);
loctra_p->SetTr2(tra2dd2);
numloctra = LocalTracks->GetEntries()-1;
IDSave = numloctra;
IDCurrent = numloctra;
nmmiss = 0;
// Loop over intermediate planes u, v for verification
for(oview=0;oview<nov;oview++) {
IDCurrentPl = IDCurrent;
i1 = ovlist[oview];
nok = 0;
for(k=0;k<NumHitCMBMod[i1];k++) {
hitcmbb = (BrCombinedHit*)HitCMBMod[i1].At(k);
u3 = hitcmbb->GetPos()[0];
z3 = hitcmbb->GetPos()[2];
du3 = hitcmbb->GetDpos()[0];
th3 = (Float_t)vlist[i1];
uc1 = u1 + au1*(z3-z1);
uc2 = u2 + au2*(z3-z2);
x = (uc2*(Float_t)TMath::Sin(th1*rad)-uc1*(Float_t)TMath::Sin(th2*rad))/(Float_t)TMath::Sin((th1-th2)*rad);
y = (uc2*(Float_t)TMath::Cos(th1*rad)-uc1*(Float_t)TMath::Cos(th2*rad))/(Float_t)TMath::Sin((th1-th2)*rad);
dx2 = dx * dx;
dy2 = dy * dy;
uc = x*(Float_t)TMath::Cos(th3*rad) - y*(Float_t)TMath::Sin(th3*rad);
arg1 = dx2*(Float_t)TMath::Cos(th3*rad);
arg2 = dy2*(Float_t)TMath::Sin(th3*rad);
arg3 = du3;
duc = (Float_t)sqrt(arg1*arg1 + arg2*arg2 + arg3*arg3);
if(DebugLevel()>4) {
cout<<"uc,duc,uc-sigma*duc,u3,uc+sigma*duc ";
cout<<uc<<","<<duc<<",";
cout<<uc-asigma*duc<<","<<u3<<","<<uc+asigma*duc<<endl;
}
if(BTWN(uc-asigma*duc,u3,uc+asigma*duc)) {
nok++;
if(nok == 1) {
// Not the right place to update the #hits on this
// track??
//
// The loctra_ID need not be the same throughout the
// whole loop.
// Update the #hits information and save the current
// hits for all the local tracks
for(i=IDSave;i<=IDCurrent;i++) {
BrLocalTrack *loctra_p2 = (BrLocalTrack*)LocalTracks->At(i);
loctra_p2->IncNhit();
BrLocalTrackModule::InsertTrahit(hitcmbb,loctra_p2);
}
}
else if(nok > 1) {
// Add ghost tracks
for(i=IDSave;i<=IDCurrentPl;i++) {
BrClonesArray &loctra = *LocalTracks;
nloctra = LocalTracks->GetEntries();
loctra_p1 = new(loctra[nloctra]) BrLocalTrack(loctra_p);
BrLocalTrack* loctra_pi = (BrLocalTrack*)LocalTracks->At(i);
CopyTrahit1(loctra_pi,loctra_p1);//???
}
IDCur = IDCurrent;
// IDCurrent = numloctra;
IDCurrent = LocalTracks->GetEntries()-1;
for(i=IDCur+1;i<=IDCurrent;i++) {
BrLocalTrack *loctra_p2 = (BrLocalTrack*)LocalTracks->At(i);
BrLocalTrackModule::InsertTrahit(hitcmbb,loctra_p2);
}
}
}
} //loop over k
if(nok == 0) {
nmmiss++;
if(nmmiss>maxmiss[2]) {
for(i=IDSave;i<=IDCurrent;i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
BrLocalTrackModule::RemoveTrack(loctra_p);
}
break;
}
}
} //loop over oview
if(DebugLevel() > 2) cout<<"nmiss for loctra = "<<nmmiss<<endl;
} //loop over j (y view tracks)
} //loop over i (x view tracks)
//set track status to be 1
if(DebugLevel()>1) cout<<"NumLoctra before cleaning is "<<LocalTracks->GetEntries()<<endl;
for(i=0;i<LocalTracks->GetEntries();i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
if(loctra_p->GetStat() == 0) loctra_p->SetStat(1);
}
nloctraa = LocalTracks->GetEntries();
BrLocalTrackModule::DeleteBadTrack(); //want to delete bad loctra, so use one in TDetector
if(DebugLevel()>1) cout<<"NumLoctra after cleaning is "<<LocalTracks->GetEntries()<<endl;
FillVTracks(ncutoff);
//printf("TrackLocalElement checkpoint 7n");
//Fit the local tracks
nloctra = LocalTracks->GetEntries();
for(i=0;i<nloctra;i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
FillHitCMBSelector(loctra_p,&HitCMBSel);
FitLocalDCTrack(&HitCMBSel,loctra_p,vlist,vlistID,nview,sol,covar[0],&chisq,&ifail);
if(DebugLevel() > 0) {
printf("DC fit %f %f %f %f %fn",sol[0],sol[2],sol[1],sol[3],chisq);
}
loctra_pos[0] = sol[1];
loctra_pos[1] = sol[3];
loctra_pos[2] = (Float_t)0.0;
loctra_p->SetPos(loctra_pos);
loctra_vec[0] = sol[0];
loctra_vec[1] = sol[2];
loctra_vec[2] = (Float_t)1.0;
loctra_p->SetVec(loctra_vec);
loctra_p->SetChisq(chisq);
}
//Need to clear following so that entries are not deleted when
//going out of scope
Return:
SavedId.Clear();
for(i=0;i<maxview;i++) HitCMBMod[i].Clear();
HitCMBSel.Clear();
Tra2DMod1.Clear();
Tra2DMod2.Clear();
return(0);
}
void BrLocalTrackDC::FitLocalDCTrack(BrDataTable *HitCMBSel,BrLocalTrack* loctra_p,Int_t vlist[maxview],Int_t vlistID[5][maxview],Int_t nview,Float_t *sol,Float_t *covar,Float_t *chisq,Int_t *ifail)
{
Int_t i,j,ic;
static const Double_t rad = .017453292519943;
static const Double_t rad1 = 2.0*TMath::ASin(1.0)/180.0;
Double_t a[4];
Double_t beta[4][1];
Double_t alpha[4][4];
Double_t *alphaa[4],*betaa[4];
Float_t work[4];
Float_t du,chisqq,dchi;
Float_t csv[20],siv[20];
Int_t iv,ip;
BrCombinedHit* hitcmb_p;
Int_t NumHitCMBSel;
//void DEQINV(Int_t n,Double_t **a,Int_t idim,Float_t *r,Int_t *ifail,Int_t k,Double_t **b);
//Int_t MOD(Int_t divident,Int_t divisor);
for(i=0;i<nview;i++) {
Double_t angl = (Double_t)vlist[i] * rad;
csv[i] = (Float_t)TMath::Cos(angl);
siv[i] = (Float_t)TMath::Sin(angl);
}
for(i=0;i<4;i++) {
beta[i][0] = (Double_t)0.0;
for(j=0;j<4;j++) {
alpha[i][j] = (Double_t)0.0;
}
}
//Is there a more elegant way to do this?
for(i=0;i<4;i++) {
alphaa[i] = &alpha[i][0];
betaa [i] = &beta[i][0];
}
NumHitCMBSel = HitCMBSel->GetEntries();
//Setup matrix equation
if(NumHitCMBSel<=2) {
cout<<"Too few points for fitting; this ought not to happen"<<endl;
return;
}
for(ic=0;ic<NumHitCMBSel;ic++) {
hitcmb_p = (BrCombinedHit*)HitCMBSel->At(ic);
iv = TSonataMath::MOD(hitcmb_p->GetImod(),100)-1;
ip = hitcmb_p->GetImod() / 100-1;
iv = vlistID[ip][iv];
a[0] = csv[iv]*hitcmb_p->GetPos()[2];
a[1] = csv[iv];
a[2] = -siv[iv]*hitcmb_p->GetPos()[2];
a[3] = -siv[iv];
du = hitcmb_p->GetDpos()[0];
for(i=0;i<4;i++) {
for(j=0;j<4;j++) alpha[i][j] += a[i]*a[j]/du/du;
beta[i][0] += a[i] * hitcmb_p->GetPos()[0]/du/du;
}
}
//Solve 4x4 matrix for xc,yc and ex,ey vector
TSonataMath::DEQINV(4,alphaa,4,work,ifail,1,betaa);
//Calculate chisqr
chisqq = (Float_t)0.0;
for(ic=0;ic<NumHitCMBSel;ic++) {
hitcmb_p = (BrCombinedHit*)HitCMBSel->At(ic);
iv = TSonataMath::MOD(hitcmb_p->GetImod(),100)-1;
ip = hitcmb_p->GetImod() / 100-1;
iv = vlistID[ip][iv];
a[0] = csv[iv]*hitcmb_p->GetPos()[2];
a[1] = csv[iv];
a[2] = -siv[iv]*hitcmb_p->GetPos()[2];
a[3] = -siv[iv];
du = hitcmb_p->GetDpos()[0];
dchi = (Float_t)0.0;
for(i=0;i<4;i++) dchi += (Float_t)(a[i] * beta[i][0]);
dchi = (dchi - hitcmb_p->GetPos()[0])/du;
chisqq += dchi*dchi;
}
*chisq = chisqq / (Float_t)(NumHitCMBSel-4);
for(i=0;i<4;i++) sol[i] = (Float_t)beta[i][0];//should this really be Float_t???
}
void BrLocalTrackDC::CopyTrahit1(BrTrackHit2D* OldTr,BrTrackHit2D* NewTr)
{
//Copy trahit entries excluding the last entry in hitcmb.
Int_t itr;
Int_t NumTraHitFound;
TObjArray TraHitFound;
BrTrackHit* trahit_p;
BrCombinedHit* hitcmb_p;
Int_t NumTraHit;
//First, figure out where the last entry is
NumTraHit = TrackHits->GetEntries();
for(itr=0;itr<NumTraHit;itr++) {
trahit_p = (BrTrackHit*)TrackHits->At(itr);
if(trahit_p->GetTra2d() == OldTr) {
TraHitFound.Add(trahit_p);
}
}
NumTraHitFound = TraHitFound.GetEntries();
for(itr=0;itr<NumTraHitFound-1;itr++) {
hitcmb_p = ((BrTrackHit*)TraHitFound.At(itr))->GetHitcmb();
if(hitcmb_p != 0) {
// trahit[NumTraHit] = new traHIT();
// trahit_p = trahit[NumTraHit];
// numtrahit = NumTraHit;
// NumTraHit++;
// if(NumTraHit > maxtrahit) cout<<"1 NumTraHit too large for "<<Name<<endl;
BrClonesArray &trahit = *TrackHits;
Int_t ntrahit = TrackHits->GetEntries();
trahit_p = new(trahit[ntrahit]) BrTrackHit();
trahit_p->SetHitcmb(hitcmb_p);
trahit_p->SetTra2d(NewTr);
// Set up relation with loctra
// numtrahit = NewTr->numtrahit;
// if(numtrahit > maxloctrahit) printf("NewTr->numtrahit too largen");
// NewTr->trahit[numtrahit] = trahit_p;
// NewTr->numtrahit++;
NewTr->AddTrackHit(trahit_p);
}
}
TraHitFound.Clear(); //clear so destructor does not delete trahits
}
void BrLocalTrackDC::CopyTrahit1(BrLocalTrack* OldTr,BrLocalTrack* NewTr)
{
//Copy trahit entries excluding the last entry in hitcmb.
Int_t itr;
Int_t NumTraHitFound;
TObjArray TraHitFound;
BrTrackHit* trahit_p;
BrCombinedHit* hitcmb_p;
Int_t NumTraHit;
NumTraHitFound = 0;
//First, figure out where the last entry is
NumTraHit = TrackHits->GetEntries();
for(itr=0;itr<NumTraHit;itr++) {
trahit_p = (BrTrackHit*)TrackHits->At(itr);
if( trahit_p->GetTrID() == OldTr ) {
TraHitFound.Add(trahit_p);
}
}
NumTraHitFound = TraHitFound.GetEntries();
for(itr=0;itr<NumTraHitFound-1;itr++) {
hitcmb_p = ((BrTrackHit*)TraHitFound.At(itr))->GetHitcmb();
if(hitcmb_p != 0) {
// trahit[NumTraHit] = new traHIT();
// trahit_p = trahit[NumTraHit];
// NumTraHit++;
// if(NumTraHit > maxtrahit) cout<<"2 NumTraHit too large for "<<Name<<endl;
BrClonesArray &trahit = *TrackHits;
Int_t ntrahit = TrackHits->GetEntries();
trahit_p = new(trahit[ntrahit]) BrTrackHit();
trahit_p->SetHitcmb(hitcmb_p);
trahit_p->SetTrID(NewTr);
trahit_p->SetLoctra(NewTr);
// Set up relation with loctra
// numtrahit = NewTr->numtrahit;
// if(numtrahit>maxloctrahit) printf("NewTr->trahit too largen");
// NewTr->trahit[numtrahit] = trahit_p;
// NewTr->numtrahit++;
NewTr->AddTrackHit(trahit_p);
}
}
}
void BrLocalTrackDC::InsertTrahit(BrCombinedHit *hitcmbb,BrTrackHit2D* tra2dd)
{
BrTrackHit* trahit_p;
BrClonesArray &trahit = *TrackHits;
Int_t ntrahit = TrackHits->GetEntries();
trahit_p = new(trahit[ntrahit]) BrTrackHit();
trahit_p->SetImod(hitcmbb->GetImod());
trahit_p->SetTrID(0);
//do we really need these two or does it just arrive because
//of differences from FORTRAN version? Do it for the moment.
trahit_p->SetTra2d(tra2dd);
trahit_p->SetHitcmb(hitcmbb);
tra2dd->AddTrackHit(trahit_p);
}
void BrLocalTrackDC::RemoveTrack(BrTrackHit2D* tra2dd)
{
//This routine should be more or less debugged although
//the results don't match exactly SONATA after the first time
//through. I think for the moment that is because of the
//different way it is done here and in using ADAMO. This should
//be monitored
Int_t NumTraHit,NumHitCMB;
Int_t i,j;
//Int_t MOD(Int_t,Int_t);
BrTrackHit* trahit_p;
BrCombinedHit* hitcmb_p;
NumTraHit = TrackHits->GetEntries();
for( i=0;i<NumTraHit;i++) {
trahit_p = (BrTrackHit*)TrackHits->At(i);
if(trahit_p->GetTra2d() == tra2dd) {
// First, find if any hitcmb corresponds to this and reset stat
NumHitCMB = CombinedHits->GetEntries();
for(j=0;j<NumHitCMB;j++) {
hitcmb_p = (BrCombinedHit*)CombinedHits->At(j);
if( trahit_p->GetHitcmb() == hitcmb_p ) {
hitcmb_p->SetStat(TSonataMath::MOD(hitcmb_p->GetStat(),16) + (hitcmb_p->GetStat()/16-1)*16);
}
}
// Now Delete this track.
// delete trahit[i];
// NumTraHit--;
// Now fill in, so we don't have blanks
// for(j=i;j<NumTraHit;j++) trahit[j] = trahit[j+1];
// i--; //Need to correct so that we check all of them
TrackHits->RemoveAt(i);
}
}
TrackHits->Compress();
tra2dd->SetStat(999);
}
void BrLocalTrackDC::DeleteBadTrack()
{
int i;
Int_t NumTra2D;
BrTrackHit2D* tra2d_p;
if(DebugLevel()>5) cout<<"DeleteBadTrack: num entries = "<<TrackHits2D->GetEntries()<<endl;
NumTra2D = TrackHits2D->GetEntries();
for(i=0;i<NumTra2D;i++) {
tra2d_p = (BrTrackHit2D*)TrackHits2D->At(i);
if(DebugLevel()>5) cout<<"i = "<<i<<"GetStat()= "<<tra2d_p->GetStat()<<endl;
if( tra2d_p->GetStat() == 999 ) {
// delete tra2d[i];
// NumTra2D--;
// Now fill in, so we don't have blanks
// for(j=i;j<NumTra2D;j++) tra2d[j] = tra2d[j+1];
// i--; //Need to correct so that we check all of them
TrackHits2D->RemoveAt(i);
}
}
TrackHits2D->Compress();
}
void BrLocalTrackDC::GetTra2DMod(Int_t stat,Int_t flg,BrDataTable *Tra2DMod)
{
Int_t i;
Int_t NumTra2D;
BrTrackHit2D *tra2d_p;
Tra2DMod->Clear();
NumTra2D = TrackHits2D->GetEntries();
for(i=0;i<NumTra2D;i++) {
tra2d_p = (BrTrackHit2D*)TrackHits2D->At(i);
if(( tra2d_p->GetStat() == stat || stat == IANY ) &&
( tra2d_p->GetFlg() == flg || flg == IANY )) {
Tra2DMod->Add(tra2d_p);
}
}
}
ROOT page - Class index - Top of the page
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.