Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C (revision 5906) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C (revision 5906) @@ -0,0 +1,560 @@ +//************************************************************************ +// +// Authors : Eva Domingo, 12/2004 +// Wolfgang Wittek, 12/2004 +// +// +// Macro for estimating the DISP parameter +// --------------------------------------- +// +// DISP is the distance between the center of gravity of the shower image +// and the estimated source position, assumed to lie on the major +// axis of the shower +// +// In order to get an estimate of DISP +// - one assumes a parametrization DISP = f(p[0],p[1],p[2],... ; s1,s2,s3,...) +// where s1,s2,s3,... are measurable quantities like +// image parameters, (theta, phi), etc. +// and p[0],p[1],p[2], ... are free parameters +// +// - and determines the free parameters p[0],p[1],p[2]... from MC gamma data +// by minimizing a chi2 which measures the deviation of the estimated +// source position (using DISP) and the true source position +// +// The following classes are used : +// +// MDisp container holding the parameters p[0],p[1],p[2],... +// +// ::Calc member function calculating DISP for a given event; +// it is called by MDispCalc::Process() +// +// MDispCalc task calculating DISP with the parameters stored in MDisp +// +// MFindDisp class with several functions : +// +// ::DefineTrainMatrix \ member functions which generate the training +// ::DefineTestMatrix | and test samples (in the form of matrices) +// ::DefineTrainTestMatrix / from the MC gamma data +// +// ::FindParams member function steering the minimization +// (definition of the fcn function for Minuit, setting up the +// event loop to be executed in each minimization step, +// call to Minuit); +// for the minimization the training matrix is used +// +// ::TestParams member function testing the quality of the DISP estimate; +// for the test the test matrix is used +// +// MHDisp container for histograms which are useful for judging +// the quality of the DISP estimate +// Also task calculating the chi2 contribution of an event +// (in Fill()) +// and computing the final chi2 of the whole event sample +// (in Finalize()) +// +// MFDisp filter to select sample for the optimization +// +// +// The classes are stored in the CVS directory Mars/mtemp/mifae/library +// and Mars/mtemp/mifae/macros +// +// +//************************************************************************ + +void OptimizeDisp() +{ + //************************************************************************ + + Bool_t CMatrix = kFALSE; // Create training and test matrices + Bool_t WOptimize = kFALSE; // Do optimization using the training sample + // and write Disp parameter values + // onto the file parDispfile + // Optimize Disp with : + //TString typeOpt = "Data"; + TString typeOpt = "MC"; + Bool_t RTest = kFALSE; // Test the quality of the Disp estimate + // using the test matrix + Bool_t WDisp = kTRUE; // Make Disp plots for the data of type : + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; + + //************************************************************************ + + gLog.SetNoColors(); + + if (gRandom) + delete gRandom; + gRandom = new TRandom3(0); + + //----------------------------------------------- + //names of files to be read for generating the training and test matrices + const char *filetrain = "..."; + const char *filetest = "..."; + + //----------------------------------------------- + // path for input for Mars + TString inPathON = + "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_16/"; + TString inPathOFF = + "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_17/"; + TString inPathMC = + "~domingo/MC/"; + + // path for output from Mars + TString outPath = "/mnt/home/pcmagic04/domingo/DispOptimization/"; + + //----------------------------------------------- + // names of files for which Disp plots should be made + const char *offfile = "*OffCrabLow42*"; + const char *onfile = "*CrabLowEn42*"; + const char *mcfile = "starmc_Gammas_30-20_Period21.27ReducedMC.Zbins0-12"; + //----------------------------------------------- + + + //--------------------------------------------------------------------- + gLog << "=====================================================" << endl; + gLog << "Macro OptimizeDisp : Start " << endl; + + gLog << "" << endl; + gLog << "Macro OptimizeDisp : CMatrix, WOptimize, RTest, WDisp = " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "kTRUE" : "kFALSE") << ", " + << (WDisp ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // files to contain the matrices (generated from filenameTrain and + // filenameTest) + // + // for the training + TString fileMatrixTrain = outPath; + fileMatrixTrain += "MatrixTrainDisp"; + fileMatrixTrain += ".root"; + gLog << "" << endl; + gLog << "Files containing the Training and Test matrices :" << endl; + gLog << " fileMatrixTrain = " << fileMatrixTrain << endl; + + // for testing + TString fileMatrixTest = outPath; + fileMatrixTest += "MatrixTestDisp"; + fileMatrixTest += ".root"; + gLog << " fileMatrixTest = " << fileMatrixTest << endl; + gLog << "" << endl; + + + //--------------- + // file to contain the optimum Disp parameter values + TString parDispFile = outPath; + // parDispFile += "parDispMC.root"; + parDispFile += "parDispMC_withcuts.root"; + + gLog << "" << endl; + gLog << "File containing the optimum Disp parameter values : parDispFile = " + << parDispFile << endl; + + + // Set the filter cuts to select a sample of events for the Disp optimization + // (islandsmin,islandsmax,usedpixelsmin,usedpixelsmax,corepixelsmin, + // corepixelsmax,sizemin,sizemax,leakage1min,leakage1max,leakage2min, + // leakage2max,lengthmin,widthmin); + MFDisp *fdisp = NULL; + fdisp = new MFDisp; + fdisp->SetCuts(0,2,7,600,0,600,0.,3000.,0.,0.,0.,0.,0.,0.); + fdisp->SetName("FilterSelector2"); + + + // Create the MFindDisp object and set the file to store optimium Disp parameters + MFindDisp finddisp(fdisp); + finddisp.SetFilenameParam(parDispFile); + + + + //====================================================================== + // Create matrices and write them onto files + //====================================================================== + + if (CMatrix) + { + gLog << "-----------------------------------------------------" << endl; + gLog << "Generate the training and test samples" << endl; + + //-------------------------------------------- + // files to be read for optimizing the Disp parameters + // + // for the training + TString filenameTrain = inPathMC; + filenameTrain += mcfile; + filenameTrain += ".root"; + Int_t howManyTrain = 30000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = inPathMC; + filenameTest += mcfile; + filenameTest += ".root"; + Int_t howManyTest = 30000; + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // For the events in the matrices define requested distribution + // in some quantities; this may be a 1-dim, 2-dim or 3-dim distribution + + /* + // Define the cos(theta) distribution + TString mname("costheta"); + MBinning bincost("Binning"+mname); + bincost.SetEdges(10, 0., 1.0); + + //MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)"); + MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mname); + MH::SetBinning(&mh3.GetHist(), &bincost); + + for (Int_t i=1; i<=mh3.GetNbins(); i++) + mh3.GetHist().SetBinContent(i, 1.0); + */ + + // Define the log10(Etrue) distribution + TString mh3name("log10Etrue"); + MBinning binslogE("Binning"+mh3name); + binslogE.SetEdges(50, 1., 3.); + + MH3 mh3("log10(MMcEvt.fEnergy)"); + mh3.SetName(mh3name); + MH::SetBinning(&mh3.GetHist(), &binslogE); + + // Flat energy distribution + // for (Int_t i=1; i<=mh3.GetNbins(); i++) + // mh3.GetHist().SetBinContent(i, 1.0); + + // Power law energy distribution + // for (Int_t i=1; i<=mh3.GetNbins(); i++) + // { + // Double_t weight = pow((Double_t)i, -1.7); + // mh3.GetHist().SetBinContent(i, weight); + // } + + /* + // define the 'cos(Theta) vs. log10(Etrue)' distribution + TString mh3name("cosThetaVslog10Etrue"); + MBinning binslogE("Binning"+mh3name+"X"); + binslogE.SetEdges(18, 1.5, 2.2); + MBinning binscost("Binning"+mh3name+"Y"); + binscost.SetEdges(10, 0., 1.0); + + // MH3 mh3("log10(MMcEvt.fEnergy)", "cos((MPointingPos.fZd)/kRad2Deg)"); + MH3 mh3("log10(MMcEvt.fEnergy)", "cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mh3name); + MH::SetBinning((TH2*)&mh3.GetHist(), &binslogE, &binscost); + + for (Int_t i=1; i<=((TH2*)mh3.GetHist()).GetNbinsX(); i++) + { + // Double_t weight = pow((Double_t)i, -1.7); + for (Int_t j=1; j<=((TH2*)mh3.GetHist()).GetNbinsY(); j++) + { + // mh3.GetHist().SetBinContent(i, j, weight); + mh3.GetHist().SetBinContent(i, j, 1.); + } + } + */ + + //-------------------------- + // Training and test samples are generated from the same input files + if (filenameTrain == filenameTest) + { + if ( !finddisp.DefineTrainTestMatrix( + filenameTrain, mh3, + howManyTrain, howManyTest, + fileMatrixTrain, fileMatrixTest, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTrainTestMatrix failed" << endl; + return; + } + } + + //-------------------------- + // Training and test samples are generated from different input files + else + { + if ( !finddisp.DefineTrainMatrix(filenameTrain, mh3, + howManyTrain, fileMatrixTrain, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTrainMatrix failed" << endl; + return; + } + + if ( !finddisp.DefineTestMatrix( filenameTest, mh3, + howManyTest, fileMatrixTest, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTestMatrix failed" << endl; + return; + } + } + + gLog << "Generation of training and test samples finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + //====================================================================== + // Optimize Disp parameters using the training sample + // + // the initial values of the parameters are taken + // - from the file parDispInit (if != "") + // - or from the arrays params and steps (if their sizes are != 0) + // - or from the MDisp constructor + //====================================================================== + + if (WOptimize) + { + gLog << "-----------------------------------------------------" << endl; + gLog << "Optimize the Disp parameters over a " + << typeOpt << " sample, using the training matrix" << endl; + + // Read training matrix from file + finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest); + + //-------------------------------------------- + // file which contains the initial values for the Disp parameters; + // if parDispInit ="" + // - the initial values are taken from the arrays params and steps + // if their sizes are different from 0 + // - otherwise they are taken from the MDisp constructor + + TString parDispInit = outPath; + //parDispInit += "parDispInit.root"; + parDispInit = ""; + + //-------------------------------------------- + + TArrayD params(0); + TArrayD steps(0); + + if (parDispInit == "") + { + Double_t vparams[5] = { + // p[0], p[1], p[2], p[3], p[4], p[5] + 1.0, 0.6, -0.8, -0.8, -1.2/*, 0.5*/}; + // 0.5, 0., 0.03, 0., 0./*, 0.5*/}; + // 0.8, 0., 0., 0., 0./*, 0.5*/}; + + Double_t vsteps[5] = { + // dp[0], dp[1], dp[2], dp[3], dp[4], dp[5] + 0.01, 0.005, 0.005, 0.005, 0.01/*, 0.001*/}; + // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/}; + // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/}; + + params.Set(5, vparams); + steps.Set (5, vsteps ); + } + + + Bool_t rf; + rf = finddisp.FindParams(parDispInit, params, steps); + + if (!rf) + { + gLog << "OptimizeDisp.C : optimization of the Disp parameters failed" << endl; + return; + } + + gLog << "Optimization of Disp parameters finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + //====================================================================== + // Test the Disp parameters on the test sample + //====================================================================== + + if (RTest) + { + // Read test matrix from file + finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest); + + gLog << "-----------------------------------------------------" << endl; + gLog << "Test the Disp parameters using on the test matrix" << endl; + Bool_t rt = finddisp.TestParams(); + if (!rt) + { + gLog << "Test of the Disp parameters on the test matrix failed" + << endl; + } + gLog << "Test of the Disp parameters finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + + //====================================================================== + // Make Disp plots + //====================================================================== + + if (WDisp) + { + gLog << "" << endl; + gLog << "-----------------------------------------------------" << endl; + gLog << "Make plots for Disp for data of the type " << typeInput << endl; + + //-------------------------------------------- + // type of data to be read (ON, OFF or MC) + if (typeInput == "ON") + TString file(onfile); + else if (typeInput == "OFF") + TString file(offfile); + else if (typeInput == "MC") + TString file(mcfile); + + // name of input root file + TString filenameData = inPathMC; + filenameData += mcfile; + filenameData += ".root"; + gLog << "Input file '" << filenameData << endl; + + //---------------------------------------------------- + // read in optimum Disp parameter values + + TFile inparam(parDispFile); + MDisp dispin; + dispin.Read("MDisp"); + inparam.Close(); + + gLog << "Disp parameter values were read in from file '" + << parDispFile << "'" << endl; + + TArrayD dispPar; + dispPar = dispin.GetParameters(); + + TArrayD dispStep; + dispStep = dispin.GetStepsizes(); + + gLog << "optimum parameter values for calculating Disp : " << endl; + for (Int_t i=0; i