/**
* \file TrkLevel2.cpp
* \author Elena Vannuccini
*/
#include <TrkLevel2.h>
#include <iostream>
#include <math.h>
using namespace std;
//......................................
// F77 routines
//......................................
extern "C" {
void dotrack_(int*, double*, double*, double*, double*, int*);
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
// int readb_(const char*);
int readb_();
void mini2_(int*,int*,int*);
void guess_();
}
//--------------------------------------
//
//
//--------------------------------------
TrkTrack::TrkTrack(){
seqno = -1;
image = -1;
chi2 = 0;
nstep = 0;
for(int it1=0;it1<5;it1++){
al[it1] = 0;
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
};
for(int ip=0;ip<6;ip++){
xgood[ip] = 0;
ygood[ip] = 0;
xm[ip] = 0;
ym[ip] = 0;
zm[ip] = 0;
resx[ip] = 0;
resy[ip] = 0;
xv[ip] = 0;
yv[ip] = 0;
zv[ip] = 0;
axv[ip] = 0;
ayv[ip] = 0;
dedx_x[ip] = 0;
dedx_y[ip] = 0;
// clx[ip] = 0;
// cly[ip] = 0;
};
clx = new TRefArray(6,0);
cly = new TRefArray(6,0);
};
//--------------------------------------
//
//
//--------------------------------------
TrkTrack::TrkTrack(const TrkTrack& t){
seqno = t.seqno;
image = t.image;
chi2 = t.chi2;
nstep = t.nstep;
for(int it1=0;it1<5;it1++){
al[it1] = t.al[it1];
for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
};
for(int ip=0;ip<6;ip++){
xgood[ip] = t.xgood[ip];
ygood[ip] = t.ygood[ip];
xm[ip] = t.xm[ip];
ym[ip] = t.ym[ip];
zm[ip] = t.zm[ip];
resx[ip] = t.resx[ip];
resy[ip] = t.resy[ip];
xv[ip] = t.xv[ip];
yv[ip] = t.yv[ip];
zv[ip] = t.zv[ip];
axv[ip] = t.axv[ip];
ayv[ip] = t.ayv[ip];
dedx_x[ip] = t.dedx_x[ip];
dedx_y[ip] = t.dedx_y[ip];
// clx[ip] = 0;//<<<<pointer
// cly[ip] = 0;//<<<<pointer
};
clx = new TRefArray(*(t.clx));
cly = new TRefArray(*(t.cly));
};
//--------------------------------------
//
//
//--------------------------------------
/**
* Evaluates the trajectory in the apparatus associated to the track.
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned.
* @param t pointer to an object of the class Trajectory,
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
* @return error flag.
*/
int TrkTrack::DoTrack(Trajectory* t){
double *dxout = new double[t->npoint];
double *dyout = new double[t->npoint];
double *dzin = new double[t->npoint];
double dal[5];
int ifail = 0;
for (int i=0; i<5; i++) dal[i] = (double)al[i];
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail);
for (int i=0; i<t->npoint; i++){
t->x[i] = (float)*dxout++;
t->y[i] = (float)*dyout++;
}
// delete [] dxout;
// delete [] dyout;
// delete [] dzin;
return ifail;
};
//--------------------------------------
//
//
//--------------------------------------
/**
* Evaluates the trajectory in the apparatus associated to the track.
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned.
* @param t pointer to an object of the class Trajectory,
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
* @return error flag.
*/
int TrkTrack::DoTrack2(Trajectory* t){
double *dxout = new double[t->npoint];
double *dyout = new double[t->npoint];
double *dthxout = new double[t->npoint];
double *dthyout = new double[t->npoint];
double *dtlout = new double[t->npoint];
double *dzin = new double[t->npoint];
double dal[5];
int ifail = 0;
for (int i=0; i<5; i++) dal[i] = (double)al[i];
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
for (int i=0; i<t->npoint; i++){
t->x[i] = (float)*dxout++;
t->y[i] = (float)*dyout++;
t->thx[i] = (float)*dthxout++;
t->thy[i] = (float)*dthyout++;
t->tl[i] = (float)*dtlout++;
}
// delete [] dxout;
// delete [] dyout;
// delete [] dzin;
return ifail;
};
//--------------------------------------
//
//
//--------------------------------------
//float TrkTrack::BdL(){
//};
//--------------------------------------
//
//
//--------------------------------------
Float_t TrkTrack::GetRigidity(){
Float_t rig=0;
if(chi2>0)rig=1./al[4];
if(rig<0) rig=-rig;
return rig;
};
//
Float_t TrkTrack::GetDeflection(){
Float_t def=0;
if(chi2>0)def=al[4];
return def;
};
//
Float_t TrkTrack::GetDEDX(){
Float_t dedx=0;
for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i];
dedx = dedx/(this->GetNX()+this->GetNY());
return dedx;
};
//--------------------------------------
//
//
//--------------------------------------
void TrkTrack::Dump(){
cout << endl << "========== Track " ;
cout << endl << "seq. n. : "<< seqno;
cout << endl << "image n. : "<< image;
cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
cout << endl << "chi^2 : "<< chi2;
cout << endl << "n.step : "<< nstep;
cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << xgood[i] ;
cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << ygood[i] ;
cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
cout << endl << "xv : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
cout << endl << "yv : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
cout << endl << "zv : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
cout << endl << "resx : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
cout << endl << "resy : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
cout << endl;
}
/**
* Set the TrkTrack position measurements
*/
void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
for(int i=0; i<6; i++) xm[i]=*xmeas++;
for(int i=0; i<6; i++) ym[i]=*ymeas++;
for(int i=0; i<6; i++) zm[i]=*zmeas++;
}
/**
* Set the TrkTrack position resolution
*/
void TrkTrack::SetResolution(double *rx, double *ry){
for(int i=0; i<6; i++) resx[i]=*rx++;
for(int i=0; i<6; i++) resy[i]=*ry++;
}
/**
* Set the TrkTrack good measurement
*/
void TrkTrack::SetGood(int *xg, int *yg){
for(int i=0; i<6; i++) xgood[i]=*xg++;
for(int i=0; i<6; i++) ygood[i]=*yg++;
}
/**
* Load the magnetic field
*/
void TrkTrack::LoadField(TString path){
strcpy(path_.path,path.Data());
path_.pathlen = path.Length();
path_.error = 0;
readb_();
};
/**
* Method to fill minimization-routine common
*/
void TrkTrack::FillMiniStruct(cMini2track& track){
for(int i=0; i<6; i++){
track.xgood[i]=xgood[i];
track.ygood[i]=ygood[i];
track.xm[i]=xm[i];
track.ym[i]=ym[i];
track.zm[i]=zm[i];
// --- temporaneo ----------------------------
// andrebbe inserita la dimensione del sensore
float segment = 100.;
track.xm_a[i]=xm[i];
track.xm_b[i]=xm[i];
track.ym_a[i]=ym[i];
track.ym_b[i]=ym[i];
if( xgood[i] && !ygood[i] ){
track.ym_a[i] = track.ym_a[i]+segment;
track.ym_b[i] = track.ym_b[i]-segment;
}else if( !xgood[i] && ygood[i]){
track.xm_a[i] = track.xm_a[i]+segment;
track.xm_b[i] = track.xm_b[i]-segment;
}
// --- temporaneo ----------------------------
track.resx[i]=resx[i];
track.resy[i]=resy[i];
}
for(int i=0; i<5; i++) track.al[i]=al[i];
track.zini = 23.5;
// ZINI = 23.5 !!! it should be the same parameter in all codes
}
/**
* Method to set values from minimization-routine common
*/
void TrkTrack::SetFromMiniStruct(cMini2track *track){
for(int i=0; i<5; i++) {
al[i]=track->al[i];
for(int j=0; j<5; j++) coval[i][j]=track->cov[i][j];
}
chi2 = track->chi2;
nstep = track->nstep;
for(int i=0; i<6; i++){
xv[i] = track->xv[i];
yv[i] = track->yv[i];
zv[i] = track->zv[i];
xm[i] = track->xm[i];
ym[i] = track->ym[i];
zm[i] = track->zm[i];
axv[i] = track->axv[i];
ayv[i] = track->ayv[i];
}
}
/**
* Tracking method. It calls F77 mini routine.
*/
void TrkTrack::Fit(double pfixed, int& fail, int iprint){
float al_ini[] = {0.,0.,0.,0.,0.};
extern cMini2track track_;
fail = 0;
FillMiniStruct(track_);
// if fit variables have been reset, evaluate the initial guess
if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
// --------------------- free momentum
if(pfixed==0.) {
track_.pfixed=0.;
}
// --------------------- fixed momentum
if(pfixed!=0.) {
al[4]=1./pfixed;
track_.pfixed=pfixed;
}
// store temporarily the initial guess
for(int i=0; i<5; i++) al_ini[i]=track_.al[i];
// ------------------------------------------
// call mini routine
int istep=0;
int ifail=0;
mini2_(&istep,&ifail, &iprint);
if(ifail!=0) {
if(iprint)cout << "ERROR: ifail= " << ifail << endl;
fail = 1;
}
// ------------------------------------------
SetFromMiniStruct(&track_);
// cout << endl << "eta ===> " << track_.al[4] << endl;
// for(int i=0; i<5; i++) al[i]=track_.al[i];
// chi2=track_.chi2;
// nstep=track_.nstep;
// for(int i=0; i<6; i++) xv[i]=track_.xv[i];
// for(int i=0; i<6; i++) yv[i]=track_.yv[i];
// for(int i=0; i<6; i++) zv[i]=track_.zv[i];
// for(int i=0; i<6; i++) axv[i]=track_.axv[i];
// for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
// for(int i=0; i<5; i++) {
// for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
// }
if(fail){
if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl;
for(int i=0; i<5; i++) al[i]=al_ini[i];
}
};
/*
* Reset the fit parameters
*/
void TrkTrack::FitReset(){
for(int i=0; i<5; i++) al[i]=-9999.;
chi2=0.;
nstep=0;
for(int i=0; i<6; i++) xv[i]=0.;
for(int i=0; i<6; i++) yv[i]=0.;
for(int i=0; i<6; i++) zv[i]=0.;
for(int i=0; i<6; i++) axv[i]=0.;
for(int i=0; i<6; i++) ayv[i]=0.;
for(int i=0; i<5; i++) {
for(int j=0; j<5; j++) coval[i][j]=0.;
}
}
//--------------------------------------
//
//
//--------------------------------------
void TrkTrack::Clear(){
seqno = -1;
image = -1;
chi2 = 0;
nstep = 0;
for(int it1=0;it1<5;it1++){
al[it1] = 0;
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
};
for(int ip=0;ip<6;ip++){
xgood[ip] = 0;
ygood[ip] = 0;
xm[ip] = 0;
ym[ip] = 0;
zm[ip] = 0;
resx[ip] = 0;
resy[ip] = 0;
xv[ip] = 0;
yv[ip] = 0;
zv[ip] = 0;
axv[ip] = 0;
ayv[ip] = 0;
dedx_x[ip] = 0;
dedx_y[ip] = 0;
// clx[ip] = 0;
// cly[ip] = 0;
};
clx->Clear();
cly->Clear();
};
//--------------------------------------
//
//
//--------------------------------------
void TrkTrack::Delete(){
Clear();
clx->Delete();
cly->Delete();
};
//--------------------------------------
//
//
//--------------------------------------
//--------------------------------------
//
//
//--------------------------------------
TrkSinglet::TrkSinglet(){
plane = 0;
coord[0] = 0;
coord[1] = 0;
sgnl = 0;
cls = 0;
};
//--------------------------------------
//
//
//--------------------------------------
TrkSinglet::TrkSinglet(const TrkSinglet& s){
plane = s.plane;
coord[0] = s.coord[0];
coord[1] = s.coord[1];
sgnl = s.sgnl;
// cls = 0;//<<<<pointer
cls = TRef(s.cls);
};
//--------------------------------------
//
//
//--------------------------------------
void TrkSinglet::Dump(){
int i=0;
cout << endl << "========== Singlet " ;
cout << endl << "plane : " << plane;
cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++;
cout << endl << "sgnl : " << sgnl;
}
//--------------------------------------
//
//
//--------------------------------------
TrkLevel2::TrkLevel2(){
// good2 = -1;
for(Int_t i=0; i<12 ; i++){
// crc[i] = -1;
good[i] = -1;
};
Track = new TClonesArray("TrkTrack");
SingletX = new TClonesArray("TrkSinglet");
SingletY = new TClonesArray("TrkSinglet");
// PhysicalTrack = new TClonesArray("TrkTrack");
//sostituire con TRefArray... appena ho capito come si usa
}
//--------------------------------------
//
//
//--------------------------------------
void TrkLevel2::Dump(){
TClonesArray &t = *Track;
TClonesArray &sx = *SingletX;
TClonesArray &sy = *SingletY;
//
cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
cout << endl << "good : "; for(int i=0; i<12; i++) cout << good[i]<<" ";
cout << endl << "ntrk() : " << this->ntrk() ;
cout << endl << "nclsx() : " << this->nclsx();
cout << endl << "nclsy() : " << this->nclsy();
for(int i=0; i<this->ntrk(); i++) ((TrkTrack *)t[i])->Dump();
for(int i=0; i<this->nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
for(int i=0; i<this->nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
*/
void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){
// temporary objects:
TrkSinglet* t_singlet = new TrkSinglet();
TrkTrack* t_track = new TrkTrack();
// **** general variables ****
// good2 = l2->good2;
for(Int_t i=0; i<12 ; i++){
// crc[i] = l2->crc[i];
good[i] = l2->good[i];
};
// *** TRACKS ***
TClonesArray &t = *Track;
for(int i=0; i<l2->ntrk; i++){
t_track->seqno = i;// NBNBNBNB deve sempre essere = i
t_track->image = l2->image[i]-1;
// cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
t_track->chi2 = l2->chi2_nt[i];
t_track->nstep = l2->nstep_nt[i];
for(int it1=0;it1<5;it1++){
t_track->al[it1] = l2->al_nt[i][it1];
for(int it2=0;it2<5;it2++)
t_track->coval[it1][it2] = l2->coval[i][it2][it1];
};
for(int ip=0;ip<6;ip++){
t_track->xgood[ip] = l2->xgood_nt[i][ip];
t_track->ygood[ip] = l2->ygood_nt[i][ip];
t_track->xm[ip] = l2->xm_nt[i][ip];
t_track->ym[ip] = l2->ym_nt[i][ip];
t_track->zm[ip] = l2->zm_nt[i][ip];
t_track->resx[ip] = l2->resx_nt[i][ip];
t_track->resy[ip] = l2->resy_nt[i][ip];
t_track->xv[ip] = l2->xv_nt[i][ip];
t_track->yv[ip] = l2->yv_nt[i][ip];
t_track->zv[ip] = l2->zv_nt[i][ip];
t_track->axv[ip] = l2->axv_nt[i][ip];
t_track->ayv[ip] = l2->ayv_nt[i][ip];
t_track->dedx_x[ip] = l2->dedx_x[i][ip];
t_track->dedx_y[ip] = l2->dedx_y[i][ip];
// t_track->clx[ip] = 0;
// t_track->cly[ip] = 0;
};
new(t[i]) TrkTrack(*t_track);
t_track->Clear();
};
// *** SINGLETS ***
TClonesArray &sx = *SingletX;
for(int i=0; i<l2->nclsx; i++){
t_singlet->plane = l2->planex[i];
t_singlet->coord[0] = l2->xs[i][0];
t_singlet->coord[1] = l2->xs[i][1];
t_singlet->sgnl = l2->signlxs[i];
// t_singlet->cls = 0;
new(sx[i]) TrkSinglet(*t_singlet);
t_singlet->Clear();
}
TClonesArray &sy = *SingletY;
for(int i=0; i<l2->nclsy; i++){
t_singlet->plane = l2->planey[i];
t_singlet->coord[0] = l2->ys[i][0];
t_singlet->coord[1] = l2->ys[i][1];
t_singlet->sgnl = l2->signlys[i];
// t_singlet->cls = 0;
new(sy[i]) TrkSinglet(*t_singlet);
t_singlet->Clear();
};
delete t_track;
delete t_singlet;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
* Ref to Level1 data (clusters) is also set.
*/
void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){
// temporary objects:
TrkSinglet* t_singlet = new TrkSinglet();
TrkTrack* t_track = new TrkTrack();
// general variables
// good2 = l2->good2;
for(Int_t i=0; i<12 ; i++){
// crc[i] = l2->crc[i];
good[i] = l2->good[i];
};
// *** TRACKS ***
TClonesArray &t = *Track;
for(int i=0; i<l2->ntrk; i++){
t_track->seqno = i;// NBNBNBNB deve sempre essere = i
t_track->image = l2->image[i]-1;
// cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
t_track->chi2 = l2->chi2_nt[i];
t_track->nstep = l2->nstep_nt[i];
for(int it1=0;it1<5;it1++){
t_track->al[it1] = l2->al_nt[i][it1];
for(int it2=0;it2<5;it2++)
t_track->coval[it1][it2] = l2->coval[i][it2][it1];
};
for(int ip=0;ip<6;ip++){
t_track->xgood[ip] = l2->xgood_nt[i][ip];
t_track->ygood[ip] = l2->ygood_nt[i][ip];
t_track->xm[ip] = l2->xm_nt[i][ip];
t_track->ym[ip] = l2->ym_nt[i][ip];
t_track->zm[ip] = l2->zm_nt[i][ip];
t_track->resx[ip] = l2->resx_nt[i][ip];
t_track->resy[ip] = l2->resy_nt[i][ip];
t_track->xv[ip] = l2->xv_nt[i][ip];
t_track->yv[ip] = l2->yv_nt[i][ip];
t_track->zv[ip] = l2->zv_nt[i][ip];
t_track->axv[ip] = l2->axv_nt[i][ip];
t_track->ayv[ip] = l2->ayv_nt[i][ip];
t_track->dedx_x[ip] = l2->dedx_x[i][ip];
t_track->dedx_y[ip] = l2->dedx_y[i][ip];
// cout << "traccia "<<i<<" -- "<< ip << " "<< l2->cltrx[i][ip] <<" "<< l2->cltry[i][ip] <<" "<< t_track->xgood[ip] << t_track->ygood[ip]<<endl;
//-----------------------------------------------------
// t_track->clx[ip] = l1->GetCluster(l2->cltrx[i][ip]-1);
// t_track->cly[ip] = l1->GetCluster(l2->cltry[i][ip]-1);
if(t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
if(t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
// if(t_track->ygood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltry[i][ip] "<<l2->cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<<l1->GetCluster(l2->cltry[i][ip]-1)<<endl;
// if(t_track->xgood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltrx[i][ip] "<<l2->cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<<l1->GetCluster(l2->cltrx[i][ip]-1)<<endl;
//-----------------------------------------------------
};
new(t[i]) TrkTrack(*t_track);
t_track->Clear();
};
// *** SINGLETS ***
TClonesArray &sx = *SingletX;
for(int i=0; i<l2->nclsx; i++){
t_singlet->plane = l2->planex[i];
t_singlet->coord[0] = l2->xs[i][0];
t_singlet->coord[1] = l2->xs[i][1];
t_singlet->sgnl = l2->signlxs[i];
//-----------------------------------------------------
// cout << "singolo x "<<i<<" -- "<< l2->clsx[i] <<endl;
t_singlet->cls = l1->GetCluster(l2->clsx[i]-1);
// cout<<" i "<<i<<" l2->clsx[i] "<<l2->clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<<l1->GetCluster(l2->clsx[i]-1)<<endl;
//-----------------------------------------------------
new(sx[i]) TrkSinglet(*t_singlet);
t_singlet->Clear();
}
TClonesArray &sy = *SingletY;
for(int i=0; i<l2->nclsy; i++){
t_singlet->plane = l2->planey[i];
t_singlet->coord[0] = l2->ys[i][0];
t_singlet->coord[1] = l2->ys[i][1];
t_singlet->sgnl = l2->signlys[i];
//-----------------------------------------------------
// cout << "singolo y "<<i<<" -- "<< l2->clsy[i] <<endl;
t_singlet->cls = l1->GetCluster(l2->clsy[i]-1);
// cout<<" i "<<i<<" l2->clsy[i] "<<l2->clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<<l1->GetCluster(l2->clsy[i]-1)<<endl;
//-----------------------------------------------------
new(sy[i]) TrkSinglet(*t_singlet);
t_singlet->Clear();
};
delete t_track;
delete t_singlet;
}
/**
* Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common).
*/
void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const {
// general variables
// l2->good2 = good2 ;
for(Int_t i=0; i<12 ; i++){
// l2->crc[i] = crc[i];
l2->good[i] = good[i];
};
// *** TRACKS ***
l2->ntrk = Track->GetEntries();
for(Int_t i=0;i<l2->ntrk;i++){
l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image;
l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2;
l2->nstep_nt[i] = ((TrkTrack *)Track->At(i))->nstep;
for(int it1=0;it1<5;it1++){
l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1];
for(int it2=0;it2<5;it2++)
l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2];
};
for(int ip=0;ip<6;ip++){
l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->xgood[ip];
l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->ygood[ip];
l2->xm_nt[i][ip] = ((TrkTrack *)Track->At(i))->xm[ip];
l2->ym_nt[i][ip] = ((TrkTrack *)Track->At(i))->ym[ip];
l2->zm_nt[i][ip] = ((TrkTrack *)Track->At(i))->zm[ip];
l2->resx_nt[i][ip] = ((TrkTrack *)Track->At(i))->resx[ip];
l2->resy_nt[i][ip] = ((TrkTrack *)Track->At(i))->resy[ip];
l2->xv_nt[i][ip] = ((TrkTrack *)Track->At(i))->xv[ip];
l2->yv_nt[i][ip] = ((TrkTrack *)Track->At(i))->yv[ip];
l2->zv_nt[i][ip] = ((TrkTrack *)Track->At(i))->zv[ip];
l2->axv_nt[i][ip] = ((TrkTrack *)Track->At(i))->axv[ip];
l2->ayv_nt[i][ip] = ((TrkTrack *)Track->At(i))->ayv[ip];
l2->dedx_x[i][ip] = ((TrkTrack *)Track->At(i))->dedx_x[ip];
l2->dedx_y[i][ip] = ((TrkTrack *)Track->At(i))->dedx_y[ip];
};
}
// *** SINGLETS ***
l2->nclsx = SingletX->GetEntries();
for(Int_t i=0;i<l2->nclsx;i++){
l2->planex[i] = ((TrkSinglet *)SingletX->At(i))->plane;
l2->xs[i][0] = ((TrkSinglet *)SingletX->At(i))->coord[0];
l2->xs[i][1] = ((TrkSinglet *)SingletX->At(i))->coord[1];
l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl;
}
l2->nclsy = SingletY->GetEntries();
for(Int_t i=0;i<l2->nclsy;i++){
l2->planey[i] = ((TrkSinglet *)SingletY->At(i))->plane;
l2->ys[i][0] = ((TrkSinglet *)SingletY->At(i))->coord[0];
l2->ys[i][1] = ((TrkSinglet *)SingletY->At(i))->coord[1];
l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl;
}
}
//--------------------------------------
//
//
//--------------------------------------
void TrkLevel2::Clear(){
// good2 = -1;
for(Int_t i=0; i<12 ; i++){
// crc[i] = -1;
good[i] = -1;
};
/* Track->RemoveAll();
SingletX->RemoveAll();
SingletY->RemoveAll();*/
// modify to avoid memory leakage
Track->Clear();
SingletX->Clear();
SingletY->Clear();
}
//--------------------------------------
//
//
//--------------------------------------
void TrkLevel2::Delete(){
Clear();
Track->Delete();
SingletX->Delete();
SingletY->Delete();
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Sort physical tracks and stores them in a TObjectArray, ordering by increasing chi**2 value (in case of track image, it selects the one with lower chi**2). The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it).
* This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used.
*/
TRefArray *TrkLevel2::GetTracks_NFitSorted(){
TRefArray *sorted = new TRefArray();
TClonesArray &t = *Track;
// TClonesArray &ts = *PhysicalTrack;
int N = ntrk();
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
// int m[50]; for(int i=0; i<N; i++)m[i]=1;
int indo=0;
int indi=0;
while(N != 0){
int nfit =0;
float chi2ref = numeric_limits<float>::max();
// first loop to search maximum num. of fit points
for(int i=0; i < ntrk(); i++){
if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
nfit = ((TrkTrack *)t[i])->GetNtot();
}
}
//second loop to search minimum chi2 among selected
for(int i=0; i<this->ntrk(); i++){
Float_t chi2 = ((TrkTrack *)t[i])->chi2;
if(chi2 < 0) chi2 = chi2*1000;
if( chi2 < chi2ref
&& ((TrkTrack *)t[i])->GetNtot() == nfit
&& m[i]==1){
chi2ref = ((TrkTrack *)t[i])->chi2;
indi = i;
};
};
if( ((TrkTrack *)t[indi])->HasImage() ){
m[((TrkTrack *)t[indi])->image] = 0;
N--;
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
};
sorted->Add( (TrkTrack*)t[indi] );
m[indi] = 0;
// cout << "SORTED "<< indo << " "<< indi << " "<< N << endl;
N--;
indo++;
}
m.clear();
// cout << "GetTracks_NFitSorted(it): Done"<< endl;
return sorted;
// return PhysicalTrack;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Retrieves the is-th stored track.
* @param it Track number, ranging from 0 to ntrk().
* Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine.
*/
TrkTrack *TrkLevel2::GetStoredTrack(int is){
if(is >= this->ntrk()){
cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl;
cout << " Stored tracks ntrk() = "<< this->ntrk() << endl;
return 0;
}
TClonesArray &t = *(Track);
TrkTrack *track = (TrkTrack*)t[is];
return track;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Retrieves the is-th stored X singlet.
* @param it Singlet number, ranging from 0 to nclsx().
*/
TrkSinglet *TrkLevel2::GetSingletX(int is){
if(is >= this->nclsx()){
cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
cout << " Stored x-singlets nclsx() = "<< this->nclsx() << endl;
return 0;
}
TClonesArray &t = *(SingletX);
TrkSinglet *singlet = (TrkSinglet*)t[is];
return singlet;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Retrieves the is-th stored Y singlet.
* @param it Singlet number, ranging from 0 to nclsx().
*/
TrkSinglet *TrkLevel2::GetSingletY(int is){
if(is >= this->nclsy()){
cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
cout << " Stored y-singlets nclsy() = "<< this->nclsx() << endl;
return 0;
}
TClonesArray &t = *(SingletY);
TrkSinglet *singlet = (TrkSinglet*)t[is];
return singlet;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Retrieves the it-th "physical" track, sorted by the method GetNTracks().
* @param it Track number, ranging from 0 to GetNTracks().
*/
TrkTrack *TrkLevel2::GetTrack(int it){
if(it >= this->GetNTracks()){
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
return 0;
}
TRefArray *sorted = GetTracks(); //TEMPORANEO
TrkTrack *track = (TrkTrack*)sorted->At(it);
sorted->Delete();
return track;
}
/**
* Give the number of "physical" tracks, sorted by the method GetTracks().
*/
Int_t TrkLevel2::GetNTracks(){
Float_t ntot=0;
TClonesArray &t = *Track;
for(int i=0; i<ntrk(); i++) {
if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
else ntot+=0.5;
}
return (Int_t)ntot;
};
//--------------------------------------
//
//
//--------------------------------------
/**
* Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks().
* @param it Track number, ranging from 0 to GetNTracks().
*/
TrkTrack *TrkLevel2::GetTrackImage(int it){
if(it >= this->GetNTracks()){
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
return 0;
}
TRefArray* sorted = GetTracks(); //TEMPORANEO
TrkTrack *track = (TrkTrack*)sorted->At(it);
if(!track->HasImage()){
cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;
return 0;
}
TrkTrack *image = (TrkTrack*)(*Track)[track->image];
sorted->Delete();
return image;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Loads the magnetic field.
* @param s Path of the magnetic-field files.
*/
void TrkLevel2::LoadField(TString path){
//
strcpy(path_.path,path.Data());
path_.pathlen = path.Length();
path_.error = 0;
readb_();
//
};
//--------------------------------------
//
//
//--------------------------------------
/**
* Get tracker-plane (mechanical) z-coordinate
* @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM)
*/
Float_t TrkLevel2::GetZTrk(Int_t plane_id){
switch(plane_id){
case 1: return ZTRK1;
case 2: return ZTRK2;
case 3: return ZTRK3;
case 4: return ZTRK4;
case 5: return ZTRK5;
case 6: return ZTRK6;
default: return 0.;
};
};
//--------------------------------------
//
//
//--------------------------------------
/**
* Trajectory default constructor.
* (By default is created with z-coordinates inside the tracking volume)
*/
Trajectory::Trajectory(){
npoint = 10;
x = new float[npoint];
y = new float[npoint];
z = new float[npoint];
thx = new float[npoint];
thy = new float[npoint];
tl = new float[npoint];
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
for(int i=0; i<npoint; i++){
x[i] = 0;
y[i] = 0;
z[i] = (ZTRK1) - i*dz;
thx[i] = 0;
thy[i] = 0;
tl[i] = 0;
}
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Trajectory constructor.
* (By default is created with z-coordinates inside the tracking volume)
* \param n Number of points
*/
Trajectory::Trajectory(int n){
if(n<=0){
cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
n=10;
}
npoint = n;
x = new float[npoint];
y = new float[npoint];
z = new float[npoint];
thx = new float[npoint];
thy = new float[npoint];
tl = new float[npoint];
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
for(int i=0; i<npoint; i++){
x[i] = 0;
y[i] = 0;
z[i] = (ZTRK1) - i*dz;
thx[i] = 0;
thy[i] = 0;
tl[i] = 0;
}
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Trajectory constructor.
* \param n Number of points
* \param pz Pointer to float array, defining z coordinates
*/
Trajectory::Trajectory(int n, float* zin){
npoint = 10;
if(n>0)npoint = n;
x = new float[npoint];
y = new float[npoint];
z = new float[npoint];
thx = new float[npoint];
thy = new float[npoint];
tl = new float[npoint];
int i=0;
do{
x[i] = 0;
y[i] = 0;
z[i] = zin[i];
thx[i] = 0;
thy[i] = 0;
tl[i] = 0;
i++;
}while(zin[i-1] > zin[i] && i < npoint);
npoint=i;
if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Dump the trajectory coordinates.
*/
void Trajectory::Dump(){
cout <<endl<< "Trajectory ========== "<<endl;
for (int i=0; i<npoint; i++){
cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
cout <<" -- " << thx[i] <<" "<< thy[i] ;
cout <<" -- " << tl[i] << endl;
};
}
//--------------------------------------
//
//
//--------------------------------------
/**
* Get trajectory length between two points
* @param ifirst first point (default 0)
* @param ilast last point (default npoint)
*/
float Trajectory::GetLength(int ifirst, int ilast){
if( ifirst<0 ) ifirst = 0;
if( ilast>=npoint) ilast = npoint-1;
float l=0;
for(int i=ifirst;i<=ilast;i++){
l=l+tl[i];
};
if(z[ilast] > ZINI)l=l-tl[ilast];
if(z[ifirst] < ZINI) l=l-tl[ifirst];
return l;
}
/**
* Evaluates the trajectory in the apparatus associated to the track.
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned.
* @param t pointer to an object of the class Trajectory,
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
* @return error flag.
*/
int Trajectory::DoTrack2(float* al){
double *dxout = new double[npoint];
double *dyout = new double[npoint];
double *dthxout = new double[npoint];
double *dthyout = new double[npoint];
double *dtlout = new double[npoint];
double *dzin = new double[npoint];
double dal[5];
int ifail = 0;
for (int i=0; i<5; i++) dal[i] = (double)al[i];
for (int i=0; i<npoint; i++) dzin[i] = (double)z[i];
dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
for (int i=0; i<npoint; i++){
x[i] = (float)*dxout++;
y[i] = (float)*dyout++;
thx[i] = (float)*dthxout++;
thy[i] = (float)*dthyout++;
tl[i] = (float)*dtlout++;
}
return ifail;
};
ClassImp(TrkLevel2);
ClassImp(TrkSinglet);
ClassImp(TrkTrack);
ClassImp(Trajectory);