00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "Validation/RecoVertex/interface/V0Validator.h"
00022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00023
00024 typedef std::vector<TrackingVertex> TrackingVertexCollection;
00025 typedef edm::Ref<TrackingVertexCollection> TrackingVertexRef;
00026 typedef edm::RefVector<edm::HepMCProduct, HepMC::GenVertex > GenVertexRefVector;
00027 typedef edm::RefVector<edm::HepMCProduct, HepMC::GenParticle > GenParticleRefVector;
00028
00029 const double piMass = 0.13957018;
00030 const double piMassSquared = piMass*piMass;
00031 const double protonMass = 0.93827203;
00032 const double protonMassSquared = protonMass*protonMass;
00033
00034
00035
00036 V0Validator::V0Validator(const edm::ParameterSet& iConfig) :
00037 theDQMRootFileName(iConfig.getParameter<std::string>("DQMRootFileName")),
00038 k0sCollectionTag(iConfig.getParameter<edm::InputTag>("kShortCollection")),
00039 lamCollectionTag(iConfig.getParameter<edm::InputTag>("lambdaCollection")),
00040 dirName(iConfig.getParameter<std::string>("dirName")) {
00041 genLam = genK0s = realLamFoundEff = realK0sFoundEff = lamCandFound =
00042 k0sCandFound = noTPforK0sCand = noTPforLamCand = realK0sFound = realLamFound = 0;
00043 theDQMstore = edm::Service<DQMStore>().operator->();
00044 }
00045
00046
00047 V0Validator::~V0Validator() {
00048
00049 }
00050
00051
00052
00053
00054
00055
00056 void V0Validator::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
00057
00058
00059
00060
00061
00062 theDQMstore->cd();
00063 std::string subDirName = dirName + "/EffFakes";
00064 theDQMstore->setCurrentFolder(subDirName.c_str());
00065
00066 ksEffVsR = theDQMstore->book1D("K0sEffVsR",
00067 "K^{0}_{S} Efficiency vs #rho", 40, 0., 40.);
00068 ksEffVsEta = theDQMstore->book1D("K0sEffVsEta",
00069 "K^{0}_{S} Efficiency vs #eta", 40, -2.5, 2.5);
00070 ksEffVsPt = theDQMstore->book1D("K0sEffVsPt",
00071 "K^{0}_{S} Efficiency vs p_{T}", 70, 0., 20.);;
00072
00073 ksTkEffVsR = theDQMstore->book1D("K0sTkEffVsR",
00074 "K^{0}_{S} Tracking Efficiency vs #rho", 40, 0., 40.);
00075 ksTkEffVsEta = theDQMstore->book1D("K0sTkEffVsEta",
00076 "K^{0}_{S} Tracking Efficiency vs #eta", 40, -2.5, 2.5);
00077 ksTkEffVsPt = theDQMstore->book1D("K0sTkEffVsPt",
00078 "K^{0}_{S} Tracking Efficiency vs p_{T}", 70, 0., 20.);
00079
00080 ksEffVsR_num = theDQMstore->book1D("K0sEffVsR_num",
00081 "K^{0}_{S} Efficiency vs #rho", 40, 0., 40.);
00082 ksEffVsEta_num = theDQMstore->book1D("K0sEffVsEta_num",
00083 "K^{0}_{S} Efficiency vs #eta", 40, -2.5, 2.5);
00084 ksEffVsPt_num = theDQMstore->book1D("K0sEffVsPt_num",
00085 "K^{0}_{S} Efficiency vs p_{T}", 70, 0., 20.);;
00086
00087 ksTkEffVsR_num = theDQMstore->book1D("K0sTkEffVsR_num",
00088 "K^{0}_{S} Tracking Efficiency vs #rho", 40, 0., 40.);
00089 ksTkEffVsEta_num = theDQMstore->book1D("K0sTkEffVsEta_num",
00090 "K^{0}_{S} Tracking Efficiency vs #eta", 40, -2.5, 2.5);
00091 ksTkEffVsPt_num = theDQMstore->book1D("K0sTkEffVsPt_num",
00092 "K^{0}_{S} Tracking Efficiency vs p_{T}", 70, 0., 20.);;
00093
00094
00095 ksEffVsR_denom = theDQMstore->book1D("K0sEffVsR_denom",
00096 "K^{0}_{S} Efficiency vs #rho", 40, 0., 40.);
00097 ksEffVsEta_denom = theDQMstore->book1D("K0sEffVsEta_denom",
00098 "K^{0}_{S} Efficiency vs #eta", 40, -2.5, 2.5);
00099 ksEffVsPt_denom = theDQMstore->book1D("K0sEffVsPt_denom",
00100 "K^{0}_{S} Efficiency vs p_{T}", 70, 0., 20.);;
00101
00102
00103 lamEffVsR = theDQMstore->book1D("LamEffVsR",
00104 "#Lambda^{0} Efficiency vs #rho", 40, 0., 40.);
00105 lamEffVsEta = theDQMstore->book1D("LamEffVsEta",
00106 "#Lambda^{0} Efficiency vs #eta", 40, -2.5, 2.5);
00107 lamEffVsPt = theDQMstore->book1D("LamEffVsPt",
00108 "#Lambda^{0} Efficiency vs p_{T}", 70, 0., 20.);
00109
00110
00111 lamTkEffVsR = theDQMstore->book1D("LamTkEffVsR",
00112 "#Lambda^{0} TrackingEfficiency vs #rho", 40, 0., 40.);
00113 lamTkEffVsEta = theDQMstore->book1D("LamTkEffVsEta",
00114 "#Lambda^{0} Tracking Efficiency vs #eta", 40, -2.5, 2.5);
00115 lamTkEffVsPt = theDQMstore->book1D("LamTkEffVsPt",
00116 "#Lambda^{0} Tracking Efficiency vs p_{T}", 70, 0., 20.);
00117
00118 lamEffVsR_num = theDQMstore->book1D("LamEffVsR_num",
00119 "#Lambda^{0} Efficiency vs #rho", 40, 0., 40.);
00120 lamEffVsEta_num = theDQMstore->book1D("LamEffVsEta_num",
00121 "#Lambda^{0} Efficiency vs #eta", 40, -2.5, 2.5);
00122 lamEffVsPt_num = theDQMstore->book1D("LamEffVsPt_num",
00123 "#Lambda^{0} Efficiency vs p_{T}", 70, 0., 20.);
00124
00125
00126 lamTkEffVsR_num = theDQMstore->book1D("LamTkEffVsR_num",
00127 "#Lambda^{0} TrackingEfficiency vs #rho", 40, 0., 40.);
00128 lamTkEffVsEta_num = theDQMstore->book1D("LamTkEffVsEta_num",
00129 "#Lambda^{0} Tracking Efficiency vs #eta", 40, -2.5, 2.5);
00130 lamTkEffVsPt_num = theDQMstore->book1D("LamTkEffVsPt_num",
00131 "#Lambda^{0} Tracking Efficiency vs p_{T}", 70, 0., 20.);
00132
00133
00134 lamEffVsR_denom = theDQMstore->book1D("LamEffVsR_denom",
00135 "#Lambda^{0} Efficiency vs #rho", 40, 0., 40.);
00136 lamEffVsEta_denom = theDQMstore->book1D("LamEffVsEta_denom",
00137 "#Lambda^{0} Efficiency vs #eta", 40, -2.5, 2.5);
00138 lamEffVsPt_denom = theDQMstore->book1D("LamEffVsPt_denom",
00139 "#Lambda^{0} Efficiency vs p_{T}", 70, 0., 20.);
00140
00141
00142
00143
00144
00145
00146 ksFakeVsR = theDQMstore->book1D("K0sFakeVsR",
00147 "K^{0}_{S} Fake Rate vs #rho", 40, 0., 40.);
00148 ksFakeVsEta = theDQMstore->book1D("K0sFakeVsEta",
00149 "K^{0}_{S} Fake Rate vs #eta", 40, -2.5, 2.5);
00150 ksFakeVsPt = theDQMstore->book1D("K0sFakeVsPt",
00151 "K^{0}_{S} Fake Rate vs p_{T}", 70, 0., 20.);
00152 ksTkFakeVsR = theDQMstore->book1D("K0sTkFakeVsR",
00153 "K^{0}_{S} Tracking Fake Rate vs #rho", 40, 0., 40.);
00154 ksTkFakeVsEta = theDQMstore->book1D("K0sTkFakeVsEta",
00155 "K^{0}_{S} Tracking Fake Rate vs #eta", 40, -2.5, 2.5);
00156 ksTkFakeVsPt = theDQMstore->book1D("K0sTkFakeVsPt",
00157 "K^{0}_{S} Tracking Fake Rate vs p_{T}", 70, 0., 20.);
00158
00159 ksFakeVsR_num = theDQMstore->book1D("K0sFakeVsR_num",
00160 "K^{0}_{S} Fake Rate vs #rho", 40, 0., 40.);
00161 ksFakeVsEta_num = theDQMstore->book1D("K0sFakeVsEta_num",
00162 "K^{0}_{S} Fake Rate vs #eta", 40, -2.5, 2.5);
00163 ksFakeVsPt_num = theDQMstore->book1D("K0sFakeVsPt_num",
00164 "K^{0}_{S} Fake Rate vs p_{T}", 70, 0., 20.);
00165 ksTkFakeVsR_num = theDQMstore->book1D("K0sTkFakeVsR_num",
00166 "K^{0}_{S} Tracking Fake Rate vs #rho", 40, 0., 40.);
00167 ksTkFakeVsEta_num = theDQMstore->book1D("K0sTkFakeVsEta_num",
00168 "K^{0}_{S} Tracking Fake Rate vs #eta", 40, -2.5, 2.5);
00169 ksTkFakeVsPt_num = theDQMstore->book1D("K0sTkFakeVsPt_num",
00170 "K^{0}_{S} Tracking Fake Rate vs p_{T}", 70, 0., 20.);
00171
00172 ksFakeVsR_denom = theDQMstore->book1D("K0sFakeVsR_denom",
00173 "K^{0}_{S} Fake Rate vs #rho", 40, 0., 40.);
00174 ksFakeVsEta_denom = theDQMstore->book1D("K0sFakeVsEta_denom",
00175 "K^{0}_{S} Fake Rate vs #eta", 40, -2.5, 2.5);
00176 ksFakeVsPt_denom = theDQMstore->book1D("K0sFakeVsPt_denom",
00177 "K^{0}_{S} Fake Rate vs p_{T}", 70, 0., 20.);
00178
00179 lamFakeVsR = theDQMstore->book1D("LamFakeVsR",
00180 "#Lambda^{0} Fake Rate vs #rho", 40, 0., 40.);
00181 lamFakeVsEta = theDQMstore->book1D("LamFakeVsEta",
00182 "#Lambda^{0} Fake Rate vs #eta", 40, -2.5, 2.5);
00183 lamFakeVsPt = theDQMstore->book1D("LamFakeVsPt",
00184 "#Lambda^{0} Fake Rate vs p_{T}", 70, 0., 20.);
00185 lamTkFakeVsR = theDQMstore->book1D("LamTkFakeVsR",
00186 "#Lambda^{0} Tracking Fake Rate vs #rho", 40, 0., 40.);
00187 lamTkFakeVsEta = theDQMstore->book1D("LamTkFakeVsEta",
00188 "#Lambda^{0} Tracking Fake Rate vs #eta", 40, -2.5, 2.5);
00189 lamTkFakeVsPt = theDQMstore->book1D("LamTkFakeVsPt",
00190 "#Lambda^{0} Tracking Fake Rate vs p_{T}", 70, 0., 20.);
00191
00192 lamFakeVsR_num = theDQMstore->book1D("LamFakeVsR_num",
00193 "#Lambda^{0} Fake Rate vs #rho", 40, 0., 40.);
00194 lamFakeVsEta_num = theDQMstore->book1D("LamFakeVsEta_num",
00195 "#Lambda^{0} Fake Rate vs #eta", 40, -2.5, 2.5);
00196 lamFakeVsPt_num = theDQMstore->book1D("LamFakeVsPt_num",
00197 "#Lambda^{0} Fake Rate vs p_{T}", 70, 0., 20.);
00198 lamTkFakeVsR_num = theDQMstore->book1D("LamTkFakeVsR_num",
00199 "#Lambda^{0} Tracking Fake Rate vs #rho", 40, 0., 40.);
00200 lamTkFakeVsEta_num = theDQMstore->book1D("LamTkFakeVsEta_num",
00201 "#Lambda^{0} Tracking Fake Rate vs #eta", 40, -2.5, 2.5);
00202 lamTkFakeVsPt_num = theDQMstore->book1D("LamTkFakeVsPt_num",
00203 "#Lambda^{0} Tracking Fake Rate vs p_{T}", 70, 0., 20.);
00204
00205 lamFakeVsR_denom = theDQMstore->book1D("LamFakeVsR_denom",
00206 "#Lambda^{0} Fake Rate vs #rho", 40, 0., 40.);
00207 lamFakeVsEta_denom = theDQMstore->book1D("LamFakeVsEta_denom",
00208 "#Lambda^{0} Fake Rate vs #eta", 40, -2.5, 2.5);
00209 lamFakeVsPt_denom = theDQMstore->book1D("LamFakeVsPt_denom",
00210 "#Lambda^{0} Fake Rate vs p_{T}", 70, 0., 20.);
00211
00212 theDQMstore->cd();
00213 subDirName = dirName + "/Other";
00214 theDQMstore->setCurrentFolder(subDirName.c_str());
00215
00216 nKs = theDQMstore->book1D("nK0s",
00217 "Number of K^{0}_{S} found per event", 60, 0., 60.);
00218 nLam = theDQMstore->book1D("nLam",
00219 "Number of #Lambda^{0} found per event", 60, 0., 60.);
00220
00221 ksXResolution = theDQMstore->book1D("ksXResolution",
00222 "Resolution of V0 decay vertex X coordinate", 50, 0., 50.);
00223 ksYResolution = theDQMstore->book1D("ksYResolution",
00224 "Resolution of V0 decay vertex Y coordinate", 50, 0., 50.);
00225 ksZResolution = theDQMstore->book1D("ksZResolution",
00226 "Resolution of V0 decay vertex Z coordinate", 50, 0., 50.);
00227 lamXResolution = theDQMstore->book1D("lamXResolution",
00228 "Resolution of V0 decay vertex X coordinate", 50, 0., 50.);
00229 lamYResolution = theDQMstore->book1D("lamYResolution",
00230 "Resolution of V0 decay vertex Y coordinate", 50, 0., 50.);
00231 lamZResolution = theDQMstore->book1D("lamZResolution",
00232 "Resolution of V0 decay vertex Z coordinate", 50, 0., 50.);
00233 ksAbsoluteDistResolution = theDQMstore->book1D("ksRResolution",
00234 "Resolution of absolute distance from primary vertex to V0 vertex",
00235 100, 0., 50.);
00236 lamAbsoluteDistResolution = theDQMstore->book1D("lamRResolution",
00237 "Resolution of absolute distance from primary vertex to V0 vertex",
00238 100, 0., 50.);
00239
00240 ksCandStatus = theDQMstore->book1D("ksCandStatus",
00241 "Fake type by cand status",
00242 10, 0., 10.);
00243 lamCandStatus = theDQMstore->book1D("ksCandStatus",
00244 "Fake type by cand status",
00245 10, 0., 10.);
00246
00247 double minKsMass = 0.49767 - 0.07;
00248 double maxKsMass = 0.49767 + 0.07;
00249 double minLamMass = 1.1156 - 0.05;
00250 double maxLamMass = 1.1156 + 0.05;
00251 int ksMassNbins = 100;
00252 double ksMassXmin = minKsMass;
00253 double ksMassXmax = maxKsMass;
00254 int lamMassNbins = 100;
00255 double lamMassXmin = minLamMass;
00256 double lamMassXmax = maxLamMass;
00257
00258 fakeKsMass = theDQMstore->book1D("ksMassFake",
00259 "Mass of fake K0S",
00260 ksMassNbins, minKsMass, maxKsMass);
00261 goodKsMass = theDQMstore->book1D("ksMassGood",
00262 "Mass of good reco K0S",
00263 ksMassNbins, minKsMass, maxKsMass);
00264 fakeLamMass = theDQMstore->book1D("lamMassFake",
00265 "Mass of fake Lambda",
00266 lamMassNbins, minLamMass, maxLamMass);
00267 goodLamMass = theDQMstore->book1D("lamMassGood",
00268 "Mass of good Lambda",
00269 lamMassNbins, minLamMass, maxLamMass);
00270
00271 ksMassAll = theDQMstore->book1D("ksMassAll",
00272 "Invariant mass of all K0S",
00273 ksMassNbins, ksMassXmin, ksMassXmax);
00274 lamMassAll = theDQMstore->book1D("lamMassAll",
00275 "Invariant mass of all #Lambda^{0}",
00276 lamMassNbins, lamMassXmin, lamMassXmax);
00277
00278 ksFakeDauRadDist = theDQMstore->book1D("radDistFakeKs",
00279 "Production radius of daughter particle of Ks fake",
00280 100, 0., 15.);
00281 lamFakeDauRadDist = theDQMstore->book1D("radDistFakeLam",
00282 "Production radius of daughter particle of Lam fake",
00283 100, 0., 15.);
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 }
00450
00451 void V0Validator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00452
00453 using std::cout;
00454 using std::endl;
00455 using namespace edm;
00456 using namespace std;
00457
00458
00459
00460 ESHandle<MagneticField> bFieldHandle;
00461 iSetup.get<IdealMagneticFieldRecord>().get(bFieldHandle);
00462 ESHandle<GlobalTrackingGeometry> globTkGeomHandle;
00463 iSetup.get<GlobalTrackingGeometryRecord>().get(globTkGeomHandle);
00464
00465
00466
00467
00468
00469 Handle<reco::RecoToSimCollection > recotosimCollectionH;
00470 iEvent.getByLabel("trackingParticleRecoTrackAsssociation", recotosimCollectionH);
00471
00472
00473 Handle<reco::SimToRecoCollection> simtorecoCollectionH;
00474 iEvent.getByLabel("trackingParticleRecoTrackAsssociation", simtorecoCollectionH);
00475
00476
00477 edm::Handle<TrackingParticleCollection> TPCollectionEff ;
00478 iEvent.getByLabel("mergedtruth", "MergedTrackTruth", TPCollectionEff);
00479 const TrackingParticleCollection tPCeff = *( TPCollectionEff.product() );
00480
00481 edm::ESHandle<TrackAssociatorBase> associatorByHits;
00482 iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits", associatorByHits);
00483
00484
00485
00486
00487
00488
00489
00490
00491 Handle< View<reco::Track> > trackCollectionH;
00492 iEvent.getByLabel("generalTracks", trackCollectionH);
00493
00494 Handle<SimTrackContainer> simTrackCollection;
00495 iEvent.getByLabel("g4SimHits", simTrackCollection);
00496 const SimTrackContainer simTC = *(simTrackCollection.product());
00497
00498 Handle<SimVertexContainer> simVertexCollection;
00499 iEvent.getByLabel("g4SimHits", simVertexCollection);
00500 const SimVertexContainer simVC = *(simVertexCollection.product());
00501
00502
00503
00504 edm::Handle<TrackingParticleCollection> TPCollectionH ;
00505 iEvent.getByLabel("mergedtruth", "MergedTrackTruth", TPCollectionH);
00506 const View<reco::Track> tC = *( trackCollectionH.product() );
00507
00508
00509
00510
00511
00512
00513 edm::Handle< std::vector<reco::Vertex> > primaryVtxCollectionH;
00514 iEvent.getByLabel("offlinePrimaryVertices", primaryVtxCollectionH);
00515 const reco::VertexCollection primaryVertexCollection = *(primaryVtxCollectionH.product());
00516
00517 reco::Vertex* thePrimary = 0;
00518 std::vector<reco::Vertex>::const_iterator iVtxPH = primaryVtxCollectionH->begin();
00519 for(std::vector<reco::Vertex>::const_iterator iVtx = primaryVtxCollectionH->begin();
00520 iVtx < primaryVtxCollectionH->end();
00521 iVtx++) {
00522 if(primaryVtxCollectionH->size() > 1) {
00523 if(iVtx->tracksSize() > iVtxPH->tracksSize()) {
00524 iVtxPH = iVtx;
00525 }
00526 }
00527 else iVtxPH = iVtx;
00528 }
00529 thePrimary = new reco::Vertex(*iVtxPH);
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 edm::Handle<reco::VertexCompositeCandidateCollection> k0sCollection;
00542 edm::Handle<reco::VertexCompositeCandidateCollection> lambdaCollection;
00543
00544
00545 iEvent.getByLabel(k0sCollectionTag, k0sCollection);
00546 iEvent.getByLabel(lamCollectionTag, lambdaCollection);
00547
00548
00549 std::vector< pair<TrackingParticleRef, TrackingParticleRef> > trueK0s;
00550 std::vector< pair<TrackingParticleRef, TrackingParticleRef> > trueLams;
00551 std::vector<double> trueKsMasses;
00552 std::vector<double> trueLamMasses;
00553
00555
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00580
00582
00583
00584
00585 double numK0sFound = 0.;
00586 double mass = 0.;
00587 std::vector<double> radDist;
00588
00589
00590 if ( k0sCollection->size() > 0 ) {
00591
00592
00593 vector<reco::TrackRef> theDaughterTracks;
00594 for( reco::VertexCompositeCandidateCollection::const_iterator iK0s = k0sCollection->begin();
00595 iK0s != k0sCollection->end();
00596 iK0s++) {
00597
00598
00599 ksMassAll->Fill( iK0s->mass() );
00600
00601 K0sCandpT = (sqrt( iK0s->momentum().perp2() ));
00602 K0sCandEta = iK0s->momentum().eta();
00603 K0sCandR = (sqrt( iK0s->vertex().perp2() ));
00604 K0sCandStatus = 0;
00605
00606 mass = iK0s->mass();
00607
00608
00609 theDaughterTracks.push_back( (*(dynamic_cast<const reco::RecoChargedCandidate *> (iK0s->daughter(0)) )).track() );
00610 theDaughterTracks.push_back( (*(dynamic_cast<const reco::RecoChargedCandidate *> (iK0s->daughter(1)) )).track() );
00611
00612
00613 for (int itrack = 0; itrack < 2; itrack++) {
00614 K0sPiCandStatus[itrack] = 0;
00615 }
00616
00617 std::vector< std::pair<TrackingParticleRef, double> > tp;
00618 TrackingParticleRef tpref;
00619 TrackingParticleRef firstDauTP;
00620 TrackingVertexRef k0sVtx;
00621
00622
00623
00624 for(View<reco::Track>::size_type i=0; i<theDaughterTracks.size(); ++i){
00625
00626 RefToBase<reco::Track> track( theDaughterTracks.at(i) );
00627
00628
00629 if(recotosimCollectionH->find(track) != recotosimCollectionH->end()) {
00630
00631 tp = (*recotosimCollectionH)[track];
00632 if (tp.size() != 0) {
00633 K0sPiCandStatus[i] = 1;
00634 tpref = tp.begin()->first;
00635
00636
00637 if( simtorecoCollectionH->find(tpref) == simtorecoCollectionH->end() ) {
00638 K0sPiCandStatus[i] = 3;
00639 }
00640
00641 TrackingVertexRef parentVertex = tpref->parentVertex();
00642 if(parentVertex.isNonnull()) radDist.push_back(parentVertex->position().R());
00643
00644 if( parentVertex.isNonnull() ) {
00645 if( k0sVtx.isNonnull() ) {
00646 if( k0sVtx->position() == parentVertex->position() ) {
00647 if( parentVertex->nDaughterTracks() == 2 ) {
00648 if( parentVertex->nSourceTracks() == 0 ) {
00649
00650 K0sCandStatus = 6;
00651 }
00652
00653 for( TrackingVertex::tp_iterator iTP = parentVertex->sourceTracks_begin();
00654 iTP != parentVertex->sourceTracks_end(); iTP++) {
00655 if( (*iTP)->pdgId() == 310 ) {
00656
00657 K0sCandStatus = 1;
00658 realK0sFound++;
00659 numK0sFound += 1.;
00660 std::pair<TrackingParticleRef, TrackingParticleRef> pair(firstDauTP, tpref);
00661
00662 trueK0s.push_back(pair);
00663 trueKsMasses.push_back(mass);
00664 }
00665 else {
00666 K0sCandStatus = 2;
00667 if( (*iTP)->pdgId() == 3122 ) {
00668 K0sCandStatus = 7;
00669 }
00670 }
00671 }
00672 }
00673 else {
00674
00675 K0sCandStatus = 3;
00676 }
00677 }
00678 else {
00679
00680 K0sCandStatus = 4;
00681 }
00682 }
00683 else {
00684
00685 k0sVtx = parentVertex;
00686 firstDauTP = tpref;
00687 }
00688 }
00689 }
00690 }
00691 else {
00692
00693 K0sPiCandStatus[i] = 2;
00694 noTPforK0sCand++;
00695 K0sCandStatus = 5;
00696 theDaughterTracks.clear();
00697 }
00698 }
00699
00700 theDaughterTracks.clear();
00701
00702 if( K0sCandStatus > 1 ) {
00703
00704 ksFakeVsR_num->Fill(K0sCandR);
00705 ksFakeVsEta_num->Fill(K0sCandEta);
00706 ksFakeVsPt_num->Fill(K0sCandpT);
00707 ksCandStatus->Fill((float) K0sCandStatus);
00708 fakeKsMass->Fill(mass);
00709 for( unsigned int ndx = 0; ndx < radDist.size(); ndx++ ) {
00710 ksFakeDauRadDist->Fill(radDist[ndx]);
00711 }
00712 }
00713 if( K0sCandStatus == 5 ) {
00714 ksTkFakeVsR_num->Fill(K0sCandR);
00715 ksTkFakeVsEta_num->Fill(K0sCandEta);
00716 ksTkFakeVsPt_num->Fill(K0sCandpT);
00717 }
00718 ksFakeVsR_denom->Fill(K0sCandR);
00719 ksFakeVsEta_denom->Fill(K0sCandEta);
00720 ksFakeVsPt_denom->Fill(K0sCandpT);
00721 }
00722 }
00723
00724
00725
00726 nKs->Fill( (float) numK0sFound );
00727 numK0sFound = 0.;
00728
00729
00730
00731 double numLamFound = 0.;
00732 mass = 0.;
00733 radDist.clear();
00734
00735 if ( lambdaCollection->size() > 0 ) {
00736
00737
00738 vector<reco::TrackRef> theDaughterTracks;
00739 for( reco::VertexCompositeCandidateCollection::const_iterator iLam = lambdaCollection->begin();
00740 iLam != lambdaCollection->end();
00741 iLam++) {
00742
00743 lamMassAll->Fill( iLam->mass() );
00744
00745 LamCandpT = (sqrt( iLam->momentum().perp2() ));
00746 LamCandEta = iLam->momentum().eta();
00747 LamCandR = (sqrt( iLam->vertex().perp2() ));
00748 LamCandStatus = 0;
00749 mass = iLam->mass();
00750
00751
00752 theDaughterTracks.push_back( (*(dynamic_cast<const reco::RecoChargedCandidate *> (iLam->daughter(0)) )).track() );
00753 theDaughterTracks.push_back( (*(dynamic_cast<const reco::RecoChargedCandidate *> (iLam->daughter(1)) )).track() );
00754
00755 for (int itrack = 0; itrack < 2; itrack++) {
00756 LamPiCandStatus[itrack] = 0;
00757 }
00758
00759 std::vector< std::pair<TrackingParticleRef, double> > tp;
00760 TrackingParticleRef tpref;
00761 TrackingParticleRef firstDauTP;
00762 TrackingVertexRef LamVtx;
00763
00764 for(View<reco::Track>::size_type i=0; i<theDaughterTracks.size(); ++i){
00765
00766
00767 RefToBase<reco::Track> track( theDaughterTracks.at(i) );
00768
00769
00770 if(recotosimCollectionH->find(track) != recotosimCollectionH->end()) {
00771
00772 tp = (*recotosimCollectionH)[track];
00773 if (tp.size() != 0) {
00774 LamPiCandStatus[i] = 1;
00775 tpref = tp.begin()->first;
00776
00777
00778 if( simtorecoCollectionH->find(tpref) == simtorecoCollectionH->end() ) {
00779 LamPiCandStatus[i] = 3;
00780 }
00781 TrackingVertexRef parentVertex = tpref->parentVertex();
00782 if( parentVertex.isNonnull() ) radDist.push_back(parentVertex->position().R());
00783
00784 if( parentVertex.isNonnull() ) {
00785 if( LamVtx.isNonnull() ) {
00786 if( LamVtx->position() == parentVertex->position() ) {
00787 if( parentVertex->nDaughterTracks() == 2 ) {
00788 if( parentVertex->nSourceTracks() == 0 ) {
00789
00790 LamCandStatus = 6;
00791 }
00792
00793 for( TrackingVertex::tp_iterator iTP = parentVertex->sourceTracks_begin();
00794 iTP != parentVertex->sourceTracks_end(); ++iTP) {
00795 if( abs((*iTP)->pdgId()) == 3122 ) {
00796 LamCandStatus = 1;
00797 realLamFound++;
00798 numLamFound += 1.;
00799 std::pair<TrackingParticleRef, TrackingParticleRef> pair(firstDauTP, tpref);
00800
00801 trueLams.push_back(pair);
00802 trueLamMasses.push_back(mass);
00803 }
00804 else {
00805 LamCandStatus = 2;
00806 if( abs((*iTP)->pdgId() ) == 310 ) {
00807 LamCandStatus = 7;
00808 }
00809 }
00810
00811
00812
00813 }
00814 }
00815 else {
00816
00817 LamCandStatus = 3;
00818 }
00819 }
00820 else {
00821
00822 LamCandStatus = 4;
00823 }
00824 }
00825 else {
00826
00827 LamVtx = parentVertex;
00828 firstDauTP = tpref;
00829 }
00830 }
00831 }
00832 }
00833 else {
00834 LamPiCandStatus[i] = 2;
00835 noTPforLamCand++;
00836 LamCandStatus = 5;
00837 theDaughterTracks.clear();
00838 }
00839 }
00840 theDaughterTracks.clear();
00841
00842
00843 if( LamCandStatus > 1 ) {
00844
00845
00846 lamFakeVsR_num->Fill(LamCandR);
00847
00848 lamFakeVsEta_num->Fill(LamCandEta);
00849
00850 lamFakeVsPt_num->Fill(LamCandpT);
00851
00852 lamCandStatus->Fill((float) LamCandStatus);
00853
00854 fakeLamMass->Fill(mass);
00855
00856 for( unsigned int ndx = 0; ndx < radDist.size(); ndx++ ) {
00857 lamFakeDauRadDist->Fill(radDist[ndx]);
00858 }
00859 }
00860
00861 if( K0sCandStatus == 5 ) {
00862 lamTkFakeVsR_num->Fill(LamCandR);
00863 lamTkFakeVsEta_num->Fill(LamCandEta);
00864 lamTkFakeVsPt_num->Fill(LamCandpT);
00865 }
00866
00867 lamFakeVsR_denom->Fill(LamCandR);
00868 lamFakeVsEta_denom->Fill(LamCandEta);
00869 lamFakeVsPt_denom->Fill(LamCandpT);
00870 }
00871 }
00872
00873 nLam->Fill( (double) numLamFound );
00874 numLamFound = 0.;
00875
00876
00878
00880
00881
00882
00883
00884 for(TrackingParticleCollection::size_type i = 0; i < tPCeff.size(); i++) {
00885 TrackingParticleRef tpr1(TPCollectionEff, i);
00886 TrackingParticle* itp1 = const_cast<TrackingParticle*>(tpr1.get());
00887 if( (itp1->pdgId() == 211 || itp1->pdgId() == 2212)
00888 && itp1->parentVertex().isNonnull()
00889 && abs(itp1->momentum().eta()) < 2.4
00890 && sqrt( itp1->momentum().perp2() ) > 0.9) {
00891 bool isLambda = false;
00892 if( itp1->pdgId() == 2212 ) isLambda = true;
00893 if( itp1->parentVertex()->nDaughterTracks() == 2 ) {
00894
00895 TrackingVertexRef piCand1Vertex = itp1->parentVertex();
00896 for(TrackingVertex::tp_iterator iTP1 = piCand1Vertex->sourceTracks_begin();
00897 iTP1 != piCand1Vertex->sourceTracks_end(); iTP1++) {
00898 if( abs((*iTP1)->pdgId()) == 3122 ) {
00899
00900
00901 for(TrackingParticleCollection::size_type j=0;
00902 j < tPCeff.size();
00903 j++) {
00904 TrackingParticleRef tpr2(TPCollectionEff, j);
00905 TrackingParticle* itp2 = const_cast<TrackingParticle*>(tpr2.get());
00906 int particle2pdgId;
00907 if (isLambda) particle2pdgId = -211;
00908 else particle2pdgId = -2212;
00909 if( itp2->pdgId() == particle2pdgId
00910 && itp2->parentVertex().isNonnull()
00911 && abs(itp2->momentum().eta()) < 2.4
00912 && sqrt(itp2->momentum().perp2()) > 0.9) {
00913 if(itp2->parentVertex() == itp1->parentVertex()) {
00914
00915 TrackingVertexRef piCand2Vertex = itp2->parentVertex();
00916 for (TrackingVertex::tp_iterator iTP2 = piCand2Vertex->sourceTracks_begin();
00917 iTP2 != piCand2Vertex->sourceTracks_end();
00918 ++iTP2) {
00919 LamGenEta = LamGenpT = LamGenR = 0.;
00920 LamGenStatus = 0;
00921 for(int ifill = 0;
00922 ifill < 2;
00923 ifill++) {
00924
00925 }
00926 if( abs((*iTP2)->pdgId()) == 3122 ) {
00927
00928
00929 LamGenpT = sqrt((*iTP2)->momentum().perp2());
00930 LamGenEta = (*iTP2)->momentum().eta();
00931 LamGenR = sqrt(itp2->vertex().perp2());
00932 genLam++;
00933 if(trueLams.size() > 0) {
00934 int loop_1 = 0;
00935 for(std::vector< pair<TrackingParticleRef, TrackingParticleRef> >::const_iterator iEffCheck = trueLams.begin();
00936 iEffCheck != trueLams.end();
00937 iEffCheck++) {
00938
00939 if( itp1->parentVertex() == iEffCheck->first->parentVertex()
00940 && itp2->parentVertex() == iEffCheck->second->parentVertex() ) {
00941 realLamFoundEff++;
00942
00943 LamGenStatus = 1;
00944
00945 goodLamMass->Fill(trueLamMasses[loop_1]);
00946
00947 break;
00948 }
00949 else {
00950
00951 LamGenStatus = 2;
00952 }
00953 loop_1++;
00954 }
00955 }
00956 else {
00957
00958 LamGenStatus = 2;
00959 }
00960 std::vector< std::pair<RefToBase<reco::Track>, double> > rt1;
00961 std::vector< std::pair<RefToBase<reco::Track>, double> > rt2;
00962
00963
00964 if( simtorecoCollectionH->find(tpr1) != simtorecoCollectionH->end() ) {
00965
00966 rt1 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) (*simtorecoCollectionH)[tpr1];
00967 if(rt1.size() != 0) {
00968 LamPiEff[0] = 1;
00969 edm::RefToBase<reco::Track> t1 = rt1.begin()->first;
00970 }
00971 }
00972 else {
00973 LamPiEff[0] = 2;
00974 }
00975
00976 if( (simtorecoCollectionH->find(tpr2) != simtorecoCollectionH->end()) ) {
00977
00978 rt2 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) (*simtorecoCollectionH)[tpr2];
00979 if(rt2.size() != 0) {
00980 LamPiEff[1] = 1;
00981 edm::RefToBase<reco::Track> t2 = rt2.begin()->first;
00982 }
00983 }
00984 else {
00985 LamPiEff[1] = 2;
00986 }
00987
00988 if( LamGenStatus == 1
00989 && (LamPiEff[0] == 2 || LamPiEff[1] == 2) ) {
00990
00991 LamGenStatus = 4;
00992 realLamFoundEff--;
00993 }
00994 if( LamGenStatus == 2
00995 && (LamPiEff[0] == 2 || LamPiEff[1] == 2) ) {
00996
00997 LamGenStatus = 3;
00998 }
00999
01000
01001 if(LamGenR > 0.) {
01002 if(LamGenStatus == 1) {
01003 lamEffVsR_num->Fill(LamGenR);
01004 }
01005 if((double) LamGenStatus < 2.5) {
01006 lamTkEffVsR_num->Fill(LamGenR);
01007 }
01008 lamEffVsR_denom->Fill(LamGenR);
01009 }
01010 if(abs(LamGenEta) > 0.) {
01011 if(LamGenStatus == 1) {
01012 lamEffVsEta_num->Fill(LamGenEta);
01013 }
01014 if((double) LamGenStatus < 2.5) {
01015 lamTkEffVsEta_num->Fill(LamGenEta);
01016 }
01017 lamEffVsEta_denom->Fill(LamGenEta);
01018 }
01019 if(LamGenpT > 0.) {
01020 if(LamGenStatus == 1) {
01021 lamEffVsPt_num->Fill(LamGenpT);
01022 }
01023 if((double) LamGenStatus < 2.5) {
01024 lamTkEffVsPt_num->Fill(LamGenpT);
01025 }
01026 lamEffVsPt_denom->Fill(LamGenpT);
01027 }
01028 }
01029 }
01030 }
01031 }
01032 }
01033 }
01034 }
01035 }
01036 }
01037 }
01038
01039
01040
01041
01042 for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
01043 TrackingParticleRef tpr1(TPCollectionEff, i);
01044 TrackingParticle* itp1=const_cast<TrackingParticle*>(tpr1.get());
01045
01046 if ( itp1->pdgId() == 211
01047 && itp1->parentVertex().isNonnull()
01048 && abs(itp1->momentum().eta()) < 2.4
01049 && sqrt(itp1->momentum().perp2()) > 0.9) {
01050 if ( itp1->parentVertex()->nDaughterTracks() == 2 ) {
01051 TrackingVertexRef piCand1Vertex = itp1->parentVertex();
01052
01053 for (TrackingVertex::tp_iterator iTP1 = piCand1Vertex->sourceTracks_begin();
01054 iTP1 != piCand1Vertex->sourceTracks_end(); ++iTP1) {
01055
01056 if ( (*iTP1)->pdgId()==310 ) {
01057
01058 for (TrackingParticleCollection::size_type j=0; j<tPCeff.size(); j++){
01059 TrackingParticleRef tpr2(TPCollectionEff, j);
01060 TrackingParticle* itp2=const_cast<TrackingParticle*>(tpr2.get());
01061
01062 if ( itp2->pdgId() == -211 && itp2->parentVertex().isNonnull()
01063 && abs(itp2->momentum().eta()) < 2.4
01064 && sqrt(itp2->momentum().perp2()) > 0.9) {
01065
01066 if ( itp2->parentVertex() == itp1->parentVertex() ) {
01067 TrackingVertexRef piCand2Vertex = itp2->parentVertex();
01068 for (TrackingVertex::tp_iterator iTP2 = piCand2Vertex->sourceTracks_begin();
01069 iTP2 != piCand2Vertex->sourceTracks_end(); ++iTP2) {
01070
01071 K0sGenEta = K0sGenpT = K0sGenR = 0.;
01072 K0sGenStatus = 0;
01073 if( (*iTP2)->pdgId() == 310 ) {
01074 K0sGenpT = sqrt( (*iTP2)->momentum().perp2() );
01075 K0sGenEta = (*iTP2)->momentum().eta();
01076 K0sGenR = sqrt(itp2->vertex().perp2());
01077 genK0s++;
01078 int loop_2 = 0;
01079 if( trueK0s.size() > 0 ) {
01080 for( std::vector< pair<TrackingParticleRef, TrackingParticleRef> >::const_iterator iEffCheck = trueK0s.begin();
01081 iEffCheck != trueK0s.end();
01082 iEffCheck++) {
01083
01084 if (itp1->parentVertex()==iEffCheck->first->parentVertex() &&
01085 itp2->parentVertex()==iEffCheck->second->parentVertex()) {
01086 realK0sFoundEff++;
01087 K0sGenStatus = 1;
01088
01089 goodKsMass->Fill(trueKsMasses[loop_2]);
01090
01091 break;
01092 }
01093 else {
01094 K0sGenStatus = 2;
01095 }
01096 }
01097 }
01098 else {
01099 K0sGenStatus = 2;
01100 }
01101
01102
01103
01104
01105 std::vector<std::pair<RefToBase<reco::Track>, double> > rt1;
01106 std::vector<std::pair<RefToBase<reco::Track>, double> > rt2;
01107
01108
01109 if( simtorecoCollectionH->find(tpr1) != simtorecoCollectionH->end() ) {
01110 rt1 = (std::vector< std::pair<RefToBase<reco::Track>, double> >) (*simtorecoCollectionH)[tpr1];
01111
01112 if(rt1.size() != 0) {
01113
01114 K0sPiEff[0] = 1;
01115 edm::RefToBase<reco::Track> t1 = rt1.begin()->first;
01116 }
01117 }
01118 else {
01119
01120 K0sPiEff[0] = 2;
01121 }
01122
01123
01124 if( simtorecoCollectionH->find(tpr2) != simtorecoCollectionH->end() ) {
01125 rt2 = (std::vector< std::pair<RefToBase<reco::Track>, double> >) (*simtorecoCollectionH)[tpr2];
01126
01127 if(rt2.size() != 0) {
01128
01129 K0sPiEff[1] = 1;
01130 edm::RefToBase<reco::Track> t2 = rt2.begin()->first;
01131 }
01132 }
01133 else {
01134 K0sPiEff[1] = 2;
01135 }
01136
01137 if(K0sGenStatus == 1
01138 && (K0sPiEff[0] == 2 || K0sPiEff[1] == 2)) {
01139 K0sGenStatus = 4;
01140 realK0sFoundEff--;
01141 }
01142 if(K0sGenStatus == 2
01143 && (K0sPiEff[0] == 2 || K0sPiEff[1] == 2)) {
01144 K0sGenStatus = 3;
01145 }
01146 if(K0sPiEff[0] == 1 && K0sPiEff[1] == 1) {
01147 k0sTracksFound++;
01148 }
01149
01150 if(K0sGenR > 0.) {
01151 if(K0sGenStatus == 1) {
01152 ksEffVsR_num->Fill(K0sGenR);
01153 }
01154 if((double) K0sGenStatus < 2.5) {
01155 ksTkEffVsR_num->Fill(K0sGenR);
01156 }
01157 ksEffVsR_denom->Fill(K0sGenR);
01158 }
01159 if(abs(K0sGenEta) > 0.) {
01160 if(K0sGenStatus == 1) {
01161 ksEffVsEta_num->Fill(K0sGenEta);
01162 }
01163 if((double) K0sGenStatus < 2.5) {
01164 ksTkEffVsEta_num->Fill(K0sGenEta);
01165 }
01166 ksEffVsEta_denom->Fill(K0sGenEta);
01167 }
01168 if(K0sGenpT > 0.) {
01169 if(K0sGenStatus == 1) {
01170 ksEffVsPt_num->Fill(K0sGenpT);
01171 }
01172 if((double) K0sGenStatus < 2.5) {
01173 ksTkEffVsPt_num->Fill(K0sGenpT);
01174 }
01175 ksEffVsPt_denom->Fill(K0sGenpT);
01176 }
01177 }
01178 }
01179 }
01180 }
01181 }
01182 }
01183 }
01184 }
01185 }
01186 }
01187
01188 delete thePrimary;
01189 }
01190
01191 void V0Validator::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
01192
01193 if(theDQMRootFileName.size() && theDQMstore) {
01194 theDQMstore->save(theDQMRootFileName);
01195 }
01196 }
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326