Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2363) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2364) @@ -252,9 +252,11 @@ // - update the input files with the hadroness (ON1.root or MC1.root) - Bool_t JobB_SC_UP = kFALSE; + Bool_t JobB_SC_UP = kTRUE; 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 ? + Bool_t WOptimize = kFALSE; // do optimization using the training sample + // and write supercuts parameter values + // onto the file parSCfile + Bool_t RTest = kFALSE; // test the supercuts using the test matrix + Bool_t WSC = kTRUE; // update input root file ? @@ -281,18 +283,6 @@ - // Job E_XX : - // - select g/h separation method XX - // - read MC root file -   // - calculate eff. collection area - // - optimize energy estimator - // - read ON root file - // - apply final cuts - // - write root file for ON data after final cuts - - Bool_t JobE_XX = kFALSE; - Bool_t WXX = kFALSE; // write out root file ? - - - // Job F_XX : extended version of E_XX (including flux plots) + + // Job E_XX : extended version of E_XX (including flux plots) // - select g/h separation method XX // - read MC root file @@ -304,9 +294,8 @@ // - write root file for ON data after final cuts - 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 ? + + Bool_t JobE_XX = kFALSE; + Bool_t WEX = kFALSE; // write out root file ? + Bool_t WRobert = kFALSE; // write out Robert's file ? @@ -1550,9 +1539,10 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = " + gLog << "Macro CT1Analysis : JobB_SC_UP, RMatrix, WOptimize, RTest, WSC = " << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " - << (WParSC ? "kTRUE" : "kFALSE") << ", " + << (RMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "kTRUE" : "kFALSE") << ", " << (WSC ? "kTRUE" : "kFALSE") << endl; - @@ -1580,6 +1570,6 @@ // file to be updated (either ON or MC) - //TString typeInput = "ON"; - TString typeInput = "MC"; + TString typeInput = "ON"; + //TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -1652,5 +1642,5 @@ //-------------------------- // create matrices and write them onto files - if (!RMatrix) + if (!RMatrix && (WOptimize || RTest) ) { MH3 &mh3 = *(new MH3("MHillas.fSize")); @@ -1705,6 +1695,9 @@ -if (WParSC) +if (WOptimize) { + gLog << "CT1Analysis.C : optimize the supercuts using the training matrix" + << endl; + TArrayD params(0); TArrayD steps(0); @@ -1806,10 +1799,21 @@ return; } + } //-------------------------------------- // test the supercuts on the test sample // + + if (RTest) + { + gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl; Bool_t rt = findsuper.TestParams(); - } + if (!rt) + { + gLog << "CT1Analysis.C : test of supercuts on the test matrix failed" + << endl; + } + + } @@ -1906,4 +1910,6 @@ //MWriteRootFile write(outNameImage, "UPDATE"); //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE"); + + /* MWriteRootFile write(outNameImage, "RECREATE"); @@ -1921,5 +1927,5 @@ write.AddContainer("HadRF", "Events"); write.AddContainer(hadSCName, "Events"); - + */ //----------------------------------------------------------------- @@ -1995,5 +2001,5 @@ tliston.AddToList(&fillhadsc); - tliston.AddToList(&write); + //tliston.AddToList(&write); tliston.AddToList(&contfinalgh); @@ -2462,4 +2468,5 @@ + //--------------------------------------------------------------------- // Job E_XX @@ -2481,7 +2488,7 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_XX, WXX = " + gLog << "Macro CT1Analysis : JobE_XX, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " - << (WXX ? "kTRUE" : "kFALSE") << endl; + << (WEX ? "kTRUE" : "kFALSE") << endl; @@ -2561,4 +2568,13 @@ gLog << "filenameDataout = " << filenameDataout << endl; + //------------------------------ + // name of file containing histograms for flux calculastion + TString filenameResults(outPath); + filenameResults += typeData; + filenameResults += "Results_"; + filenameResults += XX; + filenameResults += extout; + gLog << "filenameResults = " << filenameResults << endl; + //==================================================================== @@ -2580,6 +2596,6 @@ MContinue conthadrons(&cuthadrons); - //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); - //MHMcCT1CollectionArea* collarea; + MHMcCT1CollectionArea collarea; + collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); @@ -2591,5 +2607,5 @@ parlist.AddToList(&tasklist); InitBinnings(&parlist); - //parlist.AddToList(collarea); + parlist.AddToList(&collarea); //******************************** @@ -2618,8 +2634,8 @@ // and display the histograms // - MHMcCT1CollectionArea *collarea = - (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); - collarea->CalcEfficiency(); - collarea->DrawClone("lego"); + //MHMcCT1CollectionArea *collarea = + // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); + collarea.CalcEfficiency(); + collarea.DrawClone(); // save binnings for call to CT1EEst @@ -2643,7 +2659,7 @@ TFile f(collareaName, "RECREATE"); - collarea->GetHist()->Write(); - collarea->GetHAll()->Write(); - collarea->GetHSel()->Write(); + collarea.GetHist()->Write(); + collarea.GetHAll()->Write(); + collarea.GetHSel()->Write(); f.Close(); @@ -2660,5 +2676,5 @@ TString fImgParName = "MNewImagePar"; - + //=========================================================== // @@ -2677,5 +2693,5 @@ - if (WXX) + if (WEX) { //----------------------------------------------------------- @@ -2697,6 +2713,6 @@ gLog << "Parameters of energy estimator were read in" << endl; - TArrayD parA(5); - TArrayD parB(7); + TArrayD parA(mcest.GetNumCoeffA()); + TArrayD parB(mcest.GetNumCoeffB()); for (Int_t i=0; iDrawClone(); - pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); - pliston.FindObject("hadSC", "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 E_XX" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - - - //--------------------------------------------------------------------- - // Job F_XX - //========= - - // - select g/h separation method XX - // - read MC_XX2.root file -   // - calculate eff. collection area - // - read ON_XX2.root file - // - apply final cuts - // - calculate flux - // - write root file for ON data after final cuts (ON_XX3.root)) - - - if (JobF_XX) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobF_XX, WFX = " - << (JobF_XX ? "kTRUE" : "kFALSE") << ", " - << (WFX ? "kTRUE" : "kFALSE") << endl; - - - // type of data to be analysed - TString typeData = "ON"; - //TString typeData = "OFF"; - //TString typeData = "MC"; - gLog << "typeData = " << typeData << endl; - - TString typeMC = "MC"; - TString ext = "3.root"; - TString extout = "4.root"; - - //------------------------------ - // selection of g/h separation method - // and definition of final selections - - //TString XX("NN"); - TString XX("SC"); - //TString XX("RF"); - TString fhadronnessName("Had"); - fhadronnessName += XX; - gLog << "fhadronnessName = " << fhadronnessName << endl; - - // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - gLog << "Maximum values of hadronness, |ALPHA| and DIST = " - << maxhadronness << ", " << maxalpha << ", " - << maxdist << endl; - - //------------------------------ - // name of MC file to be used for optimizing the energy estimator - TString filenameOpt(outPath); - filenameOpt += typeMC; - filenameOpt += ext; - gLog << "filenameOpt = " << filenameOpt << endl; - - //------------------------------ - // name of file containing the parameters of the energy estimator - TString energyParName(outPath); - energyParName += "energyest_"; - energyParName += XX; - energyParName += ".root"; - gLog << "energyParName = " << energyParName << endl; - - //------------------------------ - // name of MC file to be used for calculating the eff. collection areas - TString filenameArea(outPath); - filenameArea += typeMC; - filenameArea += ext; - gLog << "filenameArea = " << filenameArea << endl; - - //------------------------------ - // name of file containing the eff. collection areas - TString collareaName(outPath); - collareaName += "area_"; - collareaName += XX; - collareaName += ".root"; - gLog << "collareaName = " << collareaName << endl; - - //------------------------------ - // name of data file to be analysed - TString filenameData(outPath); - filenameData += typeData; - filenameData += ext; - gLog << "filenameData = " << filenameData << endl; - - //------------------------------ - // name of output data file (after the final cuts) - TString filenameDataout(outPath); - filenameDataout += typeData; - filenameDataout += "_"; - filenameDataout += XX; - filenameDataout += extout; - gLog << "filenameDataout = " << filenameDataout << endl; - - //------------------------------ - // name of file containing histograms for flux calculastion - TString filenameResults(outPath); - filenameResults += typeData; - filenameResults += "Results_"; - filenameResults += XX; - filenameResults += extout; - gLog << "filenameResults = " << filenameResults << endl; - - - //==================================================================== - gLog << "-----------------------------------------------" << endl; - gLog << "Start calculation of effective collection areas" << endl; - MParList parlist; - MTaskList tasklist; - - //--------------------------------------- - // Setup the tasks to be executed - // - MReadMarsFile reader("Events", filenameArea); - reader.DisableAutoScheme(); - - MFCT1SelFinal cuthadrons; - cuthadrons.SetHadronnessName(fhadronnessName); - cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); - - MContinue conthadrons(&cuthadrons); - - MHMcCT1CollectionArea collarea; - collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); - - MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); - filler.SetName("CollectionArea"); - - //******************************** - // entries in MParList - - parlist.AddToList(&tasklist); - InitBinnings(&parlist); - parlist.AddToList(&collarea); - - //******************************** - // entries in MTaskList - - tasklist.AddToList(&reader); - tasklist.AddToList(&conthadrons); - tasklist.AddToList(&filler); - - //******************************** - - //----------------------------------------- - // Execute event loop - // - MEvtLoop magic; - magic.SetParList(&parlist); - - MProgressBar bar; - magic.SetProgressBar(&bar); - if (!magic.Eventloop()) - return; - - tasklist.PrintStatistics(0, kTRUE); - - // Calculate effective collection areas - // and display the histograms - // - //MHMcCT1CollectionArea *collarea = - // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); - collarea.CalcEfficiency(); - collarea.DrawClone(); - - // save binnings for call to CT1EEst - MBinning *binsE = (MBinning*)parlist.FindObject("BinningE"); - if (!binsE) - { - gLog << "Object 'BinningE' not found in MParList" << endl; - return; - } - MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta"); - if (!binsTheta) - { - gLog << "Object 'BinningTheta' not found in MParList" << endl; - return; - } - - - //--------------------------------------------- - // 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"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - - //=========================================================== - // - // Optimization of energy estimator - // - gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl; - - TString inpath(""); - TString outpath(""); - Int_t howMany = 2000; - CT1EEst(inpath, filenameOpt, outpath, energyParName, - fHilName, fHilNameSrc, fhadronnessName, - howMany, maxhadronness, maxalpha, maxdist, - binsE, binsTheta); - gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; - - - if (WFX) - { - //----------------------------------------------------------- - // - // Read in parameters of energy estimator ("MMcEnergyEst") - // and migration matrix ("MHMcEnergyMigration") - // - gLog << "================================================================" - << endl; - gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '" - << energyParName << "'" << endl; - TFile enparam(energyParName); - enparam.ls(); - - //MMcEnergyEst mcest("MMcEnergyEst"); - //mcest.Read("MMcEnergyEst"); - MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst")); - - gLog << "Parameters of energy estimator were read in" << endl; - - TArrayD parA(mcest.GetNumCoeffA()); - TArrayD parB(mcest.GetNumCoeffB()); - for (Int_t i=0; iFindObject("MHMcEnergyMigration")); - - gLog << "Migration matrix was read in" << endl; - - //************************************************************************* - // - // Analyse the data - // - gLog << "============================================================" - << endl; - gLog << "Analyse the data" << 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(); - - //....................................................................... - - gLog << "CT1Analysis.C : write root file '" << filenameDataout - << "'" << endl; - - //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); - /* - MWriteRootFile write(filenameDataout, "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("HadRF", "Events"); - - write.AddContainer("MEnergyEst", "Events"); - */ - - //----------------------------------------------------------------- - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - MFCT1SelFinal selfinalgh(fHilNameSrc); - selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(fhadronnessName); - selfinalgh.SetName("SelFinalgh"); - MContinue contfinalgh(&selfinalgh); - contfinalgh.SetName("ContSelFinalgh"); - - //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); - //fillhadnn.SetName("HhadNN"); - MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); - fillhadsc.SetName("HhadSC"); - MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); - fillhadrf.SetName("HhadRF"); - - //--------------------------- - // calculate estimated energy - - 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); - //eeston.SetCoeffA(parA); - //eeston.SetCoeffB(parB); - - - //--------------------------- - - - 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"); - - //--------------------------- - // new from Robert - - MFillH hfill6("MHTimeDiffTheta", "MMcEvt"); - hfill6.SetName("HTimeDiffTheta"); - - MFillH hfill6a("MHTimeDiffTime", "MMcEvt"); - hfill6a.SetName("HTimeDiffTime"); - - MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc); - hfill7.SetName("HAlphaEnergyTheta"); - - MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc); - hfill7a.SetName("HAlphaEnergyTime"); - - MFillH hfill7b("MHThetabarTime", fHilNameSrc); - hfill7b.SetName("HThetabarTime"); - - MFillH hfill7c("MHEnergyTime", "MMcEvt"); - hfill7c.SetName("HEnergyTime"); - - - //--------------------------- - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(fhadronnessName); - selfinal.SetName("SelFinal"); - MContinue contfinal(&selfinal); - contfinal.SetName("ContSelFinal"); - - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - // robert tliston.AddToList(&hfill6); //timediff @@ -3309,5 +2885,5 @@ tliston.AddToList(&eeston); - //tliston.AddToList(&write); + tliston.AddToList(&write); //tliston.AddToList(&fillhadnn);