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