Changeset 534d945 in git for validation
- Timestamp:
- Apr 21, 2017, 2:36:33 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 78d1846
- Parents:
- e605a99
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
validation/Validation.cpp
re605a99 r534d945 253 253 if(eta > etamax || eta < etamin ) continue; 254 254 255 //cout<<"b parton: "<<pt<<endl;256 255 if (particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax ) 256 //if (TMath::Abs(particle->PID) == pdgID && (particle->Status>20 && particle->Status <30) && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax ) 257 257 { 258 258 // Loop over all reco object in event … … 261 261 recoObj = (T*)branchReco->At(j); 262 262 recoMomentum = recoObj->P4(); 263 // this is simply to avoid warnings from initial state particle264 // having infite rapidity ...265 263 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue; 266 264 … … 317 315 } 318 316 317 319 318 template<typename T> 320 319 TH1D* GetEffEta(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader) … … 363 362 364 363 if (particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax ) 364 //if (TMath::Abs(particle->PID) == pdgID && (particle->Status>20 && particle->Status <30) && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax ) 365 365 { 366 366 // Loop over all reco object in event … … 422 422 423 423 424 //------------------------------------------------------------------------------ 425 426 template<typename T> 427 TH1D* GetJetEffPt(TClonesArray *branchJet, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader) 428 { 429 430 cout << "** Computing Efficiency of reconstructing "<< branchJet->GetName() << " with PID " << pdgID << endl; 431 432 Long64_t allEntries = treeReader->GetEntries(); 433 434 T *recoObj; 435 436 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum; 437 438 Float_t pt, eta; 439 Long64_t entry; 440 441 Int_t j; 442 443 TH1D *histGenPt = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax)); 444 TH1D *histRecoPt = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax)); 445 446 histGenPt->SetDirectory(0); 447 histRecoPt->SetDirectory(0); 448 449 BinLogX(histGenPt); 450 BinLogX(histRecoPt); 451 452 // Loop over all events 453 for(entry = 0; entry < allEntries; ++entry) 454 { 455 // Load selected branches with data from specified event 456 treeReader->ReadEntry(entry); 457 458 // Loop over all reco object in event 459 for(j = 0; j < branchJet->GetEntriesFast(); ++j) 460 { 461 recoObj = (T*)branchJet->At(j); 462 recoMomentum = recoObj->P4(); 463 pt = recoMomentum.Pt(); 464 eta = TMath::Abs(recoMomentum.Eta()); 465 Jet *jet = (Jet *)recoObj; 466 467 if(eta > etamax || eta < etamin ) continue; 468 if(pt < ptmin || pt > ptmax) continue; 469 470 Int_t flavor = jet->Flavor; 471 if(flavor == 21) flavor = 0; 472 473 if(TMath::Abs(pdgID) == 1) 474 { 475 476 if(flavor < 4) 477 { 478 histGenPt->Fill(pt); 479 if( jet->BTag & (1 << 0) ) histRecoPt->Fill(pt); 480 } 481 } 482 if(TMath::Abs(pdgID) == 4) 483 { 484 if(flavor == 4) 485 { 486 histGenPt->Fill(pt); 487 if( jet->BTag & (1 << 0) ) histRecoPt->Fill(pt); 488 } 489 } 490 if(TMath::Abs(pdgID) == 5) 491 { 492 if(flavor == 5) 493 { 494 histGenPt->Fill(pt); 495 if( jet->BTag & (1 << 0) ) histRecoPt->Fill(pt); 496 } 497 } 498 } 499 500 } 501 502 histRecoPt->Sumw2(); 503 histGenPt->Sumw2(); 504 505 histRecoPt->Divide(histGenPt); 506 histRecoPt->Scale(100.); 507 508 return histRecoPt; 509 } 510 511 // ------------------------------------------------------------------------------------------------------------------------------------------------------ 512 513 template<typename T> 514 TH1D* GetJetEffEta(TClonesArray *branchJet, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader) 515 { 516 517 cout << "** Computing Efficiency of reconstructing "<< branchJet->GetName() << " with PID " << pdgID << endl; 518 519 Long64_t allEntries = treeReader->GetEntries(); 520 521 T *recoObj; 522 523 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum; 524 525 Float_t pt, eta; 526 Long64_t entry; 527 528 Int_t j; 529 530 TH1D *histGenEta = new TH1D(name+" gen spectra Eta",name+" gen spectra", Nbins, etamin, etamax); 531 TH1D *histRecoEta = new TH1D(name+" reco spectra Eta",name+" reco spectra", Nbins, etamin, etamax); 532 533 histGenEta->SetDirectory(0); 534 histRecoEta->SetDirectory(0); 535 // Loop over all events 536 for(entry = 0; entry < allEntries; ++entry) 537 { 538 // Load selected branches with data from specified event 539 treeReader->ReadEntry(entry); 540 541 // Loop over all reco object in event 542 for(j = 0; j < branchJet->GetEntriesFast(); ++j) 543 { 544 recoObj = (T*)branchJet->At(j); 545 recoMomentum = recoObj->P4(); 546 pt = recoMomentum.Pt(); 547 eta = recoMomentum.Eta(); 548 Jet *jet = (Jet *)recoObj; 549 550 if(eta > etamax || eta < etamin ) continue; 551 if(pt < ptmin || pt > ptmax) continue; 552 553 Int_t flavor = jet->Flavor; 554 if(flavor == 21) flavor = 0; 555 556 if(TMath::Abs(pdgID) == 1) 557 { 558 if(flavor == 1 || flavor == 21) 559 { 560 histGenEta->Fill(eta); 561 if( jet->BTag & (1 << 0) ) histRecoEta->Fill(eta); 562 } 563 } 564 if(TMath::Abs(pdgID) == 4) 565 { 566 if(flavor == 4) 567 { 568 histGenEta->Fill(eta); 569 if( jet->BTag & (1 << 0) ) histRecoEta->Fill(eta); 570 } 571 } 572 if(TMath::Abs(pdgID) == 5) 573 { 574 if(flavor == 5) 575 { 576 histGenEta->Fill(eta); 577 if( jet->BTag & (1 << 0) ) histRecoEta->Fill(eta); 578 } 579 } 580 } 581 582 } 583 584 585 histRecoEta->Sumw2(); 586 histGenEta->Sumw2(); 587 588 histRecoEta->Divide(histGenEta); 589 histRecoEta->Scale(100.); 590 591 return histRecoEta; 592 } 593 594 595 // ----------------------------------------------------------------------------------------------------------------------------------------------------- 424 596 425 597 template<typename T> … … 2468 2640 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.pdf","pdf"); 2469 2641 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.png","png"); 2470 2471 2642 2472 2643 ///////////////////////////////////////// 2473 2644 // B-jets Efficiency/ mistag rates /// … … 2488 2659 { 2489 2660 2490 h_recbjet_eff_pt = GetEffPt<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet); 2661 h_recbjet_eff_pt = GetJetEffPt<Jet>(branchPFBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet); 2662 //h_recbjet_eff_pt = GetEffPt<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet); 2491 2663 gr_recbjet_eff_pt[k] = TGraphErrors(h_recbjet_eff_pt); 2492 2664 … … 2504 2676 for (k = 0; k < ptVals.size(); k++) 2505 2677 { 2506 h_recbjet_eff_eta = GetEffEta<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet); 2678 h_recbjet_eff_eta = GetJetEffEta<Jet>(branchPFBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet); 2679 //h_recbjet_eff_eta = GetEffEta<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet); 2507 2680 gr_recbjet_eff_eta[k] = TGraphErrors(h_recbjet_eff_eta); 2508 2681 … … 2551 2724 { 2552 2725 2553 h_recbjet_cmis_pt = GetEffPt<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet); 2726 h_recbjet_cmis_pt = GetJetEffPt<Jet>(branchPFCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet); 2727 //h_recbjet_cmis_pt = GetEffPt<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet); 2554 2728 gr_recbjet_cmis_pt[k] = TGraphErrors(h_recbjet_cmis_pt); 2555 2729 … … 2567 2741 for (k = 0; k < ptVals.size(); k++) 2568 2742 { 2569 h_recbjet_cmis_eta = GetEffEta<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet); 2743 h_recbjet_cmis_eta = GetJetEffEta<Jet>(branchPFCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet); 2744 //h_recbjet_cmis_eta = GetEffEta<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet); 2570 2745 gr_recbjet_cmis_eta[k] = TGraphErrors(h_recbjet_cmis_eta); 2571 2746 … … 2614 2789 { 2615 2790 2616 h_recbjet_lmis_pt = GetEffPt<Jet>(branchJet, branchParticleJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet); 2791 h_recbjet_lmis_pt = GetJetEffPt<Jet>(branchJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet); 2792 //h_recbjet_lmis_pt = GetEffPt<Jet>(branchJet, branchParticleJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet); 2617 2793 gr_recbjet_lmis_pt[k] = TGraphErrors(h_recbjet_lmis_pt); 2618 2794 … … 2630 2806 for (k = 0; k < ptVals.size(); k++) 2631 2807 { 2632 h_recbjet_lmis_eta = GetEffEta<Jet>(branchJet, branchParticleJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet); 2808 h_recbjet_lmis_eta = GetJetEffEta<Jet>(branchJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet); 2809 //h_recbjet_lmis_eta = GetEffEta<Jet>(branchJet, branchParticleJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet); 2633 2810 gr_recbjet_lmis_eta[k] = TGraphErrors(h_recbjet_lmis_eta); 2634 2811 … … 2643 2820 mg_recbjet_lmis_pt->Draw("APE"); 2644 2821 2645 DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 1 0.0, "p_{T} [GeV]", "light - mistag rate (%)", true, false);2822 DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 1.0, "p_{T} [GeV]", "light - mistag rate (%)", true, false); 2646 2823 2647 2824 leg_recbjet_lmis_pt->Draw(); … … 2655 2832 2656 2833 mg_recbjet_lmis_eta->Draw("APE"); 2657 DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 1 0.0, " #eta ", "light - mistag rate (%)", false, false);2834 DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 1.0, " #eta ", "light - mistag rate (%)", false, false); 2658 2835 leg_recbjet_lmis_eta->Draw(); 2659 2836 pave->Draw(); … … 2837 3014 2838 3015 fout->Write(); 2839 2840 2841 3016 2842 3017
Note:
See TracChangeset
for help on using the changeset viewer.