Index: /trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- /trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2435) +++ /trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2436) @@ -210,6 +210,6 @@ //----------------------------------------------- - TEnv env("macros/CT1env.rc"); - Bool_t printEnv = kFALSE; + //TEnv env("macros/CT1env.rc"); + //Bool_t printEnv = kFALSE; //************************************************************************ @@ -219,7 +219,7 @@ // - write root file for ON data (ON1.root); - Bool_t JobA_ON = kTRUE; - Bool_t WHistON = kTRUE; // write out histogram sigmabar vs. Theta ? - Bool_t WON1 = kTRUE; // write out root file ON1.root ? + Bool_t JobA_ON = kFALSE; + Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? + Bool_t WON1 = kFALSE; // write out root file ON1.root ? @@ -234,11 +234,13 @@ // Job B_RF_UP : read ON1.root (or MC1.root) file - // - depending on RLook : create (or read in) hadron and gamma matrix - // - depending on RTree : create (or read in) trees + // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas + // - if CTrainRF = TRUE : create training matrices for hadrons and gammas + // - if RTree = TRUE : read in trees, otherwise create trees // - calculate hadroness for method of RANDOM FOREST - // - update the input files with the hadroness (ON1.root or MC1.root) + // - update the input files with the hadroness (==>ON2.root or MC2.root) Bool_t JobB_RF_UP = kFALSE; - Bool_t RLookRF = kFALSE; // read in look-alike 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 WRF = kFALSE; // update input root file ? @@ -247,8 +249,8 @@ - // Job B_SC_UP : read ON1.root (or MC1.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 (ON1.root or MC1.root) + // - update the input files with the hadroness (==>ON3.root or MC3.root) Bool_t JobB_SC_UP = kFALSE; @@ -263,10 +265,9 @@ // 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 @@ -276,5 +277,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 @@ -296,5 +297,6 @@ Bool_t JobE_XX = kFALSE; - Bool_t WEX = kFALSE; // write out root file ? + Bool_t OEEst = kFALSE; // optimize energy estimator + Bool_t WEX = kFALSE; // update root file ? Bool_t WRobert = kFALSE; // write out Robert's file ? @@ -492,5 +494,5 @@ MEvtLoop evtloop; evtloop.SetParList(&pliston); - evtloop.ReadEnv(env, "", printEnv); + //evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); //if (WON1) @@ -822,5 +824,5 @@ MEvtLoop evtloop; evtloop.SetParList(&plist); - evtloop.ReadEnv(env, "", printEnv); + //evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); //if (WMC1) @@ -863,5 +865,5 @@ - // - create (or read in) the matrices of look-alike events for gammas + // - create (or read in) the matrices of training events for gammas // and hadrons // - create (or read in) the trees @@ -877,7 +879,8 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " + gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = " << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " - << (RLookRF ? "kTRUE" : "kFALSE") << ", " + << (RTrainRF ? "kTRUE" : "kFALSE") << ", " + << (CTrainRF ? "kTRUE" : "kFALSE") << ", " << (RTree ? "kTRUE" : "kFALSE") << ", " << (WRF ? "kTRUE" : "kFALSE") << endl; @@ -890,10 +893,21 @@ Int_t NdSize = 1; + // cut in RF hadronness + TString hadRFName = "HadRF"; + Float_t maxhadronness = 0.23; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + //-------------------------------------------- // file to be updated (either ON or MC) - TString typeInput = "ON"; - //TString typeInput = "MC"; + //TString typeInput = "ON"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -913,5 +927,5 @@ //-------------------------------------------- - // files to be read for generating the look-alike events + // files to be read for generating the matrices of training events // "hadrons" : TString filenameON = outPath; @@ -932,5 +946,5 @@ //-------------------------------------------- - // files of look-alike events + // files of training events TString outNameGammas = outPath; @@ -951,7 +965,20 @@ matrixg.EnableGraphicalOutput(); + matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrixg.AddColumn("MSigmabar.fSigmabar"); + matrixg.AddColumn("log10(MHillas.fSize)"); + matrixg.AddColumn("MHillasSrc.fDist"); + matrixg.AddColumn("MHillas.fWidth"); + matrixg.AddColumn("MHillas.fLength"); + matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + matrixg.AddColumn("MNewImagePar.fConc"); + matrixg.AddColumn("MNewImagePar.fLeakage1"); + MHMatrix matrixh("MatrixHadrons"); matrixh.EnableGraphicalOutput(); + matrixh.AddColumns(matrixg.GetColumns()); + //-------------------------------------------- // file of trees of the random forest @@ -962,6 +989,6 @@ //************************************************************************* - // read in matrices of look-alike events -if (RLookRF) + // read in matrices of training events +if (RTrainRF) { const char* mtxName = "MatrixGammas"; @@ -1005,10 +1032,10 @@ //************************************************************************* - // create matrices of look-alike events -if (!RLookRF) + // create matrices of training events +if (CTrainRF) { gLog << "" << endl; gLog << "========================================================" << endl; - gLog << " Create matrices of look-alike events" << endl; + gLog << " Create matrices of training events" << endl; gLog << " Gammas :" << endl; @@ -1085,5 +1112,5 @@ MEvtLoop evtloopg; evtloopg.SetParList(&plistg); - evtloopg.ReadEnv(env, "", printEnv); + //evtloopg.ReadEnv(env, "", printEnv); evtloopg.SetProgressBar(&matrixbar); @@ -1141,5 +1168,5 @@ MEvtLoop evtlooph; evtlooph.SetParList(&plisth); - evtlooph.ReadEnv(env, "", printEnv); + //evtlooph.ReadEnv(env, "", printEnv); evtlooph.SetProgressBar(&matrixbar); @@ -1155,9 +1182,9 @@ - // write out look-alike events + // write out training events gLog << "" << endl; gLog << "========================================================" << endl; - gLog << "Write out look=alike events" << endl; + gLog << "Write out training events" << endl; @@ -1171,5 +1198,5 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " + gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file " << outNameGammas << endl; @@ -1183,9 +1210,9 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " + gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file " << outNameHadrons << endl; } - //********** end of creating matrices of look-alike events *********** + //********** end of creating matrices of training events *********** @@ -1301,9 +1328,7 @@ - //----------------------------------------------------------------- // Update the input files with the RF hadronness // - if (WRF) { @@ -1328,13 +1353,9 @@ read.DisableAutoScheme(); - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; //....................................................................... // calculate hadronnes for method of RANDOM FOREST - TString hadRFName = "HadRF"; + MRanForestCalc rfcalc; rfcalc.SetHadronnessName(hadRFName); @@ -1360,12 +1381,5 @@ //----------------------------------------------------------------- - // 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); @@ -1482,9 +1496,10 @@ 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 = 4; + Int_t degree = 2; Double_t significance = -99.0; Bool_t drawpoly = kTRUE; @@ -1563,5 +1578,5 @@ TString parSCfile = outPath; - parSCfile += "parSC_2209a"; + parSCfile += "parSC_2310a"; gLog << "parSCfile = " << parSCfile << endl; @@ -1570,6 +1585,6 @@ // file to be updated (either ON or MC) - TString typeInput = "ON"; - //TString typeInput = "MC"; + //TString typeInput = "ON"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -1612,6 +1627,6 @@ //-------------------------------------------- - // files to contain the matrices (genertated from filenameTrain and - // filenameTest) + // files to contain the matrices (generated from filenameTrain and + // filenameTest) // // for the training @@ -1642,5 +1657,5 @@ //-------------------------- // create matrices and write them onto files - if (!RMatrix && (WOptimize || RTest) ) + if (!RMatrix) { MH3 &mh3 = *(new MH3("MHillas.fSize")); @@ -1911,5 +1926,5 @@ //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE"); - /* + MWriteRootFile write(outNameImage, "RECREATE"); @@ -1927,5 +1942,5 @@ write.AddContainer("HadRF", "Events"); write.AddContainer(hadSCName, "Events"); - */ + //----------------------------------------------------------------- @@ -2001,5 +2016,5 @@ tliston.AddToList(&fillhadsc); - //tliston.AddToList(&write); + tliston.AddToList(&write); tliston.AddToList(&contfinalgh); @@ -2292,7 +2307,6 @@ // and definition of final selections - //TString XX("NN"); - TString XX("SC"); - //TString XX("RF"); + //TString XX("SC"); + TString XX("RF"); TString fhadronnessName("Had"); fhadronnessName += XX; @@ -2300,5 +2314,5 @@ // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.30; + Float_t maxhadronness = 0.233; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; @@ -2355,6 +2369,4 @@ contfinalgh.SetName("ContSelFinalgh"); - MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); - fillhadnn.SetName("HhadNN"); MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); fillhadsc.SetName("HhadSC"); @@ -2362,4 +2374,15 @@ 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); @@ -2391,5 +2414,6 @@ pliston.AddToList(&tliston); InitBinnings(&pliston); - + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); //***************************** @@ -2400,8 +2424,8 @@ tliston.AddToList(&contfinalgh); - tliston.AddToList(&fillhadnn); tliston.AddToList(&fillhadsc); tliston.AddToList(&fillhadrf); + tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); @@ -2434,5 +2458,4 @@ // - pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); @@ -2446,13 +2469,38 @@ //------------------------------------------- - // fit alpha distribution to get the number of excess events - // - - MHHillasSrc* hillasSrc = - (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); - TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot - MHOnSubtraction onsub; - onsub.Calc(alphaHist, kTRUE, 13.1); + 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); + //------------------------------------------- @@ -2488,6 +2536,7 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_XX, WEX = " + gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (OEEst ? "kTRUE" : "kFALSE") << ", " << (WEX ? "kTRUE" : "kFALSE") << endl; @@ -2507,5 +2556,4 @@ // and definition of final selections - //TString XX("NN"); TString XX("SC"); //TString XX("RF"); @@ -2515,5 +2563,5 @@ // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.40; + Float_t maxhadronness = 0.4; Float_t maxalpha = 20.0; Float_t maxdist = 10.0; @@ -2676,6 +2724,8 @@ TString fImgParName = "MNewImagePar"; - - //=========================================================== + + if (OEEst) + { + //=========================================================== // // Optimization of energy estimator @@ -2691,8 +2741,8 @@ binsE, binsTheta); gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl; - - - if (WEX) - { + } + + if (WEX) + { //----------------------------------------------------------- // @@ -2706,10 +2756,20 @@ TFile enparam(energyParName); enparam.ls(); - - //MMcEnergyEst mcest("MMcEnergyEst"); - //mcest.Read("MMcEnergyEst"); - MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst")); - + 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()); @@ -2720,13 +2780,4 @@ parB[i] = mcest.GetCoeff( i+parA.GetSize() ); - 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; - //************************************************************************* // @@ -2758,4 +2809,5 @@ //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + MWriteRootFile write(filenameDataout, "RECREATE"); @@ -2916,5 +2968,4 @@ Int_t maxevents = -1; - //Int_t maxevents = 10000; if ( !evtloop.Eventloop(maxevents) ) return; @@ -2952,5 +3003,5 @@ gLog << "before DeleteBinnings" << endl; - //DeleteBinnings(&pliston); + DeleteBinnings(&pliston); gLog << "before Robert's code" << endl; @@ -3036,5 +3087,5 @@ } - gLog << "Macro CT1Analysis : End of Job F_XX" << endl; + gLog << "Macro CT1Analysis : End of Job E_XX" << endl; gLog << "=======================================================" << endl; }