Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7129) +++ trunk/MagicSoft/Mars/Changelog (revision 7130) @@ -23,4 +23,52 @@ 2005/06/03 Thomas Bretz + + * mjobs/hilocalib_sp1.root, mjobs/hilocalib_sp1_mc.root: + - updated + + * callisto.rc: + - MJPedestalY2.MaxEvents: 2000 replaced by 5000 as in + callisto_Dec04Jan05.txt + + * manalysis/MMultiDimDistCalc.h: + - changes to layout + + * mbadpixels/MBadPixelsCalc.cc: + - improved sanity checks + + * mbase/MEvtLoop.cc: + - fixed a bug which could cause Eventloop to crahs if + parlist was not initialized + + * mdata/MDataArray.[h,cc]: + - added copy constructor + + * mhbase/MFillH.cc: + - made sure that no constructor can crash due to NULL pointers + + * mhbase/MHMatrix.[h,cc]: + - first check in AddColumn if the column is available. Afterwards + check whether it can be added + - added new interfaced to single rows + + * mhflux/MMcSpectrumWeight.cc: + - slight change to screen output + + * mjobs/MJPedestal.cc: + - slight change to screen output + + * mpedestal/MPedCalcPedRun.cc: + - fixed a bug which caused MC files not to work treat them now + as pedestal files (always) + + * mranforest/MRFEnergyEst.[h,cc]: + - improved the code and the interface + + * mranforest/MRanForestGrow.[h,cc]: + - derives now from MRead to be able to use the bar in the + display + + * mtools/MTFillMatrix.[h,cc]: + - allow to fill a single matrix with all events * datacenter/macros/plotdb.C: Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 7129) +++ trunk/MagicSoft/Mars/NEWS (revision 7130) @@ -59,4 +59,10 @@ - callisto: In the resource file callisto_Dec04Jan05.rc MJPedestalY2.ExtractWinRight has been reduced from 4.0 to 2.0 + + - callisto: new Hi-/Lo-Gain intercalibration constants + hilocalib_sp1.root and hilocalib_sp1_mc.root + + - callisto: changed default for MJPedestalY2.MaxEvents + from 2000 to 5000 like in callisto_Dec04Jan05.txt - callisto: in MCalibrationChargeCalc the limit fgPheErrLowerLimit Index: trunk/MagicSoft/Mars/callisto.rc =================================================================== --- trunk/MagicSoft/Mars/callisto.rc (revision 7129) +++ trunk/MagicSoft/Mars/callisto.rc (revision 7130) @@ -317,5 +317,5 @@ #MJPedestalY3.UseData: Yes MJPedestalY1.MaxEvents: 500 -MJPedestalY2.MaxEvents: 2000 +MJPedestalY2.MaxEvents: 5000 MJPedestalY3.MaxEvents: 500 Index: trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h (revision 7129) +++ trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h (revision 7130) @@ -19,12 +19,12 @@ TString fHadronnessName; // Name of container storing hadronness - MHMatrix *fMGammas; //! Gammas describing matrix - MHMatrix *fMHadrons; //! Hadrons (non gammas) describing matrix + MHMatrix *fMGammas; //! Gammas describing matrix + MHMatrix *fMHadrons; //! Hadrons (non gammas) describing matrix MParameterD *fHadronness; //! Output container for calculated hadroness - MDataArray *fData; //! Used to store the MDataChains to get the event values + MDataArray *fData; //! Used to store the MDataChains to get the event values - void StreamPrimitive(ofstream &out) const; + void StreamPrimitive(ofstream &out) const; Int_t PreProcess(MParList *plist); Int_t Process(); Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 7129) +++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 7130) @@ -139,16 +139,16 @@ Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const { - if (!fGeomCam || !fPedPhotCam || !fBadPixels) - { - *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - One of the necessary container are not initialized..." << endl; - return kFALSE; - } - - const Bool_t checklo = fPedestalLevelVarianceLo>0; - const Bool_t checkhi = fPedestalLevelVarianceHi>0; + const Bool_t checklo = fPedestalLevelVarianceLo>0; + const Bool_t checkhi = fPedestalLevelVarianceHi>0; if (fPedestalLevel<=0 && !checklo && !checkhi) return kTRUE; + if (!fGeomCam || !fPedPhotCam || !fBadPixels) + { + *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - One of the necessary container are not initialized..." << endl; + return kFALSE; + } + const Int_t entries = fPedPhotCam->GetSize(); @@ -175,14 +175,21 @@ //if no pixel has a minimum signal, return + Int_t counter=0; for (int i=0; iGetNumExecutions()>0) + if (fTaskList && fTaskList->GetNumExecutions()>0) switch (printstat) { Index: trunk/MagicSoft/Mars/mdata/MDataArray.cc =================================================================== --- trunk/MagicSoft/Mars/mdata/MDataArray.cc (revision 7129) +++ trunk/MagicSoft/Mars/mdata/MDataArray.cc (revision 7130) @@ -51,5 +51,5 @@ // -------------------------------------------------------------------------- // -// Constructor +// Default Constructor // MDataArray::MDataArray(const char *name, const char *title) @@ -60,4 +60,22 @@ gROOT->GetListOfCleanups()->Add(&fList); fList.SetBit(kMustCleanup); +} + +// -------------------------------------------------------------------------- +// +// Copy Constructor +// +MDataArray::MDataArray(MDataArray &a, const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + gROOT->GetListOfCleanups()->Add(&fList); + fList.SetBit(kMustCleanup); + + TIter Next(&a.fList); + MData *data = NULL; + while ((data=(MData*)Next())) + AddEntry(data->GetRule()); } Index: trunk/MagicSoft/Mars/mdata/MDataArray.h =================================================================== --- trunk/MagicSoft/Mars/mdata/MDataArray.h (revision 7129) +++ trunk/MagicSoft/Mars/mdata/MDataArray.h (revision 7130) @@ -28,4 +28,5 @@ public: MDataArray(const char *name=NULL, const char *title=NULL); + MDataArray(MDataArray &a, const char *name=NULL, const char *title=NULL); void AddEntry(const TString rule); Index: trunk/MagicSoft/Mars/mhbase/MFillH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MFillH.cc (revision 7129) +++ trunk/MagicSoft/Mars/mhbase/MFillH.cc (revision 7130) @@ -203,11 +203,13 @@ fHName = hist; fParContainer = par; - fParContainerName = par->GetName(); + if (par) + fParContainerName = par->GetName(); AddToBranchList(Form("%s.*", (const char*)ExtractName(hist))); - AddToBranchList(Form("%s.*", par->GetName())); + if (par) + AddToBranchList(Form("%s.*", par->GetName())); if (!title) - fTitle = Form("Fill %s from %s", fName.Data(), par->GetDescriptor().Data()); + fTitle = Form("Fill %s from %s", fName.Data(), par?par->GetDescriptor().Data():"NULL"); } @@ -226,8 +228,10 @@ fH = hist; - fHName = hist->GetName(); + if (hist) + fHName = hist->GetName(); fParContainerName = par; - AddToBranchList(fH->GetDataMember()); + if (fH) + AddToBranchList(fH->GetDataMember()); if (par) AddToBranchList(Form("%s.*", (const char*)ExtractName(par))); @@ -236,5 +240,5 @@ return; - fTitle = Form("Fill %s", hist->GetDescriptor().Data()); + fTitle = Form("Fill %s", hist ? hist->GetDescriptor().Data() : "NULL"); if (!par) return; @@ -257,13 +261,19 @@ fH = hist; - fHName = hist->GetName(); + if (hist) + fHName = hist->GetName(); fParContainer = par; - fParContainerName = par->GetName(); - - AddToBranchList(fH->GetDataMember()); - AddToBranchList(Form("%s.*", par->GetName())); + if (par) + fParContainerName = par->GetName(); + + if (fH) + AddToBranchList(fH->GetDataMember()); + if (par) + AddToBranchList(Form("%s.*", par->GetName())); if (!title) - fTitle = Form("Fill %s from %s", hist->GetDescriptor().Data(), par->GetDescriptor().Data()); + fTitle = Form("Fill %s from %s", + hist?hist->GetDescriptor().Data():"NULL", + par ? par->GetDescriptor().Data():"NULL"); } Index: trunk/MagicSoft/Mars/mhbase/MHMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MHMatrix.cc (revision 7129) +++ trunk/MagicSoft/Mars/mhbase/MHMatrix.cc (revision 7130) @@ -143,4 +143,8 @@ Int_t MHMatrix::AddColumn(const char *rule) { + const Int_t idx = fData ? fData->FindRule(rule) : -1; + if (idx>=0) + return idx; + if (IsValid(fM)) { @@ -160,8 +164,4 @@ SetBit(kIsOwner); } - - const Int_t idx = fData->FindRule(rule); - if (idx>=0) - return idx; fData->AddEntry(rule); @@ -1262,2 +1262,16 @@ return kTRUE; } + +// -------------------------------------------------------------------------- +// +// Returns the row pointed to by fNumRow into TVector v +// +void MHMatrix::GetRow(TVector &v) const +{ + Int_t ncols = fM.GetNcols(); + + v.ResizeTo(ncols); + + while (ncols--) + v(ncols) = fM(fRow, ncols); +} Index: trunk/MagicSoft/Mars/mhbase/MHMatrix.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MHMatrix.h (revision 7129) +++ trunk/MagicSoft/Mars/mhbase/MHMatrix.h (revision 7130) @@ -73,4 +73,6 @@ void AddColumns(MDataArray *mat); + //MDataArray *GetColumns() { return fData; } + const MDataArray *GetColumns() const { return fData; } MDataArray *GetColumns() { return fData; } @@ -102,4 +104,10 @@ Double_t operator[](Int_t col) { return fM(fRow, col); } + void GetRow(TVector &v) const; + void operator>>(TVector &v) const + { + GetRow(v); + } + Bool_t Fill(MParList *plist, MTask *read, MFilter *filter=NULL); Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7129) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7130) @@ -328,5 +328,5 @@ *fLog << " Simulated energy range: " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl; *fLog << " Simulated spectral slope: " << fOldSlope << endl; - *fLog << " New slope: " << fNewSlope << endl; + *fLog << " New spectral slope: " << fNewSlope << endl; *fLog << " User normalization: " << fNorm << endl; *fLog << " Old Spectrum: " << GetFormulaSpecOldX() << " (I=" << GetSpecOldIntegral() << ")" << endl; Index: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 7129) +++ trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 7130) @@ -1122,5 +1122,5 @@ { fillpul.SetFilter(&fcalib); - tlist.AddToList(&decode); + tlist.AddToList(&decode); tlist.AddToList(&fcalib); tlist.AddToList(&fillpul); @@ -1255,5 +1255,5 @@ if (!calc.CheckPedestalRms(fBadPixels, fPedestalCamOut)) { - *fLog << err << "ERROR - Checking PedestalRMS via MBadPixelsCalc failed...." << endl; + *fLog << err << "ERROR - MBadPixelsCalc::CheckPedestalRms failed...." << endl; return kFALSE; } Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 7129) +++ trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 7130) @@ -242,28 +242,29 @@ // is not the first anymore // - if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal) - { + switch (fRunHeader->GetRunType()) + { + case MRawRunHeader::kRTPedestal: + case MRawRunHeader::kRTMonteCarlo: fIsFirstPedRun = kFALSE; fIsNotPedRun = kFALSE; return kTRUE; - } - - if (fRunHeader->GetRunType()==MRawRunHeader::kRTCalibration) - { - TString proj(fRunHeader->GetProjectName()); - proj.ToLower(); - - // test if a continuous light run has been declared as calibration... - if (proj.Contains("cl")) - { - fIsFirstPedRun = kFALSE; - fIsNotPedRun = kFALSE; - return kTRUE; - } - } - - + + case MRawRunHeader::kRTCalibration: + { + TString proj(fRunHeader->GetProjectName()); + proj.ToLower(); + + // test if a continuous light run has been declared as calibration... + if (proj.Contains("cl")) + { + fIsFirstPedRun = kFALSE; + fIsNotPedRun = kFALSE; + return kTRUE; + } + } + } fIsNotPedRun = kTRUE; + // // If this is the first call to ReInit (before reading the first file) @@ -464,5 +465,5 @@ } -//------------------------------------------------------------- +//----------------------------------------------------------------------- // // Return if the pedestal bit was set from the calibration trigger box. @@ -473,8 +474,8 @@ Bool_t MPedCalcPedRun::IsPedBitSet() { + if (fRunHeader->GetRunNumber()GetRunNumber() < gkFirstRunWithFinalBits) return kFALSE; Index: trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc (revision 7129) +++ trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc (revision 7130) @@ -57,5 +57,4 @@ #include "MParameters.h" - ClassImp(MRFEnergyEst); @@ -67,10 +66,10 @@ // -------------------------------------------------------------------------- // -// -MRFEnergyEst::MRFEnergyEst(const char *name, const char *title):fNumTrees(-1) -{ - // - // set the name and title of this object - // +// Default constructor. Set the name and title of this object +// +MRFEnergyEst::MRFEnergyEst(const char *name, const char *title) + : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fData(0), fEnergyEst(0), + fTestMatrix(0) +{ fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); @@ -79,78 +78,51 @@ // -------------------------------------------------------------------------- // -// Delete the data chains -// -MRFEnergyEst::~MRFEnergyEst() -{ - // delete fData; -} - -// -------------------------------------------------------------------------- -Int_t MRFEnergyEst::Train() -{ - if(!fMatrixTrain) - { - *fLog << err << dbginf << "fMatrixTrain not found... aborting." << endl; - return kFALSE; - } - - if(!fMatrixTrain->GetColumns()) - { - *fLog << err << dbginf << "fMatrixTrain does not contain rules... aborting." << endl; - return kFALSE; - } - - const Int_t ncols = (fMatrixTrain->GetM()).GetNcols(); - const Int_t nrows = (fMatrixTrain->GetM()).GetNrows(); - - cout<<"ncols="<GetColumns(); - MData &data=(*rules)[i]; - - if(iGetM()); - const Int_t nbins = fEnergyGrid.GetSize()-1; - if(nbins<=0) - { - *fLog << err << dbginf << "Energy grid not vaild... aborting." << endl; - return kFALSE; - } + TFile fileRF(fFileName, "recreate"); + if (!fileRF.IsOpen()) + { + *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl; + return kFALSE; + } + + MDataArray usedrules; + for(Int_t i=0; ipow(10.,fEnergyGrid[ie]) && energy<=pow(10.,fEnergyGrid[ie+1])) - { - TMatrixRow(mat1, irow1) = TMatrixRow(*mptr,j); - irow1++; - }else{ - TMatrixRow(mat0, irow0) = TMatrixRow(*mptr,j); - irow0++; - } + const Double_t energy = m(j,ncols-1); + + if (energy>grid[ie] && energy<=grid[ie+1]) + TMatrixFRow(mat1, irow1++) = TMatrixFRow_const(m,j); + else + TMatrixFRow(mat0, irow0++) = TMatrixFRow_const(m,j); } - if(irow1*irow0==0)continue; - - *fLog << inf << dbginf << "Training RF for energy bin "<0; + + if (!valid) + *fLog << warn << "WARNING - Skipping"; + else + *fLog << inf << "Training RF for"; + + *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ")" << endl; + + if (!valid) + continue; + + gLog.SetNullOutput(kTRUE); // last column must be removed (true energy col.) - mat1.ResizeTo(irow1,ncols-1); - mat0.ResizeTo(irow0,ncols-1); + mat1.ResizeTo(irow1, ncols-1); + mat0.ResizeTo(irow0, ncols-1); // create MHMatrix as input for RF - MHMatrix matrix1(mat1,"MatrixHadrons"); - MHMatrix matrix0(mat0,"MatrixGammas"); - - MDataArray *rules1=matrix1.GetColumns(); - MDataArray *rules0=matrix0.GetColumns(); - // rules of new matrices be re-set - if(rules1)delete rules1; rules1=new MDataArray(); - if(rules0)delete rules0; rules0=new MDataArray(); - - for(Int_t i=0;iGetColumns(); - //MData &data=(*rules)[i]; - MData &data=used_rules[i]; - rules1->AddEntry(data.GetRule()); - rules0->AddEntry(data.GetRule()); - } + MHMatrix matrix1(mat1, "MatrixHadrons"); + MHMatrix matrix0(mat0, "MatrixGammas"); + + matrix1.AddColumns(&usedrules); + matrix0.AddColumns(&usedrules); // training of RF + MTaskList tlist; MParList plist; - MTaskList tlist; plist.AddToList(&tlist); plist.AddToList(&matrix0); @@ -212,21 +180,22 @@ MEvtLoop evtloop; + evtloop.SetDisplay(fDisplay); evtloop.SetParList(&plist); - //gLog.SetNullOutput(kTRUE); - - if (!evtloop.Eventloop()) return kFALSE; - - //gLog.SetNullOutput(kFALSE); + gLog.SetNullOutput(kFALSE); + + if (!evtloop.Eventloop()) + return kFALSE; // save whole forest MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest"); - forest->SetTitle(Form("%f",0.5*(fEnergyGrid[ie]+fEnergyGrid[ie+1]))); - forest->Write(Form("%d",numbins)); - numbins++; + const TString title = Form("%f", TMath::Sqrt(grid[ie]*grid[ie+1])); + //const TString title = Form("%f", (grid[ie]+grid[ie+1])/2); + forest->SetTitle(title); + forest->Write(title); } // save rules - used_rules.Write("rules"); + usedrules.Write("rules"); fileRF.Close(); @@ -234,25 +203,19 @@ return kTRUE; } - -Int_t MRFEnergyEst::Test() -{ - if(!fMatrixTest) - { - *fLog << err << dbginf << "fMatrixTest not found... aborting." << endl; - return kFALSE; - } - - const Int_t ncols = (fMatrixTest->GetM()).GetNcols(); - const Int_t nrows = (fMatrixTest->GetM()).GetNrows(); - - if(ncols<=0 || nrows <=0) - { - *fLog << err << dbginf << "No. of columns or no. of rows of fMatrixTrain equal 0 ... aborting." << endl; - return kFALSE; - } - - TMatrix *mptr=(TMatrix*)&(fMatrixTest->GetM()); - - if(!ReadForests()) +/* +Int_t MRFEnergyEst::Test(const MHMatrix &matrixtest) +{ + gLog.Separator("MRFEnergyEst - Test"); + + const Int_t ncols = matrixtest.GetM().GetNcols(); + const Int_t nrows = matrixtest.GetM().GetNrows(); + + if (ncols<=0 || nrows <=0) + { + *fLog << err << dbginf << "No. of columns or no. of rows of matrixtrain equal 0 ... aborting." << endl; + return kFALSE; + } + + if (!ReadForests()) { *fLog << err << dbginf << "Reading RFs failed... aborting." << endl; @@ -262,85 +225,80 @@ const Int_t nbins=fEForests.GetSize(); - Double_t e_low = 100; - Double_t e_up = 0; + Double_t elow = FLT_MAX; + Double_t eup = -FLT_MAX; for(Int_t j=0;jGetTitle()),e_low); - e_up = TMath::Max(atof((fEForests[j])->GetTitle()),e_up); - } - - TH1F hres("hres","",1000,-10,10); - TH2F henergy("henergy","",100,e_low,e_up,100,e_low,e_up); - + elow = TMath::Min(atof(fEForests[j]->GetTitle()), elow); + eup = TMath::Max(atof(fEForests[j]->GetTitle()), eup); + } + + TH1F hres("hres", "", 1000, -10, 10); + TH2F henergy("henergy", "",100, elow, eup,100, elow, eup); + + const TMatrix &m=matrixtest.GetM(); for(Int_t i=0;iCalcHadroness(v); - Double_t e = atof((fEForests[j])->GetTitle()); - /*if(h>=hmax) - { - hmax=h; - e_est=pow(10.,e); - }*/ - hsum+=h; - e_est+=h*e; + const TVector v = TMatrixFRow_const(m,i); + + const Double_t h = ((MRanForest*)fEForests[j])->CalcHadroness(v); + const Double_t e = atof(fEForests[j]->GetTitle()); + + hsum += h; + eest += h*e; } - e_est/=hsum; - e_est=pow(10.,e_est); - - if(e_true>80.) hres.Fill((e_est-e_true)/e_true); - henergy.Fill(log10(e_true),log10(e_est)); - } - - if(TestBit(kEnableGraphicalOutput)) - { - if(gStyle) gStyle->SetOptFit(1); - // show results - TCanvas *c=MH::MakeDefCanvas("c","",300,900); - - c->Divide(1,3); - c->cd(1); - henergy.SetTitle("Estimated vs true energy"); - henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])"); - henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])"); - henergy.DrawCopy(); - - c->cd(2); - - TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S"); - hptr->SetTitle("Estimated vs true energy - profile"); - hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])"); - hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])"); - hptr->DrawCopy(); - - c->cd(3); - hres.Fit("gaus"); - hres.SetTitle("Energy resolution for E_{true}>80Gev"); - hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})"); - hres.GetYaxis()->SetTitle("counts"); - hres.DrawCopy(); - - - c->GetPad(1)->SetGrid(); - c->GetPad(2)->SetGrid(); - c->GetPad(3)->SetGrid(); - - } - - return kTRUE; -} - + eest /= hsum; + eest = pow(10., eest); + + //if (etrue>80.) + // hres.Fill((eest-etrue)/etrue); + + henergy.Fill(log10(etrue),log10(eest)); + } + + if(gStyle) + gStyle->SetOptFit(1); + + // show results + TCanvas *c=MH::MakeDefCanvas("c","",300,900); + + c->Divide(1,3); + c->cd(1); + henergy.SetTitle("Estimated vs true energy"); + henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])"); + henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])"); + henergy.DrawCopy(); + + c->cd(2); + TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S"); + hptr->SetTitle("Estimated vs true energy - profile"); + hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])"); + hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])"); + hptr->DrawCopy(); + + c->cd(3); + hres.Fit("gaus"); + hres.SetTitle("Energy resolution for E_{true}>80Gev"); + hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})"); + hres.GetYaxis()->SetTitle("counts"); + hres.DrawCopy(); + + c->GetPad(1)->SetGrid(); + c->GetPad(2)->SetGrid(); + c->GetPad(3)->SetGrid(); + + return kTRUE; +} +*/ Int_t MRFEnergyEst::ReadForests(MParList *plist) { - TFile fileRF(fRFfileName,"read"); - if(!fileRF.IsOpen()) + TFile fileRF(fFileName,"read"); + if (!fileRF.IsOpen()) { *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl; @@ -348,28 +306,28 @@ } - TList *list=(TList*)fileRF.GetListOfKeys(); - const Int_t n=list->GetSize()-1;// subtract 1 since 1 key belongs to MDataArray - - fEForests.Expand(n); - - MRanForest forest; - for(Int_t i=0;iFindObject(Form("%d",i))->GetTitle(); - const char *bin =list->FindObject(Form("%d",i))->GetName(); - - curforest->SetTitle(energy); - curforest->SetName(bin); - - fEForests[i]=curforest; - } - - if(plist) - { - fData=(MDataArray*)plist->FindCreateObj("MDataArray"); + fEForests.Delete(); + + TIter Next(fileRF.GetListOfKeys()); + TObject *o=0; + while ((o=Next())) + { + MRanForest *forest; + fileRF.GetObject(o->GetName(), forest); + if (!forest) + continue; + + forest->SetTitle(o->GetTitle()); + forest->SetBit(kCanDelete); + + fEForests.Add(forest); + + } + + if (plist) + { + fData = (MDataArray*)plist->FindCreateObj("MDataArray"); fData->Read("rules"); } + fileRF.Close(); @@ -381,18 +339,18 @@ fEnergyEst = (MParameterD*)plist->FindCreateObj("MParameterD", "MEnergyEst"); if (!fEnergyEst) - { - *fLog << err << dbginf << "MEnergyEst [MParameterD] not found... aborting." << endl; - return kFALSE; - } - - if(!ReadForests(plist)) - { - *fLog << err << dbginf << "Reading RFs failed... aborting." << endl; - return kFALSE; - } + return kFALSE; + + if (!ReadForests(plist)) + { + *fLog << err << "Reading RFs failed... aborting." << endl; + return kFALSE; + } + + if (fTestMatrix) + return kTRUE; if (!fData) { - *fLog << err << dbginf << "MDataArray not found... aborting." << endl; + *fLog << err << "MDataArray not found... aborting." << endl; return kFALSE; } @@ -400,5 +358,5 @@ if (!fData->PreProcess(plist)) { - *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl; + *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl; return kFALSE; } @@ -413,26 +371,24 @@ { TVector event; - *fData >> event; + if (fTestMatrix) + *fTestMatrix >> event; + else + *fData >> event; Double_t eest = 0; - //Double_t hmax = 0; - Double_t hsum = 0; - - for(Int_t j=0;jCalcHadroness(event); - Double_t e = atof((fEForests[j])->GetTitle()); - /*if(h>=hmax) - { - hmax=h; - e_est=pow(10.,e); - }*/ - hsum+=h; - eest+=h*e; - } - eest/=hsum; - eest=pow(10.,eest); - - fEnergyEst->SetVal(eest); + Double_t hsum = 0; + + TIter Next(&fEForests); + MRanForest *rf = 0; + while ((rf=(MRanForest*)Next())) + { + const Double_t h = rf->CalcHadroness(event); + const Double_t e = atof(rf->GetTitle()); + + hsum += h; + eest += h*e; + } + + fEnergyEst->SetVal(eest/hsum); fEnergyEst->SetReadyToSave(); Index: trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.h (revision 7129) +++ trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.h (revision 7130) @@ -6,11 +6,9 @@ #endif -#ifndef ROOT_TArrayD -#include -#endif - #ifndef ROOT_TObjArray #include #endif + +class TArrayD; class MHMatrix; @@ -21,17 +19,15 @@ { private: - Int_t fNumTrees; - Int_t fNumTry; - Int_t fNdSize; + Int_t fNumTrees; // Training parameters + Int_t fNumTry; // Training parameters + Int_t fNdSize; // Training parameters - TString fRFfileName; - MHMatrix *fMatrixTrain; - MHMatrix *fMatrixTest; - TArrayD fEnergyGrid; + TString fFileName; + TObjArray fEForests; MDataArray *fData; //! Used to store the MDataChains to get the event values - TObjArray fEForests; + MParameterD *fEnergyEst; //! Used to storeestimated energy - MParameterD *fEnergyEst; + MHMatrix *fTestMatrix; Int_t PreProcess(MParList *plist); @@ -42,18 +38,15 @@ public: MRFEnergyEst(const char *name=NULL, const char *title=NULL); - ~MRFEnergyEst(); - void SetMatrixTrain(MHMatrix *mat) { fMatrixTrain = mat; } - void SetMatrixTest( MHMatrix *mat) { fMatrixTest = mat; } - void SetFile(TString str) { fRFfileName = str; } + void SetFileName(TString str) { fFileName = str; } - void SetLogEnergyGrid(TArrayD &egrid) { fEnergyGrid = egrid ; } + void SetNumTrees(UShort_t n=-1) { fNumTrees = n; } + void SetNdSize(UShort_t n=-1) { fNdSize = n; } + void SetNumTry(UShort_t n=-1) { fNumTry = n; } - void SetNumTrees(UShort_t n=100) { fNumTrees = n; } - void SetNdSize(UShort_t n=1) { fNdSize = n; } - void SetNumTry(UShort_t n) { fNumTry = n; } + void SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; } + void InitMapping(MHMatrix *m=0) { fTestMatrix=m; } - Int_t Train(); - Int_t Test(); + Int_t Train(const MHMatrix &n, const TArrayD &grid); ClassDef(MRFEnergyEst, 0) // Task Index: trunk/MagicSoft/Mars/mranforest/MRanForestGrow.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestGrow.cc (revision 7129) +++ trunk/MagicSoft/Mars/mranforest/MRanForestGrow.cc (revision 7130) @@ -46,5 +46,5 @@ using namespace std; -static const TString gsDefName = "MRanForestGrow"; +static const TString gsDefName = "MRead"; static const TString gsDefTitle = "Tree Classification Loop 1/2"; @@ -55,5 +55,4 @@ // MRanForestGrow::MRanForestGrow(const char *name, const char *title) - : fNumTrees(100),fNumTry(3),fNdSize(1) { // @@ -62,4 +61,8 @@ fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); + + SetNumTrees(); + SetNumTry(); + SetNdSize(); } @@ -124,5 +127,6 @@ Int_t MRanForestGrow::Process() { - Bool_t not_last=fRanForest->GrowForest(); + const Bool_t not_last=fRanForest->GrowForest(); + fRanTree->SetReadyToSave(); Index: trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h (revision 7129) +++ trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h (revision 7130) @@ -2,14 +2,6 @@ #define MARS_MRanForestGrow -///////////////////////////////////////////////////////////////////////////// -// // -// MRanForestGrow // -// // -// Task to grow (train) a random forest // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MTask -#include "MTask.h" +#ifndef MARS_MRead +#include "MRead.h" #endif @@ -19,5 +11,5 @@ class MRanTree; -class MRanForestGrow : public MTask +class MRanForestGrow : public MRead { private: @@ -35,10 +27,14 @@ Int_t PostProcess(); + UInt_t GetEntries() { return fNumTrees; } + TString GetFileName() const { return "MRanForestGrow"; } + TString GetFullFileName() const { return "MRanForestGrow"; } + public: MRanForestGrow(const char *name=NULL, const char *title=NULL); - void SetNumTrees(Int_t n){ fNumTrees=n;} - void SetNumTry(Int_t n) { fNumTry=n; } - void SetNdSize(Int_t n) { fNdSize=n; } + void SetNumTrees(Int_t n=-1) { fNumTrees=n>0?n:100; } + void SetNumTry(Int_t n=-1) { fNumTry =n>0?n: 3; } + void SetNdSize(Int_t n=-1) { fNdSize =n>0?n: 1; } ClassDef(MRanForestGrow, 0) // Task to grow a random forest Index: trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc (revision 7129) +++ trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc (revision 7130) @@ -109,5 +109,5 @@ Bool_t MTFillMatrix::CheckResult(MHMatrix *m, Int_t num) const { - if (!m) + if (!m || num==0) return kTRUE; @@ -190,5 +190,6 @@ fLog->Separator(GetDescriptor()); *fLog << "Fill " << fDestMatrix1->GetDescriptor() << " with " << fNumDestEvents1 << endl; - *fLog << "Fill " << fDestMatrix2->GetDescriptor() << " with " << fNumDestEvents2 << endl; + if (fDestMatrix2) + *fLog << "Fill " << fDestMatrix2->GetDescriptor() << " with " << fNumDestEvents2 << endl; *fLog << "Distribution choosen "; if (fReference && fReference->GetHist().GetEntries()>0) @@ -210,23 +211,26 @@ // FIXME: Merge MFEventSelector and MFEventSelector2 MFilter *selector=0; - if (fReference) - { - // Case of a reference/nominal distribution - // The events must be read before selection - MFEventSelector2 *sel = new MFEventSelector2(*fReference); - sel->SetNumMax(fNumDestEvents1+fNumDestEvents2); - sel->SetInverted(); - - selector = sel; - } - else - { - // Case of a random distribution - // The events can be selected before reading - MFEventSelector *sel = new MFEventSelector; - sel->SetNumSelectEvts(fNumDestEvents1+fNumDestEvents2); - fReader->SetSelector(sel); - - selector = sel; + if (fNumDestEvents1>0 || fNumDestEvents2>0) + { + if (fReference) + { + // Case of a reference/nominal distribution + // The events must be read before selection + MFEventSelector2 *sel = new MFEventSelector2(*fReference); + sel->SetNumMax(fNumDestEvents1+fNumDestEvents2); + sel->SetInverted(); + + selector = sel; + } + else + { + // Case of a random distribution + // The events can be selected before reading + MFEventSelector *sel = new MFEventSelector; + sel->SetNumSelectEvts(fNumDestEvents1+fNumDestEvents2); + fReader->SetSelector(sel); + + selector = sel; + } } @@ -256,15 +260,18 @@ MFillH fill1(fDestMatrix1); MFillH fill2(fDestMatrix2); - fill1.SetFilter(&split); - fill2.SetFilter(&invsplit); + if (selector) + { + fill1.SetFilter(&split); + fill2.SetFilter(&invsplit); + } // entries in MTaskList tlist.AddToList(fReader); // Read events - if (fReference) + if (fReference && selector) tlist.AddToList(&cont); // select a sample of events tlist.AddToList(&invsplit); // process invsplit (which implicitly processes split) - if (fDestMatrix1 && fNumDestEvents1>0) + if (fDestMatrix1) tlist.AddToList(&fill1); // fill matrix 1 - if (fDestMatrix2 && fNumDestEvents2>0) + if (fDestMatrix2) tlist.AddToList(&fill2); // fill matrix 2 if (fWriteFile1) @@ -289,9 +296,6 @@ const Bool_t rc = evtloop.Eventloop(); - // Print execution statistics of the tasklist - if (rc) - tlist.PrintStatistics(); - - delete selector; + if (selector) + delete selector; if (!rc)