12 |
extern "C" { |
extern "C" { |
13 |
void dotrack_(int*, double*, double*, double*, double*, int*); |
void dotrack_(int*, double*, double*, double*, double*, int*); |
14 |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
15 |
|
void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*); |
16 |
void mini2_(int*,int*,int*); |
void mini2_(int*,int*,int*); |
17 |
void guess_(); |
void guess_(); |
18 |
void gufld_(float*, float*); |
void gufld_(float*, float*); |
178 |
// |
// |
179 |
//-------------------------------------- |
//-------------------------------------- |
180 |
/** |
/** |
|
* 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. |
|
181 |
* |
* |
182 |
* >>> OBSOLETE !!! use TrkTrack::DoTrack2(Trajectory* t) instead |
* >>> OBSOLETE !!! use TrkTrack::DoTrack(Trajectory* t) instead |
183 |
* |
* |
184 |
*/ |
*/ |
185 |
int TrkTrack::DoTrack(Trajectory* t){ |
int TrkTrack::DoTrack2(Trajectory* t){ |
|
|
|
|
cout << " int TrkTrack::DoTrack(Trajectory* t) --->> OBSOLETE !!! "<<endl; |
|
|
cout << " use int TrkTrack::DoTrack2(Trajectory* t)"<<endl; |
|
|
|
|
|
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]; |
|
186 |
|
|
187 |
TrkParams::Load(1); |
cout << endl; |
188 |
if( !TrkParams::IsLoaded(1) ){ |
cout << " int TrkTrack::DoTrack2(Trajectory* t) --->> NB NB !! this method is going to be eliminated !!! "<<endl; |
189 |
cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; |
190 |
return 0; |
cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; |
191 |
} |
cout << endl; |
|
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++; |
|
|
} |
|
192 |
|
|
193 |
// delete [] dxout; |
return DoTrack(t); |
|
// delete [] dyout; |
|
|
// delete [] dzin; |
|
194 |
|
|
|
return ifail; |
|
195 |
}; |
}; |
196 |
//-------------------------------------- |
//-------------------------------------- |
197 |
// |
// |
198 |
// |
// |
199 |
//-------------------------------------- |
//-------------------------------------- |
200 |
/** |
/** |
201 |
* Evaluates the trajectory in the apparatus associated to the track. |
* Evaluates the trajectory in the apparatus associated to the track state-vector. |
202 |
* 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. |
* It integrates the equations of motion in the magnetic field. |
203 |
* @param t pointer to an object of the class Trajectory, |
* @param t pointer to an object of the class Trajectory, |
204 |
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
* which z coordinates should be previously assigned. |
205 |
* @return error flag. |
* @return error flag. |
206 |
*/ |
*/ |
207 |
int TrkTrack::DoTrack2(Trajectory* t){ |
int TrkTrack::DoTrack(Trajectory* t){ |
208 |
|
|
209 |
double *dxout = new double[t->npoint]; |
double *dxout = new double[t->npoint]; |
210 |
double *dyout = new double[t->npoint]; |
double *dyout = new double[t->npoint]; |
221 |
|
|
222 |
TrkParams::Load(1); |
TrkParams::Load(1); |
223 |
if( !TrkParams::IsLoaded(1) ){ |
if( !TrkParams::IsLoaded(1) ){ |
224 |
cout << "int TrkTrack::DoTrack2(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
225 |
return 0; |
return 0; |
226 |
} |
} |
227 |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
228 |
|
|
229 |
for (int i=0; i<t->npoint; i++){ |
for (int i=0; i<t->npoint; i++){ |
230 |
t->x[i] = (float)*dxout++; |
t->x[i] = (float)*(dxout+i); |
231 |
t->y[i] = (float)*dyout++; |
t->y[i] = (float)*(dyout+i); |
232 |
t->thx[i] = (float)*dthxout++; |
t->thx[i] = (float)*(dthxout+i); |
233 |
t->thy[i] = (float)*dthyout++; |
t->thy[i] = (float)*(dthyout+i); |
234 |
t->tl[i] = (float)*dtlout++; |
t->tl[i] = (float)*(dtlout+i); |
235 |
} |
} |
236 |
|
|
237 |
// delete [] dxout; |
delete [] dxout; |
238 |
// delete [] dyout; |
delete [] dyout; |
239 |
// delete [] dzin; |
delete [] dzin; |
240 |
|
delete [] dthxout; |
241 |
|
delete [] dthyout; |
242 |
|
delete [] dtlout; |
243 |
|
|
244 |
return ifail; |
return ifail; |
245 |
}; |
}; |
374 |
return (last_plane-first_plane+1); |
return (last_plane-first_plane+1); |
375 |
} |
} |
376 |
/** |
/** |
377 |
|
* Returns the track "lever-arm" on the x+y view, defined as the distance (in planes) between |
378 |
|
* the upper and lower x,y (couple) measurements (the maximum value of lever-arm is 6). |
379 |
|
*/ |
380 |
|
Int_t TrkTrack::GetLeverArmXY(){ |
381 |
|
int first_plane = -1; |
382 |
|
int last_plane = -1; |
383 |
|
for(Int_t ip=0; ip<6; ip++){ |
384 |
|
if( XGood(ip)*YGood(ip) && first_plane == -1 )first_plane = ip; |
385 |
|
if( XGood(ip)*YGood(ip) && first_plane != -1 )last_plane = ip; |
386 |
|
} |
387 |
|
if( first_plane == -1 || last_plane == -1){ |
388 |
|
cout<< "Int_t TrkTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl; |
389 |
|
return 0; |
390 |
|
} |
391 |
|
return (last_plane-first_plane+1); |
392 |
|
} |
393 |
|
/** |
394 |
* Returns the reduced chi-square of track x-projection |
* Returns the reduced chi-square of track x-projection |
395 |
*/ |
*/ |
396 |
Float_t TrkTrack::GetChi2X(){ |
Float_t TrkTrack::GetChi2X(){ |
601 |
4.52043, |
4.52043, |
602 |
4.29926}; |
4.29926}; |
603 |
int index; |
int index; |
604 |
float fact; |
float fact=0.; |
605 |
for(int i=0; i<6; i++) { |
for(int i=0; i<6; i++) { |
606 |
index = int((fabs(axv[i])+1.)/2.); |
index = int((fabs(axv[i])+1.)/2.); |
607 |
if(index>10) index=10; |
if(index>10) index=10; |
648 |
|
|
649 |
TrkParams::Set(path,1); |
TrkParams::Set(path,1); |
650 |
TrkParams::Load(1); |
TrkParams::Load(1); |
651 |
|
if( !TrkParams::IsLoaded(1) ){ |
652 |
|
cout << "void TrkTrack::LoadField(TString path) --- ERROR --- m.field not loaded"<<endl; |
653 |
|
} |
654 |
|
|
655 |
}; |
}; |
656 |
|
|
672 |
track.zm[i]=zm[i]; |
track.zm[i]=zm[i]; |
673 |
|
|
674 |
// --- temporaneo ---------------------------- |
// --- temporaneo ---------------------------- |
675 |
// andrebbe inserita la dimensione del sensore |
// float segment = 100.; |
676 |
float segment = 100.; |
// track.xm_a[i]=xm[i]; |
677 |
track.xm_a[i]=xm[i]; |
// track.xm_b[i]=xm[i]; |
678 |
track.xm_b[i]=xm[i]; |
// track.ym_a[i]=ym[i]; |
679 |
track.ym_a[i]=ym[i]; |
// track.ym_b[i]=ym[i]; |
680 |
track.ym_b[i]=ym[i]; |
// if( XGood(i) && !YGood(i) ){ |
681 |
if( XGood(i) && !YGood(i) ){ |
// track.ym_a[i] = track.ym_a[i]+segment; |
682 |
track.ym_a[i] = track.ym_a[i]+segment; |
// track.ym_b[i] = track.ym_b[i]-segment; |
683 |
track.ym_b[i] = track.ym_b[i]-segment; |
// }else if( !XGood(i) && YGood(i)){ |
684 |
}else if( !XGood(i) && YGood(i)){ |
// track.xm_a[i] = track.xm_a[i]+segment; |
685 |
track.xm_a[i] = track.xm_a[i]+segment; |
// track.xm_b[i] = track.xm_b[i]-segment; |
686 |
track.xm_b[i] = track.xm_b[i]-segment; |
// } |
|
} |
|
687 |
// --- temporaneo ---------------------------- |
// --- temporaneo ---------------------------- |
688 |
|
|
689 |
|
if( XGood(i) || YGood(i) ){ |
690 |
|
double segment = 2.;//cm |
691 |
|
// NB: i parametri di allineamento hanno una notazione particolare!!! |
692 |
|
// sensor = 0 (hybrid side), 1 |
693 |
|
// ladder = 0-2 (increasing x) |
694 |
|
// plane = 0-5 (from bottom to top!!!) |
695 |
|
int is = (int)GetSensor(i); if(i==5)is=1-is; |
696 |
|
int ip = 5-i; |
697 |
|
int il = (int)GetLadder(i); |
698 |
|
|
699 |
|
double omega = 0.; |
700 |
|
double beta = 0.; |
701 |
|
double gamma = 0.; |
702 |
|
if( |
703 |
|
(is < 0 || is > 1 || ip < 0 || ip > 5 || il < 0 || il > 2) && |
704 |
|
true){ |
705 |
|
// se il piano risulta colpito, ladder e sensore devono essere |
706 |
|
// assegnati correttamente |
707 |
|
cout << " void TrkTrack::FillMiniStruct(cMini2track&) --- WARNING --- sensor not defined, cannot read alignment parameters "<<endl; |
708 |
|
cout << " is ip il = "<<is<<" "<<ip<<" "<<il<<endl; |
709 |
|
}else{ |
710 |
|
omega = alignparameters_.omega[is][il][ip]; |
711 |
|
beta = alignparameters_.beta[is][il][ip]; |
712 |
|
gamma = alignparameters_.gamma[is][il][ip]; |
713 |
|
} |
714 |
|
|
715 |
|
if( XGood(i) && !YGood(i) ){ |
716 |
|
track.xm_a[i] = xm[i] - omega * segment; |
717 |
|
track.ym_a[i] = ym[i] + segment; |
718 |
|
// track.zm_a[i] = zm[i] + beta * segment;//not used yet |
719 |
|
track.xm_b[i] = xm[i] + omega * segment; |
720 |
|
track.ym_b[i] = ym[i] - segment; |
721 |
|
// track.zm_b[i] = zm[i] - beta * segment;//not used yet |
722 |
|
}else if( !XGood(i) && YGood(i) ){ |
723 |
|
track.xm_a[i] = xm[i] + segment; |
724 |
|
track.ym_a[i] = ym[i] + omega * segment; |
725 |
|
// track.zm_a[i] = zm[i] - gamma * segment;//not used yet |
726 |
|
track.xm_b[i] = xm[i] - segment; |
727 |
|
track.ym_b[i] = ym[i] - omega * segment; |
728 |
|
// track.zm_b[i] = zm[i] + gamma * segment;//not used yet |
729 |
|
} |
730 |
|
} |
731 |
|
|
732 |
track.resx[i]=resx[i]; |
track.resx[i]=resx[i]; |
733 |
track.resy[i]=resy[i]; |
track.resy[i]=resy[i]; |
734 |
track.tailx[i]=tailx[i]; |
track.tailx[i]=tailx[i]; |
735 |
track.taily[i]=taily[i]; |
track.taily[i]=taily[i]; |
792 |
|
|
793 |
// cout << "void TrkTrack::GetClusterositions() "<<endl; |
// cout << "void TrkTrack::GetClusterositions() "<<endl; |
794 |
|
|
795 |
TrkParams::Load( ); |
TrkParams::Load(1); |
796 |
if( !TrkParams::IsLoaded() )return false; |
if( !TrkParams::IsLoaded(1) ){ |
797 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- m.field not loaded "<<endl; |
798 |
|
return false; |
799 |
|
} |
800 |
|
TrkParams::Load(4); |
801 |
|
if( !TrkParams::IsLoaded(4) ){ |
802 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- p.f.a. par. not loaded "<<endl; |
803 |
|
return false; |
804 |
|
} |
805 |
|
TrkParams::Load(5); |
806 |
|
if( !TrkParams::IsLoaded(5) ){ |
807 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- alignment par. not loaded "<<endl; |
808 |
|
return false; |
809 |
|
} |
810 |
|
|
811 |
for(int ip=0; ip<6; ip++){ |
for(int ip=0; ip<6; ip++){ |
812 |
// cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;; |
// cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;; |
813 |
int icx = GetClusterX_ID(ip)+1; |
int icx = GetClusterX_ID(ip)+1; |
859 |
*/ |
*/ |
860 |
void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){ |
void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){ |
861 |
|
|
862 |
float al_ini[] = {0.,0.,0.,0.,0.}; |
TrkParams::Load(1); |
863 |
|
if( !TrkParams::IsLoaded(1) ){ |
864 |
|
cout << "void TrkTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl; |
865 |
|
return; |
866 |
|
} |
867 |
|
TrkParams::Load(5); |
868 |
|
if( !TrkParams::IsLoaded(5) ){ |
869 |
|
cout << "void TrkTrack::Fit(double,int&,int,int) ---ERROR--- align.param. not loaded "<<endl; |
870 |
|
return; |
871 |
|
} |
872 |
|
|
873 |
TrkParams::Load( ); |
float al_ini[] = {0.,0.,0.,0.,0.}; |
|
if( !TrkParams::IsLoaded() )return; |
|
874 |
|
|
875 |
extern cMini2track track_; |
extern cMini2track track_; |
876 |
fail = 0; |
fail = 0; |
906 |
|
|
907 |
// ------------------------------------------ |
// ------------------------------------------ |
908 |
// call mini routine |
// call mini routine |
909 |
// TrkParams::Load(1); |
// ------------------------------------------ |
|
// if( !TrkParams::IsLoaded(1) ){ |
|
|
// cout << "void TrkTrack::Fit(double pfixed, int& fail, int iprint) --- ERROR --- m.field not loaded"<<endl; |
|
|
// return; |
|
|
// } |
|
910 |
int istep=0; |
int istep=0; |
911 |
int ifail=0; |
int ifail=0; |
912 |
mini2_(&istep,&ifail, &iprint); |
mini2_(&istep,&ifail, &iprint); |
1052 |
}; |
}; |
1053 |
|
|
1054 |
/** |
/** |
1055 |
* Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane. |
* Method to retrieve the ladder (0-2, increasing x) traversed by the track on this plane. |
1056 |
* If no ladder is traversed (dead area) the metod retuns -1. |
* If no ladder is traversed (dead area) the metod retuns -1. |
1057 |
* @param ip Tracker plane (0-5) |
* @param ip Tracker plane (0-5) |
1058 |
*/ |
*/ |
1075 |
/** |
/** |
1076 |
* \brief Method to include a x-cluster to the track. |
* \brief Method to include a x-cluster to the track. |
1077 |
* @param ip Tracker plane (0-5) |
* @param ip Tracker plane (0-5) |
1078 |
* @param clid Cluster ID (0,1,...) |
* @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) |
1079 |
* @param is Sensor (0-1, increasing y) |
* @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) |
1080 |
|
* @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) |
1081 |
|
* @param bad True if the cluster contains bad strips |
1082 |
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
1083 |
*/ |
*/ |
1084 |
void TrkTrack::SetXGood(int ip, int clid, int is){ |
void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){ |
1085 |
int il=0; //ladder (temporary) |
// int il=0; //ladder (temporary) |
1086 |
bool bad=false; //ladder (temporary) |
// bool bad=false; //ladder (temporary) |
1087 |
xgood[ip]=il*100000000+is*10000000+clid; |
if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1) |
1088 |
|
cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; |
1089 |
|
xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid; |
1090 |
if(bad)xgood[ip]=-xgood[ip]; |
if(bad)xgood[ip]=-xgood[ip]; |
1091 |
}; |
}; |
1092 |
/** |
/** |
1093 |
* \brief Method to include a y-cluster to the track. |
* \brief Method to include a y-cluster to the track. |
1094 |
* @param ip Tracker plane (0-5) |
* @param ip Tracker plane (0-5) |
1095 |
* @param clid Cluster ID (0,1,...) |
* @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) |
1096 |
* @param is Sensor (0-1) |
* @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) |
1097 |
|
* @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) |
1098 |
|
* @param bad True if the cluster contains bad strips |
1099 |
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
1100 |
*/ |
*/ |
1101 |
void TrkTrack::SetYGood(int ip, int clid, int is){ |
void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){ |
1102 |
int il=0; //ladder (temporary) |
// int il=0; //ladder (temporary) |
1103 |
bool bad=false; //ladder (temporary) |
// bool bad=false; //ladder (temporary) |
1104 |
ygood[ip]=il*100000000+is*10000000+clid; |
if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1) |
1105 |
|
cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; |
1106 |
|
ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid; |
1107 |
if(bad)ygood[ip]=-ygood[ip]; |
if(bad)ygood[ip]=-ygood[ip]; |
1108 |
}; |
}; |
1109 |
|
|
1206 |
vto = iv+1; |
vto = iv+1; |
1207 |
} |
} |
1208 |
for(int i=pfrom; i<pto; i++) |
for(int i=pfrom; i<pto; i++) |
1209 |
for(int j=0; j<vto; j++) |
for(int j=vfrom; j<vto; j++){ |
1210 |
if(GetDEDX(i,j)>max)max=GetDEDX(i,j); |
if(j==0 && XGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); |
1211 |
|
if(j==1 && YGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); |
1212 |
|
} |
1213 |
return max; |
return max; |
1214 |
|
|
1215 |
}; |
}; |
1232 |
vto = iv+1; |
vto = iv+1; |
1233 |
} |
} |
1234 |
for(int i=pfrom; i<pto; i++) |
for(int i=pfrom; i<pto; i++) |
1235 |
for(int j=0; j<vto; j++) |
for(int j=vfrom; j<vto; j++){ |
1236 |
if(GetDEDX(i,j)<min)min=GetDEDX(i,j); |
if(j==0 && XGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); |
1237 |
|
if(j==1 && YGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); |
1238 |
|
} |
1239 |
return min; |
return min; |
1240 |
|
|
1241 |
}; |
}; |
1242 |
|
|
1243 |
/** |
/** |
1244 |
* \brief Give the maximum spatial residual release |
* \brief Give the maximum spatial residual |
1245 |
*/ |
*/ |
1246 |
Float_t TrkTrack::GetResidual_max(int ip, int iv){ |
Float_t TrkTrack::GetResidual_max(int ip, int iv){ |
1247 |
Float_t max=0; |
Float_t max=0; |
1258 |
vto = iv+1; |
vto = iv+1; |
1259 |
} |
} |
1260 |
for(int i=pfrom; i<pto; i++){ |
for(int i=pfrom; i<pto; i++){ |
1261 |
for(int j=0; j<vto; j++){ |
for(int j=vfrom; j<vto; j++){ |
1262 |
if(j==0 && XGood(i) && fabs(xm[i]-xv[i])>fabs(max))max=xv[i]-xm[i]; |
if(j==0 && XGood(i) && fabs(xm[i]-xv[i])>fabs(max))max=xm[i]-xv[i]; |
1263 |
if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=yv[i]-ym[i]; |
if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=ym[i]-yv[i]; |
1264 |
} |
} |
1265 |
} |
} |
1266 |
return max; |
return max; |
1267 |
|
|
1268 |
}; |
}; |
1269 |
|
/** |
1270 |
|
* \brief Give the anerage spatial residual |
1271 |
|
*/ |
1272 |
|
Float_t TrkTrack::GetResidual_av(int ip, int iv){ |
1273 |
|
// |
1274 |
|
//Sum$((xm>-50)*(xm-xv)/resx)/sqrt(TrkTrack.GetNX()*TrkTrack.GetChi2X())<0.3 |
1275 |
|
|
1276 |
|
Float_t av = 0.; |
1277 |
|
int nav = 0; |
1278 |
|
// |
1279 |
|
int pfrom = 0; |
1280 |
|
int pto = 6; |
1281 |
|
int vfrom = 0; |
1282 |
|
int vto = 2; |
1283 |
|
if(ip>=0&&ip<6){ |
1284 |
|
pfrom = ip; |
1285 |
|
pto = ip+1; |
1286 |
|
} |
1287 |
|
if(iv>=0&&iv<2){ |
1288 |
|
vfrom = iv; |
1289 |
|
vto = iv+1; |
1290 |
|
} |
1291 |
|
for(int i=pfrom; i<pto; i++){ |
1292 |
|
for(int j=vfrom; j<vto; j++){ |
1293 |
|
nav++; |
1294 |
|
if(j==0 && XGood(i)) av += (xm[i]-xv[i])/resx[i]; |
1295 |
|
if(j==1 && YGood(i)) av += (ym[i]-yv[i])/resy[i]; |
1296 |
|
} |
1297 |
|
} |
1298 |
|
if(nav==0)return -100.; |
1299 |
|
return av/nav; |
1300 |
|
|
1301 |
|
}; |
1302 |
|
|
1303 |
|
|
1304 |
/** |
/** |
1681 |
int ngf = TrkParams::nGF; |
int ngf = TrkParams::nGF; |
1682 |
float *zgf = TrkParams::zGF; |
float *zgf = TrkParams::zGF; |
1683 |
Trajectory tgf = Trajectory(ngf,zgf); |
Trajectory tgf = Trajectory(ngf,zgf); |
1684 |
tgf.DoTrack2(t_track->al);//<<<< integrate the trajectory |
tgf.DoTrack(t_track->al);//<<<< integrate the trajectory |
1685 |
for(int ip=0; ip<ngf; ip++){ |
for(int ip=0; ip<ngf; ip++){ |
1686 |
t_track->xGF[ip] = tgf.x[ip]; |
t_track->xGF[ip] = tgf.x[ip]; |
1687 |
t_track->yGF[ip] = tgf.y[ip]; |
t_track->yGF[ip] = tgf.y[ip]; |
1839 |
|
|
1840 |
if(!Track)return 0; |
if(!Track)return 0; |
1841 |
|
|
1842 |
TRefArray *sorted = new TRefArray(); |
// TRefArray *sorted = new TRefArray(); |
1843 |
|
TRefArray *sorted = NULL; |
1844 |
|
|
1845 |
TClonesArray &t = *Track; |
TClonesArray &t = *Track; |
1846 |
// TClonesArray &ts = *PhysicalTrack; |
// TClonesArray &ts = *PhysicalTrack; |
1878 |
|
|
1879 |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
1880 |
}; |
}; |
1881 |
|
if(!sorted)sorted = new TRefArray( TProcessID::GetProcessWithUID(t[indi])); |
1882 |
sorted->Add( (TrkTrack*)t[indi] ); |
sorted->Add( (TrkTrack*)t[indi] ); |
1883 |
|
|
1884 |
m[indi] = 0; |
m[indi] = 0; |
2049 |
|
|
2050 |
TrkParams::Set(path,1); |
TrkParams::Set(path,1); |
2051 |
TrkParams::Load(1); |
TrkParams::Load(1); |
2052 |
|
if( !TrkParams::IsLoaded(1) ){ |
2053 |
|
cout << "void TrkLevel2::LoadField(TString path) --- ERROR --- m.field not loaded"<<endl; |
2054 |
|
} |
2055 |
|
|
2056 |
// |
// |
2057 |
}; |
}; |
2239 |
|
|
2240 |
/** |
/** |
2241 |
* Evaluates the trajectory in the apparatus associated to the track. |
* Evaluates the trajectory in the apparatus associated to the track. |
2242 |
* 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. |
* It integrates the equations of motion in the magnetic field. |
2243 |
* @param t pointer to an object of the class Trajectory, |
* @param al Track state-vector (X0,Y0,sin(theta),phi,deflection). |
2244 |
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
* @param zini z-coordinate of the reference plane (Z0). |
2245 |
* @return error flag. |
* @return error flag. |
2246 |
*/ |
* |
2247 |
int Trajectory::DoTrack2(float* al){ |
* This method is needed when you want to integrate the particle trajectory |
2248 |
|
* starting from a track state-vector relative to an arbitrary reference plane. |
2249 |
double *dxout = new double[npoint]; |
* The default reference plane, used by the tracker routines, is at zini=23.5. |
2250 |
double *dyout = new double[npoint]; |
* If you give as input the track state-vector from a TrkTrack object, |
2251 |
double *dthxout = new double[npoint]; |
* you can use Trajectory::DoTrack(float* al) instead. |
2252 |
double *dthyout = new double[npoint]; |
*/ |
2253 |
double *dtlout = new double[npoint]; |
int Trajectory::DoTrack(float* al, float zini){ |
2254 |
double *dzin = new double[npoint]; |
|
2255 |
|
// double *dxout = new double[npoint]; |
2256 |
|
// double *dyout = new double[npoint]; |
2257 |
|
// double *dthxout = new double[npoint]; |
2258 |
|
// double *dthyout = new double[npoint]; |
2259 |
|
// double *dtlout = new double[npoint]; |
2260 |
|
// double *dzin = new double[npoint]; |
2261 |
|
|
2262 |
|
double *dxout; |
2263 |
|
double *dyout; |
2264 |
|
double *dthxout; |
2265 |
|
double *dthyout; |
2266 |
|
double *dtlout; |
2267 |
|
double *dzin; |
2268 |
|
|
2269 |
|
dxout = (double*) malloc(npoint*sizeof(double)); |
2270 |
|
dyout = (double*) malloc(npoint*sizeof(double)); |
2271 |
|
dthxout = (double*) malloc(npoint*sizeof(double)); |
2272 |
|
dthyout = (double*) malloc(npoint*sizeof(double)); |
2273 |
|
dtlout = (double*) malloc(npoint*sizeof(double)); |
2274 |
|
dzin = (double*) malloc(npoint*sizeof(double)); |
2275 |
|
|
2276 |
double dal[5]; |
double dal[5]; |
2277 |
|
|
2278 |
int ifail = 0; |
double dzini = (double)zini; |
2279 |
|
|
2280 |
|
int ifail = 0; |
2281 |
|
|
2282 |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
2283 |
for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; |
for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; |
2284 |
|
|
2285 |
TrkParams::Load(1); |
TrkParams::Load(1); |
2286 |
if( !TrkParams::IsLoaded(1) ){ |
if( !TrkParams::IsLoaded(1) ){ |
2287 |
cout << "int Trajectory::DoTrack2(float* al) --- ERROR --- m.field not loaded"<<endl; |
cout << "int Trajectory::DoTrack(float* al) --- ERROR --- m.field not loaded"<<endl; |
2288 |
return 0; |
return 0; |
2289 |
} |
} |
2290 |
dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
// dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
2291 |
|
dotrack3_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&dzini,&ifail); |
2292 |
|
|
2293 |
for (int i=0; i<npoint; i++){ |
for (int i=0; i<npoint; i++){ |
2294 |
x[i] = (float)*dxout++; |
x[i] = (float)*(dxout+i); |
2295 |
y[i] = (float)*dyout++; |
y[i] = (float)*(dyout+i); |
2296 |
thx[i] = (float)*dthxout++; |
thx[i] = (float)*(dthxout+i); |
2297 |
thy[i] = (float)*dthyout++; |
thy[i] = (float)*(dthyout+i); |
2298 |
tl[i] = (float)*dtlout++; |
tl[i] = (float)*(dtlout+i); |
2299 |
} |
} |
2300 |
|
|
2301 |
|
if(dxout) free( dxout ); |
2302 |
|
if(dyout) free( dyout ); |
2303 |
|
if(dthxout)free( dthxout ); |
2304 |
|
if(dthyout)free( dthyout ); |
2305 |
|
if(dtlout) free( dtlout ); |
2306 |
|
if(dzin) free( dzin ); |
2307 |
|
|
2308 |
|
// delete [] dxout; |
2309 |
|
// delete [] dyout; |
2310 |
|
// delete [] dthxout; |
2311 |
|
// delete [] dthyout; |
2312 |
|
// delete [] dtlout; |
2313 |
|
// delete [] dzin; |
2314 |
|
|
2315 |
|
|
2316 |
return ifail; |
return ifail; |
2317 |
}; |
}; |
2318 |
|
|
2319 |
|
/** |
2320 |
|
* |
2321 |
|
* >>> OBSOLETE !!! use Trajectory::DoTrack(float* al, float zini) instead |
2322 |
|
* |
2323 |
|
*/ |
2324 |
|
int Trajectory::DoTrack2(float* al, float zini){ |
2325 |
|
|
2326 |
|
cout << endl; |
2327 |
|
cout << " int Trajectory::DoTrack2(float* al, float zini) --->> NB NB !! this method is going to be eliminated !!! "<<endl; |
2328 |
|
cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; |
2329 |
|
cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; |
2330 |
|
cout << endl; |
2331 |
|
|
2332 |
|
return DoTrack(al,zini); |
2333 |
|
|
2334 |
|
}; |
2335 |
|
|
2336 |
|
|
2337 |
|
|
2338 |
ClassImp(TrkLevel2); |
ClassImp(TrkLevel2); |
2339 |
ClassImp(TrkSinglet); |
ClassImp(TrkSinglet); |
2340 |
ClassImp(TrkTrack); |
ClassImp(TrkTrack); |