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; } //--------------------------------------------------------------------- }