--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp 2006/09/05 12:52:20 1.9 +++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp 2006/11/07 15:55:11 1.15 @@ -4,6 +4,7 @@ */ #include #include +#include using namespace std; //...................................... // F77 routines @@ -11,7 +12,9 @@ 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_(const char*); + void mini2_(int*,int*,int*); + void guess_(); } //-------------------------------------- // @@ -21,9 +24,10 @@ 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; + al[it1] = 0; + for(int it2=0;it2<5;it2++)coval[it1][it2] = 0; }; for(int ip=0;ip<6;ip++){ xgood[ip] = 0; @@ -54,9 +58,10 @@ 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]; + 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]; @@ -194,16 +199,230 @@ //-------------------------------------- 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 s){ + readb_(s.Data()); +}; +/** + * 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; +// extern cMini2fitinfo fit_info_; +// extern void mini_2_(int*,int*); +// track_.xm[0]=1.0; + + cout << "fitting routine "< " << 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){ + cout << " >>>> fit failed >>>> drawing initial par"<Clear(); + cly->Clear(); +}; +//-------------------------------------- +// +// +//-------------------------------------- +void TrkTrack::Delete(){ + Clear(); + clx->Delete(); + cly->Delete(); +}; + //-------------------------------------- +// +// +//-------------------------------------- + //-------------------------------------- // // @@ -243,10 +462,11 @@ // //-------------------------------------- TrkLevel2::TrkLevel2(){ - good2 = -1; +// good2 = -1; for(Int_t i=0; i<12 ; i++){ - crc[i] = -1; - }; +// crc[i] = -1; + good[i] = -1; + }; Track = new TClonesArray("TrkTrack"); SingletX = new TClonesArray("TrkSinglet"); SingletY = new TClonesArray("TrkSinglet"); @@ -259,13 +479,13 @@ // //-------------------------------------- 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 << "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(); @@ -281,25 +501,25 @@ * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). */ void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){ - // -// Track = new TClonesArray("TrkTrack"); -// SingletX = new TClonesArray("TrkSinglet"); -// SingletY = new TClonesArray("TrkSinglet"); -// temporary objects: + + // temporary objects: TrkSinglet* t_singlet = new TrkSinglet(); TrkTrack* t_track = new TrkTrack(); -// general variables - good2 = l2->good2; + + // **** general variables **** +// good2 = l2->good2; for(Int_t i=0; i<12 ; i++){ - crc[i] = l2->crc[i]; - }; -// *** TRACKS *** +// crc[i] = l2->crc[i]; + good[i] = l2->good[i]; + }; + // *** TRACKS *** TClonesArray &t = *Track; for(int i=0; intrk; 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++) @@ -365,9 +585,10 @@ TrkSinglet* t_singlet = new TrkSinglet(); TrkTrack* t_track = new TrkTrack(); // general variables - good2 = l2->good2; +// good2 = l2->good2; for(Int_t i=0; i<12 ; i++){ - crc[i] = l2->crc[i]; +// crc[i] = l2->crc[i]; + good[i] = l2->good[i]; }; // *** TRACKS *** TClonesArray &t = *Track; @@ -376,6 +597,7 @@ 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++) @@ -402,6 +624,8 @@ // 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 "<cltry[i][ip] "<cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<GetCluster(l2->cltry[i][ip]-1)<xgood[ip])cout<<" i "<cltrx[i][ip] "<cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<GetCluster(l2->cltrx[i][ip]-1)<clsx[i] <cls = l1->GetCluster(l2->clsx[i]-1); +// cout<<" i "<clsx[i] "<clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<GetCluster(l2->clsx[i]-1)<Clear(); @@ -430,6 +655,7 @@ //----------------------------------------------------- // cout << "singolo y "<clsy[i] <cls = l1->GetCluster(l2->clsy[i]-1); +// cout<<" i "<clsy[i] "<clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<GetCluster(l2->clsy[i]-1)<Clear(); @@ -445,9 +671,10 @@ 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 *** @@ -455,7 +682,8 @@ 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->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]; @@ -499,10 +727,11 @@ // //-------------------------------------- void TrkLevel2::Clear(){ - good2 = -1; +// good2 = -1; for(Int_t i=0; i<12 ; i++){ - crc[i] = -1; - }; +// crc[i] = -1; + good[i] = -1; + }; /* Track->RemoveAll(); SingletX->RemoveAll(); SingletY->RemoveAll();*/ @@ -515,101 +744,23 @@ // // //-------------------------------------- +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. */ -// TClonesArray *TrkLevel2::GetTracks(){ -// TClonesArray *sorted = GetTracks_NFitSorted(); -// return sorted; -// }; - -/*TClonesArray *TrkLevel2::GetTracks_Chi2Sorted(){ - - TClonesArray *sorted = new TClonesArray("TrkTrack"); - TClonesArray &t = *Track; - TClonesArray &ts = *sorted; - int N=this->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; - } - } - if( ((TrkTrack *)t[indi])->image != -1 ){ - m[((TrkTrack *)t[indi])->image] = 0; - N--; - } - new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); - m[indi] = 0; - N--; - indo++; - } - return sorted; -}*/ -/*TClonesArray *TrkLevel2::GetTracks_NFitSorted(){ - - TClonesArray *sorted = new TClonesArray("TrkTrack"); - TClonesArray &t = *Track; - TClonesArray &ts = *sorted; -// TClonesArray &ts = *PhysicalTrack; - int N = ntrk(); - vector 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 <<" "<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; @@ -731,21 +877,7 @@ * 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()){ @@ -766,7 +898,7 @@ Float_t ntot=0; TClonesArray &t = *Track; - for(int i=0; iGetImageSeqNo() == -1 ) ntot+=1.; else ntot+=0.5; } @@ -781,31 +913,6 @@ * 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; - } - - 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; - 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()){ @@ -984,6 +1091,7 @@ } + ClassImp(TrkLevel2); ClassImp(TrkSinglet); ClassImp(TrkTrack);