rootplot.cxx Source File

00001 
00002 #include "rootplot.h"
00003 
00004 
00005 rootplot::rootplot(std::vector<char*> argv): NUM_BINS(30),LOOP(30000),
00006 TIME(0.01), ENERGY_MIN(0.01*1000.), ENERGY_MAX(100.0*1000.)
00007 {
00008     
00009     int argc = argv.size();
00010     static const char * default_arg="default";
00011     static const char * default_graph="log";
00012     
00013     int num_bins = NUM_BINS;  // Initialize to default number of bins
00014     int loop = LOOP; // Initialize to default number of iterations
00015     int num_sources = 0;
00016     int num_longsources = 0;
00017     int current_arg = 0;
00018     int longtime=1;
00019     int longtimemax=20;
00020     double time=TIME;  //time to use for flux and rate functions
00021     double energy_min = ENERGY_MIN;
00022     double energy_max = ENERGY_MAX;
00023     bool use_trueflux = false;
00024     bool use_flux = false;
00025     bool use_flux_min = false;
00026     bool use_flux_max = false;
00027     bool write_to_file = false;
00028     bool sum_mode = false;
00029     bool longterm=false;
00030     double flux_min;
00031     double flux_max;
00032     std::string arg_name(default_arg);
00033     std::string output_file_name;
00034     std::string ylabel_flux   = "Flux (particles/m^2/s/MeV/sr)";
00035     std::string ylabel_eflux  = "E*Flux (particles/m^2/s/sr)";
00036     std::string ylabel_sflux  = "Flux (particles/m^2/s/MeV)"; // Flux integrated over all solid angles
00037     std::string ylabel_seflux = "E*Flux (particles/m^2/s)";   // E*Flux integrated over all solid angles
00038     
00039     std::vector<std::string> sources;
00040     std::vector<int> longsources;
00041     
00042     flux_load();
00043     
00044     FluxMgr fm(sources); 
00045     
00046     
00047     // Process Command Line Arguments
00048     
00049     std::cout << "------------------------------------------------------" << std::endl;
00050     std::cout << " Flux test program: type 'rootplot -help' for help" << std::endl;
00051     std::cout << ( ( argc == 1)?  " No command line args, using defaults"
00052         :  "") << std::endl;
00053     
00054     
00055     
00056     while(current_arg < argc)
00057     {
00058         arg_name = argv[current_arg];
00059         if("-help" == arg_name || "help" == arg_name) 
00060         { 
00061             help(); 
00062             return;// 0; 
00063         }
00064         else if("-bins" == arg_name) 
00065             num_bins = atoi(argv[++current_arg]);
00066         else if("-events" == arg_name) 
00067             loop = atoi(argv[++current_arg]);
00068         else if("-min" == arg_name || "-energy_min" == arg_name) 
00069             energy_min = atof(argv[++current_arg]);
00070         else if("-max" == arg_name || "-energy_max" == arg_name) 
00071             energy_max = atof(argv[++current_arg]);
00072         else if("-flux_min" == arg_name)
00073         {
00074             use_flux_min = true;
00075             flux_min = atof(argv[++current_arg]);
00076         }
00077         else if("-flux_max" == arg_name)
00078         {
00079             use_flux_max = true;
00080             flux_max = atof(argv[++current_arg]);
00081         }
00082         else if("-flux" == arg_name)
00083             use_flux = true;
00084         else if("-trueflux" == arg_name)
00085             use_trueflux = true;
00086         else if("-file" == arg_name)
00087         {
00088             write_to_file = true;
00089             output_file_name = argv[++current_arg];
00090         }
00091         else if("-sum" == arg_name)
00092             sum_mode = true;
00093         else if("-list" == arg_name) 
00094         { 
00095             listSources(fm.sourceList());
00096             listSpectra(); 
00097             return; 
00098         }
00099         else if("-graph" == arg_name) 
00100             default_graph = argv[++current_arg];
00101         else if("-liny" == arg_name) default_graph = "semilogx";
00102         else if("-longsrc" == arg_name) {
00103             sources.push_back(argv[++current_arg]);
00104             //put the number of this source into the list for reference
00105             longsources.push_back(num_sources);
00106             num_sources++;
00107         }
00108         else if("-time" == arg_name) {
00109             time = atof(argv[++current_arg]);
00110             std::cout<<" TIME = "<< time << std::endl;
00111         }
00112         else if('-' == arg_name[0]) {std::cerr << "Unrecognized option "<< arg_name << ", -help for help" << std::endl;}
00113         else
00114         {
00115             sources.push_back(arg_name);
00116             num_sources++;
00117         }
00118         current_arg++;
00119     }
00120     
00121     
00122     // Use default source if no source was specified in arguments
00123     if(0 == sources.size())
00124     {
00125         sources.push_back(default_arg);
00126         num_sources++;
00127     }
00128     
00129     rootEnergyHist energy_hist(num_bins,energy_min,energy_max);
00130     rootAngleHist angle_hist(num_bins);
00131     
00132     // Process all the sources
00133     for(int i = 0; i < num_sources; i++)
00134     {
00135         int j;
00136         
00137         //decide whether or not this run should be longterm
00138         longterm = false;
00139         std::vector<int>::iterator longiter;
00140         for( longiter=longsources.begin(); longiter!=longsources.end() ;longiter++){
00141             if(*longiter==i) longterm=true;
00142         }
00143         
00144         if(longterm)
00145             fm.setExpansion(-1.);
00146         
00147         
00148         if((false == sum_mode && false==longterm)||(true==longterm && (longtime==1))) 
00149         {
00150             // Reset for new source
00151             energy_hist.reset();
00152             angle_hist.reset();
00153         }
00154         
00155         EventSource *e = fm.source(sources[i]);
00156         
00157         if(longterm){
00158             fm.pass(2.);
00159             time+=2.;
00160         }
00161         
00162         std::pair<double,double> loc=fm.location();
00163         std::cout << loc.first << "   " << loc.second << std::endl;
00164         //        std::cout << "orbit angle=" << GPS::instance()-> << "orbit phase=" << << std::endl;
00165         
00166         if( 0==e ) {std::cerr << "Source \"" << sources[i] << "\" not found: -list for a list" << std::endl;
00167         return;}
00168         
00169         energy_hist.setGraphType(default_graph);
00170         energy_hist.setTitle( sources[i] );
00171         
00172         energy_hist.setXLabel("Kinetic Energy (GeV)");
00173         
00174         if(true == use_flux)
00175         {
00176             energy_hist.setFluxMode();
00177             if(false == use_trueflux)
00178                 energy_hist.setYLabel(ylabel_sflux);
00179             else
00180                 energy_hist.setYLabel(ylabel_flux);
00181         }
00182         else
00183         {
00184             if(false == use_trueflux)
00185                 energy_hist.setYLabel(ylabel_seflux);
00186             else
00187                 energy_hist.setYLabel(ylabel_eflux);
00188         }
00189         
00190         if(true == use_flux_min)
00191             energy_hist.setFluxMin(flux_min);
00192         
00193         if(true == use_flux_max)
00194             energy_hist.setFluxMax(flux_max);
00195         
00196         angle_hist.setGraphType(default_graph);
00197         angle_hist.setTitle( sources[i] );
00198         angle_hist.setPhiXLabel("Angle (degrees)");
00199         angle_hist.setPhiYLabel("Particles");
00200         angle_hist.setThetaXLabel("Cos(Theta)");
00201         angle_hist.setThetaYLabel("Particles");
00202         
00203         double scale_factor;
00204         
00205         if(false == use_trueflux)
00206         {
00207             double rate = e->rate(time)/e->totalArea();
00208             scale_factor = rate/loop*num_bins/log(energy_max/energy_min);
00209             std::cout << "RATE IS: " << rate << std::endl;
00210         }
00211         else
00212         {
00213             double flux = e->flux(time);
00214             scale_factor = flux/loop*num_bins/log(energy_max/energy_min);
00215             std::cout << "FLUX IS: " << 1.0*flux << std::endl;
00216         }
00217         
00218         for(j = 0; j < loop; j++) 
00219         {
00220             FluxSource *f = e->event(time);
00221             double energy = f->energy();
00222             Vector dir = f->launchDir();
00223             double cos_theta = dir.z();
00224             
00225             double phi = atan2(dir.y(),dir.x());
00226             if(phi < 0)
00227                 phi = 2*M_PI + phi;
00228             
00229             energy_hist.store(energy);
00230             angle_hist.storeTheta(cos_theta);
00231             angle_hist.storePhi(phi);
00232             
00233             if(j % 1000 == 0) std::cerr << "\r" << j << ": " << energy << "...";
00234         }
00235         
00236         std::cerr << "\n";
00237         // Scale factor must be applied before the draw member function can be called
00238         energy_hist.apply(scale_factor);
00239         angle_hist.apply(scale_factor);
00240         
00241         for(j = 0; j < num_bins; j++) 
00242         {
00243             std::cout << energy_min*pow(10,((j + 0.5)/num_bins) * energy_hist.retrieveRange() ) 
00244                 << "   \t" << energy_hist.retrieveFlux(j) << "\t";
00245             std::cout << std::endl;         
00246         }
00247         
00248         
00249         if(false == sum_mode && false==longterm)
00250         {
00251             angle_hist.draw(1,"begin",i,num_sources);
00252             energy_hist.draw(1,"end",i,num_sources);
00253             
00254             delete e;
00255             std::cout << "Normal method" << std::endl;        
00256         }
00257         else if(true == sum_mode && i == num_sources - 1)
00258         {
00259             angle_hist.draw(1,"begin",0,1);
00260             energy_hist.draw(1,"end",0,1);
00261             
00262             delete e;
00263             std::cout << "Sum Mode method" << std::endl;
00264         }
00265         else if(longterm==true && longtime<=longtimemax)
00266         {
00267             longtime++;
00268             i--;
00269             std::cout << "Longterm run number " << longtime << std::endl; 
00270         }
00271         else
00272         {       
00273             if(true == write_to_file)
00274             {
00275                 std::ofstream output_file;
00276                 output_file_name += ".txt";
00277                 output_file.open(output_file_name.c_str());
00278                 
00279                 energy_hist.writeFile(1./(longtime),output_file);
00280                 output_file.close();
00281             }
00282             angle_hist.draw(1./double(longtime),"begin",i,num_sources);
00283             energy_hist.draw(1./double(longtime),"end",i,num_sources);
00284             
00285             delete e;
00286             longtime=1; 
00287             
00288         }
00289         
00290    } // for all sources   
00291    
00292    //return 0;
00293 }