--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp 2006/08/04 08:18:06 1.8 +++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp 2014/10/15 08:45:51 1.59 @@ -4,6 +4,7 @@ */ #include #include +#include using namespace std; //...................................... // F77 routines @@ -11,20 +12,29 @@ 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*); + void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*); + void mini2_(int*,int*,int*); + void guess_(); + void gufld_(float*, float*); + float risxeta2_(float *); + float risxeta3_(float *); + float risxeta4_(float *); + float risyeta2_(float *); } + //-------------------------------------- // // //-------------------------------------- TrkTrack::TrkTrack(){ +// cout << "TrkTrack::TrkTrack()" << endl; 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 it2=0;it2<5;it2++)coval[it1][it2] = 0; }; for(int ip=0;ip<6;ip++){ xgood[ip] = 0; @@ -34,6 +44,8 @@ zm[ip] = 0; resx[ip] = 0; resy[ip] = 0; + tailx[ip] = 0; + taily[ip] = 0; xv[ip] = 0; yv[ip] = 0; zv[ip] = 0; @@ -41,7 +53,28 @@ ayv[ip] = 0; dedx_x[ip] = 0; dedx_y[ip] = 0; - }; + multmaxx[ip] = 0; + multmaxy[ip] = 0; + seedx[ip] = 0; + seedy[ip] = 0; + xpu[ip] = 0; + ypu[ip] = 0; + + }; + +// TrkParams::SetTrackingMode(); +// TrkParams::SetPrecisionFactor(); +// TrkParams::SetStepMin(); + TrkParams::SetMiniDefault(); + TrkParams::SetPFA(); + + int ngf = TrkParams::nGF; + for(int i=0; i>> OBSOLETE !!! use TrkTrack::DoTrack(Trajectory* t) instead + * + */ +int TrkTrack::DoTrack2(Trajectory* t){ + + cout << endl; + cout << " int TrkTrack::DoTrack2(Trajectory* t) --->> NB NB !! this method is going to be eliminated !!! "<>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<npoint]; - double *dyout = new double[t->npoint]; - double *dzin = new double[t->npoint]; + 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; @@ -96,16 +219,27 @@ for (int i=0; i<5; i++) dal[i] = (double)al[i]; for (int i=0; inpoint; i++) dzin[i] = (double)t->z[i]; - dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail); + TrkParams::Load(1); + if( !TrkParams::IsLoaded(1) ){ + cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); for (int i=0; inpoint; i++){ - t->x[i] = (float)*dxout++; - t->y[i] = (float)*dyout++; + t->x[i] = (float)*(dxout+i); + t->y[i] = (float)*(dyout+i); + t->thx[i] = (float)*(dthxout+i); + t->thy[i] = (float)*(dthyout+i); + t->tl[i] = (float)*(dtlout+i); } -// delete [] dxout; -// delete [] dyout; -// delete [] dzin; + delete [] dxout; + delete [] dyout; + delete [] dzin; + delete [] dthxout; + delete [] dthyout; + delete [] dtlout; return ifail; }; @@ -113,109 +247,1286 @@ // // //-------------------------------------- +//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; +}; +// /** - * 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. + * Method to retrieve the dE/dx measured on a tracker view. + * @param ip plane (0-5) + * @param iv view (0=x 1=y) */ -int TrkTrack::DoTrack2(Trajectory* t){ +Float_t TrkTrack::GetDEDX(int ip, int iv){ + if(iv==0 && ip>=0 && ip<6)return fabs(dedx_x[ip]); + else if(iv==1 && ip>=0 && ip<6)return fabs(dedx_y[ip]); + else { + cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<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]; +/** + * Method to evaluate the dE/dx averaged over all planes. + */ +Float_t TrkTrack::GetDEDX(){ + Float_t dedx=0; + for(Int_t ip=0; ip<6; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip); + dedx = dedx/(GetNX()+GetNY()); + return dedx; +}; +/** + * Returns 1 if the cluster on a tracker view includes bad strips + * (at least one bad strip among the four strip used by p.f.a.) + * @param ip plane (0-5) + * @param iv view (0=x 1=y) + */ +Bool_t TrkTrack::IsBad(int ip,int iv){ + if(iv==0 && ip>=0 && ip<6)return (xgood[ip]<0) ; + else if(iv==1 && ip>=0 && ip<6)return (ygood[ip]<0) ; + else { + cout << "TrkTrack::IsBad(int ip, int iv) -- wrong input parameters "<=0 && ip<6)return (dedx_x[ip]<0) ; + else if(iv==1 && ip>=0 && ip<6)return (dedx_y[ip]<0) ; + else { + cout << "TrkTrack::IsSaturated(int ip, int iv) -- wrong input parameters "<3)chiq=chiq/(GetNX()-3); + else chiq=0; + // if(chiq==0)cout << " Float_t TrkTrack::GetChi2X() -- WARNING -- value not defined "<2)chiq=chiq/(GetNY()-2); + else chiq=0; + // if(chiq==0)cout << " Float_t TrkTrack::GetChi2Y() -- WARNING -- value not defined "<3)lnl=lnl/(GetNX()-3); + else lnl=0; + if(lnl==0){ + cout << " Float_t TrkTrack::GetLnLX() -- WARNING -- value not defined "<2)lnl=lnl/(GetNY()-2); + else lnl=0; + if(lnl==0){ + cout << " Float_t TrkTrack::GetLnLY() -- WARNING -- value not defined "<5){ + cout << "Float_t TrkTrack::GetEffectiveAngle(int "< wrong input"< wrong input"<10) index=10; + tailx[i]=tx[index]; + if(flag==1) { + if(fabs(axv[i])<=10.) fact = resx[i]/risxeta2_(&(axv[i])); + if(fabs(axv[i])>10.&&fabs(axv[i])<=15.) fact = resx[i]/risxeta3_(&(axv[i])); + if(fabs(axv[i])>15.) fact = resx[i]/risxeta4_(&(axv[i])); + } else fact = 1.; + resx[i] = sx[index]*fact; + } + for(int i=0; i<6; i++) { + index = int((fabs(ayv[i])+1.)/2.); + if(index>10) index=10; + taily[i]=ty[index]; + if(flag==1) fact = resy[i]/risyeta2_(&(ayv[i])); + else fact = 1.; + resy[i] = sy[index]*fact; + } +} +/** + * 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_(); + +// TrkParams::SetTrackingMode(); +// TrkParams::SetPrecisionFactor(); +// TrkParams::SetStepMin(); + TrkParams::SetMiniDefault(); + + TrkParams::Set(path,1); + TrkParams::Load(1); + if( !TrkParams::IsLoaded(1) ){ + cout << "void TrkTrack::LoadField(TString path) --- ERROR --- m.field not loaded"< 1 || ip < 0 || ip > 5 || il < 0 || il > 2) && + true){ + // se il piano risulta colpito, ladder e sensore devono essere + // assegnati correttamente + cout << " void TrkTrack::FillMiniStruct(cMini2track&) --- WARNING --- sensor not defined, cannot read alignment parameters "<al[i]); + for(int j=0; j<5; j++) coval[i][j]= (float) (track->cov[i][j]); + } + chi2 = (float) (track->chi2); + nstep = (float) (track->nstep); + for(int i=0; i<6; i++){ + xv[i] = (float) (track->xv[i]); + yv[i] = (float) (track->yv[i]); + zv[i] = (float) (track->zv[i]); + xm[i] = (float) (track->xm[i]); + ym[i] = (float) (track->ym[i]); + zm[i] = (float) (track->zm[i]); + axv[i] = (float) (track->axv[i]); + ayv[i] = (float) (track->ayv[i]); + resx[i] = (float) (track->resx[i]); //Elena 10th + resy[i] = (float) (track->resy[i]); + } + +} +/** + * \brief Method to re-evaluate coordinates of clusters associated with a track. + * + * The method can be applied only after recovering level1 information + * (either by reprocessing single events from level0 or from + * the TrkLevel1 branch, if present); it calls F77 subroutines that + * read the level1 common and fill the minimization-routine common. + * Some clusters can be excluded or added by means of the methods: + * + * TrkTrack::ResetXGood(int ip) + * TrkTrack::ResetYGood(int ip) + * TrkTrack::SetXGood(int ip, int cid, int is) + * TrkTrack::SetYGood(int ip, int cid, int is) + * + * NB! The method TrkTrack::SetGood(int *xg, int *yg) set the plane-mask (0-1) + * for the minimization-routine common. It deletes the cluster information + * (at least for the moment...) thus cannot be applied before + * TrkTrack::EvaluateClusterPositions(). + * + * Different p.f.a. can be applied by calling (once) the method: + * + * TrkParams::SetPFA(0); //Set ETA p.f.a. + * + * @see TrkParams::SetPFA(int) + */ +Bool_t TrkTrack::EvaluateClusterPositions(){ + +// cout << "void TrkTrack::GetClusterositions() "<GetTrkLevel0() )return false; + * event->GetTrkLevel0()->ProcessEvent(); // re-processing level0->level1 + * int fail=0; + * event->GetTrkLevel2()->GetTrack(0)->Fit(0.,fail,0,1); + * + * @see EvaluateClusterPositions() + * + * The fitting procedure can be varied by changing the tracking mode, + * the fit-precision factor, the minimum number of step, etc. + * @see SetTrackingMode(int) + * @see SetPrecisionFactor(double) + * @see SetStepMin(int) + * @see SetDeltaB(int,double) + */ +void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){ + + TrkParams::Load(1); + if( !TrkParams::IsLoaded(1) ){ + cout << "void TrkTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<>>> fit failed "<0 define an inner fiducial volume) + */ +Bool_t TrkTrack::IsInsideCavity(float toll){ + +// float xmagntop, ymagntop, xmagnbottom, ymagnbottom; +// xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.); +// ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.); +// xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.); +// ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.); +// if( xmagntop>XMAGNLOW && xmagntopYMAGNLOW && ymagntopXMAGNLOW && xmagnbottomYMAGNLOW && ymagnbottom= TrkParams::xGF_max[i] - toll || + yGF[i] <= TrkParams::yGF_min[i] + toll || + yGF[i] >= TrkParams::yGF_max[i] - toll || + false){ + + return false; + } + } + return true; + + +} +/** + * Returns true if the track is inside the nominal acceptance, which is defined + * by the intersection among magnet cavity, silicon-plane sensitive area and + * ToF sensitive area (nominal values from the official document used to + * calculate the geometrical factor) + * @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume) + */ +// Bool_t TrkTrack::IsInsideAcceptance(){ + +// int ngf = TrkParams::nGF; +// for(int i=0; i= TrkParams::xGF_max[i] || +// yGF[i] <= TrkParams::yGF_min[i] || +// yGF[i] >= TrkParams::yGF_max[i] || +// false)return false; +// } +// return true; + +// } +Bool_t TrkTrack::IsInsideAcceptance(float toll){ + + + int ngf = TrkParams::nGF; + for(int i=0; i= TrkParams::xGF_max[i] - toll || + yGF[i] <= TrkParams::yGF_min[i] + toll || + yGF[i] >= TrkParams::yGF_max[i] - toll || + false){ + + return false; + } + } + return true; +} + +/** + * Returns true if the track is inside one of the surfaces which define the + * geometrical acceptance. + * @param surf tag of the surface (possible values are: S11 S12 S21 S22 T1 + * CUF T2 T3 T4 T5 CLF T6 S31 S32). + * @param toll Tolerance around the nominal surface (toll>0 define an inner + * fiducial surface) +*/ +Bool_t TrkTrack::IsInsideGFSurface(const char* surf, float toll){ + + + int ngf = TrkParams::nGF; + bool SURFOK = false; + for(int i=0; i TrkParams::xGF_min[i] + toll && + xGF[i] < TrkParams::xGF_max[i] - toll && + yGF[i] > TrkParams::yGF_min[i] + toll && + yGF[i] < TrkParams::yGF_max[i] - toll && + true)return true; + } + } + if( !SURFOK )cout << " Bool_t TrkTrack::IsInsideGFSurface(char* surf, float toll) --> suface "<5||clid<1||il<-1||il>2||is<-1||is>1) + cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<5||clid<1||il<-1||il>2||is<-1||is>1) + cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<, evaluated from the first to the last hit x view. + */ +Float_t TrkTrack::GetXav(){ + + int first_plane = -1; + int last_plane = -1; + for(Int_t ip=0; ip<6; ip++){ + if( XGood(ip) && first_plane == -1 )first_plane = ip; + if( XGood(ip) && first_plane != -1 )last_plane = ip; + } + if( first_plane == -1 || last_plane == -1){ + return -100; + } + if( last_plane-first_plane+1 ==0 )return -100; + + Float_t av = 0; + for(int ip=first_plane; ip<=last_plane; ip++)av+=xv[ip]; + + return (av/(last_plane-first_plane+1)); +} +/** + * \brief Average Y + * Average value of , evaluated from the first to the last hit x view. + */ +Float_t TrkTrack::GetYav(){ + + int first_plane = -1; + int last_plane = -1; + for(Int_t ip=0; ip<6; ip++){ + if( XGood(ip) && first_plane == -1 )first_plane = ip; + if( XGood(ip) && first_plane != -1 )last_plane = ip; + } + if( first_plane == -1 || last_plane == -1){ + return -100; + } + if( last_plane-first_plane+1 ==0 )return -100; + + Float_t av = 0; + for(int ip=first_plane; ip<=last_plane; ip++)av+=yv[ip]; + + return (av/(last_plane-first_plane+1)); +} +/** + * \brief Average Z + * Average value of , evaluated from the first to the last hit x view. + */ +Float_t TrkTrack::GetZav(){ + + int first_plane = -1; + int last_plane = -1; + for(Int_t ip=0; ip<6; ip++){ + if( XGood(ip) && first_plane == -1 )first_plane = ip; + if( XGood(ip) && first_plane != -1 )last_plane = ip; + } + if( first_plane == -1 || last_plane == -1){ + return -100; + } + if( last_plane-first_plane+1 ==0 )return -100; + + Float_t av = 0; + for(int ip=first_plane; ip<=last_plane; ip++)av+=zv[ip]; + + return (av/(last_plane-first_plane+1)); +} + +/** + * \brief Number of column traversed + */ +Int_t TrkTrack::GetNColumns(){ + int sensors[] = {0,0,0,0,0,0}; + for(int ip=0; ip<6; ip++){ + int sensorid = GetLadder(ip)+3*GetSensor(ip); + if(XGood(ip)||YGood(ip)) + if(sensorid>=0 && sensorid<6)sensors[sensorid]=1; + } + int nsensors=0; + for(int is=0; is<6; is++)nsensors += sensors[is]; + return nsensors; +}; +/** + * \brief Give the maximum energy release + */ +Float_t TrkTrack::GetDEDX_max(int ip, int iv){ + Float_t max=0; + int pfrom = 0; + int pto = 6; + int vfrom = 0; + int vto = 2; + if(ip>=0&&ip<6){ + pfrom = ip; + pto = ip+1; + } + if(iv>=0&&iv<2){ + vfrom = iv; + vto = iv+1; + } + for(int i=pfrom; imax)max=GetDEDX(i,j); + if(j==1 && YGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); + } + return max; + +}; + +/** + * \brief Give the minimum energy release + */ +Float_t TrkTrack::GetDEDX_min(int ip, int iv){ + Float_t min=100000000; + int pfrom = 0; + int pto = 6; + int vfrom = 0; + int vto = 2; + if(ip>=0&&ip<6){ + pfrom = ip; + pto = ip+1; + } + if(iv>=0&&iv<2){ + vfrom = iv; + vto = iv+1; + } + for(int i=pfrom; i=0&&ip<6){ + pfrom = ip; + pto = ip+1; + } + if(iv>=0&&iv<2){ + vfrom = iv; + vto = iv+1; + } + for(int i=pfrom; ifabs(max))max=xm[i]-xv[i]; + if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=ym[i]-yv[i]; + } + } + return max; + +}; +/** + * \brief Give the anerage spatial residual + */ +Float_t TrkTrack::GetResidual_av(int ip, int iv){ + // +//Sum$((xm>-50)*(xm-xv)/resx)/sqrt(TrkTrack.GetNX()*TrkTrack.GetChi2X())<0.3 + + Float_t av = 0.; + int nav = 0; + // + int pfrom = 0; + int pto = 6; + int vfrom = 0; + int vto = 2; + if(ip>=0&&ip<6){ + pfrom = ip; + pto = ip+1; + } + if(iv>=0&&iv<2){ + vfrom = iv; + vto = iv+1; + } + for(int i=pfrom; inpoint; i++) dzin[i] = (double)t->z[i]; +}; - dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); - - for (int i=0; inpoint; 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; +/** + * \brief Give the maximum multiplicity on the x view + */ +Int_t TrkTrack::GetClusterX_Multiplicity_max(){ + int max=0; + for(int ip=0; ip<6; ip++) + if(GetClusterX_Multiplicity(ip)>max)max=GetClusterX_Multiplicity(ip); + return max; +}; +/** + * \brief Give the minimum multiplicity on the x view + */ +Int_t TrkTrack::GetClusterX_Multiplicity_min(){ + int min=50; + for(int ip=0; ip<6; ip++) + if(GetClusterX_Multiplicity(ip)max)max=GetClusterY_Multiplicity(ip); + return max; +}; +/** + * \brief Give the minimum multiplicity on the x view + */ +Int_t TrkTrack::GetClusterY_Multiplicity_min(){ + int min=50; + for(int ip=0; ip<6; ip++) + if(GetClusterY_Multiplicity(ip)Clear(); +// if(cly)cly->Clear(); +// clx.Clear(); +// cly.Clear(); +}; //-------------------------------------- // // //-------------------------------------- -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::Delete(){ +// cout << "TrkTrack::Delete()"<GetClass()->GetName(),"TrkTrack")==0){ + if(Track)Track->Clear("C"); + Track = track; + } } //-------------------------------------- // // //-------------------------------------- void TrkLevel2::Dump(){ - TClonesArray &t = *Track; - TClonesArray &sx = *SingletX; - TClonesArray &sy = *SingletY; - + + // cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; - cout << endl << "good2 : " << good2; - cout << endl << "crc : "; for(int i=0; i<12; i++) cout << crc[i]; - cout << endl << "ntrk() : " << this->ntrk() ; - cout << endl << "nclsx() : " << this->nclsx(); - cout << endl << "nclsy() : " << this->nclsy(); - for(int i=0; intrk(); i++) ((TrkTrack *)t[i])->Dump(); - for(int i=0; inclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); - for(int i=0; inclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); + cout << endl << "good : "; for(int i=0; i<12; i++) cout << hex <<" 0x"<< good[i]<Dump(); + } +// if(SingletX){ +// TClonesArray &sx = *SingletX; +// for(int i=0; iDump(); +// } +// if(SingletY){ +// TClonesArray &sy = *SingletY; +// for(int i=0; iDump(); +// } + cout << endl; +} +/** + * \brief Dump processing status + */ +void TrkLevel2::StatusDump(int view){ + cout << "DSP n. "<= 12)return false; + return !(good[view]&flagmask); + +}; + + +//-------------------------------------- +// +// +//-------------------------------------- +/** + * The method returns false if the viking-chip was masked + * either apriori ,on the basis of the mask read from the DB, + * or run-by-run, on the basis of the calibration parameters) + * @param iv Tracker view (0-11) + * @param ivk Viking-chip number (0-23) + */ +Bool_t TrkLevel2::GetVKMask(int iv, int ivk){ + Int_t whichbit = (Int_t)pow(2,ivk); + return (whichbit&VKmask[iv])!=0; } +/** + * The method returns false if the viking-chip was masked + * for this event due to common-noise computation failure. + * @param iv Tracker view (0-11) + * @param ivk Viking-chip number (0-23) + */ +Bool_t TrkLevel2::GetVKFlag(int iv, int ivk){ + Int_t whichbit = (Int_t)pow(2,ivk); + return (whichbit&VKflag[iv])!=0; +} +/** + * The method returns true if the viking-chip was masked, either + * forced (see TrkLevel2::GetVKMask(int,int)) or + * for this event only (TrkLevel2::GetVKFlag(int,int)). + * @param iv Tracker view (0-11) + * @param ivk Viking-chip number (0-23) + */ +Bool_t TrkLevel2::IsMaskedVK(int iv, int ivk){ + return !(GetVKMask(iv,ivk)&&GetVKFlag(iv,ivk) ); +}; + //-------------------------------------- // // //-------------------------------------- /** - * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). + * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). + * Ref to Level1 data (clusters) is also set. If l1==NULL no references are set. + * (NB It make sense to set references only if events are stored in a tree that contains also the Level1 branch) */ -void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){ - // -// Track = new TClonesArray("TrkTrack"); -// SingletX = new TClonesArray("TrkSinglet"); -// SingletY = new TClonesArray("TrkSinglet"); +void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){ + +// cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<good2; +// ----------------- for(Int_t i=0; i<12 ; i++){ - crc[i] = l2->crc[i]; + good[i] = l2->good[i]; + VKmask[i]=0; + VKflag[i]=0; + for(Int_t ii=0; ii<24 ; ii++){ + Int_t setbit = (Int_t)pow(2,ii); + if( l2->vkflag[ii][i]!=-1 )VKmask[i]=VKmask[i]|setbit; + if( l2->vkflag[ii][i]!=0 )VKflag[i]=VKflag[i]|setbit; + }; }; +// -------------- // *** TRACKS *** +// -------------- + if(!Track) Track = new TClonesArray("TrkTrack"); TClonesArray &t = *Track; + for(int i=0; intrk; i++){ - t_track->seqno = i; + t_track->seqno = i;// NBNBNBNB deve sempre essere = i t_track->image = l2->image[i]-1; -// cout << "track "<seqno << t_track->image<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]; + // --------------------------------- + // new implementation of xgood/ygood + // --------------------------------- + t_track->xgood[ip] = l2->cltrx[i][ip]; //cluster ID + t_track->ygood[ip] = l2->cltry[i][ip]; //cluster ID + t_track->xgood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor + t_track->ygood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor + if(l2->xbad[i][ip]>0)t_track->xgood[ip]=-t_track->xgood[ip]; + if(l2->ybad[i][ip]>0)t_track->ygood[ip]=-t_track->ygood[ip]; +// if(l2->xbad[i][ip]>0 || l2->ybad[i][ip]>0){ +// if(l2->dedx_x[i][ip]<0 || l2->dedx_y[i][ip]<0){ +// cout << ip << " - "<< l2->cltrx[i][ip] << " "<cltry[i][ip]<<" "<ls[i][ip]<xgood[ip]<<" "<ygood[ip]<GetClusterX_ID(ip)<<" "<GetClusterY_ID(ip)<<" "<GetLadder(ip)<<" "<GetSensor(ip)<BadClusterX(ip)<<" "<BadClusterY(ip)<SaturatedClusterX(ip)<<" "<SaturatedClusterY(ip)<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->tailx[ip] = l2->tailx[i][ip]; + t_track->taily[ip] = l2->taily[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]; @@ -310,32 +1758,72 @@ 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->multmaxx[ip] = l2->multmaxx[i][ip]; + t_track->multmaxy[ip] = l2->multmaxy[i][ip]; + t_track->seedx[ip] = l2->seedx[i][ip]; + t_track->seedy[ip] = l2->seedy[i][ip]; + t_track->xpu[ip] = l2->xpu[i][ip]; + t_track->ypu[ip] = l2->ypu[i][ip]; + //----------------------------------------------------- + //----------------------------------------------------- + //----------------------------------------------------- + //----------------------------------------------------- }; + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // evaluated coordinates (to define GF) + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + int ngf = TrkParams::nGF; + float *zgf = TrkParams::zGF; + Trajectory tgf = Trajectory(ngf,zgf); + tgf.DoTrack(t_track->al);//<<<< integrate the trajectory + for(int ip=0; ipxGF[ip] = tgf.x[ip]; + t_track->yGF[ip] = tgf.y[ip]; + } + +// if(t_track->IsSaturated())t_track->Dump(); new(t[i]) TrkTrack(*t_track); t_track->Clear(); - }; + };//end loop over track + +// ---------------- // *** SINGLETS *** +// ---------------- + if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); TClonesArray &sx = *SingletX; for(int i=0; inclsx; 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->multmax = l2->multmaxsx[i]; + if(l2->sxbad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; + //----------------------------------------------------- +// if(l1) t_singlet->cls = l1->GetCluster(l2->clsx[i]-1); + //----------------------------------------------------- new(sx[i]) TrkSinglet(*t_singlet); t_singlet->Clear(); } + if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); TClonesArray &sy = *SingletY; for(int i=0; inclsy; 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->multmax = l2->multmaxsy[i]; + if(l2->sybad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; + //----------------------------------------------------- +// if(l1) t_singlet->cls = l1->GetCluster(l2->clsy[i]-1); + //----------------------------------------------------- new(sy[i]) TrkSinglet(*t_singlet); t_singlet->Clear(); - }; + }; + + - delete t_track; - delete t_singlet; + delete t_track; + delete t_singlet; } /** * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common). @@ -344,53 +1832,63 @@ void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const { // general variables - l2->good2 = good2 ; +// l2->good2 = good2 ; for(Int_t i=0; i<12 ; i++){ - l2->crc[i] = crc[i]; +// l2->crc[i] = crc[i]; + l2->good[i] = good[i]; }; // *** TRACKS *** - l2->ntrk = Track->GetEntries(); - for(Int_t i=0;intrk;i++){ - l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image; - l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2; - 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]; - }; + if(Track){ + l2->ntrk = Track->GetEntries(); + for(Int_t i=0;intrk;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->tailx[i][ip] = ((TrkTrack *)Track->At(i))->tailx[ip]; + l2->taily[i][ip] = ((TrkTrack *)Track->At(i))->taily[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;inclsx;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;inclsy;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; + if(SingletX){ + l2->nclsx = SingletX->GetEntries(); + for(Int_t i=0;inclsx;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; + } + } + + if(SingletY){ + l2->nclsy = SingletY->GetEntries(); + for(Int_t i=0;inclsy;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; + } } } //-------------------------------------- @@ -398,17 +1896,30 @@ // //-------------------------------------- void TrkLevel2::Clear(){ - good2 = -1; for(Int_t i=0; i<12 ; i++){ - crc[i] = -1; + good[i] = -1; + VKflag[i] = 0; + VKmask[i] = 0; }; -/* Track->RemoveAll(); - SingletX->RemoveAll(); - SingletY->RemoveAll();*/ - // modify to avoid memory leakage - Track->Clear(); - SingletX->Clear(); - SingletY->Clear(); +// if(Track)Track->Clear("C"); +// if(SingletX)SingletX->Clear("C"); +// if(SingletY)SingletY->Clear("C"); + if(Track)Track->Delete(); + if(SingletX)SingletX->Delete(); + if(SingletY)SingletY->Delete(); +} +// //-------------------------------------- +// // +// // +// //-------------------------------------- +void TrkLevel2::Delete(){ + +// cout << "void TrkLevel2::Delete()"<ntrk(); + // TRefArray *sorted = new TRefArray(); + TRefArray *sorted = NULL; + + TClonesArray &t = *Track; +// TClonesArray &ts = *PhysicalTrack; + int N = ntrk(); vector m(N); for(int i=0; intrk(); i++){ - if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){ - chi2ref = ((TrkTrack *)t[i])->chi2; - indi = i; + while(N > 0){ +// while(N != 0){ + int nfit =0; + float chi2ref = numeric_limits::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(); } } - if( ((TrkTrack *)t[indi])->image != -1 ){ + //second loop to search minimum chi2 among selected + for(int i=0; ichi2; + 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--; - } - new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); + +// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<Add( (TrkTrack*)t[indi] ); + m[indi] = 0; +// cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "< m(N); for(int i=0; i::max(); - // first loop to search maximum num. of fit points - for(int i=0; i < ntrk(); i++){ -// if(N==ntrk())cout << "** "<GetImageSeqNo()<< " " <<((TrkTrack *)t[i])->GetNtot() << " " <<((TrkTrack *)t[i])->chi2 << endl; - 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; intrk(); i++){ - if( ((TrkTrack *)t[i])->chi2 < chi2ref - && ((TrkTrack *)t[i])->GetNtot()== nfit - && m[i]==1){ - chi2ref = ((TrkTrack *)t[i])->chi2; - indi = i; - // cout << "2** "<HasImage() ){ - m[((TrkTrack *)t[indi])->image] = 0; - N--; - - // Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->GetNtot(); - // Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->chi2; - - // cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "< m(N); for(int i=0; i::max(); - // first loop to search maximum num. of fit points - for(int i=0; i < ntrk(); i++){ -// if(N==ntrk())cout << "** "<GetImageSeqNo()<< " " <<((TrkTrack *)t[i])->GetNtot() << " " <<((TrkTrack *)t[i])->chi2 << endl; - 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; intrk(); i++){ - if( ((TrkTrack *)t[i])->chi2 < chi2ref - && ((TrkTrack *)t[i])->GetNtot()== nfit - && m[i]==1){ - chi2ref = ((TrkTrack *)t[i])->chi2; - indi = i; - // cout << "2** "<HasImage() ){ - m[((TrkTrack *)t[indi])->image] = 0; - N--; - - // Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->GetNtot(); - // Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->chi2; - - // cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<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; } //-------------------------------------- // @@ -575,10 +1998,13 @@ TrkTrack *TrkLevel2::GetStoredTrack(int is){ if(is >= this->ntrk()){ - cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl; - cout << " Stored tracks ntrk() = "<< this->ntrk() << endl; + cout << "TrkTrack *TrkLevel2::GetStoredTrack(int) >> Track "<< is << "doen not exits! " << endl; + cout << "Stored tracks ntrk() = "<< this->ntrk() << endl; return 0; } + if(!Track){ + cout << "TrkTrack *TrkLevel2::GetStoredTrack(int is) >> (TClonesArray*) Track ==0 "<= this->nclsx()){ - cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl; - cout << " Stored x-singlets nclsx() = "<< this->nclsx() << endl; + cout << "TrkSinglet *TrkLevel2::GetSingletX(int) >> Singlet "<< is << "doen not exits! " << endl; + cout << "Stored x-singlets nclsx() = "<< this->nclsx() << endl; return 0; } + if(!SingletX)return 0; TClonesArray &t = *(SingletX); TrkSinglet *singlet = (TrkSinglet*)t[is]; return singlet; @@ -613,10 +2040,11 @@ 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; + cout << "TrkSinglet *TrkLevel2::GetSingletY(int) >> Singlet "<< is << "doen not exits! " << endl; + cout << "Stored y-singlets nclsx() = "<< this->nclsx() << endl; return 0; } + if(!SingletY)return 0; TClonesArray &t = *(SingletY); TrkSinglet *singlet = (TrkSinglet*)t[is]; return singlet; @@ -629,32 +2057,20 @@ * 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; - } - TClonesArray* sorted = GetTracks(); - TClonesArray &ts = *sorted; - -// TrkTrack *track = (TrkTrack*)ts[it]; - TrkTrack *track = new TrkTrack( *(TrkTrack*)ts[it] ); //copia - sorted->Delete("all");////TEMPORANEO - return track; -}*/ + TrkTrack *TrkLevel2::GetTrack(int it){ if(it >= this->GetNTracks()){ - cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; - cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; + cout << "TrkTrack *TrkLevel2::GetTrack(int) >> Track "<< it << "does not exits! " << endl; + cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; return 0; } TRefArray *sorted = GetTracks(); //TEMPORANEO + if(!sorted)return 0; TrkTrack *track = (TrkTrack*)sorted->At(it); - sorted->Delete(); + sorted->Clear(); + delete sorted; return track; } /** @@ -663,8 +2079,9 @@ Int_t TrkLevel2::GetNTracks(){ Float_t ntot=0; + if(!Track)return 0; TClonesArray &t = *Track; - for(int i=0; iGetImageSeqNo() == -1 ) ntot+=1.; else ntot+=0.5; } @@ -679,51 +2096,29 @@ * 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){ +TrkTrack *TrkLevel2::GetTrackImage(int it){ if(it >= this->GetNTracks()){ - cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; - cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; + cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not exits! " << endl; + cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; return 0; } - TClonesArray* sorted = GetTracks(); - TClonesArray &ts = *sorted; - TrkTrack *track = (TrkTrack*)ts[it]; -// TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it]; - - if(!track->HasImage()){ - cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; + TRefArray* sorted = GetTracks(); //TEMPORANEO + if(!sorted)return 0; + TrkTrack *track = (TrkTrack*)sorted->At(it); + + if(!track->HasImage()){ + cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not have image! " << endl; return 0; } + if(!Track)return 0; TrkTrack *image = (TrkTrack*)(*Track)[track->image]; -// GetTracks()->Delete(); ////TEMPORANEO - sorted->Delete(); ////TEMPORANEO - - return image; - -}*/ -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(); + delete sorted; - sorted->Delete(); - - return image; + return image; } //-------------------------------------- @@ -734,9 +2129,53 @@ * Loads the magnetic field. * @param s Path of the magnetic-field files. */ -void TrkLevel2::LoadField(TString s){ - readb_(s.Data()); +void TrkLevel2::LoadField(TString path){ +// +// strcpy(path_.path,path.Data()); +// path_.pathlen = path.Length(); +// path_.error = 0; +// readb_(); + +// TrkParams::SetTrackingMode(); +// TrkParams::SetPrecisionFactor(); +// TrkParams::SetStepMin(); + TrkParams::SetMiniDefault(); + + TrkParams::Set(path,1); + TrkParams::Load(1); + if( !TrkParams::IsLoaded(1) ){ + cout << "void TrkLevel2::LoadField(TString path) --- ERROR --- m.field not loaded"< zin[i] && i < npoint); npoint=i; - if(npoint != n)cout << "NB! Trajectory created with "<>> OBSOLETE !!! use Trajectory::DoTrack(float* al, float zini) instead + * + */ +int Trajectory::DoTrack2(float* al, float zini){ + + cout << endl; + cout << " int Trajectory::DoTrack2(float* al, float zini) --->> NB NB !! this method is going to be eliminated !!! "<>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<