scatter_exp.C Source File

00001 
00002 #ifdef USE_MPI
00003 #include "mpi++.h"
00004 #else
00005 #include "localmpi++.h"
00006 #endif
00007 
00008 #include "scatter_exp.h"
00009 
00010 void scatter_exp::initialize_to_zero() {
00011 
00012         
00013   //  initial_form = new char[255];
00014   //  final_form = new char[255];
00015 
00016   int i;
00017   for(i=0; i< N_RHO_ZONE_MAX; i++) 
00018     n_hits[i] = n_hit[i] = 0;
00019   for(i=0; i<N_STEP_BIN; i++)
00020     step_counter[i] = 0;
00021   for(i=0; i<N_OSC_BIN; i++)
00022     osc_counter[i] = 0;
00023 
00024   n_coll = 0;
00025   n_star = 0;
00026 
00027   energy_error = 0;
00028 
00029   final_bound = false;
00030   resonance = false;
00031   stop = false;
00032   sd = not_specified; 
00033   n_zone = 0;
00034   form_changes = 0;
00035   n_initial = 0;
00036   n_final = 0;
00037   id_scenario = 0;
00038   n_found = 0;
00039   sigma = 0;
00040   sigma_err_sq = 0; 
00041 }
00042 
00043 real scatter_exp::get_specific_sigma_err_sq(scatter_discriptor d) {
00044 
00045     if(d == sd) 
00046       return sigma_err_sq;
00047     else 
00048       return 0;
00049 }
00050 
00051 real scatter_exp::get_specific_sigma(scatter_discriptor d) {
00052 
00053     if(d == sd) 
00054       return sigma;
00055     else 
00056       return 0;
00057 }
00058 
00059 int scatter_exp::get_specific_counts(scatter_discriptor d) {
00060 
00061     if(d == sd) 
00062       return n_found;
00063     else 
00064       return 0;
00065 }
00066 
00067 int scatter_exp::get_specific_counts(scatter_discriptor d, int i) {
00068 
00069     if(d == sd) 
00070       return n_hits[i];
00071     else 
00072       return 0;
00073 }
00074 
00075 bool scatter_exp::identical_to(scatter_exp *ha) {
00076 
00077   bool the_same = false;
00078 
00079   if (!strcmp(&final_form[0], ha->get_final_form()) && 
00080       resonance == ha->get_resonance() &&
00081       sd == ha->get_scatter_discriptor()) {
00082 
00083     the_same = true;
00084     set_nzone(ha->get_nzone());
00085   }
00086 
00087   return the_same;
00088 }
00089 
00090 bool scatter_exp::operator == (scatter_exp& ha) const {
00091 
00092   bool the_same = false;
00093   if(strcmp(&initial_form[0], ha.get_initial_form())==0 && 
00094      strcmp(&final_form[0], ha.get_final_form())==0 &&
00095      resonance == ha.get_resonance() &&
00096     sd == ha.get_scatter_discriptor()) {
00097       
00098     the_same = true;
00099   }
00100 
00101   return the_same;
00102 }
00103 
00104 bool scatter_exp::operator != (scatter_exp& ha) const {
00105 
00106   bool different = true;
00107   if (*this == ha)
00108     different = false;
00109 
00110   return different;
00111 }
00112  
00113 local char * type_string(scatter_discriptor sd) {
00114   
00115   switch(sd) {
00116   case preservation:                     return "pres";
00117        break;
00118   case exchange:                         return "exch";
00119        break;
00120   case exchange_ionization:              return "ex_ion";
00121        break;
00122   case single_ionization:                return "s_ion";
00123        break;
00124   case multiple_ionization:              return "m_ion";
00125        break;
00126   case total_ionization:                 return "t_ion";
00127        break;
00128   case collision_binary:                 return "bcoll";
00129        break;
00130   case collision_ionization:             return "c_ion";
00131        break;
00132   case stable_triple:                    return "trip";
00133        break;
00134   case stable_higher_order:              return "ntup";
00135        break;
00136   case two_body_collision:               return "2_coll";
00137        break;
00138   case multiple_body_collision:          return "n_coll";
00139        break;
00140   case unknown:                          return "unknown";
00141        break;
00142   case error:                            return "error";
00143        break;
00144   case stopped:                          return "stopped";
00145        break;
00146   case not_specified:                    return "NAS";
00147        break;
00148   case number_of_scatter_discriptors:    return "N";
00149        break;
00150   default:                               return "???";
00151   };
00152 }
00153 
00154 //char dummy_string[74];
00155 
00156 ostream& operator<<(ostream& s, scatter_exp& hi) {
00157 
00158 #if 0
00159   s << hi.id_scenario << " " << hi.n_found << " "
00160     << hi.sd << " " <<hi.resonance << " "
00161     << hi.get_final_form() << " \t";
00162 
00163   for(int i=0; i<N_RHO_ZONE_MAX; i++) 
00164     s << hi.get_nhits(i) << " ";
00165   s << endl;
00166 #endif
00167 
00168   int i;
00169   s << hi.id_scenario;
00170   for(i=0; i<3-floor(log10(hi.id_scenario)+1); i++)
00171     s << " ";
00172   s << hi.n_found;
00173   for(i=0; i<2-floor(log10(hi.n_found+1)); i++)
00174     s << " ";
00175   s << type_string(hi.sd);
00176   for(i=0; i<8-strlen(type_string(hi.sd)); i++)
00177     s << " ";
00178   if(hi.resonance)
00179     s << "r  ";
00180   else
00181     s << "   ";
00182   int istrl = strlen(hi.get_initial_form());
00183   for(i=0; i<istrl-strlen(hi.get_final_form()); i++)
00184     s << " ";
00185   s << hi.get_final_form();
00186   //       for(i=0; i<N_RHO_ZONE_MAX; i++) 
00187   //     s << hi.get_nhits(i);
00188 
00189   
00190   int n_space = 0;
00191   for(i=0; i<N_RHO_ZONE_MAX; i++) 
00192     n_space += (int)floor(log10(hi.get_nhits(i)+1));
00193   for(i=0; i<max(0,N_RHO_ZONE_MAX-n_space-15); i++) 
00194     s << ".";
00195 
00196   for(i=0; i<N_RHO_ZONE_MAX; i++) {
00197     for(int j=0; j<1-floor(log10(hi.get_nhits(i)+1)); j++)
00198       s << " ";
00199     s << hi.get_nhits(i);
00200   }
00201 
00202   return s;
00203 
00204 }
00205 
00206 scatter_exp* read_scatter_exp(istream& s) {
00207 
00208   scatter_exp* hi = new scatter_exp();
00209     if(s.eof() || hi == NULL) 
00210       return NULL;
00211     
00212     int id;
00213     char iform[255];
00214     char fform[255];
00215     s >> id >> &iform[0] >> &fform[0]; 
00216 
00217     hi->set_id_scenario(id);
00218     hi->set_final_form(fform);
00219     hi->set_initial_form(iform);
00220 
00221     return hi;
00222 }
00223 
00224 
00225 void scatter_exp::init_scatter_exp(sdyn* b) {
00226 
00227   time = 0;
00228   form_changes = 0;
00229   resonance = false;
00230   sd = not_specified; 
00231   //  initial_form = new char[255];
00232   strcpy(&initial_form[0], get_normal_form(b));
00233   n_initial = b->n_leaves();
00234   //  final_form = NULL;
00235 }
00236 
00237 void scatter_exp::final_scatter_exp(sdyn* b) {
00238 
00239   if(form_changes>1)
00240     resonance = true;
00241 
00242   time = b->get_time();
00243   n_star = b->n_leaves();
00244   strcpy(&final_form[0], get_normal_form(b));
00245   // is the final system a binary?
00246   n_final = b->n_leaves();
00247   sd = classify_scatter();
00248 
00249 }
00250 
00251 
00252 int scatter_exp::count_character_in_string(char string[], char search) {
00253 
00254 
00255   int counts = 0;
00256   for(int i=0; i<strlen(string); i++) {
00257     if(strncmp(&string[i], &search, 1)==0)
00258       counts++;
00259   }
00260 
00261   return counts;
00262 
00263 }
00264 
00265 int scatter_exp::count_character_in_string(char string[], char search[],
00266                                            int n_char) {
00267 
00268 
00269   int counts = 0;
00270   for(int i=0; i<strlen(string); i++) {
00271     if(strncmp(&string[i], &search[0], n_char)==0)
00272       counts++;
00273   }
00274 
00275   return counts;
00276 
00277 }
00278 
00279 // returns 3 if triple, 4 if quadruple etc.
00280 int scatter_exp::count_multiplicity() {
00281 
00282   char *final = get_final_form();
00283   int n_fstr = strlen(final);
00284 
00285   int b=0, bn=0;
00286   for(int i=0; i<n_fstr; i++) {
00287     if (strncmp(&final[i], "(", 1)==0)
00288       b++;
00289     else if(strncmp(&final[i], ")", 1)==0)
00290       b--;
00291     bn = max(bn, b);
00292   }
00293   if(b!=0)
00294     cerr << "Not encough braces found in scatter_exp::count_mulitplicity()"
00295          << endl;
00296 
00297   return bn;
00298 }
00299 
00300 // At the moment quite rudimentary
00301 bool scatter_exp::check_for_exchange() {
00302 
00303   bool exchange = false;
00304 
00305   char *init = get_initial_form();
00306   char *final = get_final_form();
00307   int n_istr = strlen(init);
00308   int n_fstr = strlen(final);
00309   char inames[64];
00310   char fnames[64];
00311 
00312   int icnt=0;
00313   for(int i=0; i<n_istr; i++) {
00314     if (!(strncmp(&init[i], "(", 1)==0 ||
00315           strncmp(&init[i], ")", 1)==0 ||
00316           strncmp(&init[i], "+", 1)==0 ||
00317           strncmp(&init[i], ",", 1)==0)) {
00318       sprintf(&inames[icnt], "%c", init[i]);
00319       icnt++;
00320     }
00321   }
00322 
00323   int fcnt=0;
00324   for(int f=0; f<n_fstr; f++) {
00325     if (!(strncmp(&final[f], "(", 1)==0 ||
00326           strncmp(&final[f], ")", 1)==0 ||
00327           strncmp(&final[f], "+", 1)==0 ||
00328           strncmp(&final[f], ",", 1)==0)) {
00329       sprintf(&fnames[fcnt], "%c", final[f]);
00330       fcnt++;
00331     }
00332   }
00333 
00334   if(icnt!=fcnt) {
00335     cerr << "ERROR: icnt (" << icnt << ") != fcnt (" << fcnt << ")" << endl;
00336     cerr << "       in   bool scatter_exp::check_for_exchange()" << endl;
00337     icnt = min(icnt, fcnt);
00338   }
00339   if(strncmp(inames, fnames, icnt)!=0)
00340     exchange = true;
00341 
00342   return exchange;
00343 }
00344 
00345 
00346 scatter_discriptor scatter_exp::classify_scatter() {
00347 
00348   scatter_discriptor discriptor = preservation;
00349 
00350   if(sd == stopped) {
00351     discriptor = stopped;
00352     stop = true;
00353   }
00354   else if(abs(energy_error)>1.e-5)
00355     discriptor = error; 
00356   else if(strcmp(get_initial_form(), get_final_form())==0) 
00357     discriptor = preservation;                           // ((p1,p2),(t1,t2))
00358   else if(count_character_in_string(get_final_form(), '+')>=1 &&
00359           count_character_in_string(get_final_form(), '(')==1) {
00360     if (get_n_final()==2 && get_final_bound())
00361       discriptor = collision_binary;                    // (p1+p2,t1+t2)
00362     else
00363       discriptor = collision_ionization;                // (p1+p2,t1,t2)
00364   }
00365   else if(count_character_in_string(get_final_form(), '+')>=1 &&
00366           count_character_in_string(get_final_form(), '(')>=2)
00367     discriptor = collision_binary;                      // (p1+p2,(t1,t2))
00368   else if(count_multiplicity()==3) 
00369     discriptor = stable_triple;                         // (p1,(p2,(t1,t2)))
00370   else if(count_multiplicity()>3) 
00371     discriptor = stable_higher_order;               // (p11,(p12,(p2,(t1,t2))))
00372   else if(count_character_in_string(get_final_form(), '(')==1)
00373     discriptor = total_ionization;                      // (p1,p2,t1,t2)
00374   else if(check_for_exchange()) {
00375     if (count_character_in_string(get_initial_form(), '(') >
00376         count_character_in_string(get_final_form(), '('))
00377       discriptor = exchange_ionization;                 // (p1,t1,(p2,t2))
00378     else
00379       discriptor = exchange;                            // ((p1,t1),(p2,t2))
00380   }
00381   else if(count_character_in_string(get_initial_form(), '(') >
00382           count_character_in_string(get_final_form(), '(')+1)
00383     discriptor = multiple_ionization;                   // (p11,p12,p2,t1,t2)
00384   else if(count_character_in_string(get_initial_form(), '(') >
00385           count_character_in_string(get_final_form(), '('))
00386     discriptor = single_ionization;                     // (p1,p2,(t1,t2))
00387   else
00388     discriptor = unknown;
00389 
00390   n_coll = count_character_in_string(get_final_form(), '+');
00391 
00392   return discriptor;
00393 }
00394 
00395 #ifdef USE_MPI
00396 
00397 MPI_Datatype scatter_exp::initialize_data_structures_MPI() {
00398 
00399   // MPI definition of datatype scatter_exp
00400   int charlength = 2*255;
00401   int reallength = 4;
00402   int intlength = 10 + 2*N_RHO_ZONE_MAX + N_STEP_BIN + N_OSC_BIN;
00403   int boolength = 3;
00404   int blockcounts[3] = {charlength, reallength, intlength+boolength};
00405   MPI::Aint exp_displs[3];
00406   MPI_Datatype scatter_exp_type;
00407 
00408   MPI_Address(get_initial_form(), &exp_displs[0]);
00409   MPI_Address(get_time_ptr(), &exp_displs[1]);
00410   MPI_Address(get_scatter_discriptor_ptr(), &exp_displs[2]);
00411   MPI_Datatype nexp_types[3] = {MPI::CHAR, MPI::DOUBLE, MPI::INT};
00412   exp_displs[1] -= exp_displs[0];
00413   exp_displs[2] -= exp_displs[0];
00414   exp_displs[0] = 0;
00415 
00416   MPI_Type_struct(3, blockcounts, exp_displs, nexp_types, &scatter_exp_type);
00417   MPI_Type_commit(&scatter_exp_type);
00418 
00419   return scatter_exp_type;
00420 }
00421 
00422 #else
00423 
00424 MPI_Datatype scatter_exp::initialize_data_structures_MPI() {
00425   MPI_Datatype dummy;
00426   return dummy;
00427 }
00428 
00429 #endif
00430 
00431 #ifdef TOOLBOX
00432 
00433 void main(int argc, char **argv) {
00434 
00435     check_help();
00436     extern char *poptarg;
00437     int c;
00438     char* param_string = "c:";
00439 
00440     real cpu_time_check;
00441     while ((c = pgetopt(argc, argv, param_string)) != -1)
00442         switch(c) {
00443             case 'c': cpu_time_check = atof(poptarg);
00444                       break;
00445             case '?': params_to_usage(cerr, argv[0], param_string);
00446                       get_help();
00447         }            
00448 
00449     do {
00450  
00451       //      scatter_exp *hi = read_scatter_exp(cin);
00452       scatter_exp *hi = new scatter_exp(cin);
00453       cerr << *hi << endl; 
00454       delete hi; 
00455     }
00456     while(!cin.eof());
00457 }
00458 
00459 #endif