/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/Alignment/MuonAlignmentAlgorithms/plugins/MuonMillepedeAlgorithm.cc

Go to the documentation of this file.

        #include <fstream>

        #include "TFile.h"
        #include "TTree.h"
        #include "TKey.h"

        #include "FWCore/MessageLogger/interface/MessageLogger.h"
        #include "FWCore/ParameterSet/interface/ParameterSet.h"
        #include "FWCore/Framework/interface/MakerMacros.h"
        #include "FWCore/Framework/interface/Frameworkfwd.h"


        #include "TrackingTools/PatternTools/interface/Trajectory.h"
        #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"

        #include "Alignment/CommonAlignment/interface/Alignable.h"  
        #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"  
        #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"  
        #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
        #include "Alignment/CommonAlignment/interface/SurveyResidual.h"
        #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
        #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterSelector.h"
        #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentIORoot.h"
        #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
        #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
        #include <DataFormats/GeometrySurface/interface/LocalError.h> 
        #include "DataFormats/TrackReco/interface/Track.h"

        #include "Alignment/MuonAlignmentAlgorithms/plugins/MuonMillepedeAlgorithm.h"

        #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
         

        // Constructor ----------------------------------------------------------------

        MuonMillepedeAlgorithm::MuonMillepedeAlgorithm(const edm::ParameterSet& cfg):
          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" );


        }

        // Call at beginning of job ---------------------------------------------------

        void 
        MuonMillepedeAlgorithm::initialize( const edm::EventSetup& setup, 
                                            AlignableTracker* tracker, AlignableMuon* muon, 
                                            AlignableExtras* extras,
                                            AlignmentParameterStore* store )
        {
          
          edm::LogWarning("Alignment") << "[MuonMillepedeAlgorithm] Initializing...";

          // accessor Det->AlignableDet
          theAlignableDetAccessor = new AlignableNavigator(tracker, muon);

          // set alignmentParameterStore
          theAlignmentParameterStore=store;

          // get alignables
          theAlignables = theAlignmentParameterStore->alignables();

        }

        void MuonMillepedeAlgorithm::collect()
        {

          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();
        }


        // Call at end of job ---------------------------------------------------------

        void MuonMillepedeAlgorithm::terminate(void)
        {

          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)
        {
          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];
            }
          }
        }



        // Run the algorithm on trajectories and tracks -------------------------------

        void MuonMillepedeAlgorithm::run(const edm::EventSetup& setup, const EventInfo &eventInfo)
        {
          
          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

        }


        //Auxiliar
        void MuonMillepedeAlgorithm::printM(AlgebraicMatrix m)
        {
          //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::updateInfo(AlgebraicMatrix m_invCov,
                                                 AlgebraicMatrix m_weightRes,
                                                 AlgebraicMatrix m_res,
                                                 std::string id)
        {

          

          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]);
          }
        } 



DEFINE_EDM_PLUGIN( AlignmentAlgorithmPluginFactory, MuonMillepedeAlgorithm, "MuonMillepedeAlgorithm" );