Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 1963) +++ trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 1963) @@ -0,0 +1,2347 @@ + +#include "CT1EEst.C" + +void InitBinnings(MParList *plist) + { + cout << "InitBinnings" << endl; + + 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 *binsht = new MBinning("BinningHeadTail"); + binsht->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsht); + + 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); + + Int_t nbins = 8; + TArrayD xbins(nbins+1); + xbins[0] = 0.0; + xbins[1] = 7.5; + xbins[2] = 17.5; + xbins[3] = 23.5; + xbins[4] = 29.5; + xbins[5] = 35.5; + xbins[6] = 42.0; + xbins[7] = 55.0; + xbins[8] = 68.0; + + MBinning *binth = new MBinning("BinningTheta"); + binth->SetEdges(xbins); + plist->AddToList(binth); + + MBinning *binsdiff = new MBinning("BinningDiffsigma2"); + binsdiff->SetEdges(100, -5.0, 20.0); + plist->AddToList(binsdiff); + } + +void DeleteBinnings(MParList *plist) + { + cout << "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("BinningHeadTail"); + 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("BinningDiffsigma2"); + if (!bin) delete bin; + } + +void CT1Analysis() +{ + //----------------------------------------------- + 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/optionB/"; + + // path for output from Mars + TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/"; + + + + //************************************************************************ + + // Job A_ON : read ON data + // - generate sigmabar vs. Theta plot; + // - write root file for ON data (ON1.root); + // - generate hadron matrix for g/h separation + + Bool_t JobA_ON = kFALSE; + Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? + Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ? + Bool_t WON1 = kFALSE; // write out root file ON1.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); + // - generate gamma matrix for g/h separation) + + Bool_t JobA_MC = kFALSE; + Bool_t WLookMC = kFALSE; // write out look-alike events for gammas ? + Bool_t WMC1 = kFALSE; // write out root file MC1.root ? + + + // Job B_NN_UP : read ON1.root (or MC1.root) file + // - read hadron and gamma matrix + // - calculate hadroness for method of Nearest Neighbors + // and for the supercuts + // - update the input files with the hadronesses (ON2.root or MC2.root) + + Bool_t JobB_NN_UP = kFALSE; + Bool_t WNN = kFALSE; // write root file ? + + + // Job B_EST_UP : + // - read MC2.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 ON2.root and MC2.root files with estimated energy + // (ON_XX2.root and MC_XX2.root) + + Bool_t JobB_EST_UP = kFALSE; + Bool_t WESTUP = kFALSE; // update root file ? + + + // Job C_NN : read ON1.root and MC1.root files + // - read look-alike events for hadrons and gammas + // - calculate hadronness for method of NEAREST NEIGHBOR + // - produce Neyman-Pearson plot + + Bool_t JobC_NN = kFALSE; + + + // Job C_RF : read ON1.root and MC1.root files + // - read look-alike events for hadrons and gammas + // - calculate hadronness for method of RANDOM FOREST + // - produce Neyman-Pearson plot + + Bool_t JobC_RF = kFALSE; + + + // Job D_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 JobD_XX = kTRUE; + Bool_t WXXD = kFALSE; // write out root file ON3.root ? + + + //************************************************************************ + + + //--------------------------------------------------------------------- + // Job A_ON + //========= + // 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 look-alike events (matrix_hadrons_ON.root) + // (to be used later together with the corresponding MC gamma file + // for the g/h separation in the analysis of the ON 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_ON) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job A_ON" << endl; + cout << "" << endl; + cout << "Macro CT1Analysis : JobA_ON, WHistON, WLookON, WON1 = " + << JobA_ON << ", " << WHistON << ", " << WLookON << ", " + << WON1 << endl; + + TString typeInput = "ON"; + cout << "typeInput = " << typeInput << endl; + + // name of input root file + TString filenamein(onfile); + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "1.root"; + + + //----------------------------------------------------------- + MTaskList tliston; + MParList pliston; + MFilterList flist; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + + //------------------------------------------- + // create the tasks which should be executed + // + + MCT1ReadPreProc read(filenamein); + + MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", + "MCT1PointingCorrCalc"); + + MBlindPixelCalc blind; + + MFCT1SelBasic selbasic; + MContinue contbasic(&selbasic); + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); + + MImgCleanStd clean; + + // calculate also extended image parameters + TString fHilName = "MHillas"; + MHillasExt hext(fHilName); + + MHillasCalc hcalc(fHilName); + + TString fHilSrcName = "MHillasSrc"; + MHillasSrc hilsrc(fHilSrcName); + + MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", ""); + csrc1.SetInput(fHilName); + + TString fnewparName = "MNewImagePar"; + MNewImagePar newpar(fnewparName); + + MNewImageParCalc newimage(sourceName, fnewparName); + newimage.SetInput(fHilName); + + MFCT1SelStandard selstandard(fHilName, fHilSrcName); + MContinue contstandard(&selstandard); + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetTitle("Task to plot the source independent Hillas parameters"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetTitle("Task to plot the star map (points along main axis of shower)"); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + hfill3.SetTitle("Task to plot the extended Hillas parameters"); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); + + MFillH hfill5("MHNewImagePar", fnewparName); + hfill5.SetTitle("Task to plot the new image parameters"); + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + // prepare writing of look-alike events (for hadrons) + + const char* mtxName = "MatrixHadrons"; + Int_t howMany = 500000; + TString outName = outPath; + outName += "matrix_hadrons_"; + outName += typeInput; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Matrix for (hadrons)" << endl; + cout << "input file = " << filenamein<< endl; + cout << "matrix name = " << mtxName << endl; + cout << "max. no.of look-alike events = " << howMany << endl; + cout << "name of output root file = " << outName << endl; + cout << "" << endl; + + + // matrix limitation for look-alike events (approximate number) + MFEventSelector selector("", ""); + selector.SetNumSelectEvts(howMany); + + // + // --- setup the matrix and add it to the parlist + // + MHMatrix *pmatrix = new MHMatrix(mtxName); + MHMatrix& matrix = *pmatrix; + + matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrix.AddColumn("MSigmabar.fSigmabar"); + matrix.AddColumn("log10(MHillas.fSize)"); + matrix.AddColumn("MHillasSrc.fDist"); + matrix.AddColumn("MHillas.fWidth"); + matrix.AddColumn("MHillas.fLength"); + matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)"); + matrix.AddColumn("MHillas.fConc"); + matrix.AddColumn("MNewImagePar.fLeakage1"); + + matrix.Print("cols"); + + MFillH fillmatg(mtxName); + fillmatg.SetFilter(&flist); + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + + if (WON1) + { + MWriteRootFile write(outNameImage); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + } + + //***************************** + // entries in MFilterList + + flist.AddToList(&selector); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&source); + pliston.AddToList(&sigbarcalc); + + pliston.AddToList(&hext); + pliston.AddToList(&hilsrc); + pliston.AddToList(&newpar); + pliston.AddToList(&hhilext); + + pliston.AddToList(pmatrix); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + tliston.AddToList(&pointcorr); + tliston.AddToList(&blind); + + tliston.AddToList(&contbasic); + tliston.AddToList(&sigbarcalc); + tliston.AddToList(&fillsigtheta); + tliston.AddToList(&clean); + + tliston.AddToList(&hcalc); + tliston.AddToList(&csrc1); + tliston.AddToList(&newimage); + + tliston.AddToList(&flist); + tliston.AddToList(&fillmatg); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contstandard); + if (WON1) + tliston.AddToList(&write); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + //evtloop.Write(); + + //Int_t maxevents = 1000000000; + Int_t maxevents = 20000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + c = ((MHSigmaTheta*)pliston.FindObject("SigmaTheta", "MHSigmaTheta")) + ->DrawClone(); + + c = pliston.FindObject("MHHillas")->DrawClone(); + c = pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + c = pliston.FindObject("MHNewImagePar")->DrawClone(); + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + // + // ---------------------------------------------------------- + // Definition of the reference sample of look-alike events + // (this is a subsample of the original sample) + // ---------------------------------------------------------- + // + cout << "" << endl; + cout << "========================================================" << 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->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; + } + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Macro CT1Analysis : define reference sample for hadrons" << endl; + cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrix.EnableGraphicalOutput(); + Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts); + + if ( !matrixok ) + { + cout << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl; + return; + } + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + //------------------------------------------- + // write out look-alike events (for hadrons) + // + if (WLookON) + { + TFile *writejens = new TFile(outName, "RECREATE", ""); + matrix.Write(); + + cout << "" << endl; + cout << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file " + << outName << endl; + + writejens->Close(); + delete writejens; + } + + + //------------------------------------------- + // Write histograms onto a file + if (WHistON) + { + MHSigmaTheta *sigtheta = + (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); + if (!sigtheta) + { + cout << "Object 'SigmaTheta' not found" << endl; + return; + } + TH2D *fHSigTh = sigtheta->GetSigmaTheta(); + TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); + TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); + + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta_"; + outNameSigTh += typeInput; + outNameSigTh += ".root"; + + TFile outfile(outNameSigTh, "RECREATE"); + fHSigTh->Write(); + fHSigPixTh->Write(); + fHDifPixTh->Write(); + outfile.Close(); + + cout << "" << endl; + cout << "File " << outNameSigTh << " was written out" << endl; + } + + + cout << "Macro CT1Analysis : End of Job A_ON" << endl; + cout << "===================================================" << 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 look-alike events (matrix_gammas_MC.root) + // (to be used later together with the corresponding hadron file + // for the g/h separation in the analysis 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) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job A_MC" << endl; + + cout << "" << endl; + cout << "Macro CT1Analysis : JobA_MC, WLookMC, WMC1 = " + << JobA_MC << ", " << WLookMC << ", " << WMC1 << endl; + + TString typeInput = "MC"; + cout << "typeInput = " << typeInput << endl; + + + // name of input root file + TString filenamein(mcfile); + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "1.root"; + + + // use for padding sigmabar vs. Theta from ON data + TString typeHist = "ON"; + cout << "typeHist = " << typeHist << endl; + + //------------------------------------ + // Get the histograms "2D-ThetaSigmabar" + // and "3D-ThetaPixSigma" + // and "3D-ThetaPixDiff" + + + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta_"; + outNameSigTh += typeHist; + outNameSigTh += ".root"; + + + cout << "Reading in file " << outNameSigTh << endl; + + TFile *infile = new TFile(outNameSigTh); + infile->ls(); + + TH2D *fHSigmaTheta = + (TH2D*) gROOT.FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + cout << "Object '2D-ThetaSigmabar' not found on root file" << endl; + return; + } + cout << "Object '2D-ThetaSigmabar' was read in" << endl; + + TH3D *fHSigmaPixTheta = + (TH3D*) gROOT.FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + cout << "Object '3D-ThetaPixSigma' not found on root file" << endl; + return; + } + cout << "Object '3D-ThetaPixSigma' was read in" << endl; + + TH3D *fHDiffPixTheta = + (TH3D*) gROOT.FindObject("3D-ThetaPixDiff"); + if (!fHSigmaPixTheta) + { + cout << "Object '3D-ThetaPixDiff' not found on root file" << endl; + return; + } + cout << "Object '3D-ThetaPixDiff' was read in" << endl; + + //------------------------------------ + + MTaskList tlist; + MParList plist; + MFilterList flist; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + //------------------------------------------- + // create the tasks which should be executed + // + + MCT1ReadPreProc read(filenamein); + + MBlindPixelCalc blind; + + MFCT1SelBasic selbasic; + MContinue contbasic(&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)", + fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta); + padthomas.SetPadFlag(1); + + //........................................... + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); + + MImgCleanStd clean; + + // calculate also extended image parameters + TString fHilName = "MHillas"; + MHillasExt hext(fHilName); + + MHillasCalc hcalc(fHilName); + + TString fHilSrcName = "MHillasSrc"; + MHillasSrc hilsrc(fHilSrcName); + + MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", ""); + csrc1.SetInput(fHilName); + + TString fnewparName = "MNewImagePar"; + MNewImagePar newpar(fnewparName); + + MNewImageParCalc newimage(sourceName, fnewparName); + newimage.SetInput(fHilName); + + MFCT1SelStandard selstandard(fHilName, fHilSrcName); + MContinue contstandard(&selstandard); + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetTitle("Task to plot the source independent Hillas parameters"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetTitle("Task to plot the star map (points along main axis of shower)"); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + hfill3.SetTitle("Task to plot the extended Hillas parameters"); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + hfill4.SetTitle("Task to plot the source dependent Hillas parameters"); + + MFillH hfill5("MHNewImagePar", fnewparName); + hfill5.SetTitle("Task to plot the new image parameters"); + + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + // prepare writing of look-alike events (for gammas) + + const char* mtxName = "MatrixGammas"; + Int_t howMany = 500000; + + TString outName = outPath; + outName += "matrix_gammas_"; + outName += typeInput; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Matrix for (gammas)" << endl; + cout << "input file = " << filenamein<< endl; + cout << "matrix name = " << mtxName << endl; + cout << "max. no.of look-alike events = " << howMany << endl; + cout << "name of output root file = " << outName << endl; + cout << "" << endl; + + + // matrix limitation for look-alike events (approximate number) + MFEventSelector selector("", ""); + selector.SetNumSelectEvts(howMany); + + // + // --- setup the matrix and add it to the parlist + // + MHMatrix *pmatrix = new MHMatrix(mtxName); + MHMatrix& matrix = *pmatrix; + + matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrix.AddColumn("MSigmabar.fSigmabar"); + matrix.AddColumn("log10(MHillas.fSize)"); + matrix.AddColumn("MHillasSrc.fDist"); + matrix.AddColumn("MHillas.fWidth"); + matrix.AddColumn("MHillas.fLength"); + matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)"); + matrix.AddColumn("MHillas.fConc"); + matrix.AddColumn("MNewImagePar.fLeakage1"); + + matrix.Print("cols"); + + MFillH fillmatg(mtxName); + fillmatg.SetFilter(&flist); + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + + + if (WMC1) + { + MWriteRootFile write(outNameImage); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + } + + //***************************** + // entries in MFilterList + + flist.AddToList(&selector); + + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + InitBinnings(&plist); + + plist.AddToList(&source); + plist.AddToList(&sigbarcalc); + + plist.AddToList(&hext); + plist.AddToList(&hilsrc); + plist.AddToList(&newpar); + plist.AddToList(&hhilext); + + plist.AddToList(pmatrix); + + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&blind); + tlist.AddToList(&padthomas); + + tlist.AddToList(&contbasic); + tlist.AddToList(&sigbarcalc); + tlist.AddToList(&fillsigtheta); + tlist.AddToList(&clean); + + tlist.AddToList(&hcalc); + tlist.AddToList(&csrc1); + tlist.AddToList(&newimage); + + tlist.AddToList(&flist); + tlist.AddToList(&fillmatg); + + 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.SetProgressBar(&bar); + //evtloop.Write(); + + //Int_t maxevents = 1000000000; + Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tlist.PrintStatistics(0, kTRUE); + DeleteBinnings(&plist); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + c = ((MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta")) + ->DrawClone(); + c = plist.FindObject("MHHillas")->DrawClone(); + c = plist.FindObject("MHHillasExt")->DrawClone(); + c = plist.FindObject("MHHillasSrc")->DrawClone(); + c = plist.FindObject("MHNewImagePar")->DrawClone(); + c = plist.FindObject("MHStarMap")->DrawClone(); + + + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + // + // ---------------------------------------------------------- + // Definition of the reference sample of look-alike events + // (this is a subsample of the original sample) + // ---------------------------------------------------------- + // + cout << "" << endl; + cout << "========================================================" << 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); + 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; + } + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Macro CT1Analysis : define reference sample for gammas" << endl; + cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + matrix.EnableGraphicalOutput(); + Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts); + if ( !matrixok ) + { + cout << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl; + return; + } + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + //------------------------------------------- + // write out look-alike events (for gammas) + // + if (WLookMC) + { + TFile *writejens = new TFile(outName, "RECREATE", ""); + matrix.Write(); + cout << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file " + << outName << endl; + + writejens->Close(); + delete writejens; + } + + cout << "Macro CT1Analysis : End of Job A_MC" + << endl; + cout << "=========================================================" + << 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) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job B_NN_UP" << endl; + + cout << "" << endl; + cout << "Macro CT1Analysis : JobB_NN_UP, WNN = " + << JobB_NN_UP << ", " << WNN << endl; + + + //TString typeInput = "ON"; + TString typeInput = "MC"; + cout << "typeInput = " << typeInput << endl; + + // name of input root file + TString filenameData = outPath; + filenameData += typeInput; + filenameData += "1.root"; + cout << "filenameData = " << filenameData << endl; + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "2.root"; + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + outNameImage = filenameData; + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + cout << "outNameImage = " << outNameImage << endl; + + TString typeMatrixHadrons = "ON"; + cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + //************************************************************************* + // read in matrices of look-alike events + + const char* mtxName = "MatrixGammas"; + + TString outName = outPath; + outName += "matrix_gammas_"; + outName += "MC"; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Get matrix for (gammas)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outName' + // + TFile *fileg = new TFile(outName); + + MHMatrix matrixg; + matrixg.Read(mtxName); + matrixg.Print("cols"); + + delete fileg; + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + TString outName = outPath; + outName += "matrix_hadrons_"; + outName += typeMatrixHadrons; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << " Get matrix for (hadrons)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name mtxName + // + TFile *fileh = new TFile(outName); + + MHMatrix matrixh; + matrixh.Read(mtxName); + matrixh.Print("cols"); + + delete fileh; + + //************************************************************************* + + + MTaskList tliston; + MParList pliston; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + TString fHilName = "MHillas"; + TString fHilSrcName = "MHillasSrc"; + TString fnewparName = "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, fHilSrcName); + sccalc.SetHadronnessName(hadSCName); + + //....................................................................... + + if (WNN) + { + //MWriteRootFile write(outNameImage); + MWriteRootFile write(outNameImage, "UPDATE"); + + //write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MTime", "Events"); + //write.AddContainer("MMcEvt", "Events"); + //write.AddContainer("MSrcPosCam", "Events"); + //write.AddContainer("MSigmabar", "Events"); + //write.AddContainer("MHillas", "Events"); + //write.AddContainer("MHillasSrc", "Events"); + //write.AddContainer("MNewImagePar", "Events"); + + write.AddContainer(hadNNName, "Events"); + write.AddContainer(hadSCName, "Events"); + } + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + + Float_t hadcut = 0.40; + Float_t alphacut = 100.0; + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(hadcut, alphacut); + selfinalgh.SetHadronnessName(hadNNName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + + Float_t hadcut = 0.40; + Float_t alphacut = 20.0; + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(hadcut, alphacut); + selfinal.SetHadronnessName(hadNNName); + MContinue contfinal(&selfinal); + + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + MFillH hfill4("MHHillasSrc", fHilSrcName); + MFillH hfill5("MHNewImagePar", fnewparName); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&matrixg); + pliston.AddToList(&matrixh); + + pliston.AddToList(&hhilext); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&nncalc); + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + + if (WNN) + tliston.AddToList(&write); + + 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); + //evtloop.Write(); + + if ( !evtloop.Eventloop() ) + //if ( !evtloop.Eventloop(1000) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + + c = pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + c = pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + c = pliston.FindObject("MHHillas")->DrawClone(); + c = pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + c = pliston.FindObject("MHNewImagePar")->DrawClone(); + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + cout << "Macro CT1Analysis : End of Job B_NN_UP" << endl; + cout << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + //--------------------------------------------------------------------- + // Job B_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 (JobB_EST_UP) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job B_EST_UP" << endl; + + cout << "" << endl; + cout << "Macro CT1Analysis : JobB_EST_UP, WESTUP = " + << JobB_EST_UP << ", " << WESTUP << endl; + + + TString typeON = "ON"; + TString typeMC = "MC"; + TString ext = "3.root"; + + //------------------------------ + // name of MC file to be used for optimizing the energy estimator + TString filenameOpt(outPath); + filenameOpt += typeMC; + filenameOpt += ext; + cout << "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; + cout << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness and of |alpha| + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + cout << "Maximum values of hadronness and |alpha| = " + << maxhadronness << ", " << maxalpha << endl; + + // name of file containing the parameters of the energy estimator + TString energyParName(outPath); + energyParName += "energyest_"; + energyParName += XX; + energyParName += ".root"; + cout << "energyParName = " << energyParName << endl; + + + //------------------------------ + // name of ON file to be updated + TString filenameON(outPath); + filenameON += typeON; + filenameON += ext; + cout << "filenameON = " << filenameON << endl; + + // name of MC file to be updated + TString filenameMC(outPath); + filenameMC += typeMC; + filenameMC += ext; + cout << "filenameMC = " << filenameMC << endl; + + //------------------------------ + // name of updated ON file + TString filenameONup(outPath); + filenameONup += typeON; + filenameONup += "_"; + filenameONup += XX; + filenameONup += ext; + cout << "filenameONup = " << filenameONup << endl; + + // name of updated MC file + TString filenameMCup(outPath); + filenameMCup += typeMC; + filenameMCup += "_"; + filenameMCup += XX; + filenameMCup += ext; + cout << "filenameMCup = " << filenameMCup << endl; + + //----------------------------------------------------------- + + TString fHilName = "MHillas"; + TString fHilSrcName = "MHillasSrc"; + TString fnewparName = "MNewImagePar"; + + //=========================================================== + // + // Optimization of energy estimator + // + + TString inpath(""); + TString outpath(""); + Int_t howMany = 2000; + CT1EEst(inpath, filenameOpt, outpath, energyParName, + fHilName, fHilSrcName, fhadronnessName, + howMany, maxhadronness, maxalpha); + + //----------------------------------------------------------- + // + // Read in parameters of energy estimator + // + + TFile enparam(energyParName); + + TArrayD parA(5); + TArrayD parB(7); + for (Int_t i=0; iFindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "MHillas"; + TString fHilSrcName = "MHillasSrc"; + TString fnewparName = "MNewImagePar"; + + Float_t hadcut = 0.36; + Float_t alphacut = 100.0; + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(hadcut, alphacut); + selfinalgh.SetHadronnessName(hadName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadnn("MHHadronness", hadName); + + Float_t hadcut = 0.36; + Float_t alphacut = 20.0; + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(hadcut, alphacut); + selfinal.SetHadronnessName(hadName); + MContinue contfinal(&selfinal); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + MFillH hfill4("MHHillasSrc", fHilSrcName); + MFillH hfill5("MHNewImagePar", fnewparName); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&matrixg); + pliston.AddToList(&matrixh); + + pliston.AddToList(&hhilext); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&fgamma); + tliston.AddToList(&fhadron); + tliston.AddToList(&nncalc); + tliston.AddToList(&fillhadnn); + + 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); + //evtloop.Write(); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + c = pliston.FindObject("MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + c = pliston.FindObject("MHHillas")->DrawClone(); + c = pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + c = pliston.FindObject("MHNewImagePar")->DrawClone(); + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + + cout << "Macro CT1Analysis : End of Job C_NN" << endl; + cout << "===================================================" << endl; + } + + + + //--------------------------------------------------------------------- + // Job C_RF + //========= + + // read ON1 and MC1 data files + // + // - calculate hadronness for method of Random Forest + // - produce Neyman-Pearson plot + + if (JobC_RF) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job C_RF" << endl; + + cout << "" << endl; + cout << "Macro CT1Analysis : JobC_RF = " << JobC_RF << endl; + + + TString typeMC = "MC"; + cout << "typeMC = " << typeMC << endl; + + TString typeData = "ON"; + cout << "typeData = " << typeData << endl; + + TString typeMatrixHadrons = "ON"; + cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + //************************************************************************* + // read in matrices of look-alike events + + const char* mtxName = "MatrixGammas"; + + TString outName = outPath; + outName += "matrix_gammas_"; + outName += "MC"; + outName += ".root"; + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Get matrix for (gammas)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outName' + // + TFile *fileg = new TFile(outName); + + MHMatrix matrixg; + matrixg.Read(mtxName); + matrixg.Print("cols"); + + delete fileg; + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + TString outName = outPath; + outName += "matrix_hadrons_"; + outName += typeMatrixHadrons; + outName += ".root"; + + cout << "" << endl; + cout << "========================================================" << endl; + cout << " Get matrix for (hadrons)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name mtxName + // + TFile *fileh = new TFile(outName); + + MHMatrix matrixh; + matrixh.Read(mtxName); + matrixh.Print("cols"); + + delete fileh; + + //***************************************************************** + + //----------------------------------------------------------------- + // grow the trees of the random forest (event loop = tree loop) + + cout << "Macro CT1Analysis : start growing trees" << endl; + + MTaskList tlist2; + MParList plist2; + plist2.AddToList(&tlist2); + + plist2.AddToList(&matrixg); + plist2.AddToList(&matrixh); + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(100); + rfgrow2.SetNumTry(4); + rfgrow2.SetNdSize(10); + + TString outRF = outPath; + outRF += "RF.root"; + 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 + MRanForest *fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + MRanTree *fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + //----------------------------------------------------------------- + + MTaskList tliston; + MParList pliston; + + //------------------------------------------- + // create the tasks which should be executed + // + TString filenameData = outPath; + filenameData += typeData; + filenameData += "1.root"; + + TString filenameMC = outPath; + filenameMC += typeMC; + filenameMC += "1.root"; + + cout << "filenameData = " << filenameData << endl; + cout << "filenameMC = " << filenameMC << endl; + + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + MFParticleId fgamma ("MMcEvt", '=', kGAMMA); + MFParticleId fhadron("MMcEvt", '!', kGAMMA); + + //....................................................................... + // calculate hadronnes for method of Random Forest + + TString hadName = "HadRF"; + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadName); + + //....................................................................... + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "MHillas"; + TString fHilSrcName = "MHillasSrc"; + TString fnewparName = "MNewImagePar"; + + Float_t hadcut = 0.36; + Float_t alphacut = 100.0; + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(hadcut, alphacut); + selfinalgh.SetHadronnessName(hadName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadrf("MHHadronness", hadName); + MFillH fillhrf("MHRanForest"); + + Float_t hadcut = 0.36; + Float_t alphacut = 20.0; + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(hadcut, alphacut); + selfinal.SetHadronnessName(hadName); + MContinue contfinal(&selfinal); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + MFillH hfill4("MHHillasSrc", fHilSrcName); + MFillH hfill5("MHNewImagePar", fnewparName); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(fRanForest); + pliston.AddToList(fRanTree); + + pliston.AddToList(&hhilext); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&fgamma); + tliston.AddToList(&fhadron); + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillhadrf); + tliston.AddToList(&fillhrf); + + 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); + //evtloop.Write(); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + c = pliston.FindObject("MHRanForest")->DrawClone(); + c = pliston.FindObject("MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + + c = pliston.FindObject("MHHillas")->DrawClone(); + c = pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + c = pliston.FindObject("MHNewImagePar")->DrawClone(); + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + + cout << "Macro CT1Analysis : End of Job C_RF" << endl; + cout << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job D_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 (JobD_XX) + { + cout << "=====================================================" << endl; + cout << "Macro CT1Analysis : Start of Job D_XX" << endl; + + cout << "" << endl; + cout << "Macro CT1Analysis : JobD_XX, WXXD = " + << JobD_XX << ", " << WXXD << endl; + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + cout << "typeData = " << typeData << endl; + + TString typeMC = "MC"; + TString ext = "3.root"; + TString extout = "4.root"; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + TString XX("NN"); + //TString XX("SC"); + //TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + cout << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness and of |alpha| + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + cout << "Maximum values of hadronness and |alpha| = " + << maxhadronness << ", " << maxalpha << endl; + + + //------------------------------ + // name of MC file to be used for calculating the eff. collection areas + TString filenameArea(outPath); + filenameArea += typeMC; + filenameArea += "_"; + filenameArea += XX; + filenameArea += ext; + cout << "filenameArea = " << filenameArea << endl; + + //------------------------------ + // name of file containing the eff. collection areas + TString collareaName(outPath); + collareaName += "area_"; + collareaName += XX; + collareaName += ".root"; + cout << "collareaName = " << collareaName << endl; + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += "_"; + filenameData += XX; + filenameData += ext; + cout << "filenameData = " << filenameData << endl; + + //------------------------------ + // name of output data file (after the final cuts) + TString filenameDataout(outPath); + filenameDataout += typeData; + filenameDataout += "_"; + filenameDataout += XX; + filenameDataout += extout; + cout << "filenameDataout = " << filenameDataout << endl; + + + //==================================================================== + cout << "-----------------------------------------------" << endl; + cout << "Start calculation of effective collection areas" << endl; + MParList parlist; + MTaskList tasklist; + + //--------------------------------------- + // Setup the tasks to be executed + // + MReadMarsFile reader("Events", filenameArea); + + MFCT1SelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha); + + MContinue conthadrons(&cuthadrons); + + MHMcCT1CollectionArea* collarea = + new MHMcCT1CollectionArea("MHMcCT1CollectionArea","",30,2.,5.); + MMcCT1CollectionAreaCalc effi; + + //******************************** + // entries in MParList + + parlist.AddToList(&tasklist); + parlist.AddToList(collarea); + + //******************************** + // entries in MTaskList + + tasklist.AddToList(&reader); + tasklist.AddToList(&conthadrons); + tasklist.AddToList(&effi); + + //******************************** + + //----------------------------------------- + // Execute event loop + // + MEvtLoop magic; + magic.SetParList(&parlist); + + MProgressBar bar; + magic.SetProgressBar(&bar); + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(0, kTRUE); + + // Display the histograms + // + parlist.FindObject("MHMcCT1CollectionArea")->DrawClone("lego"); + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + f.Close(); + + + cout << "Calculation of effective collection areas done" << endl; + cout << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + TString hadName("Had"); + hadName += XX; + + TString fHilName = "MHillas"; + TString fHilSrcName = "MHillasSrc"; + TString fnewparName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + //....................................................................... + + if (WXXD) + { + MWriteRootFile write(filenameDataout); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + write.AddContainer("HadNN", "Events"); + write.AddContainer("HadSC", "Events"); + + write.AddContainer("MEnergyEst", "Events"); + } + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + MFCT1SelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetCuts(maxhadronness, 100.0); + selfinalgh.SetHadronnessName(fhadronnessName); + MContinue contfinalgh(&selfinalgh); + + MFillH fillhadnn("MHadNN[MHHadronness]", "HadNN"); + MFillH fillhadsc("MHadSC[MHHadronness]", "HadSC"); + //MFillH fillhadrf("MHadRF[MHHadronness]", "HadRF"); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName); + + MFillH hfill3("MHHillasExt", fHilSrcName); + MFillH hfill4("MHHillasSrc", fHilSrcName); + MFillH hfill5("MHNewImagePar", fnewparName); + + MFCT1SelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetCuts(maxhadronness, maxalpha); + selfinal.SetHadronnessName(fhadronnessName); + MContinue contfinal(&selfinal); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&hhilext); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + if (WXXD) + tliston.AddToList(&write); + + 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); + //evtloop.Write(); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + DeleteBinnings(&pliston); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + + c = pliston.FindObject("MHadNN", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + c = pliston.FindObject("MHadRF", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + c = pliston.FindObject("MHadSC", "MHHadronness"); + if (c) + { + c->DrawClone(); + c->Print(); + } + + c = pliston.FindObject("MHHillas")->DrawClone(); + c = pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + c = pliston.FindObject("MHNewImagePar")->DrawClone(); + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + cout << "Macro CT1Analysis : End of Job D_XX" << endl; + cout << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + +} + + + + + + + Index: trunk/MagicSoft/Mars/macros/CT1EEst.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1EEst.C (revision 1963) +++ trunk/MagicSoft/Mars/macros/CT1EEst.C (revision 1963) @@ -0,0 +1,618 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 09/2002 +! Abelardo Moralejo, 03/2003 +! Wolfgang Wittek, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ +// +//------------------------------------------------------------------------ +// +// fcn calculates the function to be minimized (using TMinuit::Migrad) +// +void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit(); + + MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); + + MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); + MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); + + eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); + + evtloop->Eventloop(); + + f = eval->GetChisq(); +} +//------------------------------------------------------------------------ + +void CT1EgyEst() +{ + TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + TString fileNameIn = "MC2.root"; + + TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + TString energyParName = "energyest.root"; + + TString hilName = "MHillas"; + TString hilSrcName = "MHillasSrc"; + + TString hadronnessName = "MHadronness"; + + Int_t howMany = 2200; + + Float_t maxhadronness = 0.4; + Float_t maxalpha = 20.0; + + CT1EEst(inPath, fileNameIn, + outPath, energyParName, + hilName, hilSrcName, hadronnessName, + howMany, maxhadronness, maxalpha); +} + + + +//------------------------------------------------------------------------ +// +// Arguments of CT1EEst : +// ------------------------ +// +// inPath, fileNameIn path and name of input file (MC gamma events) +// outPath, energyParName path and name of file containing the parameters +// of the energy estimator +// hilName, hilSrcName names of "MHillas" and "MHillasSrc" containers +// hadronnessName name of container holding the hadronness +// howMany no.of events to be read from input file +// maxhadronness maximum value of hadronness to be accepted +// maxalpha maximum value of |ALPHA| to be accepted +// +void CT1EEst(TString inPath, TString fileNameIn, + TString outPath, TString energyParName, + TString hilName, TString hilSrcName, TString hadronnessName, + Int_t howMany, + Float_t maxhadronness, Float_t maxalpha) +{ + cout << "================================================================" + << endl; + cout << "Start of energy estimation part" << endl; + cout << "" << endl; + cout << "CT1EEst input values : " << endl; + cout << "---------------------- " << endl; + cout << "inPath, fileNameIn = " + << inPath << ", " << fileNameIn << endl; + cout << "outPath, energyParName = " + << outPath << ", " << energyParName << endl; + cout << "hilName, hilSrcName, hadronnessName = " << hilName << ", " + << hilSrcName << ", " << hadronnessName << endl; + cout << "howMany, maxhadronness, maxalpha = " << howMany << ", " + << maxhadronness << ", " << maxalpha << endl; + cout << "" << endl; + + + // convert from [deg] to [mm] + // + MGeomCamCT1 gct1; + //maxdist /= gct1.GetConvMm2Deg(); + + + //========================================================================== + // Start of Part 1 (determination of the parameters of the energy estimator) + // + + //----------------------------------------------------------------------- + // Fill events into a MHMatrix, + // to be used later in the minimization by MINUIT + // + + MParList parlist; + + TString inputfile(outPath); + inputfile += fileNameIn; + + MFEventSelector selector; + selector.SetNumSelectEvts(howMany); + + cout << "Read events from file '" << inputfile << "'" << endl; + MReadTree read("Events"); + read.AddFile(inputfile); + read.DisableAutoScheme(); + read.SetSelector(&selector); + + MFCT1SelFinal filterhadrons; + filterhadrons.SetCuts(maxhadronness, maxalpha); + filterhadrons.SetHadronnessName(hadronnessName); + filterhadrons.SetInverted(); + + cout << "Define columns of matrix" << endl; + MHMatrix matrix; + Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); + Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact"); + + if (colenergy < 0 || colimpact < 0) + { + cout << "colenergy, colimpact = " << colenergy << ", " + << colimpact << endl; + return; + } + + MEnergyEstParam eest(hilName); + eest.Add(hilSrcName); + eest.InitMapping(&matrix); + + cout << "--------------------------------------" << endl; + cout << "Fill events into the matrix" << endl; + if ( !matrix.Fill(&parlist, &read, &filterhadrons) ) + return; + cout << "Matrix was filled with " << matrix.GetNumRows() + << " events" << endl; + + //----------------------------------------------------------------------- + // + // Setup the eventloop which will be executed in the fcn of MINUIT + // + cout << "--------------------------------------" << endl; + cout << "Setup event loop to be used in 'fcn'" << endl; + + MParList parlist; + MTaskList tasklist; + + MMatrixLoop loop(&matrix); + + + MChisqEval eval; + eval.SetY1(new MDataElement(&matrix, colenergy)); + eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); + eval.SetOwner(); + + //******************************** + // Entries in MParList + + parlist.AddToList(&tasklist); + + //******************************** + // Entries in MTaskList + + tasklist.AddToList(&loop); + tasklist.AddToList(&eest); + tasklist.AddToList(&eval); + + //******************************** + + cout << "Event loop was setup" << endl; + MEvtLoop evtloop; + evtloop.SetParList(&parlist); + + + + //.......................................................................... + + //---------- Start of minimization part -------------------- + // + // Do the minimization with MINUIT + // + // Be careful: This is not thread safe + // + cout << "--------------------------------------" << endl; + cout << "Start minimization" << endl; + + TMinuit minuit(12); + minuit.SetPrintLevel(-1); + minuit.SetFCN(fcn); + + // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg) + + if (minuit.SetErrorDef(1)) + { + cout << "SetErrorDef failed." << endl; + return; + } + + // + // Set initial values of the parameters (close enough to the final ones, taken + // from previous runs of the procedure). Parameter fA[4] is not used in the default + // energy estimation model (from D. Kranich). + // + TArrayD fA(5); + TArrayD fB(7); + + fA[0] = 21006.2; + fA[1] = -43.2648; + fA[2] = -690.403; + fA[3] = -0.0428544; + fA[4] = 1.; + fB[0] = -3391.05; + fB[1] = 136.58; + fB[2] = 0.253807; + fB[3] = 254.363; + fB[4] = 61.0386; + fB[5] = -0.0190984; + fB[6] = -421695; + + // + // Set starting values and step sizes for parameters + // + for (Int_t i=0; iWrite("EnergyParams"); + outfile.Close(); + + cout << "--------------------------------------" << endl; + cout << "Write parameters of energy estimator onto file '" + << paramout << endl; + cout << "--------------------------------------" << endl; + // + // End of Part 1 (determination of the parameters of the energy estimator) + //========================================================================= + + + //////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////// + + + //========================================================================== + // Start of Part 2 (test quality of energy estimation) + // + // + + //-------------------------------------------- + // Read the parameters of the energy estimator + // + + TString paramin(outPath); + paramin += energyParName; + + cout << "--------------------------------------" << endl; + cout << "Read parameters of energy estimator from file '" + << paramin << endl; + TFile enparam(paramin); + + //======================================================================= + // Setup the event loop + // + cout << "--------------------------------------" << endl; + cout << "Setup event loop for checking quality of energy estimation" << endl; + + + MTaskList tlist2; + MParList parlist2; + + //----------------------------------------------- + // Read events + // + + TString inputfile(inPath); + inputfile += fileNameIn; + + cout << "Read events from file '" << inputfile << "'" << endl; + MReadTree read2("Events"); + read2.AddFile(inputfile); + read2.DisableAutoScheme(); + + + //----------------------------------------------- + // Select events + // + cout << "Select events with hadronness < " << maxhadronness + << " and |alpha| < " << maxalpha << endl; + MFCT1SelFinal hcut2; + hcut2.SetCuts(maxhadronness, maxalpha); + hcut2.SetHadronnessName(hadronnessName); + + MContinue cont(&hcut2); + + + //----------------------------------------------- + // Create some histograms to display the results + // + + MH3 mh3ed ("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1.0"); + MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1.0"); + MH3 mhe ("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)"); + + mhe.SetName("HistEE"); + mh3ed.SetName("HistEnergyDelta"); + mh3ed2.SetName("HistEnergyDelta"); + + MBinning binsedx("BinningHistEnergyDeltaX"); + MBinning binsedy("BinningHistEnergyDeltaY"); + MBinning binseex("BinningHistEEX"); + MBinning binseey("BinningHistEEY"); + + binsedx.SetEdges(35, 2.0, 5.5); + binsedy.SetEdges(40,-1.5, 2.5); + + binseex.SetEdges(35, 2.0, 5.5); + binseey.SetEdges(35, 2.0, 5.5); + + MFillH hfilled(&mh3ed); + MFillH hfilled2(&mh3ed2); + MFillH hfillee(&mhe); + + + //----------------------------------------------- + // Create energy estimator task, add necessary containers and + // initialize parameters read from file: + // + + MEnergyEstParam eest2(hilName); + eest2.Add(hilSrcName); + + TArrayD parA(5); + TArrayD parB(7); + + for (Int_t i=0; iGetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); + mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)"); + mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2); + + TCanvas *c = new TCanvas("energy","CT1 Energy estimation"); + c->Divide(2,2); + c->cd(1); + energy_1->SetBottomMargin(0.12); + mhe.DrawClone("PROFXnonew"); + + mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); + mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); + mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2); + mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5); + c->cd(2); + energy_2->SetLeftMargin(0.15); + energy_2->SetBottomMargin(0.12); + mh3ed.DrawClone("PROFXnonew"); + + mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)"); + mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); + mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2); + mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5); + c->cd(3); + energy_3->SetLeftMargin(0.15); + energy_3->SetBottomMargin(0.12); + mh3ed2.DrawClone("PROFXnonew"); + + ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae"); + + c->cd(4); + energy_4->SetLeftMargin(0.1); + energy_4->SetRightMargin(0.05); + energy_4->SetBottomMargin(0.12); + deltae.GetXaxis()->SetTitleOffset(1.2); + deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}"); + deltae.Draw("e"); + + //------------------------------------------- + // + // include the histograms + // into the root file containing the parameters of the energy estimator + // + + + TString paramout(outPath); + paramout += energyParName; + + TFile out2(paramout, "UPDATE"); + ((TH2*)mh3ed.GetHist())->Write("mh3ed"); + ((TH2*)mh3ed2.GetHist())->Write("mh3ed2"); + ((TH2*)mhe.GetHist())->Write("mhe"); + deltae.Write(); + out2.Close(); + + cout << "Quality histograms were added onto the file '" << paramout << endl; + + cout << "" << endl; + cout << "End of energy estimation part" << endl; + cout << "================================================================" + << endl; + // + // End of Part 2 (test quality of energy estimation) + //========================================================================== + + + return; + +} +//============================================================================ + + + + + + + + +