/* ======================================================================== *\ ! ! * ! * 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/2009 ! ! Copyright: MAGIC Software Development, 2000-2009 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MJSimulation // // // Force reading a corsika file even if the footer (RUNE-section) is missing // by setting fForceMode to kTRUE or from the resource file by // // ForceMode: Yes // // // To switch off the simulation of the camera electronics, use: // // Camera: Off // // Note, that the border of the camera and the propertied of the cones // are still simulated (simply because it is fast). Furthermore, this // switches off the trigger for the output, i.e. all data which deposits // at least one photon in at least one pixel is written to the output file. // // // In case of a pedestal or calibration run the artificial trigger can // be "switched off" and the cosmics trrigger "switched on" by setting // fForceTrigger to kTRUE or from the resource file by // // ForceTrigger: Yes // // ///////////////////////////////////////////////////////////////////////////// #include "MJSimulation.h" #include #include // Core #include "MLog.h" #include "MLogManip.h" #include "MArgs.h" //#include "MDirIter.h" #include "MParList.h" #include "MTaskList.h" #include "MEvtLoop.h" #include "MStatusDisplay.h" // Tasks #include "MCorsikaRead.h" #include "MContinue.h" #include "MFillH.h" #include "MGeomApply.h" #include "MHillasCalc.h" #include "MImgCleanStd.h" #include "MWriteRootFile.h" #include "MSimMMCS.h" #include "MSimAbsorption.h" #include "MSimAtmosphere.h" #include "MSimReflector.h" #include "MSimPointingPos.h" #include "MSimPSF.h" #include "MSimGeomCam.h" #include "MSimSignalCam.h" #include "MSimAPD.h" #include "MSimExcessNoise.h" #include "MSimCamera.h" #include "MSimTrigger.h" #include "MSimReadout.h" #include "MSimRandomPhotons.h" #include "MSimBundlePhotons.h" #include "MSimCalibrationSignal.h" // Histograms #include "MBinning.h" #include "MHn.h" #include "MHCamera.h" #include "MHCamEvent.h" #include "MHPhotonEvent.h" // Container #include "MRawRunHeader.h" #include "MParameters.h" #include "MReflector.h" #include "MParEnv.h" #include "MPulseShape.h" #include "MGeomCam.h" #include "MPedestalCam.h" #include "MPedestalPix.h" ClassImp(MJSimulation); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // // Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE // MJSimulation::MJSimulation(const char *name, const char *title) : fForceMode(kFALSE), fCamera(kTRUE) ,fForceTrigger(kFALSE) { fName = name ? name : "MJSimulation"; fTitle = title ? title : "Standard analysis and reconstruction"; } Bool_t MJSimulation::CheckEnvLocal() { fForceMode = GetEnv("ForceMode", fForceMode); fForceTrigger = GetEnv("ForceTrigger", fForceTrigger); fCamera = GetEnv("Camera", fCamera); return kTRUE; } Bool_t MJSimulation::WriteResult(const MParList &plist) { if (fPathOut.IsNull()) { *fLog << inf << "No output path specified via SetPathOut - no output written." << endl; return kTRUE; } TObjArray cont; cont.Add(const_cast(GetEnv())); //cont.Add(const_cast(&fSequence)); cont.Add(plist.FindObject("MPulseShape")); cont.Add(plist.FindObject("Reflector")); cont.Add(plist.FindObject("MGeomCam")); cont.Add(plist.FindObject("GeomCones")); if (fDisplay) { // TString title = "-- Reflector: "; // title += fSequence.GetSequence(); // title += " --"; // fDisplay->SetTitle(title, kFALSE); fDisplay->SetTitle("Ceres", kFALSE); cont.Add(fDisplay); } // const TString oname = Form("reflector%08d.root", fSequence.GetSequence()); const TString oname = "ceres.root"; return WriteContainer(cont, oname, "RECREATE"); } void MJSimulation::SetupHist(MHn &hist) const { hist.AddHist("MCorsikaEvtHeader.fTotalEnergy"); hist.InitName("Energy"); hist.InitTitle("Energy;E [GeV]"); hist.SetLog(kTRUE, kTRUE, kFALSE); hist.AddHist("MPhotonEvent.GetNumExternal"); hist.InitName("Size"); hist.InitTitle("Size;S [#]"); hist.SetLog(kTRUE, kTRUE, kFALSE); hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100"); hist.SetDrawOption("colz"); hist.InitName("Impact;Impact;Impact"); hist.InitTitle("Impact;West <--> East [m];South <--> North [m]"); hist.SetAutoRange(); hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000"); hist.InitName("Height"); hist.InitTitle("FirstInteractionHeight;h [km]"); hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()"); hist.InitName("SkyOrigin;Az;Zd"); hist.InitTitle("Sky Origin;Az [deg];Zd [deg]"); hist.SetDrawOption("colz"); hist.SetAutoRange(); TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())"; TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())"; TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())"; TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()"; hist.AddHist(form); hist.InitName("ViewCone"); hist.InitTitle("Incident Angle;\\alpha [\\deg]"); } void MJSimulation::SetupCommonFileStructure(MWriteRootFile &write) const { // Common run headers write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MGeomCam", "RunHeaders"); write.AddContainer("MMcRunHeader", "RunHeaders"); // Common events write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); write.AddContainer("MRawEvtHeader", "Events"); write.AddContainer("MMcEvt", "Events"); } Bool_t MJSimulation::Process(const MArgs &args) { /* if (!fSequence.IsValid()) { *fLog << err << "ERROR - Sequence invalid!" << endl; return kFALSE; } */ //if (!HasWritePermission(GetPathOut())) // return kFALSE; *fLog << inf; fLog->Separator(GetDescriptor()); if (!CheckEnv()) return kFALSE; *fLog << warn << "FIXME: Monte Carlo simulation: Sequences not supported yet."; //*fLog << fSequence.GetFileName() << endl; *fLog << endl; // -------------------------------------------------------------------------------- //MDirIter iter; //if (fSequence.GetRuns(iter, MSequence::kCalibrated)<=0) // return kFALSE; // Setup Parlist MParList plist; plist.AddToList(this); // take care of fDisplay! // setup TaskList MTaskList tasks; plist.AddToList(&tasks); // -------------------------------------------------------------------------------- // FIXME: Allow empty GeomCones! MParEnv env0("Reflector"); MParEnv env1("GeomCones"); MParEnv env2("MGeomCam"); env0.SetClassName("MReflector"); env1.SetClassName("MGeomCam"); env2.SetClassName("MGeomCam"); plist.AddToList(&env0); plist.AddToList(&env1); plist.AddToList(&env2); plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm"); MPulseShape shape; plist.AddToList(&shape); // *** FIXME *** FIXME *** FIXME *** plist.FindCreateObj("MMcRunHeader"); MRawRunHeader header; header.SetValidMagicNumber(); //header.InitFadcType(3); header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTData); if (args.GetNumArguments()==1) { if (!args.GetArgumentStr(0).CompareTo("pedestal", TString::kIgnoreCase)) { header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPedestal); header.SetSourceInfo("Pedestal"); } if (!args.GetArgumentStr(0).CompareTo("calibration", TString::kIgnoreCase)) { header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTCalibration); header.SetSourceInfo("Calibration"); } if (!args.GetArgumentStr(0).CompareTo("pointrun", TString::kIgnoreCase)) header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPointRun); } // FIXME: Move to MSimPointingPos, MSimCalibrationSignal // Can we use this as input for MSimPointingPos? header.SetObservation("On", "MonteCarlo"); plist.AddToList(&header); // ++++++++ FIXME FIXME FIXME +++++++++++++ /* MPedestalCam pedcam; pedcam.Init(geomcam.GetNumPixels()); for (UInt_t i=0; i(env2.GetCont()); if (binstr.IsDefault()) binstr.SetEdgesLin(150, -shape.GetPulseWidth(), header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetPulseWidth()); if (binsd.IsDefault() && cam) binsd.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg()); if (binsd0.IsDefault() && cam) binsd0.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg()); header.Print(); // FIXME: Display acceptance curves! if (fDisplay) { TCanvas &c = fDisplay->AddTab("Info"); c.Divide(2,2); c.cd(1); gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); gPad->SetGridx(); gPad->SetGridy(); gROOT->SetSelectedPad(0); shape.DrawClone()->SetBit(kCanDelete); if (env0.GetCont() && (header.IsDataRun() || header.IsPointRun())) { c.cd(3); gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); gPad->SetGridx(); gPad->SetGridy(); gROOT->SetSelectedPad(0); static_cast(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete); } if (fCamera) { if (env1.GetCont()) { c.cd(2); gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); gPad->SetGridx(); gPad->SetGridy(); gROOT->SetSelectedPad(0); MHCamera *c = new MHCamera(static_cast(*env1.GetCont())); c->SetStats(kFALSE); c->SetBit(MHCamera::kNoLegend); c->SetBit(kCanDelete); c->Draw(); } if (cam) { c.cd(4); gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); gPad->SetGridx(); gPad->SetGridy(); gROOT->SetSelectedPad(0); MHCamera *c = new MHCamera(*cam); c->SetStats(kFALSE); c->SetBit(MHCamera::kNoLegend); c->SetBit(kCanDelete); c->Draw(); } } } //------------------------------------------- // Execute first analysis if (!evtloop.Eventloop(fMaxEvents)) { gLog << err << GetDescriptor() << ": Failed." << endl; return kFALSE; } //------------------------------------------- // FIXME: Display Spline in tab // FIXME: Display Reflector in tab // FIXME: Display Routing(?) in tab // FIXME: Display Camera(s) in tab //------------------------------------------- if (!WriteResult(plist)) return kFALSE; *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;; return kTRUE; }