Hi,
I have a problem saving histograms to a file.
operating details :
root Version 3.01/06
cc5
solaris : SunOS 5.9
but I am sure this is just a case of not understanding C++ scope rules
since I can view histograms with a TBrowser during the session but
not in a new root session.
essentially I have in my script :
===============================
//To run : Root> .L TestRoot.cxx
// Root> TestRoot()
//OR
// Root> .x TestRoot.cxx
//
// Root> .L TestRoot.cxx+
// Root> TestRoot()
//
//OR...all at once
//
// Root> .x TestRoot.cxx+ //excutes function with same name
some includes
#ifndef __CINT__
some more includes
//______________________________________________________________________
int main()
{
int n;
//TROOT test("test","Test of histogramming and I/O");
n = TestRoot();
printf("\n %d",n);
return n;
}
#endif
/***************************************************/
/***************************************************/
int TestRoot() {
Data structure, variable declarations.
/* Root File */
TFile *hfile = new TFile("test.root","RECREATE","test");
hfile->mkdir("Raw_Ge");
/* Histograms */
Char_t name[50];
Char_t label[50];
TH1S *hStatus = new TH1S("Status", "Event Stat Code", 150, -0.5, 149.5);
TH1S *hEvLen = new TH1S("EvLen", "Event Length", 150, -0.5, 149.5);
more histos :
/* Raw Ge directory */
hfile->cd("Raw_Ge");
TH1S *hRaw_Ge_E[21];
for (int i = 0; i < 21; i++) {
sprintf(name, "E_%d",i+1);
sprintf(label,"Ge%d Energy",i+1);
hRaw_Ge_E[i] = new TH1S(name,label, 8192, -0.5, 8191.5);
}
/**********************/
/* Parse Data on disc */
/**********************/
while(fread(Buf,16384,1,i_f)) {
do some stuff......
}
fclose(i_f);
hfile->Write();
//hfile->Close();
return nBlocks;
}
//---------------------------------------
Included is a full version of the script. It seems to run
as a macro, within ACLIC and runs when compiled separately.
So, what have I declared wrongly that does not enable me to save the
histograms in the root file ?
Cheers,
Karl
==========================================================================
CEA Saclay, DAPNIA/SPhN Phone : (33) 01 69 08 7553
Bat 703 - l'Orme des Merisiers Fax : (33) 01 69 08 7584
F-91191 Gif-sur-Yvette E-mail : khauschild@cea.fr
France karl_hauschild@yahoo.co.uk
WWW: http://www-dapnia.cea.fr/Sphn
//To run : Root> .L TestRoot.cxx
// Root> TestRoot()
//OR
// Root> .x TestRoot.cxx
//
// Root> .L TestRoot.cxx+
// Root> TestRoot()
//
//OR...all at once
//
// Root> .x TestRoot.cxx+ //excutes function with same name
#include <stdlib.h>
#include <stdio.h>
int TestRoot();
/* includes if not cint */
#ifndef __CINT__
#include "TROOT.h"
//#include "Class.h"
//#include "Method.h"
#include "TFile.h"
//#include "TTree.h"
#include "TRandom.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TStopwatch.h"
//#define TRUE 1
//#define FALSE 0
//______________________________________________________________________
int main()
{
int n;
//TROOT test("test","Test of histogramming and I/O");
n = TestRoot();
printf("\n %d",n);
return n;
}
#endif
/***************************************************/
/***************************************************/
int TestRoot() {
/* Data Structure */
struct Data_t {
unsigned short int id;
unsigned short int tac;
unsigned short int chan;
};
/* Used to store detector data */
Data_t Ge[21], *p_Ge;
Data_t BGO[50], *p_BGO;
Data_t CsI[50], *p_CsI;
/* Re-map CsI id's */
int new_id[40] = { /* old id */
0, 0, 0, 0, 1, /* 0, 1, 2, 3, 4 */
2, 3, 4, 0, 5, /* 5, 6, 7, 8, 9 */
6, 7, 8, 0, 9, /*10,11,12,13,14 */
10,11, 0, 0,12, /*15,16,17,18,19 */
13,14,15,16,17, /*20,21,22,23,24 */
18,19,20,21,22, /*25,26,27,28,29 */
23,24, 0, 0, 0, /*30,31,32,33,34 */
0, 25,26,27,28}; /*35,36,37,38,39 */
/* Root File */
TFile *hfile = new TFile("test.root","RECREATE","test");
hfile->mkdir("Raw_Ge");
/* Histograms */
Char_t name[50];
Char_t label[50];
TH1S *hStatus = new TH1S("Status", "Event Stat Code", 150, -0.5, 149.5);
TH1S *hEvLen = new TH1S("EvLen", "Event Length", 150, -0.5, 149.5);
TH2S *hHK = new TH2S("HK", "BGO HK", 2500,-0.5,2499.5,50,-0.5,49.5);
TH1S *hnGe = new TH1S("nGe", "Raw Ge Mult", 21, -0.5, 20.5);
TH1S *hnBGO = new TH1S("nBGO", "Raw BGO Mult", 21, -0.5, 20.5);
TH1S *hnCsI = new TH1S("nCsI", "MiniBall Mult", 51, -0.5, 50.5);
TH1S *hncGe = new TH1S("ncGe", "Clean Ge Mult", 21, -0.5, 20.5);
TH1S *hncBGO = new TH1S("ncBGO", "Clean BGO Mult", 21, -0.5, 20.5);
TH1S *hncCsI = new TH1S("ncCsI", "Clean CsI Mult", 51, -0.5, 50.5);
TH1S *hHitPat_Ge = new TH1S("HitPat_Ge", "Hit Pattern", 21, -0.5, 20.5);
TH1S *hHitPat_BGO = new TH1S("HitPat_BGO", "Hit Pattern", 100, -0.5, 99.5);
/* Raw Ge directory */
hfile->cd("Raw_Ge");
TH1S *hRaw_Ge_E[21];
TH1S *hRaw_Ge_T[21];
for (int i = 0; i < 21; i++) {
sprintf(name, "E_%d",i+1);
sprintf(label,"Ge%d Energy",i+1);
hRaw_Ge_E[i] = new TH1S(name,label, 8192, -0.5, 8191.5);
sprintf(name, "T_%d",i+1);
sprintf(label,"Ge%d Time",i+1);
hRaw_Ge_T[i] = new TH1S(name,label, 2048, -0.5, 2047.5);
}
/* Ge-totals */
hfile->cd();
TH1S *hRaw_Ge_E_tot = new TH1S("Raw_Ge_E_tot","All Ge's",
8192, -0.5, 8191.5);
TH1S *hRaw_Ge_T_tot = new TH1S("Raw_Ge_T_tot","All Ge's",
2048, -0.5, 2047.5);
/* Rate monitor */
TStopwatch timer;
timer.Start();
double told = 0;
double tnew = 0;
int printev = 1000;
/* Data File */
FILE *i_f;
i_f=fopen("run2", "rb");
// Buffer / pointers
Char_t Buf[16384];
Char_t *pBuf;
Char_t tmp;
UShort_t *current_word;
UShort_t *cword; //current word
UShort_t *end_of_Buf;
UShort_t *end_of_data;
UShort_t *end_of_block;
UShort_t *start_of_event;
//Data varaibles
UShort_t H, K, KE2, ncGe, T_Isomer;
UShort_t nBGO, nCsI;
//Random Number generator
gRandom = new TRandom(27071970);
float random, E, fy;
// variables
int nBlocks = 0;
int n, i, j, t, id, chan, energy, tac, temp;
int good_event = 0;
int bad_event = 0;
int nevents = 0;
/**********************/
/* Parse Data on disc */
/**********************/
while(fread(Buf,16384,1,i_f)) {
nBlocks++;
if (nBlocks > 1)
return nBlocks;
// Reset pointer at start of block
pBuf = &Buf[0];
current_word = (UShort_t *)Buf;
end_of_block = current_word + 8191; /* 8192 words -> 16384 chars */
/* Find the last good data word */
end_of_data = end_of_block;
while( *end_of_data == 0xFFFF)
end_of_data--;
/*********************************************/
/* Process the buffer to last good data word */
/*********************************************/
for (nevents = 0; current_word < end_of_data; ) {
//printf("\ngot here");
hStatus->Fill(100);
nevents++;
start_of_event = current_word;
/* Check event length */
for (i = 0; current_word <= end_of_data; i++)
if (*(current_word+i) == 0xFFFF) break;
hEvLen->Fill(i);
/* minimum length is 10 words */
if ( i < 10) { /* run 1,3 */
/* if ( i < 13) { /* run 2 */
hStatus->Fill(101);
current_word = start_of_event;
goto Bad_Event;
}
/* extract header data */
H = *current_word++; /* BGO ball sum energy */
K = *current_word++; /* BGO ball multiplicity */
KE2 = *current_word++; /* "nothing important" */
ncGe = *current_word++; /* clean HPGe multiplicity */
hnGe->Fill(ncGe);
//check above header
if((H == 0xFFFF)||(K == 0xFFFF) ||(KE2 == 0xFFFF)||(ncGe == 0xFFFF) ||
((H == 0)&&(K == 0))||(ncGe == 0) ) {
printf("\n %d",nevents);
printf("\n%4x %4x %4x %4x",H,K,KE2,ncGe);
hStatus->Fill(103);
current_word = start_of_event;
goto Bad_Event;
}
hHK->Fill(H,K);
/* extract Ge Data */
for (temp = i = 0; i < ncGe; i++) {
id = *current_word++; /* HPGe ID */
energy = *current_word++; /* HPGe Energy */
tac = *current_word++; /* HPGe Time */
if ((id >= 0) && (id <= 20) &&
(tac > 0) && (energy > 0) ) {
hHitPat_Ge->Fill(id);
hRaw_Ge_E[id]->Fill(energy);
hRaw_Ge_T[id]->Fill(tac);
hRaw_Ge_E_tot->Fill(energy);
hRaw_Ge_T_tot->Fill(tac);
Ge[temp].id = id;
Ge[temp].chan = energy;
Ge[temp].tac = tac;
temp++;
} else {
// probably dropped the Ge id word
printf("\n block %d event %d\n",nBlocks, nevents);
for(j = -5; j < 20; j++) {
printf("%6x",*(start_of_event+j));
}
current_word = start_of_event;
goto Bad_Event;
}
}
/* re-assign Ge mult */
ncGe = temp;
hncGe->Fill(ncGe);
T_Isomer = *current_word++; /* "nothing important" */
/* extract BGO data */
nBGO = *current_word++; /* BGO ball multiplicity in time gate */
hnBGO->Fill(nBGO);
if (nBGO > 100) { /* screwed up somewhere */
printf("\n nBgo %d",nBGO);
current_word = start_of_event;
hStatus->Fill(105);
goto Bad_Event;
}
/* extract nBGO */
for (temp = i = 0; i < nBGO; i++) {
id = *current_word++; /* BGO ID */
chan = *current_word++; /* BGO Energy */
tac = *current_word++; /* BGO Time */
hHitPat_BGO->Fill(id);
if ((id >= 0) && (id <= 20) ) {
BGO[temp].id = id;
BGO[temp].chan = chan;
BGO[temp].tac = tac;
temp++;
} else {
/* probably dropped the Ge id word */
printf("\n Bad BGO Id %d : event %d\n",BGO[i].id,nevents);
for(j = -1; j < 20; j++) {
printf("%6x",*(start_of_event+j));
}
current_word = start_of_event;
hStatus->Fill(106);
goto Bad_Event;
}
}
nBGO = temp;
hncBGO->Fill(nBGO);
/* extract particle data */
nCsI = *current_word++; /* MiniBall multiplicity */
if (nCsI == 0xFFFF) /* did not write CsI data */
goto Bad_Event;
hnCsI->Fill(nCsI);
for (i = 0; i < nCsI; i++) {
CsI[i].id = *current_word++; /* Particle ID */
CsI[i].chan = *current_word++; /* Particle Energy */
CsI[i].tac = *current_word++; /* Particle Zero Crossing Time */
//inc_hitpat_CsI(CsI[i].id);
}
if (*current_word != 0xFFFF) {
goto Bad_Event;
}
good_event++;
goto Next_Event;
Bad_Event :
bad_event++;
/* Find the event separator */
if ( *current_word != 0xFFFF) {
while ( (*current_word != 0xFFFF) && (current_word <= end_of_data) )
current_word++;
}
Next_Event :
/* move to start of next event */
if (current_word == end_of_data) {
break;
} else {
/* move to start of next event */
current_word++;
}
/* end of buffer ? FFFF padding */
if (*current_word == 0xFFFF) {
//printf("\n padding \n");
break;
}
}
//printf("\n nevents %d",nevents);
// Timing
if (nBlocks%printev == 0) {
tnew = timer.RealTime();
printf("Block:%d, rtime=%.2f s -> %.2f MBytes/s \n",
nBlocks,tnew-told,(0.016*printev)/(tnew-told));
told=tnew;
timer.Continue();
//hfile->Write();
}
}
fclose(i_f);
hfile->Write();
//hfile->Close();
return nBlocks;
}
This archive was generated by hypermail 2b29 : Tue Jan 01 2002 - 17:51:03 MET