[ROOT] Problem with global variable

From: Perfetto Francesco (Francesco.Perfetto@na.infn.it)
Date: Wed Jul 07 2004 - 12:12:04 MEST


Hi All,

I have the following problem:
I have a program (compiled) in which there is a main (Likelihood.cpp)
that call the likelihood.cc in which there are any function.

Likelihood.cpp:

#include <iostream.h>
#include "TROOT.h"
#include "TObject.h"

void llhood(Double_t ALFA);

int main(int argc, char* argv[])
{
  Double_t ALFA[] = {0,-.45,-.07};//,.04,-.07,.07,-.12,.12,-.22,.22,-.45,.45,-.78,.97,1,.8};

  for(Int_t i=0; i<3; i++){
    llhood(ALFA[i]);
  }

  return 0;

}

likelihood.cc:

#define res_class_cxx
#include "res_class.h"
#include <iostream.h>
#include "TMath.h"
#include "TFormula.h"
#include "TH1.h"
#include "TF1.h"
#include "TROOT.h"
#include "TLorentzVector.h"
#include "TVectorD.h"
#include "TObject.h"
#include "TFile.h"
#include "TTree.h"
#include "TMinuit.h"
#include "TGraph.h"
#include "TCut.h"
#include "TChain.h"


Double_t func(Double_t *xv,Double_t *par)
{
  Double_t z = xv[0];
  return (1. + par[0]*z);
}

const Int_t nparams = 1;
const Int_t TreGaus = 10;     //  <<<<<-----------RICORDA DI CAMBIARE IL
NUMERO DI 3G
const Int_t TreGaus = 10;     //  <<<<<-----------RICORDA DI CAMBIARE IL
NUMERO DI 3G
const Int_t TreGaus = 10;     //  <<<<<-----------RICORDA DI CAMBIARE IL
NUMERO DI 3G
const Double_t PI = 3.1415926;

//apro il file in cui salvo gli istogrammi per controllare la bonta' del
fit
TFile *f2 = new
TFile("/home/perfetto/AUTO/tmpLLHOOD/llhood_camp2_10x60slices.root","RECREATE");

// istogrammi calcolati con il tree t che  utilizzo per il folding
(riempire hist.ris)
TH1D *gen  = new TH1D("gen" ,"gen" ,nslices,0,1);
TH1D *rgen = new TH1D("rgen","rgen",nslices,0,1);

// istogrammi che calcolo con il tree t_eff che utilizzo per calcolare
l'efficienza (nel folding)
TH1D *gen_eff  = new TH1D("gen_eff" ,"gen_eff" ,nslices,0,1);
TH1D *rgen_eff = new TH1D("rgen_eff","rgen_eff",nslices,0,1);
TH1D *eff  = new TH1D("eff" ,"eff" ,nslices,0,1);

// istogramma che calcolo con il tree t_eff che non utilizzo
TH1D *rec_eff  = new TH1D("rec_eff" ,"rec_eff" ,2*nslices,0,2);
// istogramma che calcolo con il tree t che non utilizzo
TH1D *rec  = new TH1D("rec" ,"rec" ,2*nslices,0,2);

// istogramma della variabile z_rec pesata per la funzione teorica
TH1D *rec_alfa  = new TH1D("rec_alfa" ,"rec_alfa" ,2*nslices,0,2);

// istogramma che contiene gli eventi dopo il folding
TH1D *ris  = new TH1D("ris" ,"ris" ,2*nslices,0,2);

//definizione della f teorica che utilizzo nel folding
TF1 *f_th = new TF1("f_th",func,0,2,nparams);



//per sapere quante Fermi Dirac ci sono tra le slices e quante ne uso per
il fit

Int_t last_FD = FermiDirac;
Int_t usa_FD = 1;


// definizione della Fermi Dirac e della tripla gaussiana
Double_t norm_FD[] = {225.69};
Double_t parFD[4][FermiDirac];
Double_t par3G[9][TreGaus];

void riempi()
{
...
}

void fcnllhood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par,
Int_t iflag)
{
  // Calcolo -ln(llhood)da minimizzare.
  // calcolo gli eventi convoluti (dopo il folding)

  // il numero di bin e' pari al numero di slices
  Int_t nbins = nslices;
  // max e min bin nel fit
  Int_t maxbinfitted = 16;  // <<<<<<--------RICORDA DI CAMBIARE I VALORI
  Int_t minbinfitted = 4;   //
...
	//eventi puro spazio delle fasi
	Double_t c_gen = gen->GetBinContent(i)/n_ministeps;
	// efficienza
 	Double_t c_eff = eff->GetBinContent(i);
	//funzione teorica
	Double_t fvalue = f_th->Eval(zp);
	//eventi convoluti
	c_prod += (c_gen * c_eff * c_gz)*fvalue;
      }

    }
    ris->SetBinContent(j,c_prod);

  }
...
    //Double_t n_misj = rec->GetBinContent(j);
    Double_t n_misj = rec_alfa->GetBinContent(j);

    llh += nj*n_misj;
  }
  // funzione da minimizzare
  f = -2.*llh;
}

void llhood(Double_t ALFA)
{

  gen ->Sumw2();
  rgen->Sumw2();
  eff ->Sumw2();
  rec ->Sumw2();
  rec_alfa->Sumw2();
  ris ->Sumw2();

  // funzione che riempie le variabili con i parametri dai fit alle slices
  riempi();


  //creo una Chain di files MC con cui faccio le efficienze
  TChain *t_eff = new TChain("/data2/perfetto/bunch_1_minuit.root/res");
  t_eff->Add("/data2/perfetto/bunch_1_minuit.root/res");
  //creo una Chain di files MC con cui faccio il fit
  TChain *t = new TChain("/data2/perfetto/bunch_4_minuit.root/res");
  t->Add("/data2/perfetto/bunch_4_minuit.root/res");
...
  res_class *resc = new res_class(t);
  res_class *resc_eff = new res_class(t_eff);

  // per conoscere il numero di entries del tree t
  Int_t nentries = Int_t(t->GetEntries());
  Int_t nbytes = 0, nb = 0;
  cout<<nentries<<endl;
  // Loop per riempire tutti gli istogrammi tranne le efficienze
  for (Int_t jentry=0; jentry<nentries; jentry++) {
    Int_t ientry = resc->LoadTree(jentry);
    nb = resc->GetEntry(jentry);
    nbytes += nb;

    Bool_t bc1 = resc->eta_fit_pchi2 > 0.01;
    Bool_t bc2 = resc->chi_2_ord[0] < 5;
    Bool_t bc3 = (resc->chi_2_ord[1]-resc->chi_2_ord[0]) > 3;
    Bool_t bc4 = resc->eta_filter_flag==-1;
    Bool_t bc5 = resc->p4_ph_recA__fE[0]>320 &&
resc->p4_ph_recA__fE[0]<400;
    Bool_t bc6 = resc->isecl==1;
    Bool_t bcene = resc->eta_ini_par[5]>10 && resc->eta_ini_par[10]>10 &&
resc->eta_ini_par[15]>10 && resc->eta_ini_par[20]>10 &&
resc->eta_ini_par[25]>10 && resc->eta_ini_par[30]>10;
    Bool_t bcang = resc->angolo_ini_ord[0] > 18;

    Bool_t bct = bc1&&bc2&&bc3&&bc4&&bc5&&bc6&&bcene&&bcang;

    gen->Fill(resc->daliz_gen);
    Double_t weight = (1.+ ALFA*(resc->daliz_gen));
    if(bct){
    Double_t zalfa = resc->daliz_rec_fit;
    rec_alfa->Fill(zalfa,weight);

    rec->Fill(resc->daliz_rec_fit);
    rgen->Fill(resc->daliz_gen);

    }

  }

  Int_t nentries_eff = Int_t(t_eff->GetEntries());
  Int_t nbytes_eff = 0, nb_eff = 0;
  cout<<nentries_eff<<endl;
  // Loop per riempire tutti gli istogrammi per le efficienze
  for (Int_t jentry=0; jentry<nentries_eff; jentry++) {
    Int_t ientry = resc_eff->LoadTree(jentry);
    nb_eff = resc_eff->GetEntry(jentry);
    nbytes_eff += nb_eff;

    Bool_t bc1 = resc_eff->eta_fit_pchi2 > 0.01;
    Bool_t bc2 = resc_eff->chi_2_ord[0] < 5;
    Bool_t bc3 = (resc_eff->chi_2_ord[1]-resc_eff->chi_2_ord[0]) > 3;
    Bool_t bc4 = resc_eff->eta_filter_flag==-1;
    Bool_t bc5 = resc_eff->p4_ph_recA__fE[0]>320 &&
resc_eff->p4_ph_recA__fE[0]<400;
    Bool_t bc6 = resc_eff->isecl==1;
    Bool_t bcene = resc_eff->eta_ini_par[5]>10 &&
resc_eff->eta_ini_par[10]>10 && resc_eff->eta_ini_par[15]>10 &&
resc_eff->eta_ini_par[20]>10 && resc_eff->eta_ini_par[25]>10 &&
resc_eff->eta_ini_par[30]>10;
    Bool_t bcang = resc_eff->angolo_ini_ord[0] > 18;

    Bool_t bct = bc1&&bc2&&bc3&&bc4&&bc5&&bc6&&bcene&&bcang;

    gen_eff->Fill(resc_eff->daliz_gen);
    if(bct){
      rec_eff->Fill(resc_eff->daliz_rec_fit);
      rgen_eff->Fill(resc_eff->daliz_gen);
    }

  }
  // calcolo l'efficienza
  eff->Divide(rgen_eff,gen_eff,1,1,"B");


  // USO DI MINUIT.

  // Inizializza minuit con nparams parametri.
  TMinuit *minuit = new TMinuit(nparams);

  // Fissa la fcn da minimizzare.
  minuit->SetFCN(fcnllhood);

  Double_t arglist[100];
  Int_t ierflag = 0;

  // Settiamo i parametri e lo step

  Double_t start[nparams] = {0.35};
  Double_t step[nparams]  = {0.01};


  minuit->mnparm(0,"alpha",start[0],step[0],0,0,ierflag);

  // arglist = lista degli argomenti del comando che chiamo con mnexcm.

  arglist[0] = 1;
  minuit->mnexcm("CALL FCN",arglist,1,ierflag);

  arglist[0] = 0;
  minuit->mnexcm("MIGRAD",arglist,0,ierflag);
  minuit->mnexcm("MINOS" ,arglist,0,ierflag);
  minuit->mnexcm("SET PRINT",arglist,1,ierflag);

  minuit->mnscan();

  // Stampa dei risultati.
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  minuit->mnprin(4,amin);
  cout<<"icstat = "<<icstat<<endl;
  Double_t alfa,err_alfa;
  minuit->GetParameter(0,alfa,err_alfa);
  cout<<"alpha = "<<alfa<<endl;
  cout<<"err_alpha = "<<err_alfa<<endl;

  f2->Write();

  gen -> Delete();
  rgen -> Delete();
  gen_eff -> Delete();
  rgen_eff -> Delete();
  eff -> Delete();
  rec_eff -> Delete();
  rec -> Delete();
  rec_alfa -> Delete();
  ris -> Delete();
  f_th -> Delete();
}

The problem is that when run the program (compiling with gmake OK) the
first value of ALFA (0) don't give any problem, but when keep the
following value if I don't delate the histos (how do now at the end of
program) the program refill the histos again on the sames (I give a bad
value of fit) But if I delate the histos when run the program I have at
second value the following error:

lxkloe2:~/AUTO/tmpLLHOOD> Linux/Likelihood.exe > pippo2.txt

 *** Break *** segmentation violation
 Generating stack trace...
 0x0804943a in main + 0x2e from Linux/Likelihood.exe
 0x42015704 in __libc_start_main + 0xe4 from /lib/tls/libc.so.6
 0x0804933d in TFile::TFile[in-charge](char const*, char const*, char
const*, int) + 0x31 from Linux/Likelihood.exe
Segnale di annullamento (creato file core)

Can someone help me?
 Thanks
Francesco.

I use root 4.00/04 on redhat 9.0



This archive was generated by hypermail 2b29 : Sun Jan 02 2005 - 05:50:08 MET