Implements edm::EDAnalyzer.
Definition at line 148 of file ElectronPixelSeedAnalyzer.cc.
References alongMomentum, deltaR(), DetId::det(), lat::endl(), relval_parameters_module::energy, eta, PV3DBase< T, PVType, FrameType >::eta(), funct::exp(), edm::Event::getByLabel(), edm::Event::getByType(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), histeMC_, histeoverp_, histeta_, histetaclu_, histetaMC_, histetclu_, histnbclus_, histnbseeds_, histnrseeds_, histpt_, histptMC_, histq_, edm::Event::id(), inputCollection_, int, TrajectoryStateOnSurface::isValid(), it, TrajectoryStateOnSurface::localPosition(), LogDebug, lp, PV3DBase< T, PVType, FrameType >::mag(), mcEnergy, mcEta, mcPhi, mcPt, mcQ, FreeTrajectoryState::momentum(), oppositeToMomentum, p, pDD, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, pi, FreeTrajectoryState::position(), funct::pow(), edm::Handle< T >::product(), PropagatorWithMaterial::propagate(), r, range, DetId::rawId(), seedDphi1, seedDphi2, seedDrz1, seedDrz2, seedEta, seedLayer1, seedLayer2, seedMomentum, seedPhi, seedPhi1, seedPhi2, seedPt, seedQ, seedRz1, seedRz2, seedSide1, seedSide2, seedSubdet1, seedSubdet2, funct::sin(), funct::sqrt(), DetId::subdetId(), superclusterEnergy, superclusterEt, superclusterEta, superclusterPhi, GeomDet::surface(), t, theMagField, Surface::toGlobal(), GeomDet::toGlobal(), transformer_, TrajectoryStateTransform::transientState(), tree_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
00149 {
00150
00151
00152 typedef edm::OwnVector<TrackingRecHit> recHitContainer;
00153 typedef recHitContainer::const_iterator const_iterator;
00154 typedef std::pair<const_iterator,const_iterator> range;
00155
00156
00157
00158 edm::Handle<reco::BeamSpot> theBeamSpot;
00159 e.getByType(theBeamSpot);
00160
00161
00162
00163 edm::Handle<ElectronPixelSeedCollection> elSeeds;
00164 e.getByLabel(inputCollection_,elSeeds);
00165 edm::LogInfo("")<<"\n\n =================> Treating event "<<e.id()<<" Number of seeds "<<elSeeds.product()->size();
00166 int is=0;
00167
00168 FTSFromVertexToPointFactory myFTS;
00169 float mass=.000511;
00170 PropagatorWithMaterial* prop1stLayer = new PropagatorWithMaterial(oppositeToMomentum,mass,&(*theMagField));
00171 PropagatorWithMaterial* prop2ndLayer = new PropagatorWithMaterial(alongMomentum,mass,&(*theMagField));
00172
00173 float dphi1=0., dphi2=0., drz1=0., drz2=0.;
00174 float phi1=0., phi2=0., rz1=0., rz2=0.;
00175
00176 for( ElectronPixelSeedCollection::const_iterator MyS= (*elSeeds).begin(); MyS != (*elSeeds).end(); ++MyS) {
00177
00178 LogDebug("") <<"\nSeed nr "<<is<<": ";
00179 range r=(*MyS).recHits();
00180 LogDebug("")<<" Number of RecHits= "<<(*MyS).nHits();
00181 const GeomDet *det1=0;const GeomDet *det2=0;
00182
00183 TrajectorySeed::const_iterator it=r.first;
00184 DetId id1 = (*it).geographicalId();
00185 det1 = pDD->idToDet(id1);
00186 LogDebug("") <<" First hit local x,y,z "<<(*it).localPosition()<<" det "<<id1.det()<<" subdet "<<id1.subdetId();
00187 LogDebug("") <<" First hit global "<<det1->toGlobal((*it).localPosition());
00188
00189
00190 it++;
00191 DetId id2 = (*it).geographicalId();
00192 det2 = pDD->idToDet(id2);
00193 LogDebug("") <<" Second hit local x,y,z "<<(*it).localPosition()<<" det "<<id2.det()<<" subdet "<<id2.subdetId();
00194 LogDebug("") <<" Second hit global "<<det2->toGlobal((*it).localPosition());
00195
00196
00197
00198
00199 const GeomDet *det=0;
00200 for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
00201 TrajectoryStateOnSurface t= transformer_.transientState((*MyS).startingState(), &(det->surface()), &(*theMagField));
00202
00203
00204
00205 LogDebug("")<<" ElectronPixelSeed outermost state position: "<<t.globalPosition();
00206 LogDebug("")<<" ElectronPixelSeed outermost state momentum: "<<t.globalMomentum();
00207 edm::Ref<SuperClusterCollection> theClus=(*MyS).superCluster();
00208 LogDebug("")<<" ElectronPixelSeed superCluster energy: "<<theClus->energy()<<", position: "<<theClus->position();
00209 LogDebug("")<<" ElectronPixelSeed outermost state Pt: "<<t.globalMomentum().perp();
00210 LogDebug("")<<" ElectronPixelSeed supercluster Et: "<<theClus->energy()*sin(2.*atan(exp(-theClus->position().eta())));
00211 LogDebug("")<<" ElectronPixelSeed outermost momentum direction eta: "<<t.globalMomentum().eta();
00212 LogDebug("")<<" ElectronPixelSeed supercluster eta: "<<theClus->position().eta();
00213 LogDebug("")<<" ElectronPixelSeed seed charge: "<<(*MyS).getCharge();
00214 LogDebug("")<<" ElectronPixelSeed E/p: "<<theClus->energy()/t.globalMomentum().mag();
00215
00216
00217
00218
00219
00220 int charge = int((*MyS).getCharge());
00221 GlobalPoint xmeas(theClus->position().x(),theClus->position().y(),theClus->position().z());
00222 GlobalPoint vprim(theBeamSpot->position().x(),theBeamSpot->position().y(),theBeamSpot->position().z());
00223 float energy = theClus->energy();
00224
00225 FreeTrajectoryState fts = myFTS(&(*theMagField),xmeas, vprim,
00226 energy, charge);
00227
00228
00229
00230 PerpendicularBoundPlaneBuilder bpb;
00231 TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));
00232
00233
00234
00235
00236
00237 it=r.first;
00238 DetId id=(*it).geographicalId();
00239 const GeomDet *geomdet=pDD->idToDet((*it).geographicalId());
00240 LocalPoint lp=(*it).localPosition();
00241 GlobalPoint hitPos=geomdet->surface().toGlobal(lp);
00242
00243 TrajectoryStateOnSurface tsos1;
00244 tsos1 = prop1stLayer->propagate(tsos,geomdet->surface()) ;
00245
00246 if (tsos1.isValid()) {
00247
00248 std::pair<bool,double> est;
00249
00250
00251 float SCl_phi = xmeas.phi();
00252 float localDphi = SCl_phi-hitPos.phi();
00253 if(localDphi>CLHEP::pi)localDphi-=(2*CLHEP::pi);
00254 if(localDphi<-CLHEP::pi)localDphi+=(2*CLHEP::pi);
00255 if(fabs(localDphi)>2.5)continue;
00256
00257 phi1 = hitPos.phi();
00258 dphi1 = hitPos.phi() - tsos1.globalPosition().phi();
00259 rz1 = hitPos.perp();
00260 drz1 = hitPos.perp() - tsos1.globalPosition().perp();
00261 if (id.subdetId()%2==1) {
00262 drz1 = hitPos.z() - tsos1.globalPosition().z();
00263 rz1 = hitPos.z();
00264 }
00265
00266
00267 it++;
00268 DetId id2=(*it).geographicalId();
00269 const GeomDet *geomdet2=pDD->idToDet((*it).geographicalId());
00270 TrajectoryStateOnSurface tsos2;
00271
00272
00273 double pxHit1z = hitPos.z();
00274 double pxHit1x = hitPos.x();
00275 double pxHit1y = hitPos.y();
00276 double r1diff = (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) + (pxHit1y-vprim.y())*(pxHit1y-vprim.y());
00277 r1diff=sqrt(r1diff);
00278 double r2diff = (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) + (xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y);
00279 r2diff=sqrt(r2diff);
00280 double zVertexPred = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;
00281
00282 GlobalPoint vertexPred(vprim.x(),vprim.y(),zVertexPred);
00283
00284 FreeTrajectoryState fts2 = myFTS(&(*theMagField),hitPos,vertexPred,energy, charge);
00285 tsos2 = prop2ndLayer->propagate(fts2,geomdet2->surface()) ;
00286
00287 if (tsos2.isValid()) {
00288 LocalPoint lp2=(*it).localPosition();
00289 GlobalPoint hitPos2=geomdet2->surface().toGlobal(lp2);
00290 phi2 = hitPos2.phi();
00291 dphi2 = hitPos2.phi() - tsos2.globalPosition().phi();
00292 rz2 = hitPos2.perp();
00293 drz2 = hitPos2.perp() - tsos2.globalPosition().perp();
00294 if (id2.subdetId()%2==1) {
00295 rz2 = hitPos2.z();
00296 drz2 = hitPos2.z() - tsos2.globalPosition().z();
00297 }
00298 }
00299
00300 }
00301
00302
00303
00304 histpt_->Fill(t.globalMomentum().perp());
00305 histetclu_->Fill(theClus->energy()*sin(2.*atan(exp(-theClus->position().eta()))));
00306 histeta_->Fill(t.globalMomentum().eta());
00307 histetaclu_->Fill(theClus->position().eta());
00308 histq_->Fill((*MyS).getCharge());
00309 histeoverp_->Fill(theClus->energy()/t.globalMomentum().mag());
00310
00311 if (is<10) {
00312 superclusterEnergy[is] = theClus->energy();
00313 superclusterEta[is] = theClus->position().eta();
00314 superclusterPhi[is] = theClus->position().phi();
00315 superclusterEt[is] = theClus->energy()*sin(2.*atan(exp(-theClus->position().eta())));
00316 seedMomentum[is] = t.globalMomentum().mag();
00317 seedEta[is] = t.globalMomentum().eta();
00318 seedPhi[is] = t.globalMomentum().phi();
00319 seedPt[is] = t.globalMomentum().perp();
00320 seedQ[is] = (*MyS).getCharge();
00321 seedSubdet1[is] = id1.subdetId();
00322 switch (seedSubdet1[is]) {
00323 case 1:
00324 {
00325 PXBDetId pxbid1( id1.rawId() );
00326 seedLayer1[is] = pxbid1.layer();
00327 seedSide1[is] = 0;
00328 break;
00329 }
00330 case 2:
00331 {
00332 PXFDetId pxfid1( id1.rawId() );
00333 seedLayer1[is] = pxfid1.disk();
00334 seedSide1[is] = pxfid1.side();
00335 break;
00336 }
00337 case 3:
00338 {
00339 TIBDetId tibid1( id1.rawId() );
00340 seedLayer1[is] = tibid1.layer();
00341 seedSide1[is] = 0;
00342 break;
00343 }
00344 case 4:
00345 {
00346 TIDDetId tidid1( id1.rawId() );
00347 seedLayer1[is] = tidid1.wheel();
00348 seedSide1[is] = tidid1.side();
00349 break;
00350 }
00351 case 5:
00352 {
00353 TOBDetId tobid1( id1.rawId() );
00354 seedLayer1[is] = tobid1.layer();
00355 seedSide1[is] = 0;
00356 break;
00357 }
00358 case 6:
00359 {
00360 TECDetId tecid1( id1.rawId() );
00361 seedLayer1[is] = tecid1.wheel();
00362 seedSide1[is] = tecid1.side();
00363 break;
00364 }
00365 }
00366 seedPhi1[is] = phi1;
00367 seedRz1[is] = rz1;
00368 seedDphi1[is] = dphi1;
00369 seedDrz1[is] = drz1;
00370 seedSubdet2[is] = id2.subdetId();
00371 switch (seedSubdet2[is]) {
00372 case 1:
00373 {
00374 PXBDetId pxbid2( id2.rawId() );
00375 seedLayer2[is] = pxbid2.layer();
00376 seedSide2[is] = 0;
00377 break;
00378 }
00379 case 2:
00380 {
00381 PXFDetId pxfid2( id2.rawId() );
00382 seedLayer2[is] = pxfid2.disk();
00383 seedSide2[is] = pxfid2.side();
00384 break;
00385 }
00386 case 3:
00387 {
00388 TIBDetId tibid2( id2.rawId() );
00389 seedLayer2[is] = tibid2.layer();
00390 seedSide2[is] = 0;
00391 break;
00392 }
00393 case 4:
00394 {
00395 TIDDetId tidid2( id2.rawId() );
00396 seedLayer2[is] = tidid2.wheel();
00397 seedSide2[is] = tidid2.side();
00398 break;
00399 }
00400 case 5:
00401 {
00402 TOBDetId tobid2( id2.rawId() );
00403 seedLayer2[is] = tobid2.layer();
00404 seedSide2[is] = 0;
00405 break;
00406 }
00407 case 6:
00408 {
00409 TECDetId tecid2( id2.rawId() );
00410 seedLayer2[is] = tecid2.wheel();
00411 seedSide2[is] = tecid2.side();
00412 break;
00413 }
00414 }
00415 seedDphi2[is] = dphi2;
00416 seedDrz2[is] = drz2;
00417 seedPhi2[is] = phi2;
00418 seedRz2[is] = rz2;
00419 }
00420
00421 is++;
00422
00423 }
00424
00425 histnbseeds_->Fill(elSeeds.product()->size());
00426
00427
00428
00429 edm::Handle<SuperClusterCollection> clusters;
00430
00431 e.getByLabel("correctedHybridSuperClusters", clusters);
00432 histnbclus_->Fill(clusters.product()->size());
00433 if (clusters.product()->size()>0) histnrseeds_->Fill(elSeeds.product()->size());
00434
00435
00436
00437 edm::Handle<edm::HepMCProduct> HepMCEvt;
00438
00439
00440 e.getByLabel("source", "", HepMCEvt);
00441
00442 const HepMC::GenEvent* MCEvt = HepMCEvt->GetEvent();
00443 HepMC::GenParticle* genPc=0;
00444 HepMC::FourVector pAssSim;
00445 int ip=0;
00446 for (HepMC::GenEvent::particle_const_iterator partIter = MCEvt->particles_begin();
00447 partIter != MCEvt->particles_end(); ++partIter) {
00448
00449 for (HepMC::GenEvent::vertex_const_iterator vertIter = MCEvt->vertices_begin();
00450 vertIter != MCEvt->vertices_end(); ++vertIter) {
00451
00452
00453 HepMC::GenVertex * creation = (*partIter)->production_vertex();
00454
00455 HepMC::FourVector momentum = (*partIter)->momentum();
00456
00457 int id = (*partIter)->pdg_id();
00458 LogDebug("") << "MC particle id " << id << ", creationVertex " << (*creation) << " cm, initialMomentum " << momentum.mag() << " GeV/c" << std::endl;
00459
00460 if (id == 11 || id == -11) {
00461
00462
00463 HepMC::GenParticle* mother = 0;
00464 if ( (*partIter)->production_vertex() ) {
00465 if ( (*partIter)->production_vertex()->particles_begin(HepMC::parents) !=
00466 (*partIter)->production_vertex()->particles_end(HepMC::parents))
00467 mother = *((*partIter)->production_vertex()->particles_begin(HepMC::parents));
00468 }
00469 if ( ((mother == 0) || ((mother != 0) && (mother->pdg_id() == 23))
00470 || ((mother != 0) && (mother->pdg_id() == 32))
00471 || ((mother != 0) && (fabs(mother->pdg_id()) == 24)))) {
00472 genPc=(*partIter);
00473 pAssSim = genPc->momentum();
00474
00475 if (pAssSim.perp()> 100. || fabs(pAssSim.eta())> 2.5) continue;
00476
00477
00478 bool okSeedFound = false;
00479 double seedOkRatio = 999999.;
00480
00481
00482 reco::ElectronPixelSeed bestElectronSeed;
00483 for( ElectronPixelSeedCollection::const_iterator gsfIter= (*elSeeds).begin(); gsfIter != (*elSeeds).end(); ++gsfIter) {
00484
00485 range r=gsfIter->recHits();
00486 const GeomDet *det=0;
00487 for (TrackingRecHitCollection::const_iterator rhits=r.first; rhits!=r.second; rhits++) det = pDD->idToDet(((*rhits)).geographicalId());
00488 TrajectoryStateOnSurface t= transformer_.transientState(gsfIter->startingState(), &(det->surface()), &(*theMagField));
00489
00490 float eta = t.globalMomentum().eta();
00491 float phi = t.globalMomentum().phi();
00492 float p = t.globalMomentum().mag();
00493 double deltaR = sqrt(pow((eta-pAssSim.eta()),2) + pow((phi-pAssSim.phi()),2));
00494 if ( deltaR < 0.05 ){
00495
00496
00497 double tmpSeedRatio = p/pAssSim.t();
00498 if ( fabs(tmpSeedRatio-1) < fabs(seedOkRatio-1) ) {
00499 seedOkRatio = tmpSeedRatio;
00500 bestElectronSeed=*gsfIter;
00501 okSeedFound = true;
00502 }
00503
00504 }
00505 }
00506
00507
00508 if (okSeedFound){
00509
00510 histptMC_->Fill(pAssSim.perp());
00511 histetaMC_->Fill(pAssSim.eta());
00512 histeMC_->Fill(pAssSim.rho());
00513 if (ip<10) {
00514 mcEnergy[ip] = pAssSim.rho();
00515 mcEta[ip] = pAssSim.eta();
00516 mcPhi[ip] = pAssSim.phi();
00517 mcPt[ip] = pAssSim.perp();
00518 mcQ[ip] = ((id == 11) ? -1.: +1.);
00519 }
00520 }
00521 }
00522
00523 }
00524
00525 }
00526
00527 ip++;
00528
00529 }
00530
00531 tree_->Fill();
00532
00533 }
1.5.4