00001 #include "AtlfastAlgs/PionMatrixManager.h"
00002 #include "PathResolver/PathResolver.h"
00003 #include <iostream>
00004
00005 namespace Atlfast
00006 {
00007 using std::abs;
00008 using std::ifstream;
00009
00010
00011
00012
00013 PionMatrixManager::PionMatrixManager( string paramFile, int randSeed, MsgStream log ) :
00014 m_file( paramFile ),
00015 m_randSeed( randSeed )
00016 {
00017 m_log = &log;
00018 *m_log << MSG::INFO
00019 << "Constructing PionMatrixManager with parameter file "
00020 << m_file
00021 << endreq;
00022 m_file = PathResolver::find_file( m_file, "DATAPATH" );
00023 m_correlatedData = new CorrelatedData(randSeed);
00024 this->initialise();
00025 *m_log << MSG::INFO << "Constructed PionMatrixManager" << endreq;
00026 }
00027
00028
00029
00030
00031
00032 PionMatrixManager::~PionMatrixManager()
00033 {
00034 delete m_correlatedData;
00035 map<BinID, IBinData*>::iterator iter = m_binData.begin();
00036 map<BinID, IBinData*>::iterator end = m_binData.end();
00037 for ( ; iter != end; ++iter )
00038 {
00039 delete ( iter->second );
00040 }
00041 }
00042
00043
00044
00045
00046 void PionMatrixManager::initialise()
00047 {
00048
00049 ifstream input;
00050 input.open( m_file.c_str() );
00051
00052 if (input)
00053 {
00054 *m_log << MSG::INFO << "PionMatrixManager: File " << m_file << " open." << endreq;
00055
00056 double pTMin, etaMin, etaMax, rtBoundary;
00057
00058
00059 input >> pTMin;
00060 input >> etaMin;
00061 input >> etaMax;
00062 input >> m_nEtaBins;
00063 input >> m_nRTBins;
00064
00065
00066 double etaStep = (etaMax - etaMin) / double(m_nEtaBins);
00067 for ( int i = 0; i < m_nEtaBins; i++ )
00068 {
00069 m_etaBoundaries.push_back( etaMin + etaStep/2. + double(i)*etaStep );
00070 }
00071
00072 for ( int i = 0; i < m_nRTBins + 1; i++ )
00073 {
00074 input >> rtBoundary;
00075 m_rTBoundaries.push_back(rtBoundary);
00076 }
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 int iBin = 0;
00098
00099 for ( int rt = 0; rt < m_nRTBins; rt++ )
00100 {
00101
00102 vector<double> empty;
00103
00104 vector< vector<double> > coreC0( m_nEtaBins, empty );
00105 vector< vector<double> > coreC1( m_nEtaBins, empty );
00106 vector< vector<double> > coreC2( m_nEtaBins, empty );
00107 vector< vector<double> > coreC3( m_nEtaBins, empty );
00108 vector< vector<double> > coreC4( m_nEtaBins, empty );
00109
00110 vector< vector<double> > tailsC0( m_nEtaBins, empty );
00111 vector< vector<double> > tailsC1( m_nEtaBins, empty );
00112 vector< vector<double> > tailsC2( m_nEtaBins, empty );
00113 vector< vector<double> > tailsC3( m_nEtaBins, empty );
00114 vector< vector<double> > tailsC4( m_nEtaBins, empty );
00115
00116 vector< vector<double> > fractionsC0( m_nEtaBins, empty );
00117 vector< vector<double> > fractionsC1( m_nEtaBins, empty );
00118 vector< vector<double> > fractionsC2( m_nEtaBins, empty );
00119 vector< vector<double> > fractionsC3( m_nEtaBins, empty );
00120 vector< vector<double> > fractionsC4( m_nEtaBins, empty );
00121
00122 vector< vector<double> > correlationsC0( m_nEtaBins, empty );
00123 vector< vector<double> > correlationsC1( m_nEtaBins, empty );
00124 vector< vector<double> > correlationsC2( m_nEtaBins, empty );
00125 vector< vector<double> > correlationsC3( m_nEtaBins, empty );
00126 vector< vector<double> > correlationsC4( m_nEtaBins, empty );
00127
00128
00129
00130 fillVector( input, coreC0, 5 );
00131 fillVector( input, coreC1, 5 );
00132 fillVector( input, coreC2, 5 );
00133 fillVector( input, coreC3, 5 );
00134 fillVector( input, coreC4, 5 );
00135
00136
00137 fillVector( input, tailsC0, 5 );
00138 fillVector( input, tailsC1, 5 );
00139 fillVector( input, tailsC2, 5 );
00140 fillVector( input, tailsC3, 5 );
00141 fillVector( input, tailsC4, 5 );
00142
00143
00144 fillVector( input, fractionsC0, 5 );
00145 fillVector( input, fractionsC1, 5 );
00146 fillVector( input, fractionsC2, 5 );
00147 fillVector( input, fractionsC3, 5 );
00148 fillVector( input, fractionsC4, 5 );
00149
00150
00151 fillVector( input, correlationsC0, 4 );
00152 fillVector( input, correlationsC1, 4 );
00153 fillVector( input, correlationsC2, 4 );
00154 fillVector( input, correlationsC3, 4 );
00155 fillVector( input, correlationsC4, 4 );
00156
00157
00158
00159
00160 for ( int i = 0; i < m_nEtaBins - 1; ++i )
00161 {
00162 double etaLow = m_etaBoundaries[i];
00163 double etaHigh = m_etaBoundaries[i+1];
00164
00165
00166 BinID rTEtaBin( iBin, m_rTBoundaries[rt], m_rTBoundaries[rt+1], etaLow, etaHigh );
00167
00168
00169 vector<ParameterResolutions*> core;
00170 vector<ParameterResolutions*> tails;
00171 vector<ParameterResolutions*> fractions;
00172 vector<ParameterResolutions*> correlations;
00173 for ( int param = 0; param < 5; param++ )
00174 {
00175 vector<BinID> coreCoeffBins;
00176 coreCoeffBins.push_back( BinID( 0, coreC0[i][param], coreC0[i+1][param] ) );
00177 coreCoeffBins.push_back( BinID( 0, coreC1[i][param], coreC1[i+1][param] ) );
00178 coreCoeffBins.push_back( BinID( 0, coreC2[i][param], coreC2[i+1][param] ) );
00179 coreCoeffBins.push_back( BinID( 0, coreC3[i][param], coreC3[i+1][param] ) );
00180 coreCoeffBins.push_back( BinID( 0, coreC4[i][param], coreC4[i+1][param] ) );
00181 core.push_back( new ParameterResolutions( coreCoeffBins, etaLow, etaHigh ) );
00182
00183 vector<BinID> tailsCoeffBins;
00184 tailsCoeffBins.push_back( BinID( 0, tailsC0[i][param], tailsC0[i+1][param] ) );
00185 tailsCoeffBins.push_back( BinID( 0, tailsC1[i][param], tailsC1[i+1][param] ) );
00186 tailsCoeffBins.push_back( BinID( 0, tailsC2[i][param], tailsC2[i+1][param] ) );
00187 tailsCoeffBins.push_back( BinID( 0, tailsC3[i][param], tailsC3[i+1][param] ) );
00188 tailsCoeffBins.push_back( BinID( 0, tailsC4[i][param], tailsC4[i+1][param] ) );
00189 tails.push_back( new ParameterResolutions( tailsCoeffBins, etaLow, etaHigh ) );
00190
00191 vector<BinID> fractionsCoeffBins;
00192 fractionsCoeffBins.push_back( BinID( 0, fractionsC0[i][param], fractionsC0[i+1][param] ) );
00193 fractionsCoeffBins.push_back( BinID( 0, fractionsC1[i][param], fractionsC1[i+1][param] ) );
00194 fractionsCoeffBins.push_back( BinID( 0, fractionsC2[i][param], fractionsC2[i+1][param] ) );
00195 fractionsCoeffBins.push_back( BinID( 0, fractionsC3[i][param], fractionsC3[i+1][param] ) );
00196 fractionsCoeffBins.push_back( BinID( 0, fractionsC4[i][param], fractionsC4[i+1][param] ) );
00197 fractions.push_back( new ParameterResolutions( fractionsCoeffBins, etaLow, etaHigh ) );
00198 }
00199
00200 for ( int param = 0; param < 4; param++ )
00201 {
00202 vector<BinID> correlationsCoeffBins;
00203 correlationsCoeffBins.push_back( BinID( 0, correlationsC0[i][param], correlationsC0[i+1][param] ) );
00204 correlationsCoeffBins.push_back( BinID( 0, correlationsC1[i][param], correlationsC1[i+1][param] ) );
00205 correlationsCoeffBins.push_back( BinID( 0, correlationsC2[i][param], correlationsC2[i+1][param] ) );
00206 correlationsCoeffBins.push_back( BinID( 0, correlationsC3[i][param], correlationsC3[i+1][param] ) );
00207 correlationsCoeffBins.push_back( BinID( 0, correlationsC4[i][param], correlationsC4[i+1][param] ) );
00208 correlations.push_back( new ParameterResolutions( correlationsCoeffBins, etaLow, etaHigh ) );
00209 }
00210
00211
00212
00213 PionBinData* binData = new PionBinData( rTEtaBin, core, tails, fractions, correlations, m_randSeed );
00214
00215
00216 m_binData[rTEtaBin] = binData;
00217
00218
00219 iBin++;
00220
00221 }
00222 }
00223
00224
00225 input.close();
00226
00227 }
00228 else
00229 {
00230 *m_log << MSG::INFO
00231 << "PionMatrixManager: no data file ( " << m_file << " )!!!!"
00232 << endreq;
00233 }
00234
00235 makeHeader();
00236
00237 }
00238
00239 void PionMatrixManager::makeHeader()
00240 {
00241 HeaderPrinter hp( "Atlfast PionTrack Smearer:", *m_log );
00242 hp.add( "Total number of Bins ", m_nRTBins * m_nEtaBins );
00243 hp.add( "rT Bin Min ", m_rTBoundaries.front() );
00244 hp.add( "rT Bin Max ", m_rTBoundaries.back() );
00245 hp.add( "Number of rT Bins ", m_nRTBins );
00246 hp.add( "Eta Bin Min ", m_etaBoundaries.front() );
00247 hp.add( "Eta Bin Max ", m_etaBoundaries.back() );
00248 hp.add( "Number of Eta Bins ", m_nEtaBins );
00249 hp.print();
00250 }
00251
00252
00253
00254
00255
00256 void PionMatrixManager::fillVector( ifstream& input,
00257 vector< vector<double> >& myVector,
00258 int M = 5 )
00259 {
00260 double res;
00261 for ( int param = 0; param < M; param++ )
00262 {
00263 vector< vector<double> >::iterator binIter = myVector.begin();
00264 for ( ; binIter != myVector.end(); binIter++ )
00265 {
00266 input >> res;
00267 binIter->push_back(res);
00268 }
00269
00270 }
00271 }
00272
00273
00274
00275
00276
00277
00278 vector<double> PionMatrixManager::getVariables( const TrackTrajectory& track,
00279 HepSymMatrix& returnSigma ) const
00280 {
00281 HepSymMatrix sigma;
00282 vector<double> variables;
00283
00284 IBinData* binData = getBinData(track);
00285
00286 sigma = binData->getMatrix(track);
00287
00288
00289 variables = m_correlatedData->generate( m_correlatedData->root(sigma) );
00290 returnSigma = sigma;
00291 return variables;
00292 }
00293
00294
00295
00296
00297
00298 IBinData* PionMatrixManager::getBinData( const TrackTrajectory& traj ) const
00299 {
00300 TrackParameters track = traj.parameters();
00301
00302 vector<double> rTEta;
00303 double rT = abs( traj.radius() );
00304 double eta = abs( track.eta() );
00305
00306
00307 double rTLow = ( (m_binData.begin())->first ).low(0);
00308 double rTHigh = ( (m_binData.rbegin())->first ).high(0);
00309 double etaLow = ( (m_binData.begin())->first ).low(1);
00310 double etaHigh = ( (m_binData.rbegin())->first ).high(1);
00311
00312 if ( rT < rTLow ) rT = rTLow;
00313 if ( rT > rTHigh ) rT = rTHigh;
00314 if ( eta < etaLow ) eta = etaLow;
00315 if ( eta > etaHigh ) eta = etaHigh;
00316
00317
00318 rTEta.push_back(rT);
00319 rTEta.push_back(eta);
00320
00321 map<BinID, IBinData*>::const_iterator binIter = m_binData.begin();
00322 map<BinID, IBinData*>::const_iterator binEnd = m_binData.end();
00323
00324 for ( ; binIter != binEnd; ++binIter )
00325 {
00326 if ( binIter->first.isInBin(rTEta) )
00327 {
00328 return binIter->second;
00329 }
00330 }
00331
00332 *m_log << MSG::WARNING
00333 << "WARNING: PionMatrixManager - No bin; rT " << rT << ", eta " << eta
00334 << endreq;
00335
00336 return ( m_binData.begin() )->second;
00337 }
00338
00339
00340 }