FluxAlg class Reference

00154 :Algorithm(name, pSvcLocator) , m_sequence(0), m_initialTime(0)
00155 , m_SAAreject(0), m_insideSAA(false)
00156 {
00157     
00158     // declare properties with setProperties calls
00159     declareProperty("source_name",  m_source_name="default");
00160     declareProperty("sources",     m_source_list);
00161     declareProperty("area",        m_area=6.0); // target area in m^2
00162     declareProperty("backoff",     m_backoff=2.0); //backoff distance in m
00163 
00164     declareProperty("rocking_angle", m_rocking_angle=0); // set non-zero to enable rocking
00165 
00166     declareProperty("PointingHistory",  m_pointingHistory); // doublet, filename and launch date
00167 
00168 // deprecate these
00169     declareProperty("pointing_info_tree_name",  m_root_tree="");
00170     declareProperty("save_pointing_info",  m_save_tuple=false);
00171 
00172     declareProperty("AvoidSAA",   m_avoidSAA=false);
00173     declareProperty("Prescale",   m_prescale=1);
00174     declareProperty("source_info", m_source_info_filename="source_info.txt");
00175     declareProperty("misalignment", m_misalignmentRotation);
00176     declareProperty("alignment", m_alignmentRotation);
00177     declareProperty("pointingDirection", m_pointingDirection);
00178     declareProperty("zenithTheta", m_zenithTheta=-99);
00179     declareProperty("FilterCone",  m_filterCone);
00180     declareProperty("sourceListFile", m_sourceListFile="");
00181     declareProperty("abortOnException", m_abort_on_exception=false);
00182 
00183 }

00360 {
00361     astro::EarthCoordinate pos = GPS::instance()->earthpos();
00362     m_insideSAA = pos.insideSAA();
00363 
00364     return 0; // can't be void in observer pattern
00365 }

00078 {return m_currentRate;}

00370 {
00371     StatusCode  sc = StatusCode::SUCCESS;
00372     MsgStream   log( msgSvc(), name() );
00373     //
00374     // Purpose: have the flux service create parameters of an incoming particle 
00375     // if nothing has changed, then use the existing m_flux,
00376     // but if the "current" IFlux is not the same as the one we have now,
00377     // then change our m_flux pointer to be the new one.
00378     // Output:  a staturCode to ensure the function executed properly.
00379     m_flux = m_fluxSvc->currentFlux();
00380 
00381     // check the current random number seed
00382     int seed = CLHEP::HepRandom::getTheSeed();
00383     log << MSG::DEBUG << "random seed: " << seed << endreq;
00384 
00385     std::string particleName;
00386     if( m_avoidSAA) m_SAAreject--;
00387     int count = m_prescale;
00388     do{ // loop if we are rejecting particles generated during SAA
00389         // also do prescale here
00390         try {
00391             bool valid =m_flux->generate();
00392             if( !valid) {
00393                 log << MSG::ERROR << "Ran out of valid sources, aborting" << endreq;
00394                 return StatusCode::FAILURE;
00395             }
00396         } catch( const std::exception & e) {
00397             std::cout << "FluxAlg caught exception, aborting this event " << e.what() 
00398                 << "\n\tprocessing source " << m_flux->particleName() 
00399                 << "\n\tcurrent time: " <<GPS::instance()->time() << std::endl;
00400             if( m_abort_on_exception ){
00401                 return StatusCode::FAILURE;
00402             }
00403             setFilterPassed( false); // should go to clocks.
00404             return StatusCode::SUCCESS;
00405             
00406         }
00407         particleName = m_flux->particleName();
00408 
00409         //if it's a clock then ExposureAlg will take care of it, and no othe algorithms should care about it.
00410         if(particleName == "TimeTick" || particleName == "Clock"){
00411             m_pointing_info.set();
00412 
00413             setFilterPassed( false );
00414             return sc;
00415         }
00416         if(m_insideSAA && m_avoidSAA.value() ){
00417             double time(GPS::instance()->time()), 
00418                 endtime( m_fluxSvc->endruntime() );
00419             if( time >endtime ){
00420                 log << MSG::INFO << "Ran out of time while in SAA"<< endreq;
00421                 setFilterPassed( false );
00422                 break;  //return sc;
00423             }
00424                 
00425             ++m_SAAreject;
00426         }
00427         else break;
00428     } while(m_insideSAA && m_avoidSAA.value() || --count>0);
00429 
00430     Hep3Vector p = m_flux->launchPoint();
00431     Hep3Vector d = m_flux->launchDir();
00432 
00433     double ke = m_flux->energy(); // kinetic energy in MeV
00434 
00435     //here's where we get the particleID and mass for later.
00436     // Note that the Gaudi particle table now only has p+ and n0 for proton and neutron: 
00437     if( particleName=="p" || particleName=="proton") particleName="p+";
00438     if( particleName=="neutron") particleName="n0";
00439     ParticleProperty* prop = m_partSvc->find(particleName);
00440 
00441 //    if( prop==0 && particleName=="He" ){
00442 //        // If He didn't work (mystery!) try alpha instead
00443 //        prop = m_partSvc->find("alpha");
00444 //    }
00445     if( prop==0) {
00446         log << MSG::ERROR << "Particle name " << particleName << " not found by particle properties" << endreq;
00447         return StatusCode::FAILURE;
00448     }
00449 
00450     int partID = prop->jetsetID(); // same as stdhep id
00451 
00452     log << MSG::DEBUG ;
00453     if( log.isActive()){
00454         log<< particleName << ", flux("<<m_flux->name() << ") "
00455         << "(" << m_flux->energy()
00456         << " MeV), Launch: " 
00457         << "(" << p.x() <<", "<< p.y() <<", "<<p.z()<<")" 
00458         << " mm, Dir " 
00459         << "(" << d.x() <<", "<< d.y() <<", "<<d.z()<<")" ;
00460     }
00461     log << endreq;
00462 
00463 
00464     // Here the TDS is prepared to receive hits vectors
00465     // Check for the MC branch - it will be created if it is not available
00466     Event::MCEvent* mch = 0;
00467 
00468     SmartDataPtr<Event::MCEvent> mcheader(eventSvc(), EventModel::MC::Event);
00469     if (mcheader == 0) {
00470         sc=eventSvc()->registerObject(EventModel::MC::Event , mch= new Event::MCEvent);
00471         mch->initialize(0,0,m_sequence, m_flux->time(), m_flux->name());
00472         if(sc.isFailure()) {
00473             log << MSG::WARNING << EventModel::MC::Event  <<" could not be registered on data store" << endreq;
00474             delete mch;
00475             return sc;
00476         }
00477 
00478     }else {
00479         mch = mcheader;
00480     }
00481 
00482     mch->initialize(mch->getRunNumber(), m_flux->numSource(), mch->getSequence(), m_flux->time());
00483     mch->setSourceName(m_flux->name());
00484 
00485     Event::McParticleCol* pcol = new Event::McParticleCol;
00486     sc = eventSvc()->registerObject(EventModel::MC::McParticleCol, pcol);
00487     if( sc.isFailure()) {
00488 
00489         log << MSG::ERROR << "Could not Register "<< EventModel::MC::McParticleCol << endreq;
00490 
00491         return sc;
00492     }
00493     Event::McParticle * parent= new Event::McParticle;
00494     pcol->push_back(parent);
00495 
00496     double mass = prop->mass() , 
00497         energy = (ke+mass),
00498         momentum=sqrt(energy*energy - mass*mass); 
00499     CLHEP::HepLorentzVector pin(d*momentum,energy);
00500 
00501     // This parent particle decay at the start in the first particle, 
00502     // so initial momentum and final one are the same
00503     parent->initialize(parent, partID, 
00504         Event::McParticle::PRIMARY,
00505         pin,p, m_flux->name());
00506     parent->finalize(pin, p);
00507 
00508     // get the event header to set the time
00509     Event::EventHeader* h = 0; 
00510 
00511     SmartDataPtr<Event::EventHeader> header(eventSvc(), EventModel::EventHeader);
00512     if(0==header) {
00513         // not already there: try to register instead
00514         //sc = eventSvc()->registerObject(EventModel::EventHeader, h=new Event::EventHeader);
00515         sc = eventSvc()->registerObject(EventModel::EventHeader, EventModel::EventHeader, h=new Event::EventHeader);
00516         if( sc.isFailure()) {
00517             log << MSG::WARNING << " could not find or register the event header" << endreq;
00518         }
00519     } else{  h = header;
00520     }
00521 
00522     m_currentTime=m_flux->time();
00523 
00524     // is this proper here?
00525     m_pointing_info.set();
00526     
00527     // put pointing stuff into the root tree
00528     if( m_rootTupleSvc!=0 && !m_root_tree.value().empty()){
00529         m_rootTupleSvc->storeRowFlag(this->m_root_tree.value(), m_save_tuple);
00530     }
00531 
00532     if( m_initialTime==0) m_initialTime=m_currentTime;
00533     h->setTime(m_currentTime);
00534     m_run = h->run();  // save
00535     int numEvents = ++m_sequence;
00536     m_counts[m_flux->numSource()]++; // update count
00537     m_flux_names[m_flux->numSource()]= m_flux->name(); // save (or resave!) name 
00538 
00539     m_currentRate=numEvents/(m_currentTime-m_initialTime);
00540     return StatusCode::SUCCESS;
00541 }

00550                             {
00551     StatusCode  sc = StatusCode::SUCCESS;
00552     static bool done = false;
00553     if( done  ) return sc;
00554     done=true;
00555 
00556     if( m_rootTupleSvc!=0 ){
00557         // create the jobinfo tuple: copy to statics
00558         run = m_run;
00559         sequence=m_sequence + (m_SAAreject>0? m_SAAreject: 0 );
00560         initialTime=m_initialTime;
00561         currentTime=m_currentTime;
00562         m_rootTupleSvc->addItem("jobinfo", "run", &run);
00563         m_rootTupleSvc->addItem("jobinfo", "generated", &sequence);
00564         m_rootTupleSvc->addItem("jobinfo", "start", &initialTime);
00565         m_rootTupleSvc->addItem("jobinfo", "stop",  &currentTime);
00566     }
00567     MsgStream log(msgSvc(), name());
00568     log << MSG::INFO << "Computed Rate: "<< currentRate() << " Hz" ;
00569     summary(log.stream(), "\n\t\t\t");
00570     log  << endreq;
00571 
00572     if( !m_source_info_filename.value().empty() ){
00573         std::ofstream infofile(m_source_info_filename.value().c_str());
00574         summary( infofile, std::string("\n") );
00575     }
00576 
00577     
00578     if( m_avoidSAA && m_SAAreject>0 ){
00579         log << "\t\tRejected by SAA: " << m_SAAreject << endreq;
00580             log << "\t\t(note that this may invalidate the rate calculation)" << endreq;
00581     }
00582     return sc;
00583 }

00186                               {
00187     StatusCode  sc = StatusCode::SUCCESS;
00188     MsgStream log(msgSvc(), name());
00189     GPS* gps = GPS::instance();
00190 
00191     // Use the Job options service to set the Algorithm's parameters
00192     setProperties();
00193     // set target area for random point generation, and the backoff distance
00194     EventSource::totalArea(m_area);
00195     EventSource::s_backoff =m_backoff*1e3;// convert from m to mm
00196 
00197     // set pointing mode and associated rocking angle 
00198     if ( service("FluxSvc", m_fluxSvc).isFailure() ){
00199         log << MSG::ERROR << "Couldn't find the FluxSvc!" << endreq;
00200         return StatusCode::FAILURE;
00201     }
00202 
00203     if( m_pointingDirection.value().size()==2 ) {
00204         // if direction set, ignore everything else!
00205         double ra(m_pointingDirection.value()[0]), dec(m_pointingDirection.value()[1]);
00206         m_fluxSvc->setPointingDirection( astro::SkyDir(ra, dec));
00207         log << MSG::INFO << "set to point at ra,dec= " << ra << ", "<<dec << endreq;
00208 
00209     }else if( m_zenithTheta>=0) {
00210         // if set from defalt, set to this value
00211         m_fluxSvc->setRockType(astro::GPS::EXPLICIT, m_zenithTheta);
00212         log << MSG::INFO << "set to zenith angle " << m_zenithTheta << " degrees" << endreq;
00213 
00214     }else {
00215         if(m_rocking_angle >0 ){
00216             //output to record the pointing settings
00217             //then this line sets the rocking type, as well as the rocking angle.
00218             m_fluxSvc->setRockType(GPS::ONEPERORBIT ,m_rocking_angle);
00219             log << MSG::INFO << "Once per orbit rocking Angle: " << m_rocking_angle << " degrees" << endreq;
00220         }
00221         
00222     
00223         //set the input file to be used as the pointing database, if used
00224         if(! m_pointingHistory.value().empty()){
00225             std::string filename(m_pointingHistory.value()[0]);
00226             facilities::Util::expandEnvVar(&filename);
00227             double offset = 0;
00228             bool horizontalflag(false);
00229             if( m_pointingHistory.value().size()>1){
00230                 std::string field(m_pointingHistory.value()[1]);
00231                 if(! field.empty() ) { // allow null string
00232                     facilities::Timestamp jt(m_pointingHistory.value()[1]);
00233                     offset = (astro::JulianDate(jt.getJulian())-astro::JulianDate::missionStart())*astro::JulianDate::secondsPerDay;
00234                 }
00235             }
00236 
00237             if( m_pointingHistory.value().size()>2){
00238                 std::string field(m_pointingHistory.value()[2]);
00239                 horizontalflag =! field.empty();
00240             }
00241             log << MSG::INFO << "Loading Pointing History File : " << filename 
00242                 << " with MET offset "<< offset <<  endreq;
00243             if( horizontalflag){
00244                 log << MSG::INFO << "Will override x-direction to be horizontal"<<endreq;
00245             }
00246 
00247             gps->setPointingHistoryFile(filename, offset, horizontalflag);
00248         }
00249     }
00250     double current_time = gps->time(); // preserve time to protect against Pulsar, etc.
00251 
00252     // -------------- source name processing ------------------
00253     std::vector<std::string> sources(m_source_list.value()); // list from job options
00254 
00255     // parse file with source library entries (consistent with obssim)
00256     if( !m_sourceListFile.value().empty() ) {
00257 
00258         std::string fname(m_sourceListFile.value().c_str() );
00259         facilities::Util::expandEnvVar(&fname);
00260         std::ifstream file(fname.c_str());
00261         if( !file.is_open() ){
00262             throw std::invalid_argument("File not found: " + fname);
00263         }
00264 
00265         while( ! file.eof()){
00266             std::string line; std::getline(file, line);
00267             if( line.empty() || line[0]=='#' ) continue; 
00268             sources.push_back(line);
00269         }
00270     }
00271 
00272     // now, are there any in the list?
00273 
00274     if( !sources.empty()){
00275         log << MSG::INFO << "loading sources " << endreq;
00276         for(std::vector<std::string>::const_iterator it= sources.begin(); it!=sources.end(); ++it){
00277             log << MSG::INFO << "\t" << (*it) << endreq;
00278         }
00279         sc =  m_fluxSvc->compositeSource(sources, m_flux);
00280         if( sc.isFailure()) {
00281             log << MSG::ERROR << "Could not find one of the sources" << endreq;
00282             return sc;
00283         }
00284 
00285 
00286     }else{
00287         // no, try a single source for compatibility
00288 
00289         log << MSG::INFO << "loading source " << m_source_name << endreq;
00290 
00291         sc =  m_fluxSvc->source(m_source_name, m_flux);
00292         if( sc.isFailure()) {
00293             log << MSG::ERROR << "Could not find flux " << m_source_name << endreq;
00294             return sc;
00295         }
00296     }
00297 
00298     gps->time(current_time);  // restore time if it was modified
00299     gps->synch();             // and first notification of attached observers
00300     std::string title(m_flux->title()); if(title.length()>100) title = title.substr(0,100)+"...";
00301     log << MSG::INFO << "Source title: " << title << endreq;
00302     log << MSG::INFO << "        area: " << m_flux->targetArea() << endreq;
00303     log << MSG::INFO << "        rate: " << m_flux->rate() << endreq;
00304     if( m_prescale>1) {
00305         log << MSG::INFO << "    prescale: "<< m_prescale << endreq;
00306     }
00307 
00308     // check for (mis) alignment
00309     int alignment_pars ( m_alignmentRotation.value().size() );
00310     if (alignment_pars >0){
00311 
00312         double qx(m_alignmentRotation.value()[0]),qy(0),qz(0);
00313         if(alignment_pars >1) qy = m_alignmentRotation.value()[1];
00314         if(alignment_pars >2) qz = m_alignmentRotation.value()[2];
00315         m_fluxSvc->setAlignmentRotation(qx, qy, qz, false);
00316     }
00317 
00318     alignment_pars = ( m_misalignmentRotation.value().size() );
00319     if (alignment_pars >0){
00320 
00321         double qx(m_misalignmentRotation.value()[0]),qy(0),qz(0);
00322         if(alignment_pars >1) qy = m_misalignmentRotation.value()[1];
00323         if(alignment_pars >2) qz = m_misalignmentRotation.value()[2];
00324         m_fluxSvc->setAlignmentRotation(qx, qy, qz, true);
00325     }
00326 
00327     // check for filter cone
00328     if( m_filterCone.value().size()==3) {
00329 
00330         m_fluxSvc->setFilterCone(m_filterCone);
00331         
00332     }
00333 
00334     if ( service("ParticlePropertySvc", m_partSvc).isFailure() ){
00335         log << MSG::ERROR << "Couldn't find the ParticlePropertySvc!" << endreq;
00336         return StatusCode::FAILURE;
00337     }
00338 
00339     // get a pointer to RootTupleSvc, use only if available 
00340     if( (service("RootTupleSvc", m_rootTupleSvc, true) ). isFailure() ) {
00341         log << MSG::WARNING << " RootTupleSvc is not available, will not write Pt tuple" << endreq;
00342         m_rootTupleSvc=0;
00343     }else if( !m_root_tree.value().empty() ) {
00344         
00345         m_pointing_info.setPtTuple(m_rootTupleSvc, m_root_tree.value());
00346     }
00347 
00348 
00349     // attach an observer to be notified when orbital position changes
00350         // set callback to be notified when the position changes
00351     m_observer.setAdapter( new ActionAdapter<FluxAlg>
00352         (this, &FluxAlg::askGPS) );
00353 
00354     m_fluxSvc->attachGpsObserver(&m_observer);
00355     gps->notifyObservers();
00356     return sc;
00357 }

00586 {
00587     log << indent << " Source ID   Source Name                  counts";
00588     for(std::map<int,int>::const_iterator im=m_counts.begin(); im !=m_counts.end(); ++im) {
00589         log << indent
00590             << std::setw(10) <<im->first 
00591             << "   "  << std::setw(25) << std::left<< m_flux_names[im->first]
00592             << std::setw(10)<< std::right << im->second;
00593     }
00594     log << std::endl;
00595 }


Public Member Functions
	FluxAlg (const std::string &name, ISvcLocator *pSvcLocator)
	ctor
double	currentRate ()
StatusCode	initialize ()
	set parameters and attach to various perhaps useful services.
StatusCode	execute ()
	process an event
StatusCode	finalize ()
	clean up, summarize
Private Member Functions
int	askGPS ()
void	summary (std::ostream &log, std::string indent)
	StringProperty (m_sourceListFile)
	file name for list of sources (obssim compatible)
Private Attributes
double	m_currentRate
StringProperty	m_source_name
	a single source
StringArrayProperty	m_source_list
	allow specification of a list of sources, which will be combined
IFluxSvc *	m_fluxSvc
IFlux *	m_flux
unsigned int	m_run
unsigned int	m_sequence
IDataProviderSvc *	m_eds
IParticlePropertySvc *	m_partSvc
DoubleProperty	m_area
DoubleProperty	m_rocking_angle
DoubleProperty	m_rocking_angle_z
std::map< int, int >	m_counts
double	m_initialTime
	for measuring the total number generated per code.
double	m_currentTime
int	m_SAAreject
PointingInfo	m_pointing_info
StringArrayProperty	m_pointingHistory
	history file name and launch date
StringProperty	m_root_tree
BooleanProperty	m_save_tuple
BooleanProperty	m_avoidSAA
INTupleWriterSvc *	m_rootTupleSvc
ObserverAdapter< FluxAlg >	m_observer
bool	m_insideSAA
IntegerProperty	m_prescale
StringProperty	m_source_info_filename
std::map< int, std::string >	m_flux_names
DoubleArrayProperty	m_misalignmentRotation
DoubleArrayProperty	m_alignmentRotation
DoubleArrayProperty	m_pointingDirection
	(ra, dec) for pointing
DoubleProperty	m_backoff
	backoff distance
DoubleProperty	m_zenithTheta
	set for zenith
DoubleArrayProperty	m_filterCone
	set parameters of a cone
BooleanProperty	m_abort_on_exception
	what to do if an exception

FluxAlg Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation