CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/RecoVertex/src/TrackParameterAnalyzer.cc

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