Reference Manual

Go to the documentation of this file.
00001 
00014 #include "RecoVertex/BeamSpotProducer/interface/BeamFitter.h"
00015 #include "FWCore/ServiceRegistry/interface/Service.h"
00016 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00017 #include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
00018 
00019 #include "FWCore/Utilities/interface/InputTag.h"
00020 #include "DataFormats/Common/interface/View.h"
00021 
00022 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00023 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00024 #include "DataFormats/TrackReco/interface/Track.h"
00025 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00026 #include "DataFormats/TrackReco/interface/HitPattern.h"
00027 #include "DataFormats/VertexReco/interface/Vertex.h"
00028 
00029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00030 
00031 // ----------------------------------------------------------------------
00032 // Useful function:
00033 // ----------------------------------------------------------------------
00034 const char * BeamFitter::formatBTime(const std::time_t & t)  {
00035   struct std::tm * ptm;
00036   ptm = gmtime(&t);
00037   static char ts[] = "yyyy.mn.dd hh:mm:ss zzz ";
00038   // This value should be taken directly from edm::Event::time(), which
00039   // returns GMT (to be confirmed)
00040   strftime(ts,sizeof(ts),"%Y.%m.%d %H:%M:%S GMT",ptm);
00041   return ts;
00042 }
00043 // Update the string representations of the time
00044 void BeamFitter::updateBTime() {
00045   const char* fbeginTime = formatBTime(freftime[0]);
00046   sprintf(fbeginTimeOfFit,"%s",fbeginTime);
00047   const char* fendTime = formatBTime(freftime[1]);
00048   sprintf(fendTimeOfFit,"%s",fendTime);
00049 }
00050 
00051 
00052 BeamFitter::BeamFitter(const edm::ParameterSet& iConfig): fPVTree_(0)
00053 {
00054 
00055   debug_             = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("Debug");
00056   tracksLabel_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<edm::InputTag>("TrackCollection");
00057   vertexLabel_       = iConfig.getUntrackedParameter<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
00058   writeTxt_          = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteAscii");
00059   outputTxt_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("AsciiFileName");
00060   appendRunTxt_      = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("AppendRunToFileName");
00061   writeDIPTxt_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPAscii");
00062   // Specify whether we want to write the DIP file even if the fit is failed.
00063   writeDIPBadFit_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPOnBadFit", true);
00064   outputDIPTxt_      = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("DIPFileName");
00065   saveNtuple_        = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveNtuple");
00066   saveBeamFit_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveFitResults");
00067   savePVVertices_    = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SavePVVertices");
00068   isMuon_            = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("IsMuonCollection");
00069 
00070   trk_MinpT_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MinimumPt");
00071   trk_MaxEta_        = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumEta");
00072   trk_MaxIP_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumImpactParameter");
00073   trk_MaxZ_          = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
00074   trk_MinNTotLayers_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumTotalLayers");
00075   trk_MinNPixLayers_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumPixelLayers");
00076   trk_MaxNormChi2_   = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumNormChi2");
00077   trk_Algorithm_     = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::vector<std::string> >("TrackAlgorithm");
00078   trk_Quality_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::vector<std::string> >("TrackQuality");
00079   min_Ntrks_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
00080   convergence_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("FractionOfFittedTrks");
00081   inputBeamWidth_    = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("InputBeamWidth",-1.);
00082 
00083   for (unsigned int j=0;j<trk_Algorithm_.size();j++)
00084     algorithm_.push_back(reco::TrackBase::algoByName(trk_Algorithm_[j]));
00085   for (unsigned int j=0;j<trk_Quality_.size();j++)
00086     quality_.push_back(reco::TrackBase::qualityByName(trk_Quality_[j]));
00087 
00088   if (saveNtuple_ || saveBeamFit_ || savePVVertices_) {
00089     outputfilename_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("OutputFileName");
00090     file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
00091   }
00092   if (saveNtuple_) {
00093     ftree_ = new TTree("mytree","mytree");
00094     ftree_->AutoSave();
00095 
00096     ftree_->Branch("pt",&fpt,"fpt/D");
00097     ftree_->Branch("d0",&fd0,"fd0/D");
00098     ftree_->Branch("d0bs",&fd0bs,"fd0bs/D");
00099     ftree_->Branch("sigmad0",&fsigmad0,"fsigmad0/D");
00100     ftree_->Branch("phi0",&fphi0,"fphi0/D");
00101     ftree_->Branch("z0",&fz0,"fz0/D");
00102     ftree_->Branch("sigmaz0",&fsigmaz0,"fsigmaz0/D");
00103     ftree_->Branch("theta",&ftheta,"ftheta/D");
00104     ftree_->Branch("eta",&feta,"feta/D");
00105     ftree_->Branch("charge",&fcharge,"fcharge/I");
00106     ftree_->Branch("normchi2",&fnormchi2,"fnormchi2/D");
00107     ftree_->Branch("nTotLayerMeas",&fnTotLayerMeas,"fnTotLayerMeas/i");
00108     ftree_->Branch("nStripLayerMeas",&fnStripLayerMeas,"fnStripLayerMeas/i");
00109     ftree_->Branch("nPixelLayerMeas",&fnPixelLayerMeas,"fnPixelLayerMeas/i");
00110     ftree_->Branch("nTIBLayerMeas",&fnTIBLayerMeas,"fnTIBLayerMeas/i");
00111     ftree_->Branch("nTOBLayerMeas",&fnTOBLayerMeas,"fnTOBLayerMeas/i");
00112     ftree_->Branch("nTIDLayerMeas",&fnTIDLayerMeas,"fnTIDLayerMeas/i");
00113     ftree_->Branch("nTECLayerMeas",&fnTECLayerMeas,"fnTECLayerMeas/i");
00114     ftree_->Branch("nPXBLayerMeas",&fnPXBLayerMeas,"fnPXBLayerMeas/i");
00115     ftree_->Branch("nPXFLayerMeas",&fnPXFLayerMeas,"fnPXFLayerMeas/i");
00116     ftree_->Branch("cov",&fcov,"fcov[7][7]/D");
00117     ftree_->Branch("vx",&fvx,"fvx/D");
00118     ftree_->Branch("vy",&fvy,"fvy/D");
00119     ftree_->Branch("quality",&fquality,"fquality/O");
00120     ftree_->Branch("algo",&falgo,"falgo/O");
00121     ftree_->Branch("run",&frun,"frun/i");
00122     ftree_->Branch("lumi",&flumi,"flumi/i");
00123     ftree_->Branch("pvValid",&fpvValid,"fpvValid/O");
00124     ftree_->Branch("pvx", &fpvx, "fpvx/D");
00125     ftree_->Branch("pvy", &fpvy, "fpvy/D");
00126     ftree_->Branch("pvz", &fpvz, "fpvz/D");
00127   }
00128   if (saveBeamFit_){
00129     ftreeFit_ = new TTree("fitResults","fitResults");
00130     ftreeFit_->AutoSave();
00131     ftreeFit_->Branch("run",&frunFit,"frunFit/i");
00132     ftreeFit_->Branch("beginLumi",&fbeginLumiOfFit,"fbeginLumiOfFit/i");
00133     ftreeFit_->Branch("endLumi",&fendLumiOfFit,"fendLumiOfFit/i");
00134     ftreeFit_->Branch("beginTime",fbeginTimeOfFit,"fbeginTimeOfFit/C");
00135     ftreeFit_->Branch("endTime",fendTimeOfFit,"fendTimeOfFit/C");
00136     ftreeFit_->Branch("x",&fx,"fx/D");
00137     ftreeFit_->Branch("y",&fy,"fy/D");
00138     ftreeFit_->Branch("z",&fz,"fz/D");
00139     ftreeFit_->Branch("sigmaZ",&fsigmaZ,"fsigmaZ/D");
00140     ftreeFit_->Branch("dxdz",&fdxdz,"fdxdz/D");
00141     ftreeFit_->Branch("dydz",&fdydz,"fdydz/D");
00142     ftreeFit_->Branch("xErr",&fxErr,"fxErr/D");
00143     ftreeFit_->Branch("yErr",&fyErr,"fyErr/D");
00144     ftreeFit_->Branch("zErr",&fzErr,"fzErr/D");
00145     ftreeFit_->Branch("sigmaZErr",&fsigmaZErr,"fsigmaZErr/D");
00146     ftreeFit_->Branch("dxdzErr",&fdxdzErr,"fdxdzErr/D");
00147     ftreeFit_->Branch("dydzErr",&fdydzErr,"fdydzErr/D");
00148     ftreeFit_->Branch("widthX",&fwidthX,"fwidthX/D");
00149     ftreeFit_->Branch("widthY",&fwidthY,"fwidthY/D");
00150     ftreeFit_->Branch("widthXErr",&fwidthXErr,"fwidthXErr/D");
00151     ftreeFit_->Branch("widthYErr",&fwidthYErr,"fwidthYErr/D");
00152   }
00153 
00154   fBSvector.clear();
00155   ftotal_tracks = 0;
00156   fnTotLayerMeas = fnPixelLayerMeas = fnStripLayerMeas = fnTIBLayerMeas = 0;
00157   fnTIDLayerMeas = fnTOBLayerMeas = fnTECLayerMeas = fnPXBLayerMeas = fnPXFLayerMeas = 0;
00158   frun = flumi = -1;
00159   frunFit = fbeginLumiOfFit = fendLumiOfFit = -1;
00160   fquality = falgo = true;
00161   fpvValid = true;
00162   fpvx = fpvy = fpvz = 0;
00163   fitted_ = false;
00164   resetRefTime();
00165 
00166   //debug histograms
00167   h1ntrks = new TH1F("h1ntrks","number of tracks per event",50,0,50);
00168   h1vz_event = new TH1F("h1vz_event","track Vz", 50, -30, 30 );
00169   h1cutFlow = new TH1F("h1cutFlow","Cut flow table of track selection", 9, 0, 9);
00170   h1cutFlow->GetXaxis()->SetBinLabel(1,"No cut");
00171   h1cutFlow->GetXaxis()->SetBinLabel(2,"Traker hits");
00172   h1cutFlow->GetXaxis()->SetBinLabel(3,"Pixel hits");
00173   h1cutFlow->GetXaxis()->SetBinLabel(4,"norm. #chi^{2}");
00174   h1cutFlow->GetXaxis()->SetBinLabel(5,"algo");
00175   h1cutFlow->GetXaxis()->SetBinLabel(6,"quality");
00176   h1cutFlow->GetXaxis()->SetBinLabel(7,"d_{0}");
00177   h1cutFlow->GetXaxis()->SetBinLabel(8,"z_{0}");
00178   h1cutFlow->GetXaxis()->SetBinLabel(9,"p_{T}");
00179   resetCutFlow();
00180 
00181   // Primary vertex fitter
00182   MyPVFitter = new PVFitter(iConfig);
00183   MyPVFitter->resetAll();
00184   if (savePVVertices_){
00185     fPVTree_ = new TTree("PrimaryVertices","PrimaryVertices");
00186     MyPVFitter->setTree(fPVTree_);
00187   }
00188 
00189   // check filename
00190   ffilename_changed = false;
00191   if (writeDIPTxt_)
00192     fasciiDIP.open(outputDIPTxt_.c_str());
00193 }
00194 
00195 BeamFitter::~BeamFitter() {
00196 
00197   if (saveNtuple_) {
00198     file_->cd();
00199     if (fitted_ && h1z) h1z->Write();
00200     h1ntrks->Write();
00201     h1vz_event->Write();
00202     if (h1cutFlow) h1cutFlow->Write();
00203     ftree_->Write();
00204   }
00205   if (saveBeamFit_){
00206     file_->cd();
00207     ftreeFit_->Write();
00208   }
00209   if (savePVVertices_){
00210     file_->cd();
00211     fPVTree_->Write();
00212   }
00213 
00214 
00215   if (saveNtuple_ || saveBeamFit_ || savePVVertices_){
00216     file_->Close();
00217     delete file_;
00218   }
00219   delete MyPVFitter;
00220 }
00221 
00222 
00223 void BeamFitter::readEvent(const edm::Event& iEvent)
00224 {
00225 
00226 
00227   frun = iEvent.id().run();
00228   const edm::TimeValue_t ftimestamp = iEvent.time().value();
00229   const std::time_t ftmptime = ftimestamp >> 32;
00230 
00231   if (fbeginLumiOfFit == -1) freftime[0] = freftime[1] = ftmptime;
00232   if (freftime[0] == 0 || ftmptime < freftime[0]) freftime[0] = ftmptime;
00233   if (freftime[1] == 0 || ftmptime > freftime[1]) freftime[1] = ftmptime;
00234   // Update the human-readable string versions of the time
00235   updateBTime();
00236 
00237   flumi = iEvent.luminosityBlock();
00238   frunFit = frun;
00239   if (fbeginLumiOfFit == -1 || fbeginLumiOfFit > flumi) fbeginLumiOfFit = flumi;
00240   if (fendLumiOfFit == -1 || fendLumiOfFit < flumi) fendLumiOfFit = flumi;
00241 
00242   edm::Handle<reco::TrackCollection> TrackCollection;
00243   iEvent.getByLabel(tracksLabel_, TrackCollection);
00244 
00245   //------ Primary Vertices
00246   edm::Handle< edm::View<reco::Vertex> > PVCollection;
00247   bool hasPVs = false;
00248   edm::View<reco::Vertex> pv;
00249   if ( iEvent.getByLabel(vertexLabel_, PVCollection ) ) {
00250       pv = *PVCollection;
00251       hasPVs = true;
00252   }
00253   //------
00254 
00255   //------ Beam Spot in current event
00256   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00257   const reco::BeamSpot *refBS =  0;
00258   if ( iEvent.getByLabel("offlineBeamSpot",recoBeamSpotHandle) )
00259       refBS = recoBeamSpotHandle.product();
00260   //-------
00261 
00262   const reco::TrackCollection *tracks = TrackCollection.product();
00263 
00264   double eventZ = 0;
00265   double averageZ = 0;
00266 
00267   for ( reco::TrackCollection::const_iterator track = tracks->begin();
00268         track != tracks->end();
00269         ++track ) {
00270 
00271     if ( ! isMuon_) {
00272 
00273       const reco::HitPattern& trkHP = track->hitPattern();
00274 
00275       fnPixelLayerMeas = trkHP.pixelLayersWithMeasurement();
00276       fnStripLayerMeas = trkHP.stripLayersWithMeasurement();
00277       fnTotLayerMeas = trkHP.trackerLayersWithMeasurement();
00278       fnPXBLayerMeas = trkHP.pixelBarrelLayersWithMeasurement();
00279       fnPXFLayerMeas = trkHP.pixelEndcapLayersWithMeasurement();
00280       fnTIBLayerMeas = trkHP.stripTIBLayersWithMeasurement();
00281       fnTIDLayerMeas = trkHP.stripTIDLayersWithMeasurement();
00282       fnTOBLayerMeas = trkHP.stripTOBLayersWithMeasurement();
00283       fnTECLayerMeas = trkHP.stripTECLayersWithMeasurement();
00284     } else {
00285 
00286       fnTotLayerMeas = track->numberOfValidHits();
00287 
00288     }
00289 
00290     fpt = track->pt();
00291     feta = track->eta();
00292     fphi0 = track->phi();
00293     fcharge = track->charge();
00294     fnormchi2 = track->normalizedChi2();
00295     fd0 = track->d0();
00296     if (refBS) fd0bs = -1*track->dxy(refBS->position());
00297     else fd0bs = 0.;
00298 
00299     fsigmad0 = track->d0Error();
00300     fz0 = track->dz();
00301     fsigmaz0 = track->dzError();
00302     ftheta = track->theta();
00303     fvx = track->vx();
00304     fvy = track->vy();
00305 
00306     for (int i=0; i<5; ++i) {
00307       for (int j=0; j<5; ++j) {
00308         fcov[i][j] = track->covariance(i,j);
00309       }
00310     }
00311 
00312     fquality = true;
00313     falgo = true;
00314 
00315     if (! isMuon_ ) {
00316       if (quality_.size()!=0) {
00317           fquality = false;
00318           for (unsigned int i = 0; i<quality_.size();++i) {
00319               if(debug_) edm::LogInfo("BeamFitter") << "quality_[" << i << "] = " << track->qualityName(quality_[i]) << std::endl;
00320               if (track->quality(quality_[i])) {
00321                   fquality = true;
00322                   break;
00323               }
00324           }
00325       }
00326 
00327 
00328       // Track algorithm
00329 
00330       if (algorithm_.size()!=0) {
00331         if (std::find(algorithm_.begin(),algorithm_.end(),track->algo())==algorithm_.end())
00332           falgo = false;
00333       }
00334 
00335     }
00336 
00337     // check if we have a valid PV
00338     fpvValid = false;
00339 
00340     if ( hasPVs ) {
00341 
00342         for ( size_t ipv=0; ipv != pv.size(); ++ipv ) {
00343 
00344             if (! pv[ipv].isFake() ) fpvValid = true;
00345 
00346             if ( ipv==0 && !pv[0].isFake() ) { fpvx = pv[0].x(); fpvy = pv[0].y(); fpvz = pv[0].z(); } // fix this later
00347 
00348 
00349         }
00350 
00351     }
00352 
00353 
00354     if (saveNtuple_) ftree_->Fill();
00355     ftotal_tracks++;
00356 
00357     countPass[0] = ftotal_tracks;
00358     // Track selection
00359     if (fnTotLayerMeas >= trk_MinNTotLayers_) { countPass[1] += 1;
00360       if (fnPixelLayerMeas >= trk_MinNPixLayers_) { countPass[2] += 1;
00361         if (fnormchi2 < trk_MaxNormChi2_) { countPass[3] += 1;
00362           if (falgo) {countPass[4] += 1;
00363             if (fquality) { countPass[5] += 1;
00364               if (std::abs( fd0 ) < trk_MaxIP_) { countPass[6] += 1;
00365                 if (std::abs( fz0 ) < trk_MaxZ_){ countPass[7] += 1;
00366                   if (fpt > trk_MinpT_) {
00367                     countPass[8] += 1;
00368                     if (std::abs( feta ) < trk_MaxEta_
00369                         //&& fpvValid
00370                         ) {
00371                       if (debug_) {
00372                   edm::LogInfo("BeamFitter") << "Selected track quality = " << track->qualityMask()
00373                                              << "; track algorithm = " << track->algoName() << "= TrackAlgorithm: " << track->algo() << std::endl;
00374                       }
00375                       BSTrkParameters BSTrk(fz0,fsigmaz0,fd0,fsigmad0,fphi0,fpt,0.,0.);
00376                       BSTrk.setVx(fvx);
00377                       BSTrk.setVy(fvy);
00378                       fBSvector.push_back(BSTrk);
00379                       averageZ += fz0;
00380                     }
00381                   }
00382                 }
00383               }
00384             }
00385           }
00386         }
00387       }
00388     }// track selection
00389 
00390   }// tracks
00391 
00392   averageZ = averageZ/(float)(fBSvector.size());
00393 
00394   for( std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
00395 
00396     eventZ += fabs( iparam->z0() - averageZ );
00397 
00398   }
00399 
00400   h1ntrks->Fill( fBSvector.size() );
00401   h1vz_event->Fill( eventZ/(float)(fBSvector.size() ) ) ;
00402   for (unsigned int i=0; i < sizeof(countPass)/sizeof(countPass[0]); i++)
00403     h1cutFlow->SetBinContent(i+1,countPass[i]);
00404 
00405   MyPVFitter->readEvent(iEvent);
00406 
00407 }
00408 
00409 bool BeamFitter::runPVandTrkFitter() {
00410 // run both PV and track fitters
00411     bool fit_ok = false;
00412     bool pv_fit_ok = false;
00413     reco::BeamSpot bspotPV;
00414     reco::BeamSpot bspotTrk;
00415 
00416     // First run PV fitter
00417     if ( MyPVFitter->IsFitPerBunchCrossing() ){
00418       edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[0] = " << freftime[0]
00419                                  << "; address =  " << &freftime[0]
00420                                  << " = " << fbeginTimeOfFit << std::endl;
00421       edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[1] = " << freftime[1]
00422                                  << "; address =  " << &freftime[1]
00423                                  << " = " << fendTimeOfFit << std::endl;
00424 
00425       if ( MyPVFitter->runBXFitter() ) {
00426         fbspotPVMap = MyPVFitter->getBeamSpotMap();
00427         pv_fit_ok = true;
00428       }
00429       if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
00430       if(writeDIPTxt_ && (pv_fit_ok || writeDIPBadFit_)) {
00431           dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP
00432       }
00433       return pv_fit_ok;
00434     }
00435 
00436     if ( MyPVFitter->runFitter() ) {
00437 
00438         bspotPV = MyPVFitter->getBeamSpot();
00439 
00440         // take beam width from PV fit and pass it to track fitter
00441         // assume circular transverse beam width
00442         inputBeamWidth_ = sqrt( pow(bspotPV.BeamWidthX(),2) + pow(bspotPV.BeamWidthY(),2) )/sqrt(2);
00443         pv_fit_ok = true;
00444 
00445     } else {
00446         // problems with PV fit
00447         bspotPV.setType(reco::BeamSpot::Unknown);
00448         bspotTrk.setType(reco::BeamSpot::Unknown); //propagate error to trk beam spot
00449     }
00450 
00451     if ( runFitterNoTxt() ) {
00452 
00453         bspotTrk = fbeamspot;
00454         fit_ok = true;
00455     } else {
00456       // beamfit failed, flagged as empty beam spot
00457       bspotTrk.setType(reco::BeamSpot::Fake);
00458       fit_ok = false;
00459     }
00460 
00461     // combined results into one single beam spot
00462 
00463     // to pass errors I have to create another beam spot object
00464     reco::BeamSpot::CovarianceMatrix matrix;
00465     for (int j = 0 ; j < 7 ; ++j) {
00466         for(int k = j ; k < 7 ; ++k) {
00467             matrix(j,k) = bspotTrk.covariance(j,k);
00468         }
00469     }
00470     // change beam width error to one from PV
00471     if (pv_fit_ok && fit_ok ) {
00472       matrix(6,6) = MyPVFitter->getWidthXerr() * MyPVFitter->getWidthXerr();
00473 
00474       // get Z and sigmaZ from PV fit
00475       matrix(2,2) = bspotPV.covariance(2,2);
00476       matrix(3,3) = bspotPV.covariance(3,3);
00477       reco::BeamSpot tmpbs(reco::BeamSpot::Point(bspotTrk.x0(), bspotTrk.y0(),
00478                                                  bspotPV.z0() ),
00479                            bspotPV.sigmaZ() ,
00480                            bspotTrk.dxdz(),
00481                            bspotTrk.dydz(),
00482                            bspotPV.BeamWidthX(),
00483                            matrix,
00484                            bspotPV.type() );
00485       tmpbs.setBeamWidthY( bspotPV.BeamWidthY() );
00486       // overwrite beam spot result
00487       fbeamspot = tmpbs;
00488     }
00489     if (pv_fit_ok && fit_ok) {
00490       fbeamspot.setType(bspotPV.type());
00491     }
00492     else if(!pv_fit_ok && fit_ok){
00493       fbeamspot.setType(reco::BeamSpot::Unknown);
00494     }
00495     else if(pv_fit_ok && !fit_ok){
00496       fbeamspot.setType(reco::BeamSpot::Unknown);
00497     }
00498     else if(!pv_fit_ok && !fit_ok){
00499       fbeamspot.setType(reco::BeamSpot::Fake);
00500     }
00501 
00502     if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
00503     if(writeDIPTxt_ && ((fit_ok && pv_fit_ok) || writeDIPBadFit_)) {
00504         dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP
00505         for(size_t i= 0; i < 7; i++)ForDIPPV_.push_back(0.0);
00506     }
00507 
00508     return fit_ok && pv_fit_ok;
00509 }
00510 
00511 bool BeamFitter::runFitterNoTxt() {
00512   edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[0] = " << freftime[0]
00513                              << "; address =  " << &freftime[0]
00514                              << " = " << fbeginTimeOfFit << std::endl;
00515   edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[1] = " << freftime[1]
00516                              << "; address =  " << &freftime[1]
00517                              << " = " << fendTimeOfFit << std::endl;
00518 
00519   if (fbeginLumiOfFit == -1 || fendLumiOfFit == -1) {
00520       edm::LogWarning("BeamFitter") << "No event read! No Fitting!" << std::endl;
00521     return false;
00522   }
00523 
00524   bool fit_ok = false;
00525   // default fit to extract beam spot info
00526   if(fBSvector.size() > 1 ){
00527 
00528       edm::LogInfo("BeamFitter") << "Calculating beam spot..." << std::endl
00529                                  << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks << std::endl;
00530 
00531     BSFitter *myalgo = new BSFitter( fBSvector );
00532     myalgo->SetMaximumZ( trk_MaxZ_ );
00533     myalgo->SetConvergence( convergence_ );
00534     myalgo->SetMinimumNTrks( min_Ntrks_ );
00535     if (inputBeamWidth_ > 0 ) myalgo->SetInputBeamWidth( inputBeamWidth_ );
00536 
00537     fbeamspot = myalgo->Fit();
00538 
00539 
00540     // retrieve histogram for Vz
00541     h1z = (TH1F*) myalgo->GetVzHisto();
00542 
00543     delete myalgo;
00544     if ( fbeamspot.type() != 0 ) {// save all results except for Fake (all 0.)
00545       fit_ok = true;
00546       if (saveBeamFit_){
00547         fx = fbeamspot.x0();
00548         fy = fbeamspot.y0();
00549         fz = fbeamspot.z0();
00550         fsigmaZ = fbeamspot.sigmaZ();
00551         fdxdz = fbeamspot.dxdz();
00552         fdydz = fbeamspot.dydz();
00553         fwidthX = fbeamspot.BeamWidthX();
00554         fwidthY = fbeamspot.BeamWidthY();
00555         fxErr = fbeamspot.x0Error();
00556         fyErr = fbeamspot.y0Error();
00557         fzErr = fbeamspot.z0Error();
00558         fsigmaZErr = fbeamspot.sigmaZ0Error();
00559         fdxdzErr = fbeamspot.dxdzError();
00560         fdydzErr = fbeamspot.dydzError();
00561         fwidthXErr = fbeamspot.BeamWidthXError();
00562         fwidthYErr = fbeamspot.BeamWidthYError();
00563         ftreeFit_->Fill();
00564       }
00565     }
00566   }
00567   else{ // tracks <= 1
00568     reco::BeamSpot tmpbs;
00569     fbeamspot = tmpbs;
00570     fbeamspot.setType(reco::BeamSpot::Fake);
00571     edm::LogInfo("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;
00572 
00573   }
00574   fitted_ = true;
00575   return fit_ok;
00576 
00577 }
00578 
00579 bool BeamFitter::runFitter() {
00580 
00581     bool fit_ok = runFitterNoTxt();
00582 
00583     if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
00584     if(writeDIPTxt_ && (fit_ok || writeDIPBadFit_)) {
00585       dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP
00586     }
00587     return fit_ok;
00588 }
00589 
00590 bool BeamFitter::runBeamWidthFitter() {
00591   bool widthfit_ok = false;
00592   // default fit to extract beam spot info
00593   if(fBSvector.size() > 1 ){
00594 
00595       edm::LogInfo("BeamFitter") << "Calculating beam spot positions('d0-phi0' method) and width using llh Fit"<< std::endl
00596                                  << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks << std::endl;
00597 
00598         BSFitter *myalgo = new BSFitter( fBSvector );
00599         myalgo->SetMaximumZ( trk_MaxZ_ );
00600         myalgo->SetConvergence( convergence_ );
00601         myalgo->SetMinimumNTrks(min_Ntrks_);
00602         if (inputBeamWidth_ > 0 ) myalgo->SetInputBeamWidth( inputBeamWidth_ );
00603 
00604 
00605    myalgo->SetFitVariable(std::string("d*z"));
00606    myalgo->SetFitType(std::string("likelihood"));
00607    fbeamWidthFit = myalgo->Fit();
00608 
00609    //Add to .txt file
00610    if(writeTxt_   ) dumpBWTxtFile(outputTxt_);
00611 
00612    delete myalgo;
00613 
00614    // not fake
00615    if ( fbeamspot.type() != 0 )
00616        widthfit_ok = true;
00617   }
00618   else{
00619     fbeamspot.setType(reco::BeamSpot::Fake);
00620     edm::LogWarning("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;
00621   }
00622   return widthfit_ok;
00623 }
00624 
00625 void BeamFitter::dumpBWTxtFile(std::string & fileName ){
00626     std::ofstream outFile;
00627     outFile.open(fileName.c_str(),std::ios::app);
00628     outFile<<"-------------------------------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
00629     outFile<<"Beam width(in cm) from Log-likelihood fit (Here we assume a symmetric beam(SigmaX=SigmaY)!)"<<std::endl;
00630     outFile<<"   "<<std::endl;
00631     outFile << "BeamWidth =  " <<fbeamWidthFit.BeamWidthX() <<" +/- "<<fbeamWidthFit.BeamWidthXError() << std::endl;
00632     outFile.close();
00633 }
00634 
00635 void BeamFitter::dumpTxtFile(std::string & fileName, bool append){
00636   std::ofstream outFile;
00637 
00638   std::string tmpname = outputTxt_;
00639   char index[15];
00640   if (appendRunTxt_ && !ffilename_changed ) {
00641       sprintf(index,"%s%i","_Run", frun );
00642       tmpname.insert(outputTxt_.length()-4,index);
00643       fileName = tmpname;
00644       ffilename_changed = true;
00645   }
00646 
00647   if (!append)
00648     outFile.open(fileName.c_str());
00649   else
00650     outFile.open(fileName.c_str(),std::ios::app);
00651 
00652   if ( MyPVFitter->IsFitPerBunchCrossing() ) {
00653 
00654     for (std::map<int,reco::BeamSpot>::const_iterator abspot = fbspotPVMap.begin(); abspot!= fbspotPVMap.end(); ++abspot) {
00655       reco::BeamSpot beamspottmp = abspot->second;
00656       int bx = abspot->first;
00657 
00658       outFile << "Runnumber " << frun << " bx " << bx << std::endl;
00659       outFile << "BeginTimeOfFit " << fbeginTimeOfFit << " " << freftime[0] << std::endl;
00660       outFile << "EndTimeOfFit " << fendTimeOfFit << " " << freftime[1] << std::endl;
00661       outFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
00662       outFile << "Type " << beamspottmp.type() << std::endl;
00663       outFile << "X0 " << beamspottmp.x0() << std::endl;
00664       outFile << "Y0 " << beamspottmp.y0() << std::endl;
00665       outFile << "Z0 " << beamspottmp.z0() << std::endl;
00666       outFile << "sigmaZ0 " << beamspottmp.sigmaZ() << std::endl;
00667       outFile << "dxdz " << beamspottmp.dxdz() << std::endl;
00668       outFile << "dydz " << beamspottmp.dydz() << std::endl;
00669       outFile << "BeamWidthX " << beamspottmp.BeamWidthX() << std::endl;
00670       outFile << "BeamWidthY " << beamspottmp.BeamWidthY() << std::endl;
00671       for (int i = 0; i<6; ++i) {
00672         outFile << "Cov("<<i<<",j) ";
00673         for (int j=0; j<7; ++j) {
00674           outFile << beamspottmp.covariance(i,j) << " ";
00675         }
00676         outFile << std::endl;
00677       }
00678       outFile << "Cov(6,j) 0 0 0 0 0 0 " << beamspottmp.covariance(6,6) << std::endl;
00679       //}
00680       outFile << "EmittanceX " << beamspottmp.emittanceX() << std::endl;
00681       outFile << "EmittanceY " << beamspottmp.emittanceY() << std::endl;
00682       outFile << "BetaStar " << beamspottmp.betaStar() << std::endl;
00683 
00684     }
00685   }//if bx results needed
00686   else {
00687     outFile << "Runnumber " << frun << std::endl;
00688     outFile << "BeginTimeOfFit " << fbeginTimeOfFit << " " << freftime[0] << std::endl;
00689     outFile << "EndTimeOfFit " << fendTimeOfFit << " " << freftime[1] << std::endl;
00690     outFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
00691     outFile << "Type " << fbeamspot.type() << std::endl;
00692     outFile << "X0 " << fbeamspot.x0() << std::endl;
00693     outFile << "Y0 " << fbeamspot.y0() << std::endl;
00694     outFile << "Z0 " << fbeamspot.z0() << std::endl;
00695     outFile << "sigmaZ0 " << fbeamspot.sigmaZ() << std::endl;
00696     outFile << "dxdz " << fbeamspot.dxdz() << std::endl;
00697     outFile << "dydz " << fbeamspot.dydz() << std::endl;
00698     //  if (inputBeamWidth_ > 0 ) {
00699     //    outFile << "BeamWidthX " << inputBeamWidth_ << std::endl;
00700     //    outFile << "BeamWidthY " << inputBeamWidth_ << std::endl;
00701     //  } else {
00702     outFile << "BeamWidthX " << fbeamspot.BeamWidthX() << std::endl;
00703     outFile << "BeamWidthY " << fbeamspot.BeamWidthY() << std::endl;
00704     //  }
00705 
00706     for (int i = 0; i<6; ++i) {
00707       outFile << "Cov("<<i<<",j) ";
00708       for (int j=0; j<7; ++j) {
00709         outFile << fbeamspot.covariance(i,j) << " ";
00710       }
00711       outFile << std::endl;
00712     }
00713 
00714     // beam width error
00715     //if (inputBeamWidth_ > 0 ) {
00716     //  outFile << "Cov(6,j) 0 0 0 0 0 0 " << "1e-4" << std::endl;
00717     //} else {
00718     outFile << "Cov(6,j) 0 0 0 0 0 0 " << fbeamspot.covariance(6,6) << std::endl;
00719     //}
00720     outFile << "EmittanceX " << fbeamspot.emittanceX() << std::endl;
00721     outFile << "EmittanceY " << fbeamspot.emittanceY() << std::endl;
00722     outFile << "BetaStar " << fbeamspot.betaStar() << std::endl;
00723 
00724     //write here Pv info for DIP only: This added only if append is false, which happen for DIP only :)
00725   if(!append){
00726     outFile << "events "<< (int)ForDIPPV_[0] << std::endl;
00727     outFile << "meanPV "<< ForDIPPV_[1] << std::endl;
00728     outFile << "meanErrPV "<< ForDIPPV_[2] << std::endl;
00729     outFile << "rmsPV "<< ForDIPPV_[3] << std::endl;
00730     outFile << "rmsErrPV "<< ForDIPPV_[4] << std::endl;
00731     outFile << "maxPV "<< (int)ForDIPPV_[5] << std::endl;
00732     outFile << "nPV "<< (int)ForDIPPV_[6] << std::endl;
00733    }//writeDIPPVInfo_
00734   }//else end  here
00735 
00736   outFile.close();
00737 }
00738 
00739 void BeamFitter::write2DB(){
00740   BeamSpotObjects *pBSObjects = new BeamSpotObjects();
00741 
00742   pBSObjects->SetPosition(fbeamspot.position().X(),fbeamspot.position().Y(),fbeamspot.position().Z());
00743   //std::cout << " wrote: x= " << fbeamspot.position().X() << " y= "<< fbeamspot.position().Y() << " z= " << fbeamspot.position().Z() << std::endl;
00744   pBSObjects->SetSigmaZ(fbeamspot.sigmaZ());
00745   pBSObjects->Setdxdz(fbeamspot.dxdz());
00746   pBSObjects->Setdydz(fbeamspot.dydz());
00747   //if (inputBeamWidth_ > 0 ) {
00748   //  std::cout << " beam width value forced to be " << inputBeamWidth_ << std::endl;
00749   //  pBSObjects->SetBeamWidthX(inputBeamWidth_);
00750   //  pBSObjects->SetBeamWidthY(inputBeamWidth_);
00751   //} else {
00752     // need to fix this
00753     //std::cout << " using default value, 15e-4, for beam width!!!"<<std::endl;
00754   pBSObjects->SetBeamWidthX(fbeamspot.BeamWidthX() );
00755   pBSObjects->SetBeamWidthY(fbeamspot.BeamWidthY() );
00756   //}
00757 
00758   for (int i = 0; i<7; ++i) {
00759     for (int j=0; j<7; ++j) {
00760       pBSObjects->SetCovariance(i,j,fbeamspot.covariance(i,j));
00761     }
00762   }
00763   edm::Service<cond::service::PoolDBOutputService> poolDbService;
00764   if( poolDbService.isAvailable() ) {
00765     std::cout << "poolDBService available"<<std::endl;
00766     if ( poolDbService->isNewTagRequest( "BeamSpotObjectsRcd" ) ) {
00767       std::cout << "new tag requested" << std::endl;
00768       poolDbService->createNewIOV<BeamSpotObjects>( pBSObjects, poolDbService->beginOfTime(),poolDbService->endOfTime(),
00769                                                     "BeamSpotObjectsRcd"  );
00770     }
00771     else {
00772       std::cout << "no new tag requested" << std::endl;
00773       poolDbService->appendSinceTime<BeamSpotObjects>( pBSObjects, poolDbService->currentTime(),
00774                                                        "BeamSpotObjectsRcd" );
00775     }
00776   }
00777 }
00778 
00779 void BeamFitter::runAllFitter() {
00780   if(fBSvector.size()!=0){
00781     BSFitter *myalgo = new BSFitter( fBSvector );
00782     fbeamspot = myalgo->Fit_d0phi();
00783 
00784     // iterative
00785     if(debug_)
00786       std::cout << " d0-phi Iterative:" << std::endl;
00787     BSFitter *myitealgo = new BSFitter( fBSvector );
00788     myitealgo->Setd0Cut_d0phi(4.0);
00789     reco::BeamSpot beam_ite = myitealgo->Fit_ited0phi();
00790     if (debug_){
00791       std::cout << beam_ite << std::endl;
00792       std::cout << "\n Now run tests of the different fits\n";
00793     }
00794     // from here are just tests
00795     std::string fit_type = "chi2";
00796     myalgo->SetFitVariable(std::string("z"));
00797     myalgo->SetFitType(std::string("chi2"));
00798     reco::BeamSpot beam_fit_z_chi2 = myalgo->Fit();
00799     if (debug_){
00800       std::cout << " z Chi2 Fit ONLY:" << std::endl;
00801       std::cout << beam_fit_z_chi2 << std::endl;
00802     }
00803 
00804     fit_type = "combined";
00805     myalgo->SetFitVariable("z");
00806     myalgo->SetFitType("combined");
00807     reco::BeamSpot beam_fit_z_lh = myalgo->Fit();
00808     if (debug_){
00809       std::cout << " z Log-Likelihood Fit ONLY:" << std::endl;
00810       std::cout << beam_fit_z_lh << std::endl;
00811     }
00812 
00813     myalgo->SetFitVariable("d");
00814     myalgo->SetFitType("d0phi");
00815     reco::BeamSpot beam_fit_dphi = myalgo->Fit();
00816     if (debug_){
00817       std::cout << " d0-phi0 Fit ONLY:" << std::endl;
00818       std::cout << beam_fit_dphi << std::endl;
00819     }
00820     /*
00821     myalgo->SetFitVariable(std::string("d*z"));
00822     myalgo->SetFitType(std::string("likelihood"));
00823     reco::BeamSpot beam_fit_dz_lh = myalgo->Fit();
00824     if (debug_){
00825       std::cout << " Log-Likelihood Fit:" << std::endl;
00826       std::cout << beam_fit_dz_lh << std::endl;
00827     }
00828 
00829     myalgo->SetFitVariable(std::string("d*z"));
00830     myalgo->SetFitType(std::string("resolution"));
00831     reco::BeamSpot beam_fit_dresz_lh = myalgo->Fit();
00832     if(debug_){
00833       std::cout << " IP Resolution Fit" << std::endl;
00834       std::cout << beam_fit_dresz_lh << std::endl;
00835 
00836       std::cout << "c0 = " << myalgo->GetResPar0() << " +- " << myalgo->GetResPar0Err() << std::endl;
00837       std::cout << "c1 = " << myalgo->GetResPar1() << " +- " << myalgo->GetResPar1Err() << std::endl;
00838     }
00839     */
00840   }
00841   else
00842     if (debug_) std::cout << "No good track selected! No beam fit!" << std::endl;
00843 }
00844
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/RecoVertex/BeamSpotProducer/src/BeamFitter.cc
