Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2356) +++ trunk/MagicSoft/Mars/Changelog (revision 2357) @@ -1,3 +1,51 @@ -*-*- END OF LINE -*-*- + + + 2003/09/24: Wolfgang Wittek + + * mfilter/MFEventSelector2.[h,cc] + - execution statistics added + + * mhist/MHFindSignificance.cc + - add fHist->UseCurrentStyle() + to get the y-axis + labels drawn + + * mhist/MHMatrix.h + - replace Int_t fNumRow //! + by Int_t fNumRow // + because otherwise fNumRow is not defined when MHMatrix object is read in + after it had been written out + + * mhist/MHCT1Supercuts.cc + - change title of object + + * manalysis/MMinuitInterface.cc + - add arguments maxcalls and tolerance to SIMPLEX call + + * manalysis/MCT1SupercutsCalc.[h,cc] + - add variables asymmetry, conc, leakage + + * manalysis/MCT1Supercuts.[h,cc] + - add variables asymmetry, conc, leakage + - add TArrayD fStepsizes (initial step sizes for the parameters) + + * manalysis/MCT1FindSupercuts.cc + - replace MGeomCamCT1Daniel by MGeomCamCT1 + - arguments 'parSCinit', 'params' and 'steps' added in FindParams() ; + parSCinit is the name of the file containing the initial + values of the parameters + 'params' and 'steps' are the initial values in case parSCinit == "" + - add member functions GetMatrixTrain() and GetMatrixTest() + - remove member function WriteMatrix() + because it didn't work; now the matrices are written out in + DefineTrainMatrix(), DefineTestMatri() and DefineTrainTestMatrix() + + * macros/CT1EgyEst.C + - don't use Daniel's energy estimator + + * mmontecarlo/MMcEnergyEst.cc + - extend printout of comments + + 2003/09/17: Abelardo Moralejo @@ -51,5 +99,4 @@ - added 'Time Spectra of Cosmics' button - added size argument to instatiation of MHFadcCam - Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2356) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2357) @@ -22,4 +22,5 @@ #include "MWriteRootFile.h" +//#include "TH3D.h" @@ -77,12 +78,12 @@ MBinning *binth = new MBinning("BinningTheta"); - //Double_t yedge[9] = - // {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; // this is Daniel's binning in theta - Double_t yedge[8] = - {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99}; - //TArrayD yed(9,yedge); + //Double_t yedge[8] = + // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99}; + // this is our binning + Double_t yedge[9] = + {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; TArrayD yed; - yed.Set(8,yedge); + yed.Set(9,yedge); binth->SetEdges(yed); plist->AddToList(binth); @@ -94,5 +95,6 @@ zedge[8-i] = cos(yedge[i]/kRad2Deg); } - TArrayD zed(9,zedge); + TArrayD zed; + zed.Set(9,zedge); bincosth->SetEdges(zed); plist->AddToList(bincosth); @@ -116,4 +118,5 @@ } + void DeleteBinnings(MParList *plist) { @@ -122,4 +125,10 @@ TObject *bin; + //-------------------------------------------- + bin = plist->FindObject("BinningE"); + if (bin) delete bin; + + //-------------------------------------------- + bin = plist->FindObject("BinningSize"); if (bin) delete bin; @@ -159,5 +168,6 @@ if (bin) delete bin; - //robert + + // robert ---------------------------------------------- bin = plist->FindObject("BinningAlphaFlux"); if (bin) delete bin; @@ -170,7 +180,11 @@ } + //************************************************************************ void CT1Analysis() { + + gLog << "Entry CT1Analysis()" << endl; + gLog.SetNoColors(); @@ -189,8 +203,8 @@ // path for input for Mars - TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/"; + TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/"; // path for output from Mars - TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/"; + TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/"; //----------------------------------------------- @@ -217,27 +231,4 @@ Bool_t JobA_MC = kFALSE; Bool_t WMC1 = kFALSE; // write out root file MC1.root ? - - - // Job B_NN_UP : read ON1.root (or MC1.root) file - // - depending on RlookNN : create (or read in) hadron and gamma matrix - // - calculate hadroness for method of NEAREST NEIGHBORS - // and for the SUPERCUTS - // - update the input files with the hadronesses (ON1.root or MC1.root) - - Bool_t JobB_NN_UP = kFALSE; - Bool_t RLookNN = kFALSE; // read in look-alike events - Bool_t WNN = kFALSE; // update input root file ? - - - // Job B_SC_UP : read ON1.root (or MC1.root) file - // - depending on WParSC : create (or read in) supercuts parameter values - // - calculate hadroness for the SUPERCUTS - // - update the input files with the hadroness (ON1.root or MC1.root) - - Bool_t JobB_SC_UP = kTRUE; - Bool_t RMatrix = kFALSE; // read matrices from file - Bool_t WParSC = kTRUE; // do optimization and write supercuts - // parameter values onto a file - Bool_t WSC = kFALSE; // update input root file ? @@ -252,4 +243,18 @@ Bool_t RTree = kFALSE; // read in trees Bool_t WRF = kFALSE; // update input root file ? + + + + + // Job B_SC_UP : read ON1.root (or MC1.root) file + // - depending on WParSC : create (or read in) supercuts parameter values + // - calculate hadroness for the SUPERCUTS + // - update the input files with the hadroness (ON1.root or MC1.root) + + Bool_t JobB_SC_UP = kFALSE; + Bool_t RMatrix = kFALSE; // read training and test matrices from file + Bool_t WParSC = kFALSE; // do optimization and write supercuts + // parameter values onto a file + Bool_t WSC = kFALSE; // update input root file ? @@ -299,6 +304,9 @@ // - write root file for ON data after final cuts - Bool_t JobF_XX = kFALSE; - Bool_t WFX = kFALSE; // write out root file ? + gLog << "before setting switches for JobF_XX" << endl; + + Bool_t JobF_XX = kTRUE; + Bool_t WFX = kTRUE; // write out root file ? + Bool_t WRobert = kTRUE; // write out Robert's file ? @@ -327,5 +335,7 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobA_ON, WHistON, WON1 = " - << JobA_ON << ", " << WHistON << ", " << WON1 << endl; + << (JobA_ON ? "kTRUE" : "kFALSE") << ", " + << (WHistON ? "kTRUE" : "kFALSE") << ", " + << (WON1 ? "kTRUE" : "kFALSE") << endl; @@ -424,7 +434,6 @@ contstandard.SetName("SelStandard"); - if (WON1) - { - MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + //MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + MWriteRootFile write(outNameImage, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -438,5 +447,5 @@ write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); - } + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ @@ -500,5 +509,4 @@ Int_t maxevents = -1; - //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; @@ -587,5 +595,6 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobA_MC, WMC1 = " - << JobA_MC << ", " << WMC1 << endl; + << (JobA_MC ? "kTRUE" : "kFALSE") << ", " + << (WMC1 ? "kTRUE" : "kFALSE") << endl; @@ -688,6 +697,11 @@ MCT1ReadPreProc read(filenamein); + MBlindPixelCalc blindbeforepad; + blindbeforepad.SetUseBlindPixels(); + blindbeforepad.SetName("BlindBeforePadding"); + MBlindPixelCalc blind; blind.SetUseBlindPixels(); + blind.SetName("BlindAfterPadding"); MFCT1SelBasic selbasic; @@ -758,7 +772,6 @@ - if (WMC1) - { - MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + //MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + MWriteRootFile write(outNameImage, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -772,5 +785,5 @@ write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); - } + @@ -788,4 +801,5 @@ tlist.AddToList(&read); + tlist.AddToList(&blindbeforepad); tlist.AddToList(&padthomas); tlist.AddToList(&blind); @@ -855,28 +869,35 @@ } - - //--------------------------------------------------------------------- - // Job B_NN_UP + // Job B_RF_UP //============ - // - read in the matrices of look-alike events for gammas and hadrons - - // then read ON1.root (or MC1.root) file - // - // - calculate the hadroness for the method of nearest neighbors - // - // - update input root file, including the hadroness - - - if (JobB_NN_UP) + + // - create (or read in) the matrices of look-alike events for gammas + // and hadrons + // - create (or read in) the trees + // - then read ON1.root (or MC1.root) file + // - calculate the hadroness for the method of RANDOM FOREST + // - update input root file with the hadroness + + + if (JobB_RF_UP) { gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; + gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " - << JobB_NN_UP << ", " << RLookNN << ", " << WNN << endl; - + gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " + << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " + << (RLookRF ? "kTRUE" : "kFALSE") << ", " + << (RTree ? "kTRUE" : "kFALSE") << ", " + << (WRF ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // parameters for the random forest + Int_t NumTrees = 100; + Int_t NumTry = 3; + Int_t NdSize = 1; @@ -898,5 +919,4 @@ outNameImage += typeInput; outNameImage += "2.root"; - //TString outNameImage = filenameData; @@ -909,5 +929,5 @@ filenameON += "ON"; filenameON += "1.root"; - Int_t howManyHadrons = 500000; + Int_t howManyHadrons = 35000; gLog << "filenameON = " << filenameON << ", howManyHadrons = " << howManyHadrons << endl; @@ -918,5 +938,5 @@ filenameMC += "MC"; filenameMC += "1.root"; - Int_t howManyGammas = 50000; + Int_t howManyGammas = 35000; gLog << "filenameMC = " << filenameMC << ", howManyGammas = " << howManyGammas << endl; @@ -926,5 +946,5 @@ TString outNameGammas = outPath; - outNameGammas += "matrix_gammas_"; + outNameGammas += "RFmatrix_gammas_"; outNameGammas += "MC"; outNameGammas += ".root"; @@ -934,5 +954,5 @@ TString outNameHadrons = outPath; - outNameHadrons += "matrix_hadrons_"; + outNameHadrons += "RFmatrix_hadrons_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; @@ -940,9 +960,19 @@ MHMatrix matrixg("MatrixGammas"); + matrixg.EnableGraphicalOutput(); + MHMatrix matrixh("MatrixHadrons"); + matrixh.EnableGraphicalOutput(); + + //-------------------------------------------- + // file of trees of the random forest + + TString outRF = outPath; + outRF += "RF.root"; + //************************************************************************* // read in matrices of look-alike events -if (RLookNN) +if (RLookRF) { const char* mtxName = "MatrixGammas"; @@ -961,5 +991,5 @@ matrixg.Read(mtxName); - matrixg.Print("cols"); + matrixg.Print("SizeCols"); @@ -981,5 +1011,5 @@ matrixh.Read(mtxName); - matrixh.Print("cols"); + matrixh.Print("SizeCols"); } @@ -987,5 +1017,5 @@ //************************************************************************* // create matrices of look-alike events -if (!RLookNN) +if (!RLookRF) { gLog << "" << endl; @@ -1011,13 +1041,23 @@ MFParticleId fgamma("MMcEvt", '=', kGAMMA); fgamma.SetName("gammaID"); + MFParticleId fhadrons("MMcEvt", '!', kGAMMA); fhadrons.SetName("hadronID)"); - MFEventSelector selectorg; - selectorg.SetNumSelectEvts(howManyGammas); + TString mgname("costhg"); + MBinning bing("Binning"+mgname); + bing.SetEdges(10, 0., 1.0); + + MH3 gref("cos(MMcEvt.fTelescopeTheta)"); + gref.SetName(mgname); + MH::SetBinning(&gref.GetHist(), &bing); + for (Int_t i=1; i<=gref.GetNbins(); i++) + gref.GetHist().SetBinContent(i, 1.0); + + + MFEventSelector2 selectorg(gref); + selectorg.SetNumMax(howManyGammas); selectorg.SetName("selectGammas"); - MFEventSelector selectorh; - selectorh.SetNumSelectEvts(howManyHadrons); - selectorh.SetName("selectHadrons"); + MFillH fillmatg("MatrixGammas"); @@ -1029,5 +1069,5 @@ fillmath.SetName("fillHadrons"); - + //***************************** fill gammas *** // entries in MFilterList @@ -1065,8 +1105,26 @@ tlistg.PrintStatistics(0, kTRUE); - + //***************************** fill hadrons *** gLog << " Hadrons :" << endl; + + TString mhname("costhh"); + MBinning binh("Binning"+mhname); + binh.SetEdges(10, 0., 1.0); + + MH3 href("cos(MMcEvt.fTelescopeTheta)"); + href.SetName(mhname); + MH::SetBinning(&href.GetHist(), &binh); + for (Int_t i=1; i<=href.GetNbins(); i++) + href.GetHist().SetBinContent(i, 1.0); + + MFEventSelector2 selectorh(href); + //selectorh.SetNumMax(howManyHadrons); + // select as many hadrons as gammas + selectorh.SetNumMax(matrixg.GetM().GetNrows()); + selectorh.SetName("selectHadrons"); + + // entries in MFilterList @@ -1102,80 +1160,9 @@ tlisth.PrintStatistics(0, kTRUE); + + + //***************************************************** - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - // - // ---------------------------------------------------------- - // Definition of the reference sample of look-alike events - // (this is a subsample of the original sample) - // ---------------------------------------------------------- - // - gLog << "" << endl; - gLog << "========================================================" << endl; - // Select a maximum of nmaxevts events from the sample of look-alike - // events. They will form the reference sample. - Int_t nmaxevts = 2000; - - // target distribution for the variable in column refcolumn (start at 0); - //Int_t refcolumn = 0; - //Int_t nbins = 5; - //Float_t frombin = 0.5; - //Float_t tobin = 1.0; - //TH1F *thsh = new TH1F("thsh","target distribution", - // nbins, frombin, tobin); - //thsh->SetDirectory(NULL); - //thsh->SetXTitle("cos( \\Theta)"); - //thsh->SetYTitle("Counts"); - //Float_t dbin = (tobin-frombin)/nbins; - //Float_t lbin = frombin +dbin*0.5; - //for (Int_t j=1; j<=nbins; j++) - //{ - // thsh->Fill(lbin,1.0); - // lbin += dbin; - //} - - Int_t refcolumn = 0; - MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); - TH1F *thsh = new TH1F(); - thsh->SetNameTitle("thsh","target distribution"); - MH::SetBinning(thsh, binscostheta); - Int_t nbins = thsh->GetNbinsX(); - Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; - for (Int_t j=1; j<=nbins; j++) - { - thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixg.EnableGraphicalOutput(); - Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; - return; - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixh.EnableGraphicalOutput(); - matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); - delete thsh; - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; - return; - } - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // write out look-alike events @@ -1189,5 +1176,5 @@ // "gammas" gLog << "Gammas :" << endl; - matrixg.Print("cols"); + matrixg.Print("SizeCols"); TFile writeg(outNameGammas, "RECREATE", ""); @@ -1201,5 +1188,5 @@ // "hadrons" gLog << "Hadrons :" << endl; - matrixh.Print("cols"); + matrixh.Print("SizeCols"); TFile writeh(outNameHadrons, "RECREATE", ""); @@ -1214,14 +1201,126 @@ + MRanForest *fRanForest; + MRanTree *fRanTree; //----------------------------------------------------------------- - // Update the input files with the NN and SC hadronness - // - - if (WNN) + // read in the trees of the random forest + if (RTree) + { + MParList plisttr; + MTaskList tlisttr; + plisttr.AddToList(&tlisttr); + + MReadTree readtr("TREE", outRF); + readtr.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(NumTrees); + + // list of tasks for the loop over the trees + + tlisttr.AddToList(&readtr); + tlisttr.AddToList(&rffill); + + //------------------- + // Execute tree loop + // + MEvtLoop evtlooptr; + evtlooptr.SetParList(&plisttr); + if (!evtlooptr.Eventloop()) + return; + + tlisttr.PrintStatistics(0, kTRUE); + + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + + //----------------------------------------------------------------- + // grow the trees of the random forest (event loop = tree loop) + + if (!RTree) + { + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : start growing trees" << endl; + + MTaskList tlist2; + MParList plist2; + plist2.AddToList(&tlist2); + + plist2.AddToList(&matrixg); + plist2.AddToList(&matrixh); + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(NumTrees); + rfgrow2.SetNumTry(NumTry); + rfgrow2.SetNdSize(NdSize); + + MWriteRootFile rfwrite2(outRF); + rfwrite2.AddContainer("MRanTree", "TREE"); + + // list of tasks for the loop over the trees + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + + //------------------- + // Execute tree loop + // + MEvtLoop treeloop; + treeloop.SetParList(&plist2); + + if ( !treeloop.Eventloop() ) + return; + + tlist2.PrintStatistics(0, kTRUE); + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + // end of growing the trees of the random forest + //----------------------------------------------------------------- + + + + + //----------------------------------------------------------------- + // Update the input files with the RF hadronness + // + + if (WRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update input file '" << filenameData - << "' with the NN and SC hadronness" << endl; + << "' with the RF hadronness" << endl; MTaskList tliston; @@ -1246,19 +1345,9 @@ //....................................................................... - // calculation of hadroness for method of Nearest Neighbors - - TString hadNNName = "HadNN"; - MMultiDimDistCalc nncalc; - nncalc.SetUseNumRows(25); - nncalc.SetUseKernelMethod(kFALSE); - nncalc.SetHadronnessName(hadNNName); - - //....................................................................... - // calculation of hadroness for the supercuts - // (=0.25 if fullfilled, =0.75 otherwise) - - TString hadSCName = "HadSC"; - MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); - sccalc.SetHadronnessName(hadSCName); + // calculate hadronnes for method of RANDOM FOREST + + TString hadRFName = "HadRF"; + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadRFName); //....................................................................... @@ -1278,6 +1367,5 @@ write.AddContainer("MNewImagePar", "Events"); - write.AddContainer("HadNN", "Events"); - write.AddContainer("HadSC", "Events"); + write.AddContainer(hadRFName, "Events"); @@ -1287,4 +1375,5 @@ (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + Float_t maxhadronness = 0.40; Float_t maxalpha = 20.0; @@ -1293,21 +1382,30 @@ MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadNNName); + selfinalgh.SetHadronnessName(hadRFName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); - MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); - fillhadnn.SetName("HhadNN"); - MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); - fillhadsc.SetName("HhadSC"); + MFillH fillranfor("MHRanForest"); + fillranfor.SetName("HRanForest"); + + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadNNName); + selfinal.SetHadronnessName(hadRFName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); MFillH hfill1("MHHillas", fHilName); @@ -1332,7 +1430,9 @@ InitBinnings(&pliston); - pliston.AddToList(&matrixg); - pliston.AddToList(&matrixh); - + pliston.AddToList(fRanForest); + pliston.AddToList(fRanTree); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); //***************************** @@ -1341,9 +1441,12 @@ tliston.AddToList(&read); - tliston.AddToList(&nncalc); - tliston.AddToList(&sccalc); - tliston.AddToList(&fillhadnn); - tliston.AddToList(&fillhadsc); - + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillranfor); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); @@ -1352,7 +1455,4 @@ tliston.AddToList(&hfill5); - tliston.AddToList(&write); - - tliston.AddToList(&contfinalgh); tliston.AddToList(&contfinal); @@ -1373,9 +1473,10 @@ tliston.PrintStatistics(0, kTRUE); + //------------------------------------------- // Display the histograms // - pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); - pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + pliston.FindObject("MHRanForest")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); @@ -1385,13 +1486,49 @@ pliston.FindObject("MHStarMap")->DrawClone(); + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); + TH1 &alphaHist = absalpha->GetHist(); + alphaHist.SetXTitle("|alpha| [\\circ]"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 4; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << filenameData << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + DeleteBinnings(&pliston); } - - - gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl; + gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; gLog << "=======================================================" << endl; } //--------------------------------------------------------------------- + + @@ -1413,13 +1550,36 @@ gLog << "" << endl; - //gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = " - // << JobB_SC_UP << ", " << WParSC << ", " << WSC << endl; - + gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = " + << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " + << (WParSC ? "kTRUE" : "kFALSE") << ", " + << (WSC ? "kTRUE" : "kFALSE") << endl; + + + + //-------------------------------------------- + // file which contains the initial parameter values for the supercuts + // if parSCinit ="" the initial values are taken from the constructor of + // MCT1Supercuts + + TString parSCinit = outPath; + //parSCinit += "parSC_1709d"; + parSCinit = ""; + + gLog << "parSCinit = " << parSCinit << endl; + + //--------------- + // file onto which the optimal parameter values for the supercuts + // are written + + TString parSCfile = outPath; + parSCfile += "parSC_2209a"; + + gLog << "parSCfile = " << parSCfile << endl; //-------------------------------------------- // file to be updated (either ON or MC) - TString typeInput = "ON"; - //TString typeInput = "MC"; + //TString typeInput = "ON"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -1434,5 +1594,4 @@ outNameImage += typeInput; outNameImage += "3.root"; - outNameImage += "xxxxx"; @@ -1447,5 +1606,5 @@ TString filenameTrain = outPath; filenameTrain += "ON"; - filenameTrain += "2.root"; + filenameTrain += "1.root"; Int_t howManyTrain = 800000; gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " @@ -1455,6 +1614,6 @@ TString filenameTest = outPath; filenameTest += "ON"; - filenameTest += "2.root"; - Int_t howManyTest = 100000; + filenameTest += "1.root"; + Int_t howManyTest = 800000; gLog << "filenameTest = " << filenameTest << ", howManyTest = " @@ -1479,19 +1638,11 @@ - //-------------------------------------------- - // file which contains the optimum parameter values for the supercuts - - TString parSCfile = outPath; - parSCfile += "parSC"; - - gLog << "parSCfile = " << parSCfile << endl; - //--------------------------------------------------------------------- - // optimize supercuts using the training sample - // and test them using the test sample - -if (WParSC) - { + // Training and test matrices : + // - either create them and write them onto a file + // - or read them from a file + + MCT1FindSupercuts findsuper; findsuper.SetFilenameParam(parSCfile); @@ -1500,5 +1651,5 @@ //-------------------------- - // create matrices + // create matrices and write them onto files if (!RMatrix) { @@ -1509,6 +1660,6 @@ { if ( !findsuper.DefineTrainTestMatrix(filenameTrain, - howManyTrain, mh3, - howManyTest, mh3) ) + howManyTrain, mh3, howManyTest, mh3, + fileMatrixTrain, fileMatrixTest) ) { *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl; @@ -1519,5 +1670,6 @@ else { - if ( !findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3) ) + if ( !findsuper.DefineTrainMatrix(filenameTrain, + howManyTrain, mh3, fileMatrixTrain) ) { *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl; @@ -1525,5 +1677,6 @@ } - if ( !findsuper.DefineTestMatrix( filenameTest, howManyTest, mh3) ) + if ( !findsuper.DefineTestMatrix( filenameTest, + howManyTest, mh3, fileMatrixTest) ) { *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl; @@ -1531,24 +1684,124 @@ } } - - // writing doesn't work ??? - //findsuper.WriteMatrix(fileMatrixTrain, fileMatrixTest); - } + } //-------------------------- // read matrices from files - // NOT YET TESTED - //if (RMatrix) - // findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); + // + + if (RMatrix) + findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); //-------------------------- - //-------------------------- - // optimize the supercuts - Bool_t rf = findsuper.FindParams(); + + //--------------------------------------------------------------------- + // optimize supercuts using the training sample + // + // the initial values are taken + // - from the file parSCinit (if != "") + // - or from the arrays params and steps (if their sizes are != 0) + // - or from the MCT1Supercuts constructor + + +if (WParSC) + { + TArrayD params(0); + TArrayD steps(0); + + if (parSCinit == "") + { + Double_t vparams[104] = { + // LengthUp + 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362, + 0.007388, -0.013463, + // LengthLo + 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993, + 0.000080, -0.007157, + // WidthUp + 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353, + 0.020711, -0.016703, + // WidthLo + 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101, + 0.006126, -0.002849, + // DistUp + 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909, + 0.189697, 0.0, + // DistLo + 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374, + -0.001750, 0.0, + // AsymUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + Double_t vsteps[104] = { + // LengthUp + 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002, + 0.0008, 0.002, + // LengthLo + 0.02, 0.003, 0.05, 0.006, 0.002, 0.3, + 0.0001, 0.0008, + // WidthUp + 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006, + 0.002, 0.002, + // WidthLo + 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002, + 0.0007, 0.003, + // DistUp + 0.2, 0.0, 0.3, 0.02, 0.0, 0.03, + 0.02, 0.0 + // DistLo + 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005, + 0.0002, 0.0 + // AsymUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + params.Set(104, vparams); + steps.Set (104, vsteps ); + } + + Bool_t rf; + rf = findsuper.FindParams(parSCinit, params, steps); if (!rf) { - *fLog << "CT1Analysis.C : optimization of supercuts failed" << endl; + gLog << "CT1Analysis.C : optimization of supercuts failed" << endl; return; } @@ -1565,5 +1818,6 @@ // - + if (WSC) + { gLog << "" << endl; gLog << "========================================================" << endl; @@ -1580,5 +1834,5 @@ inparam.Close(); - gLog << "Optimum parameter values for supercuts were read in from file '" + gLog << "Parameter values for supercuts were read in from file '" << parSCfile << "'" << endl; @@ -1586,8 +1840,18 @@ supercutsPar = scin.GetParameters(); - gLog << "Optimum parameter values for supercuts : " << endl; - for (Int_t i=0; i<72; i++) + TArrayD supercutsStep; + supercutsStep = scin.GetStepsizes(); + + gLog << "Parameter values for supercuts : " << endl; + for (Int_t i=0; iGetHist(); + MFillH alpha(&alphaabs); alpha.SetName("FillAlphaAbs"); @@ -1716,7 +1981,9 @@ pliston.AddToList(&tliston); InitBinnings(&pliston); + + pliston.AddToList(&supercuts); + pliston.AddToList(&binsalphaabs); pliston.AddToList(&alphaabs); - pliston.AddToList(&supercuts); //***************************** @@ -1727,4 +1994,6 @@ tliston.AddToList(&sccalc); tliston.AddToList(&fillhadsc); + + tliston.AddToList(&write); tliston.AddToList(&contfinalgh); @@ -1735,7 +2004,4 @@ tliston.AddToList(&hfill4); tliston.AddToList(&hfill5); - - if (WSC) - tliston.AddToList(&write); tliston.AddToList(&contfinal); @@ -1776,9 +2042,10 @@ TH1 &alphaHist = absalpha->GetHist(); alphaHist.SetXTitle("|alpha| [\\circ]"); + alphaHist.SetName("alpha-macro"); Double_t alphasig = 13.1; Double_t alphamin = 30.0; Double_t alphamax = 90.0; - Int_t degree = 4; + Int_t degree = 2; Double_t significance = -99.0; Bool_t drawpoly = kTRUE; @@ -1788,5 +2055,5 @@ MHFindSignificance findsig; findsig.SetRebin(kTRUE); - findsig.SetReduceDegree(kTRUE); + findsig.SetReduceDegree(kFALSE); findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, @@ -1805,4 +2072,6 @@ DeleteBinnings(&pliston); + } + gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl; @@ -1810,661 +2079,4 @@ } //--------------------------------------------------------------------- - - - //--------------------------------------------------------------------- - // Job B_RF_UP - //============ - - - // - create (or read in) the matrices of look-alike events for gammas - // and hadrons - // - create (or read in) the trees - // - then read ON1.root (or MC1.root) file - // - calculate the hadroness for the method of RANDOM FOREST - // - update input root file with the hadroness - - - if (JobB_RF_UP) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " - << JobB_RF_UP << ", " << RLookRF << ", " << RTree << ", " - << WRF << endl; - - - - //-------------------------------------------- - // parameters for the random forest - Int_t NumTrees = 100; - Int_t NumTry = 3; - Int_t NdSize = 3; - - - //-------------------------------------------- - // file to be updated (either ON or MC) - - TString typeInput = "ON"; - //TString typeInput = "MC"; - gLog << "typeInput = " << typeInput << endl; - - // name of input root file - TString filenameData = outPath; - filenameData += typeInput; - filenameData += "2.root"; - gLog << "filenameData = " << filenameData << endl; - - // name of output root file - TString outNameImage = outPath; - outNameImage += typeInput; - outNameImage += "3.root"; - //TString outNameImage = filenameData; - - gLog << "outNameImage = " << outNameImage << endl; - - //-------------------------------------------- - // files to be read for generating the look-alike events - // "hadrons" : - TString filenameON = outPath; - filenameON += "ON"; - filenameON += "1.root"; - Int_t howManyHadrons = 1000000; - gLog << "filenameON = " << filenameON << ", howManyHadrons = " - << howManyHadrons << endl; - - - // "gammas" : - TString filenameMC = outPath; - filenameMC += "MC"; - filenameMC += "1.root"; - Int_t howManyGammas = 50000; - gLog << "filenameMC = " << filenameMC << ", howManyGammas = " - << howManyGammas << endl; - - //-------------------------------------------- - // files of look-alike events - - TString outNameGammas = outPath; - outNameGammas += "RFmatrix_gammas_"; - outNameGammas += "MC"; - outNameGammas += ".root"; - - TString typeMatrixHadrons = "ON"; - gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; - - TString outNameHadrons = outPath; - outNameHadrons += "RFmatrix_hadrons_"; - outNameHadrons += typeMatrixHadrons; - outNameHadrons += ".root"; - - - MHMatrix matrixg("MatrixGammas"); - MHMatrix matrixh("MatrixHadrons"); - - //-------------------------------------------- - // file of trees of the random forest - - TString outRF = outPath; - outRF += "RF.root"; - - - //************************************************************************* - // read in matrices of look-alike events -if (RLookRF) - { - const char* mtxName = "MatrixGammas"; - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Get matrix for (gammas)" << endl; - gLog << "matrix name = " << mtxName << endl; - gLog << "name of root file = " << outNameGammas << endl; - gLog << "" << endl; - - - // read in the object with the name 'mtxName' from file 'outNameGammas' - // - TFile fileg(outNameGammas); - - matrixg.Read(mtxName); - matrixg.Print("cols"); - - - //***************************************************************** - - const char* mtxName = "MatrixHadrons"; - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << " Get matrix for (hadrons)" << endl; - gLog << "matrix name = " << mtxName << endl; - gLog << "name of root file = " << outNameHadrons << endl; - gLog << "" << endl; - - - // read in the object with the name 'mtxName' from file 'outNameHadrons' - // - TFile fileh(outNameHadrons); - - matrixh.Read(mtxName); - matrixh.Print("cols"); - } - - - //************************************************************************* - // create matrices of look-alike events -if (!RLookRF) - { - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << " Create matrices of look-alike events" << endl; - gLog << " Gammas :" << endl; - - - MParList plistg; - MTaskList tlistg; - MFilterList flistg; - - MParList plisth; - MTaskList tlisth; - MFilterList flisth; - - MReadMarsFile readg("Events", filenameMC); - readg.DisableAutoScheme(); - - MReadMarsFile readh("Events", filenameON); - readh.DisableAutoScheme(); - - MFParticleId fgamma("MMcEvt", '=', kGAMMA); - fgamma.SetName("gammaID"); - MFParticleId fhadrons("MMcEvt", '!', kGAMMA); - fhadrons.SetName("hadronID)"); - - MFEventSelector selectorg; - selectorg.SetNumSelectEvts(howManyGammas); - selectorg.SetName("selectGammas"); - MFEventSelector selectorh; - selectorh.SetNumSelectEvts(howManyHadrons); - selectorh.SetName("selectHadrons"); - - MFillH fillmatg("MatrixGammas"); - fillmatg.SetFilter(&flistg); - fillmatg.SetName("fillGammas"); - - MFillH fillmath("MatrixHadrons"); - fillmath.SetFilter(&flisth); - fillmath.SetName("fillHadrons"); - - - //***************************** fill gammas *** - // entries in MFilterList - - flistg.AddToList(&fgamma); - flistg.AddToList(&selectorg); - - //***************************** - // entries in MParList - - plistg.AddToList(&tlistg); - InitBinnings(&plistg); - - plistg.AddToList(&matrixg); - - //***************************** - // entries in MTaskList - - tlistg.AddToList(&readg); - tlistg.AddToList(&flistg); - tlistg.AddToList(&fillmatg); - - //***************************** - - MProgressBar matrixbar; - MEvtLoop evtloopg; - evtloopg.SetParList(&plistg); - evtloopg.ReadEnv(env, "", printEnv); - evtloopg.SetProgressBar(&matrixbar); - - Int_t maxevents = -1; - if (!evtloopg.Eventloop(maxevents)) - return; - - tlistg.PrintStatistics(0, kTRUE); - - - //***************************** fill hadrons *** - - gLog << " Hadrons :" << endl; - // entries in MFilterList - - flisth.AddToList(&fhadrons); - flisth.AddToList(&selectorh); - - //***************************** - // entries in MParList - - plisth.AddToList(&tlisth); - InitBinnings(&plisth); - - plisth.AddToList(&matrixh); - - //***************************** - // entries in MTaskList - - tlisth.AddToList(&readh); - tlisth.AddToList(&flisth); - tlisth.AddToList(&fillmath); - - //***************************** - - MProgressBar matrixbar; - MEvtLoop evtlooph; - evtlooph.SetParList(&plisth); - evtlooph.ReadEnv(env, "", printEnv); - evtlooph.SetProgressBar(&matrixbar); - - Int_t maxevents = -1; - if (!evtlooph.Eventloop(maxevents)) - return; - - tlisth.PrintStatistics(0, kTRUE); - //***************************************************** - - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - // - // ---------------------------------------------------------- - // Definition of the reference sample of look-alike events - // (this is a subsample of the original sample) - // ---------------------------------------------------------- - // - gLog << "" << endl; - gLog << "========================================================" << endl; - // Select a maximum of nmaxevts events from the sample of look-alike - // events. They will form the reference sample. - Int_t nmaxevts = 10000; - - // target distribution for the variable in column refcolumn (start at 0); - //Int_t refcolumn = 0; - //Int_t nbins = 5; - //Float_t frombin = 0.5; - //Float_t tobin = 1.0; - //TH1F *thsh = new TH1F("thsh","target distribution", - // nbins, frombin, tobin); - //thsh->SetDirectory(NULL); - //thsh->SetXTitle("cos( \\Theta)"); - //thsh->SetYTitle("Counts"); - //Float_t dbin = (tobin-frombin)/nbins; - //Float_t lbin = frombin +dbin*0.5; - //for (Int_t j=1; j<=nbins; j++) - //{ - // thsh->Fill(lbin,1.0); - // lbin += dbin; - //} - - Int_t refcolumn = 0; - MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); - TH1F *thsh = new TH1F(); - thsh->SetNameTitle("thsh","target distribution"); - MH::SetBinning(thsh, binscostheta); - Int_t nbins = thsh->GetNbinsX(); - Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; - for (Int_t j=1; j<=nbins; j++) - { - thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixg.EnableGraphicalOutput(); - Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; - return; - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixh.EnableGraphicalOutput(); - matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); - delete thsh; - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; - return; - } - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - - // write out look-alike events - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Write out look=alike events" << endl; - - - //------------------------------------------- - // "gammas" - gLog << "Gammas :" << endl; - matrixg.Print("cols"); - - TFile writeg(outNameGammas, "RECREATE", ""); - matrixg.Write(); - - gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " - << outNameGammas << endl; - - //------------------------------------------- - // "hadrons" - gLog << "Hadrons :" << endl; - matrixh.Print("cols"); - - TFile writeh(outNameHadrons, "RECREATE", ""); - matrixh.Write(); - - gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " - << outNameHadrons << endl; - - } - //********** end of creating matrices of look-alike events *********** - - - MRanForest *fRanForest; - MRanTree *fRanTree; - //----------------------------------------------------------------- - // read in the trees of the random forest - if (RTree) - { - MParList plisttr; - MTaskList tlisttr; - plisttr.AddToList(&tlisttr); - - MReadTree readtr("TREE", outRF); - readtr.DisableAutoScheme(); - - MRanForestFill rffill; - rffill.SetNumTrees(NumTrees); - - // list of tasks for the loop over the trees - - tlisttr.AddToList(&readtr); - tlisttr.AddToList(&rffill); - - //------------------- - // Execute tree loop - // - MEvtLoop evtlooptr; - evtlooptr.SetParList(&plisttr); - if (!evtlooptr.Eventloop()) - return; - - tlisttr.PrintStatistics(0, kTRUE); - - - // get adresses of objects which are used in the next eventloop - fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); - if (!fRanForest) - { - *fLog << err << dbginf << "MRanForest not found... aborting." << endl; - return kFALSE; - } - - fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); - if (!fRanTree) - { - *fLog << err << dbginf << "MRanTree not found... aborting." << endl; - return kFALSE; - } - - } - - //----------------------------------------------------------------- - // grow the trees of the random forest (event loop = tree loop) - - if (!RTree) - { - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : start growing trees" << endl; - - MTaskList tlist2; - MParList plist2; - plist2.AddToList(&tlist2); - - plist2.AddToList(&matrixg); - plist2.AddToList(&matrixh); - - MRanForestGrow rfgrow2; - rfgrow2.SetNumTrees(NumTrees); - rfgrow2.SetNumTry(NumTry); - rfgrow2.SetNdSize(NdSize); - - MWriteRootFile rfwrite2(outRF); - rfwrite2.AddContainer("MRanTree", "TREE"); - - // list of tasks for the loop over the trees - - tlist2.AddToList(&rfgrow2); - tlist2.AddToList(&rfwrite2); - - //------------------- - // Execute tree loop - // - MEvtLoop treeloop; - treeloop.SetParList(&plist2); - - if ( !treeloop.Eventloop() ) - return; - - tlist2.PrintStatistics(0, kTRUE); - - // get adresses of objects which are used in the next eventloop - fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); - if (!fRanForest) - { - *fLog << err << dbginf << "MRanForest not found... aborting." << endl; - return kFALSE; - } - - fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); - if (!fRanTree) - { - *fLog << err << dbginf << "MRanTree not found... aborting." << endl; - return kFALSE; - } - - } - // end of growing the trees of the random forest - //----------------------------------------------------------------- - - - - - //----------------------------------------------------------------- - // Update the input files with the RF hadronness - // - - if (WRF) - { - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Update input file '" << filenameData - << "' with the RF hadronness" << endl; - - MTaskList tliston; - MParList pliston; - - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //....................................................................... - // calculate hadronnes for method of RANDOM FOREST - - TString hadRFName = "HadRF"; - MRanForestCalc rfcalc; - rfcalc.SetHadronnessName(hadRFName); - - //....................................................................... - - //MWriteRootFile write(outNameImage, "UPDATE"); - MWriteRootFile write(outNameImage, "RECREATE"); - - write.AddContainer("MRawRunHeader", "RunHeaders"); - write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("ThetaOrig", "Events"); - write.AddContainer("MSrcPosCam", "Events"); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("MHillas", "Events"); - write.AddContainer("MHillasExt", "Events"); - write.AddContainer("MHillasSrc", "Events"); - write.AddContainer("MNewImagePar", "Events"); - - write.AddContainer("HadNN", "Events"); - write.AddContainer("HadSC", "Events"); - - write.AddContainer(hadRFName, "Events"); - - - //----------------------------------------------------------------- - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - - MFCT1SelFinal selfinalgh(fHilNameSrc); - selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadRFName); - selfinalgh.SetName("SelFinalgh"); - MContinue contfinalgh(&selfinalgh); - contfinalgh.SetName("ContSelFinalgh"); - - MFillH fillranfor("MHRanForest"); - fillranfor.SetName("HRanForest"); - - MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); - fillhadrf.SetName("HhadRF"); - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadRFName); - selfinal.SetName("SelFinal"); - MContinue contfinal(&selfinal); - contfinal.SetName("ContSelFinal"); - - - MFillH hfill1("MHHillas", fHilName); - hfill1.SetName("HHillas"); - - MFillH hfill2("MHStarMap", fHilName); - hfill2.SetName("HStarMap"); - - MFillH hfill3("MHHillasExt", fHilNameSrc); - hfill3.SetName("HHillasExt"); - - MFillH hfill4("MHHillasSrc", fHilNameSrc); - hfill4.SetName("HHillasSrc"); - - MFillH hfill5("MHNewImagePar", fImgParName); - hfill5.SetName("HNewImagePar"); - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - pliston.AddToList(fRanForest); - pliston.AddToList(fRanTree); - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - - tliston.AddToList(&rfcalc); - tliston.AddToList(&fillranfor); - tliston.AddToList(&fillhadrf); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - tliston.AddToList(&hfill5); - - tliston.AddToList(&write); - - tliston.AddToList(&contfinalgh); - tliston.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - pliston.FindObject("MHRanForest")->DrawClone(); - pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); - - pliston.FindObject("MHHillas")->DrawClone(); - pliston.FindObject("MHHillasExt")->DrawClone(); - pliston.FindObject("MHHillasSrc")->DrawClone(); - pliston.FindObject("MHNewImagePar")->DrawClone(); - pliston.FindObject("MHStarMap")->DrawClone(); - - DeleteBinnings(&pliston); - } - - gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - @@ -2484,5 +2096,6 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobC = " << JobC << endl; + gLog << "Macro CT1Analysis : JobC = " + << (JobC ? "kTRUE" : "kFALSE") << endl; @@ -2522,5 +2135,5 @@ // names of hadronness containers - TString hadNNName = "HadNN"; + //TString hadNNName = "HadNN"; TString hadSCName = "HadSC"; TString hadRFName = "HadRF"; @@ -2540,11 +2153,11 @@ MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadNNName); + selfinalgh.SetHadronnessName(hadSCName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); - MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); - fillhadnn.SetName("HhadNN"); + //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + //fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); fillhadsc.SetName("HhadSC"); @@ -2554,5 +2167,5 @@ MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadNNName); + selfinal.SetHadronnessName(hadSCName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); @@ -2588,5 +2201,5 @@ tliston.AddToList(&read); - tliston.AddToList(&fillhadnn); + //tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); @@ -2623,5 +2236,5 @@ // - pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); @@ -2657,5 +2270,6 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobD = " - << JobD << endl; + << (JobD ? "kTRUE" : "kFALSE") << endl; + // type of data to be analysed @@ -2868,5 +2482,7 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobE_XX, WXX = " - << JobE_XX << ", " << WXX << endl; + << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (WXX ? "kTRUE" : "kFALSE") << endl; + // type of data to be analysed @@ -3060,4 +2676,7 @@ gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; + + if (WXX) + { //----------------------------------------------------------- // @@ -3119,7 +2738,7 @@ //....................................................................... - if (WXX) - { - MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + + //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + MWriteRootFile write(filenameDataout, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -3139,5 +2758,5 @@ write.AddContainer("MEnergyEst", "Events"); - } + //----------------------------------------------------------------- @@ -3265,4 +2884,6 @@ DeleteBinnings(&pliston); + + } gLog << "Macro CT1Analysis : End of Job E_XX" << endl; @@ -3292,5 +2913,7 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobF_XX, WFX = " - << JobF_XX << ", " << WFX << endl; + << (JobF_XX ? "kTRUE" : "kFALSE") << ", " + << (WFX ? "kTRUE" : "kFALSE") << endl; + // type of data to be analysed @@ -3493,4 +3116,7 @@ gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; + + if (WFX) + { //----------------------------------------------------------- // @@ -3552,10 +3178,9 @@ //....................................................................... - //if (WFX) - //{ gLog << "CT1Analysis.C : write root file '" << filenameDataout << "'" << endl; //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + /* MWriteRootFile write(filenameDataout, "RECREATE"); @@ -3571,10 +3196,10 @@ write.AddContainer("MNewImagePar", "Events"); - write.AddContainer("HadNN", "Events"); + //write.AddContainer("HadNN", "Events"); write.AddContainer("HadSC", "Events"); write.AddContainer("HadRF", "Events"); write.AddContainer("MEnergyEst", "Events"); - //} + */ //----------------------------------------------------------------- @@ -3590,6 +3215,6 @@ contfinalgh.SetName("ContSelFinalgh"); - MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); - fillhadnn.SetName("HhadNN"); + //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); + //fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); fillhadsc.SetName("HhadSC"); @@ -3600,15 +3225,15 @@ // calculate estimated energy - //MEnergyEstParam eeston(fHilName); - //eeston.Add(fHilNameSrc); - - //eeston.SetCoeffA(parA); - //eeston.SetCoeffB(parB); + MEnergyEstParam eeston(fHilName); + eeston.Add(fHilNameSrc); + + eeston.SetCoeffA(parA); + eeston.SetCoeffB(parB); //--------------------------- // calculate estimated energy using Daniel's parameters - MEnergyEstParamDanielMkn421 eeston(fHilName); - eeston.Add(fHilNameSrc); + //MEnergyEstParamDanielMkn421 eeston(fHilName); + //eeston.Add(fHilNameSrc); //eeston.SetCoeffA(parA); //eeston.SetCoeffB(parB); @@ -3684,8 +3309,7 @@ tliston.AddToList(&eeston); - if (WFX) - tliston.AddToList(&write); - - tliston.AddToList(&fillhadnn); + //tliston.AddToList(&write); + + //tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); @@ -3727,19 +3351,38 @@ // - pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + + gLog << "before hadRF" << endl; pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + gLog << "before hadSC" << endl; pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + gLog << "before MHHillas" << endl; pliston.FindObject("MHHillas")->DrawClone(); + + gLog << "before MHHillasExt" << endl; pliston.FindObject("MHHillasExt")->DrawClone(); + + gLog << "before MHHillasSrc" << endl; pliston.FindObject("MHHillasSrc")->DrawClone(); + + gLog << "before MHNewImagePar" << endl; pliston.FindObject("MHNewImagePar")->DrawClone(); + + gLog << "before MHStarMap" << endl; pliston.FindObject("MHStarMap")->DrawClone(); + gLog << "before DeleteBinnings" << endl; + //DeleteBinnings(&pliston); + gLog << "before Robert's code" << endl; + //rwagner write all relevant histograms onto a file + if (WRobert) + { gLog << "=======================================================" << endl; gLog << "Write results onto file '" << filenameResults << "'" << endl; @@ -3751,5 +3394,7 @@ TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); alphaHist->Write(); - + gLog << "Alpha plot has been written out" << endl; + + MHAlphaEnergyTheta* aetH = (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); @@ -3757,4 +3402,5 @@ aetHist->SetName("aetHist"); aetHist->Write(); + gLog << "AlphaEnergyTheta plot has been written out" << endl; MHAlphaEnergyTime* aetH2 = @@ -3764,4 +3410,5 @@ // aetHist2->DrawClone(); aetHist2->Write(); + gLog << "AlphaEnergyTime plot has been written out" << endl; MHThetabarTime* tbt = @@ -3770,4 +3417,5 @@ tbtHist->SetName("tbtHist"); tbtHist->Write(); + gLog << "ThetabarTime plot has been written out" << endl; MHEnergyTime* ent = @@ -3776,4 +3424,5 @@ entHist->SetName("entHist"); entHist->Write(); + gLog << "EnergyTime plot has been written out" << endl; MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); @@ -3782,4 +3431,6 @@ timeHist->SetTitle("Time diffs"); timeHist->Write(); + gLog << "TimeDiffTheta plot has been written out" << endl; + MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); @@ -3788,4 +3439,5 @@ timeHist2->SetTitle("Time diffs"); timeHist2->Write(); + gLog << "TimeDiffTime plot has been written out" << endl; //rwagner write also collareas to same file @@ -3793,12 +3445,18 @@ collarea->GetHAll()->Write(); collarea->GetHSel()->Write(); - + gLog << "Effective collection areas have been written out" << endl; + //rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data //rwagner (ON/OFF/MC), MJDmin, MJDmax to this file + gLog << "before closing outfile" << endl; + //outfile.Close(); - + gLog << "Results were written onto file '" << filenameResults + << "'" << endl; gLog << "=======================================================" << endl; - + } + + } gLog << "Macro CT1Analysis : End of Job F_XX" << endl; Index: trunk/MagicSoft/Mars/macros/CT1EgyEst.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2356) +++ trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2357) @@ -33,5 +33,5 @@ #include "MHMatrix.h" #include "MEnergyEstParam.h" -#include "MEnergyEstParamDanielMkn421.h" +//#include "MEnergyEstParamDanielMkn421.h" #include "MMatrixLoop.h" #include "MChisqEval.h" @@ -206,13 +206,13 @@ // - //MEnergyEstParam eest2(hilName); + MEnergyEstParam eest2(hilName); + eest2.Add(hilSrcName); + + eest2.SetCoeffA(parA); + eest2.SetCoeffB(parB); + + // estimate energy using Daniel's parameters + //MEnergyEstParamDanielMkn421 eest2(hilName); //eest2.Add(hilSrcName); - - //eest2.SetCoeffA(parA); - //eest2.SetCoeffB(parB); - - // estimate energy using Daniel's parameters - MEnergyEstParamDanielMkn421 eest2(hilName); - eest2.Add(hilSrcName); Index: trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2357) @@ -57,5 +57,6 @@ #include "MFEventSelector2.h" #include "MFillH.h" -#include "MGeomCamCT1Daniel.h" +//#include "MGeomCamCT1Daniel.h" +#include "MGeomCamCT1.h" #include "MH3.h" #include "MHCT1Supercuts.h" @@ -95,4 +96,6 @@ Double_t *par, Int_t iflag) { + //cout << "entry fcnSupercuts" << endl; + //------------------------------------------------------------- // save pointer to the MINUIT object for optimizing the supercuts @@ -106,4 +109,5 @@ MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); + //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts"); @@ -130,4 +134,6 @@ // for testing //TArrayD checkparameters = super->GetParameters(); + //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " + // << fNparx << endl; //gLog << "fcnsupercuts : i, par, checkparameters =" << endl; //for (Int_t i=0; iEventloop(); + //tasklistfcn->PrintStatistics(0, kTRUE); + MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); if (!alpha) @@ -147,5 +155,5 @@ TH1 &alphaHist = alpha->GetHist(); - //alphaHist.SetXTitle("|\\alpha| [\\circ]"); + alphaHist.SetName("alpha-fcnSupercuts"); //------------------------------------------- @@ -161,5 +169,5 @@ const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; - const Int_t degree = 4; + const Int_t degree = 2; Bool_t drawpoly; @@ -240,7 +248,8 @@ //--------------------------- // camera geometry is needed for conversion mm ==> degree - fCam = new MGeomCamCT1Daniel; - - // matrices to contain the the training/test samples + //fCam = new MGeomCamCT1Daniel; + fCam = new MGeomCamCT1; + + // matrices to contain the training/test samples fMatrixTrain = new MHMatrix("MatrixTrain"); fMatrixTest = new MHMatrix("MatrixTest"); @@ -277,5 +286,6 @@ // Bool_t MCT1FindSupercuts::DefineTrainMatrix(const TString &nametrain, - const Int_t howmanytrain, MH3 &hreftrain) + const Int_t howmanytrain, MH3 &hreftrain, + const TString &filetrain) { if (nametrain.IsNull() || howmanytrain <= 0) @@ -343,4 +353,18 @@ *fLog << "=============================================" << endl; + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + return kTRUE; } @@ -356,5 +380,6 @@ // Bool_t MCT1FindSupercuts::DefineTestMatrix(const TString &nametest, - const Int_t howmanytest, MH3 &hreftest) + const Int_t howmanytest, MH3 &hreftest, + const TString &filetest) { if (nametest.IsNull() || howmanytest<=0) @@ -422,4 +447,18 @@ *fLog << "=============================================" << endl; + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + return kTRUE; } @@ -434,5 +473,6 @@ const TString &name, const Int_t howmanytrain, MH3 &hreftrain, - const Int_t howmanytest, MH3 &hreftest) + const Int_t howmanytest, MH3 &hreftest, + const TString &filetrain, const TString &filetest) { *fLog << "=============================================" << endl; @@ -511,5 +551,6 @@ evtloop.SetProgressBar(&bar); - if (!evtloop.Eventloop()) + Int_t maxev = -1; + if (!evtloop.Eventloop(maxev)) return kFALSE; @@ -520,5 +561,5 @@ if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { - *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events (" + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtrain << ") is incompatible with the no.of requested events (" @@ -530,5 +571,5 @@ if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) { - *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events (" + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtest << ") is incompatible with the no.of requested events (" @@ -541,4 +582,30 @@ + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + return kTRUE; } @@ -548,6 +615,4 @@ // Read training and test matrices from files // -// -// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ // @@ -563,4 +628,5 @@ *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '" << filetrain << "'" << endl; + filetr.Close(); @@ -574,45 +640,9 @@ *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '" << filetest << "'" << endl; - + filete.Close(); return kTRUE; } -// -------------------------------------------------------------------------- -// -// Write training and test matrices onto files -// -// -// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ -// -// -Bool_t MCT1FindSupercuts::WriteMatrix(const TString &filetrain, const TString &filetest) -{ - //-------------------------- - // write out training matrix - - TFile filetr(filetrain, "RECREATE"); - - *fLog << "nach TFile" << endl; - - fMatrixTrain->Write(); - - *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '" - << filetrain << "'" << endl; - - - //-------------------------- - // write out test matrix - - TFile filete(filetest, "RECREATE"); - fMatrixTest->Print("SizeCols"); - fMatrixTest->Write(); - - *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '" - << filetest << "'" << endl; - - - return kTRUE; -} //------------------------------------------------------------------------ @@ -645,13 +675,18 @@ // put the supercuts hadronness // -// - for the minimization, the starting values of the parameters are taken as -// MCT1Supercuts::GetParameters(fVinit) +// - for the minimization, the starting values of the parameters are taken +// - from file parSCinit (if it is != "") +// - or from the arrays params and/or steps // //---------------------------------------------------------------------- -Bool_t MCT1FindSupercuts::FindParams() +Bool_t MCT1FindSupercuts::FindParams(TString parSCinit, + TArrayD ¶ms, TArrayD &steps) { // Setup the event loop which will be executed in the // fcnSupercuts function of MINUIT // + // parSCinit is the name of the file containing the initial values + // of the parameters; + // if parSCinit = "" 'params' and 'steps' are taken as initial values *fLog << "=============================================" << endl; @@ -685,8 +720,43 @@ MMatrixLoop loop(fMatrixTrain); + //-------------------------------- // create container for the supercut parameters - // and set them to their default values + // and set them to their initial values MCT1Supercuts super; + // take initial values from file parSCinit + if (parSCinit != "") + { + TFile inparam(parSCinit); + super.Read("MCT1Supercuts"); + inparam.Close(); + *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file " + << parSCinit << endl; + } + + // take initial values from 'params' and/or 'steps' + else if (params.GetSize() != 0 || steps.GetSize() != 0 ) + { + if (params.GetSize() != 0) + { + *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'" + << endl; + super.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'" + << endl; + super.SetStepsizes(steps); + } + } + else + { + *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor" + << endl; + } + + + //-------------------------------- // calculate supercuts hadronness fCalcHadTrain->SetHadronnessName(fHadronnessName); @@ -761,6 +831,7 @@ // - // get initial values of parameters from MCT1Supercuts + // get initial values of parameters fVinit = super.GetParameters(); + fStep = super.GetStepsizes(); TString name[fVinit.GetSize()]; @@ -776,16 +847,45 @@ name[i] = "p"; name[i] += i+1; - fStep[i] = TMath::Abs(fVinit[i]/10.0); + //fStep[i] = TMath::Abs(fVinit[i]/10.0); fLimlo[i] = -100.0; fLimup[i] = 100.0; fFix[i] = 0; - - // vary only first 48 parameters - if (i >= 48) - { - fStep[i] = 0.0; - fFix[i] = 1; - } - } + } + + // these parameters make no sense, fix them at 0.0 + fVinit[33] = 0.0; + fStep[33] = 0.0; + fFix[33] = 1; + + fVinit[36] = 0.0; + fStep[36] = 0.0; + fFix[36] = 1; + + fVinit[39] = 0.0; + fStep[39] = 0.0; + fFix[39] = 1; + + fVinit[41] = 0.0; + fStep[41] = 0.0; + fFix[41] = 1; + + fVinit[44] = 0.0; + fStep[44] = 0.0; + fFix[44] = 1; + + fVinit[47] = 0.0; + fStep[47] = 0.0; + fFix[47] = 1; + + // vary only first 48 parameters + //for (UInt_t i=0; i= 48) + // { + // fStep[i] = 0.0; + // fFix[i] = 1; + // } + //} + // ------------------------------------------- @@ -910,4 +1010,8 @@ MH3 alphabefore("MatrixTest[7]"); alphabefore.SetName(mh3NameB); + + TH1 &alphahistb = alphabefore.GetHist(); + alphahistb.SetName("alphaBefore-TestParams"); + MFillH fillalphabefore(&alphabefore); fillalphabefore.SetName("FillAlphaBefore"); @@ -924,4 +1028,8 @@ MH3 alphaafter("MatrixTest[7]"); alphaafter.SetName(mh3NameA); + + TH1 &alphahista = alphaafter.GetHist(); + alphahista.SetName("alphaAfter-TestParams"); + MFillH fillalphaafter(&alphaafter); fillalphaafter.SetName("FillAlphaAfter"); @@ -978,4 +1086,5 @@ TH1 &alphaHist2 = alphaAfter->GetHist(); alphaHist2.SetXTitle("|\\alpha| [\\circ]"); + alphaHist2.SetName("alpha-TestParams)"); TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); @@ -999,5 +1108,5 @@ const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; - const Int_t degree = 4; + const Int_t degree = 2; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kFALSE; @@ -1025,2 +1134,3 @@ return kTRUE; } + Index: trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2357) @@ -89,17 +89,20 @@ Bool_t DefineTrainMatrix(const TString &name, const Int_t howmany, - MH3 &href); + MH3 &href, const TString &filetrain); Bool_t DefineTestMatrix(const TString &name, const Int_t howmany, - MH3 &href); + MH3 &href, const TString &filetest); Bool_t DefineTrainTestMatrix(const TString &name, - const Int_t howmanytrain, MH3 &hreftrain, - const Int_t howmanytest, MH3 &hreftest); + const Int_t howmanytrain, MH3 &hreftrain, + const Int_t howmanytest, MH3 &hreftest, + const TString &filetrain, const TString &filetest); + + MHMatrix *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); - Bool_t WriteMatrix(const TString &filetrain, const TString &filetest); - Bool_t FindParams(); + Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); Bool_t TestParams(); Index: trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2357) @@ -50,10 +50,12 @@ // MCT1Supercuts::MCT1Supercuts(const char *name, const char *title) - : fParameters(72), - fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), + : fParameters(104), fStepsizes(104), + fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24), - fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), - fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), - fAlphaUp(fParameters.GetArray()+64) + fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), + fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), + fConcUp(fParameters.GetArray()+64), fConcLo(fParameters.GetArray()+72), + fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88), + fAlphaUp(fParameters.GetArray()+96) { fName = name ? name : "MCT1Supercuts"; @@ -71,4 +73,7 @@ void MCT1Supercuts::InitParameters() { + //--------------------------------------------------- + // these are Daniel's original values for Mkn 421 + fLengthUp[0] = 0.315585; fLengthUp[1] = 0.001455; @@ -80,4 +85,13 @@ fLengthUp[7] = -0.013463; + fLengthLo[0] = 0.151530; + fLengthLo[1] = 0.028323; + fLengthLo[2] = 0.510707; + fLengthLo[3] = 0.053089; + fLengthLo[4] = 0.013708; + fLengthLo[5] = 2.357993; + fLengthLo[6] = 0.000080; + fLengthLo[7] = -0.007157; + fWidthUp[0] = 0.145412; fWidthUp[1] = -0.001771; @@ -89,4 +103,13 @@ fWidthUp[7] = -0.016703; + fWidthLo[0] = 0.089187; + fWidthLo[1] = -0.006430; + fWidthLo[2] = 0.074442; + fWidthLo[3] = 0.003738; + fWidthLo[4] = -0.004256; + fWidthLo[5] = -0.014101; + fWidthLo[6] = 0.006126; + fWidthLo[7] = -0.002849; + fDistUp[0] = 1.787943; fDistUp[1] = 0; @@ -98,22 +121,4 @@ fDistUp[7] = 0; - fLengthLo[0] = 0.151530; - fLengthLo[1] = 0.028323; - fLengthLo[2] = 0.510707; - fLengthLo[3] = 0.053089; - fLengthLo[4] = 0.013708; - fLengthLo[5] = 2.357993; - fLengthLo[6] = 0.000080; - fLengthLo[7] = -0.007157; - - fWidthLo[0] = 0.089187; - fWidthLo[1] = -0.006430; - fWidthLo[2] = 0.074442; - fWidthLo[3] = 0.003738; - fWidthLo[4] = -0.004256; - fWidthLo[5] = -0.014101; - fWidthLo[6] = 0.006126; - fWidthLo[7] = -0.002849; - fDistLo[0] = 0.589406; fDistLo[1] = 0; @@ -124,22 +129,61 @@ fDistLo[6] = -0.001750; fDistLo[7] = 0; - - fAsymUp[0] = 0.061267; - fAsymUp[1] = 0.014462; - fAsymUp[2] = 0.014327; - fAsymUp[3] = 0.014540; - fAsymUp[4] = 0.013391; - fAsymUp[5] = 0.012319; - fAsymUp[6] = 0.010444; - fAsymUp[7] = 0.008328; - - fAsymLo[0] = -0.012055; - fAsymLo[1] = 0.009157; - fAsymLo[2] = 0.005441; - fAsymLo[3] = 0.000399; - fAsymLo[4] = 0.001433; - fAsymLo[5] = -0.002050; - fAsymLo[6] = -0.000104; - fAsymLo[7] = -0.001188; + + + // dummy values + + fAsymUp[0] = 1.e10; + fAsymUp[1] = 0.0; + fAsymUp[2] = 0.0; + fAsymUp[3] = 0.0; + fAsymUp[4] = 0.0; + fAsymUp[5] = 0.0; + fAsymUp[6] = 0.0; + fAsymUp[7] = 0.0; + + fAsymLo[0] = -1.e10; + fAsymLo[1] = 0.0; + fAsymLo[2] = 0.0; + fAsymLo[3] = 0.0; + fAsymLo[4] = 0.0; + fAsymLo[5] = 0.0; + fAsymLo[6] = 0.0; + fAsymLo[7] = 0.0; + + fConcUp[0] = 1.e10; + fConcUp[1] = 0.0; + fConcUp[2] = 0.0; + fConcUp[3] = 0.0; + fConcUp[4] = 0.0; + fConcUp[5] = 0.0; + fConcUp[6] = 0.0; + fConcUp[7] = 0.0; + + fConcLo[0] = -1.e10; + fConcLo[1] = 0.0; + fConcLo[2] = 0.0; + fConcLo[3] = 0.0; + fConcLo[4] = 0.0; + fConcLo[5] = 0.0; + fConcLo[6] = 0.0; + fConcLo[7] = 0.0; + + fLeakage1Up[0] = 1.e10; + fLeakage1Up[1] = 0.0; + fLeakage1Up[2] = 0.0; + fLeakage1Up[3] = 0.0; + fLeakage1Up[4] = 0.0; + fLeakage1Up[5] = 0.0; + fLeakage1Up[6] = 0.0; + fLeakage1Up[7] = 0.0; + + fLeakage1Lo[0] = -1.e10; + fLeakage1Lo[1] = 0.0; + fLeakage1Lo[2] = 0.0; + fLeakage1Lo[3] = 0.0; + fLeakage1Lo[4] = 0.0; + fLeakage1Lo[5] = 0.0; + fLeakage1Lo[6] = 0.0; + fLeakage1Lo[7] = 0.0; fAlphaUp[0] = 13.123440; @@ -151,4 +195,139 @@ fAlphaUp[6] = 0; fAlphaUp[7] = 0; + + //--------------------------------------------------- + // fStepsizes + // if == 0.0 the parameter will be fixed in the minimization + // != 0.0 initial step sizes for the parameters + + // LengthUp + fStepsizes[0] = 0.03; + fStepsizes[1] = 0.0002; + fStepsizes[2] = 0.02; + fStepsizes[3] = 0.0006; + fStepsizes[4] = 0.0002; + fStepsizes[5] = 0.002; + fStepsizes[6] = 0.0008; + fStepsizes[7] = 0.002; + + // LengthLo + fStepsizes[8] = 0.02; + fStepsizes[9] = 0.003; + fStepsizes[10] = 0.05; + fStepsizes[11] = 0.006; + fStepsizes[12] = 0.002; + fStepsizes[13] = 0.3; + fStepsizes[14] = 0.0001; + fStepsizes[15] = 0.0008; + + // WidthUp + fStepsizes[16] = 0.02; + fStepsizes[17] = 0.0002; + fStepsizes[18] = 0.006; + fStepsizes[19] = 0.003; + fStepsizes[20] = 0.002; + fStepsizes[21] = 0.006; + fStepsizes[22] = 0.002; + fStepsizes[23] = 0.002; + + // WidthLo + fStepsizes[24] = 0.009; + fStepsizes[25] = 0.0007; + fStepsizes[26] = 0.008; + fStepsizes[27] = 0.0004; + fStepsizes[28] = 0.0005; + fStepsizes[29] = 0.002; + fStepsizes[30] = 0.0007; + fStepsizes[31] = 0.003; + + // DistUp + fStepsizes[32] = 0.2; + fStepsizes[33] = 0.0; + fStepsizes[34] = 0.3; + fStepsizes[35] = 0.02; + fStepsizes[36] = 0.0; + fStepsizes[37] = 0.03; + fStepsizes[38] = 0.02; + fStepsizes[39] = 0.0; + + // DistLo + fStepsizes[40] = 0.06; + fStepsizes[41] = 0.0; + fStepsizes[42] = 0.009; + fStepsizes[43] = 0.0008; + fStepsizes[44] = 0.0; + fStepsizes[45] = 0.005; + fStepsizes[46] = 0.0002; + fStepsizes[47] = 0.0; + + // AsymUp + fStepsizes[48] = 0.0; + fStepsizes[49] = 0.0; + fStepsizes[50] = 0.0; + fStepsizes[51] = 0.0; + fStepsizes[52] = 0.0; + fStepsizes[53] = 0.0; + fStepsizes[54] = 0.0; + fStepsizes[55] = 0.0; + + // AsymLo + fStepsizes[56] = 0.0; + fStepsizes[57] = 0.0; + fStepsizes[58] = 0.0; + fStepsizes[59] = 0.0; + fStepsizes[60] = 0.0; + fStepsizes[61] = 0.0; + fStepsizes[62] = 0.0; + fStepsizes[63] = 0.0; + + // ConcUp + fStepsizes[64] = 0.0; + fStepsizes[65] = 0.0; + fStepsizes[66] = 0.0; + fStepsizes[67] = 0.0; + fStepsizes[68] = 0.0; + fStepsizes[69] = 0.0; + fStepsizes[70] = 0.0; + fStepsizes[71] = 0.0; + + // ConcLo + fStepsizes[72] = 0.0; + fStepsizes[73] = 0.0; + fStepsizes[74] = 0.0; + fStepsizes[75] = 0.0; + fStepsizes[76] = 0.0; + fStepsizes[77] = 0.0; + fStepsizes[78] = 0.0; + fStepsizes[79] = 0.0; + + // Leakage1Up + fStepsizes[80] = 0.0; + fStepsizes[81] = 0.0; + fStepsizes[82] = 0.0; + fStepsizes[83] = 0.0; + fStepsizes[84] = 0.0; + fStepsizes[85] = 0.0; + fStepsizes[86] = 0.0; + fStepsizes[87] = 0.0; + + // Leakage1Lo + fStepsizes[88] = 0.0; + fStepsizes[89] = 0.0; + fStepsizes[90] = 0.0; + fStepsizes[91] = 0.0; + fStepsizes[92] = 0.0; + fStepsizes[93] = 0.0; + fStepsizes[94] = 0.0; + fStepsizes[95] = 0.0; + + // AlphaUp + fStepsizes[96] = 0.0; + fStepsizes[97] = 0.0; + fStepsizes[98] = 0.0; + fStepsizes[99] = 0.0; + fStepsizes[100] = 0.0; + fStepsizes[101] = 0.0; + fStepsizes[102] = 0.0; + fStepsizes[103] = 0.0; } @@ -173,13 +352,37 @@ } - - - - - - - - - - - +// -------------------------------------------------------------------------- +// +// Set the step sizes from the array 'd' +// +// +Bool_t MCT1Supercuts::SetStepsizes(const TArrayD &d) +{ + if (d.GetSize() != fStepsizes.GetSize()) + { + *fLog << err << "Sizes of d and of fStepsizes are different : " + << d.GetSize() << ", " << fStepsizes.GetSize() << endl; + return kFALSE; + } + + fStepsizes = d; + + return kTRUE; +} + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2357) @@ -14,4 +14,5 @@ private: TArrayD fParameters; // supercut parameters + TArrayD fStepsizes; // step sizes of supercut parameters Double_t *fLengthUp; //! @@ -23,5 +24,12 @@ Double_t *fAsymUp; //! Double_t *fAsymLo; //! + + Double_t *fConcUp; //! + Double_t *fConcLo; //! + Double_t *fLeakage1Up; //! + Double_t *fLeakage1Lo; //! + Double_t *fAlphaUp; //! + public: @@ -31,5 +39,8 @@ Bool_t SetParameters(const TArrayD &d); + Bool_t SetStepsizes(const TArrayD &d); + const TArrayD &GetParameters() const { return fParameters; } + const TArrayD &GetStepsizes() const { return fStepsizes; } const Double_t *GetLengthUp() const { return fLengthUp; } @@ -41,4 +52,11 @@ const Double_t *GetAsymUp() const { return fAsymUp; } const Double_t *GetAsymLo() const { return fAsymLo; } + + const Double_t *GetConcUp() const { return fConcUp; } + const Double_t *GetConcLo() const { return fConcLo; } + + const Double_t *GetLeakage1Up() const { return fLeakage1Up; } + const Double_t *GetLeakage1Lo() const { return fLeakage1Lo; } + const Double_t *GetAlphaUp() const { return fAlphaUp; } @@ -47,2 +65,5 @@ #endif + + + Index: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2357) @@ -44,4 +44,5 @@ #include "MHillasExt.h" #include "MHillasSrc.h" +#include "MNewImagePar.h" #include "MMcEvt.hxx" #include "MCerPhotEvt.h" @@ -65,7 +66,9 @@ MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, - const char *hilsrcname, const char *name, const char *title) + const char *hilsrcname, + const char *name, const char *title) : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname), - fSuperName("MCT1Supercuts") + fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), + fSuperName("MCT1Supercuts") { fName = name ? name : "MCT1SupercutsCalc"; @@ -101,5 +104,4 @@ return kFALSE; } - if (fMatrix) @@ -114,8 +116,22 @@ } + fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt"); + if (!fHilExt) + { + *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); if (!fHilSrc) { *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar"); + if (!fNewPar) + { + *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl; return kFALSE; } @@ -172,5 +188,13 @@ Double_t MCT1SupercutsCalc::GetVal(Int_t i) const { - return (*fMatrix)[fMap[i]]; + + Double_t val = (*fMatrix)[fMap[i]]; + + + //*fLog << "MCT1SupercutsCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + + return val; } @@ -187,5 +211,5 @@ { if (fMatrix) - return; + return; fMatrix = mat; @@ -199,10 +223,12 @@ fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); + fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); } // --------------------------------------------------------------------------- // -// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich) -// optimized for mkn 421 2001 data +// Evaluate dynamical supercuts // // set hadronness to 0.25 if cuts are fullfilled @@ -222,4 +248,12 @@ const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist(); + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), + fHilSrc->GetCosDeltaAlpha()); + const Double_t asym0 = fMatrix ? GetVal(8) : help; + const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); + const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + const Double_t newdist = dist0 * fMm2Deg; @@ -236,20 +270,67 @@ const Double_t length = length0 * fMm2Deg; const Double_t width = width0 * fMm2Deg; - - if (newdist < 1.05 && + const Double_t asym = asym0 * fMm2Deg; + + /* + *fLog << "newdist, length, width, asym, dist, conc, leakage = " + << newdist << ", " << length << ", " << width << ", " + << asym << ", " << dist << ", " << conc << ", " << leakage + << endl; + + *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", " + << endl; + */ + + + if ( + //dist < 1.05 && + //newdist < 1.05 && + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && - dist < 1.05 && + length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) && length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) && + width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) && width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) && - //asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && - //asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + + asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && - dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) ) - fHadronness->SetHadronness(0.25); + dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) ) + + fHadronness->SetHadronness(0.25); else - fHadronness->SetHadronness(0.75); + fHadronness->SetHadronness(0.75); + + //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl; fHadronness->SetReadyToSave(); @@ -257,2 +338,11 @@ return kTRUE; } + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2357) @@ -13,4 +13,6 @@ class MHillas; class MHillasSrc; +class MHillasExt; +class MNewImagePar; class MMcEvt; class MCerPhotEvt; @@ -25,4 +27,6 @@ MHillas *fHil; MHillasSrc *fHilSrc; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; MMcEvt *fMcEvt; MHadronness *fHadronness; //! output container for hadronness @@ -32,9 +36,11 @@ TString fHilName; TString fHilSrcName; + TString fHilExtName; + TString fNewParName; TString fSuperName; // name of container for supercut parameters Double_t fMm2Deg; //! - Int_t fMap[8]; //! + Int_t fMap[11]; //! MHMatrix *fMatrix; //! @@ -63,2 +69,14 @@ #endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2356) +++ trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2357) @@ -103,5 +103,5 @@ // 3 maximum output, showing progress of minimizations // - minuit.SetPrintLevel(1); + minuit.SetPrintLevel(-1); //.............................................. @@ -175,7 +175,8 @@ //.............................................. + // This doesn't seem to have any effect // Set maximum number of iterations (default = 500) //Int_t maxiter = 100000; - minuit.SetMaxIterations(10); + //minuit.SetMaxIterations(maxiter); //.............................................. @@ -211,9 +212,13 @@ if (nulloutput) fLog->SetNullOutput(kTRUE); - Double_t tmp = 0; - minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize); + Int_t maxcalls = 3000; + Double_t tolerance = 0.1; + Double_t tmp[2]; + tmp[0] = maxcalls; + tmp[1] = tolerance; + minuit.mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize); if (nulloutput) fLog->SetNullOutput(kFALSE); - //*fLog << "return from SIMPLEX" << endl; + *fLog << "return from SIMPLEX" << endl; } @@ -269,5 +274,5 @@ // 2 values, errors, step sizes, internal values // 3 values, errors, step sizes, 1st derivatives - // 4 values, paraboloc errors, MINOS errors + // 4 values, parabolic errors, MINOS errors //Int_t nkode = 4; @@ -292,2 +297,11 @@ return kTRUE; } + + + + + + + + + Index: trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc (revision 2356) +++ trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc (revision 2357) @@ -188,4 +188,7 @@ } + *fLog << "-------------------------" << endl; + *fLog << "MFEventSelector2::ReadDistribution; read input file to generate the original distribution" << endl; + MEvtLoop run(GetName()); MParList plist; @@ -215,4 +218,10 @@ } + tlist.PrintStatistics(0, kTRUE); + + *fLog << "MFEventSelector2::ReadDistribution; Original distribution has " + << fHistOrig->GetHist().GetEntries() << " entries" << endl; + *fLog << "-------------------------" << endl; + return read.Rewind(); } @@ -259,7 +268,9 @@ if (fNumMax>0) { - cout << "SCALE: " << fNumMax/hn.Integral() << endl; - cout << "SCALE: " << fNumMax << endl; - cout << "SCALE: " << hn.Integral() << endl; + *fLog << "MFEventSelector2::PrepareHistograms; SCALE : fNumMax = " + << fNumMax << ", hn.Integral() = " << hn.Integral() + << ", fNumMax/hn.Integral() = " << fNumMax/hn.Integral() + << endl; + hn.Scale(fNumMax/hn.Integral()); } @@ -329,4 +340,6 @@ Int_t MFEventSelector2::PreProcess(MParList *parlist) { + memset(fErrors, 0, sizeof(fErrors)); + MTaskList *tasklist = (MTaskList*)parlist->FindObject("MTaskList"); if (!tasklist) @@ -383,4 +396,5 @@ const Double_t valz=fDataZ.GetValue(); + // get corresponding bin number const Int_t bin = fHistNom->FindFixBin(valx, valy, valz)-1; @@ -389,4 +403,5 @@ if (bin<0) return kTRUE; + if (gRandom->Rndm()*fIs[bin]<=fNom[bin]) @@ -404,4 +419,16 @@ fIs[bin]-=1; + //---------------------- + Int_t rc; + if (fResult) + { + rc = 0; + fErrors[rc]++; + return kTRUE; + } + rc = 1; + fErrors[rc]++; + //---------------------- + return kTRUE; } @@ -413,4 +440,21 @@ Int_t MFEventSelector2::PostProcess() { + //--------------------------------- + if (GetNumExecutions() != 0) + { + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Events not selected" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Events selected!" << endl; + *fLog << endl; + } + //--------------------------------- + if (!fCanvas || !fDisplay) return kTRUE; Index: trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h (revision 2356) +++ trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h (revision 2357) @@ -42,4 +42,5 @@ Bool_t fResult; + Int_t fErrors[2]; TH1 &InitHistogram(MH3* &hist); Index: trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc (revision 2356) +++ trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc (revision 2357) @@ -61,5 +61,5 @@ // fName = name ? name : "MHCT1Supercuts"; - fTitle = title ? title : "Container for supercuts"; + fTitle = title ? title : "Container for histograms for the supercuts"; Index: trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2356) +++ trunk/MagicSoft/Mars/mhist/MHFindSignificance.cc (revision 2357) @@ -353,5 +353,6 @@ fHistOrig = fhist; - fHist = (TH1*)fHistOrig->Clone("|alpha| plot"); + fHist = (TH1*)fHistOrig->Clone(); + fHist->SetName(fhist->GetName()); if ( !fHist ) { @@ -362,7 +363,8 @@ fHist->Sumw2(); - fHist->SetNameTitle("Alpha", "alpha plot"); + //fHist->SetNameTitle("Alpha", "alpha plot"); fHist->SetXTitle("|alpha| [\\circ]"); fHist->SetYTitle("Counts"); + fHist->UseCurrentStyle(); fAlphamin = alphamin; @@ -502,10 +504,11 @@ fHistOrig = fhist; - fHist = (TH1*)fHistOrig->Clone("|alpha| plot"); + fHist = (TH1*)fHistOrig->Clone(); + fHist->SetName(fhist->GetName()); fHist->Sumw2(); - fHist->SetNameTitle("alpha", "alpha plot"); + //fHist->SetNameTitle("alpha", "alpha plot"); fHist->SetXTitle("|alpha| [\\circ]"); fHist->SetYTitle("Counts"); - + fHist->UseCurrentStyle(); fAlphamin = alphamin; @@ -754,7 +757,7 @@ nrebin += 1; - //if (fHist) - delete fHist; - fHist = NULL; + TString histname = fHist->GetName(); + delete fHist; + fHist = NULL; *fLog << "MHFindSignificance::FitPolynomial; rebin the |alpha| plot, grouping " @@ -766,6 +769,6 @@ fHist = new TH1F; fHist->Sumw2(); - fHist->SetNameTitle("Rebinned", "Rebinned alpha plot"); - + fHist->SetNameTitle(histname, histname); + fHist->UseCurrentStyle(); // do rebinning such that x0 remains a lower bin edge @@ -911,5 +914,5 @@ } - *fLog << "FitPolynomial : before CallMinuit()" << endl; + //*fLog << "FitPolynomial : before CallMinuit()" << endl; MMinuitInterface inter; @@ -918,5 +921,5 @@ kFALSE); - *fLog << "FitPolynomial : after CallMinuit()" << endl; + //*fLog << "FitPolynomial : after CallMinuit()" << endl; if (rc != 0) @@ -1909,4 +1912,5 @@ // get errors + /* Double_t eplus; Double_t eminus; @@ -1943,6 +1947,8 @@ } } + */ //---------------------------------------- + /* *fLog << "Covariance matrix :" << endl; for (Int_t j=0; j<=fDegree; j++) @@ -1966,4 +1972,5 @@ *fLog << endl; } + */ *fLog << "---------------------------" << endl; Index: trunk/MagicSoft/Mars/mhist/MHMatrix.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 2356) +++ trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 2357) @@ -30,5 +30,5 @@ static const TString gsDefTitle; //! Default Title - Int_t fNumRow; //! Number of dimensions of histogram + Int_t fNumRow; // Number of dimensions of histogram Int_t fRow; //! Present row TMatrix fM; // Matrix to be filled Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 2356) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 2357) @@ -219,5 +219,5 @@ // *fLog << inf << "--------------------------------------" << endl; - *fLog << inf << "Start minimization" << endl; + *fLog << inf << "Start minimization for the energy estimator" << endl; gMinuit = new TMinuit(12); @@ -316,5 +316,5 @@ // eest.Print(); eest.StopMapping(); - *fLog << inf << "Minimization finished" << endl; + *fLog << inf << "Minimization for the energy estimator finished" << endl; }