Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2300) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2301) @@ -1,4 +1,26 @@ #include "CT1EgyEst.C" + +#include "MBinning.h" +#include "MBlindPixelCalc.h" +#include "MContinue.h" +#include "MCT1PointingCorrCalc.h" + +#include "MFCT1SelBasic.h" +#include "MFCT1SelStandard.h" +#include "MFCT1SelFinal.h" +#include "MFillH.h" + +#include "MHillasCalc.h" +#include "MHillasSrcCalc.h" +#include "MImgCleanStd.h" + +#include "MParList.h" +#include "MSigmabarCalc.h" +#include "MSrcPosCam.h" +#include "MTaskList.h" +#include "MWriteRootFile.h" + + void InitBinnings(MParList *plist) @@ -6,8 +28,13 @@ gLog << "InitBinnings" << endl; + //-------------------------------------------- MBinning *binse = new MBinning("BinningE"); //binse->SetEdgesLog(30, 1.0e2, 1.0e5); - binse->SetEdges(30, 2, 5); + + //This is Daniel's binning in energy: + binse->SetEdgesLog(14, 296.296, 86497.6); plist->AddToList(binse); + + //-------------------------------------------- MBinning *binssize = new MBinning("BinningSize"); @@ -51,5 +78,5 @@ MBinning *binth = new MBinning("BinningTheta"); Double_t yedge[9] = - {7.5, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; + {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; TArrayD yed(9,yedge); binth->SetEdges(yed); @@ -69,4 +96,17 @@ binsdiff->SetEdges(100, -5.0, 20.0); plist->AddToList(binsdiff); + + // robert ---------------------------------------------- + MBinning *binsalphaf = new MBinning("BinningAlphaFlux"); + binsalphaf->SetEdges(100, -100, 100); + plist->AddToList(binsalphaf); + + MBinning *binsdifftime = new MBinning("BinningTimeDiff"); + binsdifftime->SetEdges(50, 0, 10); + plist->AddToList(binsdifftime); + + MBinning *binstime = new MBinning("BinningTime"); + binstime->SetEdges(50, 44500, 61000); + plist->AddToList(binstime); } @@ -112,4 +152,14 @@ bin = plist->FindObject("BinningDiffsigma2"); + if (bin) delete bin; + + //robert + bin = plist->FindObject("BinningAlphaFlux"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTimeDiff"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTime"); if (bin) delete bin; } @@ -173,4 +223,16 @@ 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 ? @@ -209,16 +271,29 @@ - // Job E_XX : // - select g/h separation method XX - // - read MC_XX2.root file + // - read MC root file   // - calculate eff. collection area - // - read ON_XX2.root file + // - 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) + // - select g/h separation method XX + // - read MC root file +   // - calculate eff. collection area + // - optimize energy estimator + // - read ON root file // - apply final cuts // - calculate flux - // - write root file for ON data after final cuts (ON3.root)) - - Bool_t JobE_XX = kTRUE; - Bool_t WXX = kFALSE; // write out root file ON3.root ? + // - write root file for ON data after final cuts + + Bool_t JobF_XX = kFALSE; + Bool_t WFX = kFALSE; // write out root file ? @@ -344,6 +419,7 @@ contstandard.SetName("SelStandard"); - - MWriteRootFile write(outNameImage); + if (WON1) + { + MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -357,5 +433,5 @@ write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); - + } //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ @@ -677,6 +753,7 @@ - - MWriteRootFile write(outNameImage); + if (WMC1) + { + MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -690,5 +767,5 @@ write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MNewImagePar", "Events"); - + } @@ -1315,4 +1392,407 @@ //--------------------------------------------------------------------- + // Job B_SC_UP + //============ + + // - create (or read in) optimum supercuts parameter values + // + // - calculate the hadroness for the supercuts + // + // - update input root file, including the hadroness + + + if (JobB_SC_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = " + << JobB_SC_UP << ", " << WParSC << ", " << WSC << endl; + + + + //-------------------------------------------- + // 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"; + outNameImage += "xxxxx"; + + + //TString outNameImage = filenameData; + + gLog << "outNameImage = " << outNameImage << endl; + + //-------------------------------------------- + // files to be read for optimizing the supercuts + // + // for the training + TString filenameTrain = outPath; + filenameTrain += "ON"; + filenameTrain += "2.root"; + Int_t howManyTrain = 800000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = outPath; + filenameTest += "ON"; + filenameTest += "2.root"; + Int_t howManyTest = 100000; + + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // files to contain the matrices (genertated from filenameTrain and + // filenameTest) + // + // for the training + TString fileMatrixTrain = outPath; + fileMatrixTrain += "MatrixTrainSC"; + fileMatrixTrain += ".root"; + gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; + + // for testing + TString fileMatrixTest = outPath; + fileMatrixTest += "MatrixTestSC"; + fileMatrixTest += ".root"; + gLog << "fileMatrixTest = " << fileMatrixTest << endl; + + + //-------------------------------------------- + // 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) + { + MCT1FindSupercuts findsuper; + findsuper.SetFilenameParam(parSCfile); + findsuper.SetHadronnessName("HadSC"); + + + //-------------------------- + // create matrices + if (!RMatrix) + { + MH3 &mh3 = *(new MH3("MHillas.fSize")); + mh3.SetName("Target distribution for SIZE"); + + if (filenameTrain == filenameTest) + { + findsuper.DefineTrainTestMatrix(filenameTrain, + howManyTrain, mh3, + howManyTest, mh3); + } + else + { + findsuper.DefineTrainMatrix(filenameTrain, howManyTrain, mh3); + findsuper.DefineTestMatrix( filenameTest, howManyTest, mh3); + } + + // writing doesn't work ??? + //findsuper.WriteMatrix(fileMatrixTrain, fileMatrixTest); + } + + //-------------------------- + // read matrices from files + // NOT YET TESTED + //if (RMatrix) + // findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); + //-------------------------- + + + //-------------------------- + // optimize the supercuts + Bool_t rf = findsuper.FindParams(); + + if (!rf) + { + *fLog << "CT1Analysis.C : optimization of supercuts failed" << endl; + return; + } + + //-------------------------------------- + // test the supercuts on the test sample + // + Bool_t rt = findsuper.TestParams(); + } + + + //----------------------------------------------------------------- + // Update the input files with the SC hadronness + // + + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update input file '" << filenameData + << "' with the SC hadronness" << endl; + + + //---------------------------------------------------- + // read in optimum parameter values for the supercuts + + TFile inparam(parSCfile); + MCT1Supercuts scin; + scin.Read("MCT1Supercuts"); + inparam.Close(); + + gLog << "Optimum parameter values for supercuts were read in from file '" + << parSCfile << "'" << endl; + + TArrayD supercutsPar(72); + scin.GetParams(72, supercutsPar); + + gLog << "Optimum parameter values for supercuts : " << endl; + for (Int_t i=0; i<72; i++) + { + gLog << supercutsPar[i] << ", "; + } + gLog << endl; + + + //---------------------------------------------------- + MTaskList tliston; + MParList pliston; + + // set the parameters of the supercuts + MCT1Supercuts supercuts; + supercuts.SetParams(72, supercutsPar); + gLog << "parameter values for the supercuts used for updating the input file ' " + << filenameData << "'" << endl; + supercuts.GetParams(72, supercutsPar); + for (Int_t i=0; i<72; i++) + { + gLog << supercutsPar[i] << ", "; + } + gLog << endl; + + + // 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"; + + + //....................................................................... + // calculation of hadroness for the supercuts + // (=0.25 if fullfilled, =0.75 otherwise) + + TString hadSCName = "HadSC"; + MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); + sccalc.SetHadronnessName(hadSCName); + + + //....................................................................... + + + if (WSC) + { + //MWriteRootFile write(outNameImage, "UPDATE"); + MWriteRootFile write = new MWriteRootFile(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"); + } + + + //----------------------------------------------------------------- + // 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(hadSCName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadSCName); + 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); + 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(&binsalphaabs); + pliston.AddToList(&alphaabs); + pliston.AddToList(&supercuts); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + if (WSC) + tliston.AddToList(&write); + + 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("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + 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_SC_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- // Job B_RF_UP //============ @@ -2336,5 +2816,5 @@ MHOnSubtraction onsub; - onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); + onsub.Calc(alphaHist, kTRUE, 13.1); //------------------------------------------- @@ -2387,6 +2867,6 @@ //TString XX("NN"); - //TString XX("SC"); - TString XX("RF"); + TString XX("SC"); + //TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; @@ -2412,5 +2892,5 @@ TString energyParName(outPath); energyParName += "energyest_"; - //energyParName += XX; + energyParName += XX; energyParName += ".root"; gLog << "energyParName = " << energyParName << endl; @@ -2420,6 +2900,4 @@ TString filenameArea(outPath); filenameArea += typeMC; - //filenameArea += "_"; - //filenameArea += XX; filenameArea += ext; gLog << "filenameArea = " << filenameArea << endl; @@ -2429,5 +2907,5 @@ TString collareaName(outPath); collareaName += "area_"; - //collareaName += XX; + collareaName += XX; collareaName += ".root"; gLog << "collareaName = " << collareaName << endl; @@ -2437,6 +2915,4 @@ TString filenameData(outPath); filenameData += typeData; - //filenameData += "_"; - //filenameData += XX; filenameData += ext; gLog << "filenameData = " << filenameData << endl; @@ -2446,6 +2922,6 @@ TString filenameDataout(outPath); filenameDataout += typeData; - //filenameDataout += "_"; - //filenameDataout += XX; + filenameDataout += "_"; + filenameDataout += XX; filenameDataout += extout; gLog << "filenameDataout = " << filenameDataout << endl; @@ -2538,4 +3014,5 @@ f.Close(); + gLog << "Collection area plots written onto file " << collareaName << endl; gLog << "Calculation of effective collection areas done" << endl; @@ -2554,4 +3031,5 @@ // Optimization of energy estimator // + gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl; TString inpath(""); @@ -2562,4 +3040,5 @@ howMany, maxhadronness, maxalpha, maxdist, binsE, binsTheta); + gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; //----------------------------------------------------------- @@ -2622,6 +3101,7 @@ //....................................................................... - - MWriteRootFile write(filenameDataout); + if (WXX) + { + MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -2641,5 +3121,5 @@ write.AddContainer("MEnergyEst", "Events"); - + } //----------------------------------------------------------------- @@ -2670,4 +3150,11 @@ eeston.SetCoeffA(parA); eeston.SetCoeffB(parB); + + //--------------------------- + // calculate estimated energy using Daniel's parameters + + //MEnergyEstParam eeston(fHilName); + //eeston.Add(fHilNameSrc); + //--------------------------- @@ -2766,4 +3253,536 @@ //--------------------------------------------------------------------- + + //--------------------------------------------------------------------- + // 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, WXX = " + << JobF_XX << ", " << WFX << 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; + + //----------------------------------------------------------- + // + // 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(); + + //....................................................................... + + if (WFX) + { + MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + + 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 + tliston.AddToList(&hfill6a); //timediff + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&eeston); + + if (WFX) + tliston.AddToList(&write); + + tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + //robert + tliston.AddToList(&hfill7); + tliston.AddToList(&hfill7a); + tliston.AddToList(&hfill7b); + tliston.AddToList(&hfill7c); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + 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); + + +//rwagner write all relevant histograms onto a file + + gLog << "=======================================================" << endl; + gLog << "Write results onto file '" << filenameResults << "'" << endl; + + TFile outfile(filenameResults,"recreate"); + + MHHillasSrc* hillasSrc = + (MHHillasSrc*)(pliston->FindObject("MHHillasSrc")); + TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + alphaHist->Write(); + + MHAlphaEnergyTheta* aetH = + (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); + TH3D* aetHist = (TH3D*)(aetH->GetHist()); + aetHist->SetName("aetHist"); + aetHist->Write(); + + MHAlphaEnergyTime* aetH2 = + (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime")); + TH3D* aetHist2 = (TH3D*)(aetH2->GetHist()); + aetHist2->SetName("aetimeHist"); +// aetHist2->DrawClone(); + aetHist2->Write(); + + MHThetabarTime* tbt = + (MHThetabarTime*)(pliston->FindObject("MHThetabarTime")); + TProfile* tbtHist = (TProfile*)(tbt->GetHist()); + tbtHist->SetName("tbtHist"); + tbtHist->Write(); + + MHEnergyTime* ent = + (MHEnergyTime*)(pliston->FindObject("MHEnergyTime")); + TH2D* entHist = (TH2D*)(ent->GetHist()); + entHist->SetName("entHist"); + entHist->Write(); + + MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); + TH2D* timeHist = (TH2D*)(time->GetHist()); + timeHist->SetName("MHTimeDiffTheta"); + timeHist->SetTitle("Time diffs"); + timeHist->Write(); + + MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); + TH2D* timeHist2 = (TH2D*)(time2->GetHist()); + timeHist2->SetName("MHTimeDiffTime"); + timeHist2->SetTitle("Time diffs"); + timeHist2->Write(); + +//rwagner write also collareas to same file + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + +//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 + + outfile.Close(); + + gLog << "=======================================================" << endl; + + + gLog << "Macro CT1Analysis : End of Job F_XX" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + } @@ -2781,2 +3800,15 @@ + + + + + + + + + + + + +