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/PionMatrixManager.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 Used by tracksmearer to provide smear matrices .
011    *
012    * corresponding to given track trajectories. It reads a flat
013    * file containing smear matrix data and creates a BinData object for
014    * every eta/rT bin.
015    */
016 
017 
018   //-----------------------------------------------------------
019   // PUBLIC: Constructor
020   //-----------------------------------------------------------
021   PionMatrixManager::PionMatrixManager( string paramFile, int randSeed, MsgStream log ) : 
022        m_file( paramFile ), 
023        m_randSeed( randSeed )
024   {
025     m_log = &log;
026     *m_log << MSG::INFO 
027            << "Constructing PionMatrixManager with parameter file "
028            << m_file 
029            << endreq;
030     m_file = PathResolver::find_file( m_file, "DATAPATH" );
031     m_correlatedData = new CorrelatedData(randSeed);
032     this->initialise();
033     *m_log << MSG::INFO << "Constructed PionMatrixManager" << endreq;
034   }
035   
036   
037   //-----------------------------------------------------------
038   // PUBLIC: Destructor
039   //-----------------------------------------------------------
040   PionMatrixManager::~PionMatrixManager()
041   {
042     delete m_correlatedData;
043     map<BinID, IBinData*>::iterator iter = m_binData.begin();
044     map<BinID, IBinData*>::iterator end = m_binData.end();
045     for ( ; iter != end; ++iter )
046     {
047       delete ( iter->second ); 
048     }
049   }
050   
051   //------------------------------------------------------------
052   // PRIVATE: initialise : read file and construct data bins
053   //------------------------------------------------------------
054   void PionMatrixManager::initialise()
055   {
056     // open file
057     ifstream input;
058     input.open( m_file.c_str() );
059     
060     if (input) 
061     {
062       *m_log << MSG::INFO << "PionMatrixManager: File " << m_file << " open." << endreq;
063       
064       double pTMin, etaMin, etaMax, rtBoundary;
065       
066       // read some parameters
067       input >> pTMin;
068       input >> etaMin;
069       input >> etaMax;
070       input >> m_nEtaBins;
071       input >> m_nRTBins;
072       
073       // construct vector<binboundaries>
074       double etaStep = (etaMax - etaMin) / double(m_nEtaBins);
075       for ( int i = 0; i < m_nEtaBins; i++ )
076       { 
077         m_etaBoundaries.push_back( etaMin + etaStep/2. + double(i)*etaStep );
078       }
079       
080       for ( int i = 0; i < m_nRTBins + 1; i++ )
081       { 
082         input >> rtBoundary;
083         m_rTBoundaries.push_back(rtBoundary);
084       }
085       
086       // parameters are stored in rT bin blocks----
087       // with parameters in format:
088       // core  resolutions: 
089       //               C0(d0, z0, phi0, cot(theta0), q/pT), 
090       //               C1(d0, z0, phi0, cot(theta0), q/pT), etc.
091       // tails resolutions: 
092       //               C0(d0, z0, phi0, cot(theta0), q/pT), 
093       //               C1(d0, z0, phi0, cot(theta0), q/pT), etc.
094       // core fractions: 
095       //               C0(d0, z0, phi0, cot(theta0), q/pT), 
096       //               C1(d0, z0, phi0, cot(theta0), q/pT), etc.
097       // correlation coefficients rho: 
098       //               C0(rho(d0,phi0), rho(d0,q/pT), rho(phi0,q/pT), rho(z0,cot(theta0))), 
099       //               C1(rho(d0,phi0), rho(d0,q/pT), rho(phi0,q/pT), rho(z0,cot(theta0))), etc.
100       //
101       // coefficients C0,C1,C2,C3,C4 define pT-dependence of the respective quantities
102       // according to f(pT) = C0 + C1/sqrt(pT) + C2/pT + C3/pT/sqrt(pT) + C4/pT^2
103         
104       //start bin id ints at zero
105       int iBin = 0; 
106       
107       for ( int rt = 0; rt < m_nRTBins; rt++ )
108       {
109         // make vectors to hold all resolution parameters for this rT bin
110         vector<double> empty;
111         
112         vector< vector<double> >  coreC0( m_nEtaBins, empty ); 
113         vector< vector<double> >  coreC1( m_nEtaBins, empty );
114         vector< vector<double> >  coreC2( m_nEtaBins, empty );
115         vector< vector<double> >  coreC3( m_nEtaBins, empty );
116         vector< vector<double> >  coreC4( m_nEtaBins, empty );
117         
118         vector< vector<double> >  tailsC0( m_nEtaBins, empty );
119         vector< vector<double> >  tailsC1( m_nEtaBins, empty );
120         vector< vector<double> >  tailsC2( m_nEtaBins, empty );
121         vector< vector<double> >  tailsC3( m_nEtaBins, empty );
122         vector< vector<double> >  tailsC4( m_nEtaBins, empty );
123         
124         vector< vector<double> >  fractionsC0( m_nEtaBins, empty );
125         vector< vector<double> >  fractionsC1( m_nEtaBins, empty );
126         vector< vector<double> >  fractionsC2( m_nEtaBins, empty );
127         vector< vector<double> >  fractionsC3( m_nEtaBins, empty );
128         vector< vector<double> >  fractionsC4( m_nEtaBins, empty );
129 
130         vector< vector<double> >  correlationsC0( m_nEtaBins, empty );
131         vector< vector<double> >  correlationsC1( m_nEtaBins, empty );
132         vector< vector<double> >  correlationsC2( m_nEtaBins, empty );
133         vector< vector<double> >  correlationsC3( m_nEtaBins, empty );
134         vector< vector<double> >  correlationsC4( m_nEtaBins, empty );
135         
136         // read eta bin number of values for each Ci, i=0,1,2,3,4, and parameter
137         // read core data
138         fillVector( input, coreC0, 5 );
139         fillVector( input, coreC1, 5 );
140         fillVector( input, coreC2, 5 );
141         fillVector( input, coreC3, 5 );
142         fillVector( input, coreC4, 5 );
143 
144         // read tail data
145         fillVector( input, tailsC0, 5 );
146         fillVector( input, tailsC1, 5 );
147         fillVector( input, tailsC2, 5 );
148         fillVector( input, tailsC3, 5 );
149         fillVector( input, tailsC4, 5 );
150 
151         // read fractions data
152         fillVector( input, fractionsC0, 5 );
153         fillVector( input, fractionsC1, 5 );
154         fillVector( input, fractionsC2, 5 );
155         fillVector( input, fractionsC3, 5 );
156         fillVector( input, fractionsC4, 5 );
157 
158         // read correlations data
159         fillVector( input, correlationsC0, 4 );
160         fillVector( input, correlationsC1, 4 );
161         fillVector( input, correlationsC2, 4 );
162         fillVector( input, correlationsC3, 4 );
163         fillVector( input, correlationsC4, 4 );
164         
165         // DATA READING FINISHED FOR THIS RT BIN
166         
167         // got all values, now make rT/eta bins and resolution objects
168         for ( int i = 0; i < m_nEtaBins - 1; ++i ) 
169         {         
170           double etaLow = m_etaBoundaries[i];
171           double etaHigh = m_etaBoundaries[i+1];
172           
173           // make bin id with rT and eta boundaries
174           BinID rTEtaBin( iBin, m_rTBoundaries[rt], m_rTBoundaries[rt+1], etaLow, etaHigh );
175           
176           // make parameter resolutions for each parameter
177           vector<ParameterResolutions*> core;
178           vector<ParameterResolutions*> tails;
179           vector<ParameterResolutions*> fractions;
180           vector<ParameterResolutions*> correlations;
181           for ( int param = 0; param < 5; param++ ) 
182           {
183             vector<BinID> coreCoeffBins;
184             coreCoeffBins.push_back( BinID( 0, coreC0[i][param], coreC0[i+1][param] ) );
185             coreCoeffBins.push_back( BinID( 0, coreC1[i][param], coreC1[i+1][param] ) );
186             coreCoeffBins.push_back( BinID( 0, coreC2[i][param], coreC2[i+1][param] ) );
187             coreCoeffBins.push_back( BinID( 0, coreC3[i][param], coreC3[i+1][param] ) );
188             coreCoeffBins.push_back( BinID( 0, coreC4[i][param], coreC4[i+1][param] ) );
189             core.push_back( new ParameterResolutions( coreCoeffBins, etaLow, etaHigh ) );
190 
191             vector<BinID> tailsCoeffBins;
192             tailsCoeffBins.push_back( BinID( 0, tailsC0[i][param], tailsC0[i+1][param] ) );
193             tailsCoeffBins.push_back( BinID( 0, tailsC1[i][param], tailsC1[i+1][param] ) );
194             tailsCoeffBins.push_back( BinID( 0, tailsC2[i][param], tailsC2[i+1][param] ) );
195             tailsCoeffBins.push_back( BinID( 0, tailsC3[i][param], tailsC3[i+1][param] ) );
196             tailsCoeffBins.push_back( BinID( 0, tailsC4[i][param], tailsC4[i+1][param] ) );
197             tails.push_back( new ParameterResolutions( tailsCoeffBins, etaLow, etaHigh ) );
198 
199             vector<BinID> fractionsCoeffBins;
200             fractionsCoeffBins.push_back( BinID( 0, fractionsC0[i][param], fractionsC0[i+1][param] ) );
201             fractionsCoeffBins.push_back( BinID( 0, fractionsC1[i][param], fractionsC1[i+1][param] ) );
202             fractionsCoeffBins.push_back( BinID( 0, fractionsC2[i][param], fractionsC2[i+1][param] ) );
203             fractionsCoeffBins.push_back( BinID( 0, fractionsC3[i][param], fractionsC3[i+1][param] ) );
204             fractionsCoeffBins.push_back( BinID( 0, fractionsC4[i][param], fractionsC4[i+1][param] ) );
205             fractions.push_back( new ParameterResolutions( fractionsCoeffBins, etaLow, etaHigh ) );
206           }
207           
208           for ( int param = 0; param < 4; param++ ) 
209           {
210             vector<BinID> correlationsCoeffBins;
211             correlationsCoeffBins.push_back( BinID( 0, correlationsC0[i][param], correlationsC0[i+1][param] ) );
212             correlationsCoeffBins.push_back( BinID( 0, correlationsC1[i][param], correlationsC1[i+1][param] ) );
213             correlationsCoeffBins.push_back( BinID( 0, correlationsC2[i][param], correlationsC2[i+1][param] ) );
214             correlationsCoeffBins.push_back( BinID( 0, correlationsC3[i][param], correlationsC3[i+1][param] ) );
215             correlationsCoeffBins.push_back( BinID( 0, correlationsC4[i][param], correlationsC4[i+1][param] ) );
216             correlations.push_back( new ParameterResolutions( correlationsCoeffBins, etaLow, etaHigh ) );
217           }
218           
219           // Now make a bin data with edges rT and eta, containing the 
220           // correlation data
221           PionBinData* binData = new PionBinData( rTEtaBin, core, tails, fractions, correlations, m_randSeed );
222           
223           // enter bin data into map
224           m_binData[rTEtaBin] = binData;
225           
226           // increment bin ID
227           iBin++;
228           
229         } // end of eta bin loop
230       } // end of rT bin loop
231 
232       // close file
233       input.close();
234       
235     }
236     else  // no input file 
237     { 
238       *m_log << MSG::INFO 
239              << "PionMatrixManager: no data file ( " << m_file << " )!!!!" 
240              << endreq;
241     }
242     
243     makeHeader();
244     
245   }
246   
247   void PionMatrixManager::makeHeader()
248   {
249     HeaderPrinter hp( "Atlfast PionTrack Smearer:", *m_log );
250     hp.add( "Total number of Bins     ", m_nRTBins * m_nEtaBins );
251     hp.add( "rT Bin Min               ", m_rTBoundaries.front() );        
252     hp.add( "rT Bin Max               ", m_rTBoundaries.back() );
253     hp.add( "Number of rT Bins        ", m_nRTBins );
254     hp.add( "Eta Bin Min              ", m_etaBoundaries.front() );        
255     hp.add( "Eta Bin Max              ", m_etaBoundaries.back() );
256     hp.add( "Number of Eta Bins       ", m_nEtaBins );
257     hp.print();
258   }
259   
260   //-----------------------------------------------------------
261   // PRIVATE: fillVector
262   //             reads n x M parameters into member variable
263   //-----------------------------------------------------------
264   void PionMatrixManager::fillVector( ifstream& input, 
265                                       vector< vector<double> >& myVector, 
266                                       int M = 5 )
267   {
268     double res;
269     for ( int param = 0; param < M; param++ ) 
270     { 
271       vector< vector<double> >::iterator binIter = myVector.begin();
272       for ( ; binIter != myVector.end(); binIter++ )
273       {
274         input >> res;
275         binIter->push_back(res);
276       } 
277       
278     }
279   }
280   
281  
282   //-----------------------------------------------------------
283   // PUBLIC: getMatrix : returns the sigma matrix corresponding
284   //                     to a given track trajectory
285   //-----------------------------------------------------------
286   vector<double> PionMatrixManager::getVariables( const TrackTrajectory& track, 
287                                                   HepSymMatrix& returnSigma ) const
288   {
289     HepSymMatrix sigma;
290     vector<double> variables;
291 
292     IBinData* binData = getBinData(track);
293 
294     sigma = binData->getMatrix(track);
295     
296     // do Dcorset and Dcorgen to get smear parameters
297     variables = m_correlatedData->generate( m_correlatedData->root(sigma) );
298     returnSigma = sigma;
299     return variables;
300   }
301   
302   //-----------------------------------------------------------
303   // Private: getBinData : returns the IBinData of bin corresponding
304   //                     to a given track trajectory
305   //-----------------------------------------------------------
306   IBinData* PionMatrixManager::getBinData( const TrackTrajectory& traj ) const
307   {
308     TrackParameters track = traj.parameters();
309 
310     vector<double> rTEta;
311     double rT = abs( traj.radius() );
312     double eta = abs( track.eta() );
313 
314     // validity check
315     double rTLow = ( (m_binData.begin())->first ).low(0);
316     double rTHigh = ( (m_binData.rbegin())->first ).high(0);
317     double etaLow = ( (m_binData.begin())->first ).low(1);
318     double etaHigh = ( (m_binData.rbegin())->first ).high(1);
319 
320     if ( rT < rTLow )  rT = rTLow;
321     if ( rT > rTHigh ) rT = rTHigh;  
322     if ( eta < etaLow )  eta = etaLow;
323     if ( eta > etaHigh ) eta = etaHigh;  
324 
325     // find BinID
326     rTEta.push_back(rT);
327     rTEta.push_back(eta);
328 
329     map<BinID, IBinData*>::const_iterator binIter = m_binData.begin();
330     map<BinID, IBinData*>::const_iterator binEnd  = m_binData.end();
331 
332     for ( ; binIter != binEnd; ++binIter )
333     {
334       if ( binIter->first.isInBin(rTEta) ) 
335       {
336         return binIter->second;
337       }
338     }
339     // OOPS! couldn't fin bin
340     *m_log << MSG::WARNING 
341            << "WARNING: PionMatrixManager - No bin; rT " << rT << ", eta " << eta 
342            << endreq;
343     
344     return ( m_binData.begin() )->second;
345   }
346   
347   
348 } // namespace
349
[ 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!