Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2436) +++ trunk/MagicSoft/Mars/Changelog (revision 2437) @@ -1,3 +1,21 @@ -*-*- END OF LINE -*-*- + + 2003/10/28: Wolfgang Wittek + + * manalysis/MCT1PadONOFF.cc + - replace GetMeanRms() by GetPedestalRms() + - replace SetMeanRms() by SetPedestalRms() + - reactivate code which was commented out by tgb + (no compilation errors on Alpha OSF) + + * manalysis/AnalysisLinkDef.h + /Makefile + - put back MCT1PadONOFF + + * macros/CT1Analysis.C + /ONOFFCT1Analysis.C + - current versions of macros for the analysis of CT1 data + + 2003/10/26: Oscar Blanch Bigas @@ -51,5 +69,4 @@ - removed FillLevels - obsolete - added SetLevels instead - Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2436) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2437) @@ -1081,7 +1081,4 @@ fillmatg.SetName("fillGammas"); - MFillH fillmath("MatrixHadrons"); - fillmath.SetFilter(&flisth); - fillmath.SetName("fillHadrons"); @@ -1142,4 +1139,8 @@ selectorh.SetName("selectHadrons"); + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&flisth); + fillmath.SetName("fillHadrons"); + // entries in MFilterList @@ -1182,9 +1183,9 @@ - // write out training events + // write out matrices of training events gLog << "" << endl; gLog << "========================================================" << endl; - gLog << "Write out training events" << endl; + gLog << "Write out matrices of training events" << endl; @@ -1242,4 +1243,5 @@ // MEvtLoop evtlooptr; + evtlooptr.SetName("ReadRFTrees"); evtlooptr.SetParList(&plisttr); if (!evtlooptr.Eventloop()) @@ -1300,4 +1302,5 @@ // MEvtLoop treeloop; + treeloop.SetName("GrowRFTrees"); treeloop.SetParList(&plist2); Index: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2436) +++ trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2437) @@ -1,4 +1,5 @@ #include "CT1EgyEst.C" + void InitBinnings(MParList *plist) @@ -142,13 +143,13 @@ // path for input for Mars - TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/"; + TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/"; // path for output from Mars - TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/"; + TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/"; //----------------------------------------------- - TEnv env("macros/CT1env.rc"); - Bool_t printEnv = kFALSE; + //TEnv env("macros/CT1env.rc"); + //Bool_t printEnv = kFALSE; //************************************************************************ @@ -167,25 +168,16 @@ // (or OFF1.root or MC1.root)? - // Job B_NN_UP : read ON1.root (OFF1.root or MC1.root) file - // - depending on RlookNN : create (or read in) hadron and gamma matrix - // - calculate hadroness for method of NEAREST NEIGHBORS - // - update the input files with the hadroness (ON1.root, OFF1.root - // or MC1.root) - - Bool_t JobB_NN_UP = kFALSE; - Bool_t RLookNN = kFALSE; // read in look-alike events - Bool_t WNN = kFALSE; // update input root file ? - // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file - // - depending on RLook : create (or read in) hadron and gamma matrix - // - depending on RTree : create (or read in) trees + // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas + // - if CTrainRF = TRUE : create matrices of training events + // - if RTree = TRUE : read in trees, otherwise create trees // - calculate hadroness for method of RANDOM FOREST - // and for the SUPERCUTS - // - update the input files with the hadronesses (ON1.root, OFF1.root - // or MC1.root) + // - update the input files with the hadronesses (ON2.root, OFF2.root + // or MC2.root) Bool_t JobB_RF_UP = kFALSE; - Bool_t RLookRF = kFALSE; // read in look-alike events + Bool_t RTrainRF = kFALSE; // read in matrices of training events + Bool_t CTrainRF = kFALSE; // create matrices of training events Bool_t RTree = kFALSE; // read in trees Bool_t WRF = kFALSE; // update input root file ? @@ -265,6 +257,9 @@ gLog << "" << endl; gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = " - << JobA << ", " << WPad << ", " << RPad << ", " << Wout << endl; - + << (JobA ? "kTRUE" : "kFALSE") << ", " + << (WPad ? "kTRUE" : "kFALSE") << ", " + << (RPad ? "kTRUE" : "kFALSE") << ", " + << (Wout ? "kTRUE" : "kFALSE") << endl; + //-------------------------------------------------- @@ -334,7 +329,7 @@ MCT1ReadPreProc readON(fileON); - MFCT1SelBasic selthetaon; - selthetaon.SetCuts(-100.0, 29.5, 35.5); - MContinue contthetaon(&selthetaon); + //MFCT1SelBasic selthetaon; + //selthetaon.SetCuts(-100.0, 29.5, 35.5); + //MContinue contthetaon(&selthetaon); MBlindPixelCalc blindon; @@ -346,4 +341,5 @@ MHBlindPixels blindON("BlindPixelsON"); MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels"); + fillblindON.SetName("FillBlind"); MSigmabarCalc sigbarcalcon; @@ -351,5 +347,6 @@ MHSigmaTheta sigthON("SigmaThetaON"); MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); - + fillsigthetaON.SetName("FillSigTheta"); + //***************************** // entries in MParList @@ -380,5 +377,5 @@ Int_t maxeventson = -1; - //Int_t maxeventson = 20000; + //Int_t maxeventson = 10000; if ( !evtloopon.Eventloop(maxeventson) ) return; @@ -420,8 +417,6 @@ MHBlindPixels blindOFF("BlindPixelsOFF"); - MHBlindPixels *pblindOFF = &blindOFF; - gLog << "pblindOFF = " << pblindOFF << endl; - MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels"); + fillblindOFF.SetName("FillBlindOFF"); MSigmabarCalc sigbarcalcoff; @@ -429,5 +424,6 @@ MHSigmaTheta sigthOFF("SigmaThetaOFF"); MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); - + fillsigthetaOFF.SetName("FillSigThetaOFF"); + //***************************** // entries in MParList @@ -657,5 +653,5 @@ MEvtLoop evtloop; evtloop.SetParList(&pliston); - evtloop.ReadEnv(env, "", printEnv); + //evtloop.ReadEnv(env, "", printEnv); evtloop.SetProgressBar(&bar); // evtloop.Write(); @@ -694,26 +690,49 @@ + + //--------------------------------------------------------------------- - // Job B_NN_UP + // Job B_RF_UP //============ - // - read in the matrices of look-alike events for gammas and hadrons - - // then read ON1.root (or MC1.root) file - // - // - calculate the hadroness for the method of nearest neighbors - // - // - update input root file, including the hadroness - - - if (JobB_NN_UP) + + // - create (or read in) the matrices of training events for gammas + // and hadrons + // - create (or read in) the trees + // - then read ON1.root (or MC1.root) file + // - calculate the hadroness for the method of RANDOM FOREST + // - update input root file with the hadroness + + + if (JobB_RF_UP) { gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; + gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " - << JobB_NN_UP << ", " << RLookNN << ", " << WNN << endl; - + gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = " + << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " + << (RTrainRF ? "kTRUE" : "kFALSE") << ", " + << (CTrainRF ? "kTRUE" : "kFALSE") << ", " + << (RTree ? "kTRUE" : "kFALSE") << ", " + << (WRF ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // parameters for the random forest + Int_t NumTrees = 100; + Int_t NumTry = 3; + Int_t NdSize = 1; + + + TString hadRFName = "HadRF"; + Float_t maxhadronness = 0.23; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; @@ -721,7 +740,7 @@ // file to be updated (ON, OFF or MC) - //TString typeInput = "ON"; + TString typeInput = "ON"; //TString typeInput = "OFF"; - TString typeInput = "MC"; + //TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -736,5 +755,4 @@ outNameImage += typeInput; outNameImage += "2.root"; - //TString outNameImage = filenameData; @@ -742,10 +760,10 @@ //-------------------------------------------- - // files to be read for generating the look-alike events + // files to be read for generating the matrices of training events // "hadrons" : TString filenameHad = outPath; filenameHad += "OFF"; filenameHad += "1.root"; - Int_t howManyHadrons = 500000; + Int_t howManyHadrons = 1000000; gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " << howManyHadrons << endl; @@ -761,8 +779,8 @@ //-------------------------------------------- - // files of look-alike events + // files of training events TString outNameGammas = outPath; - outNameGammas += "matrix_gammas_"; + outNameGammas += "RFmatrix_gammas_"; outNameGammas += "MC"; outNameGammas += ".root"; @@ -772,5 +790,5 @@ TString outNameHadrons = outPath; - outNameHadrons += "matrix_hadrons_"; + outNameHadrons += "RFmatrix_hadrons_"; outNameHadrons += typeMatrixHadrons; outNameHadrons += ".root"; @@ -778,9 +796,32 @@ MHMatrix matrixg("MatrixGammas"); + matrixg.EnableGraphicalOutput(); + + matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrixg.AddColumn("MSigmabar.fSigmabar"); + matrixg.AddColumn("log10(MHillas.fSize)"); + matrixg.AddColumn("MHillasSrc.fDist"); + matrixg.AddColumn("MHillas.fWidth"); + matrixg.AddColumn("MHillas.fLength"); + matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + matrixg.AddColumn("MNewImagePar.fConc"); + matrixg.AddColumn("MNewImagePar.fLeakage1"); + MHMatrix matrixh("MatrixHadrons"); + matrixh.EnableGraphicalOutput(); + + matrixh.AddColumns(matrixg.GetColumns()); + + //-------------------------------------------- + // file of trees of the random forest + + TString outRF = outPath; + outRF += "RF.root"; + //************************************************************************* - // read in matrices of look-alike events -if (RLookNN) + // read in matrices of training events +if (RTrainRF) { const char* mtxName = "MatrixGammas"; @@ -799,5 +840,5 @@ matrixg.Read(mtxName); - matrixg.Print("cols"); + matrixg.Print("SizeCols"); @@ -819,15 +860,15 @@ matrixh.Read(mtxName); - matrixh.Print("cols"); + matrixh.Print("SizeCols"); } //************************************************************************* - // create matrices of look-alike events -if (!RLookNN) + // create matrices of training events +if (CTrainRF) { gLog << "" << endl; gLog << "========================================================" << endl; - gLog << " Create matrices of look-alike events" << endl; + gLog << " Create matrices of training events" << endl; gLog << " Gammas :" << endl; @@ -852,18 +893,21 @@ fhadrons.SetName("hadronID)"); - MFEventSelector selectorg; - selectorg.SetNumSelectEvts(howManyGammas); + TString mgname("costhg"); + MBinning bing("Binning"+mgname); + bing.SetEdges(10, 0., 1.0); + + MH3 gref("cos(MMcEvt.fTelescopeTheta)"); + gref.SetName(mgname); + MH::SetBinning(&gref.GetHist(), &bing); + for (Int_t i=1; i<=gref.GetNbins(); i++) + gref.GetHist().SetBinContent(i, 1.0); + + MFEventSelector2 selectorg(gref); + selectorg.SetNumMax(howManyGammas); selectorg.SetName("selectGammas"); - MFEventSelector selectorh; - selectorh.SetNumSelectEvts(howManyHadrons); - selectorh.SetName("selectHadrons"); MFillH fillmatg("MatrixGammas"); fillmatg.SetFilter(&flistg); fillmatg.SetName("fillGammas"); - - MFillH fillmath("MatrixHadrons"); - fillmath.SetFilter(&flisth); - fillmath.SetName("fillHadrons"); @@ -894,5 +938,5 @@ MEvtLoop evtloopg; evtloopg.SetParList(&plistg); - evtloopg.ReadEnv(env, "", printEnv); + //evtloopg.ReadEnv(env, "", printEnv); evtloopg.SetProgressBar(&matrixbar); @@ -907,4 +951,26 @@ gLog << " Hadrons :" << endl; + + TString mhname("costhh"); + MBinning binh("Binning"+mhname); + binh.SetEdges(10, 0., 1.0); + + MH3 href("cos(MMcEvt.fTelescopeTheta)"); + href.SetName(mhname); + MH::SetBinning(&href.GetHist(), &binh); + for (Int_t i=1; i<=href.GetNbins(); i++) + href.GetHist().SetBinContent(i, 1.0); + + MFEventSelector2 selectorh(href); + //selectorh.SetNumMax(howManyHadrons); + // select as many hadrons as gammas + selectorh.SetNumMax(matrixg.GetM().GetNrows()); + selectorh.SetName("selectHadrons"); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&flisth); + fillmath.SetName("fillHadrons"); + + // entries in MFilterList @@ -932,5 +998,5 @@ MEvtLoop evtlooph; evtlooph.SetParList(&plisth); - evtlooph.ReadEnv(env, "", printEnv); + //evtlooph.ReadEnv(env, "", printEnv); evtlooph.SetProgressBar(&matrixbar); @@ -942,84 +1008,10 @@ //***************************************************** - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - // - // ---------------------------------------------------------- - // Definition of the reference sample of look-alike events - // (this is a subsample of the original sample) - // ---------------------------------------------------------- - // + + // write out matrices of training events events + gLog << "" << endl; gLog << "========================================================" << endl; - // Select a maximum of nmaxevts events from the sample of look-alike - // events. They will form the reference sample. - Int_t nmaxevts = 2000; - - // target distribution for the variable in column refcolumn (start at 0); - //Int_t refcolumn = 0; - //Int_t nbins = 5; - //Float_t frombin = 0.5; - //Float_t tobin = 1.0; - //TH1F *thsh = new TH1F("thsh","target distribution", - // nbins, frombin, tobin); - //thsh->SetDirectory(NULL); - //thsh->SetXTitle("cos( \\Theta)"); - //thsh->SetYTitle("Counts"); - //Float_t dbin = (tobin-frombin)/nbins; - //Float_t lbin = frombin +dbin*0.5; - //for (Int_t j=1; j<=nbins; j++) - //{ - // thsh->Fill(lbin,1.0); - // lbin += dbin; - //} - - Int_t refcolumn = 0; - MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); - TH1F *thsh = new TH1F(); - thsh->SetNameTitle("thsh","target distribution"); - MH::SetBinning(thsh, binscostheta); - Int_t nbins = thsh->GetNbinsX(); - Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; - for (Int_t j=1; j<=nbins; j++) - { - thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixg.EnableGraphicalOutput(); - Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; - return; - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixh.EnableGraphicalOutput(); - matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); - delete thsh; - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; - return; - } - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - - // write out look-alike events - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Write out look=alike events" << endl; + gLog << "Write out matrices of training events" << endl; @@ -1027,5 +1019,5 @@ // "gammas" gLog << "Gammas :" << endl; - matrixg.Print("cols"); + matrixg.Print("SizeCols"); TFile writeg(outNameGammas, "RECREATE", ""); @@ -1033,5 +1025,5 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " + gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file " << outNameGammas << endl; @@ -1039,5 +1031,5 @@ // "hadrons" gLog << "Hadrons :" << endl; - matrixh.Print("cols"); + matrixh.Print("SizeCols"); TFile writeh(outNameHadrons, "RECREATE", ""); @@ -1045,21 +1037,135 @@ gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " + gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file " << outNameHadrons << endl; } - //********** end of creating matrices of look-alike events *********** - - + //********** end of creating matrices of training events *********** + + + MRanForest *fRanForest; + MRanTree *fRanTree; //----------------------------------------------------------------- - // Update the input files with the NN and SC hadronness - // - - if (WNN) + // read in the trees of the random forest + if (RTree) + { + MParList plisttr; + MTaskList tlisttr; + plisttr.AddToList(&tlisttr); + + MReadTree readtr("TREE", outRF); + readtr.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(NumTrees); + + // list of tasks for the loop over the trees + + tlisttr.AddToList(&readtr); + tlisttr.AddToList(&rffill); + + //------------------- + // Execute tree loop + // + MEvtLoop evtlooptr; + evtlooptr.SetName("ReadRFTrees"); + evtlooptr.SetParList(&plisttr); + if (!evtlooptr.Eventloop()) + return; + + tlisttr.PrintStatistics(0, kTRUE); + + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + + //----------------------------------------------------------------- + // grow the trees of the random forest (event loop = tree loop) + + if (!RTree) + { + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : start growing trees" << endl; + + MTaskList tlist2; + MParList plist2; + plist2.AddToList(&tlist2); + + plist2.AddToList(&matrixg); + plist2.AddToList(&matrixh); + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(NumTrees); + rfgrow2.SetNumTry(NumTry); + rfgrow2.SetNdSize(NdSize); + + MWriteRootFile rfwrite2(outRF); + rfwrite2.AddContainer("MRanTree", "TREE"); + + // list of tasks for the loop over the trees + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + + //------------------- + // Execute tree loop + // + MEvtLoop treeloop; + treeloop.SetName("GrowRFTrees"); + treeloop.SetParList(&plist2); + + if ( !treeloop.Eventloop() ) + return; + + tlist2.PrintStatistics(0, kTRUE); + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + // end of growing the trees of the random forest + //----------------------------------------------------------------- + + + + + //----------------------------------------------------------------- + // Update the input files with the RF hadronness + // + + if (WRF) { gLog << "" << endl; gLog << "========================================================" << endl; gLog << "Update input file '" << filenameData - << "' with the NN and SC hadronness" << endl; + << "' with the RF hadronness" << endl; MTaskList tliston; @@ -1078,17 +1184,12 @@ read.DisableAutoScheme(); - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; //....................................................................... - // calculation of hadroness for method of Nearest Neighbors - - TString hadNNName = "HadNN"; - MMultiDimDistCalc nncalc; - nncalc.SetUseNumRows(25); - nncalc.SetUseKernelMethod(kFALSE); - nncalc.SetHadronnessName(hadNNName); + // calculate hadronnes for method of RANDOM FOREST + + + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadRFName); + //....................................................................... @@ -1108,8 +1209,5 @@ write.AddContainer("MNewImagePar", "Events"); - write.AddContainer("HadRF", "Events"); - write.AddContainer("HadSC", "Events"); - write.AddContainer("HadNN", "Events"); - + write.AddContainer(hadRFName, "Events"); //----------------------------------------------------------------- @@ -1118,24 +1216,34 @@ (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; + MFCT1SelFinal selfinalgh(fHilNameSrc); selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadNNName); + selfinalgh.SetHadronnessName(hadRFName); selfinalgh.SetName("SelFinalgh"); MContinue contfinalgh(&selfinalgh); contfinalgh.SetName("ContSelFinalgh"); - MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); - fillhadnn.SetName("HhadNN"); + MFillH fillranfor("MHRanForest"); + fillranfor.SetName("HRanForest"); + + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); MFCT1SelFinal selfinal(fHilNameSrc); selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadNNName); + selfinal.SetHadronnessName(hadRFName); selfinal.SetName("SelFinal"); MContinue contfinal(&selfinal); contfinal.SetName("ContSelFinal"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); @@ -1161,6 +1269,9 @@ InitBinnings(&pliston); - pliston.AddToList(&matrixg); - pliston.AddToList(&matrixh); + pliston.AddToList(fRanForest); + pliston.AddToList(fRanTree); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); @@ -1170,7 +1281,12 @@ tliston.AddToList(&read); - tliston.AddToList(&nncalc); - tliston.AddToList(&fillhadnn); - + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillranfor); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); @@ -1179,7 +1295,4 @@ tliston.AddToList(&hfill5); - tliston.AddToList(&write); - - tliston.AddToList(&contfinalgh); tliston.AddToList(&contfinal); @@ -1200,8 +1313,10 @@ tliston.PrintStatistics(0, kTRUE); + //------------------------------------------- // Display the histograms // - pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("MHRanForest")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); pliston.FindObject("MHHillas")->DrawClone(); @@ -1211,10 +1326,46 @@ pliston.FindObject("MHStarMap")->DrawClone(); + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); + TH1 &alphaHist = absalpha->GetHist(); + alphaHist.SetXTitle("|alpha| [\\circ]"); + alphaHist.SetName("alpha-macro"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << filenameData << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + DeleteBinnings(&pliston); } - - - gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl; + gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; gLog << "=======================================================" << endl; } @@ -1222,498 +1373,37 @@ + //--------------------------------------------------------------------- - // Job B_RF_UP - //============ - - - // - create (or read in) the matrices of look-alike events for gammas - // and hadrons - // - create (or read in) the trees - // - then read ON1.root (or MC1.root) file - // - calculate the hadroness for the method of RANDOM FOREST - // - update input root file with the hadroness - - - if (JobB_RF_UP) + // Job C + //====== + + // - read ON1 and MC1 data files + // which should have been updated to contain the hadronnesses + // for the method of NEAREST NEIGHBORS and for the SUOERCUTS + // - produce Neyman-Pearson plots + + if (JobC) { gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; + gLog << "Macro CT1Analysis : Start of Job C" << endl; gLog << "" << endl; - gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " - << JobB_RF_UP << ", " << RLookRF << ", " << RTree << ", " - << WRF << endl; - - - - //-------------------------------------------- - // parameters for the random forest - Int_t NumTrees = 100; - Int_t NumTry = 3; - Int_t NdSize = 3; - - - //-------------------------------------------- - // file to be updated (ON, OFF or MC) - - TString typeInput = "ON"; - //TString typeInput = "OFF"; - //TString typeInput = "MC"; - gLog << "typeInput = " << typeInput << endl; - - // name of input root file + gLog << "Macro CT1Analysis : JobC = " << JobC << endl; + + + // name of input data file TString filenameData = outPath; - filenameData += typeInput; - filenameData += "1.root"; - gLog << "filenameData = " << filenameData << endl; - - // name of output root file - TString outNameImage = outPath; - outNameImage += typeInput; - outNameImage += "2.root"; - //TString outNameImage = filenameData; - - gLog << "outNameImage = " << outNameImage << endl; - - //-------------------------------------------- - // files to be read for generating the look-alike events - // "hadrons" : - TString filenameHad = outPath; - filenameHad += "OFF"; - filenameHad += "1.root"; - Int_t howManyHadrons = 1000000; - gLog << "filenameHad = " << filenameHad << ", howManyHadrons = " - << howManyHadrons << endl; - - - // "gammas" : + filenameData += "OFF"; + filenameData += "2.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of input MC file TString filenameMC = outPath; filenameMC += "MC"; - filenameMC += "1.root"; - Int_t howManyGammas = 50000; - gLog << "filenameMC = " << filenameMC << ", howManyGammas = " - << howManyGammas << endl; - - //-------------------------------------------- - // files of look-alike events - - TString outNameGammas = outPath; - outNameGammas += "RFmatrix_gammas_"; - outNameGammas += "MC"; - outNameGammas += ".root"; - - TString typeMatrixHadrons = "OFF"; - gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; - - TString outNameHadrons = outPath; - outNameHadrons += "RFmatrix_hadrons_"; - outNameHadrons += typeMatrixHadrons; - outNameHadrons += ".root"; - - - MHMatrix matrixg("MatrixGammas"); - MHMatrix matrixh("MatrixHadrons"); - - //-------------------------------------------- - // file of trees of the random forest - - TString outRF = outPath; - outRF += "RF.root"; - - - //************************************************************************* - // read in matrices of look-alike events -if (RLookRF) - { - const char* mtxName = "MatrixGammas"; - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Get matrix for (gammas)" << endl; - gLog << "matrix name = " << mtxName << endl; - gLog << "name of root file = " << outNameGammas << endl; - gLog << "" << endl; - - - // read in the object with the name 'mtxName' from file 'outNameGammas' - // - TFile fileg(outNameGammas); - - matrixg.Read(mtxName); - matrixg.Print("cols"); - - - //***************************************************************** - - const char* mtxName = "MatrixHadrons"; - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << " Get matrix for (hadrons)" << endl; - gLog << "matrix name = " << mtxName << endl; - gLog << "name of root file = " << outNameHadrons << endl; - gLog << "" << endl; - - - // read in the object with the name 'mtxName' from file 'outNameHadrons' - // - TFile fileh(outNameHadrons); - - matrixh.Read(mtxName); - matrixh.Print("cols"); - } - - - //************************************************************************* - // create matrices of look-alike events -if (!RLookRF) - { - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << " Create matrices of look-alike events" << endl; - gLog << " Gammas :" << endl; - - - MParList plistg; - MTaskList tlistg; - MFilterList flistg; - - MParList plisth; - MTaskList tlisth; - MFilterList flisth; - - MReadMarsFile readg("Events", filenameMC); - readg.DisableAutoScheme(); - - MReadMarsFile readh("Events", filenameHad); - readh.DisableAutoScheme(); - - MFParticleId fgamma("MMcEvt", '=', kGAMMA); - fgamma.SetName("gammaID"); - MFParticleId fhadrons("MMcEvt", '!', kGAMMA); - fhadrons.SetName("hadronID)"); - - MFEventSelector selectorg; - selectorg.SetNumSelectEvts(howManyGammas); - selectorg.SetName("selectGammas"); - MFEventSelector selectorh; - selectorh.SetNumSelectEvts(howManyHadrons); - selectorh.SetName("selectHadrons"); - - MFillH fillmatg("MatrixGammas"); - fillmatg.SetFilter(&flistg); - fillmatg.SetName("fillGammas"); - - MFillH fillmath("MatrixHadrons"); - fillmath.SetFilter(&flisth); - fillmath.SetName("fillHadrons"); - - - //***************************** fill gammas *** - // entries in MFilterList - - flistg.AddToList(&fgamma); - flistg.AddToList(&selectorg); - - //***************************** - // entries in MParList - - plistg.AddToList(&tlistg); - InitBinnings(&plistg); - - plistg.AddToList(&matrixg); - - //***************************** - // entries in MTaskList - - tlistg.AddToList(&readg); - tlistg.AddToList(&flistg); - tlistg.AddToList(&fillmatg); - - //***************************** - - MProgressBar matrixbar; - MEvtLoop evtloopg; - evtloopg.SetParList(&plistg); - evtloopg.ReadEnv(env, "", printEnv); - evtloopg.SetProgressBar(&matrixbar); - - Int_t maxevents = -1; - if (!evtloopg.Eventloop(maxevents)) - return; - - tlistg.PrintStatistics(0, kTRUE); - - - //***************************** fill hadrons *** - - gLog << " Hadrons :" << endl; - // entries in MFilterList - - flisth.AddToList(&fhadrons); - flisth.AddToList(&selectorh); - - //***************************** - // entries in MParList - - plisth.AddToList(&tlisth); - InitBinnings(&plisth); - - plisth.AddToList(&matrixh); - - //***************************** - // entries in MTaskList - - tlisth.AddToList(&readh); - tlisth.AddToList(&flisth); - tlisth.AddToList(&fillmath); - - //***************************** - - MProgressBar matrixbar; - MEvtLoop evtlooph; - evtlooph.SetParList(&plisth); - evtlooph.ReadEnv(env, "", printEnv); - evtlooph.SetProgressBar(&matrixbar); - - Int_t maxevents = -1; - if (!evtlooph.Eventloop(maxevents)) - return; - - tlisth.PrintStatistics(0, kTRUE); - //***************************************************** - - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - // - // ---------------------------------------------------------- - // Definition of the reference sample of look-alike events - // (this is a subsample of the original sample) - // ---------------------------------------------------------- - // - gLog << "" << endl; - gLog << "========================================================" << endl; - // Select a maximum of nmaxevts events from the sample of look-alike - // events. They will form the reference sample. - Int_t nmaxevts = 10000; - - // target distribution for the variable in column refcolumn (start at 0); - //Int_t refcolumn = 0; - //Int_t nbins = 5; - //Float_t frombin = 0.5; - //Float_t tobin = 1.0; - //TH1F *thsh = new TH1F("thsh","target distribution", - // nbins, frombin, tobin); - //thsh->SetDirectory(NULL); - //thsh->SetXTitle("cos( \\Theta)"); - //thsh->SetYTitle("Counts"); - //Float_t dbin = (tobin-frombin)/nbins; - //Float_t lbin = frombin +dbin*0.5; - //for (Int_t j=1; j<=nbins; j++) - //{ - // thsh->Fill(lbin,1.0); - // lbin += dbin; - //} - - Int_t refcolumn = 0; - MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); - TH1F *thsh = new TH1F(); - thsh->SetNameTitle("thsh","target distribution"); - MH::SetBinning(thsh, binscostheta); - Int_t nbins = thsh->GetNbinsX(); - Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; - for (Int_t j=1; j<=nbins; j++) - { - thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixg.EnableGraphicalOutput(); - Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; - return; - } - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; - gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - matrixh.EnableGraphicalOutput(); - matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); - delete thsh; - - if ( !matrixok ) - { - gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; - return; - } - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - - // write out look-alike events - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Write out look=alike events" << endl; - - - //------------------------------------------- - // "gammas" - gLog << "Gammas :" << endl; - matrixg.Print("cols"); - - TFile writeg(outNameGammas, "RECREATE", ""); - matrixg.Write(); - - gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " - << outNameGammas << endl; - - //------------------------------------------- - // "hadrons" - gLog << "Hadrons :" << endl; - matrixh.Print("cols"); - - TFile writeh(outNameHadrons, "RECREATE", ""); - matrixh.Write(); - - gLog << "" << endl; - gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " - << outNameHadrons << endl; - - } - //********** end of creating matrices of look-alike events *********** - - - MRanForest *fRanForest; - MRanTree *fRanTree; + filenameMC += "2.root"; + gLog << "filenameMC = " << filenameMC << endl; + + //----------------------------------------------------------------- - // read in the trees of the random forest - if (RTree) - { - MParList plisttr; - MTaskList tlisttr; - plisttr.AddToList(&tlisttr); - - MReadTree readtr("TREE", outRF); - readtr.DisableAutoScheme(); - - MRanForestFill rffill; - rffill.SetNumTrees(NumTrees); - - // list of tasks for the loop over the trees - - tlisttr.AddToList(&readtr); - tlisttr.AddToList(&rffill); - - //------------------- - // Execute tree loop - // - MEvtLoop evtlooptr; - evtlooptr.SetParList(&plisttr); - if (!evtlooptr.Eventloop()) - return; - - tlisttr.PrintStatistics(0, kTRUE); - - - // get adresses of objects which are used in the next eventloop - fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); - if (!fRanForest) - { - gLog << err << dbginf << "MRanForest not found... aborting." << endl; - return kFALSE; - } - - fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); - if (!fRanTree) - { - gLog << err << dbginf << "MRanTree not found... aborting." << endl; - return kFALSE; - } - - } - - //----------------------------------------------------------------- - // grow the trees of the random forest (event loop = tree loop) - - if (!RTree) - { - - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Macro CT1Analysis : start growing trees" << endl; - - MTaskList tlist2; - MParList plist2; - plist2.AddToList(&tlist2); - - plist2.AddToList(&matrixg); - plist2.AddToList(&matrixh); - - MRanForestGrow rfgrow2; - rfgrow2.SetNumTrees(NumTrees); - rfgrow2.SetNumTry(NumTry); - rfgrow2.SetNdSize(NdSize); - - MWriteRootFile rfwrite2(outRF); - rfwrite2.AddContainer("MRanTree", "TREE"); - - // list of tasks for the loop over the trees - - tlist2.AddToList(&rfgrow2); - tlist2.AddToList(&rfwrite2); - - //------------------- - // Execute tree loop - // - MEvtLoop treeloop; - treeloop.SetParList(&plist2); - - if ( !treeloop.Eventloop() ) - return; - - tlist2.PrintStatistics(0, kTRUE); - - // get adresses of objects which are used in the next eventloop - fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); - if (!fRanForest) - { - gLog << err << dbginf << "MRanForest not found... aborting." << endl; - return kFALSE; - } - - fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); - if (!fRanTree) - { - gLog << err << dbginf << "MRanTree not found... aborting." << endl; - return kFALSE; - } - - } - // end of growing the trees of the random forest - //----------------------------------------------------------------- - - - - - //----------------------------------------------------------------- - // Update the input files with the RF hadronness - // - - if (WRF) - { - gLog << "" << endl; - gLog << "========================================================" << endl; - gLog << "Update input file '" << filenameData - << "' with the RF hadronness" << endl; MTaskList tliston; @@ -1729,6 +1419,18 @@ // - MReadMarsFile read("Events", filenameData); + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); read.DisableAutoScheme(); + + + //....................................................................... + // names of hadronness containers + + TString hadNNName = "HadNN"; + TString hadSCName = "HadSC"; + TString hadRFName = "HadRF"; + + //....................................................................... + TString fHilName = "MHillas"; @@ -1737,24 +1439,504 @@ TString fImgParName = "MNewImagePar"; + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadRFName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + //fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadRFName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + //tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "Macro CT1Analysis : End of Job C" << endl; + gLog << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job D + //====== + + // - select g/h separation method XX + // - read ON2 (or MC2) root file + // - apply cuts in hadronness + // - make plots + + + if (JobD) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job D" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobD = " + << JobD << endl; + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString ext = "2.root"; + + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("NN"); + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.30; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); + //fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + //tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + + //------------------------------------------- + // Display the histograms + // + + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + //------------------------------------------- + // fit alpha distribution to get the number of excess events + // + + MHHillasSrc* hillasSrc = + (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); + TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + + MHOnSubtraction onsub; + onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); + //------------------------------------------- + + DeleteBinnings(&pliston); + + gLog << "Macro CT1Analysis : End of Job D" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job E_EST_UP + //============ + + // - read MC2.root file + // - select g/h separation method XX + // - optimize energy estimator for events passing final cuts + // - write parameters of energy estimator onto file "energyest_XX.root" + // + // - read ON2.root and MC2.root files + // - update input root file with the estimated energy + // (ON_XX2.root, MC_XX2.root) + + + if (JobE_EST_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " + << JobE_EST_UP << ", " << WESTUP << endl; + + + TString typeON = "ON"; + TString typeOFF = "OFF"; + TString typeMC = "MC"; + TString ext = "2.root"; + TString extout = "3.root"; + + //------------------------------ + // name of MC file to be used for optimizing the energy estimator + TString filenameOpt(outPath); + filenameOpt += typeMC; + filenameOpt += ext; + gLog << "filenameOpt = " << filenameOpt << endl; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("NN"); + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |alpha| and dist + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + // name of file containing the parameters of the energy estimator + TString energyParName(outPath); + energyParName += "energyest_"; + energyParName += XX; + energyParName += ".root"; + gLog << "energyParName = " << energyParName << endl; + + + //------------------------------ + // name of ON file to be updated + TString filenameON(outPath); + filenameON += typeON; + filenameON += ext; + gLog << "filenameON = " << filenameON << endl; + + // name of OFF file to be updated + TString filenameOFF(outPath); + filenameOFF += typeOFF; + filenameOFF += ext; + gLog << "filenameOFF = " << filenameOFF << endl; + + // name of MC file to be updated + TString filenameMC(outPath); + filenameMC += typeMC; + filenameMC += ext; + gLog << "filenameMC = " << filenameMC << endl; + + //------------------------------ + // name of updated ON file + TString filenameONup(outPath); + filenameONup += typeON; + filenameONup += "_"; + filenameONup += XX; + filenameONup += extout; + gLog << "filenameONup = " << filenameONup << endl; + + // name of updated OFF file + TString filenameOFFup(outPath); + filenameOFFup += typeOFF; + filenameOFFup += "_"; + filenameOFFup += XX; + filenameOFFup += extout; + gLog << "filenameOFFup = " << filenameOFFup << endl; + + // name of updated MC file + TString filenameMCup(outPath); + filenameMCup += typeMC; + filenameMCup += "_"; + filenameMCup += XX; + filenameMCup += extout; + gLog << "filenameMCup = " << filenameMCup << endl; + + //----------------------------------------------------------- + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //=========================================================== + // + // Optimization of energy estimator + // + + TString inpath(""); + TString outpath(""); + Int_t howMany = 2000; + CT1EEst(inpath, filenameOpt, outpath, energyParName, + fHilName, fHilNameSrc, fhadronnessName, + howMany, maxhadronness, maxalpha, maxdist); + + //----------------------------------------------------------- + // + // Read in parameters of energy estimator + // + gLog << "================================================================" + << endl; + gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" + << energyParName << "'" << endl; + TFile enparam(energyParName); + MMcEnergyEst mcest("MMcEnergyEst"); + mcest.Read("MMcEnergyEst"); + enparam.Close(); + + TArrayD parA(5); + TArrayD parB(7); + for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameOFF); + read.DisableAutoScheme(); + + //--------------------------- + // calculate estimated energy + + MEnergyEstParam eest2(fHilName); + eest2.Add(fHilNameSrc); + + eest2.SetCoeffA(parA); + eest2.SetCoeffB(parB); + //....................................................................... - // calculate hadronnes for method of RANDOM FOREST - - TString hadRFName = "HadRF"; - MRanForestCalc rfcalc; - rfcalc.SetHadronnessName(hadRFName); - - //....................................................................... - // calculation of hadroness for the supercuts - // (=0.25 if fullfilled, =0.75 otherwise) - - TString hadSCName = "HadSC"; - MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); - sccalc.SetHadronnessName(hadSCName); - - - //....................................................................... - - //MWriteRootFile write(outNameImage, "UPDATE"); - MWriteRootFile write(outNameImage, "RECREATE"); + + MWriteRootFile write(filenameOFFup); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -1769,678 +1951,74 @@ write.AddContainer("MNewImagePar", "Events"); + //write.AddContainer("HadNN", "Events"); + write.AddContainer("HadSC", "Events"); write.AddContainer("HadRF", "Events"); - write.AddContainer("HadSC", "Events"); + + write.AddContainer("MEnergyEst", "Events"); //----------------------------------------------------------------- + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + MContinue contfinal(&selfinal); + + + //***************************** + // entries in MParList + + plistoff.AddToList(&tlistoff); + InitBinnings(&plistoff); + + + //***************************** + // entries in MTaskList + + tlistoff.AddToList(&read); + tlistoff.AddToList(&eest2); + tlistoff.AddToList(&write); + tlistoff.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plistoff); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tlistoff.PrintStatistics(0, kTRUE); + DeleteBinnings(&plistoff); + + //--------------------------------------------------- + + //--------------------------------------------------- + // Update ON data + // + gLog << "============================================================" + << endl; + gLog << "Macro CT1Analysis.C : update file '" << filenameON + << "'" << endl; + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes MGeomCam *fGeom = (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - - MFCT1SelFinal selfinalgh(fHilNameSrc); - selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadRFName); - selfinalgh.SetName("SelFinalgh"); - MContinue contfinalgh(&selfinalgh); - contfinalgh.SetName("ContSelFinalgh"); - - MFillH fillranfor("MHRanForest"); - fillranfor.SetName("HRanForest"); - - MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); - fillhadrf.SetName("HhadRF"); - MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); - fillhadsc.SetName("HhadSC"); - - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadRFName); - selfinal.SetName("SelFinal"); - MContinue contfinal(&selfinal); - contfinal.SetName("ContSelFinal"); - - - MFillH hfill1("MHHillas", fHilName); - hfill1.SetName("HHillas"); - - MFillH hfill2("MHStarMap", fHilName); - hfill2.SetName("HStarMap"); - - MFillH hfill3("MHHillasExt", fHilNameSrc); - hfill3.SetName("HHillasExt"); - - MFillH hfill4("MHHillasSrc", fHilNameSrc); - hfill4.SetName("HHillasSrc"); - - MFillH hfill5("MHNewImagePar", fImgParName); - hfill5.SetName("HNewImagePar"); - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - pliston.AddToList(fRanForest); - pliston.AddToList(fRanTree); - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - - tliston.AddToList(&rfcalc); - tliston.AddToList(&sccalc); - tliston.AddToList(&fillranfor); - tliston.AddToList(&fillhadrf); - tliston.AddToList(&fillhadsc); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - tliston.AddToList(&hfill5); - - tliston.AddToList(&write); - - tliston.AddToList(&contfinalgh); - tliston.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - pliston.FindObject("MHRanForest")->DrawClone(); - pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); - pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); - - pliston.FindObject("MHHillas")->DrawClone(); - pliston.FindObject("MHHillasExt")->DrawClone(); - pliston.FindObject("MHHillasSrc")->DrawClone(); - pliston.FindObject("MHNewImagePar")->DrawClone(); - pliston.FindObject("MHStarMap")->DrawClone(); - - DeleteBinnings(&pliston); - } - - gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - - - - //--------------------------------------------------------------------- - // Job C - //====== - - // - read ON1 and MC1 data files - // which should have been updated to contain the hadronnesses - // for the method of NEAREST NEIGHBORS and for the SUOERCUTS - // - produce Neyman-Pearson plots - - if (JobC) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job C" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobC = " << JobC << endl; - - - // name of input data file - TString filenameData = outPath; - filenameData += "OFF"; - filenameData += "2.root"; - gLog << "filenameData = " << filenameData << endl; - - // name of input MC file - TString filenameMC = outPath; - filenameMC += "MC"; - filenameMC += "2.root"; - gLog << "filenameMC = " << filenameMC << endl; - - - //----------------------------------------------------------------- - - MTaskList tliston; - MParList pliston; - - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - //------------------------------------------- // create the tasks which should be executed // - MReadMarsFile read("Events", filenameMC); - read.AddFile(filenameData); - read.DisableAutoScheme(); - - - //....................................................................... - // names of hadronness containers - - TString hadNNName = "HadNN"; - TString hadSCName = "HadSC"; - TString hadRFName = "HadRF"; - - //....................................................................... - - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - - MFCT1SelFinal selfinalgh(fHilNameSrc); - selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(hadRFName); - selfinalgh.SetName("SelFinalgh"); - MContinue contfinalgh(&selfinalgh); - contfinalgh.SetName("ContSelFinalgh"); - - //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); - //fillhadnn.SetName("HhadNN"); - MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); - fillhadsc.SetName("HhadSC"); - MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); - fillhadrf.SetName("HhadRF"); - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(hadRFName); - selfinal.SetName("SelFinal"); - MContinue contfinal(&selfinal); - contfinal.SetName("ContSelFinal"); - - - MFillH hfill1("MHHillas", fHilName); - hfill1.SetName("HHillas"); - - MFillH hfill2("MHStarMap", fHilName); - hfill2.SetName("HStarMap"); - - MFillH hfill3("MHHillasExt", fHilNameSrc); - hfill3.SetName("HHillasExt"); - - MFillH hfill4("MHHillasSrc", fHilNameSrc); - hfill4.SetName("HHillasSrc"); - - MFillH hfill5("MHNewImagePar", fImgParName); - hfill5.SetName("HNewImagePar"); - - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - - //tliston.AddToList(&fillhadnn); - tliston.AddToList(&fillhadsc); - tliston.AddToList(&fillhadrf); - - tliston.AddToList(&contfinalgh); - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - tliston.AddToList(&hfill5); - - tliston.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - //Int_t maxevents = 35000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - - //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); - pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); - pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); - - pliston.FindObject("MHHillas")->DrawClone(); - pliston.FindObject("MHHillasExt")->DrawClone(); - pliston.FindObject("MHHillasSrc")->DrawClone(); - pliston.FindObject("MHNewImagePar")->DrawClone(); - pliston.FindObject("MHStarMap")->DrawClone(); - - DeleteBinnings(&pliston); - - gLog << "Macro CT1Analysis : End of Job C" << endl; - gLog << "===================================================" << endl; - } - - - //--------------------------------------------------------------------- - // Job D - //====== - - // - select g/h separation method XX - // - read ON2 (or MC2) root file - // - apply cuts in hadronness - // - make plots - - - if (JobD) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job D" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobD = " - << JobD << endl; - - // type of data to be analysed - TString typeData = "ON"; - //TString typeData = "OFF"; - //TString typeData = "MC"; - gLog << "typeData = " << typeData << endl; - - TString ext = "2.root"; - - - //------------------------------ - // selection of g/h separation method - // and definition of final selections - - //TString XX("NN"); - //TString XX("SC"); - TString XX("RF"); - TString fhadronnessName("Had"); - fhadronnessName += XX; - gLog << "fhadronnessName = " << fhadronnessName << endl; - - // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.30; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - gLog << "Maximum values of hadronness, |ALPHA| and DIST = " - << maxhadronness << ", " << maxalpha << ", " - << maxdist << endl; - - - //------------------------------ - // name of data file to be analysed - TString filenameData(outPath); - filenameData += typeData; - filenameData += ext; - gLog << "filenameData = " << filenameData << endl; - - - - //************************************************************************* - // - // Analyse the data - // - - MTaskList tliston; - MParList pliston; - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - - //----------------------------------------------------------------- - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - MFCT1SelFinal selfinalgh(fHilNameSrc); - selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); - selfinalgh.SetHadronnessName(fhadronnessName); - selfinalgh.SetName("SelFinalgh"); - MContinue contfinalgh(&selfinalgh); - contfinalgh.SetName("ContSelFinalgh"); - - //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); - //fillhadnn.SetName("HhadNN"); - MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); - fillhadsc.SetName("HhadSC"); - MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); - fillhadrf.SetName("HhadRF"); - - - MFillH hfill1("MHHillas", fHilName); - hfill1.SetName("HHillas"); - - MFillH hfill2("MHStarMap", fHilName); - hfill2.SetName("HStarMap"); - - MFillH hfill3("MHHillasExt", fHilNameSrc); - hfill3.SetName("HHillasExt"); - - MFillH hfill4("MHHillasSrc", fHilNameSrc); - hfill4.SetName("HHillasSrc"); - - MFillH hfill5("MHNewImagePar", fImgParName); - hfill5.SetName("HNewImagePar"); - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(fhadronnessName); - selfinal.SetName("SelFinal"); - MContinue contfinal(&selfinal); - contfinal.SetName("ContSelFinal"); - - - //***************************** - // entries in MParList - - pliston.AddToList(&tliston); - InitBinnings(&pliston); - - - //***************************** - // entries in MTaskList - - tliston.AddToList(&read); - - tliston.AddToList(&contfinalgh); - - //tliston.AddToList(&fillhadnn); - tliston.AddToList(&fillhadsc); - tliston.AddToList(&fillhadrf); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - tliston.AddToList(&hfill5); - - tliston.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - //Int_t maxevents = 10000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - - //------------------------------------------- - // Display the histograms - // - - //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); - pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); - pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); - - pliston.FindObject("MHHillas")->DrawClone(); - pliston.FindObject("MHHillasExt")->DrawClone(); - pliston.FindObject("MHHillasSrc")->DrawClone(); - pliston.FindObject("MHNewImagePar")->DrawClone(); - pliston.FindObject("MHStarMap")->DrawClone(); - - //------------------------------------------- - // fit alpha distribution to get the number of excess events - // - - MHHillasSrc* hillasSrc = - (MHHillasSrc*)(pliston.FindObject("MHHillasSrc")); - TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); - - MHOnSubtraction onsub; - onsub.Calc(alphaHist, &pliston, kTRUE, 13.1); - //------------------------------------------- - - DeleteBinnings(&pliston); - - gLog << "Macro CT1Analysis : End of Job D" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - - - //--------------------------------------------------------------------- - // Job E_EST_UP - //============ - - // - read MC2.root file - // - select g/h separation method XX - // - optimize energy estimator for events passing final cuts - // - write parameters of energy estimator onto file "energyest_XX.root" - // - // - read ON2.root and MC2.root files - // - update input root file with the estimated energy - // (ON_XX2.root, MC_XX2.root) - - - if (JobE_EST_UP) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " - << JobE_EST_UP << ", " << WESTUP << endl; - - - TString typeON = "ON"; - TString typeOFF = "OFF"; - TString typeMC = "MC"; - TString ext = "2.root"; - TString extout = "3.root"; - - //------------------------------ - // name of MC file to be used for optimizing the energy estimator - TString filenameOpt(outPath); - filenameOpt += typeMC; - filenameOpt += ext; - gLog << "filenameOpt = " << filenameOpt << endl; - - //------------------------------ - // selection of g/h separation method - // and definition of final selections - - //TString XX("NN"); - //TString XX("SC"); - TString XX("RF"); - TString fhadronnessName("Had"); - fhadronnessName += XX; - gLog << "fhadronnessName = " << fhadronnessName << endl; - - // maximum values of the hadronness, |alpha| and dist - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - gLog << "Maximum values of hadronness, |ALPHA| and DIST = " - << maxhadronness << ", " << maxalpha << ", " - << maxdist << endl; - - // name of file containing the parameters of the energy estimator - TString energyParName(outPath); - energyParName += "energyest_"; - energyParName += XX; - energyParName += ".root"; - gLog << "energyParName = " << energyParName << endl; - - - //------------------------------ - // name of ON file to be updated - TString filenameON(outPath); - filenameON += typeON; - filenameON += ext; - gLog << "filenameON = " << filenameON << endl; - - // name of OFF file to be updated - TString filenameOFF(outPath); - filenameOFF += typeOFF; - filenameOFF += ext; - gLog << "filenameOFF = " << filenameOFF << endl; - - // name of MC file to be updated - TString filenameMC(outPath); - filenameMC += typeMC; - filenameMC += ext; - gLog << "filenameMC = " << filenameMC << endl; - - //------------------------------ - // name of updated ON file - TString filenameONup(outPath); - filenameONup += typeON; - filenameONup += "_"; - filenameONup += XX; - filenameONup += extout; - gLog << "filenameONup = " << filenameONup << endl; - - // name of updated OFF file - TString filenameOFFup(outPath); - filenameOFFup += typeOFF; - filenameOFFup += "_"; - filenameOFFup += XX; - filenameOFFup += extout; - gLog << "filenameOFFup = " << filenameOFFup << endl; - - // name of updated MC file - TString filenameMCup(outPath); - filenameMCup += typeMC; - filenameMCup += "_"; - filenameMCup += XX; - filenameMCup += extout; - gLog << "filenameMCup = " << filenameMCup << endl; - - //----------------------------------------------------------- - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //=========================================================== - // - // Optimization of energy estimator - // - - TString inpath(""); - TString outpath(""); - Int_t howMany = 2000; - CT1EEst(inpath, filenameOpt, outpath, energyParName, - fHilName, fHilNameSrc, fhadronnessName, - howMany, maxhadronness, maxalpha, maxdist); - - //----------------------------------------------------------- - // - // Read in parameters of energy estimator - // - gLog << "================================================================" - << endl; - gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" - << energyParName << "'" << endl; - TFile enparam(energyParName); - MMcEnergyEst mcest("MMcEnergyEst"); - mcest.Read("MMcEnergyEst"); - enparam.Close(); - - TArrayD parA(5); - TArrayD parB(7); - for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameOFF); + MReadMarsFile read("Events", filenameON); read.DisableAutoScheme(); @@ -2456,5 +2034,5 @@ //....................................................................... - MWriteRootFile write(filenameOFFup); + MWriteRootFile write(filenameONup); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -2486,6 +2064,6 @@ // entries in MParList - plistoff.AddToList(&tlistoff); - InitBinnings(&plistoff); + pliston.AddToList(&tliston); + InitBinnings(&pliston); @@ -2493,8 +2071,8 @@ // entries in MTaskList - tlistoff.AddToList(&read); - tlistoff.AddToList(&eest2); - tlistoff.AddToList(&write); - tlistoff.AddToList(&contfinal); + tliston.AddToList(&read); + tliston.AddToList(&eest2); + tliston.AddToList(&write); + tliston.AddToList(&contfinal); //***************************** @@ -2505,5 +2083,5 @@ MProgressBar bar; MEvtLoop evtloop; - evtloop.SetParList(&plistoff); + evtloop.SetParList(&pliston); evtloop.SetProgressBar(&bar); @@ -2513,24 +2091,19 @@ return; - tlistoff.PrintStatistics(0, kTRUE); - DeleteBinnings(&plistoff); + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); //--------------------------------------------------- //--------------------------------------------------- - // Update ON data + // Update MC data // gLog << "============================================================" << endl; - gLog << "Macro CT1Analysis.C : update file '" << filenameON + gLog << "Macro CT1Analysis.C : update file '" << filenameMC << "'" << endl; - MTaskList tliston; - MParList pliston; - - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + MTaskList tlistmc; + MParList plistmc; //------------------------------------------- @@ -2538,5 +2111,5 @@ // - MReadMarsFile read("Events", filenameON); + MReadMarsFile read("Events", filenameMC); read.DisableAutoScheme(); @@ -2552,5 +2125,5 @@ //....................................................................... - MWriteRootFile write(filenameONup); + MWriteRootFile write(filenameMCup); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -2582,6 +2155,6 @@ // entries in MParList - pliston.AddToList(&tliston); - InitBinnings(&pliston); + plistmc.AddToList(&tlistmc); + InitBinnings(&plistmc); @@ -2589,8 +2162,8 @@ // entries in MTaskList - tliston.AddToList(&read); - tliston.AddToList(&eest2); - tliston.AddToList(&write); - tliston.AddToList(&contfinal); + tlistmc.AddToList(&read); + tlistmc.AddToList(&eest2); + tlistmc.AddToList(&write); + tlistmc.AddToList(&contfinal); //***************************** @@ -2601,5 +2174,5 @@ MProgressBar bar; MEvtLoop evtloop; - evtloop.SetParList(&pliston); + evtloop.SetParList(&plistmc); evtloop.SetProgressBar(&bar); @@ -2609,19 +2182,203 @@ return; - tliston.PrintStatistics(0, kTRUE); - DeleteBinnings(&pliston); - - //--------------------------------------------------- - - //--------------------------------------------------- - // Update MC data - // - gLog << "============================================================" - << endl; - gLog << "Macro CT1Analysis.C : update file '" << filenameMC - << "'" << endl; - - MTaskList tlistmc; - MParList plistmc; + tlistmc.PrintStatistics(0, kTRUE); + DeleteBinnings(&plistmc); + + + //========== end update ============================================ + } + + enparam.Close(); + + gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job F_XX + //========= + + // - select g/h separation method XX + // - read MC_XX2.root file +   // - calculate eff. collection area + // - read ON_XX2.root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts (ON_XX3.root)) + + + if (JobF_XX) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobF_XX, WXX = " + << JobF_XX << ", " << WXX << endl; + + // type of data to be analysed + //TString typeData = "ON"; + //TString typeData = "OFF"; + TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString typeMC = "MC"; + TString ext = "3.root"; + TString extout = "4.root"; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("NN"); + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of MC file to be used for calculating the eff. collection areas + TString filenameArea(outPath); + filenameArea += typeMC; + filenameArea += "_"; + filenameArea += XX; + filenameArea += ext; + gLog << "filenameArea = " << filenameArea << endl; + + //------------------------------ + // name of file containing the eff. collection areas + TString collareaName(outPath); + collareaName += "area_"; + collareaName += XX; + collareaName += ".root"; + gLog << "collareaName = " << collareaName << endl; + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += "_"; + filenameData += XX; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + //------------------------------ + // name of output data file (after the final cuts) + TString filenameDataout(outPath); + filenameDataout += typeData; + filenameDataout += "_"; + filenameDataout += XX; + filenameDataout += extout; + gLog << "filenameDataout = " << filenameDataout << endl; + + + //==================================================================== + gLog << "-----------------------------------------------" << endl; + gLog << "Start calculation of effective collection areas" << endl; + MParList parlist; + MTaskList tasklist; + + //--------------------------------------- + // Setup the tasks to be executed + // + MReadMarsFile reader("Events", filenameArea); + reader.DisableAutoScheme(); + + MFCT1SelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); + + MContinue conthadrons(&cuthadrons); + + //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); + //MHMcCT1CollectionArea* collarea; + + MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); + filler.SetName("CollectionArea"); + + //******************************** + // entries in MParList + + parlist.AddToList(&tasklist); + InitBinnings(&parlist); + //parlist.AddToList(collarea); + + //******************************** + // entries in MTaskList + + tasklist.AddToList(&reader); + tasklist.AddToList(&conthadrons); + tasklist.AddToList(&filler); + + //******************************** + + //----------------------------------------- + // Execute event loop + // + MEvtLoop magic; + magic.SetParList(&parlist); + + MProgressBar bar; + magic.SetProgressBar(&bar); + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(0, kTRUE); + + // Calculate effective collection areas + // and display the histograms + // + + MHMcCT1CollectionArea *collarea = + (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); + + collarea->CalcEfficiency(); + collarea->DrawClone("lego"); + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + f.Close(); + + //delete collarea; + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; //------------------------------------------- @@ -2629,19 +2386,11 @@ // - MReadMarsFile read("Events", filenameMC); + MReadMarsFile read("Events", filenameData); read.DisableAutoScheme(); - //--------------------------- - // calculate estimated energy - - MEnergyEstParam eest2(fHilName); - eest2.Add(fHilNameSrc); - - eest2.SetCoeffA(parA); - eest2.SetCoeffB(parB); - //....................................................................... - MWriteRootFile write(filenameMCup); + + MWriteRootFile write(filenameDataout); write.AddContainer("MRawRunHeader", "RunHeaders"); @@ -2662,271 +2411,4 @@ write.AddContainer("MEnergyEst", "Events"); - //----------------------------------------------------------------- - - MFCT1SelFinal selfinal(fHilNameSrc); - selfinal.SetCuts(maxhadronness, maxalpha, maxdist); - selfinal.SetHadronnessName(fhadronnessName); - MContinue contfinal(&selfinal); - - - //***************************** - // entries in MParList - - plistmc.AddToList(&tlistmc); - InitBinnings(&plistmc); - - - //***************************** - // entries in MTaskList - - tlistmc.AddToList(&read); - tlistmc.AddToList(&eest2); - tlistmc.AddToList(&write); - tlistmc.AddToList(&contfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&plistmc); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = -1; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tlistmc.PrintStatistics(0, kTRUE); - DeleteBinnings(&plistmc); - - - //========== end update ============================================ - } - - enparam.Close(); - - gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl; - gLog << "=======================================================" << endl; - } - //--------------------------------------------------------------------- - - - //--------------------------------------------------------------------- - // Job F_XX - //========= - - // - select g/h separation method XX - // - read MC_XX2.root file -   // - calculate eff. collection area - // - read ON_XX2.root file - // - apply final cuts - // - calculate flux - // - write root file for ON data after final cuts (ON_XX3.root)) - - - if (JobF_XX) - { - gLog << "=====================================================" << endl; - gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; - - gLog << "" << endl; - gLog << "Macro CT1Analysis : JobF_XX, WXX = " - << JobF_XX << ", " << WXX << endl; - - // type of data to be analysed - //TString typeData = "ON"; - //TString typeData = "OFF"; - TString typeData = "MC"; - gLog << "typeData = " << typeData << endl; - - TString typeMC = "MC"; - TString ext = "3.root"; - TString extout = "4.root"; - - //------------------------------ - // selection of g/h separation method - // and definition of final selections - - //TString XX("NN"); - //TString XX("SC"); - TString XX("RF"); - TString fhadronnessName("Had"); - fhadronnessName += XX; - gLog << "fhadronnessName = " << fhadronnessName << endl; - - // maximum values of the hadronness, |ALPHA| and DIST - Float_t maxhadronness = 0.40; - Float_t maxalpha = 20.0; - Float_t maxdist = 10.0; - gLog << "Maximum values of hadronness, |ALPHA| and DIST = " - << maxhadronness << ", " << maxalpha << ", " - << maxdist << endl; - - - //------------------------------ - // name of MC file to be used for calculating the eff. collection areas - TString filenameArea(outPath); - filenameArea += typeMC; - filenameArea += "_"; - filenameArea += XX; - filenameArea += ext; - gLog << "filenameArea = " << filenameArea << endl; - - //------------------------------ - // name of file containing the eff. collection areas - TString collareaName(outPath); - collareaName += "area_"; - collareaName += XX; - collareaName += ".root"; - gLog << "collareaName = " << collareaName << endl; - - //------------------------------ - // name of data file to be analysed - TString filenameData(outPath); - filenameData += typeData; - filenameData += "_"; - filenameData += XX; - filenameData += ext; - gLog << "filenameData = " << filenameData << endl; - - //------------------------------ - // name of output data file (after the final cuts) - TString filenameDataout(outPath); - filenameDataout += typeData; - filenameDataout += "_"; - filenameDataout += XX; - filenameDataout += extout; - gLog << "filenameDataout = " << filenameDataout << endl; - - - //==================================================================== - gLog << "-----------------------------------------------" << endl; - gLog << "Start calculation of effective collection areas" << endl; - MParList parlist; - MTaskList tasklist; - - //--------------------------------------- - // Setup the tasks to be executed - // - MReadMarsFile reader("Events", filenameArea); - reader.DisableAutoScheme(); - - MFCT1SelFinal cuthadrons; - cuthadrons.SetHadronnessName(fhadronnessName); - cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); - - MContinue conthadrons(&cuthadrons); - - //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); - //MHMcCT1CollectionArea* collarea; - - MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); - filler.SetName("CollectionArea"); - - //******************************** - // entries in MParList - - parlist.AddToList(&tasklist); - InitBinnings(&parlist); - //parlist.AddToList(collarea); - - //******************************** - // entries in MTaskList - - tasklist.AddToList(&reader); - tasklist.AddToList(&conthadrons); - tasklist.AddToList(&filler); - - //******************************** - - //----------------------------------------- - // Execute event loop - // - MEvtLoop magic; - magic.SetParList(&parlist); - - MProgressBar bar; - magic.SetProgressBar(&bar); - if (!magic.Eventloop()) - return; - - tasklist.PrintStatistics(0, kTRUE); - - // Calculate effective collection areas - // and display the histograms - // - - MHMcCT1CollectionArea *collarea = - (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); - - collarea->CalcEfficiency(); - collarea->DrawClone("lego"); - - //--------------------------------------------- - // Write histograms to a file - // - - TFile f(collareaName, "RECREATE"); - collarea->GetHist()->Write(); - collarea->GetHAll()->Write(); - collarea->GetHSel()->Write(); - f.Close(); - - //delete collarea; - - gLog << "Calculation of effective collection areas done" << endl; - gLog << "-----------------------------------------------" << endl; - //------------------------------------------------------------------ - - - //************************************************************************* - // - // Analyse the data - // - - MTaskList tliston; - MParList pliston; - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - TString fHilName = "MHillas"; - TString fHilNameExt = "MHillasExt"; - TString fHilNameSrc = "MHillasSrc"; - TString fImgParName = "MNewImagePar"; - - //------------------------------------------- - // create the tasks which should be executed - // - - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - //....................................................................... - - - MWriteRootFile write(filenameDataout); - - write.AddContainer("MRawRunHeader", "RunHeaders"); - write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("ThetaOrig", "Events"); - write.AddContainer("MSrcPosCam", "Events"); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("MHillas", "Events"); - write.AddContainer("MHillasExt", "Events"); - write.AddContainer("MHillasSrc", "Events"); - write.AddContainer("MNewImagePar", "Events"); - - //write.AddContainer("HadNN", "Events"); - write.AddContainer("HadSC", "Events"); - write.AddContainer("HadRF", "Events"); - - write.AddContainer("MEnergyEst", "Events"); - //----------------------------------------------------------------- Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2436) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2437) @@ -42,5 +42,5 @@ #pragma link C++ class MCT1PadSchweizer+; -//#pragma link C++ class MCT1PadONOFF+; +#pragma link C++ class MCT1PadONOFF+; #pragma link C++ class MCT1PointingCorrCalc+; @@ -59,4 +59,5 @@ #pragma link C++ class MMcPedestalRead+; + #endif @@ -70,2 +71,4 @@ + + Index: trunk/MagicSoft/Mars/manalysis/MCT1PadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1PadONOFF.cc (revision 2436) +++ trunk/MagicSoft/Mars/manalysis/MCT1PadONOFF.cc (revision 2437) @@ -200,5 +200,5 @@ { edgesx[i] = ax->GetBinLowEdge(i+1); - //*fLog << "i, theta low edge = " << i << ", " << edgesx[i] << endl; + // *fLog << "i, theta low edge = " << i << ", " << edgesx[i] << endl; } MBinning binth; @@ -213,5 +213,5 @@ { edgesy[i] = ay->GetBinLowEdge(i+1); - //*fLog << "i, sigmabar low edge = " << i << ", " << edgesy[i] << endl; + // *fLog << "i, sigmabar low edge = " << i << ", " << edgesy[i] << endl; } MBinning binsg; @@ -989,5 +989,5 @@ Int_t MCT1PadONOFF::Process() { - //*fLog << "Entry MCT1PadONOFF::Process();" << endl; + // *fLog << "Entry MCT1PadONOFF::Process();" << endl; //------------------------------------------------ @@ -1014,5 +1014,5 @@ continue; - fEvt->AddPixel(i, 0.0, (*fPed)[i].GetMeanRms()); + fEvt->AddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); } @@ -1026,5 +1026,5 @@ //fSigmabar->Calc(*fCam, *fPed, *fEvt); - //*fLog << "before padding : " << endl; + // *fLog << "before padding : " << endl; //fSigmabar->Print(""); @@ -1055,5 +1055,5 @@ if (sigbarold > 0) { - //*fLog << "MCT1PadONOFF::Process(); sigmabar of event to be padded is > 0 : " + // *fLog << "MCT1PadONOFF::Process(); sigmabar of event to be padded is > 0 : " // << sigbarold << ". Stop event loop " << endl; // input data should have sigmabar = 0; stop event loop @@ -1066,10 +1066,10 @@ const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - //*fLog << "theta = " << theta << endl; + // *fLog << "theta = " << theta << endl; Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) { - //*fLog << "MCT1PadONOFF::Process(); binNumber out of range : theta, binTheta = " + // *fLog << "MCT1PadONOFF::Process(); binNumber out of range : theta, binTheta = " // << theta << ", " << binTheta << "; Skip event " << endl; // event cannot be padded; skip event @@ -1095,5 +1095,5 @@ nSigma = hn->GetBinContent(binSigma); - //*fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = " + // *fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = " // << theta << ", " << sigbarold << ", " << binTheta << ", " // << binSigma << ", " << nTheta << ", " << nSigma << endl; @@ -1111,5 +1111,5 @@ nSigma = hn->GetBinContent(binSigma); - //*fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = " + // *fLog << "Theta, sigbarold, binTheta, binSigma, nTheta, nSigma = " // << theta << ", " << sigbarold << ", " << binTheta << ", " // << binSigma << ", " << nTheta << ", " << nSigma << endl; @@ -1141,5 +1141,5 @@ if ( nblind->Integral() == 0.0 ) { - //*fLog << "MCT1PadONOFF::Process(); projection for Theta bin " + // *fLog << "MCT1PadONOFF::Process(); projection for Theta bin " // << binTheta << " has no entries; Skip event " << endl; // event cannot be padded; skip event @@ -1156,6 +1156,6 @@ delete nblind; -#warn Code commented out due no compilation errors on Alpha OSF (tgb) -/* + //warn Code commented out due no compilation errors on Alpha OSF (tgb) + // throw the Id of numBlind different pixels in this event TH1D *hblind; @@ -1185,5 +1185,5 @@ fBlinds->SetPixelBlind(idBlind); - //*fLog << "idBlind = " << idBlind << endl; + // *fLog << "idBlind = " << idBlind << endl; } fBlinds->SetReadyToSave(); @@ -1192,5 +1192,5 @@ } -*/ + //****************************************************************** // has the event to be padded ? @@ -1199,5 +1199,5 @@ Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold); - //*fLog << "binSig, sigbarold = " << binSig << ", " << sigbarold << endl; + // *fLog << "binSig, sigbarold = " << binSig << ", " << sigbarold << endl; Double_t prob; @@ -1219,5 +1219,5 @@ prob = 0.0; - //*fLog << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + // *fLog << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma // << ", " << prob << endl; } @@ -1238,5 +1238,5 @@ delete hpad; // event should not be padded - //*fLog << "event is not padded" << endl; + // *fLog << "event is not padded" << endl; rc = 8; @@ -1245,5 +1245,5 @@ } // event should be padded - //*fLog << "event is padded" << endl; + // *fLog << "event is padded" << endl; @@ -1266,5 +1266,5 @@ sigmabar = hpad->GetRandom(); - //*fLog << "sigmabar = " << sigmabar << endl; + // *fLog << "sigmabar = " << sigmabar << endl; delete hpad; @@ -1299,5 +1299,5 @@ { sigmabar = hsigma->GetRandom(); - //*fLog << "Theta, bin number = " << theta << ", " << binTheta + // *fLog << "Theta, bin number = " << theta << ", " << binTheta // << ", sigmabar = " << sigmabar << endl } @@ -1309,5 +1309,5 @@ //------------------------------------------- - //*fLog << "MCT1PadONOFF::Process(); sigbarold, sigmabar = " + // *fLog << "MCT1PadONOFF::Process(); sigbarold, sigmabar = " // << sigbarold << ", "<< sigmabar << endl; @@ -1367,5 +1367,5 @@ elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); - //*fLog << "elNoise2 = " << elNoise2 << endl; + // *fLog << "elNoise2 = " << elNoise2 << endl; lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photons] @@ -1405,5 +1405,5 @@ MPedestalPix &ppix = (*fPed)[j]; - Double_t oldsigma = ppix.GetMeanRms(); + Double_t oldsigma = ppix.GetPedestalRms(); Double_t oldsigma2 = oldsigma*oldsigma; @@ -1463,5 +1463,5 @@ if (!ok) { - //*fLog << "theta, j, count, sigmabar, diff = " << theta << ", " + // *fLog << "theta, j, count, sigmabar, diff = " << theta << ", " // << j << ", " << count << ", " << sigmabar << ", " // << diff << endl; @@ -1524,5 +1524,5 @@ if (!ok) { - //*fLog << "theta, j, count, sigmabar, sig = " << theta << ", " + // *fLog << "theta, j, count, sigmabar, sig = " << theta << ", " // << j << ", " << count << ", " << sigmabar << ", " // << sig << endl; @@ -1602,5 +1602,5 @@ Double_t newsigma = sqrt( oldsigma2 + addSig2 ); - ppix.SetMeanRms( newsigma ); + ppix.SetPedestalRms( newsigma ); fHSigmaPedestal->Fill( oldsigma, newsigma ); @@ -1611,5 +1611,5 @@ //fSigmabar->Calc(*fCam, *fPed, *fEvt); - //*fLog << "after padding : " << endl; + // *fLog << "after padding : " << endl; //fSigmabar->Print(""); Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2436) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2437) @@ -66,6 +66,7 @@ MMinuitInterface.cc \ MFiltercutsCalc.cc \ + MCT1PadONOFF.cc \ MMcPedestalRead.cc -# MCT1PadONOFF.cc \ + SRCS = $(SRCFILES)