// ESAF : Euso Simulation and Analysis Framework
// $Id: EusoHistoFactory.cc,v 1.10 2005/02/12 18:17:40 thea Exp $
// created Apr, 18 2003
#include "Riostream.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TText.h"
#include "TPad.h"
#include "TStyle.h"
#include "EusoHistoFactory.hh"
#include "EEvent.hh"
#include "EDetPhoton.hh"
#include "EFee.hh"
#include "EAFee.hh"
#include "Config.hh"
ClassImp(EusoHistoFactory)
EusoHistoFactory *EusoHistoFactory::fMe = NULL;
// constructor
EusoHistoFactory::EusoHistoFactory() {
SetCanvas( NULL );
}
// destructor
EusoHistoFactory::~EusoHistoFactory() {
Close();
}
// close canvas if owner
void EusoHistoFactory::Close() {
if ( Canvas() )
delete Canvas();
SetCanvas( NULL );
}
// instance
EusoHistoFactory* EusoHistoFactory::Get() {
if ( fMe == NULL )
fMe = new EusoHistoFactory();
return fMe;
}
// standard interface to event histograms
void EusoHistoFactory::DoHisto(EusoStandardHistoId h, EEvent* pEv, TCanvas*c) {
if (!pEv) return;
TCanvas *canvas;
if ( c ) {
canvas = c;
SetCanvas( NULL );
}
else {
canvas = new TCanvas( "EusoHisto", "Euso Histogram" );
SetCanvas( canvas );
}
switch ( h ) {
case HS_GTUTIME_1: DoGtuPhotonsHisto(pEv,canvas); break;
case HS_PHOTON_PHICOSTH_1: DoPhiCosThetaHisto(pEv,canvas); break;
case HS_PHOTON_XY_1: DoXYHisto(pEv, canvas); break;
case HS_PHOTON_WLXY_1: DoWavelengthXYHisto(pEv, canvas); break;
default: break;
}
}
// global histogram on a tree
void EusoHistoFactory::DoGlobalHisto(EusoStandardHistoId h,TTree* tree,TCanvas* c) {
if (!tree) return;
TCanvas *canvas;
if ( c ) {
canvas = c;
SetCanvas( NULL );
}
else {
canvas = new TCanvas( "EusoGlobalHisto", "Euso Global Histogram" );
SetCanvas( canvas );
}
switch ( h ) {
case HG_TRIGGER_EFF_1: DoTriggerEffic1(tree,canvas); break;
default: break;
}
}
// Time distribution of photons at various stage of Euso detector
void EusoHistoFactory::DoGtuPhotonsHisto(EEvent* pEv,TCanvas *canvas) {
cout << "DoGtuPhotonsOnHito called" << endl;
EDetector* pDet;
if ( !(pDet = pEv->GetDetector()) ) {
cout << "No EDetector in this event, exiting." << endl;
return;
}
canvas->Clear();
canvas->cd();
TPad *pad1 = new TPad("Pad1","Pad1",0.05,0.50,0.45,0.95,21);
TPad *pad2 = new TPad("Pad2","Pad2",0.50,0.50,0.95,0.95,21);
TPad *pad3 = new TPad("Pad3","Pad3",0.05,0.05,0.45,0.45,21);
TPad *pad4 = new TPad("Pad4","Pad4",0.50,0.05,0.95,0.45,21);
pad1->Draw();
pad2->Draw();
pad3->Draw();
pad4->Draw();
double GtuLength = 2500.; //Config::Get()->GetCF("Electronics","Macrocell")->GetNum("GtuTimeLength");
Int_t Gtus = (Int_t) ((pDet->GetLastTime() - pDet->GetFirstTime() ) / GtuLength);
gROOT->Delete( "EusoHF_GtuPho*" );
TH1F *th0 = new TH1F( "EusoHF_GtuPho0", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th1 = new TH1F( "EusoHF_GtuPho1", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th2 = new TH1F( "EusoHF_GtuPho2", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
for( Int_t i=0; i < pDet->GetNumPhotons(); i++ ) {
EDetPhoton& ph = *(pDet->GetPhoton( i ));
Float_t time = ( ph.GetTime() - pDet->GetFirstTime() ) / GtuLength;
th0->Fill( time );
if ( ph.CrossedIFS() ) th1->Fill( time );
if ( ph.CrossedIFS() && ph.MadeSignal() ) th2->Fill( time );
}
pad1->cd();
// draw histograms with colors and filling
th0->SetStats(0);
th0->SetLineColor(9);
th0->SetFillColor(0);
th0->SetXTitle("GTU number");
th0->SetYTitle("Number of Photons");
th0->Draw();
th1->SetStats(0);
th1->SetLineColor(2);
th1->SetFillColor(0);
th1->Draw("same");
th2->SetStats(0);
th2->SetLineColor(8);
th2->SetFillColor(8);
th2->Draw("same");
// write meaning and statistics
TText *t = new TText();
char txt[200];
t->SetTextFont(62);
t->SetTextColor(9);
t->SetTextSize(0.04);
sprintf(txt,"Ph. on Pupil:%d",(int)th0->GetEntries());
t->DrawText(60.,th0->GetMaximum()*0.94,txt);
t->SetTextColor(2);
sprintf(txt,"Ph. on FS:%d",(int)th1->GetEntries());
t->DrawText(60.,th0->GetMaximum()*0.86,txt);
t->SetTextColor(8);
sprintf(txt,"Ph. detected:%d",(int)th2->GetEntries());
t->DrawText(60.,th0->GetMaximum()*0.78,txt);
canvas->Update();
// second pad
pad2->cd();
TH1F *th3 = new TH1F( "EusoHF_GtuPho3", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th4 = new TH1F( "EusoHF_GtuPho4", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th5 = new TH1F( "EusoHF_GtuPho5", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th6 = new TH1F( "EusoHF_GtuPho6", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
for( Int_t i=0; i < pDet->GetNumPhotons(); i++ ) {
EDetPhoton& ph = *(pDet->GetPhoton( i ));
Float_t time = ( ph.GetTime() - pDet->GetFirstTime() ) / GtuLength;
if ( ph.CrossedIFS() ) th3->Fill( time );
if ( ph.CrossedIFS() && ph.MadeCount() ) th5->Fill( time );
if ( ph.CrossedIFS() && ph.MadeFastOR() ) th6->Fill( time );
}
th3->SetStats(0);
th3->SetLineColor(9);
th3->SetFillColor(0);
th3->SetXTitle("GTU number");
th3->SetYTitle("Number of Photons");
th3->Draw();
th4->SetStats(0);
th4->SetLineColor(6);
th4->SetFillColor(0);
th4->Draw("same");
th5->SetStats(0);
th5->SetLineColor(2);
th5->SetFillColor(2);
th5->Draw("same");
th6->SetStats(0);
th6->SetLineColor(8);
th6->SetFillColor(8);
th6->Draw("same");
TText *t2 = new TText();
t2->SetTextFont(62);
t2->SetTextColor(9);
t2->SetTextSize(0.04);
sprintf(txt,"Photons on FS:%d",(int)th3->GetEntries());
t2->DrawText(65.,th3->GetMaximum()*0.94,txt);
t2->SetTextColor(6);
sprintf(txt,"Pmt signals:%d",(int)th4->GetEntries());
t2->DrawText(65.,th3->GetMaximum()*0.86,txt);
t2->SetTextColor(2);
sprintf(txt,"Frontend hits:%d",(int)th5->GetEntries());
t2->DrawText(65.,th3->GetMaximum()*0.78,txt);
t2->SetTextColor(8);
sprintf(txt,"Fast OR counts:%d",(int)th6->GetEntries());
t2->DrawText(65.,th3->GetMaximum()*0.70,txt);
canvas->Update();
// third pad: charge
pad3->cd();
TH1F *th7 = new TH1F( "EusoHF_GtuPho7", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th8 = new TH1F( "EusoHF_GtuPho8", "Photons vs GTU", Gtus, 0., (Float_t) Gtus);
TH1F *th9 = new TH1F( "EusoHF_GtuPho9", "Charge vs GTU", Gtus, 0., (Float_t) Gtus);
for( Int_t i=0; i < pDet->GetNumPhotons(); i++ ) {
EDetPhoton& ph = *(pDet->GetPhoton( i ));
Float_t time = ( ph.GetTime() - pDet->GetFirstTime() ) / GtuLength;
if ( ph.CrossedIFS() && ph.MadeCount() ) th7->Fill( time );
if ( ph.CrossedIFS() && ph.MadeFastOR() ) th8->Fill( time );
}
Float_t total_afee_charge=0.;
for(Int_t i=0; i < pDet->GetNumAFee(); i++ ) {
EAFee *afee = pDet->GetAFee( i );
Int_t gtu = afee->GetGtu();
Float_t charge = afee->GetDyCharge();
Int_t numpe = (Int_t)(charge*2.e13); // convert into NOMINAL photoelectrons, gain 3.E5
if ( gtu > 0 && gtu < Gtus) {
Int_t old = (Int_t)(th3->GetBinContent( gtu ));
th9->SetBinContent( gtu, old + numpe );
total_afee_charge += numpe;
}
}
th7->SetStats(0);
th7->SetLineColor(9);
th7->SetFillColor(0);
th7->SetXTitle("GTU number");
th7->SetYTitle("Number of Photons");
th7->Draw();
th8->SetStats(0);
th8->SetLineColor(6);
th8->SetFillColor(0);
th8->Draw("same");
th9->SetStats(0);
th9->SetLineColor(8);
th9->SetFillColor(8);
th9->Draw("same");
TText *t3 = new TText();
t3->SetTextFont(62);
t3->SetTextColor(9);
t3->SetTextSize(0.04);
sprintf(txt,"Ph. detected:%d",(Int_t)th7->GetEntries());
t3->DrawText(65.,th7->GetMaximum()*0.94,txt);
t3->SetTextColor(6);
sprintf(txt,"Ph. in DFEE:%d",(Int_t)th8->GetEntries());
t3->DrawText(65.,th7->GetMaximum()*0.86,txt);
t3->SetTextColor(2);
sprintf(txt,"Charge (pC,dynode):%f",total_afee_charge);
t3->DrawText(65.,th7->GetMaximum()*0.78,txt);
canvas->Update();
}
// photons on focal surface
void EusoHistoFactory::DoPhiCosThetaHisto(EEvent* ev, TCanvas *canvas){
cout << "DoPhiCosThetaHisto called" << endl;
EDetector* pDet;
if ( !(pDet = ev->GetDetector()) ) {
cout << "No EDetector in this event, exiting." << endl;
return;
}
canvas->Clear();
canvas->Divide(2,2);
Float_t th0CosThMin(10),th0CosThMax(-10);
Float_t th0PhiMin(10),th0PhiMax(-10);
for ( Int_t i=0; i < pDet->GetNumPhotons(); i++){
EDetPhoton& ph=*(pDet->GetPhoton(i));
if (TMath::Cos(ph.GetTheta()) > th0CosThMax)
th0CosThMax = TMath::Cos(ph.GetTheta());
if (ph.GetPhi() > th0PhiMax)
th0PhiMax = ph.GetPhi();
if (TMath::Cos(ph.GetTheta()) < th0CosThMin)
th0CosThMin = TMath::Cos(ph.GetTheta());
if (ph.GetPhi() < th0PhiMin)
th0PhiMin = ph.GetPhi();
}
Float_t th0Delta = th0CosThMax-th0CosThMin;
th0CosThMin-=0.1*th0Delta;
th0CosThMax+=0.1*th0Delta;
th0Delta = th0PhiMax-th0PhiMin;
th0PhiMin-=0.1*th0Delta;
th0PhiMax+=0.1*th0Delta;
gDirectory->Delete("EusoHF_PhiCosTh*");
TH2F *th0 = new TH2F( "EusoHF_PhiCosTh0", "Cos(Theta) vs Phi", 100,
th0CosThMin, th0CosThMax, 100, th0PhiMin, th0PhiMax);
TH2F *th1 = new TH2F( "EusoHF_PhiCosTh1", "Cos(Theta) vs Phi", 100,
th0CosThMin, th0CosThMax, 100, th0PhiMin, th0PhiMax);
TH2F *th2 = new TH2F( "EusoHF_PhiCosTh2", "Cos(Theta) vs Phi", 100,
th0CosThMin, th0CosThMax, 100, th0PhiMin, th0PhiMax);
TH2F *th3 = new TH2F( "EusoHF_PhiCosTh3", "Cos(Theta) vs Phi", 100,
th0CosThMin, th0CosThMax, 100, th0PhiMin, th0PhiMax);
for ( Int_t i=0; i < pDet->GetNumPhotons(); i++){
EDetPhoton& ph = *(pDet->GetPhoton( i ));
th0->Fill(TMath::Cos(ph.GetTheta()),ph.GetPhi());
if (ph.MadeSignal() && ph.CrossedIFS()) th1->Fill(TMath::Cos(ph.GetTheta()),ph.GetPhi());
if (ph.CrossedIFS()) th2->Fill(TMath::Cos(ph.GetTheta()),ph.GetPhi());
if (ph.MadeFastOR()) th3->Fill(TMath::Cos(ph.GetTheta()),ph.GetPhi());
}
th0->SetStats(0);
th1->SetStats(0);
th2->SetStats(0);
th3->SetStats(0);
canvas->cd(1);
th0->Draw("colz");
canvas->cd(3);
th1->Draw("colz");
canvas->cd(2);
th2->Draw("colz");
canvas->cd(4);
th3->Draw("colz");
canvas->Update();
}
//photons on XY plane
void EusoHistoFactory::DoXYHisto(EEvent* ev, TCanvas *c){
cout << "DoXY called" << endl;
EDetector* pDet;
if ( !(pDet = ev->GetDetector()) ) {
cout << "No EDetector in this event, exiting." << endl;
return;
}
c->Clear();
c->Divide(2,2);
//Get the number of photons
Int_t numph=pDet->GetNumPhotons();
cout << "Total number of photons: " << numph << endl;
gDirectory->Delete("EusoHF_XY*");
Float_t x[numph];
Float_t y[numph];
Float_t xsh[numph];
Float_t ysh[numph];
Float_t xshmax=0, xshmin=0, yshmax=0, yshmin=0;
Int_t nlist[numph];
//read coordinates x,y for each photon
for (Int_t i=0; i<numph; i++) {
EDetPhoton& ph=*(pDet->GetPhoton(i));
x[i]=ph.GetPos()[0];
y[i]=ph.GetPos()[1];
}
Int_t j=0;
//roughly determinate which photons belong to the shower
for (Int_t i=0; i<numph; i++){
Float_t deltax=TMath::Abs(x[i]-x[i+1]);
if (deltax < 1.5) {
Float_t deltay=TMath::Abs(y[i]-y[i+1]);
if (deltay < 1.5) {
xsh[j]=x[i];
ysh[j]=y[i];
xsh[j+1]=x[i+1];
ysh[j+1]=y[i+1];
j++;
}
}
}
j++;
//calculate x,y of the shower minimum and maximum
for (Int_t i=0; i<j; i++) {
if (i==0) {
xshmin=xsh[0];
xshmax=xsh[0];
yshmin=ysh[0];
yshmax=ysh[0];
}
if (i>0) {
if (xsh[i]<xshmin) {xshmin=xsh[i];}
if (xsh[i]>xshmax) {xshmax=xsh[i];}
if (ysh[i]<yshmin) {yshmin=ysh[i];}
if (ysh[i]>yshmax) {yshmax=ysh[i];}
}
}
Float_t deltashx = xshmax - xshmin;
Float_t deltashy = yshmax - yshmin;
xshmin -= deltashx / 10.0;
xshmax += deltashx / 10.0;
yshmin -= deltashy / 10.0;
yshmax += deltashy / 10.0;
cout << "Creating istograms..." << endl;
cout << "X-Min: " << xshmin << endl;
cout << "X-Max: " << xshmax << endl;
cout << "Y-Min: " << yshmin << endl;
cout << "Y-Max: " << yshmax << endl;
//Select photons in the interval (xshmin,xshmax);(yshmin,yshmax)
Int_t l=0;
for (Int_t i=0; i<numph; i++) {
if ((((x[i]>xshmin && x[i]<xshmax) && y[i]>yshmin) && y[i]<yshmax)) {
nlist[l]=i;
l++;
}
}
l++;
//Create istograms and fill them
TH2F *th1=new TH2F("EusoHF_XY1","Hits on IFS",1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th2=new TH2F("EusoHF_XY2","MadeSignal",1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th3=new TH2F("EusoHF_XY3","MadeCount",1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th4=new TH2F("EusoHF_XY4","MadeFastOR",1000,xshmin,xshmax,1000,yshmin,yshmax);
for (Int_t i=0; i<l; i++) {
Int_t k=nlist[i];
EDetPhoton& ph=*(pDet->GetPhoton(k));
if (ph.CrossedIFS()) th1->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.MadeSignal()) th2->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.MadeCount()) th3->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.MadeFastOR()) th4->Fill(ph.GetPos()[0],ph.GetPos()[1]);
}
gStyle->SetOptStat(10);
th1->SetStats(kTRUE);
th1->SetXTitle("x");
th1->SetYTitle("y");
th2->SetStats(kTRUE);
th2->SetXTitle("x");
th2->SetYTitle("y");
th3->SetStats(kTRUE);
th3->SetXTitle("x");
th3->SetYTitle("y");
th4->SetStats(kTRUE);
th4->SetXTitle("x");
th4->SetYTitle("y");
c->cd(1);
th1->Draw();
c->cd(2);
th2->Draw();
c->cd(3);
th3->Draw();
c->cd(4);
th4->Draw();
c->Update();
cout << "DoXY ended" << endl;
}
//photons on the IFS in XY plane at various wavelenghts
void EusoHistoFactory::DoWavelengthXYHisto(EEvent* ev, TCanvas *c){
cout << "DoWaveLength called" << endl;
EDetector* pDet;
if ( !(pDet = ev->GetDetector()) ) {
cout << "No EDetector in this event, exiting." << endl;
return;
}
c->Clear();
c->Divide(2,2);
//Get the number of photons
Int_t numph=pDet->GetNumPhotons();
cout << "Total number of photons: " << numph << endl;
gDirectory->Delete("EusoHF_WLXY*");
Float_t x[numph];
Float_t y[numph];
Float_t xsh[numph];
Float_t ysh[numph];
Float_t xshmax=0, xshmin=0, yshmax=0, yshmin=0;
Int_t nlist[numph];
//read coordinates x,y for each photon
for (Int_t i=0; i<numph; i++) {
EDetPhoton& ph=*(pDet->GetPhoton(i));
x[i]=ph.GetPos()[0];
y[i]=ph.GetPos()[1];
}
Int_t j=0;
//determinate which photons belong to the shower
for (Int_t i=0; i<numph; i++){
Float_t deltax=TMath::Abs(x[i]-x[i+1]);
if (deltax < 1.5) {
Float_t deltay=TMath::Abs(y[i]-y[i+1]);
if (deltay < 1.5) {
xsh[j]=x[i];
ysh[j]=y[i];
xsh[j+1]=x[i+1];
ysh[j+1]=y[i+1];
j++;
}
}
}
j++;
//calculate x,y of the shower minimum and maximum
for (Int_t i=0; i<j; i++) {
if (i==0) {
xshmin=xsh[0];
xshmax=xsh[0];
yshmin=ysh[0];
yshmax=ysh[0];
}
if (i>0) {
if (xsh[i]<xshmin) {xshmin=xsh[i];}
if (xsh[i]>xshmax) {xshmax=xsh[i];}
if (ysh[i]<yshmin) {yshmin=ysh[i];}
if (ysh[i]>yshmax) {yshmax=ysh[i];}
}
}
Float_t deltashx = xshmax - xshmin;
Float_t deltashy = yshmax - yshmin;
xshmin -= deltashx / 10.0;
xshmax += deltashx / 10.0;
yshmin -= deltashy / 10.0;
yshmax += deltashy / 10.0;
cout << "Creating istograms..." << endl;
cout << "X-Min: " << xshmin << endl;
cout << "X-Max: " << xshmax << endl;
cout << "Y-Min: " << yshmin << endl;
cout << "Y-Max: " << yshmax << endl;
//Select photons in the interval (xshmin,xshmax);(yshmin,yshmax)
Int_t l=0;
for (Int_t i=0; i<numph; i++) {
if ((((x[i]>xshmin && x[i]<xshmax) && y[i]>yshmin) && y[i]<yshmax)) {
nlist[l]=i;
l++;
}
}
l++;
//Create istograms and fill them
TH2F *th1=new TH2F("EusoHF_WLXY1","Hits on IFS, all wavelengths", 1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th2=new TH2F("EusoHF_WLXY2","Hits on IFS, #lambda < 350 nm", 1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th3=new TH2F("EusoHF_WLXY3","Hits on IFS, 350 nm #leq #lambda #leq 370 nm", 1000,xshmin,xshmax,1000,yshmin,yshmax);
TH2F *th4=new TH2F("EusoHF_WLXY4","Hits on IFS, #lambda > 370 nm", 1000,xshmin,xshmax,1000,yshmin,yshmax);
for (Int_t i=0; i<l; i++) {
Int_t k=nlist[i];
EDetPhoton& ph=*(pDet->GetPhoton(k));
if (ph.CrossedIFS()) th1->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.CrossedIFS() && ph.GetLambda()<0.00035) th2->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.CrossedIFS() && ph.GetLambda()>=0.00035 && ph.GetLambda()<=0.00037) th3->Fill(ph.GetPos()[0],ph.GetPos()[1]);
if (ph.CrossedIFS() && ph.GetLambda()>0.00037) th4->Fill(ph.GetPos()[0],ph.GetPos()[1]);
}
gStyle->SetOptStat(10);
th1->SetStats(kTRUE);
th1->SetXTitle("x");
th1->SetYTitle("y");
th2->SetStats(kTRUE);
th2->SetXTitle("x");
th2->SetYTitle("y");
th3->SetStats(kTRUE);
th3->SetXTitle("x");
th3->SetYTitle("y");
th4->SetStats(kTRUE);
th4->SetXTitle("x");
th4->SetYTitle("y");
c->cd(1);
th1->Draw();
c->cd(2);
th2->Draw();
c->cd(3);
th3->Draw();
c->cd(4);
th4->Draw();
c->Update();
cout << "DoWaveLenght ended" << endl;
}
// global histograms
void EusoHistoFactory::DoTriggerEffic1(TTree* tree, TCanvas* canvas) {
canvas->cd();
}