Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:05 

0001 /**_________________________________________________________________
0002    class:   BeamFitter.cc
0003    package: RecoVertex/BeamSpotProducer
0004 
0005 
0006 
0007    author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
0008            Geng-Yuan Jeng, UC Riverside (Geng-Yuan.Jeng@cern.ch)
0009 
0010 
0011 ________________________________________________________________**/
0012 
0013 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0014 #include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
0015 
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "FWCore/Framework/interface/ConsumesCollector.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 
0021 #include "RecoVertex/BeamSpotProducer/interface/BeamFitter.h"
0022 #include "RecoVertex/BeamSpotProducer/interface/BeamSpotWrite2Txt.h"
0023 
0024 #include "DataFormats/Common/interface/View.h"
0025 #include "DataFormats/TrackReco/interface/HitPattern.h"
0026 #include "DataFormats/VertexReco/interface/Vertex.h"
0027 
0028 // Update the string representations of the time
0029 void BeamFitter::updateBTime() {
0030   char ts[] = "yyyy.mn.dd hh:mm:ss zzz ";
0031   char *fbeginTime = ts;
0032   strftime(fbeginTime, sizeof(ts), "%Y.%m.%d %H:%M:%S GMT", gmtime(&freftime[0]));
0033   sprintf(fbeginTimeOfFit, "%s", fbeginTime);
0034   char *fendTime = ts;
0035   strftime(fendTime, sizeof(ts), "%Y.%m.%d %H:%M:%S GMT", gmtime(&freftime[1]));
0036   sprintf(fendTimeOfFit, "%s", fendTime);
0037 }
0038 
0039 BeamFitter::BeamFitter(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iColl) : fPVTree_(nullptr) {
0040   debug_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("Debug");
0041   tracksToken_ = iColl.consumes<reco::TrackCollection>(
0042       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<edm::InputTag>("TrackCollection"));
0043   vertexToken_ = iColl.consumes<edm::View<reco::Vertex> >(
0044       iConfig.getUntrackedParameter<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices")));
0045   beamSpotToken_ = iColl.consumes<reco::BeamSpot>(
0046       iConfig.getUntrackedParameter<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot")));
0047   writeTxt_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteAscii");
0048   outputTxt_ =
0049       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("AsciiFileName");
0050   appendRunTxt_ =
0051       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("AppendRunToFileName");
0052   writeDIPTxt_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPAscii");
0053   // Specify whether we want to write the DIP file even if the fit is failed.
0054   writeDIPBadFit_ =
0055       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPOnBadFit", true);
0056   outputDIPTxt_ =
0057       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("DIPFileName");
0058   saveNtuple_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveNtuple");
0059   saveBeamFit_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveFitResults");
0060   savePVVertices_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SavePVVertices");
0061   isMuon_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("IsMuonCollection");
0062 
0063   trk_MinpT_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MinimumPt");
0064   trk_MaxEta_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumEta");
0065   trk_MaxIP_ =
0066       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumImpactParameter");
0067   trk_MaxZ_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
0068   trk_MinNTotLayers_ =
0069       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumTotalLayers");
0070   trk_MinNPixLayers_ =
0071       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumPixelLayers");
0072   trk_MaxNormChi2_ =
0073       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumNormChi2");
0074   trk_Algorithm_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter")
0075                        .getUntrackedParameter<std::vector<std::string> >("TrackAlgorithm");
0076   trk_Quality_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter")
0077                      .getUntrackedParameter<std::vector<std::string> >("TrackQuality");
0078   min_Ntrks_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
0079   convergence_ =
0080       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("FractionOfFittedTrks");
0081   inputBeamWidth_ =
0082       iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("InputBeamWidth", -1.);
0083 
0084   for (unsigned int j = 0; j < trk_Algorithm_.size(); j++)
0085     algorithm_.push_back(reco::TrackBase::algoByName(trk_Algorithm_[j]));
0086   for (unsigned int j = 0; j < trk_Quality_.size(); j++)
0087     quality_.push_back(reco::TrackBase::qualityByName(trk_Quality_[j]));
0088 
0089   if (saveNtuple_ || saveBeamFit_ || savePVVertices_) {
0090     outputfilename_ =
0091         iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("OutputFileName");
0092     file_ = TFile::Open(outputfilename_.c_str(), "RECREATE");
0093   }
0094   if (saveNtuple_) {
0095     ftree_ = new TTree("mytree", "mytree");
0096     ftree_->AutoSave();
0097 
0098     ftree_->Branch("pt", &fpt, "fpt/D");
0099     ftree_->Branch("d0", &fd0, "fd0/D");
0100     ftree_->Branch("d0bs", &fd0bs, "fd0bs/D");
0101     ftree_->Branch("sigmad0", &fsigmad0, "fsigmad0/D");
0102     ftree_->Branch("phi0", &fphi0, "fphi0/D");
0103     ftree_->Branch("z0", &fz0, "fz0/D");
0104     ftree_->Branch("sigmaz0", &fsigmaz0, "fsigmaz0/D");
0105     ftree_->Branch("theta", &ftheta, "ftheta/D");
0106     ftree_->Branch("eta", &feta, "feta/D");
0107     ftree_->Branch("charge", &fcharge, "fcharge/I");
0108     ftree_->Branch("normchi2", &fnormchi2, "fnormchi2/D");
0109     ftree_->Branch("nTotLayerMeas", &fnTotLayerMeas, "fnTotLayerMeas/i");
0110     ftree_->Branch("nStripLayerMeas", &fnStripLayerMeas, "fnStripLayerMeas/i");
0111     ftree_->Branch("nPixelLayerMeas", &fnPixelLayerMeas, "fnPixelLayerMeas/i");
0112     ftree_->Branch("nTIBLayerMeas", &fnTIBLayerMeas, "fnTIBLayerMeas/i");
0113     ftree_->Branch("nTOBLayerMeas", &fnTOBLayerMeas, "fnTOBLayerMeas/i");
0114     ftree_->Branch("nTIDLayerMeas", &fnTIDLayerMeas, "fnTIDLayerMeas/i");
0115     ftree_->Branch("nTECLayerMeas", &fnTECLayerMeas, "fnTECLayerMeas/i");
0116     ftree_->Branch("nPXBLayerMeas", &fnPXBLayerMeas, "fnPXBLayerMeas/i");
0117     ftree_->Branch("nPXFLayerMeas", &fnPXFLayerMeas, "fnPXFLayerMeas/i");
0118     ftree_->Branch("cov", &fcov, "fcov[7][7]/D");
0119     ftree_->Branch("vx", &fvx, "fvx/D");
0120     ftree_->Branch("vy", &fvy, "fvy/D");
0121     ftree_->Branch("quality", &fquality, "fquality/O");
0122     ftree_->Branch("algo", &falgo, "falgo/O");
0123     ftree_->Branch("run", &frun, "frun/i");
0124     ftree_->Branch("lumi", &flumi, "flumi/i");
0125     ftree_->Branch("pvValid", &fpvValid, "fpvValid/O");
0126     ftree_->Branch("pvx", &fpvx, "fpvx/D");
0127     ftree_->Branch("pvy", &fpvy, "fpvy/D");
0128     ftree_->Branch("pvz", &fpvz, "fpvz/D");
0129   }
0130   if (saveBeamFit_) {
0131     ftreeFit_ = new TTree("fitResults", "fitResults");
0132     ftreeFit_->AutoSave();
0133     ftreeFit_->Branch("run", &frunFit, "frunFit/i");
0134     ftreeFit_->Branch("beginLumi", &fbeginLumiOfFit, "fbeginLumiOfFit/i");
0135     ftreeFit_->Branch("endLumi", &fendLumiOfFit, "fendLumiOfFit/i");
0136     ftreeFit_->Branch("beginTime", fbeginTimeOfFit, "fbeginTimeOfFit/C");
0137     ftreeFit_->Branch("endTime", fendTimeOfFit, "fendTimeOfFit/C");
0138     ftreeFit_->Branch("x", &fx, "fx/D");
0139     ftreeFit_->Branch("y", &fy, "fy/D");
0140     ftreeFit_->Branch("z", &fz, "fz/D");
0141     ftreeFit_->Branch("sigmaZ", &fsigmaZ, "fsigmaZ/D");
0142     ftreeFit_->Branch("dxdz", &fdxdz, "fdxdz/D");
0143     ftreeFit_->Branch("dydz", &fdydz, "fdydz/D");
0144     ftreeFit_->Branch("xErr", &fxErr, "fxErr/D");
0145     ftreeFit_->Branch("yErr", &fyErr, "fyErr/D");
0146     ftreeFit_->Branch("zErr", &fzErr, "fzErr/D");
0147     ftreeFit_->Branch("sigmaZErr", &fsigmaZErr, "fsigmaZErr/D");
0148     ftreeFit_->Branch("dxdzErr", &fdxdzErr, "fdxdzErr/D");
0149     ftreeFit_->Branch("dydzErr", &fdydzErr, "fdydzErr/D");
0150     ftreeFit_->Branch("widthX", &fwidthX, "fwidthX/D");
0151     ftreeFit_->Branch("widthY", &fwidthY, "fwidthY/D");
0152     ftreeFit_->Branch("widthXErr", &fwidthXErr, "fwidthXErr/D");
0153     ftreeFit_->Branch("widthYErr", &fwidthYErr, "fwidthYErr/D");
0154   }
0155 
0156   fBSvector.clear();
0157   ftotal_tracks = 0;
0158   fnTotLayerMeas = fnPixelLayerMeas = fnStripLayerMeas = fnTIBLayerMeas = 0;
0159   fnTIDLayerMeas = fnTOBLayerMeas = fnTECLayerMeas = fnPXBLayerMeas = fnPXFLayerMeas = 0;
0160   frun = flumi = -1;
0161   frunFit = fbeginLumiOfFit = fendLumiOfFit = -1;
0162   fquality = falgo = true;
0163   fpvValid = true;
0164   fpvx = fpvy = fpvz = 0;
0165   fitted_ = false;
0166   resetRefTime();
0167 
0168   //debug histograms
0169   h1ntrks = new TH1F("h1ntrks", "number of tracks per event", 50, 0, 50);
0170   h1vz_event = new TH1F("h1vz_event", "track Vz", 50, -30, 30);
0171   h1cutFlow = new TH1F("h1cutFlow", "Cut flow table of track selection", 9, 0, 9);
0172   h1cutFlow->GetXaxis()->SetBinLabel(1, "No cut");
0173   h1cutFlow->GetXaxis()->SetBinLabel(2, "Traker hits");
0174   h1cutFlow->GetXaxis()->SetBinLabel(3, "Pixel hits");
0175   h1cutFlow->GetXaxis()->SetBinLabel(4, "norm. #chi^{2}");
0176   h1cutFlow->GetXaxis()->SetBinLabel(5, "algo");
0177   h1cutFlow->GetXaxis()->SetBinLabel(6, "quality");
0178   h1cutFlow->GetXaxis()->SetBinLabel(7, "d_{0}");
0179   h1cutFlow->GetXaxis()->SetBinLabel(8, "z_{0}");
0180   h1cutFlow->GetXaxis()->SetBinLabel(9, "p_{T}");
0181   resetCutFlow();
0182 
0183   // Primary vertex fitter
0184   MyPVFitter = new PVFitter(iConfig, iColl);
0185   MyPVFitter->resetAll();
0186   if (savePVVertices_) {
0187     fPVTree_ = new TTree("PrimaryVertices", "PrimaryVertices");
0188     MyPVFitter->setTree(fPVTree_);
0189   }
0190 
0191   // check filename
0192   ffilename_changed = false;
0193   if (writeDIPTxt_)
0194     fasciiDIP.open(outputDIPTxt_.c_str());
0195 }
0196 
0197 BeamFitter::~BeamFitter() {
0198   if (saveNtuple_) {
0199     file_->cd();
0200     if (fitted_ && h1z)
0201       h1z->Write();
0202     h1ntrks->Write();
0203     h1vz_event->Write();
0204     if (h1cutFlow)
0205       h1cutFlow->Write();
0206     ftree_->Write();
0207   }
0208   if (saveBeamFit_) {
0209     file_->cd();
0210     ftreeFit_->Write();
0211   }
0212   if (savePVVertices_) {
0213     file_->cd();
0214     fPVTree_->Write();
0215   }
0216 
0217   if (saveNtuple_ || saveBeamFit_ || savePVVertices_) {
0218     file_->Close();
0219     delete file_;
0220   }
0221   delete MyPVFitter;
0222 }
0223 
0224 void BeamFitter::fillDescription(edm::ParameterSetDescription &iDesc) {
0225   edm::ParameterSetDescription beamFitter;
0226 
0227   beamFitter.addUntracked<bool>("Debug");
0228   beamFitter.addUntracked<edm::InputTag>("TrackCollection");
0229   iDesc.addUntracked<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0230   iDesc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0231   beamFitter.addUntracked<bool>("WriteAscii");
0232   beamFitter.addUntracked<std::string>("AsciiFileName");
0233   beamFitter.addUntracked<bool>("AppendRunToFileName");
0234   beamFitter.addUntracked<bool>("WriteDIPAscii");
0235   // Specify whether we want to write the DIP file even if the fit is failed.
0236   beamFitter.addUntracked<bool>("WriteDIPOnBadFit", true);
0237   beamFitter.addUntracked<std::string>("DIPFileName");
0238   beamFitter.addUntracked<bool>("SaveNtuple");
0239   beamFitter.addUntracked<bool>("SaveFitResults");
0240   beamFitter.addUntracked<bool>("SavePVVertices");
0241   beamFitter.addUntracked<bool>("IsMuonCollection");
0242 
0243   beamFitter.addUntracked<double>("MinimumPt");
0244   beamFitter.addUntracked<double>("MaximumEta");
0245   beamFitter.addUntracked<double>("MaximumImpactParameter");
0246   beamFitter.addUntracked<double>("MaximumZ");
0247   beamFitter.addUntracked<int>("MinimumTotalLayers");
0248   beamFitter.addUntracked<int>("MinimumPixelLayers");
0249   beamFitter.addUntracked<double>("MaximumNormChi2");
0250   beamFitter.addUntracked<std::vector<std::string> >("TrackAlgorithm");
0251   beamFitter.addUntracked<std::vector<std::string> >("TrackQuality");
0252   beamFitter.addUntracked<int>("MinimumInputTracks");
0253   beamFitter.addUntracked<double>("FractionOfFittedTrks");
0254   beamFitter.addUntracked<double>("InputBeamWidth", -1.);
0255 
0256   beamFitter.addUntracked<std::string>("OutputFileName", "");
0257 
0258   iDesc.add<edm::ParameterSetDescription>("BeamFitter", beamFitter);
0259 }
0260 
0261 void BeamFitter::readEvent(const edm::Event &iEvent) {
0262   frun = iEvent.id().run();
0263   const edm::TimeValue_t ftimestamp = iEvent.time().value();
0264   const std::time_t ftmptime = ftimestamp >> 32;
0265 
0266   if (fbeginLumiOfFit == -1)
0267     freftime[0] = freftime[1] = ftmptime;
0268   if (freftime[0] == 0 || ftmptime < freftime[0])
0269     freftime[0] = ftmptime;
0270   if (freftime[1] == 0 || ftmptime > freftime[1])
0271     freftime[1] = ftmptime;
0272   // Update the human-readable string versions of the time
0273   updateBTime();
0274 
0275   flumi = iEvent.luminosityBlock();
0276   frunFit = frun;
0277   if (fbeginLumiOfFit == -1 || fbeginLumiOfFit > flumi)
0278     fbeginLumiOfFit = flumi;
0279   if (fendLumiOfFit == -1 || fendLumiOfFit < flumi)
0280     fendLumiOfFit = flumi;
0281 
0282   edm::Handle<reco::TrackCollection> TrackCollection;
0283   iEvent.getByToken(tracksToken_, TrackCollection);
0284 
0285   //------ Primary Vertices
0286   edm::Handle<edm::View<reco::Vertex> > PVCollection;
0287   bool hasPVs = false;
0288   edm::View<reco::Vertex> pv;
0289   if (iEvent.getByToken(vertexToken_, PVCollection)) {
0290     pv = *PVCollection;
0291     hasPVs = true;
0292   }
0293   //------
0294 
0295   //------ Beam Spot in current event
0296   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0297   const reco::BeamSpot *refBS = nullptr;
0298   if (iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle))
0299     refBS = recoBeamSpotHandle.product();
0300   //-------
0301 
0302   const reco::TrackCollection *tracks = TrackCollection.product();
0303 
0304   double eventZ = 0;
0305   double averageZ = 0;
0306 
0307   for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
0308     if (!isMuon_) {
0309       const reco::HitPattern &trkHP = track->hitPattern();
0310 
0311       fnPixelLayerMeas = trkHP.pixelLayersWithMeasurement();
0312       fnStripLayerMeas = trkHP.stripLayersWithMeasurement();
0313       fnTotLayerMeas = trkHP.trackerLayersWithMeasurement();
0314       fnPXBLayerMeas = trkHP.pixelBarrelLayersWithMeasurement();
0315       fnPXFLayerMeas = trkHP.pixelEndcapLayersWithMeasurement();
0316       fnTIBLayerMeas = trkHP.stripTIBLayersWithMeasurement();
0317       fnTIDLayerMeas = trkHP.stripTIDLayersWithMeasurement();
0318       fnTOBLayerMeas = trkHP.stripTOBLayersWithMeasurement();
0319       fnTECLayerMeas = trkHP.stripTECLayersWithMeasurement();
0320     } else {
0321       fnTotLayerMeas = track->numberOfValidHits();
0322     }
0323 
0324     fpt = track->pt();
0325     feta = track->eta();
0326     fphi0 = track->phi();
0327     fcharge = track->charge();
0328     fnormchi2 = track->normalizedChi2();
0329     fd0 = track->d0();
0330     if (refBS)
0331       fd0bs = -1 * track->dxy(refBS->position());
0332     else
0333       fd0bs = 0.;
0334 
0335     fsigmad0 = track->d0Error();
0336     fz0 = track->dz();
0337     fsigmaz0 = track->dzError();
0338     ftheta = track->theta();
0339     fvx = track->vx();
0340     fvy = track->vy();
0341 
0342     for (int i = 0; i < 5; ++i) {
0343       for (int j = 0; j < 5; ++j) {
0344         fcov[i][j] = track->covariance(i, j);
0345       }
0346     }
0347 
0348     fquality = true;
0349     falgo = true;
0350 
0351     if (!isMuon_) {
0352       if (!quality_.empty()) {
0353         fquality = false;
0354         for (unsigned int i = 0; i < quality_.size(); ++i) {
0355           if (debug_)
0356             edm::LogInfo("BeamFitter") << "quality_[" << i << "] = " << track->qualityName(quality_[i]) << std::endl;
0357           if (track->quality(quality_[i])) {
0358             fquality = true;
0359             break;
0360           }
0361         }
0362       }
0363 
0364       // Track algorithm
0365 
0366       if (!algorithm_.empty()) {
0367         if (std::find(algorithm_.begin(), algorithm_.end(), track->algo()) == algorithm_.end())
0368           falgo = false;
0369       }
0370     }
0371 
0372     // check if we have a valid PV
0373     fpvValid = false;
0374 
0375     if (hasPVs) {
0376       for (size_t ipv = 0; ipv != pv.size(); ++ipv) {
0377         if (!pv[ipv].isFake())
0378           fpvValid = true;
0379 
0380         if (ipv == 0 && !pv[0].isFake()) {
0381           fpvx = pv[0].x();
0382           fpvy = pv[0].y();
0383           fpvz = pv[0].z();
0384         }  // fix this later
0385       }
0386     }
0387 
0388     if (saveNtuple_)
0389       ftree_->Fill();
0390     ftotal_tracks++;
0391 
0392     countPass[0] = ftotal_tracks;
0393     // Track selection
0394     if (fnTotLayerMeas >= trk_MinNTotLayers_) {
0395       countPass[1] += 1;
0396       if (fnPixelLayerMeas >= trk_MinNPixLayers_) {
0397         countPass[2] += 1;
0398         if (fnormchi2 < trk_MaxNormChi2_) {
0399           countPass[3] += 1;
0400           if (falgo) {
0401             countPass[4] += 1;
0402             if (fquality) {
0403               countPass[5] += 1;
0404               if (std::abs(fd0) < trk_MaxIP_) {
0405                 countPass[6] += 1;
0406                 if (std::abs(fz0) < trk_MaxZ_) {
0407                   countPass[7] += 1;
0408                   if (fpt > trk_MinpT_) {
0409                     countPass[8] += 1;
0410                     if (std::abs(feta) < trk_MaxEta_
0411                         //&& fpvValid
0412                     ) {
0413                       if (debug_) {
0414                         edm::LogInfo("BeamFitter") << "Selected track quality = " << track->qualityMask()
0415                                                    << "; track algorithm = " << track->algoName()
0416                                                    << "= TrackAlgorithm: " << track->algo() << std::endl;
0417                       }
0418                       BSTrkParameters BSTrk(fz0, fsigmaz0, fd0, fsigmad0, fphi0, fpt, 0., 0.);
0419                       BSTrk.setVx(fvx);
0420                       BSTrk.setVy(fvy);
0421                       fBSvector.push_back(BSTrk);
0422                       averageZ += fz0;
0423                     }
0424                   }
0425                 }
0426               }
0427             }
0428           }
0429         }
0430       }
0431     }  // track selection
0432 
0433   }  // tracks
0434 
0435   averageZ = averageZ / (float)(fBSvector.size());
0436 
0437   for (std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
0438     eventZ += fabs(iparam->z0() - averageZ);
0439   }
0440 
0441   h1ntrks->Fill(fBSvector.size());
0442   h1vz_event->Fill(eventZ / (float)(fBSvector.size()));
0443   for (unsigned int i = 0; i < sizeof(countPass) / sizeof(countPass[0]); i++)
0444     h1cutFlow->SetBinContent(i + 1, countPass[i]);
0445 
0446   MyPVFitter->readEvent(iEvent);
0447 }
0448 
0449 bool BeamFitter::runPVandTrkFitter() {
0450   // run both PV and track fitters
0451   bool fit_ok = false;
0452   bool pv_fit_ok = false;
0453   reco::BeamSpot bspotPV;
0454   reco::BeamSpot bspotTrk;
0455 
0456   // First run PV fitter
0457   if (MyPVFitter->IsFitPerBunchCrossing()) {
0458     edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[0] = " << freftime[0]
0459                                << "; address =  " << &freftime[0] << " = " << fbeginTimeOfFit << std::endl;
0460     edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[1] = " << freftime[1]
0461                                << "; address =  " << &freftime[1] << " = " << fendTimeOfFit << std::endl;
0462 
0463     if (MyPVFitter->runBXFitter()) {
0464       fbspotPVMap = MyPVFitter->getBeamSpotMap();
0465       pv_fit_ok = true;
0466     }
0467     if (writeTxt_)
0468       dumpTxtFile(outputTxt_, true);  // all reaults
0469     if (writeDIPTxt_ && (pv_fit_ok || writeDIPBadFit_)) {
0470       dumpTxtFile(outputDIPTxt_, false);  // for DQM/DIP
0471     }
0472     return pv_fit_ok;
0473   }
0474 
0475   if (MyPVFitter->runFitter()) {
0476     bspotPV = MyPVFitter->getBeamSpot();
0477 
0478     // take beam width from PV fit and pass it to track fitter
0479     // assume circular transverse beam width
0480     inputBeamWidth_ = sqrt(pow(bspotPV.BeamWidthX(), 2) + pow(bspotPV.BeamWidthY(), 2)) / sqrt(2);
0481     pv_fit_ok = true;
0482 
0483   } else {
0484     // problems with PV fit
0485     bspotPV.setType(reco::BeamSpot::Unknown);
0486     bspotTrk.setType(reco::BeamSpot::Unknown);  //propagate error to trk beam spot
0487   }
0488 
0489   if (runFitterNoTxt()) {
0490     bspotTrk = fbeamspot;
0491     fit_ok = true;
0492   } else {
0493     // beamfit failed, flagged as empty beam spot
0494     bspotTrk.setType(reco::BeamSpot::Fake);
0495     fit_ok = false;
0496   }
0497 
0498   // combined results into one single beam spot
0499 
0500   // to pass errors I have to create another beam spot object
0501   reco::BeamSpot::CovarianceMatrix matrix;
0502   for (int j = 0; j < 7; ++j) {
0503     for (int k = j; k < 7; ++k) {
0504       matrix(j, k) = bspotTrk.covariance(j, k);
0505     }
0506   }
0507   // change beam width error to one from PV
0508   if (pv_fit_ok && fit_ok) {
0509     matrix(6, 6) = MyPVFitter->getWidthXerr() * MyPVFitter->getWidthXerr();
0510 
0511     // get Z and sigmaZ from PV fit
0512     matrix(2, 2) = bspotPV.covariance(2, 2);
0513     matrix(3, 3) = bspotPV.covariance(3, 3);
0514     reco::BeamSpot tmpbs(reco::BeamSpot::Point(bspotTrk.x0(), bspotTrk.y0(), bspotPV.z0()),
0515                          bspotPV.sigmaZ(),
0516                          bspotTrk.dxdz(),
0517                          bspotTrk.dydz(),
0518                          bspotPV.BeamWidthX(),
0519                          matrix,
0520                          bspotPV.type());
0521     tmpbs.setBeamWidthY(bspotPV.BeamWidthY());
0522     // overwrite beam spot result
0523     fbeamspot = tmpbs;
0524   }
0525   if (pv_fit_ok && fit_ok) {
0526     fbeamspot.setType(bspotPV.type());
0527   } else if (!pv_fit_ok && fit_ok) {
0528     fbeamspot.setType(reco::BeamSpot::Unknown);
0529   } else if (pv_fit_ok && !fit_ok) {
0530     fbeamspot.setType(reco::BeamSpot::Unknown);
0531   } else if (!pv_fit_ok && !fit_ok) {
0532     fbeamspot.setType(reco::BeamSpot::Fake);
0533   }
0534 
0535   if (writeTxt_)
0536     dumpTxtFile(outputTxt_, true);  // all reaults
0537   if (writeDIPTxt_ && ((fit_ok && pv_fit_ok) || writeDIPBadFit_)) {
0538     dumpTxtFile(outputDIPTxt_, false);  // for DQM/DIP
0539     for (size_t i = 0; i < 7; i++)
0540       ForDIPPV_.push_back(0.0);
0541   }
0542 
0543   return fit_ok && pv_fit_ok;
0544 }
0545 
0546 bool BeamFitter::runFitterNoTxt() {
0547   edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[0] = " << freftime[0]
0548                              << "; address =  " << &freftime[0] << " = " << fbeginTimeOfFit << std::endl;
0549   edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[1] = " << freftime[1]
0550                              << "; address =  " << &freftime[1] << " = " << fendTimeOfFit << std::endl;
0551 
0552   if (fbeginLumiOfFit == -1 || fendLumiOfFit == -1) {
0553     edm::LogWarning("BeamFitter") << "No event read! No Fitting!" << std::endl;
0554     return false;
0555   }
0556 
0557   bool fit_ok = false;
0558   // default fit to extract beam spot info
0559   if (fBSvector.size() > 1) {
0560     edm::LogInfo("BeamFitter") << "Calculating beam spot..." << std::endl
0561                                << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks
0562                                << std::endl;
0563 
0564     BSFitter *myalgo = new BSFitter(fBSvector);
0565     myalgo->SetMaximumZ(trk_MaxZ_);
0566     myalgo->SetConvergence(convergence_);
0567     myalgo->SetMinimumNTrks(min_Ntrks_);
0568     if (inputBeamWidth_ > 0)
0569       myalgo->SetInputBeamWidth(inputBeamWidth_);
0570 
0571     fbeamspot = myalgo->Fit();
0572 
0573     // retrieve histogram for Vz
0574     h1z = (TH1F *)myalgo->GetVzHisto();
0575 
0576     delete myalgo;
0577     if (fbeamspot.type() > 0) {  // save all results except for Fake and Unknown (all 0.)
0578       fit_ok = true;
0579       if (saveBeamFit_) {
0580         fx = fbeamspot.x0();
0581         fy = fbeamspot.y0();
0582         fz = fbeamspot.z0();
0583         fsigmaZ = fbeamspot.sigmaZ();
0584         fdxdz = fbeamspot.dxdz();
0585         fdydz = fbeamspot.dydz();
0586         fwidthX = fbeamspot.BeamWidthX();
0587         fwidthY = fbeamspot.BeamWidthY();
0588         fxErr = fbeamspot.x0Error();
0589         fyErr = fbeamspot.y0Error();
0590         fzErr = fbeamspot.z0Error();
0591         fsigmaZErr = fbeamspot.sigmaZ0Error();
0592         fdxdzErr = fbeamspot.dxdzError();
0593         fdydzErr = fbeamspot.dydzError();
0594         fwidthXErr = fbeamspot.BeamWidthXError();
0595         fwidthYErr = fbeamspot.BeamWidthYError();
0596         ftreeFit_->Fill();
0597       }
0598     }
0599   } else {  // tracks <= 1
0600     reco::BeamSpot tmpbs;
0601     fbeamspot = tmpbs;
0602     fbeamspot.setType(reco::BeamSpot::Fake);
0603     edm::LogInfo("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;
0604   }
0605   fitted_ = true;
0606   return fit_ok;
0607 }
0608 
0609 bool BeamFitter::runFitter() {
0610   bool fit_ok = runFitterNoTxt();
0611 
0612   if (writeTxt_)
0613     dumpTxtFile(outputTxt_, true);  // all reaults
0614   if (writeDIPTxt_ && (fit_ok || writeDIPBadFit_)) {
0615     dumpTxtFile(outputDIPTxt_, false);  // for DQM/DIP
0616   }
0617   return fit_ok;
0618 }
0619 
0620 bool BeamFitter::runBeamWidthFitter() {
0621   bool widthfit_ok = false;
0622   // default fit to extract beam spot info
0623   if (fBSvector.size() > 1) {
0624     edm::LogInfo("BeamFitter") << "Calculating beam spot positions('d0-phi0' method) and width using llh Fit"
0625                                << std::endl
0626                                << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks
0627                                << std::endl;
0628 
0629     BSFitter *myalgo = new BSFitter(fBSvector);
0630     myalgo->SetMaximumZ(trk_MaxZ_);
0631     myalgo->SetConvergence(convergence_);
0632     myalgo->SetMinimumNTrks(min_Ntrks_);
0633     if (inputBeamWidth_ > 0)
0634       myalgo->SetInputBeamWidth(inputBeamWidth_);
0635 
0636     myalgo->SetFitVariable(std::string("d*z"));
0637     myalgo->SetFitType(std::string("likelihood"));
0638     fbeamWidthFit = myalgo->Fit();
0639 
0640     //Add to .txt file
0641     if (writeTxt_)
0642       dumpBWTxtFile(outputTxt_);
0643 
0644     delete myalgo;
0645 
0646     // not fake
0647     if (fbeamspot.type() != 0)
0648       widthfit_ok = true;
0649   } else {
0650     fbeamspot.setType(reco::BeamSpot::Fake);
0651     edm::LogWarning("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;
0652   }
0653   return widthfit_ok;
0654 }
0655 
0656 void BeamFitter::dumpBWTxtFile(std::string &fileName) {
0657   std::ofstream outFile;
0658   outFile.open(fileName.c_str(), std::ios::app);
0659   outFile << "---------------------------------------------------------------------------------------------------------"
0660              "----------------------------------------------------"
0661           << std::endl;
0662   outFile << "Beam width(in cm) from Log-likelihood fit (Here we assume a symmetric beam(SigmaX=SigmaY)!)" << std::endl;
0663   outFile << "   " << std::endl;
0664   outFile << "BeamWidth =  " << fbeamWidthFit.BeamWidthX() << " +/- " << fbeamWidthFit.BeamWidthXError() << std::endl;
0665   outFile.close();
0666 }
0667 
0668 void BeamFitter::dumpTxtFile(std::string &fileName, bool append) {
0669   std::ofstream outFile;
0670 
0671   std::string tmpname = outputTxt_;
0672   char index[15];
0673   if (appendRunTxt_ && writeTxt_ && !ffilename_changed) {
0674     sprintf(index, "%s%i", "_Run", frun);
0675     tmpname.insert(outputTxt_.length() - 4, index);
0676     fileName = tmpname;
0677     ffilename_changed = true;
0678   }
0679 
0680   if (!append)
0681     outFile.open(fileName.c_str());
0682   else
0683     outFile.open(fileName.c_str(), std::ios::app);
0684 
0685   if (MyPVFitter->IsFitPerBunchCrossing()) {
0686     for (std::map<int, reco::BeamSpot>::const_iterator abspot = fbspotPVMap.begin(); abspot != fbspotPVMap.end();
0687          ++abspot) {
0688       reco::BeamSpot beamspottmp = abspot->second;
0689       int bx = abspot->first;
0690 
0691       outFile << "Runnumber " << frun << " bx " << bx << std::endl;
0692       outFile << "BeginTimeOfFit " << fbeginTimeOfFit << " " << freftime[0] << std::endl;
0693       outFile << "EndTimeOfFit " << fendTimeOfFit << " " << freftime[1] << std::endl;
0694       outFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
0695       outFile << "Type " << beamspottmp.type() << std::endl;
0696       outFile << "X0 " << beamspottmp.x0() << std::endl;
0697       outFile << "Y0 " << beamspottmp.y0() << std::endl;
0698       outFile << "Z0 " << beamspottmp.z0() << std::endl;
0699       outFile << "sigmaZ0 " << beamspottmp.sigmaZ() << std::endl;
0700       outFile << "dxdz " << beamspottmp.dxdz() << std::endl;
0701       outFile << "dydz " << beamspottmp.dydz() << std::endl;
0702       outFile << "BeamWidthX " << beamspottmp.BeamWidthX() << std::endl;
0703       outFile << "BeamWidthY " << beamspottmp.BeamWidthY() << std::endl;
0704       for (int i = 0; i < 6; ++i) {
0705         outFile << "Cov(" << i << ",j) ";
0706         for (int j = 0; j < 7; ++j) {
0707           outFile << beamspottmp.covariance(i, j) << " ";
0708         }
0709         outFile << std::endl;
0710       }
0711       outFile << "Cov(6,j) 0 0 0 0 0 0 " << beamspottmp.covariance(6, 6) << std::endl;
0712       //}
0713       outFile << "EmittanceX " << beamspottmp.emittanceX() << std::endl;
0714       outFile << "EmittanceY " << beamspottmp.emittanceY() << std::endl;
0715       outFile << "BetaStar " << beamspottmp.betaStar() << std::endl;
0716     }
0717   }  //if bx results needed
0718   else {
0719     beamspot::BeamSpotContainer currentBS;
0720 
0721     currentBS.beamspot = fbeamspot;
0722     currentBS.run = frun;
0723     std::copy(fbeginTimeOfFit, fbeginTimeOfFit + 32, currentBS.beginTimeOfFit);
0724     std::copy(fendTimeOfFit, fendTimeOfFit + 32, currentBS.endTimeOfFit);
0725     currentBS.beginLumiOfFit = fbeginLumiOfFit;
0726     currentBS.endLumiOfFit = fendLumiOfFit;
0727     std::copy(freftime, freftime + 2, currentBS.reftime);
0728 
0729     beamspot::dumpBeamSpotTxt(outFile, currentBS);
0730 
0731     //write here Pv info for DIP only: This added only if append is false, which happen for DIP only :)
0732     if (!append) {
0733       outFile << "events " << (int)ForDIPPV_[0] << std::endl;
0734       outFile << "meanPV " << ForDIPPV_[1] << std::endl;
0735       outFile << "meanErrPV " << ForDIPPV_[2] << std::endl;
0736       outFile << "rmsPV " << ForDIPPV_[3] << std::endl;
0737       outFile << "rmsErrPV " << ForDIPPV_[4] << std::endl;
0738       outFile << "maxPV " << (int)ForDIPPV_[5] << std::endl;
0739       outFile << "nPV " << (int)ForDIPPV_[6] << std::endl;
0740     }  //writeDIPPVInfo_
0741   }    //else end  here
0742 
0743   outFile.close();
0744 }
0745 
0746 void BeamFitter::write2DB() {
0747   BeamSpotObjects pBSObjects;
0748 
0749   pBSObjects.setPosition(fbeamspot.position().X(), fbeamspot.position().Y(), fbeamspot.position().Z());
0750   //std::cout << " wrote: x= " << fbeamspot.position().X() << " y= "<< fbeamspot.position().Y() << " z= " << fbeamspot.position().Z() << std::endl;
0751   pBSObjects.setSigmaZ(fbeamspot.sigmaZ());
0752   pBSObjects.setdxdz(fbeamspot.dxdz());
0753   pBSObjects.setdydz(fbeamspot.dydz());
0754   //if (inputBeamWidth_ > 0 ) {
0755   //  std::cout << " beam width value forced to be " << inputBeamWidth_ << std::endl;
0756   //  pBSObjects->SetBeamWidthX(inputBeamWidth_);
0757   //  pBSObjects->SetBeamWidthY(inputBeamWidth_);
0758   //} else {
0759   // need to fix this
0760   //std::cout << " using default value, 15e-4, for beam width!!!"<<std::endl;
0761   pBSObjects.setBeamWidthX(fbeamspot.BeamWidthX());
0762   pBSObjects.setBeamWidthY(fbeamspot.BeamWidthY());
0763   //}
0764 
0765   for (int i = 0; i < 7; ++i) {
0766     for (int j = 0; j < 7; ++j) {
0767       pBSObjects.setCovariance(i, j, fbeamspot.covariance(i, j));
0768     }
0769   }
0770   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0771   if (poolDbService.isAvailable()) {
0772     std::cout << "poolDBService available" << std::endl;
0773     if (poolDbService->isNewTagRequest("BeamSpotObjectsRcd")) {
0774       std::cout << "new tag requested" << std::endl;
0775       poolDbService->createOneIOV<BeamSpotObjects>(pBSObjects, poolDbService->beginOfTime(), "BeamSpotObjectsRcd");
0776     } else {
0777       std::cout << "no new tag requested" << std::endl;
0778       poolDbService->appendOneIOV<BeamSpotObjects>(pBSObjects, poolDbService->currentTime(), "BeamSpotObjectsRcd");
0779     }
0780   }
0781 }
0782 
0783 void BeamFitter::runAllFitter() {
0784   if (!fBSvector.empty()) {
0785     BSFitter *myalgo = new BSFitter(fBSvector);
0786     fbeamspot = myalgo->Fit_d0phi();
0787 
0788     // iterative
0789     if (debug_)
0790       std::cout << " d0-phi Iterative:" << std::endl;
0791     BSFitter *myitealgo = new BSFitter(fBSvector);
0792     myitealgo->Setd0Cut_d0phi(4.0);
0793     reco::BeamSpot beam_ite = myitealgo->Fit_ited0phi();
0794     if (debug_) {
0795       std::cout << beam_ite << std::endl;
0796       std::cout << "\n Now run tests of the different fits\n";
0797     }
0798     // from here are just tests
0799     std::string fit_type = "chi2";
0800     myalgo->SetFitVariable(std::string("z"));
0801     myalgo->SetFitType(std::string("chi2"));
0802     reco::BeamSpot beam_fit_z_chi2 = myalgo->Fit();
0803     if (debug_) {
0804       std::cout << " z Chi2 Fit ONLY:" << std::endl;
0805       std::cout << beam_fit_z_chi2 << std::endl;
0806     }
0807 
0808     fit_type = "combined";
0809     myalgo->SetFitVariable("z");
0810     myalgo->SetFitType("combined");
0811     reco::BeamSpot beam_fit_z_lh = myalgo->Fit();
0812     if (debug_) {
0813       std::cout << " z Log-Likelihood Fit ONLY:" << std::endl;
0814       std::cout << beam_fit_z_lh << std::endl;
0815     }
0816 
0817     myalgo->SetFitVariable("d");
0818     myalgo->SetFitType("d0phi");
0819     reco::BeamSpot beam_fit_dphi = myalgo->Fit();
0820     if (debug_) {
0821       std::cout << " d0-phi0 Fit ONLY:" << std::endl;
0822       std::cout << beam_fit_dphi << std::endl;
0823     }
0824     /*
0825     myalgo->SetFitVariable(std::string("d*z"));
0826     myalgo->SetFitType(std::string("likelihood"));
0827     reco::BeamSpot beam_fit_dz_lh = myalgo->Fit();
0828     if (debug_){
0829       std::cout << " Log-Likelihood Fit:" << std::endl;
0830       std::cout << beam_fit_dz_lh << std::endl;
0831     }
0832 
0833     myalgo->SetFitVariable(std::string("d*z"));
0834     myalgo->SetFitType(std::string("resolution"));
0835     reco::BeamSpot beam_fit_dresz_lh = myalgo->Fit();
0836     if(debug_){
0837       std::cout << " IP Resolution Fit" << std::endl;
0838       std::cout << beam_fit_dresz_lh << std::endl;
0839 
0840       std::cout << "c0 = " << myalgo->GetResPar0() << " +- " << myalgo->GetResPar0Err() << std::endl;
0841       std::cout << "c1 = " << myalgo->GetResPar1() << " +- " << myalgo->GetResPar1Err() << std::endl;
0842     }
0843     */
0844   } else if (debug_)
0845     std::cout << "No good track selected! No beam fit!" << std::endl;
0846 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!