Fork me on GitHub

Changeset 6563d4d in git


Ignore:
Timestamp:
Aug 23, 2016, 2:37:13 PM (8 years ago)
Author:
Alexandre Mertens <alexandre.mertens@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
20ca0cd, df315dc
Parents:
9f8d4e7
Message:

added met plot

File:
1 edited

Legend:

Unmodified
Added
Removed
  • examples/Validation.C

    r9f8d4e7 r6563d4d  
    4949    double ptmin;
    5050    double ptmax;
     51    double xmin;
     52    double xmax;
    5153    TString obj;
    5254
    5355    resolPlot();
    5456    resolPlot(double ptdown, double ptup, TString object);
    55     void set(double ptdown, double ptup, TString object);
     57    void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
    5658    void print(){std::cout << ptmin << std::endl;}
    5759};
     
    6769}
    6870
    69 void resolPlot::set(double ptdown, double ptup, TString object){
     71void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax){
    7072    ptmin = ptdown;
    7173    ptmax = ptup;
    7274    obj = object;
    7375
    74     cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 500,  0, 2);
    75     fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 500, 0.4, 0.4);
    76 
    77 }
    78 
    79 void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj)
     76    cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 500,  xmin, xmax);
     77    fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 500,  xmin, xmax);
     78
     79}
     80
     81void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
    8082{
    8183    double width;
     
    8991        ptdown = TMath::Power(10,ptmin + i * width );
    9092        ptup = TMath::Power(10,ptmin + (i+1) * width );
    91         ptemp.set(ptdown, ptup, obj);
     93        ptemp.set(ptdown, ptup, obj, xmin, xmax);
    9294        histos->push_back(ptemp);
    9395    }
     
    362364  }
    363365}
     366
     367void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
     368{
     369
     370  Long64_t allEntries = treeReader->GetEntries();
     371
     372  cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
     373
     374  MissingET *met;
     375  ScalarHT *scalarHT;
     376
     377  Long64_t entry;
     378
     379  Int_t bin;
     380  Double_t ht;
     381
     382  // Loop over all events
     383  for(entry = 0; entry < allEntries; ++entry)
     384  {
     385    // Load selected branches with data from specified event
     386    treeReader->ReadEntry(entry);
     387
     388    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
     389
     390    if (branchJet->GetEntriesFast() > 1)
     391    {
     392      met = (MissingET*) branchMet->At(0);
     393      scalarHT = (ScalarHT*) branchScalarHT->At(0);
     394      ht = scalarHT->HT;
     395
     396      for (bin = 0; bin < Nbins; bin++)
     397      {
     398        if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
     399        {
     400          histos->at(bin).cenResolHist->Fill(met->P4().Px());
     401          histos->at(bin).fwdResolHist->Fill(met->P4().Py());
     402        }
     403      }
     404    }
     405  }
     406}
     407
     408
    364409std::pair<Double_t, Double_t> GausFit(TH1* hist)
    365410{
     
    374419
    375420
    376 TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central)
     421TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
    377422{
    378423    Int_t bin;
     
    389434            std::cout << " entries : " << histos->at(bin).cenResolHist->GetEntries() << std::endl;
    390435            std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
    391             gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
    392             gr.SetPointError(count,0, sigvalues.second);
     436            if (rms == true)
     437            {
     438              gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
     439              gr.SetPointError(count,0, sigvalues.second); // to correct
     440            }
     441            else
     442            {
     443              gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
     444              gr.SetPointError(count,0, sigvalues.second);
     445            }
    393446            count++;
    394447        }
     
    551604  TClonesArray *branchPFJet = treeReader->UseBranch("Jet");
    552605  TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet");
     606  TClonesArray *branchScalarHT = treeReader->UseBranch("ScalarHT");
     607  TClonesArray *branchMet = treeReader->UseBranch("MissingET");
    553608
    554609#ifdef ELECTRON
     
    811866  gr_pfjets.SetName("pfJet");
    812867
    813 
    814   // PFJets Energy Resolution
     868  // CaloJets Energy Resolution
    815869  std::vector<resolPlot> plots_calojets;
    816870  HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
     
    819873  gr_calojets.SetName("caloJet");
    820874
     875  // MET Resolution vs HT
     876  std::vector<resolPlot> plots_met;
     877  HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -200, 200);
     878  GetMetres( &plots_met, branchScalarHT, branchMet, branchPFJet, treeReader);
     879  TGraphErrors gr_met = EresGraph(&plots_met, true);
     880  gr_calojets.SetName("MET");
     881
    821882  // Canvas
    822883
     
    882943  C_btag1->SaveAs(btagEff+".eps");
    883944  C_jet->SaveAs(jetRes+".eps");
     945
     946  TString metRes = "MetRes";
     947  TCanvas *C_met = new TCanvas(metRes,metRes, 1000, 500);
     948
     949  TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
     950  TLegend *leg_met = new TLegend(0.52,0.7,0.9,0.9);
     951
     952  addGraph(mg_met, &gr_met, leg_met, 3);
     953
     954  mg_met->Draw("ACX");
     955  leg_met->Draw();
     956
     957  DrawAxis(mg_met, leg_met, 100);
     958
     959  C_met->SaveAs(metRes+".eps");
     960
    884961
    885962#endif
     
    9341011      plots_pfjets.at(bin).cenResolHist->Write();
    9351012      plots_calojets.at(bin).cenResolHist->Write();
     1013      plots_met.at(bin).cenResolHist->Write();
    9361014  }
    9371015  histos_btag.first->Write();
     
    9421020  C_tautag1->Write();
    9431021  C_jet->Write();
     1022  C_met->Write();
    9441023#endif
    9451024
Note: See TracChangeset for help on using the changeset viewer.