#include <Validation/RecoVertex/src/TrackParameterAnalyzer.cc>

Inheritance diagram for TrackParameterAnalyzer:

Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)

virtual void	beginJob ()

virtual void	endJob ()

	TrackParameterAnalyzer (const edm::ParameterSet &)

	~TrackParameterAnalyzer ()

Public Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzer ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Member Functions
bool	match (const ParameterVector &a, const ParameterVector &b)

Private Attributes
edm::EDGetTokenT < edm::SimTrackContainer >	edmSimTrackContainerToken_

edm::EDGetTokenT < edm::SimVertexContainer >	edmSimVertexContainerToken_

TH1 *	h1_Beff_

TH1 *	h1_par0_

TH1 *	h1_par1_

TH1 *	h1_par2_

TH1 *	h1_par3_

TH1 *	h1_par4_

TH1 *	h1_pull0_

TH1 *	h1_pull1_

TH1 *	h1_pull2_

TH1 *	h1_pull3_

TH1 *	h1_pull4_

TH1 *	h1_res0_

TH1 *	h1_res1_

TH1 *	h1_res2_

TH1 *	h1_res3_

TH1 *	h1_res4_

TH2 *	h2_dvsphi_

std::string	outputFile_

edm::EDGetTokenT < reco::TrackCollection >	recoTrackCollectionToken_

TFile *	rootFile_

double	simUnit_

bool	verbose_

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 47 of file TrackParameterAnalyzer.h.

Constructor & Destructor Documentation

TrackParameterAnalyzer::TrackParameterAnalyzer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 36 of file TrackParameterAnalyzer.cc.

References edm::getReleaseVersion(), outputFile_, rootFile_, simUnit_, and AlCaHLTBitMon_QueryRunRegistry::string.

   : edmSimVertexContainerToken_( consumes<edm::SimVertexContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
   , edmSimTrackContainerToken_( consumes<edm::SimTrackContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
   , recoTrackCollectionToken_( consumes<reco::TrackCollection>( edm::InputTag( iConfig.getUntrackedParameter<std::string>( "recoTrackProducer" ) ) ) )
   , outputFile_( iConfig.getUntrackedParameter<std::string>( "outputFile" ) )
   , simUnit_( 1.0 ) //  starting from  CMSSW_1_2_x, I think
   , verbose_( iConfig.getUntrackedParameter<bool>( "verbose", false ) ) {
    //now do whatever initialization is needed
   // open output file to store histograms}
   TString tversion(edm::getReleaseVersion());
   tversion = tversion.Remove(0,1);
   tversion = tversion.Remove(tversion.Length()-1,tversion.Length());
   outputFile_  = std::string(tversion)+"_"+outputFile_;
   rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE"); 
   if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
     simUnit_=0.1;  // for use in  CMSSW_1_1_1 tutorial
   }
 }

TrackParameterAnalyzer::~TrackParameterAnalyzer ( )

Definition at line 56 of file TrackParameterAnalyzer.cc.

References rootFile_.

 {
    // do anything here that needs to be done at destruction time
    // (e.g. close files, deallocate resources etc.)
   delete rootFile_;
 }

Member Function Documentation

void TrackParameterAnalyzer::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

virtual

Implements edm::EDAnalyzer.

Definition at line 128 of file TrackParameterAnalyzer.cc.

References EnergyCorrector::c, funct::cos(), gather_cfg::cout, edmSimTrackContainerToken_, edmSimVertexContainerToken_, edm::Event::getByToken(), h1_Beff_, h1_par0_, h1_par1_, h1_par2_, h1_par3_, h1_par4_, h1_pull0_, h1_pull1_, h1_pull2_, h1_pull3_, h1_pull4_, h1_res0_, h1_res1_, h1_res2_, h1_res3_, h1_res4_, h2_dvsphi_, reco::TrackBase::i_dsz, reco::TrackBase::i_dxy, reco::TrackBase::i_lambda, reco::TrackBase::i_phi, reco::TrackBase::i_qoverp, kappa, M_PI, match(), AlCaHLTBitMon_ParallelJobs::p, funct::pow(), class-composition::Q, lumiQueryAPI::q, recoTrackCollectionToken_, alignCSCRings::s, simUnit_, funct::sin(), mathSSE::sqrt(), lumiQTWidget::t, findQualityFiles::v, and verbose_.

 {
    using CLHEP::HepLorentzVector;
 
    const double fBfield=3.8;
   
    edm::Handle<edm::SimVertexContainer> simVtcs;
    iEvent.getByToken( edmSimVertexContainerToken_, simVtcs );
    if(verbose_){
      std::cout << "SimVertex " << simVtcs->size() << std::endl;
      for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
      v!=simVtcs->end(); ++v){
        std::cout << "simvtx "
          << std::setw(10) << std::setprecision(4)
          << v->position().x() << " "
          << v->position().y() << " "
          << v->position().z() << " "
          << v->parentIndex() << " "
          << v->noParent() << " "
          << std::endl;
      }
    }
    
    // get the simulated tracks, extract perigee parameters
    edm::Handle<edm::SimTrackContainer> simTrks;
    iEvent.getByToken( edmSimTrackContainerToken_, simTrks );
    
    if(verbose_){std::cout << "simtrks " << simTrks->size() << std::endl;}
    std::vector<ParameterVector > tsim;
    for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
        t!=simTrks->end(); ++t){
      if (t->noVertex()){
        std::cout << "simtrk  has no vertex" << std::endl;
        return;
      }else{
        // get the vertex position
        HepLorentzVector v((*simVtcs)[t->vertIndex()].position().x(),
                           (*simVtcs)[t->vertIndex()].position().y(),
                           (*simVtcs)[t->vertIndex()].position().z(),
                           (*simVtcs)[t->vertIndex()].position().e());
        int pdgCode=t->type();
 
        if( pdgCode==-99 ){
      // such entries cause crashes, no idea what they are
      std::cout << "funny particle skipped  , code="  << pdgCode << std::endl;
        }else{
      double Q=0;
      if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==321)){Q=-1;}
      else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)){Q=1;}
      else {
        std::cout << pdgCode << " " <<std::endl;
      }
      HepLorentzVector p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
      if(verbose_){
        std::cout << "simtrk "
                << " gen=" << std::setw( 4) << t->genpartIndex()
                << " vtx=" << std::setw( 4) << t->vertIndex() 
                << " pdg=" << std::setw( 5) << t->type() 
                << " Q="   << std::setw( 3) << Q 
                << " pt="  << std::setw(6) << p.perp()
                << " vx="  << std::setw(6) << v.x() 
                << " vy="  << std::setw(6) << v.y() 
                << " vz="  << std::setw(6) << v.z()
                << std::endl;
      }
      if ( (Q != 0) && (p.perp()>0.1) ){
        double x0=v.x()*simUnit_;
        double y0=v.y()*simUnit_;
        double z0=v.z()*simUnit_;
        double kappa=-Q*0.002998*fBfield/p.perp();
        double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
        double q=sqrt(1.-2.*kappa*D0);
        double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
        double s1;
        if (fabs(kappa*s0)>0.001){
          s1=asin(kappa*s0)/kappa;
        }else{
          double ks02=(kappa*s0)*(kappa*s0);
          s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
        }
        ParameterVector par;
        par[reco::TrackBase::i_qoverp] = Q/sqrt(p.perp2()+p.pz()*p.pz());
        par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
        par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
        par[reco::TrackBase::i_dxy] = 2.*D0/(1.+q);
        par[reco::TrackBase::i_dsz] = z0*sin(p.theta())-s1*cos(p.theta());
        tsim.push_back(par);
      }
        }
      }// has vertex
    }//for loop
    
    // simtrack parameters are in now tsim
    // loop over tracks and try to match them to simulated tracks
 
 
    edm::Handle<reco::TrackCollection> recTracks;
    iEvent.getByToken( recoTrackCollectionToken_, recTracks );
 
    for(reco::TrackCollection::const_iterator t=recTracks->begin();
        t!=recTracks->end(); ++t){
      reco::TrackBase::ParameterVector  p = t->parameters();
      reco::TrackBase::CovarianceMatrix c = t->covariance();
      if(verbose_){
        std::cout << "reco pars= " << p << std::endl;
      }
      for(std::vector<ParameterVector>::const_iterator s=tsim.begin();
      s!=tsim.end(); ++s){
        if (match(*s,p)){
      h1_pull0_->Fill((p(0)-(*s)(0))/sqrt(c(0,0)));
      h1_pull1_->Fill((p(1)-(*s)(1))/sqrt(c(1,1)));
      h1_pull2_->Fill((p(2)-(*s)(2))/sqrt(c(2,2)));
      h1_pull3_->Fill((p(3)-(*s)(3))/sqrt(c(3,3)));
      h1_pull4_->Fill((p(4)-(*s)(4))/sqrt(c(4,4)));
 
      h1_res0_->Fill(p(0)-(*s)(0));
      h1_res1_->Fill(p(1)-(*s)(1));
      h1_res2_->Fill(p(2)-(*s)(2));
      h1_res3_->Fill(p(3)-(*s)(3));
      h1_res4_->Fill(p(4)-(*s)(4));
 
      h1_Beff_->Fill(p(0)/(*s)(0)*fBfield);
      h2_dvsphi_->Fill(p(2),p(3));
      h1_par0_->Fill(p(0));
      h1_par1_->Fill(p(1));
      h1_par2_->Fill(p(2));
      h1_par3_->Fill(p(3));
      h1_par4_->Fill(p(4));
       }
      }
    }
 
 }

void TrackParameterAnalyzer::beginJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 68 of file TrackParameterAnalyzer.cc.

References gather_cfg::cout, h1_Beff_, h1_par0_, h1_par1_, h1_par2_, h1_par3_, h1_par4_, h1_pull0_, h1_pull1_, h1_pull2_, h1_pull3_, h1_pull4_, h1_res0_, h1_res1_, h1_res2_, h1_res3_, h1_res4_, h2_dvsphi_, M_PI, rootFile_, and simUnit_.

                                      {
   std::cout << " TrackParameterAnalyzer::beginJob  conversion from sim units to rec units is " << simUnit_ << std::endl;
 
   rootFile_->cd();
   h1_pull0_ = new TH1F("pull0","pull q/p",100,-25.,25.);
   h1_pull1_ = new TH1F("pull1","pull lambda",100,-25.,25.);
   h1_pull2_ = new TH1F("pull2","pull phi  ",100,-25.,25.);
   h1_pull3_ = new TH1F("pull3","pull dca  ",100,-25.,25.);
   h1_pull4_ = new TH1F("pull4","pull zdca ",100,-25.,25.);
 
   h1_res0_ = new TH1F("res0","res q/p",100,-0.1,0.1);
   h1_res1_ = new TH1F("res1","res lambda",100,-0.1,0.1);
   h1_res2_ = new TH1F("res2","res phi  ",100,-0.1,0.1);
   h1_res3_ = new TH1F("res3","res dca  ",100,-0.1,0.1);
   h1_res4_ = new TH1F("res4","res zdca ",100,-0.1,0.1);
 
   h1_Beff_  = new TH1F("Beff", "Beff",2000,-10.,10.);
   h2_dvsphi_ = new TH2F("dvsphi","dvsphi",360,-M_PI,M_PI,100,-0.1,0.1);
   h1_par0_ = new TH1F("par0","q/p",100,-0.1,0.1);
   h1_par1_ = new TH1F("par1","lambda",100, -M_PI/2.,M_PI/2.);
   h1_par2_ = new TH1F("par2","phi  ",100,-M_PI,M_PI);
   h1_par3_ = new TH1F("par3","dca  ",100,-0.1,0.1);
   h1_par4_ = new TH1F("par4","zdca ",1000,-10.,10.);
 }

void TrackParameterAnalyzer::endJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 94 of file TrackParameterAnalyzer.cc.

References h1_Beff_, h1_par0_, h1_par1_, h1_par2_, h1_par3_, h1_par4_, h1_pull0_, h1_pull1_, h1_pull2_, h1_pull3_, h1_pull4_, h1_res0_, h1_res1_, h1_res2_, h1_res3_, h1_res4_, h2_dvsphi_, and rootFile_.

                                     {
   rootFile_->cd();
   h1_pull0_->Write();
   h1_pull1_->Write();
   h1_pull2_->Write();
   h1_pull3_->Write();
   h1_pull4_->Write();
 
   h1_res0_->Write();
   h1_res1_->Write();
   h1_res2_->Write();
   h1_res3_->Write();
   h1_res4_->Write();
 
   h1_Beff_->Write();
   h2_dvsphi_->Write();
   h1_par0_->Write();
   h1_par1_->Write();
   h1_par2_->Write();
   h1_par3_->Write();
   h1_par4_->Write();
 }

bool TrackParameterAnalyzer::match	(	const ParameterVector &	a,
		const ParameterVector &	b
	)

private

Definition at line 118 of file TrackParameterAnalyzer.cc.

References a, b, and M_PI.

Referenced by analyze().

                                              {
   double dtheta=a(1)-b(1);
   double dphi  =a(2)-b(2);
   if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
   return ( (fabs(dtheta)<0.02) && (fabs(dphi)<0.04) );
 }