#include "BrLocalTrackTPC.h"
#include "BrDetectorParamsTPC.h"
#include "BrDetectorVolume.h"
#include "BrEventNode.h"
#include "btwn.h"
#include "BrTPCCluster.h"
#include "TSonataMath.h"
ClassImp(BrLocalTrackTPC)
//////////////////////////////////////////////////////////////////////////////////
// //
// The BrLocalTrackTPC class provides an interface to the BRAHMS TPC //
// Local Tracking //
// //
//////////////////////////////////////////////////////////////////////////////////
BrLocalTrackTPC::BrLocalTrackTPC():BrLocalTrackModule()
{
// Default Constructor.
h_TPCpsig = NULL;
h_TPCtsig = NULL;
h_TPCpsig_tsig = NULL;
}
BrLocalTrackTPC::BrLocalTrackTPC(Text_t *Name,Char_t *Title)
:BrLocalTrackModule(Name,Title)
{
fParams_p = NULL;
fVolume_p = NULL;
Char_t HistName[32];
sprintf(HistName,"TPCpsig%sx0",GetName());
cout<<"Creating histogram "<<HistName<<endl;
h_TPCpsig = new TH1F(HistName,HistName,50,-.5,2.5);
sprintf(HistName,"TPCtsig%sx0",GetName());
cout<<"Creating histogram "<<HistName<<endl;
h_TPCtsig = new TH1F(HistName,"tsig",50,-.5,2.5);
sprintf(HistName,"TPCpsig_tsig%sx0",GetName());
cout<<"Creating histogram "<<HistName<<endl;
h_TPCpsig_tsig = new TH2F(HistName,"psig vs tsig",50,-.5,2.5,50,-.5,2.5);
}
BrLocalTrackTPC::~BrLocalTrackTPC()
{
// 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;
if( fVolume_p != 0)
delete fVolume_p;
if(h_TPCpsig) delete h_TPCpsig;
if(h_TPCtsig) delete h_TPCtsig;
if(h_TPCpsig_tsig) delete h_TPCpsig_tsig;
}
void BrLocalTrackTPC::SetDetectorParamsTPC(BrDetectorParamsTPC* 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 BrDetectorParamsTPC("Dummy","TPCINT");
*(BrDetectorParamsTPC*)fParams_p = *par;
// The conversion cast is needed because fParams_p is of
// type BrDetectorParams.
}
void BrLocalTrackTPC::SetDetectorParamsTPC(const BrDetectorParamsTPC& 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 (BrDetectorParamsTPC*) fParams_p;
fParams_p = new BrDetectorParamsTPC("Dummy","ParamsTPC");
*(BrDetectorParamsTPC*)fParams_p = par;
// The conversion cast is needed because fParams_p is of
// type BrDetectorParams.
}
void BrLocalTrackTPC::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 BrLocalTrackTPC::Event(BrEventNode* InputTable,BrEventNode* OutputTable)
{
//The routine that is executed event by event.
//First local tables are cleared, then the different subroutines are executed
//The data path is:
//TPCPadCluster:
//Get DigitizedTPC from the Input Table
//DigitizedTPC table has a list of BrDigTPC
//Generate TPCCluster intermediate table with entries of BrTPCClusters
//From that information, generate list of CombinedHits
//
//TPCPadTrack
//Take CombinedHits and generate LocalTracks in TPCPadRecon
//Use LocalTracks to generate VirtualTracks
//
//FillLocalTracks:
//Use VirtualTracks generated to minimize on chisq and get list of detector tracks
Int_t status;
//These next three need to be loaded in. Are they done in LoadCuts??
static Float_t dx = (Float_t)1.5;
static Float_t dy = (Float_t)0.6;
static Int_t maxmiss = 2;
//LoadCuts();????? Appears never to come there in FORTRAN Sonata
CombinedHits->Clear();
LocalTracks->Clear();
TrackHits->Clear();
VirtualTracks->Clear();
TPCPadCluster(InputTable);
status = TPCPadTrack(0,maxmiss,dx,dy);
if(!status) FillLocalTracks(OutputTable);
}
Int_t BrLocalTrackTPC::TPCPadTrack(Int_t idet,Int_t maxmiss,Float_t dx,Float_t dy)
{
//This routine calls TPCPadRecon in a loop over the number of pads a track is
//allowed to miss before throwing away
//Then virtual tracks are filled.
//Then loop over the number of local tracks previously found, generating
//the chisqr
//
//Local tables needed:
//CombinedHits
//Local tables generated:
//LocalTracks
Int_t i,nrow;
Float_t sol[4];
Float_t loctra_pos[3];
Float_t loctra_vec[3];
Float_t covar[4][4];
Float_t chisq;
Int_t ifail;
Int_t status;
Int_t NumLoctra;
BrDetectorParamsTPC *par;
BrDetectorVolume *vol;
BrDataTable HitCMBMod;
BrLocalTrack *loctra_p;
//Get stuff from the input table
Char_t TableName[32];
sprintf(TableName,"%s Combined Hitsx0",GetName());
//Clear before starting new ones.
TrackHits->Clear();
LocalTracks->Clear();
VirtualTracks->Clear();
par = GetDetectorParamsTPC();
vol = GetDetectorVolume();
nrow = par->GetNumberOfRows();
status = TPCPadRecon(&HitCMBMod,nrow,1,-1,maxmiss,dx,dy);
if( !status ) {
if(maxmiss > 0 ) {
for(i=1;i<=maxmiss;i++) {
status = TPCPadRecon(&HitCMBMod,nrow-i,i,-1,maxmiss-i+1,dx,dy);
if( status ) return(status);
}
}
}
else return(status);
FillVTracks(nrow-2);
NumLoctra = LocalTracks->GetEntries();
for(i=0;i<NumLoctra;i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
FitLocalTPCTrack(loctra_p,sol,covar[0],&chisq,&ifail);
printf("PTPC 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_vec[0] = sol[0];
loctra_vec[1] = sol[2];
loctra_vec[2] = (Float_t)1.0;
loctra_p->SetPos(loctra_pos);
loctra_p->SetVec(loctra_vec);
loctra_p->SetChisq(chisq);
}
HitCMBMod.Clear();
return(0);
}
void BrLocalTrackTPC::FitLocalTPCTrack(BrLocalTrack* loctraa,Float_t *sol,Float_t *covar,Float_t *chisq,Int_t *ifail)
{
//FitLocalTPCTrack takes the input Local Track, extracts the combined hits that
//make up this track and does a fit to get a chisqr
//Input arguments:
//BrLocalTrack *loctraa : the local track to be fit
//Output arguments:
//Float_t *sol : output parameters from the fit
//Float_t *covar : covar
//Input tables:
//None
//Output tables:
//None
Int_t i,j,ic;
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,dchi1,dchi2;
Int_t NumHitCMBSel;
BrCombinedHit *hitcmb_p;
NumHitCMBSel = loctraa->GetTrackEntries();
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];
}
//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 = loctraa->GetTrackHitAt(ic)->GetHitcmb();
a[0] = hitcmb_p->GetPos()[2];
a[1] = 1;
a[2] = hitcmb_p->GetPos()[2];
a[3] = 1;
du = (Float_t)0.02;
for(i=0;i<2;i++) {
for(j=0;j<2;j++) alpha[i][j] += a[i]*a[j]/du/du;
beta[i][0] += a[i] * hitcmb_p->GetPos()[0]/du/du;
}
for(i=2;i<4;i++) {
for(j=2;j<4;j++) alpha[i][j] += a[i]*a[j]/du/du;
beta[i][0] += a[i] * hitcmb_p->GetPos()[1]/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 = loctraa->GetTrackHitAt(ic)->GetHitcmb();
a[0] = hitcmb_p->GetPos()[2];
a[1] = 1;
a[2] = hitcmb_p->GetPos()[2];
a[3] = 1;
du = (Float_t)0.02;
dchi = (Float_t)0.0;
dchi1 = (Float_t)(beta[1][0] + beta[0][0]*a[0] - hitcmb_p->GetPos()[0]);
dchi2 = (Float_t)(beta[3][0] + beta[2][0]*a[2] - hitcmb_p->GetPos()[1]);
dchi = (dchi1*dchi1 + dchi2*dchi2)/du;
chisqq += dchi*dchi;
}
*chisq = chisqq;
for(i=0;i<4;i++) sol[i] = (Float_t)beta[i][0];//should this really be Float_t???
}
Int_t BrLocalTrackTPC::TPCPadRecon(BrDataTable *HitCMBMod,Int_t istart,Int_t ilast,Int_t idir,Int_t maxmiss,Float_t dx,Float_t dy)
{
//Routine to generate trackhits.
//First loop over combined hits selected on istart
//Generate local tracks
//Loop over local tracks
//Input parameters
//BrDataTable *HitCMBMod
//Int_t istart
//Int_t ilast
//Int_t idir
//Int_t maxmiss : Max number of pads that can be missed before throwing away
//Float_t dx
//Float_t dy
//Input Tables
//Output Tables
//LocalTracks
//BrTrackHit generated in InsertTrahit
Int_t i,j,jj,imod,irow,next_row;
Int_t ncand;
Int_t iclist[maxcand];
Float_t proj[2];
Float_t dZ,ax,ay;
Float_t temp_pos[3];
Float_t temp_vec[3];
Float_t loctra_pos[3];
Float_t loctra_vec[3];
Int_t temp_nhit,jcmb;
Int_t NumLoctra,NumLoctraInt,NumHitCMBMod;
BrDetectorParamsTPC *par;
BrDetectorVolume *vol;
BrCombinedHit* hitcmb_p;
BrLocalTrack* loctra_p;
BrLocalTrack* IDsave;
par = GetDetectorParamsTPC();
vol = GetDetectorVolume();
imod = 100 + istart;
HitCMBMod->Clear();
GetHitCMBMod(imod,HitCMBMod);
NumHitCMBMod = HitCMBMod->GetEntries();
for(i=0;i<NumHitCMBMod;i++) {
hitcmb_p = (BrCombinedHit*)HitCMBMod->At(i);
if(hitcmb_p->GetStat()/16 == 0) {
hitcmb_p->SetStatBit(1<<4);
BrClonesArray &loctra = *LocalTracks;
Int_t nloctra = LocalTracks->GetEntries();
loctra_p = new(loctra[nloctra]) BrLocalTrack();
loctra_p->SetPos(hitcmb_p->GetPos());
loctra_p->SetNhit(1);
loctra_p->SetFlg(istart);
loctra_p->SetStat(0);
InsertTrahit(hitcmb_p,loctra_p);
}
}
//This part needs to be checked in order to have arbitrary
//direction. fortran equivalent for i = istart,ilast,idir
//where istart might be greater than ilast and idir=-1 or
//istart is less than ilast and idir = +1. for now we do it
//this way because everything I have found in Sonata calling
//is idir=-1
//for(irow=istart;irow>=ilast-idir;irow +=idir) {
for(irow=istart;irow>=ilast-idir;irow +=idir) {
next_row = irow + idir;
NumLoctra = LocalTracks->GetEntries();
NumLoctraInt = NumLoctra;
for(i=0;i<NumLoctraInt;i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
if(loctra_p->GetStat() == 0) {
IDsave = loctra_p;
dZ = par->GetRowDistance()*(next_row - loctra_p->GetFlg());
if( loctra_p->GetNhit() == 1 ) {
proj[0] = loctra_p->GetPos()[0];
proj[1] = loctra_p->GetPos()[1];
}
else {
ax = loctra_p->GetVec()[0] / loctra_p->GetVec()[2];
ay = loctra_p->GetVec()[1] / loctra_p->GetVec()[2];
proj[0] = loctra_p->GetPos()[0] + ax * dZ;
proj[1] = loctra_p->GetPos()[1] + ay * dZ;
}
TPCPadFindHit(HitCMBMod,next_row,proj,
dx/(Float_t)sqrt((Double_t)loctra_p->GetNhit()),
dy/(Float_t)sqrt((Double_t)loctra_p->GetNhit()),
&ncand,maxcand,iclist);
if(ncand == 0) {
if( (next_row - istart)/idir+1-loctra_p->GetNhit()>maxmiss) {
RemoveTrack(loctra_p);
}
}
else if(ncand == 1) {
// We now have one loctra which can have one or more hits in
// the present row.
// If ncand is 1, ie only one hit, the algorithm is simple
// otherwise we should save the present local track in a
// local copy before updating the information. Otherwise
// the nhit information and position is wrong.
j = iclist[0];
hitcmb_p = (BrCombinedHit*)HitCMBMod->At(j);
hitcmb_p->SetStat(TSonataMath::MOD(hitcmb_p->GetStat(),16)
+ (hitcmb_p->GetStat()/16+1)*16);
loctra_vec[0] = hitcmb_p->GetPos()[0] - loctra_p->GetPos()[0];
loctra_vec[1] = hitcmb_p->GetPos()[1] - loctra_p->GetPos()[1];
loctra_vec[2] = hitcmb_p->GetPos()[2] - loctra_p->GetPos()[2];
loctra_p->AddVec(loctra_vec);
// Update position to be on the last plane
// We take a weighted average on proj[0], proj[1] and
// then find the actual position.
temp_pos[0] = (hitcmb_p->GetPos()[0] + loctra_p->GetNhit()*proj[0])/(loctra_p->GetNhit()+1);
temp_pos[1] = (hitcmb_p->GetPos()[1] + loctra_p->GetNhit()*proj[1])/(loctra_p->GetNhit()+1);
temp_pos[2] = hitcmb_p->GetPos()[2];
loctra_p->SetPos(temp_pos);
loctra_p->IncNhit();
loctra_p->SetFlg(next_row);
loctra_p->SetStat(0);
InsertTrahit(hitcmb_p,loctra_p);
}
else {
// There is more than one hit
// Store present value in local vector
IDsave = loctra_p;
for(jj=0;jj<3;jj++) {
temp_pos[jj] = loctra_p->GetPos()[jj];
temp_vec[jj] = loctra_p->GetVec()[jj];
}
temp_nhit = loctra_p->GetNhit();
for(jcmb=0;jcmb<ncand;jcmb++) {
j = iclist[jcmb];
hitcmb_p = (BrCombinedHit*)HitCMBMod->At(j);
hitcmb_p->SetStat(TSonataMath::MOD(hitcmb_p->GetStat(),16)
+ (hitcmb_p->GetStat()/16+1)*16);
if(jcmb == ncand-1) {
// Calculate the new vector as average over previous
// given values.
// Remember, loctra_p was defined on the last pass
// when jcmb was ncand-2
for(jj=0;jj<3;jj++) loctra_vec[jj] = hitcmb_p->GetPos()[jj] - temp_pos[jj] + temp_vec[jj];
loctra_p->SetVec(loctra_vec);
// Update position to be last plane
loctra_pos[0] = (hitcmb_p->GetPos()[0] + temp_nhit*proj[0])/(temp_nhit+1);
loctra_pos[1] = (hitcmb_p->GetPos()[1] + temp_nhit*proj[1])/(temp_nhit+1);
loctra_pos[2] = hitcmb_p->GetPos()[2];
loctra_p->SetPos(loctra_pos);
loctra_p->SetNhit(temp_nhit + 1);
loctra_p->SetFlg(next_row);
loctra_p->SetStat(0);
InsertTrahit(hitcmb_p,loctra_p);
}
else {
BrClonesArray &loctra = *LocalTracks;
Int_t nloctra = LocalTracks->GetEntries();
loctra_p = new(loctra[nloctra]) BrLocalTrack();
for(jj=0;jj<3;jj++) loctra_vec[jj] = hitcmb_p->GetPos()[jj] - temp_pos[jj] + temp_vec[jj];
loctra_p->SetVec(loctra_vec);
loctra_pos[0] = (hitcmb_p->GetPos()[0] + temp_nhit*proj[0])/(temp_nhit+1);
loctra_pos[1] = (hitcmb_p->GetPos()[1] + temp_nhit*proj[1])/(temp_nhit+1);
loctra_pos[2] = hitcmb_p->GetPos()[2];
loctra_p->SetPos(loctra_pos);
loctra_p->SetNhit(temp_nhit+1);
loctra_p->SetFlg(next_row);
loctra_p->SetStat(0);
InsertTrahit(hitcmb_p,loctra_p);
CopyTrahit(IDsave,loctra_p,hitcmb_p->GetImod());
} //jcmb == ncand
} //loop over jcmb (=>j)
} //ncand>=1
} //loctra->stat == 0
} //loctra loop
} //irow loop
NumLoctra = LocalTracks->GetEntries();
for(i=0;i<NumLoctra;i++) {
loctra_p = (BrLocalTrack*)LocalTracks->At(i);
if(loctra_p->GetStat() == 0) loctra_p->SetStat(1);
}
DeleteBadTrack();
return(0);
}
void BrLocalTrackTPC::TPCPadFindHit(BrDataTable *HitCMBMod,Int_t irow,Float_t *proj,Float_t dx,Float_t dy,Int_t *ncand,Int_t maxcandd,Int_t *iclist)
{
//This routine selects CombinedHits on irow
//Loop over selected CombinedHits and take the ones in the list
//where (projx - dx < hitcmb_x < projx + dx) and (projy - dy < hitcmb_y < projy + dy)
//
//Input Arguments:
//BrDataTable *HitCMBMod : Data table to hold CombinedHits selected on padrow
//Int_t irow : the padrow of the TPC
//Float_t *proj :
//Float_t dx : x-spacing to be accepted as a hit
//Float_t dy : y-spacing to be accepted as a hit
//Int_t iclist: list of hitcmb's found
//Int_t ncand: number of hitcmb's found
Int_t i,imod,ncandd;
BrCombinedHit* hitcmb_p;
Int_t NumHitCMBMod;
ncandd = 0;
imod = 100 + irow;
HitCMBMod->Clear();
GetHitCMBMod(imod,HitCMBMod);
NumHitCMBMod = HitCMBMod->GetEntries();
for(i=0;i<NumHitCMBMod;i++) {
hitcmb_p = (BrCombinedHit*)HitCMBMod->At(i);
if(BTWN(proj[0]-dx,hitcmb_p->GetPos()[0],proj[0]+dx)) {
if(BTWN(proj[1]-dy,hitcmb_p->GetPos()[1],proj[1]+dy)) {
if(ncandd<maxcandd) {
iclist[ncandd] = i;
ncandd++;
}
}
}
}
*ncand = ncandd;
}
void BrLocalTrackTPC::TPCPadCluster(BrEventNode* InputTable)
{
//Input Table
//DigitizedTPC : from calling routine
//Intermediate Table
//TPCCluster : has entries of BrTPCCluster to be used in ScanMultiHitClusters
Int_t irow,nrow;
Int_t padadc[maxpad][maxtime];
BrDigTPC* padid [maxpad][maxtime];
Int_t nsegment[maxpad];
Int_t first_segment[maxpad];
Int_t last_segment [maxpad];
Int_t ihit,ihitid,first_pad,last_pad,iseg;
Int_t ip,it,ir,is,cid,cid2;
Int_t ok,in_seg;
Int_t it1,it1a,it2,it2a,is1,lseg,ic,npix;
Int_t next_segment;
Float_t tmin,tmax,pmin,pmax,tsum,tav,psum,pav,t1,t2,dsum;
Float_t psig; //RMS in pad direction
Float_t tsig; //RMS in time direction
//Lower and upper limits of psig and tsig. These are presently hardwired, later
//they should be read in.
static const Float_t psig_lowlimit = (Float_t)0.25;
static const Float_t psig_highlimit = (Float_t)0.90;
static const Float_t tsig_lowlimit = (Float_t)0.50;
static const Float_t tsig_highlimit = (Float_t)1.25;
Int_t pdadc;
static const Float_t adc_threshold = (Float_t)1.;
Int_t NumDigTPCMod;
#define maxseg 1000
struct {
Int_t time1;
Int_t time2;
Int_t row;
Int_t clid;
Int_t next;
} segment[maxseg];
BrTPCCluster* clus_p;
BrDataTable DigTPCMod;
BrDigTPC *digtpc_p;
BrDetectorParamsTPC *par;
BrDetectorVolume *vol;
Float_t detZ;
BrDataTable *DigitizedTPC;
BrDataTable TPCCluster;
//Get stuff from the input table
Char_t TableName[32];
sprintf(TableName,"DigitizedTPC %sx0",GetName());
DigitizedTPC = InputTable->GetDataTable(TableName);
if(DigitizedTPC == 0) {
if(DebugLevel() > 0 ) cout<<"No DigitizedTPC in this event for "<<GetName()<<endl;
return;
}
par = GetDetectorParamsTPC();
vol = GetDetectorVolume();
detZ = (Float_t)0.0; //For now with TPC for test
TPCCluster.Clear(); //Be sure it is clear before starting
nrow = par->GetNumberOfRows();
ihitid = 0;
for(irow=0;irow<nrow;irow++) {
ihitid++;
TSonataMath::vzero((Int_t *)&padadc[0][0],maxpad*maxtime);
TSonataMath::vzero((Int_t *)&padid[0][0] ,maxpad*maxtime);
DigTPCMod.Clear();
GetDigTPCMod(irow,IANY,IANY,DigitizedTPC,&DigTPCMod);
NumDigTPCMod = DigTPCMod.GetEntries();
for( ihit=0;ihit<NumDigTPCMod;ihit++) {
digtpc_p = (BrDigTPC*)DigTPCMod.At(ihit);
ip = digtpc_p->GetPadno();
it = (Int_t)digtpc_p->GetTime();
if(digtpc_p->GetAdc() >= adc_threshold ) {
if( ip > maxpad ) {
cout<<"Pad outside limits, it is "<<ip<<endl;
break; //Same as goto 100 in FORTRAN version
}
if( it > maxtime ) {
cout<<"Time outside limits, it is "<<it<<endl;
break; //Same as goto 100 in FORTRAN version
}
padadc[ip][it] = (Int_t)digtpc_p->GetAdc();
padid [ip][it] = digtpc_p;
}
}
last_pad = maxpad;
first_pad = 1;
// The array padadc contains the ADC values of the Hits
// Create segments for each row which contains the first and last
// points pr padrow which has valid time hits
iseg = 0;
for( ir=0;ir<maxpad;ir++) {
nsegment[ir] = 0;
last_segment[ir] = -1;
it = 0;
ok = kTRUE;
in_seg = kFALSE;
while(ok) {
if(padadc[ir][it] > 0 ) {
if(in_seg) {
segment[iseg-1].time2 = it;
}
else {
iseg++;
if(nsegment[ir] == 0) first_segment[ir] = iseg;
nsegment[ir]++;
segment[iseg-1].row = ir;
segment[iseg-1].clid = 0;
segment[iseg-1].next = NULL_SEGMENT;
segment[iseg-1].time1 = it;
segment[iseg-1].time2 = it;
in_seg = kTRUE;
last_segment[ir] = iseg;
}
}
else {
in_seg = kFALSE;
}
it++;
if(it>=maxtime) ok=kFALSE;
}
}
// Segments are now created for this ROW plane. Now make clusters
// from the segments. First clear the cluster struct.
InitCluster();
for(ir=first_pad;ir<last_pad;ir++) {
for(is=first_segment[ir]-1;is<last_segment[ir];is++) {
if(segment[is].clid == 0) {
// Can all of this be done with classes and pointers?
cid = GetFreeCluster();
cluster[cid].head = is;
cluster[cid].tail = is;
segment[is].next = NULL_SEGMENT;
segment[is].clid = cid;
}
else {
cid = segment[is].clid;
}
// If this is not the last row, look for overlapping segments in
// the following padrow. These segments may already have been
// used. In this case the previous used segments will be added to
// the chain of the ongoing list and their clusterid changed.
// eg. in row ir the two segments likely will have been assigned to
// two different clusters, but really belong to the same.
// row ir : |-------| |-------|
// row ir+1: |-------|
if( ir < last_pad ) {
it1 = segment[is].time1;
it2 = segment[is].time2;
for(is1=first_segment[ir+1]-1;is1<last_segment[ir+1];is1++) {
it1a = segment[is1].time1;
it2a = segment[is1].time2;
if(((it1<=it1a) && (it1a<=it2)) || ((it1a <=it1) && (it1 <= it2a))) {
if(segment[is1].clid == 0 ) {
segment[is1].next = cluster[cid].head;
segment[is1].clid = cid;
cluster[cid].head = is1;
}
else if( segment[is1].clid != cid ) {
// Already in a cluster
cid2 = segment[is1].clid;
lseg = cluster[cid].tail;
// Update segment cluster information
next_segment = cluster[cid2].head;
while(next_segment!=NULL_SEGMENT) {
segment[next_segment].clid = cid;
next_segment = segment[next_segment].next;
}
segment[lseg].next = cluster[cid2].head;
cluster[cid].tail = cluster[cid2].tail;
cluster[cid2].head = NULL_CLUSTER;
cluster[cid2].tail = NULL_CLUSTER;
}
}
}
} //not last pad
} //segments
} //pads
// Calculate cluster parameters
for(ic=1;ic<max_cluster;ic++) {
if(cluster[ic].head != NULL_SEGMENT) {
npix = 0;
tmax = (Float_t)0.0;
tmin = (Float_t)maxtime;
pmin = (Float_t)last_pad;
pmax = (Float_t)0.0;
tsum = (Float_t)0.0;
tav = (Float_t)0.0;
tsig = (Float_t)0.0;
psum = (Float_t)0.0;
pav = (Float_t)0.0;
psig = (Float_t)0.0;
// Loop over segments belonging to this cluster
next_segment = cluster[ic].head;
while(next_segment != NULL_SEGMENT) {
ir = segment[next_segment].row;
it1 = segment[next_segment].time1;
it2 = segment[next_segment].time2;
t1 = (Float_t)it1;
t2 = (Float_t)it2;
pmin = TMath::Min(pmin,(Float_t)ir);
pmax = TMath::Max(pmax,(Float_t)ir);
tmin = TMath::Min(tmin,t1);
tmax = TMath::Max(tmax,t2);
dsum = (Float_t)0.0;
for(it=it1;it<=it2;it++) {
pdadc = padadc[ir][it];
tsum += pdadc;
dsum += pdadc;
tav += pdadc * (Float_t)it;
tsig += pdadc * (Float_t)it * (Float_t)it;
npix++;
}
psum += dsum;
pav += dsum * (Float_t)ir;
psig += dsum * (Float_t)ir * (Float_t)ir;
next_segment = segment[next_segment].next;
}
// Analyze this segment
pav = pav / psum;
psig = psig / psum - pav * pav;
psig = (Float_t)sqrt(psig);
tav = tav / tsum;
tsig = tsig / tsum - tav * tav;
tsig = (Float_t)sqrt(tsig);
// Create clusters and relations to digtpc (for future salvage)
clus_p = new BrTPCCluster();
TPCCluster.Add(clus_p);
clus_p->SetGvID(0); //This needs to be addressed
clus_p->SetRow(irow);
clus_p->SetPav(pav);
clus_p->SetPsig(psig);
clus_p->SetTav(tav);
clus_p->SetTsig(tsig);
clus_p->SetDedx(tsum);
clus_p->SetNpix(npix);
clus_p->SetSubcl(0);
// Asociate a quality to signal
// 3 == noise or on edge of detector
// 2 == likely multihit
// 1 == good single hit
// The cut parameters are presently hard wired; this must be changed
// if ( psig < 0.25 || tsig < 1.4 ) clus_p->SetStat(3);
// else if( psig < 0.90 && tsig < 2.4 ) clus_p->SetStat(1);
// else clus_p->SetStat(2);
if (psig < psig_lowlimit || tsig < tsig_lowlimit) clus_p->SetStat(3);
else if (psig <psig_highlimit && tsig <tsig_highlimit) clus_p->SetStat(1);
else clus_p->SetStat(2);
h_TPCpsig->Fill(psig);
h_TPCtsig->Fill(tsig);
h_TPCpsig_tsig->Fill(psig,tsig);
// Insert relation for digtpc which contributes to this cluster
// by looping over segments etc for the cluster again.
next_segment = cluster[ic].head;
while(next_segment != NULL_SEGMENT ) {
ir = segment[next_segment].row;
it1 = segment[next_segment].time1;
it2 = segment[next_segment].time2;
for(it=it1;it<=it2;it++) {
padid[ir][it]->SetClus(clus_p);
}
next_segment = segment[next_segment].next;
}
} //cluster[ic].head !=0
} //clusters
} //rows
//Analyze clusters; look for non-single hit clusters;
//use the prescribe method
ScanMultiHitClusters(DigitizedTPC,&TPCCluster);
//Fill the single clusters into HITCMB
TPCClusterToHitCMB(&TPCCluster);
DigTPCMod.Clear();
//TPCCluster.Clear();
TPCCluster.Delete(); //delete clusters as we do not need them anymore
}
Int_t BrLocalTrackTPC::GetFreeCluster()
{
Int_t i = 1;
while(i<max_cluster) {
if( cluster[i].head == NULL_CLUSTER ) return(i);
i++;
}
return(-1); //Error return; need to be sure and handle in calling routine
}
void BrLocalTrackTPC::InitCluster()
{
for(Int_t i = 0;i<max_cluster;i++) {
cluster[i].head = NULL_CLUSTER;
cluster[i].tail = NULL_CLUSTER;
}
}
void BrLocalTrackTPC::ScanMultiHitClusters(BrDataTable* DigitizedTPC,BrDataTable* TPCCluster)
{
struct {
Int_t ix;
Int_t iy;
} list[40][2];
Int_t nrow,irow,ip,it,ih,ic,ok,ix,iy,ixm,iym;
BrTPCCluster* cluster_ID;
Int_t sub_cluster,next_item,present_item,present,nextp,item;
Int_t j,ixc,iyc,ics,npix;
Float_t vmax,tmin,tmax,pmin,pmax,tsum,tav,tsig,psum,pav,psig;
static const Float_t adc_threshold = (Float_t)1.;
static const Int_t dirx[4] = { 1, 0,-1, 0};
static const Int_t diry[4] = { 0, 1, 0,-1};
Int_t NumDigTPCMod;
BrDataTable DigTPCMod;
BrDigTPC* digtpc_p;
Int_t NumClusTPCMod;
BrDataTable ClusTPCMod;
BrTPCCluster* clus_p;
BrDetectorParamsTPC *par;
par = GetDetectorParamsTPC();
nrow = par->GetNumberOfRows();
for(irow=0;irow<nrow;irow++) {
for(ip=0;ip<maxpad;ip++) {
for(it=0;it<maxtime;it++) {
pixels[ip][it].ID = 0;
}
}
DigTPCMod.Clear();
GetDigTPCMod(irow,IANY,IANY,DigitizedTPC,&DigTPCMod);
ClusTPCMod.Clear();
GetClusTPCMod(irow,2,TPCCluster,&ClusTPCMod);
NumDigTPCMod = DigTPCMod.GetEntries();
for(ih=0;ih<NumDigTPCMod;ih++) {
digtpc_p = (BrDigTPC*)DigTPCMod.At(ih);
ip = digtpc_p->GetPadno();
it = digtpc_p->GetTime();
if(digtpc_p->GetAdc() > adc_threshold) {
if( ip > maxpad ) {
cout<<"Pad outside limits, it is "<<ip<<endl;
break;
}
if( it > maxtime ) {
cout<<"Time outside limits, it is "<<it<<endl;
break;
}
pixels[ip][it].adcval = (Int_t)digtpc_p->GetAdc();
pixels[ip][it].ID = digtpc_p->GetClus();
pixels[ip][it].subID = 0;
pixels[ip][it].digtpc = digtpc_p;
}
}
NumClusTPCMod = ClusTPCMod.GetEntries();
for( ic=0;ic<NumClusTPCMod;ic++ ) {
cluster_ID = (BrTPCCluster*)ClusTPCMod.At(ic);
ok = GetMaxPixel(cluster_ID,&ixm,&iym);
sub_cluster = 0;
while( ok ) {
next_item = 0;
sub_cluster++;
present_item = 1;
present = 0;
nextp = 1;
list[0][0].ix = ixm;
list[0][0].iy = iym;
// Mark starting pixel as belonging to sub-cluster
pixels[ixm][iym].subID = sub_cluster;
while( present_item > 0 ) {
for(item=0;item<present_item;item++) {
ixc = list[item][present].ix;
iyc = list[item][present].iy;
vmax = (Float_t)pixels[ixc][iyc].adcval;
for( j=0;j<4;j++ ) {
ix = ixc + dirx[j];
iy = iyc + diry[j];
if( ix >= 0 && ix < maxpad &&
iy >= 0 && iy < maxtime ) {
if( pixels[ix][iy].ID == cluster_ID &&
pixels[ix][iy].subID == 0 &&
pixels[ix][iy].adcval <= vmax ) {
list[next_item][nextp].ix = ix;
list[next_item][nextp].iy = iy;
pixels[ix][iy].subID = sub_cluster;
next_item++;
}
}
}
}
// New scan - for next set of rows.
present_item = next_item;
if(nextp == 0) {
nextp = 1;
present = 0;
}
else {
nextp = 0;
present = 1;
}
next_item = 0;
}
// One search has been completed. Find the maximum value for the remaining
// subclusters
ok = GetMaxPixel(cluster_ID,&ixm,&iym);
} //Loop over sub_clusters
for(ics=1;ics<=sub_cluster;ics++) {
npix = 0;
tmax = (Float_t)0.0;
tmin = (Float_t)maxtime;
pmin = (Float_t)maxpad;
pmax = (Float_t)0.0;
tsum = (Float_t)0.0;
tav = (Float_t)0.0;
tsig = (Float_t)0.0;
psum = (Float_t)0.0;
pav = (Float_t)0.0;
psig = (Float_t)0.0;
for(ix=0;ix<maxpad;ix++) {
for(iy=0;iy<maxtime;iy++) {
if( pixels[ix][iy].ID == cluster_ID &&
pixels[ix][iy].subID == ics ) {
tsum += pixels[ix][iy].adcval;
tav += pixels[ix][iy].adcval * (Float_t)iy;
pav += pixels[ix][iy].adcval * (Float_t)ix;
tsig += pixels[ix][iy].adcval * (Float_t)(iy*iy);
psig += pixels[ix][iy].adcval * (Float_t)(ix*ix);
npix++;
}
}
}
pav = pav / tsum;
psig = psig / tsum - pav*pav;
psig = (Float_t)sqrt(psig);
tav = tav / tsum;
tsig = tsig / tsum - tav*tav;
tsig = (Float_t)sqrt(tsig);
// Add this info as a new cluster to the present set
clus_p = new BrTPCCluster();
TPCCluster->Add(clus_p);
clus_p->SetGvID(0); //This needs to be addressed
clus_p->SetRow(irow);
clus_p->SetPav(pav);
clus_p->SetPsig(psig);
clus_p->SetTav(tav);
clus_p->SetTsig(tsig);
clus_p->SetDedx(tsum);
clus_p->SetNpix(npix);
clus_p->SetSubcl(cluster_ID);
clus_p->SetStat(4);
// Insert the relations to digtpc by looping over the pixel
// information.
for(ix=0;ix<maxpad;ix++) {
for(iy=0;iy<maxtime;iy++) {
if(pixels[ix][iy].ID == cluster_ID &&
pixels[ix][iy].subID == ics ) {
BrDigTPC *digtpcins = pixels[ix][iy].digtpc;
digtpcins->SetClus(clus_p);
clus_p->SetDigtpc(digtpcins);
}
}
}
}
}
}
DigTPCMod.Clear();
ClusTPCMod.Clear();
}
Int_t BrLocalTrackTPC::GetMaxPixel(BrTPCCluster* cluster_ID,Int_t *ixm,Int_t *iym)
{
Int_t ix,iy;
Float_t vmax;
//This routine finds the maximum value and the pixel for a given
//cluster. Only pixels which have not been assigned to a SUBCLUSTER
//are considered.
*ixm = 0;
*iym = 0;
vmax = (Float_t)0.0;
for(ix=0;ix<maxpad;ix++) {
for(iy=0;iy<maxtime;iy++) {
if(pixels[ix][iy].ID == cluster_ID) {
if(pixels[ix][iy].subID == 0 && pixels[ix][iy].adcval>vmax) {
vmax = (Float_t)pixels[ix][iy].adcval;
*ixm = ix;
*iym = iy;
}
}
}
}
if(*ixm == 0) return(kFALSE);
return(kTRUE);
}
void BrLocalTrackTPC::TPCClusterToHitCMB(BrDataTable* TPCCluster)
{
//Generate BrCombinedHit elements from information in TPCCluster
Float_t xbase,ybase,xb,cm_to_timeb,detZ;
Int_t irow,ic;
Float_t hitcmb_pos[3],hitcmb_dpos[3];
Int_t NumClusTPCMod;
BrDataTable ClusTPCMod;
BrTPCCluster* clus_p;
BrDetectorParamsTPC *par;
BrDetectorVolume *vol;
BrCombinedHit *hitcmb_p;
par = GetDetectorParamsTPC();
vol = GetDetectorVolume();
if( vol->GetCut() == 0 ) {
ybase = -vol->GetSize()[1] / (Float_t)2.;
xbase = -vol->GetSize()[0] / (Float_t)2.;
}
else {
ybase = -vol->GetFsize()[1] / (Float_t)2.;
xbase = -vol->GetFsize()[0] / (Float_t)2.;
}
xb = -par->GetPadsprow() * par->GetPadDistance();
xbase = TMath::Max(xb,xbase);
cm_to_timeb = (Float_t)1.0/(par->GetDriftv() * par->GetTimeBucket());
//detZ = DetectorElementPosition[idet];
detZ = (Float_t)0.0;
for(irow=0;irow<par->GetNumberOfRows();irow++) {
// First, pick up all the single clusters
ClusTPCMod.Clear();
GetClusTPCMod(irow,1,TPCCluster,&ClusTPCMod);
NumClusTPCMod = ClusTPCMod.GetEntries();
for(ic=0;ic<NumClusTPCMod;ic++) {
clus_p = (BrTPCCluster*)ClusTPCMod.At(ic);
BrClonesArray &hitcmb = *CombinedHits;
Int_t nhitcmb = CombinedHits->GetEntries();
hitcmb_p = new(hitcmb[nhitcmb]) BrCombinedHit();
hitcmb_pos[0] = (clus_p->GetPav() - (Float_t)0.5)*par->GetPadDistance() + xbase;
hitcmb_pos[1] = (clus_p->GetTav() + (Float_t)0.5)/cm_to_timeb + ybase;
hitcmb_pos[2] = detZ + par->GetRowDistance()*(irow+1-(par->GetNumberOfRows()+1)/(Float_t)2.0);
hitcmb_dpos[0] = clus_p->GetPsig()*par->GetPadDistance()/(Float_t)10.;
hitcmb_dpos[1] = (clus_p->GetTsig()/cm_to_timeb)/(Float_t)10.;
hitcmb_p->SetPos(hitcmb_pos);
hitcmb_p->SetDpos(hitcmb_dpos);
hitcmb_p->SetImod(100 + irow + 1);
hitcmb_p->SetNhit(1);
hitcmb_p->SetStat(1);
hitcmb_p->SetType(3);
}
// Now pick up the multi-hit clusters
// they are stored as stat=4
// The primary multihit clusters are stored as stat=2
ClusTPCMod.Clear();
GetClusTPCMod(irow,4,TPCCluster,&ClusTPCMod);
NumClusTPCMod = ClusTPCMod.GetEntries();
for(ic=0;ic<NumClusTPCMod;ic++) {
clus_p = (BrTPCCluster*)ClusTPCMod.At(ic);
BrClonesArray &hitcmb = *CombinedHits;
Int_t nhitcmb = CombinedHits->GetEntries();
hitcmb_p = new(hitcmb[nhitcmb]) BrCombinedHit();
hitcmb_pos[0] = (clus_p->GetPav() - (Float_t)0.5)*par->GetPadDistance() + xbase;
hitcmb_pos[1] = (clus_p->GetTav() + (Float_t)0.5)/cm_to_timeb + ybase;
hitcmb_pos[2] = detZ + par->GetRowDistance()*(irow+1-(par->GetNumberOfRows()+1)/(Float_t)2.0);
hitcmb_dpos[0] = clus_p->GetPsig()*par->GetPadDistance()/(Float_t)10.;
hitcmb_dpos[1] = (clus_p->GetTsig()/cm_to_timeb)/(Float_t)10.;
hitcmb_p->SetPos(hitcmb_pos);
hitcmb_p->SetDpos(hitcmb_dpos);
hitcmb_p->SetImod(100 + irow + 1);
hitcmb_p->SetNhit(1);
hitcmb_p->SetStat(1);
hitcmb_p->SetType(3);
}
}
}
void BrLocalTrackTPC::GetDigTPCMod(Int_t irow,Int_t padno,Int_t time,BrDataTable* DigTPC,BrDataTable* DigTPCMod)
{
//Generate a list of BrDigTPC elements out of DigTPC table selected on irow, padno and time
BrDigTPC *digtpc_p;
Int_t NumDigTPC;
Int_t i;
NumDigTPC = DigTPC->GetEntries();
for(i=0;i<NumDigTPC;i++) {
digtpc_p = (BrDigTPC*)DigTPC->At(i);
if(( digtpc_p->GetIrow() == irow || irow == IANY ) &&
( digtpc_p->GetPadno() == padno || padno == IANY ) &&
( digtpc_p->GetTime() == time || time == IANY )) {
DigTPCMod->Add(digtpc_p);
}
}
//Now sort them in order;
DigTPCMod->Sort();
}
void BrLocalTrackTPC::GetClusTPCMod(Int_t irow,Int_t stat,BrDataTable *TPCCluster,BrDataTable *ClusTPCMod)
{
//Get a list of BrTPCCluster elements out of TPCCluster table selected on irow and stat
BrTPCCluster *clus_p;
Int_t i,NumClus;
NumClus = TPCCluster->GetEntries();
for(i=0;i<NumClus;i++) {
clus_p = (BrTPCCluster*)TPCCluster->At(i);
if(( clus_p->GetRow() == irow || irow == IANY ) &&
( clus_p->GetStat() == stat || stat == IANY ) ) {
ClusTPCMod->Add(clus_p);
}
}
//ClusTPCMod->Sort();//probabally not needed here
}
void BrLocalTrackTPC::FillHitCMBSelectorr( BrLocalTrack* loctra,BrDataTable *HitCMBSel)
{
int ic;
//Find the hitcmb which are related through the loctra part.
//i.e. the non-mainview planes.
Int_t NumTraHit;
BrTrackHit *trahit_p;
HitCMBSel->Clear();
NumTraHit = TrackHits->GetEntries();
for(ic=0;ic<NumTraHit;ic++) {
trahit_p = (BrTrackHit*)TrackHits->At(ic);
if(loctra == trahit_p->GetLoctra() ) {
if( trahit_p->GetHitcmb() ) {
HitCMBSel->Add(trahit_p->GetHitcmb());
}
}
}
//Next statement needs to be activated when Tra2D is activated
Int_t isilly = 0;
if( isilly == 0 )return;
/*
for(ic=0;ic<NumTraHit;ic++) {
trahit_p = (BrTrackHit*)TrackHits.At(ic);
if(loctra->GetTr1 == trahit_p->GetTra2d() ) {
if( trahit_p->GetHitcmb() ) {
HitCMBSel->Add(trahit_p->GetHitcmb());
}
}
}
for(ic=0;ic<NumTraHit;ic++) {
trahit_p = (BrTrackHit*)TrackHits.At(ic);
if(loctra_p->GetTr2() == trahit_p->GetTra2d() ) {
if( trahit_p->GetHitcmb() ) {
HitCMBSel->Add(trahit_p->GetHitcmb());
}
}
}
*/
//*NumHitCMBSel = numhitcmbsel;
}
void BrLocalTrackTPC::CopyTrahit(BrLocalTrack* OldTr,BrLocalTrack* NewTr,Int_t Imod)
{
//Loop over TrackHits and select BrTrackHit entries Whose TrID matches OldTr
//If Imod of that BrTrackHit also matches the input parameter:
// 1. Get the combined hit associated with that BrTrackHit
// 2. Set its status to appropriate value
// 3. Make a new BrTrackHit and add to TrackHits
// 4. Add this BrTrackHit to NewTr
//
//Arguments:
//BrLocalTrack *OldTr :
//BrLocalTrack *NewTr :
//Int_t Imod :
//
//Input table
//TrackHits (modified in this routine)
Int_t itr;
Int_t NumTraHit;
BrCombinedHit* hitcmb_p;
BrTrackHit* trahit_p;
BrTrackHit* trahit_p1;
NumTraHit = TrackHits->GetEntries();
for(itr=0;itr<NumTraHit;itr++) {
trahit_p1 = (BrTrackHit*)TrackHits->At(itr);
if( trahit_p1->GetTrID() == OldTr ) {
if( trahit_p1->GetImod() != Imod ) {
hitcmb_p = trahit_p1->GetHitcmb();
hitcmb_p->SetStat(TSonataMath::MOD(hitcmb_p->GetStat(),16)+(hitcmb_p->GetStat()/16+1)*16);
BrClonesArray &trahit = *TrackHits;
Int_t ntrahit = TrackHits->GetEntries();
trahit_p = new(trahit[ntrahit]) BrTrackHit();
// This is the new information
trahit_p->SetTrID(NewTr);
trahit_p->SetLoctra(NewTr);
trahit_p->SetHitcmb(hitcmb_p);
// Set up relation with loctra
NewTr->AddTrackHit(trahit_p);
// Copy the rest of the stuff
// trahit_p->gvolu = trahit[itr]->gvolu;
// trahit_p->SetImod(trahit_p1->GetImod());
// trahit_p->SetTra2d(trahit_p1->GetTra2d());
}
}
}
}
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.