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 *filename; Bool_t Data = kFALSE; // loop over data ? Bool_t WLook = kFALSE; // write out look-alike events for hadrons ? Bool_t WHist = kFALSE; // write out histograms for padding ? Bool_t MC = kTRUE; // loop over MC gamma data ? Bool_t RHist = kTRUE; // read in histograms for padding ? Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? cout << "" << endl; cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ", " << WHist << endl; cout << " RHist, MC = " << RHist << ", " << MC << endl; cout << "" << endl; //--------------------------------------------------------------------- // start of section for ON data // (in the CT1 analysis the ON data are also used as a sample representing // the hadrons for the g/h separation) // read ON data file // - to produce the 2D-histogram "sigmabar versus Theta" // (to be used for the padding of the MC gamma data) // - to write a file of look-alike events (for hadrons) // (to be used later together with the corresponding MC gamma file // for the g/h separation in the analysis of the ON data) // - to write a file of ON events (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 (Data) { cout << "=====================================================" << endl; cout << "Macro AnalyseCT1 : Start of section for ON data" << 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; readname = "ReadCT1data"; MCT1ReadPreProc read(filename, readname); 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(fHillas, fHillasSrc); 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 = 30000; char* outName = "matrix_hadrons.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("", "", 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("Hillas.fWidth"); matrix.AddColumn("Hillas.fLength"); matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); matrix.AddColumn("abs(Hillas.fAsym)"); matrix.AddColumn("abs(Hillas.fM3Long)"); matrix.AddColumn("abs(Hillas.fM3Trans)"); matrix.AddColumn("abs(HillasSrc.fHeadTail)"); matrix.AddColumn("Hillas.fConc"); matrix.AddColumn("Hillas.fConc1"); matrix.AddColumn("HillasSrc.fDist"); matrix.AddColumn("log10(Hillas.fSize)"); matrix.AddColumn("HillasSrc.fAlpha"); matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); //matrix.AddColumn("MMcEvt.fTheta"); pliston.AddToList(pmatrix); MFillH fillmatg(mtxName); fillmatg.SetFilter(&flist); //+++++++++++++++++++++++++++++++++++++++++++++++++++ MSelFinal selfinal(fHillas, fHillasSrc); //***************************** // 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); tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); tliston.AddToList(&hfill3); tliston.AddToList(&hfill4); tliston.AddToList(&flist); tliston.AddToList(&fillmatg); //tliston.AddToList(&selfinal); //***************************** //------------------------------------------- // Execute event loop // MEvtLoop evtloop; evtloop.SetParList(&pliston); Int_t maxevents = 1000000000; if ( !evtloop.Eventloop(maxevents) ) return; tliston.PrintStatistics(0, kTRUE); //------------------------------------------- // Display the histograms // pliston.FindObject("MHHillas")->DrawClone(); pliston.FindObject("MHHillasSrc")->DrawClone(); //pliston.FindObject("MHHillasExt")->DrawClone(); pliston.FindObject(nxt)->DrawClone(); 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; // set refcolumn negative if distribution is not to be changed Int_t refcolumn = 12; 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 (WLook) { 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 (WHist) { MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); if (!sigtheta) { *fLog << "Object 'SigmaTheta' not found" << endl; return; } TH2D *fHSigTh = sigtheta->GetSigmaTheta(); TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); TFile outfile("SigmaTheta.root", "RECREATE"); fHSigTh->Write(); fHSigPixTh->Write(); fHDifPixTh->Write(); outfile.Close(); cout << "File 'SigmaTheta.root' was written out" << endl; } cout << "Macro AnalyseCT1 : End of section for ON data" << endl; cout << "===================================================" << endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- // start of section for MC gamma data // 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 (for gammas) // (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 // (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) // - to write a file of padded MC gamma events (after the final cuts) // (to be used for the optimization of the energy estimator // and for calculating the effective collection areas) if (MC) { cout << "=============================================================" << endl; cout << "Macro : AnalyseCT1 : Start of section for MC gamma data" << endl; MTaskList tlist; MParList plist; plist.AddToList(&tlist); //------------------------------------ // Get the histograms "2D-ThetaSigmabar" // and "3D-ThetaPixSigma" // and "3D-ThetaPixDiff" if (RHist) { cout << "Reading in file 'SigmaTheta.root' " << endl; TFile *infile = new TFile("SigmaTheta.root"); infile->ls(); TH2D *fHSigmaTheta = (TH2D*) gROOT.FindObject("2D-ThetaSigmabar"); if (!fHSigmaTheta) { *fLog << "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) { *fLog << "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) { *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; return; } cout << "Object '3D-ThetaPixDiff' was read in" << endl; } else { MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); if (!sigtheta) { cout << "Object with name 'SigmaTheta' not found" << endl; return; } TH2D *fHSigmaTheta = sigtheta->GetSigmaTheta(); TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta(); TH3D *fHDiffPixTheta = sigtheta->GetDiffPixTheta(); } //------------------------------------ //:::::::::::::::::::::::::::::::::::::::::: 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(fHillas, fHillasSrc); 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 = 30000; char* outName = "matrix_gammas.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("Hillas.fWidth"); matrix.AddColumn("Hillas.fLength"); matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); matrix.AddColumn("abs(Hillas.fAsym)"); matrix.AddColumn("abs(Hillas.fM3Long)"); matrix.AddColumn("abs(Hillas.fM3Trans)"); matrix.AddColumn("abs(HillasSrc.fHeadTail)"); matrix.AddColumn("Hillas.fConc"); matrix.AddColumn("Hillas.fConc1"); matrix.AddColumn("HillasSrc.fDist"); matrix.AddColumn("log10(Hillas.fSize)"); matrix.AddColumn("HillasSrc.fAlpha"); matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); //matrix.AddColumn("MMcEvt.fTheta"); plist.AddToList(pmatrix); MFillH fillmatg(mtxName); fillmatg.SetFilter(&flist); //+++++++++++++++++++++++++++++++++++++++++++++++++++ MSelFinal selfinal(fHillas, fHillasSrc); //***************************** // 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); tlist.AddToList(&hfill1); tlist.AddToList(&hfill2); tlist.AddToList(&hfill3); tlist.AddToList(&hfill4); tlist.AddToList(&flist); tlist.AddToList(&fillmatg); //tlist.AddToList(&selfinal); //***************************** //------------------------------------------- // Execute event loop // MEvtLoop evtloop; evtloop.SetParList(&plist); Int_t maxevents = 1000000000; 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 = 12; 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 section for MC gamma data" << endl; cout << "=========================================================" << endl; } }