Index: unk/MagicSoft/Mars/macros/AnalyseCT1.C =================================================================== --- /trunk/MagicSoft/Mars/macros/AnalyseCT1.C (revision 2685) +++ (revision ) @@ -1,2050 +1,0 @@ -void AnalyseCT1() -{ - //----------------------------------------------- - //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; - //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; - //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6"; - //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc"; - //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc"; - - //----------------------------------------------- - const char *offfile = "~/Mars200203/CT1data/offdata.preproc"; - - const char *onfile = "~/Mars200203/CT1data/mkn421_on.preproc"; - //const char *onfile = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6"; - - //const char *mcfile = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5 - //const char *mcfile = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6"; - const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc"; - - //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6"; - //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303"; - - //----------------------------------------------- - //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc"; - //const char *onfile = "~magican/home/ct1test/PreprocData/mkn421_on.preproc"; - //const char *mcfile = "~magican/home/ct1test/PreprocData/mc_gammas.preproc"; - - - const char *filename; - - //************************************************************************ - // Job 1 : read OFF data - // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1); - // generate hadron matrix for g/h separation) - - Bool_t Job1 = kFALSE; - Bool_t WHistOFF = kTRUE; // write out histogram sigmabar vs. Theta ? - Bool_t WLookOFF = kTRUE; // write out look-alike events for hadrons ? - Bool_t WOFF1 = kTRUE; // write out root file OFF1 ? - - - // Job 2 : read ON data - // (generate sigmabar vs. Theta plot; write root file for ON data (ON1)) - - Bool_t Job2 = 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 ? - - // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data - // (do padding; write root file for MC gammas (MC1); - // generate gamma matrix for g/h separation) - - Bool_t Job3 = kFALSE; - Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? - Bool_t WMC1 = kTRUE; // write out root file MC1 ? - - - - // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix - // (produce the Neyman-Pearson plot) - Bool_t Job4 = kTRUE; - - - // Job 5 : read MC1 file, read hadron and gamma matrix - // (do g/h separation; write root file for MC gammas after final cuts (MC2)) - Bool_t Job5 = kFALSE; - Bool_t WMC2 = kFALSE; // write out root file MC2 ? - - - // Job 6 : read ON data, read hadron and gamma matrix - // (do g/h separation; write root file for ON data after final cuts (ON2)) - Bool_t Job6 = kFALSE; - Bool_t WON2 = kTRUE; // write out root file ON2 ? - //************************************************************************ - - - - - //--------------------------------------------------------------------- - // Job 1 - //====== - - // read OFF data file - - // - produce the 2D-histogram "sigmabar versus Theta" for OFF data - // (to be used for the padding of the MC gamma data) - - // - write a file of look-alike events (hadron matrix) - // (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 OFF events (OFF1) - // (after the standard cuts, before the g/h separation) - // (to be used together with the corresponding MC gamma file - // for the optimization of the g/h separation) - - if (Job1) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 1" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = " - << Job1 << ", " << WHistOFF << ", " << WLookOFF << ", " - << WOFF1 << endl; - - - TString typeData = "OFF"; - cout << "typeData = " << typeData << endl; - - MTaskList tliston; - MParList pliston; - pliston.AddToList(&tliston); - - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - pliston.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - pliston.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - pliston.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - pliston.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - pliston.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - pliston.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - pliston.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - pliston.AddToList(&binsdiff); - //:::::::::::::::::::::::::::::::::::::::::: - - - //------------------------------------------- - // create the tasks which should be executed - // - - filename = offfile; - RName = "ReadCT1data"; - MCT1ReadPreProc read(filename, RName); - - MSelBasic selbasic; - - MBlindPixelCalc blind; - blind.SetUseInterpolation(); - blind.SetUseBlindPixels(); - // blind.SetUseCetralPixel(); - - // create an object of class MSigmabar, - // because class MHSigmaTheta will look for it - MSigmabar sigbar; - pliston.AddToList(&sigbar); - MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); - fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); - - MImgCleanStd clean; - - // calculate also extended image parameters - TString fHilName = "Hillas"; - MHillasExt hext(fHilName); - pliston.AddToList(&hext); - MHillasExt *fHillas = &hext; - - // name of output container for MHillasCalc - // = name of MHillas object - MHillasCalc hcalc(fHilName); - - // name of output container for MHillasSrcCalc - // = name of MHillasSrc object - TString fHilSrcName = "HillasSrc"; - MHillasSrc hilsrc(fHilSrcName); - MHillasSrc *fHillasSrc = &hilsrc; - pliston.AddToList(&hilsrc); - MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); - - MSelStandard selstand(fHilName, fHilSrcName); - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - pliston.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - // prepare writing of look-alike events (for hadrons) - - const char* mtxName = "MatrixHadrons"; - Int_t howMany = 50000; - TString outName = "matrix_hadrons_"; - outName += typeData; - outName += ".root"; - - cout << "" << endl; - cout << "========================================================" << endl; - cout << "Matrix for (hadrons)" << endl; - cout << "input file = " << filename<< 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("", "", RName); - selector.SetNumSelectEvts(howMany); - MFilterList flist; - flist.AddToList(&selector); - - // - // --- 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(Hillas.fSize)"); - matrix.AddColumn("HillasSrc.fDist"); - matrix.AddColumn("Hillas.fWidth"); - matrix.AddColumn("Hillas.fLength"); - matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); - matrix.AddColumn("abs(Hillas.fM3Long)"); - matrix.AddColumn("Hillas.fConc"); - matrix.AddColumn("Hillas.fConc1"); - - pliston.AddToList(pmatrix); - - MFillH fillmatg(mtxName); - fillmatg.SetFilter(&flist); - - - if (WOFF1) - { - TString outNameImage = typeData; - outNameImage += "1.root"; - MWriteRootFile write(outNameImage); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("Hillas", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("HillasSrc", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MMcRunHeader","RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - } - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - - - //***************************** - // tasks to be executed - - tliston.AddToList(&read); - - tliston.AddToList(&selbasic); - tliston.AddToList(&blind); - tliston.AddToList(&fillsigtheta); - tliston.AddToList(&clean); - - tliston.AddToList(&hcalc); - tliston.AddToList(&csrc1); - - tliston.AddToList(&selstand); - if (WOFF1) - tliston.AddToList(&write); - - tliston.AddToList(&flist); - tliston.AddToList(&fillmatg); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - //------------------------------------------- - // Display the histograms - // - TObject *c; - c = pliston.FindObject("MHHillas")->DrawClone(); - - c = pliston.FindObject("MHHillasSrc")->DrawClone(); - - //pliston.FindObject("MHHillasExt")->DrawClone(); - c = pliston.FindObject(nxt)->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 1); - // set refcolumn negative if distribution is not to be changed - Int_t refcolumn = 1; - Int_t nbins = 10; - 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 AnalyseCT1 : define reference sample for hadrons" << endl; - cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) - { - cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl; - return; - } - //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - - //------------------------------------------- - // write out look-alike events (for hadrons) - // - if (WLookOFF) - { - TFile *writejens = new TFile(outName, "RECREATE", ""); - matrix.Write(); - cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file " - << outName << endl; - - writejens->Close(); - delete writejens; - } - - //------------------------------------------- - // Write histograms onto a file - if (WHistOFF) - { - 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 = "SigmaTheta_"; - outNameSigTh += typeData; - outNameSigTh += ".root"; - - TFile outfile(outNameSigTh, "RECREATE"); - fHSigTh->Write(); - fHSigPixTh->Write(); - fHDifPixTh->Write(); - outfile.Close(); - - cout << "File " << outNameSigTh << " was written out" << endl; - } - - - cout << "Macro AnalyseCT1 : End of Job 1" << endl; - cout << "===================================================" << endl; - } - //--------------------------------------------------------------------- - - //--------------------------------------------------------------------- - // Job 2 - //====== - // read ON data file - - // - produce the 2D-histogram "sigmabar versus Theta" for ON data - // (to be used for the padding of the MC gamma data) - - // - write a file of look-alike events (hadron matrix) - // (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) - // (after the standard cuts, before the g/h separation) - // (to be used together with the corresponding MC gamma file - // for the optimization of the g/h separation) - - if (Job2) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 2" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = " - << Job2 << ", " << WHistON << ", " << WLookON << ", " - << WON1 << endl; - - TString typeData = "ON"; - cout << "typeData = " << typeData << endl; - - MTaskList tliston; - MParList pliston; - pliston.AddToList(&tliston); - - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - pliston.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - pliston.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - pliston.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - pliston.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - pliston.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - pliston.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - pliston.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - pliston.AddToList(&binsdiff); - //:::::::::::::::::::::::::::::::::::::::::: - - - //------------------------------------------- - // create the tasks which should be executed - // - - filename = onfile; - RName = "ReadCT1data"; - MCT1ReadPreProc read(filename, RName); - - MSelBasic selbasic; - - MBlindPixelCalc blind; - blind.SetUseInterpolation(); - blind.SetUseBlindPixels(); - // blind.SetUseCetralPixel(); - - // create an object of class MSigmabar, - // because class MHSigmaTheta will look for it - MSigmabar sigbar; - pliston.AddToList(&sigbar); - MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); - fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); - - MImgCleanStd clean; - - // calculate also extended image parameters - TString fHilName = "Hillas"; - MHillasExt hext(fHilName); - pliston.AddToList(&hext); - MHillasExt *fHillas = &hext; - - // name of output container for MHillasCalc - // = name of MHillas object - MHillasCalc hcalc(fHilName); - - // name of output container for MHillasSrcCalc - // = name of MHillasSrc object - TString fHilSrcName = "HillasSrc"; - MHillasSrc hilsrc(fHilSrcName); - MHillasSrc *fHillasSrc = &hilsrc; - pliston.AddToList(&hilsrc); - MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); - - MSelStandard selstand(fHilName, fHilSrcName); - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - pliston.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - // prepare writing of look-alike events (for hadrons) - - const char* mtxName = "MatrixHadrons"; - Int_t howMany = 50000; - TString outName = "matrix_hadrons_"; - outName += typeData; - outName += ".root"; - - - cout << "" << endl; - cout << "========================================================" << endl; - cout << "Matrix for (hadrons)" << endl; - cout << "input file = " << filename<< 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("", "", RName); - selector.SetNumSelectEvts(howMany); - MFilterList flist; - flist.AddToList(&selector); - - // - // --- 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(Hillas.fSize)"); - matrix.AddColumn("HillasSrc.fDist"); - matrix.AddColumn("Hillas.fWidth"); - matrix.AddColumn("Hillas.fLength"); - matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); - matrix.AddColumn("abs(Hillas.fM3Long)"); - matrix.AddColumn("Hillas.fConc"); - matrix.AddColumn("Hillas.fConc1"); - - pliston.AddToList(pmatrix); - - MFillH fillmatg(mtxName); - fillmatg.SetFilter(&flist); - - - if (WON1) - { - TString outNameImage = typeData; - outNameImage += "1.root"; - MWriteRootFile write(outNameImage); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("Hillas", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("HillasSrc", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MMcRunHeader","RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - } - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - - - //***************************** - // tasks to be executed - - tliston.AddToList(&read); - - tliston.AddToList(&selbasic); - tliston.AddToList(&blind); - tliston.AddToList(&fillsigtheta); - tliston.AddToList(&clean); - - tliston.AddToList(&hcalc); - tliston.AddToList(&csrc1); - - tliston.AddToList(&selstand); - if (WON1) - tliston.AddToList(&write); - - tliston.AddToList(&flist); - tliston.AddToList(&fillmatg); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - //------------------------------------------- - // Display the histograms - // - TObject *c; - c = pliston.FindObject("MHHillas")->DrawClone(); - - c = pliston.FindObject("MHHillasSrc")->DrawClone(); - - //pliston.FindObject("MHHillasExt")->DrawClone(); - c = pliston.FindObject(nxt)->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 1); - // set refcolumn negative if distribution is not to be changed - Int_t refcolumn = 1; - Int_t nbins = 10; - 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 AnalyseCT1 : define reference sample for hadrons" << endl; - cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) - { - cout << "Macro AnalyseCT1 : 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 << "Macro AnalyseCT1 : 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 = "SigmaTheta_"; - outNameSigTh += typeData; - outNameSigTh += ".root"; - - TFile outfile(outNameSigTh, "RECREATE"); - fHSigTh->Write(); - fHSigPixTh->Write(); - fHDifPixTh->Write(); - outfile.Close(); - - cout << "File " << outNameSigTh << " was written out" << endl; - } - - - cout << "Macro AnalyseCT1 : End of Job 2" << endl; - cout << "===================================================" << endl; - } - - - //--------------------------------------------------------------------- - // Job 3 - //====== - - // read MC gamma data - - // - to pad them - // (using the 2D-histogram "sigmabar versus Theta" of the OFF data) - - // - to write a file of look-alike events (gammas matrix) - // (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) - // (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 (Job3) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 3" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = " - << Job3 << ", " << WLookMC << ", " << WMC1 << endl; - - - TString typeMC = "MC"; - cout << "typeMC = " << typeMC << endl; - - // use for padding sigmabar vs. Theta from ON or OFF data - TString typeData = "ON"; - //TString typeData = "OFF"; - cout << "typeData = " << typeData << endl; - - //------------------------------------ - // Get the histograms "2D-ThetaSigmabar" - // and "3D-ThetaPixSigma" - // and "3D-ThetaPixDiff" - - - TString outNameSigTh = "SigmaTheta_"; - outNameSigTh += typeData; - 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; - - plist.AddToList(&tlist); - - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - plist.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - plist.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - plist.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - plist.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - plist.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - plist.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - plist.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - plist.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - plist.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - plist.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - plist.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - plist.AddToList(&binsdiff); - - //:::::::::::::::::::::::::::::::::::::::::: - - //------------------------------------------- - // create the tasks which should be executed - // - filename = mcfile; - readname = "ReadCT1MC"; - MCT1ReadPreProc read(filename, readname); - - MSelBasic selbasic; - - MBlindPixelCalc blind; - blind.SetUseInterpolation(); - blind.SetUseBlindPixels(); - // blind.SetUseCetralPixel(); - - // 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); - - //........................................... - - MImgCleanStd clean; //(0, 0); - - // calculate also extended image parameters - TString fHilName = "Hillas"; - MHillasExt hext(fHilName); - plist.AddToList(&hext); - MHillasExt *fHillas = &hext; - - // name of output container for MHillasCalc - // = name of MHillas object - MHillasCalc hcalc(fHilName); - - // name of output container for MHillasSrcCalc - // = name of MHillasSrc object - TString fHilSrcName = "HillasSrc"; - MHillasSrc hilsrc(fHilSrcName); - MHillasSrc *fHillasSrc = &hilsrc; - plist.AddToList(&hilsrc); - MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); - - MSelStandard selstand(fHilName, fHilSrcName); - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - plist.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - // prepare writing of look-alike events (for gammas) - - const char* mtxName = "MatrixGammas"; - Int_t howMany = 50000; - - TString outName = "matrix_gammas_"; - outName += typeMC; - outName += "_"; - outName += typeData; - outName += ".root"; - - - cout << "" << endl; - cout << "========================================================" << endl; - cout << "Matrix for (gammas)" << endl; - cout << "input file = " << filename<< 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("", "", readname); - selector.SetNumSelectEvts(howMany); - MFilterList flist; - flist.AddToList(&selector); - - // - // --- 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(Hillas.fSize)"); - matrix.AddColumn("HillasSrc.fDist"); - matrix.AddColumn("Hillas.fWidth"); - matrix.AddColumn("Hillas.fLength"); - matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); - matrix.AddColumn("abs(Hillas.fM3Long)"); - matrix.AddColumn("Hillas.fConc"); - matrix.AddColumn("Hillas.fConc1"); - - plist.AddToList(pmatrix); - - MFillH fillmatg(mtxName); - fillmatg.SetFilter(&flist); - - - if (WMC1) - { - TString outNameImage = typeMC; - outNameImage += "_"; - outNameImage += typeData; - outNameImage += "1.root"; - MWriteRootFile write(outNameImage); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("Hillas", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("HillasSrc", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MMcRunHeader","RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - } - - - //+++++++++++++++++++++++++++++++++++++++++++++++++++ - - - //***************************** - // tasks to be executed - - tlist.AddToList(&read); - - tlist.AddToList(&selbasic); - tlist.AddToList(&blind); - tlist.AddToList(&padthomas); - tlist.AddToList(&clean); - - tlist.AddToList(&hcalc); - tlist.AddToList(&csrc1); - - tlist.AddToList(&selstand); - if (WMC1) - tlist.AddToList(&write); - - tlist.AddToList(&flist); - tlist.AddToList(&fillmatg); - - tlist.AddToList(&hfill1); - tlist.AddToList(&hfill2); - tlist.AddToList(&hfill3); - tlist.AddToList(&hfill4); - - - //***************************** - - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&plist); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 100; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tlist.PrintStatistics(0, kTRUE); - - //------------------------------------------- - // Display the histograms - // - plist.FindObject("MHHillas")->DrawClone(); - plist.FindObject("MHHillasSrc")->DrawClone(); - - //plist.FindObject("MHHillasExt")->DrawClone(); - plist.FindObject(nxt)->DrawClone(); - - 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; - // set refcolumn negative if distribution is not to be changed - Int_t refcolumn = 1; - Int_t nbins = 10; - 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 AnalyseCT1 : define reference sample for gammas" << endl; - cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = " - << refcolumn << ", " << nmaxevts << endl; - - if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) - { - cout << "Macro AnalyseCT1 : 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 AnalyseCT1 : matrix of look-alike events for gammas written onto file " - << outName << endl; - - writejens->Close(); - delete writejens; - } - - cout << "Macro AnalyseCT1 : End of Job 3" - << endl; - cout << "=========================================================" - << endl; - } - - - - //--------------------------------------------------------------------- - // Job 4 - //====== - - // read OFF1 and MC1 data files - // - // - produce Neyman-Pearson plot - - if (Job4) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 4" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job4 = " << Job4 << endl; - - - TString typeMC = "MC"; - cout << "typeMC = " << typeMC << endl; - - // use hadron matrix from ON or OFF data - TString typeData = "ON"; - //TString typeData = "OFF"; - cout << "typeData = " << typeData << endl; - - //************************************************************************* - // read in matrices of look-alike events - - const char* mtxName = "MatrixGammas"; - - - TString outName = "matrix_gammas_"; - outName += typeMC; - outName += "_"; - outName += typeData; - 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); - pmatrixg = &matrixg; - - delete fileg; - - - //***************************************************************** - - const char* mtxName = "MatrixHadrons"; - - TString outName = "matrix_hadrons_"; - outName += typeData; - 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); - pmatrixh = &matrixh; - - delete fileh; - - //***************************************************************** - - MTaskList tliston; - MParList pliston; - pliston.AddToList(&tliston); - - pliston.AddToList(pmatrixg); - pliston.AddToList(pmatrixh); - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - pliston.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - pliston.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - pliston.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - pliston.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - pliston.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - pliston.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - pliston.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - pliston.AddToList(&binsdiff); - //:::::::::::::::::::::::::::::::::::::::::: - - - //------------------------------------------- - // create the tasks which should be executed - // - - TString filenameData = typeData; - filenameData += "1.root"; - - TString filenameMC = typeMC; - filenameMC += "_"; - filenameMC += typeData; - filenameMC += "1.root"; - - - cout << "filenameData = " << filenameData << endl; - cout << "filenameMC = " << filenameMC << endl; - - MReadMarsFile read("Events", filenameMC); - read.AddFile(filenameData); - read.DisableAutoScheme(); - - //....................................................................... - // g/h separation - - MMultiDimDistCalc ghsep; - ghsep.SetUseNumRows(25); - ghsep.SetUseKernelMethod(kFALSE); - - //....................................................................... - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - TString fHilName = "Hillas"; - TString fHilSrcName = "HillasSrc"; - - Float_t hadcut = 0.2; - MSelFinal selfinalgh(fHilName, fHilSrcName); - selfinalgh.SetHadronnessCut(hadcut); - selfinalgh.SetAlphaCut(100.0); - - MFillH fillh("MHHadronness"); - - MSelFinal selfinal(fHilName, fHilSrcName); - selfinal.SetHadronnessCut(hadcut); - selfinal.SetAlphaCut(20.0); - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - pliston.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - - //***************************** - // tasks to be executed - - tliston.AddToList(&read); - - tliston.AddToList(&ghsep); - tliston.AddToList(&fillh); - - tliston.AddToList(&selfinalgh); - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - - tliston.AddToList(&selfinal); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - TObject *c; - - c = pliston.FindObject("MHHadronness")->DrawClone(); - c->Print(); - - //c = pliston.FindObject("MHHillas")->DrawClone(); - - c = pliston.FindObject("MHHillasSrc")->DrawClone(); - - //pliston.FindObject("MHHillasExt")->DrawClone(); - //c = pliston.FindObject(nxt)->DrawClone(); - - c = pliston.FindObject("MHStarMap")->DrawClone(); - - - - cout << "Macro AnalyseCT1 : End of Job 4" << endl; - cout << "===================================================" << endl; - } - - - //--------------------------------------------------------------------- - // Job 5 - //====== - - // - read in the matrices of look-alike events for gammas and hadrons - - // then read MC1 data file - // - // - perform the g/h separation - // - apply the final cuts - // - // - write a file of MC gamma events (MC2) - // (after the g/h separation and after the final cuts) - - if (Job5) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 5" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job5, WMC2 = " - << Job5 << ", " << WMC2 << endl; - - TString typeInput = "MC"; - cout << "typeInput = " << typeInput << endl; - - TString typeMC = "MC"; - cout << "typeMC = " << typeMC << endl; - - // use hadron matrix from ON or OFF data - TString typeData = "ON"; - //TString typeData = "OFF"; - cout << "typeData = " << typeData << endl; - - //************************************************************************* - // read in matrices of look-alike events - - const char* mtxName = "MatrixGammas"; - - TString outName = "matrix_gammas_"; - outName += typeMC; - outName += "_"; - outName += typeData; - 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); - pmatrixg = &matrixg; - - delete fileg; - - //***************************************************************** - - const char* mtxName = "MatrixHadrons"; - - TString outName = "matrix_hadrons_"; - outName += typeData; - 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); - pmatrixh = &matrixh; - - delete fileh; - - //************************************************************************* - - - MTaskList tliston; - MParList pliston; - pliston.AddToList(&tliston); - - pliston.AddToList(pmatrixg); - pliston.AddToList(pmatrixh); - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - pliston.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - pliston.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - pliston.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - pliston.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - pliston.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - pliston.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - pliston.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - pliston.AddToList(&binsdiff); - //:::::::::::::::::::::::::::::::::::::::::: - - - //------------------------------------------- - // create the tasks which should be executed - // - - TString filenameMC = typeInput; - filenameMC += "_"; - filenameMC += typeData; - filenameMC += "1.root"; - - readname = "ReadCT1MCdata"; - MReadMarsFile read("Events", filenameMC); - read.DisableAutoScheme(); - - - //....................................................................... - // g/h separation - - MMultiDimDistCalc ghsep; - ghsep.SetUseNumRows(25); - ghsep.SetUseKernelMethod(kFALSE); - - - - //....................................................................... - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - TString fHilName = "Hillas"; - TString fHilSrcName = "HillasSrc"; - - Float_t hadcut = 0.2; - MSelFinal selfinalgh(fHilName, fHilSrcName); - selfinalgh.SetHadronnessCut(hadcut); - selfinalgh.SetAlphaCut(100.0); - - MFillH fillh("MHHadronness"); - - MSelFinal selfinal(fHilName, fHilSrcName); - selfinal.SetHadronnessCut(hadcut); - selfinal.SetAlphaCut(20.0); - - if (WMC2) - { - TString outNameImage = typeInput; - outNameImage += "_"; - outNameImage += typeData; - outNameImage += "2.root"; - MWriteRootFile write(outNameImage); - write.AddContainer("MHadronness", "Events"); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("Hillas", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("HillasSrc", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MMcRunHeader","RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - } - - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - pliston.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - - - //***************************** - // tasks to be executed - - tliston.AddToList(&read); - - tliston.AddToList(&ghsep); - tliston.AddToList(&fillh); - - tliston.AddToList(&selfinalgh); - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - - tliston.AddToList(&selfinal); - if (WMC2) - tliston.AddToList(&write); - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - TObject *c; - - c = pliston.FindObject("MHHadronness")->DrawClone(); - c->Print(); - - c = pliston.FindObject("MHHillas")->DrawClone(); - - c = pliston.FindObject("MHHillasSrc")->DrawClone(); - - //pliston.FindObject("MHHillasExt")->DrawClone(); - c = pliston.FindObject(nxt)->DrawClone(); - - c = pliston.FindObject("MHStarMap")->DrawClone(); - - - - cout << "Macro AnalyseCT1 : End of Job 5" << endl; - cout << "===================================================" << endl; - } - - - //--------------------------------------------------------------------- - // Job 6 - //====== - - // - read in the matrices of look-alike events for gammas and hadrons - - // then read ON1 data file - // - // - perform the g/h separation - // - apply the final cuts - // - // - write a file of ON events (ON2) - // (after the g/h separation and after the final cuts) - // (to be used for the flux calculation) - // - // - do the flux calculation - - if (Job6) - { - cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of Job 6" << endl; - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Job6, WON2 = " - << Job6 << ", " << WON2 << endl; - - TString typeInput = "ON"; - cout << "typeInput = " << typeInput << endl; - - TString typeMC = "MC"; - cout << "typeMC = " << typeMC << endl; - - // use hadron matrix from ON or OFF data - TString typeData = "ON"; - //TString typeData = "OFF"; - cout << "typeData = " << typeData << endl; - - - //************************************************************************* - // read in matrices of look-alike events - - const char* mtxName = "MatrixGammas"; - - TString outName = "matrix_gammas_"; - outName += typeMC; - outName += "_"; - outName += typeData; - 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); - pmatrixg = &matrixg; - - delete fileg; - - //***************************************************************** - - const char* mtxName = "MatrixHadrons"; - - TString outName = "matrix_hadrons_"; - outName += typeData; - 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); - pmatrixh = &matrixh; - - delete fileh; - - //************************************************************************* - - - MTaskList tliston; - MParList pliston; - pliston.AddToList(&tliston); - - pliston.AddToList(pmatrixg); - pliston.AddToList(pmatrixh); - - //:::::::::::::::::::::::::::::::::::::::::: - - MBinning binssize("BinningSize"); - binssize.SetEdgesLog(50, 10, 1.0e5); - pliston.AddToList(&binssize); - - MBinning binsdistc("BinningDist"); - binsdistc.SetEdges(50, 0, 1.4); - pliston.AddToList(&binsdistc); - - MBinning binswidth("BinningWidth"); - binswidth.SetEdges(50, 0, 1.0); - pliston.AddToList(&binswidth); - - MBinning binslength("BinningLength"); - binslength.SetEdges(50, 0, 1.0); - pliston.AddToList(&binslength); - - MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(100, -100, 100); - pliston.AddToList(&binsalpha); - - MBinning binsasym("BinningAsym"); - binsasym.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsasym); - - MBinning binsht("BinningHeadTail"); - binsht.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsht); - - MBinning binsm3l("BinningM3Long"); - binsm3l.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3l); - - MBinning binsm3t("BinningM3Trans"); - binsm3t.SetEdges(50, -1.5, 1.5); - pliston.AddToList(&binsm3t); - - - //..... - MBinning binsb("BinningSigmabar"); - binsb.SetEdges( 100, 0.0, 5.0); - pliston.AddToList(&binsb); - - MBinning binth("BinningTheta"); - binth.SetEdges( 70, -0.5, 69.5); - pliston.AddToList(&binth); - - MBinning binsdiff("BinningDiffsigma2"); - binsdiff.SetEdges(100, -5.0, 20.0); - pliston.AddToList(&binsdiff); - //:::::::::::::::::::::::::::::::::::::::::: - - - //------------------------------------------- - // create the tasks which should be executed - // - - TString filenameData = typeInput; - filenameData += "1.root"; - - readname = "ReadCT1ONdata"; - MReadMarsFile read("Events", filenameData); - read.DisableAutoScheme(); - - - //....................................................................... - // g/h separation - - MMultiDimDistCalc ghsep; - ghsep.SetUseNumRows(25); - ghsep.SetUseKernelMethod(kFALSE); - - - - //....................................................................... - - // geometry is needed in MHHillas... classes - MGeomCam *fGeom = - (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); - - TString fHilName = "Hillas"; - TString fHilSrcName = "HillasSrc"; - - Float_t hadcut = 0.2; - MSelFinal selfinalgh(fHilName, fHilSrcName); - selfinalgh.SetHadronnessCut(hadcut); - selfinalgh.SetAlphaCut(100.0); - - MFillH fillh("MHHadronness"); - - MSelFinal selfinal(fHilName, fHilSrcName); - selfinal.SetHadronnessCut(hadcut); - selfinal.SetAlphaCut(20.0); - - if (WON2) - { - TString outNameImage = typeInput; - outNameImage += "_"; - outNameImage += typeData; - outNameImage += "2.root"; - MWriteRootFile write(outNameImage); - write.AddContainer("MHadronness", "Events"); - write.AddContainer("MSigmabar", "Events"); - write.AddContainer("Hillas", "Events"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("HillasSrc", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MMcRunHeader","RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - } - - - MFillH hfill1("MHHillas", fHilName); - MFillH hfill2("MHStarMap", fHilName); - - TString nxt = "HillasExt"; - MHHillasExt hhilext(nxt, "", fHilName); - pliston.AddToList(&hhilext); - TString namext = nxt; - namext += "[MHHillasExt]"; - MFillH hfill3(namext, fHilSrcName); - - MFillH hfill4("MHHillasSrc", fHilSrcName); - - - - //***************************** - // tasks to be executed - - tliston.AddToList(&read); - - tliston.AddToList(&ghsep); - tliston.AddToList(&fillh); - - tliston.AddToList(&selfinalgh); - - tliston.AddToList(&hfill1); - tliston.AddToList(&hfill2); - tliston.AddToList(&hfill3); - tliston.AddToList(&hfill4); - - tliston.AddToList(&selfinal); - if (WON2) - tliston.AddToList(&write); - - - - //***************************** - - //------------------------------------------- - // Execute event loop - // - MProgressBar bar; - MEvtLoop evtloop; - evtloop.SetParList(&pliston); - evtloop.SetProgressBar(&bar); - - Int_t maxevents = 1000000000; - //Int_t maxevents = 1000; - if ( !evtloop.Eventloop(maxevents) ) - return; - - tliston.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // Display the histograms - // - TObject *c; - - c = pliston.FindObject("MHHadronness")->DrawClone(); - c->Print(); - - //c = pliston.FindObject("MHHillas")->DrawClone(); - - c = pliston.FindObject("MHHillasSrc")->DrawClone(); - - ////pliston.FindObject("MHHillasExt")->DrawClone(); - //c = pliston.FindObject(nxt)->DrawClone(); - - c = pliston.FindObject("MHStarMap")->DrawClone(); - - - - cout << "Macro AnalyseCT1 : End of Job 6" << endl; - cout << "=======================================================" << endl; - } - //--------------------------------------------------------------------- -} - - - - - - - - - -