Skip to content
Snippets Groups Projects
Commit 812ae99c authored by Vardan Tadevosyan's avatar Vardan Tadevosyan Committed by Stephen A. Wood
Browse files

Add a script for event-by-event comparison of energies of the best

track in Engine and the golden track in hcana.

      Add directory examples/gtrk_e holding the code.
parent c4736336
No related branches found
No related tags found
No related merge requests found
#define hcana_class_cxx
#define engine_class_cxx
#include <TTree.h>
#include <TFile.h>
#include <TROOT.h>
#include <iostream>
#include <TMath.h>
#include <TLegend.h>
#include "hcana_class.h"
#include "engine_class.h"
using namespace std;
void comp_gtrk_e() {
// Compare Ecalo of the golden track from Engine and hcana.
// hcana tree
hcana_class* hcana_tree = new hcana_class(0);
Long64_t hcana_nentries = hcana_tree->fChain->GetEntriesFast();
cout << "hcana_nentries = " << hcana_nentries << endl;
// hcana_tree->fChain->SetBranchStatus("*",0); // disable all branches
// hcana_tree->fChain->SetBranchStatus("*H.gold.e",1); // activate branchname
// hcana_tree->fChain->SetBranchStatus("*H.gold.p",1);
// hcana_tree->fChain->SetBranchStatus("*H.gold.index",1);
// hcana_tree->fChain->SetBranchStatus("*H.tr.n",1);
// hcana_tree->fChain->SetBranchStatus("*g.evtyp",1);
// hcana_tree->fChain->SetBranchStatus("*fEvtHdr.fEvtType",1);
// hcana_tree->fChain->SetBranchStatus("*fEvtHdr.fEvtNum",1);
// Engine tree
engine_class* engine_tree = new engine_class(0);
Long64_t engine_nentries = engine_tree->fChain->GetEntriesFast();
cout << "engine_nentries = " << engine_nentries << endl;
// engine_tree->fChain->SetBranchStatus("*",0);
// engine_tree->fChain->SetBranchStatus("*hstrk_et",1);
// engine_tree->fChain->SetBranchStatus("*hsp",1);
// engine_tree->fChain->SetBranchStatus("*ev_type",1);
// engine_tree->fChain->SetBranchStatus("*eventid",1);
// engine_tree->fChain->SetBranchStatus("*hstrk_in",1);
// engine_tree->fChain->SetBranchStatus("*hntracks",1);
// Histograms.
TH1F* he_hcana = new TH1F("e_hcana", "Golden track Ecalo",
100, 0.011, 1.8);
he_hcana->SetLineColor(kRed);
TH1F* he_engine = (TH1F*) he_hcana->Clone("e_engine");
he_engine->SetLineColor(kBlue);
TH1F* ee_diff = new TH1F("ee_diff",
"Event by event hcana - Engine Ecalo difference", 200,-1.,1.);
TH1F* pp_diff = new TH1F("pp_diff","Event by event P difference",
200,-1.,1.);
TH2F* pvse_diff = new TH2F("pvse_diff","Event by event P vs Ecalo difference",
200,-1.,1., 200,-100.,100.);
pvse_diff->SetMarkerStyle(3);
TH2F* evse = new TH2F("evse","Event by event hcana vs Engine Ecalo",
210,-0.1,2., 210,-0.1,2.);
evse->SetMarkerStyle(3);
// Loop over entries.
Long64_t nentries = TMath::Min(hcana_nentries, engine_nentries);
cout << "nentries = " << nentries << endl;
Long64_t jentry_hcana = 0;
Long64_t jentry_engine = 0;
for (Long64_t jentry=0; jentry<nentries; jentry++) {
hcana_tree->fChain->GetEntry(jentry_hcana++);
engine_tree->fChain->GetEntry(jentry_engine++);
while (hcana_tree->fEvtHdr_fEvtNum > engine_tree->eventid &&
jentry_engine < engine_nentries) {
engine_tree->fChain->GetEntry(jentry_engine++);
//cout<< hcana_tree->fEvtHdr_fEvtNum << " " << engine_tree->eventid << endl;
}
while (hcana_tree->fEvtHdr_fEvtNum < engine_tree->eventid &&
jentry_hcana < hcana_nentries) {
hcana_tree->fChain->GetEntry(jentry_hcana++);
//cout<< hcana_tree->fEvtHdr_fEvtNum << " " << engine_tree->eventid << endl;
}
if (hcana_tree->fEvtHdr_fEvtType==1 && engine_tree->ev_type==1
&& hcana_tree->H_gold_index>-1 && engine_tree->hstrk_in>0
&& hcana_tree->H_gold_index == engine_tree->hstrk_in - 1) {
he_hcana->Fill(hcana_tree->H_gold_e, 1.);
//cout << hcana_tree->H_gold_e << " " << hcana_tree->g_evtyp << " "
// << hcana_tree->fEvtHdr_fEvtType << endl;
he_engine->Fill(engine_tree->hstrk_et, 1.);
// cout << engine_tree->hsshtrk << " " << engine_tree->ev_type << " "
// << engine_tree->eventid << endl;
ee_diff->Fill(hcana_tree->H_gold_e - engine_tree->hstrk_et, 1.);
pp_diff->Fill(hcana_tree->H_gold_p - engine_tree->hsp, 1.);
pvse_diff->Fill(hcana_tree->H_gold_e - engine_tree->hstrk_et,
hcana_tree->H_gold_p - engine_tree->hsp, 1.);
evse->Fill(engine_tree->hstrk_et, hcana_tree->H_gold_e, 1.);
//Debug print out.
// if (hcana_tree->H_gold_e>0.6 && hcana_tree->H_gold_e<0.95 &&
// engine_tree->hstrk_et<0.1) {
// if (hcana_tree->H_gold_e>0.6 && hcana_tree->H_gold_e<0.95 &&
// TMath::Abs(engine_tree->hstrk_et - hcana_tree->H_gold_e) < 0.0001 ) {
// if (engine_tree->hstrk_et>0.6 && engine_tree->hstrk_et<0.95 &&
// hcana_tree->H_gold_e<0.1) {
if (TMath::Abs(engine_tree->hstrk_et - hcana_tree->H_gold_e) < 0.65 &&
TMath::Abs(engine_tree->hstrk_et - hcana_tree->H_gold_e) > 0.05 ) {
cout << "===========================================================\n";
cout << "Engine event " << engine_tree->eventid
<< ": hstrk_et = " << engine_tree->hstrk_et
<< ", hstrk_index = " << engine_tree->hstrk_in
<< ", hntracks = " << engine_tree->hntracks << endl;
cout << " hsxfp = " << engine_tree->hsxfp
<< ", hsyfp = " << engine_tree->hsyfp << endl;
cout << " hsxpfp= " << engine_tree->hsxpfp
<< ", hsypfp= " << engine_tree->hsypfp << endl;
cout << " hsp = " << engine_tree->hsp << endl;
cout << " At calorimeter:"
<< " x = " << engine_tree->hsxfp+338.69*engine_tree->hsxpfp
<< " y = " << engine_tree->hsyfp+338.69*engine_tree->hsypfp
<< endl;
cout << "-----------------------------------------------------------\n";
cout << "hcana event " << hcana_tree->fEvtHdr_fEvtNum
<< ": H_gold_e = " << hcana_tree->H_gold_e
<< " H_gold_index = " << hcana_tree->H_gold_index
<< " H_tr_n = " << hcana_tree->H_tr_n << endl;
cout << " H_gold_fp_x = " << hcana_tree->H_gold_fp_x
<< " H_gold_fp_y = " << hcana_tree->H_gold_fp_y << endl;
cout << " H_gold_fp_th= " << hcana_tree->H_gold_fp_th
<< " H_gold_fp_ph= " << hcana_tree->H_gold_fp_ph << endl;
cout << " H_gold_p = " << hcana_tree->H_gold_p << endl;
cout << " At calorimeter:"
<< " x = " << hcana_tree->H_gold_fp_x+338.69*hcana_tree->H_gold_fp_th
<< " y = " << hcana_tree->H_gold_fp_y+338.69*hcana_tree->H_gold_fp_ph
<< endl;
}
}
} //over entries
// Plot.
TCanvas* c1 = new TCanvas("c1","Golden Track E from hcana",600,750);
c1->Divide(1,2);
c1->cd(1);
he_hcana->Draw();
he_engine->Draw("same");
TLegend* leg = new TLegend(0.25, 0.70, 0.40, 0.82);
leg->AddEntry(he_hcana,"hcana","l");
leg->AddEntry(he_engine,"Engine","l");
leg->Draw();
// Difference between the histograms.
c1->cd(2);
TH1F* dif = he_hcana->Clone("difference");
dif->Reset();
dif->Add(he_hcana,he_engine,1.,-1.);
dif->SetTitle("hcana -- Engine difference");
dif->SetFillColor(kRed);
dif->SetLineColor(kRed);
dif->SetLineWidth(1);
dif->SetFillStyle(1111);
dif->Draw();
TCanvas* c2 = new TCanvas("c2","Event by event hcana - Engine Ecalo diff.",
600,750);
c2->Divide(1,2);
c2->cd(1);
gPad->SetLogy();
ee_diff->Draw();
c2->cd(2);
// gPad->SetLogy();
// pp_diff->Draw();
// pvse_diff->Draw();
evse->Draw();
}
#define engine_class_cxx
#include "engine_class.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void engine_class::Loop()
{
// In a ROOT session, you can do:
// Root > .L engine_class.C
// Root > engine_class t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
}
}
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Mon Sep 29 17:54:32 2014 by ROOT version 5.34/06
// from TTree h9010/
// found on file: hms52949_rz.root
//////////////////////////////////////////////////////////
#ifndef engine_class_h
#define engine_class_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
// Header file for the classes stored in the TTree if any.
// Fixed size dimensions of array or collections stored in the TTree if any.
class engine_class {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
Int_t fCurrent; //!current Tree number in a TChain
// Declaration of leaf types
Float_t hcer_npe;
Float_t hsp;
Float_t hse;
Float_t charge;
Float_t hsdelta;
Float_t hstheta;
Float_t hsphi;
Float_t w;
Float_t hszbeam;
Float_t hsdedx1;
Float_t hsbeta;
Float_t hsshtrk;
Float_t hsprtrk;
Float_t hsxfp;
Float_t hsyfp;
Float_t hsxpfp;
Float_t hsypfp;
Float_t hsytar;
Float_t hsxptar;
Float_t hsyptar;
Float_t hstart;
Float_t eventid;
Float_t ev_type;
Float_t gfrx_raw;
Float_t gfry_raw;
Float_t gbeam_x;
Float_t gbeam_y;
Float_t bpma_x;
Float_t bpma_y;
Float_t bpmb_x;
Float_t bpmb_y;
Float_t bpmc_x;
Float_t bpmc_y;
Float_t hseloss;
Float_t hntracks;
Float_t hcal_et;
Float_t hgoodsc;
Float_t hcal_e1;
Float_t hcal_e2;
Float_t hcal_e3;
Float_t hcal_e4;
Float_t hstrk_et;
Float_t hstrk_in;
// List of branches
TBranch *b_hcer_npe; //!
TBranch *b_hsp; //!
TBranch *b_hse; //!
TBranch *b_charge; //!
TBranch *b_hsdelta; //!
TBranch *b_hstheta; //!
TBranch *b_hsphi; //!
TBranch *b_w; //!
TBranch *b_hszbeam; //!
TBranch *b_hsdedx1; //!
TBranch *b_hsbeta; //!
TBranch *b_hsshtrk; //!
TBranch *b_hsprtrk; //!
TBranch *b_hsxfp; //!
TBranch *b_hsyfp; //!
TBranch *b_hsxpfp; //!
TBranch *b_hsypfp; //!
TBranch *b_hsytar; //!
TBranch *b_hsxptar; //!
TBranch *b_hsyptar; //!
TBranch *b_hstart; //!
TBranch *b_eventid; //!
TBranch *b_ev_type; //!
TBranch *b_gfrx_raw; //!
TBranch *b_gfry_raw; //!
TBranch *b_gbeam_x; //!
TBranch *b_gbeam_y; //!
TBranch *b_bpma_x; //!
TBranch *b_bpma_y; //!
TBranch *b_bpmb_x; //!
TBranch *b_bpmb_y; //!
TBranch *b_bpmc_x; //!
TBranch *b_bpmc_y; //!
TBranch *b_hseloss; //!
TBranch *b_hntracks; //!
TBranch *b_hcal_et; //!
TBranch *b_hgoodsc; //!
TBranch *b_hcal_e1; //!
TBranch *b_hcal_e2; //!
TBranch *b_hcal_e3; //!
TBranch *b_hcal_e4; //!
TBranch *b_hstrk_et; //!
TBranch *b_hstrk_in; //!
engine_class(TTree *tree=0);
virtual ~engine_class();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
virtual Long64_t LoadTree(Long64_t entry);
virtual void Init(TTree *tree);
virtual void Loop();
virtual Bool_t Notify();
virtual void Show(Long64_t entry = -1);
};
#endif
#ifdef engine_class_cxx
engine_class::engine_class(TTree *tree) : fChain(0)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
if (tree == 0) {
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("hms52949_rz.root");
if (!f || !f->IsOpen()) {
f = new TFile("hms52949_rz.root");
}
f->GetObject("h9010",tree);
}
Init(tree);
}
engine_class::~engine_class()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t engine_class::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t engine_class::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->GetTreeNumber() != fCurrent) {
fCurrent = fChain->GetTreeNumber();
Notify();
}
return centry;
}
void engine_class::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("hcer_npe", &hcer_npe, &b_hcer_npe);
fChain->SetBranchAddress("hsp", &hsp, &b_hsp);
fChain->SetBranchAddress("hse", &hse, &b_hse);
fChain->SetBranchAddress("charge", &charge, &b_charge);
fChain->SetBranchAddress("hsdelta", &hsdelta, &b_hsdelta);
fChain->SetBranchAddress("hstheta", &hstheta, &b_hstheta);
fChain->SetBranchAddress("hsphi", &hsphi, &b_hsphi);
fChain->SetBranchAddress("w", &w, &b_w);
fChain->SetBranchAddress("hszbeam", &hszbeam, &b_hszbeam);
fChain->SetBranchAddress("hsdedx1", &hsdedx1, &b_hsdedx1);
fChain->SetBranchAddress("hsbeta", &hsbeta, &b_hsbeta);
fChain->SetBranchAddress("hsshtrk", &hsshtrk, &b_hsshtrk);
fChain->SetBranchAddress("hsprtrk", &hsprtrk, &b_hsprtrk);
fChain->SetBranchAddress("hsxfp", &hsxfp, &b_hsxfp);
fChain->SetBranchAddress("hsyfp", &hsyfp, &b_hsyfp);
fChain->SetBranchAddress("hsxpfp", &hsxpfp, &b_hsxpfp);
fChain->SetBranchAddress("hsypfp", &hsypfp, &b_hsypfp);
fChain->SetBranchAddress("hsytar", &hsytar, &b_hsytar);
fChain->SetBranchAddress("hsxptar", &hsxptar, &b_hsxptar);
fChain->SetBranchAddress("hsyptar", &hsyptar, &b_hsyptar);
fChain->SetBranchAddress("hstart", &hstart, &b_hstart);
fChain->SetBranchAddress("eventid", &eventid, &b_eventid);
fChain->SetBranchAddress("ev_type", &ev_type, &b_ev_type);
fChain->SetBranchAddress("gfrx_raw", &gfrx_raw, &b_gfrx_raw);
fChain->SetBranchAddress("gfry_raw", &gfry_raw, &b_gfry_raw);
fChain->SetBranchAddress("gbeam_x", &gbeam_x, &b_gbeam_x);
fChain->SetBranchAddress("gbeam_y", &gbeam_y, &b_gbeam_y);
fChain->SetBranchAddress("bpma_x", &bpma_x, &b_bpma_x);
fChain->SetBranchAddress("bpma_y", &bpma_y, &b_bpma_y);
fChain->SetBranchAddress("bpmb_x", &bpmb_x, &b_bpmb_x);
fChain->SetBranchAddress("bpmb_y", &bpmb_y, &b_bpmb_y);
fChain->SetBranchAddress("bpmc_x", &bpmc_x, &b_bpmc_x);
fChain->SetBranchAddress("bpmc_y", &bpmc_y, &b_bpmc_y);
fChain->SetBranchAddress("hseloss", &hseloss, &b_hseloss);
fChain->SetBranchAddress("hntracks", &hntracks, &b_hntracks);
fChain->SetBranchAddress("hcal_et", &hcal_et, &b_hcal_et);
fChain->SetBranchAddress("hgoodsc", &hgoodsc, &b_hgoodsc);
fChain->SetBranchAddress("hcal_e1", &hcal_e1, &b_hcal_e1);
fChain->SetBranchAddress("hcal_e2", &hcal_e2, &b_hcal_e2);
fChain->SetBranchAddress("hcal_e3", &hcal_e3, &b_hcal_e3);
fChain->SetBranchAddress("hcal_e4", &hcal_e4, &b_hcal_e4);
fChain->SetBranchAddress("hstrk_et", &hstrk_et, &b_hstrk_et);
fChain->SetBranchAddress("hstrk_in", &hstrk_in, &b_hstrk_in);
Notify();
}
Bool_t engine_class::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
void engine_class::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
Int_t engine_class::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
#endif // #ifdef engine_class_cxx
#define hcana_class_cxx
#include "hcana_class.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void hcana_class::Loop()
{
// In a ROOT session, you can do:
// Root > .L hcana_class.C
// Root > hcana_class t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
}
}
This diff is collapsed.
../hms52949_rz.root
\ No newline at end of file
../hodtest_52949.root
\ No newline at end of file
#include <TTree.h>
#include <TFile.h>
#include <TROOT.h>
#include <iostream>
using namespace std;
void make_class(string root_file="hodtest_52949.root",
string tree_name="T;1",
string class_name="hcana_class") {
TFile* r_file = new TFile(root_file.c_str());
// TTree *T = (TTree*)gROOT->FindObject("T;1");
TTree *T = (TTree*)gROOT->FindObject(tree_name.c_str());
T->MakeClass(class_name.c_str());
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment