Hello ROOTpeople
I think I have found a bug in TClonesArray, but probably I just still
haven't
understood how this thing works...
I have attached two small classes which I have been using to test
TClonesArray-based event classes. The behaviour I am trying to get from
these classes is
1. several different "events" can co-exist, including events of
classes derived from some base "event" class, so no static
TClonesArray a la $ROOTSYS/test/Event.cxx
2. from one event to the next, the number of objects in the
TClonesArray (i.e. the multiplicity of nuclei in the event) can
vary from 1 to 100 (rough order of magnitude). The TClonesArray
therefore has to change size and adapt to the current event, in
order to optimize the size of the resulting TTree.
After some work (about 6 months ... :-D ) the working solution I have come
up with (see MyEvt.cxx) is, roughly speaking, the following:
* initialise TClonesArray with an arbitrary number of slots:
TClonesArray("MyNuc",50) (BTW, does the number of slots - 50 here
- have any importance whatsoever ?)
* before looping over events, call ExpandCreate(1) in order to
create 1 object in the array
* for each new particle in each event, call ExpandCreate(mult) with
"mult" the new multiplicity if "mult" is larger than the
multiplicity of any previous event, otherwise just use the
existing array
* once the event is set up (all nuclei created) call
ExpandCreate(mult) with "mult" the final multiplicity in order to
set the size of the array i.e. for writing in a TTree. If the
current event is smaller than the previous one, this will
"liberate" the extra nuclei, otherwise it should have no effect
* before starting a new event, call Clear("C") to wipe all nuclei
(liberate any allocated memory if necessary) and then call
ExpandCreate(1) again in order to have 1 nucleus in the array
ready for a new event.
This seems to work fine (see Main.cxx for example of use), the low
compression factor of the TTree branch for the events tells me that
empty slots were not written (i.e. for low multiplicity events), as does
the fact that if I use the TTree to plot, e.g., data.fParticles.GetZ()
then I do not see a huge peak at Z=0 corresponding to the empty slots.
However I am a little worried by the following behaviour. In an
interactive session I load libPhysics.so and do ".L yNuc.cxx+". Notice
that MyNuc has a static data member to count the number of existing
objects. So if I do:
root [3] MyNuc a
root [4] a.GetNbObj()
I get the answer:
(Int_t)1
Now I create a TClonesArray and initialise it with 5 MyNuc objects:
root [5] TClonesArray tca("MyNuc",50)
root [6] tca.ExpandCreate(5)
If I now type
root [7] a.GetNbObj()
(Int_t)6
I see I have 6 MyNuc objects, i.e. "a" plus the 5 in the array. OK. Here
are the objects:
root [8] (MyNuc*)tca.At(0)
(class MyNuc*)0x894da88
root [9] (MyNuc*)tca.At(1)
(class MyNuc*)0x89f6668
root [10] (MyNuc*)tca.At(2)
(class MyNuc*)0x89f6488
root [11] (MyNuc*)tca.At(3)
(class MyNuc*)0x89f6598
root [12] (MyNuc*)tca.At(4)
(class MyNuc*)0x89f66f0
root [13] (MyNuc*)tca.At(5)
(class MyNuc*)0x0
(the last line was just to make sure that there are only 5 objects in
the array ;-) ).
Now let's suppose that the next event has only 1 nucleus. I do:
root [14] tca.Clear("C")
and I see that all the slots are empty:
root [15] (MyNuc*)tca.At(0)
(class MyNuc*)0x0
root [16] (MyNuc*)tca.At(4)
(class MyNuc*)0x0
but the objects still exist (this took me a long time to understand):
root [17] a.GetNbObj()
(Int_t)6
Now I set the size of the event to 1:
root [18] tca.ExpandCreate(1)
root [19] a.GetNbObj()
(Int_t)6
root [20] (MyNuc*)tca.At(0)
(class MyNuc*)0x894da88
root [21] (MyNuc*)tca.At(1)
(class MyNuc*)0x0
I see that no new events have been created (good !), that the first slot
of the array has been filled,
and that the object in the slot is the same one that was occupying this
slot in the previous event (good again !).
Now let's suppose that the next event is multiplicity 3:
root [22] tca.Clear("C")
root [23] tca.ExpandCreate(3)
root [24] (MyNuc*)tca.At(0)
(class MyNuc*)0x894da88
root [25] (MyNuc*)tca.At(1)
(class MyNuc*)0x89f66f0
root [26] (MyNuc*)tca.At(2)
(class MyNuc*)0x89f6598
root [27] a.GetNbObj()
(Int_t)8
Two new objects have been created (the difference between the old
multiplicity - 1 - and the new - 3 - ?)!! However, the slots in the
array are filled with the old objects, although not in the same order as
before: slot "0" has not changed, but slots "1" and "2" now contain the
objects from the old slots "4" and "3" respectively.... I was
expecting/hoping that no new objects would be created (no need), and
that the 3 slots in this event would be the same as the first three
slots of the event with multiplicity 5 (this isn't essential, but it
seems a little strange to mix them up like this).
If I try to get back to the initial event size, i.e. multiplicity 5, the
following occurs:
root [29] tca.Clear("C")
root [30] tca.ExpandCreate(5)
root [31] a.GetNbObj()
(Int_t)10
root [32] (MyNuc*)tca.At(0)
(class MyNuc*)0x894da88
root [33] (MyNuc*)tca.At(1)
(class MyNuc*)0x89f66f0
root [34] (MyNuc*)tca.At(2)
(class MyNuc*)0x89f6598
root [35] (MyNuc*)tca.At(3)
(class MyNuc*)0x89f6488
root [36] (MyNuc*)tca.At(4)
(class MyNuc*)0x89f6668
root [37] (MyNuc*)tca.At(5)
(class MyNuc*)0x0
Once again, the fact of using ExpandCreate(n) followed by
ExpandCreate(n+i) creates i objects, even if (n+i) is smaller or the
same as the number of objects created in a previous call. Either this
behaviour is wrong or I haven't understood. However, there is no sign of
the new objects in the array slots : all slots are filled with
pre-existing objects (albeit in a different order).
Now, if I execute the MainTest() example in Main.cxx I do not see the
number of objects increasing up to infinity (or segmentation fault,
whichever is sooner ;-) ) - rather the number of objects simply
increases up to the size of the largest event treated and no greater.
This is the "correct" behaviour - at least, it's the one I want. I then
realised that in MyEvt I use not "ExpandCreate" but "ExpandCreateFast",
and indeed if I repeat all the previous instructions replacing
ExpandCreate by ExpandCreateFast there is no more memory leak.
Is this difference normal ?
Best regards
John
--
ganil logo <http://www.ganil.fr>
John D. Frankland <mailto:frankland@ganil.fr>
Beam Coordinator
GANIL
B.P. 55027
14076 CAEN Cedex 05
*tel:* +33 (0)231454628
*fax:* +33 (0)231454665
#ifndef __MyNuc_H
#define __MyNuc_H
#include "TLorentzVector.h"
class MyNuc : public TLorentzVector
{
UChar_t fZ;
UChar_t fA;
static Int_t fNb_Objects;
public:
MyNuc();
MyNuc(const MyNuc &);
MyNuc(int z, int a);
virtual ~MyNuc();
void SetZ(int z){fZ=(UChar_t)z;}
void SetA(int a){fA=(UChar_t)a;}
Int_t GetZ() const;
Int_t GetA() const;
Int_t GetNbObj();
void SetPxPyPz(Float_t px, Float_t py, Float_t pz);
void Clear(Option_t *opt="");
ClassDef(MyNuc,1) // Just for tests
};
#endif
#include "MyNuc.h"
ClassImp(MyNuc)
//__________________________________________________________________
//MyNuc
//For testing TClonesarray behaviour
//
Int_t MyNuc::fNb_Objects=0;
MyNuc::MyNuc()
{
fZ = fA = 0;
fNb_Objects++;
}
MyNuc::MyNuc(int z, int a)
{
fZ = (UChar_t)z;
fA = (UChar_t)a;
fNb_Objects++;
}
MyNuc::MyNuc(const MyNuc &obj)
{
fZ = obj.GetZ();
fA = obj.GetA();
}
MyNuc::~MyNuc()
{
fNb_Objects--;
}
Int_t MyNuc::GetZ() const
{
return (Int_t)fZ;
}
Int_t MyNuc::GetA() const
{
return (Int_t)fA;
}
Int_t MyNuc::GetNbObj()
{
return fNb_Objects;
}
void MyNuc::Clear(const Option_t *opt)
{
TLorentzVector::Clear(opt);
fZ=fA=0;
}
void MyNuc::SetPxPyPz(Float_t px, Float_t py, Float_t pz)
{
Float_t mass = 931.5*GetA();
TLorentzVector::SetXYZM(px,py,pz,mass);
}
#ifndef __MyEvt_H
#define __MyEvt_H
#define DEFAULT_EVENT_SIZE 50
#include "TObject.h"
#include "TClonesArray.h"
#include "MyNuc.h"
class MyEvt : public TObject
{
Int_t fSizeArray;// used to manage TClonesArray
Int_t fMult;//multiplicity of event
TClonesArray *fParticles;//-> array of nuclei in event
public:
MyEvt();
MyEvt(Int_t mult);
~MyEvt();
void SetSizeEvent(Int_t mult);
MyNuc *GetNextParticle();
MyNuc *GetParticle(Int_t npart) const;
Int_t GetMult(){ return fMult; }
void Clear(const Option_t *opt="");
ClassDef(MyEvt,1);
};
#endif
#include "MyEvt.h"
#include <iostream>
using std::cout;
using std::endl;
ClassImp(MyEvt)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MyEvt::MyEvt()
{
fMult = 0;
fSizeArray = 0;
fParticles = 0;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MyEvt::MyEvt(Int_t mult)
{
fMult = 0;
SetSizeEvent(mult);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MyEvt::~MyEvt()
{
if(fParticles){
fParticles->Delete();
delete fParticles;
fParticles=0;
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void MyEvt::SetSizeEvent(Int_t mult)
{
//Create TClonesArray if necessary.
//Set size of TClonesArray to mult ready for filling.
if(!fParticles){
fParticles = new TClonesArray("MyNuc", DEFAULT_EVENT_SIZE);
}
fParticles->ExpandCreateFast(mult);
fSizeArray = mult;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MyNuc *MyEvt::GetParticle(Int_t npart)const
{
//
//Access to event member with index npart (1<=npart<=fMult)
//
if(fParticles){
if(npart<1 || npart>fMult)
{
Warning("GetParticle","Particle index %d out of range (1-%d)",npart,fMult);
return 0;
}
return (MyNuc*)fParticles->At(npart-1);
}
else
{
Warning("GetParticle", "No particles in event");
}
return 0;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
MyNuc* MyEvt::GetNextParticle()
{
//Only use when filling a new event, not for reading the event.
//Increases multiplicity by one and returns pointer to next particle in array
//Increases size of TClonesArray if necessary
fMult++;
MyNuc *tmp;
if(fMult<=fSizeArray)
{
tmp = GetParticle(fMult);
}
else
{
//increase size of TClonesArray
SetSizeEvent(fMult);
tmp = GetParticle(fMult);
}
if( tmp )
{
return tmp;
}
else
{
Warning("GetNextParticle","TClonesArray is too small, fMult=%d",fMult);
fMult--;
return 0;
}
}
void MyEvt::Clear(Option_t *opt)
{
fParticles->Clear("C"); // call Clear for all nuclei in event
SetSizeEvent(1); // reset TClonesArray to contain 1 nucleus
fMult = 0; // reset event multiplicity to 0
}
#include <iostream>
using std::cout;
using std::endl;
#include "TFile.h"
#include "TTree.h"
#include "TRandom.h"
#include "TMath.h"
#include "MyNuc.h"
#include "MyEvt.h"
void MainTest(){
MyNuc a;
cout << "NbObj=" << a.GetNbObj() << endl;
int event, nobj;
MyEvt *mevt = new MyEvt();
cout << "NbObj=" << a.GetNbObj() << endl;
TFile f("test.root","recreate");
TTree t("tree","Test Tree");
t.Branch("data", "MyEvt", &mevt, 64000, 0);
t.Branch("numevt", &event, "event/I");
t.Branch("numobj", &nobj, "nobj/I");
mevt->SetSizeEvent(1);
cout << "NbObj=" << a.GetNbObj() << endl;
for(event=1; event<200000; event++)
{
int mult = (int)gRandom->Gaus(20., 15.);
mult = TMath::Abs(mult) + 1;
for(int part=1; part<=mult; part++)
{
int z = (int)(gRandom->Uniform(0.,50.) + 1);
int a = 2 * z;
float px = gRandom->Gaus(0., 270.);
float py = gRandom->Gaus(0., 270.);
float pz = gRandom->Gaus(0., 270.);
MyNuc *tmp = mevt->GetNextParticle();
if(!tmp){
mevt->Error("GetNextParticle", "Out of memory ?");
exit(1);
}
tmp->SetZ(z);
tmp->SetA(a);
tmp->SetPxPyPz(px, py, pz);
}
mevt->SetSizeEvent(mult);
nobj = a.GetNbObj();
t.Fill();
mevt->Clear();
if(event && !(event%1000)) {
cout << "Event " << event;
cout << " ---- NbObj=" << a.GetNbObj() << endl;
}
}
t.Write();
f.Close();
}
This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:13 MET