CMSSW Reference Manual

00001 #include "Validation/RecoVertex/interface/TrackParameterAnalyzer.h"
00002 #include <string>
00003 #include <vector>
00004 #include "FWCore/Utilities/interface/GetReleaseVersion.h"
00005 
00006 //
00007 //
00008 // constants, enums and typedefs
00009 //
00010 const double fBfield=4.06;
00011 
00012 //
00013 // static data member definitions
00014 //
00015 
00016 //
00017 // constructors and destructor
00018 //
00019 TrackParameterAnalyzer::TrackParameterAnalyzer(const edm::ParameterSet& iConfig) :
00020   simG4_( iConfig.getParameter<edm::InputTag>( "simG4" ) ) {
00021    //now do whatever initialization is needed
00022 
00023   recoTrackProducer_   = iConfig.getUntrackedParameter<std::string>("recoTrackProducer");
00024   // open output file to store histograms}
00025   outputFile_   = iConfig.getUntrackedParameter<std::string>("outputFile");
00026   TString tversion(edm::getReleaseVersion());
00027   tversion = tversion.Remove(0,1);
00028   tversion = tversion.Remove(tversion.Length()-1,tversion.Length());
00029   outputFile_  = std::string(tversion)+"_"+outputFile_;
00030 
00031   rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE"); 
00032   verbose_= iConfig.getUntrackedParameter<bool>("verbose", false);
00033   simUnit_=1.0;  //  starting from  CMSSW_1_2_x, I think
00034   if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
00035     simUnit_=0.1;  // for use in  CMSSW_1_1_1 tutorial
00036   }
00037 }
00038 
00039 
00040 TrackParameterAnalyzer::~TrackParameterAnalyzer()
00041 {
00042  
00043    // do anything here that needs to be done at destruction time
00044    // (e.g. close files, deallocate resources etc.)
00045   delete rootFile_;
00046 }
00047 
00048 
00049 
00050 //
00051 // member functions
00052 //
00053 void TrackParameterAnalyzer::beginJob(edm::EventSetup const&){
00054   std::cout << " TrackParameterAnalyzer::beginJob  conversion from sim units to rec units is " << simUnit_ << std::endl;
00055 
00056   rootFile_->cd();
00057   h1_pull0_ = new TH1F("pull0","pull q/p",100,-25.,25.);
00058   h1_pull1_ = new TH1F("pull1","pull lambda",100,-25.,25.);
00059   h1_pull2_ = new TH1F("pull2","pull phi  ",100,-25.,25.);
00060   h1_pull3_ = new TH1F("pull3","pull dca  ",100,-25.,25.);
00061   h1_pull4_ = new TH1F("pull4","pull zdca ",100,-25.,25.);
00062 
00063   h1_res0_ = new TH1F("res0","res q/p",100,-0.1,0.1);
00064   h1_res1_ = new TH1F("res1","res lambda",100,-0.1,0.1);
00065   h1_res2_ = new TH1F("res2","res phi  ",100,-0.1,0.1);
00066   h1_res3_ = new TH1F("res3","res dca  ",100,-0.1,0.1);
00067   h1_res4_ = new TH1F("res4","res zdca ",100,-0.1,0.1);
00068 
00069   h1_Beff_  = new TH1F("Beff", "Beff",2000,-10.,10.);
00070   h2_dvsphi_ = new TH2F("dvsphi","dvsphi",360,-M_PI,M_PI,100,-0.1,0.1);
00071   h1_par0_ = new TH1F("par0","q/p",100,-0.1,0.1);
00072   h1_par1_ = new TH1F("par1","lambda",100, -M_PI/2.,M_PI/2.);
00073   h1_par2_ = new TH1F("par2","phi  ",100,-M_PI,M_PI);
00074   h1_par3_ = new TH1F("par3","dca  ",100,-0.1,0.1);
00075   h1_par4_ = new TH1F("par4","zdca ",1000,-10.,10.);
00076 }
00077 
00078 
00079 void TrackParameterAnalyzer::endJob() {
00080   rootFile_->cd();
00081   h1_pull0_->Write();
00082   h1_pull1_->Write();
00083   h1_pull2_->Write();
00084   h1_pull3_->Write();
00085   h1_pull4_->Write();
00086 
00087   h1_res0_->Write();
00088   h1_res1_->Write();
00089   h1_res2_->Write();
00090   h1_res3_->Write();
00091   h1_res4_->Write();
00092 
00093   h1_Beff_->Write();
00094   h2_dvsphi_->Write();
00095   h1_par0_->Write();
00096   h1_par1_->Write();
00097   h1_par2_->Write();
00098   h1_par3_->Write();
00099   h1_par4_->Write();
00100 }
00101 
00102 // helper function
00103 bool TrackParameterAnalyzer::match(const ParameterVector  &a, 
00104                                    const ParameterVector  &b){
00105   double dtheta=a(1)-b(1);
00106   double dphi  =a(2)-b(2);
00107   if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
00108   return ( (fabs(dtheta)<0.02) && (fabs(dphi)<0.04) );
00109 }
00110 
00111 // ------------ method called to produce the data  ------------
00112 void
00113 TrackParameterAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00114 {
00115    using namespace edm;
00116   
00117    Handle<edm::SimVertexContainer> simVtcs;
00118    iEvent.getByLabel( simG4_, simVtcs);
00119    if(verbose_){
00120      std::cout << "SimVertex " << simVtcs->size() << std::endl;
00121      for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
00122          v!=simVtcs->end(); ++v){
00123        std::cout << "simvtx "
00124                  << std::setw(10) << std::setprecision(4)
00125                  << v->position().x() << " "
00126                  << v->position().y() << " "
00127                  << v->position().z() << " "
00128                  << v->parentIndex() << " "
00129                  << v->noParent() << " "
00130                  << std::endl;
00131      }
00132    }
00133    
00134    // get the simulated tracks, extract perigee parameters
00135    Handle<SimTrackContainer> simTrks;
00136    iEvent.getByLabel( simG4_, simTrks);
00137    
00138    if(verbose_){std::cout << "simtrks " << simTrks->size() << std::endl;}
00139    std::vector<ParameterVector > tsim;
00140    for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
00141        t!=simTrks->end(); ++t){
00142      if (t->noVertex()){
00143        std::cout << "simtrk  has no vertex" << std::endl;
00144        return;
00145      }else{
00146        // get the vertex position
00147        HepLorentzVector v((*simVtcs)[t->vertIndex()].position().x(),
00148                           (*simVtcs)[t->vertIndex()].position().y(),
00149                           (*simVtcs)[t->vertIndex()].position().z(),
00150                           (*simVtcs)[t->vertIndex()].position().e());
00151        int pdgCode=t->type();
00152 
00153        if( pdgCode==-99 ){
00154          // such entries cause crashes, no idea what they are
00155          std::cout << "funny particle skipped  , code="  << pdgCode << std::endl;
00156        }else{
00157          double Q=0; //double Q=HepPDT::theTable().getParticleData(pdgCode)->charge();
00158          if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==321)){Q=-1;}
00159          else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)){Q=1;}
00160          else {
00161            std::cout << pdgCode << " " <<std::endl;
00162          }
00163          HepLorentzVector p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
00164          if(verbose_){
00165            std::cout << "simtrk "
00166                        << " gen=" << std::setw( 4) << t->genpartIndex()
00167                        << " vtx=" << std::setw( 4) << t->vertIndex() 
00168                        << " pdg=" << std::setw( 5) << t->type() 
00169                        << " Q="   << std::setw( 3) << Q 
00170                        << " pt="  << std::setw(6) << p.perp()
00171                        << " vx="  << std::setw(6) << v.x() 
00172                        << " vy="  << std::setw(6) << v.y() 
00173                        << " vz="  << std::setw(6) << v.z()
00174                        << std::endl;
00175          }
00176          if ( (Q != 0) && (p.perp()>0.1) ){
00177            double x0=v.x()*simUnit_;
00178            double y0=v.y()*simUnit_;
00179            double z0=v.z()*simUnit_;
00180            double kappa=-Q*0.002998*fBfield/p.perp();
00181            double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
00182            double q=sqrt(1.-2.*kappa*D0);
00183            double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
00184            double s1;
00185            if (fabs(kappa*s0)>0.001){
00186              s1=asin(kappa*s0)/kappa;
00187            }else{
00188              double ks02=(kappa*s0)*(kappa*s0);
00189              s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
00190            }
00191            ParameterVector par;
00192            /*
00193              par[reco::TrackBase::i_transverseCurvature] = kappa;
00194              par[reco::TrackBase::i_theta] = p.theta();
00195              par[reco::TrackBase::i_phi0] = p.phi()-asin(kappa*s0);
00196              par[reco::TrackBase::i_d0] = 2.*D0/(1.+q);
00197              par[reco::TrackBase::i_dz] = z0-s1/tan(p.theta()); 
00198              */
00199            // new improved parametrization
00200            par[reco::TrackBase::i_qoverp] = Q/sqrt(p.perp2()+p.pz()*p.pz());
00201            par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
00202            par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
00203            par[reco::TrackBase::i_dxy] = 2.*D0/(1.+q);
00204            par[reco::TrackBase::i_dsz] = z0-s1/tan(p.theta()); 
00205            tsim.push_back(par);
00206          }
00207        }
00208      }// has vertex
00209    }//for loop
00210    
00211    // simtrack parameters are in now tsim
00212    // loop over tracks and try to match them to simulated tracks
00213 
00214 
00215    Handle<reco::TrackCollection> recTracks;
00216    iEvent.getByLabel(recoTrackProducer_, recTracks);
00217 
00218    for(reco::TrackCollection::const_iterator t=recTracks->begin();
00219        t!=recTracks->end(); ++t){
00220      reco::TrackBase::ParameterVector  p = t->parameters();
00221      reco::TrackBase::CovarianceMatrix c = t->covariance();
00222      if(verbose_){
00223        std::cout << "reco pars= " << p << std::endl;
00224        //std::cout << "z0=" << p(4) << std::endl;
00225      }
00226      for(std::vector<ParameterVector>::const_iterator s=tsim.begin();
00227          s!=tsim.end(); ++s){
00228        if (match(*s,p)){
00229          //if(verbose_){ std::cout << "match found" << std::endl;}
00230          h1_pull0_->Fill((p(0)-(*s)(0))/sqrt(c(0,0)));
00231          h1_pull1_->Fill((p(1)-(*s)(1))/sqrt(c(1,1)));
00232          h1_pull2_->Fill((p(2)-(*s)(2))/sqrt(c(2,2)));
00233          h1_pull3_->Fill((p(3)-(*s)(3))/sqrt(c(3,3)));
00234          h1_pull4_->Fill((p(4)-(*s)(4))/sqrt(c(4,4)));
00235 
00236          h1_res0_->Fill(p(0)-(*s)(0));
00237          h1_res1_->Fill(p(1)-(*s)(1));
00238          h1_res2_->Fill(p(2)-(*s)(2));
00239          h1_res3_->Fill(p(3)-(*s)(3));
00240          h1_res4_->Fill(p(4)-(*s)(4));
00241 
00242          h1_Beff_->Fill(p(0)/(*s)(0)*fBfield);
00243          h2_dvsphi_->Fill(p(2),p(3));
00244          h1_par0_->Fill(p(0));
00245          h1_par1_->Fill(p(1));
00246          h1_par2_->Fill(p(2));
00247          h1_par3_->Fill(p(3));
00248          h1_par4_->Fill(p(4));
00249       }
00250      }
00251    }
00252 
00253 
00254 
00255 
00256 }
00257
TrackParameterAnalyzer.cc
