#include <MuonMillepedeAlgorithm.h>

Inheritance diagram for MuonMillepedeAlgorithm:

Public Member Functions
void	collect ()

void	initialize (const edm::EventSetup &setup, AlignableTracker tracker, AlignableMuon muon, AlignableExtras extras, AlignmentParameterStore store)
	Call at beginning of job. More...

	MuonMillepedeAlgorithm (const edm::ParameterSet &cfg)
	Constructor. More...

void	run (const edm::EventSetup &setup, const EventInfo &eventInfo)
	Run the algorithm. More...

void	terminate (const edm::EventSetup &setup)
	Call at end of job. More...

void	toTMat (AlgebraicMatrix , TMatrixD )

void	updateInfo (const AlgebraicMatrix &, const AlgebraicMatrix &, const AlgebraicMatrix &, std::string)

	~MuonMillepedeAlgorithm ()
	Destructor. More...

Public Member Functions inherited from AlignmentAlgorithmBase
virtual bool	addCalibrations (const Calibrations &)

	AlignmentAlgorithmBase (const edm::ParameterSet &)
	Constructor. More...

virtual void	beginLuminosityBlock (const edm::EventSetup &setup)
	called at begin of luminosity block (no lumi block info passed yet) More...

virtual void	beginRun (const edm::EventSetup &setup)
	called at begin of run More...

virtual void	endLuminosityBlock (const edm::EventSetup &setup)
	called at end of luminosity block (no lumi block info passed yet) More...

virtual void	endRun (const EndRunInfo &runInfo, const edm::EventSetup &setup)
	called at end of run - order of arguments like in EDProducer etc. More...

virtual bool	processesEvents ()
	Returns whether algorithm proccesses events in current configuration. More...

virtual bool	setParametersForRunRange (const RunRange &rr)

virtual void	startNewLoop ()

virtual bool	supportsCalibrations ()

virtual void	terminate ()
	Called at end of job (must be implemented in derived class) More...

virtual	~AlignmentAlgorithmBase ()
	Destructor. More...

Private Member Functions
void	printM (const AlgebraicMatrix &)

Private Attributes
double	chi2nCut

std::string	collec_f

int	collec_number

std::string	collec_path

edm::Service< TFileService >	fs

std::map< std::string, TH1D * >	histoMap

bool	isCollectionJob

std::map< std::string, AlgebraicMatrix * >	map_invCov

std::map< std::string, AlgebraicMatrix * >	map_N

std::map< std::string, AlgebraicMatrix * >	map_weightRes

std::string	outputCollName

double	ptCut

AlignableNavigator *	theAlignableDetAccessor

std::vector< Alignable * >	theAlignables

AlignmentParameterStore *	theAlignmentParameterStore

bool	verbose

Additional Inherited Members
Public Types inherited from AlignmentAlgorithmBase
typedef std::pair< const Trajectory , const reco::Track >	ConstTrajTrackPair

typedef std::vector < ConstTrajTrackPair >	ConstTrajTrackPairCollection

typedef cond::RealTimeType < cond::runnumber >::type	RunNumber

typedef std::pair< RunNumber, RunNumber >	RunRange

Detailed Description

Definition at line 20 of file MuonMillepedeAlgorithm.h.

Constructor & Destructor Documentation

MuonMillepedeAlgorithm::MuonMillepedeAlgorithm ( const edm::ParameterSet & cfg )

Constructor.

Definition at line 36 of file MuonMillepedeAlgorithm.cc.

References chi2nCut, collec_f, collec_number, collec_path, edm::ParameterSet::getParameter(), isCollectionJob, outputCollName, ptCut, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                             :
       AlignmentAlgorithmBase( cfg )
     {
 
       // parse parameters
 
       edm::LogWarning("Alignment") << "[MuonMillepedeAlgorithm] constructed.";
 
       collec_f = cfg.getParameter<std::string>( "CollectionFile" );
 
       isCollectionJob = cfg.getParameter<bool>( "isCollectionJob" );
      
       collec_path = cfg.getParameter<std::string>( "collectionPath" );
       
       collec_number = cfg.getParameter<int>( "collectionNumber" );
         
       outputCollName =  cfg.getParameter<std::string>( "outputCollName" );
      
       ptCut = cfg.getParameter<double>( "ptCut" );
       
       chi2nCut = cfg.getParameter<double>( "chi2nCut" );
 
     }

MuonMillepedeAlgorithm::~MuonMillepedeAlgorithm ( )

inline

Destructor.

Definition at line 29 of file MuonMillepedeAlgorithm.h.

29 {};

Member Function Documentation

void MuonMillepedeAlgorithm::collect ( )

Definition at line 84 of file MuonMillepedeAlgorithm.cc.

References collec_f, collec_number, collec_path, cmsHarvester::index, gen::n, summarizeEdmComparisonLogfiles::objectName, outputCollName, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by terminate().

     {
 
       std::map<std::string, TMatrixD *> map;
 
       for(int c_job = 0; c_job < collec_number; ++c_job)
       {
         char name_f[40];
         snprintf(name_f, sizeof(name_f), "%s_%d/%s", collec_path.c_str(), c_job, collec_f.c_str());
         TFile file_it(name_f);
 
         if(file_it.IsZombie())
           continue;
 
         TList *m_list = file_it.GetListOfKeys();
         if(m_list == 0) {
           return;
         }
         TKey *index = (TKey *)m_list->First();
         if(index == 0) {
         }
         if( index != 0 )
         {
           do
           {
 
         std::string objectName(index->GetName());
         TMatrixD *mat = (TMatrixD *)index->ReadObj();
         std::map<std::string, TMatrixD *>::iterator node = map.find(objectName);
         if(node == map.end())
         {
           TMatrixD *n_mat = new TMatrixD(mat->GetNrows(), mat->GetNcols());
           map.insert(make_pair(objectName, n_mat));
         }
         *(map[objectName]) += *mat;
         index = (TKey*)m_list->After(index);
           } while(index != 0);
         }
         file_it.Close();    
       }
       
 
       TFile theFile2(outputCollName.c_str(), "recreate");
       theFile2.cd();
 
       std::map<std::string, TMatrixD *>::iterator m_it = map.begin();
       for(; m_it != map.end(); ++m_it)
       {  
         if(m_it->first.find("_invCov") != std::string::npos)
         {
            std::string id_s = m_it->first.substr(0, m_it->first.find("_invCov"));
            std::string id_w = id_s + "_weightRes";
            std::string id_n = id_s + "_N";
            std::string cov = id_s + "_cov";
            std::string sol = id_s + "_sol";
 
            //Covariance calculation 
            TMatrixD invMat( m_it->second->GetNrows(), m_it->second->GetNcols()); 
            invMat = *(m_it->second);
            invMat.Invert();
            //weighted residuals
            TMatrixD weightMat( m_it->second->GetNcols(), 1);
            weightMat = *(map[id_w]);
            //Solution of the linear system
            TMatrixD solution( m_it->second->GetNrows(), 1);
            solution = invMat * weightMat;
            //Number of Tracks
            TMatrixD n(1,1);
            n = *(map[id_n]);     
            
            invMat.Write(cov.c_str());
            n.Write(id_n.c_str());
            solution.Write(sol.c_str());
         }
       }
       theFile2.Write();
       theFile2.Close();
     }

void MuonMillepedeAlgorithm::initialize	(	const edm::EventSetup &	setup,
		AlignableTracker *	tracker,
		AlignableMuon *	muon,
		AlignableExtras *	extras,
		AlignmentParameterStore *	store
	)

virtual

Call at beginning of job.

Implements AlignmentAlgorithmBase.

Definition at line 65 of file MuonMillepedeAlgorithm.cc.

References AlignmentParameterStore::alignables(), theAlignableDetAccessor, theAlignables, and theAlignmentParameterStore.

     {
       
       edm::LogWarning("Alignment") << "[MuonMillepedeAlgorithm] Initializing...";
 
       // accessor Det->AlignableDet
       theAlignableDetAccessor = new AlignableNavigator(tracker, muon);
 
       // set alignmentParameterStore
       theAlignmentParameterStore=store;
 
       // get alignables
       theAlignables = theAlignmentParameterStore->alignables();
 
     }

void MuonMillepedeAlgorithm::printM ( const AlgebraicMatrix & m )

private

Definition at line 480 of file MuonMillepedeAlgorithm.cc.

     {
       //for(int i = 0; i < m.num_row(); ++i)
       // {
       //  for(int j = 0; j < m.num_col(); ++j)
       //  {
       //    std::cout << m[i][j] << " ";
       //  }
       //  std::cout << std::endl;
       //}
     }

void MuonMillepedeAlgorithm::run	(	const edm::EventSetup &	setup,
		const EventInfo &	eventInfo
	)

virtual

Run the algorithm.

Implements AlignmentAlgorithmBase.

Definition at line 232 of file MuonMillepedeAlgorithm.cc.

References CompositeAlignmentParameters::alignableFromAlignableDet(), AlignableNavigator::alignableFromGeomDet(), AlignableNavigator::alignablesFromHits(), Alignable::alignmentParameters(), chi2n, chi2nCut, AlignableNavigator::detAndSubdetInMap(), eta, reco::TrackBase::eta(), TrajectoryMeasurement::forwardPredictedState(), GeomDet::geographicalId(), TrackingRecHit::geographicalId(), isCollectionJob, TrajectoryStateOnSurface::isValid(), TrackingRecHit::isValid(), Trajectory::measurements(), mergeVDriftHistosByStation::name, reco::TrackBase::normalizedChi2(), reco::TrackBase::numberOfValidHits(), phi, reco::TrackBase::phi(), EnergyCorrector::pt, reco::TrackBase::pt(), ptCut, TrajectoryMeasurement::recHit(), AlignmentParameters::selectedDerivatives(), corrVsCorr::selection, AlignmentParameterStore::selectParameters(), AlCaHLTBitMon_QueryRunRegistry::string, theAlignableDetAccessor, theAlignmentParameterStore, testEve_cfg::tracks, AlignmentAlgorithmBase::EventInfo::trajTrackPairs(), and updateInfo().

Referenced by Types.EventID::cppID(), and Types.LuminosityBlockID::cppID().

         {
       
       if( isCollectionJob )
       {
         return;
       }
 
 
       // loop over tracks  
       //int t_counter = 0;
       const ConstTrajTrackPairCollection &tracks = eventInfo.trajTrackPairs();
       for( ConstTrajTrackPairCollection::const_iterator it=tracks.begin();
            it!=tracks.end();it++) {
 
         const Trajectory* traj = (*it).first;
         const reco::Track* track = (*it).second;
 
         float pt    = track->pt();
         float eta   = track->eta();
         float phi   = track->phi();
         float chi2n = track->normalizedChi2();
         //int   ndof = track->ndof();
         int   nhit  = track->numberOfValidHits();
 
         if (0) edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: " << pt <<","<< eta <<","<< phi <<","<< chi2n << ","<<nhit;
 
         //Accept or not accept the track
         if( pt > ptCut && chi2n < chi2nCut ) 
         {
 
 
           std::vector<const TransientTrackingRecHit*> hitvec;
           std::vector<TrajectoryStateOnSurface> tsosvec;
           
           std::vector<TrajectoryMeasurement> measurements = traj->measurements();
          
           //In this loop the measurements and hits are extracted and put on two vectors 
           for (std::vector<TrajectoryMeasurement>::iterator im=measurements.begin();
              im!=measurements.end(); im++)
           {
         TrajectoryMeasurement meas = *im;
         const TransientTrackingRecHit* hit = &(*meas.recHit());
         //We are not very strict at this point
         if (hit->isValid()  && theAlignableDetAccessor->detAndSubdetInMap( hit->geographicalId() ))
         {
           //***Forward
           TrajectoryStateOnSurface tsos = meas.forwardPredictedState();
           if (tsos.isValid())
           {
             hitvec.push_back(hit);
             tsosvec.push_back(tsos);
           }
         }
           }
         
           // transform RecHit vector to AlignableDet vector
           std::vector <AlignableDetOrUnitPtr> alidetvec = 
         theAlignableDetAccessor->alignablesFromHits(hitvec);
 
           // get concatenated alignment parameters for list of alignables
           CompositeAlignmentParameters aap = 
         theAlignmentParameterStore->selectParameters(alidetvec);
 
           std::vector<TrajectoryStateOnSurface>::const_iterator itsos=tsosvec.begin();
           std::vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin();
 
           
           //int ch_counter = 0;
 
           while (itsos != tsosvec.end()) 
           {
         // get AlignableDet for this hit
         const GeomDet* det=(*ihit)->det();
         AlignableDetOrUnitPtr alidet = 
           theAlignableDetAccessor->alignableFromGeomDet(det);
           
         // get relevant Alignable
         Alignable* ali=aap.alignableFromAlignableDet(alidet);
          
         //To be sure that the ali is not null and that it's a DT segment 
         if ( ali!=0 && (*ihit)->geographicalId().subdetId() == 1)
         {
           DTChamberId m_Chamber(det->geographicalId());
           //Station 4 does not contain Theta SL 
           if((*ihit)->dimension() == 4 || ((*ihit)->dimension() == 2 && m_Chamber.station() == 4))
           //if((*ihit)->dimension() == 4)
           {
             edm::LogInfo("Alignment") << "Entrando";
           
             AlignmentParameters* params = ali->alignmentParameters();
           
             edm::LogInfo("Alignment") << "Entrando";
             //(dx/dz,dy/dz,x,y) for a 4DSegment
             AlgebraicVector ihit4D = (*ihit)->parameters();
 
             //The innerMostState always contains the Z
             //(q/pt,dx/dz,dy/dz,x,y)
             AlgebraicVector5 alivec = (*itsos).localParameters().mixedFormatVector();
             
             //The covariance matrix follows the sequence
             //(q/pt,dx/dz,dy/dz,x,y) but we reorder to
             //(x,y,dx/dz,dy/dz)
             AlgebraicSymMatrix55 rawCovMat = (*itsos).localError().matrix();
             AlgebraicMatrix CovMat(4,4);
             int m_index[] = {2,3,0,1};
             for(int c_ei = 0; c_ei < 4; ++c_ei)
             {
               for(int c_ej = 0; c_ej < 4; ++c_ej)
               {
             CovMat[m_index[c_ei]][m_index[c_ej]] = rawCovMat(c_ei+1,c_ej+1);
               }
             }
             
             int inv_check;
             //printM(CovMat);
             CovMat.invert(inv_check);
             if (inv_check != 0) { 
               edm::LogError("Alignment") << "Covariance Matrix inversion failed"; 
               return; 
             } 
 
 
 
             //The order is changed to:
             // (x,0,dx/dz,0) MB4 Chamber
             // (x,y,dx/dz,dy/dz) Not MB4 Chamber
             AlgebraicMatrix residuals(4,1);
             if(m_Chamber.station() == 4)
             {
               //Filling Residuals
               residuals[0][0] = ihit4D[2]-alivec[3];
               residuals[1][0] = 0.0;
               residuals[2][0] = ihit4D[0]-alivec[1];
               residuals[3][0] = 0.0;
               //The error in the Theta coord is set to infinite
               CovMat[1][0] = 0.0; CovMat[1][1] = 0.0; CovMat[1][2] = 0.0;
               CovMat[1][3] = 0.0; CovMat[0][1] = 0.0; CovMat[2][1] = 0.0; 
               CovMat[3][1] = 0.0; CovMat[3][0] = 0.0; CovMat[3][2] = 0.0;
               CovMat[3][3] = 0.0; CovMat[0][3] = 0.0; CovMat[2][3] = 0.0; 
             }
             else
             {
               //Filling Residuals
               residuals[0][0] = ihit4D[2]-alivec[3];
               residuals[1][0] = ihit4D[3]-alivec[4];
               residuals[2][0] = ihit4D[0]-alivec[1];
               residuals[3][0] = ihit4D[1]-alivec[2];
             }
       
             // Derivatives 
             AlgebraicMatrix derivsAux = params->selectedDerivatives(*itsos,alidet);
             
             std::vector<bool> mb4_mask;
             std::vector<bool> selection;
             //To be checked 
             mb4_mask.push_back(true); mb4_mask.push_back(false); mb4_mask.push_back(true); 
             mb4_mask.push_back(true); mb4_mask.push_back(true); mb4_mask.push_back(false); 
             selection.push_back(true); selection.push_back(true); selection.push_back(false); 
             selection.push_back(false); selection.push_back(false); selection.push_back(false); 
             int nAlignParam = 0;
             if(m_Chamber.station() == 4)
             {  
               for(int icor = 0; icor < 6; ++icor)
               {
             if(mb4_mask[icor] && selection[icor]) nAlignParam++;                      
               }
             }
             else
             {
               nAlignParam = derivsAux.num_row();
             }
             
             AlgebraicMatrix derivs(nAlignParam, 4);
             if(m_Chamber.station() == 4)
             { 
               int der_c = 0;
               for(int icor = 0; icor < 6; ++icor)
               {
             if(mb4_mask[icor] && selection[icor])
             {
                for(int ccor = 0; ccor < 4; ++ccor)
                {
                  derivs[der_c][ccor] = derivsAux[icor][ccor];
                  ++der_c;
                }
             } 
               }
             }
             else
             {
               derivs = derivsAux;
             }
      
 
             //for(int co = 0; co < derivs.num_row(); ++co)
             //{
             //  for(int ci = 0; ci < derivs.num_col(); ++ci)
             //  {
             //     edm::LogInfo("Alignment") << "Derivatives: " << co << " " << ci << " " << derivs[co][ci] << " ";
             //  }
             //} 
       
             AlgebraicMatrix derivsT = derivs.T();
             AlgebraicMatrix invCov = derivs*CovMat*derivsT;
             AlgebraicMatrix weightRes = derivs*CovMat*residuals;
 
             //this->printM(derivs);
             //this->printM(CovMat);
             //this->printM(residuals);
             
             char name[40];
             snprintf(name, sizeof(name),
                  "Chamber_%d_%d_%d", m_Chamber.wheel(), m_Chamber.station(), m_Chamber.sector());
             std::string chamId(name);
             //MB4 need a special treatment
             /*AlgebraicMatrix invCovMB4(2,2);
             AlgebraicMatrix weightResMB4(2,1); 
             if( m_Chamber.station() == 4 )
             { 
               int m_index_2[] = {0,2};
               for(int c_i = 0; c_i < 2; ++c_i)
               {
             weightResMB4[c_i][0] = weightRes[m_index_2[c_i]][0];
             for(int c_j = 0; c_j < 2; ++c_j)
             {
               invCovMB4[c_i][c_j] = invCov[m_index_2[c_i]][m_index_2[c_j]];
             }  
               }
               this->updateInfo(invCovMB4, weightResMB4, residuals, chamId); 
             }
             else
             {
               this->updateInfo(invCov, weightRes, residuals, chamId); 
             }*/
             this->updateInfo(invCov, weightRes, residuals, chamId); 
           }
         }
         itsos++;
         ihit++;
           }
         }
       } // end of track loop
 
     }

void MuonMillepedeAlgorithm::terminate ( const edm::EventSetup & setup )

virtual

Call at end of job.

Implements AlignmentAlgorithmBase.

Definition at line 166 of file MuonMillepedeAlgorithm.cc.

References Alignable::alignmentParameters(), AlignmentParameterStore::applyParameters(), collec_f, collect(), isCollectionJob, map_invCov, map_N, map_weightRes, AlignmentParameters::setValid(), theAlignables, theAlignmentParameterStore, interactiveExample::theFile, and toTMat().

     {
 
       if( isCollectionJob )
       {
          this->collect();
          return;
       }
 
 
       edm::LogWarning("Alignment") << "[MuonMillepedeAlgorithm] Terminating";
 
       // iterate over alignment parameters
       for(std::vector<Alignable*>::const_iterator
         it=theAlignables.begin(); it!=theAlignables.end(); it++) {
         Alignable* ali=(*it);
         // Alignment parameters
         // AlignmentParameters* par = ali->alignmentParameters();
         edm::LogInfo("Alignment") << "now apply params";
         theAlignmentParameterStore->applyParameters(ali);
         // set these parameters 'valid'
         ali->alignmentParameters()->setValid(true);
 
       }
 
       edm::LogWarning("Alignment") << "[MuonMillepedeAlgorithm] Writing aligned parameters to file: " << theAlignables.size();
 
       TFile *theFile = new TFile(collec_f.c_str(), "recreate");
       theFile->cd();
       std::map<std::string, AlgebraicMatrix *>::iterator invCov_it = map_invCov.begin();
       std::map<std::string, AlgebraicMatrix *>::iterator weightRes_it = map_weightRes.begin();
       std::map<std::string, AlgebraicMatrix *>::iterator n_it = map_N.begin();
       for(; n_it != map_N.end(); ++invCov_it, ++weightRes_it, ++n_it)
       {
         TMatrixD tmat_invcov(0,0);
         this->toTMat(invCov_it->second, &tmat_invcov);
         TMatrixD tmat_weightres(0,0);
         this->toTMat(weightRes_it->second, &tmat_weightres);
         TMatrixD tmat_n(0,0);
         this->toTMat(n_it->second, &tmat_n);
         
         tmat_invcov.Write(invCov_it->first.c_str());
         tmat_weightres.Write(weightRes_it->first.c_str());
         tmat_n.Write(n_it->first.c_str());
       }
 
       theFile->Write();
       theFile->Close();
 
     }

void MuonMillepedeAlgorithm::toTMat	(	AlgebraicMatrix *	am_mat,
		TMatrixD *	tmat_mat
	)

Definition at line 218 of file MuonMillepedeAlgorithm.cc.

Referenced by terminate().

     {
       tmat_mat->ResizeTo(am_mat->num_row(), am_mat->num_col());
       for(int c_i = 0; c_i < am_mat->num_row(); ++c_i) {
         for(int c_j = 0; c_j < am_mat->num_col(); ++c_j) {
           (*tmat_mat)(c_i, c_j) = (*am_mat)[c_i][c_j];
         }
       }
     }

void MuonMillepedeAlgorithm::updateInfo	(	const AlgebraicMatrix &	m_invCov,
		const AlgebraicMatrix &	m_weightRes,
		const AlgebraicMatrix &	m_res,
		std::string	id
	)

Definition at line 492 of file MuonMillepedeAlgorithm.cc.

References fs, timingPdfMaker::histo, histoMap, cutBasedElectronHLTPreselecition_Summer16_V1_cff::idName, TFileService::make(), map_invCov, map_N, map_weightRes, mergeVDriftHistosByStation::name, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by run().

     {
 
       
 
       std::string id_invCov = id + "_invCov";
       std::string id_weightRes = id + "_weightRes";
       std::string id_n = id + "_N";
       
       edm::LogInfo("Alignment") << "Entrando";
 
       std::map<std::string, AlgebraicMatrix *>::iterator node = map_invCov.find(id_invCov);
       if(node == map_invCov.end())
       {
         AlgebraicMatrix *f_invCov = new AlgebraicMatrix(m_invCov.num_row(), m_invCov.num_col());
         AlgebraicMatrix *f_weightRes = new AlgebraicMatrix(m_weightRes.num_row(), m_weightRes.num_col());
         AlgebraicMatrix *f_n = new AlgebraicMatrix(1,1);
 
         map_invCov.insert(make_pair(id_invCov, f_invCov));
         map_weightRes.insert(make_pair(id_weightRes, f_weightRes));
         map_N.insert(make_pair(id_n, f_n));
 
         for(int iCount = 0; iCount < 4; ++iCount)
         {
           char name[40];
           snprintf(name, sizeof(name), "%s_var_%d", id.c_str(), iCount);
           std::string idName(name);
           float range = 5.0;
           //if( iCount == 0 || iCount == 1 ) {
           //  range = 0.01;
           //}
           TH1D *histo = fs->make<TH1D>(idName.c_str(), idName.c_str(), 200, -range, range );
           histoMap.insert(make_pair(idName, histo));
         }
       }
       
       *map_invCov[id_invCov] = *map_invCov[id_invCov] + m_invCov;
       *map_weightRes[id_weightRes] = *map_weightRes[id_weightRes] + m_weightRes;
       (*map_N[id_n])[0][0]++;
 
       for(int iCount = 0; iCount < 4; ++iCount)
       {
         char name[40];
         snprintf(name, sizeof(name), "%s_var_%d", id.c_str(), iCount);
         std::string idName(name); 
         histoMap[idName]->Fill(m_res[iCount][0]);
       }
     } 