/** * \file TrkLevel2.cpp * \author Elena Vannuccini */ #include #include #include 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*); 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 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; tailx[ip] = 0; taily[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; 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 *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; inpoint; i++) dzin[i] = (double)t->z[i]; 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+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 [] dthxout; delete [] dthyout; delete [] dtlout; 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; }; // /** * Method to retrieve the dE/dx measured on a tracker view. * @param ip plane (0-5) * @param iv view (0=x 1=y) */ 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 "<=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; imax)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(); }; //-------------------------------------- // // //-------------------------------------- void TrkTrack::Delete(){ // cout << "TrkTrack::Delete()"<GetClass()->GetName(),"TrkTrack")==0){ if(Track)Track->Clear("C"); Track = track; } } //-------------------------------------- // // //-------------------------------------- void TrkLevel2::Dump(){ // cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; 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). * 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, TrkLevel1 *l1){ // cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<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;// NBNBNBNB deve sempre essere = i t_track->image = l2->image[i]-1; 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++){ // --------------------------------- // 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]; 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->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; } /** * 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 *** 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 *** 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; } } } //-------------------------------------- // // //-------------------------------------- void TrkLevel2::Clear(){ for(Int_t i=0; i<12 ; i++){ good[i] = -1; VKflag[i] = 0; VKmask[i] = 0; }; // 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()"< m(N); for(int i=0; i 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(); } } //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--; // cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<Add( (TrkTrack*)t[indi] ); m[indi] = 0; // cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "<= this->ntrk()){ 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 << "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; } //-------------------------------------- // // //-------------------------------------- /** * 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 << "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; } //-------------------------------------- // // //-------------------------------------- /** * 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 << "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->Clear(); delete sorted; return track; } /** * Give the number of "physical" tracks, sorted by the method GetTracks(). */ 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; } 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 << "TrkTrack *TrkLevel2::GetTrackImage(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); 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]; sorted->Delete(); delete sorted; 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_(); // 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"<>> 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; i0)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 "<> " << 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. * @param al Track state-vector (X0,Y0,sin(theta),phi,deflection). * @param zini z-coordinate of the reference plane (Z0). * @return error flag. * * This method is needed when you want to integrate the particle trajectory * starting from a track state-vector relative to an arbitrary reference plane. * The default reference plane, used by the tracker routines, is at zini=23.5. * If you give as input the track state-vector from a TrkTrack object, * you can use Trajectory::DoTrack(float* al) instead. */ int Trajectory::DoTrack(float* al, float zini){ // 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 *dxout; double *dyout; double *dthxout; double *dthyout; double *dtlout; double *dzin; dxout = (double*) malloc(npoint*sizeof(double)); dyout = (double*) malloc(npoint*sizeof(double)); dthxout = (double*) malloc(npoint*sizeof(double)); dthyout = (double*) malloc(npoint*sizeof(double)); dtlout = (double*) malloc(npoint*sizeof(double)); dzin = (double*) malloc(npoint*sizeof(double)); double dal[5]; double dzini = (double)zini; int ifail = 0; for (int i=0; i<5; i++) dal[i] = (double)al[i]; for (int i=0; i>> 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) <<<<"<