Fork me on GitHub

Changeset 534d945 in git for validation/Validation.cpp


Ignore:
Timestamp:
Apr 21, 2017, 2:36:33 PM (8 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
78d1846
Parents:
e605a99
Message:

fixed btag eff plots

File:
1 edited

Legend:

Unmodified
Added
Removed

  • validation/Validation.cpp

    re605a99 r534d945  
    253253      if(eta > etamax || eta < etamin ) continue;
    254254
    255       //cout<<"b parton: "<<pt<<endl;
    256255      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 )
    257257      {
    258258        // Loop over all reco object in event
    … …  
    261261          recoObj = (T*)branchReco->At(j);
    262262          recoMomentum = recoObj->P4();
    263           // this is simply to avoid warnings from initial state particle
    264           // having infite rapidity ...
    265263        //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
    266264
    … …  
    317315}
    318316
     317
    319318template<typename T>
    320319TH1D* GetEffEta(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
    … …  
    363362
    364363      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 )
    365365      {
    366366        // Loop over all reco object in event
    … …  
    422422
    423423
     424//------------------------------------------------------------------------------
     425
     426template<typename T>
     427TH1D* 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
     513template<typename T>
     514TH1D* 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// -----------------------------------------------------------------------------------------------------------------------------------------------------
    424596
    425597template<typename T>
    … …  
    24682640    c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.pdf","pdf");
    24692641    c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.png","png");
    2470 
    2471 
     2642   
    24722643    /////////////////////////////////////////
    24732644    // B-jets  Efficiency/ mistag rates   ///
    … …  
    24882659    {
    24892660
    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);
    24912663       gr_recbjet_eff_pt[k] = TGraphErrors(h_recbjet_eff_pt);
    24922664
    … …  
    25042676    for (k = 0; k < ptVals.size(); k++)
    25052677    {
    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);
    25072680       gr_recbjet_eff_eta[k] = TGraphErrors(h_recbjet_eff_eta);
    25082681
    … …  
    25512724    {
    25522725
    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);
    25542728       gr_recbjet_cmis_pt[k] = TGraphErrors(h_recbjet_cmis_pt);
    25552729
    … …  
    25672741    for (k = 0; k < ptVals.size(); k++)
    25682742    {
    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);
    25702745       gr_recbjet_cmis_eta[k] = TGraphErrors(h_recbjet_cmis_eta);
    25712746
    … …  
    26142789    {
    26152790
    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);
    26172793       gr_recbjet_lmis_pt[k] = TGraphErrors(h_recbjet_lmis_pt);
    26182794
    … …  
    26302806    for (k = 0; k < ptVals.size(); k++)
    26312807    {
    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);
    26332810       gr_recbjet_lmis_eta[k] = TGraphErrors(h_recbjet_lmis_eta);
    26342811
    … …  
    26432820    mg_recbjet_lmis_pt->Draw("APE");
    26442821   
    2645     DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 10.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);
    26462823   
    26472824    leg_recbjet_lmis_pt->Draw();
    … …  
    26552832
    26562833    mg_recbjet_lmis_eta->Draw("APE");
    2657     DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 10.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);
    26582835    leg_recbjet_lmis_eta->Draw();
    26592836    pave->Draw();
    … …  
    28373014
    28383015  fout->Write();
    2839 
    2840 
    28413016
    28423017
Note: See TracChangeset for help on using the changeset viewer.