//* Nilanga Liyanage //* Jefferson Lab //*-*-*-*-*-*-*-*-*-*-*-*Hall A database optimization *-*-*-*-*-*-*-*-* //*-* ============================ * //*-* * //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* #include "TROOT.h" #include "TMinuit.h" #include "TMath.h" #include "TStopwatch.h" #include "TFile.h" #include "TString.h" #include #include #include #include #include "opt_parameter.h" #include "spec.h" #include "Database.h" #include "Event.h" #include "Peak.h" using namespace std; extern void write_test_kumac(double beamx, double beamy); extern void compare_peak( int itr, double* sum_x, double* sum_y, double* sum_phi, double* sum_theta, double* sum_chi, int* sum_event, double beamx, double beamy ); void fcnk0(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag); double mat_mult(int iarm, int i_evt, int ii_line, int j_powx); double mat_mult_old(int iarm, int i_evt, int ii_line, int j_powx); void fp_vertex(int iarm, int i_evt); void usage(void); void getargs( int argc, char **argv ); database db; static data dat; static spec sp[3]; static int nvary; static bool lTest = false; // Test mode static bool lESPACE = true; // Use ESPACE databases by defaults (for now) static bool lPhi = true; // Phi angle mode (default) static bool lTheta = false; // Theta angle mode static bool lCalcXtg = false; // Calculate x_tg internally static bool lSiv = true; // Use sieve in ytg calibration static const char* spectro_names[] = { "R", "L", 0 }; // Helper classes for reading/writing matrix elements from/to databases class DatabaseIO { public: enum EMode { kRead, kWrite }; int read() { return do_io(kRead); } int write() { return do_io(kWrite); } protected: virtual int do_io( EMode mode ) = 0; }; class cpp_db : public DatabaseIO { protected: virtual int do_io( EMode mode ); }; class espace_db : public DatabaseIO { protected: virtual int do_io( EMode mode ); }; //_____________________________________________________________________________ int main(int argc, char **argv) { static const char rcsid[] = "$Id: Optimize.cpp,v 1.16 2004/02/25 22:16:46 ole Exp $"; TROOT minexam("optimize",rcsid); Double_t p0=0; Double_t p1=1; Double_t p3=3; double step; const double resolution=0.00005; TStopwatch timer; Int_t ierflg = 0; DatabaseIO* dbio; timer.Start(); getargs(argc,argv); cout << "fit index " << db.fit_indx << endl; // Initialize and read input parameters if(db.init() != 0) { cerr << "Error reading input file " << db.infil << endl; exit(1); } // Define database interface if( lESPACE ) dbio = new espace_db; else dbio = new cpp_db; // Read matrix elements if( dbio->read() != 0) { cerr << "Error reading input database file" << endl; exit(2); } // // copy coefficients of the spectrometer constants: // // FIXME: necessary? where do we really store what? for(int arm=1; arm<=2; arm++) { for(int i=0; i<4; i++) { db.sp[arm].Gamma[i] = sp[arm].gamma[i]; } db.sp[arm].lines = sp[arm].noptics; } // // In case of dp optimization calculate the expected // dp_kin values for the peaks here. They are constants // for each peak and need to be calculated only once // if( db.fit_indx == 2 ) db.dpk_peak(); // Print what we've read, and count nvary nvary = 0; for(int arm=1; arm<=2; arm++) { printf("Arm %i\n",arm); for(int i=0; i 0); } } // // Except for the offset optimization case, the optimization can NOT // fix offsets in the target variables. Add an offset term to the // variables so that the matrix elements do not get distorted in // case the offsets are not perfect. This offset is not written into // the espace database. // // if(db.fit_indx != 0) nvary++; int ivary=0; // Initialize Minuit with a maximum of nvary TMinuit *gMinuit = new TMinuit(nvary); gMinuit->mninit(5,6,7); gMinuit->SetFCN(fcnk0); if(db.fit_indx==4) { gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg); goto L99; } for(int arm=1; arm<=2; arm++) { char fname[7]; if(arm==1) fname[0]='e'; else fname[0]='h'; for(int i=0; imnparm(ivary, fname, Opt.op_coeff[j],step, 0, 0,ierflg); ivary++; if (ierflg) { Printf(" UNABLE TO DEFINE PARAMETER NO. %d",ivary); return ierflg; } } } } if(db.fit_indx != 0) { gMinuit->mnparm(ivary,"offset", 0.002, 0.000001, 0, 0,ierflg); ivary++; if (ierflg) { Printf(" UNABLE TO DEFINE OFFSET PARAMETER"); return ierflg; } } gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg); // Quit here if test mode if(lTest) goto L99; Printf(" calling migrad."); gMinuit->mnexcm("MIGRAD", &p0 ,0,ierflg); gMinuit->mnexcm("CALL FCN", &p3 , 1,ierflg); // // Write the output database files // if( dbio->write() != 0) { cerr << "Error writing output database file" << endl; } else cout << "copy of output database successfully written" << endl; L99: delete gMinuit; delete dbio; cout<<"Time at the end of job = "<.dat // and db_R..dat, where is the value of // the "db_in" field in the optimization input file. // // In write mode, the files db_R/L..dat are copied to // db_R/L..dat, replacing the optics matrix elements // of the infiles with the optimized ones stored in memory. // is the value of the "db_out" field in the input file. // // The database files are expected in the basic ASCII file format // of the standard VDC class. // // Configuration tags are not supported at this time. If present, // they will be skipped. This routine always reads or writes the first // configuration. const int LEN = 200; const char** theSpectro = spectro_names; if(db.spectro<1||db.spectro>2){ cout<<"arm number error, 1 or 2"<2 ) break; const char* const* pelmt = elmt_names; int ielmt = 0; while( *pelmt && strcmp(w,*pelmt) ) { pelmt++; ielmt++; } if( !*pelmt ) break; // Here, we're not interested in the pathlength or forward elements if( (*pelmt)[0] == 'L' || len == 2 ) { if( writing ) fprintf(outfile,"%s",buff); continue; } // The good matrix elements have already been written at this point, // so quit here if( writing ) continue; // This looks like a good line, let's read it // Assume the following fixed format: // Optics& ropt = Spec.optics[Spec.noptics]; int* k = ropt.index; double* d = ropt.op_coeff; int n; if( (n = sscanf(buff, "%s %d %d %d %lf %lf %lf %lf %lf %lf %lf %d", w, k+1, k+2, k+3, d, d+1, d+2, d+3, d+4, d+5, d+6, k )) != 12 ) { // Fail on apparently incomplete lines if( n < 11 ) { cerr << "ERROR: not enough data for matrix element: " << endl << buff << "in database file " << filename << " " << "(need 12 fields = 4 codes + 7 coeffcients + varycode)" << endl; return 5; } // Hmmm ... only 11 elements - varycode is missing - old format // database. Print a more informative warning, then die cerr << "ERROR: varycode missing for matrix element: " << endl << buff << "in database file " << filename << " " << "(probably old database)" << endl << "Add varycode (int) right after coefficients." << endl; return 6; } strcpy(ropt.name, Form("%s%d%d%d",w,k[1],k[2],k[3])); Spec.noptics++; } // If writing then copy the rest of the file if( writing ) { while( !feof(file) && !ferror(file) && !ferror(outfile)) { fprintf(outfile,"%s",buff); fgets(buff,LEN,file); } } // Next file fclose(file); if( writing ) fclose(outfile); theSpectro+=3-db.spectro; //arm++; } return 0; } //_____________________________________________________________________________ int espace_db::do_io( EMode mode ) { // Read or write matrix elements from/to ESPACE database // Names of data sections const TString SPECR("specr"), SPECL("specl"), END("+++"); const int NR = 520, NL = 530; const int RIGHT = 1, LEFT = 2; const int LEN = 200; char buff[LEN], w[LEN]; bool writing = (mode == kWrite); // Open files FILE* file = fopen(db.db_in,"r"); if( !file ) return 1; FILE* outfile = 0; if( writing ) { outfile = fopen(db.db_out,"w"); if( !outfile ) { fclose(file); return 2; } } bool doing_matrix = false; int nlines = 0, arm = 0; // Read lines from the input file while( fgets(buff, LEN, file) ) { // When writing, copy lines to output except for the matrix elements if( writing && !doing_matrix ) fprintf(outfile,"%s",buff); // Data section? if( !arm && sscanf(buff,"%s",w) == 1 ) { if( SPECR == w ) arm = RIGHT; else if( SPECL == w ) arm = LEFT; } if( !arm ) continue; // Read fixed number of lines until we reach the matrix elements if( !doing_matrix ) { nlines++; if( (arm==RIGHT && nlines < NR) || (arm==LEFT && nlines < NL) ) continue; doing_matrix = true; continue; //get next line since this one already written in write mode } // Read gamma, skip other spectrometer parameters spec& Spec = sp[arm]; double* g = Spec.gamma; if( writing ) fprintf(outfile,"%10.4E %10.4E %10.4E %10.4E\n", g[0], g[1], g[2], g[3]); else sscanf(buff,"%lf %lf %lf %lf", g, g+1, g+2, g+3); fgets(buff,LEN,file); if( writing ) fprintf(outfile,"%s",buff); fgets(buff,LEN,file); if( writing ) fprintf(outfile,"%s",buff); // When writing, dump the data from memory if( writing ) { for(int i=0; i] " << "" << endl; cerr << " -t Test. Set everything up, but do not actually call MIGRAD." << endl; cerr << " -e Use ESPACE database format (default for now)." << endl; cerr << " -c Use C++ Analyzer database format." << endl; cerr << " -x Override x_tg from event data with a value calculated " << "from beam positions." << endl; cerr << " -P Phi angle optimization (default). Negates -T." << endl; cerr << " -T Theta angle optimization. Negates -P" << endl; cerr << " -B Optimize both Theta and Phi." << endl; cerr << " -A Optimize Ytg without sieve." < Use instead of opt_xxx.dat." << endl; cerr << " -f Use instead of the dat_file specified " << "in the input file." << endl; cerr << " Variable to optimize. One of: " << "offset, y_tg, dp, angles, emiss, test" < 1) { const char *opt = *++argv; if (*opt == '-') { while (*++opt != '\0') { switch (*opt) { case 't': lTest = true; break; case 'e': lESPACE = true; break; case 'c': lESPACE = false; break; case 'x': lCalcXtg = true; break; case 'T': lPhi = false; lTheta = true; break; case 'P': lPhi = true; lTheta = false; break; case 'B': lPhi = true; lTheta = true; break; case 'A': lSiv = false; break; case 'i': if (!*++opt) { if (argc-- < 1) usage(); opt = *++argv; } got_infile = true; db.infil = opt; opt = "?"; break; default: usage(); } } } else { // Only allow one type of optimization at a time if( gotvar ) usage(); // Parse argument not starting with '-' const Optvar* item = optvar; while( item->c ) { if( item->c == *opt ) { db.fit_indx = item->indx; if( !got_infile ) db.infil = item->infil; gotvar = true; break; } item++; } // Quit if unknown argument if( !gotvar ) break; } } if( !gotvar ) usage(); return; }