Portals
Pythia
Pythia Portal (prototype) at CACR Envoy:
https://envoy.cacr.caltech.edu:8443/clarens/web/pythia/pythia.html
Screen shot (click for large version)
ROOT
Example analysis output histogram from ROOT, QCD
background and 140Gev Higgs diphoton events, fitted to Gaussian and polynomial
background.
Code for analysis:
This function creates the chain of .root files for input, and calls the analysis loop
{
TBrowser b;
TChain* ch = new TChain("EGAMMA");
ch->Add("../eg05/eg05_Hgg_IVBfusion_m140_1.root");
//ch->Add("../eg05/eg05_Hgg_IVBfusion_m140_*.root");
//ch->Add("../QCD/eg05_jets_merged.root");
gROOT->ProcessLine(".L Analyse_Eg05.C");
Analyse_Eg05 t(ch);
t.Loop("QCD_histograms.root");
}
This contains the analysis loop function for all events in the chain.
#define Analyse_Eg05_cxx
#include "Analyse_Eg05.h"
#include
#include
#include
#define MY_PI 3.141592654
Int_t called = 0;
// Quadratic background function
Double_t background(Double_t *x, Double_t *par) {
return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}
// Lorenzian Peak function
Double_t lorentzianPeak(Double_t *x, Double_t *par) {
return (0.5*par[0]*par[1]/TMath::Pi()) /
TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
}
// Gaussian
Double_t gaussianPeak(Double_t *x, Double_t *par) {
Double_t A = par[0];
Double_t mu = par[1];
Double_t sigma = par[2];
A /= sigma;
A /= TMath::Sqrt(2.*MY_PI);
Double_t B = x[0]-mu;
Double_t Gauss = A*TMath::Exp(-B*B/(2*sigma*sigma));
called++;
return Gauss;
}
// Sum of background and peak function
Double_t fitFunction(Double_t *x, Double_t *par) {
return background(x,par) + gaussianPeak(x,&par[3]);
}
// Vladimir's energy corection function 1
Float_t Corfac_S25( Double_t eta)
{
// nore ES25 avec correction HoE,barrel
Float_t xcor_S25=1.;
if (fabs(eta)<0.25) { xcor_S25=1/0.963; } //0.970
else if (fabs(eta)<0.5) { xcor_S25=1/0.965; } //0.970
else if (fabs(eta)<0.75){ xcor_S25=1/0.961; } //0.969
else if (fabs(eta)<1.) { xcor_S25=1/0.960; }
else if (fabs(eta)<1.25){ xcor_S25=1/0.959; }
else if (fabs(eta)<1.5) { xcor_S25=1/0.957; }
else if (fabs(eta)<1.75){ xcor_S25=1/0.961; }
else if (fabs(eta)<2.) { xcor_S25=1/0.961; }
else if (fabs(eta)<2.25){ xcor_S25=1/0.962; }
else { xcor_S25=1/0.961; }
return xcor_S25;
}
// Vladimir's energy corection function 2
Float_t Corfac_Esc( Double_t eta)
{
// norme E super_clusters avec correction HoE,barrel
Float_t xcor_Esc=1.;
if (fabs(eta)<0.25) { xcor_Esc=1/0.971; }
else if (fabs(eta)<0.5) { xcor_Esc=1/0.972; }
else if (fabs(eta)<0.75){ xcor_Esc=1/0.970; }
else if (fabs(eta)<1.) { xcor_Esc=1/0.969; }
else if (fabs(eta)<1.25){ xcor_Esc=1/0.968; }
else if (fabs(eta)<1.5) { xcor_Esc=1/0.965; }
else if (fabs(eta)<1.75){ xcor_Esc=1/0.976; }
else if (fabs(eta)<2.) { xcor_Esc=1/0.980; }
else if (fabs(eta)<2.25){ xcor_Esc=1/0.981; }
else { xcor_Esc=1/0.979; }
return xcor_Esc;
}
void Analyse_Eg05::Loop(char histofile[])
{
if (fChain == 0) return;
// Photon four vectors
Double_t photon[2][4]; // 0 = Px, 1 = Py, 2 = Pz, 3 = E
Double_t correction[2][2]; // Energy corrections for each photon
// create array of histos
TObjArray hlist(0);
TH1F *hSCEnergies = new TH1F("hSCEnergies","SuperCluster Energy (GeV/c^2)",100,0.0,200.0);
TH1F *hSCET = new TH1F("hSCET","SuperCluster Et (GeV/c^2)",100,0.0,200.0);
TH1F *hSCEta = new TH1F("hSCEta","SuperCluster Eta",100,-5.0,5.0);
TH1F *hSCPhi = new TH1F("hSCPhi","SuperCluster Phi (Degrees)",100,0.0,2.0*MY_PI);
TH1F *hSCRadius = new TH1F("hSCRadius","SuperCluster Radius (cm?)",100,0.0,500.0);
TH1F *hSCTheta = new TH1F("hSCTheta","SuperCluster Theta (Degrees)",100,-2.0*MY_PI,2.0*MY_PI);
TH2F *hC9b25 = new TH2F("hC9b25","Seed Cluster S9 vs S25 Energy",100,1.0,1.5,100,1.0,1.5);
TH1F *hRecMass = new TH1F("hRecMass","Invariant Mass (Uncorrected) of Supercluster Pairs (GeV/c^2)",100,100.,250.);
TH1F *hRecMass1 = new TH1F("hRecMass1","Invariant Mass (Correction S25) of Supercluster Pairs (GeV/c^2)",100,100.,250.);
TH1F *hRecMass2 = new TH1F("hRecMass2","Invariant Mass (Correction Esc) of Supercluster Pairs (GeV/c^2)",100,100.,250.);
hlist.Add(hSCEnergies);
hlist.Add(hSCET);
hlist.Add(hSCEta);
hlist.Add(hSCPhi);
hlist.Add(hSCRadius);
hlist.Add(hSCTheta);
hlist.Add(hC9b25);
hlist.Add(hRecMass);
hlist.Add(hRecMass1);
hlist.Add(hRecMass2);
ofstream cfile;
cfile.open("Signal1.data", ios::out);
cfile << "ETtot,Eg1,Etg1,Eg2,Etg2,deta,dphi,MHiggs" << endl;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0, nclusters = 0;
for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
nclusters += negsc;
if(negsc < 2) continue; // discard events with less than two Egamma superclusters
cfile << L1glbl_L1TotEt << " ";
Float_t eta1;
Float_t phi1;
// Extract the measured parameters of each photon, and store as a fourvector
for(Int_t k=0;k<2;k++) {
Float_t E = egsc_scenergy[k];
Float_t ET = egsc_scet[k];
Float_t eta = egsc_sceta[k];
Float_t phi = egsc_scphi[k];
Float_t theta = egsc_sctheta[k];
Float_t radius = egsc_scradius[k];
cfile << E << " " << ET << " ";
if(k == 0) {
eta1 = eta;
phi1 = phi;
} else if (k == 1) {
cfile << fabs(eta1-eta) << " " << fabs(phi1-phi);
}
hSCEnergies->Fill(E);
hSCET->Fill(ET);
hSCEta->Fill(eta);
hSCPhi->Fill(phi);
hSCTheta->Fill(theta);
hSCRadius->Fill(radius);
correction[k][0] = Corfac_S25(eta); // this correction is the best
correction[k][1] = Corfac_Esc(eta);
Double_t pt = E*sin(theta); // transverse momentum
photon[k][0] = pt*cos(phi); // px
photon[k][1] = pt*sin(phi); // py
photon[k][2] = E*cos(theta); // pz
photon[k][3] = E; // E
// Get the Supercluster's seed cluster
Int_t seed = egsc_sciseed[k];
// Get the seed cluster's S9, S25
Float_t S9 = 1.0e-9 * egc_S9[seed];
Float_t S25 = 1.0e-9 * egc_S25[seed];
//cout << "seed " << seed << " S9 " << S9 << " S25 " << S25 << endl;
hC9b25->Fill(S9,S25);
}
// Calculate the invariant mass of the photon pair, using the fourvectors
Double_t P12[4]; //sum of 4-momenta
P12[0] = photon[0][0]+photon[1][0];
P12[1] = photon[0][1]+photon[1][1];
P12[2] = photon[0][2]+photon[1][2];
P12[3] = photon[0][3]+photon[1][3];
Double_t recmass = sqrt(P12[3]*P12[3]-P12[0]*P12[0]-P12[1]*P12[1]-P12[2]*P12[2]);
hRecMass->Fill(recmass);
// Calculate the invariant mass using the E25 correction for energy
P12[0] = correction[0][0]*photon[0][0] + correction[1][0]*photon[1][0];
P12[1] = correction[0][0]*photon[0][1] + correction[1][0]*photon[1][1];
P12[2] = correction[0][0]*photon[0][2] + correction[1][0]*photon[1][2];
P12[3] = correction[0][0]*photon[0][3] + correction[1][0]*photon[1][3];
Double_t recmass1 = sqrt(P12[3]*P12[3]-P12[0]*P12[0]-P12[1]*P12[1]-P12[2]*P12[2]);
hRecMass1->Fill(recmass1);
cfile << " " << recmass1 << endl;
// Calculate the invariant mass using the Esc correction for energy
P12[0] = correction[0][1]*photon[0][0] + correction[1][1]*photon[1][0];
P12[1] = correction[0][1]*photon[0][1] + correction[1][1]*photon[1][1];
P12[2] = correction[0][1]*photon[0][2] + correction[1][1]*photon[1][2];
P12[3] = correction[0][1]*photon[0][3] + correction[1][1]*photon[1][3];
Double_t recmass2 = sqrt(P12[3]*P12[3]-P12[0]*P12[0]-P12[1]*P12[1]-P12[2]*P12[2]);
hRecMass2->Fill(recmass2);
}
cout << "Total bytes " << nb << endl;
cout << "Total clusters " << nclusters << endl;
Double_t par[6];
// create a TF1 with the range from 0 to 3 and 6 parameters
TF1 *fitFcn = new TF1("fitFcn",fitFunction,100,250,6);
TF1 *backFcn = new TF1("backFcn",background,100,250,3);
TF1 *signalFcn = new TF1("signalFcn",gaussianPeak,100,250,3);
fitFcn->SetNpx(500); // Number of points used to draw the function
fitFcn->SetLineWidth(2);
fitFcn->SetLineColor(kMagenta);
backFcn->SetNpx(500);
backFcn->SetLineColor(kRed);
backFcn->SetLineWidth(2);
signalFcn->SetLineColor(kBlue);
signalFcn->SetNpx(500);
fitFcn->SetParLimits(3,100.0,10000.0); // Limits on Gauss A
fitFcn->SetParLimits(4,120.0,160.0); // Limits on Gauss mu
fitFcn->SetParLimits(5,1.0,20.); // Limits on Gauss sigma
fitFcn->SetParameters(1,1,1,1000.0,140.0,5.0);
// V+ means produce some output from Minuit, and add the fit result to the histogram
// "ep" are drawing options - error bars with point
//hRecMass1->Fit("fitFcn","V+","ep");
hRecMass1->Fit("fitFcn","V+","ep");
// writes the fit results into the par array
fitFcn->GetParameters(par);
//fitFcn->Draw("same");
backFcn->SetParameters(par);
backFcn->Draw("same");
//signalFcn->SetParameters(&par[3]);
//signalFcn->Draw("same");
TString sDetails = "Fit: M=" + TString(Form("%5.1f s=%5.2f GeV/c^2",par[4],par[5]));
// draw the legend
TLegend *legend=new TLegend(0.6,0.8,1.0,0.6);
legend->SetTextFont(72);
legend->SetTextSize(0.03);
legend->AddEntry(hRecMass,"Data","lpe");
legend->AddEntry(backFcn,"Background","l");
//legend->AddEntry(signalFcn,"Signal Fit","l");
legend->AddEntry(fitFcn,sDetails,"l");
legend->Draw();
TString histos(histofile);
if(histos != "") {
TFile fout(histofile,"recreate");
hlist.Write();
fout.Close();
}
cfile.flush();
cfile.close();
}
