Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 5691) +++ trunk/MagicSoft/Mars/Makefile (revision 5692) @@ -58,4 +58,5 @@ mimage \ mmontecarlo \ + mhft \ mmc \ mraw \ Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 5691) +++ trunk/MagicSoft/Mars/NEWS (revision 5692) @@ -20,4 +20,8 @@ MGeomCamMagic to a enhanced geometry MGeomCamMagicXT having only small pixels + + - added possibility to write data to memory (TTree) using MReadTree + + - added possibility to read a TTree stored only in memory by MReadTree Index: trunk/MagicSoft/Mars/mbase/MLog.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MLog.cc (revision 5691) +++ trunk/MagicSoft/Mars/mbase/MLog.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 12/2000 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -116,4 +116,8 @@ using namespace std; +#undef DEBUG +//#define DEBUG + + // root 3.02: // check for TObjectWarning, TObject::Info, gErrorIgnoreLevel @@ -216,5 +220,22 @@ DeallocateFile(); +#ifdef DEBUG + TIter Next(fPlugins); + TObject *o=0; + while ((o=Next())) + { + cout << "Delete: " << o->GetName() << std::flush; + cout << " [" << o->ClassName() << "]" << endl; + delete o; + } + + cout << "Delete: fPlugins " << fPlugins << "..." << std::flush; +#endif + delete fPlugins; +#ifdef DEBUG + cout << "done." << endl; +#endif + #ifdef _REENTRANT delete fMuxStream; Index: trunk/MagicSoft/Mars/mbase/MParContainer.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MParContainer.cc (revision 5691) +++ trunk/MagicSoft/Mars/mbase/MParContainer.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 12/2000 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -109,5 +109,21 @@ { #ifdef DEBUG + if (fName.IsNull() || fName==(TString)"MTime") + return; + *fLog << all << "Deleting " << GetDescriptor() << endl; + if (TestBit(kMustCleanup) && gROOT && gROOT->MustClean()) + { + *fLog << "Recursive Remove..." << flush; + if (TestBit(kCanDelete)) + *fLog << "kCanDelete..." << flush; + TIter Next(gROOT->GetListOfCleanups()); + TObject *o=0; + while ((o=Next())) + *fLog << dbg << o->GetName() << " [" << o->ClassName() << "]" << endl; + *fLog << dbg << "Removing..." << flush; + gROOT->GetListOfCleanups()->RecursiveRemove(this); + *fLog << "Removed." << endl; + } #endif } Index: trunk/MagicSoft/Mars/mbase/MTaskList.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MTaskList.cc (revision 5691) +++ trunk/MagicSoft/Mars/mbase/MTaskList.cc (revision 5692) @@ -879,2 +879,51 @@ } + +// -------------------------------------------------------------------------- +// +// Find an object with the same name in the list and replace it with +// the new one. If the kIsOwner flag is set and the object was not +// created automatically, the object is deleted. +// +Bool_t MTaskList::Replace(MTask *task) +{ + // + // check if the object (you want to add) exists + // + if (!task) + return kFALSE; + + if (task==this) + { + *fLog << warn << "WARNING - You cannot add a tasklist to itself. This" << endl; + *fLog << " would create infinite recursions...ignored." << endl; + return kFALSE; + + } + + MTask *obj = (MTask*)FindObject(task->GetName()); + if (!obj) + { + *fLog << warn << "No object with the same name '"; + *fLog << task->GetName() << "' in list... adding." << endl; + return AddToList(task); + } + + if (task==obj) + return kTRUE; + + *fLog << inf << "Adding " << task->GetName() << " to " << GetName() << " for " << obj->GetStreamId() << "... " << flush; + task->SetStreamId(obj->GetStreamId()); + task->SetBit(kMustCleanup); + fTasks->AddAfter((TObject*)obj, task); + *fLog << "Done." << endl; + + RemoveFromList(obj); + + if (TestBit(kIsOwner)) + delete obj; + + //*fLog << inf << "MTask '" << task->GetName() << "' found and replaced..." << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mbase/MTaskList.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MTaskList.h (revision 5691) +++ trunk/MagicSoft/Mars/mbase/MTaskList.h (revision 5692) @@ -49,4 +49,5 @@ void SetSerialNumber(Byte_t num); + Bool_t Replace(MTask *obj); Bool_t RemoveFromList(MTask *task); Index: trunk/MagicSoft/Mars/mdata/MDataArray.cc =================================================================== --- trunk/MagicSoft/Mars/mdata/MDataArray.cc (revision 5691) +++ trunk/MagicSoft/Mars/mdata/MDataArray.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 08/2002 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -86,4 +86,30 @@ // -------------------------------------------------------------------------- // +// Try to find a MData which has the same rule (GetRule()). +// Be carefull: This may already fail if the rule is A*B and you are +// searching for B*A - FIXME: A more intelligent comaparton (based on +// IsEqual()) should be developed. +// +// Returns the index in the array, -1 if rule was not found. +// +Int_t MDataArray::FindRule(const char *rule) const +{ + const MDataChain data(rule); + if (!data.IsValid()) + return -1; + + const TString r(data.GetRule()); + + TIter Next(&fList); + const MData *d = NULL; + while ((d=(MData*)Next())) + if (d->GetRule()==r) + return fList.IndexOf(d); + + return -1; +} + +// -------------------------------------------------------------------------- +// // Return the i-th entry // @@ -131,9 +157,11 @@ Int_t n=0; - TIter Next(&fList); - MData *data = NULL; - while ((data=(MData*)Next())) - { - *fLog << all << " " << fName << "[" << setw(3) << n++ << "] = " << flush; + const Int_t w = 1+(Int_t)TMath::Log10(fList.GetEntries()); + + TIter Next(&fList); + MData *data = NULL; + while ((data=(MData*)Next())) + { + *fLog << all << " " << fName << "[" << setw(w) << n++ << "] = " << flush; data->Print(); *fLog << endl; Index: trunk/MagicSoft/Mars/mdata/MDataArray.h =================================================================== --- trunk/MagicSoft/Mars/mdata/MDataArray.h (revision 5691) +++ trunk/MagicSoft/Mars/mdata/MDataArray.h (revision 5692) @@ -32,4 +32,6 @@ void AddEntry(MData *data); + Int_t FindRule(const char *rule) const; + MData &operator[](Int_t i) const; Double_t operator()(Int_t i) const; Index: trunk/MagicSoft/Mars/mdata/MDataChain.cc =================================================================== --- trunk/MagicSoft/Mars/mdata/MDataChain.cc (revision 5691) +++ trunk/MagicSoft/Mars/mdata/MDataChain.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 04/2002 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -211,5 +211,5 @@ return; - *fLog << inf << "Trying to resolve rule... " << flush; + *fLog << inf << "Parsing rule... " << flush; if (!(fMember=ParseString(rule, 1))) { Index: trunk/MagicSoft/Mars/mfbase/MF.cc =================================================================== --- trunk/MagicSoft/Mars/mfbase/MF.cc (revision 5691) +++ trunk/MagicSoft/Mars/mfbase/MF.cc (revision 5692) @@ -16,7 +16,7 @@ ! ! -! Author(s): Thomas Bretz 01/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 +! Author(s): Thomas Bretz 01/2002 +! +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -122,5 +122,5 @@ fTitle = title ? title : gsDefTitle.Data(); - *fLog << inf << "Trying to resolve filter rule... " << flush; + *fLog << inf << "Parsing filter rule... " << flush; if (!(fF=ParseString(text, 1))) { Index: trunk/MagicSoft/Mars/mfileio/MReadTree.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 5691) +++ trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 12/2000 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -24,28 +24,33 @@ ///////////////////////////////////////////////////////////////////////////// -// // -// MReadTree // -// // -// This tasks opens all branches in a specified tree and creates the // -// corresponding parameter containers if not already existing in the // -// parameter list. // -// // -// The Process function reads one events from the tree. To go through the // -// events of one tree make sure that the event number is increased from // -// outside. It makes also possible to go back by decreasing the number. // -// // -// If you don't want to start reading the first event you have to call // -// MReadTree::SetEventNum after instantiating your MReadTree-object. // -// // -// To make reading much faster (up to a factor of 10 to 20) you can // -// ensure that only the data you are really processing is enabled by // -// calling MReadTree::UseLeaf. // -// // -// If the chain switches from one file to another file all // -// TObject::Notify() functions are called of TObject objects which were // -// added to the Notifier list view MReadTree::AddNotify. If MReadTree // -// is the owner (viw MReadTree::SetOwner) all this objects are deleted // -// by the destructor of MReadTree // -// // +// +// MReadTree +// +// This tasks opens all branches in a specified tree and creates the +// corresponding parameter containers if not already existing in the +// parameter list. +// +// The Process function reads one events from the tree. To go through the +// events of one tree make sure that the event number is increased from +// outside. It makes also possible to go back by decreasing the number. +// +// If you don't want to start reading the first event you have to call +// MReadTree::SetEventNum after instantiating your MReadTree-object. +// +// To make reading much faster (up to a factor of 10 to 20) you can +// ensure that only the data you are really processing is enabled by +// calling MReadTree::UseLeaf. +// +// If the chain switches from one file to another file all +// TObject::Notify() functions are called of TObject objects which were +// added to the Notifier list view MReadTree::AddNotify. If MReadTree +// is the owner (viw MReadTree::SetOwner) all this objects are deleted +// by the destructor of MReadTree +// +// +// ToDo: +// ----- +// - Auto Scheme and Branch choosing doesn't work for memory trees +// ///////////////////////////////////////////////////////////////////////////// #include "MReadTree.h" @@ -73,7 +78,9 @@ // -------------------------------------------------------------------------- // -// Default constructor. Don't use it. -// -MReadTree::MReadTree() +// Default constructor. Use this constructor (ONLY) if you want to +// access a memory tree (a TTree not stored in a file) created by +// MWriteRootFile or manually. +// +MReadTree::MReadTree(TTree *tree) : fNumEntry(0), fNumEntries(0), fBranchChoosing(kFALSE), fAutoEnable(kTRUE) { @@ -81,8 +88,25 @@ fTitle = "Task to loop over all events in one single tree"; - fVetoList = NULL; - fNotify = NULL; + fVetoList = new TList; + fVetoList->SetOwner(); + + fNotify = new TList; fChain = NULL; + + fTree = tree; + SetBit(kChainWasChanged); + + if (fTree) + fAutoEnable = kFALSE; + /* + if (!fTree) + return; + + TIter Next(fTree->GetListOfBranches()); + TBranch *b=0; + while ((b=(TBranch*)Next())) + b->ResetAddress(); + */ } @@ -114,4 +138,5 @@ // fChain = new MChain(tname); + fTree = fChain; // root 3.02: @@ -122,6 +147,4 @@ if (fname) AddFile(fname); -// if (fChain->Add(fname)>0) -// SetBit(kChainWasChanged); } @@ -137,19 +160,34 @@ // creates a memory leak! // - TIter Next(fChain->GetStatus()); - - TChainElement *element = NULL; - while ((element=(TChainElement*)Next())) - if (element->GetBaddress()) - delete (MParContainer**)element->GetBaddress(); - - // - // Delete the chain and the veto list - // + if (fChain) // FIXME: MEMORY LEAK for fTree!=0 + { + TIter Next(fChain->GetStatus()); + + TChainElement *element = NULL; + while ((element=(TChainElement*)Next())) + if (element->GetBaddress()) + delete (MParContainer**)element->GetBaddress(); + + // + // Delete the chain and the veto list + // #if ROOT_VERSION_CODE < ROOT_VERSION(3,03,00) - if (fChain->GetFile()) - delete fChain->GetFile(); + if (fChain->GetFile()) + delete fChain->GetFile(); #endif - delete fChain; + delete fChain; + } + + /* This is AND MUST be done in PostProcess. See there for more details + else + { + TIter Next(fTree->GetListOfBranches()); + + TBranch *element = NULL; + while ((element=(TBranch*)Next())) + if (element->GetAddress()) + delete (MParContainer**)element->GetAddress(); + } + */ delete fNotify; @@ -179,22 +217,37 @@ Bool_t MReadTree::CheckBranchSize() { - TArrayI entries(fChain->GetStatus()->GetSize()); + //if (!fChain) // >FIXME: fTree!=0 + // return kTRUE; + + TArrayI entries(fChain ? fChain->GetStatus()->GetSize() : fTree->GetListOfBranches()->GetSize()); Int_t num=0; // Loop over all branches which have a corresponding container - TIter Next(fChain->GetStatus()); - - TChainElement *element = NULL; - while ((element=(TChainElement*)Next())) - { - // Get branch name and find pointer to corresponding branch - const TString name = element->GetName(); - const TBranch *b = fChain->FindBranch(name); - - // Skip element without corresponding branches (like "*") - if (!b) - continue; - - entries[num++] = (Int_t)b->GetEntries(); + /* + if (fChain) + { + TIter Next(fChain->GetStatus()); + + TChainElement *element = NULL; + while ((element=(TChainElement*)Next())) + { + // Get branch name and find pointer to corresponding branch + const TString name = element->GetName(); + const TBranch *b = fChain->FindBranch(name); + + // Skip element without corresponding branches (like "*") + if (!b) + continue; + + entries[num++] = (Int_t)b->GetEntries(); + } + } + else */ + { + TIter Next(fTree->GetListOfBranches()); + + TBranch *element = NULL; + while ((element=(TBranch*)Next())) + entries[num++] = (Int_t)element->GetEntries(); } @@ -209,5 +262,5 @@ *fLog << " Due to several circumstances (such at a bug in MReadTree or wrong" << endl; *fLog << " usage of the file UPDATE mode) you may have produced a file in which" << endl; - *fLog << " at least two branches in the same tree (" << fChain->GetName() << ") have different" << endl; + *fLog << " at least two branches in the same tree (" << fTree->GetName() << ") have different" << endl; *fLog << " number of entries. Sorry, but this file (" << GetFileName() << ")" << endl; *fLog << " is unusable." << endl; @@ -278,4 +331,12 @@ Int_t MReadTree::AddFile(const char *fname, Int_t entries) { + if (!fChain) + { + *fLog << err << "MReadTree::AddFile - ERROR: You cannot add a file, because MReadTree" << endl; + *fLog << " handles a memory based tree or its default" << endl; + *fLog << " constructor was called." << endl; + return 0; + } + #if ROOT_VERSION_CODE < ROOT_VERSION(3,03,01) // @@ -353,4 +414,12 @@ Int_t MReadTree::AddFiles(const MReadTree &read) { + if (!fChain) + { + *fLog << err << "MReadTree::AddFiles - ERROR: You cannot add a file, because MReadTree" << endl; + *fLog << " handles a memory based tree or its default" << endl; + *fLog << " constructor was called." << endl; + return 0; + } + const Int_t rc = fChain->Add(read.fChain); @@ -376,5 +445,6 @@ void MReadTree::SortFiles() { - fChain->GetListOfFiles()->Sort(); + if (fChain) + fChain->GetListOfFiles()->Sort(); } @@ -392,5 +462,5 @@ *fLog << inf << GetDescriptor() << ": Branch choosing enabled (only enabled branches are read)." << endl; - fChain->SetBranchStatus("*", kFALSE); + fTree->SetBranchStatus("*", kFALSE); // *CHANGED-fChain-to-fTree* fBranchChoosing = kTRUE; } @@ -405,4 +475,10 @@ void MReadTree::EnableBranch(const char *name) { + if (!fChain) + { + *fLog << warn << "MReadTree::EnableBranch - WARNING: EnableBranch doesn't work with memory based trees... ignored." << endl; + return; + } + if (fChain->GetListOfFiles()->GetEntries()==0) { @@ -424,5 +500,5 @@ void MReadTree::SetBranchStatus(const char *name, Bool_t status) { - fChain->SetBranchStatus(name, status); + fTree->SetBranchStatus(name, status); // *CHANGED-fChain-to-fTree* *fLog << inf << (status ? "Enabled" : "Disabled"); @@ -449,5 +525,5 @@ bn.Remove(bn.Length()-1); - if (fChain->GetBranch(bn)) + if (fTree->GetBranch(bn)) // *CHANGED-fChain-to-fTree* SetBranchStatus(name, status); @@ -461,5 +537,5 @@ return; - if (fChain->GetBranch(dot+1)) + if (fTree->GetBranch(dot+1)) // *CHANGED-fChain-to-fTree* SetBranchStatus(dot+1, status); } @@ -575,6 +651,7 @@ if (TestBit(kChainWasChanged)) { - *fLog << inf << "Scanning chain " << fChain->GetName() << "... " << flush; - fNumEntries = (UInt_t)fChain->GetEntries(); + // *CHANGED-fChain-to-fTree* + *fLog << inf << "Scanning chain " << fTree->GetName() << "... " << flush; + fNumEntries = (UInt_t)fTree->GetEntries(); *fLog << fNumEntries << " events found." << endl; ResetBit(kChainWasChanged); @@ -624,5 +701,5 @@ // tasks are not preprocessed. // - fChain->SetNotify(NULL); + fTree->SetNotify(NULL); //*CHANGED-fChain-to-fTree* // @@ -633,5 +710,6 @@ // it crashes. ResetTree makes sure, that the tree number is -1 // - fChain->ResetTree(); + if (fChain) // *CHANGED-fChain-to-fTree* + fChain->ResetTree(); // Maybe this would be enough, but the above might be safer... //if (fChain->GetTreeNumber()==0) @@ -641,5 +719,5 @@ // check for files and for the tree! // - if (!fChain->GetFile()) + if (fChain && !fChain->GetFile()) // *CHANGED-fChain-to-fTree* { *fLog << err << GetDescriptor() << ": No file or no tree with name " << fChain->GetName() << " in file." << endl; @@ -665,5 +743,5 @@ // create the Iterator to loop over all branches // - TIter Next(fChain->GetListOfBranches()); + TIter Next(fTree->GetListOfBranches()); // *CHANGED-fChain-to-fTree* TBranch *branch=NULL; @@ -728,7 +806,23 @@ // If we created one already, delete it. // - TChainElement *element = (TChainElement*)fChain->GetStatus()->FindObject(bname); - if (element) - delete (MParContainer**)element->GetBaddress(); + // *CHANGED-fChain-to-fTree* + if (fChain) + { + TChainElement *element = (TChainElement*)fChain->GetStatus()->FindObject(bname); + if (element) + delete (MParContainer**)element->GetBaddress(); + } + /* This can't be done here for memory trees - see + PostProcess for more details. + else + { + TBranch *branch = (TBranch*)fTree->GetBranch(bname); + if (branch) + { + *fLog << bname << " " << (void*)branch->GetAddress() << " " << pcont << endl; + delete (MParContainer**)branch->GetAddress(); + } + } + */ // @@ -736,5 +830,5 @@ // the actual branch should be stored - enable branch. // - fChain->SetBranchAddress(bname, pcont); + fTree->SetBranchAddress(bname, pcont); // *CHANGED-fChain-to-fTree* *fLog << inf << "Master branch address '" << bname << "' ["; @@ -762,8 +856,14 @@ // PreProcess-function. // - fChain->ResetTree(); - fChain->SetNotify(this); - - return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; + if (fChain) + fChain->ResetTree(); // *CHANGED-fChain-to-fTree* + + fTree->SetNotify(this); // *CHANGED-fChain-to-fTree* + + const Int_t rc = GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; + if (rc!=kTRUE || fChain) + return rc; + + return Notify(); } @@ -775,4 +875,7 @@ void MReadTree::SetReadyToSave(Bool_t flag) { + if (!fChain) + return; + TIter Next(fChain->GetStatus()); @@ -870,8 +973,8 @@ } - const Bool_t rc = fChain->GetEntry(fNumEntry++) != 0; + const Bool_t rc = fTree->GetEntry(fNumEntry++) != 0; // *CHANGED-fChain-to-fTree* if (fTaskList) - fTaskList->SetStreamId(fChain->GetName()); + fTaskList->SetStreamId(fTree->GetName()); // *CHANGED-fChain-to-fTree* if (rc) @@ -887,4 +990,23 @@ Int_t MReadTree::PostProcess() { + // In the case of a memory tree I don't know how we can + // make a decision in PreProcess between a self allocated + // memory address or a pending address set long before. + // So we delete the stuff in PostProcess and not the destructor + // (which might seg faullt if PreProcess has never been called) + if (!fChain) + { + TIter Next(fTree->GetListOfBranches()); + TBranch *b=0; + while ((b=(TBranch*)Next())) + { + if (b->GetAddress()) + { + delete b->GetAddress(); + b->ResetAddress(); + } + } + } + return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE; } @@ -896,5 +1018,5 @@ Bool_t MReadTree::GetEvent() { - Bool_t rc = fChain->GetEntry(fNumEntry) != 0; + Bool_t rc = fTree->GetEntry(fNumEntry) != 0; // *CHANGED-fChain-to-fTree* if (rc) @@ -981,5 +1103,5 @@ TString MReadTree::GetFileName() const { - const TFile *file = fChain->GetFile(); + const TFile *file = fChain ? fChain->GetFile() : fTree->GetCurrentFile(); if (!file) @@ -997,5 +1119,5 @@ Int_t MReadTree::GetFileIndex() const { - return fChain->GetTreeNumber(); + return fChain ? fChain->GetTreeNumber() : 0; /* const TString filename = fChain->GetFile()->GetName(); Index: trunk/MagicSoft/Mars/mfileio/MReadTree.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadTree.h (revision 5691) +++ trunk/MagicSoft/Mars/mfileio/MReadTree.h (revision 5692) @@ -8,4 +8,5 @@ class MChain; class TBranch; +class TTree; class MTaskList; @@ -15,5 +16,6 @@ private: - MChain *fChain; // Pointer to tree + MChain *fChain; // Pointer to tree==fTree (only if fTree inherits from MChain) + TTree *fTree; // Pointer to tree UInt_t fNumEntry; // Number of actual entry in chain @@ -23,6 +25,6 @@ Bool_t fAutoEnable; // Flag for auto enabeling scheme - TList *fVetoList; // List of Branches which are not allowed to get enabled - TList *fNotify; // List of TObjects to notify when switching files + TList *fVetoList; //-> List of Branches which are not allowed to get enabled + TList *fNotify; //-> List of TObjects to notify when switching files MTaskList *fTaskList; // Tasklist to set StreamId @@ -47,5 +49,5 @@ public: - MReadTree(); + MReadTree(TTree *tree=0); MReadTree(const char *treename, const char *filename=NULL, const char *name=NULL, const char *title=NULL); ~MReadTree(); Index: trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc (revision 5691) +++ trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc (revision 5692) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 6/2001 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -37,4 +37,24 @@ // ReInit()) For more details se the corresponding constructor. // +// Memory based trees +// ------------------ +// It is possible to store the data into memory (TTrees) instead of +// writing the data into a file. To do this either call the default +// constructor or specify 'memory' as option in the constructor. +// +// Afterwards the tree can be found using gROOT->FindObject("treename") +// +// Currently(!) the tree is not deleted at all! Please make sure to +// delete it if it is not used anymore - otherwise you could wast a LOT +// of memory. Please consider that this behaviour might change in the +// future. +// +// Such trees are usefull if you want to use more basic root-tools +// like TMultiLayerPerceptron or TEventList. +// +// If you want to process such memory based Trees using Mars make sure, +// that you don't need data from the RunHeader tree because you can +// only use MReadTree but not MReadMarsFile with such a tree. +// ///////////////////////////////////////////////////////////////////////////// #include "MWriteRootFile.h" @@ -72,5 +92,5 @@ // - // Believing the root user guide, TTree instanced are owned by the + // Believing the root user guide, TTree instances are owned by the // directory (file) in which they are. This means we don't have to // care about their destruction. @@ -87,10 +107,4 @@ { Init(); - - // - // Set the Arrays the owner of its entries. This means, that the - // destructor of the arrays will delete all its entries. - // - fBranches.SetOwner(); } @@ -133,11 +147,23 @@ // MWriteRootFile::MWriteRootFile(const char *fname, - const Option_t *opt, + const Option_t *option, const char *ftitle, const Int_t comp, const char *name, - const char *title) + const char *title) : fOut(NULL) { Init(name, title); + + TString opt(option); + opt.ToLower(); + + // + // Check if we are writing to memory + // + if (opt.Contains("memory")) + { + fSplitRule = fname; + return; + } // @@ -172,16 +198,19 @@ Print(); - // - // If the file is still open (no error) write the keys. This is necessary - // for appearance of the all trees and branches. - // - if (IsFileOpen()) - fOut->Write(); - - // - // Delete the file. This'll also close the file (if open) - // - delete fOut; - fOut = 0; + if (fOut) + { + // + // If the file is still open (no error) write the keys. This is necessary + // for appearance of the all trees and branches. + // + if (IsFileOpen()) + fOut->Write(); + + // + // Delete the file. This'll also close the file (if open) + // + delete fOut; + fOut = 0; + } // @@ -257,4 +286,7 @@ void MWriteRootFile::AddContainer(const char *cname, const char *tname, Bool_t must) { + if (!fOut && !tname) + tname = fSplitRule; + TIter Next(&fBranches); TObject *o=0; @@ -262,4 +294,5 @@ if (TString(o->GetName())==TString(tname) && TString(o->GetTitle())==TString(cname)) { + *fLog << warn; *fLog << "WARNING - Container '" << cname << "' in Tree '" << tname << "' already scheduled... ignored." << endl; return; @@ -290,4 +323,7 @@ Bool_t must) { + if (!fOut && !tname) + tname = fSplitRule; + TIter Next(&fBranches); TObject *o=0; @@ -296,4 +332,5 @@ static_cast(o)->GetContainer()==cont) { + *fLog << warn; *fLog << "WARNING - Container " << cont << " in Tree '" << tname << "' already scheduled... ignored." << endl; return; @@ -378,7 +415,27 @@ // - // Check if the tree is already existing (part of the file) - // - TTree *tree = (TTree*)fOut->Get(tname); + // Check if the tree is already existing (part of the file or memory) + // + TTree *tree = fOut ? (TTree*)fOut->Get(tname) : dynamic_cast(gROOT->FindObject(tname)); + if (!fOut && tree) + { + if (tree->GetCurrentFile()) + { + *fLog << err; + *fLog << "ERROR - You are trying to write data into a memory stored root tree," << endl; + *fLog << " because you either called the default constructor or have" << endl; + *fLog << " instantiated MWriteRootFile using the write option 'memory'." << endl; + *fLog << " This tree '" << tname << "' is already existing in" << endl; + *fLog << " memory (gROOT->FindObject) and is already belonging to a" << endl; + *fLog << " file (" << tree->GetCurrentFile()->GetName() << ")." << endl; + *fLog << " This can - for example - happen if you are reading from a" << endl; + *fLog << " tree with the same name. The easiest solution in this case" << endl; + *fLog << " is to change the name of the tree you want to write to." << endl; + *fLog << endl; + return kFALSE; + } + *fLog << inf << "Tree '" << tname << "' found already in Memory... using." << endl; + } + if (!tree) { @@ -389,7 +446,10 @@ // TDirectory *save = gDirectory; - fOut->cd(); - - tree = new TTree(tname, ttitle); + if (fOut) + fOut->cd(); + else + gROOT->cd(); + + tree = new TTree(tname, ttitle, fOut ? 99 : 1); fTrees.AddLast(tree); @@ -567,4 +627,10 @@ // + // If we are writing into memory we don't split into seperate files + // + if (!fOut) + return kTRUE; + + // // For more information see TTree:ChangeFile() // @@ -774,8 +840,12 @@ // -------------------------------------------------------------------------- // -// return open state of the root file. +// return open state of the root file. If the file is 'memory' kTRUE is +// returned. // Bool_t MWriteRootFile::IsFileOpen() const { + if (!fOut) + return kTRUE; + const char *n = fOut->GetName(); return n==0 || *n==0 ? kTRUE : fOut->IsOpen(); @@ -784,8 +854,12 @@ // -------------------------------------------------------------------------- // -// return name of the root-file +// return name of the root-file. If the file is "memory" "" is +// returned. // const char *MWriteRootFile::GetFileName() const { + if (!fOut) + return ""; + const char *n = fOut->GetName(); return n==0 || *n==0 ? "" : n; @@ -794,9 +868,9 @@ // -------------------------------------------------------------------------- // -// cd into file. See TFile::cd() +// cd into file. See TFile::cd(). If the file is "memory" kTRUE is returned. // Bool_t MWriteRootFile::cd(const char *path) { - return fOut->cd(path); + return fOut ? fOut->cd(path) : kTRUE; } @@ -809,18 +883,20 @@ void MWriteRootFile::StreamPrimitive(ofstream &out) const { - out << " MWriteRootFile " << GetUniqueName() << "(\""; - out << fOut->GetName() << "\", \""; - out << fOut->GetOption() << "\", \""; - out << fOut->GetTitle() << "\", "; - out << fOut->GetCompressionLevel(); - - if (fName!=gsDefName || fTitle!=gsDefTitle) - { - out << ", \"" << fName << "\""; - if (fTitle!=gsDefTitle) - out << ", \"" << fTitle << "\""; - } - out << ");" << endl; - + out << " MWriteRootFile " << GetUniqueName(); + if (fOut) + { + out << "(\""; + out << fOut->GetName() << "\", \""; + out << fOut->GetOption() << "\", \""; + out << fOut->GetTitle() << "\", "; + out << fOut->GetCompressionLevel(); + out << ")"; + } + out << ";" << endl; + + if (fName!=gsDefName) + out << " " << GetUniqueName() << ".SetName(\"" << fName << "\");" << endl; + if (fTitle!=gsDefTitle) + out << " " << GetUniqueName() << ".SetTitle(\"" << fTitle << "\");" << endl; MRootFileBranch *entry; Index: trunk/MagicSoft/Mars/mhbase/MBinning.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 5691) +++ trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 5692) @@ -38,13 +38,20 @@ MBinning(const char *name=NULL, const char *title=NULL); MBinning(Int_t nbins, Axis_t lo, Axis_t hi, const char *name=0, const char *opt="", const char *title=NULL); + MBinning(const MBinning &bins) { SetEdges(bins); } + + void Copy(TObject &named) const + { + MBinning &bins = (MBinning&)named; + bins.SetEdges(*this); + } void SetEdges(const TArrayD &arr) { fEdges = arr; - fType = kIsUserArray; + fType = kIsUserArray; } void SetEdges(const TAxis &axe); - void SetEdges(const MBinning &bins) { SetEdges(fEdges); } + void SetEdges(const MBinning &bins) { SetEdges(bins.fEdges); fType = bins.fType; } void SetEdges(const TH1 &h, const Char_t axis='x'); void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up); Index: trunk/MagicSoft/Mars/mhbase/MHMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MHMatrix.cc (revision 5691) +++ trunk/MagicSoft/Mars/mhbase/MHMatrix.cc (revision 5692) @@ -160,4 +160,8 @@ } + const Int_t idx = fData->FindRule(rule); + if (idx>=0) + return idx; + fData->AddEntry(rule); return fData->GetNumEntries()-1; Index: trunk/MagicSoft/Mars/mhbase/MHMatrix.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MHMatrix.h (revision 5691) +++ trunk/MagicSoft/Mars/mhbase/MHMatrix.h (revision 5692) @@ -115,5 +115,5 @@ void ReduceRows(UInt_t num); - ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events + ClassDef(MHMatrix, 1) // Multidimensional Matrix (TMatrix) to store events }; Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5691) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5692) @@ -41,4 +41,7 @@ #include #include +#include + +#include #include @@ -51,4 +54,18 @@ using namespace std; + +void MAlphaFitter::Clear(Option_t *o) +{ + fSignificance=0; + fEventsExcess=0; + fEventsSignal=0; + fEventsBackground=0; + + fChiSqSignal=0; + fChiSqBg=0; + fIntegralMax=0; + + fCoefficients.Reset(); +} // -------------------------------------------------------------------------- @@ -84,4 +101,8 @@ Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint) { + Clear(); + if (h.GetEntries()==0) + return kFALSE; + Double_t sigmax=fSigMax; Double_t bgmin =fBgMin; @@ -175,5 +196,4 @@ //fSignificance = MMath::SignificanceLiMaSigned(s, b); - const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999); @@ -190,10 +210,90 @@ } +Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Bool_t paint) +{ + /* + Clear(); + if (hon.GetEntries()==0) + return kFALSE; + */ + + TH1D h(hon); + h.Add(&hof, -1); + + MAlphaFitter fit(*this); + fit.SetPolynomOrder(1); + + if (!fit.Fit(h, paint)) + return kFALSE; + + fChiSqSignal = fit.GetChiSqSignal(); + fChiSqBg = fit.GetChiSqBg(); + fCoefficients = fit.GetCoefficients(); + + const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999); + + fIntegralMax = hon.GetBinLowEdge(bin+1); + fEventsBackground = hof.Integral(0, bin); + fEventsSignal = hon.Integral(0, bin); + fEventsExcess = fEventsSignal-fEventsBackground; + fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground); + + if (fEventsExcess<0) + fEventsExcess=0; +/* + TF1 func("", "gaus(0)+pol0(3)", 0., 90.); + + const Double_t A = fEventsSignal/bin; + const Double_t dA = TMath::Abs(A); + func.SetParLimits(0, -dA*4, dA*4); + func.SetParLimits(2, 0, 90); + func.SetParLimits(3, -dA, dA); + + func.SetParameter(0, A); + func.FixParameter(1, 0); + func.SetParameter(2, fSigMax*0.75); + func.SetParameter(3, 0); + + // options : N do not store the function, do not draw + // I use integral of function in bin rather than value at bin center + // R use the range specified in the function range + // Q quiet mode + TH1D h(hon); + h.Add(&hof, -1); + h.Fit(&func, "NQI", "", 0, 90); + + fChiSqSignal = func.GetChisquare()/func.GetNDF(); + + const Int_t bin1 = h.GetXaxis()->FindFixBin(func.GetParameter(2)*2); + + fChiSqBg = 0; + for (int i=bin1; i<=h.GetNbinsX(); i++) + { + const Float_t val = h.GetBinContent(i); + fChiSqBg = val*val; + } + if (fChiSqBg>0) + fChiSqBg /= h.GetNbinsX()+1-bin1; + + fCoefficients.Set(func.GetNpar(), func.GetParameters()); + + // ------------------------------------ + if (paint) + { + func.SetLineColor(kBlue); + func.SetLineWidth(2); + func.Paint("same"); + } + // ------------------------------------ + */ + return kTRUE; +} + void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const { - TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}=%.1f \\chi_{s}^{2}=%.1f)", + TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}/ndf=%.1f \\chi_{s}^{2}/ndf=%.1f c_{0}=%.1f)", fSignificance, GetGausSigma(), (int)fEventsExcess, fIntegralMax, - fChiSqBg, fChiSqSignal)); + fChiSqBg, fChiSqSignal, fCoefficients[3])); text.SetBit(TLatex::kTextNDC); @@ -210,5 +310,9 @@ f.fBgMin = fBgMin; f.fBgMax = fBgMax; + f.fScaleMin = fScaleMin; + f.fScaleMax = fScaleMax; f.fPolynomOrder = fPolynomOrder; + f.fScaleMode = fScaleMode; + f.fScaleUser = fScaleUser; f.fCoefficients.Set(fCoefficients.GetSize()); f.fCoefficients.Reset(); @@ -227,3 +331,158 @@ *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl; *fLog << " ...polynom order " << fPolynomOrder << endl; -} + *fLog << " ...scale mode: "; + switch (fScaleMode) + { + case kNone: *fLog << "none."; break; + case kEntries: *fLog << "entries."; break; + case kIntegral: *fLog << "integral."; break; + case kOffRegion: *fLog << "off region."; break; + case kLeastSquare: *fLog << "least square."; break; + case kUserScale: *fLog << "user def (" << fScaleUser << ")"; break; + } + *fLog << endl; + + if (TString(o).Contains("result")) + { + *fLog << "Result:" << endl; + *fLog << " - Significance " << fSignificance << endl; + *fLog << " - Excess Events " << fEventsExcess << endl; + *fLog << " - Signal Events " << fEventsSignal << endl; + *fLog << " - Background Events " << fEventsBackground << endl; + *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl; + *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl; + *fLog << " - Signal integrated " << fIntegralMax << "°" << endl; + } +} + +Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint) +{ + const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000))); + TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E"); + h->SetDirectory(0); + + const Bool_t rc = Fit(*h, paint); + delete h; + return rc; +} + +Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint) +{ + const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000))); + TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E"); + h->SetDirectory(0); + + const Bool_t rc = Fit(*h, paint); + delete h; + return rc; +} + +Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint) +{ + const TString name(Form("TempAlpha%06d", gRandom->Integer(1000000))); + TH1D *h = hon.ProjectionZ(name, -1, 9999, -1, 9999, "E"); + h->SetDirectory(0); + + const Bool_t rc = Fit(*h, paint); + delete h; + return rc; +} + +Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint) +{ + const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000))); + const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000))); + TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E"); + TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E"); + h1->SetDirectory(0); + h0->SetDirectory(0); + + Scale(*h0, *h1); + + const Bool_t rc = Fit(*h1, *h0, paint); + + delete h0; + delete h1; + + return rc; +} + +Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint) +{ + const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000))); + const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000))); + + TH1D *h1 = hon.ProjectionZ(name1, bin, bin, -1, 9999, "E"); + TH1D *h0 = hof.ProjectionZ(name0, bin, bin, -1, 9999, "E"); + h1->SetDirectory(0); + h0->SetDirectory(0); + + Scale(*h0, *h1); + + const Bool_t rc = Fit(*h1, *h0, paint); + + delete h0; + delete h1; + + return rc; +} + +Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint) +{ + const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000))); + const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000))); + + TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E"); + TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E"); + h1->SetDirectory(0); + h0->SetDirectory(0); + + Scale(*h0, *h1); + + const Bool_t rc = Fit(*h1, *h0, paint); + + delete h0; + delete h1; + + return rc; +} + +void MAlphaFitter::Scale(TH1D &of, const TH1D &on) const +{ + Float_t scaleon = 1; + Float_t scaleof = 1; + switch (fScaleMode) + { + case kNone: + return; + + case kEntries: + scaleon = on.GetEntries(); + scaleof = of.GetEntries(); + break; + + case kIntegral: + scaleon = on.Integral(); + scaleof = of.Integral(); + break; + + case kOffRegion: + { + const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin); + const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax); + scaleon = on.Integral(min, max); + scaleof = of.Integral(min, max); + } + break; + + case kUserScale: + scaleon = fScaleUser; + break; + + default: + return; + } + + if (scaleof!=0) + of.Scale(scaleon/scaleof); +} Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5691) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5692) @@ -15,7 +15,17 @@ class TH1D; +class TH3D; class MAlphaFitter : public MParContainer { +public: + enum ScaleMode_t { + kNone, + kEntries, + kIntegral, + kOffRegion, + kLeastSquare, + kUserScale + }; private: Float_t fSigInt; @@ -23,4 +33,6 @@ Float_t fBgMin; Float_t fBgMax; + Float_t fScaleMin; + Float_t fScaleMax; Int_t fPolynomOrder; @@ -38,7 +50,10 @@ TF1 *fFunc; + ScaleMode_t fScaleMode; + Double_t fScaleUser; + public: // Implementing the function yourself is only about 5% faster - MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)) + MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kEntries), fScaleUser(1) { fName = name ? name : "MAlphaFitter"; @@ -58,5 +73,6 @@ } - void Print(Option_t *o=0) const; + void Clear(Option_t *o=""); + void Print(Option_t *o="") const; void Copy(TObject &o) const; /* @@ -67,8 +83,12 @@ */ + void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; } + void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; } void SetSignalIntegralMax(Float_t s) { fSigInt = s; } void SetSignalFitMax(Float_t s) { fSigMax = s; } void SetBackgroundFitMin(Float_t s) { fBgMin = s; } void SetBackgroundFitMax(Float_t s) { fBgMax = s; } + void SetScaleMin(Float_t s) { fScaleMin = s; } + void SetScaleMax(Float_t s) { fScaleMax = s; } void SetPolynomOrder(Int_t s) { fPolynomOrder = s; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", s)); gROOT->GetListOfFunctions()->Remove(fFunc); @@ -93,4 +113,32 @@ Bool_t Fit(TH1D &h, Bool_t paint=kFALSE); + Bool_t Fit(TH1D &on, TH1D &off, Bool_t paint=kFALSE); + Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE) + { + return off ? Fit(on, *off, paint) : Fit(on, paint); + } + + Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE); + Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE); + Bool_t FitTheta(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE); + + Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE); + Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE); + Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE); + + Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE) + { + return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint); + } + Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE) + { + return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint); + } + Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE) + { + return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint); + } + + void Scale(TH1D &off, const TH1D &on) const; ClassDef(MAlphaFitter, 1) Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5691) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5692) @@ -75,6 +75,9 @@ // MHAlpha::MHAlpha(const char *name, const char *title) - : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0), - fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0) + : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0), + fPointPos(0), fTimeEffOn(0), fTime(0), + fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE), + //fEnergyMin(-1), fEnergyMax(-1), fSizeMin(-1), fSizeMax(-1), + fMatrix(0) { // @@ -100,7 +103,7 @@ - fHEnergy.SetName("Energy"); - fHEnergy.SetTitle(" N_{exc} vs. E_{est} "); - fHEnergy.SetXTitle("E_{est} [GeV]"); + //fHEnergy.SetName("Energy"); + //fHEnergy.SetTitle(" N_{exc} vs. E_{est} "); + //fHEnergy.SetXTitle("E_{est} [GeV]"); fHEnergy.SetYTitle("N_{exc}"); fHEnergy.SetDirectory(NULL); @@ -184,11 +187,9 @@ for (int i=1; i<=n; i++) { - TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":""); - if (fit.Fit(*h)) + if (fit.FitEnergy(fHAlpha, fOffData, i)) { fHEnergy.SetBinContent(i, fit.GetEventsExcess()); fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2); } - delete h; } } @@ -203,11 +204,9 @@ for (int i=1; i<=n; i++) { - TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":""); - if (fit.Fit(*h)) + if (fit.FitTheta(fHAlpha, fOffData, i)) { fHTheta.SetBinContent(i, fit.GetEventsExcess()); fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2); } - delete h; } } @@ -219,7 +218,56 @@ fHTheta.Reset(); + if (fName!=(TString)"MHAlphaOff") + { + MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff"); + if (!hoff) + *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl; + else + { + *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl; + SetOffData(*hoff); + } + } + + fHillas = 0; + /* + if (fSizeMin>=0 || fSizeMax>=0) + { + fHillas = (MHillas*)pl->FindObject("MHillas"); + if (!fHillas) + { + *fLog << warn << "Size cut set, but MHillas not found... abort." << endl; + return kFALSE; + } + } + */ fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst"); if (!fEnergy) - *fLog << warn << "MEnergyEst not found... ignored." << endl; + { /* + if (fEnergyMin>=0 || fEnergyMax>=0) + { + *fLog << warn << "Energy cut set, but MEnergyEst not found... abort." << endl; + return kFALSE; + } */ + + *fLog << warn << "MEnergyEst not found... " << flush; + + if (!fHillas) + fHillas = (MHillas*)pl->FindObject("MHillas"); + if (fHillas) + *fLog << "using SIZE instead!" << endl; + else + *fLog << "ignored." << endl; + + fHEnergy.SetName("Size"); + fHEnergy.SetTitle(" N_{exc} vs. Size "); + fHEnergy.SetXTitle("Size [\\gamma]"); + } + else + { + fHEnergy.SetName("Energy"); + fHEnergy.SetTitle(" N_{exc} vs. E_{est} "); + fHEnergy.SetXTitle("E_{est} [GeV]"); + } fPointPos = (MPointingPos*)pl->FindObject("MPointingPos"); @@ -230,4 +278,6 @@ if (!fTimeEffOn) *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl; + else + *fTimeEffOn = MTime(); // FIXME: How to do it different? fTime = (MTime*)pl->FindObject("MTime"); @@ -246,23 +296,26 @@ MBinning binst, binse, binsa; - binst.SetEdges(fHAlpha, 'x'); - binse.SetEdges(fHAlpha, 'y'); - binsa.SetEdges(fHAlpha, 'z'); - - MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning"); - if (fPointPos && bins) - binst.SetEdges(*bins); - if (!fPointPos) - binst.SetEdges(1, 0, 90); - - bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning"); - if (fEnergy && bins) - binse.SetEdges(*bins); - if (!fEnergy) - binse.SetEdges(1, 10, 100000); - - bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning"); - if (bins) - binsa.SetEdges(*bins); + binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x'); + binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y'); + binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z'); + + if (!fOffData) + { + MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning"); + if (fPointPos && bins) + binst.SetEdges(*bins); + if (!fPointPos) + binst.SetEdges(1, 0, 60); + + bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning"); + if ((fEnergy||fHillas) && bins) + binse.SetEdges(*bins); + if (!fEnergy && !fHillas) + binse.SetEdges(1, 10, 100000); + + bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning"); + if (bins) + binsa.SetEdges(*bins); + } binse.Apply(fHEnergy); @@ -321,6 +374,26 @@ } + // Check new 'last time' + MTime *time = final ? fTime : fTimeEffOn; + + if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax()) + { + *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl; + *fLog << "than upper edge of histogram... skipped." << endl; + *fLog << "This should not happen. Maybe you started you eventloop with" << endl; + *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl; + rebin++; + return; + } + + // Fit time histogram MAlphaFitter fit(fFit); - if (!fit.Fit(fHAlphaTime)) + + TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0; + const Bool_t rc = fit.Fit(fHAlphaTime, h); + if (h) + delete h; + + if (!rc) return; @@ -334,6 +407,4 @@ // Prepare Histogram // - - MTime *time = final ? fTime : fTimeEffOn; if (final) time->Plus1ns(); @@ -363,10 +434,21 @@ { Double_t alpha, energy, theta; + Double_t size=-1; if (fMatrix) { alpha = (*fMatrix)[fMap[0]]; - energy = 1000; - theta = 0; + energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]]; + size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]]; + //<0 ? 1000 : (*fMatrix)[fMap[1]]; + theta = 0; + + if (energy<0) + energy=size; + if (size<0) + size=energy; + + if (energy<0 && size<0) + energy = size = 1000; } else @@ -380,12 +462,27 @@ alpha = hil->GetAlpha(); - energy = fEnergy ? fEnergy->GetEnergy() : 1000; + if (fHillas) + size = fHillas->GetSize(); + energy = fEnergy ? fEnergy->GetEnergy() : (fHillas?fHillas->GetSize():1000); theta = fPointPos ? fPointPos->GetZd() : 0; } + // enhance histogram if necessary UpdateAlphaTime(); + // Fill histograms fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w); - fHAlphaTime.Fill(TMath::Abs(alpha), w); + + // Check cuts + /* + if ( (fEnergyMin>=0 && energy=0 && energy>fEnergyMax) || + (fSizeMin>=0 && size =0 && size >fSizeMin) ) + return kTRUE; + */ + + if (!fSkipHistTime) + fHAlphaTime.Fill(TMath::Abs(alpha), w); return kTRUE; @@ -441,5 +538,21 @@ padsave->cd(1); - fHAlpha.ProjectionZ(fNameProjAlpha); + TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha"); + if (fOffData) + { + TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff"); + fFit.Scale(*hoff, *hon); + + hon->SetMaximum(); + hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05); + + if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff"))) + { + h0->Reset(); + h0->Add(hoff, hon, -1); + const Float_t min = h0->GetMinimum()*1.05; + hon->SetMinimum(min<0 ? min : 0); + } + } FitEnergyBins(); FitThetaBins(); @@ -447,9 +560,10 @@ if (o==(TString)"alpha") - if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha))) + if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha"))) { // Do not store the 'final' result in fFit MAlphaFitter fit(fFit); - fit.Fit(*h0, kTRUE); + TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff"); + fit.Fit(*h0, hoff, kTRUE); fit.PaintResult(); } @@ -459,8 +573,31 @@ if (o==(TString)"theta") + { + TH1 *h = (TH1*)gPad->FindObject("Alpha_x"); + if (h) + { + TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x"); + h2->SetDirectory(0); + h2->Scale(fHTheta.Integral()/h2->Integral()); + h->Reset(); + h->Add(h2); + delete h2; + } PaintText(fHTheta.Integral(), 0); + } if (o==(TString)"energy") { + TH1 *h = (TH1*)gPad->FindObject("Alpha_y"); + if (h) + { + TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y"); + h2->SetDirectory(0); + h2->Scale(fHEnergy.Integral()/h2->Integral()); + h->Reset(); + h->Add(h2); + delete h2; + } + if (fHEnergy.GetMaximum()>1) { @@ -469,5 +606,4 @@ } FitEnergySpec(kTRUE); - } @@ -494,6 +630,8 @@ pad->cd(1); gPad->SetBorderMode(0); - h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E"); + + h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E"); h->SetBit(TH1::kNoTitle); + h->SetStats(kTRUE); h->SetXTitle("\\alpha [\\circ]"); h->SetYTitle("Counts"); @@ -501,5 +639,35 @@ h->SetMarkerStyle(kFullDotMedium); h->SetBit(kCanDelete); - h->Draw(); + h->Draw(""); + + if (fOffData) + { + h->SetMarkerColor(kGreen); + + h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E"); + h->SetBit(TH1::kNoTitle|TH1::kNoStats); + h->SetXTitle("\\alpha [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->SetMarkerColor(kRed); + h->Draw("same"); + + h = (TH1D*)h->Clone("ProjAlphaOnOff"); + h->SetBit(TH1::kNoTitle); + h->SetXTitle("\\alpha [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->SetMarkerColor(kBlue); + h->Draw("same"); + + TLine lin; + lin.SetLineStyle(kDashed); + lin.DrawLine(0, 0, 90, 0); + } + // After the Alpha-projection has been drawn. Fit the histogram // and paint the result into this pad @@ -511,5 +679,17 @@ gPad->SetBorderMode(0); fHEnergy.Draw(); + AppendPad("energy"); + + h = (TH1D*)fHAlpha.Project3D("y"); + h->SetBit(TH1::kNoTitle|TH1::kNoStats); + h->SetXTitle("E [GeV]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->SetMarkerColor(kCyan); + h->SetLineColor(kCyan); + h->Draw("Psame"); } else @@ -531,20 +711,107 @@ gPad->SetBorderMode(0); fHTheta.Draw(); + AppendPad("theta"); + + h = (TH1D*)fHAlpha.Project3D("x"); + h->SetBit(TH1::kNoTitle|TH1::kNoStats); + h->SetXTitle("\\theta [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->SetMarkerColor(kCyan); + h->SetLineColor(kCyan); + h->Draw("Psame"); } else delete pad->GetPad(4); - -} +} + +void MHAlpha::DrawAll() +{ + // FIXME: Do in Paint if special option given! + TCanvas *c = new TCanvas; + Int_t n = fHAlpha.GetNbinsY(); + Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1); + c->Divide(nc, nc, 0, 0); + + // Do not store the 'final' result in fFit + MAlphaFitter fit(fFit); + + for (int i=1; i<=fHAlpha.GetNbinsY(); i++) + { + c->cd(i); + + TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E"); + hon->SetBit(TH1::kNoTitle); + hon->SetStats(kTRUE); + hon->SetXTitle("\\alpha [\\circ]"); + hon->SetYTitle("Counts"); + hon->SetDirectory(NULL); + hon->SetMarkerStyle(kFullDotMedium); + hon->SetBit(kCanDelete); + hon->Draw(""); + + TH1D *hof = 0; + + if (fOffData) + { + hon->SetMarkerColor(kGreen); + + hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E"); + hof->SetBit(TH1::kNoTitle|TH1::kNoStats); + hof->SetXTitle("\\alpha [\\circ]"); + hof->SetYTitle("Counts"); + hof->SetDirectory(NULL); + hof->SetMarkerStyle(kFullDotMedium); + hof->SetBit(kCanDelete); + hof->SetMarkerColor(kRed); + hof->Draw("same"); + + fit.Scale(*hof, *hon); + + hon->SetMaximum(); + hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05); + + TH1D *diff = new TH1D(*hon); + diff->Add(hof, -1); + diff->SetBit(TH1::kNoTitle); + diff->SetXTitle("\\alpha [\\circ]"); + diff->SetYTitle("Counts"); + diff->SetDirectory(NULL); + diff->SetMarkerStyle(kFullDotMedium); + diff->SetBit(kCanDelete); + diff->SetMarkerColor(kBlue); + diff->Draw("same"); + + TLine lin; + lin.SetLineStyle(kDashed); + lin.DrawLine(0, 0, 90, 0); + + const Float_t min = diff->GetMinimum()*1.05; + hon->SetMinimum(min<0 ? min : 0); + } + + if (hof ? fit.Fit(*hon, *hof) : fit.Fit(*hon)) + *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl; + /* + if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE)) + { + fHEnergy.SetBinContent(i, fit.GetEventsExcess()); + fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2); + }*/ + } + +} + Bool_t MHAlpha::Finalize() { - // Store the final result in fFit - TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E"); - h->SetDirectory(0); - fFit.Print(); - Bool_t rc = fFit.Fit(*h); - delete h; - if (!rc) + //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E"); + //h->SetDirectory(0); + //Bool_t rc = fFit.Fit(*h); + //delete h; + if (!fFit.FitAlpha(fHAlpha, fOffData)) { *fLog << warn << "Histogram empty." << endl; @@ -552,11 +819,26 @@ } + // Store the final result in fFit + fFit.Print("result"); + if (fResult) fResult->SetVal(fFit.GetSignificance()); - FitEnergyBins(); - FitThetaBins(); - UpdateAlphaTime(kTRUE); - MH::RemoveFirstBin(fHTime); + if (!fSkipHistEnergy) + { + *fLog << inf << "Processing energy bins..." << endl; + FitEnergyBins(); + } + if (!fSkipHistTheta) + { + *fLog << inf << "Processing theta bins..." << endl; + FitThetaBins(); + } + if (!fSkipHistTime) + { + *fLog << inf << "Processing time bins..." << endl; + UpdateAlphaTime(kTRUE); + MH::RemoveFirstBin(fHTime); + } return kTRUE; @@ -572,5 +854,7 @@ // will take the values from the matrix instead of the containers. // -void MHAlpha::InitMapping(MHMatrix *mat) +// It takes fSkipHist* into account! +// +void MHAlpha::InitMapping(MHMatrix *mat, Int_t type) { if (fMatrix) @@ -580,5 +864,26 @@ fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha"); - //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy"); + fMap[1] = -1; + fMap[2] = -1; + fMap[3] = -1; + fMap[4] = -1; + + if (!fSkipHistEnergy) + if (type==0) + { + fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy"); + fMap[2] = -1; + } + else + { + fMap[1] = -1; + fMap[2] = fMatrix->AddColumn("MHillas.fSize"); + } + + if (!fSkipHistTheta) + fMap[3] = fMatrix->AddColumn("MPointingPos.fZd"); + + // if (!fSkipHistTime) + // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime"); } Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 5691) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 5692) @@ -21,4 +21,5 @@ class MParameterD; class MEnergyEst; +class MHillas; class MHMatrix; class MPointingPos; @@ -84,4 +85,6 @@ { private: + const TH3D *fOffData; + MAlphaFitter fFit; // SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT) @@ -94,4 +97,5 @@ MParameterD *fResult; //! MEnergyEst *fEnergy; //! + MHillas *fHillas; //! MPointingPos *fPointPos; //! @@ -100,8 +104,17 @@ MTime fLastTime; //! Last fTimeEffOn - const TString fNameProjAlpha; //! This should make sure, that gROOT doen't confuse the projection with something else + //Float_t fEnergyMin; + //Float_t fEnergyMax; + //Float_t fSizeMin; + //Float_t fSizeMax; + + Bool_t fSkipHistTime; + Bool_t fSkipHistTheta; + Bool_t fSkipHistEnergy; + + //const TString fNameProjAlpha; //! This should make sure, that gROOT doen't confuse the projection with something else MHMatrix *fMatrix; //! - Int_t fMap[2]; //! + Int_t fMap[5]; //! void FitEnergySpec(Bool_t paint=kFALSE); @@ -115,4 +128,6 @@ void PaintText(Double_t val, Double_t error) const; + Int_t DistancetoPrimitive(Int_t px, Int_t py) { return 0; } + public: MHAlpha(const char *name=NULL, const char *title=NULL); @@ -125,8 +140,32 @@ const MAlphaFitter &GetAlphaFitter() const { return fFit; } + void SetOffData(const MHAlpha &h) + { + fOffData = &h.fHAlpha; + } +/* + void SetSizeCuts(Float_t min, Float_t max) { fSizeMin=min; fSizeMax=max; } + void SetSizeMin(Float_t min) { fSizeMin=min; } + void SetSizeMax(Float_t max) { fSizeMax=max; } + void SetEnergyCuts(Float_t min, Float_t max) { fEnergyMin=min; fEnergyMax=max; } + void SetEnergyMin(Float_t min) { fEnergyMin=min; } + void SetEnergyMax(Float_t max) { fEnergyMax=max; } + + void SetCuts(const MHAlpha &h) { + fSizeMin = h.fSizeMin; fEnergyMin = h.fEnergyMin; + fSizeMax = h.fSizeMax; fEnergyMax = h.fEnergyMax; + } + */ + + void SkipHistTime(Bool_t b=kTRUE) { fSkipHistTime=b; } + void SkipHistTheta(Bool_t b=kTRUE) { fSkipHistTheta=b; } + void SkipHistEnergy(Bool_t b=kTRUE) { fSkipHistEnergy=b; } + void Paint(Option_t *opt=""); void Draw(Option_t *option=""); - void InitMapping(MHMatrix *mat); + void DrawAll(); //*MENU* + + void InitMapping(MHMatrix *mat, Int_t type=0); void StopMapping(); Index: trunk/MagicSoft/Mars/mjobs/MSequence.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MSequence.cc (revision 5691) +++ trunk/MagicSoft/Mars/mjobs/MSequence.cc (revision 5692) @@ -145,4 +145,6 @@ d += fNight.GetStringFmt("%Y_%m_%d"); } + else + gSystem->ExpandPathName(d); for (int i=0; iExpandPathName(fname); + + fTitle = Form("Sequence contained in file %s", expname); + + TEnv env(expname); + delete [] expname; TString str;