File indexing completed on 2021-02-14 12:45:58
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.h"
0021
0022
0023
0024 #include "DataFormats/TrackReco/interface/Track.h"
0025 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0026 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0027 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0028
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0033 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0034
0035
0036
0037
0038 ValidationMisalignedTracker::ValidationMisalignedTracker(const edm::ParameterSet& iConfig)
0039 : geomToken_(esConsumes()), magFieldToken_(esConsumes()) {
0040
0041 mzmu = 0., recmzmu = 0., ptzmu = 0., recptzmu = 0., etazmu = 0., recetazmu = 0., thetazmu = 0., recthetazmu = 0.,
0042 phizmu = 0., recphizmu = 0.;
0043 recenezmu = 0., enezmu = 0., pLzmu = 0., recpLzmu = 0., yzmu = 0., recyzmu = 0., mxptmu = 0., recmxptmu = 0.,
0044 minptmu = 0., recminptmu = 0.;
0045
0046
0047 flag = 0, flagrec = 0, count = 0, countrec = 0;
0048 nAssoc = 0;
0049
0050 for (int i = 0; i < 2; i++) {
0051 countpart[i] = 0;
0052 countpartrec[i] = 0;
0053 for (int j = 0; j < 2; j++) {
0054 ene[i][j] = 0.;
0055 p[i][j] = 0.;
0056 px[i][j] = 0.;
0057 py[i][j] = 0.;
0058 pz[i][j] = 0.;
0059 ptmu[i][j] = 0.;
0060 recene[i][j] = 0.;
0061 recp[i][j] = 0.;
0062 recpx[i][j] = 0.;
0063 recpy[i][j] = 0.;
0064 recpz[i][j] = 0.;
0065 recptmu[i][j] = 0.;
0066 }
0067 }
0068
0069 eventCount_ = 0;
0070
0071 selection_eff = iConfig.getUntrackedParameter<bool>("selection_eff", false);
0072 selection_fake = iConfig.getUntrackedParameter<bool>("selection_fake", true);
0073 ZmassSelection_ = iConfig.getUntrackedParameter<bool>("ZmassSelection", false);
0074 simobject = iConfig.getUntrackedParameter<std::string>("simobject", "g4SimHits");
0075 trackassociator = iConfig.getUntrackedParameter<std::string>("TrackAssociator", "ByHits");
0076 associators = iConfig.getParameter<std::vector<std::string> >("associators");
0077 label = iConfig.getParameter<std::vector<edm::InputTag> >("label");
0078 label_tp_effic = iConfig.getParameter<edm::InputTag>("label_tp_effic");
0079 label_tp_fake = iConfig.getParameter<edm::InputTag>("label_tp_fake");
0080
0081 rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile", "myroot.root");
0082 file_ = new TFile(rootfile_.c_str(), "RECREATE");
0083
0084
0085 tree_eff = new TTree("EffTracks", "Efficiency Tracks Tree");
0086
0087 tree_eff->Branch("Run", &irun, "irun/i");
0088 tree_eff->Branch("Event", &ievt, "ievt/i");
0089
0090
0091 tree_eff->Branch("TrackID", &trackType, "trackType/i");
0092 tree_eff->Branch("pt", &pt, "pt/F");
0093 tree_eff->Branch("eta", &eta, "eta/F");
0094 tree_eff->Branch("CotTheta", &cottheta, "cottheta/F");
0095 tree_eff->Branch("phi", &phi, "phi/F");
0096 tree_eff->Branch("d0", &d0, "d0/F");
0097 tree_eff->Branch("z0", &z0, "z0/F");
0098 tree_eff->Branch("nhit", &nhit, "nhit/i");
0099
0100
0101 tree_eff->Branch("recpt", &recpt, "recpt/F");
0102 tree_eff->Branch("receta", &receta, "receta/F");
0103 tree_eff->Branch("CotRecTheta", &reccottheta, "reccottheta/F");
0104 tree_eff->Branch("recphi", &recphi, "recphi/F");
0105 tree_eff->Branch("recd0", &recd0, "recd0/F");
0106 tree_eff->Branch("recz0", &recz0, "recz0/F");
0107 tree_eff->Branch("nAssoc", &nAssoc, "nAssoc/i");
0108 tree_eff->Branch("recnhit", &recnhit, "recnhit/i");
0109 tree_eff->Branch("CHISQ", &recchiq, "recchiq/F");
0110
0111 tree_eff->Branch("reseta", &reseta, "reseta/F");
0112 tree_eff->Branch("respt", &respt, "respt/F");
0113 tree_eff->Branch("resd0", &resd0, "resd0/F");
0114 tree_eff->Branch("resz0", &resz0, "resz0/F");
0115 tree_eff->Branch("resphi", &resphi, "resphi/F");
0116 tree_eff->Branch("rescottheta", &rescottheta, "rescottheta/F");
0117 tree_eff->Branch("eff", &eff, "eff/F");
0118
0119
0120 tree_eff->Branch("mzmu", &mzmu, "mzmu/F");
0121 tree_eff->Branch("ptzmu", &ptzmu, "ptzmu/F");
0122 tree_eff->Branch("pLzmu", &pLzmu, "pLzmu/F");
0123 tree_eff->Branch("enezmu", &enezmu, "enezmu/F");
0124 tree_eff->Branch("etazmu", &etazmu, "etazmu/F");
0125 tree_eff->Branch("thetazmu", &thetazmu, "thetazmu/F");
0126 tree_eff->Branch("phizmu", &phizmu, "phizmu/F");
0127 tree_eff->Branch("yzmu", &yzmu, "yzmu/F");
0128 tree_eff->Branch("mxptmu", &mxptmu, "mxptmu/F");
0129 tree_eff->Branch("minptmu", &minptmu, "minptmu/F");
0130
0131 tree_eff->Branch("recmzmu", &recmzmu, "recmzmu/F");
0132 tree_eff->Branch("recptzmu", &recptzmu, "recptzmu/F");
0133 tree_eff->Branch("recpLzmu", &recpLzmu, "recpLzmu/F");
0134 tree_eff->Branch("recenezmu", &recenezmu, "recenezmu/F");
0135 tree_eff->Branch("recetazmu", &recetazmu, "recetazmu/F");
0136 tree_eff->Branch("recthetazmu", &recthetazmu, "recthetazmu/F");
0137 tree_eff->Branch("recphizmu", &recphizmu, "recphizmu/F");
0138 tree_eff->Branch("recyzmu", &recyzmu, "recyzmu/F");
0139 tree_eff->Branch("recmxptmu", &recmxptmu, "recmxptmu/F");
0140 tree_eff->Branch("recminptmu", &recminptmu, "recminptmu/F");
0141
0142
0143
0144 tree_eff->Branch("chi2Associator", &recchiq, "recchiq/F");
0145
0146
0147
0148 tree_fake = new TTree("FakeTracks", "Fake Rate Tracks Tree");
0149
0150 tree_fake->Branch("Run", &irun, "irun/i");
0151 tree_fake->Branch("Event", &ievt, "ievt/i");
0152
0153
0154 tree_fake->Branch("fakeTrackID", &faketrackType, "faketrackType/i");
0155 tree_fake->Branch("fakept", &fakept, "fakept/F");
0156 tree_fake->Branch("fakeeta", &fakeeta, "fakeeta/F");
0157 tree_fake->Branch("fakeCotTheta", &fakecottheta, "fakecottheta/F");
0158 tree_fake->Branch("fakephi", &fakephi, "fakephi/F");
0159 tree_fake->Branch("faked0", &faked0, "faked0/F");
0160 tree_fake->Branch("fakez0", &fakez0, "fakez0/F");
0161 tree_fake->Branch("fakenhit", &fakenhit, "fakenhit/i");
0162
0163
0164 tree_fake->Branch("fakerecpt", &fakerecpt, "fakerecpt/F");
0165 tree_fake->Branch("fakereceta", &fakereceta, "fakereceta/F");
0166 tree_fake->Branch("fakeCotRecTheta", &fakereccottheta, "fakereccottheta/F");
0167 tree_fake->Branch("fakerecphi", &fakerecphi, "fakerecphi/F");
0168 tree_fake->Branch("fakerecd0", &fakerecd0, "fakerecd0/F");
0169 tree_fake->Branch("fakerecz0", &fakerecz0, "fakerecz0/F");
0170 tree_fake->Branch("fakenAssoc", &fakenAssoc, "fakenAssoc/i");
0171 tree_fake->Branch("fakerecnhit", &fakerecnhit, "fakerecnhit/i");
0172 tree_fake->Branch("fakeCHISQ", &fakerecchiq, "fakerecchiq/F");
0173
0174 tree_fake->Branch("fakereseta", &fakereseta, "fakereseta/F");
0175 tree_fake->Branch("fakerespt", &fakerespt, "fakerespt/F");
0176 tree_fake->Branch("fakeresd0", &fakeresd0, "fakeresd0/F");
0177 tree_fake->Branch("fakeresz0", &fakeresz0, "fakeresz0/F");
0178 tree_fake->Branch("fakeresphi", &fakeresphi, "fakeresphi/F");
0179 tree_fake->Branch("fakerescottheta", &fakerescottheta, "fakerescottheta/F");
0180 tree_fake->Branch("fake", &fake, "fake/F");
0181
0182
0183 tree_fake->Branch("fakemzmu", &fakemzmu, "fakemzmu/F");
0184 tree_fake->Branch("fakeptzmu", &fakeptzmu, "fakeptzmu/F");
0185 tree_fake->Branch("fakepLzmu", &fakepLzmu, "fakepLzmu/F");
0186 tree_fake->Branch("fakeenezmu", &fakeenezmu, "fakeenezmu/F");
0187 tree_fake->Branch("fakeetazmu", &fakeetazmu, "fakeetazmu/F");
0188 tree_fake->Branch("fakethetazmu", &fakethetazmu, "fakethetazmu/F");
0189 tree_fake->Branch("fakephizmu", &fakephizmu, "fakephizmu/F");
0190 tree_fake->Branch("fakeyzmu", &fakeyzmu, "fakeyzmu/F");
0191 tree_fake->Branch("fakemxptmu", &fakemxptmu, "fakemxptmu/F");
0192 tree_fake->Branch("fakeminptmu", &fakeminptmu, "fakeminptmu/F");
0193
0194 tree_fake->Branch("fakerecmzmu", &fakerecmzmu, "fakerecmzmu/F");
0195 tree_fake->Branch("fakerecptzmu", &fakerecptzmu, "fakerecptzmu/F");
0196 tree_fake->Branch("fakerecpLzmu", &fakerecpLzmu, "fakerecpLzmu/F");
0197 tree_fake->Branch("fakerecenezmu", &fakerecenezmu, "fakerecenezmu/F");
0198 tree_fake->Branch("fakerecetazmu", &fakerecetazmu, "fakerecetazmu/F");
0199 tree_fake->Branch("fakerecthetazmu", &fakerecthetazmu, "fakerecthetazmu/F");
0200 tree_fake->Branch("fakerecphizmu", &fakerecphizmu, "fakerecphizmu/F");
0201 tree_fake->Branch("fakerecyzmu", &fakerecyzmu, "fakerecyzmu/F");
0202 tree_fake->Branch("fakerecmxptmu", &fakerecmxptmu, "fakerecmxptmu/F");
0203 tree_fake->Branch("fakerecminptmu", &fakerecminptmu, "fakerecminptmu/F");
0204
0205 tree_fake->Branch("fakechi2Associator", &fakerecchiq, "fakerecchiq/F");
0206 }
0207
0208 ValidationMisalignedTracker::~ValidationMisalignedTracker() {
0209 std::cout << "ValidationMisalignedTracker::endJob Processed " << eventCount_ << " events" << std::endl;
0210
0211
0212 file_->Write();
0213
0214
0215 file_->Close();
0216 tree_eff = nullptr;
0217 tree_fake = nullptr;
0218 }
0219
0220
0221
0222
0223
0224
0225 void ValidationMisalignedTracker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0226 std::vector<const reco::TrackToTrackingParticleAssociator*> associatore;
0227
0228 {
0229 edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
0230 for (unsigned int w = 0; w < associators.size(); w++) {
0231 iEvent.getByLabel(associators[w], theAssociator);
0232 associatore.push_back(theAssociator.product());
0233 }
0234 }
0235
0236 edm::LogInfo("Tracker Misalignment Validation") << "\n Starting!";
0237
0238
0239 skip = false;
0240 std::vector<int> indmu;
0241
0242 if (selection_eff && ZmassSelection_) {
0243 edm::Handle<edm::HepMCProduct> evt;
0244 iEvent.getByLabel("source", evt);
0245 bool accepted = false;
0246 bool foundmuons = false;
0247 HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0248
0249 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0250 if (!accepted && ((*p)->pdg_id() == 23) && (*p)->status() == 3) {
0251 accepted = true;
0252 for (HepMC::GenVertex::particle_iterator aDaughter = (*p)->end_vertex()->particles_begin(HepMC::descendants);
0253 aDaughter != (*p)->end_vertex()->particles_end(HepMC::descendants);
0254 aDaughter++) {
0255 if (abs((*aDaughter)->pdg_id()) == 13) {
0256 foundmuons = true;
0257 if ((*aDaughter)->status() != 1) {
0258 for (HepMC::GenVertex::particle_iterator byaDaughter =
0259 (*aDaughter)->end_vertex()->particles_begin(HepMC::descendants);
0260 byaDaughter != (*aDaughter)->end_vertex()->particles_end(HepMC::descendants);
0261 byaDaughter++) {
0262 if ((*byaDaughter)->status() == 1 && abs((*byaDaughter)->pdg_id()) == 13) {
0263 indmu.push_back((*byaDaughter)->barcode());
0264 std::cout << "Stable muon from Z with charge " << (*byaDaughter)->pdg_id() << " and index "
0265 << (*byaDaughter)->barcode() << std::endl;
0266 }
0267 }
0268 } else {
0269 indmu.push_back((*aDaughter)->barcode());
0270 std::cout << "Stable muon from Z with charge " << (*aDaughter)->pdg_id() << " and index "
0271 << (*aDaughter)->barcode() << std::endl;
0272 }
0273 }
0274 }
0275 if (!foundmuons) {
0276 std::cout << "No muons from Z ...skip event" << std::endl;
0277 skip = true;
0278 }
0279 }
0280 }
0281 if (!accepted) {
0282 std::cout << "No Z particles in the event ...skip event" << std::endl;
0283 skip = true;
0284 }
0285 } else {
0286 skip = false;
0287 }
0288
0289
0290
0291
0292 const TrackerGeometry* trackerGeometry = &iSetup.getData(geomToken_);
0293 auto testGeomDet = trackerGeometry->detsTOB().front();
0294 std::cout << testGeomDet->position() << std::endl;
0295
0296
0297 irun = iEvent.id().run();
0298 ievt = iEvent.id().event();
0299
0300
0301 int countpart[2] = {0, 0}, countpartrec[2] = {0, 0}, flag = 0, flagrec = 0, count = 0, countrec = 0;
0302
0303 float ene[2][2], px[2][2], py[2][2], pz[2][2], ptmu[2][2];
0304 float recene[2][2], recp[2][2], recpx[2][2], recpy[2][2], recpz[2][2], recptmu[2][2];
0305
0306 for (int i = 0; i < 2; i++) {
0307 for (int j = 0; j < 2; j++) {
0308 ene[i][j] = 0.;
0309 px[i][j] = 0.;
0310 py[i][j] = 0.;
0311 pz[i][j] = 0.;
0312 ptmu[i][j] = 0.;
0313 recene[i][j] = 0.;
0314 recp[i][j] = 0.;
0315 recpx[i][j] = 0.;
0316 recpy[i][j] = 0.;
0317 recpz[i][j] = 0.;
0318 recptmu[i][j] = 0.;
0319 }
0320 }
0321
0322 edm::Handle<TrackingParticleCollection> TPCollectionHeff;
0323 iEvent.getByLabel(label_tp_effic, TPCollectionHeff);
0324 const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
0325
0326 edm::Handle<TrackingParticleCollection> TPCollectionHfake;
0327 iEvent.getByLabel(label_tp_fake, TPCollectionHfake);
0328 const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
0329
0330 int w = 0;
0331 for (unsigned int ww = 0; ww < associators.size(); ww++) {
0332
0333
0334
0335
0336 edm::InputTag algo = label[0];
0337
0338 edm::Handle<edm::View<reco::Track> > trackCollection;
0339 iEvent.getByLabel(algo, trackCollection);
0340 const edm::View<reco::Track> tC = *(trackCollection.product());
0341
0342
0343 LogTrace("TrackValidator") << "Calling associateRecoToSim method"
0344 << "\n";
0345 reco::RecoToSimCollection recSimColl = associatore[ww]->associateRecoToSim(trackCollection, TPCollectionHfake);
0346
0347 LogTrace("TrackValidator") << "Calling associateSimToReco method"
0348 << "\n";
0349 reco::SimToRecoCollection simRecColl = associatore[ww]->associateSimToReco(trackCollection, TPCollectionHeff);
0350
0351
0352
0353
0354
0355 if (selection_eff && !skip) {
0356 std::cout << "Computing Efficiency" << std::endl;
0357
0358 edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles (before cuts): " << tPCeff.size() << "\n";
0359 int ats = 0;
0360 int st = 0;
0361 for (TrackingParticleCollection::size_type i = 0; i < tPCeff.size(); i++) {
0362
0363 eta = 0., theta = 0., phi = 0., pt = 0., cottheta = 0., costheta = 0.;
0364 d0 = 0., z0 = 0.;
0365 nhit = 0;
0366 receta = 0., rectheta = 0., recphi = 0., recpt = 0., reccottheta = 0., recd0 = 0., recz0 = 0.;
0367 respt = 0., resd0 = 0., resz0 = 0., reseta = 0., resphi = 0., rescottheta = 0.;
0368 recchiq = 0.;
0369 recnhit = 0;
0370 trackType = 0;
0371 eff = 0;
0372
0373
0374 TrackingParticleRef tp(TPCollectionHeff, i);
0375 if (tp->charge() == 0)
0376 continue;
0377 st++;
0378
0379
0380
0381
0382 const SimTrack* simulatedTrack = &(*tp->g4Track_begin());
0383
0384 const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0385 FreeTrajectoryState ftsAtProduction(
0386 GlobalPoint(tp->vertex().x(), tp->vertex().y(), tp->vertex().z()),
0387 GlobalVector(
0388 simulatedTrack->momentum().x(), simulatedTrack->momentum().y(), simulatedTrack->momentum().z()),
0389 TrackCharge(tp->charge()),
0390 theMF);
0391 TSCPBuilderNoMaterial tscpBuilder;
0392 TrajectoryStateClosestToPoint tsAtClosestApproach =
0393 tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));
0394 GlobalPoint v = tsAtClosestApproach.theState().position();
0395 GlobalVector p = tsAtClosestApproach.theState().momentum();
0396
0397
0398
0399
0400 double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0401 double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0402 d0 = float(-dxySim);
0403 z0 = float(dszSim * p.mag() / p.perp());
0404
0405 if (abs(simulatedTrack->type()) == 13 && simulatedTrack->genpartIndex() != -1) {
0406 std::cout << " TRACCIA SIM DI MUONI " << std::endl;
0407 std::cout << "Gen part " << simulatedTrack->genpartIndex() << std::endl;
0408 trackType = simulatedTrack->type();
0409 theta = simulatedTrack->momentum().theta();
0410 costheta = cos(theta);
0411 cottheta = 1. / tan(theta);
0412
0413 eta = simulatedTrack->momentum().eta();
0414 phi = simulatedTrack->momentum().phi();
0415 pt = simulatedTrack->momentum().pt();
0416 nhit = tp->matchedHit();
0417
0418 std::cout << "3) Before assoc: SimTrack of type = " << simulatedTrack->type() << " ,at eta = " << eta
0419 << " ,with pt at vertex = " << simulatedTrack->momentum().pt() << " GeV/c"
0420 << " ,d0 =" << d0 << " ,z0 =" << z0 << " ,nhit=" << nhit << std::endl;
0421
0422 if (ZmassSelection_) {
0423 if (abs(trackType) == 13 &&
0424 (simulatedTrack->genpartIndex() == indmu[0] || simulatedTrack->genpartIndex() == indmu[1])) {
0425 std::cout << " TRACK sim of muons from Z " << std::endl;
0426 flag = 0;
0427 count = countpart[0];
0428 countpart[0]++;
0429 } else if (abs(trackType) == 11) {
0430
0431 flag = 1;
0432 count = countpart[1];
0433 countpart[1]++;
0434 }
0435
0436 px[flag][count] = simulatedTrack->momentum().x();
0437 py[flag][count] = simulatedTrack->momentum().y();
0438 pz[flag][count] = simulatedTrack->momentum().z();
0439 ptmu[flag][count] = simulatedTrack->momentum().pt();
0440 ene[flag][count] = simulatedTrack->momentum().e();
0441 }
0442
0443 std::vector<std::pair<edm::RefToBase<reco::Track>, double> > rt;
0444 if (simRecColl.find(tp) != simRecColl.end()) {
0445 rt = simRecColl[tp];
0446 if (!rt.empty()) {
0447 edm::RefToBase<reco::Track> t = rt.begin()->first;
0448 ats++;
0449
0450
0451
0452
0453
0454
0455
0456
0457 edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st << " with pt=" << t->pt()
0458 << " associated with quality:" << rt.begin()->second << "\n";
0459 std::cout << "Reconstructed Track:" << t->pt() << std::endl;
0460 std::cout << "\tpT: " << t->pt() << std::endl;
0461 std::cout << "\timpact parameter:d0: " << t->d0() << std::endl;
0462 std::cout << "\timpact parameter:z0: " << t->dz() << std::endl;
0463 std::cout << "\tAzimuthal angle of point of closest approach:" << t->phi() << std::endl;
0464 std::cout << "\tcharge: " << t->charge() << std::endl;
0465 std::cout << "\teta: " << t->eta() << std::endl;
0466 std::cout << "\tnormalizedChi2: " << t->normalizedChi2() << std::endl;
0467
0468 recnhit = t->numberOfValidHits();
0469 recchiq = t->normalizedChi2();
0470 rectheta = t->theta();
0471 reccottheta = 1. / tan(rectheta);
0472
0473 receta = t->momentum().eta();
0474
0475 recphi = t->phi();
0476 recpt = t->pt();
0477 ptused.push_back(recpt);
0478 recd0 = t->d0();
0479 recz0 = t->dz();
0480
0481 std::cout << "5) After call to associator: the best match has " << recnhit << " hits, Chi2 = " << recchiq
0482 << ", pt at vertex = " << recpt << " GeV/c, "
0483 << ", recd0 = " << recd0 << ", recz0= " << recz0 << std::endl;
0484
0485 respt = recpt - pt;
0486 resd0 = recd0 - d0;
0487 resz0 = recz0 - z0;
0488 reseta = receta - eta;
0489 resphi = recphi - phi;
0490 rescottheta = reccottheta - cottheta;
0491 eff = 1;
0492
0493 std::cout << "6) Transverse momentum residual=" << respt << " ,d0 residual=" << resd0
0494 << " ,z0 residual=" << resz0 << " with eff=" << eff << std::endl;
0495
0496 if (ZmassSelection_) {
0497 if (abs(trackType) == 13) {
0498 std::cout << " TRACCIA RECO DI MUONI " << std::endl;
0499 flagrec = 0;
0500 countrec = countpartrec[0];
0501 countpartrec[0]++;
0502 } else if (abs(trackType) == 11) {
0503 std::cout << " TRACCIA RECO DI ELETTRONI " << std::endl;
0504 flagrec = 1;
0505 countrec = countpartrec[1];
0506 countpartrec[1]++;
0507 }
0508
0509 recp[flagrec][countrec] = sqrt(t->momentum().mag2());
0510 recpx[flagrec][countrec] = t->momentum().x();
0511 recpy[flagrec][countrec] = t->momentum().y();
0512 recpz[flagrec][countrec] = t->momentum().z();
0513 recptmu[flagrec][countrec] =
0514 sqrt((t->momentum().x() * t->momentum().x()) + (t->momentum().y() * t->momentum().y()));
0515 if (abs(trackType) == 13)
0516 recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.105 * 0.105);
0517 if (abs(trackType) == 11)
0518 recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.0005 * 0.0005);
0519 }
0520
0521 std::cout << "7) Transverse momentum reconstructed =" << recpt << " at eta= " << receta
0522 << " and phi= " << recphi << std::endl;
0523 }
0524 } else {
0525 edm::LogVerbatim("TrackValidator")
0526 << "TrackingParticle #" << st << " with pt=" << sqrt(tp->momentum().perp2())
0527 << " NOT associated to any reco::Track"
0528 << "\n";
0529 receta = -100.;
0530 recphi = -100.;
0531 recpt = -100.;
0532 recd0 = -100.;
0533 recz0 = -100;
0534 respt = -100.;
0535 resd0 = -100.;
0536 resz0 = -100.;
0537 resphi = -100.;
0538 reseta = -100.;
0539 rescottheta = -100.;
0540 recnhit = 100;
0541 recchiq = -100;
0542 eff = 0;
0543 flagrec = 100;
0544 }
0545
0546 std::cout << "Eff=" << eff << std::endl;
0547
0548
0549
0550 std::cout << "Flag is" << flag << std::endl;
0551 std::cout << "RecFlag is" << flagrec << std::endl;
0552
0553 if (countpart[0] == 2 && flag == 0) {
0554 mzmu =
0555 sqrt((ene[0][0] + ene[0][1]) * (ene[0][0] + ene[0][1]) - (px[0][0] + px[0][1]) * (px[0][0] + px[0][1]) -
0556 (py[0][0] + py[0][1]) * (py[0][0] + py[0][1]) - (pz[0][0] + pz[0][1]) * (pz[0][0] + pz[0][1]));
0557 std::cout << "Mzmu " << mzmu << std::endl;
0558 ptzmu = sqrt((px[0][0] + px[0][1]) * (px[0][0] + px[0][1]) + (py[0][0] + py[0][1]) * (py[0][0] + py[0][1]));
0559
0560 pLzmu = pz[0][0] + pz[0][1];
0561 enezmu = ene[0][0] + ene[0][1];
0562 phizmu = atan2((py[0][0] + py[0][1]), (px[0][0] + px[0][1]));
0563 thetazmu = atan2(ptzmu, (pz[0][0] + pz[0][1]));
0564 etazmu = -log(tan(thetazmu * 3.14 / 360.));
0565 yzmu = 0.5 * log((enezmu + pLzmu) / (enezmu - pLzmu));
0566 mxptmu = std::max(ptmu[0][0], ptmu[0][1]);
0567 minptmu = std::min(ptmu[0][0], ptmu[0][1]);
0568 } else {
0569 mzmu = -100.;
0570 ptzmu = -100.;
0571 pLzmu = -100.;
0572 enezmu = -100.;
0573 etazmu = -100.;
0574 phizmu = -100.;
0575 thetazmu = -100.;
0576 yzmu = -100.;
0577 mxptmu = -100.;
0578 minptmu = -100.;
0579 }
0580
0581
0582 if (countpartrec[0] == 2 && flagrec == 0) {
0583 recmzmu = sqrt((recene[0][0] + recene[0][1]) * (recene[0][0] + recene[0][1]) -
0584 (recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) -
0585 (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]) -
0586 (recpz[0][0] + recpz[0][1]) * (recpz[0][0] + recpz[0][1]));
0587 std::cout << "RecMzmu " << recmzmu << std::endl;
0588 recptzmu = sqrt((recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) +
0589 (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]));
0590
0591 recpLzmu = recpz[0][0] + recpz[0][1];
0592 recenezmu = recene[0][0] + recene[0][1];
0593 recphizmu = atan2((recpy[0][0] + recpy[0][1]), (recpx[0][0] + recpx[0][1]));
0594 recthetazmu = atan2(recptzmu, (recpz[0][0] + recpz[0][1]));
0595 recetazmu = -log(tan(recthetazmu * 3.14 / 360.));
0596 recyzmu = 0.5 * log((recenezmu + recpLzmu) / (recenezmu - recpLzmu));
0597 recmxptmu = std::max(recptmu[0][0], recptmu[0][1]);
0598 recminptmu = std::min(recptmu[0][0], recptmu[0][1]);
0599 } else {
0600 recmzmu = -100.;
0601 recptzmu = -100.;
0602 recpLzmu = -100.;
0603 recenezmu = -100.;
0604 recetazmu = -100.;
0605 recphizmu = -100.;
0606 recthetazmu = -100.;
0607 recyzmu = -100.;
0608 recmxptmu = -100;
0609 recminptmu = -100.;
0610 }
0611
0612 tree_eff->Fill();
0613
0614 }
0615 }
0616 }
0617
0618
0619
0620
0621 if (selection_fake) {
0622 std::cout << "Computing Fake Rate" << std::endl;
0623
0624 fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0625 faked0 = 0., fakez0 = 0.;
0626 fakenhit = 0;
0627 fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0628 fakerecz0 = 0.;
0629 fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0630 fakerecchiq = 0.;
0631 fakerecnhit = 0;
0632 faketrackType = 0;
0633 fake = 0;
0634
0635
0636 int rT = 0;
0637 for (reco::TrackCollection::size_type i = 0; i < tC.size(); ++i) {
0638 edm::RefToBase<reco::Track> track(trackCollection, i);
0639 rT++;
0640
0641 fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0642 faked0 = 0., fakez0 = 0.;
0643 fakenhit = 0;
0644 fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0645 fakerecz0 = 0.;
0646 fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0647 fakerecchiq = 0.;
0648 fakerecnhit = 0;
0649 faketrackType = 0;
0650 fake = 0;
0651
0652 fakerecnhit = track->numberOfValidHits();
0653 fakerecchiq = track->normalizedChi2();
0654 fakerectheta = track->theta();
0655 fakereccottheta = 1. / tan(rectheta);
0656
0657 fakereceta = track->momentum().eta();
0658
0659 fakerecphi = track->phi();
0660 fakerecpt = track->pt();
0661 fakerecd0 = track->d0();
0662 fakerecz0 = track->dz();
0663
0664 std::cout << "1) Before assoc: TkRecTrack at eta = " << fakereceta << std::endl;
0665 std::cout << "Track number " << i << std::endl;
0666 std::cout << "\tPT: " << track->pt() << std::endl;
0667 std::cout << "\timpact parameter:d0: " << track->d0() << std::endl;
0668 std::cout << "\timpact parameter:z0: " << track->dz() << std::endl;
0669 std::cout << "\tAzimuthal angle of point of closest approach:" << track->phi() << std::endl;
0670 std::cout << "\tcharge: " << track->charge() << std::endl;
0671 std::cout << "\teta: " << track->eta() << std::endl;
0672 std::cout << "\tnormalizedChi2: " << track->normalizedChi2() << std::endl;
0673
0674 std::vector<std::pair<TrackingParticleRef, double> > tp;
0675
0676
0677 if (recSimColl.find(track) != recSimColl.end()) {
0678 tp = recSimColl[track];
0679 if (!tp.empty()) {
0680 edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0681 << " associated with quality:" << tp.begin()->second << "\n";
0682
0683 TrackingParticleRef tpr = tp.begin()->first;
0684 const SimTrack* fakeassocTrack = &(*tpr->g4Track_begin());
0685
0686 const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0687 FreeTrajectoryState ftsAtProduction(
0688 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0689 GlobalVector(
0690 fakeassocTrack->momentum().x(), fakeassocTrack->momentum().y(), fakeassocTrack->momentum().z()),
0691 TrackCharge(tpr->charge()),
0692 theMF);
0693 TSCPBuilderNoMaterial tscpBuilder;
0694 TrajectoryStateClosestToPoint tsAtClosestApproach =
0695 tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));
0696 GlobalPoint v = tsAtClosestApproach.theState().position();
0697 GlobalVector p = tsAtClosestApproach.theState().momentum();
0698
0699
0700
0701
0702 double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0703 double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0704 faked0 = float(-dxySim);
0705 fakez0 = float(dszSim * p.mag() / p.perp());
0706
0707 faketrackType = fakeassocTrack->type();
0708 faketheta = fakeassocTrack->momentum().theta();
0709 fakecottheta = 1. / tan(faketheta);
0710 fakeeta = fakeassocTrack->momentum().eta();
0711 fakephi = fakeassocTrack->momentum().phi();
0712 fakept = fakeassocTrack->momentum().pt();
0713 fakenhit = tpr->matchedHit();
0714
0715 std::cout << "4) After call to associator: the best SimTrack match is of type" << fakeassocTrack->type()
0716 << " ,at eta = " << fakeeta << " and phi = " << fakephi << " ,with pt at vertex = " << fakept
0717 << " GeV/c"
0718 << " ,d0 global = " << faked0 << " ,z0 = " << fakez0 << std::endl;
0719 fake = 1;
0720
0721 fakerespt = fakerecpt - fakept;
0722 fakeresd0 = fakerecd0 - faked0;
0723 fakeresz0 = fakerecz0 - fakez0;
0724 fakereseta = -log(tan(fakerectheta / 2.)) - (-log(tan(faketheta / 2.)));
0725 fakeresphi = fakerecphi - fakephi;
0726 fakerescottheta = fakereccottheta - fakecottheta;
0727 }
0728 } else {
0729 edm::LogVerbatim("TrackValidator")
0730 << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any TrackingParticle"
0731 << "\n";
0732
0733 fakeeta = -100.;
0734 faketheta = -100;
0735 fakephi = -100.;
0736 fakept = -100.;
0737 faked0 = -100.;
0738 fakez0 = -100;
0739 fakerespt = -100.;
0740 fakeresd0 = -100.;
0741 fakeresz0 = -100.;
0742 fakeresphi = -100.;
0743 fakereseta = -100.;
0744 fakerescottheta = -100.;
0745 fakenhit = 100;
0746 fake = 0;
0747 }
0748
0749 tree_fake->Fill();
0750 }
0751
0752 }
0753
0754 w++;
0755
0756 }
0757 }
0758
0759
0760 void ValidationMisalignedTracker::endJob() { std::cout << "\t Misalignment analysis completed \n" << std::endl; }
0761
0762 DEFINE_FWK_MODULE(ValidationMisalignedTracker);