[ROOT] Scope problem : save histograms

From: K. Hauschild (karlhaus@hep.saclay.cea.fr)
Date: Thu Oct 18 2001 - 18:57:16 MEST


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