Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2662) +++ trunk/MagicSoft/Mars/Changelog (revision 2663) @@ -4,4 +4,52 @@ -*-*- END OF LINE -*-*- + + 2003/12/12: Sabrina Stark + + * mfilter/MFSelBasic.[h,cc] + - new class for MAGIC analysis, corresponding to MFCT1SelBasic + + * mfilter/MFSelStandard.[h,cc] + - new class for MAGIC analysis, corresponding to MFCT1SelStandard + + * mfilter/MFSelFinal.[h,cc] + - new class for MAGIC analysis, corresponding to MFCT1SelFinal + + * mfilter/Makefile + /FilterLinkDef.h + - include MFSelBasic, MFSelStandard, MFSelFinal + + * manalysis/MPadONOFF.[h,cc] + - new class for MAGIC analysis, corresponding to MPadONOFF + + * manalysis/MSupercuts.[h,cc] + - new class for MAGIC analysis, corresponding to MCT1Supercuts + + * manalysis/MSupercutsCalc.[h,cc] + - new class for MAGIC analysis, corresponding to MCT1SupercutsCalc + + * manalysis/MFindSupercuts.[h,cc] + - new class for MAGIC analysis, corresponding to MCT1FindSupercuts + + * manalysis/MPedestalWorkaround.[h,cc] + - temporary task for MAGIC analysis to copy pedestal values and + RMS from MPedPhotCam to MPedestalCam + + * manalysis/Makefile + /AnalysisLinkDef.h + - include MPadONOFF, MSupercuts, MSupercutsCalc, MFindSupercuts, + MPedestalWorkaround + + * mhist/MHSupercuts.[h,cc] + - new class for MAGIC analysis, corresponding to MHCT1Supercuts + + * mhist/Makefile + /HistLinkDef.h + - include MHSupercuts + + * macros/ONOFFAnalysis.C: + - adaption of ONOFFCT1Analysis to apply on Magic data + + 2003/12/12: Markus Gaug Index: trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 2663) +++ trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 2663) @@ -0,0 +1,3342 @@ + + +//#include "MagicEgyEst.C" + + +void InitBinnings(MParList *plist) +{ + gLog << "InitBinnings" << endl; + + //-------------------------------------------- + MBinning *binse = new MBinning("BinningE"); + //binse->SetEdgesLog(30, 1.0e2, 1.0e5); + + //This is Daniel's binning in energy: + binse->SetEdgesLog(14, 296.296, 86497.6); + plist->AddToList(binse); + + //-------------------------------------------- + + MBinning *binssize = new MBinning("BinningSize"); + binssize->SetEdgesLog(50, 10, 1.0e5); + plist->AddToList(binssize); + + MBinning *binsdistc = new MBinning("BinningDist"); + binsdistc->SetEdges(50, 0, 1.4); + plist->AddToList(binsdistc); + + MBinning *binswidth = new MBinning("BinningWidth"); + binswidth->SetEdges(50, 0, 1.0); + plist->AddToList(binswidth); + + MBinning *binslength = new MBinning("BinningLength"); + binslength->SetEdges(50, 0, 1.0); + plist->AddToList(binslength); + + MBinning *binsalpha = new MBinning("BinningAlpha"); + binsalpha->SetEdges(100, -100, 100); + plist->AddToList(binsalpha); + + MBinning *binsasym = new MBinning("BinningAsym"); + binsasym->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsasym); + + MBinning *binsm3l = new MBinning("BinningM3Long"); + binsm3l->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsm3l); + + MBinning *binsm3t = new MBinning("BinningM3Trans"); + binsm3t->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsm3t); + + + //..... + MBinning *binsb = new MBinning("BinningSigmabar"); + binsb->SetEdges( 100, 0.0, 5.0); + plist->AddToList(binsb); + + MBinning *binth = new MBinning("BinningTheta"); + // this is Daniel's binning in theta + //Double_t yedge[8] = + // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99}; + // this is our binning + Double_t yedge[9] = + {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; + TArrayD yed; + yed.Set(9,yedge); + binth->SetEdges(yed); + plist->AddToList(binth); + + MBinning *bincosth = new MBinning("BinningCosTheta"); + Double_t zedge[9]; + for (Int_t i=0; i<9; i++) + { + zedge[8-i] = cos(yedge[i]/kRad2Deg); + } + TArrayD zed; + zed.Set(9,zedge); + bincosth->SetEdges(zed); + plist->AddToList(bincosth); + + MBinning *binsdiff = new MBinning("BinningDiffsigma2"); + binsdiff->SetEdges(100, -5.0, 20.0); + plist->AddToList(binsdiff); + + // robert ---------------------------------------------- + MBinning *binsalphaf = new MBinning("BinningAlphaFlux"); + binsalphaf->SetEdges(100, -100, 100); + plist->AddToList(binsalphaf); + + MBinning *binsdifftime = new MBinning("BinningTimeDiff"); + binsdifftime->SetEdges(50, 0, 10); + plist->AddToList(binsdifftime); + + MBinning *binstime = new MBinning("BinningTime"); + binstime->SetEdges(50, 44500, 61000); + plist->AddToList(binstime); +} + + +void DeleteBinnings(MParList *plist) +{ + gLog << "DeleteBinnings" << endl; + + TObject *bin; + + //-------------------------------------------- + bin = plist->FindObject("BinningE"); + if (bin) delete bin; + + //-------------------------------------------- + + bin = plist->FindObject("BinningSize"); + if (bin) delete bin; + + bin = plist->FindObject("BinningDist"); + if (bin) delete bin; + + bin = plist->FindObject("BinningWidth"); + if (bin) delete bin; + + bin = plist->FindObject("BinningLength"); + if (bin) delete bin; + + bin = plist->FindObject("BinningAlpha"); + if (bin) delete bin; + + bin = plist->FindObject("BinningAsym"); + if (bin) delete bin; + + bin = plist->FindObject("BinningM3Long"); + if (bin) delete bin; + + bin = plist->FindObject("BinningM3Trans"); + if (bin) delete bin; + + //..... + bin = plist->FindObject("BinningSigmabar"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTheta"); + if (bin) delete bin; + + bin = plist->FindObject("BinningCosTheta"); + if (bin) delete bin; + + bin = plist->FindObject("BinningDiffsigma2"); + if (bin) delete bin; + + + // robert ---------------------------------------------- + bin = plist->FindObject("BinningAlphaFlux"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTimeDiff"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTime"); + if (bin) delete bin; +} + + + +//************************************************************************ +void ONOFFAnalysis() +{ + gLog.SetNoColors(); + + if (gRandom) + delete gRandom; + gRandom = new TRandom3(0); + + //----------------------------------------------- + //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; + const char *offfile = "/data/MAGIC/rootdata/2003_11_29/20031128_032*_D_OffCrab1_E.root"; + + //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc"; + // const char *onfile = "~magican/ct1test/wittek/mkn421_00-01"; + const char *onfile = "/data/MAGIC/rootdata/2003_11_29/20031128_031*_D_Crab-Nebula_E.root"; + + const char *mcfile = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root"; + //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_4/gh/0/0/G_M0_00_0_550*.root"; + //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_5/gh/0/0/G_M0_00_0_550*.root"; + //----------------------------------------------- + + // path for input for Mars + TString inPath = "/data/MAGIC/scratch/wittek/"; + + // path for output from Mars + TString outPath = "/data/MAGIC/scratch/wittek/"; + + //----------------------------------------------- + + //TEnv env("macros/CT1env.rc"); + //Bool_t printEnv = kFALSE; + + //************************************************************************ + + // Job A : + // - produce MHSigmaTheta plots for ON and OFF data + // - write out (or read in) these MHSigmaTheta plots + // - read ON (or OFF or MC) data + // - pad the events; + // - write root file for ON (or OFF or MC) data (ON1.root, ...); + + Bool_t JobA = kTRUE; + Bool_t GPad = kFALSE; // generate padding histograms? + Bool_t WPad = kFALSE; // write out padding histograms ? + Bool_t RPad = kFALSE; // read in padding histograms ? + Bool_t Wout = kTRUE; // write out root file ON1.root + // (or OFF1.root or MC1.root)? + + + // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file + // - if CTrainRF = TRUE : create matrices of training events + // and root files of training and test events + // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas + // - if RTree = TRUE : read in trees, otherwise create trees + // - calculate hadroness for method of RANDOM FOREST + // - update the input files with the hadronesses (ON2.root, OFF2.root + // or MC2.root) + + Bool_t JobB_RF_UP = kFALSE; + Bool_t CTrainRF = kFALSE; // create matrices of training events + // and root files of training and test events + Bool_t RTrainRF = kFALSE; // read in matrices of training events + Bool_t RTree = kFALSE; // read in trees (otherwise grow them) + Bool_t WRF = kFALSE; // update input root file ? + + + // Job B_SC_UP : read ON2.root (or MC2.root) file + // - depending on WParSC : create (or read in) supercuts parameter values + // - calculate hadroness for the SUPERCUTS + // - update the input files with the hadroness (==>ON3.root or MC3.root) + + Bool_t JobB_SC_UP = kTRUE; + Bool_t CMatrix = kTRUE; // create training and test matrices + Bool_t RMatrix = kFALSE; // read training and test matrices from file + Bool_t WOptimize = kFALSE; // do optimization using the training sample + // and write supercuts parameter values + // onto the file parSCfile + Bool_t RTest = kFALSE; // test the supercuts using the test matrix + Bool_t WSC = kFALSE; // update input root file ? + + + // Job C: + // - read ON3.root and MC3.root files + // which should have been updated to contain the hadronnesses + // for the method of + // RF + // SUPERCUTS and + // - produce Neyman-Pearson plots + + Bool_t JobC = kFALSE; + + + // Job D : + // - select g/h separation method XX + // - read ON3 (or MC3) root file + // - apply cuts in hadronness + // - make plots + + Bool_t JobD = kFALSE; + + + + // Job E_XX : extended version of E_XX (including flux plots) + // - select g/h separation method XX + // - read MC root file +   // - calculate eff. collection area + // - optimize energy estimator + // - read ON root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts + + + Bool_t JobE_XX = kFALSE; + Bool_t CCollArea= kFALSE; // calculate eff. collection areas + Bool_t OEEst = kFALSE; // optimize energy estimator + Bool_t WEX = kFALSE; // update root file ? + Bool_t WRobert = kFALSE; // write out Robert's file ? + + + + //************************************************************************ + + + //--------------------------------------------------------------------- + // Job A + //========= + // read ON data file + + // - produce the 2D-histogram "sigmabar versus Theta" + // (SigmaTheta_ON.root) for ON data + // (to be used for the padding of the MC gamma data) + + // - write a file of ON events (ON1.root) + // (after the standard cuts, before the g/h separation) + // (to be used together with the corresponding MC gamma file (MC1.root) + // for the optimization of the g/h separation) + + + if (JobA) + { + gLog << "=====================================================" << endl; + gLog << "Macro MagicAnalysis : Start of Job A" << endl; + gLog << "" << endl; + gLog << "Macro MagicAnalysis : JobA, WPad, RPad, Wout = " + << (JobA ? "kTRUE" : "kFALSE") << ", " + << (WPad ? "kTRUE" : "kFALSE") << ", " + << (RPad ? "kTRUE" : "kFALSE") << ", " + << (Wout ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------------- + // names of ON and OFF files to be read + // for generating the histograms to be used in the padding + TString fileON = onfile; + TString fileOFF = offfile; + gLog << "fileON, fileOFF = " << fileON << ", " << fileOFF << endl; + + // name of file to conatin the histograms for the padding + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta"; + outNameSigTh += ".root"; + + //-------------------------------------------------- + // type of data to be padded + TString typeInput = "ON"; + //TString typeInput = "OFF"; + //TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + + // name of input root file + if (typeInput == "ON") + TString filenamein(onfile); + else if (typeInput == "OFF") + TString filenamein(offfile); + else if (typeInput == "MC") + TString filenamein(mcfile); + gLog << "data to be padded : " << filenamein << endl; + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "1.root"; + gLog << "padded data to be written onto : " << outNameImage << endl; + + //-------------------------------------------------- + + //************************************************************ + // generate histograms to be used in the padding + // + // read ON and OFF data files + // generate (or read in) the padding histograms for ON and OFF data + // and merge these histograms + + MPadONOFF pad; + pad.SetName("MPadONOFF"); + pad.SetPadFlag(1); + pad.SetDataType(typeInput); + + // generate the padding histograms + if (GPad) + { + gLog << "=====================================================" << endl; + gLog << "Start generating the padding histograms" << endl; + //----------------------------------------- + // ON events + + gLog << "-----------" << endl; + gLog << "ON events :" << endl; + gLog << "-----------" << endl; + + MTaskList tliston; + MParList pliston; + + MReadMarsFile readON("Events", fileON); + read.DisableAutoScheme(); + MCT1ReadPreProc readON(fileON); + + //MFSelBasic selthetaon; + //selthetaon.SetCuts(-100.0, 29.5, 35.5); + //MContinue contthetaon(&selthetaon); + + MBlindPixelCalc blindon; + blindon.SetUseBlindPixels(); + + MFSelBasic selbasicon; + MContinue contbasicon(&selbasicon); + + MHBlindPixels blindON("BlindPixelsON"); + MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels"); + fillblindON.SetName("FillBlind"); + + MSigmabarCalc sigbarcalcon; + + MHSigmaTheta sigthON("SigmaThetaON"); + MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); + fillsigthetaON.SetName("FillSigTheta"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + pliston.AddToList(&blindON); + pliston.AddToList(&sigthON); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&readON); + //tliston.AddToList(&contthetaon); + + tliston.AddToList(&blindon); + + tliston.AddToList(&contbasicon); + tliston.AddToList(&fillblindON); + tliston.AddToList(&sigbarcalcon); + tliston.AddToList(&fillsigthetaON); + + MProgressBar baron; + MEvtLoop evtloopon; + evtloopon.SetParList(&pliston); + evtloopon.SetProgressBar(&baron); + + Int_t maxeventson = -1; + //Int_t maxeventson = 10000; + if ( !evtloopon.Eventloop(maxeventson) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + blindON.DrawClone(); + sigthON.DrawClone(); + + // save the histograms for the padding + TH2D *sigthon = sigthON.GetSigmaTheta(); + TH3D *sigpixthon = sigthON.GetSigmaPixTheta(); + TH3D *diffpixthon = sigthON.GetDiffPixTheta(); + + TH2D *blindidthon = blindON.GetBlindId(); + TH2D *blindnthon = blindON.GetBlindN(); + + //----------------------------------------- + // OFF events + + gLog << "------------" << endl; + gLog << "OFF events :" << endl; + gLog << "------------" << endl; + + MTaskList tlistoff; + MParList plistoff; + + MReadMarsFile readOFF("Events", fileOFF); + read.DisableAutoScheme(); + // MCT1ReadPreProc readOFF(fileOFF); + + MFSelBasic selthetaoff; + selthetaoff.SetCuts(-100.0, 29.5, 35.5); + MContinue contthetaoff(&selthetaoff); + + MBlindPixelCalc blindoff; + blindoff.SetUseBlindPixels(); + + MFSelBasic selbasicoff; + MContinue contbasicoff(&selbasicoff); + + MHBlindPixels blindOFF("BlindPixelsOFF"); + MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels"); + fillblindOFF.SetName("FillBlindOFF"); + + MSigmabarCalc sigbarcalcoff; + + MHSigmaTheta sigthOFF("SigmaThetaOFF"); + MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); + fillsigthetaOFF.SetName("FillSigThetaOFF"); + + //***************************** + // entries in MParList + + plistoff.AddToList(&tlistoff); + InitBinnings(&plistoff); + plistoff.AddToList(&blindOFF); + plistoff.AddToList(&sigthOFF); + + + //***************************** + // entries in MTaskList + + tlistoff.AddToList(&readOFF); + //tlistoff.AddToList(&contthetaoff); + + tlistoff.AddToList(&blindoff); + + tlistoff.AddToList(&contbasicoff); + tlistoff.AddToList(&fillblindOFF); + tlistoff.AddToList(&sigbarcalcoff); + tlistoff.AddToList(&fillsigthetaOFF); + + MProgressBar baroff; + MEvtLoop evtloopoff; + evtloopoff.SetParList(&plistoff); + evtloopoff.SetProgressBar(&baroff); + + Int_t maxeventsoff = -1; + //Int_t maxeventsoff = 20000; + if ( !evtloopoff.Eventloop(maxeventsoff) ) + return; + + tlistoff.PrintStatistics(0, kTRUE); + + blindOFF.DrawClone(); + sigthOFF.DrawClone(); + + // save the histograms for the padding + TH2D *sigthoff = sigthOFF.GetSigmaTheta(); + TH3D *sigpixthoff = sigthOFF.GetSigmaPixTheta(); + TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta(); + + TH2D *blindidthoff = blindOFF.GetBlindId(); + TH2D *blindnthoff = blindOFF.GetBlindN(); + + + //----------------------------------------- + + gLog << "End of generating the padding histograms" << endl; + gLog << "=====================================================" << endl; + + pad.MergeHistograms(sigthon, sigthoff, + sigpixthon, sigpixthoff, + diffpixthon, diffpixthoff, + blindidthon, blindidthoff, + blindnthon, blindnthoff); + + if (WPad) + { + // write the padding histograms onto a file --------- + pad.WriteTargetDist(outNameSigTh); + } + } + + // read the padding histograms --------------------------- + if (RPad) + { + pad.ReadTargetDist(outNameSigTh); + } + + + //************************************************************ + + if (Wout) + { + gLog << "=====================================================" << endl; + gLog << "Start the padding" << endl; + + //----------------------------------------------------------- + MTaskList tliston; + MParList pliston; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenamein); + read.DisableAutoScheme(); + + MGeomApply apply; + MMcPedestalCopy pcopy; + MMcPedestalNSBAdd pnsb; + + MPedestalWorkaround waround; + + // a way to find out whether one is dealing with MC : + MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC + f1.SetName("Select MC"); + MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data + f2.SetName("Select Data"); + + MCerPhotCalc ncalc; + ncalc.SetFilter(&f1); + MCerPhotAnal2 nanal; + nanal.SetFilter(&f2); + + //if (typeInput == "ON") + //{ + // MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", + // "MCT1PointingCorrCalc"); + //} + + MBlindPixelCalc blindbeforepad; + //blindbeforepad.SetUseBlindPixels(); + blindbeforepad.SetName("BlindBeforePadding"); + + MBlindPixelCalc blind; + //blind.SetUseBlindPixels(); + blind.SetName("BlindAfterPadding"); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels"); + fillblind.SetName("HBlind"); + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetName("HSigmaTheta"); + + MImgCleanStd clean; + + + // calculation of image parameters --------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + MHillasCalc hcalc; + hcalc.SetNameHillas(fHilName); + hcalc.SetNameHillasExt(fHilNameExt); + hcalc.SetNameNewImgPar(fImgParName); + + MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc); + hsrccalc.SetInput(fHilName); + + 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"); + // -------------------------------------------------- + + MFSelStandard selstandard(fHilNameSrc); + selstandard.SetHillasName(fHilName); + selstandard.SetImgParName(fImgParName); + selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); + MContinue contstandard(&selstandard); + contstandard.SetName("SelStandard"); + + + MWriteRootFile write(outNameImage); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + //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"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&source); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&f1); + tliston.AddToList(&f2); + tliston.AddToList(&apply); + tliston.AddToList(&pcopy); + //tliston.AddToList(&waround); + + tliston.AddToList(&pnsb); + tliston.AddToList(&ncalc); + tliston.AddToList(&nanal); + + tliston.AddToList(&blindbeforepad); + // tliston.AddToList(&pad); + // if (typeInput == "ON") + // tliston.AddToList(&pointcorr); + tliston.AddToList(&blind); + tliston.AddToList(&contbasic); + + tliston.AddToList(&fillblind); + tliston.AddToList(&sigbarcalc); + tliston.AddToList(&fillsigtheta); + tliston.AddToList(&clean); + + tliston.AddToList(&hcalc); + tliston.AddToList(&hsrccalc); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contstandard); + tliston.AddToList(&write); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + //evtloop.ReadEnv(env, "", printEnv); + evtloop.SetProgressBar(&bar); + // evtloop.Write(); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + + pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone(); + pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "End of padding" << endl; + gLog << "=====================================================" << endl; + } + + + gLog << "Macro MagicAnalysis : End of Job A" << endl; + gLog << "===================================================" << endl; + } + + + + + + //--------------------------------------------------------------------- + // Job B_RF_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 MagicAnalysis : Start of Job B_RF_UP" << endl; + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : 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"; + + + TString extin = "1.root"; + TString extout = "2.root"; + + //-------------------------------------------- + // for the analysis using ON data only set typeMatrixHadrons = "ON" + // ON and OFF data = "OFF" + TString typeMatrixHadrons = "OFF"; + gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + + // 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 + TString NameData = outPath; + NameData += typeInput; + TString inNameData(NameData); + inNameData += extin; + gLog << "inNameData = " << inNameData << endl; + + // name of output root file + TString outNameData(NameData); + outNameData += extout; + gLog << "outNameData = " << outNameData << endl; + + //-------------------------------------------- + // files to be read for generating + // - the matrices of training events + // - and the root files of training and test events + + + // "hadrons" : + TString filenameHad = outPath; + filenameHad += typeMatrixHadrons; + filenameHad += extin; + Int_t howManyHadronsTrain = 12000; + Int_t howManyHadronsTest = 12000; + gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = " + << howManyHadronsTrain << ", howManyHadronsTest = " + << howManyHadronsTest << endl; + + + // "gammas" : + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += extin; + Int_t howManyGammasTrain = 12000; + Int_t howManyGammasTest = 12000; + gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = " + << howManyGammasTrain << ", howManyGammasTest = " + << howManyGammasTest << endl; + + //-------------------------------------------- + // files for the matrices of training events + + TString NameGammas = outPath; + NameGammas += "RFmatrix_gammas_Train_"; + NameGammas += "MC"; + NameGammas += extin; + + TString NameHadrons = outPath; + NameHadrons += "RFmatrix_hadrons_Train_"; + NameHadrons += typeMatrixHadrons; + NameHadrons += extin; + + + //-------------------------------------------- + // root files for the training events + + TString NameGammasTrain = outPath; + NameGammasTrain += "RF_gammas_Train_"; + NameGammasTrain += "MC"; + TString inNameGammasTrain(NameGammasTrain); + inNameGammasTrain += extin; + TString outNameGammasTrain(NameGammasTrain); + outNameGammasTrain += extout; + + + TString NameHadronsTrain = outPath; + NameHadronsTrain += "RF_hadrons_Train_"; + NameHadronsTrain += typeMatrixHadrons; + TString inNameHadronsTrain(NameHadronsTrain); + inNameHadronsTrain += extin; + TString outNameHadronsTrain(NameHadronsTrain); + outNameHadronsTrain += extout; + + + //-------------------------------------------- + // root files for the test events + + TString NameGammasTest = outPath; + NameGammasTest += "RF_gammas_Test_"; + NameGammasTest += "MC"; + TString inNameGammasTest(NameGammasTest); + inNameGammasTest += extin; + TString outNameGammasTest(NameGammasTest); + outNameGammasTest += extout; + + TString NameHadronsTest = outPath; + NameHadronsTest += "RF_hadrons_Test_"; + NameHadronsTest += typeMatrixHadrons; + TString inNameHadronsTest(NameHadronsTest); + inNameHadronsTest += extin; + TString outNameHadronsTest(NameHadronsTest); + outNameHadronsTest += extout; + + //-------------------------------------------------------------------- + + + 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 training events +if (RTrainRF) + { + const char* mtxName = "MatrixGammas"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Get matrix for (gammas)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << NameGammas << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'NameGammas' + // + TFile fileg(NameGammas); + + matrixg.Read(mtxName); + matrixg.Print("SizeCols"); + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Get matrix for (hadrons)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << NameHadrons << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'NameHadrons' + // + TFile fileh(NameHadrons); + + matrixh.Read(mtxName); + matrixh.Print("SizeCols"); + } + + + //************************************************************************* + // create matrices of training events + // and root files of training and test events + +if (CTrainRF) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Create matrices of training events and root files of training and test events" + << endl; + gLog << " Gammas :" << endl; + gLog << "---------" << endl; + + MParList plistg; + MTaskList tlistg; + + MReadMarsFile readg("Events", filenameMC); + readg.DisableAutoScheme(); + + 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(howManyGammasTrain+howManyGammasTest); + selectorg.SetName("selectGammasTrainTest"); + selectorg.SetInverted(); + selectorg.SetUseOrigDistribution(kTRUE); + + MContinue contg(&selectorg); + contg.SetName("ContGammas"); + + Double_t probg = ( (Double_t) howManyGammasTrain ) + / ( (Double_t)(howManyGammasTrain+howManyGammasTest) ); + MFRandomSplit splitg(probg); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&splitg); + fillmatg.SetName("fillGammas"); + + //----------------------- + // for writing the root files of training and test events + // for gammas + + MWriteRootFile writetraing(inNameGammasTrain, "RECREATE"); + writetraing.SetName("WriteGammasTrain"); + writetraing.SetFilter(&splitg); + + writetraing.AddContainer("MRawRunHeader", "RunHeaders"); + writetraing.AddContainer("MTime", "Events"); + writetraing.AddContainer("MMcEvt", "Events"); + writetraing.AddContainer("ThetaOrig", "Events"); + writetraing.AddContainer("MSrcPosCam", "Events"); + writetraing.AddContainer("MSigmabar", "Events"); + writetraing.AddContainer("MHillas", "Events"); + writetraing.AddContainer("MHillasExt", "Events"); + writetraing.AddContainer("MHillasSrc", "Events"); + writetraing.AddContainer("MNewImagePar", "Events"); + + MContinue contgtrain(&splitg); + contgtrain.SetName("ContGammaTrain"); + + MWriteRootFile writetestg(inNameGammasTest, "RECREATE"); + writetestg.SetName("WriteGammasTest"); + + writetestg.AddContainer("MRawRunHeader", "RunHeaders"); + writetestg.AddContainer("MTime", "Events"); + writetestg.AddContainer("MMcEvt", "Events"); + writetestg.AddContainer("ThetaOrig", "Events"); + writetestg.AddContainer("MSrcPosCam", "Events"); + writetestg.AddContainer("MSigmabar", "Events"); + writetestg.AddContainer("MHillas", "Events"); + writetestg.AddContainer("MHillasExt", "Events"); + writetestg.AddContainer("MHillasSrc", "Events"); + writetestg.AddContainer("MNewImagePar", "Events"); + + //----------------------- + + //***************************** fill gammas *** + // entries in MParList + + plistg.AddToList(&tlistg); + InitBinnings(&plistg); + + plistg.AddToList(&matrixg); + + //***************************** + // entries in MTaskList + + tlistg.AddToList(&readg); + tlistg.AddToList(&contg); + + tlistg.AddToList(&splitg); + tlistg.AddToList(&fillmatg); + tlistg.AddToList(&writetraing); + tlistg.AddToList(&contgtrain); + + tlistg.AddToList(&writetestg); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtloopg; + evtloopg.SetName("FillGammaMatrix"); + evtloopg.SetParList(&plistg); + //evtloopg.ReadEnv(env, "", printEnv); + evtloopg.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtloopg.Eventloop(maxevents)) + return; + + tlistg.PrintStatistics(0, kTRUE); + + matrixg.Print("SizeCols"); + Int_t generatedgTrain = matrixg.GetM().GetNrows(); + if ( fabs(generatedgTrain-howManyGammasTrain) > + 3.0*sqrt(howManyGammasTrain) ) + { + gLog << "ONOFFAnalysis.C : no.of generated gamma training events (" + << generatedgTrain << ") is incompatible with the no.of requested events (" + << howManyGammasTrain << ")" << endl; + } + + + Int_t generatedgTest = writetestg.GetNumExecutions(); + if ( fabs(generatedgTest-howManyGammasTest) > + 3.0*sqrt(howManyGammasTest) ) + { + gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" + << generatedgTest << ") is incompatible with the no.of requested events (" + << howManyGammasTest << ")" << endl; + } + + //***************************** fill hadrons *** + gLog << "---------------------------------------------------------------" + << endl; + gLog << " Hadrons :" << endl; + gLog << "----------" << endl; + + MParList plisth; + MTaskList tlisth; + + MReadMarsFile readh("Events", filenameHad); + readh.DisableAutoScheme(); + + 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); + + //use the original distribution from the gammas + MH3 &href = *(selectorg.GetHistOrig()); + + MFEventSelector2 selectorh(href); + selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest); + selectorh.SetName("selectHadronsTrainTest"); + selectorh.SetInverted(); + + MContinue conth(&selectorh); + conth.SetName("ContHadrons"); + + Double_t probh = ( (Double_t) howManyHadronsTrain ) + / ( (Double_t)(howManyHadronsTrain+howManyHadronsTest) ); + MFRandomSplit splith(probh); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&splith); + fillmath.SetName("fillHadrons"); + + //----------------------- + // for writing the root files of training and test events + // for hadrons + + MWriteRootFile writetrainh(inNameHadronsTrain, "RECREATE"); + writetrainh.SetName("WriteHadronsTrain"); + writetrainh.SetFilter(&splith); + + writetrainh.AddContainer("MRawRunHeader", "RunHeaders"); + writetrainh.AddContainer("MTime", "Events"); + writetrainh.AddContainer("MMcEvt", "Events"); + writetrainh.AddContainer("ThetaOrig", "Events"); + writetrainh.AddContainer("MSrcPosCam", "Events"); + writetrainh.AddContainer("MSigmabar", "Events"); + writetrainh.AddContainer("MHillas", "Events"); + writetrainh.AddContainer("MHillasExt", "Events"); + writetrainh.AddContainer("MHillasSrc", "Events"); + writetrainh.AddContainer("MNewImagePar", "Events"); + + MContinue conthtrain(&splith); + + MWriteRootFile writetesth(inNameHadronsTest, "RECREATE"); + writetesth.SetName("WriteHadronsTest"); + + writetesth.AddContainer("MRawRunHeader", "RunHeaders"); + writetesth.AddContainer("MTime", "Events"); + writetesth.AddContainer("MMcEvt", "Events"); + writetesth.AddContainer("ThetaOrig", "Events"); + writetesth.AddContainer("MSrcPosCam", "Events"); + writetesth.AddContainer("MSigmabar", "Events"); + writetesth.AddContainer("MHillas", "Events"); + writetesth.AddContainer("MHillasExt", "Events"); + writetesth.AddContainer("MHillasSrc", "Events"); + writetesth.AddContainer("MNewImagePar", "Events"); + + + //***************************** + // entries in MParList + + plisth.AddToList(&tlisth); + InitBinnings(&plisth); + + plisth.AddToList(&matrixh); + + //***************************** + // entries in MTaskList + + tlisth.AddToList(&readh); + tlisth.AddToList(&conth); + + tlisth.AddToList(&splith); + tlisth.AddToList(&fillmath); + tlisth.AddToList(&writetrainh); + tlisth.AddToList(&conthtrain); + + tlisth.AddToList(&writetesth); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtlooph; + evtlooph.SetName("FillHadronMatrix"); + evtlooph.SetParList(&plisth); + //evtlooph.ReadEnv(env, "", printEnv); + evtlooph.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtlooph.Eventloop(maxevents)) + return; + + tlisth.PrintStatistics(0, kTRUE); + + matrixh.Print("SizeCols"); + Int_t generatedhTrain = matrixh.GetM().GetNrows(); + if ( fabs(generatedhTrain-howManyHadronsTrain) > + 3.0*sqrt(howManyHadronsTrain) ) + { + gLog << "ONOFFAnalysis.C : no.of generated hadron training events (" + << generatedhTrain << ") is incompatible with the no.of requested events (" + << howManyHadronsTrain << ")" << endl; + } + + + Int_t generatedhTest = writetesth.GetNumExecutions(); + if ( fabs(generatedhTest-howManyHadronsTest) > + 3.0*sqrt(howManyHadronsTest) ) + { + gLog << "ONOFFAnalysis.C : no.of generated gamma test events (" + << generatedhTest << ") is incompatible with the no.of requested events (" + << howManyHadronsTest << ")" << endl; + } + + + //***************************************************** + + + // write out matrices of training events + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Write out matrices of training events" << endl; + + + //------------------------------------------- + // "gammas" + gLog << "Gammas :" << endl; + matrixg.Print("SizeCols"); + + TFile writeg(NameGammas, "RECREATE", ""); + matrixg.Write(); + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : matrix of training events for gammas written onto file " + << NameGammas << endl; + + //------------------------------------------- + // "hadrons" + gLog << "Hadrons :" << endl; + matrixh.Print("SizeCols"); + + TFile writeh(NameHadrons, "RECREATE", ""); + matrixh.Write(); + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : matrix of training events for hadrons written onto file " + << NameHadrons << endl; + + } + //********** end of creating matrices of training events *********** + + + MRanForest *fRanForest; + MRanTree *fRanTree; + //----------------------------------------------------------------- + // 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); + + gLog << "ONOFFAnalysis : RF trees were read in from file " + << outRF << endl; + + // 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 MagicAnalysis : 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"); + + MFillH fillh2("MHRanForestGini"); + + // list of tasks for the loop over the trees + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); + + //------------------- + // Execute tree loop + // + MEvtLoop treeloop; + treeloop.SetName("GrowRFTrees"); + treeloop.SetParList(&plist2); + + if ( !treeloop.Eventloop() ) + return; + + tlist2.PrintStatistics(0, kTRUE); + + plist2.FindObject("MHRanForestGini")->DrawClone(); + + + // 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 root files with the RF hadronness + // + + if (WRF) + { + //TString fileName(inNameHadronsTrain); + //TString outName(outNameHadronsTrain); + + //TString fileName(inNameHadronsTest); + //TString outName(outNameHadronsTest); + + //TString fileName(inNameGammasTrain); + //TString outName(outNameGammasTrain); + + //TString fileName(inNameGammasTest); + //TString outName(outNameGammasTest); + + TString fileName(inNameData); + TString outName(outNameData); + + + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update root file '" << fileName + << "' with the RF hadronness; ==> " << outName << endl; + + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", fileName); + read.DisableAutoScheme(); + + + //....................................................................... + // calculate hadronnes for method of RANDOM FOREST + + + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadRFName); + + + //....................................................................... + + //MWriteRootFile write(outName, "UPDATE"); + MWriteRootFile write(outName, "RECREATE"); + + 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(hadRFName, "Events"); + + //----------------------------------------------------------------- + + + MFSelFinal 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"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + 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"); + + + 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); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillranfor); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + 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.SetName("UpdateRootFile"); + 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("hadRF", "MHHadronness")->Print(); + + 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 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 '" << fileName << "' : " << 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 MagicAnalysis : End of Job B_RF_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job B_SC_UP + //============ + + // - create (or read in) optimum supercuts parameter values + // + // - calculate the hadroness for the supercuts + // + // - update input root file, including the hadroness + + + if (JobB_SC_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro MagicAnalysis : Start of Job B_SC_UP" << endl; + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = " + << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " + << (RMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "kTRUE" : "kFALSE") << ", " + << (WSC ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // file which contains the initial parameter values for the supercuts + // if parSCinit ="" the initial values are taken from the constructor of + // MSupercuts + + TString parSCinit = outPath; + //parSCinit += "parSC_1709d"; + parSCinit = ""; + + gLog << "parSCinit = " << parSCinit << endl; + + //--------------- + // file onto which the optimal parameter values for the supercuts + // are written + + TString parSCfile = outPath; + parSCfile += "parSC_2310a"; + + gLog << "parSCfile = " << parSCfile << endl; + + //-------------------------------------------- + // file to be updated (either ON or MC) + + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + // name of input root file + TString filenameData = outPath; + filenameData += typeInput; + filenameData += "2.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "3.root"; + + + //TString outNameImage = filenameData; + + gLog << "outNameImage = " << outNameImage << endl; + + //-------------------------------------------- + // files to be read for optimizing the supercuts + // + // for the training + TString filenameTrain = outPath; + filenameTrain += "ON"; + filenameTrain += "1.root"; + Int_t howManyTrain = 800000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = outPath; + filenameTest += "ON"; + filenameTest += "1.root"; + Int_t howManyTest = 800000; + + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // files to contain the matrices (generated from filenameTrain and + // filenameTest) + // + // for the training + TString fileMatrixTrain = outPath; + fileMatrixTrain += "MatrixTrainSC"; + fileMatrixTrain += ".root"; + gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; + + // for testing + TString fileMatrixTest = outPath; + fileMatrixTest += "MatrixTestSC"; + fileMatrixTest += ".root"; + gLog << "fileMatrixTest = " << fileMatrixTest << endl; + + + + //--------------------------------------------------------------------- + // Training and test matrices : + // - either create them and write them onto a file + // - or read them from a file + + + MFindSupercuts findsuper; + findsuper.SetFilenameParam(parSCfile); + findsuper.SetHadronnessName("HadSC"); + findsuper.SetUseOrigDistribution(kTRUE); + + //-------------------------- + // create matrices and write them onto files + if (CMatrix) + { + TString mname("costheta"); + MBinning bin("Binning"+mname); + bin.SetEdges(10, 0., 1.0); + + MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mname); + MH::SetBinning(&mh3.GetHist(), &bin); + for (Int_t i=1; i<=mh3.GetNbins(); i++) + mh3.GetHist().SetBinContent(i, 1.0); + + + if (filenameTrain == filenameTest) + { + if ( !findsuper.DefineTrainTestMatrix( + filenameTrain, mh3, + howManyTrain, howManyTest, + fileMatrixTrain, fileMatrixTest) ) + { + *fLog << "MagicAnalysis.C : DefineTrainTestMatrix failed" << endl; + return; + } + + } + else + { + if ( !findsuper.DefineTrainMatrix(filenameTrain, mh3, + howManyTrain, fileMatrixTrain) ) + { + *fLog << "MagicAnalysis.C : DefineTrainMatrix failed" << endl; + return; + } + + if ( !findsuper.DefineTestMatrix( filenameTest, mh3, + howManyTest, fileMatrixTest) ) + { + *fLog << "MagicAnalysis.C : DefineTestMatrix failed" << endl; + return; + } + } + } + + //-------------------------- + // read matrices from files + // + + if (RMatrix) + findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); + //-------------------------- + + + + //--------------------------------------------------------------------- + // optimize supercuts using the training sample + // + // the initial values are taken + // - from the file parSCinit (if != "") + // - or from the arrays params and steps (if their sizes are != 0) + // - or from the MSupercuts constructor + + +if (WOptimize) + { + gLog << "MagicAnalysis.C : optimize the supercuts using the training matrix" + << endl; + + TArrayD params(0); + TArrayD steps(0); + + if (parSCinit == "") + { + Double_t vparams[104] = { + // LengthUp + 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362, + 0.007388, -0.013463, + // LengthLo + 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993, + 0.000080, -0.007157, + // WidthUp + 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353, + 0.020711, -0.016703, + // WidthLo + 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101, + 0.006126, -0.002849, + // DistUp + 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909, + 0.189697, 0.0, + // DistLo + 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374, + -0.001750, 0.0, + // AsymUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + Double_t vsteps[104] = { + // LengthUp + 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002, + 0.0008, 0.002, + // LengthLo + 0.02, 0.003, 0.05, 0.006, 0.002, 0.3, + 0.0001, 0.0008, + // WidthUp + 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006, + 0.002, 0.002, + // WidthLo + 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002, + 0.0007, 0.003, + // DistUp + 0.2, 0.0, 0.3, 0.02, 0.0, 0.03, + 0.02, 0.0 + // DistLo + 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005, + 0.0002, 0.0 + // AsymUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + params.Set(104, vparams); + steps.Set (104, vsteps ); + } + + Bool_t rf; + rf = findsuper.FindParams(parSCinit, params, steps); + + if (!rf) + { + gLog << "MagicAnalysis.C : optimization of supercuts failed" << endl; + return; + } + } + + //-------------------------------------- + // test the supercuts on the test sample + // + + if (RTest) + { + gLog << "MagicAnalysis.C : test the supercuts on the test matrix" << endl; + Bool_t rt = findsuper.TestParams(); + if (!rt) + { + gLog << "MagicAnalysis.C : test of supercuts on the test matrix failed" + << endl; + } + + } + + + //----------------------------------------------------------------- + // Update the input files with the SC hadronness + // + + if (WSC) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update input file '" << filenameData + << "' with the SC hadronness" << endl; + + + //---------------------------------------------------- + // read in optimum parameter values for the supercuts + + TFile inparam(parSCfile); + MSupercuts scin; + scin.Read("MSupercuts"); + inparam.Close(); + + gLog << "Parameter values for supercuts were read in from file '" + << parSCfile << "'" << endl; + + TArrayD supercutsPar; + supercutsPar = scin.GetParameters(); + + TArrayD supercutsStep; + supercutsStep = scin.GetStepsizes(); + + gLog << "Parameter values for supercuts : " << endl; + for (Int_t i=0; iFindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + + //....................................................................... + // calculation of hadroness for the supercuts + // (=0.25 if fullfilled, =0.75 otherwise) + + TString hadSCName = "HadSC"; + MSupercutsCalc sccalc(fHilName, fHilNameSrc); + sccalc.SetHadronnessName(hadSCName); + + + //....................................................................... + + + //MWriteRootFile write(outNameImage, "UPDATE"); + //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE"); + + + MWriteRootFile write(outNameImage, "RECREATE"); + + 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("HadRF", "Events"); + write.AddContainer(hadSCName, "Events"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadSCName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadSCName); + 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); + + TH1 &alphahist = alphaabs->GetHist(); + + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); + + 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(&supercuts); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadsc); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + 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; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + 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 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 MagicAnalysis : End of Job B_SC_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 Random Forest and for the SUPERCUTS + // - produce Neyman-Pearson plots + + if (JobC) + { + gLog << "=====================================================" << endl; + gLog << "Macro MagicAnalysis : Start of Job C" << endl; + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : JobC = " + << (JobC ? "kTRUE" : "kFALSE") << endl; + + + + TString ext("2.root"); + TString extout("2.root"); + + TString typeHadrons("OFF"); + TString typeGammas("MC"); + + //-------------------------------------------- + // name of input data file + TString NameData = outPath; + NameData += typeHadrons; + TString inNameData(NameData); + inNameData += ext; + gLog << "inNameData = " << inNameData << endl; + + // name of input MC file + TString NameMC = outPath; + NameMC += typeGammas; + TString inNameMC(NameMC); + inNameMC += ext; + gLog << "inNameMC = " << inNameMC << endl; + + + //-------------------------------------------- + // root files for the training events + + + + TString NameGammasTrain = outPath; + NameGammasTrain += "RF_gammas_Train_"; + NameGammasTrain += typeGammas; + TString outNameGammasTrain(NameGammasTrain); + outNameGammasTrain += extout; + + + TString NameHadronsTrain = outPath; + NameHadronsTrain += "RF_hadrons_Train_"; + NameHadronsTrain += typeHadrons; + TString outNameHadronsTrain(NameHadronsTrain); + outNameHadronsTrain += extout; + + + //-------------------------------------------- + // root files for the test events + + TString NameGammasTest = outPath; + NameGammasTest += "RF_gammas_Test_"; + NameGammasTest += typeGammas; + TString outNameGammasTest(NameGammasTest); + outNameGammasTest += extout; + + TString NameHadronsTest = outPath; + NameHadronsTest += "RF_hadrons_Test_"; + NameHadronsTest += typeHadrons; + TString outNameHadronsTest(NameHadronsTest); + outNameHadronsTest += extout; + + //-------------------------------------------------------------------- + + //TString filenameData(inNameData); + //TString filenameMC(inNameMC); + + //TString filenameData(outNameHadronsTrain); + //TString filenameMC(outNameGammasTrain); + + TString filenameData(outNameHadronsTest); + TString filenameMC(outNameGammasTest); + + gLog << "filenameData = " << filenameData << endl; + gLog << "filenameMC = " << filenameMC << endl; + + //----------------------------------------------------------------- + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + + //....................................................................... + // names of hadronness containers + + 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; + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadRFName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + + //MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + //fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFSelFinal 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(&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("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 MagicAnalysis : 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 MagicAnalysis : Start of Job D" << endl; + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : JobD = " + << (JobD ? "kTRUE" : "kFALSE") << endl; + + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString ext = "3.root"; + + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //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.233; + 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("MGeomCamMagic", "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("MGeomCamMagic", "MGeomCam"); + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + + TH1 &alphahist = alphaabs->GetHist(); + + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); + + MFillH hfill1("MHHillas", fHilName); + 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"); + + MFSelFinal 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); + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&alpha); + 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("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 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-JobD"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << filenameData << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + + DeleteBinnings(&pliston); + + gLog << "Macro MagicAnalysis : End of Job D" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + + + //--------------------------------------------------------------------- + // Job E_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 (JobE_XX) + { + gLog << "=====================================================" << endl; + gLog << "Macro MagicAnalysis : Start of Job E_XX" << endl; + + gLog << "" << endl; + gLog << "Macro MagicAnalysis : JobE_XX, CCollArea, OEEst, WEX = " + << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (CCollArea?"kTRUE" : "kFALSE") << ", " + << (OEEst ? "kTRUE" : "kFALSE") << ", " + << (WEX ? "kTRUE" : "kFALSE") << 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("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.23; + 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 optimizing the energy estimator + TString filenameOpt(outPath); + filenameOpt += typeMC; + filenameOpt += ext; + gLog << "filenameOpt = " << filenameOpt << 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 MC file to be used for calculating the eff. collection areas + TString filenameArea(outPath); + filenameArea += typeMC; + 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 += 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; + + //------------------------------ + // name of file containing histograms for flux calculastion + TString filenameResults(outPath); + filenameResults += typeData; + filenameResults += "Results_"; + filenameResults += XX; + filenameResults += extout; + gLog << "filenameResults = " << filenameResults << endl; + + + //==================================================================== + + MHMcCollectionArea collarea; + collarea.SetEaxis(MHMcCollectionArea::kLinear); + + MParList parlist; + InitBinnings(&parlist); + + if (CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Start calculation of effective collection areas" << endl; + + + MTaskList tasklist; + + //--------------------------------------- + // Setup the tasks to be executed + // + MReadMarsFile reader("Events", filenameArea); + reader.DisableAutoScheme(); + + MFSelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); + + MContinue conthadrons(&cuthadrons); + + + MFillH filler("MHMcCollectionArea", "MMcEvt"); + filler.SetName("CollectionArea"); + + //******************************** + // entries in MParList + + parlist.AddToList(&tasklist); + + 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 + // + //MHMcCollectionArea *collarea = + // (MHMcCollectionArea*)parlist.FindObject("MHMcCollectionArea"); + collarea.CalcEfficiency(); + collarea.DrawClone(); + + + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + //collarea.GetHist()->Write(); + //collarea.GetHAll()->Write(); + //collarea.GetHSel()->Write(); + collarea.Write(); + + f.Close(); + + gLog << "Collection area plots written onto file " << collareaName << endl; + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + } + + if (!CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Read in effective collection areas from file " + << collareaName << endl; + + TFile collfile(collareaName); + collfile.ls(); + collarea.Read("MHMcCollectionArea"); + collarea.DrawClone(); + + gLog << "Effective collection areas were read in from file " + << collareaName << endl; + gLog << "-----------------------------------------------" << endl; + } + + + // save binnings for call to MagicEEst + MBinning *binsE = (MBinning*)parlist.FindObject("BinningE"); + if (!binsE) + { + gLog << "Object 'BinningE' not found in MParList" << endl; + return; + } + MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta"); + if (!binsTheta) + { + gLog << "Object 'BinningTheta' not found in MParList" << endl; + return; + } + + //------------------------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + + if (OEEst) + { + //=========================================================== + // + // Optimization of energy estimator + // + gLog << "Macro MagicAnalysis.C : calling MagicEEst" << endl; + + TString inpath(""); + TString outpath(""); + Int_t howMany = 2000; + MagicEEst(inpath, filenameOpt, outpath, energyParName, + fHilName, fHilNameSrc, fhadronnessName, + howMany, maxhadronness, maxalpha, maxdist, + binsE, binsTheta); + gLog << "Macro MagicAnalysis.C : returning from MagicEEst" << endl; + } + + if (WEX) + { + //----------------------------------------------------------- + // + // Read in parameters of energy estimator ("MMcEnergyEst") + // and migration matrix ("MHMcEnergyMigration") + // + gLog << "================================================================" + << endl; + gLog << "Macro MagicAnalysis.C : read parameters of energy estimator and migration matrix from file '" + << energyParName << "'" << endl; + TFile enparam(energyParName); + enparam.ls(); + MMcEnergyEst mcest("MMcEnergyEst"); + mcest.Read("MMcEnergyEst"); + + //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst")); + gLog << "Parameters of energy estimator were read in" << endl; + + + gLog << "Read in Migration matrix" << endl; + + MHMcEnergyMigration mighiston("MHMcEnergyMigration"); + mighiston.Read("MHMcEnergyMigration"); + //MHMcEnergyMigration &mighiston = + // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration")); + + gLog << "Migration matrix was read in" << endl; + + + TArrayD parA(mcest.GetNumCoeffA()); + TArrayD parB(mcest.GetNumCoeffB()); + for (Int_t i=0; iFindCreateObj("MGeomCamMagic", "MGeomCam"); + + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + //....................................................................... + + gLog << "MagicAnalysis.C : write root file '" << filenameDataout + << "'" << endl; + + //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + + + MWriteRootFile write(filenameDataout, "RECREATE"); + + 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"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + MFSelFinal 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"); + + //--------------------------- + // calculate estimated energy + + MEnergyEstParam eeston(fHilName); + eeston.Add(fHilNameSrc); + + eeston.SetCoeffA(parA); + eeston.SetCoeffB(parB); + + //--------------------------- + // calculate estimated energy using Daniel's parameters + + //MEnergyEstParamDanielMkn421 eeston(fHilName); + //eeston.Add(fHilNameSrc); + //eeston.SetCoeffA(parA); + //eeston.SetCoeffB(parB); + + + //--------------------------- + + + 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"); + + //--------------------------- + // new from Robert + + MFillH hfill6("MHTimeDiffTheta", "MMcEvt"); + hfill6.SetName("HTimeDiffTheta"); + + MFillH hfill6a("MHTimeDiffTime", "MMcEvt"); + hfill6a.SetName("HTimeDiffTime"); + + MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc); + hfill7.SetName("HAlphaEnergyTheta"); + + MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc); + hfill7a.SetName("HAlphaEnergyTime"); + + MFillH hfill7b("MHThetabarTime", fHilNameSrc); + hfill7b.SetName("HThetabarTime"); + + MFillH hfill7c("MHEnergyTime", "MMcEvt"); + hfill7c.SetName("HEnergyTime"); + + + //--------------------------- + + MFSelFinal 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); + + // robert + tliston.AddToList(&hfill6); //timediff + tliston.AddToList(&hfill6a); //timediff + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&eeston); + + tliston.AddToList(&write); + + //tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + //robert + tliston.AddToList(&hfill7); + tliston.AddToList(&hfill7a); + tliston.AddToList(&hfill7b); + tliston.AddToList(&hfill7c); + + 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("hadNN", "MHHadronness")->DrawClone(); + + gLog << "before hadRF" << endl; + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + gLog << "before hadSC" << endl; + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + gLog << "before MHHillas" << endl; + pliston.FindObject("MHHillas")->DrawClone(); + + gLog << "before MHHillasExt" << endl; + pliston.FindObject("MHHillasExt")->DrawClone(); + + gLog << "before MHHillasSrc" << endl; + pliston.FindObject("MHHillasSrc")->DrawClone(); + + gLog << "before MHNewImagePar" << endl; + pliston.FindObject("MHNewImagePar")->DrawClone(); + + gLog << "before MHStarMap" << endl; + pliston.FindObject("MHStarMap")->DrawClone(); + + gLog << "before DeleteBinnings" << endl; + + DeleteBinnings(&pliston); + + gLog << "before Robert's code" << endl; + + +//rwagner write all relevant histograms onto a file + + if (WRobert) + { + gLog << "=======================================================" << endl; + gLog << "Write results onto file '" << filenameResults << "'" << endl; + + TFile outfile(filenameResults,"recreate"); + + MHHillasSrc* hillasSrc = + (MHHillasSrc*)(pliston->FindObject("MHHillasSrc")); + TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + alphaHist->Write(); + gLog << "Alpha plot has been written out" << endl; + + + MHAlphaEnergyTheta* aetH = + (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); + TH3D* aetHist = (TH3D*)(aetH->GetHist()); + aetHist->SetName("aetHist"); + aetHist->Write(); + gLog << "AlphaEnergyTheta plot has been written out" << endl; + + MHAlphaEnergyTime* aetH2 = + (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime")); + TH3D* aetHist2 = (TH3D*)(aetH2->GetHist()); + aetHist2->SetName("aetimeHist"); +// aetHist2->DrawClone(); + aetHist2->Write(); + gLog << "AlphaEnergyTime plot has been written out" << endl; + + MHThetabarTime* tbt = + (MHThetabarTime*)(pliston->FindObject("MHThetabarTime")); + TProfile* tbtHist = (TProfile*)(tbt->GetHist()); + tbtHist->SetName("tbtHist"); + tbtHist->Write(); + gLog << "ThetabarTime plot has been written out" << endl; + + MHEnergyTime* ent = + (MHEnergyTime*)(pliston->FindObject("MHEnergyTime")); + TH2D* entHist = (TH2D*)(ent->GetHist()); + entHist->SetName("entHist"); + entHist->Write(); + gLog << "EnergyTime plot has been written out" << endl; + + MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); + TH2D* timeHist = (TH2D*)(time->GetHist()); + timeHist->SetName("MHTimeDiffTheta"); + timeHist->SetTitle("Time diffs"); + timeHist->Write(); + gLog << "TimeDiffTheta plot has been written out" << endl; + + + MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); + TH2D* timeHist2 = (TH2D*)(time2->GetHist()); + timeHist2->SetName("MHTimeDiffTime"); + timeHist2->SetTitle("Time diffs"); + timeHist2->Write(); + gLog << "TimeDiffTime plot has been written out" << endl; + +//rwagner write also collareas to same file + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + gLog << "Effective collection areas have been written out" << endl; + +//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data +//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file + + gLog << "before closing outfile" << endl; + + //outfile.Close(); + gLog << "Results were written onto file '" << filenameResults + << "'" << endl; + gLog << "=======================================================" << endl; + } + + } + + gLog << "Macro Analysis : End of Job E_XX" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2662) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2663) @@ -48,4 +48,5 @@ #pragma link C++ class MCT1PadSchweizer+; #pragma link C++ class MCT1PadONOFF+; +#pragma link C++ class MPadONOFF+; #pragma link C++ class MCT1PointingCorrCalc+; @@ -60,4 +61,7 @@ #pragma link C++ class MCT1SupercutsCalc+; #pragma link C++ class MCT1FindSupercuts+; +#pragma link C++ class MSupercuts+; +#pragma link C++ class MSupercutsCalc+; +#pragma link C++ class MFindSupercuts+; #pragma link C++ class MMinuitInterface+; #pragma link C++ class MFiltercutsCalc+; @@ -68,4 +72,5 @@ #pragma link C++ class MCalibrationCam+; #pragma link C++ class MCalibrationCalc+; +#pragma link C++ class MPedestalWorkaround+; #pragma link C++ class MExtractedSignalCam+; Index: trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc (revision 2663) @@ -0,0 +1,1173 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 6/2003 +! Author(s): Wolfgang Wittek, 7/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MFindSupercuts // +// // +// Class for otimizing the parameters of the supercuts // +// // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MFindSupercuts.h" + +#include // fabs + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MContinue.h" +#include "MSupercuts.h" +#include "MSupercutsCalc.h" +#include "MDataElement.h" +#include "MDataMember.h" + +#include "MEvtLoop.h" +#include "MFSelFinal.h" +#include "MF.h" +#include "MFEventSelector.h" +#include "MFEventSelector2.h" +#include "MFillH.h" +//#include "MGeomCamCT1Daniel.h" +#include "MFRandomSplit.h" +#include "MGeomCamMagic.h" +#include "MH3.h" +#include "MHSupercuts.h" +#include "MHFindSignificance.h" +#include "MHMatrix.h" +#include "MHOnSubtraction.h" + +#include "MLog.h" +#include "MLogManip.h" +#include "MMatrixLoop.h" +#include "MMinuitInterface.h" +#include "MParList.h" +#include "MProgressBar.h" +#include "MReadMarsFile.h" +#include "MReadTree.h" +#include "MTaskList.h" + + +ClassImp(MFindSupercuts); + +using namespace std; + + +//------------------------------------------------------------------------ +// +// fcnSupercuts +// +// - calculates the quantity to be minimized (using TMinuit) +// +// - the quantity to be minimized is (-1)*significance of the gamma signal +// in the alpha distribution (after cuts) +// +// - the parameters to be varied in the minimization are the cut parameters +// (par) +// +static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f, + Double_t *par, Int_t iflag) +{ + //cout << "entry fcnSupercuts" << endl; + + //------------------------------------------------------------- + // save pointer to the MINUIT object for optimizing the supercuts + // because it will be overwritten + // when fitting the alpha distribution in MHFindSignificance + TMinuit *savePointer = gMinuit; + //------------------------------------------------------------- + + + MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); + + MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); + //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); + + MSupercuts *super = (MSupercuts*)plistfcn->FindObject("MSupercuts"); + if (!super) + { + gLog << "fcnSupercuts : MSupercuts object '" << "MSupercuts" + << "' not found... aborting" << endl; + return; + } + + // + // transfer current parameter values to MSupercuts + // + // Attention : npar is the number of variable parameters + // not the total number of parameters + // + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + super->SetParameters(TArrayD(fNparx, par)); + + //$$$$$$$$$$$$$$$$$$$$$ + // for testing + //TArrayD checkparameters = super->GetParameters(); + //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " + // << fNparx << endl; + //gLog << "fcnsupercuts : i, par, checkparameters =" << endl; + //for (Int_t i=0; iEventloop(); + + //tasklistfcn->PrintStatistics(0, kTRUE); + + MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); + if (!alpha) + return; + + TH1 &alphaHist = alpha->GetHist(); + alphaHist.SetName("alpha-fcnSupercuts"); + + //------------------------------------------- + // set Minuit pointer to zero in order not to destroy the TMinuit + // object for optimizing the supercuts + gMinuit = NULL; + + //================================================================= + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + const Double_t alphasig = 20.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 2; + + Bool_t drawpoly; + Bool_t fitgauss; + if (iflag == 3) + { + drawpoly = kTRUE; + fitgauss = kTRUE; + } + else + { + drawpoly = kFALSE; + fitgauss = kFALSE; + } + //drawpoly = kFALSE; + //fitgauss = kFALSE; + + const Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + + //================================================================= + + // reset gMinuit to the MINUIT object for optimizing the supercuts + gMinuit = savePointer; + //------------------------------------------- + + if (!rc) + { + gLog << "fcnSupercuts : FindSigma() failed" << endl; + f = 1.e10; + return; + } + + // plot some quantities during the optimization + MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts"); + if (plotsuper) + plotsuper->Fill(&findsig); + + //------------------------ + // get significance + const Double_t significance = findsig.GetSignificance(); + f = significance>0 ? -significance : 0; + + + //------------------------ + // optimize signal/background ratio + //Double_t ratio = findsig.GetNbg()>0.0 ? + // findsig.GetNex()/findsig.GetNbg() : 0.0; + //f = -ratio; + + //------------------------------------------- + // final calculations + //if (iflag == 3) + //{ + // + //} + + //------------------------------------------------------------- +} + + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFindSupercuts::MFindSupercuts(const char *name, const char *title) +: fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL) +{ + fName = name ? name : "MFindSupercuts"; + fTitle = title ? title : "Optimizer of the supercuts"; + + fUseOrigDistribution = kFALSE; + + //--------------------------- + // camera geometry is needed for conversion mm ==> degree + //fCam = new MGeomCamCT1Daniel; + fCam = new MGeomCamMagic; + + // matrices to contain the training/test samples + fMatrixTrain = new MHMatrix("MatrixTrain"); + fMatrixTest = new MHMatrix("MatrixTest"); + + // objects of MSupercutsCalc to which these matrices are attached + fCalcHadTrain = new MSupercutsCalc("SupercutsCalcTrain"); + fCalcHadTest = new MSupercutsCalc("SupercutsCalcTest"); + + // Define columns of matrices + fCalcHadTrain->InitMapping(fMatrixTrain); + fCalcHadTest->InitMapping(fMatrixTest); +} + +// -------------------------------------------------------------------------- +// +// Default destructor. +// +MFindSupercuts::~MFindSupercuts() +{ + delete fCam; + delete fMatrixTrain; + delete fMatrixTest; + delete fCalcHadTrain; + delete fCalcHadTest; +} + + +// -------------------------------------------------------------------------- +// +// +// +void MFindSupercuts::SetUseOrigDistribution(Bool_t b) +{ + fUseOrigDistribution = b; + + *fLog << "MFindSupercuts : when defining training and test matrices use the original distribution" + << endl; +} + +// -------------------------------------------------------------------------- +// +// Define the matrix 'fMatrixTrain' for the training sample +// +// alltogether 'howmanytrain' events are read from file 'nametrain'; +// the events are selected according to a target distribution 'hreftrain' +// +// +Bool_t MFindSupercuts::DefineTrainMatrix( + const TString &nametrain, MH3 &hreftrain, + const Int_t howmanytrain, const TString &filetrain) +{ + if (nametrain.IsNull() || howmanytrain <= 0) + return kFALSE; + + *fLog << "=============================================" << endl; + *fLog << "fill training matrix from file '" << nametrain + << "', select " << howmanytrain + << " events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << hreftrain.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", nametrain); + read.DisableAutoScheme(); + + MFEventSelector2 seltrain(hreftrain); + seltrain.SetNumMax(howmanytrain); + seltrain.SetName("selectTrain"); + if (fUseOrigDistribution) + seltrain.SetUseOrigDistribution(kTRUE); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&seltrain); + filltrain.SetName("fillMatrixTrain"); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&seltrain); + tlist.AddToList(&filltrain); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrain"); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindSupercuts::DefineTrainMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + *fLog << "training matrix was filled" << endl; + *fLog << "=============================================" << endl; + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindSupercuts::DefineTrainMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Define the matrix for the test sample +// +// alltogether 'howmanytest' events are read from file 'nametest' +// +// the events are selected according to a target distribution 'hreftest' +// +// +Bool_t MFindSupercuts::DefineTestMatrix( + const TString &nametest, MH3 &hreftest, + const Int_t howmanytest, const TString &filetest) +{ + if (nametest.IsNull() || howmanytest<=0) + return kFALSE; + + *fLog << "=============================================" << endl; + *fLog << "fill test matrix from file '" << nametest + << "', select " << howmanytest + << " events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << hreftest.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", nametest); + read.DisableAutoScheme(); + + MFEventSelector2 seltest(hreftest); + seltest.SetNumMax(howmanytest); + seltest.SetName("selectTest"); + if (fUseOrigDistribution) + seltest.SetUseOrigDistribution(kTRUE); + + MFillH filltest(fMatrixTest); + filltest.SetFilter(&seltest); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&seltest); + tlist.AddToList(&filltest); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTest"); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTest->Print("SizeCols"); + const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindSupercuts::DefineTestMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + *fLog << "test matrix was filled" << endl; + *fLog << "=============================================" << endl; + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindSupercuts::DefineTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Define the matrices for the training and test sample respectively +// +// +// +Bool_t MFindSupercuts::DefineTrainTestMatrix( + const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest) +{ + *fLog << "=============================================" << endl; + *fLog << "fill training and test matrix from file '" << name + << "', select " << howmanytrain + << " training and " << howmanytest << " test events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << href.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", name); + read.DisableAutoScheme(); + + MFEventSelector2 selector(href); + selector.SetNumMax(howmanytrain+howmanytest); + selector.SetName("selectTrainTest"); + selector.SetInverted(); + if (fUseOrigDistribution) + selector.SetUseOrigDistribution(kTRUE); + + MContinue cont(&selector); + cont.SetName("ContTrainTest"); + + Double_t prob = ( (Double_t) howmanytrain ) + / ( (Double_t)(howmanytrain+howmanytest) ); + MFRandomSplit split(prob); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&split); + filltrain.SetName("fillMatrixTrain"); + + + // consider this event as candidate for a test event + // only if event was not accepted as a training event + + MContinue conttrain(&split); + conttrain.SetName("ContTrain"); + + MFillH filltest(fMatrixTest); + filltest.SetName("fillMatrixTest"); + + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&cont); + + tlist.AddToList(&split); + tlist.AddToList(&filltrain); + tlist.AddToList(&conttrain); + + tlist.AddToList(&filltest); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrainTest"); + evtloop.SetProgressBar(&bar); + + Int_t maxev = -1; + if (!evtloop.Eventloop(maxev)) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtrain + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + fMatrixTest->Print("SizeCols"); + const Int_t generatedtest = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtest + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + + *fLog << "training and test matrix were filled" << endl; + *fLog << "=============================================" << endl; + + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read training and test matrices from files +// +// + +Bool_t MFindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest) +{ + //-------------------------- + // read in training matrix + + TFile filetr(filetrain); + fMatrixTrain->Read("MatrixTrain"); + fMatrixTrain->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrix; Training matrix was read in from file '" + << filetrain << "'" << endl; + filetr.Close(); + + + //-------------------------- + // read in test matrix + + TFile filete(filetest); + fMatrixTest->Read("MatrixTest"); + fMatrixTest->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrix; Test matrix was read in from file '" + << filetest << "'" << endl; + filete.Close(); + + return kTRUE; +} + + +//------------------------------------------------------------------------ +// +// Steering program for optimizing the supercuts +// --------------------------------------------- +// +// the criterion for the 'optimum' is +// +// - a maximum significance of the gamma signal in the alpha plot, +// in which the supercuts have been applied +// +// The various steps are : +// +// - setup the event loop to be executed for each call to fcnSupercuts +// - call TMinuit to do the minimization of (-significance) : +// the fcnSupercuts function calculates the significance +// for the current cut values +// for this - the alpha plot is produced in the event loop, +// in which the cuts have been applied +// - the significance of the gamma signal in the alpha plot +// is calculated +// +// Needed as input : (to be set by the Set functions) +// +// - fFilenameParam name of file to which optimum values of the +// parameters are written +// +// - fHadronnessName name of container where MSupercutsCalc will +// put the supercuts hadronness +// +// - for the minimization, the starting values of the parameters are taken +// - from file parSCinit (if it is != "") +// - or from the arrays params and/or steps +// +//---------------------------------------------------------------------- +Bool_t MFindSupercuts::FindParams(TString parSCinit, + TArrayD ¶ms, TArrayD &steps) +{ + // Setup the event loop which will be executed in the + // fcnSupercuts function of MINUIT + // + // parSCinit is the name of the file containing the initial values + // of the parameters; + // if parSCinit = "" 'params' and 'steps' are taken as initial values + + *fLog << "=============================================" << endl; + *fLog << "Setup event loop for fcnSupercuts" << endl; + + if (fHadronnessName.IsNull()) + { + *fLog << "MFindSupercuts::FindParams; hadronness name is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain == NULL) + { + *fLog << "MFindSupercuts::FindParams; training matrix is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain->GetM().GetNrows() <= 0) + { + *fLog << "MFindSupercuts::FindParams; training matrix has no entries" + << endl; + return kFALSE; + } + + MParList parlistfcn; + MTaskList tasklistfcn; + + // loop over rows of matrix + MMatrixLoop loop(fMatrixTrain); + + //-------------------------------- + // create container for the supercut parameters + // and set them to their initial values + MSupercuts super; + + // take initial values from file parSCinit + if (parSCinit != "") + { + TFile inparam(parSCinit); + super.Read("MSupercuts"); + inparam.Close(); + *fLog << "MFindSupercuts::FindParams; initial values of parameters are taken from file " + << parSCinit << endl; + } + + // take initial values from 'params' and/or 'steps' + else if (params.GetSize() != 0 || steps.GetSize() != 0 ) + { + if (params.GetSize() != 0) + { + *fLog << "MFindSupercuts::FindParams; initial values of parameters are taken from 'params'" + << endl; + super.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MFindSupercuts::FindParams; initial step sizes are taken from 'steps'" + << endl; + super.SetStepsizes(steps); + } + } + else + { + *fLog << "MFindSupercuts::FindParams; initial values and step sizes are taken from the MSupercuts constructor" + << endl; + } + + + //-------------------------------- + // calculate supercuts hadronness + fCalcHadTrain->SetHadronnessName(fHadronnessName); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue supercuts(&scfilter); + + // plot |alpha| + const TString mh3Name = "AlphaFcn"; + MBinning binsalpha("Binning"+mh3Name); + binsalpha.SetEdges(54, -12.0, 96.0); + + *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl; + + // |alpha| is assumed to be in column 7 of the matrix + MH3 alpha("MatrixTrain[7]"); + alpha.SetName(mh3Name); + + MFillH fillalpha(&alpha); + + // book histograms to be filled during the optimization + // ! not in the event loop ! + MHSupercuts plotsuper; + plotsuper.SetupFill(&parlistfcn); + + //****************************** + // entries in MParList (extension of old MParList) + + parlistfcn.AddToList(&tasklistfcn); + parlistfcn.AddToList(&super); + parlistfcn.AddToList(fCam); + parlistfcn.AddToList(fMatrixTrain); + + parlistfcn.AddToList(&binsalpha); + parlistfcn.AddToList(&alpha); + + parlistfcn.AddToList(&plotsuper); + + //****************************** + // entries in MTaskList + + tasklistfcn.AddToList(&loop); + tasklistfcn.AddToList(fCalcHadTrain); + tasklistfcn.AddToList(&supercuts); + tasklistfcn.AddToList(&fillalpha); + + + //****************************** + + MEvtLoop evtloopfcn("EvtLoopFCN"); + evtloopfcn.SetParList(&parlistfcn); + *fLog << "Event loop for fcnSupercuts has been setup" << endl; + + // address of evtloopfcn is used in CallMinuit() + + + //----------------------------------------------------------------------- + // + //---------- Start of minimization part -------------------- + // + // Do the minimization with MINUIT + // + // Be careful: This is not thread safe + // + *fLog << "========================================================" << endl; + *fLog << "Start minimization for supercuts" << endl; + + + // ------------------------------------------- + // prepare call to MINUIT + // + + // get initial values of parameters + fVinit = super.GetParameters(); + fStep = super.GetStepsizes(); + + TString name[fVinit.GetSize()]; + fStep.Set(fVinit.GetSize()); + fLimlo.Set(fVinit.GetSize()); + fLimup.Set(fVinit.GetSize()); + fFix.Set(fVinit.GetSize()); + + fNpar = fVinit.GetSize(); + + for (UInt_t i=0; i= 48) + // { + // fStep[i] = 0.0; + // fFix[i] = 1; + // } + //} + + + // ------------------------------------------- + // call MINUIT + + TStopwatch clock; + clock.Start(); + + *fLog << "before calling CallMinuit" << endl; + + MMinuitInterface inter; + Bool_t rc = inter.CallMinuit(fcnSupercuts, name, + fVinit, fStep, fLimlo, fLimup, fFix, + &evtloopfcn, "SIMPLEX", kFALSE); + + *fLog << "after calling CallMinuit" << endl; + + *fLog << "Time spent for the minimization in MINUIT : " << endl;; + clock.Stop(); + clock.Print(); + + plotsuper.DrawClone(); + + if (!rc) + return kFALSE; + + *fLog << "Minimization for supercuts finished" << endl; + *fLog << "========================================================" << endl; + + + // ----------------------------------------------------------------- + // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put + // into MSupercuts + + // write optimum parameter values onto file filenameParam + + TFile outparam(fFilenameParam, "RECREATE"); + super.Write(); + outparam.Close(); + + *fLog << "Optimum parameter values for supercuts were written onto file '" + << fFilenameParam << "' :" << endl; + + const TArrayD &check = super.GetParameters(); + for (Int_t i=0; iGetM().GetNrows() <= 0) + { + *fLog << "MFindSupercuts::TestParams; test matrix has no entries" + << endl; + return kFALSE; + } + + // ------------- TEST the supercuts ------------------------------ + // + *fLog << "Test the supercuts on the test sample" << endl; + + // ----------------------------------------------------------------- + // read optimum parameter values from file filenameParam + // into array 'supercutsPar' + + TFile inparam(fFilenameParam); + MSupercuts scin; + scin.Read("MSupercuts"); + inparam.Close(); + + *fLog << "Optimum parameter values for supercuts were read from file '"; + *fLog << fFilenameParam << "' :" << endl; + + const TArrayD &supercutsPar = scin.GetParameters(); + for (Int_t i=0; iSetHadronnessName(fHadronnessName); + + + //------------------------------------------- + + MMatrixLoop loop(fMatrixTest); + + // plot alpha before applying the supercuts + const TString mh3NameB = "AlphaBefore"; + MBinning binsalphabef("Binning"+mh3NameB); + binsalphabef.SetEdges(54, -12.0, 96.0); + + // |alpha| is assumed to be in column 7 of the matrix + MH3 alphabefore("MatrixTest[7]"); + alphabefore.SetName(mh3NameB); + + TH1 &alphahistb = alphabefore.GetHist(); + alphahistb.SetName("alphaBefore-TestParams"); + + MFillH fillalphabefore(&alphabefore); + fillalphabefore.SetName("FillAlphaBefore"); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue applysupercuts(&scfilter); + + // plot alpha after applying the supercuts + const TString mh3NameA = "AlphaAfter"; + MBinning binsalphaaft("Binning"+mh3NameA); + binsalphaaft.SetEdges(54, -12.0, 96.0); + + MH3 alphaafter("MatrixTest[7]"); + alphaafter.SetName(mh3NameA); + + TH1 &alphahista = alphaafter.GetHist(); + alphahista.SetName("alphaAfter-TestParams"); + + MFillH fillalphaafter(&alphaafter); + fillalphaafter.SetName("FillAlphaAfter"); + + //****************************** + // entries in MParList + + parlist2.AddToList(&tasklist2); + + parlist2.AddToList(&supercuts); + + parlist2.AddToList(fCam); + parlist2.AddToList(fMatrixTest); + + parlist2.AddToList(&binsalphabef); + parlist2.AddToList(&alphabefore); + + parlist2.AddToList(&binsalphaaft); + parlist2.AddToList(&alphaafter); + + //****************************** + // entries in MTaskList + + tasklist2.AddToList(&loop); + tasklist2.AddToList(&fillalphabefore); + + tasklist2.AddToList(fCalcHadTest); + tasklist2.AddToList(&applysupercuts); + + tasklist2.AddToList(&fillalphaafter); + + //****************************** + + MProgressBar bar2; + MEvtLoop evtloop2; + evtloop2.SetParList(&parlist2); + evtloop2.SetName("EvtLoopTestParams"); + evtloop2.SetProgressBar(&bar2); + Int_t maxevents2 = -1; + if (!evtloop2.Eventloop(maxevents2)) + return kFALSE; + + tasklist2.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // draw the alpha plots + + MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3"); + TH1 &alphaHist1 = alphaBefore->GetHist(); + alphaHist1.SetXTitle("|\\alpha| [\\circ]"); + + MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); + TH1 &alphaHist2 = alphaAfter->GetHist(); + alphaHist2.SetXTitle("|\\alpha| [\\circ]"); + alphaHist2.SetName("alpha-TestParams"); + + TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); + c->Divide(2,1); + + gROOT->SetSelectedPad(NULL); + + c->cd(1); + alphaHist1.DrawCopy(); + + c->cd(2); + alphaHist2.DrawCopy(); + + + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + const Double_t alphasig = 20.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 2; + const Bool_t drawpoly = kTRUE; + const Bool_t fitgauss = kTRUE; + const Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + + const Double_t significance = findsig.GetSignificance(); + const Double_t alphasi = findsig.GetAlphasi(); + + *fLog << "Significance of gamma signal after supercuts in test sample : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + //------------------------------------------- + + + *fLog << "Test of supercuts part finished" << endl; + *fLog << "======================================================" << endl; + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h (revision 2663) @@ -0,0 +1,134 @@ +#ifndef MARS_MFindSupercuts +#define MARS_MFindSupercuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif + +class MFilter; +class MEvtLoop; +class MH3; +class MSupercutsCalc; +class MGeomCam; +class MHMatrix; +/* +#include "MFilter.h" +#include "MEvtLoop.h" +#include "MH3.h" +#include "MSupercutsCalc.h" +#include "MGeomCam.h" +#include "MHMatrix.h" +*/ + +class MFindSupercuts : public MParContainer +{ +private: + + TString fFilenameTrain; + TString fFilenameTest; + + Int_t fHowManyTrain; + Int_t fHowManyTest; + + Bool_t fUseOrigDistribution; + + TString fFilenameParam; + + TString fHadronnessName; + + MSupercutsCalc *fCalcHadTrain; + MSupercutsCalc *fCalcHadTest; + + MHMatrix *fMatrixTrain; + MHMatrix *fMatrixTest; + MGeomCam *fCam; + + MEvtLoop *fObjectFit; + + MFilter *fMatrixFilter; + + // to comunicate with MINUIT ----------------- + // attention : dimensions must agree with those in + // MMinuitInterface::CallMinuit() + //char fParName [80][100]; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + + UInt_t fNpar; + + TString fMethod; + + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; + //-------------------------------------------- + + +public: + MFindSupercuts(const char *name=NULL, const char *title=NULL); + ~MFindSupercuts(); + + void SetFilenameTraining(const TString &name, const Int_t howmany) + {fFilenameTrain = name; fHowManyTrain = howmany; } + + void SetFilenameTest(const TString &name, const Int_t howmany) + {fFilenameTest = name; fHowManyTest = howmany; } + + void SetFilenameParam(const TString &name) {fFilenameParam = name;} + void SetHadronnessName(const TString &name) {fHadronnessName = name;} + + void SetMatrixFilter(MFilter *filter) {fMatrixFilter = filter;} + + Bool_t DefineTrainMatrix(const TString &name, MH3 &href, + const Int_t howmany, const TString &filetrain); + + Bool_t DefineTestMatrix(const TString &name, MH3 &href, + const Int_t howmany, const TString &filetest); + + Bool_t DefineTrainTestMatrix(const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest); + + void SetUseOrigDistribution(Bool_t b); + + MHMatrix *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } + + Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); + + Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); + Bool_t TestParams(); + + ClassDef(MFindSupercuts, 1) // Class for optimization of the Supercuts +}; + +#endif + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2663) @@ -0,0 +1,1874 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 06/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPadONOFF +// +// This task applies padding such that for a given pixel and for a given +// Theta bin the resulting distribution of the pedestal sigma is identical +// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta. +// +// The number of photons, its error and the pedestal sigmas are altered. +// On average, the number of photons added is zero. +// +// The formulas used can be found in Thomas Schweizer's Thesis, +// Section 2.2.1 +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// THE FOLLOWING OPTION HAS NOT BEEN TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +// +// The padding has to be done before the image cleaning because the +// image cleaning depends on the pedestal sigmas. +// +// For random numbers gRandom is used. +// +///////////////////////////////////////////////////////////////////////////// +#include "MPadONOFF.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MSigmabar.h" +#include "MMcEvt.hxx" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MBlindPixels.h" + +#include "MRead.h" +#include "MFilterList.h" +#include "MTaskList.h" +#include "MBlindPixelCalc.h" +#include "MHBlindPixels.h" +#include "MFillH.h" +#include "MHSigmaTheta.h" +#include "MEvtLoop.h" + +ClassImp(MPadONOFF); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPadONOFF::MPadONOFF(const char *name, const char *title) +{ + fName = name ? name : "MPadONOFF"; + fTitle = title ? title : "Task for the ON-OFF padding"; + + fPadFlag = 1; + *fLog << all << "fPadFlag = " << fPadFlag << endl; + + fType = ""; + + fHSigmaTheta = NULL; + fHSigmaThetaON = NULL; + fHSigmaThetaOFF = NULL; + + fHSigmaPixTheta = NULL; + fHSigmaPixThetaON = NULL; + fHSigmaPixThetaOFF = NULL; + + fHDiffPixTheta = NULL; + fHDiffPixThetaON = NULL; + fHDiffPixThetaOFF = NULL; + + fHgON = NULL; + fHgOFF = NULL; + + fHBlindPixIdTheta = NULL; + fHBlindPixNTheta = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPadONOFF::~MPadONOFF() +{ + if (fHBlindPixIdTheta != NULL) delete fHBlindPixIdTheta; + if (fHBlindPixNTheta != NULL) delete fHBlindPixNTheta; + + if (fHSigmaTheta != NULL) delete fHSigmaTheta; + if (fHSigmaPixTheta != NULL) delete fHSigmaPixTheta; + if (fHDiffPixTheta != NULL) delete fHDiffPixTheta; + + if (fHSigmaPedestal != NULL) delete fHSigmaPedestal; + if (fHPhotons != NULL) delete fHPhotons; + if (fHNSB != NULL) delete fHNSB; + + if (fInfile != NULL) delete fInfile; +} + +// -------------------------------------------------------------------------- +// +// Merge the distributions of ON and OFF data +// +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) +// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) +// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) +// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) +// +// and define the target distributions for the padding +// +Bool_t MPadONOFF::MergeHistograms(TH2D *sigthon, TH2D *sigthoff, + TH3D *sigpixthon, TH3D *sigpixthoff, + TH3D *diffpixthon, TH3D *diffpixthoff, + TH2D *blindidthon, TH2D *blindidthoff, + TH2D *blindnthon, TH2D *blindnthoff) +{ + *fLog << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "Merge the ON and OFF histograms to obtain the target distributions for the padding" + << endl; + + //------------------------------------------------------------- + // merge the distributions of ON and OFF events + // to obtain the target distribution fHSigmaTheta + // + + Double_t eps = 1.e-10; + + fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + // get binning for fHgON and fHgOFF from sigthon + // binning in Theta + TAxis *ax = sigthon->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon->GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); + fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); + + //............ loop over Theta bins ........................... + + // hista is the normalized 1D histogram of sigmabar for ON data + // histb is the normalized 1D histogram of sigmabar for OFF data + // histc is the difference ON-OFF + + // at the beginning, histap is a copy of hista + // at the end, it will be the 1D histogram for ON data after padding + + // at the beginning, histbp is a copy of histb + // at the end, it will be the 1D histogram for OFF data after padding + + // at the beginning, histcp is a copy of histc + // at the end, it should be zero + + *fLog << all << "MergeHistograms; bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + Stat_t suma = hista->Integral(); + hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + + TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); + Stat_t sumb = histb->Integral(); + histb->Scale(1./sumb); + + TH1D *histbp = new TH1D( (const TH1D&)*histb ); + + TH1D *histc = new TH1D( (const TH1D&)*hista ); + histc->Add(histb, -1.0); + + TH1D *histcp = new TH1D( (const TH1D&)*histc ); + + + // calculate matrix g + + // fHgON[k][j] tells which fraction of ON events should be padded + // from sigma_k to sigma_j + + + // fHgOFF[k][j] tells which fraction of OFF events should be padded + // from sigma_k to sigma_j + + //-------- start j loop ------------------------------------------------ + // loop over bins in histc, starting from the end + Double_t v, s, w, t, x, u, a, b, c, arg; + + for (Int_t j=nbinssig; j >= 1; j--) + { + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) continue; + if (v >= 0.0) + s = 1.0; + else + s = -1.0; + + //................ start k loop ...................................... + // look for a bin k which may compensate the content of bin j + for (Int_t k=j-1; k >= 1; k--) + { + w = histcp->GetBinContent(k); + if ( fabs(w) < eps ) continue; + if (w >= 0.0) + t = 1.0; + else + t = -1.0; + + if (s==t) continue; + + x = v + w; + if (x >= 0.0) + u = 1.0; + else + u = -1.0; + + if (u == s) + { + arg = -w; + + if (arg <=0.0) + { + fHgON->SetBinContent (l, k, j, -arg); + fHgOFF->SetBinContent(l, k, j, 0.0); + } + else + { + fHgON->SetBinContent (l, k, j, 0.0); + fHgOFF->SetBinContent(l, k, j, arg); + } + + + *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j + << ", " << arg << endl; + + // g(k-1, j-1) = arg; + //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " + // << arg << endl; + + //...................................... + // this is for checking the procedure + if (arg < 0.0) + { + a = histap->GetBinContent(k); + histap->SetBinContent(k, a+arg); + a = histap->GetBinContent(j); + histap->SetBinContent(j, a-arg); + } + else + { + b = histbp->GetBinContent(k); + histbp->SetBinContent(k, b-arg); + b = histbp->GetBinContent(j); + histbp->SetBinContent(j, b+arg); + } + //...................................... + + histcp->SetBinContent(k, 0.0); + histcp->SetBinContent(j, x); + + //...................................... + // redefine v + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) break; + if (v >= 0.0) + s = 1.0; + else + s = -1.0; + //...................................... + + continue; + } + + arg = v; + + if (arg <=0.0) + { + fHgON->SetBinContent (l, k, j, -arg); + fHgOFF->SetBinContent(l, k, j, 0.0); + } + else + { + fHgON->SetBinContent (l, k, j, 0.0); + fHgOFF->SetBinContent(l, k, j, arg); + } + + *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j + << ", " << arg << endl; + + //g(k-1, j-1) = arg; + //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " + // << arg << endl; + + //...................................... + // this is for checking the procedure + if (arg < 0.0) + { + a = histap->GetBinContent(k); + histap->SetBinContent(k, a+arg); + a = histap->GetBinContent(j); + histap->SetBinContent(j, a-arg); + } + else + { + b = histbp->GetBinContent(k); + + histbp->SetBinContent(k, b-arg); + b = histbp->GetBinContent(j); + histbp->SetBinContent(j, b+arg); + } + //...................................... + + histcp->SetBinContent(k, x); + histcp->SetBinContent(j, 0.0); + + break; + } + //................ end k loop ...................................... + } + //-------- end j loop ------------------------------------------------ + + // check results for this Theta bin + for (Int_t j=1; j<=nbinssig; j++) + { + a = histap->GetBinContent(j); + b = histbp->GetBinContent(j); + c = histcp->GetBinContent(j); + + if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) + *fLog << err << "Read; inconsistency in results; a, b, c = " + << a << ", " << b << ", " << c << endl; + } + + + // fill target distribution SigmaTheta + // for this Theta bin + // + for (Int_t j=1; j<=nbinssig; j++) + { + a = histap->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, a); + } + + delete hista; + delete histb; + delete histc; + + delete histap; + delete histbp; + delete histcp; + } + //............ end of loop over Theta bins .................... + + + // target distributions for MC + // SigmaPixTheta and DiffPixTheta + // BlindPixIdTheta and BlindPixNTheta + // are calculated as averages of the ON and OFF distributions + + // Note : in all histograms the pedestal RMS is assumed to be given in units of "number of photons" + + fHSigmaThetaON = sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (target)"); + + fHSigmaThetaOFF = sigthoff; + fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (target)"); + + + fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&)*sigpixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHSigmaPixThetaON = sigpixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (target)"); + + fHSigmaPixThetaOFF = sigpixthoff; + fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (target)"); + + fHDiffPixThetaOFF = diffpixthoff; + fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (target)"); + + + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/OFF events to be padded to a sigmabar + // of the OFF/ON events : fracON, fracOFF + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + TAxis *ax; + TAxis *ay; + TAxis *az; + + Double_t entON; + Double_t entOFF; + + Double_t normON; + Double_t normOFF; + + // set ranges for 2D-projections of 3D-histograms + ax = sigpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = sigpixthoff->GetXaxis(); + ax->SetRange(j, j); + + ax = diffpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthoff->GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histOFF; + + TH1D *hist1; + TH1D *hist1OFF; + + // weights for ON and OFF distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracOFF); + Double_t facOFF = 0.5 * (1. + fracON - fracOFF); + + *fLog << all << fracON << ", " << fracOFF << ", " + << facON << ", " << facOFF << endl; + //------------------------------------------- + ay = blindnthon->GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = (TH1D*)blindnthon->ProjectionY ("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = (TH1D*)blindnthoff->ProjectionY("", j, j, ""); + + // number of events in this theta bin + entON = hist1->Integral(); + entOFF = hist1OFF->Integral(); + + *fLog << all << "BlindPixNTheta : j, entON, entOFF = " << j << ", " + << entON << ", " << entOFF << endl; + + normON = entON<=0.0 ? 0.0 : 1.0/entON; + normOFF = entOFF<=0.0 ? 0.0 : 1.0/entOFF; + + hist1->Scale(normON); + hist1OFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist1->Scale(facON); + hist1->Add(hist1OFF, facOFF); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------- + ay = blindidthon->GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = (TH1D*)blindidthon->ProjectionY ("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = (TH1D*)blindidthoff->ProjectionY("", j, j, ""); + + hist1->Scale(normON); + hist1OFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist1->Scale(facON); + hist1->Add(hist1OFF, facOFF); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------- + ay = sigpixthon->GetYaxis(); + Int_t nbinspix = ay->GetNbins(); + + az = sigpixthon->GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigpixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)sigpixthoff->Project3D("zy"); + + hist->Scale(normON); + histOFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspix; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------- + ay = diffpixthon->GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + az = diffpixthon->GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)diffpixthoff->Project3D("zy"); + + hist->Scale(normON); + histOFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------- + } + + + *fLog << all << "The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaTheta->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOFF->SetDirectory(NULL); + + fHSigmaPixTheta->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + fHSigmaPixThetaOFF->SetDirectory(NULL); + + fHDiffPixTheta->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetDirectory(NULL); + + + fHgON->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Read target distributions from a file +// +// +Bool_t MPadONOFF::ReadTargetDist(const char* namefilein) +{ + *fLog << all << "Read padding histograms from file " << namefilein << endl; + + fInfile = new TFile(namefilein); + fInfile->ls(); + + //------------------------------------ + + fHSigmaTheta = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + *fLog << all << "Object '2D-ThetaSigmabar' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabar' was read in" << endl; + + fHSigmaThetaON = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); + if (!fHSigmaThetaON) + { + *fLog << all << "Object '2D-ThetaSigmabarON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabarON' was read in" << endl; + + fHSigmaThetaOFF = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); + if (!fHSigmaThetaOFF) + { + *fLog << "Object '2D-ThetaSigmabarOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabarOFF' was read in" << endl; + + fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + *fLog << all << "Object '3D-ThetaPixSigma' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixSigma' was read in" << endl; + + fHSigmaPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); + if (!fHSigmaPixThetaON) + { + *fLog << all << "Object '3D-ThetaPixSigmaON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixSigmaON' was read in" << endl; + + fHSigmaPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); + if (!fHSigmaPixThetaOFF) + { + *fLog << all << "Object '3D-ThetaPixSigmaOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixSigmaOFF' was read in" << endl; + + fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + *fLog << all << "Object '3D-ThetaPixDiff' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixDiff' was read in" << endl; + + fHDiffPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); + if (!fHDiffPixThetaON) + { + *fLog << all << "Object '3D-ThetaPixDiffON' not found on root file" + << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixDiffON' was read in" << endl; + + fHDiffPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); + if (!fHDiffPixThetaOFF) + { + *fLog << all << "Object '3D-ThetaPixDiffOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixDiffOFF' was read in" << endl; + + fHgON = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); + if (!fHgON) + { + *fLog << all << "Object '3D-PaddingMatrixON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-PaddingMatrixON' was read in" << endl; + + fHgOFF = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); + if (!fHgOFF) + { + *fLog << all << "Object '3D-PaddingMatrixOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-PaddingMatrixOFF' was read in" << endl; + + + fHBlindPixIdTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindId"); + if (!fHBlindPixIdTheta) + { + *fLog << all << "Object '2D-ThetaBlindId' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaBlindId' was read in" << endl; + + fHBlindPixNTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindN"); + if (!fHBlindPixNTheta) + { + *fLog << all << "Object '2D-ThetaBlindN' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaBlindN' was read in" << endl; + + + //------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Write target distributions onto a file +// +// +Bool_t MPadONOFF::WriteTargetDist(const char* namefileout) +{ + *fLog << all << "Write padding histograms onto file " << namefileout << endl; + + TFile outfile(namefileout, "RECREATE"); + + fHBlindPixNTheta->Write(); + fHBlindPixIdTheta->Write(); + + fHSigmaTheta->Write(); + fHSigmaThetaON->Write(); + fHSigmaThetaOFF->Write(); + + fHSigmaPixTheta->Write(); + fHSigmaPixThetaON->Write(); + fHSigmaPixThetaOFF->Write(); + + fHDiffPixTheta->Write(); + fHDiffPixThetaON->Write(); + fHDiffPixThetaOFF->Write(); + + fHgON->Write(); + fHgOFF->Write(); + + *fLog << all << "WriteTargetDist; target histograms written onto file " + << namefileout << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set type of data to be padded +// +// this is not necessary if the type of the data can be recognized +// directly from the input +// +// +void MPadONOFF::SetDataType(const char* type) +{ + fType = type; + *fLog << all << "SetDataType(); type of data to be padded : " + << fType << endl; + + return; +} + + +// -------------------------------------------------------------------------- +// +// Set the option for the padding +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// NOT YET TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +void MPadONOFF::SetPadFlag(Int_t padflag) +{ + fPadFlag = padflag; + *fLog << all << "SetPadFlag(); choose option " << fPadFlag << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers and prepare the histograms +// +Int_t MPadONOFF::PreProcess(MParList *pList) +{ + if ( !fHSigmaTheta || !fHSigmaThetaON || !fHSigmaThetaOFF || + !fHSigmaPixTheta || !fHSigmaPixThetaON || !fHSigmaPixThetaOFF || + !fHDiffPixTheta || !fHDiffPixThetaON || !fHDiffPixThetaOFF || + !fHBlindPixIdTheta || !fHBlindPixNTheta || + !fHgON || !fHgOFF) + { + *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << err << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << err << "MSigmabar not found... aborting." << endl; + return kFALSE; + } + + fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlinds) + { + *fLog << err << "MBlindPixels not found... aborting." << endl; + return kFALSE; + } + + if (fType !="ON" && fType !="OFF" && fType !="MC") + { + *fLog << err << "Type of data to be padded not defined... aborting." << endl; + return kFALSE; + } + + + //-------------------------------------------------------------------- + // histograms for checking the padding + // + // plot pedestal sigmas + fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", + 100, 0.0, 5.0, 100, 0.0, 5.0); + fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); + fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); + + // plot no.of photons (before vs. after padding) + fHPhotons = new TH2D("Photons","Photons: after vs.before padding", + 100, -10.0, 90.0, 100, -10, 90); + fHPhotons->SetXTitle("no.of photons before padding"); + fHPhotons->SetYTitle("no.of photons after padding"); + + // plot of added NSB + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB->SetXTitle("no.of added NSB photons"); + fHNSB->SetYTitle("no.of pixels"); + + + //-------------------------------------------------------------------- + + fIter = 20; + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the Padding +// idealy the events to be padded should have been generated without NSB +// +Int_t MPadONOFF::Process() +{ + *fLog << all << "Entry Process();" << endl; + + // this is the conversion factor from the number of photons + // to the number of photo electrons + // later to be taken from MCalibration + // fPEpePhoton = PW * LG * QE * 1D + Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + + //------------------------------------------------ + // add pixels to MCerPhotEvt which are not yet in; + // set their number of photons equal to zero + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetNumPixels(); + for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); + } + + + + //----------------------------------------- + Int_t rc=0; + Int_t rw=0; + + const UInt_t npix = fEvt->GetNumPixels(); + + //fSigmabar->Calc(*fCam, *fPed, *fEvt); + // *fLog << all << "before padding : " << endl; + //fSigmabar->Print(""); + + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + // to simulate the situation that before the padding the NSB and + // electronic noise are zero : set Sigma = 0 for all pixels + //for (UInt_t i=0; ioperator[](i); + // Int_t j = pix.GetPixId(); + + // MPedestalPix &ppix = fPed->operator[](j); + // ppix.SetMeanRms(0.0); + //} + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + //------------------------------------------- + // Calculate average sigma of the event + // + Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "before padding : " << endl; + fSigmabar->Print(""); + Double_t sigbarold = sigbarold_phot * fPEperPhoton; + + + Double_t sigbarold2 = sigbarold*sigbarold; + + + // for MC data : expect sigmabar to be zero before the padding + if (fType == "MC") + { + if (sigbarold_phot > 0) + { + *fLog << warn << "Process(); sigmabar of event to be padded is > 0 : " + << sigbarold_phot << ". Stop event loop " << endl; + // input data should have sigmabar = 0; stop event loop + + rc = 1; + fErrors[rc]++; + return kCONTINUE; + } + } + + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + *fLog << all << "theta = " << theta << endl; + + Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) + { + *fLog << warn << "Process(); binNumber out of range : theta, binTheta = " + << theta << ", " << binTheta << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + //------------------------------------------- + // get number of events in this theta bin + // and number of events in this sigbarold_phot bin + + Double_t nTheta; + Double_t nSigma; + if (fType == "ON") + { + TH1D *hn; + + hn = fHSigmaThetaON->ProjectionY("", binTheta, binTheta, ""); + nTheta = hn->Integral(); + Int_t binSigma = hn->FindBin(sigbarold_phot); + nSigma = hn->GetBinContent(binSigma); + + *fLog << all + << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + << theta << ", " << sigbarold_phot << ", " << binTheta << ", " + << binSigma << ", " << nTheta << ", " << nSigma << endl; + + delete hn; + } + + else if (fType == "OFF") + { + TH1D *hn; + + hn = fHSigmaThetaOFF->ProjectionY("", binTheta, binTheta, ""); + nTheta = hn->Integral(); + Int_t binSigma = hn->FindBin(sigbarold_phot); + nSigma = hn->GetBinContent(binSigma); + + *fLog << all + << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + << theta << ", " << sigbarold_phot << ", " << binTheta << ", " + << binSigma << ", " << nTheta << ", " << nSigma << endl; + + delete hn; + } + + else + { + nTheta = 0.0; + nSigma = 0.0; + } + + //------------------------------------------- + // for the current theta, + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // do this only for MC data + + + if (fType == "MC") + { + + // numBlind is the number of blind pixels in this event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, ""); + if ( nblind->Integral() == 0.0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + } + delete nblind; + + //warn Code commented out due no compilation errors on Alpha OSF (tgb) + + // throw the Id of numBlind different pixels in this event + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, ""); + if ( hblind->Integral() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + *fLog << all << "idBlind = " << idBlind << endl; + } + fBlinds->SetReadyToSave(); + + delete hblind; + + } + + //****************************************************************** + // has the event to be padded ? + // if yes : to which sigmabar should it be padded ? + // + + Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); + *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " + << sigbarold_phot << endl; + + Double_t prob; + TH1D *hpad = NULL; + if (fType == "ON") + { + hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); + + //Int_t nb = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nb; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHgON = " << i << ", " + // << hpad->GetBinContent(i) << endl; + //} + + if (nSigma != 0.0) + prob = hpad->Integral() * nTheta/nSigma; + else + prob = 0.0; + + *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + << ", " << prob << endl; + } + else if (fType == "OFF") + { + hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, ""); + if (nSigma != 0.0) + prob = hpad->Integral() * nTheta/nSigma; + else + prob = 0.0; + } + else + prob = 1.0; + + if ( fType != "MC" && + (prob <= 0.0 || gRandom->Uniform() > prob) ) + { + delete hpad; + // event should not be padded + *fLog << all << "event is not padded" << endl; + + rc = 8; + fErrors[rc]++; + return kTRUE; + } + // event should be padded + *fLog << all << "event is padded" << endl; + + + //------------------------------------------- + // for the current theta, generate a sigmabar + // + // for ON/OFF data according to the matrix fHgON/OFF + // for MC data according to the histogram fHSigmaTheta + // + Double_t sigmabar_phot = 0; + Double_t sigmabar = 0; + if (fType == "ON" || fType == "OFF") + { + //Int_t nbinsx = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nbinsx; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i) + // << endl; + //} + + sigmabar_phot = hpad->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + + *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; + + delete hpad; + } + + else if (fType == "MC") + { + Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta); + + if (binTheta != bincheck) + { + cout << "The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + TH1D *hsigma; + + hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, ""); + if ( hsigma->Integral() == 0.0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete hsigma; + + rc = 3; + fErrors[rc]++; + return kCONTINUE; + } + else + { + sigmabar_phot = hsigma->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + + *fLog << all << "Theta, bin number = " << theta << ", " << binTheta + << ", sigmabar_phot = " << sigmabar_phot << endl; + } + delete hsigma; + } + + const Double_t sigmabar2 = sigmabar*sigmabar; + + //------------------------------------------- + + *fLog << all << "Process(); sigbarold, sigmabar = " + << sigbarold << ", "<< sigmabar << endl; + + // Skip event if target sigmabar is <= sigbarold + if (sigmabar <= sigbarold) + { + *fLog << all << "Process(); target sigmabar is less than sigbarold : " + << sigmabar << ", " << sigbarold << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; + return kCONTINUE; + } + + + //------------------------------------------- + // + // Calculate average number of NSB photo electrons to be added (lambdabar) + // from the value of sigmabar, + // - making assumptions about the average electronic noise (elNoise2) and + // - using a fixed value (F2excess) for the excess noise factor + + Double_t elNoise2; // [photons] + Double_t F2excess = 1.3; + Double_t lambdabar; // [photons] + + + + Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta); + if (binTheta != bincheck) + { + cout << "The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. + // The average electronic noise (to be added) has to be well above this RMS, + // otherwise the electronic noise of an individual pixel (elNoise2Pix) + // may become negative + + TH1D *hnoise; + if (fType == "MC") + hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "ON") + hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "OFF") + hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else + { + *fLog << err << "Process; illegal data type... aborting" << endl; + return kERROR; + } + + Double_t RMS_phot = hnoise->GetRMS(1); + Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; + delete hnoise; + + elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); + *fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photo electrons] + + // This value of lambdabar is the same for all pixels; + // note that lambdabar is normalized to the area of pixel 0 + + //---------- start loop over pixels --------------------------------- + // do the padding for each pixel + // + // pad only pixels - which are used (before image cleaning) + // + Double_t sig_phot = 0; + Double_t sig = 0; + + Double_t sigma2 = 0; + + Double_t diff_phot = 0; + Double_t diff = 0; + + Double_t addSig2_phot= 0; + Double_t addSig2 = 0; + + Double_t elNoise2Pix = 0; + + + for (UInt_t i=0; iGetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; + Double_t oldsigma_phot = ppix.GetPedestalRms(); + Double_t oldsigma = oldsigma_phot * fPEperPhoton; + Double_t oldsigma2 = oldsigma*oldsigma; + + //--------------------------------- + // throw the Sigma for this pixel + // + Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); + + Int_t count; + Bool_t ok; + + TH1D *hdiff; + TH1D *hsig; + + switch (fPadFlag) + { + case 1 : + // throw the Sigma for this pixel from the distribution fHDiffPixTheta + + if (fType == "MC") + hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else if (fType == "ON") + hdiff = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else + hdiff = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + + if ( hdiff->Integral() == 0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hdiff; + + rc = 5; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + diff = diff_phot * fPEperPhoton * fPEperPhoton; + + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2/ratioArea + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, diff = " << theta + << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; + + rw = 1; + fWarnings[rw]++; + } + else + { + rw = 0; + fWarnings[rw]++; + } + + delete hdiff; + sigma2 = diff + sigmabar2; + break; + + case 2 : + // throw the Sigma for this pixel from the distribution fHSigmaPixTheta + + if (fType == "MC") + hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else if (fType == "ON") + hsig = fHSigmaPixThetaOFF->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else + hsig = fHSigmaPixThetaON->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + + if ( hsig->Integral() == 0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hsig; + + rc = 6; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + sig = sig_phot * fPEperPhoton; + + sigma2 = sig*sig/ratioArea; + // the following condition ensures that elNoise2Pix > 0.0 + + if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, sig = " << theta + << ", " + << j << ", " << count << ", " << sigmabar << ", " + << sig << endl; + sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; + + rw = 1; + fWarnings[rw]++; + } + else + { + rw = 0; + fWarnings[rw]++; + } + + delete hsig; + break; + } + + //--------------------------------- + // get the additional sigma^2 for this pixel (due to the padding) + + addSig2 = sigma2*ratioArea - oldsigma2; + addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; + + + //--------------------------------- + // throw actual number of additional NSB photo electrons (NSB) + // and its RMS (sigmaNSB) + + Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); + Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; + Double_t sigmaNSB0; + + if (arg >= 0.0) + { + sigmaNSB0 = sqrt( arg ); + } + else + { + sigmaNSB0 = 0.0000001; + if (arg < -1.e-10) + { + *fLog << err << "Process(); argument of sqrt < 0 : " + << arg << endl; + } + } + + + //--------------------------------- + // smear NSB0 according to sigmaNSB0 + // and subtract lambdabar because of AC coupling + + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; + Double_t NSB_phot = NSB / fPEperPhoton; + + //--------------------------------- + + // add additional NSB to the number of photons + Double_t oldphotons_phot = pix.GetNumPhotons(); + Double_t oldphotons = oldphotons_phot * fPEperPhoton; + Double_t newphotons = oldphotons + NSB; + Double_t newphotons_phot = newphotons / fPEperPhoton;; + pix.SetNumPhotons( newphotons_phot ); + + + fHNSB->Fill( NSB_phot/ratioArea ); + fHPhotons->Fill( oldphotons_phot/ratioArea, + newphotons_phot/ratioArea ); + + + // error: add sigma of padded noise quadratically + Double_t olderror_phot = pix.GetErrorPhot(); + Double_t newerror_phot = + sqrt( olderror_phot*olderror_phot + addSig2_phot ); + pix.SetErrorPhot( newerror_phot ); + + + Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + Double_t newsigma_phot = newsigma / fPEperPhoton;; + ppix.SetPedestalRms( newsigma_phot ); + + fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); + } + //---------- end of loop over pixels --------------------------------- + + // Calculate sigmabar again and crosscheck + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "after padding : " << endl; + fSigmabar->Print(""); + + rc = 0; + fErrors[rc]++; + return kTRUE; + //****************************************************************** +} + +// -------------------------------------------------------------------------- +// +// +Int_t MPadONOFF::PostProcess() +{ + if (GetNumExecutions() != 0) + { + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + + int temp; + if (fWarnings[0] != 0.0) + temp = (int)(fWarnings[1]*100/fWarnings[0]); + else + temp = 100; + + *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) + << temp + << "%) Warning : iteration to find acceptable sigma failed" + << ", fIter = " << fIter << endl; + + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) + << (int)(fErrors[2]*100/GetNumExecutions()) + << "%) Evts skipped due to: Zenith angle out of range" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) + << (int)(fErrors[3]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigmabar" << endl; + + *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) + << (int)(fErrors[4]*100/GetNumExecutions()) + << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; + + *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) + << (int)(fErrors[5]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; + + *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) + << (int)(fErrors[6]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma" << endl; + + *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) + << (int)(fErrors[7]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Blind pixels" << endl; + + *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) + << (int)(fErrors[8]*100/GetNumExecutions()) + << "%) Evts didn't have to be padded" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts were successfully padded!" << endl; + *fLog << endl; + + } + + //--------------------------------------------------------------- + TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200)); + c.Divide(3, 4); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + fHNSB->SetDirectory(NULL); + fHNSB->DrawCopy(); + fHNSB->SetBit(kCanDelete); + + c.cd(2); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->DrawCopy(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(3); + fHPhotons->SetDirectory(NULL); + fHPhotons->DrawCopy(); + fHPhotons->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(4); + fHSigmaTheta->SetDirectory(NULL); + fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta"); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 + + c.cd(5); + TH2D *l1; + l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); + l1->SetDirectory(NULL); + l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l1->SetXTitle("\\Theta [\\circ]"); + l1->SetYTitle("Sigma^2-Sigmabar^2"); + + l1->DrawCopy("box"); + l1->SetBit(kCanDelete);; + + c.cd(8); + TH2D *l2; + l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); + l2->SetDirectory(NULL); + l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l2->SetXTitle("pixel"); + l2->SetYTitle("Sigma^2-Sigmabar^2"); + + l2->DrawCopy("box"); + l2->SetBit(kCanDelete);; + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma + + c.cd(6); + TH2D *k1; + k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); + k1->SetDirectory(NULL); + k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)"); + k1->SetXTitle("\\Theta [\\circ]"); + k1->SetYTitle("Sigma"); + + k1->DrawCopy("box"); + k1->SetBit(kCanDelete); + + c.cd(9); + TH2D *k2; + k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); + k2->SetDirectory(NULL); + k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)"); + k2->SetXTitle("pixel"); + k2->SetYTitle("Sigma"); + + k2->DrawCopy("box"); + k2->SetBit(kCanDelete);; + + + //-------------------------------------------------------------------- + + c.cd(11); + TH2D *m1; + m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); + m1->SetDirectory(NULL); + m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetXTitle("Sigmabar-old"); + m1->SetYTitle("Sigmabar-new"); + + m1->DrawCopy("box"); + m1->SetBit(kCanDelete);; + + c.cd(12); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetXTitle("Sigmabar-old"); + m2->SetYTitle("Sigmabar-new"); + + m2->DrawCopy("box"); + m2->SetBit(kCanDelete);; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2663) @@ -0,0 +1,105 @@ +#ifndef MARS_MPadONOFF +#define MARS_MPadONOFF + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class TH3D; + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MSigmabar; +class MParList; +class MBlindPixels; +class MRead; +class MFilterList; + + +class MPadONOFF : public MTask +{ +private: + MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + MBlindPixels *fBlinds; + + TString fType; // type of data to be padded + TFile *fInfile; // input file containing padding histograms + + Int_t fPadFlag; + Int_t fIter; + + Int_t fErrors[9]; + Int_t fWarnings[2]; + + + // plots used for the padding + // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" + TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) + + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) + + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) + + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + + TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data + TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data + + // plots for checking the padding + TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after + // versus before padding + TH2D *fHPhotons; // 2D-histogram : no.of photons after + // versus before padding + TH1D *fHNSB; // 1D-histogram : additional NSB + + +public: + MPadONOFF(const char *name=NULL, const char *title=NULL); + ~MPadONOFF(); + + Bool_t MergeHistograms(TH2D *sigthon, TH2D *sigthoff, + TH3D *sigpixthon, TH3D *sigpixthoff, + TH3D *diffpixthon, TH3D *diffpixthoff, + TH2D *blindidthon, TH2D *blindidthoff, + TH2D *blindnthon, TH2D *blindnthoff); + + Bool_t ReadTargetDist(const char *filein); + Bool_t WriteTargetDist(const char *fileout); + + void SetDataType(const char *type); // type of data to be padded + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + void SetPadFlag(Int_t padflag); + + ClassDef(MPadONOFF, 0) // task for the ON-OFF padding +}; + +#endif + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2663) @@ -0,0 +1,114 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Sabrina Stark 12/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MPedestalWorkaround // +// // +// This class contains the function to store the pedestal value and the // +// RMS of the pedestal in units of photons. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MPedestalWorkaround.h" + +#include +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MPedPhotCam.h" + +using namespace std; + +ClassImp(MPedestalWorkaround); + +MPedestalWorkaround::MPedestalWorkaround(const char *name, const char *title) +{ + fName = name ? name : "MPedestalWorkaround"; + fTitle = title ? title : "Storage of pedestal values and RMS in units of photons"; + +} + +// ------------------------------------------------------------------------ +// + +Int_t MPedestalWorkaround::PreProcess(MParList *pList) +{ + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + fPedPhot = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); + if (!fPedPhot) + { + *fLog << err << "MPedPhotCam not found... aborting." << endl; + return kFALSE; + } + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +Int_t MPedestalWorkaround::Process() +{ + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetPixelContent( val, i, *fCam, type); + valout = (*fPed)[i].GetPedestal(); + (*fPed)[i].SetPedestal(val); + *fLog << "i, val, valout : " << i <<", "<< val<<", " << valout << endl; + type = 1; + fPedPhot->GetPixelContent( val, i, *fCam, type); + valout = (*fPed)[i].GetPedestalRms(); + (*fPed)[i].SetPedestalRms(val); + *fLog << "RMS : i, val, valout : " << i <<", "<< val<<", " << valout << endl; + } + + + return kTRUE; +} + +Int_t MPedestalWorkaround::PostProcess() +{ + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2663) @@ -0,0 +1,44 @@ +#ifndef MARS_MPedestalWorkaround +#define MARS_MPedestalWorkaround + +///////////////////////////////////////////////////////////////////////////// +// // +// MPedestalWorkaround // +// // +// Temporary copy of pedestal values and RMS into MPedestalCam from // +// MPedPhotCam in order to have the pedestal RMS in units of number of // +// photons // +// // +///////////////////////////////////////////////////////////////////////////// +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MPedestalCam; +class MPedPhotCam; +class MGeomCam; +class MParList; + + +class MPedestalWorkaround : public MTask +{ +private: + + MPedestalCam *fPed; // + MPedPhotCam *fPedPhot; // + MGeomCam *fCam; // + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + MPedestalWorkaround(const char *name=NULL, const char *title=NULL); + + ClassDef(MPedestalWorkaround, 0) // Task to copy pedestal RMS +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercuts.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercuts.cc (revision 2663) @@ -0,0 +1,388 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 08/2003 +! Author(s): Thomas Bretz, 08/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSupercuts // +// // +// this is the container for the parameters of the supercuts // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSupercuts.h" + +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSupercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// constructor +// +MSupercuts::MSupercuts(const char *name, const char *title) + : fParameters(104), fStepsizes(104), + fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), + fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24), + fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), + fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), + fConcUp(fParameters.GetArray()+64), fConcLo(fParameters.GetArray()+72), + fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88), + fAlphaUp(fParameters.GetArray()+96) +{ + fName = name ? name : "MSupercuts"; + fTitle = title ? title : "Container for the supercut parameters"; + + // set supercut parameters to their default values + InitParameters(); +} + + +// -------------------------------------------------------------------------- +// +// set default values for the supercut parameters +// +void MSupercuts::InitParameters() +{ + //--------------------------------------------------- + // these are Daniel's original values for Mkn 421 + + fLengthUp[0] = 0.315585; + fLengthUp[1] = 0.001455; + fLengthUp[2] = 0.203198; + fLengthUp[3] = 0.005532; + fLengthUp[4] = -0.001670; + fLengthUp[5] = -0.020362; + fLengthUp[6] = 0.007388; + fLengthUp[7] = -0.013463; + + fLengthLo[0] = 0.151530; + fLengthLo[1] = 0.028323; + fLengthLo[2] = 0.510707; + fLengthLo[3] = 0.053089; + fLengthLo[4] = 0.013708; + fLengthLo[5] = 2.357993; + fLengthLo[6] = 0.000080; + fLengthLo[7] = -0.007157; + + fWidthUp[0] = 0.145412; + fWidthUp[1] = -0.001771; + fWidthUp[2] = 0.054462; + fWidthUp[3] = 0.022280; + fWidthUp[4] = -0.009893; + fWidthUp[5] = 0.056353; + fWidthUp[6] = 0.020711; + fWidthUp[7] = -0.016703; + + fWidthLo[0] = 0.089187; + fWidthLo[1] = -0.006430; + fWidthLo[2] = 0.074442; + fWidthLo[3] = 0.003738; + fWidthLo[4] = -0.004256; + fWidthLo[5] = -0.014101; + fWidthLo[6] = 0.006126; + fWidthLo[7] = -0.002849; + + fDistUp[0] = 1.787943; + fDistUp[1] = 0; + fDistUp[2] = 2.942310; + fDistUp[3] = 0.199815; + fDistUp[4] = 0; + fDistUp[5] = 0.249909; + fDistUp[6] = 0.189697; + fDistUp[7] = 0; + + fDistLo[0] = 0.589406; + fDistLo[1] = 0; + fDistLo[2] = -0.083964; + fDistLo[3] = -0.007975; + fDistLo[4] = 0; + fDistLo[5] = 0.045374; + fDistLo[6] = -0.001750; + fDistLo[7] = 0; + + + // dummy values + + fAsymUp[0] = 1.e10; + fAsymUp[1] = 0.0; + fAsymUp[2] = 0.0; + fAsymUp[3] = 0.0; + fAsymUp[4] = 0.0; + fAsymUp[5] = 0.0; + fAsymUp[6] = 0.0; + fAsymUp[7] = 0.0; + + fAsymLo[0] = -1.e10; + fAsymLo[1] = 0.0; + fAsymLo[2] = 0.0; + fAsymLo[3] = 0.0; + fAsymLo[4] = 0.0; + fAsymLo[5] = 0.0; + fAsymLo[6] = 0.0; + fAsymLo[7] = 0.0; + + fConcUp[0] = 1.e10; + fConcUp[1] = 0.0; + fConcUp[2] = 0.0; + fConcUp[3] = 0.0; + fConcUp[4] = 0.0; + fConcUp[5] = 0.0; + fConcUp[6] = 0.0; + fConcUp[7] = 0.0; + + fConcLo[0] = -1.e10; + fConcLo[1] = 0.0; + fConcLo[2] = 0.0; + fConcLo[3] = 0.0; + fConcLo[4] = 0.0; + fConcLo[5] = 0.0; + fConcLo[6] = 0.0; + fConcLo[7] = 0.0; + + fLeakage1Up[0] = 1.e10; + fLeakage1Up[1] = 0.0; + fLeakage1Up[2] = 0.0; + fLeakage1Up[3] = 0.0; + fLeakage1Up[4] = 0.0; + fLeakage1Up[5] = 0.0; + fLeakage1Up[6] = 0.0; + fLeakage1Up[7] = 0.0; + + fLeakage1Lo[0] = -1.e10; + fLeakage1Lo[1] = 0.0; + fLeakage1Lo[2] = 0.0; + fLeakage1Lo[3] = 0.0; + fLeakage1Lo[4] = 0.0; + fLeakage1Lo[5] = 0.0; + fLeakage1Lo[6] = 0.0; + fLeakage1Lo[7] = 0.0; + + fAlphaUp[0] = 13.123440; + fAlphaUp[1] = 0; + fAlphaUp[2] = 0; + fAlphaUp[3] = 0; + fAlphaUp[4] = 0; + fAlphaUp[5] = 0; + fAlphaUp[6] = 0; + fAlphaUp[7] = 0; + + //--------------------------------------------------- + // fStepsizes + // if == 0.0 the parameter will be fixed in the minimization + // != 0.0 initial step sizes for the parameters + + // LengthUp + fStepsizes[0] = 0.03; + fStepsizes[1] = 0.0002; + fStepsizes[2] = 0.02; + fStepsizes[3] = 0.0006; + fStepsizes[4] = 0.0002; + fStepsizes[5] = 0.002; + fStepsizes[6] = 0.0008; + fStepsizes[7] = 0.002; + + // LengthLo + fStepsizes[8] = 0.02; + fStepsizes[9] = 0.003; + fStepsizes[10] = 0.05; + fStepsizes[11] = 0.006; + fStepsizes[12] = 0.002; + fStepsizes[13] = 0.3; + fStepsizes[14] = 0.0001; + fStepsizes[15] = 0.0008; + + // WidthUp + fStepsizes[16] = 0.02; + fStepsizes[17] = 0.0002; + fStepsizes[18] = 0.006; + fStepsizes[19] = 0.003; + fStepsizes[20] = 0.002; + fStepsizes[21] = 0.006; + fStepsizes[22] = 0.002; + fStepsizes[23] = 0.002; + + // WidthLo + fStepsizes[24] = 0.009; + fStepsizes[25] = 0.0007; + fStepsizes[26] = 0.008; + fStepsizes[27] = 0.0004; + fStepsizes[28] = 0.0005; + fStepsizes[29] = 0.002; + fStepsizes[30] = 0.0007; + fStepsizes[31] = 0.003; + + // DistUp + fStepsizes[32] = 0.2; + fStepsizes[33] = 0.0; + fStepsizes[34] = 0.3; + fStepsizes[35] = 0.02; + fStepsizes[36] = 0.0; + fStepsizes[37] = 0.03; + fStepsizes[38] = 0.02; + fStepsizes[39] = 0.0; + + // DistLo + fStepsizes[40] = 0.06; + fStepsizes[41] = 0.0; + fStepsizes[42] = 0.009; + fStepsizes[43] = 0.0008; + fStepsizes[44] = 0.0; + fStepsizes[45] = 0.005; + fStepsizes[46] = 0.0002; + fStepsizes[47] = 0.0; + + // AsymUp + fStepsizes[48] = 0.0; + fStepsizes[49] = 0.0; + fStepsizes[50] = 0.0; + fStepsizes[51] = 0.0; + fStepsizes[52] = 0.0; + fStepsizes[53] = 0.0; + fStepsizes[54] = 0.0; + fStepsizes[55] = 0.0; + + // AsymLo + fStepsizes[56] = 0.0; + fStepsizes[57] = 0.0; + fStepsizes[58] = 0.0; + fStepsizes[59] = 0.0; + fStepsizes[60] = 0.0; + fStepsizes[61] = 0.0; + fStepsizes[62] = 0.0; + fStepsizes[63] = 0.0; + + // ConcUp + fStepsizes[64] = 0.0; + fStepsizes[65] = 0.0; + fStepsizes[66] = 0.0; + fStepsizes[67] = 0.0; + fStepsizes[68] = 0.0; + fStepsizes[69] = 0.0; + fStepsizes[70] = 0.0; + fStepsizes[71] = 0.0; + + // ConcLo + fStepsizes[72] = 0.0; + fStepsizes[73] = 0.0; + fStepsizes[74] = 0.0; + fStepsizes[75] = 0.0; + fStepsizes[76] = 0.0; + fStepsizes[77] = 0.0; + fStepsizes[78] = 0.0; + fStepsizes[79] = 0.0; + + // Leakage1Up + fStepsizes[80] = 0.0; + fStepsizes[81] = 0.0; + fStepsizes[82] = 0.0; + fStepsizes[83] = 0.0; + fStepsizes[84] = 0.0; + fStepsizes[85] = 0.0; + fStepsizes[86] = 0.0; + fStepsizes[87] = 0.0; + + // Leakage1Lo + fStepsizes[88] = 0.0; + fStepsizes[89] = 0.0; + fStepsizes[90] = 0.0; + fStepsizes[91] = 0.0; + fStepsizes[92] = 0.0; + fStepsizes[93] = 0.0; + fStepsizes[94] = 0.0; + fStepsizes[95] = 0.0; + + // AlphaUp + fStepsizes[96] = 0.0; + fStepsizes[97] = 0.0; + fStepsizes[98] = 0.0; + fStepsizes[99] = 0.0; + fStepsizes[100] = 0.0; + fStepsizes[101] = 0.0; + fStepsizes[102] = 0.0; + fStepsizes[103] = 0.0; +} + + +// -------------------------------------------------------------------------- +// +// Set the parameter values from the array 'd' +// +// +Bool_t MSupercuts::SetParameters(const TArrayD &d) +{ + if (d.GetSize() != fParameters.GetSize()) + { + *fLog << err << "Sizes of d and of fParameters are different : " + << d.GetSize() << ", " << fParameters.GetSize() << endl; + return kFALSE; + } + + fParameters = d; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the step sizes from the array 'd' +// +// +Bool_t MSupercuts::SetStepsizes(const TArrayD &d) +{ + if (d.GetSize() != fStepsizes.GetSize()) + { + *fLog << err << "Sizes of d and of fStepsizes are different : " + << d.GetSize() << ", " << fStepsizes.GetSize() << endl; + return kFALSE; + } + + fStepsizes = d; + + return kTRUE; +} + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercuts.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercuts.h (revision 2663) @@ -0,0 +1,69 @@ +#ifndef MARS_MSupercuts +#define MARS_MSupercuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MSupercuts : public MParContainer +{ +private: + TArrayD fParameters; // supercut parameters + TArrayD fStepsizes; // step sizes of supercut parameters + + Double_t *fLengthUp; //! + Double_t *fLengthLo; //! + Double_t *fWidthUp; //! + Double_t *fWidthLo; //! + Double_t *fDistUp; //! + Double_t *fDistLo; //! + Double_t *fAsymUp; //! + Double_t *fAsymLo; //! + + Double_t *fConcUp; //! + Double_t *fConcLo; //! + Double_t *fLeakage1Up; //! + Double_t *fLeakage1Lo; //! + + Double_t *fAlphaUp; //! + + +public: + MSupercuts(const char *name=NULL, const char *title=NULL); + + void InitParameters(); + + Bool_t SetParameters(const TArrayD &d); + Bool_t SetStepsizes(const TArrayD &d); + + const TArrayD &GetParameters() const { return fParameters; } + const TArrayD &GetStepsizes() const { return fStepsizes; } + + const Double_t *GetLengthUp() const { return fLengthUp; } + const Double_t *GetLengthLo() const { return fLengthLo; } + const Double_t *GetWidthUp() const { return fWidthUp; } + const Double_t *GetWidthLo() const { return fWidthLo; } + const Double_t *GetDistUp() const { return fDistUp; } + const Double_t *GetDistLo() const { return fDistLo; } + const Double_t *GetAsymUp() const { return fAsymUp; } + const Double_t *GetAsymLo() const { return fAsymLo; } + + const Double_t *GetConcUp() const { return fConcUp; } + const Double_t *GetConcLo() const { return fConcLo; } + + const Double_t *GetLeakage1Up() const { return fLeakage1Up; } + const Double_t *GetLeakage1Lo() const { return fLeakage1Lo; } + + const Double_t *GetAlphaUp() const { return fAlphaUp; } + + ClassDef(MSupercuts, 1) // A container for the Supercut parameters +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc (revision 2663) @@ -0,0 +1,348 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSupercutsCalc // +// // +// this class calculates the hadronness for the supercuts // +// the parameters of the supercuts are taken // +// from the container MSupercuts // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSupercutsCalc.h" + +#include +#include + +#include "TFile.h" +#include "TArrayD.h" + +#include "MParList.h" +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" +#include "MMcEvt.hxx" +#include "MCerPhotEvt.h" +#include "MGeomCam.h" +#include "MHadronness.h" +#include "MHMatrix.h" +#include "MSupercuts.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSupercutsCalc); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// constructor +// + +MSupercutsCalc::MSupercutsCalc(const char *hilname, + const char *hilsrcname, + const char *name, const char *title) + : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname), + fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), + fSuperName("MSupercuts") +{ + fName = name ? name : "MSupercutsCalc"; + fTitle = title ? title : "Class to evaluate the Supercuts"; + + fMatrix = NULL; +} + +// -------------------------------------------------------------------------- +// +Int_t MSupercutsCalc::PreProcess(MParList *pList) +{ + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); + if (!fHadronness) + { + *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl; + return kFALSE; + } + + fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts"); + if (!fSuper) + { + *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl; + return kFALSE; + } + + if (fMatrix) + return kTRUE; + + //----------------------------------------------------------- + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); + if (!fHil) + { + *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt"); + if (!fHilExt) + { + *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar"); + if (!fNewPar) + { + *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculation of upper and lower limits +// +Double_t MSupercutsCalc::CtsMCut(const Double_t* a, Double_t ls, Double_t ct, + Double_t ls2, Double_t dd2) const +{ + // define cut-function + // + // dNOMLOGSIZE = 4.1 (=log(60.0) + // dNOMCOSZA = 1.0 + // + // a: array of cut parameters + // ls: log(SIZE) - dNOMLOGSIZE + // ls2: ls^2 + // ct: Cos(ZA.) - dNOMCOSZA + // dd2: DIST^2 + const Double_t limit = + a[0] + a[1] * dd2 + a[2] * ct + + ls * (a[3] + a[4] * dd2 + a[5] * ct) + + ls2 * (a[6] + a[7] * dd2); + + //*fLog << "MSupercutsCalc::CtsMCut; *a = " + // << *a << ", " << *(a+1) << ", " << *(a+2) << ", " + // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", " + // << *(a+6) << ", " << *(a+7) << endl; + + //*fLog << "MSupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls + // << ", " << ls2 << ", " << ct << ", " << dd2 << ", " + // << limit << endl; + + return limit; +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MSupercutsCalc::GetVal(Int_t i) const +{ + + Double_t val = (*fMatrix)[fMap[i]]; + + + //*fLog << "MSupercutsCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + + return val; +} + +// -------------------------------------------------------------------------- +// +// You can use this function if you want to use a MHMatrix instead of the +// given containers. This function adds all necessary columns to the +// given matrix. Afterward you should fill the matrix with the corresponding +// data (eg from a file by using MHMatrix::Fill). If you now loop +// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process +// will take the values from the matrix instead of the containers. +// +void MSupercutsCalc::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); + fMap[1] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[2] = fMatrix->AddColumn("MHillas.fLength"); + fMap[3] = fMatrix->AddColumn("MHillas.fSize"); + fMap[4] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[5] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); + fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); + fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); +} + +// --------------------------------------------------------------------------- +// +// Evaluate dynamical supercuts +// +// set hadronness to 0.25 if cuts are fullfilled +// 0.75 otherwise +// +Int_t MSupercutsCalc::Process() +{ + const Double_t kNomLogSize = 4.1; + const Double_t kNomCosZA = 1.0; + + const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta(); + const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength(); + const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize(); + const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX(); + const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY(); + const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist(); + + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), + fHilSrc->GetCosDeltaAlpha()); + const Double_t asym0 = fMatrix ? GetVal(8) : help; + const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); + const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + + const Double_t newdist = dist0 * fMm2Deg; + + const Double_t dist2 = meanx*meanx + meany*meany; + const Double_t dist = sqrt(dist2) * fMm2Deg; + const Double_t dd2 = dist*dist; + + + const Double_t dmls = log(size) - kNomLogSize; + const Double_t dmls2 = dmls * dmls; + + const Double_t dmcza = cos(theta) - kNomCosZA; + + const Double_t length = length0 * fMm2Deg; + const Double_t width = width0 * fMm2Deg; + const Double_t asym = asym0 * fMm2Deg; + + /* + *fLog << "newdist, length, width, asym, dist, conc, leakage = " + << newdist << ", " << length << ", " << width << ", " + << asym << ", " << dist << ", " << conc << ", " << leakage + << endl; + + *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", " + << endl; + */ + + + if ( + //dist < 1.05 && + //newdist < 1.05 && + + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) && + length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) && + + width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) && + width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) && + + asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + + dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) ) + + fHadronness->SetHadronness(0.25); + else + fHadronness->SetHadronness(0.75); + + //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl; + + fHadronness->SetReadyToSave(); + + return kTRUE; +} + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h (revision 2663) @@ -0,0 +1,82 @@ +#ifndef MARS_MSupercutsCalc +#define MARS_MSupercutsCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MParList; +class MHillas; +class MHillasSrc; +class MHillasExt; +class MNewImagePar; +class MMcEvt; +class MCerPhotEvt; +class MGeomCam; +class MHadronness; +class MHMatrix; +class MSupercuts; + +class MSupercutsCalc : public MTask +{ +private: + MHillas *fHil; + MHillasSrc *fHilSrc; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; + MMcEvt *fMcEvt; + MHadronness *fHadronness; //! output container for hadronness + MSupercuts *fSuper; // container for supercut parameters + + TString fHadronnessName; // name of container to store hadronness + TString fHilName; + TString fHilSrcName; + TString fHilExtName; + TString fNewParName; + TString fSuperName; // name of container for supercut parameters + + Double_t fMm2Deg; //! + + Int_t fMap[11]; //! + MHMatrix *fMatrix; //! + + Int_t PreProcess(MParList *pList); + Int_t Process(); + + Double_t GetVal(Int_t i) const; + + Double_t CtsMCut(const Double_t* a, Double_t ls, Double_t ct, + Double_t ls2, Double_t dd2) const; + +public: + MSupercutsCalc(const char *hilname="MHillas", + const char *hilsrcname="MHillasSrc", + const char *name=NULL, const char *title=NULL); + + void SetHadronnessName(const TString name) { fHadronnessName = name; } + TString GetHadronnessName() const { return fHadronnessName; } + + void InitMapping(MHMatrix *mat); + void StopMapping() { InitMapping(NULL); } + + + ClassDef(MSupercutsCalc, 0) // A class to evaluate the Supercuts +}; + +#endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2662) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2663) @@ -68,7 +68,11 @@ MCT1SupercutsCalc.cc \ MCT1FindSupercuts.cc \ + MSupercuts.cc \ + MSupercutsCalc.cc \ + MFindSupercuts.cc \ MMinuitInterface.cc \ MFiltercutsCalc.cc \ MCT1PadONOFF.cc \ + MPadONOFF.cc \ MCalibrationCalc.cc \ MCalibrationPix.cc \ @@ -76,4 +80,5 @@ MCalibrationPINDiode.cc \ MCalibrationCam.cc \ + MPedestalWorkaround.cc \ MExtractedSignalCam.cc \ MExtractedSignalPix.cc \ @@ -82,4 +87,5 @@ MArrivalTimeCalc.cc \ MMcCalibrationCalc.cc +>>>>>>> 1.65 SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 2662) +++ trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 2663) @@ -22,4 +22,7 @@ #pragma link C++ class MFCT1SelStandard+; #pragma link C++ class MFCT1SelFinal+; +#pragma link C++ class MFSelBasic+; +#pragma link C++ class MFSelStandard+; +#pragma link C++ class MFSelFinal+; #pragma link C++ class MFEnergySlope+; Index: trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc (revision 2663) @@ -0,0 +1,271 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 04/2003 +! Author(s): Thomas Bretz, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MFSelBasic +// +// This is a class to evaluate basic cuts +// +// to be called after the calibration (when the number of photons is +// available for all pixels) +// +// The basic cuts are : +// +// remove bad runs +// thetamin < theta < thetamax +// software trigger fullfilled (with minimum no.of photons = minphotons) +// +// +///////////////////////////////////////////////////////////////////////////// + +#include "MFSelBasic.h" + +#include "MParList.h" + +#include "MMcEvt.hxx" + +#include "MCerPhotEvt.h" +#include "MRawRunHeader.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MPedestalPix.h" +#include "MPedestalCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MFSelBasic); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFSelBasic::MFSelBasic(const char *name, const char *title) +{ + fName = name ? name : "MFSelBasic"; + fTitle = title ? title : "Filter to evaluate basic cuts"; + + // default values of cuts + SetCuts(13.0, 0.0, 60.0); +} + +// -------------------------------------------------------------------------- +// +// Set the cut values +// +// +void MFSelBasic::SetCuts(Float_t minphotons, + Float_t thetamin, Float_t thetamax) +{ + fMinPhotons = minphotons; + fThetaMin = thetamin; + fThetaMax = thetamax; + + *fLog << inf << "MFSelBasic cut values : fMinPhotons, fThetaMin, fThetaMax = "; + *fLog << fMinPhotons <<", " << fThetaMin << ", " << fThetaMax << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers +// +// +Int_t MFSelBasic::PreProcess(MParList *pList) +{ + fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRawRun) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } +/* + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam missing in Parameter List... aborting." << endl; + return kFALSE; + } +*/ + memset(fCut, 0, sizeof(fCut)); + + return kTRUE; +} + +Int_t MFSelBasic::Set(Int_t rc) +{ + fCut[rc]++; + fResult=kTRUE; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Evaluate basic cuts +// +// bad events : fResult = kTRUE; +// good events : fResult = kFALSE; +// +Int_t MFSelBasic::Process() +{ + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + + fResult = kFALSE; + + // remove bad runs for MC gammas + if (fMcEvt->GetEnergy() == 0.0 && fMcEvt->GetImpact() == 0.0) + { + if (fRawRun->GetRunNumber() == 601 || + fRawRun->GetRunNumber() == 613 || + fRawRun->GetRunNumber() == 614 ) + return Set(1); + } + + if (thetafThetaMax) + return Set(3); + + if (!SwTrigger()) + return Set(4); + + fCut[0]++; + + return kTRUE; +} +// -------------------------------------------------------------------------- +// +// Software trigger +// +// require 2 neighboring pixels (which are not in the outermost ring), +// each having at least 'fMinPhotons' photons +// +// +Bool_t MFSelBasic::SwTrigger() +{ + const Int_t entries = fEvt->GetNumPixels(); + + for (Int_t i=0; iIsPixelUsed(id1) ) + continue; + + const MGeomPix &gpix1 = (*fCam)[id1]; + if ( gpix1.IsInOutermostRing() ) + continue; + + const MCerPhotPix &pix1 = *fEvt->GetPixById(id1); + + const Double_t photons1 = pix1.GetNumPhotons(); + if (photons1 >= fMinPhotons) + return kTRUE; + } + } + return kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the Basic selections. +// +Int_t MFSelBasic::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + + *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; + *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: bad run " << endl; + + *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; + *fLog << (int)(fCut[2]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: Zenith angle < " << fThetaMin << endl; + + *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; + *fLog << (int)(fCut[3]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: Zenith angle > " << fThetaMax << endl; + + *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ; + *fLog << (int)(fCut[4]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: Software trigger not fullfilled" ; + *fLog << " (with fMinPhotons = " << fMinPhotons << ")" << endl; + + *fLog << " " << fCut[0] << " (" << (int)(fCut[0]*100/GetNumExecutions()) ; + *fLog << "%) Evts survived Basic selections!" << endl; + *fLog << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mfilter/MFSelBasic.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelBasic.h (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelBasic.h (revision 2663) @@ -0,0 +1,68 @@ +#ifndef MARS_MFSelBasic +#define MARS_MFSelBasic + +///////////////////////////////////////////////////////////////////////////// +// // +// MFSelBasic // +// // +// Class to evaluate basic cuts // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MMcEvt; +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MRawRunHeader; + +class MFSelBasic : public MFilter +{ +private: + const MMcEvt *fMcEvt; + const MGeomCam *fCam; // Camera Geometry + const MCerPhotEvt *fEvt; // Cerenkov Photon Event + //const MPedestalCam *fPed; // Pedestal information + const MRawRunHeader *fRawRun; + + Float_t fMinPhotons; + Float_t fThetaMin; + Float_t fThetaMax; + + Int_t fCut[5]; + + Bool_t fResult; + + Int_t Set(Int_t rc); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + Bool_t IsExpressionTrue() const { return fResult; } + + Bool_t SwTrigger(); + +public: + MFSelBasic(const char *name=NULL, const char *title=NULL); + + void SetCuts(Float_t minphotons, Float_t thetamin, Float_t thetamax); + + ClassDef(MFSelBasic, 0) // Class to evaluate basic cuts +}; + +#endif + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc (revision 2663) @@ -0,0 +1,193 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MFSelFinal +// +// WHAT ARE THE FINAL CUTS? +// +// This is a class to evaluate the Final Cuts +// (these cuts define the final sample of gammas; +// relevant for the calculation of the effective collection areas) +// +// to be called after the calculation of the hadroness +// +///////////////////////////////////////////////////////////////////////////// +#include "MFSelFinal.h" + +#include + +#include "MParList.h" + +#include "MMcEvt.hxx" + +#include "MCerPhotEvt.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MHadronness.h" + +ClassImp(MFSelFinal); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFSelFinal::MFSelFinal(const char *hilsrcname, + const char *name, const char *title) + : fHilSrcName(hilsrcname), fHadronnessName("MHadronness") +{ + fName = name ? name : "MFSelFinal"; + fTitle = title ? title : "Class to evaluate the Final Cuts"; + + // default values of cuts + SetCuts(1.0, 100.0, 10.); +} + +// -------------------------------------------------------------------------- +// +// Set cut values +// +// +void MFSelFinal::SetCuts(Float_t hadmax, Float_t alphamax, Float_t distmax) +{ + fHadronnessMax = hadmax; + fAlphaMax = alphamax; + fDistMax = distmax; + + *fLog << inf << "MFSelFinal cut values : fHadronnessMax, fAlphaMax, fDistMax = "; + *fLog << fHadronnessMax << ", " << fAlphaMax << ", " << fDistMax << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers +// +Int_t MFSelFinal::PreProcess(MParList *pList) +{ + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; + return kFALSE; + } + + fHadronness = (MHadronness*)pList->FindObject(fHadronnessName, "MHadronness"); + if (!fHadronness) + { + *fLog << dbginf << "MHadronness not found... aborting." << endl; + return kFALSE; + } + + MGeomCam *cam = (MGeomCam*)pList->FindCreateObj("MGeomCamMagic","MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + memset(fCut, 0, sizeof(fCut)); + + return kTRUE; +} + +Int_t MFSelFinal::Set(Int_t rc) +{ + fCut[rc]++; + fResult=kTRUE; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Evaluate final cuts +// +// if cuts are fulfilled set fResult = kTRUE +// +Int_t MFSelFinal::Process() +{ + const Double_t modalpha = fabs( fHilSrc->GetAlpha() ); + const Double_t h = fHadronness->GetHadronness(); + + fResult = kFALSE; + + if (h>fHadronnessMax) + return Set(1); + + if (modalpha>fAlphaMax) + return Set(2); + + if (fMm2Deg*fHilSrc->GetDist()>fDistMax) + return Set(3); + + fCut[0]++; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the Final selections. +// +Int_t MFSelFinal::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; + *fLog << (int)(fCut[1]*100/GetNumExecutions()); + *fLog << "%) Evts skipped due to: hadronness > "<< fHadronnessMax; + *fLog << " (hadronness from '" << fHadronnessName << "')" << endl; + + *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; + *fLog << (int)(fCut[2]*100/GetNumExecutions()); + *fLog << "%) Evts skipped due to: |ALPHA| > " << fAlphaMax; + *fLog << " [degrees]" << endl; + + *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; + *fLog << (int)(fCut[3]*100/GetNumExecutions()); + *fLog << "%) Evts skipped due to: DIST > " << fDistMax; + *fLog << " [degrees]" << endl; + + *fLog << " " << fCut[0] << " (" ; + *fLog << (int)(fCut[0]*100/GetNumExecutions()); + *fLog << "%) Evts survived Final selections!" << endl; + *fLog << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mfilter/MFSelFinal.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelFinal.h (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelFinal.h (revision 2663) @@ -0,0 +1,70 @@ +#ifndef MARS_MFSelFinal +#define MARS_MFSelFinal + +///////////////////////////////////////////////////////////////////////////// +// // +// MSelFinal // +// // +// Class to evaluate final cuts // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MHillas; +class MHillasSrc; +class MHadronness; + +class MFSelFinal : public MFilter +{ +private: + MHillasSrc *fHilSrc; + MHadronness *fHadronness; + + TString fHilName; + TString fHilSrcName; + TString fHadronnessName; + + Float_t fHadronnessMax; + Float_t fAlphaMax; + Float_t fDistMax; + + Double_t fMm2Deg; // conversion mm to degrees in camera + + Int_t fCut[4]; + + Bool_t fResult; + + Int_t Set(Int_t rc); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + Bool_t IsExpressionTrue() const { return fResult; } + +public: + MFSelFinal(const char *HilSrcName="MHillasSrc", + const char *name=NULL, const char *title=NULL); + + void SetHadronnessName(const TString name) { fHadronnessName = name; } + + void SetCuts(Float_t hadmax, Float_t alphamax, Float_t distmax); + + ClassDef(MFSelFinal, 0) // Class to evaluate final cuts +}; + +#endif + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc (revision 2663) @@ -0,0 +1,221 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 04/2003 +! Author(s): Thomas Bretz, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MFSelStandard +// +// This is a class to evaluate the Standard Cuts +// +// WHAT ARE THE STANDARD CUTS? // +// +// to be called after the calculation of the image parameters +// before the g/h separation +// +///////////////////////////////////////////////////////////////////////////// + +#include "MFSelStandard.h" + +#include "MParList.h" + +#include "MMcEvt.hxx" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MCerPhotEvt.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHillas.h" +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" + +ClassImp(MFSelStandard); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFSelStandard::MFSelStandard(const char *hilsrcname, + const char *name, const char *title) + : fHilName("MHillas"), fHilSrcName(hilsrcname), fImgParName("MNewImagePar") +{ + fName = name ? name : "MFSelStandard"; + fTitle = title ? title : "Class to evaluate the Standard Cuts"; + + // default values of cuts + SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); +} + +// -------------------------------------------------------------------------- +// +// Set the values for the cuts +// +// +void MFSelStandard::SetCuts(Float_t usedpixelsmax, Float_t corepixelsmin, + Float_t sizemin, Float_t distmin, Float_t distmax, + Float_t lengthmin, Float_t widthmin) +{ + fUsedPixelsMax = usedpixelsmax; + fCorePixelsMin = corepixelsmin; + fSizeMin = sizemin; + fDistMin = distmin; + fDistMax = distmax; + fLengthMin = lengthmin; + fWidthMin = widthmin; + + *fLog << inf << "MFSelStandard cut values : fUsedPixelsMax, fCorePixelsMin = "; + *fLog << fUsedPixelsMax << ", " << fCorePixelsMin << endl; + *fLog << inf << " fSizeMin, fDistMin, fDistMax = " << fSizeMin ; + *fLog << ", " << fDistMin << ", " << fDistMax << endl; + *fLog << inf << " fLengthMin, fWidthMin = " << fLengthMin ; + *fLog << ", " << fWidthMin << endl; +} + +// -------------------------------------------------------------------------- +// +// MISSING +// +Int_t MFSelStandard::PreProcess(MParList *pList) +{ + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); + if (!fHil) + { + *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewImgPar = (MNewImagePar*)pList->FindObject(fImgParName, "MNewImagePar"); + if (!fNewImgPar) + { + *fLog << err << fImgParName << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + memset(fCut, 0, sizeof(fCut)); + + return kTRUE; +} + +Bool_t MFSelStandard::Set(Int_t rc) +{ + fResult = kTRUE; + fCut[rc]++; + return kTRUE; +} +// -------------------------------------------------------------------------- +// +// Evaluate standard cuts +// +// if selections are fulfilled set fResult = kTRUE; +// +// +Int_t MFSelStandard::Process() +{ + const Double_t length = fHil->GetLength() * fMm2Deg; + const Double_t width = fHil->GetWidth() * fMm2Deg; + const Double_t dist = fHilSrc->GetDist()* fMm2Deg; + //const Double_t delta = fHil->GetDelta() * kRad2Deg; + const Double_t size = fHil->GetSize(); + const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels(); + const Int_t numcorepixels = fNewImgPar->GetNumCorePixels(); + + fResult = kFALSE; + + if (numusedpixels>=fUsedPixelsMax || numcorepixels<=fCorePixelsMin) + return Set(1); + + if (size<=fSizeMin ) + return Set(2); + + if (distfDistMax) + return Set(3); + + if (length<=fLengthMin || width<=fWidthMin) + return Set(4); + + fCut[0]++; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the Standard selections. +// +Int_t MFSelStandard::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fCut[1] << " (" << setw(3); + *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: Used pixels >= " << fUsedPixelsMax ; + *fLog << " or Core pixels <= " << fCorePixelsMin << endl; + + *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; + *fLog << (int)(fCut[2]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: SIZE <= " << fSizeMin << endl; + + *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; + *fLog << (int)(fCut[3]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: DIST < " << fDistMin; + *fLog << " or DIST > " << fDistMax << endl; + + *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ; + *fLog << (int)(fCut[4]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin; + *fLog << " or WIDTH <= " << fWidthMin << endl; + + *fLog << " " << fCut[0] << " (" ; + *fLog << (int)(fCut[0]*100/GetNumExecutions()) ; + *fLog << "%) Evts survived Standard selections!" << endl; + *fLog << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mfilter/MFSelStandard.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelStandard.h (revision 2663) +++ trunk/MagicSoft/Mars/mfilter/MFSelStandard.h (revision 2663) @@ -0,0 +1,78 @@ +#ifndef MARS_MFSelStandard +#define MARS_MFSelStandard + +///////////////////////////////////////////////////////////////////////////// +// // +// MFSelStandard // +// // +// Class to evaluate standard cuts // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MHillas; +class MHillasSrc; +class MNewImagePar; + +class MFSelStandard : public MFilter +{ +private: + MHillas *fHil; + MHillasSrc *fHilSrc; + MNewImagePar *fNewImgPar; + + TString fHilName; + TString fHilSrcName; + TString fImgParName; + + Float_t fUsedPixelsMax; + Float_t fCorePixelsMin; + Float_t fSizeMin; + Float_t fDistMin; + Float_t fDistMax; + Float_t fLengthMin; + Float_t fWidthMin; + + Double_t fMm2Deg; // conversion mm to degrees in camera + + Int_t fCut[5]; + + Bool_t fResult; + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + Bool_t IsExpressionTrue() const { return fResult; } + + Bool_t Set(Int_t rc); + +public: + MFSelStandard(const char *HilSrcName="MHillasSrc", + const char *name=NULL, const char *title=NULL); + + void SetHillasName(const char *name) { fHilName = name; } + void SetImgParName(const char *name) { fImgParName = name; } + + void SetCuts(Float_t usedpixelsmax, Float_t corepixelsmin, + Float_t sizemin, Float_t distmin, Float_t distmax, + Float_t lengthmin, Float_t widthmin); + + ClassDef(MFSelStandard, 0) // Class to evaluate standard cuts +}; + +#endif + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- trunk/MagicSoft/Mars/mfilter/Makefile (revision 2662) +++ trunk/MagicSoft/Mars/mfilter/Makefile (revision 2663) @@ -48,4 +48,7 @@ MFCT1SelStandard.cc \ MFCT1SelFinal.cc \ + MFSelBasic.cc \ + MFSelStandard.cc \ + MFSelFinal.cc \ MFEnergySlope.cc \ MFRandomSplit.cc Index: trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2662) +++ trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2663) @@ -44,4 +44,5 @@ #pragma link C++ class MHFindSignificance+; #pragma link C++ class MHCT1Supercuts+; +#pragma link C++ class MHSupercuts+; #pragma link C++ class MHCompProb+; Index: trunk/MagicSoft/Mars/mhist/MHSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSupercuts.cc (revision 2663) +++ trunk/MagicSoft/Mars/mhist/MHSupercuts.cc (revision 2663) @@ -0,0 +1,225 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek 2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////// +// +// MHSupercuts +// +// This class contains histograms for the supercuts +// +// the histograms are filled during the optimization of the supercuts +// +/////////////////////////////////////////////////////////////////////// +#include "MHSupercuts.h" + +#include + +#include +#include +#include +#include + +#include "MParList.h" +#include "MHFindSignificance.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHSupercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup four histograms for Alpha, and Dist +// +MHSupercuts::MHSupercuts(const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHSupercuts"; + fTitle = title ? title : "Container for histograms for the supercuts"; + + + fDegree = new TH1F("Degree", "Degree of polynomial", 5, -0.5, 4.5); + fProb = new TH1F("Prob", "chi2 probability", 40, 0, 1.0); + fNdf = new TH1F("NDF", "NDF of polynomial fit", 60, -0.5, 59.5); + fGamma = new TH1F("gamma", "gamma", 40, 0.0, 8.0); + fNexNon = new TH1F("NexNon", "Nex / Non", 50, 0.0, 1.0); + fSigLiMa= new TH1F("Significance", "significance of gamma signal", + 60, 0.0, 120.0); + fSigtoBackg= new TH2F("SigtoBackg", "Significance vs signal/backg ratio", + 50, 0.0, 10.0, 60, 0.0, 120.0); + fSigDegree = new TH2F("SigDegree", "Significance vs Degree of polynomial", + 5, -0.5, 4.5, 60, 0.0, 120.0); + fSigNbins = new TH2F("SigNbins", "Significance vs number of bins", + 40, -0.5, 79.5, 60, 0.0, 120.0); + + fDegree->SetDirectory(NULL); + fProb->SetDirectory(NULL); + fNdf->SetDirectory(NULL); + fGamma->SetDirectory(NULL); + fNexNon->SetDirectory(NULL); + fSigLiMa->SetDirectory(NULL); + fSigtoBackg->SetDirectory(NULL); + fSigDegree->SetDirectory(NULL); + fSigNbins->SetDirectory(NULL); + + fDegree->SetXTitle("order of polynomial"); + fProb->SetXTitle("chi2 probability of polynomial fit"); + fNdf->SetXTitle("NDF of polynomial fit"); + fGamma->SetXTitle("gamma"); + fNexNon->SetXTitle("Nex / Non"); + fSigLiMa->SetXTitle("significance"); + + fSigtoBackg->SetXTitle("signa./background ratio"); + fSigtoBackg->SetYTitle("significance"); + + fSigDegree->SetXTitle("order of polynomial"); + fSigDegree->SetYTitle("significance"); + + fSigNbins->SetXTitle("number of bins"); + fSigNbins->SetYTitle("significance"); + + fDegree->SetYTitle("Counts"); + fProb->SetYTitle("Counts"); + fNdf->SetYTitle("Counts"); + fGamma->SetYTitle("Counts"); + fNexNon->SetYTitle("Counts"); + fSigLiMa->SetYTitle("Counts"); + + fSigtoBackg->SetZTitle("Counts"); + fSigDegree->SetZTitle("Counts"); + fSigNbins->SetZTitle("Counts"); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHSupercuts::~MHSupercuts() +{ + delete fDegree; + delete fProb; + delete fNdf; + delete fGamma; + delete fNexNon; + delete fSigLiMa; + + delete fSigtoBackg; + delete fSigDegree; + delete fSigNbins; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms from the 'MHFindSignificance' container +// +Bool_t MHSupercuts::Fill(const MParContainer *par, const Stat_t w) +{ + if (!par) + { + *fLog << err << "MHSupercuts::Fill: Pointer (!=NULL) expected." << endl; + return kFALSE; + } + + MHFindSignificance &h = *(MHFindSignificance*)par; + + fDegree ->Fill(h.GetDegree( ), w); + fProb ->Fill(h.GetProb(), w); + fNdf ->Fill(h.GetNdf(), w); + fGamma ->Fill(h.GetGamma(), w); + + Double_t ratio = h.GetNon()>0.0 ? h.GetNex()/h.GetNon() : 0.0; + fNexNon ->Fill(ratio, w); + + fSigLiMa ->Fill(h.GetSigLiMa(), w); + + Double_t sigtobackg = h.GetNbg()!=0.0 ? h.GetNex() / h.GetNbg() : 0.0; + fSigtoBackg->Fill(sigtobackg, h.GetSigLiMa(), w); + + fSigDegree ->Fill(h.GetDegree(), h.GetSigLiMa(), w); + fSigNbins ->Fill(h.GetMbins(), h.GetSigLiMa(), w); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the two histograms into it. +// Be careful: The histograms belongs to this object and won't get deleted +// together with the canvas. +// +void MHSupercuts::Draw(Option_t *) +{ + //TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + //pad->SetBorderMode(0); + //AppendPad(""); + + TCanvas *pad = new TCanvas("Supercuts", "Supercut plots", 900, 900); + gROOT->SetSelectedPad(NULL); + + pad->Divide(3, 3); + + pad->cd(1); + gPad->SetBorderMode(0); + fDegree->Draw(); + + pad->cd(2); + gPad->SetBorderMode(0); + fProb->Draw(); + + pad->cd(3); + gPad->SetBorderMode(0); + fNdf->Draw(); + + pad->cd(4); + gPad->SetBorderMode(0); + fGamma->Draw(); + + pad->cd(5); + gPad->SetBorderMode(0); + fNexNon->Draw(); + + pad->cd(6); + gPad->SetBorderMode(0); + fSigLiMa->Draw(); + + pad->cd(7); + gPad->SetBorderMode(0); + fSigtoBackg->Draw(); + + pad->cd(8); + gPad->SetBorderMode(0); + fSigDegree->Draw(); + + pad->cd(9); + gPad->SetBorderMode(0); + fSigNbins->Draw(); + + pad->Modified(); + pad->Update(); +} Index: trunk/MagicSoft/Mars/mhist/MHSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSupercuts.h (revision 2663) +++ trunk/MagicSoft/Mars/mhist/MHSupercuts.h (revision 2663) @@ -0,0 +1,42 @@ +#ifndef MARS_MHSupercuts +#define MARS_MHSupercuts + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1F; +class TH2F; + +class MHSupercuts : public MH +{ +private: + TH1F *fDegree; // order of polynomial for background fit + TH1F *fProb; // chi2 probability of polynomial fit + TH1F *fNdf; // NDF of polynomial fit + TH1F *fGamma; // Nbg = gamma * Noff + TH1F *fNexNon; // no.of excess events / no.of events in signal region + TH1F *fSigLiMa; // significance of gamma signal + + TH2F *fSigtoBackg; // significance vs signal to background ratio (Nex/Nbg) + TH2F *fSigDegree; // significance vs order of polynomial + TH2F *fSigNbins; // significance vs number of bins + + +public: + MHSupercuts(const char *name=NULL, const char *title=NULL); + ~MHSupercuts(); + + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + void Draw(Option_t *opt=NULL); + + ClassDef(MHSupercuts, 1) // Container which holds histograms for the supercuts +}; + +#endif + + + + + Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 2662) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 2663) @@ -68,5 +68,6 @@ MHCalibrationPixel.cc \ MHCalibrationBlindPixel.cc \ - MHCalibrationPINDiode.cc + MHCalibrationPINDiode.cc \ + MHSupercuts.cc # MHCurrents.cc \