00001 00069 #include <iostream> 00070 #include "treedo.h" 00071 00072 using namespace std; 00073 00074 extern EUppLow operator++(enum EUppLow &rs, int ); 00075 extern ERegion operator++(enum ERegion &rs, int ); 00076 extern Etype operator++(enum Etype &rs, int ); 00077 extern Edir operator++(enum Edir &rs, int ); 00078 /*extern Eorientation operator++(enum Eorientation &rs, int );*/ 00079 00080 extern TreeLine *trelin; 00081 extern int trelinanz; 00082 extern treeregion *rcTreeRegion[2][3][4][4]; 00083 extern Det *rcDETRegion[2][3][4]; 00084 extern Options opt; 00085 int parttrackanz; 00086 //________________________________________________________________________ 00087 /*---------------------------------------------------------------------------*\ 00088 00089 Bcheck() - this function compares the momentum look-up energy with the 00090 shooting method result. This function is for debugging 00091 purposes only. 00092 00093 inputs: (1) double E - the momentum lookup energy 00094 (2) PartTrack *f - pointer to the front partial track for the 00095 track 00096 (3) PartTrack *b - pointer to the back partial track of the 00097 track 00098 (4) double TVertex - transverse position of the vertex for the 00099 track 00100 (5) double ZVertex - Z position of the vertex for the track 00101 00102 outputs: an ASCII ntuple file. 00103 00104 \*---------------------------------------------------------------------------*/ 00105 double rcShootP( double Pest,PartTrack *front, PartTrack *back, double accuracy){ 00106 cerr << "This function is only a stub" << endl; 00107 return -1000; 00108 } 00109 void treedo::BCheck( double E, PartTrack *f, PartTrack *b, double TVertex, double ZVertex) 00110 { 00111 double Es = rcShootP(0.0,f,b,0.005); 00112 //extern physic phys_carlo; 00113 static FILE *test = 0; 00114 00115 if( !test ) 00116 test = fopen("magcheck","w"); 00117 if( test && E > 0.0) { 00118 fprintf(test,"%f %f %f %f %f %f\n", 00119 /*phys_carlo.E,*/E,Es,f->mx,f->my, TVertex, ZVertex); 00120 } 00121 } 00122 //________________________________________________________________________ 00123 /*---------------------------------------------------------------------------*\ 00124 00125 rcLinkUsedTracks() - this function links the track cross-linked list to 00126 a single-linked list to ease the access for later 00127 functions. The new link is done via the usenext 00128 pointer. 00129 00130 inputs: (1) Track *track - 00131 (2) int upplow - 00132 00133 returns: the head to the linked list. 00134 00135 \*---------------------------------------------------------------------------*/ 00136 00137 00138 Track * treedo::rcLinkUsedTracks( Track *track, int upplow ) 00139 { 00140 Track *ret = 0, *usedwalk = 0; 00141 Track *trackwalk, *ytrack; 00142 00143 /* loop over all found tracks */ 00144 for(trackwalk = track; trackwalk; trackwalk = trackwalk->next ) { 00145 /* and the possible y-tracks */ 00146 for( ytrack = trackwalk; ytrack; ytrack = ytrack->ynext ) { 00147 //Statist[ytrack->method].TracksGenerated[upplow] ++; 00148 if(!ytrack->Used) 00149 continue; 00150 //Statist[ytrack->method].TracksUsed[upplow] ++; 00151 if( !ret ) /* return the first used track */ 00152 ret = usedwalk = ytrack; 00153 else { 00154 usedwalk->usednext = ytrack; /* and link them together */ 00155 usedwalk = ytrack; 00156 } 00157 } 00158 } 00159 return ret; 00160 } 00161 //________________________________________________________________________ 00162 /*---------------------------------------------------------------------------*\ 00163 00164 rcTreeDo() - this function is the main tracking function. It 00165 performes tracking in projections (u/v/x) to form treelines, 00166 it combines projection tracks to tracks in space and bridges 00167 tracks in space before and after the magnet to form recon- 00168 structed particle tracks. The function calculates the momentum 00169 of the tracks and afterwards reiterates the track parameters 00170 using the now know track position and momentum for iterating 00171 the hit positions in the tracking chambers. 00172 00173 inputs: (1) int iEventNo - the event number 00174 00175 outputs: (1) Event *rcTreeDo() - pointer to the structure with all 00176 of the reconstruction information for 00177 this event. 00178 00179 \*---------------------------------------------------------------------------*/ 00180 Event * treedo::rcTreeDo(int iEventNo){ 00181 cerr << "Initialize" << endl; 00182 enum EUppLow upplow; /* Loop counter for TOP and BOTTOM detectors */ 00183 // enum Emethod method; /* Loop counter over reconstruction method */ 00184 // enum EFrontBack frback; /* Loop counter over FRONT or BACK regions */ 00185 enum Edir dir; /* Loop counter over wire pitches (U,V,X) */ 00186 enum ERegion region; 00187 enum Etype type; 00188 /* enum Eorientation orient;*/ 00189 00190 int i,k; 00191 int dlayer = 0; /* number of detector planes in the search */ 00192 double A[3][2]; /* conversion between xy and uv */ 00193 Event *event; /* area for storing the reconstruction info */ 00194 PartTrack *area = 0; 00195 Det *rd/*, *rnd*/; /* pointers for moving through the linked 00196 lists of detector id and hit information */ 00197 TreeLine *trelin1,*trelin2; 00198 00199 treesearch TreeSearch; 00200 treecombine TreeCombine; 00201 treesort TreeSort; 00202 treematch TreeMatch; 00203 00204 /*int charge; 00205 int found_front = 0; 00206 double ZVertex, bending, theta, phi, ESpec; 00207 int outbounds = 0;*/ 00208 int levelmax; 00209 static int init = 0; 00210 static char *channel[TLAYERS]; 00211 static int *hashchannel[TLAYERS]; 00212 int rcSETiLevels0 = 4; 00213 const int numWiresr3 = 282; 00214 static char *channelr3[numWiresr3]; 00215 static int *hashchannelr3[numWiresr3]; 00216 00217 event = (Event *)calloc(1,sizeof(Event)); 00218 assert(event); 00219 // Initialize stuff if first event 00220 if( !init ) { 00221 init++; /* don't initialize again */ 00222 levelmax = rcSETiLevels0; /* determine the bin-division depth 00223 of the tree-search */ 00224 for( i = 0; i < TLAYERS; i++ ) { /* reserve space for the bit patterns */ 00225 channel[i] = (char*)malloc( 1UL << levelmax ); 00226 hashchannel[i] = (int*)malloc( (sizeof(int)*(1UL << (levelmax-1) ))); 00227 assert( channel[i] && hashchannel[i] ); 00228 } 00229 for( i = 0; i < numWiresr3; i++){ /* reserve space for region 3 bit patterns */ 00230 channelr3[i] = (char*)malloc( 1UL << (opt.levels[0][2][0])); 00231 hashchannelr3[i] = (int*)malloc((sizeof(int)*(1UL << (opt.levels[0][2][0]-1)))); 00232 assert(channelr3[i] && hashchannelr3[i]); 00233 } 00234 00235 } 00236 00237 00238 00239 if(!init){ 00240 //do stuff 00241 init++; 00242 } 00243 //loop through all detectors/regions/wire directions 00244 for(upplow = w_upper;upplow <=w_upper/*w_lower*/;upplow++){ 00245 for(region = r3/*r1*/;region<=r3;region++){ 00246 for(type=d_drift;type<=d_drift/*d_cerenkov*/;type++){ 00247 parttrackanz = 0;/* no partial tracks found yet */ 00248 for(dir=u_dir;dir<=v_dir/*y_dir*/;dir++){ 00249 /* ---- 1st: check that the direction is tree-searchable ---- */ 00250 if(rcTreeRegion[upplow][region-1][type][dir]->searchable == false){ 00251 (event->treeline[upplow][region-1][type][dir])=0; 00252 printf("%i %i %i %i\n",upplow,region,type,dir); 00253 continue;//'searchable' is set by tree.cc 00254 } 00255 00256 /* ---- 2nd: do some variable initialization for the found treeline 00257 linked list ---- */ 00258 00259 trelin = 0; /* clear pointer to linked list of treelines */ 00260 trelinanz = 0; /* clear the "treelines found" counter */ 00261 00262 /* ---- 3rd: create the bit patterns for the hits ---- */ 00263 00264 //_______________________________________ 00265 if(region == r3){ 00266 dlayer = 0; /* set "number of detectors" to zero */ 00267 int decrease; 00268 trelin1 = 0;trelin2 = 0; 00269 for(k=0,rd = rcDETRegion[upplow][region-1][dir],decrease = 0; 00270 rd; rd = rd->nextsame,decrease+=282,k++ ) { /* over the like-pitched planes in a region */ 00271 trelin = 0;//clear the linked list for the 2nd layer 00272 Hit *H,*Hnext; 00273 A[dir][0] = rd->rCos; /* cos(angle of wire pitch) */ 00274 A[dir][1] = rd->rSin; /* sin(angle of wire pitch) */ 00275 H = rd->hitbydet; 00276 cerr << "Set Pattern Hits" << endl; 00277 for( i = 0; i < numWiresr3; i++){ 00278 memset(channelr3[i],0,1UL<<(opt.levels[0][2][0])); 00279 memset(hashchannelr3[i], 0, sizeof(int)* 00280 (1UL<<(opt.levels[0][2][0]-1))); 00281 } 00282 while(H && H->wire < 282+decrease){ 00283 TreeSearch.TsSetPoint(rcTreeRegion[upplow][region-1][type][dir]->rWidth,H,channelr3[H->wire-decrease],hashchannelr3[H->wire-decrease],1U<<(opt.levels[upplow][region-1][type]-1)); 00284 //for(int i=0;i<8;i++)cerr << hashchannelr3[H->wire-decrease][i] << " ";cerr << endl; 00285 //cerr << H->wire-decrease << ":"; 00286 if(opt.showEventPattern){ 00287 cerr << "w" << H->wire << ":"; 00288 for(int i=0;i<(1UL<<(opt.levels[upplow][region-1][type]))-1;i++){ 00289 00290 cerr << "<" ; 00291 cerr << channelr3[H->wire-decrease][i] ; 00292 cerr << ">"; 00293 } 00294 cerr << endl; 00295 } 00296 Hnext = H->nextdet; 00297 H = Hnext; 00298 } 00299 if(1)cerr << "Search for Matching Patterns " << endl; 00300 00301 TreeSearch.TsSearch(&(rcTreeRegion[upplow][region-1][type][dir]->node), 00302 channelr3, hashchannelr3, opt.levels[upplow][region-1][type],numWiresr3); 00303 00304 00305 if(1)cerr << "Sort Patterns " << endl; 00306 if( rcTreeRegion[upplow][region-1][type][dir] ){//<- This if statement may be done wrong 00307 if(k==0){ 00308 TreeCombine.TlTreeLineSort(trelin,upplow,region,type,dir, 00309 1UL<<(opt.levels[upplow][region-1][type]-1),0,dlayer); 00310 trelin1 = trelin; 00311 } 00312 else if(k==1){ 00313 TreeCombine.TlTreeLineSort(trelin,upplow,region,type,dir, 00314 1UL<<(opt.levels[upplow][region-1][type]-1),0,dlayer); 00315 trelin2 = trelin; 00316 } 00317 00318 } 00319 dlayer++; 00320 00321 } 00322 if(1)cerr << "Match Region 3 Segments" << endl; 00323 00324 trelin = TreeMatch.MatchR3(trelin1, 00325 trelin2,upplow,region,dir); 00326 event->treeline[upplow][region-1][type][dir] = trelin; 00327 tlayers = TLAYERS; /* remember the number of tree-detector */ 00328 tlaym1 = tlayers - 1; /* remember tlayers-1 for convenience */ 00329 } 00330 else{ 00331 cerr << "Error, no code for this type of detector here\n"; 00332 return event; 00333 } 00334 00335 00336 00337 //_______________________________________ 00338 /* 00339 00340 00341 tlayer = 0; /* set "number of tree-detectors" to zero * 00342 dlayer = 0; /* set "number of detectors" to zero * 00343 for( i= 0;i<numWiresr3;i++){//this is redundant 00344 memset(channelr3[i],0,1UL<<(opt.levels[upplow][region-1][type])); 00345 memset(hashchannelr3[i], 0, sizeof(int)*(1UL<<(opt.levels[upplow][region-1][type]-1))); 00346 } 00347 00348 for(rd = rcDETRegion[upplow][region-1][dir]; 00349 rd; rd = rd->nextsame ) { /* over the like-pitched planes in a region * 00350 if(rd->type == d_drift && region == r3){ 00351 Hit *H,*Hnext; 00352 A[dir][0] = rd->rCos; /* cos(angle of wire pitch) * 00353 A[dir][1] = rd->rSin; /* sin(angle of wire pitch) * 00354 H = rd->hitbydet; 00355 while(H){ 00356 TreeSearch.TsSetPoint(rcTreeRegion[upplow][region-1][type][dir]->rWidth,H,channelr3[H->wire],hashchannelr3[H->wire],1U<<(opt.levels[upplow][region-1][type]-1)); 00357 //for(int i=0;i<8;i++)cerr << hashchannelr3[H->wire][i] << " ";cerr << endl; 00358 //cerr << H->wire << ":"; 00359 //for(int i=0;i<15;i++)cerr << "<" << channelr3[H->wire][i] << ">";cerr << endl; 00360 00361 Hnext = H->nextdet; 00362 H = Hnext; 00363 00364 } 00365 } 00366 00367 else{ 00368 cerr << "Error, no code for this type of detector here\n"; 00369 return event; 00370 } 00371 tlayer++; 00372 dlayer++; 00373 }/* end of loop over the planes in the region * 00374 //cerr << "tlayer = " << tlayer << endl; 00375 tlayers = tlayer; /* remember the number of tree-detector * 00376 tlaym1 = tlayer - 1; /* remember tlayers-1 for convenience * 00377 ________________________________________________*/ 00378 /* ---- 4th: now the hit patterns for the planes with wires of a 00379 the same pitch have been constructed for a region of 00380 the detector, so apply the treesearch algorithm to 00381 these hit patterns to find the possible candidates for 00382 tracks (i.e. the treelines). The found treelines are 00383 stored in a linked list and the pointer to this list 00384 is put into rcTreeRegion. ---- */ 00385 00386 00387 00388 if(region == r3){ 00389 //for(int ii=0;ii<10;ii++){ 00390 //for(int i=0;i<15;i++)cerr << "<" << channelr3[ii][i] << ">";cerr << endl;} 00391 //cerr << "---" << endl; 00392 //TreeSearch.TsSearch(&(rcTreeRegion[upplow][region-1][type][dir]->node), 00393 // channelr3, hashchannelr3, opt.levels[upplow][region-1][type],numWiresr3); 00394 } 00395 else{ 00396 //TreeSearch.TsSearch(&(rcTreeRegion[upplow][region-1][type][dir]->node), 00397 // channel, hashchannel, opt.levels[upplow][region-1][type],0); 00398 00399 } 00400 00401 if( trelinanz > TREELINABORT ) { /* Too many treelines */ 00402 /* so drop the event */ 00403 return 0; 00404 } 00405 00406 /* ---- 5th: sort and search for ghosts in the resulting lines ---- */ 00407 /* 00408 if( rcTreeRegion[upplow][region-1][type][dir] ){//<- This if statement may be done wrong 00409 TreeCombine.TlTreeLineSort(trelin,upplow,region,type,dir, 00410 1UL<<(opt.levels[upplow][region-1][type]-1),tlayer,dlayer); 00411 } 00412 event->treeline[upplow][region-1][type][dir] = trelin; 00413 */ 00414 }/* end of loop over the three wire-pitch directions */ 00415 //cerr << "event : " << event->treeline[upplow][region-1][type][u_dir]->numhits << endl; 00416 /* ---- TASK 2: Combine the treelines into partial tracks ---- */ 00417 if( rcTreeRegion[upplow][region-1][type] && tlayers) {//<- This if statement may be done wrong 00418 area = TreeCombine.TlTreeCombine(event->treeline[upplow][region-1][type], 00419 1L<<(opt.levels[upplow][region-1][type]-1), 00420 upplow, 00421 region, 00422 type, 00423 tlayers, 00424 dlayer); 00425 } 00426 00427 if( parttrackanz > 600 ) { /* Adamo can't stand that */ 00428 cerr << "Too many partial tracks " << endl; 00429 return 0; 00430 } 00431 /* ---- TASK 3: Sort out the Partial Tracks ---- */ 00432 TreeSort.rcPartConnSort(area); 00433 /* ---- TASK 4: Hook up the partial track info to the event info ---- */ 00434 event->parttrack[upplow][region][type] = area; 00435 }}/* end of loop over the detectors */ 00436 /* end of loop over the regions */ 00437 00438 /* ============================== 00439 * Correlate front and back 00440 * tracks from x, y and y' infor- 00441 * mation 00442 * ============================== */ 00443 00444 00445 00446 00447 } 00448 return event; 00449 } 00450 //________________________________________________________________________
1.4.6