Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2151) +++ trunk/MagicSoft/Mars/Changelog (revision 2152) @@ -1,3 +1,23 @@ -*-*- END OF LINE -*-*- + + 2003/06/03: Wolfgang Wittek + + * mfileio/MCT1ReadPreProc.[h,cc] + - reset blind pixels for each event + (because they may have been changed by the padding) + + * macros/ONOFFCT1Analysis.C + - will be the macro for the CT1 analysis using ON and OFF data + + * manalysis/MPadONOFF.[h,cc] + - new class + - class for the padding of ON/OFF data + + * manalysis/MPadSchweizer.[h,cc] + - remove fBlinds->Clear() because the resetting of the + blind pixels is now done in MCT1ReadPreProc + + + 2003/06/02: Thomas Bretz Index: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2152) +++ trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C (revision 2152) @@ -0,0 +1,3191 @@ + +#include "CT1EgyEst.C" + +void InitBinnings(MParList *plist) +{ + gLog << "InitBinnings" << endl; + + MBinning *binse = new MBinning("BinningE"); + //binse->SetEdgesLog(30, 1.0e2, 1.0e5); + binse->SetEdges(30, 2, 5); + 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"); + Double_t yedge[9] = + {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; + TArrayD yed(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(9,zedge); + bincosth->SetEdges(zed); + plist->AddToList(bincosth); + + MBinning *binsdiff = new MBinning("BinningDiffsigma2"); + binsdiff->SetEdges(100, -5.0, 20.0); + plist->AddToList(binsdiff); +} + +void DeleteBinnings(MParList *plist) +{ + gLog << "DeleteBinnings" << endl; + + TObject *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; +} + +//************************************************************************ +void ONOFFCT1Analysis() +{ + gLog.SetNoColors(); + + if (gRandom) + delete gRandom; + gRandom = new TRandom3(0); + + //----------------------------------------------- + const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; + + //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc"; + const char *onfile = "~magican/ct1test/wittek/mkn421_00-01"; + + const char *mcfile = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc"; + //----------------------------------------------- + + // path for input for Mars + TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + + // path for output from Mars + TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + + //----------------------------------------------- + + 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 WPad = kFALSE; // write out padding histograms ? + Bool_t RPad = kTRUE; // read in padding histograms ? + Bool_t Wout = kTRUE; // write out root file ON1.root + // (or OFF1.root or MC1.root)? + + + // Job A_MC : read MC gamma data, + // - read sigmabar vs. Theta plot from ON data + // - do padding; + // - write root file for MC gammas (MC1.root); + + Bool_t JobA_MC = kFALSE; + Bool_t WMC1 = kFALSE; // write out root file MC1.root ? + + + // Job B_NN_UP : read ON1.root (or MC1.root) file + // - depending on RlookNN : create (or read in) hadron and gamma matrix + // - calculate hadroness for method of NEAREST NEIGHBORS + // and for the SUPERCUTS + // - update the input files with the hadronesses (ON1.root or MC1.root) + + Bool_t JobB_NN_UP = kFALSE; + Bool_t RLookNN = kFALSE; // read in look-alike events + Bool_t WNN = kFALSE; // update input root file ? + + + // Job B_RF_UP : read ON1.root (or MC1.root) file + // - depending on RLook : create (or read in) hadron and gamma matrix + // - depending on RTree : create (or read in) trees + // - calculate hadroness for method of RANDOM FOREST + // - update the input files with the hadroness (ON1.root or MC1.root) + + Bool_t JobB_RF_UP = kFALSE; + Bool_t RLookRF = kFALSE; // read in look-alike events + Bool_t RTree = kFALSE; // read in trees + Bool_t WRF = kFALSE; // update input root file ? + + + + // Job C: + // - read ON1.root and MC1.root files + // which should have been updated to contain the hadronnesses + // for the method of + // NEAREST NEIGHBORS + // SUPERCUTS and + // RF + // - produce Neyman-Pearson plots + + Bool_t JobC = kFALSE; + + + // Job D : + // - select g/h separation method XX + // - read ON2 (or MC2) root file + // - apply cuts in hadronness + // - make plots + + Bool_t JobD = kFALSE; + + + + // Job E_EST_UP : + // - read MC1.root file + // - select g/h separation method XX + // - optimize energy estimation for events passing the final cuts + // - write parameters of energy estimator onto file + // - update ON1.root and MC1.root files with estimated energy + // (ON_XX1.root and MC_XX1.root) + + Bool_t JobE_EST_UP = kFALSE; + Bool_t WESTUP = kFALSE; // update root files ? + + + + // Job F_XX : + // - select g/h separation method XX + // - read MC_XX2.root file +   // - calculate eff. collection area + // - read ON_XX2.root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts (ON3.root)) + + Bool_t JobF_XX = kFALSE; + Bool_t WXX = kFALSE; // write out root file ON3.root ? + + + //************************************************************************ + + + //--------------------------------------------------------------------- + // 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 CT1Analysis : Start of Job A" << endl; + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = " + << JobA << ", " << WPad << ", " << RPad << ", " << Wout << 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 (!RPad) + { + gLog << "=====================================================" << endl; + gLog << "Start generating the padding histograms" << endl; + //----------------------------------------- + // ON events + + gLog << "-----------" << endl; + gLog << "ON events :" << endl; + gLog << "-----------" << endl; + + MTaskList tliston; + MParList pliston; + + MCT1ReadPreProc readON(fileON); + + MBlindPixelCalc blindon; + blindon.SetUseBlindPixels(); + + MFCT1SelBasic selbasicon; + MContinue contbasicon(&selbasicon); + + MHBlindPixels blindON("BlindPixelsON"); + MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels"); + + MSigmabarCalc sigbarcalcon; + + MHSigmaTheta sigthON("SigmaThetaON"); + MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + pliston.AddToList(&blindON); + pliston.AddToList(&sigthON); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&readON); + + 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 = 1000; + 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; + + MCT1ReadPreProc readOFF(fileOFF); + + MBlindPixelCalc blindoff; + blindoff.SetUseBlindPixels(); + + MFCT1SelBasic selbasicoff; + MContinue contbasicoff(&selbasicoff); + + MHBlindPixels blindOFF("BlindPixelsOFF"); + MHBlindPixels *pblindOFF = &blindOFF; + gLog << "pblindOFF = " << pblindOFF << endl; + + MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels"); + + MSigmabarCalc sigbarcalcoff; + + MHSigmaTheta sigthOFF("SigmaThetaOFF"); + MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); + + //***************************** + // entries in MParList + + plistoff.AddToList(&tlistoff); + InitBinnings(&plistoff); + plistoff.AddToList(&blindOFF); + plistoff.AddToList(&sigthOFF); + + + //***************************** + // entries in MTaskList + + tlistoff.AddToList(&readOFF); + + 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 = 1000; + 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; + + gLog << "Merge the ON and OFF histograms and define the histograms for the padding :" << endl; + pad.MergeHistograms(sigthon, sigthoff, + sigpixthon, sigpixthoff, + diffpixthon, diffpixthoff, + blindidthon, blindidthoff, + blindnthon, blindnthoff); + gLog << "The histograms for the padding have been defined" << endl; + gLog << "=====================================================" << endl; + + + + 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("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MCT1ReadPreProc read(filenamein); + + if (typeInput == "ON") + { + MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", + "MCT1PointingCorrCalc"); + } + + MBlindPixelCalc blind; + blind.SetUseBlindPixels(); + + MFCT1SelBasic 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"); + // -------------------------------------------------- + + MFCT1SelStandard 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("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(&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 = 10000; + 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 CT1Analysis : End of Job A" << endl; + gLog << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job A_MC + //========= + + // read MC gamma data + + // - to pad them + // (using the 2D-histogram "sigmabar versus Theta" + // (SigmaTheta_ON.root) of the ON data) + + // - to write a file of padded MC gamma events (MC1.root) + // (after the standard cuts, before the g/h separation) + // (to be used together with the corresponding hadron file + // for the optimization of the g/h separation) + + + if (JobA_MC) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job A_MC" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobA_MC, WMC1 = " + << JobA_MC << ", " << WMC1 << endl; + + + // name of input root file + TString filenamein(mcfile); + + // name of output root file + TString outNameImage = outPath; + outNameImage += "MC"; + outNameImage += "1.root"; + + //------------------------------------------------ + // use for padding sigmabar vs. Theta from ON data + TString typeHist = "ON"; + gLog << "typeHist = " << typeHist << endl; + + // name of file containing the histograms for the padding + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta_"; + outNameSigTh += typeHist; + outNameSigTh += ".root"; + + + //------------------------------------ + // Get the histograms "2D-ThetaSigmabar" + // and "3D-ThetaPixSigma" + // and "3D-ThetaPixDiff" + // and "2D-IdBlindPixels" + // and "2D-NBlindPixels" + + + gLog << "Reading in file " << outNameSigTh << endl; + + TFile *infile = new TFile(outNameSigTh); + infile->ls(); + + TH2D *fHSigmaTheta = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + gLog << "Object '2D-ThetaSigmabar' not found on root file" << endl; + return; + } + gLog << "Object '2D-ThetaSigmabar' was read in" << endl; + + TH3D *fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + gLog << "Object '3D-ThetaPixSigma' not found on root file" << endl; + return; + } + gLog << "Object '3D-ThetaPixSigma' was read in" << endl; + + TH3D *fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + gLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; + return; + } + gLog << "Object '3D-ThetaPixDiff' was read in" << endl; + + + TH2D *fHIdBlindPixels = + (TH2D*) gROOT->FindObject("2D-IdBlindPixels"); + if (!fHIdBlindPixels) + { + gLog << "Object '2D-IdBlindPixels' not found on root file" << endl; + return; + } + gLog << "Object '2D-IdBlindPixels' was read in" << endl; + + TH2D *fHNBlindPixels = + (TH2D*) gROOT->FindObject("2D-NBlindPixels"); + if (!fHNBlindPixels) + { + gLog << "Object '2D-NBlindPixels' not found on root file" << endl; + return; + } + gLog << "Object '2D-NBlindPixels' was read in" << endl; + + //------------------------------------ + + MTaskList tlist; + MParList plist; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)plist->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MCT1ReadPreProc read(filenamein); + + MBlindPixelCalc blind; + blind.SetUseBlindPixels(); + + MFCT1SelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + + // There are 2 options for Thomas Schweizer's padding + // fPadFlag = 1 get Sigmabar from fHSigmaTheta + // and Sigma from fHDiffPixTheta + // fPadFlag = 2 get Sigma from fHSigmaPixTheta + + MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)"); + padthomas.SetHistograms(fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta, + fHIdBlindPixels, fHNBlindPixels); + padthomas.SetPadFlag(1); + + MFillH fillblind("MCBlindPixels[MHBlindPixels]", "MBlindPixels"); + fillblind.SetName("HBlind"); + + + //........................................... + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("MCSigmaTheta[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"); + // -------------------------------------------------- + + MFCT1SelStandard 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("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 + + plist.AddToList(&tlist); + InitBinnings(&plist); + + plist.AddToList(&source); + + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&padthomas); + tlist.AddToList(&blind); + + tlist.AddToList(&contbasic); + tlist.AddToList(&fillblind); + tlist.AddToList(&sigbarcalc); + tlist.AddToList(&fillsigtheta); + tlist.AddToList(&clean); + + tlist.AddToList(&hcalc); + tlist.AddToList(&hsrccalc); + + tlist.AddToList(&hfill1); + tlist.AddToList(&hfill2); + tlist.AddToList(&hfill3); + tlist.AddToList(&hfill4); + tlist.AddToList(&hfill5); + + tlist.AddToList(&contstandard); + if (WMC1) + tlist.AddToList(&write); + + //***************************** + + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.ReadEnv(env, "", printEnv); + evtloop.SetProgressBar(&bar); + //if (WMC1) + // evtloop.Write(); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tlist.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + plist.FindObject("MCSigmaTheta", "MHSigmaTheta")->DrawClone(); + plist.FindObject("MCBlindPixels", "MHBlindPixels")->DrawClone(); + + plist.FindObject("MHHillas")->DrawClone(); + plist.FindObject("MHHillasExt")->DrawClone(); + plist.FindObject("MHHillasSrc")->DrawClone(); + plist.FindObject("MHNewImagePar")->DrawClone(); + plist.FindObject("MHStarMap")->DrawClone(); + + + + DeleteBinnings(&plist); + + gLog << "Macro CT1Analysis : End of Job A_MC" + << endl; + gLog << "=========================================================" + << endl; + } + + + + //--------------------------------------------------------------------- + // Job B_NN_UP + //============ + + // - read in the matrices of look-alike events for gammas and hadrons + + // then read ON1.root (or MC1.root) file + // + // - calculate the hadroness for the method of nearest neighbors + // + // - update input root file, including the hadroness + + + if (JobB_NN_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobB_NN_UP, RLookNN, WNN = " + << JobB_NN_UP << ", " << RLookNN << ", " << WNN << endl; + + + + //-------------------------------------------- + // file to be updated (either ON or MC) + + TString typeInput = "ON"; + //TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + // name of input root file + TString filenameData = outPath; + filenameData += typeInput; + filenameData += "1.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "2.root"; + + //TString outNameImage = filenameData; + + gLog << "outNameImage = " << outNameImage << endl; + + //-------------------------------------------- + // files to be read for generating the look-alike events + // "hadrons" : + TString filenameON = outPath; + filenameON += "ON"; + filenameON += "1.root"; + Int_t howManyHadrons = 500000; + gLog << "filenameON = " << filenameON << ", howManyHadrons = " + << howManyHadrons << endl; + + + // "gammas" : + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += "1.root"; + Int_t howManyGammas = 50000; + gLog << "filenameMC = " << filenameMC << ", howManyGammas = " + << howManyGammas << endl; + + //-------------------------------------------- + // files of look-alike events + + TString outNameGammas = outPath; + outNameGammas += "matrix_gammas_"; + outNameGammas += "MC"; + outNameGammas += ".root"; + + TString typeMatrixHadrons = "ON"; + gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + TString outNameHadrons = outPath; + outNameHadrons += "matrix_hadrons_"; + outNameHadrons += typeMatrixHadrons; + outNameHadrons += ".root"; + + + MHMatrix matrixg("MatrixGammas"); + MHMatrix matrixh("MatrixHadrons"); + + //************************************************************************* + // read in matrices of look-alike events +if (RLookNN) + { + const char* mtxName = "MatrixGammas"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Get matrix for (gammas)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << outNameGammas << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outNameGammas' + // + TFile fileg(outNameGammas); + + matrixg.Read(mtxName); + matrixg.Print("cols"); + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Get matrix for (hadrons)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << outNameHadrons << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outNameHadrons' + // + TFile fileh(outNameHadrons); + + matrixh.Read(mtxName); + matrixh.Print("cols"); + } + + + //************************************************************************* + // create matrices of look-alike events +if (!RLookNN) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Create matrices of look-alike events" << endl; + gLog << " Gammas :" << endl; + + + MParList plistg; + MTaskList tlistg; + MFilterList flistg; + + MParList plisth; + MTaskList tlisth; + MFilterList flisth; + + MReadMarsFile readg("Events", filenameMC); + readg.DisableAutoScheme(); + + MReadMarsFile readh("Events", filenameON); + readh.DisableAutoScheme(); + + MFParticleId fgamma("MMcEvt", '=', kGAMMA); + fgamma.SetName("gammaID"); + MFParticleId fhadrons("MMcEvt", '!', kGAMMA); + fhadrons.SetName("hadronID)"); + + MFEventSelector selectorg; + selectorg.SetNumSelectEvts(howManyGammas); + selectorg.SetName("selectGammas"); + MFEventSelector selectorh; + selectorh.SetNumSelectEvts(howManyHadrons); + selectorh.SetName("selectHadrons"); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&flistg); + fillmatg.SetName("fillGammas"); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&flisth); + fillmath.SetName("fillHadrons"); + + + //***************************** fill gammas *** + // entries in MFilterList + + flistg.AddToList(&fgamma); + flistg.AddToList(&selectorg); + + //***************************** + // entries in MParList + + plistg.AddToList(&tlistg); + InitBinnings(&plistg); + + plistg.AddToList(&matrixg); + + //***************************** + // entries in MTaskList + + tlistg.AddToList(&readg); + tlistg.AddToList(&flistg); + tlistg.AddToList(&fillmatg); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtloopg; + evtloopg.SetParList(&plistg); + evtloopg.ReadEnv(env, "", printEnv); + evtloopg.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtloopg.Eventloop(maxevents)) + return; + + tlistg.PrintStatistics(0, kTRUE); + + + //***************************** fill hadrons *** + + gLog << " Hadrons :" << endl; + // entries in MFilterList + + flisth.AddToList(&fhadrons); + flisth.AddToList(&selectorh); + + //***************************** + // entries in MParList + + plisth.AddToList(&tlisth); + InitBinnings(&plisth); + + plisth.AddToList(&matrixh); + + //***************************** + // entries in MTaskList + + tlisth.AddToList(&readh); + tlisth.AddToList(&flisth); + tlisth.AddToList(&fillmath); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtlooph; + evtlooph.SetParList(&plisth); + evtlooph.ReadEnv(env, "", printEnv); + evtlooph.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtlooph.Eventloop(maxevents)) + return; + + tlisth.PrintStatistics(0, kTRUE); + //***************************************************** + + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + // + // ---------------------------------------------------------- + // Definition of the reference sample of look-alike events + // (this is a subsample of the original sample) + // ---------------------------------------------------------- + // + gLog << "" << endl; + gLog << "========================================================" << endl; + // Select a maximum of nmaxevts events from the sample of look-alike + // events. They will form the reference sample. + Int_t nmaxevts = 2000; + + // target distribution for the variable in column refcolumn (start at 0); + //Int_t refcolumn = 0; + //Int_t nbins = 5; + //Float_t frombin = 0.5; + //Float_t tobin = 1.0; + //TH1F *thsh = new TH1F("thsh","target distribution", + // nbins, frombin, tobin); + //thsh->SetDirectory(NULL); + //thsh->SetXTitle("cos( \\Theta)"); + //thsh->SetYTitle("Counts"); + //Float_t dbin = (tobin-frombin)/nbins; + //Float_t lbin = frombin +dbin*0.5; + //for (Int_t j=1; j<=nbins; j++) + //{ + // thsh->Fill(lbin,1.0); + // lbin += dbin; + //} + + Int_t refcolumn = 0; + MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); + TH1F *thsh = new TH1F(); + thsh->SetNameTitle("thsh","target distribution"); + MH::SetBinning(thsh, binscostheta); + Int_t nbins = thsh->GetNbinsX(); + Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; + for (Int_t j=1; j<=nbins; j++) + { + thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); + } + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; + gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrixg.EnableGraphicalOutput(); + Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); + + if ( !matrixok ) + { + gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; + return; + } + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; + gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrixh.EnableGraphicalOutput(); + matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); + delete thsh; + + if ( !matrixok ) + { + gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; + return; + } + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + // write out look-alike events + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Write out look=alike events" << endl; + + + //------------------------------------------- + // "gammas" + gLog << "Gammas :" << endl; + matrixg.Print("cols"); + + TFile writeg(outNameGammas, "RECREATE", ""); + matrixg.Write(); + + gLog << "" << endl; + gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " + << outNameGammas << endl; + + //------------------------------------------- + // "hadrons" + gLog << "Hadrons :" << endl; + matrixh.Print("cols"); + + TFile writeh(outNameHadrons, "RECREATE", ""); + matrixh.Write(); + + gLog << "" << endl; + gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " + << outNameHadrons << endl; + + } + //********** end of creating matrices of look-alike events *********** + + + //----------------------------------------------------------------- + // Update the input files with the NN and SC hadronness + // + + if (WNN) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update input file '" << filenameData + << "' with the NN and SC hadronness" << endl; + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //....................................................................... + // calculation of hadroness for method of Nearest Neighbors + + TString hadNNName = "HadNN"; + MMultiDimDistCalc nncalc; + nncalc.SetUseNumRows(25); + nncalc.SetUseKernelMethod(kFALSE); + nncalc.SetHadronnessName(hadNNName); + + //....................................................................... + // calculation of hadroness for the supercuts + // (=0.25 if fullfilled, =0.75 otherwise) + + TString hadSCName = "HadSC"; + MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc); + sccalc.SetHadronnessName(hadSCName); + + //....................................................................... + + //MWriteRootFile write(outNameImage, "UPDATE"); + MWriteRootFile write(outNameImage, "RECREATE"); + + 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"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadNNName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadNNName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&matrixg); + pliston.AddToList(&matrixh); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&nncalc); + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&write); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + //------------------------------------------- + // Display the histograms + // + pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + } + + + + gLog << "Macro CT1Analysis : End of Job B_NN_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job B_RF_UP + //============ + + + // - create (or read in) the matrices of look-alike events for gammas + // and hadrons + // - create (or read in) the trees + // - then read ON1.root (or MC1.root) file + // - calculate the hadroness for the method of RANDOM FOREST + // - update input root file with the hadroness + + + if (JobB_RF_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " + << JobB_RF_UP << ", " << RLookRF << ", " << RTree << ", " + << WRF << endl; + + + + //-------------------------------------------- + // parameters for the random forest + Int_t NumTrees = 100; + Int_t NumTry = 3; + Int_t NdSize = 3; + + + //-------------------------------------------- + // file to be updated (either ON or MC) + + TString typeInput = "ON"; + //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 generating the look-alike events + // "hadrons" : + TString filenameON = outPath; + filenameON += "ON"; + filenameON += "1.root"; + Int_t howManyHadrons = 1000000; + gLog << "filenameON = " << filenameON << ", howManyHadrons = " + << howManyHadrons << endl; + + + // "gammas" : + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += "1.root"; + Int_t howManyGammas = 50000; + gLog << "filenameMC = " << filenameMC << ", howManyGammas = " + << howManyGammas << endl; + + //-------------------------------------------- + // files of look-alike events + + TString outNameGammas = outPath; + outNameGammas += "RFmatrix_gammas_"; + outNameGammas += "MC"; + outNameGammas += ".root"; + + TString typeMatrixHadrons = "ON"; + gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + TString outNameHadrons = outPath; + outNameHadrons += "RFmatrix_hadrons_"; + outNameHadrons += typeMatrixHadrons; + outNameHadrons += ".root"; + + + MHMatrix matrixg("MatrixGammas"); + MHMatrix matrixh("MatrixHadrons"); + + //-------------------------------------------- + // file of trees of the random forest + + TString outRF = outPath; + outRF += "RF.root"; + + + //************************************************************************* + // read in matrices of look-alike events +if (RLookRF) + { + const char* mtxName = "MatrixGammas"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Get matrix for (gammas)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << outNameGammas << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outNameGammas' + // + TFile fileg(outNameGammas); + + matrixg.Read(mtxName); + matrixg.Print("cols"); + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Get matrix for (hadrons)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << outNameHadrons << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outNameHadrons' + // + TFile fileh(outNameHadrons); + + matrixh.Read(mtxName); + matrixh.Print("cols"); + } + + + //************************************************************************* + // create matrices of look-alike events +if (!RLookRF) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Create matrices of look-alike events" << endl; + gLog << " Gammas :" << endl; + + + MParList plistg; + MTaskList tlistg; + MFilterList flistg; + + MParList plisth; + MTaskList tlisth; + MFilterList flisth; + + MReadMarsFile readg("Events", filenameMC); + readg.DisableAutoScheme(); + + MReadMarsFile readh("Events", filenameON); + readh.DisableAutoScheme(); + + MFParticleId fgamma("MMcEvt", '=', kGAMMA); + fgamma.SetName("gammaID"); + MFParticleId fhadrons("MMcEvt", '!', kGAMMA); + fhadrons.SetName("hadronID)"); + + MFEventSelector selectorg; + selectorg.SetNumSelectEvts(howManyGammas); + selectorg.SetName("selectGammas"); + MFEventSelector selectorh; + selectorh.SetNumSelectEvts(howManyHadrons); + selectorh.SetName("selectHadrons"); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&flistg); + fillmatg.SetName("fillGammas"); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&flisth); + fillmath.SetName("fillHadrons"); + + + //***************************** fill gammas *** + // entries in MFilterList + + flistg.AddToList(&fgamma); + flistg.AddToList(&selectorg); + + //***************************** + // entries in MParList + + plistg.AddToList(&tlistg); + InitBinnings(&plistg); + + plistg.AddToList(&matrixg); + + //***************************** + // entries in MTaskList + + tlistg.AddToList(&readg); + tlistg.AddToList(&flistg); + tlistg.AddToList(&fillmatg); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtloopg; + evtloopg.SetParList(&plistg); + evtloopg.ReadEnv(env, "", printEnv); + evtloopg.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtloopg.Eventloop(maxevents)) + return; + + tlistg.PrintStatistics(0, kTRUE); + + + //***************************** fill hadrons *** + + gLog << " Hadrons :" << endl; + // entries in MFilterList + + flisth.AddToList(&fhadrons); + flisth.AddToList(&selectorh); + + //***************************** + // entries in MParList + + plisth.AddToList(&tlisth); + InitBinnings(&plisth); + + plisth.AddToList(&matrixh); + + //***************************** + // entries in MTaskList + + tlisth.AddToList(&readh); + tlisth.AddToList(&flisth); + tlisth.AddToList(&fillmath); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtlooph; + evtlooph.SetParList(&plisth); + evtlooph.ReadEnv(env, "", printEnv); + evtlooph.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtlooph.Eventloop(maxevents)) + return; + + tlisth.PrintStatistics(0, kTRUE); + //***************************************************** + + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + // + // ---------------------------------------------------------- + // Definition of the reference sample of look-alike events + // (this is a subsample of the original sample) + // ---------------------------------------------------------- + // + gLog << "" << endl; + gLog << "========================================================" << endl; + // Select a maximum of nmaxevts events from the sample of look-alike + // events. They will form the reference sample. + Int_t nmaxevts = 10000; + + // target distribution for the variable in column refcolumn (start at 0); + //Int_t refcolumn = 0; + //Int_t nbins = 5; + //Float_t frombin = 0.5; + //Float_t tobin = 1.0; + //TH1F *thsh = new TH1F("thsh","target distribution", + // nbins, frombin, tobin); + //thsh->SetDirectory(NULL); + //thsh->SetXTitle("cos( \\Theta)"); + //thsh->SetYTitle("Counts"); + //Float_t dbin = (tobin-frombin)/nbins; + //Float_t lbin = frombin +dbin*0.5; + //for (Int_t j=1; j<=nbins; j++) + //{ + // thsh->Fill(lbin,1.0); + // lbin += dbin; + //} + + Int_t refcolumn = 0; + MBinning *binscostheta = (MBinning*)plistg->FindObject("BinningCosTheta", "MBinning"); + TH1F *thsh = new TH1F(); + thsh->SetNameTitle("thsh","target distribution"); + MH::SetBinning(thsh, binscostheta); + Int_t nbins = thsh->GetNbinsX(); + Double_t cont[8] = {1500, 1500, 1500, 3250, 3250, 3900, 3900, 3900}; + for (Int_t j=1; j<=nbins; j++) + { + thsh->Fill(thsh->GetBinCenter(j), cont[j-1]); + } + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : define reference sample for gammas" << endl; + gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrixg.EnableGraphicalOutput(); + Bool_t matrixok = matrixg.DefRefMatrix(refcolumn, *thsh, nmaxevts); + + if ( !matrixok ) + { + gLog << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; + return; + } + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : define reference sample for hadrons" << endl; + gLog << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrixh.EnableGraphicalOutput(); + matrixok = matrixh.DefRefMatrix(refcolumn, *thsh, nmaxevts); + delete thsh; + + if ( !matrixok ) + { + gLog << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; + return; + } + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + // write out look-alike events + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Write out look=alike events" << endl; + + + //------------------------------------------- + // "gammas" + gLog << "Gammas :" << endl; + matrixg.Print("cols"); + + TFile writeg(outNameGammas, "RECREATE", ""); + matrixg.Write(); + + gLog << "" << endl; + gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " + << outNameGammas << endl; + + //------------------------------------------- + // "hadrons" + gLog << "Hadrons :" << endl; + matrixh.Print("cols"); + + TFile writeh(outNameHadrons, "RECREATE", ""); + matrixh.Write(); + + gLog << "" << endl; + gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " + << outNameHadrons << endl; + + } + //********** end of creating matrices of look-alike events *********** + + + MRanForest *fRanForest; + MRanTree *fRanTree; + //----------------------------------------------------------------- + // read in the trees of the random forest + if (RTree) + { + MParList plisttr; + MTaskList tlisttr; + plisttr.AddToList(&tlisttr); + + MReadTree readtr("TREE", outRF); + readtr.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(NumTrees); + + // list of tasks for the loop over the trees + + tlisttr.AddToList(&readtr); + tlisttr.AddToList(&rffill); + + //------------------- + // Execute tree loop + // + MEvtLoop evtlooptr; + evtlooptr.SetParList(&plisttr); + if (!evtlooptr.Eventloop()) + return; + + tlisttr.PrintStatistics(0, kTRUE); + + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + + //----------------------------------------------------------------- + // grow the trees of the random forest (event loop = tree loop) + + if (!RTree) + { + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro CT1Analysis : start growing trees" << endl; + + MTaskList tlist2; + MParList plist2; + plist2.AddToList(&tlist2); + + plist2.AddToList(&matrixg); + plist2.AddToList(&matrixh); + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(NumTrees); + rfgrow2.SetNumTry(NumTry); + rfgrow2.SetNdSize(NdSize); + + MWriteRootFile rfwrite2(outRF); + rfwrite2.AddContainer("MRanTree", "TREE"); + + // list of tasks for the loop over the trees + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + + //------------------- + // Execute tree loop + // + MEvtLoop treeloop; + treeloop.SetParList(&plist2); + + if ( !treeloop.Eventloop() ) + return; + + tlist2.PrintStatistics(0, kTRUE); + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + // end of growing the trees of the random forest + //----------------------------------------------------------------- + + + + + //----------------------------------------------------------------- + // Update the input files with the RF hadronness + // + + if (WRF) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update input file '" << filenameData + << "' with the RF hadronness" << endl; + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //....................................................................... + // calculate hadronnes for method of RANDOM FOREST + + TString hadRFName = "HadRF"; + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadRFName); + + //....................................................................... + + //MWriteRootFile write(outNameImage, "UPDATE"); + 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("HadNN", "Events"); + write.AddContainer("HadSC", "Events"); + + write.AddContainer(hadRFName, "Events"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadRFName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillranfor("MHRanForest"); + fillranfor.SetName("HRanForest"); + + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadRFName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(fRanForest); + pliston.AddToList(fRanTree); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillranfor); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&write); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + pliston.FindObject("MHRanForest")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + } + + gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + //--------------------------------------------------------------------- + // Job C + //====== + + // - read ON1 and MC1 data files + // which should have been updated to contain the hadronnesses + // for the method of NEAREST NEIGHBORS and for the SUOERCUTS + // - produce Neyman-Pearson plots + + if (JobC) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job C" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobC = " << JobC << endl; + + + // name of input data file + TString filenameData = outPath; + filenameData += "ON"; + filenameData += "3.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of input MC file + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += "3.root"; + gLog << "filenameMC = " << filenameMC << endl; + + + //----------------------------------------------------------------- + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + + //....................................................................... + // names of hadronness containers + + TString hadNNName = "HadNN"; + TString hadSCName = "HadSC"; + TString hadRFName = "HadRF"; + + //....................................................................... + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadNNName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadNNName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "Macro CT1Analysis : End of Job C" << endl; + gLog << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job D + //====== + + // - select g/h separation method XX + // - read ON2 (or MC2) root file + // - apply cuts in hadronness + // - make plots + + + if (JobD) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job D" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobD = " + << JobD << endl; + + // type of data to be analysed + //TString typeData = "ON"; + //TString typeData = "OFF"; + TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString ext = "3.root"; + + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("NN"); + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.20; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); + fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&contfinal); + + tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "Macro CT1Analysis : End of Job D" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job E_EST_UP + //============ + + // - read MC2.root file + // - select g/h separation method XX + // - optimize energy estimator for events passing final cuts + // - write parameters of energy estimator onto file "energyest_XX.root" + // + // - read ON2.root and MC2.root files + // - update input root file with the estimated energy + // (ON_XX2.root, MC_XX2.root) + + + if (JobE_EST_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = " + << JobE_EST_UP << ", " << WESTUP << endl; + + + TString typeON = "ON"; + TString typeMC = "MC"; + TString ext = "3.root"; + TString extout = "4.root"; + + //------------------------------ + // name of MC file to be used for optimizing the energy estimator + TString filenameOpt(outPath); + filenameOpt += typeMC; + filenameOpt += ext; + gLog << "filenameOpt = " << filenameOpt << endl; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + TString XX("NN"); + //TString XX("SC"); + //TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |alpha| and dist + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + // name of file containing the parameters of the energy estimator + TString energyParName(outPath); + energyParName += "energyest_"; + energyParName += XX; + energyParName += ".root"; + gLog << "energyParName = " << energyParName << endl; + + + //------------------------------ + // name of ON file to be updated + TString filenameON(outPath); + filenameON += typeON; + filenameON += ext; + gLog << "filenameON = " << filenameON << endl; + + // name of MC file to be updated + TString filenameMC(outPath); + filenameMC += typeMC; + filenameMC += ext; + gLog << "filenameMC = " << filenameMC << endl; + + //------------------------------ + // name of updated ON file + TString filenameONup(outPath); + filenameONup += typeON; + filenameONup += "_"; + filenameONup += XX; + filenameONup += extout; + gLog << "filenameONup = " << filenameONup << endl; + + // name of updated MC file + TString filenameMCup(outPath); + filenameMCup += typeMC; + filenameMCup += "_"; + filenameMCup += XX; + filenameMCup += extout; + gLog << "filenameMCup = " << filenameMCup << endl; + + //----------------------------------------------------------- + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //=========================================================== + // + // Optimization of energy estimator + // + + TString inpath(""); + TString outpath(""); + Int_t howMany = 2000; + CT1EEst(inpath, filenameOpt, outpath, energyParName, + fHilName, fHilNameSrc, fhadronnessName, + howMany, maxhadronness, maxalpha, maxdist); + + //----------------------------------------------------------- + // + // Read in parameters of energy estimator + // + gLog << "================================================================" + << endl; + gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '" + << energyParName << "'" << endl; + TFile enparam(energyParName); + MMcEnergyEst mcest("MMcEnergyEst"); + mcest.Read("MMcEnergyEst"); + enparam.Close(); + + TArrayD parA(5); + TArrayD parB(7); + for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameON); + read.DisableAutoScheme(); + + //--------------------------- + // calculate estimated energy + + MEnergyEstParam eest2(fHilName); + eest2.Add(fHilNameSrc); + + eest2.SetCoeffA(parA); + eest2.SetCoeffB(parB); + + //....................................................................... + + MWriteRootFile write(filenameONup); + + 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"); + + //----------------------------------------------------------------- + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + MContinue contfinal(&selfinal); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + tliston.AddToList(&eest2); + tliston.AddToList(&write); + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //--------------------------------------------------- + //--------------------------------------------------- + // Update MC data + // + gLog << "============================================================" + << endl; + gLog << "Macro CT1Analysis.C : update file '" << filenameMC + << "'" << endl; + + MTaskList tlistmc; + MParList plistmc; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameMC); + read.DisableAutoScheme(); + + //--------------------------- + // calculate estimated energy + + MEnergyEstParam eest2(fHilName); + eest2.Add(fHilNameSrc); + + eest2.SetCoeffA(parA); + eest2.SetCoeffB(parB); + + //....................................................................... + + MWriteRootFile write(filenameMCup); + + 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"); + + //----------------------------------------------------------------- + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + MContinue contfinal(&selfinal); + + + //***************************** + // entries in MParList + + plistmc.AddToList(&tlistmc); + InitBinnings(&plistmc); + + + //***************************** + // entries in MTaskList + + tlistmc.AddToList(&read); + tlistmc.AddToList(&eest2); + tlistmc.AddToList(&write); + tlistmc.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plistmc); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tlistmc.PrintStatistics(0, kTRUE); + DeleteBinnings(&plistmc); + + + //========== end update ============================================ + } + + enparam.Close(); + + gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job F_XX + //========= + + // - select g/h separation method XX + // - read MC_XX2.root file +   // - calculate eff. collection area + // - read ON_XX2.root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts (ON_XX3.root)) + + + if (JobF_XX) + { + gLog << "=====================================================" << endl; + gLog << "Macro CT1Analysis : Start of Job F_XX" << endl; + + gLog << "" << endl; + gLog << "Macro CT1Analysis : JobF_XX, WXX = " + << JobF_XX << ", " << WXX << endl; + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString typeMC = "MC"; + TString ext = "2.root"; + TString extout = "3.root"; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("NN"); + TString XX("SC"); + //TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of MC file to be used for calculating the eff. collection areas + TString filenameArea(outPath); + filenameArea += typeMC; + filenameArea += "_"; + filenameArea += XX; + filenameArea += ext; + gLog << "filenameArea = " << filenameArea << endl; + + //------------------------------ + // name of file containing the eff. collection areas + TString collareaName(outPath); + collareaName += "area_"; + collareaName += XX; + collareaName += ".root"; + gLog << "collareaName = " << collareaName << endl; + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += "_"; + filenameData += XX; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + //------------------------------ + // name of output data file (after the final cuts) + TString filenameDataout(outPath); + filenameDataout += typeData; + filenameDataout += "_"; + filenameDataout += XX; + filenameDataout += extout; + gLog << "filenameDataout = " << filenameDataout << endl; + + + //==================================================================== + gLog << "-----------------------------------------------" << endl; + gLog << "Start calculation of effective collection areas" << endl; + MParList parlist; + MTaskList tasklist; + + //--------------------------------------- + // Setup the tasks to be executed + // + MReadMarsFile reader("Events", filenameArea); + reader.DisableAutoScheme(); + + MFCT1SelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); + + MContinue conthadrons(&cuthadrons); + + MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea(); + + MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); + filler.SetName("CollectionArea"); + + //******************************** + // entries in MParList + + parlist.AddToList(&tasklist); + InitBinnings(&parlist); + parlist.AddToList(collarea); + + //******************************** + // entries in MTaskList + + tasklist.AddToList(&reader); + tasklist.AddToList(&conthadrons); + tasklist.AddToList(&filler); + + //******************************** + + //----------------------------------------- + // Execute event loop + // + MEvtLoop magic; + magic.SetParList(&parlist); + + MProgressBar bar; + magic.SetProgressBar(&bar); + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(0, kTRUE); + + // Calculate effective collection areas + // and display the histograms + // + //MHMcCT1CollectionArea *collarea = parlist.FindObject("MHMcCT1CollectionArea"); + collarea->CalcEfficiency(); + collarea->DrawClone("lego"); + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + f.Close(); + + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + //....................................................................... + + + MWriteRootFile write(filenameDataout); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + write.AddContainer("HadNN", "Events"); + write.AddContainer("HadSC", "Events"); + write.AddContainer("HadRF", "Events"); + + write.AddContainer("MEnergyEst", "Events"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + MFCT1SelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); + fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + MFCT1SelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + if (WXX) + 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); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "Macro CT1Analysis : End of Job F_XX" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + +} + + + Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2151) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2152) @@ -41,4 +41,5 @@ #pragma link C++ class MPadding+; #pragma link C++ class MPadSchweizer+; +#pragma link C++ class MPadONOFF+; #pragma link C++ class MCT1PointingCorrCalc+; @@ -54,2 +55,3 @@ #endif + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2152) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2152) @@ -0,0 +1,1368 @@ +/* ======================================================================== *\ +! +! * +! * 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. +// +// 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. +// +// This implementation has been tested for CT1 data. For MAGIC some +// modifications are necessary. +// +///////////////////////////////////////////////////////////////////////////// +#include "MPadONOFF.h" + +#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); + + +// -------------------------------------------------------------------------- +// +// 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 << "MPadONOFF: fPadFlag = " << fPadFlag << endl; + + fType = ""; + + fHSigmaTheta = NULL; + fHSigmaPixTheta = NULL; + fHDiffPixTheta = NULL; + + fHgON = NULL; + fHgOFF = NULL; + + fHBlindPixIdTheta = NULL; + fHBlindPixNTheta = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPadONOFF::~MPadONOFF() +{ + if (fHSigmaTheta != NULL) delete fHSigmaTheta; + + 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) +{ + //------------------------------------------------------------- + // merge the distributions of ON and OFF events + // to obtain the target distribution fHSigmaTheta + // + + fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + Int_t nbinstheta = sigthon->GetNbinsX(); + Int_t n = sigthon->GetNbinsY(); + Double_t eps = 1.e-10; + + fHgON = new TH3D; + fHgON->SetBins(nbinstheta, 0.5, 0.5+(Double_t)nbinstheta, + n, 0.5, 0.5+(Double_t)n, + n, 0.5, 0.5+(Double_t)n); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + fHgOFF->SetBins(nbinstheta, 0.5, 0.5+(Double_t)nbinstheta, + n, 0.5, 0.5+(Double_t)n, + n, 0.5, 0.5+(Double_t)n); + 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 + + 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 + + // fHg[k][j] (if <0.0) tells how many ON events in bin k should be padded + // from sigma_k to sigma_j + + + // fHg[k][j] (if >0.0) tells how many OFF events in bin k 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=n; 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 << "k, j, arg = " << 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 << "k, j, arg = " << 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<=n; 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 << "MPadONOFF::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<=n; 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 + // SigmaPixTheta and DiffPixTheta + // BlindPixIdTheta and BlindPixNTheta + // are taken from the ON data + + fHSigmaPixTheta = sigpixthon; + fHDiffPixTheta = diffpixthon; + + fHBlindPixIdTheta = blindidthon; + fHBlindPixNTheta = blindnthon; + + //-------------------------------------------- + + fHSigmaTheta->SetDirectory(NULL); + fHSigmaPixTheta->SetDirectory(NULL); + fHDiffPixTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixNTheta->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 << "Read padding histograms from file " << namefilein << endl; + + fInfile = new TFile(namefilein); + fInfile->ls(); + + //------------------------------------ + + fHSigmaTheta = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '2D-ThetaSigmabar' was read in" << endl; + + *fLog << "ReadTargetDist : fHSigmaTheta = " << fHSigmaTheta << endl; + + fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixSigma' was read in" << endl; + + fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixDiff' was read in" << endl; + + fHgON = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); + if (!fHgON) + { + *fLog << "Object '3D-PaddingMatrixON' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '3D-PaddingMatrixON' was read in" << endl; + + fHgOFF = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); + if (!fHgOFF) + { + *fLog << "Object '3D-PaddingMatrixOFF' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '3D-PaddingMatrixOFF' was read in" << endl; + + + fHBlindPixIdTheta = + (TH2D*) gROOT->FindObject("2D-IdBlindPixels"); + if (!fHBlindPixIdTheta) + { + *fLog << "Object '2D-IdBlindPixels' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '2D-IdBlindPixels' was read in" << endl; + + fHBlindPixNTheta = + (TH2D*) gROOT->FindObject("2D-NBlindPixels"); + if (!fHBlindPixNTheta) + { + *fLog << "Object '2D-NBlindPixels' not found on root file" << endl; + return kFALSE; + } + *fLog << "Object '2D-NBlindPixels' was read in" << endl; + + *fLog << "ReadTargetDist : fHBlindPixNTheta = " << fHBlindPixNTheta << endl; + + //------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Write target distributions onto a file +// +// +Bool_t MPadONOFF::WriteTargetDist(const char* namefileout) +{ + *fLog << "Write padding histograms onto file " << namefileout << endl; + + TFile outfile(namefileout, "RECREATE"); + + fHBlindPixNTheta->Write(); + fHBlindPixIdTheta->Write(); + + fHSigmaTheta->Write(); + fHSigmaPixTheta->Write(); + fHDiffPixTheta->Write(); + + fHgON->Write(); + fHgOFF->Write(); + + *fLog << "MPadONOFF::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 << "MPadONOFF::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. +// +// 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 << "MPadONOFF::SetPadFlag(); choose option " << fPadFlag << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers and prepare the histograms +// +Bool_t MPadONOFF::PreProcess(MParList *pList) +{ + if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta || + !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"); + + + //-------------------------------------------------------------------- + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the Padding +// idealy the events to be padded should have been generated without NSB +// +Bool_t MPadONOFF::Process() +{ + //*fLog << "Entry MPadONOFF::Process();" << endl; + + //------------------------------------------------ + // 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].GetMeanRms()); + } + + + + //----------------------------------------- + Int_t rc=0; + + const UInt_t npix = fEvt->GetNumPixels(); + + //fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << "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 = fSigmabar->Calc(*fCam, *fPed, *fEvt); + Double_t sigbarold2 = sigbarold*sigbarold; + //fSigmabar->Print(""); + + // for MC data : expect sigmabar to be zero before the padding + if (fType == "MC") + { + if (sigbarold > 0) + { + //*fLog << "MPadONOFF::Process(); sigmabar of event to be padded is > 0 : " + // << sigbarold << ". 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 << "theta = " << theta << endl; + + Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) + { + //*fLog << "MPadONOFF::Process(); binNumber out of range : theta, binTheta = " + // << theta << ", " << binTheta << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + //------------------------------------------- + // for the current theta, + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // + + + + // numBlind is the number of blind pixels in this event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, ""); + if ( nblind->GetEntries() == 0.0 ) + { + *fLog << "MPadONOFF::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; + + + // 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, ""); + //fBlinds->Clear(); + if ( hblind->GetEntries() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + //*fLog << "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); + + Double_t prob; + TH1D *hpad = NULL; + if (fType == "ON") + { + hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); + prob = hpad->Integral(); + } + else if (fType == "OFF") + { + hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, ""); + prob = hpad->Integral(); + } + else + prob = 1.0; + + if ( fType != "MC" && + (prob <= 0.0 || gRandom->Uniform() > prob) ) + { + delete hpad; + // event should not be padded + rc = 8; + fErrors[rc]++; + return kTRUE; + } + // event should be padded + + //------------------------------------------- + // for the current theta, generate a sigmabar + // + // for ON/OFF data according to the matrix fHg + // for MC data according to the histogram fHSigmaTheta + // + Double_t sigmabar=0; + if (fType == "ON" || fType == "OFF") + { + sigmabar = hpad->GetRandom(); + 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->GetEntries() == 0.0 ) + { + *fLog << "MPadONOFF::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 = hsigma->GetRandom(); + //*fLog << "Theta, bin number = " << theta << ", " << binTheta + // << ", sigmabar = " << sigmabar << endl + } + delete hsigma; + } + + const Double_t sigmabar2 = sigmabar*sigmabar; + + //------------------------------------------- + + //*fLog << "MPadONOFF::Process(); sigbarold, sigmabar = " + // << sigbarold << ", "<< sigmabar << endl; + + // Skip event if target sigmabar is <= sigbarold + if (sigmabar <= sigbarold) + { + *fLog << "MPadONOFF::Process(); target sigmabar is less than sigbarold : " + << sigmabar << ", " << sigbarold << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; + return kCONTINUE; + } + + + //------------------------------------------- + // + // Calculate average number of NSB photons 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 = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + 0, 9999, ""); + Double_t RMS = hnoise->GetRMS(1); + delete hnoise; + + elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); + //*fLog << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photons] + + // 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 = 0; + Double_t sigma2 = 0; + Double_t diff = 0; + Double_t addSig2 = 0; + Double_t elNoise2Pix = 0; + + + for (UInt_t i=0; iGetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; + Double_t oldsigma = ppix.GetMeanRms(); + 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 + + hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + if ( hdiff->GetEntries() == 0 ) + { + *fLog << "MPadONOFF::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; m<20; m++) + { + count += 1; + diff = hdiff->GetRandom(); + // the following condition ensures that elNoise2Pix > 0.0 + + if ( (diff + sigmabar2 - oldsigma2/ratioArea + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + + *fLog << "theta, j, count, sigmabar, diff = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; + } + delete hdiff; + sigma2 = diff + sigmabar2; + break; + + case 2 : + // throw the Sigma for this pixel from the distribution fHSigmaPixTheta + + hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + if ( hsig->GetEntries() == 0 ) + { + *fLog << "MPadONOFF::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; m<20; m++) + { + count += 1; + + sig = hsig->GetRandom(); + 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 << "theta, j, count, sigmabar, sig = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << sig << endl; + sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; + } + delete hsig; + break; + } + + //--------------------------------- + // get the additional sigma^2 for this pixel (due to the padding) + + addSig2 = sigma2*ratioArea - oldsigma2; + + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; + + + //--------------------------------- + // throw actual number of additional NSB photons (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) + { + sigmaNSB0 = sqrt( arg ); + } + else + { + *fLog << "MPadONOFF::Process(); argument of sqrt < 0 : " + << arg << endl; + sigmaNSB0 = 0.0000001; + } + + + //--------------------------------- + // smear NSB0 according to sigmaNSB0 + // and subtract lambdabar because of AC coupling + + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; + + //--------------------------------- + + // add additional NSB to the number of photons + Double_t oldphotons = pix.GetNumPhotons(); + Double_t newphotons = oldphotons + NSB; + pix.SetNumPhotons( newphotons ); + + + fHNSB->Fill( NSB/sqrt(ratioArea) ); + fHPhotons->Fill( oldphotons/sqrt(ratioArea), newphotons/sqrt(ratioArea) ); + + + // error: add sigma of padded noise quadratically + Double_t olderror = pix.GetErrorPhot(); + Double_t newerror = sqrt( olderror*olderror + addSig2 ); + pix.SetErrorPhot( newerror ); + + + Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + ppix.SetMeanRms( newsigma ); + + fHSigmaPedestal->Fill( oldsigma, newsigma ); + } + //---------- end of loop over pixels --------------------------------- + + // Calculate sigmabar again and crosscheck + + //fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << "after padding : " << endl; + //fSigmabar->Print(""); + + rc = 0; + fErrors[rc]++; + return kTRUE; + //****************************************************************** +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MPadONOFF::PostProcess() +{ + if (GetNumExecutions() != 0) + { + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *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("(Input) 2D : Sigmabar, \\Theta"); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + // draw the 3D histogram (input): Theta, pixel, Sigma^2-Sigmabar^2 + + c.cd(5); + TH2D *l1; + l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); + l1->SetDirectory(NULL); + l1->SetTitle("(Input) 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("(Input) 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 (input): Theta, pixel, Sigma + + c.cd(6); + TH2D *k1; + k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); + k1->SetDirectory(NULL); + k1->SetTitle("(Input) 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("(Input) 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("(Input) Sigmabarnew vs. Sigmabarold (ON, all \\Theta)"); + m1->SetXTitle("pixel"); + m1->SetYTitle("Sigma"); + + m1->DrawCopy("box"); + m1->SetBit(kCanDelete);; + + c.cd(12); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->SetTitle("(Input) Sigmabarnew vs. Sigmabarold (OFF, all \\Theta)"); + m2->SetXTitle("pixel"); + m2->SetYTitle("Sigma"); + + m2->DrawCopy("box"); + m2->SetBit(kCanDelete);; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2152) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2152) @@ -0,0 +1,91 @@ +#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 fRunType; + Int_t fGroup; + + Int_t fErrors[9]; + + // plots used for the padding + 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) + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHDiffPixTheta; // 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 + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + void SetPadFlag(Int_t padflag); + + ClassDef(MPadONOFF, 0) // task for the ON-OFF padding +}; + +#endif + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2151) +++ trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2152) @@ -392,5 +392,4 @@ hblind = fHBlindPixIdTheta->ProjectionY("", binPix, binPix, ""); - fBlinds->Clear(); if ( hblind->GetEntries() > 0.0 ) for (UInt_t i=0; iSetSelectedPad(NULL); @@ -801,4 +800,22 @@ fHSigmaTheta->DrawCopy(); fHSigmaTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2151) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2152) @@ -23,5 +23,7 @@ # INCLUDES = -I. -I../mbase -I../mmc -I../mraw -I../mgeom -I../mfilter \ - -I../mdata -I../mhist -I../mgui -I../mimage -I../mhistmc + -I../mdata -I../mhist -I../mgui -I../mimage -I../mhistmc \ + -I../mfileio + #------------------------------------------------------------------------------ @@ -54,4 +56,5 @@ MPadding.cc \ MPadSchweizer.cc \ + MPadONOFF.cc \ MCT1PointingCorrCalc.cc \ MParameters.cc \ Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 2151) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 2152) @@ -234,4 +234,5 @@ // Boolean bdontusepix[iMAXNUMPIX]; // bdontusepix is set true if the pixel should not be used in image analysis, otherwise it is true; + fBlinds->Clear(); *fLog << "Don't use pixels: "; @@ -245,5 +246,4 @@ *fLog << "Exclude pixels: "; - // Boolean bexcludepix[iMAXNUMPIX]; for (int i=0; iIsBlind(i) ) + { + fBlnd[i] = 1; + } + else + { + fBlnd[i] = 0; + } + } fBlinds->SetReadyToSave(); @@ -881,4 +895,13 @@ */ + // reset blind pixels for this event + fBlinds->Clear(); + for (int i=0; iSetPixelBlind(i); + } + + // reset pedestal RMS for this event for (Int_t i=0; iread((Byte_t*)&event, sizeof(struct eventrecord)); Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 2151) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 2152) @@ -4,4 +4,8 @@ #ifndef ROOT_TArrayF #include +#endif + +#ifndef ROOT_TArrayI +#include #endif @@ -58,4 +62,5 @@ TArrayF fPedRMS; + TArrayI fBlnd; Bool_t OpenNextFile();