| 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 | //    int  readb_(const char*); | void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*); | 
|  | int  readb_(); |  | 
| 16 | void mini2_(int*,int*,int*); | void mini2_(int*,int*,int*); | 
| 17 | void guess_(); | void guess_(); | 
| 18 |  | void gufld_(float*, float*); | 
| 19 |  | float risxeta2_(float *); | 
| 20 |  | float risxeta3_(float *); | 
| 21 |  | float risxeta4_(float *); | 
| 22 |  | float risyeta2_(float *); | 
| 23 | } | } | 
| 24 |  |  | 
| 25 | //-------------------------------------- | //-------------------------------------- | 
| 26 | // | // | 
| 27 | // | // | 
| 44 | zm[ip]     = 0; | zm[ip]     = 0; | 
| 45 | resx[ip]   = 0; | resx[ip]   = 0; | 
| 46 | resy[ip]   = 0; | resy[ip]   = 0; | 
| 47 |  | tailx[ip]   = 0; | 
| 48 |  | taily[ip]   = 0; | 
| 49 | xv[ip]     = 0; | xv[ip]     = 0; | 
| 50 | yv[ip]     = 0; | yv[ip]     = 0; | 
| 51 | zv[ip]     = 0; | zv[ip]     = 0; | 
| 53 | ayv[ip]    = 0; | ayv[ip]    = 0; | 
| 54 | dedx_x[ip] = 0; | dedx_x[ip] = 0; | 
| 55 | dedx_y[ip] = 0; | dedx_y[ip] = 0; | 
| 56 |  | multmaxx[ip] = 0; | 
| 57 |  | multmaxy[ip] = 0; | 
| 58 |  | seedx[ip]  = 0; | 
| 59 |  | seedy[ip]  = 0; | 
| 60 |  | xpu[ip]    = 0; | 
| 61 |  | ypu[ip]    = 0; | 
| 62 |  |  | 
| 63 | }; | }; | 
| 64 | clx = 0; |  | 
| 65 | cly = 0; | //     TrkParams::SetTrackingMode(); | 
| 66 |  | //     TrkParams::SetPrecisionFactor(); | 
| 67 |  | //     TrkParams::SetStepMin(); | 
| 68 |  | TrkParams::SetMiniDefault(); | 
| 69 |  | TrkParams::SetPFA(); | 
| 70 |  |  | 
| 71 |  | int ngf = TrkParams::nGF; | 
| 72 |  | for(int i=0; i<ngf; i++){ | 
| 73 |  | xGF[i] = 0.; | 
| 74 |  | yGF[i] = 0.; | 
| 75 |  | } | 
| 76 |  |  | 
| 77 |  |  | 
| 78 | }; | }; | 
| 79 | //-------------------------------------- | //-------------------------------------- | 
| 80 | // | // | 
| 97 | zm[ip]     = t.zm[ip]; | zm[ip]     = t.zm[ip]; | 
| 98 | resx[ip]   = t.resx[ip]; | resx[ip]   = t.resx[ip]; | 
| 99 | resy[ip]   = t.resy[ip]; | resy[ip]   = t.resy[ip]; | 
| 100 |  | tailx[ip]  = t.tailx[ip]; | 
| 101 |  | taily[ip]  = t.taily[ip]; | 
| 102 | xv[ip]     = t.xv[ip]; | xv[ip]     = t.xv[ip]; | 
| 103 | yv[ip]     = t.yv[ip]; | yv[ip]     = t.yv[ip]; | 
| 104 | zv[ip]     = t.zv[ip]; | zv[ip]     = t.zv[ip]; | 
| 106 | ayv[ip]    = t.ayv[ip]; | ayv[ip]    = t.ayv[ip]; | 
| 107 | dedx_x[ip] = t.dedx_x[ip]; | dedx_x[ip] = t.dedx_x[ip]; | 
| 108 | dedx_y[ip] = t.dedx_y[ip]; | dedx_y[ip] = t.dedx_y[ip]; | 
| 109 |  | multmaxx[ip] = t.multmaxx[ip]; | 
| 110 |  | multmaxy[ip] = t.multmaxy[ip]; | 
| 111 |  | seedx[ip]    = t.seedx[ip]; | 
| 112 |  | seedy[ip]    = t.seedy[ip]; | 
| 113 |  | xpu[ip]      = t.xpu[ip]; | 
| 114 |  | ypu[ip]      = t.ypu[ip]; | 
| 115 | }; | }; | 
| 116 | clx = new TRefArray(*(t.clx)); |  | 
| 117 | cly = new TRefArray(*(t.cly)); | //     TrkParams::SetTrackingMode(); | 
| 118 |  | //     TrkParams::SetPrecisionFactor(); | 
| 119 |  | //     TrkParams::SetStepMin(); | 
| 120 |  | TrkParams::SetMiniDefault(); | 
| 121 |  | TrkParams::SetPFA(); | 
| 122 |  |  | 
| 123 |  | int ngf = TrkParams::nGF; | 
| 124 |  | for(int i=0; i<ngf; i++){ | 
| 125 |  | xGF[i] = t.xGF[i]; | 
| 126 |  | yGF[i] = t.yGF[i]; | 
| 127 |  | } | 
| 128 | }; | }; | 
| 129 | //-------------------------------------- | //-------------------------------------- | 
| 130 | // | // | 
| 148 | t.zm[ip]     = zm[ip]; | t.zm[ip]     = zm[ip]; | 
| 149 | t.resx[ip]   = resx[ip]; | t.resx[ip]   = resx[ip]; | 
| 150 | t.resy[ip]   = resy[ip]; | t.resy[ip]   = resy[ip]; | 
| 151 |  | t.tailx[ip]  = tailx[ip]; | 
| 152 |  | t.taily[ip]  = taily[ip]; | 
| 153 | t.xv[ip]     = xv[ip]; | t.xv[ip]     = xv[ip]; | 
| 154 | t.yv[ip]     = yv[ip]; | t.yv[ip]     = yv[ip]; | 
| 155 | t.zv[ip]     = zv[ip]; | t.zv[ip]     = zv[ip]; | 
| 157 | t.ayv[ip]    = ayv[ip]; | t.ayv[ip]    = ayv[ip]; | 
| 158 | t.dedx_x[ip] = dedx_x[ip]; | t.dedx_x[ip] = dedx_x[ip]; | 
| 159 | t.dedx_y[ip] = dedx_y[ip]; | t.dedx_y[ip] = dedx_y[ip]; | 
| 160 |  | t.multmaxx[ip] = multmaxx[ip]; | 
| 161 |  | t.multmaxy[ip] = multmaxy[ip]; | 
| 162 |  | t.seedx[ip]    = seedx[ip]; | 
| 163 |  | t.seedy[ip]    = seedy[ip]; | 
| 164 |  | t.xpu[ip]      = xpu[ip]; | 
| 165 |  | t.ypu[ip]      = ypu[ip]; | 
| 166 |  |  | 
| 167 | }; | }; | 
| 168 |  | int ngf = TrkParams::nGF; | 
| 169 |  | for(int i=0; i<ngf; i++){ | 
| 170 |  | t.xGF[i] = xGF[i]; | 
| 171 |  | t.yGF[i] = yGF[i]; | 
| 172 |  | } | 
| 173 |  |  | 
| 174 |  |  | 
| 175 | }; | }; | 
| 176 | //-------------------------------------- | //-------------------------------------- | 
| 178 | // | // | 
| 179 | //-------------------------------------- | //-------------------------------------- | 
| 180 | /** | /** | 
| 181 | * Evaluates the trajectory in the apparatus associated to the track. | * | 
| 182 | * 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. | * >>> OBSOLETE !!! use TrkTrack::DoTrack(Trajectory* t) instead | 
| 183 | * @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. |  | 
| 184 | */ | */ | 
| 185 | int TrkTrack::DoTrack(Trajectory* t){ | int TrkTrack::DoTrack2(Trajectory* t){ | 
|  |  |  | 
|  | double *dxout = new double[t->npoint]; |  | 
|  | double *dyout = new double[t->npoint]; |  | 
|  | double *dzin = new double[t->npoint]; |  | 
|  | double dal[5]; |  | 
|  |  |  | 
|  | int ifail = 0; |  | 
|  |  |  | 
|  | for (int i=0; i<5; i++)         dal[i]  = (double)al[i]; |  | 
|  | for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |  | 
| 186 |  |  | 
| 187 | dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail); | cout << endl; | 
| 188 |  | cout << " int TrkTrack::DoTrack2(Trajectory* t) --->> NB NB !! this method is going to be eliminated !!! "<<endl; | 
| 189 | for (int i=0; i<t->npoint; i++){ | cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; | 
| 190 | t->x[i] = (float)*dxout++; | cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; | 
| 191 | t->y[i] = (float)*dyout++; | cout << endl; | 
|  | } |  | 
| 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]; | 
| 219 | for (int i=0; i<5; i++)         dal[i]  = (double)al[i]; | for (int i=0; i<5; i++)         dal[i]  = (double)al[i]; | 
| 220 | for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; | for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; | 
| 221 |  |  | 
| 222 |  | TrkParams::Load(1); | 
| 223 |  | if( !TrkParams::IsLoaded(1) ){ | 
| 224 |  | cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; | 
| 225 |  | 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 | }; | }; | 
| 266 | return def; | return def; | 
| 267 | }; | }; | 
| 268 | // | // | 
| 269 |  | /** | 
| 270 |  | * Method to retrieve the dE/dx measured on a tracker view. | 
| 271 |  | * @param ip plane (0-5) | 
| 272 |  | * @param iv view (0=x 1=y) | 
| 273 |  | */ | 
| 274 |  | Float_t TrkTrack::GetDEDX(int ip, int iv){ | 
| 275 |  | if(iv==0 && ip>=0 && ip<6)return fabs(dedx_x[ip]); | 
| 276 |  | else if(iv==1 && ip>=0 && ip<6)return fabs(dedx_y[ip]); | 
| 277 |  | else { | 
| 278 |  | cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; | 
| 279 |  | return 0.; | 
| 280 |  | } | 
| 281 |  | } | 
| 282 |  | /** | 
| 283 |  | * Method to evaluate the dE/dx measured on a tracker plane. | 
| 284 |  | * The two measurements on x- and y-view are averaged. | 
| 285 |  | * @param ip plane (0-5) | 
| 286 |  | */ | 
| 287 |  | Float_t TrkTrack::GetDEDX(int ip){ | 
| 288 |  | if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0; | 
| 289 |  | return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip)); | 
| 290 |  | }; | 
| 291 |  |  | 
| 292 |  | /** | 
| 293 |  | * Method to evaluate the dE/dx averaged over all planes. | 
| 294 |  | */ | 
| 295 | Float_t TrkTrack::GetDEDX(){ | Float_t TrkTrack::GetDEDX(){ | 
| 296 | Float_t dedx=0; | Float_t dedx=0; | 
| 297 | for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i]; | for(Int_t ip=0; ip<6; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip); | 
| 298 | dedx = dedx/(this->GetNX()+this->GetNY()); | dedx = dedx/(GetNX()+GetNY()); | 
| 299 | return dedx; | return dedx; | 
| 300 |  | }; | 
| 301 |  | /** | 
| 302 |  | * Returns 1 if the cluster on a tracker view includes bad strips | 
| 303 |  | * (at least one bad strip among the four strip used by p.f.a.) | 
| 304 |  | * @param ip plane (0-5) | 
| 305 |  | * @param iv view (0=x 1=y) | 
| 306 |  | */ | 
| 307 |  | Bool_t TrkTrack::IsBad(int ip,int iv){ | 
| 308 |  | if(iv==0 && ip>=0 && ip<6)return (xgood[ip]<0) ; | 
| 309 |  | else if(iv==1 && ip>=0 && ip<6)return (ygood[ip]<0) ; | 
| 310 |  | else { | 
| 311 |  | cout << "TrkTrack::IsBad(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; | 
| 312 |  | return 0.; | 
| 313 |  | } | 
| 314 |  | }; | 
| 315 |  | /** | 
| 316 |  | * Returns 1 if the signal on a tracker view is saturated. | 
| 317 |  | * @param ip plane (0-5) | 
| 318 |  | * @param iv view (0=x 1=y) | 
| 319 |  | */ | 
| 320 |  | Bool_t TrkTrack::IsSaturated(int ip,int iv){ | 
| 321 |  | if(iv==0 && ip>=0 && ip<6)return (dedx_x[ip]<0) ; | 
| 322 |  | else if(iv==1 && ip>=0 && ip<6)return (dedx_y[ip]<0) ; | 
| 323 |  | else { | 
| 324 |  | cout << "TrkTrack::IsSaturated(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; | 
| 325 |  | return 0.; | 
| 326 |  | } | 
| 327 |  | }; | 
| 328 |  | /** | 
| 329 |  | * Returns 1 if either the x or the y signal on a tracker plane is saturated. | 
| 330 |  | * @param ip plane (0-5) | 
| 331 |  | */ | 
| 332 |  | Bool_t TrkTrack::IsSaturated(int ip){ | 
| 333 |  | return (IsSaturated(ip,0)||IsSaturated(ip,1)); | 
| 334 |  | }; | 
| 335 |  | /** | 
| 336 |  | * Returns 1 if there is at least a saturated signal along the track. | 
| 337 |  | */ | 
| 338 |  | Bool_t TrkTrack::IsSaturated(){ | 
| 339 |  | for(int ip=0; ip<6; ip++)for(int iv=0; iv<2; iv++)if(IsSaturated(ip,iv))return true; | 
| 340 |  | return false; | 
| 341 |  | } | 
| 342 |  | /** | 
| 343 |  | * Returns the track "lever-arm" on the x view, defined as the distance (in planes) between | 
| 344 |  | * the upper and lower x measurements (the maximum value of lever-arm is 6). | 
| 345 |  | */ | 
| 346 |  | Int_t TrkTrack::GetLeverArmX(){ | 
| 347 |  | int first_plane = -1; | 
| 348 |  | int last_plane  = -1; | 
| 349 |  | for(Int_t ip=0; ip<6; ip++){ | 
| 350 |  | if( XGood(ip) && first_plane == -1 )first_plane = ip; | 
| 351 |  | if( XGood(ip) && first_plane != -1 )last_plane = ip; | 
| 352 |  | } | 
| 353 |  | if( first_plane == -1 || last_plane == -1){ | 
| 354 |  | cout<< "Int_t TrkTrack::GetLeverArmX() -- XGood(ip) always false ??? "<<endl; | 
| 355 |  | return 0; | 
| 356 |  | } | 
| 357 |  | return (last_plane-first_plane+1); | 
| 358 |  | } | 
| 359 |  | /** | 
| 360 |  | * Returns the track "lever-arm" on the y view, defined as the distance (in planes) between | 
| 361 |  | * the upper and lower y measurements (the maximum value of lever-arm is 6). | 
| 362 |  | */ | 
| 363 |  | Int_t TrkTrack::GetLeverArmY(){ | 
| 364 |  | int first_plane = -1; | 
| 365 |  | int last_plane  = -1; | 
| 366 |  | for(Int_t ip=0; ip<6; ip++){ | 
| 367 |  | if( YGood(ip) && first_plane == -1 )first_plane = ip; | 
| 368 |  | if( YGood(ip) && first_plane != -1 )last_plane = ip; | 
| 369 |  | } | 
| 370 |  | if( first_plane == -1 || last_plane == -1){ | 
| 371 |  | cout<< "Int_t TrkTrack::GetLeverArmY() -- YGood(ip) always false ??? "<<endl; | 
| 372 |  | return 0; | 
| 373 |  | } | 
| 374 |  | 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 | 
| 395 |  | */ | 
| 396 |  | Float_t  TrkTrack::GetChi2X(){ | 
| 397 |  | float chiq=0; | 
| 398 |  | for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.); | 
| 399 |  | if(GetNX()>3)chiq=chiq/(GetNX()-3); | 
| 400 |  | else chiq=0; | 
| 401 |  | if(chiq==0)cout << " Float_t  TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl; | 
| 402 |  | return chiq; | 
| 403 |  | } | 
| 404 |  | /** | 
| 405 |  | * Returns the reduced chi-square of track y-projection | 
| 406 |  | */ | 
| 407 |  | Float_t  TrkTrack::GetChi2Y(){ | 
| 408 |  | float chiq=0; | 
| 409 |  | for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.); | 
| 410 |  | if(GetNY()>2)chiq=chiq/(GetNY()-2); | 
| 411 |  | else chiq=0; | 
| 412 |  | if(chiq==0)cout << " Float_t  TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl; | 
| 413 |  | return chiq; | 
| 414 |  | } | 
| 415 |  | /** | 
| 416 |  | * Returns the logarythm of the likeliwood-function of  track x-projection | 
| 417 |  | */ | 
| 418 |  | Float_t TrkTrack::GetLnLX(){ | 
| 419 |  | float lnl=0; | 
| 420 |  | for(int ip=0; ip<6; ip++) | 
| 421 |  | if( XGood(ip) && tailx[ip]!=0 ) | 
| 422 |  | lnl += (tailx[ip]+1.) * log( (tailx[ip]*pow(resx[ip],2.) + pow(xv[ip]-xm[ip],2.)) / (tailx[ip]*pow(resx[ip],2)) ); | 
| 423 |  | if(GetNX()>3)lnl=lnl/(GetNX()-3); | 
| 424 |  | else lnl=0; | 
| 425 |  | if(lnl==0){ | 
| 426 |  | cout << " Float_t  TrkTrack::GetLnLX() -- WARNING -- value not defined "<<lnl<<endl; | 
| 427 |  | Dump(); | 
| 428 |  | } | 
| 429 |  | return lnl; | 
| 430 |  |  | 
| 431 |  | } | 
| 432 |  | /** | 
| 433 |  | * Returns the logarythm of the likeliwood-function of  track y-projection | 
| 434 |  | */ | 
| 435 |  | Float_t TrkTrack::GetLnLY(){ | 
| 436 |  | float lnl=0; | 
| 437 |  | for(int ip=0; ip<6; ip++) | 
| 438 |  | if( YGood(ip) && taily[ip]!=0 ) | 
| 439 |  | lnl += (taily[ip]+1.) * log( (taily[ip]*pow(resy[ip],2.) + pow(yv[ip]-ym[ip],2.)) / (taily[ip]*pow(resy[ip],2)) ); | 
| 440 |  | if(GetNY()>2)lnl=lnl/(GetNY()-2); | 
| 441 |  | else lnl=0; | 
| 442 |  | if(lnl==0){ | 
| 443 |  | cout << " Float_t  TrkTrack::GetLnLY() -- WARNING -- value not defined "<<lnl<<endl; | 
| 444 |  | Dump(); | 
| 445 |  | } | 
| 446 |  | return lnl; | 
| 447 |  |  | 
| 448 |  | } | 
| 449 |  | /** | 
| 450 |  | * Returns the effective angle, relative to the sensor, on each plane. | 
| 451 |  | * @param ip plane (0-5) | 
| 452 |  | * @param iv view (0=x 1=y) | 
| 453 |  | */ | 
| 454 |  | Float_t TrkTrack::GetEffectiveAngle(int ip, int iv){ | 
| 455 |  |  | 
| 456 |  | if(ip<0 || ip>5){ | 
| 457 |  | cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl; | 
| 458 |  | return 0.; | 
| 459 |  | } | 
| 460 |  |  | 
| 461 |  | float v[3]={xv[ip],yv[ip],zv[ip]}; | 
| 462 |  | //----------------------------------------- | 
| 463 |  | // effective angle (relative to the sensor) | 
| 464 |  | //----------------------------------------- | 
| 465 |  | float axv_geo  = axv[ip]; | 
| 466 |  | float muhall_h = 297.61; //cm**2/Vs | 
| 467 |  | float BY = TrkParams::GetBY(v); | 
| 468 |  | float axv_eff = 0; | 
| 469 |  | if(ip==5) axv_geo = -1*axv_geo; | 
| 470 |  | if(ip==5) BY      = -1*BY; | 
| 471 |  | axv_eff = 180.*atan( tan(axv_geo*acos(-1.)/180.) + muhall_h * BY * 0.0001)/acos(-1.); | 
| 472 |  | //----------------------------------------- | 
| 473 |  | // effective angle (relative to the sensor) | 
| 474 |  | //----------------------------------------- | 
| 475 |  | float ayv_geo = ayv[ip]; | 
| 476 |  | float muhall_e = 1258.18; //cm**2/Vs | 
| 477 |  | float BX = TrkParams::GetBX(v); | 
| 478 |  | float ayv_eff = 0; | 
| 479 |  | ayv_eff = 180.*atan( tan(ayv_geo*acos(-1.)/180.) + muhall_e * BX * 0.0001)/acos(-1.); | 
| 480 |  |  | 
| 481 |  | if     (iv==0)return axv_eff; | 
| 482 |  | else if(iv==1)return ayv_eff; | 
| 483 |  | else{ | 
| 484 |  | cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl; | 
| 485 |  | return 0.; | 
| 486 |  | } | 
| 487 |  |  | 
| 488 | }; | }; | 
| 489 |  |  | 
| 490 | //-------------------------------------- | //-------------------------------------- | 
| 498 | cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " "; | cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " "; | 
| 499 | cout << endl << "chi^2    : "<< chi2; | cout << endl << "chi^2    : "<< chi2; | 
| 500 | cout << endl << "n.step   : "<< nstep; | cout << endl << "n.step   : "<< nstep; | 
| 501 | cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ; | cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << XGood(i) ; | 
| 502 | cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ; | cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << YGood(i) ; | 
| 503 | cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " "; | cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " "; | 
| 504 | cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " "; | cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " "; | 
| 505 | cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " "; | cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " "; | 
| 508 | cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " "; | cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " "; | 
| 509 | cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " "; | cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " "; | 
| 510 | cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " "; | cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " "; | 
| 511 |  | cout << endl << "tailx    : "; for(int i=0; i<6; i++)cout << tailx[i] << " "; | 
| 512 |  | cout << endl << "taily    : "; for(int i=0; i<6; i++)cout << taily[i] << " "; | 
| 513 | cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" "; | cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" "; | 
| 514 | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" "; | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" "; | 
| 515 | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" "; | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" "; | 
| 517 | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" "; | cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" "; | 
| 518 | cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; | cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; | 
| 519 | cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; | cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; | 
| 520 |  | cout << endl << "maxs x   : "; for(int i=0; i<6; i++)cout << GetClusterX_MaxStrip(i) << " "; | 
| 521 |  | cout << endl << "maxs y   : "; for(int i=0; i<6; i++)cout << GetClusterY_MaxStrip(i) << " "; | 
| 522 |  | cout << endl << "mult x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Multiplicity(i) << " "; | 
| 523 |  | cout << endl << "mult y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Multiplicity(i) << " "; | 
| 524 |  | cout << endl << "seed x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Seed(i) << " "; | 
| 525 |  | cout << endl << "seed y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Seed(i) << " "; | 
| 526 |  | cout << endl << "xpu      : "; for(int i=0; i<6; i++)cout << xpu[i] << " "; | 
| 527 |  | cout << endl << "ypu      : "; for(int i=0; i<6; i++)cout << ypu[i] << " "; | 
| 528 |  |  | 
| 529 | cout << endl; | cout << endl; | 
| 530 | } | } | 
| 531 | /** | /** | 
| 544 | for(int i=0; i<6; i++) resy[i]=*ry++; | for(int i=0; i<6; i++) resy[i]=*ry++; | 
| 545 | } | } | 
| 546 | /** | /** | 
| 547 |  | * Set the TrkTrack tails position resolution | 
| 548 |  | */ | 
| 549 |  | void TrkTrack::SetTail(double *tx, double *ty, double factor){ | 
| 550 |  | for(int i=0; i<6; i++) tailx[i]=factor*(*tx++); | 
| 551 |  | for(int i=0; i<6; i++) taily[i]=factor*(*ty++); | 
| 552 |  | } | 
| 553 |  | /** | 
| 554 |  | * Set the TrkTrack Student parameter (resx,resy,tailx,taily) | 
| 555 |  | * from previous gausian fit | 
| 556 |  | *@param flag =0 standard, =1 with noise correction | 
| 557 |  | */ | 
| 558 |  | void TrkTrack::SetStudentParam(int flag){ | 
| 559 |  | float sx[11]={0.000128242, | 
| 560 |  | 0.000136942, | 
| 561 |  | 0.000162718, | 
| 562 |  | 0.000202644, | 
| 563 |  | 0.00025597, | 
| 564 |  | 0.000317456, | 
| 565 |  | 0.000349048, | 
| 566 |  | 0.000384638, | 
| 567 |  | 0.000457295, | 
| 568 |  | 0.000512319, | 
| 569 |  | 0.000538573}; | 
| 570 |  | float tx[11]={1.79402, | 
| 571 |  | 2.04876, | 
| 572 |  | 2.88376, | 
| 573 |  | 3.3, | 
| 574 |  | 3.14084, | 
| 575 |  | 4.07686, | 
| 576 |  | 4.44736, | 
| 577 |  | 3.5179, | 
| 578 |  | 3.38697, | 
| 579 |  | 3.45739, | 
| 580 |  | 3.18627}; | 
| 581 |  | float sy[11]={0.000483075, | 
| 582 |  | 0.000466925, | 
| 583 |  | 0.000431658, | 
| 584 |  | 0.000428317, | 
| 585 |  | 0.000433854, | 
| 586 |  | 0.000444044, | 
| 587 |  | 0.000482098, | 
| 588 |  | 0.000537579, | 
| 589 |  | 0.000636279, | 
| 590 |  | 0.000741998, | 
| 591 |  | 0.000864261}; | 
| 592 |  | float ty[11]={0.997032, | 
| 593 |  | 1.11147, | 
| 594 |  | 1.18526, | 
| 595 |  | 1.61404, | 
| 596 |  | 2.21908, | 
| 597 |  | 3.08959, | 
| 598 |  | 4.48833, | 
| 599 |  | 4.42687, | 
| 600 |  | 4.65253, | 
| 601 |  | 4.52043, | 
| 602 |  | 4.29926}; | 
| 603 |  | int index; | 
| 604 |  | float fact=0.; | 
| 605 |  | for(int i=0; i<6; i++) { | 
| 606 |  | index = int((fabs(axv[i])+1.)/2.); | 
| 607 |  | if(index>10) index=10; | 
| 608 |  | tailx[i]=tx[index]; | 
| 609 |  | if(flag==1) { | 
| 610 |  | if(fabs(axv[i])<=10.) fact = resx[i]/risxeta2_(&(axv[i])); | 
| 611 |  | if(fabs(axv[i])>10.&&fabs(axv[i])<=15.) fact = resx[i]/risxeta3_(&(axv[i])); | 
| 612 |  | if(fabs(axv[i])>15.) fact = resx[i]/risxeta4_(&(axv[i])); | 
| 613 |  | } else fact = 1.; | 
| 614 |  | resx[i] = sx[index]*fact; | 
| 615 |  | } | 
| 616 |  | for(int i=0; i<6; i++) { | 
| 617 |  | index = int((fabs(ayv[i])+1.)/2.); | 
| 618 |  | if(index>10) index=10; | 
| 619 |  | taily[i]=ty[index]; | 
| 620 |  | if(flag==1) fact = resy[i]/risyeta2_(&(ayv[i])); | 
| 621 |  | else fact = 1.; | 
| 622 |  | resy[i] = sy[index]*fact; | 
| 623 |  | } | 
| 624 |  | } | 
| 625 |  | /** | 
| 626 | * Set the TrkTrack good measurement | * Set the TrkTrack good measurement | 
| 627 | */ | */ | 
| 628 | void TrkTrack::SetGood(int *xg, int *yg){ | void TrkTrack::SetGood(int *xg, int *yg){ | 
| 629 |  |  | 
| 630 | for(int i=0; i<6; i++) xgood[i]=*xg++; | for(int i=0; i<6; i++) xgood[i]=*xg++; | 
| 631 | for(int i=0; i<6; i++) ygood[i]=*yg++; | for(int i=0; i<6; i++) ygood[i]=*yg++; | 
| 632 | } | } | 
| 636 | */ | */ | 
| 637 | void TrkTrack::LoadField(TString path){ | void TrkTrack::LoadField(TString path){ | 
| 638 |  |  | 
| 639 | strcpy(path_.path,path.Data()); | //     strcpy(path_.path,path.Data()); | 
| 640 | path_.pathlen = path.Length(); | //     path_.pathlen = path.Length(); | 
| 641 | path_.error   = 0; | //     path_.error   = 0; | 
| 642 | readb_(); | //     readb_(); | 
| 643 |  |  | 
| 644 |  | //     TrkParams::SetTrackingMode(); | 
| 645 |  | //     TrkParams::SetPrecisionFactor(); | 
| 646 |  | //     TrkParams::SetStepMin(); | 
| 647 |  | TrkParams::SetMiniDefault(); | 
| 648 |  |  | 
| 649 |  | TrkParams::Set(path,1); | 
| 650 |  | 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 |  |  | 
| 657 |  |  | 
| 658 | /** | /** | 
| 659 | * Method to fill minimization-routine common | * Method to fill minimization-routine common | 
| 660 | */ | */ | 
| 662 |  |  | 
| 663 | for(int i=0; i<6; i++){ | for(int i=0; i<6; i++){ | 
| 664 |  |  | 
| 665 | track.xgood[i]=xgood[i]; | //      cout << i<<" - "<<xgood[i]<<" "<<XGood(i)<<endl; | 
| 666 | track.ygood[i]=ygood[i]; | //      cout << i<<" - "<<ygood[i]<<" "<<YGood(i)<<endl; | 
| 667 |  | track.xgood[i]=XGood(i); | 
| 668 |  | track.ygood[i]=YGood(i); | 
| 669 |  |  | 
| 670 | track.xm[i]=xm[i]; | track.xm[i]=xm[i]; | 
| 671 | track.ym[i]=ym[i]; | track.ym[i]=ym[i]; | 
| 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.;// EM GCC 4.7 | 
| 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];// EM GCC 4.7 unused | 
| 712 |  | //              gamma   = alignparameters_.gamma[is][il][ip];// EM GCC 4.7 unused | 
| 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]; | 
| 735 |  | track.taily[i]=taily[i]; | 
| 736 | } | } | 
| 737 |  |  | 
| 738 | for(int i=0; i<5; i++) track.al[i]=al[i]; | for(int i=0; i<5; i++) track.al[i]=al[i]; | 
| 761 | axv[i] = track->axv[i]; | axv[i] = track->axv[i]; | 
| 762 | ayv[i] = track->ayv[i]; | ayv[i] = track->ayv[i]; | 
| 763 | } | } | 
|  |  |  | 
| 764 |  |  | 
| 765 | } | } | 
| 766 | /** | /** | 
| 767 | * Tracking method. It calls F77 mini routine. | * \brief Method to re-evaluate coordinates of clusters associated with a track. | 
| 768 |  | * | 
| 769 |  | * The method can be applied only after recovering level1 information | 
| 770 |  | * (either by reprocessing single events from level0 or from | 
| 771 |  | * the TrkLevel1 branch, if present); it calls F77 subroutines that | 
| 772 |  | * read the level1 common and fill the minimization-routine common. | 
| 773 |  | * Some clusters can be excluded or added by means of the methods: | 
| 774 |  | * | 
| 775 |  | * TrkTrack::ResetXGood(int ip) | 
| 776 |  | * TrkTrack::ResetYGood(int ip) | 
| 777 |  | * TrkTrack::SetXGood(int ip, int cid, int is) | 
| 778 |  | * TrkTrack::SetYGood(int ip, int cid, int is) | 
| 779 |  | * | 
| 780 |  | * NB! The method TrkTrack::SetGood(int *xg, int *yg) set the plane-mask (0-1) | 
| 781 |  | * for the minimization-routine common. It deletes the cluster information | 
| 782 |  | * (at least for the moment...) thus cannot be applied before | 
| 783 |  | * TrkTrack::EvaluateClusterPositions(). | 
| 784 |  | * | 
| 785 |  | * Different p.f.a. can be applied by calling (once) the method: | 
| 786 |  | * | 
| 787 |  | * TrkParams::SetPFA(0); //Set ETA p.f.a. | 
| 788 |  | * | 
| 789 |  | * @see TrkParams::SetPFA(int) | 
| 790 | */ | */ | 
| 791 | void TrkTrack::Fit(double pfixed, int& fail, int iprint){ | Bool_t TrkTrack::EvaluateClusterPositions(){ | 
| 792 |  |  | 
| 793 |  | //     cout << "void TrkTrack::GetClusterositions() "<<endl; | 
| 794 |  |  | 
| 795 |  | TrkParams::Load(1); | 
| 796 |  | 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++){ | 
| 812 |  | //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;; | 
| 813 |  | int icx = GetClusterX_ID(ip)+1; | 
| 814 |  | int icy = GetClusterY_ID(ip)+1; | 
| 815 |  | int sensor = GetSensor(ip)+1;//<< convenzione "Paolo" | 
| 816 |  | if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena" | 
| 817 |  | int ladder = GetLadder(ip)+1; | 
| 818 |  | float ax = axv[ip]; | 
| 819 |  | float ay = ayv[ip]; | 
| 820 |  | float v[3]; | 
| 821 |  | v[0]=xv[ip]; | 
| 822 |  | v[1]=yv[ip]; | 
| 823 |  | v[2]=zv[ip]; | 
| 824 |  | float bfx = 10*TrkParams::GetBX(v);//Tesla | 
| 825 |  | float bfy = 10*TrkParams::GetBY(v);//Tesla | 
| 826 |  | int ipp=ip+1; | 
| 827 |  | xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy); | 
| 828 |  | if(icx<0 || icy<0)return false; | 
| 829 |  | } | 
| 830 |  | return true; | 
| 831 |  | } | 
| 832 |  | /** | 
| 833 |  | * \brief Tracking method. It calls F77 mini routine. | 
| 834 |  | * | 
| 835 |  | * @param pfixed Particle momentum. If pfixed=0 the momentum | 
| 836 |  | * is left as a free parameter, otherwise it is fixed to the input value. | 
| 837 |  | * @param fail Output flag (!=0 if the fit failed). | 
| 838 |  | * @param iprint Flag to set debug mode ( 0 = no output; 1 = verbose; 2 = debug). | 
| 839 |  | * @param froml1 Flag to re-evaluate positions (see TrkTrack::GetClusterPositions()). | 
| 840 |  | * | 
| 841 |  | * The option to re-evaluate positions can be used only after recovering | 
| 842 |  | * level1 information, eg. by reprocessing the single event. | 
| 843 |  | * | 
| 844 |  | * Example: | 
| 845 |  | * | 
| 846 |  | * if( !event->GetTrkLevel0() )return false; | 
| 847 |  | * event->GetTrkLevel0()->ProcessEvent(); // re-processing level0->level1 | 
| 848 |  | * int fail=0; | 
| 849 |  | * event->GetTrkLevel2()->GetTrack(0)->Fit(0.,fail,0,1); | 
| 850 |  | * | 
| 851 |  | * @see EvaluateClusterPositions() | 
| 852 |  | * | 
| 853 |  | * The fitting procedure can be varied by changing the tracking mode, | 
| 854 |  | * the fit-precision factor, the minimum number of step, etc. | 
| 855 |  | * @see SetTrackingMode(int) | 
| 856 |  | * @see SetPrecisionFactor(double) | 
| 857 |  | * @see SetStepMin(int) | 
| 858 |  | * @see SetDeltaB(int,double) | 
| 859 |  | */ | 
| 860 |  | void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){ | 
| 861 |  |  | 
| 862 |  | 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 | float al_ini[] = {0.,0.,0.,0.,0.}; | float al_ini[] = {0.,0.,0.,0.,0.}; | 
| 874 |  |  | 
| 875 | extern cMini2track track_; | extern cMini2track track_; | 
| 876 | fail = 0; | fail = 0; | 
|  | FillMiniStruct(track_); |  | 
| 877 |  |  | 
| 878 |  | FillMiniStruct(track_); | 
| 879 |  |  | 
| 880 |  | if(froml1!=0){ | 
| 881 |  | if( !EvaluateClusterPositions() ){ | 
| 882 |  | cout << "void TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<") --- ERROR evaluating cluster positions "<<endl; | 
| 883 |  | FillMiniStruct(track_) ; | 
| 884 |  | fail = 1; | 
| 885 |  | return; | 
| 886 |  | } | 
| 887 |  | }else{ | 
| 888 |  | FillMiniStruct(track_); | 
| 889 |  | } | 
| 890 |  |  | 
| 891 | // if fit variables have been reset, evaluate the initial guess | // if fit variables have been reset, evaluate the initial guess | 
| 892 | if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_(); | if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_(); | 
| 893 |  |  | 
| 906 |  |  | 
| 907 | //  ------------------------------------------ | //  ------------------------------------------ | 
| 908 | //  call mini routine | //  call mini routine | 
| 909 |  | //  ------------------------------------------ | 
| 910 | int istep=0; | int istep=0; | 
| 911 | int ifail=0; | int ifail=0; | 
| 912 | mini2_(&istep,&ifail, &iprint); | mini2_(&istep,&ifail, &iprint); | 
| 917 | //  ------------------------------------------ | //  ------------------------------------------ | 
| 918 |  |  | 
| 919 | SetFromMiniStruct(&track_); | SetFromMiniStruct(&track_); | 
|  | //    cout << endl << "eta ===> " << track_.al[4] << endl; |  | 
|  |  |  | 
|  | //     for(int i=0; i<5; i++) al[i]=track_.al[i]; |  | 
|  | //     chi2=track_.chi2; |  | 
|  | //     nstep=track_.nstep; |  | 
|  | //     for(int i=0; i<6; i++) xv[i]=track_.xv[i]; |  | 
|  | //     for(int i=0; i<6; i++) yv[i]=track_.yv[i]; |  | 
|  | //     for(int i=0; i<6; i++) zv[i]=track_.zv[i]; |  | 
|  | //     for(int i=0; i<6; i++) axv[i]=track_.axv[i]; |  | 
|  | //     for(int i=0; i<6; i++) ayv[i]=track_.ayv[i]; |  | 
|  | //     for(int i=0; i<5; i++) { |  | 
|  | //      for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j]; |  | 
|  | //     } |  | 
| 920 |  |  | 
| 921 | if(fail){ | if(fail){ | 
| 922 | if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl; | if(iprint)cout << " >>>> fit failed "<<endl; | 
| 923 | for(int i=0; i<5; i++) al[i]=al_ini[i]; | for(int i=0; i<5; i++) al[i]=al_ini[i]; | 
| 924 | } | } | 
| 925 |  |  | 
| 926 | }; | }; | 
| 927 | /* | /** | 
| 928 | * Reset the fit parameters | * Reset the fit parameters | 
| 929 | */ | */ | 
| 930 | void TrkTrack::FitReset(){ | void TrkTrack::FitReset(){ | 
| 931 | for(int i=0; i<5; i++) al[i]=-9999.; | for(int i=0; i<5; i++) al[i]=-9999.; | 
| 932 | chi2=0.; | chi2=0.; | 
| 933 | nstep=0; | nstep=0; | 
| 934 | for(int i=0; i<6; i++) xv[i]=0.; | //     for(int i=0; i<6; i++) xv[i]=0.; | 
| 935 | for(int i=0; i<6; i++) yv[i]=0.; | //     for(int i=0; i<6; i++) yv[i]=0.; | 
| 936 | for(int i=0; i<6; i++) zv[i]=0.; | //     for(int i=0; i<6; i++) zv[i]=0.; | 
| 937 | for(int i=0; i<6; i++) axv[i]=0.; | //     for(int i=0; i<6; i++) axv[i]=0.; | 
| 938 | for(int i=0; i<6; i++) ayv[i]=0.; | //     for(int i=0; i<6; i++) ayv[i]=0.; | 
| 939 | for(int i=0; i<5; i++) { | for(int i=0; i<5; i++) { | 
| 940 | for(int j=0; j<5; j++) coval[i][j]=0.; | for(int j=0; j<5; j++) coval[i][j]=0.; | 
| 941 | } | } | 
| 942 | } | } | 
| 943 |  | /** | 
| 944 |  | * Set the tracking mode | 
| 945 |  | */ | 
| 946 |  | void TrkTrack::SetTrackingMode(int trackmode){ | 
| 947 |  | extern cMini2track track_; | 
| 948 |  | track_.trackmode = trackmode; | 
| 949 |  | } | 
| 950 |  | /** | 
| 951 |  | * Set the factor scale for tracking precision | 
| 952 |  | */ | 
| 953 |  | void TrkTrack::SetPrecisionFactor(double fact){ | 
| 954 |  | extern cMini2track track_; | 
| 955 |  | track_.fact = fact; | 
| 956 |  | } | 
| 957 |  | /** | 
| 958 |  | * Set the minimum number of steps for tracking precision | 
| 959 |  | */ | 
| 960 |  | void TrkTrack::SetStepMin(int istepmin){ | 
| 961 |  | extern cMini2track track_; | 
| 962 |  | track_.istepmin = istepmin; | 
| 963 |  | } | 
| 964 |  | /** | 
| 965 |  | * Set deltaB parameters (id=0,1). By default they are set to zero. | 
| 966 |  | */ | 
| 967 |  | void TrkTrack::SetDeltaB(int id, double db){ | 
| 968 |  | if(id!=0 && id!=1)cout << "void TrkTrack::SetDeltaB(int id,double db) -- wrong input parameters: "<<id<<" "<<db<<endl; | 
| 969 |  | TrkParams::SetDeltaB(id,db); | 
| 970 |  | } | 
| 971 |  |  | 
| 972 |  | /** | 
| 973 |  | * Returns true if the track is inside the magnet cavity. | 
| 974 |  | * @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume) | 
| 975 |  | */ | 
| 976 |  | Bool_t TrkTrack::IsInsideCavity(float toll){ | 
| 977 |  |  | 
| 978 |  | //     float xmagntop, ymagntop, xmagnbottom, ymagnbottom; | 
| 979 |  | //     xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.); | 
| 980 |  | //     ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.); | 
| 981 |  | //     xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.); | 
| 982 |  | //     ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.); | 
| 983 |  | //     if( xmagntop>XMAGNLOW && xmagntop<XMAGNHIGH && | 
| 984 |  | //      ymagntop>YMAGNLOW && ymagntop<YMAGNHIGH && | 
| 985 |  | //      xmagnbottom>XMAGNLOW && xmagnbottom<XMAGNHIGH && | 
| 986 |  | //      ymagnbottom>YMAGNLOW && ymagnbottom<YMAGNHIGH ) return(true); | 
| 987 |  | //     else return(false); | 
| 988 |  |  | 
| 989 |  | int ngf = TrkParams::nGF; | 
| 990 |  | for(int i=0; i<ngf; i++){ | 
| 991 |  | // | 
| 992 |  | //      cout << endl << TrkParams::GF_element[i]; | 
| 993 |  | if( | 
| 994 |  | TrkParams::GF_element[i].CompareTo("CUF") && | 
| 995 |  | TrkParams::GF_element[i].CompareTo("T2")  && | 
| 996 |  | TrkParams::GF_element[i].CompareTo("T3")  && | 
| 997 |  | TrkParams::GF_element[i].CompareTo("T4")  && | 
| 998 |  | TrkParams::GF_element[i].CompareTo("T5")  && | 
| 999 |  | TrkParams::GF_element[i].CompareTo("CLF") && | 
| 1000 |  | true)continue; | 
| 1001 |  | // apply condition only within the cavity | 
| 1002 |  | //      cout << " -- "<<xGF[i]<<" "<<yGF[i]; | 
| 1003 |  | if( | 
| 1004 |  | xGF[i] <= TrkParams::xGF_min[i] + toll || | 
| 1005 |  | xGF[i] >= TrkParams::xGF_max[i] - toll || | 
| 1006 |  | yGF[i] <= TrkParams::yGF_min[i] + toll || | 
| 1007 |  | yGF[i] >= TrkParams::yGF_max[i] - toll || | 
| 1008 |  | false){ | 
| 1009 |  |  | 
| 1010 |  | return false; | 
| 1011 |  | } | 
| 1012 |  | } | 
| 1013 |  | return true; | 
| 1014 |  |  | 
| 1015 |  |  | 
| 1016 |  | } | 
| 1017 |  | /** | 
| 1018 |  | * Returns true if the track is inside the nominal acceptance, which is defined | 
| 1019 |  | * by the intersection among magnet cavity, silicon-plane sensitive area and | 
| 1020 |  | * ToF sensitive area (nominal values from the official document used to | 
| 1021 |  | * calculate the geometrical factor) | 
| 1022 |  | * @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume) | 
| 1023 |  | */ | 
| 1024 |  | // Bool_t TrkTrack::IsInsideAcceptance(){ | 
| 1025 |  |  | 
| 1026 |  | //     int ngf = TrkParams::nGF; | 
| 1027 |  | //     for(int i=0; i<ngf; i++){ | 
| 1028 |  | //      if( | 
| 1029 |  | //          xGF[i] <= TrkParams::xGF_min[i] || | 
| 1030 |  | //          xGF[i] >= TrkParams::xGF_max[i] || | 
| 1031 |  | //          yGF[i] <= TrkParams::yGF_min[i] || | 
| 1032 |  | //          yGF[i] >= TrkParams::yGF_max[i] || | 
| 1033 |  | //          false)return false; | 
| 1034 |  | //     } | 
| 1035 |  | //     return true; | 
| 1036 |  |  | 
| 1037 |  | // } | 
| 1038 |  | Bool_t TrkTrack::IsInsideAcceptance(float toll){ | 
| 1039 |  |  | 
| 1040 |  |  | 
| 1041 |  | int ngf = TrkParams::nGF; | 
| 1042 |  | for(int i=0; i<ngf; i++){ | 
| 1043 |  | // | 
| 1044 |  | //      cout << endl << TrkParams::GF_element[i]; | 
| 1045 |  | if( | 
| 1046 |  | TrkParams::GF_element[i].CompareTo("S11") && | 
| 1047 |  | TrkParams::GF_element[i].CompareTo("S12") && | 
| 1048 |  | TrkParams::GF_element[i].CompareTo("S21") && | 
| 1049 |  | TrkParams::GF_element[i].CompareTo("S22") && | 
| 1050 |  | TrkParams::GF_element[i].CompareTo("T1")  && | 
| 1051 |  | TrkParams::GF_element[i].CompareTo("CUF") && | 
| 1052 |  | TrkParams::GF_element[i].CompareTo("T2")  && | 
| 1053 |  | TrkParams::GF_element[i].CompareTo("T3")  && | 
| 1054 |  | TrkParams::GF_element[i].CompareTo("T4")  && | 
| 1055 |  | TrkParams::GF_element[i].CompareTo("T5")  && | 
| 1056 |  | TrkParams::GF_element[i].CompareTo("CLF") && | 
| 1057 |  | TrkParams::GF_element[i].CompareTo("T6")  && | 
| 1058 |  | TrkParams::GF_element[i].CompareTo("S31") && | 
| 1059 |  | TrkParams::GF_element[i].CompareTo("S32") && | 
| 1060 |  | true)continue; | 
| 1061 |  | // apply condition only within the cavity | 
| 1062 |  | //      cout << " -- "<<xGF[i]<<" "<<yGF[i]; | 
| 1063 |  | if( | 
| 1064 |  | xGF[i] <= TrkParams::xGF_min[i] + toll || | 
| 1065 |  | xGF[i] >= TrkParams::xGF_max[i] - toll || | 
| 1066 |  | yGF[i] <= TrkParams::yGF_min[i] + toll || | 
| 1067 |  | yGF[i] >= TrkParams::yGF_max[i] - toll || | 
| 1068 |  | false){ | 
| 1069 |  |  | 
| 1070 |  | return false; | 
| 1071 |  | } | 
| 1072 |  | } | 
| 1073 |  | return true; | 
| 1074 |  | } | 
| 1075 |  |  | 
| 1076 |  | /** | 
| 1077 |  | * Returns true if the track is inside one of the surfaces which define the | 
| 1078 |  | * geometrical acceptance. | 
| 1079 |  | * @param surf tag of the surface (possible values are: S11 S12 S21 S22 T1 | 
| 1080 |  | * CUF T2 T3 T4 T5 CLF T6 S31 S32). | 
| 1081 |  | * @param toll  Tolerance around the nominal surface (toll>0 define an inner | 
| 1082 |  | * fiducial surface) | 
| 1083 |  | */ | 
| 1084 |  | Bool_t TrkTrack::IsInsideGFSurface(const char* surf, float toll){ | 
| 1085 |  |  | 
| 1086 |  |  | 
| 1087 |  | int ngf = TrkParams::nGF; | 
| 1088 |  | bool SURFOK = false; | 
| 1089 |  | for(int i=0; i<ngf; i++){ | 
| 1090 |  | if(  !TrkParams::GF_element[i].CompareTo(surf)  ){ | 
| 1091 |  | SURFOK=true; | 
| 1092 |  | if( | 
| 1093 |  | xGF[i] > TrkParams::xGF_min[i] + toll && | 
| 1094 |  | xGF[i] < TrkParams::xGF_max[i] - toll && | 
| 1095 |  | yGF[i] > TrkParams::yGF_min[i] + toll && | 
| 1096 |  | yGF[i] < TrkParams::yGF_max[i] - toll && | 
| 1097 |  | true)return true; | 
| 1098 |  | } | 
| 1099 |  | } | 
| 1100 |  | if( !SURFOK )cout << " Bool_t TrkTrack::IsInsideGFSurface(char* surf, float toll) --> suface "<<surf<<" not defined "<<endl; | 
| 1101 |  | return false; | 
| 1102 |  |  | 
| 1103 |  | } | 
| 1104 |  |  | 
| 1105 |  | /** | 
| 1106 |  | * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track. | 
| 1107 |  | * If no cluster is associated, ID=-1. | 
| 1108 |  | * @param ip Tracker plane (0-5) | 
| 1109 |  | */ | 
| 1110 |  | Int_t TrkTrack::GetClusterX_ID(int ip){ | 
| 1111 |  | return ((Int_t)fabs(xgood[ip]))%10000000-1; | 
| 1112 |  | }; | 
| 1113 |  | /** | 
| 1114 |  | * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track. | 
| 1115 |  | * If no cluster is associated, ID=-1. | 
| 1116 |  | * @param ip Tracker plane (0-5) | 
| 1117 |  | */ | 
| 1118 |  | Int_t TrkTrack::GetClusterY_ID(int ip){ | 
| 1119 |  | return ((Int_t)fabs(ygood[ip]))%10000000-1; | 
| 1120 |  | }; | 
| 1121 |  |  | 
| 1122 |  | /** | 
| 1123 |  | * Method to retrieve the ladder (0-2, increasing x) traversed by the track on this plane. | 
| 1124 |  | * If no ladder is traversed (dead area) the metod retuns -1. | 
| 1125 |  | * @param ip Tracker plane (0-5) | 
| 1126 |  | */ | 
| 1127 |  | Int_t TrkTrack::GetLadder(int ip){ | 
| 1128 |  | if(XGood(ip))return (Int_t)fabs(xgood[ip]/100000000)-1; | 
| 1129 |  | if(YGood(ip))return (Int_t)fabs(ygood[ip]/100000000)-1; | 
| 1130 |  | return -1; | 
| 1131 |  | }; | 
| 1132 |  | /** | 
| 1133 |  | * Method to retrieve the sensor (0-1, increasing y) traversed by the track on this plane. | 
| 1134 |  | * If no sensor is traversed (dead area) the metod retuns -1. | 
| 1135 |  | * @param ip Tracker plane (0-5) | 
| 1136 |  | */ | 
| 1137 |  | Int_t TrkTrack::GetSensor(int ip){ | 
| 1138 |  | if(XGood(ip))return (Int_t)((Int_t)fabs(xgood[ip]/10000000)%10)-1; | 
| 1139 |  | if(YGood(ip))return (Int_t)((Int_t)fabs(ygood[ip]/10000000)%10)-1; | 
| 1140 |  | return -1; | 
| 1141 |  | }; | 
| 1142 |  |  | 
| 1143 |  | /** | 
| 1144 |  | * \brief Method to include a x-cluster to the track. | 
| 1145 |  | * @param ip Tracker plane (0-5) | 
| 1146 |  | * @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) | 
| 1147 |  | * @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) | 
| 1148 |  | * @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) | 
| 1149 |  | * @param bad True if the cluster contains bad strips | 
| 1150 |  | * @see Fit(double pfixed, int& fail, int iprint, int froml1) | 
| 1151 |  | */ | 
| 1152 |  | void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){ | 
| 1153 |  | //    int il=0;       //ladder (temporary) | 
| 1154 |  | //    bool bad=false; //ladder (temporary) | 
| 1155 |  | if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1) | 
| 1156 |  | cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; | 
| 1157 |  | xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid; | 
| 1158 |  | if(bad)xgood[ip]=-xgood[ip]; | 
| 1159 |  | }; | 
| 1160 |  | /** | 
| 1161 |  | * \brief Method to include a y-cluster to the track. | 
| 1162 |  | * @param ip Tracker plane (0-5) | 
| 1163 |  | * @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) | 
| 1164 |  | * @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) | 
| 1165 |  | * @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) | 
| 1166 |  | * @param bad True if the cluster contains bad strips | 
| 1167 |  | * @see Fit(double pfixed, int& fail, int iprint, int froml1) | 
| 1168 |  | */ | 
| 1169 |  | void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){ | 
| 1170 |  | //    int il=0;       //ladder (temporary) | 
| 1171 |  | //    bool bad=false; //ladder (temporary) | 
| 1172 |  | if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1) | 
| 1173 |  | cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; | 
| 1174 |  | ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid; | 
| 1175 |  | if(bad)ygood[ip]=-ygood[ip]; | 
| 1176 |  | }; | 
| 1177 |  |  | 
| 1178 |  | /** | 
| 1179 |  | * \brief Average X | 
| 1180 |  | * Average value of <xv>, evaluated from the first to the last hit x view. | 
| 1181 |  | */ | 
| 1182 |  | Float_t TrkTrack::GetXav(){ | 
| 1183 |  |  | 
| 1184 |  | int first_plane = -1; | 
| 1185 |  | int last_plane  = -1; | 
| 1186 |  | for(Int_t ip=0; ip<6; ip++){ | 
| 1187 |  | if( XGood(ip) && first_plane == -1 )first_plane = ip; | 
| 1188 |  | if( XGood(ip) && first_plane != -1 )last_plane = ip; | 
| 1189 |  | } | 
| 1190 |  | if( first_plane == -1 || last_plane == -1){ | 
| 1191 |  | return -100; | 
| 1192 |  | } | 
| 1193 |  | if( last_plane-first_plane+1 ==0 )return -100; | 
| 1194 |  |  | 
| 1195 |  | Float_t av = 0; | 
| 1196 |  | for(int ip=first_plane; ip<=last_plane; ip++)av+=xv[ip]; | 
| 1197 |  |  | 
| 1198 |  | return (av/(last_plane-first_plane+1)); | 
| 1199 |  | } | 
| 1200 |  | /** | 
| 1201 |  | * \brief Average Y | 
| 1202 |  | * Average value of <yv>, evaluated from the first to the last hit x view. | 
| 1203 |  | */ | 
| 1204 |  | Float_t TrkTrack::GetYav(){ | 
| 1205 |  |  | 
| 1206 |  | int first_plane = -1; | 
| 1207 |  | int last_plane  = -1; | 
| 1208 |  | for(Int_t ip=0; ip<6; ip++){ | 
| 1209 |  | if( XGood(ip) && first_plane == -1 )first_plane = ip; | 
| 1210 |  | if( XGood(ip) && first_plane != -1 )last_plane = ip; | 
| 1211 |  | } | 
| 1212 |  | if( first_plane == -1 || last_plane == -1){ | 
| 1213 |  | return -100; | 
| 1214 |  | } | 
| 1215 |  | if( last_plane-first_plane+1 ==0 )return -100; | 
| 1216 |  |  | 
| 1217 |  | Float_t av = 0; | 
| 1218 |  | for(int ip=first_plane; ip<=last_plane; ip++)av+=yv[ip]; | 
| 1219 |  |  | 
| 1220 |  | return (av/(last_plane-first_plane+1)); | 
| 1221 |  | } | 
| 1222 |  | /** | 
| 1223 |  | * \brief Average Z | 
| 1224 |  | * Average value of <zv>, evaluated from the first to the last hit x view. | 
| 1225 |  | */ | 
| 1226 |  | Float_t TrkTrack::GetZav(){ | 
| 1227 |  |  | 
| 1228 |  | int first_plane = -1; | 
| 1229 |  | int last_plane  = -1; | 
| 1230 |  | for(Int_t ip=0; ip<6; ip++){ | 
| 1231 |  | if( XGood(ip) && first_plane == -1 )first_plane = ip; | 
| 1232 |  | if( XGood(ip) && first_plane != -1 )last_plane = ip; | 
| 1233 |  | } | 
| 1234 |  | if( first_plane == -1 || last_plane == -1){ | 
| 1235 |  | return -100; | 
| 1236 |  | } | 
| 1237 |  | if( last_plane-first_plane+1 ==0 )return -100; | 
| 1238 |  |  | 
| 1239 |  | Float_t av = 0; | 
| 1240 |  | for(int ip=first_plane; ip<=last_plane; ip++)av+=zv[ip]; | 
| 1241 |  |  | 
| 1242 |  | return (av/(last_plane-first_plane+1)); | 
| 1243 |  | } | 
| 1244 |  |  | 
| 1245 |  | /** | 
| 1246 |  | * \brief Number of column traversed | 
| 1247 |  | */ | 
| 1248 |  | Int_t TrkTrack::GetNColumns(){ | 
| 1249 |  | int sensors[] = {0,0,0,0,0,0}; | 
| 1250 |  | for(int ip=0; ip<6; ip++){ | 
| 1251 |  | int sensorid = GetLadder(ip)+3*GetSensor(ip); | 
| 1252 |  | if(XGood(ip)||YGood(ip)) | 
| 1253 |  | if(sensorid>=0 && sensorid<6)sensors[sensorid]=1; | 
| 1254 |  | } | 
| 1255 |  | int nsensors=0; | 
| 1256 |  | for(int is=0; is<6; is++)nsensors += sensors[is]; | 
| 1257 |  | return nsensors; | 
| 1258 |  | }; | 
| 1259 |  | /** | 
| 1260 |  | * \brief Give the maximum energy release | 
| 1261 |  | */ | 
| 1262 |  | Float_t TrkTrack::GetDEDX_max(int ip, int iv){ | 
| 1263 |  | Float_t max=0; | 
| 1264 |  | int pfrom = 0; | 
| 1265 |  | int pto   = 6; | 
| 1266 |  | int vfrom = 0; | 
| 1267 |  | int vto   = 2; | 
| 1268 |  | if(ip>=0&&ip<6){ | 
| 1269 |  | pfrom = ip; | 
| 1270 |  | pto   = ip+1; | 
| 1271 |  | } | 
| 1272 |  | if(iv>=0&&iv<2){ | 
| 1273 |  | vfrom = iv; | 
| 1274 |  | vto   = iv+1; | 
| 1275 |  | } | 
| 1276 |  | for(int i=pfrom; i<pto; i++) | 
| 1277 |  | for(int j=vfrom; j<vto; j++){ | 
| 1278 |  | if(j==0 && XGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); | 
| 1279 |  | if(j==1 && YGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); | 
| 1280 |  | } | 
| 1281 |  | return max; | 
| 1282 |  |  | 
| 1283 |  | }; | 
| 1284 |  |  | 
| 1285 |  | /** | 
| 1286 |  | * \brief Give the minimum energy release | 
| 1287 |  | */ | 
| 1288 |  | Float_t TrkTrack::GetDEDX_min(int ip, int iv){ | 
| 1289 |  | Float_t min=100000000; | 
| 1290 |  | int pfrom = 0; | 
| 1291 |  | int pto   = 6; | 
| 1292 |  | int vfrom = 0; | 
| 1293 |  | int vto   = 2; | 
| 1294 |  | if(ip>=0&&ip<6){ | 
| 1295 |  | pfrom = ip; | 
| 1296 |  | pto   = ip+1; | 
| 1297 |  | } | 
| 1298 |  | if(iv>=0&&iv<2){ | 
| 1299 |  | vfrom = iv; | 
| 1300 |  | vto   = iv+1; | 
| 1301 |  | } | 
| 1302 |  | for(int i=pfrom; i<pto; i++) | 
| 1303 |  | for(int j=vfrom; j<vto; j++){ | 
| 1304 |  | if(j==0 && XGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); | 
| 1305 |  | if(j==1 && YGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); | 
| 1306 |  | } | 
| 1307 |  | return min; | 
| 1308 |  |  | 
| 1309 |  | }; | 
| 1310 |  |  | 
| 1311 |  | /** | 
| 1312 |  | * \brief Give the maximum spatial residual | 
| 1313 |  | */ | 
| 1314 |  | Float_t TrkTrack::GetResidual_max(int ip, int iv){ | 
| 1315 |  | Float_t max=0; | 
| 1316 |  | int pfrom = 0; | 
| 1317 |  | int pto   = 6; | 
| 1318 |  | int vfrom = 0; | 
| 1319 |  | int vto   = 2; | 
| 1320 |  | if(ip>=0&&ip<6){ | 
| 1321 |  | pfrom = ip; | 
| 1322 |  | pto   = ip+1; | 
| 1323 |  | } | 
| 1324 |  | if(iv>=0&&iv<2){ | 
| 1325 |  | vfrom = iv; | 
| 1326 |  | vto   = iv+1; | 
| 1327 |  | } | 
| 1328 |  | for(int i=pfrom; i<pto; i++){ | 
| 1329 |  | for(int j=vfrom; j<vto; j++){ | 
| 1330 |  | if(j==0 && XGood(i) && fabs(xm[i]-xv[i])>fabs(max))max=xm[i]-xv[i]; | 
| 1331 |  | if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=ym[i]-yv[i]; | 
| 1332 |  | } | 
| 1333 |  | } | 
| 1334 |  | return max; | 
| 1335 |  |  | 
| 1336 |  | }; | 
| 1337 |  | /** | 
| 1338 |  | * \brief Give the anerage spatial residual | 
| 1339 |  | */ | 
| 1340 |  | Float_t TrkTrack::GetResidual_av(int ip, int iv){ | 
| 1341 |  | // | 
| 1342 |  | //Sum$((xm>-50)*(xm-xv)/resx)/sqrt(TrkTrack.GetNX()*TrkTrack.GetChi2X())<0.3 | 
| 1343 |  |  | 
| 1344 |  | Float_t av  = 0.; | 
| 1345 |  | int     nav = 0; | 
| 1346 |  | // | 
| 1347 |  | int pfrom = 0; | 
| 1348 |  | int pto   = 6; | 
| 1349 |  | int vfrom = 0; | 
| 1350 |  | int vto   = 2; | 
| 1351 |  | if(ip>=0&&ip<6){ | 
| 1352 |  | pfrom = ip; | 
| 1353 |  | pto   = ip+1; | 
| 1354 |  | } | 
| 1355 |  | if(iv>=0&&iv<2){ | 
| 1356 |  | vfrom = iv; | 
| 1357 |  | vto   = iv+1; | 
| 1358 |  | } | 
| 1359 |  | for(int i=pfrom; i<pto; i++){ | 
| 1360 |  | for(int j=vfrom; j<vto; j++){ | 
| 1361 |  | nav++; | 
| 1362 |  | if(j==0 && XGood(i)) av += (xm[i]-xv[i])/resx[i]; | 
| 1363 |  | if(j==1 && YGood(i)) av += (ym[i]-yv[i])/resy[i]; | 
| 1364 |  | } | 
| 1365 |  | } | 
| 1366 |  | if(nav==0)return -100.; | 
| 1367 |  | return av/nav; | 
| 1368 |  |  | 
| 1369 |  | }; | 
| 1370 |  |  | 
| 1371 |  |  | 
| 1372 |  | /** | 
| 1373 |  | * \brief Give the maximum multiplicity on the x view | 
| 1374 |  | */ | 
| 1375 |  | Int_t TrkTrack::GetClusterX_Multiplicity_max(){ | 
| 1376 |  | int max=0; | 
| 1377 |  | for(int ip=0; ip<6; ip++) | 
| 1378 |  | if(GetClusterX_Multiplicity(ip)>max)max=GetClusterX_Multiplicity(ip); | 
| 1379 |  | return max; | 
| 1380 |  | }; | 
| 1381 |  | /** | 
| 1382 |  | * \brief Give the minimum multiplicity on the x view | 
| 1383 |  | */ | 
| 1384 |  | Int_t TrkTrack::GetClusterX_Multiplicity_min(){ | 
| 1385 |  | int min=50; | 
| 1386 |  | for(int ip=0; ip<6; ip++) | 
| 1387 |  | if(GetClusterX_Multiplicity(ip)<min)min=GetClusterX_Multiplicity(ip); | 
| 1388 |  | return min; | 
| 1389 |  | }; | 
| 1390 |  | /** | 
| 1391 |  | * \brief Give the maximum multiplicity on the x view | 
| 1392 |  | */ | 
| 1393 |  | Int_t TrkTrack::GetClusterY_Multiplicity_max(){ | 
| 1394 |  | int max=0; | 
| 1395 |  | for(int ip=0; ip<6; ip++) | 
| 1396 |  | if(GetClusterY_Multiplicity(ip)>max)max=GetClusterY_Multiplicity(ip); | 
| 1397 |  | return max; | 
| 1398 |  | }; | 
| 1399 |  | /** | 
| 1400 |  | * \brief Give the minimum multiplicity on the x view | 
| 1401 |  | */ | 
| 1402 |  | Int_t TrkTrack::GetClusterY_Multiplicity_min(){ | 
| 1403 |  | int min=50; | 
| 1404 |  | for(int ip=0; ip<6; ip++) | 
| 1405 |  | if(GetClusterY_Multiplicity(ip)<min)min=GetClusterY_Multiplicity(ip); | 
| 1406 |  | return min; | 
| 1407 |  | }; | 
| 1408 |  |  | 
| 1409 |  | /** | 
| 1410 |  | * \brief Give the minimum seed on the x view | 
| 1411 |  | */ | 
| 1412 |  | Float_t TrkTrack::GetClusterX_Seed_min(){ | 
| 1413 |  | Float_t min=100000; | 
| 1414 |  | for(int ip=0; ip<6; ip++) | 
| 1415 |  | if(XGood(ip) && GetClusterX_Seed(ip)<min)min=GetClusterX_Seed(ip); | 
| 1416 |  | return min; | 
| 1417 |  | }; | 
| 1418 |  | /** | 
| 1419 |  | * \brief Give the minimum seed on the x view | 
| 1420 |  | */ | 
| 1421 |  | Float_t TrkTrack::GetClusterY_Seed_min(){ | 
| 1422 |  | Float_t min=100000; | 
| 1423 |  | for(int ip=0; ip<6; ip++) | 
| 1424 |  | if(YGood(ip) && GetClusterY_Seed(ip)<min)min=GetClusterY_Seed(ip); | 
| 1425 |  | return min; | 
| 1426 |  | }; | 
| 1427 |  |  | 
| 1428 |  |  | 
| 1429 | //-------------------------------------- | //-------------------------------------- | 
| 1430 | // | // | 
| 1448 | zm[ip]     = 0; | zm[ip]     = 0; | 
| 1449 | resx[ip]   = 0; | resx[ip]   = 0; | 
| 1450 | resy[ip]   = 0; | resy[ip]   = 0; | 
| 1451 |  | tailx[ip]  = 0; | 
| 1452 |  | taily[ip]  = 0; | 
| 1453 | xv[ip]     = 0; | xv[ip]     = 0; | 
| 1454 | yv[ip]     = 0; | yv[ip]     = 0; | 
| 1455 | zv[ip]     = 0; | zv[ip]     = 0; | 
| 1459 | dedx_y[ip] = 0; | dedx_y[ip] = 0; | 
| 1460 |  |  | 
| 1461 | }; | }; | 
| 1462 | if(clx)clx->Clear(); | int ngf = TrkParams::nGF; | 
| 1463 | if(cly)cly->Clear(); | for(int i=0; i<ngf; i++){ | 
| 1464 |  | xGF[i] = 0.; | 
| 1465 |  | yGF[i] = 0.; | 
| 1466 |  | } | 
| 1467 |  | //     if(clx)clx->Clear(); | 
| 1468 |  | //     if(cly)cly->Clear(); | 
| 1469 |  | //    clx.Clear(); | 
| 1470 |  | //    cly.Clear(); | 
| 1471 | }; | }; | 
| 1472 | //-------------------------------------- | //-------------------------------------- | 
| 1473 | // | // | 
| 1475 | //-------------------------------------- | //-------------------------------------- | 
| 1476 | void TrkTrack::Delete(){ | void TrkTrack::Delete(){ | 
| 1477 | //    cout << "TrkTrack::Delete()"<<endl; | //    cout << "TrkTrack::Delete()"<<endl; | 
| 1478 | //    Clear(); | Clear(); | 
| 1479 | if(clx)delete clx; | //    if(clx)delete clx; | 
| 1480 | if(cly)delete cly; | //    if(cly)delete cly; | 
| 1481 | }; | }; | 
| 1482 | //-------------------------------------- | //-------------------------------------- | 
| 1483 | // | // | 
| 1490 | //-------------------------------------- | //-------------------------------------- | 
| 1491 | TrkSinglet::TrkSinglet(){ | TrkSinglet::TrkSinglet(){ | 
| 1492 | //    cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl; | //    cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl; | 
| 1493 | plane    = 0; | //     plane    = 0; | 
| 1494 | coord[0] = 0; | //     coord[0] = 0; | 
| 1495 | coord[1] = 0; | //     coord[1] = 0; | 
| 1496 | sgnl     = 0; | //     sgnl     = 0; | 
| 1497 |  | //     multmax  = 0; | 
| 1498 | //    cls      = 0; | //    cls      = 0; | 
| 1499 |  | Clear(); | 
| 1500 | }; | }; | 
| 1501 | //-------------------------------------- | //-------------------------------------- | 
| 1502 | // | // | 
| 1508 | coord[0] = s.coord[0]; | coord[0] = s.coord[0]; | 
| 1509 | coord[1] = s.coord[1]; | coord[1] = s.coord[1]; | 
| 1510 | sgnl     = s.sgnl; | sgnl     = s.sgnl; | 
| 1511 |  | multmax  = s.multmax; | 
| 1512 | //      cls      = 0;//<<<<pointer | //      cls      = 0;//<<<<pointer | 
| 1513 | cls      = TRef(s.cls); | //    cls      = TRef(s.cls); | 
| 1514 | }; | }; | 
| 1515 | //-------------------------------------- | //-------------------------------------- | 
| 1516 | // | // | 
| 1519 | void TrkSinglet::Dump(){ | void TrkSinglet::Dump(){ | 
| 1520 | int i=0; | int i=0; | 
| 1521 | cout << endl << "========== Singlet " ; | cout << endl << "========== Singlet " ; | 
| 1522 | cout << endl << "plane    : " << plane; | cout << endl << "plane        : " << plane; | 
| 1523 | cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++; | cout << endl << "coord[2]     : "; while( i<2 && cout << coord[i] << " ") i++; | 
| 1524 | cout << endl << "sgnl     : " << sgnl; | cout << endl << "sgnl         : " << sgnl; | 
| 1525 |  | cout << endl << "max.strip    : "; | 
| 1526 |  | cout << endl << "multiplicity : "; | 
| 1527 | } | } | 
| 1528 | //-------------------------------------- | //-------------------------------------- | 
| 1529 | // | // | 
| 1536 | coord[0]=-999; | coord[0]=-999; | 
| 1537 | coord[1]=-999; | coord[1]=-999; | 
| 1538 | sgnl=0; | sgnl=0; | 
| 1539 |  | multmax  = 0; | 
| 1540 |  |  | 
| 1541 | } | } | 
| 1542 | //-------------------------------------- | //-------------------------------------- | 
| 1544 | // | // | 
| 1545 | //-------------------------------------- | //-------------------------------------- | 
| 1546 | TrkLevel2::TrkLevel2(){ | TrkLevel2::TrkLevel2(){ | 
| 1547 | //    cout <<"TrkLevel2::TrkLevel2()"<<endl; | //    cout <<"TrkLevel2::TrkLevel2()"<<endl; | 
| 1548 | for(Int_t i=0; i<12 ; i++){ | for(Int_t i=0; i<12 ; i++){ | 
| 1549 | good[i] = -1; | good[i] = -1; | 
| 1550 | }; | VKmask[i] = 0; | 
| 1551 | //    Track    = new TClonesArray("TrkTrack"); | VKflag[i] = 0; | 
| 1552 | //    SingletX = new TClonesArray("TrkSinglet"); | }; | 
|  | //    SingletY = new TClonesArray("TrkSinglet"); |  | 
| 1553 | Track    = 0; | Track    = 0; | 
| 1554 | SingletX = 0; | SingletX = 0; | 
| 1555 | SingletY = 0; | SingletY = 0; | 
| 1559 | // | // | 
| 1560 | // | // | 
| 1561 | //-------------------------------------- | //-------------------------------------- | 
| 1562 |  | void TrkLevel2::Set(){ | 
| 1563 |  | if(!Track)Track    = new TClonesArray("TrkTrack"); | 
| 1564 |  | if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); | 
| 1565 |  | if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); | 
| 1566 |  | } | 
| 1567 |  | //-------------------------------------- | 
| 1568 |  | // | 
| 1569 |  | // | 
| 1570 |  | //-------------------------------------- | 
| 1571 | void TrkLevel2::Dump(){ | void TrkLevel2::Dump(){ | 
| 1572 |  |  | 
| 1573 | // | // | 
| 1574 | cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; | cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; | 
| 1575 | cout << endl << "good     : "; for(int i=0; i<12; i++) cout << good[i]<<" "; | cout << endl << "good     : "; for(int i=0; i<12; i++) cout << hex <<" 0x"<< good[i]<<dec; | 
| 1576 | cout << endl << "ntrk()   : " << this->ntrk() ; | cout << endl << "ntrk()   : " << ntrk() ; | 
| 1577 | cout << endl << "nclsx()  : " << this->nclsx(); | cout << endl << "nclsx()  : " << nclsx(); | 
| 1578 | cout << endl << "nclsy()  : " << this->nclsy(); | cout << endl << "nclsy()  : " << nclsy(); | 
|  | //     TClonesArray &t  = *Track; |  | 
|  | //     TClonesArray &sx = *SingletX; |  | 
|  | //     TClonesArray &sy = *SingletY; |  | 
|  | //     for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump(); |  | 
|  | //     for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |  | 
|  | //     for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |  | 
| 1579 | if(Track){ | if(Track){ | 
| 1580 | TClonesArray &t  = *Track; | TClonesArray &t  = *Track; | 
| 1581 | for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump(); | for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump(); | 
| 1582 | } | } | 
| 1583 | if(SingletX){ | //     if(SingletX){ | 
| 1584 | TClonesArray &sx = *SingletX; | //      TClonesArray &sx = *SingletX; | 
| 1585 | for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); | //      for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); | 
| 1586 | } | //     } | 
| 1587 | if(SingletY){ | //     if(SingletY){ | 
| 1588 | TClonesArray &sy = *SingletY; | //      TClonesArray &sy = *SingletY; | 
| 1589 | for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); | //      for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); | 
| 1590 | } | //     } | 
| 1591 |  | cout << endl; | 
| 1592 | } | } | 
| 1593 |  | /** | 
| 1594 |  | * \brief Dump processing status | 
| 1595 |  | */ | 
| 1596 |  | void TrkLevel2::StatusDump(int view){ | 
| 1597 |  | cout << "DSP n. "<<view+1<<" status: "<<hex<<good[view]<<endl; | 
| 1598 |  | }; | 
| 1599 |  | /** | 
| 1600 |  | * \brief Check event status | 
| 1601 |  | * | 
| 1602 |  | * Check the event status, according to a flag-mask given as input. | 
| 1603 |  | * Return true if the view passes the check. | 
| 1604 |  | * | 
| 1605 |  | * @param view View number (0-11) | 
| 1606 |  | * @param flagmask Mask of flags to check (eg. flagmask=0x111 no missing packet, | 
| 1607 |  | *  no crc error, no software alarm) | 
| 1608 |  | * | 
| 1609 |  | * @see TrkLevel2 class definition to know how the status flag is defined | 
| 1610 |  | * | 
| 1611 |  | */ | 
| 1612 |  | Bool_t TrkLevel2::StatusCheck(int view, int flagmask){ | 
| 1613 |  |  | 
| 1614 |  | if( view<0 || view >= 12)return false; | 
| 1615 |  | return !(good[view]&flagmask); | 
| 1616 |  |  | 
| 1617 |  | }; | 
| 1618 |  |  | 
| 1619 |  |  | 
| 1620 | //-------------------------------------- | //-------------------------------------- | 
| 1621 | // | // | 
| 1622 | // | // | 
| 1623 | //-------------------------------------- | //-------------------------------------- | 
| 1624 | /** | /** | 
| 1625 | * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). | * The method returns false if the viking-chip was masked | 
| 1626 |  | * either apriori ,on the basis of the mask read from the DB, | 
| 1627 |  | * or run-by-run, on the basis of the calibration parameters) | 
| 1628 |  | * @param iv Tracker view (0-11) | 
| 1629 |  | * @param ivk Viking-chip number (0-23) | 
| 1630 |  | */ | 
| 1631 |  | Bool_t TrkLevel2::GetVKMask(int iv, int ivk){ | 
| 1632 |  | Int_t whichbit = (Int_t)pow(2,ivk); | 
| 1633 |  | return (whichbit&VKmask[iv])!=0; | 
| 1634 |  | } | 
| 1635 |  | /** | 
| 1636 |  | * The method returns false if the viking-chip was masked | 
| 1637 |  | * for this event due to common-noise computation failure. | 
| 1638 |  | * @param iv Tracker view (0-11) | 
| 1639 |  | * @param ivk Viking-chip number (0-23) | 
| 1640 |  | */ | 
| 1641 |  | Bool_t TrkLevel2::GetVKFlag(int iv, int ivk){ | 
| 1642 |  | Int_t whichbit = (Int_t)pow(2,ivk); | 
| 1643 |  | return (whichbit&VKflag[iv])!=0; | 
| 1644 |  | } | 
| 1645 |  | /** | 
| 1646 |  | * The method returns true if the viking-chip was masked, either | 
| 1647 |  | * forced (see TrkLevel2::GetVKMask(int,int)) or | 
| 1648 |  | * for this event only (TrkLevel2::GetVKFlag(int,int)). | 
| 1649 |  | * @param iv Tracker view (0-11) | 
| 1650 |  | * @param ivk Viking-chip number (0-23) | 
| 1651 |  | */ | 
| 1652 |  | Bool_t TrkLevel2::IsMaskedVK(int iv, int ivk){ | 
| 1653 |  | return !(GetVKMask(iv,ivk)&&GetVKFlag(iv,ivk) ); | 
| 1654 |  | }; | 
| 1655 |  |  | 
| 1656 |  | //-------------------------------------- | 
| 1657 |  | // | 
| 1658 |  | // | 
| 1659 |  | //-------------------------------------- | 
| 1660 |  | /** | 
| 1661 |  | * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). | 
| 1662 |  | * Ref to Level1 data (clusters) is also set. If l1==NULL no references are set. | 
| 1663 |  | * (NB It make sense to set references only if events are stored in a tree that contains also the Level1 branch) | 
| 1664 | */ | */ | 
| 1665 | void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){ | void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){ | 
| 1666 |  |  | 
| 1667 |  | //    cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<<endl; | 
| 1668 |  | Clear(); | 
| 1669 |  |  | 
| 1670 | //  temporary objects: | //  temporary objects: | 
| 1671 | TrkSinglet* t_singlet = new TrkSinglet(); | TrkSinglet* t_singlet = new TrkSinglet(); | 
| 1672 | TrkTrack*   t_track   = new TrkTrack(); | TrkTrack*   t_track   = new TrkTrack(); | 
| 1673 |  |  | 
| 1674 | //  **** general variables **** | //  ----------------- | 
| 1675 | //    good2 = l2->good2; | //  general variables | 
| 1676 |  | //  ----------------- | 
| 1677 | for(Int_t i=0; i<12 ; i++){ | for(Int_t i=0; i<12 ; i++){ | 
| 1678 | //              crc[i] = l2->crc[i]; | good[i] = l2->good[i]; | 
| 1679 | good[i] = l2->good[i]; | VKmask[i]=0; | 
| 1680 |  | VKflag[i]=0; | 
| 1681 |  | for(Int_t ii=0; ii<24 ; ii++){ | 
| 1682 |  | Int_t setbit = (Int_t)pow(2,ii); | 
| 1683 |  | if( l2->vkflag[ii][i]!=-1 )VKmask[i]=VKmask[i]|setbit; | 
| 1684 |  | if( l2->vkflag[ii][i]!=0  )VKflag[i]=VKflag[i]|setbit; | 
| 1685 | }; | }; | 
| 1686 | //  *** TRACKS *** | }; | 
| 1687 |  | //  -------------- | 
| 1688 |  | //  *** TRACKS *** | 
| 1689 |  | //  -------------- | 
| 1690 | if(!Track) Track = new TClonesArray("TrkTrack"); | if(!Track) Track = new TClonesArray("TrkTrack"); | 
| 1691 | TClonesArray &t = *Track; | TClonesArray &t = *Track; | 
| 1692 |  |  | 
| 1693 | for(int i=0; i<l2->ntrk; i++){ | for(int i=0; i<l2->ntrk; i++){ | 
| 1694 | t_track->seqno = i;// NBNBNBNB deve sempre essere = i | t_track->seqno = i;// NBNBNBNB deve sempre essere = i | 
| 1695 | t_track->image = l2->image[i]-1; | t_track->image = l2->image[i]-1; | 
|  | //      cout << "track "<<i<<t_track->seqno << t_track->image<<endl; |  | 
| 1696 | t_track->chi2  = l2->chi2_nt[i]; | t_track->chi2  = l2->chi2_nt[i]; | 
| 1697 | t_track->nstep = l2->nstep_nt[i]; | t_track->nstep = l2->nstep_nt[i]; | 
| 1698 | for(int it1=0;it1<5;it1++){ | for(int it1=0;it1<5;it1++){ | 
| 1701 | t_track->coval[it1][it2] = l2->coval[i][it2][it1]; | t_track->coval[it1][it2] = l2->coval[i][it2][it1]; | 
| 1702 | }; | }; | 
| 1703 | for(int ip=0;ip<6;ip++){ | for(int ip=0;ip<6;ip++){ | 
| 1704 | t_track->xgood[ip]  = l2->xgood_nt[i][ip]; | // --------------------------------- | 
| 1705 | t_track->ygood[ip]  = l2->ygood_nt[i][ip]; | // new implementation of xgood/ygood | 
| 1706 |  | // --------------------------------- | 
| 1707 |  | t_track->xgood[ip]  = l2->cltrx[i][ip]; //cluster ID | 
| 1708 |  | t_track->ygood[ip]  = l2->cltry[i][ip]; //cluster ID | 
| 1709 |  | t_track->xgood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor | 
| 1710 |  | t_track->ygood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor | 
| 1711 |  | if(l2->xbad[i][ip]>0)t_track->xgood[ip]=-t_track->xgood[ip]; | 
| 1712 |  | if(l2->ybad[i][ip]>0)t_track->ygood[ip]=-t_track->ygood[ip]; | 
| 1713 |  | //          if(l2->xbad[i][ip]>0 || l2->ybad[i][ip]>0){ | 
| 1714 |  | //          if(l2->dedx_x[i][ip]<0 || l2->dedx_y[i][ip]<0){ | 
| 1715 |  | //              cout << ip << " - "<< l2->cltrx[i][ip] << " "<<l2->cltry[i][ip]<<" "<<l2->ls[i][ip]<<endl; | 
| 1716 |  | //              cout << ip << " - "<<t_track->xgood[ip]<<" "<<t_track->ygood[ip]<<endl; | 
| 1717 |  | //              cout << ip << " - "<<t_track->GetClusterX_ID(ip)<<" "<<t_track->GetClusterY_ID(ip)<<" "<<t_track->GetLadder(ip)<<" "<<t_track->GetSensor(ip)<<endl; | 
| 1718 |  | //              cout << ip << " - "<<t_track->BadClusterX(ip)<<" "<<t_track->BadClusterY(ip)<<endl; | 
| 1719 |  | //              cout << ip << " - "<<t_track->SaturatedClusterX(ip)<<" "<<t_track->SaturatedClusterY(ip)<<endl; | 
| 1720 |  | //          } | 
| 1721 | t_track->xm[ip]     = l2->xm_nt[i][ip]; | t_track->xm[ip]     = l2->xm_nt[i][ip]; | 
| 1722 | t_track->ym[ip]     = l2->ym_nt[i][ip]; | t_track->ym[ip]     = l2->ym_nt[i][ip]; | 
| 1723 | t_track->zm[ip]     = l2->zm_nt[i][ip]; | t_track->zm[ip]     = l2->zm_nt[i][ip]; | 
| 1724 | t_track->resx[ip]   = l2->resx_nt[i][ip]; | t_track->resx[ip]   = l2->resx_nt[i][ip]; | 
| 1725 | t_track->resy[ip]   = l2->resy_nt[i][ip]; | t_track->resy[ip]   = l2->resy_nt[i][ip]; | 
| 1726 |  | t_track->tailx[ip]  = l2->tailx[i][ip]; | 
| 1727 |  | t_track->taily[ip]  = l2->taily[i][ip]; | 
| 1728 | t_track->xv[ip]     = l2->xv_nt[i][ip]; | t_track->xv[ip]     = l2->xv_nt[i][ip]; | 
| 1729 | t_track->yv[ip]     = l2->yv_nt[i][ip]; | t_track->yv[ip]     = l2->yv_nt[i][ip]; | 
| 1730 | t_track->zv[ip]     = l2->zv_nt[i][ip]; | t_track->zv[ip]     = l2->zv_nt[i][ip]; | 
| 1732 | t_track->ayv[ip]    = l2->ayv_nt[i][ip]; | t_track->ayv[ip]    = l2->ayv_nt[i][ip]; | 
| 1733 | t_track->dedx_x[ip] = l2->dedx_x[i][ip]; | t_track->dedx_x[ip] = l2->dedx_x[i][ip]; | 
| 1734 | t_track->dedx_y[ip] = l2->dedx_y[i][ip]; | t_track->dedx_y[ip] = l2->dedx_y[i][ip]; | 
| 1735 | //                      t_track->clx[ip] = 0; | t_track->multmaxx[ip] = l2->multmaxx[i][ip]; | 
| 1736 | //                      t_track->cly[ip] = 0; | t_track->multmaxy[ip] = l2->multmaxy[i][ip]; | 
| 1737 |  | t_track->seedx[ip]  = l2->seedx[i][ip]; | 
| 1738 |  | t_track->seedy[ip]  = l2->seedy[i][ip]; | 
| 1739 |  | t_track->xpu[ip]    = l2->xpu[i][ip]; | 
| 1740 |  | t_track->ypu[ip]    = l2->ypu[i][ip]; | 
| 1741 |  | //----------------------------------------------------- | 
| 1742 |  | //----------------------------------------------------- | 
| 1743 |  | //----------------------------------------------------- | 
| 1744 |  | //----------------------------------------------------- | 
| 1745 | }; | }; | 
| 1746 |  | // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | 
| 1747 |  | // evaluated coordinates (to define GF) | 
| 1748 |  | // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | 
| 1749 |  | int    ngf = TrkParams::nGF; | 
| 1750 |  | float *zgf = TrkParams::zGF; | 
| 1751 |  | Trajectory tgf = Trajectory(ngf,zgf); | 
| 1752 |  | tgf.DoTrack(t_track->al);//<<<< integrate the trajectory | 
| 1753 |  | for(int ip=0; ip<ngf; ip++){ | 
| 1754 |  | t_track->xGF[ip] = tgf.x[ip]; | 
| 1755 |  | t_track->yGF[ip] = tgf.y[ip]; | 
| 1756 |  | } | 
| 1757 |  |  | 
| 1758 |  | //      if(t_track->IsSaturated())t_track->Dump(); | 
| 1759 | new(t[i]) TrkTrack(*t_track); | new(t[i]) TrkTrack(*t_track); | 
| 1760 | t_track->Clear(); | t_track->Clear(); | 
| 1761 | }; | };//end loop over track | 
| 1762 |  |  | 
| 1763 |  | //  ---------------- | 
| 1764 | //  *** SINGLETS *** | //  *** SINGLETS *** | 
| 1765 |  | //  ---------------- | 
| 1766 | if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); | if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); | 
| 1767 | TClonesArray &sx = *SingletX; | TClonesArray &sx = *SingletX; | 
| 1768 | for(int i=0; i<l2->nclsx; i++){ | for(int i=0; i<l2->nclsx; i++){ | 
| 1770 | t_singlet->coord[0] = l2->xs[i][0]; | t_singlet->coord[0] = l2->xs[i][0]; | 
| 1771 | t_singlet->coord[1] = l2->xs[i][1]; | t_singlet->coord[1] = l2->xs[i][1]; | 
| 1772 | t_singlet->sgnl     = l2->signlxs[i]; | t_singlet->sgnl     = l2->signlxs[i]; | 
| 1773 | //              t_singlet->cls      = 0; | t_singlet->multmax = l2->multmaxsx[i]; | 
| 1774 |  | if(l2->sxbad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; | 
| 1775 |  | //----------------------------------------------------- | 
| 1776 |  | //      if(l1) t_singlet->cls      = l1->GetCluster(l2->clsx[i]-1); | 
| 1777 |  | //----------------------------------------------------- | 
| 1778 | new(sx[i]) TrkSinglet(*t_singlet); | new(sx[i]) TrkSinglet(*t_singlet); | 
| 1779 | t_singlet->Clear(); | t_singlet->Clear(); | 
| 1780 | } | } | 
| 1785 | t_singlet->coord[0] = l2->ys[i][0]; | t_singlet->coord[0] = l2->ys[i][0]; | 
| 1786 | t_singlet->coord[1] = l2->ys[i][1]; | t_singlet->coord[1] = l2->ys[i][1]; | 
| 1787 | t_singlet->sgnl     = l2->signlys[i]; | t_singlet->sgnl     = l2->signlys[i]; | 
| 1788 | //              t_singlet->cls      = 0; | t_singlet->multmax  = l2->multmaxsy[i]; | 
| 1789 |  | if(l2->sybad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; | 
| 1790 |  | //----------------------------------------------------- | 
| 1791 |  | //      if(l1) t_singlet->cls      = l1->GetCluster(l2->clsy[i]-1); | 
| 1792 |  | //----------------------------------------------------- | 
| 1793 | new(sy[i]) TrkSinglet(*t_singlet); | new(sy[i]) TrkSinglet(*t_singlet); | 
| 1794 | t_singlet->Clear(); | t_singlet->Clear(); | 
| 1795 | }; | }; | 
| 1796 |  |  | 
| 1797 |  |  | 
| 1798 |  |  | 
| 1799 | delete t_track; | delete t_track; | 
| 1800 | delete t_singlet; | delete t_singlet; | 
| 1801 | } | } | 
|  | //-------------------------------------- |  | 
|  | // |  | 
|  | // |  | 
|  | //-------------------------------------- |  | 
|  | /** |  | 
|  | * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |  | 
|  | * Ref to Level1 data (clusters) is also set. |  | 
|  | */ |  | 
|  | void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){ |  | 
|  |  |  | 
|  | //  temporary objects: |  | 
|  | TrkSinglet* t_singlet = new TrkSinglet(); |  | 
|  | TrkTrack*   t_track   = new TrkTrack(); |  | 
|  | // general variables |  | 
|  | //      good2 = l2->good2; |  | 
|  | for(Int_t i=0; i<12 ; i++){ |  | 
|  | //              crc[i] = l2->crc[i]; |  | 
|  | good[i] = l2->good[i]; |  | 
|  | }; |  | 
|  | // *** TRACKS *** |  | 
|  | if(!Track) Track = new TClonesArray("TrkTrack"); |  | 
|  | TClonesArray &t = *Track; |  | 
|  | if(!t_track->clx)t_track->clx = new TRefArray(6,0); |  | 
|  | if(!t_track->cly)t_track->cly = new TRefArray(6,0); |  | 
|  | for(int i=0; i<l2->ntrk; i++){ |  | 
|  | t_track->seqno = i;// NBNBNBNB deve sempre essere = i |  | 
|  | t_track->image = l2->image[i]-1; |  | 
|  | //              cout << "track "<<i<<t_track->seqno << t_track->image<<endl; |  | 
|  | t_track->chi2  = l2->chi2_nt[i]; |  | 
|  | t_track->nstep = l2->nstep_nt[i]; |  | 
|  | for(int it1=0;it1<5;it1++){ |  | 
|  | t_track->al[it1] = l2->al_nt[i][it1]; |  | 
|  | for(int it2=0;it2<5;it2++) |  | 
|  | t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |  | 
|  | }; |  | 
|  | for(int ip=0;ip<6;ip++){ |  | 
|  | t_track->xgood[ip]  = l2->xgood_nt[i][ip]; |  | 
|  | t_track->ygood[ip]  = l2->ygood_nt[i][ip]; |  | 
|  | t_track->xm[ip]     = l2->xm_nt[i][ip]; |  | 
|  | t_track->ym[ip]     = l2->ym_nt[i][ip]; |  | 
|  | t_track->zm[ip]     = l2->zm_nt[i][ip]; |  | 
|  | t_track->resx[ip]   = l2->resx_nt[i][ip]; |  | 
|  | t_track->resy[ip]   = l2->resy_nt[i][ip]; |  | 
|  | t_track->xv[ip]     = l2->xv_nt[i][ip]; |  | 
|  | t_track->yv[ip]     = l2->yv_nt[i][ip]; |  | 
|  | t_track->zv[ip]     = l2->zv_nt[i][ip]; |  | 
|  | t_track->axv[ip]    = l2->axv_nt[i][ip]; |  | 
|  | t_track->ayv[ip]    = l2->ayv_nt[i][ip]; |  | 
|  | t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |  | 
|  | t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |  | 
|  | //                      cout << "traccia "<<i<<"  --  "<< ip << " "<< l2->cltrx[i][ip] <<" "<< l2->cltry[i][ip] <<" "<< t_track->xgood[ip] << t_track->ygood[ip]<<endl; |  | 
|  | //----------------------------------------------------- |  | 
|  | //----------------------------------------------------- |  | 
|  | if(t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip); |  | 
|  | if(t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip); |  | 
|  | //                      if(t_track->ygood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltry[i][ip] "<<l2->cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<<l1->GetCluster(l2->cltry[i][ip]-1)<<endl; |  | 
|  | //                      if(t_track->xgood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltrx[i][ip] "<<l2->cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<<l1->GetCluster(l2->cltrx[i][ip]-1)<<endl; |  | 
|  | //----------------------------------------------------- |  | 
|  | //----------------------------------------------------- |  | 
|  | }; |  | 
|  | new(t[i]) TrkTrack(*t_track); |  | 
|  | t_track->Clear(); |  | 
|  | }; |  | 
|  | // *** SINGLETS *** |  | 
|  | if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |  | 
|  | TClonesArray &sx = *SingletX; |  | 
|  | for(int i=0; i<l2->nclsx; i++){ |  | 
|  | t_singlet->plane    = l2->planex[i]; |  | 
|  | t_singlet->coord[0] = l2->xs[i][0]; |  | 
|  | t_singlet->coord[1] = l2->xs[i][1]; |  | 
|  | t_singlet->sgnl     = l2->signlxs[i]; |  | 
|  | //----------------------------------------------------- |  | 
|  | //              cout << "singolo x "<<i<<"  --  "<<  l2->clsx[i] <<endl; |  | 
|  | t_singlet->cls      = l1->GetCluster(l2->clsx[i]-1); |  | 
|  | //              cout<<" i "<<i<<" l2->clsx[i] "<<l2->clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<<l1->GetCluster(l2->clsx[i]-1)<<endl; |  | 
|  | //----------------------------------------------------- |  | 
|  | new(sx[i]) TrkSinglet(*t_singlet); |  | 
|  | t_singlet->Clear(); |  | 
|  | } |  | 
|  | if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |  | 
|  | TClonesArray &sy = *SingletY; |  | 
|  | for(int i=0; i<l2->nclsy; i++){ |  | 
|  | t_singlet->plane    = l2->planey[i]; |  | 
|  | t_singlet->coord[0] = l2->ys[i][0]; |  | 
|  | t_singlet->coord[1] = l2->ys[i][1]; |  | 
|  | t_singlet->sgnl     = l2->signlys[i]; |  | 
|  | //----------------------------------------------------- |  | 
|  | //              cout << "singolo y "<<i<<"  --  "<<  l2->clsy[i] <<endl; |  | 
|  | t_singlet->cls      = l1->GetCluster(l2->clsy[i]-1); |  | 
|  | //              cout<<" i "<<i<<" l2->clsy[i] "<<l2->clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<<l1->GetCluster(l2->clsy[i]-1)<<endl; |  | 
|  | //----------------------------------------------------- |  | 
|  | new(sy[i]) TrkSinglet(*t_singlet); |  | 
|  | t_singlet->Clear(); |  | 
|  | }; |  | 
|  |  |  | 
|  | delete t_track; |  | 
|  | delete t_singlet; |  | 
|  | } |  | 
| 1802 | /** | /** | 
| 1803 | * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common). | * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common). | 
| 1804 | */ | */ | 
| 1825 | l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2]; | l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2]; | 
| 1826 | }; | }; | 
| 1827 | for(int ip=0;ip<6;ip++){ | for(int ip=0;ip<6;ip++){ | 
| 1828 | l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->xgood[ip]; | l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->XGood(ip); | 
| 1829 | l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->ygood[ip]; | l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->YGood(ip); | 
| 1830 | l2->xm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xm[ip]; | l2->xm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xm[ip]; | 
| 1831 | l2->ym_nt[i][ip]    = ((TrkTrack *)Track->At(i))->ym[ip]; | l2->ym_nt[i][ip]    = ((TrkTrack *)Track->At(i))->ym[ip]; | 
| 1832 | l2->zm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zm[ip]; | l2->zm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zm[ip]; | 
| 1833 | l2->resx_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resx[ip]; | l2->resx_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resx[ip]; | 
| 1834 | l2->resy_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resy[ip]; | l2->resy_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resy[ip]; | 
| 1835 |  | l2->tailx[i][ip]  = ((TrkTrack *)Track->At(i))->tailx[ip]; | 
| 1836 |  | l2->taily[i][ip]  = ((TrkTrack *)Track->At(i))->taily[ip]; | 
| 1837 | l2->xv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xv[ip]; | l2->xv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xv[ip]; | 
| 1838 | l2->yv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->yv[ip]; | l2->yv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->yv[ip]; | 
| 1839 | l2->zv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zv[ip]; | l2->zv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zv[ip]; | 
| 1872 | void TrkLevel2::Clear(){ | void TrkLevel2::Clear(){ | 
| 1873 | for(Int_t i=0; i<12 ; i++){ | for(Int_t i=0; i<12 ; i++){ | 
| 1874 | good[i] = -1; | good[i] = -1; | 
| 1875 |  | VKflag[i] = 0; | 
| 1876 |  | VKmask[i] = 0; | 
| 1877 | }; | }; | 
| 1878 | //    if(Track)Track->Clear("C"); | //    if(Track)Track->Clear("C"); | 
| 1879 | //    if(SingletX)SingletX->Clear("C"); | //    if(SingletX)SingletX->Clear("C"); | 
| 1907 |  |  | 
| 1908 | if(!Track)return 0; | if(!Track)return 0; | 
| 1909 |  |  | 
| 1910 | TRefArray *sorted = new TRefArray(); | //    TRefArray *sorted = new TRefArray(); | 
| 1911 |  | TRefArray *sorted = NULL; | 
| 1912 |  |  | 
| 1913 | TClonesArray &t  = *Track; | TClonesArray &t  = *Track; | 
| 1914 | //    TClonesArray &ts = *PhysicalTrack; | //    TClonesArray &ts = *PhysicalTrack; | 
| 1918 |  |  | 
| 1919 | int indo=0; | int indo=0; | 
| 1920 | int indi=0; | int indi=0; | 
| 1921 | while(N != 0){ | while(N > 0){ | 
| 1922 |  | //    while(N != 0){ | 
| 1923 | int nfit =0; | int nfit =0; | 
| 1924 | float chi2ref = numeric_limits<float>::max(); | float chi2ref = numeric_limits<float>::max(); | 
| 1925 |  |  | 
| 1930 | } | } | 
| 1931 | } | } | 
| 1932 | //second loop to search minimum chi2 among selected | //second loop to search minimum chi2 among selected | 
| 1933 | for(int i=0; i<this->ntrk(); i++){ | for(int i=0; i<ntrk(); i++){ | 
| 1934 | Float_t chi2 = ((TrkTrack *)t[i])->chi2; | Float_t chi2 = ((TrkTrack *)t[i])->chi2; | 
| 1935 | if(chi2 < 0) chi2 = chi2*1000; | if(chi2 < 0) chi2 = -chi2*1000; | 
| 1936 | if(    chi2 < chi2ref | if(    chi2 < chi2ref | 
| 1937 | && ((TrkTrack *)t[i])->GetNtot() == nfit | && ((TrkTrack *)t[i])->GetNtot() == nfit | 
| 1938 | && m[i]==1){ | && m[i]==1){ | 
| 1944 | m[((TrkTrack *)t[indi])->image] = 0; | m[((TrkTrack *)t[indi])->image] = 0; | 
| 1945 | N--; | N--; | 
| 1946 |  |  | 
| 1947 | //      cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; | //          cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; | 
| 1948 | }; | }; | 
| 1949 |  | if(!sorted)sorted = new TRefArray( TProcessID::GetProcessWithUID(t[indi])); | 
| 1950 | sorted->Add( (TrkTrack*)t[indi] ); | sorted->Add( (TrkTrack*)t[indi] ); | 
| 1951 |  |  | 
| 1952 | m[indi] = 0; | m[indi] = 0; | 
| 1953 | //              cout << "SORTED "<< indo << " "<< indi << " "<< N << endl; | //      cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "<<chi2ref<<endl; | 
| 1954 | N--; | N--; | 
| 1955 | indo++; | indo++; | 
| 1956 | } | } | 
| 1957 | m.clear(); | m.clear(); | 
| 1958 | //      cout << "GetTracks_NFitSorted(it): Done"<< endl; | //    cout << "GetTracks_NFitSorted(it): Done"<< endl; | 
| 1959 |  |  | 
| 1960 | return sorted; | return sorted; | 
| 1961 | //    return PhysicalTrack; | //    return PhysicalTrack; | 
| 1972 | TrkTrack *TrkLevel2::GetStoredTrack(int is){ | TrkTrack *TrkLevel2::GetStoredTrack(int is){ | 
| 1973 |  |  | 
| 1974 | if(is >= this->ntrk()){ | if(is >= this->ntrk()){ | 
| 1975 | cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl; | cout << "TrkTrack *TrkLevel2::GetStoredTrack(int) >> Track "<< is << "doen not exits! " << endl; | 
| 1976 | cout << "                Stored tracks ntrk() = "<< this->ntrk() << endl; | cout << "Stored tracks ntrk() = "<< this->ntrk() << endl; | 
| 1977 | return 0; | return 0; | 
| 1978 | } | } | 
| 1979 | if(!Track){ | if(!Track){ | 
| 1994 | TrkSinglet *TrkLevel2::GetSingletX(int is){ | TrkSinglet *TrkLevel2::GetSingletX(int is){ | 
| 1995 |  |  | 
| 1996 | if(is >= this->nclsx()){ | if(is >= this->nclsx()){ | 
| 1997 | cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl; | cout << "TrkSinglet *TrkLevel2::GetSingletX(int) >> Singlet "<< is << "doen not exits! " << endl; | 
| 1998 | cout << "                Stored x-singlets nclsx() = "<< this->nclsx() << endl; | cout << "Stored x-singlets nclsx() = "<< this->nclsx() << endl; | 
| 1999 | return 0; | return 0; | 
| 2000 | } | } | 
| 2001 | if(!SingletX)return 0; | if(!SingletX)return 0; | 
| 2014 | TrkSinglet *TrkLevel2::GetSingletY(int is){ | TrkSinglet *TrkLevel2::GetSingletY(int is){ | 
| 2015 |  |  | 
| 2016 | if(is >= this->nclsy()){ | if(is >= this->nclsy()){ | 
| 2017 | cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl; | cout << "TrkSinglet *TrkLevel2::GetSingletY(int) >> Singlet "<< is << "doen not exits! " << endl; | 
| 2018 | cout << "                Stored y-singlets nclsy() = "<< this->nclsx() << endl; | cout << "Stored y-singlets nclsx() = "<< this->nclsx() << endl; | 
| 2019 | return 0; | return 0; | 
| 2020 | } | } | 
| 2021 | if(!SingletY)return 0; | if(!SingletY)return 0; | 
| 2035 | TrkTrack *TrkLevel2::GetTrack(int it){ | TrkTrack *TrkLevel2::GetTrack(int it){ | 
| 2036 |  |  | 
| 2037 | if(it >= this->GetNTracks()){ | if(it >= this->GetNTracks()){ | 
| 2038 | cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; | cout << "TrkTrack *TrkLevel2::GetTrack(int) >> Track "<< it << "does not exits! " << endl; | 
| 2039 | cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl; | cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; | 
| 2040 | return 0; | return 0; | 
| 2041 | } | } | 
| 2042 |  |  | 
| 2073 | TrkTrack *TrkLevel2::GetTrackImage(int it){ | TrkTrack *TrkLevel2::GetTrackImage(int it){ | 
| 2074 |  |  | 
| 2075 | if(it >= this->GetNTracks()){ | if(it >= this->GetNTracks()){ | 
| 2076 | cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; | cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not exits! " << endl; | 
| 2077 | cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl; | cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; | 
| 2078 | return 0; | return 0; | 
| 2079 | } | } | 
| 2080 |  |  | 
| 2083 | TrkTrack *track = (TrkTrack*)sorted->At(it); | TrkTrack *track = (TrkTrack*)sorted->At(it); | 
| 2084 |  |  | 
| 2085 | if(!track->HasImage()){ | if(!track->HasImage()){ | 
| 2086 | cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; | cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not have image! " << endl; | 
| 2087 | return 0; | return 0; | 
| 2088 | } | } | 
| 2089 | if(!Track)return 0; | if(!Track)return 0; | 
| 2105 | */ | */ | 
| 2106 | void TrkLevel2::LoadField(TString path){ | void TrkLevel2::LoadField(TString path){ | 
| 2107 | // | // | 
| 2108 | strcpy(path_.path,path.Data()); | //     strcpy(path_.path,path.Data()); | 
| 2109 | path_.pathlen = path.Length(); | //     path_.pathlen = path.Length(); | 
| 2110 | path_.error   = 0; | //     path_.error   = 0; | 
| 2111 | readb_(); | //     readb_(); | 
| 2112 |  |  | 
| 2113 |  | //     TrkParams::SetTrackingMode(); | 
| 2114 |  | //     TrkParams::SetPrecisionFactor(); | 
| 2115 |  | //     TrkParams::SetStepMin(); | 
| 2116 |  | TrkParams::SetMiniDefault(); | 
| 2117 |  |  | 
| 2118 |  | TrkParams::Set(path,1); | 
| 2119 |  | TrkParams::Load(1); | 
| 2120 |  | if( !TrkParams::IsLoaded(1) ){ | 
| 2121 |  | cout << "void TrkLevel2::LoadField(TString path) --- ERROR --- m.field not loaded"<<endl; | 
| 2122 |  | } | 
| 2123 |  |  | 
| 2124 | // | // | 
| 2125 | }; | }; | 
| 2126 |  | // /** | 
| 2127 |  | //  * Get BY (kGauss) | 
| 2128 |  | //  * @param v (x,y,z) coordinates in cm | 
| 2129 |  | //  */ | 
| 2130 |  | // float TrkLevel2::GetBX(float* v){ | 
| 2131 |  | //     float b[3]; | 
| 2132 |  | //     gufld_(v,b); | 
| 2133 |  | //     return b[0]/10.; | 
| 2134 |  | // } | 
| 2135 |  | // /** | 
| 2136 |  | //  * Get BY (kGauss) | 
| 2137 |  | //  * @param v (x,y,z) coordinates in cm | 
| 2138 |  | //  */ | 
| 2139 |  | // float TrkLevel2::GetBY(float* v){ | 
| 2140 |  | //     float b[3]; | 
| 2141 |  | //     gufld_(v,b); | 
| 2142 |  | //     return b[1]/10.; | 
| 2143 |  | // } | 
| 2144 |  | // /** | 
| 2145 |  | //  * Get BY (kGauss) | 
| 2146 |  | //  * @param v (x,y,z) coordinates in cm | 
| 2147 |  | //  */ | 
| 2148 |  | // float TrkLevel2::GetBZ(float* v){ | 
| 2149 |  | //     float b[3]; | 
| 2150 |  | //     gufld_(v,b); | 
| 2151 |  | //     return b[2]/10.; | 
| 2152 |  | // } | 
| 2153 | //-------------------------------------- | //-------------------------------------- | 
| 2154 | // | // | 
| 2155 | // | // | 
| 2307 |  |  | 
| 2308 | /** | /** | 
| 2309 | * Evaluates the trajectory in the apparatus associated to the track. | * Evaluates the trajectory in the apparatus associated to the track. | 
| 2310 | * 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. | 
| 2311 | * @param t pointer to an object of the class Trajectory, | * @param al Track state-vector (X0,Y0,sin(theta),phi,deflection). | 
| 2312 | * 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). | 
| 2313 | * @return error flag. | * @return error flag. | 
| 2314 | */ | * | 
| 2315 | int Trajectory::DoTrack2(float* al){ | * This method is needed when you want to integrate the particle trajectory | 
| 2316 |  | * starting from a track state-vector relative to an arbitrary reference plane. | 
| 2317 | double *dxout   = new double[npoint]; | * The default reference plane, used by the tracker routines, is at zini=23.5. | 
| 2318 | double *dyout   = new double[npoint]; | * If you give as input the track state-vector from a TrkTrack object, | 
| 2319 | double *dthxout = new double[npoint]; | * you can use Trajectory::DoTrack(float* al) instead. | 
| 2320 | double *dthyout = new double[npoint]; | */ | 
| 2321 | double *dtlout  = new double[npoint]; | int Trajectory::DoTrack(float* al, float zini){ | 
| 2322 | double *dzin    = new double[npoint]; |  | 
| 2323 |  | //      double *dxout   = new double[npoint]; | 
| 2324 |  | //      double *dyout   = new double[npoint]; | 
| 2325 |  | //      double *dthxout = new double[npoint]; | 
| 2326 |  | //      double *dthyout = new double[npoint]; | 
| 2327 |  | //      double *dtlout  = new double[npoint]; | 
| 2328 |  | //      double *dzin    = new double[npoint]; | 
| 2329 |  |  | 
| 2330 |  | double *dxout; | 
| 2331 |  | double *dyout; | 
| 2332 |  | double *dthxout; | 
| 2333 |  | double *dthyout; | 
| 2334 |  | double *dtlout; | 
| 2335 |  | double *dzin; | 
| 2336 |  |  | 
| 2337 |  | dxout   = (double*) malloc(npoint*sizeof(double)); | 
| 2338 |  | dyout   = (double*) malloc(npoint*sizeof(double)); | 
| 2339 |  | dthxout = (double*) malloc(npoint*sizeof(double)); | 
| 2340 |  | dthyout = (double*) malloc(npoint*sizeof(double)); | 
| 2341 |  | dtlout  = (double*) malloc(npoint*sizeof(double)); | 
| 2342 |  | dzin    = (double*) malloc(npoint*sizeof(double)); | 
| 2343 |  |  | 
| 2344 | double dal[5]; | double dal[5]; | 
| 2345 |  |  | 
| 2346 | int ifail = 0; | double dzini = (double)zini; | 
| 2347 |  |  | 
| 2348 |  | int ifail = 0; | 
| 2349 |  |  | 
| 2350 | for (int i=0; i<5; i++)      dal[i]  = (double)al[i]; | for (int i=0; i<5; i++)      dal[i]  = (double)al[i]; | 
| 2351 | for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; | for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; | 
| 2352 |  |  | 
| 2353 | dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); | TrkParams::Load(1); | 
| 2354 |  | if( !TrkParams::IsLoaded(1) ){ | 
| 2355 |  | cout << "int Trajectory::DoTrack(float* al) --- ERROR --- m.field not loaded"<<endl; | 
| 2356 |  | return 0; | 
| 2357 |  | } | 
| 2358 |  | //    dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); | 
| 2359 |  | dotrack3_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&dzini,&ifail); | 
| 2360 |  |  | 
| 2361 | for (int i=0; i<npoint; i++){ | for (int i=0; i<npoint; i++){ | 
| 2362 | x[i]   = (float)*dxout++; | x[i]   = (float)*(dxout+i); | 
| 2363 | y[i]   = (float)*dyout++; | y[i]   = (float)*(dyout+i); | 
| 2364 | thx[i] = (float)*dthxout++; | thx[i] = (float)*(dthxout+i); | 
| 2365 | thy[i] = (float)*dthyout++; | thy[i] = (float)*(dthyout+i); | 
| 2366 | tl[i]  = (float)*dtlout++; | tl[i]  = (float)*(dtlout+i); | 
| 2367 | } | } | 
| 2368 |  |  | 
| 2369 |  | if(dxout)  free( dxout ); | 
| 2370 |  | if(dyout)  free( dyout ); | 
| 2371 |  | if(dthxout)free( dthxout ); | 
| 2372 |  | if(dthyout)free( dthyout ); | 
| 2373 |  | if(dtlout) free( dtlout ); | 
| 2374 |  | if(dzin)   free( dzin ); | 
| 2375 |  |  | 
| 2376 |  | //      delete [] dxout; | 
| 2377 |  | //      delete [] dyout; | 
| 2378 |  | //      delete [] dthxout; | 
| 2379 |  | //      delete [] dthyout; | 
| 2380 |  | //      delete [] dtlout; | 
| 2381 |  | //      delete [] dzin; | 
| 2382 |  |  | 
| 2383 |  |  | 
| 2384 | return ifail; | return ifail; | 
| 2385 | }; | }; | 
| 2386 |  |  | 
| 2387 |  | /** | 
| 2388 |  | * | 
| 2389 |  | * >>> OBSOLETE !!! use Trajectory::DoTrack(float* al, float zini) instead | 
| 2390 |  | * | 
| 2391 |  | */ | 
| 2392 |  | int Trajectory::DoTrack2(float* al, float zini){ | 
| 2393 |  |  | 
| 2394 |  | cout << endl; | 
| 2395 |  | cout << " int Trajectory::DoTrack2(float* al, float zini) --->> NB NB !! this method is going to be eliminated !!! "<<endl; | 
| 2396 |  | cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; | 
| 2397 |  | cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; | 
| 2398 |  | cout << endl; | 
| 2399 |  |  | 
| 2400 |  | return DoTrack(al,zini); | 
| 2401 |  |  | 
| 2402 |  | }; | 
| 2403 |  |  | 
| 2404 |  |  | 
| 2405 |  |  | 
| 2406 | ClassImp(TrkLevel2); | ClassImp(TrkLevel2); | 
| 2407 | ClassImp(TrkSinglet); | ClassImp(TrkSinglet); | 
| 2408 | ClassImp(TrkTrack); | ClassImp(TrkTrack); |