Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2474) +++ trunk/MagicSoft/Mars/Changelog (revision 2475) @@ -1,3 +1,28 @@ -*-*- END OF LINE -*-*- + + 2003/11/05: Wolfgang Wittek + + * macros/CT1Analysis.C + /ONOFFCT1Analysis.C + - current versions of the macros for the analysis of CT1 data + + * manalysis/MMarquardt.[h,cc] + - very pleliminary version of a class performing a minimization + using the Marquardt method + + * mimage/M2dimFunctionFit.[h,cc] + - very pleliminary version of a class which fits a 2-dim function to the + shower image using the maximum likelihood method + + * mimage/M2dimFunction.[h,cc] + - very pleliminary version of a container which contains the parameters + of the 2-dim function describing the shower image + + * mimage/MH2dimFunction.[h,cc] + - very pleliminary version of a container holding the histograms + for the parameters of the 2-dim function describing the shower image + + + 2003/11/05: Marcos Lopez @@ -205,4 +230,6 @@ * macros/weights.C - Added macro showing how to transform the spectrum of the MC showers. + + Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2474) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2475) @@ -241,7 +241,7 @@ Bool_t JobB_RF_UP = kFALSE; + Bool_t CTrainRF = kFALSE; // create matrices of training events Bool_t RTrainRF = kFALSE; // read in matrices of training events - Bool_t CTrainRF = kFALSE; // create matrices of training events - Bool_t RTree = kFALSE; // read in trees + Bool_t RTree = kFALSE; // read in trees (otherwise grow them) Bool_t WRF = kFALSE; // update input root file ? @@ -255,5 +255,6 @@ Bool_t JobB_SC_UP = kFALSE; - Bool_t RMatrix = kFALSE; // read training and test matrices from file + Bool_t CMatrix = kFALSE; // create training and test matrices + Bool_t RMatrix = kFALSE; // read training and test matrices from file Bool_t WOptimize = kFALSE; // do optimization using the training sample // and write supercuts parameter values @@ -1293,8 +1294,11 @@ rfwrite2.AddContainer("MRanTree", "TREE"); + MFillH fillh2("MHRanForestGini"); + // list of tasks for the loop over the trees tlist2.AddToList(&rfgrow2); tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); //------------------- @@ -1309,4 +1313,7 @@ tlist2.PrintStatistics(0, kTRUE); + + plist2.FindObject("MHRanForestGini")->DrawClone(); + // get adresses of objects which are used in the next eventloop @@ -1557,6 +1564,7 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_SC_UP, RMatrix, WOptimize, RTest, WSC = " + gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = " << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " << (RMatrix ? "kTRUE" : "kFALSE") << ", " << (WOptimize ? "kTRUE" : "kFALSE") << ", " @@ -1660,5 +1668,5 @@ //-------------------------- // create matrices and write them onto files - if (!RMatrix) + if (CMatrix) { MH3 &mh3 = *(new MH3("MHillas.fSize")); @@ -2111,5 +2119,5 @@ // - read ON1 and MC1 data files // which should have been updated to contain the hadronnesses - // for the method of NEAREST NEIGHBORS and for the SUOERCUTS + // for the method of Random Forest and for the SUPERCUTS // - produce Neyman-Pearson plots Index: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2474) +++ trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2475) @@ -1,2 +1,3 @@ + #include "CT1EgyEst.C" @@ -7,8 +8,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,7 +57,12 @@ MBinning *binth = new MBinning("BinningTheta"); + // 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}; + // this is our binning Double_t yedge[9] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; - TArrayD yed(9,yedge); + TArrayD yed; + yed.Set(9,yedge); binth->SetEdges(yed); plist->AddToList(binth); @@ -63,5 +74,6 @@ zedge[8-i] = cos(yedge[i]/kRad2Deg); } - TArrayD zed(9,zedge); + TArrayD zed; + zed.Set(9,zedge); bincosth->SetEdges(zed); plist->AddToList(bincosth); @@ -70,5 +82,19 @@ 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); } + void DeleteBinnings(MParList *plist) @@ -76,11 +102,11 @@ gLog << "DeleteBinnings" << endl; - if (!plist) - { - gLog << "Deletebinnins : MParlist no longer existing" << endl; - return; - } - TObject *bin; + + //-------------------------------------------- + bin = plist->FindObject("BinningE"); + if (bin) delete bin; + + //-------------------------------------------- bin = plist->FindObject("BinningSize"); @@ -121,6 +147,17 @@ if (bin) delete bin; - gLog << "exit DeleteBinnings" << endl; + + // 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; } + + //************************************************************************ @@ -162,6 +199,6 @@ // - write root file for ON (or OFF or MC) data (ON1.root, ...); - Bool_t JobA = kTRUE; - Bool_t WPad = kTRUE; // write out padding histograms ? + Bool_t JobA = kFALSE; + Bool_t WPad = kFALSE; // write out padding histograms ? Bool_t RPad = kFALSE; // read in padding histograms ? Bool_t Wout = kFALSE; // write out root file ON1.root @@ -178,18 +215,31 @@ Bool_t JobB_RF_UP = kFALSE; - Bool_t RTrainRF = kFALSE; // read in matrices of training events - Bool_t CTrainRF = kFALSE; // create matrices of training events - Bool_t RTree = kFALSE; // read in trees + Bool_t CTrainRF = kFALSE; // create matrices of training events + Bool_t RTrainRF = kFALSE; // read in matrices of training events + Bool_t RTree = kFALSE; // read in trees (otherwise grow them) Bool_t WRF = kFALSE; // update input root file ? + // Job B_SC_UP : read ON2.root (or MC2.root) file + // - depending on WParSC : create (or read in) supercuts parameter values + // - calculate hadroness for the SUPERCUTS + // - update the input files with the hadroness (==>ON3.root or MC3.root) + + Bool_t JobB_SC_UP = kFALSE; + Bool_t CMatrix = kFALSE; // create training and test matrices + Bool_t RMatrix = kFALSE; // read training and test matrices from 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 = kFALSE; // update input root file ? + // Job C: - // - read ON1.root and MC1.root files + // - read ON3.root and MC3.root files // which should have been updated to contain the hadronnesses // for the method of - // NEAREST NEIGHBORS + // RF // SUPERCUTS and - // RF // - produce Neyman-Pearson plots @@ -199,5 +249,5 @@ // Job D : // - select g/h separation method XX - // - read ON2 (or MC2) root file + // - read ON3 (or MC3) root file // - apply cuts in hadronness // - make plots @@ -207,28 +257,20 @@ - // Job E_EST_UP : - // - read MC1.root file + // Job E_XX : extended version of E_XX (including flux plots) // - select g/h separation method XX - // - optimize energy estimation for events passing the final cuts - // - write parameters of energy estimator onto file - // - update ON1.root, OFF1.root and MC1.root files with estimated energy - // (ON_XX1.root, OFF_XX1.root and MC_XX1.root) - - Bool_t JobE_EST_UP = kFALSE; - Bool_t WESTUP = kFALSE; // update root files ? - - - - // Job F_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 // - calculate flux - // - write root file for ON data after final cuts (ON3.root)) - - Bool_t JobF_XX = kFALSE; - Bool_t WXX = kFALSE; // write out root file ON3.root ? + // - write root file for ON data after final cuts + + + Bool_t JobE_XX = kTRUE; + Bool_t OEEst = kFALSE; // optimize energy estimator + Bool_t WEX = kTRUE; // update root file ? + Bool_t WRobert = kTRUE; // write out Robert's file ? + @@ -277,7 +319,7 @@ //-------------------------------------------------- // type of data to be padded - TString typeInput = "ON"; + //TString typeInput = "ON"; //TString typeInput = "OFF"; - //TString typeInput = "MC"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -740,7 +782,7 @@ // file to be updated (ON, OFF or MC) - TString typeInput = "ON"; + //TString typeInput = "ON"; //TString typeInput = "OFF"; - //TString typeInput = "MC"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -937,4 +979,5 @@ MProgressBar matrixbar; MEvtLoop evtloopg; + evtloopg.SetName("FillGammaMatrix"); evtloopg.SetParList(&plistg); //evtloopg.ReadEnv(env, "", printEnv); @@ -997,4 +1040,5 @@ MProgressBar matrixbar; MEvtLoop evtlooph; + evtlooph.SetName("FillHadronMatrix"); evtlooph.SetParList(&plisth); //evtlooph.ReadEnv(env, "", printEnv); @@ -1009,5 +1053,5 @@ - // write out matrices of training events events + // write out matrices of training events gLog << "" << endl; @@ -1119,8 +1163,11 @@ rfwrite2.AddContainer("MRanTree", "TREE"); + MFillH fillh2("MHRanForestGini"); + // list of tasks for the loop over the trees tlist2.AddToList(&rfgrow2); tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); //------------------- @@ -1135,4 +1182,7 @@ tlist2.PrintStatistics(0, kTRUE); + + plist2.FindObject("MHRanForestGini")->DrawClone(); + // get adresses of objects which are used in the next eventloop @@ -1212,8 +1262,4 @@ //----------------------------------------------------------------- - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - @@ -1304,4 +1350,5 @@ MProgressBar bar; MEvtLoop evtloop; + evtloop.SetName("UpdateRootFile"); evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); @@ -1373,40 +1420,361 @@ - //--------------------------------------------------------------------- - // Job C - //====== - - // - read ON1 and MC1 data files - // which should have been updated to contain the hadronnesses - // for the method of NEAREST NEIGHBORS and for the SUOERCUTS - // - produce Neyman-Pearson plots - - if (JobC) + // 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 C" << endl; + gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl; gLog << "" << endl; - gLog << "Macro CT1Analysis : JobC = " << JobC << endl; - - - // name of input data file + gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = " + << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " + << (RMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "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_2310a"; + + gLog << "parSCfile = " << parSCfile << endl; + + //-------------------------------------------- + // file to be updated (either ON or MC) + + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + // name of input root file TString filenameData = outPath; - filenameData += "OFF"; + filenameData += typeInput; filenameData += "2.root"; - gLog << "filenameData = " << filenameData << endl; - - // name of input MC file - TString filenameMC = outPath; - filenameMC += "MC"; - filenameMC += "2.root"; - gLog << "filenameMC = " << filenameMC << endl; + 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 optimizing the supercuts + // + // for the training + TString filenameTrain = outPath; + filenameTrain += "ON"; + filenameTrain += "1.root"; + Int_t howManyTrain = 800000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = outPath; + filenameTest += "ON"; + filenameTest += "1.root"; + Int_t howManyTest = 800000; + + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // files to contain the matrices (generated 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; + + + + //--------------------------------------------------------------------- + // Training and test matrices : + // - either create them and write them onto a file + // - or read them from a file + + + MCT1FindSupercuts findsuper; + findsuper.SetFilenameParam(parSCfile); + findsuper.SetHadronnessName("HadSC"); + + + //-------------------------- + // create matrices and write them onto files + if (CMatrix) + { + MH3 &mh3 = *(new MH3("MHillas.fSize")); + mh3.SetName("Target distribution for SIZE"); + + if (filenameTrain == filenameTest) + { + if ( !findsuper.DefineTrainTestMatrix(filenameTrain, + howManyTrain, mh3, howManyTest, mh3, + fileMatrixTrain, fileMatrixTest) ) + { + *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl; + return; + } + + } + else + { + if ( !findsuper.DefineTrainMatrix(filenameTrain, + howManyTrain, mh3, fileMatrixTrain) ) + { + *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl; + return; + } + + if ( !findsuper.DefineTestMatrix( filenameTest, + howManyTest, mh3, fileMatrixTest) ) + { + *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl; + return; + } + } + } + + //-------------------------- + // read matrices from files + // + + if (RMatrix) + findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); + //-------------------------- + + + + //--------------------------------------------------------------------- + // 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 (WOptimize) + { + gLog << "CT1Analysis.C : optimize the supercuts using the training matrix" + << endl; + + 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) + { + gLog << "CT1Analysis.C : optimization of supercuts failed" << endl; + 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; + } + + } //----------------------------------------------------------------- - + // Update the input files with the SC hadronness + // + + if (WSC) + { + 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 << "Parameter values for supercuts were read in from file '" + << parSCfile << "'" << endl; + + TArrayD supercutsPar; + supercutsPar = scin.GetParameters(); + + TArrayD supercutsStep; + supercutsStep = scin.GetStepsizes(); + + gLog << "Parameter values for supercuts : " << endl; + for (Int_t i=0; iDrawClone(); - pliston.FindObject("hadSC", "MHHadronness")->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 C" << endl; - gLog << "===================================================" << endl; - } - - - //--------------------------------------------------------------------- - // Job D - //====== - - // - select g/h separation method XX - // - read ON2 (or MC2) root file - // - apply cuts in hadronness - // - make plots - - - if (JobD) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job D" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobD = " - << JobD << endl; - - // type of data to be analysed - TString typeData = "ON"; - //TString typeData = "OFF"; - //TString typeData = "MC"; - gLog << "typeData = " << typeData << endl; - - TString ext = "2.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.30; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - gLog << "Maximum values of hadronness, |ALPHA| and DIST = " - << maxhadronness << ", " << maxalpha << ", " - << maxdist << endl; - - - //------------------------------ - // name of data file to be analysed - TString filenameData(outPath); - filenameData += typeData; - filenameData += ext; - gLog << "filenameData = " << filenameData << endl; - - - - //************************************************************************* - // - // Analyse the data - // - - MTaskList tliston; - MParList pliston; - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - - //----------------------------------------------------------------- - // 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"); - - - 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"); - - 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); - - tliston.AddToList(&contfinalgh); - - //tliston.AddToList(&fillhadnn); - tliston.AddToList(&fillhadsc); - tliston.AddToList(&fillhadrf); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - tliston.AddToList(&hfill5); - - 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(); - - //------------------------------------------- - // fit alpha distribution to get the number of excess events - // - - MHHillasSrc* hillasSrc = - (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); - TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); - - MHOnSubtraction onsub; - onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); - //------------------------------------------- - - DeleteBinnings(&pliston); - - gLog << "Macro CT1Analysis : End of Job D" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - - - //--------------------------------------------------------------------- - // Job E_EST_UP - //============ - - // - read MC2.root file - // - select g/h separation method XX - // - optimize energy estimator for events passing final cuts - // - write parameters of energy estimator onto file "energyest_XX.root" - // - // - read ON2.root and MC2.root files - // - update input root file with the estimated energy - // (ON_XX2.root, MC_XX2.root) - - - if (JobE_EST_UP) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " - << JobE_EST_UP << ", " << WESTUP << endl; - - - TString typeON = "ON"; - TString typeOFF = "OFF"; - TString typeMC = "MC"; - TString ext = "2.root"; - TString extout = "3.root"; - - //------------------------------ - // name of MC file to be used for optimizing the energy estimator - TString filenameOpt(outPath); - filenameOpt += typeMC; - filenameOpt += ext; - gLog << "filenameOpt = " << filenameOpt << endl; - - //------------------------------ - // 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 file containing the parameters of the energy estimator - TString energyParName(outPath); - energyParName += "energyest_"; - energyParName += XX; - energyParName += ".root"; - gLog << "energyParName = " << energyParName << endl; - - - //------------------------------ - // name of ON file to be updated - TString filenameON(outPath); - filenameON += typeON; - filenameON += ext; - gLog << "filenameON = " << filenameON << endl; - - // name of OFF file to be updated - TString filenameOFF(outPath); - filenameOFF += typeOFF; - filenameOFF += ext; - gLog << "filenameOFF = " << filenameOFF << endl; - - // name of MC file to be updated - TString filenameMC(outPath); - filenameMC += typeMC; - filenameMC += ext; - gLog << "filenameMC = " << filenameMC << endl; - - //------------------------------ - // name of updated ON file - TString filenameONup(outPath); - filenameONup += typeON; - filenameONup += "_"; - filenameONup += XX; - filenameONup += extout; - gLog << "filenameONup = " << filenameONup << endl; - - // name of updated OFF file - TString filenameOFFup(outPath); - filenameOFFup += typeOFF; - filenameOFFup += "_"; - filenameOFFup += XX; - filenameOFFup += extout; - gLog << "filenameOFFup = " << filenameOFFup << endl; - - // name of updated MC file - TString filenameMCup(outPath); - filenameMCup += typeMC; - filenameMCup += "_"; - filenameMCup += XX; - filenameMCup += extout; - gLog << "filenameMCup = " << filenameMCup << endl; - - //----------------------------------------------------------- - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //=========================================================== - // - // Optimization of energy estimator - // - - TString inpath(""); - TString outpath(""); - Int_t howMany = 2000; - CT1EEst(inpath, filenameOpt, outpath, energyParName, - fHilName, fHilNameSrc, fhadronnessName, - howMany, maxhadronness, maxalpha, maxdist); - - //----------------------------------------------------------- - // - // Read in parameters of energy estimator - // - gLog << "================================================================" - << endl; - gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" - << energyParName << "'" << endl; - TFile enparam(energyParName); - MMcEnergyEst mcest("MMcEnergyEst"); - mcest.Read("MMcEnergyEst"); - enparam.Close(); - - TArrayD parA(5); - TArrayD parB(7); - for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameOFF); - read.DisableAutoScheme(); - - //--------------------------- - // calculate estimated energy - - MEnergyEstParam eest2(fHilName); - eest2.Add(fHilNameSrc); - - eest2.SetCoeffA(parA); - eest2.SetCoeffB(parB); - - //....................................................................... - - MWriteRootFile write(filenameOFFup); + MWriteRootFile write(outNameImage, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -1951,16 +1825,469 @@ write.AddContainer("MNewImagePar", "Events"); - //write.AddContainer("HadNN", "Events"); - write.AddContainer("HadSC", "Events"); write.AddContainer("HadRF", "Events"); - - write.AddContainer("MEnergyEst", "Events"); + write.AddContainer(hadSCName, "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); + + TH1 &alphahist = alphaabs->GetHist(); + + 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(&supercuts); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadsc); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + 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]"); + alphaHist.SetName("alpha-macro"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + 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 C + //====== + + // - read ON1 and MC1 data files + // which should have been updated to contain the hadronnesses + // for the method of Random Forest and for the SUPERCUTS + // - produce Neyman-Pearson plots + + if (JobC) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job C" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobC = " + << (JobC ? "kTRUE" : "kFALSE") << endl; + + + // name of input data file + TString filenameData = outPath; + filenameData += "OFF"; + filenameData += "3.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of input MC file + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += "3.root"; + gLog << "filenameMC = " << filenameMC << 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", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + + //....................................................................... + // names of hadronness containers + + TString hadSCName = "HadSC"; + TString hadRFName = "HadRF"; + + //....................................................................... + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + 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 fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + 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); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadSC", "MHHadronness")->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 C" << endl; + gLog << "===================================================" << endl; + } + + + + //--------------------------------------------------------------------- + // Job D + //====== + + // - select g/h separation method XX + // - read ON2 (or MC2) root file + // - apply cuts in hadronness + // - make plots + + + if (JobD) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job D" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobD = " + << (JobD ? "kTRUE" : "kFALSE") << endl; + + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString ext = "3.root"; + + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //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.233; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + + //----------------------------------------------------------------- + // 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 fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + + TH1 &alphahist = alphaabs->GetHist(); + + 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"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); @@ -1968,15 +2295,27 @@ // entries in MParList - plistoff.AddToList(&tlistoff); - InitBinnings(&plistoff); - + pliston.AddToList(&tliston); + InitBinnings(&pliston); + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); //***************************** // entries in MTaskList - tlistoff.AddToList(&read); - tlistoff.AddToList(&eest2); - tlistoff.AddToList(&write); - tlistoff.AddToList(&contfinal); + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); //***************************** @@ -1987,28 +2326,351 @@ MProgressBar bar; MEvtLoop evtloop; - evtloop.SetParList(&plistoff); + evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); Int_t maxevents = -1; - //Int_t maxevents = 1000; + //Int_t maxevents = 10000; if ( !evtloop.Eventloop(maxevents) ) return; - tlistoff.PrintStatistics(0, kTRUE); - DeleteBinnings(&plistoff); - - //--------------------------------------------------- - - //--------------------------------------------------- - // Update ON data + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + 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(); + + + //------------------------------------------- + + // 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]"); + alphaHist.SetName("alpha-JobD"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + 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 D" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + + + //--------------------------------------------------------------------- + // Job E_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 (JobE_XX) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job E_XX" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " + << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (OEEst ? "kTRUE" : "kFALSE") << ", " + << (WEX ? "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("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.23; + 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"; + + + if (OEEst) + { + //=========================================================== + // + // 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 (WEX) + { + //----------------------------------------------------------- + // + // 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; + + + gLog << "Read in Migration matrix" << endl; + + MHMcEnergyMigration mighiston("MHMcEnergyMigration"); + mighiston.Read("MHMcEnergyMigration"); + //MHMcEnergyMigration &mighiston = + // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration")); + + gLog << "Migration matrix was read in" << endl; + + + TArrayD parA(mcest.GetNumCoeffA()); + TArrayD parB(mcest.GetNumCoeffB()); + for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); + //------------------------------------------- // create the tasks which should be executed // - MReadMarsFile read("Events", filenameON); + MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); - //--------------------------- - // calculate estimated energy - - MEnergyEstParam eest2(fHilName); - eest2.Add(fHilNameSrc); - - eest2.SetCoeffA(parA); - eest2.SetCoeffB(parB); - //....................................................................... - MWriteRootFile write(filenameONup); + gLog << "CT1Analysis.C : write root file '" << filenameDataout + << "'" << endl; + + //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + + + MWriteRootFile write(filenameDataout, "RECREATE"); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -2053,362 +2713,4 @@ write.AddContainer("MEnergyEst", "Events"); - //----------------------------------------------------------------- - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(fhadronnessName); - MContinue contfinal(&selfinal); - - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - tliston.AddToList(&eest2); - tliston.AddToList(&write); - tliston.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - DeleteBinnings(&pliston); - - //--------------------------------------------------- - - //--------------------------------------------------- - // Update MC data - // - gLog << "============================================================" - << endl; - gLog << "Macro CT1Analysis.C : update file '" << filenameMC - << "'" << endl; - - MTaskList tlistmc; - MParList plistmc; - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameMC); - read.DisableAutoScheme(); - - //--------------------------- - // calculate estimated energy - - MEnergyEstParam eest2(fHilName); - eest2.Add(fHilNameSrc); - - eest2.SetCoeffA(parA); - eest2.SetCoeffB(parB); - - //....................................................................... - - MWriteRootFile write(filenameMCup); - - 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"); - - //----------------------------------------------------------------- - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(fhadronnessName); - MContinue contfinal(&selfinal); - - - //***************************** - // entries in MParList - - plistmc.AddToList(&tlistmc); - InitBinnings(&plistmc); - - - //***************************** - // entries in MTaskList - - tlistmc.AddToList(&read); - tlistmc.AddToList(&eest2); - tlistmc.AddToList(&write); - tlistmc.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&plistmc); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tlistmc.PrintStatistics(0, kTRUE); - DeleteBinnings(&plistmc); - - - //========== end update ============================================ - } - - enparam.Close(); - - gLog << "Macro CT1Analysis : End of Job E_EST_UP" << 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, WXX = " - << JobF_XX << ", " << WXX << 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 calculating the eff. collection areas - TString filenameArea(outPath); - filenameArea += typeMC; - filenameArea += "_"; - filenameArea += XX; - 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 += "_"; - filenameData += XX; - 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; - - - //==================================================================== - 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 = new MHMcCT1CollectionArea(); - //MHMcCT1CollectionArea* collarea; - - 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("lego"); - - //--------------------------------------------- - // Write histograms to a file - // - - TFile f(collareaName, "RECREATE"); - collarea->GetHist()->Write(); - collarea->GetHAll()->Write(); - collarea->GetHSel()->Write(); - f.Close(); - - //delete collarea; - - gLog << "Calculation of effective collection areas done" << endl; - gLog << "-----------------------------------------------" << endl; - //------------------------------------------------------------------ - - - //************************************************************************* - // - // Analyse the data - // - - MTaskList tliston; - MParList pliston; - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - //....................................................................... - - - MWriteRootFile write(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"); - //----------------------------------------------------------------- @@ -2431,4 +2733,24 @@ 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); @@ -2446,4 +2768,28 @@ 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); @@ -2467,8 +2813,12 @@ tliston.AddToList(&read); + // robert + tliston.AddToList(&hfill6); //timediff + tliston.AddToList(&hfill6a); //timediff + tliston.AddToList(&contfinalgh); - - if (WXX) - tliston.AddToList(&write); + tliston.AddToList(&eeston); + + tliston.AddToList(&write); //tliston.AddToList(&fillhadnn); @@ -2482,4 +2832,10 @@ tliston.AddToList(&hfill5); + //robert + tliston.AddToList(&hfill7); + tliston.AddToList(&hfill7a); + tliston.AddToList(&hfill7b); + tliston.AddToList(&hfill7c); + tliston.AddToList(&contfinal); @@ -2495,5 +2851,4 @@ Int_t maxevents = -1; - //Int_t maxevents = 10000; if ( !evtloop.Eventloop(maxevents) ) return; @@ -2507,16 +2862,113 @@ //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 << "Macro CT1Analysis : End of Job F_XX" << endl; + 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; + + TFile outfile(filenameResults,"recreate"); + + MHHillasSrc* hillasSrc = + (MHHillasSrc*)(pliston->FindObject("MHHillasSrc")); + TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + alphaHist->Write(); + gLog << "Alpha plot has been written out" << endl; + + + MHAlphaEnergyTheta* aetH = + (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); + TH3D* aetHist = (TH3D*)(aetH->GetHist()); + aetHist->SetName("aetHist"); + aetHist->Write(); + gLog << "AlphaEnergyTheta plot has been written out" << endl; + + MHAlphaEnergyTime* aetH2 = + (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime")); + TH3D* aetHist2 = (TH3D*)(aetH2->GetHist()); + aetHist2->SetName("aetimeHist"); +// aetHist2->DrawClone(); + aetHist2->Write(); + gLog << "AlphaEnergyTime plot has been written out" << endl; + + MHThetabarTime* tbt = + (MHThetabarTime*)(pliston->FindObject("MHThetabarTime")); + TProfile* tbtHist = (TProfile*)(tbt->GetHist()); + tbtHist->SetName("tbtHist"); + tbtHist->Write(); + gLog << "ThetabarTime plot has been written out" << endl; + + MHEnergyTime* ent = + (MHEnergyTime*)(pliston->FindObject("MHEnergyTime")); + TH2D* entHist = (TH2D*)(ent->GetHist()); + entHist->SetName("entHist"); + entHist->Write(); + gLog << "EnergyTime plot has been written out" << endl; + + MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); + TH2D* timeHist = (TH2D*)(time->GetHist()); + timeHist->SetName("MHTimeDiffTheta"); + timeHist->SetTitle("Time diffs"); + timeHist->Write(); + gLog << "TimeDiffTheta plot has been written out" << endl; + + + MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); + TH2D* timeHist2 = (TH2D*)(time2->GetHist()); + timeHist2->SetName("MHTimeDiffTime"); + timeHist2->SetTitle("Time diffs"); + timeHist2->Write(); + gLog << "TimeDiffTime plot has been written out" << endl; + +//rwagner write also collareas to same file + collarea->GetHist()->Write(); + 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 E_XX" << endl; gLog << "=======================================================" << endl; } @@ -2527,2 +2979,28 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MMarquardt.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMarquardt.cc (revision 2475) +++ trunk/MagicSoft/Mars/manalysis/MMarquardt.cc (revision 2475) @@ -0,0 +1,414 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek 10/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MMarquardt // +// // +// Marquardt method of solving nonlinear least-squares problems // +// // +// (see Numerical recipes (2nd ed.), W.H.Press et al., p.688 ff) // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MMarquardt.h" + +#include // fabs + +#include +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MParContainer.h" + +ClassImp(MMarquardt); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMarquardt::MMarquardt(const char *name, const char *title) +{ + fName = name ? name : "MMarquardt"; + fTitle = title ? title : "Marquardt minimization"; +} + +// ----------------------------------------------------------------------- +// +// Set - the number of parameters +// - the maximum number of steps allowed in the minimization and +// - the change in chi2 signaling convergence + +void MMarquardt::SetNpar(Int_t numpar, Int_t numstepmax, Double_t loopchi2) +{ + fNpar = numpar; + fNumStepMax = numstepmax; + fLoopChi2 = loopchi2; + + fdParam.ResizeTo(fNpar); + + fParam.ResizeTo(fNpar); + fGrad.ResizeTo(fNpar); + fCovar.ResizeTo(fNpar, fNpar); + + fmyParam.ResizeTo(fNpar); + fmyGrad.ResizeTo(fNpar); + fmyCovar.ResizeTo(fNpar, fNpar); + + fIxc.ResizeTo(fNpar); + fIxr.ResizeTo(fNpar); + fIp.ResizeTo(fNpar); +} + + +// ----------------------------------------------------------------------- +// +// do the minimization +// +// fcn is the function which calculates for a given set of parameters +// - the function L to be minimized +// - beta_k = -1/2 * dL/da_k (a kind of gradient of L) +// - alfa_kl = 1/2 * dL/(da_k da_l) (a kind of 2nd derivative of L) +// +// Vinit contains the starting values of the parameters +// + +Int_t MMarquardt::Loop( + Bool_t (*fcn)(TVectorD &, TMatrixD &, TVectorD &, Double_t &), + TVectorD &Vinit) +{ + fFunc = fcn; + + // set the initial parameter values + for (Int_t i=0; i fLoopChi2 && fNumStep < fNumStepMax); + + //------------------------------------------- + // do the final calculations + if (!rcnext) + { + *fLog << "MMarquardt::Loop; function call failed in step " << fNumStep + << endl; + return -2; + } + + if (fdChi2 > fLoopChi2) + { + *fLog << "MMarquardt::Loop; minimization has not converged, fChi2, fdChi2 = " + << fChi2 << ", " << fdChi2 << endl; + return -3; + } + + *fLog << "MMarquardt::Loop; minimization has converged, fChi2, fdChi2, fNumStep = " + << fChi2 << ", " << fdChi2 << ", " << fNumStep << endl; + + + Bool_t rccov = CovMatrix(); + if (!rccov) + { + *fLog << "MMarquardt::Loop; calculation of covariance matrix failed " + << endl; + return 1; + } + + //------------------------------------------- + // results + + *fLog << "MMarquardt::Loop; Results of fit : fChi2, fNumStep, fdChi2 =" + << fChi2 << ", " << fNumStep << ", " << fdChi2 << endl; + + for (Int_t i=0; i= h) + { + h = fabs(covar(j,k)); + ir = j; + ic = k; + } + } + else + { + if (fIp(k) > 1) return kFALSE; + } + } + } + } + + fIp(ic)++; + if (ir != ic) + { + for (l=0; l=0; l--) + { + if (fIxr(l) != fIxc(l)) + { + for (k=0; k +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// M2dimFunction +// +// Storage Container for the parameters of the 2-dim function describing +// the shower image +// +// +///////////////////////////////////////////////////////////////////////////// +#include "M2dimFunction.h" + +#include > + +#include "TMath.h" +#include "TVectorD.h" +#include "TMatrixD.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHillas.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MCerPhotEvt.h" +#include "MCerPhotPix.h" + +ClassImp(M2dimFunction); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// TwodimFunction +// +// this function calculates for a given set of parameters +// +// - the log likelihood function L to be minimized +// - beta_k = -1/2 * dL/da_k (a kind of gradient of L) +// - alfa_kl = 1/2 * dL/(da_k da_l) (a kind of second derivative of L) +// +// (see Numerical recipes (2nd ed.), W.H.Press et al., p. 687) +// +// Note : this is not a member function of M2dimFunction; +// it will be called by the minimization class MMarquardt; +// the address of this function is passed to MMarquardt +// in M2dimFunction::Fit() by the call MMarquardt::Loop() +// + +Bool_t TwodimFunction(TVectord &a, TMatrixD &alfa, TVectord &beta, Double_t &L) +{ + Int_t npar = a.GetNrows(); + Int_t npixels =; + + TMatrixD dmuda (npar, npixels); + TVector dLdmu (npixels); + TVector d2Ldmu2(npixels); + + //------------------------------------------ + // these are the parameters for which alfa, beta and L are calculated + + Double_t xbar = a(0); + Double_t ybar = a(1); + Double_t delta = a(2); + Double_t amp = a(3); + Double_t leng = a(4); + Double_t wid = a(5); + Double_t asy = a(6); + + + for (Int_t m=0; mGetMeanX(); + fVinit(1) = fhillas->GetMeanY(); + fVinit(2) = fhillas->GetDelta(); + } + else + { + fVinit(0) = ; + fVinit(1) = ; + fVinit(2) = ; + } + + fVinit(3) = + fVinit(4) = + fVinit(5) = + fVinit(6) = + + return; +} + + + +// -------------------------------------------------------------------------- +// +// Fit of the 2-dim function to the shower image +// +void M2dimFunction::Fit() +{ + fMarquardt.Loop(TwodimFunction, fVinit); + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +void M2dimFunction::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Parameters of 2-dim function (" << GetName() << ")" << endl; + *fLog << " - fXbar = " << fXbar << endl; + *fLog << " - fYbar = " << fYbar << endl; + *fLog << " - fAmp = " << fAmp << endl; + *fLog << " - fMajor = " << fMajor << endl; + *fLog << " - fMinor = " << fMinor << endl; + *fLog << " - fAsym = " << fAsym << endl; +} +//============================================================================ + + + + + + + + + + Index: trunk/MagicSoft/Mars/mimage/M2dimFunction.h =================================================================== --- trunk/MagicSoft/Mars/mimage/M2dimFunction.h (revision 2475) +++ trunk/MagicSoft/Mars/mimage/M2dimFunction.h (revision 2475) @@ -0,0 +1,66 @@ +#ifndef MARS_M2dimFunction +#define MARS_M2dimFunction + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MHillas; +class MGeomCam; +class MCerPhotEvt; + +class M2dimFunction : public MParContainer +{ +private: + + // free parameters of the 2-dim function describing the shower image + Double_t fXbar, fYbar; // position where shower image has its maximum + Double_t fDelta; // angle between shower axis and x-axis + Double_t fAmp; // amplitude of shower image + Double_t fMajor; fMinor; // parameters describing the extension of + // the shower along and perpendicular + // to the shower axis + Double_t fAsym; // parameter describing the asymmetry along the + // shower axis + + TVectorD Vinit; // initial values of the parameters + +public: + M2dimFunction(const char *name=NULL, const char *title=NULL); + + void SetVinit(MHillas *fhillas); + + Fit(); + + ClassDef(M2dimFunction, 1) // Container to hold the parameters of the + // 2-dim function describing the shower image +}; + +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.cc (revision 2475) +++ trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.cc (revision 2475) @@ -0,0 +1,145 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek 10/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// M2dimFunctionFit +// +// Task to fit a 2-dim function to the shower image +// +////////////////////////////////////////////////////////////////////////////// +#include "M2dimFunctionFit.h" + +#include + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MSrcPosCam.h" +#include "MCerPhotEvt.h" +#include "MNewImagePar.h" +#include "MNewImagePar.h" +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(M2dimFunctionFit); + +static const TString gsDefName = "M2dimFunctionFit"; +static const TString gsDefTitle = "Fit 2-dim function"; + +// ------------------------------------------------------------------------- +// +// Default constructor. +// +// twodimfunname is the name of a container of type M2dimFunction, +// in which the parameters are stored; the default is "M2dimFunction" +// +// +M2dimFunctionFit::M2dimFunctionFit(const char *twodimfunname, + const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + f2dimFunName = twodimfunname; + fHillasInput = "MHillas"; +} + +// ------------------------------------------------------------------------- +// +Int_t M2dimFunctionFit::PreProcess(MParList *pList) +{ + fHillas = (MHillas*)pList->FindObject(fHillasInput, "MHillas"); + if (!fHillas) + { + *fLog << err << dbginf << "MHillas not found... aborting." << endl; + return kFALSE; + } + + + fCerPhotEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fCerPhotEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeomCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + f2dimfun = (M2dimFunction*)pList->FindCreateObj("M2dimFunction", f2dimFunName); + if (!f2dimfun) + { + *fLog << dbginf << "M2dimFunction object '" << f2dimFunName + << "' not found... aborting." << endl; + return kFALSE; + } + + fErrors = 0; + + return kTRUE; +} + +// ------------------------------------------------------------------------- +// +// fit 2-dim function to the shower image +// + +Int_t M2dimFunctionFit::Process() +{ + Bool_t rc = f2dimfun->Fit(); + if (!rc) + { + fErrors++; + return kCONTINUE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the hillas calculation. The percentage +// is calculated with respect to the number of executions of this task. +// +/* +Bool_t M2dimFunctionFit::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << fErrors << " (" << (int)(fErrors*100/GetNumExecutions()) << "%) Evts skipped due to: fit of 2-dim function failed" << endl; + *fLog << endl; + + return kTRUE; +} +*/ + Index: trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.h =================================================================== --- trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.h (revision 2475) +++ trunk/MagicSoft/Mars/mimage/M2dimFunctionFit.h (revision 2475) @@ -0,0 +1,52 @@ +#ifndef MARS_M2dimFunctionFit +#define MARS_M2dimFunctionFit + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MHillas; +class MNewImagePar; +class MSrcPosCam; +class MGeomCam; +class MCerPhotEvt; + +class M2dimFunctionFit : public MTask +{ +private: + MGeomCam *fGeomCam; + MCerPhotEvt *fCerPhotEvt; + + MHillas *fHillas; //! Pointer to the source independent hillas parameters + MSrcPosCam *fSrcPos; //! Pointer to the source position + MNewImagePar *fNewImagePar; //! Pointer to the output container for the new image parameters + + TString fSrcName; + TString fNewParName; + TString fHillasInput; + + Int_t fErrors; + + Bool_t PreProcess(MParList *plist); + Bool_t Process(); + Bool_t PostProcess(); + +public: + M2dimFunctionFit(const char *twodimfunname, + const char *name=NULL, const char *title=NULL); + + void SetInput(TString hilname) { fHillasInput = hilname; } + + ClassDef(M2dimFunctionFit, 0) // task to fit a 2-dim function to the shower image +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/mimage/MH2dimFunction.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MH2dimFunction.cc (revision 2475) +++ trunk/MagicSoft/Mars/mimage/MH2dimFunction.cc (revision 2475) @@ -0,0 +1,205 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 10/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MH2dimFunction +// +// container of histograms for the parameters of the 2-dim function +// describing the shower image +// +//////////////////////////////////////////////////////////////////////////// +#include "MH2dimFunction.h" + +#include + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" +#include "MBinning.h" +#include "MParList.h" + +#include "MHillas.h" +#include "M2dimFunction.h" + +ClassImp(MH2dimFunction); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup histograms +// +MH2dimFunction::MH2dimFunction(const char *name, const char *title) +{ + fName = name ? name : "MH2dimFunction"; + fTitle = title ? title : "Histograms of parameters of 2-dim function"; + + fHistXbarYbar.SetName("shower-maximum"); + fHistXbarYbar.SetTitle("position of maximum of shower"); + fHistXbarYbar.SetXTitle("x-position of maximum of shower image"); + fHistXbarYbar.SetYTitle("y-position of maximum of shower image"); + fHistXbarYbar.SetDirectory(NULL); + fHistXbarYbar.SetFillStyle(4000); + fHistXbarYbar.UseCurrentStyle(); + + fHistAmp.SetName("Amplitude parameter"); + fHistAmp.SetTitle("amplitude parameter"); + fHistAmp.SetXTitle("amplitude parameter"); + fHistAmp.SetYTitle("Counts"); + fHistAmp.SetDirectory(NULL); + fHistAmp.SetFillStyle(4000); + fHistAmp.UseCurrentStyle(); + + fHistMajor.SetName("Length parameter"); + fHistMajor.SetTitle("length parameter"); + fHistMajor.SetXTitle("length parameter"); + fHistMajor.SetYTitle("Counts"); + fHistMajor.SetDirectory(NULL); + fHistMajor.SetLineColor(kBlue); + fHistMajor.SetFillStyle(4000); + fHistMajor.UseCurrentStyle(); + + fHistMinor.SetName("Width parameter"); + fHistMinor.SetTitle("width parameter"); + fHistMinor.SetXTitle("width parameter"); + fHistMinor.SetYTitle("Counts"); + fHistMinor.SetDirectory(NULL); + fHistMinor.SetLineColor(kGreen); + fHistMinor.SetFillStyle(4000); + fHistMinor.UseCurrentStyle(); + + fHistAsym.SetName("Asymmetry parameter"); + fHistAsym.SetTitle("asymmetry parameter"); + fHistAsym.SetXTitle("asymmetry parameter"); + fHistAsym.SetYTitle("Counts"); + fHistAsym.SetDirectory(NULL); + fHistAsym.SetFillStyle(4000); + fHistAsym.UseCurrentStyle(); + + MBinning bins; + + bins.SetEdges(100, 0, 1); + bins.Apply(fHistAmp); + + bins.SetEdges(100, 0, 1); + bins.Apply(fHistMajor); + + bins.SetEdges(100, 0, 1); + bins.Apply(fHistMinor); + + bins.SetEdges(100, 0, 1); + bins.Apply(fHistAsym); +} + +// -------------------------------------------------------------------------- +// +// Setup the Binning for the histograms automatically if the correct +// instances of MBinning +// +Bool_t MH2dimFunction::SetupFill(const MParList *plist) +{ + ApplyBinning(*plist, "XbarYbar", &fHistXbarYbar); + ApplyBinning(*plist, "Amp", &fHistAmp); + ApplyBinning(*plist, "Major", &fHistMajor); + ApplyBinning(*plist, "Minor", &fHistMinor); + ApplyBinning(*plist, "Asym", &fHistAsym); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Fill the histograms with data from a M2dimFunction container. +// +Bool_t MH2dimFunction::Fill(const MParContainer *par, const Stat_t w) +{ + if (!par) + { + *fLog << err << "MH2dimFunction::Fill: Pointer (!=NULL) expected." << endl; + return kFALSE; + } + + const M2dimFunction &h = *(M2dimFunction*)par; + + fHistXbarYbar.Fill(h.GetXbar(), h.GetYbar(), w); + fHistAmp.Fill(h.GetAmp(), w); + + fHistMajor.Fill(h.GetMajor(), w); + fHistMinor.Fill(h.GetMinor(), w); + + fHistAsym.Fill(h.GetAsym(), w); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the two histograms into it. +// Be careful: The histograms belongs to this object and won't get deleted +// together with the canvas. +// +void MH2dimFunction::Draw(Option_t *) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + pad->Divide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + TAxis &x = *fHistLeakage1.GetXaxis(); + x.SetRangeUser(0.0, x.GetXmax()); + MH::DrawSame(fHistLeakage1, fHistLeakage2, "Leakage1 and Leakage2"); + + pad->cd(2); + gPad->SetBorderMode(0); + MH::DrawSame(fHistCorePix, fHistUsedPix, "Number of core/used Pixels"); + + pad->cd(3); + gPad->SetBorderMode(0); + MH::DrawSame(fHistConc1, fHistConc, "Concentrations"); + + pad->cd(4); + gPad->SetBorderMode(0); + + pad->Modified(); + pad->Update(); +} + +//========================================================================== + + + + + + Index: trunk/MagicSoft/Mars/mimage/MH2dimFunction.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MH2dimFunction.h (revision 2475) +++ trunk/MagicSoft/Mars/mimage/MH2dimFunction.h (revision 2475) @@ -0,0 +1,47 @@ +#ifndef MARS_MH2dimFunction +#define MARS_MH2dimFunction + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef ROOT_TH1 +#include +#endif + +class MHillas; + +class MH2dimFunction : public MH +{ +private: + + TH1F fHistXbarYbar; // position where shower image has its maximum + TH1F fHistAmp; // amplitude of shower image + TH1F fHistMajor; // 'length' of shower image + TH1F fHistMinor; // 'width' of shower image + TH1F fHistAsym; // 'asymmetry' of shower image + + +public: + MH2dimFunction(const char *name=NULL, const char *title=NULL); + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + TH1F &GetHistXbarYbar() { return fHistXbarYbar; } + TH1F &GetHistAmp() { return fHistAmp; } + + TH1F &GetHistMajor() { return fHistMajor; } + TH1F &GetHistMinor() { return fHistMinor; } + + TH1F &GetHistAsym() { return fHistAsym; } + + void Draw(Option_t *opt=NULL); + + ClassDef(MH2dimFunction, 1) // Histograms of parameters of 2-dim function +}; + +#endif + + + +