/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz, 1/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MJStar // // Resource file entries are case sensitive! // ///////////////////////////////////////////////////////////////////////////// #include "MJStar.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MDirIter.h" #include "MParList.h" #include "MTaskList.h" #include "MEvtLoop.h" #include "MStatusDisplay.h" #include "MH3.h" #include "MHVsTime.h" #include "MHSectorVsTime.h" #include "MHCamEvent.h" #include "MHCamEventRot.h" #include "MBinning.h" #include "MReadReports.h" #include "MReadMarsFile.h" #include "MF.h" #include "MFDataMember.h" #include "MFDeltaT.h" #include "MFSoftwareTrigger.h" #include "MContinue.h" #include "MGeomApply.h" #include "MEventRateCalc.h" #include "MImgCleanStd.h" #include "MSrcPosCalc.h" #include "MSrcPosCorrect.h" #include "MHillasCalc.h" #include "MMuonSearchParCalc.h" #include "MMuonCalibParCalc.h" #include "MFillH.h" #include "MWriteRootFile.h" #include "MMuonSetup.h" #include "MObservatory.h" #include "MPointingPosCalc.h" ClassImp(MJStar); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // // Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE // MJStar::MJStar(const char *name, const char *title) : fMuonAnalysis(kTRUE) { fName = name ? name : "MJStar"; fTitle = title ? title : "Standard analysis and reconstruction"; } Bool_t MJStar::CheckEnvLocal() { DisableMuonAnalysis(!GetEnv("MuonAnalysis", fMuonAnalysis)); return kTRUE; } Bool_t MJStar::WriteResult() { if (fPathOut.IsNull()) { *fLog << inf << "No output path specified via SetPathOut - no output written." << endl; return kTRUE; } const TString oname = Form("%s/star%08d.root", (const char*)fPathOut, fSequence.GetSequence()); *fLog << inf << "Writing to file: " << oname << endl; TFile file(oname, "RECREATE"); *fLog << inf << " - MStatusDisplay..." << flush; if (fDisplay && fDisplay->Write()<=0) { *fLog << err << "Unable to write MStatusDisplay to " << oname << endl; return kFALSE; } *fLog << inf << "ok." << endl; return kTRUE; } Bool_t MJStar::Process(Bool_t ismc) { if (!fSequence.IsValid()) { *fLog << err << "ERROR - Sequence invalid!" << endl; return kFALSE; } //if (!CheckEnv()) // return kFALSE; CheckEnv(); // -------------------------------------------------------------------------------- *fLog << inf; fLog->Separator(GetDescriptor()); *fLog << "Calculate image parameters of sequence "; *fLog << fSequence.GetName() << endl; *fLog << endl; // -------------------------------------------------------------------------------- MDirIter iter; const Int_t n0 = fSequence.SetupDatRuns(iter, MSequence::kCalibrated, fPathData); const Int_t n1 = fSequence.GetNumDatRuns(); if (n0==0) { *fLog << err << "ERROR - No input files of sequence found!" << endl; return kFALSE; } if (n0!=n1) { *fLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl; if (fLog->GetDebugLevel()>4) { *fLog << dbg << "Files which are searched:" << endl; iter.Print(); } return kFALSE; } // Setup Parlist MParList plist; plist.AddToList(this); // take care of fDisplay! MObservatory obs; plist.AddToList(&obs); MMuonSetup muonsetup; plist.AddToList(&muonsetup); // Setup binnings for muon analysis MBinning bins1("BinningRadius"); MBinning bins2("BinningArcWidth"); MBinning bins3("BinningRingBroadening"); MBinning bins4("BinningSizeVsArcRadius"); MBinning bins5("BinningMuonWidth"); MBinning bins6("BinningArcPhi"); plist.AddToList(&bins1); plist.AddToList(&bins2); plist.AddToList(&bins3); plist.AddToList(&bins4); plist.AddToList(&bins5); plist.AddToList(&bins6); // Setup Tasklist MTaskList tlist; plist.AddToList(&tlist); MReadReports readreal; readreal.AddTree("Events", "MTime.", MReadReports::kMaster); readreal.AddTree("Drive", MReadReports::kRequired); readreal.AddTree("Starguider", MReadReports::kRequired); readreal.AddTree("Currents", MReadReports::kRequired); readreal.AddTree("CC"); readreal.AddFiles(iter); MReadMarsFile readmc("Events"); readmc.DisableAutoScheme(); readmc.AddFiles(iter); // ------------------ Setup general tasks ---------------- MFDeltaT fdeltat; MContinue cont(&fdeltat, "FilterDeltaT", "Filter events with wrong timing"); cont.SetInverted(); MGeomApply apply; // Only necessary to craete geometry MEventRateCalc rate; rate.SetNumEvents(1200); MFSoftwareTrigger swtrig; MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger"); contsw.SetInverted(); MImgCleanStd clean; clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm"); MSrcPosCalc poscalc; MHillasCalc hcalc; hcalc.Disable(MHillasCalc::kCalcConc); // ------------------ Setup histograms and fill tasks ---------------- MHCamEvent evt0a(0, "Cleaned", "Signal after Cleaning;;S [\\gamma]"); MHCamEvent evt0b(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction"); evt0b.SetThreshold(0); MFillH fillvs("MHRate", "MTime", "FillEventRate"); MFillH fillp1("MHPointing", "MTimeDrive", "FillDrive"); MFillH fillp2("MHPointing", "MTimeStarguider", "FillStarguider"); fillp1.SetBit(MFillH::kDoNotDisplay); fillp1.SetBit(MFillH::kCanSkip); fillp2.SetBit(MFillH::kCanSkip); MFillH fill0a(&evt0a, "MSignalCam", "FillSignalCam"); MFillH fill0b(&evt0b, "MSignalCam", "FillCntUsedPixels"); MFillH fill1("MHHillas", "MHillas", "FillHillas"); MFillH fill2("MHHillasExt", "", "FillHillasExt"); MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc"); MFillH fill4("MHImagePar", "MImagePar", "FillImagePar"); MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar"); MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime"); //fillvs.SetNameTab("Rate"); fill9.SetNameTab("EffOnTime"); // ------------------ Setup write task ---------------- // Effective on-time need its own not to be skipped by (eg) image cleaning // Muons needs its own to have a unique SetReadyToSave const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data())); MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW"); MWriteRootFile writet(2, rule, fOverwrite?"RECREATE":"NEW"); // EffectiveOnTime MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW"); // Muons writem.SetName("WriteMuons"); // Data write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); write.AddContainer("MHillasSrc", "Events"); write.AddContainer("MImagePar", "Events"); write.AddContainer("MNewImagePar", "Events"); write.AddContainer("MRawEvtHeader", "Events"); write.AddContainer("MPointingPos", "Events"); // Run Header write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MBadPixelsCam", "RunHeaders"); write.AddContainer("MGeomCam", "RunHeaders"); write.AddContainer("MObservatory", "RunHeaders"); // Muon Setup write.AddContainer("BinningRadius", "RunHeaders"); write.AddContainer("BinningArcWidth", "RunHeaders"); write.AddContainer("BinningRingBroadening", "RunHeaders"); write.AddContainer("BinningSizeVsArcRadius", "RunHeaders"); write.AddContainer("MMuonSetup", "RunHeaders"); if (ismc) { // Monte Carlo Data write.AddContainer("MMcEvt", "Events"); write.AddContainer("MMcTrig", "Events"); // Monte Carlo Run Headers write.AddContainer("MMcRunHeader", "RunHeaders"); write.AddContainer("MMcTrigHeader", "RunHeaders"); write.AddContainer("MMcFadcHeader", "RunHeaders"); write.AddContainer("MMcConfigRunHeader", "RunHeaders"); write.AddContainer("MMcCorsikaRunHeader", "RunHeaders"); } else { write.AddContainer("MTime", "Events"); // Drive write.AddContainer("MReportDrive", "Drive"); write.AddContainer("MTimeDrive", "Drive"); // Starguider write.AddContainer("MReportStarguider", "Starguider", kFALSE); write.AddContainer("MTimeStarguider", "Starguider", kFALSE); // Effective On Time writet.AddContainer("MEffectiveOnTime", "EffectiveOnTime"); writet.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime"); } // What to write in muon tree writem.AddContainer("MMuonSearchPar", "Muons"); writem.AddContainer("MMuonCalibPar", "Muons"); writem.AddContainer("MHillas", "Muons"); writem.AddContainer("MHillasExt", "Muons"); writem.AddContainer("MHillasSrc", "Muons"); writem.AddContainer("MImagePar", "Muons"); writem.AddContainer("MNewImagePar", "Muons"); writem.AddContainer("MRawEvtHeader", "Muons"); writem.AddContainer("MPointingPos", "Muons"); if (ismc) { // Monte Carlo Data writem.AddContainer("MMcEvt", "Muons"); writem.AddContainer("MMcTrig", "Muons"); } if (ismc) if (fMuonAnalysis) writem.AddCopySource("OriginalMC"); else write.AddCopySource("OriginalMC"); MTaskList tlist2("Events"); tlist2.AddToList(&apply); if (!ismc) tlist2.AddToList(&cont); tlist2.AddToList(&contsw); if (!ismc) { // Calibration events don't enter star at all. tlist2.AddToList(&rate); tlist2.AddToList(&fillvs); tlist2.AddToList(&fill9); tlist2.AddToList(&writet); } tlist2.AddToList(&clean); tlist2.AddToList(&fill0a); tlist2.AddToList(&fill0b); tlist2.AddToList(&poscalc); tlist2.AddToList(&hcalc); tlist2.AddToList(&fill1); tlist2.AddToList(&fill2); tlist2.AddToList(&fill3); tlist2.AddToList(&fill4); tlist2.AddToList(&fill5); // ----------------------- Muon Analysis ---------------------- // Filter to start muon analysis MF fmuon1("MHillas.fSize>150", "MuonPreCut"); // Filter to calculate further muon parameters MF fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&" "(MMuonSearchPar.fDeviation<45)", "MuonSearchCut"); // Filter to fill the MHMuonPar MF fmuon3("(MMuonCalibPar.fArcPhi>190) && (MMuonSearchPar.fDeviation<35) &&" "(MMuonCalibPar.fArcWidth<0.20) && (MMuonCalibPar.fArcWidth>0.04)", "MuonFinalCut"); // Filter to write Muons to Muon tree MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut"); writem.SetFilter(&fmuon4); MMuonSearchParCalc muscalc; muscalc.SetFilter(&fmuon1); MMuonCalibParCalc mcalc; mcalc.SetFilter(&fmuon2); MFillH fillmuon("MHSingleMuon", "", "FillMuon"); MFillH fillmpar("MHMuonPar", "", "FillMuonPar"); fillmuon.SetFilter(&fmuon2); fillmpar.SetFilter(&fmuon3); fillmuon.SetBit(MFillH::kDoNotDisplay); if (fMuonAnalysis) { tlist2.AddToList(&fmuon1); tlist2.AddToList(&muscalc); tlist2.AddToList(&fmuon2); tlist2.AddToList(&fillmuon); tlist2.AddToList(&mcalc); tlist2.AddToList(&fmuon3); tlist2.AddToList(&fillmpar); tlist2.AddToList(&fmuon4); tlist2.AddToList(&writem); } // ------------------------------------------------------------ // Initialize histogram MHSectorVsTime histdc, histrms; histdc.SetNameTime("MTimeCurrents"); histdc.SetTitle("Mean of all DC Currents;; [nA]"); histdc.SetMinimum(0); histdc.SetMaximum(10); histrms.SetNameTime("MTimeCurrents"); histrms.SetTitle("Mean pedestal rms of all pixels;;<\\sigma_{p}> [phe]"); histrms.SetType(5); histrms.SetMinimum(0); histrms.SetMaximum(10); /* // Define area index [0=inner, 1=outer] // TArrayI inner(1); inner[0] = 0; histdc.SetAreaIndex(inner); */ // Task to fill the histogram MFillH filldc(&histdc, "MCameraDC", "FillDC"); MFillH fillrms(&histrms, "MPedPhotFromExtractorRndm", "FillPedRms"); //MFillH filltst("MHTest", "MTime", "FillTest"); filldc.SetNameTab("Currents"); fillrms.SetNameTab("MeanRms"); //filltst.SetNameTab("Test"); MFillH fillw("MHWeather", "MTimeCC", "FillWeather"); MPointingPosCalc pcalc; // ------------------------------------------------------------ tlist.AddToList(ismc ? (MTask*)&readmc : (MTask*)&readreal); tlist.AddToList(&pcalc, ismc ? "Events" : "Drive"); //tlist.AddToList(&filltst, "Events"); tlist.AddToList(&tlist2, "Events"); if (!ismc) { tlist.AddToList(&fillw, "CC"); tlist.AddToList(&fillp1, "Drive"); tlist.AddToList(&fillp2, "Starguider"); tlist.AddToList(&filldc, "Currents"); tlist.AddToList(&fillrms, "Currents"); } tlist.AddToList(&write); // Create and setup the eventloop MEvtLoop evtloop(fName); evtloop.SetParList(&plist); evtloop.SetDisplay(fDisplay); evtloop.SetLogStream(fLog); if (!SetupEnv(evtloop)) return kFALSE; // Execute first analysis if (!evtloop.Eventloop(fMaxEvents)) { *fLog << err << GetDescriptor() << ": Failed." << endl; return kFALSE; } if (!WriteResult()) return kFALSE; *fLog << all << GetDescriptor() << ": Done." << endl; *fLog << endl << endl; return kTRUE; }