Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]   [ 15.*.* ]  
001 #include "AtlfastAlgs/BremMatrixManager.h"
002 #include "PathResolver/PathResolver.h"
003 #include <iostream>
004 
005 namespace Atlfast
006 {
007   using std::abs;
008   using std::ifstream;
009 
010 /** @brief Provides smear matrices corresponding to given track trajectories.
011    *
012    * Used by tracksmearer to provide smear matrices
013    * corresponding to given track trajectories. It reads a flat
014    * file containing smear matrix data and creates a BinData object for
015    * every eta/rT bin.
016    */
017 
018 
019   //-----------------------------------------------------------
020   // PUBLIC: Constructor
021   //-----------------------------------------------------------
022   BremMatrixManager::BremMatrixManager( string paramFile, int randSeed, MsgStream log ) : 
023        m_file( paramFile ), 
024        m_randSeed( randSeed )
025   {
026     m_log = &log;
027     *m_log << MSG::INFO 
028            << "Constructing BremMatrixManager with parameter file "
029            << m_file 
030            << endreq;
031     m_file = PathResolver::find_file( m_file, "DATAPATH" );
032     this->initialise();
033     *m_log << MSG::INFO << "Constructed BremMatrixManager" << endreq;
034   }
035   
036   
037   //-----------------------------------------------------------
038   // PUBLIC: Destructor
039   //-----------------------------------------------------------
040   BremMatrixManager::~BremMatrixManager()
041   {
042     map<BinID, BremBinData*>::iterator iter = m_binData.begin();
043     map<BinID, BremBinData*>::iterator end = m_binData.end();
044     for ( ; iter != end; ++iter )
045     {
046       delete ( iter->second ); 
047     }
048   }
049   
050   //------------------------------------------------------------
051   // PRIVATE: initialise : read file and construct data bins
052   //------------------------------------------------------------
053   void BremMatrixManager::initialise()
054   {
055     // open file
056     ifstream input;
057     input.open( m_file.c_str() );
058     
059     if (input) 
060     {
061       *m_log << MSG::INFO << "BremMatrixManager: File " << m_file << " open." << endreq;
062       
063       double pTMin, etaMin, etaMax, rtBoundary;
064       
065       // read some parameters
066       input >> pTMin;
067       input >> etaMin;
068       input >> etaMax;
069       input >> m_nEtaBins;
070       input >> m_nRTBins;
071       
072       // construct vector<binboundaries>
073       double etaStep = (etaMax - etaMin) / double(m_nEtaBins);
074       for ( int i = 0; i < m_nEtaBins; i++ )
075       { 
076         m_etaBoundaries.push_back( etaMin + etaStep/2. + double(i)*etaStep );
077       }
078       
079       for ( int i = 0; i < m_nRTBins + 1; i++ )
080       { 
081         input >> rtBoundary;
082         m_rTBoundaries.push_back(rtBoundary);
083       }
084       
085       // parameters are stored in rT bin blocks----
086       // with parameters in format:
087       // power function start points: C0(d0, phi0, q/pT), 
088       //                              C1(d0, phi0, q/pT), etc.
089       // power function slopes: C0(d0, phi0, q/pT), 
090       //                        C1(d0, phi0, q/pT), etc.
091       // Gaussian sigmas: C0(d0, phi0, q/pT), 
092       //                  C1(d0, phi0, q/pT), etc.
093       //
094       // coefficients C0,C1,C2,C3,C4 define pT-dependence of the respective quantities
095       // according to f(pT) = C0 + C1/sqrt(pT) + C2/pT + C3/pT/sqrt(pT) + C4/pT^2
096         
097       //start bin id ints at zero
098       int iBin = 0; 
099       
100       for ( int rt = 0; rt < m_nRTBins; rt++ )
101       {
102         // make vectors to hold all resolution parameters for this rT bin
103         vector<double> empty;
104         
105         vector< vector<double> >  startPointsC0( m_nEtaBins, empty ); 
106         vector< vector<double> >  startPointsC1( m_nEtaBins, empty );
107         vector< vector<double> >  startPointsC2( m_nEtaBins, empty );
108         vector< vector<double> >  startPointsC3( m_nEtaBins, empty );
109         vector< vector<double> >  startPointsC4( m_nEtaBins, empty );
110         
111         vector< vector<double> >  slopesC0( m_nEtaBins, empty );
112         vector< vector<double> >  slopesC1( m_nEtaBins, empty );
113         vector< vector<double> >  slopesC2( m_nEtaBins, empty );
114         vector< vector<double> >  slopesC3( m_nEtaBins, empty );
115         vector< vector<double> >  slopesC4( m_nEtaBins, empty );
116         
117         vector< vector<double> >  sigmasC0( m_nEtaBins, empty );
118         vector< vector<double> >  sigmasC1( m_nEtaBins, empty );
119         vector< vector<double> >  sigmasC2( m_nEtaBins, empty );
120         vector< vector<double> >  sigmasC3( m_nEtaBins, empty );
121         vector< vector<double> >  sigmasC4( m_nEtaBins, empty );
122         
123         //read eta bin number of values for each Ci, i=0,1,2,3, and parameter
124         // read start point data
125         fillVector( input, startPointsC0, 3 );
126         fillVector( input, startPointsC1, 3 );
127         fillVector( input, startPointsC2, 3 );
128         fillVector( input, startPointsC3, 3 );
129         fillVector( input, startPointsC4, 3 );
130 
131         // read slope data
132         fillVector( input, slopesC0, 3 );
133         fillVector( input, slopesC1, 3 );
134         fillVector( input, slopesC2, 3 );
135         fillVector( input, slopesC3, 3 );
136         fillVector( input, slopesC4, 3 );
137 
138         // read Gaussian sigma data
139         fillVector( input, sigmasC0, 3 );
140         fillVector( input, sigmasC1, 3 );
141         fillVector( input, sigmasC2, 3 );
142         fillVector( input, sigmasC3, 3 );
143         fillVector( input, sigmasC4, 3 );
144         
145         // DATA READING FINISHED FOR THIS RT BIN
146         
147         // got all values, now make rT/eta bins and resolution objects
148         for ( int i = 0; i < m_nEtaBins - 1; ++i ) 
149         {         
150           double etaLow = m_etaBoundaries[i];
151           double etaHigh = m_etaBoundaries[i+1];
152           
153           // make bin id with rT and eta boundaries
154           BinID rTEtaBin( iBin, m_rTBoundaries[rt], m_rTBoundaries[rt+1], etaLow, etaHigh );
155           
156           // make parameter resolutions for each parameter
157           vector<ParameterResolutions*> startPoints;
158           vector<ParameterResolutions*> slopes;
159           vector<ParameterResolutions*> sigmas;
160           for ( int param = 0; param < 3; param++ ) 
161           {
162             vector<BinID> startPointCoeffBins;
163             startPointCoeffBins.push_back( BinID( 0, startPointsC0[i][param], startPointsC0[i+1][param] ) );
164             startPointCoeffBins.push_back( BinID( 0, startPointsC1[i][param], startPointsC1[i+1][param] ) );
165             startPointCoeffBins.push_back( BinID( 0, startPointsC2[i][param], startPointsC2[i+1][param] ) );
166             startPointCoeffBins.push_back( BinID( 0, startPointsC3[i][param], startPointsC3[i+1][param] ) );
167             startPointCoeffBins.push_back( BinID( 0, startPointsC4[i][param], startPointsC4[i+1][param] ) );
168             startPoints.push_back( new ParameterResolutions( startPointCoeffBins, etaLow, etaHigh ) );
169 
170             vector<BinID> slopeCoeffBins;
171             slopeCoeffBins.push_back( BinID( 0, slopesC0[i][param], slopesC0[i+1][param] ) );
172             slopeCoeffBins.push_back( BinID( 0, slopesC1[i][param], slopesC1[i+1][param] ) );
173             slopeCoeffBins.push_back( BinID( 0, slopesC2[i][param], slopesC2[i+1][param] ) );
174             slopeCoeffBins.push_back( BinID( 0, slopesC3[i][param], slopesC3[i+1][param] ) );
175             slopeCoeffBins.push_back( BinID( 0, slopesC4[i][param], slopesC4[i+1][param] ) );
176             slopes.push_back( new ParameterResolutions( slopeCoeffBins, etaLow, etaHigh ) );
177 
178             vector<BinID> sigmaCoeffBins;
179             sigmaCoeffBins.push_back( BinID( 0, sigmasC0[i][param], sigmasC0[i+1][param] ) );
180             sigmaCoeffBins.push_back( BinID( 0, sigmasC1[i][param], sigmasC1[i+1][param] ) );
181             sigmaCoeffBins.push_back( BinID( 0, sigmasC2[i][param], sigmasC2[i+1][param] ) );
182             sigmaCoeffBins.push_back( BinID( 0, sigmasC3[i][param], sigmasC3[i+1][param] ) );
183             sigmaCoeffBins.push_back( BinID( 0, sigmasC4[i][param], sigmasC4[i+1][param] ) );
184             sigmas.push_back( new ParameterResolutions( sigmaCoeffBins, etaLow, etaHigh ) );
185           }
186           
187           // Now make a bin data with edges rT and eta, containing the 
188           // data for bremsstrahlung corrections
189           BremBinData* binData = new BremBinData( rTEtaBin, startPoints, slopes, sigmas, m_randSeed );
190           
191           // enter bin data into map
192           m_binData[rTEtaBin] = binData;
193           
194           // increment bin ID
195           iBin++;
196           
197         } // end of eta bin loop
198       } // end of rT bin loop
199 
200       // close file
201       input.close();
202       
203     }
204     else  // no input file 
205     { 
206       *m_log << MSG::INFO 
207              << "BremMatrixManager: no data file ( " << m_file << " )!!!!" 
208              << endreq;
209     }
210     
211     makeHeader();
212     
213   }
214   
215   void BremMatrixManager::makeHeader()
216   {
217     HeaderPrinter hp( "Atlfast BremMatrixManager:", *m_log );
218     hp.add( "Total number of Bins     ", m_nRTBins * m_nEtaBins );
219     hp.add( "rT Bin Min               ", m_rTBoundaries.front() );        
220     hp.add( "rT Bin Max               ", m_rTBoundaries.back() );
221     hp.add( "Number of rT Bins        ", m_nRTBins );
222     hp.add( "Eta Bin Min              ", m_etaBoundaries.front() );        
223     hp.add( "Eta Bin Max              ", m_etaBoundaries.back() );
224     hp.add( "Number of Eta Bins       ", m_nEtaBins );
225     hp.print();
226   }
227   
228   //-----------------------------------------------------------
229   // PRIVATE: fillVector
230   //             reads n x M parameters into member variable
231   //-----------------------------------------------------------
232   void BremMatrixManager::fillVector( ifstream& input, 
233                                       vector< vector<double> >& myVector, 
234                                       int M = 3 )
235   {
236     double res;
237     for ( int param = 0; param < M; param++ ) 
238     { 
239       vector< vector<double> >::iterator binIter = myVector.begin();
240       for ( ; binIter != myVector.end(); binIter++ )
241       {
242         input >> res;
243         binIter->push_back(res);
244       } 
245       
246     }
247   }
248   
249  
250   //-----------------------------------------------------------
251   // PUBLIC: getMatrix : returns the sigma matrix corresponding
252   //                     to a given track trajectory
253   //-----------------------------------------------------------
254   TrackTrajectory BremMatrixManager::getBremTrack( const TrackTrajectory& track ) const
255   {
256 
257     BremBinData* binData = getBinData(track);
258     
259     TrackTrajectory newTrack = binData->getBremTrack(track);
260 
261     return newTrack;
262   }
263   
264   
265   //-----------------------------------------------------------
266   // Private: getBinData : returns the BremBinData of bin corresponding
267   //                       to a given track trajectory
268   //-----------------------------------------------------------
269   BremBinData* BremMatrixManager::getBinData( const TrackTrajectory& traj ) const
270   {
271     TrackParameters track = traj.parameters();
272 
273     vector<double> rTEta;
274     double rT = abs( traj.radius() );
275     double eta = abs( track.eta() );
276 
277     // validity check
278     double rTLow = ( (m_binData.begin())->first ).low(0);
279     double rTHigh = ( (m_binData.rbegin())->first ).high(0);
280     double etaLow = ( (m_binData.begin())->first ).low(1);
281     double etaHigh = ( (m_binData.rbegin())->first ).high(1);
282 
283     if ( rT < rTLow )  rT = rTLow;
284     if ( rT > rTHigh ) rT = rTHigh;  
285     if ( eta < etaLow )  eta = etaLow;
286     if ( eta > etaHigh ) eta = etaHigh;  
287 
288     // find BinID
289     rTEta.push_back(rT);
290     rTEta.push_back(eta);
291 
292     map<BinID, BremBinData*>::const_iterator binIter = m_binData.begin();
293     map<BinID, BremBinData*>::const_iterator binEnd  = m_binData.end();
294 
295     for ( ; binIter != binEnd; ++binIter )
296     {
297       if ( binIter->first.isInBin(rTEta) ) 
298       {
299         return binIter->second;
300       }
301     }
302     // OOPS! couldn't fin bin
303     *m_log << MSG::WARNING 
304            << "WARNING: BremMatrixManager - No bin; rT " << rT << ", eta " << eta 
305            << endreq;
306     
307     return ( m_binData.begin() )->second;
308   }
309   
310   
311 } // namespace
312
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine. Valid HTML 4.01!