Atlfast::BremFitter class Reference

00013                                                       :m_log(log)*/ {
00014     
00015     // HOW IT WORKS ---------------
00016     //  pT is scaled according to Eta of the track
00017     // There exist 3D Eta, R and PT bins.
00018     // For a track Eta Bin, there are n bremR values.
00019     // for each bremR value there are 50 p Values.
00020     // The bremR value is chosen at random.
00021     // The p Value within this bremR value is also chosen at random
00022     // The pT scale is calculated from this random number and the pT bin edges
00023 
00024     m_correlatedData = new CorrelatedData(randSeed);
00025     m_randomEngine = new HepJamesRandom(randSeed);
00026     //read files and make variables
00027     // open file
00028     m_file = "Atlfast_XKalInput.dat";
00029     ifstream input;
00030     input.open( m_file.c_str() );
00031     if (input) {
00032       
00033       log << MSG::INFO <<"BremFitter: File "<<m_file<<" open."<<endreq;
00034       log << setprecision(3);
00035       vector<double> empty;
00036       string skip;
00037       // read some parameters
00038       readSingleParameter(log, "RMin",input,m_rMin);
00039       readSingleParameter(log, "RMax",input,m_rMax);
00040       readSingleParameter(log, "N RBins",input,m_nRBins);
00041       readSingleParameter(log, "Eta Extent",input,m_etaExtent);
00042       readSingleParameter(log, "N Eta Bins",input,m_nEtaBins);
00043       
00044       // read BremR values
00045       double bremR;
00046       vector<double> bremRStore;
00047       input >> skip;
00048       log << MSG::DEBUG << "Reading " << m_nRBins*m_nEtaBins <<" "<< skip << " values" << endreq;
00049       for (int i=0; i < (m_nRBins*m_nEtaBins ); ++i) {
00050           input >> bremR;
00051           log << MSG::DEBUG << i << ": " << bremR << endreq;
00052           bremRStore.push_back(bremR);
00053       }
00054 
00055       // make displacements
00056       readDisplacements(log, input, bremRStore);
00057 
00058       // read more parameters
00059       readSingleParameter(log, "PMin",input,m_pMin);
00060       readSingleParameter(log, "PMax",input,m_pMax);
00061       readSingleParameter(log, "N PBins",input,m_nPBins);
00062       readSingleParameter(log, "P Scale 1",input,m_p1scal);
00063       readSingleParameter(log, "P Scale 2",input,m_p2scal);
00064 
00065       double pTScaleHigh;
00066       double pTScaleLow;
00067       int pTID=1, etaID=1;
00068       double pTValue;
00069       double step = (m_pMax - m_pMin) / (m_nPBins);
00070       vector<double>::iterator  bremRIter = bremRStore.begin();
00071 
00072       // make pt distribution EtaBin objects
00073       //====================================
00074       // there are nEtaBins bins
00075       for (int iEta=0; iEta < m_nEtaBins; iEta++){
00076         BremEtaBin* etaBin = new BremEtaBin(randSeed);
00077 
00078         // each bin has nRBins number of RBins
00079         for (int iR=0; iR < m_nRBins; iR++){
00080 
00081              // skip comments in file
00082              input >> skip;
00083              input >> skip;
00084              input >> skip;
00085              input >> skip;
00086              log << MSG::DEBUG << "Reading pT values " << skip << " at radial bin ";
00087              input >> skip;
00088              log << MSG::DEBUG << skip << ", eta bin ";
00089              input >> skip;
00090              log << MSG::DEBUG << skip << endreq;
00091              
00092              //-----
00093              BremRBin* rBin = new BremRBin(randSeed);
00094              pTScaleLow =0;
00095              pTValue = m_pMin;
00096           // each rBin has nPBins number of PTBins
00097 
00098           for (int iP=0; iP < m_nPBins; iP++){
00099             input >> pTScaleHigh;
00100             log << MSG::DEBUG << setw(10) << "("<<pTScaleLow<<" - "<<pTScaleHigh<<
00101               ", "<<pTID<<", "<< pTValue<<","<<step<<")"<<endreq;
00102             rBin->addBin(TestValue(++pTID, pTScaleHigh), BremPTScaleBin( pTScaleLow, pTScaleHigh, pTValue, step ) ); 
00103             pTValue += step;
00104             pTScaleLow = pTScaleHigh;
00105           }
00106 
00107           log << MSG::DEBUG << "Adding R Bin "<<etaID<<" with Brem Radius "<<*bremRIter<<endreq;
00108           etaBin->addBin( TestValue(++etaID, *bremRIter), rBin );
00109           ++bremRIter;
00110 
00111         }
00112         m_pTDistEtaBins.push_back(etaBin);
00113       }
00114       vector<BremEtaBin*>::const_iterator iter = m_pTDistEtaBins.begin();
00115       vector<BremEtaBin*>::const_iterator end = m_pTDistEtaBins.end();
00116       log << MSG::DEBUG << "Electron Brem PT Distribution Bins" << endreq;
00117       int ieta=0;
00118       for(;iter != end; ++iter) {
00119         ++ieta;
00120         (log) << MSG::DEBUG <<"Eta "<<ieta<<" "<< (**iter) << endreq;
00121       }
00122       log << MSG::INFO << "finished reading "<<m_file<<endreq;
00123     }
00124     else { log <<MSG::WARNING << "BremFitter: no data file ( "
00125                <<m_file<<" )!!!!" << endreq;}
00126 
00127     
00128     // done reading files!!    
00129   }

00021                   {
00022       //      delete m_randomEngine; delete m_correlatedData;
00023     }

00131                                                                         {
00132 
00133     TrackParameters track = traj.parameters();
00134     int iEta;
00135     double oldQPT,newQPT, newD0, newPhi;
00136     double pScale, deltaP;
00137     //    double pStep;
00138 
00139     // get eta index and find pScale
00140     iEta = int(m_nEtaBins*abs(track.eta())/ m_etaExtent);
00141     if (iEta >= m_nEtaBins) iEta = m_nEtaBins-1;
00142 
00143     // -------- Get PSCALE Random --------------
00144     double random = m_randomEngine->flat();
00145 
00146     pScale = m_pTDistEtaBins[iEta]->calculatePScale(random);
00147 
00148     // scale pT
00149     oldQPT = track.pT() * traj.signOfCharge();
00150     newQPT = scalePT(oldQPT, pScale);
00151 
00152     // do a check
00153     if (newQPT != 0 && oldQPT != 0) deltaP = 1/newQPT - 1/oldQPT;
00154     else deltaP = 0;
00155 
00156     // generate phi and d0 displacements
00157     newD0 = scaleDPhi(iEta, random, oldQPT, track.impactParameter(),
00158                       m_d0Disp, deltaP);
00159     newPhi = scaleDPhi(iEta, random, oldQPT, track.phi(), m_phiDisp, deltaP);
00160 
00161     // make new track
00162     double newInvPtCharge = 1 / newQPT;
00163 
00164     TrackParameters newParams(track.eta(),
00165                               newPhi,
00166                               abs(newQPT), 
00167                               newD0,
00168                               track.zPerigee(),
00169                               track.cotTheta(), 
00170                               newInvPtCharge);
00171 
00172     TrackTrajectory bremTrack(newParams, traj.startPoint(), traj.curvature());
00173     return bremTrack;
00174   }

00179   {
00180     double pscale, scale = 20; //GeV scale
00181     double random = m_correlatedData->generate(1);
00182     double numerator =  pT*pT + axis.off()*axis.off();
00183     double denominator = scale*scale + axis.off()*axis.off();
00184     pscale = pow( numerator/denominator, axis.pow()/2 );
00185     double mean = axis.average(iEta,randBremR);
00186     double rms = axis.rms(iEta, randBremR);
00187     double grad = pscale * (mean + random*rms);
00188 
00189     return old + grad*deltaP;  
00190   }

00208                                                                                                  {
00209         string skip;
00210         input >> skip;
00211         input >> parameter;
00212         //m_log << MSG::INFO << "Read " << skip << " with value " << parameter << endreq;
00213         s << MSG::DEBUG << "Loading " << name <<": Read " << skip << " with value " << parameter << endreq;
00214       }

00201                                                                                                     {
00202         string skip;
00203         input >> skip;
00204         input >> parameter;
00205         //m_log << MSG::INFO << "Read " << skip << " with value " << parameter << endreq;
00206         s << MSG::DEBUG << "Loading " << name <<": Read " << skip << " with value " << parameter << endreq;
00207       }

00216                                                                                                {
00217 
00218       double dtpow,pdoff,ftpow,pfoff;
00219       string skip;
00220 
00221       readSingleParameter(log, "D0 Power",input,dtpow);
00222       readSingleParameter(log, "D0 Offset",input,pdoff);
00223       readSingleParameter(log, "Phi Power",input,ftpow);
00224       readSingleParameter(log, "Phi Offset",input,pfoff);
00225 
00226       // make displacements
00227       // 4 x n(RBins) (corresponding to 5(eta) lots of the same n(RBins) ) values
00228       double data;
00229       vector< vector< map<double, double> > > axes;
00230       
00231       // read the 4 vectors(d0 mean, rms and phi mean, rms)
00232       for (int nVectors = 0; nVectors < 4; ++nVectors) {
00233       vector<double>::iterator  bremRIter = bremRStore.begin();
00234         input >> skip;
00235         log << MSG::DEBUG << "Reading " << m_nRBins << " " << skip << " values" << endreq;
00236         vector<double> vecData;
00237         // read parameters
00238         for (int iR=0; iR < (m_nRBins) ; ++iR){
00239           input >> data;
00240           log << MSG::DEBUG << iR << ": " << data << endreq;
00241           vecData.push_back(data);
00242         }
00243         // these n(RBins) parameters are the same for each eta bin, with different bremR values
00244         // run through each eta bin and make map for each
00245         vector< map<double,double> > vecDataMap;
00246         vector<double>::iterator iData, iDataEnd = vecData.end();
00247         
00248         // make n(eta) lots of data maps
00249         for (int iE=0; iE < (m_nEtaBins) ; ++iE){
00250           map<double,double> dataMap;
00251           iData = vecData.begin();
00252           for (; iData != iDataEnd; ++iData) {
00253             dataMap[*bremRIter] = *iData;
00254             ++bremRIter;
00255           }
00256           vecDataMap.push_back(dataMap);
00257         }
00258         
00259         axes.push_back(vecDataMap);
00260       }
00261       // Done all four
00262       
00263       m_d0Disp = Displacement(pdoff, dtpow, axes[0], axes[1]);
00264       m_phiDisp = Displacement(pfoff, ftpow, axes[2], axes[3]);
00265       //========================================
00266   }

00193   { 
00194     // do some maths (soft pT scaling normalised to 20GeV)
00195     double scale  = 1 + m_p2scal*log(abs(pT)/20);
00196     pscale = m_p1scal * (1 + scale*(pscale-1));
00197     if (pscale < 0) pscale=1;
00198     return ( pscale * pT );
00199   }

Atlfast::BremFitter Class Reference

Public Methods

Private Methods

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Methods
	BremFitter (int, MsgStream &)
	~BremFitter ()
TrackTrajectory	doBremFit (const TrackTrajectory &) const
Private Methods
double	scaleDPhi (int iEta, double random, double pT, double old, const Displacement &axis, double deltaP) const
void	readSingleParameter (MsgStream &, string, ifstream &, int &parameter)
void	readSingleParameter (MsgStream &, string, ifstream &, double &parameter)
void	readDisplacements (MsgStream &, ifstream &, vector< double > bremValues)
double	scalePT (double pT, double pscale) const
Private Attributes
Displacement	m_d0Disp
Displacement	m_phiDisp
string	m_file
int	m_nRBins
int	m_nEtaBins
int	m_nPBins
double	m_rMin
double	m_rMax
double	m_pMin
double	m_pMax
double	m_p1scal
double	m_p2scal
double	m_etaExtent
vector< BremEtaBin * >	m_pTDistEtaBins
CorrelatedData *	m_correlatedData
HepJamesRandom *	m_randomEngine