Index: trunk/Mars/fact/analysis/star.C =================================================================== --- trunk/Mars/fact/analysis/star.C (revision 15231) +++ trunk/Mars/fact/analysis/star.C (revision 15231) @@ -0,0 +1,295 @@ +int star(const char *seqfile="sequences/20111205_013.seq", Double_t lvl1=7.8, Double_t lvl2=3.9, const char *inpath = "output", const char *outpath = "output") +{ + double deltat = 17.5; + // lvl1 = 2.5; + // lvl2 = 0.5; + + // The sequence file which defines the files for the analysis + MSequence seq(seqfile); + if (!seq.IsValid()) + { + gLog << "ERROR - Sequence invalid!" << endl; + return 1; + } + + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << "Calculate image parameters of sequence "; + gLog << seq.GetFileName() << endl; + gLog << endl; + + // ------------------------------------------------------ + + gLog << "Inpath: " << inpath << endl; + gLog << "Outpath: " << outpath << endl; + + const TString rule(Form("s/([0-9]+_[0-9]+)_C[.]root?$/%s\\/$1_I.root/", + MJob::Esc(outpath).Data())); + gLog << "Rule: " << rule << endl; + + MDirIter iter; + if (seq.GetRuns(iter, MSequence::kFactCal, inpath)<=0) + { + gLog << "ERROR - Sequence valid but without files." << endl; + return 2; + } + + iter.Print(); + + gLog.Separator(); + + // -------------------------------------------------------- + + MBinning bins1( 36, -90, 90, "BinningAlpha"); + MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); + MBinning bins4( 25, 0, 2.5, "BinningDist"); + + MBinning binsM1(100, 0, 5, "BinningMuonRadius"); + MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation"); + MBinning binsM3(150, 0, 60, "BinningMuonTime"); + MBinning binsM4(300, 0, 30, "BinningMuonTimeRms"); + MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi"); + MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth"); + MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log"); + MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma"); + + MBinning binsM9(100, 0, 5, "BinningRadius"); + + // --------------------------------------------------------- + + // Filter to start muon analysis + //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut"); + // Filter to calculate further muon parameters +// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&" +// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut"); + MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&" + "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4", + "MuonSearchCut"); + + // Filter to fill the MHMuonPar + MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2", + "MuonFinalCut"); + // Filter to write Muons to Muon tree + //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut"); + + // --------------------------------------------------------- + + MStatusDisplay *d = new MStatusDisplay; + + MTaskList tlist; + + MParList plist2; + plist2.AddToList(&tlist); + plist2.AddToList(&bins1); + plist2.AddToList(&bins3); + plist2.AddToList(&bins4); + plist2.AddToList(&binsM1); + plist2.AddToList(&binsM2); + plist2.AddToList(&binsM3); + plist2.AddToList(&binsM4); + plist2.AddToList(&binsM5); + plist2.AddToList(&binsM6); + plist2.AddToList(&binsM7); + plist2.AddToList(&binsM8); + plist2.AddToList(&binsM9); + + MMuonSetup muonsetup; + plist2.AddToList(&muonsetup); + + MEvtLoop loop; + loop.SetDisplay(d); + loop.SetParList(&plist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + read.AddFiles(iter); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanStd clean(lvl1, lvl2); + clean.SetMethod(MImgCleanStd::kAbsolute); + clean.SetTimeLvl1(deltat);//17.5/*1.3*/); + clean.SetTimeLvl2(deltat);//17.5/*1.3*/); + clean.SetPostCleanType(3); + + MHillasCalc hcalc; + hcalc.Disable(MHillasCalc::kCalcConc); + + MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("Rate"); + hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];"); + hrate.InitLabels(MH3::kLabelsX); + hrate.GetHist().SetMinimum(0); + + MFillH fillR(&hrate, "", "FillRate"); + + // ------------------ Setup histograms and fill tasks ---------------- + MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]"); + MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction"); + MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]"); + evtC1.SetErrorSpread(kFALSE); + evtC2.SetErrorSpread(kFALSE); + evtC2.SetThreshold(0); + + MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam"); + MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels"); + MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos"); + + MFillH fillD1("MHHillas", "MHillas", "FillHillas"); + MFillH fillD2("MHHillasExt", "", "FillHillasExt"); + MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc"); + MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar"); + MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar"); + MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize"); + + MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)"); + halpha.SetName("AvsD;Dist;Alpha"); + + MFillH filla(&halpha); + filla.SetNameTab("AvsD"); + filla.SetDrawOption("colz"); + + // ----------------------- Muon Analysis ---------------------- + //writem.SetFilter(&fmuon4); + + MHn mhn1; + mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonRadius"); + mhn1.InitTitle("MuonRadius"); + + mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonDeviation"); + mhn1.InitTitle("MuonDeviation"); + + mhn1.AddHist("MMuonSearchPar.fTime"); + mhn1.InitName("MuonTime"); + mhn1.InitTitle("MuonTime"); + + mhn1.AddHist("MMuonSearchPar.fTimeRms"); + mhn1.InitName("MuonTimeRms"); + mhn1.InitTitle("MuonTimeRms"); + + MHn mhn3; + mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation"); + mhn3.InitTitle("MuonRadius"); + mhn3.SetDrawOption("colz"); + + MHn mhn2; + mhn2.AddHist("MMuonCalibPar.fArcPhi"); + mhn2.InitName("MuonArcPhi"); + mhn2.InitTitle("MuonArcPhi"); + + mhn2.AddHist("MMuonCalibPar.fArcWidth"); + mhn2.InitName("MuonArcWidth"); + mhn2.InitTitle("MuonArcWidth"); + + mhn2.AddHist("MMuonCalibPar.fMuonSize"); + mhn2.InitName("MuonSize"); + mhn2.InitTitle("MuonSize"); + + mhn2.AddHist("MMuonCalibPar.fRelTimeSigma"); + mhn2.InitName("MuonRelTimeSigma"); + mhn2.InitTitle("MuonRelTimeSigma"); + + MHn mhn4; + mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth"); + mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth"); + mhn4.InitTitle("MuonArcPhi"); + mhn4.SetDrawOption("colz"); + + MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar"); + MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar"); + MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius"); + MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi"); + fillmhn1.SetNameTab("MuonS"); + fillmhn3.SetNameTab("DevVsRad"); + fillmhn2.SetNameTab("MuonC"); + fillmhn4.SetNameTab("WidthVsPhi"); + fillmhn2.SetFilter(&fmuon2); + fillmhn4.SetFilter(&fmuon2); + + 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); + + // --------------------------------------------------------------- + + MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); // EffectiveOnTime + write.AddContainer("MTime", "Events"); + 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("MRawRunHeader", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); + + MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0", + "MuonHistCut"); + fillmhn2.SetFilter(&fmuonhn); + + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&apply); + tlist.AddToList(&clean); + tlist.AddToList(&hcalc); + tlist.AddToList(&fillC1); + tlist.AddToList(&fillC2); + tlist.AddToList(&fillC3); + tlist.AddToList(&fillR); + tlist.AddToList(&fillD1); + tlist.AddToList(&fillD2); + tlist.AddToList(&fillD3); + tlist.AddToList(&fillD4); + tlist.AddToList(&fillD5); + tlist.AddToList(&fillD6); + tlist.AddToList(&filla); + //tlist.AddToList(&fmuon1); + tlist.AddToList(&muscalc); + tlist.AddToList(&fillmhn1); + tlist.AddToList(&fillmhn3); + tlist.AddToList(&fmuon2); + tlist.AddToList(&fillmuon); + tlist.AddToList(&mcalc); + tlist.AddToList(&fmuonhn); + tlist.AddToList(&fillmhn2); + tlist.AddToList(&fillmhn4); + tlist.AddToList(&fmuon3); + tlist.AddToList(&fillmpar); + //tlist.AddToList(&fmuon4); + //tlist.AddToList(&writem); + tlist.AddToList(&write); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + TString title = "-- Image parameters #"; + title += seq.GetSequence(); + title += " --"; + d->SetTitle(title, kFALSE); + + TString path; + path += Form("%s/20%06d_%03d-images.root", outpath, + seq.GetSequence()/1000, seq.GetSequence()%1000); + + d->SaveAs(path); + + return 0; +}