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