Athena-Atlfast: User Guide | ||
---|---|---|
Prev | Chapter 5. Writing an Analysis Algorithm | Next |
One of the most powerful aspects of the Athena-Atlfast program is the ability for the user to write their own analysis algorithms and add them into the code. This then means the simulation and analysis can be run together. The example described in this chapter is based on (but not exactly the same as) the standard Atlfast ExampleAnalysis.
The user will need to construct two files, a header (*.h) file and a source code (*.cxx) file. Thus the user should now have two files, MyAnalysis.h and MyAnalysis.cxx.
If the user has checked out the Atlfast code from CVS then the .h file needs to be put in the src/Atlfast/AtlfastCode/AtlfastCode directory and the .cxx file will need to be put in the src/Atlfast/AtlfastCode/src directory. Before any new files can be compiled the following line has to be added to the src/Atlfast/AtlfastCode/GNUmakefile.in file so that the compiler knows to include it:
src/MyAnalysis.cxx \
If the user is linking directly to the Atlfast libraries without checking out the code, then the process is a little more involved and these instructions should be followed.
Once the two files are in the correct places the actual algorithms can be constructed. The analysis algorithms header file should start with a #ifndef and end with a #endif. For example:
#ifndef __Atlfast_MyAnalysis_H
#define __Atlfast_MyAnalysis_H
.....
lots of code!
.....
#endif
The algorithm should be enclosed by namespace Atlfast braces and inherit from Algorithm, the standard Athena algorithm. The default Constructor and destructor that must be included are as follows:
MyAnalysis( const std::string& name, ISvcLocator* pSvcLocator ) ;
virtual ~MyAnalysis();
Athena algorithms have three default methods that must be included:
StatusCode initialize() ;
StatusCode execute() ;
StatusCode finalize() ;
The initialize() is run at the start of the job, the finalize() at the end, and the execute() is run once per event. There are a number of Athena files that must be included and also some Atlfast ones so that the code will compile. Once this is done the user will have the following Skeleton Analysis Header file. To go with this there is a simple Skeleton Analysis source file . Both these files are already included in the Atlfast release code as AnalysisSkeleton.h and AnalysisSkeleton.cxx and the GNUmakefile.in has got AnalysisSkeleton.cxx already included, so will not need to be changed. This means that if the user has checked the code out from CVS, they only needs to edit the two AnalysisSkeleton files to make their own analysis algorithm.
This first example of adding stuff to the analysis code is a simple method to print out information about certain particle types. ReconstructedParticles are stored within the TES as a collection object. Within the algorithm the location in the TES of the particular particles being analysed must be defined. This means that the header file must define a member variable for this location to be written to:
std::string m_particleLocation ;
A default value for this can be defined at the top of the header file before the Atlfast namespace declaration:
#define DEFAULT_particleLocation "/Event/AtlfastIsolatedElectrons"
The particle location can then be set in the constructor in the source file:
m_particleLocation = DEFAULT_particleLocation ;
If the user wishes to change the particleLocation via the StandardAtlfastOptions.txt so to code doesn't have to be recompiled each time the user wants to change the particle type, the following line can be added to the constructor:
declareProperty( "ParticleLocation", m_particleLocation ) ;
Then all the user has to do to change the particleLocation is add the following line to the StandardAtlfastOptions.txt:
MyAnalysis.ParticleLocation = "/Event/AtlfastIsolatedPhotons";
A list of particle types and their locations is given in Chapter 6: Atlfast Reconstructed Entities. Now the algorithm knows where to look for the particles, the user can print out information about them. In any method which the user wishes to print out information, a MsgStream object needs to be made. For example, it is always good to put a message in the initialize method to tell the user the object has initialised:
MsgStream log( messageService(), name() ) ;
log << MSG::DEBUG << "Initialising" << endreq;
It is in the execute method, which is run once every event, that the output method for particle information needs to be added. The ReconstructedParticles stored in a collection object are typedefined by ReconstructedParticleCollection. The user needs to make a pointer to this collection:
SmartDataPtr<ReconstructedParticleCollection>
           particles( this->eventDataService(), m_particleLocation ) ;
if( ! particles ) {
     log << MSG::DEBUG
           << "No ReconstructedParticles found in TES at "
           << m_particleLocation
           << endreq ;
}
Now the the particle collection can be iterated over and information about each particle printed out, so an iterator has to be declared:
ReconstructedParticleCollection::iterator particle ;
int noPart = 0;
for( particle = particles->begin();
       particle < particles->end();
       ++particle){
       double pT = (*particle)->pT() ;
       log << MSG::DEBUG << " Particle pT was " << pT << endreq;
       ++noPart;
}
log << MSG::DEBUG << "No. of Particles " << noPart << endreq;
This then is the first step completed in writing an analysis algorithm. The complete code for these can be found in MyAnalysis.h and MyAnalysis.cxx. In these examples the code for printing out information has been separated out into a method called dump()
So far the example has only shown how to print information about the particles to the screen. Information can also be put into HBOOK histogram plots. Firstly the name of the histogram file is set in the jobOptions.txt
HistogramPersistencySvc.OutputFile = "atlfast.hbook";
The header file must declare a histogram object as a member variable and also an int as the histogram identifier:
IHistogram1D* m_hist;
int m_histStart ;
Like the particleLocation, the histogram identifier can be set in the StandardAtlfastOptions file. A default value can be declared in the header file before the Atlfast namespace declaration:
#define DEFAULT_histStart 1000
Then this property needs to be declared in the constructor in the .cxx file:
m_histStart = DEFAULT_histStart ;
declareProperty( "HistStart", m_histStart ) ;
and finally, the m_hist has to be initialised as a null pointer in the constructor:
m_hist = NULL ;
Now the histogram has to be "booked" in the initialise() method of the algorithm. Here the number and range of bins in the histogram are set, along with the histograms title:
m_hist =
  histoSvc()->book("/stat/simple1d/",
                       m_histStart,
                       "pT",
                       160,
                       0.0,
                       160.0
                     );
Once the histogram is booked it can be filled once per event in the execute() method. In this example the pT is filled just after it was printed out to the screen using:
m_hist->fill(pT,1.0);
The second argument is the weight, in this case 1.0. The histogram does not have to be filled with the pT, it could be for example the particles Eta:
m_hist->fill((*particle)->eta(), 1.0);
A full list of all reconstructed entities methods can be found in Chapter 6: Atlfast Reconstructed Entities More than one histogram can be booked and filled in a single job, but each has to have its own unique histogram object and identifier. The user can for example declare only m_histStart and then for subsequent bookings of histograms could give them the identifier m_histStart++. For example:
m_hist_pt =
  histoSvc()->book("/stat/simple1d/",
                       m_histStart++,
                       "pT",
                       160,
                       0.0,
                       160.0
                     );
m_hist_eta =
  histoSvc()->book("/stat/simple1d/",
                       m_histStart++,
                       "Eta",
                       100,
                       -5.0,
                       5.0
                     );
Two dimensional histograms can also be filled by declaring an IHistogram2D* object. This is then booked in the same way as a 1d histogram, except both the x and y axis range and number of bins has to be declared.
m_hist_2d =
  histoSvc()->book("/stat/simple2d/",
                       m_histStart++,
                       "Eta vs pT",
                       100,
                       -5.0,
                       5.0,
                       160,
                       0.0,
                       160.0
                     );
The complete example programs can be found at MyAnalysis.h and MyAnalysis.cxx.
![]() | home | Next ![]() |