/* $Id: matching.c,v 1.25 2005/11/08 10:16:51 michael Exp $ $Source: /home/michael/cvsroot/matching_soft_v1/matching.c,v $ */ /**********************************************************/ /* */ /* Matching Module */ /* All code for the hit-matching */ /**********************************************************/ /* * * Some quick_info. * * */ #include #include #include "trigger.h" #include "matching.h" #include "scheduler.h" #define THETA_BIN_RESOLUTION_SHIFT 2 #define DELTA_THETA_BIN_RESOLUTION_SHIFT 2 #define PHI_BIN_RESOLUTION_SHIFT 4 #define DELTA_THETA_TABLE_WIDTH 64 #define PHI_TABLE_WIDTH 8 /* unsigned long PHI_DEVIATION_WINDOW; float INVARIANT_MASS_WINDOW; float INVARIANT_MASS; */ unsigned long Leptons[MATCHING_MAX_LEPTON_NUMBER+3]; unsigned long DiLeptons[14]; unsigned long scaler_tof_theta_cut = 0; unsigned long scaler_tof_phi_cut = 0; unsigned long scaler_shower_theta_cut = 0; unsigned long scaler_shower_phi_cut = 0; unsigned long found_number_of_leptons = 0; unsigned long found_number_of_dileptons = 0; // hit matching unsigned long hit_matching(void) { unsigned long sector_number,j,k; // unsigned long *curr_ptr; unsigned long RICH_theta, RICH_phi; unsigned long momentum_mapped_phi; unsigned long TOF_theta, TOF_phi; unsigned long SHOWER_theta, SHOWER_phi; unsigned long momentum; long phi_deviation, delta_theta; unsigned long theta_bin, phi_bin, delta_theta_bin; unsigned long electron_flag = 0; // bended to beam => positron unsigned long current_lepton_number = 0; unsigned long current_lepton_address = 0; unsigned long lepton_number_print_counter = 0; #ifdef ENABLE_LINEAR_MATCHING_WINDOW unsigned long phi_in_half_sector, sign_flag; unsigned long linear_matching_window; #endif found_number_of_leptons = 0; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: +Hit_matching", printf_writes_total); printf_vme(print_string); } #endif for(sector_number=0; sector_number<=5; sector_number++) { #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 3) { sprintf(print_string, ":%6d: sector: %d, RICH_hits_len: 0x%d, RICH_hits_ptr[sec]:%8.8x", printf_writes_total, sector_number, RICH_hits_length[sector_number], RICH_hits_ptr[sector_number] ); printf_vme(print_string); } #endif for(j=0; j< RICH_hits_length[sector_number] ; j++) { #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 2) { sprintf(print_string, ":%6d: RICH_hit sec: %d, num: %d", printf_writes_total, sector_number, j ); printf_vme(print_string); } #endif RICH_theta = (*(RICH_hits_ptr[sector_number]+j)) & 0xff; RICH_phi = ((*(RICH_hits_ptr[sector_number]+j)) >> 8) & 0x7ff; if( (RICH_phi&0xff) > 0x7f) { momentum_mapped_phi = RICH_phi & 0x7f; } else { momentum_mapped_phi = 0x80 - (RICH_phi&0xff); if(momentum_mapped_phi > 0x7f) { momentum_mapped_phi = 0x7f; } } theta_bin = (RICH_theta>>THETA_BIN_RESOLUTION_SHIFT); phi_bin = (momentum_mapped_phi>>PHI_BIN_RESOLUTION_SHIFT); #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 2) { sprintf(print_string, ":%6d: RICH_theta: 0x%x, RICH_phi: 0x%x, TOF_hits_length: %d", printf_writes_total, RICH_theta, RICH_phi , TOF_hits_length[sector_number] ); printf_vme(print_string); } #endif #ifdef ENABLE_LINEAR_MATCHING_WINDOW phi_in_half_sector = (RICH_phi&0xff); sign_flag = RICH_phi&0x80; // if sign flag is on, then hit is from 0 to 30 degrees if( sign_flag ) { phi_in_half_sector = RICH_phi&0x7f; } else { phi_in_half_sector = 0x7f - phi_in_half_sector; } linear_matching_window = (*PHI_DEVIATION_WINDOW) + (phi_in_half_sector >> (*G_PHI_MATCHING_SLOPE)); #endif for(k=0; k< TOF_hits_length[sector_number]; k++) { TOF_theta = (*(TOF_hits_ptr[sector_number]+k)) & 0xff; TOF_phi = ((*(TOF_hits_ptr[sector_number]+k)) >> 8) & 0x7ff; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: TOF number: %d, TOF_theta: 0x%x, TOF_phi: 0x%x", printf_writes_total, k, TOF_theta, TOF_phi ); printf_vme(print_string); } #endif phi_deviation = RICH_phi - TOF_phi; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 4) { sprintf(print_string, ":%6d: signed phi_deviation: 0x%8.8x", printf_writes_total, phi_deviation ); printf_vme(print_string); } #endif if(phi_deviation < 0) { // absolute phi_deviation=-phi_deviation; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 2) { sprintf(print_string, ":%6d: phi_deviation : 0x%x", printf_writes_total, phi_deviation ); printf_vme(print_string); } #endif #ifndef ENABLE_LINEAR_MATCHING_WINDOW if( phi_deviation < (*PHI_DEVIATION_WINDOW) ) { #else if( phi_deviation < linear_matching_window ) { #endif // found hit #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 2) { sprintf(print_string, ":%6d: found hits in phi window", printf_writes_total ); printf_vme(print_string); } #endif delta_theta = RICH_theta - TOF_theta; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 4) { sprintf(print_string, ":%6d: delta_theta: 0x%8.8x", printf_writes_total, delta_theta ); printf_vme(print_string); } #endif if(delta_theta<0) { // absolute electron_flag=1; // bended to larger angles delta_theta=-delta_theta; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) == 0x2) { sprintf(print_string, ":%6d: sign: delta_theta: 0x%8.8x", printf_writes_total, delta_theta ); printf_vme(print_string); } #endif if(delta_theta > (*DELTA_THETA_CUT)) { #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 0x1) { sprintf(print_string, ":%6d: delta_theta > CUT, no lepton", printf_writes_total ); printf_vme(print_string); } #endif scaler_tof_theta_cut++; continue; } delta_theta_bin = (delta_theta>>DELTA_THETA_BIN_RESOLUTION_SHIFT); // momentum_matrix momentum = *(MomentumMatrix + theta_bin*DELTA_THETA_TABLE_WIDTH*PHI_TABLE_WIDTH + delta_theta_bin*PHI_TABLE_WIDTH + phi_bin); /* // old format Leptons[current_lepton_address++] = ((momentum&0xff)<<24 ) | (electron_flag<<20) | (*(RICH_hits_ptr[sector_number]+j)); */ // (*RICH_hits_ptr[sector_number]+j)); // 31-24 23-22 21 20 19-12 11-4 3-0 // momentum reserved electron flag Detector Bit META Nr. RICH Nr. sector // Detector Bit: 0 => TOF // 1 => Shower // for RICH if(momentum >= 0) { Leptons[current_lepton_address++] = ((momentum&0xff)<< 24 ) | (electron_flag<<21) | (0x0 << 20) | (k<<12) | (j<<4) | (sector_number); current_lepton_number++; } if(current_lepton_number > MATCHING_MAX_LEPTON_NUMBER) { lepton_number_print_counter++; if(lepton_number_print_counter >=10) { sprintf(print_string, ":%6d:++attention: lepton number > %d: %d, return!", printf_writes_total, MATCHING_MAX_LEPTON_NUMBER, current_lepton_number ); printf_vme(print_string); lepton_number_print_counter = 0; } (*G_ACCUMULATED_LEPTONS)++; found_number_of_leptons = current_lepton_number; return 1; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) == 0x2) { sprintf(print_string, ":%6d: found lepton R-TOF: sec:%d, j:%d, TOF_hit_nr: %d, curr_lepton_nr:%d", printf_writes_total, sector_number,j, k, current_lepton_number ); printf_vme(print_string); sprintf(print_string, ":%6d: R-TOF: tb:%d, dtb:%d phi-bin:%d, momentum:%d, e_flag:%d", printf_writes_total, theta_bin, delta_theta_bin, phi_bin, momentum, electron_flag ); printf_vme(print_string); } #endif } else { // phi cut rejected event scaler_tof_phi_cut++; } } // same matching game with shower // stupid Michael! You have to reset this electron flag to 0 electron_flag = 0; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: SHOWER_hits_length[%d]: %x", printf_writes_total, sector_number, SHOWER_hits_length[sector_number] ); printf_vme(print_string); } #endif for(k=0; (k < SHOWER_hits_length[sector_number]) && (k < 100) ; k++) { SHOWER_theta = (*(SHOWER_hits_ptr[sector_number]+k)) & 0xff; SHOWER_phi = ((*(SHOWER_hits_ptr[sector_number]+k)) >> 8) & 0x7ff; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: sector: %d, SHW num: %d, SHW_theta: 0x%x, SHW_phi: 0x%x", printf_writes_total, sector_number, k, SHOWER_theta, SHOWER_phi ); printf_vme(print_string); } #endif phi_deviation = RICH_phi - SHOWER_phi; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) &0xf) > 4) { sprintf(print_string, ":%6d: signed phi_deviation: 0x%8.8x", printf_writes_total, phi_deviation ); printf_vme(print_string); } #endif if(phi_deviation < 0) { // absolute phi_deviation=-phi_deviation; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4)&0xf) > 2) { sprintf(print_string, ":%6d: phi_deviation : 0x%x", printf_writes_total, phi_deviation ); printf_vme(print_string); } #endif #ifndef ENABLE_LINEAR_MATCHING_WINDOW if( phi_deviation < (*PHI_DEVIATION_WINDOW) ) { // found hit #else if( phi_deviation < linear_matching_window ) { #endif #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 2) { sprintf(print_string, ":%6d: found hits in phi window", printf_writes_total ); printf_vme(print_string); } #endif delta_theta = RICH_theta - SHOWER_theta; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 4) { sprintf(print_string, ":%6d: delta_theta (sign): 0x%8.8x", printf_writes_total, delta_theta ); printf_vme(print_string); } #endif if(delta_theta<0) { // absolute electron_flag=1; // bended to larger angles delta_theta=-delta_theta; } if(delta_theta > (*DELTA_THETA_CUT)) { #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: delta_theta shower > CUT, no lepton", printf_writes_total ); printf_vme(print_string); } #endif continue; } delta_theta_bin = (delta_theta>>DELTA_THETA_BIN_RESOLUTION_SHIFT); // momentum_matrix momentum = *(MomentumMatrix + theta_bin*DELTA_THETA_TABLE_WIDTH*PHI_TABLE_WIDTH + delta_theta_bin*PHI_TABLE_WIDTH + phi_bin); /* Leptons[current_lepton_number++] = ((momentum&0xff)<<24 ) | (electron_flag<<20) | (*(RICH_hits_ptr[sector_number]+j+k)); */ // 31-24 23-22 21 20 19-12 11-4 3-0 // momentum reserved electron flag Detector Bit META Nr. RICH Nr. sector // Detector Bit: 0 => TOF // 1 => Shower // *** for Shower if(momentum >= 0) { Leptons[current_lepton_address++] = ((momentum&0xff)<< 24 ) | (electron_flag<<21) | (0x1 << 20) | (k<<12) | (j<<4) | (sector_number); current_lepton_number++; } if(current_lepton_number > MATCHING_MAX_LEPTON_NUMBER) { lepton_number_print_counter++; if(lepton_number_print_counter >=10) { sprintf(print_string, ":%6d:++attention: lepton number > %d: %d, return!", printf_writes_total, MATCHING_MAX_LEPTON_NUMBER, current_lepton_number ); printf_vme(print_string); lepton_number_print_counter = 0; } (*G_ACCUMULATED_LEPTONS)++; found_number_of_leptons = current_lepton_number; return 1; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: found lepton R-SH: sec:%d, j:%d, k:%d, curr_lepton:%d ", printf_writes_total, sector_number,j, k, current_lepton_number ); printf_vme(print_string); sprintf(print_string, ":%6d: R-SH: (dec) tbin:%d, dtbin:%d phi-bin:%d, momentum:%d, e_flag:%d", printf_writes_total, theta_bin, delta_theta_bin, phi_bin, momentum, electron_flag ); printf_vme(print_string); } #endif continue; } } // found all leptons } // end of loop over all RICH hits } // end of loop over all sectors (*G_ACCUMULATED_LEPTONS) += current_lepton_number; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: found %d leptons", printf_writes_total, current_lepton_number ); printf_vme(print_string); } #endif found_number_of_leptons = current_lepton_number; return 0; } // hit_matching(void) void generate_sin_table(float *ptr, unsigned long resolution) { long int i; float s; for(i=0;i<=resolution;i++) { s= 3.1415926/2.0 / (float)resolution * i; *(ptr+i) = (float) sinf(s); } } float tab_sin_8(unsigned long a) { float b, s; //unsigned long r; b= 255.0/(3.1415926/2.0); b=(float)a/b; s=sinf(b); // r = (unsigned long)(s*32768.0); return s; } float tab_sin_11(unsigned long a) { float b; //unsigned long r; /* 6*256 = 1536 1536 / 4 = 384 384*1 384 384*2 768 384*3 1152 */ b= 1536.0/(3.1415926/2.0); if(a>1151) { a=a-1280; b=(float)a/b; return (-1.0+sinf(b)); } else if(a>767){ a=a-1024; b=(float)a/b; return (-sinf(b)); } else if(a>383){ a=a-768; b=(float)a/b; return (1.0-sinf(b)); } else if(a<384){ b=(float)a/b; return sinf(b); } else { sprintf(print_string, ":%6d:++error: sin: a not in range: %d. halt", printf_writes_total, a); printf_vme(print_string); for(;;); } } unsigned long lepton_matching() { unsigned long i,j; unsigned long electron_flag_a, electron_flag_b; unsigned long curr_DiLepton_Number = 0 ; unsigned long curr_DiLepton_Address = 0 ; float x1,y1,z1, x2,y2,z2, cos_omega, invariant_mass_squared; float momentum1, momentum2; //unsigned long i_momentum1, i_momentum2; unsigned long theta1, theta2, phi1, phi2; unsigned long *help_ptr; #if 0 long delta_theta; #endif unsigned long trigger_decision = 0; unsigned long lepton_array_length; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: +Lepton_matching: nr_of_lep: %d", printf_writes_total, found_number_of_leptons); printf_vme(print_string); } #endif lepton_array_length = found_number_of_leptons; for(i=0;i>20)&0x1; for(j=i+1; j< lepton_array_length; j++) { electron_flag_b = ((Leptons[j])>>20)&0x1; // test --remove-- // electron_flag_b = 1; removed if( (electron_flag_a ^ electron_flag_b) ) { // found Di-Lepton #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: found Di-Lepton: i:%d, j:%d", printf_writes_total, i,j ); printf_vme(print_string); } #endif /// *(RICH_hits_ptr[sector_number]+j) // old theta2 = Leptons[j] & 0xff; // 31-24 23-22 21 20 19-12 11-4 3-0 // momentum reserved electron flag Detector Bit META Nr. RICH Nr. sector // Detector Bit: 0 => TOF // 1 => Shower theta1 = (*( RICH_hits_ptr[Leptons[i]&0xf] + ((Leptons[i]>>4)&0xff) )) & 0xff; theta2 = (*( RICH_hits_ptr[Leptons[j]&0xf] + ((Leptons[j]>>4)&0xff) )) & 0xff; // delta theta cut #if 0 delta_theta = theta1-theta2; if(delta_theta<0) { delta_theta=-delta_theta; } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: delta_theta: 0x%8.8x", printf_writes_total, delta_theta ); printf_vme(print_string); } #endif if(delta_theta > (*DELTA_THETA_CUT)) { #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: delta_theta > CUT, no delepton", printf_writes_total,theta1,phi1, theta2, phi2 ); printf_vme(print_string); } #endif continue; } #endif // delta theta cut // phi1 = (Leptons[i]>>8) & 0x7ff; // phi2 = (Leptons[j]>>8) & 0x7ff; phi1 = ((*( RICH_hits_ptr[Leptons[i]&0xf] + ((Leptons[i]>>4)&0xff) )) >> 8 ) & 0x7ff; phi2 = ((*( RICH_hits_ptr[Leptons[j]&0xf] + ((Leptons[j]>>4)&0xff) )) >> 8 ) & 0x7ff; momentum1 = (Leptons[i]>>24)&0xff; momentum2 = (Leptons[j]>>24)&0xff; #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: (hex) t1: %8.8x, p1: %8.8x, t2: %8.8x, p2: %8.8x", printf_writes_total,theta1,phi1, theta2, phi2 ); printf_vme(print_string); } #endif z1 = tab_sin_8(theta1); z2 = tab_sin_8(theta2); y1 = tab_sin_11(phi1) * tab_sin_8(0xff-theta1); y2 = tab_sin_11(phi2) * tab_sin_8(0xff-theta2); x1 = tab_sin_11(0x5ff-phi1) * tab_sin_8(0xff-theta1); x2 = tab_sin_11(0x5ff-phi2) * tab_sin_8(0xff-theta2); #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: x1: %8.8x, x2: %8.8x, y1: %8.8x, y2: %8.8x, z1: %8.8x, z2: %8.8x", printf_writes_total, x1, x2, y1, y2, z1, z2 ); printf_vme(print_string); } #endif cos_omega = x1*x2 + y1*y2 + z1*z2; //momentum1 invariant_mass_squared = 2 *momentum1 * momentum2 * (1 - cos_omega); #if 0 if(invariant_mass_squared < 0) { sprintf(print_string, ":%6d:++error: inv. mass < 0, m1: %8.8x, m2:%8.8x, cos_o: %8.8x", printf_writes_total, momentum1, momentum2, cos_omega); printf_vme(print_string); } #endif #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: m1: %8.8x, m2: %8.8x, cosw: %8.8x, imsqu: %8.8x", printf_writes_total, momentum1 ,momentum2, cos_omega, invariant_mass_squared ); printf_vme(print_string); } #endif if( invariant_mass_squared > (*INVARIANT_MASS_SQUARED_CUT) ) { // found dilepton in inv_mass window #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: found mass: %8.8x. higher cut", printf_writes_total, invariant_mass_squared ); printf_vme(print_string); } #endif } DiLeptons[curr_DiLepton_Address++] = (i) | (j<<8); help_ptr = (unsigned long *)&invariant_mass_squared; DiLeptons[curr_DiLepton_Address++] = *(help_ptr); curr_DiLepton_Number++; trigger_decision++; // make active // return 1; if(curr_DiLepton_Number > 6) { #if 0 sprintf(print_string, ":%6d:++attention: dilepton number > 6: %d, return!", printf_writes_total, curr_DiLepton_Number ); printf_vme(print_string); #endif (*G_ACCUMULATED_DILEPTONS)+= curr_DiLepton_Number; found_number_of_dileptons = curr_DiLepton_Number; return 1; } } } } #ifdef ENABLE_DEBUG_M if( (((*debug_register)>>4) & 0xf) > 1) { sprintf(print_string, ":%6d: found %d Dileptons below mass cut", printf_writes_total, curr_DiLepton_Number ); printf_vme(print_string); } #endif (*G_ACCUMULATED_DILEPTONS) += curr_DiLepton_Number; found_number_of_dileptons = curr_DiLepton_Number; if( found_number_of_dileptons != trigger_decision ) { sprintf(print_string, ":%6d:++internal error in number of Dileptons, a:%d , b:%d", printf_writes_total, found_number_of_dileptons, trigger_decision ); printf_vme(print_string); for(;;); } return 0; } // end of lepton_matching(void)