Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2542) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2543) @@ -286,4 +286,5 @@ + // Job E_XX : extended version of E_XX (including flux plots) // - select g/h separation method XX @@ -297,8 +298,10 @@ - Bool_t JobE_XX = kFALSE; + Bool_t JobE_XX = kTRUE; + Bool_t CCollArea= kFALSE; // calculate eff. collection areas Bool_t OEEst = kFALSE; // optimize energy estimator - Bool_t WEX = kFALSE; // update root file ? - Bool_t WRobert = kFALSE; // write out Robert's file ? + Bool_t WEX = kTRUE; // update root file ? + Bool_t WRobert = kTRUE; // write out Robert's file ? + @@ -2527,5 +2530,4 @@ - //--------------------------------------------------------------------- // Job E_XX @@ -2547,6 +2549,7 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " + gLog << "Macro CT1Analysis : JobE_XX, CCollArea, OEEst, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (CCollArea?"kTRUE" : "kFALSE") << ", " << (OEEst ? "kTRUE" : "kFALSE") << ", " << (WEX ? "kTRUE" : "kFALSE") << endl; @@ -2554,7 +2557,7 @@ // type of data to be analysed - TString typeData = "ON"; + //TString typeData = "ON"; //TString typeData = "OFF"; - //TString typeData = "MC"; + TString typeData = "MC"; gLog << "typeData = " << typeData << endl; @@ -2567,6 +2570,6 @@ // and definition of final selections - TString XX("SC"); - //TString XX("RF"); + //TString XX("SC"); + TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; @@ -2574,5 +2577,5 @@ // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.4; + Float_t maxhadronness = 0.23; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; @@ -2638,7 +2641,17 @@ //==================================================================== + + MHMcCT1CollectionArea collarea; + collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); + + MParList parlist; + InitBinnings(&parlist); + + if (CCollArea) + { gLog << "-----------------------------------------------" << endl; gLog << "Start calculation of effective collection areas" << endl; - MParList parlist; + + MTaskList tasklist; @@ -2655,6 +2668,4 @@ MContinue conthadrons(&cuthadrons); - MHMcCT1CollectionArea collarea; - collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); @@ -2665,5 +2676,5 @@ parlist.AddToList(&tasklist); - InitBinnings(&parlist); + parlist.AddToList(&collarea); @@ -2697,4 +2708,42 @@ collarea.CalcEfficiency(); collarea.DrawClone(); + + + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + //collarea.GetHist()->Write(); + //collarea.GetHAll()->Write(); + //collarea.GetHSel()->Write(); + collarea.Write(); + + f.Close(); + + gLog << "Collection area plots written onto file " << collareaName << endl; + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + } + + if (!CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Read in effective collection areas from file " + << collareaName << endl; + + TFile collfile(collareaName); + collfile.ls(); + collarea.Read("MHMcCT1CollectionArea"); + collarea.DrawClone(); + + gLog << "Effective collection areas were read in from file " + << collareaName << endl; + gLog << "-----------------------------------------------" << endl; + } + // save binnings for call to CT1EEst @@ -2712,22 +2761,5 @@ } - - //--------------------------------------------- - // Write histograms to a file - // - - TFile f(collareaName, "RECREATE"); - collarea.GetHist()->Write(); - collarea.GetHAll()->Write(); - collarea.GetHSel()->Write(); - f.Close(); - - gLog << "Collection area plots written onto file " << collareaName << endl; - - gLog << "Calculation of effective collection areas done" << endl; - gLog << "-----------------------------------------------" << endl; - //------------------------------------------------------------------ - - + //------------------------------------- TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt"; Index: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2542) +++ trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2543) @@ -207,6 +207,7 @@ // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file + // - if CTrainRF = TRUE : create matrices of training events + // and root files of training and test events // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas - // - if CTrainRF = TRUE : create matrices of training events // - if RTree = TRUE : read in trees, otherwise create trees // - calculate hadroness for method of RANDOM FOREST @@ -214,6 +215,7 @@ // or MC2.root) - Bool_t JobB_RF_UP = kFALSE; - Bool_t CTrainRF = kFALSE; // create matrices of training events + Bool_t JobB_RF_UP = kTRUE; + Bool_t CTrainRF = kTRUE; // create matrices of training events + // and root files of training and test events Bool_t RTrainRF = kFALSE; // read in matrices of training events Bool_t RTree = kFALSE; // read in trees (otherwise grow them) @@ -268,8 +270,9 @@ - Bool_t JobE_XX = kTRUE; + Bool_t JobE_XX = kFALSE; + Bool_t CCollArea= kFALSE; // calculate eff. collection areas Bool_t OEEst = kFALSE; // optimize energy estimator - Bool_t WEX = kTRUE; // update root file ? - Bool_t WRobert = kTRUE; // write out Robert's file ? + Bool_t WEX = kFALSE; // update root file ? + Bool_t WRobert = kFALSE; // write out Robert's file ? @@ -802,12 +805,20 @@ //-------------------------------------------- - // files to be read for generating the matrices of training events + // files to be read for generating + // - the matrices of training events + // - and the root files of training and test events + + TString typeMatrixHadrons = "OFF"; + gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + // "hadrons" : TString filenameHad = outPath; - filenameHad += "OFF"; + filenameHad += typeMatrixHadrons; filenameHad += "1.root"; - Int_t howManyHadrons = 1000000; - gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " - << howManyHadrons << endl; + Int_t howManyHadronsTrain = 8000; + Int_t howManyHadronsTest = 8000; + gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = " + << howManyHadronsTrain << ", howManyHadronsTest = " + << howManyHadronsTest << endl; @@ -816,23 +827,52 @@ filenameMC += "MC"; filenameMC += "1.root"; - Int_t howManyGammas = 50000; - gLog << "filenameMC = " << filenameMC << ", howManyGammas = " - << howManyGammas << endl; + Int_t howManyGammasTrain = 8000; + Int_t howManyGammasTest = 8000; + gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = " + << howManyGammasTrain << ", howManyGammasTest = " + << howManyGammasTest << endl; //-------------------------------------------- - // files of training events + // files for the matrices of training events TString outNameGammas = outPath; - outNameGammas += "RFmatrix_gammas_"; + outNameGammas += "RFmatrix_gammas_Train_"; outNameGammas += "MC"; outNameGammas += ".root"; - TString typeMatrixHadrons = "OFF"; - gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; - TString outNameHadrons = outPath; - outNameHadrons += "RFmatrix_hadrons_"; + outNameHadrons += "RFmatrix_hadrons_Train_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; + + + //-------------------------------------------- + // root files for the training events + + TString outNameGammasTrain = outPath; + outNameGammasTrain += "RF_gammas_Train_"; + outNameGammasTrain += "MC"; + outNameGammasTrain += ".root"; + + TString outNameHadronsTrain = outPath; + outNameHadronsTrain += "RF_hadrons_Train_"; + outNameHadronsTrain += typeMatrixHadrons; + outNameHadronsTrain += ".root"; + + + //-------------------------------------------- + // root files for the test events + + TString outNameGammasTest = outPath; + outNameGammasTest += "RF_gammas_Test_"; + outNameGammasTest += "MC"; + outNameGammasTest += ".root"; + + TString outNameHadronsTest = outPath; + outNameHadronsTest += "RF_hadrons_Test_"; + outNameHadronsTest += typeMatrixHadrons; + outNameHadronsTest += ".root"; + + //-------------------------------------------------------------------- @@ -908,30 +948,29 @@ //************************************************************************* // create matrices of training events + // and root files of training and test events + if (CTrainRF) { gLog << "" << endl; gLog << "========================================================" << endl; - gLog << " Create matrices of training events" << endl; + gLog << " Create matrices of training events and root files of training and test events" + << endl; gLog << " Gammas :" << endl; - + gLog << "---------" << endl; MParList plistg; MTaskList tlistg; - MFilterList flistg; - - MParList plisth; - MTaskList tlisth; - MFilterList flisth; + + MFilterList flistgtrain; + flistgtrain.SetName("GammaFilterTrain"); + MFilterList flistgtest; + flistgtest.SetName("GammaFilterTest"); MReadMarsFile readg("Events", filenameMC); readg.DisableAutoScheme(); - MReadMarsFile readh("Events", filenameHad); - readh.DisableAutoScheme(); MFParticleId fgamma("MMcEvt", '=', kGAMMA); fgamma.SetName("gammaID"); - MFParticleId fhadrons("MMcEvt", '!', kGAMMA); - fhadrons.SetName("hadronID)"); TString mgname("costhg"); @@ -946,17 +985,67 @@ MFEventSelector2 selectorg(gref); - selectorg.SetNumMax(howManyGammas); - selectorg.SetName("selectGammas"); + selectorg.SetNumMax(howManyGammasTrain+howManyGammasTest); + selectorg.SetName("selectGammasTrainTest"); + + TString strgtrain = "gRandom->Rndm()*("; + strgtrain += howManyGammasTrain; + strgtrain += "+"; + strgtrain += howManyGammasTest; + strgtrain += ")<"; + strgtrain += howManyGammasTrain; + MF selsplitg(strgtrain); MFillH fillmatg("MatrixGammas"); - fillmatg.SetFilter(&flistg); + fillmatg.SetFilter(&flistgtrain); fillmatg.SetName("fillGammas"); + //----------------------- + // for writing the root files of training and test events + // for gammas + + MWriteRootFile writetraing(outNameGammasTrain, "RECREATE"); + writetraing.SetName("WriteGammasTrain"); + writetraing.SetFilter(&flistgtrain); + + writetraing.AddContainer("MRawRunHeader", "RunHeaders"); + writetraing.AddContainer("MTime", "Events"); + writetraing.AddContainer("MMcEvt", "Events"); + writetraing.AddContainer("ThetaOrig", "Events"); + writetraing.AddContainer("MSrcPosCam", "Events"); + writetraing.AddContainer("MSigmabar", "Events"); + writetraing.AddContainer("MHillas", "Events"); + writetraing.AddContainer("MHillasExt", "Events"); + writetraing.AddContainer("MHillasSrc", "Events"); + writetraing.AddContainer("MNewImagePar", "Events"); + + MContinue contg(&flistgtrain); + + MWriteRootFile writetestg(outNameGammasTest, "RECREATE"); + writetestg.SetName("WriteGammasTest"); + writetestg.SetFilter(&flistgtest); + + writetestg.AddContainer("MRawRunHeader", "RunHeaders"); + writetestg.AddContainer("MTime", "Events"); + writetestg.AddContainer("MMcEvt", "Events"); + writetestg.AddContainer("ThetaOrig", "Events"); + writetestg.AddContainer("MSrcPosCam", "Events"); + writetestg.AddContainer("MSigmabar", "Events"); + writetestg.AddContainer("MHillas", "Events"); + writetestg.AddContainer("MHillasExt", "Events"); + writetestg.AddContainer("MHillasSrc", "Events"); + writetestg.AddContainer("MNewImagePar", "Events"); + + //----------------------- //***************************** fill gammas *** // entries in MFilterList - flistg.AddToList(&fgamma); - flistg.AddToList(&selectorg); + flistgtrain.AddToList(&fgamma); + flistgtrain.AddToList(&selectorg); + flistgtrain.AddToList(&selsplitg); + + flistgtest.AddToList(&fgamma); + flistgtest.AddToList(&selectorg); + flistgtest.AddToList(&selsplitg); //***************************** @@ -972,6 +1061,13 @@ tlistg.AddToList(&readg); - tlistg.AddToList(&flistg); + + tlistg.AddToList(&flistgtrain); tlistg.AddToList(&fillmatg); + tlistg.AddToList(&writetraing); + + tlistg.AddToList(&contg); + + tlistg.AddToList(&flistgtest); + tlistg.AddToList(&writetestg); //***************************** @@ -990,8 +1086,42 @@ tlistg.PrintStatistics(0, kTRUE); + matrixg.Print("SizeCols"); + Int_t generatedgTrain = matrixg.GetM().GetNrows(); + if ( fabs(generatedgTrain-howManyGammasTrain) > + 3.0*sqrt(howManyGammasTrain) ) + { + gLog << "ONOFFCT1Analysis.C : no.of generated gamma training events (" + << generatedgTrain << ") is incompatible with the no.of requested events (" + << howManyGammasTrain << ")" << endl; + } + + + Int_t generatedgTest = writetestg.GetNumExecutions(); + if ( fabs(generatedgTest-howManyGammasTest) > + 3.0*sqrt(howManyGammasTest) ) + { + gLog << "ONOFFCT1Analysis.C : no.of generated gamma test events (" + << generatedgTest << ") is incompatible with the no.of requested events (" + << howManyGammasTest << ")" << endl; + } //***************************** fill hadrons *** gLog << " Hadrons :" << endl; + gLog << "----------" << endl; + + MParList plisth; + MTaskList tlisth; + + MFilterList flisthtrain; + flisthtrain.SetName("HadronFilterTrain"); + MFilterList flisthtest; + flisthtest.SetName("HadronFilterTest"); + + MReadMarsFile readh("Events", filenameHad); + readh.DisableAutoScheme(); + + MFParticleId fhadrons("MMcEvt", '!', kGAMMA); + fhadrons.SetName("hadronID)"); TString mhname("costhh"); @@ -1006,18 +1136,66 @@ MFEventSelector2 selectorh(href); - //selectorh.SetNumMax(howManyHadrons); - // select as many hadrons as gammas - selectorh.SetNumMax(matrixg.GetM().GetNrows()); - selectorh.SetName("selectHadrons"); + selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest); + selectorh.SetName("selectHadronsTrainTest"); + + TString strhtrain = "gRandom->Rndm()*("; + strhtrain += howManyHadronsTrain; + strhtrain += "+"; + strhtrain += howManyHadronsTest; + strhtrain += ")<"; + strhtrain += howManyHadronsTrain; + MF selsplith(strhtrain); MFillH fillmath("MatrixHadrons"); - fillmath.SetFilter(&flisth); + fillmath.SetFilter(&flisthtrain); fillmath.SetName("fillHadrons"); + //----------------------- + // for writing the root files of training and test events + // for hadrons + + MWriteRootFile writetrainh(outNameHadronsTrain, "RECREATE"); + writetrainh.SetName("WriteHadronsTrain"); + writetrainh.SetFilter(&flisthtrain); + + writetrainh.AddContainer("MRawRunHeader", "RunHeaders"); + writetrainh.AddContainer("MTime", "Events"); + writetrainh.AddContainer("MMcEvt", "Events"); + writetrainh.AddContainer("ThetaOrig", "Events"); + writetrainh.AddContainer("MSrcPosCam", "Events"); + writetrainh.AddContainer("MSigmabar", "Events"); + writetrainh.AddContainer("MHillas", "Events"); + writetrainh.AddContainer("MHillasExt", "Events"); + writetrainh.AddContainer("MHillasSrc", "Events"); + writetrainh.AddContainer("MNewImagePar", "Events"); + + MContinue conth(&flisthtrain); + + MWriteRootFile writetesth(outNameHadronsTest, "RECREATE"); + writetesth.SetName("WriteHadronsTest"); + writetesth.SetFilter(&flisthtest); + + writetesth.AddContainer("MRawRunHeader", "RunHeaders"); + writetesth.AddContainer("MTime", "Events"); + writetesth.AddContainer("MMcEvt", "Events"); + writetesth.AddContainer("ThetaOrig", "Events"); + writetesth.AddContainer("MSrcPosCam", "Events"); + writetesth.AddContainer("MSigmabar", "Events"); + writetesth.AddContainer("MHillas", "Events"); + writetesth.AddContainer("MHillasExt", "Events"); + writetesth.AddContainer("MHillasSrc", "Events"); + writetesth.AddContainer("MNewImagePar", "Events"); + + //----------------------- // entries in MFilterList - flisth.AddToList(&fhadrons); - flisth.AddToList(&selectorh); + flisthtrain.AddToList(&fhadrons); + flisthtrain.AddToList(&selectorh); + flisthtrain.AddToList(&selsplith); + + flisthtest.AddToList(&fhadrons); + flisthtest.AddToList(&selectorh); + flisthtest.AddToList(&selsplith); //***************************** @@ -1033,6 +1211,13 @@ tlisth.AddToList(&readh); - tlisth.AddToList(&flisth); + + tlisth.AddToList(&flisthtrain); tlisth.AddToList(&fillmath); + tlisth.AddToList(&writetrainh); + + tlisth.AddToList(&conth); + + tlisth.AddToList(&flisthtest); + tlisth.AddToList(&writetesth); //***************************** @@ -1050,4 +1235,26 @@ tlisth.PrintStatistics(0, kTRUE); + + matrixh.Print("SizeCols"); + Int_t generatedhTrain = matrixh.GetM().GetNrows(); + if ( fabs(generatedhTrain-howManyHadronsTrain) > + 3.0*sqrt(howManyHadronsTrain) ) + { + gLog << "ONOFFCT1Analysis.C : no.of generated hadron training events (" + << generatedhTrain << ") is incompatible with the no.of requested events (" + << howManyHadronsTrain << ")" << endl; + } + + + Int_t generatedhTest = writetesth.GetNumExecutions(); + if ( fabs(generatedhTest-howManyHadronsTest) > + 3.0*sqrt(howManyHadronsTest) ) + { + gLog << "ONOFFCT1Analysis.C : no.of generated gamma test events (" + << generatedhTest << ") is incompatible with the no.of requested events (" + << howManyHadronsTest << ")" << endl; + } + + //***************************************************** @@ -1204,6 +1411,4 @@ // end of growing the trees of the random forest //----------------------------------------------------------------- - - @@ -2419,6 +2624,7 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " + gLog << "Macro CT1Analysis : JobE_XX, CCollArea, OEEst, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (CCollArea?"kTRUE" : "kFALSE") << ", " << (OEEst ? "kTRUE" : "kFALSE") << ", " << (WEX ? "kTRUE" : "kFALSE") << endl; @@ -2510,7 +2716,17 @@ //==================================================================== + + MHMcCT1CollectionArea collarea; + collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); + + MParList parlist; + InitBinnings(&parlist); + + if (CCollArea) + { gLog << "-----------------------------------------------" << endl; gLog << "Start calculation of effective collection areas" << endl; - MParList parlist; + + MTaskList tasklist; @@ -2527,6 +2743,4 @@ MContinue conthadrons(&cuthadrons); - MHMcCT1CollectionArea collarea; - collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); @@ -2537,5 +2751,5 @@ parlist.AddToList(&tasklist); - InitBinnings(&parlist); + parlist.AddToList(&collarea); @@ -2569,4 +2783,42 @@ collarea.CalcEfficiency(); collarea.DrawClone(); + + + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + //collarea.GetHist()->Write(); + //collarea.GetHAll()->Write(); + //collarea.GetHSel()->Write(); + collarea.Write(); + + f.Close(); + + gLog << "Collection area plots written onto file " << collareaName << endl; + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + } + + if (!CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Read in effective collection areas from file " + << collareaName << endl; + + TFile collfile(collareaName); + collfile.ls(); + collarea.Read("MHMcCT1CollectionArea"); + collarea.DrawClone(); + + gLog << "Effective collection areas were read in from file " + << collareaName << endl; + gLog << "-----------------------------------------------" << endl; + } + // save binnings for call to CT1EEst @@ -2584,22 +2836,5 @@ } - - //--------------------------------------------- - // Write histograms to a file - // - - TFile f(collareaName, "RECREATE"); - collarea.GetHist()->Write(); - collarea.GetHAll()->Write(); - collarea.GetHSel()->Write(); - f.Close(); - - gLog << "Collection area plots written onto file " << collareaName << endl; - - gLog << "Calculation of effective collection areas done" << endl; - gLog << "-----------------------------------------------" << endl; - //------------------------------------------------------------------ - - + //------------------------------------- TString fHilName = "MHillas"; TString fHilNameExt = "MHillasExt";