Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C (revision 5904) +++ (revision ) @@ -1,560 +1,0 @@ -//************************************************************************ -// -// 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