Changeset 2437 for trunk/MagicSoft/Mars/macros
- Timestamp:
- 10/28/03 07:33:49 (21 years ago)
- Location:
- trunk/MagicSoft/Mars/macros
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/macros/CT1Analysis.C
r2436 r2437 1081 1081 fillmatg.SetName("fillGammas"); 1082 1082 1083 MFillH fillmath("MatrixHadrons");1084 fillmath.SetFilter(&flisth);1085 fillmath.SetName("fillHadrons");1086 1083 1087 1084 … … 1142 1139 selectorh.SetName("selectHadrons"); 1143 1140 1141 MFillH fillmath("MatrixHadrons"); 1142 fillmath.SetFilter(&flisth); 1143 fillmath.SetName("fillHadrons"); 1144 1144 1145 1145 1146 // entries in MFilterList … … 1182 1183 1183 1184 1184 // write out training events1185 // write out matrices of training events 1185 1186 1186 1187 gLog << "" << endl; 1187 1188 gLog << "========================================================" << endl; 1188 gLog << "Write out training events" << endl;1189 gLog << "Write out matrices of training events" << endl; 1189 1190 1190 1191 … … 1242 1243 // 1243 1244 MEvtLoop evtlooptr; 1245 evtlooptr.SetName("ReadRFTrees"); 1244 1246 evtlooptr.SetParList(&plisttr); 1245 1247 if (!evtlooptr.Eventloop()) … … 1300 1302 // 1301 1303 MEvtLoop treeloop; 1304 treeloop.SetName("GrowRFTrees"); 1302 1305 treeloop.SetParList(&plist2); 1303 1306 -
trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C
r2255 r2437 1 1 2 2 #include "CT1EgyEst.C" 3 3 4 4 5 void InitBinnings(MParList *plist) … … 142 143 143 144 // path for input for Mars 144 TString inPath = "~magican/ct1test/wittek/marsoutput/option C/";145 TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/"; 145 146 146 147 // path for output from Mars 147 TString outPath = "~magican/ct1test/wittek/marsoutput/option C/";148 TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/"; 148 149 149 150 //----------------------------------------------- 150 151 151 TEnv env("macros/CT1env.rc");152 Bool_t printEnv = kFALSE;152 //TEnv env("macros/CT1env.rc"); 153 //Bool_t printEnv = kFALSE; 153 154 154 155 //************************************************************************ … … 167 168 // (or OFF1.root or MC1.root)? 168 169 169 // Job B_NN_UP : read ON1.root (OFF1.root or MC1.root) file170 // - depending on RlookNN : create (or read in) hadron and gamma matrix171 // - calculate hadroness for method of NEAREST NEIGHBORS172 // - update the input files with the hadroness (ON1.root, OFF1.root173 // or MC1.root)174 175 Bool_t JobB_NN_UP = kFALSE;176 Bool_t RLookNN = kFALSE; // read in look-alike events177 Bool_t WNN = kFALSE; // update input root file ?178 179 170 180 171 // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file 181 // - depending on RLook : create (or read in) hadron and gamma matrix 182 // - depending on RTree : create (or read in) trees 172 // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas 173 // - if CTrainRF = TRUE : create matrices of training events 174 // - if RTree = TRUE : read in trees, otherwise create trees 183 175 // - calculate hadroness for method of RANDOM FOREST 184 // and for the SUPERCUTS 185 // - update the input files with the hadronesses (ON1.root, OFF1.root 186 // or MC1.root) 176 // - update the input files with the hadronesses (ON2.root, OFF2.root 177 // or MC2.root) 187 178 188 179 Bool_t JobB_RF_UP = kFALSE; 189 Bool_t RLookRF = kFALSE; // read in look-alike events 180 Bool_t RTrainRF = kFALSE; // read in matrices of training events 181 Bool_t CTrainRF = kFALSE; // create matrices of training events 190 182 Bool_t RTree = kFALSE; // read in trees 191 183 Bool_t WRF = kFALSE; // update input root file ? … … 265 257 gLog << "" << endl; 266 258 gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = " 267 << JobA << ", " << WPad << ", " << RPad << ", " << Wout << endl; 268 259 << (JobA ? "kTRUE" : "kFALSE") << ", " 260 << (WPad ? "kTRUE" : "kFALSE") << ", " 261 << (RPad ? "kTRUE" : "kFALSE") << ", " 262 << (Wout ? "kTRUE" : "kFALSE") << endl; 263 269 264 270 265 //-------------------------------------------------- … … 334 329 MCT1ReadPreProc readON(fileON); 335 330 336 MFCT1SelBasic selthetaon;337 selthetaon.SetCuts(-100.0, 29.5, 35.5);338 MContinue contthetaon(&selthetaon);331 //MFCT1SelBasic selthetaon; 332 //selthetaon.SetCuts(-100.0, 29.5, 35.5); 333 //MContinue contthetaon(&selthetaon); 339 334 340 335 MBlindPixelCalc blindon; … … 346 341 MHBlindPixels blindON("BlindPixelsON"); 347 342 MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels"); 343 fillblindON.SetName("FillBlind"); 348 344 349 345 MSigmabarCalc sigbarcalcon; … … 351 347 MHSigmaTheta sigthON("SigmaThetaON"); 352 348 MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); 353 349 fillsigthetaON.SetName("FillSigTheta"); 350 354 351 //***************************** 355 352 // entries in MParList … … 380 377 381 378 Int_t maxeventson = -1; 382 //Int_t maxeventson = 20000;379 //Int_t maxeventson = 10000; 383 380 if ( !evtloopon.Eventloop(maxeventson) ) 384 381 return; … … 420 417 421 418 MHBlindPixels blindOFF("BlindPixelsOFF"); 422 MHBlindPixels *pblindOFF = &blindOFF;423 gLog << "pblindOFF = " << pblindOFF << endl;424 425 419 MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels"); 420 fillblindOFF.SetName("FillBlindOFF"); 426 421 427 422 MSigmabarCalc sigbarcalcoff; … … 429 424 MHSigmaTheta sigthOFF("SigmaThetaOFF"); 430 425 MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); 431 426 fillsigthetaOFF.SetName("FillSigThetaOFF"); 427 432 428 //***************************** 433 429 // entries in MParList … … 657 653 MEvtLoop evtloop; 658 654 evtloop.SetParList(&pliston); 659 evtloop.ReadEnv(env, "", printEnv);655 //evtloop.ReadEnv(env, "", printEnv); 660 656 evtloop.SetProgressBar(&bar); 661 657 // evtloop.Write(); … … 694 690 695 691 692 693 696 694 //--------------------------------------------------------------------- 697 // Job B_ NN_UP695 // Job B_RF_UP 698 696 //============ 699 697 700 // - read in the matrices of look-alike events for gammas and hadrons 701 702 // then read ON1.root (or MC1.root) file703 // 704 // - calculate the hadroness for the method of nearest neighbors705 // 706 // - update input root file , includingthe hadroness707 708 709 if (JobB_ NN_UP)698 699 // - create (or read in) the matrices of training events for gammas 700 // and hadrons 701 // - create (or read in) the trees 702 // - then read ON1.root (or MC1.root) file 703 // - calculate the hadroness for the method of RANDOM FOREST 704 // - update input root file with the hadroness 705 706 707 if (JobB_RF_UP) 710 708 { 711 709 gLog << "=====================================================" << endl; 712 gLog << "Macro CT1Analysis : Start of Job B_ NN_UP" << endl;710 gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; 713 711 714 712 gLog << "" << endl; 715 gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " 716 << JobB_NN_UP << ", " << RLookNN << ", " << WNN << endl; 717 713 gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = " 714 << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " 715 << (RTrainRF ? "kTRUE" : "kFALSE") << ", " 716 << (CTrainRF ? "kTRUE" : "kFALSE") << ", " 717 << (RTree ? "kTRUE" : "kFALSE") << ", " 718 << (WRF ? "kTRUE" : "kFALSE") << endl; 719 720 721 //-------------------------------------------- 722 // parameters for the random forest 723 Int_t NumTrees = 100; 724 Int_t NumTry = 3; 725 Int_t NdSize = 1; 726 727 728 TString hadRFName = "HadRF"; 729 Float_t maxhadronness = 0.23; 730 Float_t maxalpha = 20.0; 731 Float_t maxdist = 10.0; 732 733 TString fHilName = "MHillas"; 734 TString fHilNameExt = "MHillasExt"; 735 TString fHilNameSrc = "MHillasSrc"; 736 TString fImgParName = "MNewImagePar"; 718 737 719 738 … … 721 740 // file to be updated (ON, OFF or MC) 722 741 723 //TString typeInput = "ON";742 TString typeInput = "ON"; 724 743 //TString typeInput = "OFF"; 725 TString typeInput = "MC";744 //TString typeInput = "MC"; 726 745 gLog << "typeInput = " << typeInput << endl; 727 746 … … 736 755 outNameImage += typeInput; 737 756 outNameImage += "2.root"; 738 739 757 //TString outNameImage = filenameData; 740 758 … … 742 760 743 761 //-------------------------------------------- 744 // files to be read for generating the look-alikeevents762 // files to be read for generating the matrices of training events 745 763 // "hadrons" : 746 764 TString filenameHad = outPath; 747 765 filenameHad += "OFF"; 748 766 filenameHad += "1.root"; 749 Int_t howManyHadrons = 500000;767 Int_t howManyHadrons = 1000000; 750 768 gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " 751 769 << howManyHadrons << endl; … … 761 779 762 780 //-------------------------------------------- 763 // files of look-alikeevents781 // files of training events 764 782 765 783 TString outNameGammas = outPath; 766 outNameGammas += " matrix_gammas_";784 outNameGammas += "RFmatrix_gammas_"; 767 785 outNameGammas += "MC"; 768 786 outNameGammas += ".root"; … … 772 790 773 791 TString outNameHadrons = outPath; 774 outNameHadrons += " matrix_hadrons_";792 outNameHadrons += "RFmatrix_hadrons_"; 775 793 outNameHadrons += typeMatrixHadrons; 776 794 outNameHadrons += ".root"; … … 778 796 779 797 MHMatrix matrixg("MatrixGammas"); 798 matrixg.EnableGraphicalOutput(); 799 800 matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 801 matrixg.AddColumn("MSigmabar.fSigmabar"); 802 matrixg.AddColumn("log10(MHillas.fSize)"); 803 matrixg.AddColumn("MHillasSrc.fDist"); 804 matrixg.AddColumn("MHillas.fWidth"); 805 matrixg.AddColumn("MHillas.fLength"); 806 matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); 807 matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); 808 matrixg.AddColumn("MNewImagePar.fConc"); 809 matrixg.AddColumn("MNewImagePar.fLeakage1"); 810 780 811 MHMatrix matrixh("MatrixHadrons"); 812 matrixh.EnableGraphicalOutput(); 813 814 matrixh.AddColumns(matrixg.GetColumns()); 815 816 //-------------------------------------------- 817 // file of trees of the random forest 818 819 TString outRF = outPath; 820 outRF += "RF.root"; 821 781 822 782 823 //************************************************************************* 783 // read in matrices of look-alikeevents784 if (R LookNN)824 // read in matrices of training events 825 if (RTrainRF) 785 826 { 786 827 const char* mtxName = "MatrixGammas"; … … 799 840 800 841 matrixg.Read(mtxName); 801 matrixg.Print(" cols");842 matrixg.Print("SizeCols"); 802 843 803 844 … … 819 860 820 861 matrixh.Read(mtxName); 821 matrixh.Print(" cols");862 matrixh.Print("SizeCols"); 822 863 } 823 864 824 865 825 866 //************************************************************************* 826 // create matrices of look-alikeevents827 if ( !RLookNN)867 // create matrices of training events 868 if (CTrainRF) 828 869 { 829 870 gLog << "" << endl; 830 871 gLog << "========================================================" << endl; 831 gLog << " Create matrices of look-alikeevents" << endl;872 gLog << " Create matrices of training events" << endl; 832 873 gLog << " Gammas :" << endl; 833 874 … … 852 893 fhadrons.SetName("hadronID)"); 853 894 854 MFEventSelector selectorg; 855 selectorg.SetNumSelectEvts(howManyGammas); 895 TString mgname("costhg"); 896 MBinning bing("Binning"+mgname); 897 bing.SetEdges(10, 0., 1.0); 898 899 MH3 gref("cos(MMcEvt.fTelescopeTheta)"); 900 gref.SetName(mgname); 901 MH::SetBinning(&gref.GetHist(), &bing); 902 for (Int_t i=1; i<=gref.GetNbins(); i++) 903 gref.GetHist().SetBinContent(i, 1.0); 904 905 MFEventSelector2 selectorg(gref); 906 selectorg.SetNumMax(howManyGammas); 856 907 selectorg.SetName("selectGammas"); 857 MFEventSelector selectorh;858 selectorh.SetNumSelectEvts(howManyHadrons);859 selectorh.SetName("selectHadrons");860 908 861 909 MFillH fillmatg("MatrixGammas"); 862 910 fillmatg.SetFilter(&flistg); 863 911 fillmatg.SetName("fillGammas"); 864 865 MFillH fillmath("MatrixHadrons");866 fillmath.SetFilter(&flisth);867 fillmath.SetName("fillHadrons");868 912 869 913 … … 894 938 MEvtLoop evtloopg; 895 939 evtloopg.SetParList(&plistg); 896 evtloopg.ReadEnv(env, "", printEnv);940 //evtloopg.ReadEnv(env, "", printEnv); 897 941 evtloopg.SetProgressBar(&matrixbar); 898 942 … … 907 951 908 952 gLog << " Hadrons :" << endl; 953 954 TString mhname("costhh"); 955 MBinning binh("Binning"+mhname); 956 binh.SetEdges(10, 0., 1.0); 957 958 MH3 href("cos(MMcEvt.fTelescopeTheta)"); 959 href.SetName(mhname); 960 MH::SetBinning(&href.GetHist(), &binh); 961 for (Int_t i=1; i<=href.GetNbins(); i++) 962 href.GetHist().SetBinContent(i, 1.0); 963 964 MFEventSelector2 selectorh(href); 965 //selectorh.SetNumMax(howManyHadrons); 966 // select as many hadrons as gammas 967 selectorh.SetNumMax(matrixg.GetM().GetNrows()); 968 selectorh.SetName("selectHadrons"); 969 970 MFillH fillmath("MatrixHadrons"); 971 fillmath.SetFilter(&flisth); 972 fillmath.SetName("fillHadrons"); 973 974 909 975 // entries in MFilterList 910 976 … … 932 998 MEvtLoop evtlooph; 933 999 evtlooph.SetParList(&plisth); 934 evtlooph.ReadEnv(env, "", printEnv);1000 //evtlooph.ReadEnv(env, "", printEnv); 935 1001 evtlooph.SetProgressBar(&matrixbar); 936 1002 … … 942 1008 //***************************************************** 943 1009 944 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 945 // 946 // ---------------------------------------------------------- 947 // Definition of the reference sample of look-alike events 948 // (this is a subsample of the original sample) 949 // ---------------------------------------------------------- 950 // 1010 1011 // write out matrices of training events events 1012 951 1013 gLog << "" << endl; 952 1014 gLog << "========================================================" << endl; 953 // Select a maximum of nmaxevts events from the sample of look-alike 954 // events. They will form the reference sample. 955 Int_t nmaxevts = 2000; 956 957 // target distribution for the variable in column refcolumn (start at 0); 958 //Int_t refcolumn = 0; 959 //Int_t nbins = 5; 960 //Float_t frombin = 0.5; 961 //Float_t tobin = 1.0; 962 //TH1F *thsh = new TH1F("thsh","target distribution", 963 // nbins, frombin, tobin); 964 //thsh->SetDirectory(NULL); 965 //thsh->SetXTitle("cos( \\Theta)"); 966 //thsh->SetYTitle("Counts"); 967 //Float_t dbin = (tobin-frombin)/nbins; 968 //Float_t lbin = frombin +dbin*0.5; 969 //for (Int_t j=1; j<=nbins; j++) 970 //{ 971 // thsh->Fill(lbin,1.0); 972 // lbin += dbin; 973 //} 974 975 Int_t refcolumn = 0; 976 MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); 977 TH1F *thsh = new TH1F(); 978 thsh->SetNameTitle("thsh","target distribution"); 979 MH::SetBinning(thsh, binscostheta); 980 Int_t nbins = thsh->GetNbinsX(); 981 Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; 982 for (Int_t j=1; j<=nbins; j++) 983 { 984 thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); 985 } 986 987 gLog << "" << endl; 988 gLog << "========================================================" << endl; 989 gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; 990 gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " 991 << refcolumn << ", " << nmaxevts << endl; 992 993 matrixg.EnableGraphicalOutput(); 994 Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); 995 996 if ( !matrixok ) 997 { 998 gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; 999 return; 1000 } 1001 1002 gLog << "" << endl; 1003 gLog << "========================================================" << endl; 1004 gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; 1005 gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " 1006 << refcolumn << ", " << nmaxevts << endl; 1007 1008 matrixh.EnableGraphicalOutput(); 1009 matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); 1010 delete thsh; 1011 1012 if ( !matrixok ) 1013 { 1014 gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; 1015 return; 1016 } 1017 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1018 1019 // write out look-alike events 1020 1021 gLog << "" << endl; 1022 gLog << "========================================================" << endl; 1023 gLog << "Write out look=alike events" << endl; 1015 gLog << "Write out matrices of training events" << endl; 1024 1016 1025 1017 … … 1027 1019 // "gammas" 1028 1020 gLog << "Gammas :" << endl; 1029 matrixg.Print(" cols");1021 matrixg.Print("SizeCols"); 1030 1022 1031 1023 TFile writeg(outNameGammas, "RECREATE", ""); … … 1033 1025 1034 1026 gLog << "" << endl; 1035 gLog << "Macro CT1Analysis : matrix of look-alikeevents for gammas written onto file "1027 gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file " 1036 1028 << outNameGammas << endl; 1037 1029 … … 1039 1031 // "hadrons" 1040 1032 gLog << "Hadrons :" << endl; 1041 matrixh.Print(" cols");1033 matrixh.Print("SizeCols"); 1042 1034 1043 1035 TFile writeh(outNameHadrons, "RECREATE", ""); … … 1045 1037 1046 1038 gLog << "" << endl; 1047 gLog << "Macro CT1Analysis : matrix of look-alikeevents for hadrons written onto file "1039 gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file " 1048 1040 << outNameHadrons << endl; 1049 1041 1050 1042 } 1051 //********** end of creating matrices of look-alike events *********** 1052 1053 1043 //********** end of creating matrices of training events *********** 1044 1045 1046 MRanForest *fRanForest; 1047 MRanTree *fRanTree; 1054 1048 //----------------------------------------------------------------- 1055 // Update the input files with the NN and SC hadronness 1056 // 1057 1058 if (WNN) 1049 // read in the trees of the random forest 1050 if (RTree) 1051 { 1052 MParList plisttr; 1053 MTaskList tlisttr; 1054 plisttr.AddToList(&tlisttr); 1055 1056 MReadTree readtr("TREE", outRF); 1057 readtr.DisableAutoScheme(); 1058 1059 MRanForestFill rffill; 1060 rffill.SetNumTrees(NumTrees); 1061 1062 // list of tasks for the loop over the trees 1063 1064 tlisttr.AddToList(&readtr); 1065 tlisttr.AddToList(&rffill); 1066 1067 //------------------- 1068 // Execute tree loop 1069 // 1070 MEvtLoop evtlooptr; 1071 evtlooptr.SetName("ReadRFTrees"); 1072 evtlooptr.SetParList(&plisttr); 1073 if (!evtlooptr.Eventloop()) 1074 return; 1075 1076 tlisttr.PrintStatistics(0, kTRUE); 1077 1078 1079 // get adresses of objects which are used in the next eventloop 1080 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); 1081 if (!fRanForest) 1082 { 1083 gLog << err << dbginf << "MRanForest not found... aborting." << endl; 1084 return kFALSE; 1085 } 1086 1087 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); 1088 if (!fRanTree) 1089 { 1090 gLog << err << dbginf << "MRanTree not found... aborting." << endl; 1091 return kFALSE; 1092 } 1093 1094 } 1095 1096 //----------------------------------------------------------------- 1097 // grow the trees of the random forest (event loop = tree loop) 1098 1099 if (!RTree) 1100 { 1101 1102 gLog << "" << endl; 1103 gLog << "========================================================" << endl; 1104 gLog << "Macro CT1Analysis : start growing trees" << endl; 1105 1106 MTaskList tlist2; 1107 MParList plist2; 1108 plist2.AddToList(&tlist2); 1109 1110 plist2.AddToList(&matrixg); 1111 plist2.AddToList(&matrixh); 1112 1113 MRanForestGrow rfgrow2; 1114 rfgrow2.SetNumTrees(NumTrees); 1115 rfgrow2.SetNumTry(NumTry); 1116 rfgrow2.SetNdSize(NdSize); 1117 1118 MWriteRootFile rfwrite2(outRF); 1119 rfwrite2.AddContainer("MRanTree", "TREE"); 1120 1121 // list of tasks for the loop over the trees 1122 1123 tlist2.AddToList(&rfgrow2); 1124 tlist2.AddToList(&rfwrite2); 1125 1126 //------------------- 1127 // Execute tree loop 1128 // 1129 MEvtLoop treeloop; 1130 treeloop.SetName("GrowRFTrees"); 1131 treeloop.SetParList(&plist2); 1132 1133 if ( !treeloop.Eventloop() ) 1134 return; 1135 1136 tlist2.PrintStatistics(0, kTRUE); 1137 1138 // get adresses of objects which are used in the next eventloop 1139 fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); 1140 if (!fRanForest) 1141 { 1142 gLog << err << dbginf << "MRanForest not found... aborting." << endl; 1143 return kFALSE; 1144 } 1145 1146 fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); 1147 if (!fRanTree) 1148 { 1149 gLog << err << dbginf << "MRanTree not found... aborting." << endl; 1150 return kFALSE; 1151 } 1152 1153 } 1154 // end of growing the trees of the random forest 1155 //----------------------------------------------------------------- 1156 1157 1158 1159 1160 //----------------------------------------------------------------- 1161 // Update the input files with the RF hadronness 1162 // 1163 1164 if (WRF) 1059 1165 { 1060 1166 gLog << "" << endl; 1061 1167 gLog << "========================================================" << endl; 1062 1168 gLog << "Update input file '" << filenameData 1063 << "' with the NN and SChadronness" << endl;1169 << "' with the RF hadronness" << endl; 1064 1170 1065 1171 MTaskList tliston; … … 1078 1184 read.DisableAutoScheme(); 1079 1185 1080 TString fHilName = "MHillas";1081 TString fHilNameExt = "MHillasExt";1082 TString fHilNameSrc = "MHillasSrc";1083 TString fImgParName = "MNewImagePar";1084 1186 1085 1187 //....................................................................... 1086 // calculation of hadroness for method of Nearest Neighbors 1087 1088 TString hadNNName = "HadNN"; 1089 MMultiDimDistCalc nncalc; 1090 nncalc.SetUseNumRows(25); 1091 nncalc.SetUseKernelMethod(kFALSE); 1092 nncalc.SetHadronnessName(hadNNName); 1188 // calculate hadronnes for method of RANDOM FOREST 1189 1190 1191 MRanForestCalc rfcalc; 1192 rfcalc.SetHadronnessName(hadRFName); 1193 1093 1194 1094 1195 //....................................................................... … … 1108 1209 write.AddContainer("MNewImagePar", "Events"); 1109 1210 1110 write.AddContainer("HadRF", "Events"); 1111 write.AddContainer("HadSC", "Events"); 1112 write.AddContainer("HadNN", "Events"); 1113 1211 write.AddContainer(hadRFName, "Events"); 1114 1212 1115 1213 //----------------------------------------------------------------- … … 1118 1216 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1119 1217 1120 Float_t maxhadronness = 0.40; 1121 Float_t maxalpha = 20.0; 1122 Float_t maxdist = 10.0; 1218 1123 1219 1124 1220 MFCT1SelFinal selfinalgh(fHilNameSrc); 1125 1221 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); 1126 selfinalgh.SetHadronnessName(had NNName);1222 selfinalgh.SetHadronnessName(hadRFName); 1127 1223 selfinalgh.SetName("SelFinalgh"); 1128 1224 MContinue contfinalgh(&selfinalgh); 1129 1225 contfinalgh.SetName("ContSelFinalgh"); 1130 1226 1131 MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); 1132 fillhadnn.SetName("HhadNN"); 1227 MFillH fillranfor("MHRanForest"); 1228 fillranfor.SetName("HRanForest"); 1229 1230 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); 1231 fillhadrf.SetName("HhadRF"); 1133 1232 1134 1233 MFCT1SelFinal selfinal(fHilNameSrc); 1135 1234 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 1136 selfinal.SetHadronnessName(had NNName);1235 selfinal.SetHadronnessName(hadRFName); 1137 1236 selfinal.SetName("SelFinal"); 1138 1237 MContinue contfinal(&selfinal); 1139 1238 contfinal.SetName("ContSelFinal"); 1239 1240 TString mh3name = "abs(Alpha)"; 1241 MBinning binsalphaabs("Binning"+mh3name); 1242 binsalphaabs.SetEdges(50, -2.0, 98.0); 1243 1244 MH3 alphaabs("abs(MHillasSrc.fAlpha)"); 1245 alphaabs.SetName(mh3name); 1246 MFillH alpha(&alphaabs); 1247 alpha.SetName("FillAlphaAbs"); 1140 1248 1141 1249 … … 1161 1269 InitBinnings(&pliston); 1162 1270 1163 pliston.AddToList(&matrixg); 1164 pliston.AddToList(&matrixh); 1271 pliston.AddToList(fRanForest); 1272 pliston.AddToList(fRanTree); 1273 1274 pliston.AddToList(&binsalphaabs); 1275 pliston.AddToList(&alphaabs); 1165 1276 1166 1277 … … 1170 1281 tliston.AddToList(&read); 1171 1282 1172 tliston.AddToList(&nncalc); 1173 tliston.AddToList(&fillhadnn); 1174 1283 tliston.AddToList(&rfcalc); 1284 tliston.AddToList(&fillranfor); 1285 tliston.AddToList(&fillhadrf); 1286 1287 tliston.AddToList(&write); 1288 tliston.AddToList(&contfinalgh); 1289 1290 tliston.AddToList(&alpha); 1175 1291 tliston.AddToList(&hfill1); 1176 1292 tliston.AddToList(&hfill2); … … 1179 1295 tliston.AddToList(&hfill5); 1180 1296 1181 tliston.AddToList(&write);1182 1183 tliston.AddToList(&contfinalgh);1184 1297 tliston.AddToList(&contfinal); 1185 1298 … … 1200 1313 tliston.PrintStatistics(0, kTRUE); 1201 1314 1315 1202 1316 //------------------------------------------- 1203 1317 // Display the histograms 1204 1318 // 1205 pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); 1319 pliston.FindObject("MHRanForest")->DrawClone(); 1320 pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); 1206 1321 1207 1322 pliston.FindObject("MHHillas")->DrawClone(); … … 1211 1326 pliston.FindObject("MHStarMap")->DrawClone(); 1212 1327 1328 1329 //------------------------------------------- 1330 // fit alpha distribution to get the number of excess events and 1331 // calculate significance of gamma signal in the alpha plot 1332 1333 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); 1334 TH1 &alphaHist = absalpha->GetHist(); 1335 alphaHist.SetXTitle("|alpha| [\\circ]"); 1336 alphaHist.SetName("alpha-macro"); 1337 1338 Double_t alphasig = 13.1; 1339 Double_t alphamin = 30.0; 1340 Double_t alphamax = 90.0; 1341 Int_t degree = 2; 1342 Double_t significance = -99.0; 1343 Bool_t drawpoly = kTRUE; 1344 Bool_t fitgauss = kTRUE; 1345 Bool_t print = kTRUE; 1346 1347 MHFindSignificance findsig; 1348 findsig.SetRebin(kTRUE); 1349 findsig.SetReduceDegree(kFALSE); 1350 1351 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 1352 alphasig, drawpoly, fitgauss, print); 1353 significance = findsig.GetSignificance(); 1354 Float_t alphasi = findsig.GetAlphasi(); 1355 1356 gLog << "For file '" << filenameData << "' : " << endl; 1357 gLog << "Significance of gamma signal after supercuts : " 1358 << significance << " (for |alpha| < " << alphasi << " degrees)" 1359 << endl; 1360 1361 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); 1362 1363 //------------------------------------------- 1364 1365 1213 1366 DeleteBinnings(&pliston); 1214 1367 } 1215 1368 1216 1217 1218 gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl; 1369 gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; 1219 1370 gLog << "=======================================================" << endl; 1220 1371 } … … 1222 1373 1223 1374 1375 1224 1376 //--------------------------------------------------------------------- 1225 // Job B_RF_UP 1226 //============ 1227 1228 1229 // - create (or read in) the matrices of look-alike events for gammas 1230 // and hadrons 1231 // - create (or read in) the trees 1232 // - then read ON1.root (or MC1.root) file 1233 // - calculate the hadroness for the method of RANDOM FOREST 1234 // - update input root file with the hadroness 1235 1236 1237 if (JobB_RF_UP) 1377 // Job C 1378 //====== 1379 1380 // - read ON1 and MC1 data files 1381 // which should have been updated to contain the hadronnesses 1382 // for the method of NEAREST NEIGHBORS and for the SUOERCUTS 1383 // - produce Neyman-Pearson plots 1384 1385 if (JobC) 1238 1386 { 1239 1387 gLog << "=====================================================" << endl; 1240 gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;1388 gLog << "Macro CT1Analysis : Start of Job C" << endl; 1241 1389 1242 1390 gLog << "" << endl; 1243 gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " 1244 << JobB_RF_UP << ", " << RLookRF << ", " << RTree << ", " 1245 << WRF << endl; 1246 1247 1248 1249 //-------------------------------------------- 1250 // parameters for the random forest 1251 Int_t NumTrees = 100; 1252 Int_t NumTry = 3; 1253 Int_t NdSize = 3; 1254 1255 1256 //-------------------------------------------- 1257 // file to be updated (ON, OFF or MC) 1258 1259 TString typeInput = "ON"; 1260 //TString typeInput = "OFF"; 1261 //TString typeInput = "MC"; 1262 gLog << "typeInput = " << typeInput << endl; 1263 1264 // name of input root file 1391 gLog << "Macro CT1Analysis : JobC = " << JobC << endl; 1392 1393 1394 // name of input data file 1265 1395 TString filenameData = outPath; 1266 filenameData += typeInput; 1267 filenameData += "1.root"; 1268 gLog << "filenameData = " << filenameData << endl; 1269 1270 // name of output root file 1271 TString outNameImage = outPath; 1272 outNameImage += typeInput; 1273 outNameImage += "2.root"; 1274 //TString outNameImage = filenameData; 1275 1276 gLog << "outNameImage = " << outNameImage << endl; 1277 1278 //-------------------------------------------- 1279 // files to be read for generating the look-alike events 1280 // "hadrons" : 1281 TString filenameHad = outPath; 1282 filenameHad += "OFF"; 1283 filenameHad += "1.root"; 1284 Int_t howManyHadrons = 1000000; 1285 gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " 1286 << howManyHadrons << endl; 1287 1288 1289 // "gammas" : 1396 filenameData += "OFF"; 1397 filenameData += "2.root"; 1398 gLog << "filenameData = " << filenameData << endl; 1399 1400 // name of input MC file 1290 1401 TString filenameMC = outPath; 1291 1402 filenameMC += "MC"; 1292 filenameMC += "1.root"; 1293 Int_t howManyGammas = 50000; 1294 gLog << "filenameMC = " << filenameMC << ", howManyGammas = " 1295 << howManyGammas << endl; 1296 1297 //-------------------------------------------- 1298 // files of look-alike events 1299 1300 TString outNameGammas = outPath; 1301 outNameGammas += "RFmatrix_gammas_"; 1302 outNameGammas += "MC"; 1303 outNameGammas += ".root"; 1304 1305 TString typeMatrixHadrons = "OFF"; 1306 gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; 1307 1308 TString outNameHadrons = outPath; 1309 outNameHadrons += "RFmatrix_hadrons_"; 1310 outNameHadrons += typeMatrixHadrons; 1311 outNameHadrons += ".root"; 1312 1313 1314 MHMatrix matrixg("MatrixGammas"); 1315 MHMatrix matrixh("MatrixHadrons"); 1316 1317 //-------------------------------------------- 1318 // file of trees of the random forest 1319 1320 TString outRF = outPath; 1321 outRF += "RF.root"; 1322 1323 1324 //************************************************************************* 1325 // read in matrices of look-alike events 1326 if (RLookRF) 1327 { 1328 const char* mtxName = "MatrixGammas"; 1329 1330 gLog << "" << endl; 1331 gLog << "========================================================" << endl; 1332 gLog << "Get matrix for (gammas)" << endl; 1333 gLog << "matrix name = " << mtxName << endl; 1334 gLog << "name of root file = " << outNameGammas << endl; 1335 gLog << "" << endl; 1336 1337 1338 // read in the object with the name 'mtxName' from file 'outNameGammas' 1339 // 1340 TFile fileg(outNameGammas); 1341 1342 matrixg.Read(mtxName); 1343 matrixg.Print("cols"); 1344 1345 1346 //***************************************************************** 1347 1348 const char* mtxName = "MatrixHadrons"; 1349 1350 gLog << "" << endl; 1351 gLog << "========================================================" << endl; 1352 gLog << " Get matrix for (hadrons)" << endl; 1353 gLog << "matrix name = " << mtxName << endl; 1354 gLog << "name of root file = " << outNameHadrons << endl; 1355 gLog << "" << endl; 1356 1357 1358 // read in the object with the name 'mtxName' from file 'outNameHadrons' 1359 // 1360 TFile fileh(outNameHadrons); 1361 1362 matrixh.Read(mtxName); 1363 matrixh.Print("cols"); 1364 } 1365 1366 1367 //************************************************************************* 1368 // create matrices of look-alike events 1369 if (!RLookRF) 1370 { 1371 gLog << "" << endl; 1372 gLog << "========================================================" << endl; 1373 gLog << " Create matrices of look-alike events" << endl; 1374 gLog << " Gammas :" << endl; 1375 1376 1377 MParList plistg; 1378 MTaskList tlistg; 1379 MFilterList flistg; 1380 1381 MParList plisth; 1382 MTaskList tlisth; 1383 MFilterList flisth; 1384 1385 MReadMarsFile readg("Events", filenameMC); 1386 readg.DisableAutoScheme(); 1387 1388 MReadMarsFile readh("Events", filenameHad); 1389 readh.DisableAutoScheme(); 1390 1391 MFParticleId fgamma("MMcEvt", '=', kGAMMA); 1392 fgamma.SetName("gammaID"); 1393 MFParticleId fhadrons("MMcEvt", '!', kGAMMA); 1394 fhadrons.SetName("hadronID)"); 1395 1396 MFEventSelector selectorg; 1397 selectorg.SetNumSelectEvts(howManyGammas); 1398 selectorg.SetName("selectGammas"); 1399 MFEventSelector selectorh; 1400 selectorh.SetNumSelectEvts(howManyHadrons); 1401 selectorh.SetName("selectHadrons"); 1402 1403 MFillH fillmatg("MatrixGammas"); 1404 fillmatg.SetFilter(&flistg); 1405 fillmatg.SetName("fillGammas"); 1406 1407 MFillH fillmath("MatrixHadrons"); 1408 fillmath.SetFilter(&flisth); 1409 fillmath.SetName("fillHadrons"); 1410 1411 1412 //***************************** fill gammas *** 1413 // entries in MFilterList 1414 1415 flistg.AddToList(&fgamma); 1416 flistg.AddToList(&selectorg); 1417 1418 //***************************** 1419 // entries in MParList 1420 1421 plistg.AddToList(&tlistg); 1422 InitBinnings(&plistg); 1423 1424 plistg.AddToList(&matrixg); 1425 1426 //***************************** 1427 // entries in MTaskList 1428 1429 tlistg.AddToList(&readg); 1430 tlistg.AddToList(&flistg); 1431 tlistg.AddToList(&fillmatg); 1432 1433 //***************************** 1434 1435 MProgressBar matrixbar; 1436 MEvtLoop evtloopg; 1437 evtloopg.SetParList(&plistg); 1438 evtloopg.ReadEnv(env, "", printEnv); 1439 evtloopg.SetProgressBar(&matrixbar); 1440 1441 Int_t maxevents = -1; 1442 if (!evtloopg.Eventloop(maxevents)) 1443 return; 1444 1445 tlistg.PrintStatistics(0, kTRUE); 1446 1447 1448 //***************************** fill hadrons *** 1449 1450 gLog << " Hadrons :" << endl; 1451 // entries in MFilterList 1452 1453 flisth.AddToList(&fhadrons); 1454 flisth.AddToList(&selectorh); 1455 1456 //***************************** 1457 // entries in MParList 1458 1459 plisth.AddToList(&tlisth); 1460 InitBinnings(&plisth); 1461 1462 plisth.AddToList(&matrixh); 1463 1464 //***************************** 1465 // entries in MTaskList 1466 1467 tlisth.AddToList(&readh); 1468 tlisth.AddToList(&flisth); 1469 tlisth.AddToList(&fillmath); 1470 1471 //***************************** 1472 1473 MProgressBar matrixbar; 1474 MEvtLoop evtlooph; 1475 evtlooph.SetParList(&plisth); 1476 evtlooph.ReadEnv(env, "", printEnv); 1477 evtlooph.SetProgressBar(&matrixbar); 1478 1479 Int_t maxevents = -1; 1480 if (!evtlooph.Eventloop(maxevents)) 1481 return; 1482 1483 tlisth.PrintStatistics(0, kTRUE); 1484 //***************************************************** 1485 1486 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1487 // 1488 // ---------------------------------------------------------- 1489 // Definition of the reference sample of look-alike events 1490 // (this is a subsample of the original sample) 1491 // ---------------------------------------------------------- 1492 // 1493 gLog << "" << endl; 1494 gLog << "========================================================" << endl; 1495 // Select a maximum of nmaxevts events from the sample of look-alike 1496 // events. They will form the reference sample. 1497 Int_t nmaxevts = 10000; 1498 1499 // target distribution for the variable in column refcolumn (start at 0); 1500 //Int_t refcolumn = 0; 1501 //Int_t nbins = 5; 1502 //Float_t frombin = 0.5; 1503 //Float_t tobin = 1.0; 1504 //TH1F *thsh = new TH1F("thsh","target distribution", 1505 // nbins, frombin, tobin); 1506 //thsh->SetDirectory(NULL); 1507 //thsh->SetXTitle("cos( \\Theta)"); 1508 //thsh->SetYTitle("Counts"); 1509 //Float_t dbin = (tobin-frombin)/nbins; 1510 //Float_t lbin = frombin +dbin*0.5; 1511 //for (Int_t j=1; j<=nbins; j++) 1512 //{ 1513 // thsh->Fill(lbin,1.0); 1514 // lbin += dbin; 1515 //} 1516 1517 Int_t refcolumn = 0; 1518 MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); 1519 TH1F *thsh = new TH1F(); 1520 thsh->SetNameTitle("thsh","target distribution"); 1521 MH::SetBinning(thsh, binscostheta); 1522 Int_t nbins = thsh->GetNbinsX(); 1523 Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; 1524 for (Int_t j=1; j<=nbins; j++) 1525 { 1526 thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); 1527 } 1528 1529 gLog << "" << endl; 1530 gLog << "========================================================" << endl; 1531 gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; 1532 gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " 1533 << refcolumn << ", " << nmaxevts << endl; 1534 1535 matrixg.EnableGraphicalOutput(); 1536 Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); 1537 1538 if ( !matrixok ) 1539 { 1540 gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; 1541 return; 1542 } 1543 1544 gLog << "" << endl; 1545 gLog << "========================================================" << endl; 1546 gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; 1547 gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " 1548 << refcolumn << ", " << nmaxevts << endl; 1549 1550 matrixh.EnableGraphicalOutput(); 1551 matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); 1552 delete thsh; 1553 1554 if ( !matrixok ) 1555 { 1556 gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; 1557 return; 1558 } 1559 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1560 1561 // write out look-alike events 1562 1563 gLog << "" << endl; 1564 gLog << "========================================================" << endl; 1565 gLog << "Write out look=alike events" << endl; 1566 1567 1568 //------------------------------------------- 1569 // "gammas" 1570 gLog << "Gammas :" << endl; 1571 matrixg.Print("cols"); 1572 1573 TFile writeg(outNameGammas, "RECREATE", ""); 1574 matrixg.Write(); 1575 1576 gLog << "" << endl; 1577 gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " 1578 << outNameGammas << endl; 1579 1580 //------------------------------------------- 1581 // "hadrons" 1582 gLog << "Hadrons :" << endl; 1583 matrixh.Print("cols"); 1584 1585 TFile writeh(outNameHadrons, "RECREATE", ""); 1586 matrixh.Write(); 1587 1588 gLog << "" << endl; 1589 gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " 1590 << outNameHadrons << endl; 1591 1592 } 1593 //********** end of creating matrices of look-alike events *********** 1594 1595 1596 MRanForest *fRanForest; 1597 MRanTree *fRanTree; 1403 filenameMC += "2.root"; 1404 gLog << "filenameMC = " << filenameMC << endl; 1405 1406 1598 1407 //----------------------------------------------------------------- 1599 // read in the trees of the random forest1600 if (RTree)1601 {1602 MParList plisttr;1603 MTaskList tlisttr;1604 plisttr.AddToList(&tlisttr);1605 1606 MReadTree readtr("TREE", outRF);1607 readtr.DisableAutoScheme();1608 1609 MRanForestFill rffill;1610 rffill.SetNumTrees(NumTrees);1611 1612 // list of tasks for the loop over the trees1613 1614 tlisttr.AddToList(&readtr);1615 tlisttr.AddToList(&rffill);1616 1617 //-------------------1618 // Execute tree loop1619 //1620 MEvtLoop evtlooptr;1621 evtlooptr.SetParList(&plisttr);1622 if (!evtlooptr.Eventloop())1623 return;1624 1625 tlisttr.PrintStatistics(0, kTRUE);1626 1627 1628 // get adresses of objects which are used in the next eventloop1629 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");1630 if (!fRanForest)1631 {1632 gLog << err << dbginf << "MRanForest not found... aborting." << endl;1633 return kFALSE;1634 }1635 1636 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");1637 if (!fRanTree)1638 {1639 gLog << err << dbginf << "MRanTree not found... aborting." << endl;1640 return kFALSE;1641 }1642 1643 }1644 1645 //-----------------------------------------------------------------1646 // grow the trees of the random forest (event loop = tree loop)1647 1648 if (!RTree)1649 {1650 1651 gLog << "" << endl;1652 gLog << "========================================================" << endl;1653 gLog << "Macro CT1Analysis : start growing trees" << endl;1654 1655 MTaskList tlist2;1656 MParList plist2;1657 plist2.AddToList(&tlist2);1658 1659 plist2.AddToList(&matrixg);1660 plist2.AddToList(&matrixh);1661 1662 MRanForestGrow rfgrow2;1663 rfgrow2.SetNumTrees(NumTrees);1664 rfgrow2.SetNumTry(NumTry);1665 rfgrow2.SetNdSize(NdSize);1666 1667 MWriteRootFile rfwrite2(outRF);1668 rfwrite2.AddContainer("MRanTree", "TREE");1669 1670 // list of tasks for the loop over the trees1671 1672 tlist2.AddToList(&rfgrow2);1673 tlist2.AddToList(&rfwrite2);1674 1675 //-------------------1676 // Execute tree loop1677 //1678 MEvtLoop treeloop;1679 treeloop.SetParList(&plist2);1680 1681 if ( !treeloop.Eventloop() )1682 return;1683 1684 tlist2.PrintStatistics(0, kTRUE);1685 1686 // get adresses of objects which are used in the next eventloop1687 fRanForest = (MRanForest*)plist2->FindObject("MRanForest");1688 if (!fRanForest)1689 {1690 gLog << err << dbginf << "MRanForest not found... aborting." << endl;1691 return kFALSE;1692 }1693 1694 fRanTree = (MRanTree*)plist2->FindObject("MRanTree");1695 if (!fRanTree)1696 {1697 gLog << err << dbginf << "MRanTree not found... aborting." << endl;1698 return kFALSE;1699 }1700 1701 }1702 // end of growing the trees of the random forest1703 //-----------------------------------------------------------------1704 1705 1706 1707 1708 //-----------------------------------------------------------------1709 // Update the input files with the RF hadronness1710 //1711 1712 if (WRF)1713 {1714 gLog << "" << endl;1715 gLog << "========================================================" << endl;1716 gLog << "Update input file '" << filenameData1717 << "' with the RF hadronness" << endl;1718 1408 1719 1409 MTaskList tliston; … … 1729 1419 // 1730 1420 1731 MReadMarsFile read("Events", filenameData); 1421 MReadMarsFile read("Events", filenameMC); 1422 read.AddFile(filenameData); 1732 1423 read.DisableAutoScheme(); 1424 1425 1426 //....................................................................... 1427 // names of hadronness containers 1428 1429 TString hadNNName = "HadNN"; 1430 TString hadSCName = "HadSC"; 1431 TString hadRFName = "HadRF"; 1432 1433 //....................................................................... 1434 1733 1435 1734 1436 TString fHilName = "MHillas"; … … 1737 1439 TString fImgParName = "MNewImagePar"; 1738 1440 1441 Float_t maxhadronness = 0.40; 1442 Float_t maxalpha = 20.0; 1443 Float_t maxdist = 10.0; 1444 1445 MFCT1SelFinal selfinalgh(fHilNameSrc); 1446 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); 1447 selfinalgh.SetHadronnessName(hadRFName); 1448 selfinalgh.SetName("SelFinalgh"); 1449 MContinue contfinalgh(&selfinalgh); 1450 contfinalgh.SetName("ContSelFinalgh"); 1451 1452 //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); 1453 //fillhadnn.SetName("HhadNN"); 1454 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); 1455 fillhadsc.SetName("HhadSC"); 1456 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); 1457 fillhadrf.SetName("HhadRF"); 1458 1459 MFCT1SelFinal selfinal(fHilNameSrc); 1460 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 1461 selfinal.SetHadronnessName(hadRFName); 1462 selfinal.SetName("SelFinal"); 1463 MContinue contfinal(&selfinal); 1464 contfinal.SetName("ContSelFinal"); 1465 1466 1467 MFillH hfill1("MHHillas", fHilName); 1468 hfill1.SetName("HHillas"); 1469 1470 MFillH hfill2("MHStarMap", fHilName); 1471 hfill2.SetName("HStarMap"); 1472 1473 MFillH hfill3("MHHillasExt", fHilNameSrc); 1474 hfill3.SetName("HHillasExt"); 1475 1476 MFillH hfill4("MHHillasSrc", fHilNameSrc); 1477 hfill4.SetName("HHillasSrc"); 1478 1479 MFillH hfill5("MHNewImagePar", fImgParName); 1480 hfill5.SetName("HNewImagePar"); 1481 1482 1483 //***************************** 1484 // entries in MParList 1485 1486 pliston.AddToList(&tliston); 1487 InitBinnings(&pliston); 1488 1489 1490 //***************************** 1491 // entries in MTaskList 1492 1493 tliston.AddToList(&read); 1494 1495 //tliston.AddToList(&fillhadnn); 1496 tliston.AddToList(&fillhadsc); 1497 tliston.AddToList(&fillhadrf); 1498 1499 tliston.AddToList(&contfinalgh); 1500 tliston.AddToList(&hfill1); 1501 tliston.AddToList(&hfill2); 1502 tliston.AddToList(&hfill3); 1503 tliston.AddToList(&hfill4); 1504 tliston.AddToList(&hfill5); 1505 1506 tliston.AddToList(&contfinal); 1507 1508 //***************************** 1509 1510 //------------------------------------------- 1511 // Execute event loop 1512 // 1513 MProgressBar bar; 1514 MEvtLoop evtloop; 1515 evtloop.SetParList(&pliston); 1516 evtloop.SetProgressBar(&bar); 1517 1518 Int_t maxevents = -1; 1519 //Int_t maxevents = 35000; 1520 if ( !evtloop.Eventloop(maxevents) ) 1521 return; 1522 1523 tliston.PrintStatistics(0, kTRUE); 1524 1525 1526 //------------------------------------------- 1527 // Display the histograms 1528 // 1529 1530 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); 1531 pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); 1532 pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); 1533 1534 pliston.FindObject("MHHillas")->DrawClone(); 1535 pliston.FindObject("MHHillasExt")->DrawClone(); 1536 pliston.FindObject("MHHillasSrc")->DrawClone(); 1537 pliston.FindObject("MHNewImagePar")->DrawClone(); 1538 pliston.FindObject("MHStarMap")->DrawClone(); 1539 1540 DeleteBinnings(&pliston); 1541 1542 gLog << "Macro CT1Analysis : End of Job C" << endl; 1543 gLog << "===================================================" << endl; 1544 } 1545 1546 1547 //--------------------------------------------------------------------- 1548 // Job D 1549 //====== 1550 1551 // - select g/h separation method XX 1552 // - read ON2 (or MC2) root file 1553 // - apply cuts in hadronness 1554 // - make plots 1555 1556 1557 if (JobD) 1558 { 1559 gLog << "=====================================================" << endl; 1560 gLog << "Macro CT1Analysis : Start of Job D" << endl; 1561 1562 gLog << "" << endl; 1563 gLog << "Macro CT1Analysis : JobD = " 1564 << JobD << endl; 1565 1566 // type of data to be analysed 1567 TString typeData = "ON"; 1568 //TString typeData = "OFF"; 1569 //TString typeData = "MC"; 1570 gLog << "typeData = " << typeData << endl; 1571 1572 TString ext = "2.root"; 1573 1574 1575 //------------------------------ 1576 // selection of g/h separation method 1577 // and definition of final selections 1578 1579 //TString XX("NN"); 1580 //TString XX("SC"); 1581 TString XX("RF"); 1582 TString fhadronnessName("Had"); 1583 fhadronnessName += XX; 1584 gLog << "fhadronnessName = " << fhadronnessName << endl; 1585 1586 // maximum values of the hadronness, |ALPHA| and DIST 1587 Float_t maxhadronness = 0.30; 1588 Float_t maxalpha = 20.0; 1589 Float_t maxdist = 10.0; 1590 gLog << "Maximum values of hadronness, |ALPHA| and DIST = " 1591 << maxhadronness << ", " << maxalpha << ", " 1592 << maxdist << endl; 1593 1594 1595 //------------------------------ 1596 // name of data file to be analysed 1597 TString filenameData(outPath); 1598 filenameData += typeData; 1599 filenameData += ext; 1600 gLog << "filenameData = " << filenameData << endl; 1601 1602 1603 1604 //************************************************************************* 1605 // 1606 // Analyse the data 1607 // 1608 1609 MTaskList tliston; 1610 MParList pliston; 1611 1612 // geometry is needed in MHHillas... classes 1613 MGeomCam *fGeom = 1614 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1615 1616 1617 TString fHilName = "MHillas"; 1618 TString fHilNameExt = "MHillasExt"; 1619 TString fHilNameSrc = "MHillasSrc"; 1620 TString fImgParName = "MNewImagePar"; 1621 1622 //------------------------------------------- 1623 // create the tasks which should be executed 1624 // 1625 1626 MReadMarsFile read("Events", filenameData); 1627 read.DisableAutoScheme(); 1628 1629 1630 //----------------------------------------------------------------- 1631 // geometry is needed in MHHillas... classes 1632 MGeomCam *fGeom = 1633 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1634 1635 MFCT1SelFinal selfinalgh(fHilNameSrc); 1636 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); 1637 selfinalgh.SetHadronnessName(fhadronnessName); 1638 selfinalgh.SetName("SelFinalgh"); 1639 MContinue contfinalgh(&selfinalgh); 1640 contfinalgh.SetName("ContSelFinalgh"); 1641 1642 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); 1643 //fillhadnn.SetName("HhadNN"); 1644 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); 1645 fillhadsc.SetName("HhadSC"); 1646 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); 1647 fillhadrf.SetName("HhadRF"); 1648 1649 1650 MFillH hfill1("MHHillas", fHilName); 1651 hfill1.SetName("HHillas"); 1652 1653 MFillH hfill2("MHStarMap", fHilName); 1654 hfill2.SetName("HStarMap"); 1655 1656 MFillH hfill3("MHHillasExt", fHilNameSrc); 1657 hfill3.SetName("HHillasExt"); 1658 1659 MFillH hfill4("MHHillasSrc", fHilNameSrc); 1660 hfill4.SetName("HHillasSrc"); 1661 1662 MFillH hfill5("MHNewImagePar", fImgParName); 1663 hfill5.SetName("HNewImagePar"); 1664 1665 MFCT1SelFinal selfinal(fHilNameSrc); 1666 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 1667 selfinal.SetHadronnessName(fhadronnessName); 1668 selfinal.SetName("SelFinal"); 1669 MContinue contfinal(&selfinal); 1670 contfinal.SetName("ContSelFinal"); 1671 1672 1673 //***************************** 1674 // entries in MParList 1675 1676 pliston.AddToList(&tliston); 1677 InitBinnings(&pliston); 1678 1679 1680 //***************************** 1681 // entries in MTaskList 1682 1683 tliston.AddToList(&read); 1684 1685 tliston.AddToList(&contfinalgh); 1686 1687 //tliston.AddToList(&fillhadnn); 1688 tliston.AddToList(&fillhadsc); 1689 tliston.AddToList(&fillhadrf); 1690 1691 tliston.AddToList(&hfill1); 1692 tliston.AddToList(&hfill2); 1693 tliston.AddToList(&hfill3); 1694 tliston.AddToList(&hfill4); 1695 tliston.AddToList(&hfill5); 1696 1697 tliston.AddToList(&contfinal); 1698 1699 //***************************** 1700 1701 //------------------------------------------- 1702 // Execute event loop 1703 // 1704 MProgressBar bar; 1705 MEvtLoop evtloop; 1706 evtloop.SetParList(&pliston); 1707 evtloop.SetProgressBar(&bar); 1708 1709 Int_t maxevents = -1; 1710 //Int_t maxevents = 10000; 1711 if ( !evtloop.Eventloop(maxevents) ) 1712 return; 1713 1714 tliston.PrintStatistics(0, kTRUE); 1715 1716 1717 1718 //------------------------------------------- 1719 // Display the histograms 1720 // 1721 1722 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); 1723 pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); 1724 pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); 1725 1726 pliston.FindObject("MHHillas")->DrawClone(); 1727 pliston.FindObject("MHHillasExt")->DrawClone(); 1728 pliston.FindObject("MHHillasSrc")->DrawClone(); 1729 pliston.FindObject("MHNewImagePar")->DrawClone(); 1730 pliston.FindObject("MHStarMap")->DrawClone(); 1731 1732 //------------------------------------------- 1733 // fit alpha distribution to get the number of excess events 1734 // 1735 1736 MHHillasSrc* hillasSrc = 1737 (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); 1738 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); 1739 1740 MHOnSubtraction onsub; 1741 onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); 1742 //------------------------------------------- 1743 1744 DeleteBinnings(&pliston); 1745 1746 gLog << "Macro CT1Analysis : End of Job D" << endl; 1747 gLog << "=======================================================" << endl; 1748 } 1749 //--------------------------------------------------------------------- 1750 1751 1752 //--------------------------------------------------------------------- 1753 // Job E_EST_UP 1754 //============ 1755 1756 // - read MC2.root file 1757 // - select g/h separation method XX 1758 // - optimize energy estimator for events passing final cuts 1759 // - write parameters of energy estimator onto file "energyest_XX.root" 1760 // 1761 // - read ON2.root and MC2.root files 1762 // - update input root file with the estimated energy 1763 // (ON_XX2.root, MC_XX2.root) 1764 1765 1766 if (JobE_EST_UP) 1767 { 1768 gLog << "=====================================================" << endl; 1769 gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; 1770 1771 gLog << "" << endl; 1772 gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " 1773 << JobE_EST_UP << ", " << WESTUP << endl; 1774 1775 1776 TString typeON = "ON"; 1777 TString typeOFF = "OFF"; 1778 TString typeMC = "MC"; 1779 TString ext = "2.root"; 1780 TString extout = "3.root"; 1781 1782 //------------------------------ 1783 // name of MC file to be used for optimizing the energy estimator 1784 TString filenameOpt(outPath); 1785 filenameOpt += typeMC; 1786 filenameOpt += ext; 1787 gLog << "filenameOpt = " << filenameOpt << endl; 1788 1789 //------------------------------ 1790 // selection of g/h separation method 1791 // and definition of final selections 1792 1793 //TString XX("NN"); 1794 //TString XX("SC"); 1795 TString XX("RF"); 1796 TString fhadronnessName("Had"); 1797 fhadronnessName += XX; 1798 gLog << "fhadronnessName = " << fhadronnessName << endl; 1799 1800 // maximum values of the hadronness, |alpha| and dist 1801 Float_t maxhadronness = 0.40; 1802 Float_t maxalpha = 20.0; 1803 Float_t maxdist = 10.0; 1804 gLog << "Maximum values of hadronness, |ALPHA| and DIST = " 1805 << maxhadronness << ", " << maxalpha << ", " 1806 << maxdist << endl; 1807 1808 // name of file containing the parameters of the energy estimator 1809 TString energyParName(outPath); 1810 energyParName += "energyest_"; 1811 energyParName += XX; 1812 energyParName += ".root"; 1813 gLog << "energyParName = " << energyParName << endl; 1814 1815 1816 //------------------------------ 1817 // name of ON file to be updated 1818 TString filenameON(outPath); 1819 filenameON += typeON; 1820 filenameON += ext; 1821 gLog << "filenameON = " << filenameON << endl; 1822 1823 // name of OFF file to be updated 1824 TString filenameOFF(outPath); 1825 filenameOFF += typeOFF; 1826 filenameOFF += ext; 1827 gLog << "filenameOFF = " << filenameOFF << endl; 1828 1829 // name of MC file to be updated 1830 TString filenameMC(outPath); 1831 filenameMC += typeMC; 1832 filenameMC += ext; 1833 gLog << "filenameMC = " << filenameMC << endl; 1834 1835 //------------------------------ 1836 // name of updated ON file 1837 TString filenameONup(outPath); 1838 filenameONup += typeON; 1839 filenameONup += "_"; 1840 filenameONup += XX; 1841 filenameONup += extout; 1842 gLog << "filenameONup = " << filenameONup << endl; 1843 1844 // name of updated OFF file 1845 TString filenameOFFup(outPath); 1846 filenameOFFup += typeOFF; 1847 filenameOFFup += "_"; 1848 filenameOFFup += XX; 1849 filenameOFFup += extout; 1850 gLog << "filenameOFFup = " << filenameOFFup << endl; 1851 1852 // name of updated MC file 1853 TString filenameMCup(outPath); 1854 filenameMCup += typeMC; 1855 filenameMCup += "_"; 1856 filenameMCup += XX; 1857 filenameMCup += extout; 1858 gLog << "filenameMCup = " << filenameMCup << endl; 1859 1860 //----------------------------------------------------------- 1861 1862 TString fHilName = "MHillas"; 1863 TString fHilNameExt = "MHillasExt"; 1864 TString fHilNameSrc = "MHillasSrc"; 1865 TString fImgParName = "MNewImagePar"; 1866 1867 //=========================================================== 1868 // 1869 // Optimization of energy estimator 1870 // 1871 1872 TString inpath(""); 1873 TString outpath(""); 1874 Int_t howMany = 2000; 1875 CT1EEst(inpath, filenameOpt, outpath, energyParName, 1876 fHilName, fHilNameSrc, fhadronnessName, 1877 howMany, maxhadronness, maxalpha, maxdist); 1878 1879 //----------------------------------------------------------- 1880 // 1881 // Read in parameters of energy estimator 1882 // 1883 gLog << "================================================================" 1884 << endl; 1885 gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" 1886 << energyParName << "'" << endl; 1887 TFile enparam(energyParName); 1888 MMcEnergyEst mcest("MMcEnergyEst"); 1889 mcest.Read("MMcEnergyEst"); 1890 enparam.Close(); 1891 1892 TArrayD parA(5); 1893 TArrayD parB(7); 1894 for (Int_t i=0; i<parA.GetSize(); i++) 1895 parA[i] = mcest.GetCoeff(i); 1896 for (Int_t i=0; i<parB.GetSize(); i++) 1897 parB[i] = mcest.GetCoeff( i+parA.GetSize() ); 1898 1899 1900 if (WESTUP) 1901 { 1902 //========== start update ============================================ 1903 // 1904 // Update ON, OFF and MC root files with the estimated energy 1905 1906 //--------------------------------------------------- 1907 // Update OFF data 1908 // 1909 gLog << "============================================================" 1910 << endl; 1911 gLog << "Macro CT1Analysis.C : update file '" << filenameOFF 1912 << "'" << endl; 1913 1914 MTaskList tlistoff; 1915 MParList plistoff; 1916 1917 1918 // geometry is needed in MHHillas... classes 1919 MGeomCam *fGeom = 1920 (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1921 1922 //------------------------------------------- 1923 // create the tasks which should be executed 1924 // 1925 1926 MReadMarsFile read("Events", filenameOFF); 1927 read.DisableAutoScheme(); 1928 1929 //--------------------------- 1930 // calculate estimated energy 1931 1932 MEnergyEstParam eest2(fHilName); 1933 eest2.Add(fHilNameSrc); 1934 1935 eest2.SetCoeffA(parA); 1936 eest2.SetCoeffB(parB); 1937 1739 1938 //....................................................................... 1740 // calculate hadronnes for method of RANDOM FOREST 1741 1742 TString hadRFName = "HadRF"; 1743 MRanForestCalc rfcalc; 1744 rfcalc.SetHadronnessName(hadRFName); 1745 1746 //....................................................................... 1747 // calculation of hadroness for the supercuts 1748 // (=0.25 if fullfilled, =0.75 otherwise) 1749 1750 TString hadSCName = "HadSC"; 1751 MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); 1752 sccalc.SetHadronnessName(hadSCName); 1753 1754 1755 //....................................................................... 1756 1757 //MWriteRootFile write(outNameImage, "UPDATE"); 1758 MWriteRootFile write(outNameImage, "RECREATE"); 1939 1940 MWriteRootFile write(filenameOFFup); 1759 1941 1760 1942 write.AddContainer("MRawRunHeader", "RunHeaders"); … … 1769 1951 write.AddContainer("MNewImagePar", "Events"); 1770 1952 1953 //write.AddContainer("HadNN", "Events"); 1954 write.AddContainer("HadSC", "Events"); 1771 1955 write.AddContainer("HadRF", "Events"); 1772 write.AddContainer("HadSC", "Events"); 1956 1957 write.AddContainer("MEnergyEst", "Events"); 1773 1958 1774 1959 //----------------------------------------------------------------- 1960 1961 MFCT1SelFinal selfinal(fHilNameSrc); 1962 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 1963 selfinal.SetHadronnessName(fhadronnessName); 1964 MContinue contfinal(&selfinal); 1965 1966 1967 //***************************** 1968 // entries in MParList 1969 1970 plistoff.AddToList(&tlistoff); 1971 InitBinnings(&plistoff); 1972 1973 1974 //***************************** 1975 // entries in MTaskList 1976 1977 tlistoff.AddToList(&read); 1978 tlistoff.AddToList(&eest2); 1979 tlistoff.AddToList(&write); 1980 tlistoff.AddToList(&contfinal); 1981 1982 //***************************** 1983 1984 //------------------------------------------- 1985 // Execute event loop 1986 // 1987 MProgressBar bar; 1988 MEvtLoop evtloop; 1989 evtloop.SetParList(&plistoff); 1990 evtloop.SetProgressBar(&bar); 1991 1992 Int_t maxevents = -1; 1993 //Int_t maxevents = 1000; 1994 if ( !evtloop.Eventloop(maxevents) ) 1995 return; 1996 1997 tlistoff.PrintStatistics(0, kTRUE); 1998 DeleteBinnings(&plistoff); 1999 2000 //--------------------------------------------------- 2001 2002 //--------------------------------------------------- 2003 // Update ON data 2004 // 2005 gLog << "============================================================" 2006 << endl; 2007 gLog << "Macro CT1Analysis.C : update file '" << filenameON 2008 << "'" << endl; 2009 2010 MTaskList tliston; 2011 MParList pliston; 2012 2013 1775 2014 // geometry is needed in MHHillas... classes 1776 2015 MGeomCam *fGeom = 1777 2016 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1778 2017 1779 1780 Float_t maxhadronness = 0.40;1781 Float_t maxalpha = 20.0;1782 Float_t maxdist = 10.0;1783 1784 MFCT1SelFinal selfinalgh(fHilNameSrc);1785 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);1786 selfinalgh.SetHadronnessName(hadRFName);1787 selfinalgh.SetName("SelFinalgh");1788 MContinue contfinalgh(&selfinalgh);1789 contfinalgh.SetName("ContSelFinalgh");1790 1791 MFillH fillranfor("MHRanForest");1792 fillranfor.SetName("HRanForest");1793 1794 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);1795 fillhadrf.SetName("HhadRF");1796 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);1797 fillhadsc.SetName("HhadSC");1798 1799 1800 MFCT1SelFinal selfinal(fHilNameSrc);1801 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);1802 selfinal.SetHadronnessName(hadRFName);1803 selfinal.SetName("SelFinal");1804 MContinue contfinal(&selfinal);1805 contfinal.SetName("ContSelFinal");1806 1807 1808 MFillH hfill1("MHHillas", fHilName);1809 hfill1.SetName("HHillas");1810 1811 MFillH hfill2("MHStarMap", fHilName);1812 hfill2.SetName("HStarMap");1813 1814 MFillH hfill3("MHHillasExt", fHilNameSrc);1815 hfill3.SetName("HHillasExt");1816 1817 MFillH hfill4("MHHillasSrc", fHilNameSrc);1818 hfill4.SetName("HHillasSrc");1819 1820 MFillH hfill5("MHNewImagePar", fImgParName);1821 hfill5.SetName("HNewImagePar");1822 1823 //*****************************1824 // entries in MParList1825 1826 pliston.AddToList(&tliston);1827 InitBinnings(&pliston);1828 1829 pliston.AddToList(fRanForest);1830 pliston.AddToList(fRanTree);1831 1832 //*****************************1833 // entries in MTaskList1834 1835 tliston.AddToList(&read);1836 1837 tliston.AddToList(&rfcalc);1838 tliston.AddToList(&sccalc);1839 tliston.AddToList(&fillranfor);1840 tliston.AddToList(&fillhadrf);1841 tliston.AddToList(&fillhadsc);1842 1843 tliston.AddToList(&hfill1);1844 tliston.AddToList(&hfill2);1845 tliston.AddToList(&hfill3);1846 tliston.AddToList(&hfill4);1847 tliston.AddToList(&hfill5);1848 1849 tliston.AddToList(&write);1850 1851 tliston.AddToList(&contfinalgh);1852 tliston.AddToList(&contfinal);1853 1854 //*****************************1855 1856 //-------------------------------------------1857 // Execute event loop1858 //1859 MProgressBar bar;1860 MEvtLoop evtloop;1861 evtloop.SetParList(&pliston);1862 evtloop.SetProgressBar(&bar);1863 1864 Int_t maxevents = -1;1865 if ( !evtloop.Eventloop(maxevents) )1866 return;1867 1868 tliston.PrintStatistics(0, kTRUE);1869 1870 1871 //-------------------------------------------1872 // Display the histograms1873 //1874 pliston.FindObject("MHRanForest")->DrawClone();1875 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();1876 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();1877 1878 pliston.FindObject("MHHillas")->DrawClone();1879 pliston.FindObject("MHHillasExt")->DrawClone();1880 pliston.FindObject("MHHillasSrc")->DrawClone();1881 pliston.FindObject("MHNewImagePar")->DrawClone();1882 pliston.FindObject("MHStarMap")->DrawClone();1883 1884 DeleteBinnings(&pliston);1885 }1886 1887 gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;1888 gLog << "=======================================================" << endl;1889 }1890 //---------------------------------------------------------------------1891 1892 1893 1894 //---------------------------------------------------------------------1895 // Job C1896 //======1897 1898 // - read ON1 and MC1 data files1899 // which should have been updated to contain the hadronnesses1900 // for the method of NEAREST NEIGHBORS and for the SUOERCUTS1901 // - produce Neyman-Pearson plots1902 1903 if (JobC)1904 {1905 gLog << "=====================================================" << endl;1906 gLog << "Macro CT1Analysis : Start of Job C" << endl;1907 1908 gLog << "" << endl;1909 gLog << "Macro CT1Analysis : JobC = " << JobC << endl;1910 1911 1912 // name of input data file1913 TString filenameData = outPath;1914 filenameData += "OFF";1915 filenameData += "2.root";1916 gLog << "filenameData = " << filenameData << endl;1917 1918 // name of input MC file1919 TString filenameMC = outPath;1920 filenameMC += "MC";1921 filenameMC += "2.root";1922 gLog << "filenameMC = " << filenameMC << endl;1923 1924 1925 //-----------------------------------------------------------------1926 1927 MTaskList tliston;1928 MParList pliston;1929 1930 1931 // geometry is needed in MHHillas... classes1932 MGeomCam *fGeom =1933 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");1934 1935 2018 //------------------------------------------- 1936 2019 // create the tasks which should be executed 1937 2020 // 1938 2021 1939 MReadMarsFile read("Events", filenameMC); 1940 read.AddFile(filenameData); 1941 read.DisableAutoScheme(); 1942 1943 1944 //....................................................................... 1945 // names of hadronness containers 1946 1947 TString hadNNName = "HadNN"; 1948 TString hadSCName = "HadSC"; 1949 TString hadRFName = "HadRF"; 1950 1951 //....................................................................... 1952 1953 1954 TString fHilName = "MHillas"; 1955 TString fHilNameExt = "MHillasExt"; 1956 TString fHilNameSrc = "MHillasSrc"; 1957 TString fImgParName = "MNewImagePar"; 1958 1959 Float_t maxhadronness = 0.40; 1960 Float_t maxalpha = 20.0; 1961 Float_t maxdist = 10.0; 1962 1963 MFCT1SelFinal selfinalgh(fHilNameSrc); 1964 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); 1965 selfinalgh.SetHadronnessName(hadRFName); 1966 selfinalgh.SetName("SelFinalgh"); 1967 MContinue contfinalgh(&selfinalgh); 1968 contfinalgh.SetName("ContSelFinalgh"); 1969 1970 //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); 1971 //fillhadnn.SetName("HhadNN"); 1972 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); 1973 fillhadsc.SetName("HhadSC"); 1974 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); 1975 fillhadrf.SetName("HhadRF"); 1976 1977 MFCT1SelFinal selfinal(fHilNameSrc); 1978 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 1979 selfinal.SetHadronnessName(hadRFName); 1980 selfinal.SetName("SelFinal"); 1981 MContinue contfinal(&selfinal); 1982 contfinal.SetName("ContSelFinal"); 1983 1984 1985 MFillH hfill1("MHHillas", fHilName); 1986 hfill1.SetName("HHillas"); 1987 1988 MFillH hfill2("MHStarMap", fHilName); 1989 hfill2.SetName("HStarMap"); 1990 1991 MFillH hfill3("MHHillasExt", fHilNameSrc); 1992 hfill3.SetName("HHillasExt"); 1993 1994 MFillH hfill4("MHHillasSrc", fHilNameSrc); 1995 hfill4.SetName("HHillasSrc"); 1996 1997 MFillH hfill5("MHNewImagePar", fImgParName); 1998 hfill5.SetName("HNewImagePar"); 1999 2000 2001 //***************************** 2002 // entries in MParList 2003 2004 pliston.AddToList(&tliston); 2005 InitBinnings(&pliston); 2006 2007 2008 //***************************** 2009 // entries in MTaskList 2010 2011 tliston.AddToList(&read); 2012 2013 //tliston.AddToList(&fillhadnn); 2014 tliston.AddToList(&fillhadsc); 2015 tliston.AddToList(&fillhadrf); 2016 2017 tliston.AddToList(&contfinalgh); 2018 tliston.AddToList(&hfill1); 2019 tliston.AddToList(&hfill2); 2020 tliston.AddToList(&hfill3); 2021 tliston.AddToList(&hfill4); 2022 tliston.AddToList(&hfill5); 2023 2024 tliston.AddToList(&contfinal); 2025 2026 //***************************** 2027 2028 //------------------------------------------- 2029 // Execute event loop 2030 // 2031 MProgressBar bar; 2032 MEvtLoop evtloop; 2033 evtloop.SetParList(&pliston); 2034 evtloop.SetProgressBar(&bar); 2035 2036 Int_t maxevents = -1; 2037 //Int_t maxevents = 35000; 2038 if ( !evtloop.Eventloop(maxevents) ) 2039 return; 2040 2041 tliston.PrintStatistics(0, kTRUE); 2042 2043 2044 //------------------------------------------- 2045 // Display the histograms 2046 // 2047 2048 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); 2049 pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); 2050 pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); 2051 2052 pliston.FindObject("MHHillas")->DrawClone(); 2053 pliston.FindObject("MHHillasExt")->DrawClone(); 2054 pliston.FindObject("MHHillasSrc")->DrawClone(); 2055 pliston.FindObject("MHNewImagePar")->DrawClone(); 2056 pliston.FindObject("MHStarMap")->DrawClone(); 2057 2058 DeleteBinnings(&pliston); 2059 2060 gLog << "Macro CT1Analysis : End of Job C" << endl; 2061 gLog << "===================================================" << endl; 2062 } 2063 2064 2065 //--------------------------------------------------------------------- 2066 // Job D 2067 //====== 2068 2069 // - select g/h separation method XX 2070 // - read ON2 (or MC2) root file 2071 // - apply cuts in hadronness 2072 // - make plots 2073 2074 2075 if (JobD) 2076 { 2077 gLog << "=====================================================" << endl; 2078 gLog << "Macro CT1Analysis : Start of Job D" << endl; 2079 2080 gLog << "" << endl; 2081 gLog << "Macro CT1Analysis : JobD = " 2082 << JobD << endl; 2083 2084 // type of data to be analysed 2085 TString typeData = "ON"; 2086 //TString typeData = "OFF"; 2087 //TString typeData = "MC"; 2088 gLog << "typeData = " << typeData << endl; 2089 2090 TString ext = "2.root"; 2091 2092 2093 //------------------------------ 2094 // selection of g/h separation method 2095 // and definition of final selections 2096 2097 //TString XX("NN"); 2098 //TString XX("SC"); 2099 TString XX("RF"); 2100 TString fhadronnessName("Had"); 2101 fhadronnessName += XX; 2102 gLog << "fhadronnessName = " << fhadronnessName << endl; 2103 2104 // maximum values of the hadronness, |ALPHA| and DIST 2105 Float_t maxhadronness = 0.30; 2106 Float_t maxalpha = 20.0; 2107 Float_t maxdist = 10.0; 2108 gLog << "Maximum values of hadronness, |ALPHA| and DIST = " 2109 << maxhadronness << ", " << maxalpha << ", " 2110 << maxdist << endl; 2111 2112 2113 //------------------------------ 2114 // name of data file to be analysed 2115 TString filenameData(outPath); 2116 filenameData += typeData; 2117 filenameData += ext; 2118 gLog << "filenameData = " << filenameData << endl; 2119 2120 2121 2122 //************************************************************************* 2123 // 2124 // Analyse the data 2125 // 2126 2127 MTaskList tliston; 2128 MParList pliston; 2129 2130 // geometry is needed in MHHillas... classes 2131 MGeomCam *fGeom = 2132 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 2133 2134 2135 TString fHilName = "MHillas"; 2136 TString fHilNameExt = "MHillasExt"; 2137 TString fHilNameSrc = "MHillasSrc"; 2138 TString fImgParName = "MNewImagePar"; 2139 2140 //------------------------------------------- 2141 // create the tasks which should be executed 2142 // 2143 2144 MReadMarsFile read("Events", filenameData); 2145 read.DisableAutoScheme(); 2146 2147 2148 //----------------------------------------------------------------- 2149 // geometry is needed in MHHillas... classes 2150 MGeomCam *fGeom = 2151 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 2152 2153 MFCT1SelFinal selfinalgh(fHilNameSrc); 2154 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); 2155 selfinalgh.SetHadronnessName(fhadronnessName); 2156 selfinalgh.SetName("SelFinalgh"); 2157 MContinue contfinalgh(&selfinalgh); 2158 contfinalgh.SetName("ContSelFinalgh"); 2159 2160 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); 2161 //fillhadnn.SetName("HhadNN"); 2162 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); 2163 fillhadsc.SetName("HhadSC"); 2164 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); 2165 fillhadrf.SetName("HhadRF"); 2166 2167 2168 MFillH hfill1("MHHillas", fHilName); 2169 hfill1.SetName("HHillas"); 2170 2171 MFillH hfill2("MHStarMap", fHilName); 2172 hfill2.SetName("HStarMap"); 2173 2174 MFillH hfill3("MHHillasExt", fHilNameSrc); 2175 hfill3.SetName("HHillasExt"); 2176 2177 MFillH hfill4("MHHillasSrc", fHilNameSrc); 2178 hfill4.SetName("HHillasSrc"); 2179 2180 MFillH hfill5("MHNewImagePar", fImgParName); 2181 hfill5.SetName("HNewImagePar"); 2182 2183 MFCT1SelFinal selfinal(fHilNameSrc); 2184 selfinal.SetCuts(maxhadronness, maxalpha, maxdist); 2185 selfinal.SetHadronnessName(fhadronnessName); 2186 selfinal.SetName("SelFinal"); 2187 MContinue contfinal(&selfinal); 2188 contfinal.SetName("ContSelFinal"); 2189 2190 2191 //***************************** 2192 // entries in MParList 2193 2194 pliston.AddToList(&tliston); 2195 InitBinnings(&pliston); 2196 2197 2198 //***************************** 2199 // entries in MTaskList 2200 2201 tliston.AddToList(&read); 2202 2203 tliston.AddToList(&contfinalgh); 2204 2205 //tliston.AddToList(&fillhadnn); 2206 tliston.AddToList(&fillhadsc); 2207 tliston.AddToList(&fillhadrf); 2208 2209 tliston.AddToList(&hfill1); 2210 tliston.AddToList(&hfill2); 2211 tliston.AddToList(&hfill3); 2212 tliston.AddToList(&hfill4); 2213 tliston.AddToList(&hfill5); 2214 2215 tliston.AddToList(&contfinal); 2216 2217 //***************************** 2218 2219 //------------------------------------------- 2220 // Execute event loop 2221 // 2222 MProgressBar bar; 2223 MEvtLoop evtloop; 2224 evtloop.SetParList(&pliston); 2225 evtloop.SetProgressBar(&bar); 2226 2227 Int_t maxevents = -1; 2228 //Int_t maxevents = 10000; 2229 if ( !evtloop.Eventloop(maxevents) ) 2230 return; 2231 2232 tliston.PrintStatistics(0, kTRUE); 2233 2234 2235 2236 //------------------------------------------- 2237 // Display the histograms 2238 // 2239 2240 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); 2241 pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); 2242 pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); 2243 2244 pliston.FindObject("MHHillas")->DrawClone(); 2245 pliston.FindObject("MHHillasExt")->DrawClone(); 2246 pliston.FindObject("MHHillasSrc")->DrawClone(); 2247 pliston.FindObject("MHNewImagePar")->DrawClone(); 2248 pliston.FindObject("MHStarMap")->DrawClone(); 2249 2250 //------------------------------------------- 2251 // fit alpha distribution to get the number of excess events 2252 // 2253 2254 MHHillasSrc* hillasSrc = 2255 (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); 2256 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); 2257 2258 MHOnSubtraction onsub; 2259 onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); 2260 //------------------------------------------- 2261 2262 DeleteBinnings(&pliston); 2263 2264 gLog << "Macro CT1Analysis : End of Job D" << endl; 2265 gLog << "=======================================================" << endl; 2266 } 2267 //--------------------------------------------------------------------- 2268 2269 2270 //--------------------------------------------------------------------- 2271 // Job E_EST_UP 2272 //============ 2273 2274 // - read MC2.root file 2275 // - select g/h separation method XX 2276 // - optimize energy estimator for events passing final cuts 2277 // - write parameters of energy estimator onto file "energyest_XX.root" 2278 // 2279 // - read ON2.root and MC2.root files 2280 // - update input root file with the estimated energy 2281 // (ON_XX2.root, MC_XX2.root) 2282 2283 2284 if (JobE_EST_UP) 2285 { 2286 gLog << "=====================================================" << endl; 2287 gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; 2288 2289 gLog << "" << endl; 2290 gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " 2291 << JobE_EST_UP << ", " << WESTUP << endl; 2292 2293 2294 TString typeON = "ON"; 2295 TString typeOFF = "OFF"; 2296 TString typeMC = "MC"; 2297 TString ext = "2.root"; 2298 TString extout = "3.root"; 2299 2300 //------------------------------ 2301 // name of MC file to be used for optimizing the energy estimator 2302 TString filenameOpt(outPath); 2303 filenameOpt += typeMC; 2304 filenameOpt += ext; 2305 gLog << "filenameOpt = " << filenameOpt << endl; 2306 2307 //------------------------------ 2308 // selection of g/h separation method 2309 // and definition of final selections 2310 2311 //TString XX("NN"); 2312 //TString XX("SC"); 2313 TString XX("RF"); 2314 TString fhadronnessName("Had"); 2315 fhadronnessName += XX; 2316 gLog << "fhadronnessName = " << fhadronnessName << endl; 2317 2318 // maximum values of the hadronness, |alpha| and dist 2319 Float_t maxhadronness = 0.40; 2320 Float_t maxalpha = 20.0; 2321 Float_t maxdist = 10.0; 2322 gLog << "Maximum values of hadronness, |ALPHA| and DIST = " 2323 << maxhadronness << ", " << maxalpha << ", " 2324 << maxdist << endl; 2325 2326 // name of file containing the parameters of the energy estimator 2327 TString energyParName(outPath); 2328 energyParName += "energyest_"; 2329 energyParName += XX; 2330 energyParName += ".root"; 2331 gLog << "energyParName = " << energyParName << endl; 2332 2333 2334 //------------------------------ 2335 // name of ON file to be updated 2336 TString filenameON(outPath); 2337 filenameON += typeON; 2338 filenameON += ext; 2339 gLog << "filenameON = " << filenameON << endl; 2340 2341 // name of OFF file to be updated 2342 TString filenameOFF(outPath); 2343 filenameOFF += typeOFF; 2344 filenameOFF += ext; 2345 gLog << "filenameOFF = " << filenameOFF << endl; 2346 2347 // name of MC file to be updated 2348 TString filenameMC(outPath); 2349 filenameMC += typeMC; 2350 filenameMC += ext; 2351 gLog << "filenameMC = " << filenameMC << endl; 2352 2353 //------------------------------ 2354 // name of updated ON file 2355 TString filenameONup(outPath); 2356 filenameONup += typeON; 2357 filenameONup += "_"; 2358 filenameONup += XX; 2359 filenameONup += extout; 2360 gLog << "filenameONup = " << filenameONup << endl; 2361 2362 // name of updated OFF file 2363 TString filenameOFFup(outPath); 2364 filenameOFFup += typeOFF; 2365 filenameOFFup += "_"; 2366 filenameOFFup += XX; 2367 filenameOFFup += extout; 2368 gLog << "filenameOFFup = " << filenameOFFup << endl; 2369 2370 // name of updated MC file 2371 TString filenameMCup(outPath); 2372 filenameMCup += typeMC; 2373 filenameMCup += "_"; 2374 filenameMCup += XX; 2375 filenameMCup += extout; 2376 gLog << "filenameMCup = " << filenameMCup << endl; 2377 2378 //----------------------------------------------------------- 2379 2380 TString fHilName = "MHillas"; 2381 TString fHilNameExt = "MHillasExt"; 2382 TString fHilNameSrc = "MHillasSrc"; 2383 TString fImgParName = "MNewImagePar"; 2384 2385 //=========================================================== 2386 // 2387 // Optimization of energy estimator 2388 // 2389 2390 TString inpath(""); 2391 TString outpath(""); 2392 Int_t howMany = 2000; 2393 CT1EEst(inpath, filenameOpt, outpath, energyParName, 2394 fHilName, fHilNameSrc, fhadronnessName, 2395 howMany, maxhadronness, maxalpha, maxdist); 2396 2397 //----------------------------------------------------------- 2398 // 2399 // Read in parameters of energy estimator 2400 // 2401 gLog << "================================================================" 2402 << endl; 2403 gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" 2404 << energyParName << "'" << endl; 2405 TFile enparam(energyParName); 2406 MMcEnergyEst mcest("MMcEnergyEst"); 2407 mcest.Read("MMcEnergyEst"); 2408 enparam.Close(); 2409 2410 TArrayD parA(5); 2411 TArrayD parB(7); 2412 for (Int_t i=0; i<parA.GetSize(); i++) 2413 parA[i] = mcest.GetCoeff(i); 2414 for (Int_t i=0; i<parB.GetSize(); i++) 2415 parB[i] = mcest.GetCoeff( i+parA.GetSize() ); 2416 2417 2418 if (WESTUP) 2419 { 2420 //========== start update ============================================ 2421 // 2422 // Update ON, OFF and MC root files with the estimated energy 2423 2424 //--------------------------------------------------- 2425 // Update OFF data 2426 // 2427 gLog << "============================================================" 2428 << endl; 2429 gLog << "Macro CT1Analysis.C : update file '" << filenameOFF 2430 << "'" << endl; 2431 2432 MTaskList tlistoff; 2433 MParList plistoff; 2434 2435 2436 // geometry is needed in MHHillas... classes 2437 MGeomCam *fGeom = 2438 (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam"); 2439 2440 //------------------------------------------- 2441 // create the tasks which should be executed 2442 // 2443 2444 MReadMarsFile read("Events", filenameOFF); 2022 MReadMarsFile read("Events", filenameON); 2445 2023 read.DisableAutoScheme(); 2446 2024 … … 2456 2034 //....................................................................... 2457 2035 2458 MWriteRootFile write(filenameO FFup);2036 MWriteRootFile write(filenameONup); 2459 2037 2460 2038 write.AddContainer("MRawRunHeader", "RunHeaders"); … … 2486 2064 // entries in MParList 2487 2065 2488 plisto ff.AddToList(&tlistoff);2489 InitBinnings(&plisto ff);2066 pliston.AddToList(&tliston); 2067 InitBinnings(&pliston); 2490 2068 2491 2069 … … 2493 2071 // entries in MTaskList 2494 2072 2495 tlisto ff.AddToList(&read);2496 tlisto ff.AddToList(&eest2);2497 tlisto ff.AddToList(&write);2498 tlisto ff.AddToList(&contfinal);2073 tliston.AddToList(&read); 2074 tliston.AddToList(&eest2); 2075 tliston.AddToList(&write); 2076 tliston.AddToList(&contfinal); 2499 2077 2500 2078 //***************************** … … 2505 2083 MProgressBar bar; 2506 2084 MEvtLoop evtloop; 2507 evtloop.SetParList(&plisto ff);2085 evtloop.SetParList(&pliston); 2508 2086 evtloop.SetProgressBar(&bar); 2509 2087 … … 2513 2091 return; 2514 2092 2515 tlisto ff.PrintStatistics(0, kTRUE);2516 DeleteBinnings(&plisto ff);2093 tliston.PrintStatistics(0, kTRUE); 2094 DeleteBinnings(&pliston); 2517 2095 2518 2096 //--------------------------------------------------- 2519 2097 2520 2098 //--------------------------------------------------- 2521 // Update ONdata2099 // Update MC data 2522 2100 // 2523 2101 gLog << "============================================================" 2524 2102 << endl; 2525 gLog << "Macro CT1Analysis.C : update file '" << filename ON2103 gLog << "Macro CT1Analysis.C : update file '" << filenameMC 2526 2104 << "'" << endl; 2527 2105 2528 MTaskList tliston; 2529 MParList pliston; 2530 2531 2532 // geometry is needed in MHHillas... classes 2533 MGeomCam *fGeom = 2534 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 2106 MTaskList tlistmc; 2107 MParList plistmc; 2535 2108 2536 2109 //------------------------------------------- … … 2538 2111 // 2539 2112 2540 MReadMarsFile read("Events", filename ON);2113 MReadMarsFile read("Events", filenameMC); 2541 2114 read.DisableAutoScheme(); 2542 2115 … … 2552 2125 //....................................................................... 2553 2126 2554 MWriteRootFile write(filename ONup);2127 MWriteRootFile write(filenameMCup); 2555 2128 2556 2129 write.AddContainer("MRawRunHeader", "RunHeaders"); … … 2582 2155 // entries in MParList 2583 2156 2584 plist on.AddToList(&tliston);2585 InitBinnings(&plist on);2157 plistmc.AddToList(&tlistmc); 2158 InitBinnings(&plistmc); 2586 2159 2587 2160 … … 2589 2162 // entries in MTaskList 2590 2163 2591 tlist on.AddToList(&read);2592 tlist on.AddToList(&eest2);2593 tlist on.AddToList(&write);2594 tlist on.AddToList(&contfinal);2164 tlistmc.AddToList(&read); 2165 tlistmc.AddToList(&eest2); 2166 tlistmc.AddToList(&write); 2167 tlistmc.AddToList(&contfinal); 2595 2168 2596 2169 //***************************** … … 2601 2174 MProgressBar bar; 2602 2175 MEvtLoop evtloop; 2603 evtloop.SetParList(&plist on);2176 evtloop.SetParList(&plistmc); 2604 2177 evtloop.SetProgressBar(&bar); 2605 2178 … … 2609 2182 return; 2610 2183 2611 tliston.PrintStatistics(0, kTRUE); 2612 DeleteBinnings(&pliston); 2613 2614 //--------------------------------------------------- 2615 2616 //--------------------------------------------------- 2617 // Update MC data 2618 // 2619 gLog << "============================================================" 2620 << endl; 2621 gLog << "Macro CT1Analysis.C : update file '" << filenameMC 2622 << "'" << endl; 2623 2624 MTaskList tlistmc; 2625 MParList plistmc; 2184 tlistmc.PrintStatistics(0, kTRUE); 2185 DeleteBinnings(&plistmc); 2186 2187 2188 //========== end update ============================================ 2189 } 2190 2191 enparam.Close(); 2192 2193 gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl; 2194 gLog << "=======================================================" << endl; 2195 } 2196 //--------------------------------------------------------------------- 2197 2198 2199 //--------------------------------------------------------------------- 2200 // Job F_XX 2201 //========= 2202 2203 // - select g/h separation method XX 2204 // - read MC_XX2.root file 2205 // - calculate eff. collection area 2206 // - read ON_XX2.root file 2207 // - apply final cuts 2208 // - calculate flux 2209 // - write root file for ON data after final cuts (ON_XX3.root)) 2210 2211 2212 if (JobF_XX) 2213 { 2214 gLog << "=====================================================" << endl; 2215 gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; 2216 2217 gLog << "" << endl; 2218 gLog << "Macro CT1Analysis : JobF_XX, WXX = " 2219 << JobF_XX << ", " << WXX << endl; 2220 2221 // type of data to be analysed 2222 //TString typeData = "ON"; 2223 //TString typeData = "OFF"; 2224 TString typeData = "MC"; 2225 gLog << "typeData = " << typeData << endl; 2226 2227 TString typeMC = "MC"; 2228 TString ext = "3.root"; 2229 TString extout = "4.root"; 2230 2231 //------------------------------ 2232 // selection of g/h separation method 2233 // and definition of final selections 2234 2235 //TString XX("NN"); 2236 //TString XX("SC"); 2237 TString XX("RF"); 2238 TString fhadronnessName("Had"); 2239 fhadronnessName += XX; 2240 gLog << "fhadronnessName = " << fhadronnessName << endl; 2241 2242 // maximum values of the hadronness, |ALPHA| and DIST 2243 Float_t maxhadronness = 0.40; 2244 Float_t maxalpha = 20.0; 2245 Float_t maxdist = 10.0; 2246 gLog << "Maximum values of hadronness, |ALPHA| and DIST = " 2247 << maxhadronness << ", " << maxalpha << ", " 2248 << maxdist << endl; 2249 2250 2251 //------------------------------ 2252 // name of MC file to be used for calculating the eff. collection areas 2253 TString filenameArea(outPath); 2254 filenameArea += typeMC; 2255 filenameArea += "_"; 2256 filenameArea += XX; 2257 filenameArea += ext; 2258 gLog << "filenameArea = " << filenameArea << endl; 2259 2260 //------------------------------ 2261 // name of file containing the eff. collection areas 2262 TString collareaName(outPath); 2263 collareaName += "area_"; 2264 collareaName += XX; 2265 collareaName += ".root"; 2266 gLog << "collareaName = " << collareaName << endl; 2267 2268 //------------------------------ 2269 // name of data file to be analysed 2270 TString filenameData(outPath); 2271 filenameData += typeData; 2272 filenameData += "_"; 2273 filenameData += XX; 2274 filenameData += ext; 2275 gLog << "filenameData = " << filenameData << endl; 2276 2277 //------------------------------ 2278 // name of output data file (after the final cuts) 2279 TString filenameDataout(outPath); 2280 filenameDataout += typeData; 2281 filenameDataout += "_"; 2282 filenameDataout += XX; 2283 filenameDataout += extout; 2284 gLog << "filenameDataout = " << filenameDataout << endl; 2285 2286 2287 //==================================================================== 2288 gLog << "-----------------------------------------------" << endl; 2289 gLog << "Start calculation of effective collection areas" << endl; 2290 MParList parlist; 2291 MTaskList tasklist; 2292 2293 //--------------------------------------- 2294 // Setup the tasks to be executed 2295 // 2296 MReadMarsFile reader("Events", filenameArea); 2297 reader.DisableAutoScheme(); 2298 2299 MFCT1SelFinal cuthadrons; 2300 cuthadrons.SetHadronnessName(fhadronnessName); 2301 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); 2302 2303 MContinue conthadrons(&cuthadrons); 2304 2305 //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); 2306 //MHMcCT1CollectionArea* collarea; 2307 2308 MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); 2309 filler.SetName("CollectionArea"); 2310 2311 //******************************** 2312 // entries in MParList 2313 2314 parlist.AddToList(&tasklist); 2315 InitBinnings(&parlist); 2316 //parlist.AddToList(collarea); 2317 2318 //******************************** 2319 // entries in MTaskList 2320 2321 tasklist.AddToList(&reader); 2322 tasklist.AddToList(&conthadrons); 2323 tasklist.AddToList(&filler); 2324 2325 //******************************** 2326 2327 //----------------------------------------- 2328 // Execute event loop 2329 // 2330 MEvtLoop magic; 2331 magic.SetParList(&parlist); 2332 2333 MProgressBar bar; 2334 magic.SetProgressBar(&bar); 2335 if (!magic.Eventloop()) 2336 return; 2337 2338 tasklist.PrintStatistics(0, kTRUE); 2339 2340 // Calculate effective collection areas 2341 // and display the histograms 2342 // 2343 2344 MHMcCT1CollectionArea *collarea = 2345 (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); 2346 2347 collarea->CalcEfficiency(); 2348 collarea->DrawClone("lego"); 2349 2350 //--------------------------------------------- 2351 // Write histograms to a file 2352 // 2353 2354 TFile f(collareaName, "RECREATE"); 2355 collarea->GetHist()->Write(); 2356 collarea->GetHAll()->Write(); 2357 collarea->GetHSel()->Write(); 2358 f.Close(); 2359 2360 //delete collarea; 2361 2362 gLog << "Calculation of effective collection areas done" << endl; 2363 gLog << "-----------------------------------------------" << endl; 2364 //------------------------------------------------------------------ 2365 2366 2367 //************************************************************************* 2368 // 2369 // Analyse the data 2370 // 2371 2372 MTaskList tliston; 2373 MParList pliston; 2374 2375 // geometry is needed in MHHillas... classes 2376 MGeomCam *fGeom = 2377 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 2378 2379 TString fHilName = "MHillas"; 2380 TString fHilNameExt = "MHillasExt"; 2381 TString fHilNameSrc = "MHillasSrc"; 2382 TString fImgParName = "MNewImagePar"; 2626 2383 2627 2384 //------------------------------------------- … … 2629 2386 // 2630 2387 2631 MReadMarsFile read("Events", filename MC);2388 MReadMarsFile read("Events", filenameData); 2632 2389 read.DisableAutoScheme(); 2633 2390 2634 //---------------------------2635 // calculate estimated energy2636 2637 MEnergyEstParam eest2(fHilName);2638 eest2.Add(fHilNameSrc);2639 2640 eest2.SetCoeffA(parA);2641 eest2.SetCoeffB(parB);2642 2643 2391 //....................................................................... 2644 2392 2645 MWriteRootFile write(filenameMCup); 2393 2394 MWriteRootFile write(filenameDataout); 2646 2395 2647 2396 write.AddContainer("MRawRunHeader", "RunHeaders"); … … 2662 2411 write.AddContainer("MEnergyEst", "Events"); 2663 2412 2664 //-----------------------------------------------------------------2665 2666 MFCT1SelFinal selfinal(fHilNameSrc);2667 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);2668 selfinal.SetHadronnessName(fhadronnessName);2669 MContinue contfinal(&selfinal);2670 2671 2672 //*****************************2673 // entries in MParList2674 2675 plistmc.AddToList(&tlistmc);2676 InitBinnings(&plistmc);2677 2678 2679 //*****************************2680 // entries in MTaskList2681 2682 tlistmc.AddToList(&read);2683 tlistmc.AddToList(&eest2);2684 tlistmc.AddToList(&write);2685 tlistmc.AddToList(&contfinal);2686 2687 //*****************************2688 2689 //-------------------------------------------2690 // Execute event loop2691 //2692 MProgressBar bar;2693 MEvtLoop evtloop;2694 evtloop.SetParList(&plistmc);2695 evtloop.SetProgressBar(&bar);2696 2697 Int_t maxevents = -1;2698 //Int_t maxevents = 1000;2699 if ( !evtloop.Eventloop(maxevents) )2700 return;2701 2702 tlistmc.PrintStatistics(0, kTRUE);2703 DeleteBinnings(&plistmc);2704 2705 2706 //========== end update ============================================2707 }2708 2709 enparam.Close();2710 2711 gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;2712 gLog << "=======================================================" << endl;2713 }2714 //---------------------------------------------------------------------2715 2716 2717 //---------------------------------------------------------------------2718 // Job F_XX2719 //=========2720 2721 // - select g/h separation method XX2722 // - read MC_XX2.root file2723 // - calculate eff. collection area2724 // - read ON_XX2.root file2725 // - apply final cuts2726 // - calculate flux2727 // - write root file for ON data after final cuts (ON_XX3.root))2728 2729 2730 if (JobF_XX)2731 {2732 gLog << "=====================================================" << endl;2733 gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;2734 2735 gLog << "" << endl;2736 gLog << "Macro CT1Analysis : JobF_XX, WXX = "2737 << JobF_XX << ", " << WXX << endl;2738 2739 // type of data to be analysed2740 //TString typeData = "ON";2741 //TString typeData = "OFF";2742 TString typeData = "MC";2743 gLog << "typeData = " << typeData << endl;2744 2745 TString typeMC = "MC";2746 TString ext = "3.root";2747 TString extout = "4.root";2748 2749 //------------------------------2750 // selection of g/h separation method2751 // and definition of final selections2752 2753 //TString XX("NN");2754 //TString XX("SC");2755 TString XX("RF");2756 TString fhadronnessName("Had");2757 fhadronnessName += XX;2758 gLog << "fhadronnessName = " << fhadronnessName << endl;2759 2760 // maximum values of the hadronness, |ALPHA| and DIST2761 Float_t maxhadronness = 0.40;2762 Float_t maxalpha = 20.0;2763 Float_t maxdist = 10.0;2764 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "2765 << maxhadronness << ", " << maxalpha << ", "2766 << maxdist << endl;2767 2768 2769 //------------------------------2770 // name of MC file to be used for calculating the eff. collection areas2771 TString filenameArea(outPath);2772 filenameArea += typeMC;2773 filenameArea += "_";2774 filenameArea += XX;2775 filenameArea += ext;2776 gLog << "filenameArea = " << filenameArea << endl;2777 2778 //------------------------------2779 // name of file containing the eff. collection areas2780 TString collareaName(outPath);2781 collareaName += "area_";2782 collareaName += XX;2783 collareaName += ".root";2784 gLog << "collareaName = " << collareaName << endl;2785 2786 //------------------------------2787 // name of data file to be analysed2788 TString filenameData(outPath);2789 filenameData += typeData;2790 filenameData += "_";2791 filenameData += XX;2792 filenameData += ext;2793 gLog << "filenameData = " << filenameData << endl;2794 2795 //------------------------------2796 // name of output data file (after the final cuts)2797 TString filenameDataout(outPath);2798 filenameDataout += typeData;2799 filenameDataout += "_";2800 filenameDataout += XX;2801 filenameDataout += extout;2802 gLog << "filenameDataout = " << filenameDataout << endl;2803 2804 2805 //====================================================================2806 gLog << "-----------------------------------------------" << endl;2807 gLog << "Start calculation of effective collection areas" << endl;2808 MParList parlist;2809 MTaskList tasklist;2810 2811 //---------------------------------------2812 // Setup the tasks to be executed2813 //2814 MReadMarsFile reader("Events", filenameArea);2815 reader.DisableAutoScheme();2816 2817 MFCT1SelFinal cuthadrons;2818 cuthadrons.SetHadronnessName(fhadronnessName);2819 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);2820 2821 MContinue conthadrons(&cuthadrons);2822 2823 //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();2824 //MHMcCT1CollectionArea* collarea;2825 2826 MFillH filler("MHMcCT1CollectionArea", "MMcEvt");2827 filler.SetName("CollectionArea");2828 2829 //********************************2830 // entries in MParList2831 2832 parlist.AddToList(&tasklist);2833 InitBinnings(&parlist);2834 //parlist.AddToList(collarea);2835 2836 //********************************2837 // entries in MTaskList2838 2839 tasklist.AddToList(&reader);2840 tasklist.AddToList(&conthadrons);2841 tasklist.AddToList(&filler);2842 2843 //********************************2844 2845 //-----------------------------------------2846 // Execute event loop2847 //2848 MEvtLoop magic;2849 magic.SetParList(&parlist);2850 2851 MProgressBar bar;2852 magic.SetProgressBar(&bar);2853 if (!magic.Eventloop())2854 return;2855 2856 tasklist.PrintStatistics(0, kTRUE);2857 2858 // Calculate effective collection areas2859 // and display the histograms2860 //2861 2862 MHMcCT1CollectionArea *collarea =2863 (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");2864 2865 collarea->CalcEfficiency();2866 collarea->DrawClone("lego");2867 2868 //---------------------------------------------2869 // Write histograms to a file2870 //2871 2872 TFile f(collareaName, "RECREATE");2873 collarea->GetHist()->Write();2874 collarea->GetHAll()->Write();2875 collarea->GetHSel()->Write();2876 f.Close();2877 2878 //delete collarea;2879 2880 gLog << "Calculation of effective collection areas done" << endl;2881 gLog << "-----------------------------------------------" << endl;2882 //------------------------------------------------------------------2883 2884 2885 //*************************************************************************2886 //2887 // Analyse the data2888 //2889 2890 MTaskList tliston;2891 MParList pliston;2892 2893 // geometry is needed in MHHillas... classes2894 MGeomCam *fGeom =2895 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");2896 2897 TString fHilName = "MHillas";2898 TString fHilNameExt = "MHillasExt";2899 TString fHilNameSrc = "MHillasSrc";2900 TString fImgParName = "MNewImagePar";2901 2902 //-------------------------------------------2903 // create the tasks which should be executed2904 //2905 2906 MReadMarsFile read("Events", filenameData);2907 read.DisableAutoScheme();2908 2909 //.......................................................................2910 2911 2912 MWriteRootFile write(filenameDataout);2913 2914 write.AddContainer("MRawRunHeader", "RunHeaders");2915 write.AddContainer("MTime", "Events");2916 write.AddContainer("MMcEvt", "Events");2917 write.AddContainer("ThetaOrig", "Events");2918 write.AddContainer("MSrcPosCam", "Events");2919 write.AddContainer("MSigmabar", "Events");2920 write.AddContainer("MHillas", "Events");2921 write.AddContainer("MHillasExt", "Events");2922 write.AddContainer("MHillasSrc", "Events");2923 write.AddContainer("MNewImagePar", "Events");2924 2925 //write.AddContainer("HadNN", "Events");2926 write.AddContainer("HadSC", "Events");2927 write.AddContainer("HadRF", "Events");2928 2929 write.AddContainer("MEnergyEst", "Events");2930 2931 2413 2932 2414 //-----------------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.