Index: trunk/MagicSoft/Mars/macros/CT1Analysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1Analysis.C (revision 2072) +++ (revision ) @@ -1,2347 +1,0 @@ - -#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; - } - //--------------------------------------------------------------------- - -} - - - - - - -