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:

{
  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");
}

#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();


}

  02/16/2007 by Julian Bunn, email: Julian.Bunn@caltech.edu