#include <MuonGeometryArrange.h>

Inheritance diagram for MuonGeometryArrange:

Classes
struct	MGACollection
Public Types
typedef std::vector< Alignable * >	Alignables
typedef AlignTransform	SurveyValue
typedef Alignments	SurveyValues
Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
	Read from DB and print survey info.
	MuonGeometryArrange (const edm::ParameterSet &)
	Do nothing. Required by framework.
Private Member Functions
bool	checkChosen (Alignable *ali)
void	compare (Alignable refAli, Alignable curAli, Alignable *curAliCopy2)
void	compareGeometries (Alignable refAli, Alignable curAli, Alignable *curAliCopy2)
void	createROOTGeometry (const edm::EventSetup &iSetup)
void	endHist ()
void	fillTree (Alignable *refAli, AlgebraicVector diff)
bool	isMother (Alignable *ali)
void	makeGraph (int sizeI, float smi, float sma, float minV, float maxV, TH2F dxh, TGraph grx, const char name, const char title, const char titleg, const char axis, float xp, float yp, int numEntries)
bool	passChosen (Alignable *ali)
bool	passIdCut (uint32_t)
Private Attributes
TTree *	_alignTree
float	_alphaVal
float	_betaVal
align::StructureType	_commonMuonLevel
float	_dalphaVal
float	_dbetaVal
int	_detDim
bool	_detIdFlag
std::string	_detIdFlagFile
std::vector< uint32_t >	_detIdFlagVector
float	_dgammaVal
float	_dphiVal
float	_drotxVal
float	_drotyVal
float	_drotzVal
float	_drVal
float	_dxVal
float	_dyVal
float	_dzVal
int	_endcap
float	_etaVal
std::string	_filename
float	_gammaVal
int	_id
std::string	_inputFilename1
std::string	_inputFilename2
TFile *	_inputRootFile1
TFile *	_inputRootFile2
TTree *	_inputTree1
TTree *	_inputTree2
std::string	_inputTreename
std::string	_inputXMLCurrent
std::string	_inputXMLReference
float	_ldphiVal
float	_ldrVal
float	_ldxVal
float	_ldyVal
float	_ldzVal
int	_level
std::vector< MGACollection >	_mgacollection
int	_mid
int	_mlevel
align::PositionType	_MuonCommonCM
align::EulerAngles	_MuonCommonR
align::GlobalVector	_MuonCommonT
float	_phiVal
int	_ring
float	_rotxVal
float	_rotyVal
float	_rotzVal
float	_rVal
std::string	_setCommonMuonSystem
int	_station
int	_sublevel
float	_surLength
double	_surRot [9]
float	_surWidth
TFile *	_theFile
int	_useDetId
std::string	_weightBy
bool	_weightById
std::string	_weightByIdFile
std::vector< unsigned int >	_weightByIdVector
bool	_writeToDB
float	_xVal
float	_yVal
float	_zVal
AlignableMuon *	currentMuon
AlignableMuon *	dummyMuon
bool	firstEvent_
MuonAlignment *	inputAlign1
MuonAlignment *	inputAlign2
MuonAlignment *	inputAlign2a
Alignable *	inputGeometry1
Alignable *	inputGeometry2
edm::ParameterSet	m_params
AlignableMuon *	referenceMuon
std::vector< align::StructureType >	theLevels
const SurveyErrors *	theSurveyErrors
unsigned int	theSurveyIndex
const Alignments *	theSurveyValues

Detailed Description

Module that reads survey info from DB and prints them out.

Usage: module comparator = MuonGeometryArrange {

lots of stuff

} path p = { comparator }

Date:: 2010/01/04 17:04:08

Revision:: 1.3

Author:: Nhan Tran

Definition at line 38 of file MuonGeometryArrange.h.

Member Typedef Documentation

typedef std::vector<Alignable*> MuonGeometryArrange::Alignables

Definition at line 44 of file MuonGeometryArrange.h.

typedef AlignTransform MuonGeometryArrange::SurveyValue

Definition at line 42 of file MuonGeometryArrange.h.

typedef Alignments MuonGeometryArrange::SurveyValues

Definition at line 43 of file MuonGeometryArrange.h.

Constructor & Destructor Documentation

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

Do nothing. Required by framework.

Definition at line 43 of file MuonGeometryArrange.cc.

References _alignTree, _alphaVal, _betaVal, _dalphaVal, _dbetaVal, _detDim, _detIdFlag, _detIdFlagFile, _detIdFlagVector, _dgammaVal, _dphiVal, _drotxVal, _drotyVal, _drotzVal, _drVal, _dxVal, _dyVal, _dzVal, _endcap, _etaVal, _filename, _gammaVal, _id, _inputFilename1, _inputFilename2, _inputTreename, _inputXMLCurrent, _inputXMLReference, _ldphiVal, _ldrVal, _ldxVal, _ldyVal, _ldzVal, _level, _mgacollection, _mid, _mlevel, _phiVal, _ring, _rotxVal, _rotyVal, _rotzVal, _rVal, _station, _sublevel, _surLength, _surRot, _surWidth, _theFile, _useDetId, _weightBy, _weightById, _weightByIdFile, _weightByIdVector, _xVal, _yVal, _zVal, currentMuon, groupFilesInBlocks::fin, edm::ParameterSet::getUntrackedParameter(), errorMatrix2Lands_multiChannel::id, prof2calltree::l, referenceMuon, and theLevels.

                                                                   :
  theSurveyIndex(0), _writeToDB(false), _commonMuonLevel(align::invalid), firstEvent_(true)
{
        referenceMuon=0x0;
        currentMuon=0x0;
        // Input is XML
        _inputXMLCurrent = cfg.getUntrackedParameter<std::string> ("inputXMLCurrent");
        _inputXMLReference = cfg.getUntrackedParameter<std::string> ("inputXMLReference");

        //input is ROOT
        _inputFilename1 = cfg.getUntrackedParameter< std::string >
                                ("inputROOTFile1");
        _inputFilename2 = cfg.getUntrackedParameter< std::string >
                                ("inputROOTFile2");
        _inputTreename = cfg.getUntrackedParameter< std::string > ("treeName");
        
        //output file
        _filename = cfg.getUntrackedParameter< std::string > ("outputFile");
        
        
        const std::vector<std::string>& levels = 
            cfg.getUntrackedParameter< std::vector<std::string> > ("levels");
        
        _weightBy = cfg.getUntrackedParameter< std::string > ("weightBy");
        _detIdFlag = cfg.getUntrackedParameter< bool > ("detIdFlag");
        _detIdFlagFile = cfg.getUntrackedParameter< std::string >
                          ("detIdFlagFile");
        _weightById  = cfg.getUntrackedParameter< bool > ("weightById");
        _weightByIdFile = cfg.getUntrackedParameter< std::string >
                         ("weightByIdFile");
        _endcap  = cfg.getUntrackedParameter<int> ("endcapNumber");
        _station = cfg.getUntrackedParameter<int> ("stationNumber");
        _ring    = cfg.getUntrackedParameter<int> ("ringNumber");
        
        //setting the levels being used in the geometry comparator
        AlignableObjectId dummy;
        edm::LogInfo("MuonGeometryArrange") << "levels: " << levels.size();
        for (unsigned int l = 0; l < levels.size(); ++l){
                theLevels.push_back( dummy.nameToType(levels[l]));
                edm::LogInfo("MuonGeometryArrange") << "level: " << levels[l];
        }
        
                
        // if want to use, make id cut list
        if (_detIdFlag){
        ifstream fin;
        fin.open( _detIdFlagFile.c_str() );
        
        while (!fin.eof() && fin.good() ){
                        
                        uint32_t id;
                        fin >> id;
                        _detIdFlagVector.push_back(id);
        }
        fin.close();
        }               
        
        // turn weightByIdFile into weightByIdVector
        unsigned int lastID=999999999;
        if (_weightById){
                std::ifstream inFile;
                inFile.open( _weightByIdFile.c_str() );
                int ctr = 0;
                while ( !inFile.eof() ){
                        ctr++;
                        unsigned int listId;
                        inFile >> listId;
                        inFile.ignore(256, '\n');
                        if(listId!=lastID){
                          _weightByIdVector.push_back( listId );
                        }
                        lastID=listId;
                }
                inFile.close();
        }
        
        
        
        //root configuration
        _theFile = new TFile(_filename.c_str(),"RECREATE");
        _alignTree = new TTree("alignTree","alignTree");
        _alignTree->Branch("id", &_id, "id/I");
        _alignTree->Branch("level", &_level, "level/I");
        _alignTree->Branch("mid", &_mid, "mid/I");
        _alignTree->Branch("mlevel", &_mlevel, "mlevel/I");
        _alignTree->Branch("sublevel", &_sublevel, "sublevel/I");
        _alignTree->Branch("x", &_xVal, "x/F");
        _alignTree->Branch("y", &_yVal, "y/F");
        _alignTree->Branch("z", &_zVal, "z/F");
        _alignTree->Branch("r", &_rVal, "r/F");
        _alignTree->Branch("phi", &_phiVal, "phi/F");
        _alignTree->Branch("eta", &_etaVal, "eta/F");
        _alignTree->Branch("alpha", &_alphaVal, "alpha/F");
        _alignTree->Branch("beta", &_betaVal, "beta/F");
        _alignTree->Branch("gamma", &_gammaVal, "gamma/F");
        _alignTree->Branch("dx", &_dxVal, "dx/F");
        _alignTree->Branch("dy", &_dyVal, "dy/F");
        _alignTree->Branch("dz", &_dzVal, "dz/F");
        _alignTree->Branch("dr", &_drVal, "dr/F");
        _alignTree->Branch("dphi", &_dphiVal, "dphi/F");
        _alignTree->Branch("dalpha", &_dalphaVal, "dalpha/F");
        _alignTree->Branch("dbeta", &_dbetaVal, "dbeta/F");
        _alignTree->Branch("dgamma", &_dgammaVal, "dgamma/F");
        _alignTree->Branch("ldx", &_ldxVal, "ldx/F");
        _alignTree->Branch("ldy", &_ldyVal, "ldy/F");
        _alignTree->Branch("ldz", &_ldzVal, "ldz/F");
        _alignTree->Branch("ldr", &_ldrVal, "ldr/F");
        _alignTree->Branch("ldphi", &_ldphiVal, "ldphi/F");
        _alignTree->Branch("useDetId", &_useDetId, "useDetId/I");
        _alignTree->Branch("detDim", &_detDim, "detDim/I");
        _alignTree->Branch("rotx",&_rotxVal, "rotx/F"); 
        _alignTree->Branch("roty",&_rotyVal, "roty/F"); 
        _alignTree->Branch("rotz",&_rotzVal, "rotz/F"); 
        _alignTree->Branch("drotx",&_drotxVal, "drotx/F");      
        _alignTree->Branch("droty",&_drotyVal, "droty/F");      
        _alignTree->Branch("drotz",&_drotzVal, "drotz/F");      
        _alignTree->Branch("surW", &_surWidth, "surW/F");
        _alignTree->Branch("surL", &_surLength, "surL/F");
        _alignTree->Branch("surRot", &_surRot, "surRot[9]/D");

        _mgacollection.clear(); 
}

Member Function Documentation

void MuonGeometryArrange::analyze	(	const edm::Event &	,
		const edm::EventSetup &	iSetup
	)		`[virtual]`

Implements edm::EDAnalyzer.

Definition at line 457 of file MuonGeometryArrange.cc.

References _alignTree, _inputXMLCurrent, _inputXMLReference, _theFile, compare(), endHist(), MuonAlignment::fillGapsInSurvey(), firstEvent_, MuonAlignment::getAlignableMuon(), inputAlign1, inputAlign2, inputAlign2a, inputGeometry1, and inputGeometry2.

                                                              {
  if (firstEvent_) {
    
    // My stuff
    MuonAlignmentInputXML inputMethod1(_inputXMLCurrent);
    inputAlign1 = new MuonAlignment(iSetup, inputMethod1);
    inputAlign1->fillGapsInSurvey(0, 0);
    MuonAlignmentInputXML inputMethod2(_inputXMLReference);
    inputAlign2 = new MuonAlignment(iSetup, inputMethod2);
    inputAlign2->fillGapsInSurvey(0, 0);
    MuonAlignmentInputXML inputMethod3(_inputXMLReference);
    inputAlign2a = new MuonAlignment(iSetup, inputMethod3);
    inputAlign2a->fillGapsInSurvey(0, 0);
    
    inputGeometry1 = static_cast<Alignable*> (inputAlign1->getAlignableMuon());
    inputGeometry2 = static_cast<Alignable*> (inputAlign2->getAlignableMuon());
    Alignable* inputGeometry2Copy2 = 
      static_cast<Alignable*> (inputAlign2a->getAlignableMuon());
    
    //compare the goemetries
    compare(inputGeometry1, inputGeometry2, inputGeometry2Copy2);
    
    //write out ntuple
    //might be better to do within output module
    _theFile->cd();
    _alignTree->Write();
    endHist();
    //   _theFile->Close();
    
    firstEvent_ = false;
  }
}

void MuonGeometryArrange::beginJob ( void ) [virtual]

Read from DB and print survey info.

Reimplemented from edm::EDAnalyzer.

Definition at line 450 of file MuonGeometryArrange.cc.

References firstEvent_.

                                  {
  firstEvent_ = true;
}

bool MuonGeometryArrange::checkChosen ( Alignable * ali ) [private]

Definition at line 896 of file MuonGeometryArrange.cc.

References _endcap, _ring, _station, gather_cfg::cout, CSC(), DetId::det(), Alignable::geomDetId(), Alignable::id(), DetId::Muon, and DetId::subdetId().

Referenced by isMother(), and passChosen().

                                                   {
 // Check whether the item passed satisfies the criteria given.
  if(ali==0x0) return false;    // elementary sanity
 // Is this in the CSC section?  If not, bail.  Later may extend.
  if(ali->geomDetId().det()!=DetId::Muon ||
     ali->geomDetId().subdetId()!=MuonSubdetId::CSC) return false;
 // If it is a CSC alignable, then check that the station, etc are
 // those requested.
 // One might think of aligning more than a single ring at a time,
 // by using a vector of ring numbers.  I don't see the sense in
 // trying to align more than one station at a time for comparison.
  CSCDetId cscId(ali->geomDetId());
#ifdef jnbdebug
std::cout<<"JNB "<<ali->id()<<" "<<cscId.endcap()<<" "
<<cscId.station()<<" "<<cscId.ring()<<" "<<cscId.chamber()<<"   "
<<_endcap<<" "<<_station<<" "<<_ring
<<"\n"<<std::flush;
#endif
  if(cscId.endcap()==_endcap && cscId.station()==_station &&
     cscId.ring()==_ring) {
      return true;
  }
  return false;
}

void MuonGeometryArrange::compare	(	Alignable *	refAli,
		Alignable *	curAli,
		Alignable *	curAliCopy2
	)		`[private]`

Definition at line 492 of file MuonGeometryArrange.cc.

References compareGeometries(), Alignable::components(), and i.

Referenced by analyze().

                                       {

 // First sanity
  if(refAli==0x0){return;}      
  if(curAli==0x0){return;}
  
  const std::vector<Alignable*>& refComp = refAli->components();
  const std::vector<Alignable*>& curComp = curAli->components();
  const std::vector<Alignable*>& curComp2 = curAliCopy2->components();
  compareGeometries(refAli, curAli, curAliCopy2);

  int nComp=refComp.size();
  for(int i=0; i<nComp; i++){
    compare(refComp[i], curComp[i], curComp2[i]);
  }
  return;
}

void MuonGeometryArrange::compareGeometries	(	Alignable *	refAli,
		Alignable *	curAli,
		Alignable *	curAliCopy2
	)		`[private]`

Definition at line 512 of file MuonGeometryArrange.cc.

References _weightBy, _weightById, _weightByIdVector, align::centerOfMass(), change(), CastorDataFrameFilter_impl::check(), Alignable::components(), align::createPoints(), cond::rpcobgas::detid, diffTreeTool::diff, align::diffAlignables(), align::diffRot(), alignCSCRings::e, fillTree(), i, isMother(), PV3DBase< T, PVType, FrameType >::mag(), Alignable::mother(), align::moveAlignable(), align::readModuleList(), makeMuonMisalignmentScenario::rot, align::toAngles(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by compare().

                                                      {
 // First sanity
  if(refAli==0x0){return;}      
  if(curAli==0x0){return;}
 // Is this the Ring we want to align?  If so it will contain the
 // chambers specified in the configuration file        
  if(!isMother(refAli)) return; // Not the desired alignable object
 // But... There are granddaughters involved--and I don't want to monkey with
 // the layers of the chambers.  So, if the mother of this is also an approved
 // mother, bail.
  if(isMother(refAli->mother() )) return;
  const std::vector<Alignable*>& refComp = refAli->components();
  const std::vector<Alignable*>& curComp = curCopy->components();
  if(refComp.size()!=curComp.size()){
     return;
  }
 // GlobalVectors is a vector of GlobalVector which is a 3D vector
  align::GlobalVectors originalVectors;
  align::GlobalVectors currentVectors;
  align::GlobalVectors originalRelativeVectors;
  align::GlobalVectors currentRelativeVectors;

        
  int nComp = refComp.size();
  int nUsed = 0;
 // Use the total displacements here:
  CLHEP::Hep3Vector TotalX, TotalL;
  TotalX.set(0.,0.,0.);   TotalL.set(0., 0., 0.);
//  CLHEP::Hep3Vector* Rsubtotal, Wsubtotal, DRsubtotal, DWsubtotal;
  std::vector<CLHEP::Hep3Vector> Positions;
  std::vector<CLHEP::Hep3Vector> DelPositions;

  double xrcenter=0.;
  double yrcenter=0.;
  double zrcenter=0.;
  double xccenter=0.;
  double yccenter=0.;
  double zccenter=0.;

  bool useIt;
 // Create the "center" for the reference alignment chambers, and
 // load a vector of their centers
  for(int ich=0; ich<nComp; ich++){
    useIt=true;
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         useIt=false;
    }
    if(!useIt) continue;
    align::GlobalVectors curVs;
    align::createPoints(&curVs, refComp[ich], 
      _weightBy, _weightById, _weightByIdVector);
    align::GlobalVector pointsCM = align::centerOfMass(curVs);
    originalVectors.push_back(pointsCM);
    nUsed++;
    xrcenter+= pointsCM.x();
    yrcenter+= pointsCM.y();
    zrcenter+= pointsCM.z();
  }
  xrcenter=xrcenter/nUsed;
  yrcenter=yrcenter/nUsed;
  zrcenter=zrcenter/nUsed;

 // Create the "center" for the current alignment chambers, and
 // load a vector of their centers
   for(int ich=0; ich<nComp; ich++){
    useIt=true;
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         useIt=false;
    }
    if(!useIt)continue;
    align::GlobalVectors curVs;
    align::createPoints(&curVs, curComp[ich], 
      _weightBy, _weightById, _weightByIdVector);
    align::GlobalVector pointsCM = align::centerOfMass(curVs);
    currentVectors.push_back(pointsCM);

    xccenter+= pointsCM.x();
    yccenter+= pointsCM.y();
    zccenter+= pointsCM.z();
  }
  xccenter=xccenter/nUsed;
  yccenter=yccenter/nUsed;
  zccenter=zccenter/nUsed;


 // OK, now load the <very approximate> vectors from the ring "centers"
  align::GlobalVector CCur(xccenter, yccenter, zccenter);
  align::GlobalVector CRef(xrcenter, yrcenter, zrcenter);
  int nCompR = currentVectors.size();
  for(int ich=0; ich<nCompR; ich++){
    originalRelativeVectors.push_back(originalVectors[ich]-CRef);
    currentRelativeVectors.push_back(currentVectors[ich]-CCur);
  }

 // All right.  Now let the hacking begin.
 // First out of the gate let's try using the raw values and see what
 // diffRot does for us.

 
  align::RotationType rtype3=align::diffRot(currentRelativeVectors,
                                            originalRelativeVectors);


  align::EulerAngles angles(3);
  angles = align::toAngles(rtype3);

  for(int ich=0; ich<nComp; ich++){
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         continue;
    }
      CLHEP::Hep3Vector Rtotal, Wtotal;
      Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
      for (int i = 0; i < 100; i++){
        AlgebraicVector diff = align::diffAlignables(refComp[ich],curComp[ich],
                           _weightBy, _weightById, _weightByIdVector);
        CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
        Rtotal+=dR;
        CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
        CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
        CLHEP::HepRotation drot(dW.unit(),dW.mag());
        rot*=drot;
        Wtotal.set(rot.axis().x()*rot.delta(),
                   rot.axis().y()*rot.delta(), rot.axis().z()*rot.delta());
        align::moveAlignable(curComp[ich], diff);
        float tolerance = 1e-7;
        AlgebraicVector check = align::diffAlignables(refComp[ich],curComp[ich],
                        _weightBy, _weightById, _weightByIdVector);
        align::GlobalVector checkR(check[0],check[1],check[2]);
        align::GlobalVector checkW(check[3],check[4],check[5]);
        DetId detid(refComp[ich]->id());
        if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){
//       edm::LogInfo("CompareGeoms") << "Tolerance Exceeded!(alObjId: " 
//       << refAli->alignableObjectId()
//       << ", rawId: " << refComp[ich]->geomDetId().rawId()
//       << ", subdetId: "<< detid.subdetId() << "): " << diff;
        }
        else{
         TotalX+=Rtotal;
         break;
        }         // end of else
      }  // end of for on int i
  }     // end of for on ich

  // At this point we should have a total displacement and total L
   TotalX=TotalX/nUsed;

  // Now start again!
   AlgebraicVector change(6);
   change(1)=TotalX.x();
   change(2)=TotalX.y();
   change(3)=TotalX.z();

   change(4)=angles[0];
   change(5)=angles[1];
   change(6)=angles[2];
   align::moveAlignable(curAli, change);        // move as a chunk

  // Now get the components again.  They should be in new locations
   const std::vector<Alignable*>& curComp2 = curAli->components();

  for(int ich=0; ich<nComp; ich++){
     CLHEP::Hep3Vector Rtotal, Wtotal;
     Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
     if(_weightById){
       if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         continue;
     }
    
     for (int i = 0; i < 100; i++){
       AlgebraicVector diff = align::diffAlignables(refComp[ich],curComp2[ich],
                          _weightBy, _weightById, _weightByIdVector);
       CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
       Rtotal+=dR;
       CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
       CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
       CLHEP::HepRotation drot(dW.unit(),dW.mag());
       rot*=drot;
       Wtotal.set(rot.axis().x()*rot.delta(), rot.axis().y()*rot.delta(),
                  rot.axis().z()*rot.delta());
       align::moveAlignable(curComp2[ich], diff);
       float tolerance = 1e-7;
       AlgebraicVector check = align::diffAlignables(refComp[ich],curComp2[ich],
                       _weightBy, _weightById, _weightByIdVector);
       align::GlobalVector checkR(check[0],check[1],check[2]);
       align::GlobalVector checkW(check[3],check[4],check[5]);
       if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){}
       else{break;}
     }   // end of for on int i
     AlgebraicVector TRtot(6);
     TRtot(1) = Rtotal.x(); TRtot(2) = Rtotal.y(); TRtot(3) = Rtotal.z();
     TRtot(4) = Wtotal.x(); TRtot(5) = Wtotal.y(); TRtot(6) = Wtotal.z();
     fillTree(refComp[ich], TRtot);
  }    // end of for on ich



}

void MuonGeometryArrange::createROOTGeometry ( const edm::EventSetup & iSetup ) [private]

Definition at line 455 of file MuonGeometryArrange.cc.

{}

void MuonGeometryArrange::endHist ( ) [private]

Definition at line 166 of file MuonGeometryArrange.cc.

References _mgacollection, _theFile, i, makeGraph(), findQualityFiles::maxI, findQualityFiles::minI, run_regression::ret, pileupReCalc_HLTpaths::scale, findQualityFiles::size, mathSSE::sqrt(), x, and detailsBasic3DVector::y.

Referenced by analyze().

                                 {
 // Unpack the list and create ntuples here.

   int size=_mgacollection.size();
   if(size<=0) return;  // nothing to do here.
   float* xp = new float[size+1];
   float* yp = new float[size+1];
   int i;
   float minV, maxV;
   int minI, maxI;

   minV=99999999.; maxV=-minV;  minI=9999999; maxI=-minI;
   TGraph* grx=0x0;
   TH2F* dxh=0x0;

// for position plots:
   for(i=0; i<size; i++){
     if(_mgacollection[i].phipos<minI) minI=_mgacollection[i].phipos;
     if(_mgacollection[i].phipos>maxI) maxI=_mgacollection[i].phipos;
     xp[i]=_mgacollection[i].phipos;
   }
   if(minI>=maxI) return;       // can't do anything?
   xp[size]=xp[size-1]+1;       // wraparound point

   if(1<minI) minI=1;
   if(size>maxI) maxI=size;
   maxI++;      // allow for wraparound to show neighbors
   int sizeI=maxI+1-minI;
   float smi=minI-1;
   float sma=maxI+1;
 

// Dx plot

   for(i=0; i<size; i++){
     if(_mgacollection[i].ldx<minV) minV=_mgacollection[i].ldx;
     if(_mgacollection[i].ldx>maxV) maxV=_mgacollection[i].ldx;
     yp[i]=_mgacollection[i].ldx;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delX_vs_position", "Local #delta X vs position", 
       "GdelX_vs_position","#delta x in cm", xp, yp, size);
// Dy plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].ldy<minV) minV=_mgacollection[i].ldy;
     if(_mgacollection[i].ldy>maxV) maxV=_mgacollection[i].ldy;
     yp[i]=_mgacollection[i].ldy;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delY_vs_position", "Local #delta Y vs position", 
       "GdelY_vs_position","#delta y in cm", xp, yp, size);

// Dz plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dz<minV) minV=_mgacollection[i].dz;
     if(_mgacollection[i].dz>maxV) maxV=_mgacollection[i].dz;
     yp[i]=_mgacollection[i].dz;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delZ_vs_position", "Local #delta Z vs position", 
       "GdelZ_vs_position","#delta z in cm", xp, yp, size);

// Dphi plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dphi<minV) minV=_mgacollection[i].dphi;
     if(_mgacollection[i].dphi>maxV) maxV=_mgacollection[i].dphi;
     yp[i]=_mgacollection[i].dphi;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delphi_vs_position", "#delta #phi vs position", 
       "Gdelphi_vs_position","#delta #phi in radians", xp, yp, size);

// Dr plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dr<minV) minV=_mgacollection[i].dr;
     if(_mgacollection[i].dr>maxV) maxV=_mgacollection[i].dr;
     yp[i]=_mgacollection[i].dr;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delR_vs_position", "#delta R vs position", 
       "GdelR_vs_position","#delta R in cm", xp, yp, size);

// Drphi plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     float ttemp=_mgacollection[i].r*_mgacollection[i].dphi;
     if(ttemp<minV) minV=ttemp;
     if(ttemp>maxV) maxV=ttemp;
     yp[i]=ttemp;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delRphi_vs_position", "R #delta #phi vs position", 
       "GdelRphi_vs_position","R #delta #phi in cm", xp, yp, size);

// Dalpha plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dalpha<minV) minV=_mgacollection[i].dalpha;
     if(_mgacollection[i].dalpha>maxV) maxV=_mgacollection[i].dalpha;
     yp[i]=_mgacollection[i].dalpha;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delalpha_vs_position", "#delta #alpha vs position", 
       "Gdelalpha_vs_position","#delta #alpha in rad", xp, yp, size);

// Dbeta plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dbeta<minV) minV=_mgacollection[i].dbeta;
     if(_mgacollection[i].dbeta>maxV) maxV=_mgacollection[i].dbeta;
     yp[i]=_mgacollection[i].dbeta;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delbeta_vs_position", "#delta #beta vs position", 
       "Gdelbeta_vs_position","#delta #beta in rad", xp, yp, size);

// Dgamma plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dgamma<minV) minV=_mgacollection[i].dgamma;
     if(_mgacollection[i].dgamma>maxV) maxV=_mgacollection[i].dgamma;
     yp[i]=_mgacollection[i].dgamma;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delgamma_vs_position", "#delta #gamma vs position", 
       "Gdelgamma_vs_position","#delta #gamma in rad", xp, yp, size);

// Drotx plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].drotx<minV) minV=_mgacollection[i].drotx;
     if(_mgacollection[i].drotx>maxV) maxV=_mgacollection[i].drotx;
     yp[i]=_mgacollection[i].drotx;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotX_vs_position", "#delta rotX vs position", 
       "GdelrotX_vs_position","#delta rotX in rad", xp, yp, size);

// Droty plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].droty<minV) minV=_mgacollection[i].droty;
     if(_mgacollection[i].droty>maxV) maxV=_mgacollection[i].droty;
     yp[i]=_mgacollection[i].droty;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotY_vs_position", "#delta rotY vs position", 
       "GdelrotY_vs_position","#delta rotY in rad", xp, yp, size);

// Drotz plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].drotz<minV) minV=_mgacollection[i].drotz;
     if(_mgacollection[i].drotz>maxV) maxV=_mgacollection[i].drotz;
     yp[i]=_mgacollection[i].drotz;
   }
   yp[size]=yp[0];      // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotZ_vs_position", "#delta rotZ vs position", 
       "GdelrotZ_vs_position","#delta rotZ in rad", xp, yp, size);



// Vector plots
// First find the maximum length of sqrt(dx*dx+dy*dy):  we'll have to
// scale these for visibility
   maxV=-99999999.;
   float ttemp, rtemp;
   float maxR=-9999999.;
   for(i=0; i<size; i++){
     ttemp= sqrt(_mgacollection[i].dx*_mgacollection[i].dx+
        _mgacollection[i].dy*_mgacollection[i].dy);
     rtemp= sqrt(_mgacollection[i].x*_mgacollection[i].x+
        _mgacollection[i].y*_mgacollection[i].y);
     if(ttemp>maxV) maxV=ttemp;
     if(rtemp>maxR) maxR=rtemp;
   }
   
  // Don't try to scale rediculously small values 
   float smallestVcm=.001;      // 10 microns
   if(maxV<smallestVcm) maxV=smallestVcm;
   float scale=0.;
   float lside=1.1*maxR;
   if(lside<=0) lside=100.;
   if(maxV>0){scale=.09*lside/maxV;}    // units of pad length!
   char scalename[50];
   int ret=snprintf(scalename,50,"#delta #bar{x}   length =%f cm",maxV);
  // If ret<=0 we don't want to print the scale!
  
   if(ret>0){
     dxh=new TH2F("vecdrplot",scalename,80,-lside,lside,80,-lside,lside);
   }
   else{
     dxh=new TH2F("vecdrplot","delta #bar{x} Bad scale",80,-lside,lside,80,-lside,lside);
   }
   dxh->GetXaxis()->SetTitle("x in cm");
   dxh->GetYaxis()->SetTitle("y in cm");
   dxh->SetStats(kFALSE);
   dxh->Draw();
   TArrow* arrow;
   for(i=0; i<size; i++){
     ttemp= sqrt(_mgacollection[i].dx*_mgacollection[i].dx+
       _mgacollection[i].dy*_mgacollection[i].dy);
//     ttemp=ttemp*scale;
     float nx=_mgacollection[i].x+scale*_mgacollection[i].dx;
     float ny=_mgacollection[i].y+scale*_mgacollection[i].dy;
     arrow = new TArrow(_mgacollection[i].x, 
                   _mgacollection[i].y, nx, ny);// ttemp*.3*.05, "->");
     arrow->SetLineWidth(2); arrow->SetArrowSize(ttemp*.2*.05/maxV);
     arrow->SetLineColor(1); arrow->SetLineStyle(1);
     arrow->Paint();
     dxh->GetListOfFunctions()->Add(static_cast<TObject*>(arrow));
//     arrow->Draw();
//     arrow->Write();
   }
   dxh->Write();
 
   _theFile->Write();   
   _theFile->Close();

   delete [] yp;  delete [] xp;

}

void MuonGeometryArrange::fillTree	(	Alignable *	refAli,
		AlgebraicVector	diff
	)		`[private]`

Definition at line 716 of file MuonGeometryArrange.cc.

Referenced by compareGeometries().

                                                                         {
        
        
        _id = refAli->id();
        _level = refAli->alignableObjectId();
        //need if ali has no mother
        if (refAli->mother()){
                _mid = refAli->mother()->geomDetId().rawId();
                _mlevel = refAli->mother()->alignableObjectId();
        }
        else{
                _mid = -1;
                _mlevel = -1;
        }
        DetId detid(_id);
        _sublevel = detid.subdetId();
        int ringPhiPos=-99;
        if(detid.det()==DetId::Muon && detid.subdetId()== MuonSubdetId::CSC){
           CSCDetId cscId(refAli->geomDetId());
           ringPhiPos = cscId.chamber();
        }
        _xVal = refAli->globalPosition().x();
        _yVal = refAli->globalPosition().y();
        _zVal = refAli->globalPosition().z();
        align::GlobalVector vec(_xVal,_yVal,_zVal);
        _rVal = vec.perp();
        _phiVal = vec.phi();
        _etaVal = vec.eta();
        align::RotationType rot = refAli->globalRotation();
        align::EulerAngles eulerAngles = align::toAngles(rot);
        _rotxVal = atan2(rot.yz(), rot.zz());
        float ttt=-rot.xz();
        if(ttt>1.) ttt=1.;
        if(ttt<-1.) ttt=-1.;
        _rotyVal = asin(ttt);
        _rotzVal = atan2(rot.xy(), rot.xx());
        _alphaVal = eulerAngles[0];
        _betaVal = eulerAngles[1];
        _gammaVal = eulerAngles[2];
        _dxVal = diff[0];
        _dyVal = diff[1];
        _dzVal = diff[2];
        //getting dR and dPhi
        align::GlobalVector vRef(_xVal,_yVal,_zVal);
        align::GlobalVector vCur(_xVal - _dxVal, _yVal-_dyVal, _zVal - _dzVal);
        _drVal = vCur.perp() - vRef.perp();
        _dphiVal = vCur.phi() - vRef.phi();
        
        _dalphaVal = diff[3];
        _dbetaVal = diff[4];
        _dgammaVal = diff[5];
        _drotxVal=-999.; _drotyVal=-999.; _drotzVal=-999.; 

        align::EulerAngles deuler(3);
        deuler(1)=_dalphaVal;
        deuler(2)= _dbetaVal;
        deuler(3)= _dgammaVal;
        align::RotationType drot = align::toMatrix(deuler);
        double xx=rot.xx();
        double xy=rot.xy();
        double xz=rot.xz();
        double yx=rot.yx();
        double yy=rot.yy();
        double yz=rot.yz();
        double zx=rot.zx();
        double zy=rot.zy();
        double zz=rot.zz();
        double detrot=(zz*yy - zy*yz)*xx + (-zz*yx + zx*yz)*xy + (zy*yx - zx*yy)*xz;
        detrot=1/detrot;
        double ixx=(zz*yy - zy*yz)*detrot;
        double ixy=(-zz*xy + zy*xz)*detrot;
        double ixz=(yz*xy - yy*xz)*detrot;
        double iyx=(-zz*yx + zx*yz)*detrot;
        double iyy=(zz*xx - zx*xz)*detrot;
        double iyz=(-yz*xx + yx*xz)*detrot;
        double izx=(zy*yx - zx*yy)*detrot;
        double izy=(-zy*xx + zx*xy)*detrot;
        double izz=(yy*xx - yx*xy)*detrot;
        align::RotationType invrot(ixx,ixy,ixz, iyx,iyy,iyz, izx,izy,izz);
        align::RotationType prot = rot*drot*invrot;
//      align::RotationType prot = rot*drot;
        float protx; //, proty, protz;
        protx = atan2(prot.yz(), prot.zz());
        _drotxVal = protx;//_rotxVal-protx; //atan2(drot.yz(), drot.zz());
        ttt=-prot.xz();
        if(ttt>1.) ttt=1.;
        if(ttt<-1.) ttt=-1.;
        _drotyVal = asin(ttt);// -_rotyVal;
        _drotzVal = atan2(prot.xy(), prot.xx());// - _rotzVal;
// Above does not account for 2Pi wraparounds!
// Prior knowledge:  these are supposed to be small rotations.  Therefore:
        if(_drotxVal>3.141592656) _drotxVal=-6.2831853072+_drotxVal;
        if(_drotxVal<-3.141592656) _drotxVal=6.2831853072+_drotxVal;
        if(_drotyVal>3.141592656) _drotyVal=-6.2831853072+_drotyVal;
        if(_drotyVal<-3.141592656) _drotyVal=6.2831853072+_drotyVal;
        if(_drotzVal>3.141592656) _drotzVal=-6.2831853072+_drotzVal;
        if(_drotzVal<-3.141592656) _drotzVal=6.2831853072+_drotzVal;
        
        _ldxVal=-999.; _ldyVal=-999.; _ldxVal=-999.; 
        _ldrVal=-999.; _ldphiVal=-999; // set fake

//      if(refAli->alignableObjectId() == align::AlignableDetUnit){
         align::GlobalVector dV(_dxVal, _dyVal, _dzVal); 
         align::LocalVector pointL = refAli->surface().toLocal(dV);
         //align::LocalVector pointL = (refAli->mother())->surface().toLocal(dV);
         _ldxVal=pointL.x(); _ldyVal=pointL.y(); _ldzVal=pointL.z();
         _ldphiVal=pointL.phi(); _ldrVal=pointL.perp();
//      }
        //detIdFlag
        if (refAli->alignableObjectId() == align::AlignableDetUnit){
          if (_detIdFlag){
           if ((passIdCut(refAli->id()))||(passIdCut(refAli->mother()->id()))){
                _useDetId = 1;
           }
           else{
                _useDetId = 0;
           }
          }
        }
        // det module dimension
        if (refAli->alignableObjectId() == align::AlignableDetUnit){
           if (refAli->mother()->alignableObjectId() != align::AlignableDet){
             _detDim = 1;}
           else if (refAli->mother()->alignableObjectId() == 
                                       align::AlignableDet) {_detDim = 2;}
        }
        else _detDim = 0;
        
        
        
        _surWidth = refAli->surface().width();
        _surLength = refAli->surface().length();
        align::RotationType rt = refAli->globalRotation();
        _surRot[0] = rt.xx(); _surRot[1] = rt.xy(); _surRot[2] = rt.xz();
        _surRot[3] = rt.yx(); _surRot[4] = rt.yy(); _surRot[5] = rt.yz();
        _surRot[6] = rt.zx(); _surRot[7] = rt.zy(); _surRot[8] = rt.zz();

        MGACollection holdit;
        holdit.id=_id;  holdit.level=_level;  holdit.mid=_mid;
        holdit.mlevel=_mlevel;
        holdit.sublevel=_sublevel;
        holdit.x=_xVal;  holdit.y=_yVal;  holdit.z=_zVal;       
        holdit.r=_rVal;  holdit.phi=_phiVal;  holdit.eta=_etaVal;
        holdit.alpha=_alphaVal;  holdit.beta=_betaVal;  holdit.gamma=_gammaVal;
        holdit.dx=_dxVal;  holdit.dy=_dyVal;  holdit.dz=_dzVal; 
        holdit.dr=_drVal;  holdit.dphi=_dphiVal; 
        holdit.dalpha=_dalphaVal;  holdit.dbeta=_dbetaVal;  
        holdit.dgamma=_dgammaVal;
        holdit.useDetId=_useDetId; holdit.detDim=_detDim;
        holdit.surW=_surWidth; holdit.surL=_surLength;
        holdit.ldx=_ldxVal; holdit.ldy=_ldyVal; holdit.ldz=_ldzVal;
        holdit.ldr=_ldrVal; holdit.ldphi=_ldphiVal;
        holdit.rotx=_rotxVal; holdit.roty=_rotyVal; holdit.rotz=_rotzVal;
        holdit.drotx=_drotxVal; holdit.droty=_drotyVal; holdit.drotz=_drotzVal;
        for(int i=0; i<9; i++){holdit.surRot[i]=_surRot[i];}
        holdit.phipos=ringPhiPos;
        _mgacollection.push_back(holdit);

        
        //Fill
        _alignTree->Fill();
        
}

bool MuonGeometryArrange::isMother ( Alignable * ali ) [private]

Definition at line 881 of file MuonGeometryArrange.cc.

References checkChosen(), Alignable::components(), i, and findQualityFiles::size.

Referenced by compareGeometries().

                                                {
  // Is this the mother ring?
 if(ali==0x0) return false;     // elementary sanity
 const std::vector<Alignable*>& aliComp = ali->components();

 int size=aliComp.size();
 if(size<=0) return false;      // no subcomponents
 
 for(int i=0; i<size; i++){
   if(checkChosen(aliComp[i])) return true;     // A ring has CSC chambers
 }                                      // as subcomponents
 return false;          // 1'st layer of subcomponents weren't CSC chambers
}

void MuonGeometryArrange::makeGraph	(	int	sizeI,
		float	smi,
		float	sma,
		float	minV,
		float	maxV,
		TH2F *	dxh,
		TGraph *	grx,
		const char *	name,
		const char *	title,
		const char *	titleg,
		const char *	axis,
		float *	xp,
		float *	yp,
		int	numEntries
	)		`[private]`

Definition at line 417 of file MuonGeometryArrange.cc.

References diffTreeTool::diff.

Referenced by endHist().

                                       {

  if(minV>=maxV || smi>=sma || sizeI<=1 || xp==0x0 || yp==0x0) return;
        // out of bounds, bail
  float diff=maxV-minV;
  float over=.05*diff;
  double ylo=minV-over;
  double yhi=maxV+over;
  double dsmi, dsma;
  dsmi=smi; dsma=sma;
  dxh= new TH2F(name, title,
    sizeI+2, dsmi, dsma, 50, ylo, yhi);
  dxh->GetXaxis()->SetTitle("Position around ring");
  dxh->GetYaxis()->SetTitle(axis);
  dxh->SetStats(kFALSE);
  dxh->Draw();
  grx = new TGraph(size, xp, yp);
  grx->SetName(titleg);
  grx->SetTitle(title);
  grx->SetMarkerColor(2); grx->SetMarkerStyle(3);
  grx->GetXaxis()->SetLimits(dsmi, dsma);
  grx->GetXaxis()->SetTitle("position number");
  grx->GetYaxis()->SetLimits(ylo,yhi);
  grx->GetYaxis()->SetTitle(axis);
  grx->Draw("A*");
  grx->Write();
  return; 
}

bool MuonGeometryArrange::passChosen ( Alignable * ali ) [private]

Definition at line 922 of file MuonGeometryArrange.cc.

References checkChosen(), Alignable::components(), i, and findQualityFiles::size.

                                                    {

 // Check to see if this contains CSC components of the appropriate ring
 // Ring will contain N Alignables which represent chambers, each of which
 // in turn contains M planes.  For our purposes we don't care about the
 // planes.
 // Hmm.  Interesting question:  Do I want to try to fit the chamber as
 // such, or use the geometry?
 // I want to fit the chamber, so I'll try to use its presence as the marker.
 // What specifically identifies a chamber as a chamber, and not as a layer?
 // The fact that it has layers as sub components, or the fact that it is
 // the first item with a non-zero ID breakdown?  Pick the latter.
 //
 if(ali==0x0) return false;
 if(checkChosen(ali)) return true;      // If this is one of the desired
                                        // CSC chambers, accept it
 const std::vector<Alignable*>& aliComp = ali->components();

 int size=aliComp.size();
 if(size<=0) return false;      // no subcomponents
 
 for(int i=0; i<size; i++){
   if(checkChosen(aliComp[i])) return true;     // A ring has CSC chambers
 }                                      // as subcomponents
 return false;          // 1'st layer of subcomponents weren't CSC chambers
}

bool MuonGeometryArrange::passIdCut ( uint32_t id ) [private]

Definition at line 949 of file MuonGeometryArrange.cc.

References _detIdFlagVector, cond::rpcobgas::detid, and i.

Referenced by fillTree().

                                                {
        
        bool pass = false;
        DetId detid(id);
//      if(detid.det()==DetId::Muon && detid.subdetId()== MuonSubdetId::CSC){
//         CSCDetId cscId(refAli->geomDetId());
//         if(cscId.layer()!=1) return false;           // ONLY FIRST LAYER!
//      }
        int nEntries = _detIdFlagVector.size();
        
        for (int i = 0; i < nEntries; i++){
                if (_detIdFlagVector[i] == id) pass = true;
        }
        
        return pass;
        
}