Index: trunk/Mars/fact/analysis/mc/ganymed.C =================================================================== --- trunk/Mars/fact/analysis/mc/ganymed.C (revision 17733) +++ trunk/Mars/fact/analysis/mc/ganymed.C (revision 17733) @@ -0,0 +1,525 @@ +void SetupWriter(MWriteRootFile *write, const char *name, Bool_t ismc) const +{ + if (!write) + return; + + write->SetName(name); + write->AddContainer("MHillas", "Events"); + write->AddContainer("MHillasSrc", "Events"); + write->AddContainer("MHillasExt", "Events"); + write->AddContainer("MPointingPos", "Events"); + write->AddContainer("MHillasSrcAnti", "Events", kFALSE); + write->AddContainer("MImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar2", "Events", kFALSE); + write->AddContainer("Hadronness", "Events", kFALSE); + write->AddContainer("MSrcPosCam", "Events", kFALSE); + write->AddContainer("MSrcPosAnti", "Events", kFALSE); + write->AddContainer("ThetaSquared", "Events", kFALSE); + write->AddContainer("OpticalAxis", "Events", kFALSE); + write->AddContainer("Disp", "Events", kFALSE); + write->AddContainer("Ghostbuster", "Events", kFALSE); + write->AddContainer("MEnergyEst", "Events", kFALSE); + write->AddContainer("MTime", "Events", kFALSE); + write->AddContainer("MMcEvt", "Events", kFALSE); + write->AddContainer("MEnergyEst", "Events", kFALSE); + write->AddContainer("DataType", "Events"); + write->AddContainer("Weight", "Events"); + write->AddContainer("FileId", "Events"); + write->AddContainer("EvtNumber", "Events"); + + write->AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + write->AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write->AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); + + //if (ismc) + // write->AddCopySource("OriginalMC"); +} + +int ganymed(const char *dataset, const char *outfile, Float_t ra, Float_t dec, Bool_t ismc=false) +{ + MDirIter files; + if (files.ReadList(dataset)<=0) + return 2; + files.Sort(); + + // ============================================================= + + // Crab (05h34'32" 22d00'52") + MPointingPos source("MSourcePos"); + //source.SetSkyPosition(5.5755, 22.0144); // ra[h], dec[deg] + //source.SetSkyPosition(12.356089, 30.176897); // ra[h], dec[deg]//1218 + //source.SetSkyPosition(16.897867,39.760201);// mrk501 + //source.SetSkyPosition(11.0742, 38.2089); // ra[h], dec[deg]//mrk421 + + if (ismc) + source.SetLocalPosition(ra, dec); // RA is ZD, DEC, is AZ!!!! + else + source.SetSkyPosition(ra, dec); // ra[h], dec[deg]//mrk421 + + UInt_t fNumOffSourcePos = 5; + + // ============================================================= + + // FIXME: If fPathIn read cuts and energy estimator from file! + MContinue contq("MImagePar.fNumIslands>3.5" + "|| MNewImagePar.fNumUsedPixels<5.5" + "|| MNewImagePar.fLeakage1>0.1", + "CutQ"); + + MContinue cont0("MHillas.fSize<0" + "|| MHillasExt.fSlopeSpreadWeighted>0.68" + "|| MHillasExt.fSlopeSpreadWeighted<0.18" + "|| log10(MHillas.GetArea*MGeomCam.fConvMm2Deg^2)>(log10(MHillas.fSize)-2)*1.1-1.55" + "|| MNewImagePar.fConcCore < 0.13" + "|| MNewImagePar.fConcCOG < 0.15" + , "Cut0"); + + // Cuts before plots + MContinue cont2("0" + "|| MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg<-0.07" + "|| log10(MNewImagePar.fConc1)<-log10(MHillas.fSize)*0.57+0.13" + //"|| log10(MNewImagePar.fConc1)>-log10(MHillas.fSize)*0.57+0.5" + //"|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.42+0.29" + "|| MHillas.GetArea*MGeomCam.fConvMm2Deg^2<0.01" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.5+0.095" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg<0.4" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>1.2", + "Cut2"); + + TArrayD param(12); + // Parametrization of Disp + param[0] = 1.15136; // constant + param[8] = 0.0681437; // slope + param[9] = 2.62932; // leak + param[10] = 1.51279; // size-offset + param[11] = 0.0507821; // size-slope + // Parametrization for sign of disp (m3long, slope) + param[5] = -0.07; + param[6] = 7.2; + param[7] = 0.5; + // ThetaSq-Cut + param[1] = 0.11; // 0.215 + // Area-Cut + param[2] = 0.215468; + param[3] = 5.63973; + param[4] = 0.0836169; + + // -------------------------------------------------------------------------------- + + // Setup fitter and histograms + MAlphaFitter fit; + fit.SetBackgroundFitMin(0.09); + fit.SetBackgroundFitMax(0.25); + fit.SetPolynomOrder(1); + fit.SetSignalFunction(MAlphaFitter::kThetaSq); + fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale); + + // -------------------------------------------------------------------------------- + + MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); + MBinning binsD( 25, 0, 2.5, "BinningDist"); + MBinning binsC( 25, 0, 2.5, "BinningDistC"); + MBinning binsR( 60, 0, 300, "BinningRate"); + + // -------------------------------------------------------------------------------- + + MParList plist; + plist.AddToList(&source); + plist.AddToList(&fit); + + // Initialize default binnings + plist.AddToList(&bins3); + plist.AddToList(&binsD); + plist.AddToList(&binsC); + plist.AddToList(&binsR); + + MParameterI par("DataType"); + MParameterD weight("Weight"); + plist.AddToList(&par); + plist.AddToList(&weight); + + // Setup Tasklist + MTaskList tlist; + plist.AddToList(&tlist); + + // La Palma Magic1 + MObservatory obs; + plist.AddToList(&obs); + + // -------------------------------------------------------------------------------- + + MParameterD scale; + scale.SetVal(1./fNumOffSourcePos); + + MHThetaSq *halphaoff = plist.FindCreateObj("MHThetaSq", "HistOff"); + halphaoff->ForceUsingSize(); + + MFillH falpha(halphaoff, "", "FillHist"); + MFillH fenergy1("HistEOff [MHThetaSq]", "", "FillHistE"); + + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kNoCut); + cuts.SetVariables(param); + + // Cuts before thetasq + MContinue cont1(&cuts, "Cut1"); + + // Energy estimation + MEnergyEstimate energyest; + energyest.SetRule("MHillas.fSize^0.815*pow(10, 0.0199*(MHillasSrc.fDist*1.17e-2 - (-6.95))^2*(MHillasSrc.fDist*1.17e-2>(-6.95)) + 0.5997 *(MHillasSrc.fDist*1.17e-2 - 0.7568)^2 * (MHillasSrc.fDist*1.17e-2>0.7568) )"); + + // Cuts before writing + MContinue cont3("", "Cut3"); + contq.SetAllowEmpty(); + cont0.SetAllowEmpty(); + cont1.SetAllowEmpty(); + cont2.SetAllowEmpty(); + cont3.SetAllowEmpty(); + + cont1.SetInverted(); + + // Filter for VsSize + MFDataPhrase ftheta("ThetaSquared.fVal < ThetaSquaredCut.fVal", "CutT"); + + // ------------- Loop Off Data -------------------- + MRead *readoffdata = NULL; + if (ismc) + { + MReadMarsFile *read = new MReadMarsFile("Events"); + read->DisableAutoScheme(); + readoffdata = read; + } + else + { + MReadReports *read = new MReadReports; + read->AddTree("Events", "MTime.", MReadReports::kMaster); + read->AddTree("Drive",MReadReports::kRequired); + read->AddTree("Rates",MReadReports::kRequired); + readoffdata = read; + } + + readoffdata->AddFiles(files); + + TString fname0 = Form("%s-summary.root", outfile); + TString fname1 = Form("%s-analysis.root", outfile); + + MWriteRootFile dummy0(fname0.Data(), "RECREATE"); + MWriteRootFile dummy1(fname1.Data(), "RECREATE"); + + MWriteRootFile *write0 = &dummy0; + MWriteRootFile *write1 = &dummy1; + SetupWriter(write0, "WriteAfterCut0", ismc); + SetupWriter(write1, "WriteAfterCut3", ismc); + + MWriteRootFile write2(fname1.Data(), "RECRATE"); + write2.AddContainer("MReportRates", "Rates"); + write2.AddContainer("MTimeRates", "Rates"); + + MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber"); + setevtnum.SetNameParameter("EvtNumber"); + + MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId"); + setrunnum.SetNameParameter("FileId"); + + MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill3c("MHVsSizeOffTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta"); + MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + //MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime"); + fill1a.SetNameTab("PreCut"); + fill2a.SetNameTab("PostCut"); + fill3a.SetNameTab("VsSize"); + fill3c.SetNameTab("CutT"); + fill4a.SetNameTab("HilExt"); + fill5a.SetNameTab("HilSrc"); + fill6a.SetNameTab("ImgPar"); + fill7a.SetNameTab("NewPar"); + //fill9a.SetNameTab("EffOffT"); + + fill3c.SetFilter(&ftheta); + + //MFDataMember fbin("Bin.fVal", '>', 0); + //fill9a.SetFilter(&fbin); + + MTaskList tlist2; + tlist2.SetNumPasses(fNumOffSourcePos); + fill2a.SetWeight(&scale); + fill3a.SetWeight(&scale); + fill3c.SetWeight(&scale); + fill4a.SetWeight(&scale); + fill5a.SetWeight(&scale); + fill6a.SetWeight(&scale); + fill7a.SetWeight(&scale); + + // How to get source position from off- and on-data? + MSrcPosCalc scalc; + scalc.SetMode(MSrcPosCalc::kWobble); + scalc.SetCallback(&tlist2); + + MSrcPosCorrect scor; + + MHillasCalc hcalc; + MHillasCalc hcalc2("MHillasCalcAnti"); + hcalc.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetNameHillasSrc("MHillasSrcAnti"); + hcalc2.SetNameSrcPosCam("MSrcPosAnti"); + + MSrcPosRndm srcrndm; + + tlist2.AddToList(&scalc); + //tlist2.AddToList(&scor); + tlist2.AddToList(&srcrndm); + tlist2.AddToList(&hcalc); + tlist2.AddToList(&hcalc2); + tlist2.AddToList(&setrunnum); + tlist2.AddToList(&setevtnum); + tlist2.AddToList(&cuts); + tlist2.AddToList(&energyest); + tlist2.AddToList(write0); + tlist2.AddToList(&cont0); // Post/Pre shown cuts + tlist2.AddToList(&cont1); // Calculate ThetaSq + tlist2.AddToList(&fill2a); + tlist2.AddToList(&ftheta); // Calculate theta Cut + if (!ismc) + { + tlist2.AddToList(&falpha); + tlist2.AddToList(&fenergy1); + } + tlist2.AddToList(&fill3c); // With theta cut + tlist2.AddToList(&cont2); + tlist2.AddToList(&fill3a); + tlist2.AddToList(&fill4a); + tlist2.AddToList(&fill5a); + tlist2.AddToList(&fill6a); + tlist2.AddToList(&fill7a); + tlist2.AddToList(&cont3); + tlist2.AddToList(write1); + + MPointingPosCalc pcalc; + MPointingDevCalc devcalc; + + tlist.AddToList(readoffdata); + if (!ismc) + tlist.AddToList(&pcalc, "Drive"); + tlist.AddToList(&contq, "Events"); + tlist.AddToList(&fill1a, "Events"); + tlist.AddToList(&tlist2, "Events"); + + // by setting it here it is distributed to all consecutive tasks + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + par.SetVal(0); + weight.SetVal(-1./fNumOffSourcePos); + + MStatusDisplay *d = new MStatusDisplay; + + // Create and setup the eventloop + MEvtLoop evtloop("Off"); + evtloop.SetParList(&plist); + evtloop.SetDisplay(d); + + // Execute first analysis + if (!evtloop.Eventloop()) + return 2; + + if (!evtloop.GetDisplay()) + return 2; + + // ------------- Loop On Data -------------------- + MRead *readondata = NULL; + if (ismc) + { + MReadMarsFile *read = new MReadMarsFile("Events"); + read->DisableAutoScheme(); + readondata = read; + } + else + { + MReadReports *read = new MReadReports; + read->AddTree("Events", "MTime.", MReadReports::kMaster); + read->AddTree("Drive", MReadReports::kRequired); + read->AddTree("Rates", MReadReports::kRequired); + read->AddTree("Humidity", MReadReports::kRequired); + read->AddTree("Temperatures", MReadReports::kRequired); + readondata = read; + } + readondata->AddFiles(files); + + scalc.SetMode(MSrcPosCalc::kDefault); + scalc.SetNumRandomOffPositions(fNumOffSourcePos); + + MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta"); + MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime"); + fill1b.SetNameTab("PreCut"); + fill2b.SetNameTab("PostCut"); + fill3b.SetNameTab("VsSize"); + fill3d.SetNameTab("CutT"); + fill4b.SetNameTab("HilExt"); + fill5b.SetNameTab("HilSrc"); + fill6b.SetNameTab("ImgPar"); + fill7b.SetNameTab("NewPar"); + //fill9b.SetNameTab("EffOnT"); + fill1b.SetDrawOption("same"); + fill2b.SetDrawOption("same"); + fill3b.SetDrawOption("same"); + fill3d.SetDrawOption("same"); + fill4b.SetDrawOption("same"); + fill5b.SetDrawOption("same"); + fill6b.SetDrawOption("same"); + fill7b.SetDrawOption("same"); + //fill9b.SetFilter(&fbin); + + fill3d.SetFilter(&ftheta); + + MHThetaSq *halphaon = plist.FindCreateObj("MHThetaSq", "Hist"); + halphaon->ForceUsingSize(); + + MFillH falpha2(halphaon, "", "FillHist"); + MFillH fenergy2("HistE [MHThetaSq]", "", "FillHistE"); + + tlist2.SetNumPasses(); + + tlist.Replace(readondata); + + MHSrcPosCam hsrcpos(kTRUE); + + MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam"); + fillsrc.SetNameTab("SrcPos"); + + tlist2.AddToListBefore(&fillsrc, &hcalc); + + MH3 hvs1("MPointingPos.fZd"); + MH3 hvs2("MPointingPos.fZd"); + MHVsTime hvstM("MTimeDrive.GetMjd"); + MHVsTime hvst0("MPointingPos.fZd"); + MHVsTime hvst1("MReportRates.fTriggerRate"); + MH3 hvst2("MPointingPos.fZd", "MReportRates.fTriggerRate"); + MHVsTime hvst("MReportTemperatures.GetTempMean"); + MHVsTime hvsh("MReportHumidity.GetMean"); + hvs1.SetName("ObsTime;Theta"); + hvs2.SetName("OnTime;Theta"); + hvstM.SetName("MJD"); + hvst0.SetName("ZdTime"); + hvst1.SetName("Rate"); + hvst2.SetName("RateZd;Theta;Rate"); + hvst.SetName("CamTemp"); + hvsh.SetName("CamHum"); + hvs1.SetTitle("Observation time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]"); + hvs2.SetTitle("On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]"); + hvstM.SetTitle("MJD vs time;Time;MJD"); + hvst0.SetTitle("Zenith distance vs time;Time;Zd [\\circ]"); + hvst1.SetTitle("Camera trigger rate vs time;Time;Rate [Hz]"); + hvst2.SetTitle("Camera trigger rate vs Zd; Zd [\\circ];Rate [Hz]"); + hvst.SetTitle("Average sensor compartment temperature;Time;T [\\circ C]"); + hvsh.SetTitle("Average camera humidity;Time;Hum [%]"); + + MFillH fillvs1(&hvs1, "", "FillObsTime"); + MFillH fillvs2(&hvs2, "", "FillOnTime"); + MFillH fillvstM(&hvstM, "MTime", "FillMjdVsTime"); + MFillH fillvst0(&hvst0, "MTimeDrive", "FillZdVsTime"); + MFillH fillvst1(&hvst1, "MTimeRates", "FillRate"); + MFillH fillvst2(&hvst2, "", "FillRateVsZd"); + MFillH fillvst(&hvst, "MTimeTemperatures", "FillCamTemp"); + MFillH fillvsh(&hvsh, "MTimeHumidity", "FillCamHum"); + fillvs1.SetNameTab("ObsTime"); + fillvs2.SetNameTab("OnTime"); + fillvst1.SetNameTab("Rate"); + fillvst2.SetNameTab("RateZd"); + fillvst.SetNameTab("CamTemp"); + fillvsh.SetNameTab("CamHum"); + + hvs1.SetWeight("MReportRates.fElapsedTime"); + hvs2.SetWeight("MReportRates.fElapsedOnTime"); + + hvst1.SetMinimum(0); + hvst2.SetDrawOption("colz"); + + MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("RateVsRun"); + hrate.SetTitle("Event rate after quality cuts;File Id;Event Rate [Hz];"); + hrate.InitLabels(MH3::kLabelsX); + hrate.GetHist().SetMinimum(0); + + MFillH fillr(&hrate, "", "FillRateVsRun"); + tlist.AddToListAfter(&fillr, &contq, "Events"); + + tlist.Replace(&fill1b); + + tlist2.Replace(&fill2b); + tlist2.Replace(&fill3b); + tlist2.Replace(&fill3d); + tlist2.Replace(&fill4b); + tlist2.Replace(&fill5b); + tlist2.Replace(&fill6b); + tlist2.Replace(&fill7b); + tlist2.Replace(&falpha2); + tlist2.Replace(&fenergy2); + if (!ismc) + { + tlist.AddToList(&fillvs1, "Rates"); + tlist.AddToList(&fillvs2, "Rates"); + tlist.AddToList(&fillvst1, "Rates"); + tlist.AddToList(&fillvst2, "Rates"); + tlist.AddToList(&write2, "Rates"); + tlist.AddToList(&fillvst, "Temperatures"); + tlist.AddToList(&fillvsh, "Humidity"); + tlist.AddToList(&fillvst0, "Drive"); + tlist.AddToList(&fillvstM, "Drive"); + } + + // by setting it here it is distributed to all consecutive tasks + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + par.SetVal(1); + weight.SetVal(1); + + // Execute first analysis + if (!evtloop.Eventloop()) + return 2; + + if (!evtloop.GetDisplay()) + return 2; + + delete readondata; + delete readoffdata; + + TObjArray cont; + //cont.Add(const_cast(GetEnv())); + //cont.Add(const_cast(&fSequence)); + + //TNamed cmdline("CommandLine", fCommandLine.Data()); + //cont.Add(&cmdline); + + if (d) + { + TString title = "-- Analysis #"; + title += dataset; + title += " --"; + d->SetTitle(title, kFALSE); + + cont.Add(d); + TString path; + path += fname1; + + d->Write(); + } + +// if (!WriteContainer(cont, fname1, "RECREATE")) +// return; + + return 0; +}