Index: branches/MarsISDCBranchBasedOn17887/fact/analysis/mc/star_NEWCLEANING.C =================================================================== --- branches/MarsISDCBranchBasedOn17887/fact/analysis/mc/star_NEWCLEANING.C (revision 18000) +++ branches/MarsISDCBranchBasedOn17887/fact/analysis/mc/star_NEWCLEANING.C (revision 18000) @@ -0,0 +1,420 @@ +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.h" + +#include "../mcore/DrsCalib.h" +#include "MDrsCalibration.h" +//#include "MLogManip.h" +#include "MExtralgoSpline.h" +#include "MSequence.h" +#include "MStatusArray.h" +#include "MHCamera.h" +#include "MJob.h" +#include "MWriteRootFile.h" +#include "MHCamera.h" +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" +#include "MDirIter.h" +#include "MTaskList.h" +#include "MFDataPhrase.h" +#include "MArrayF.h" +#include "MBadPixelsTreat.h" +#include "MCalibrateDrsTimes.h" +#include "MHSectorVsTime.h" +#include "MHCamEvent.h" +#include "MExtractTimeAndChargeSpline.h" +#include "MFillH.h" +#include "MDrsCalibApply.h" +#include "MGeomApply.h" +#include "MContinue.h" +#include "MRawFitsRead.h" +#include "MEvtLoop.h" +#include "MParList.h" +#include "MStatusDisplay.h" +#include "MDrsCalibrationTime.h" +#include "MH3.h" +#include "MGeomCamFACT.h" +#include "MCalibrateFact.h" +#include "MParameters.h" +#include "MWriteAsciiFile.h" +#include "MFMagicCuts.h" + +#include "MMuonSetup.h" +#include "MReadMarsFile.h" +#include "MHillasCalc.h" +#include "MHn.h" +#include "MMuonSearchParCalc.h" +#include "MMuonCalibParCalc.h" +#include "MBinning.h" +#include "MImgCleanTime.h" + + +#endif + +using namespace std; + +int star(const TString callistofile, TString outfile, + TString displayfile, TString displaytitle, + const TString txtfile, + Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5); + + +int star(const TString infile = "callisto.root", TString outfile = "", + Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5) +{ + return star(infile, outfile, "", "", "", lvl1, lvl2, deltat); +} + +int star(const TString callistofile, TString outfile, + TString displayfile, TString displaytitle, + const TString txtfile, + Double_t lvl1, Double_t lvl2, Double_t deltat) +{ + + //Double_t lvl1=7.8, + //Double_t lvl2=3.9, + // lvl1 = 2.5; + // lvl2 = 0.5; + // ------------------------------------------------------ + + if (displaytitle.IsNull()) + displaytitle = gSystem->BaseName(callistofile); + + FileStat_t fstat; + int rc = gSystem->GetPathInfo(outfile, fstat); + bool isdir = !rc || R_ISDIR(fstat.fMode); + + char *buf = gSystem->ConcatFileName(outfile, "star.root"); + outfile = buf; + delete [] buf; + + if (displayfile.IsNull()) + { + displayfile = outfile; + displayfile.Insert(displayfile.Last('.'), "-display"); + } + else + { + if (isdir && gSystem->DirName(displayfile)==TString(".")) + { + buf = gSystem->ConcatFileName(outfile, displayfile); + displayfile = buf; + delete [] buf; + } + } + + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << endl; + gLog << "Callisto file: " << callistofile << endl; + gLog << "Output file: " << outfile << endl; + gLog << "Display file: " << displayfile << endl; + gLog << "Display title: " << displaytitle << endl; + gLog << endl; + gLog << "lvl1: " << lvl1 <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.AddFile(callistofile); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanTime clean; + clean.SetMinCount(0); + clean.SetMinSize(25); + clean.SetDeltaT(1.25*17.5*0.1111); + + 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); + + // --------------------------------------------------------------- + + MWriteAsciiFile writeascii(txtfile, "MHillas"); + writeascii.AddColumns("MHillasExt"); + writeascii.AddColumns("MHillasSrc"); + + + MWriteRootFile write(outfile, "RECREATE", "Image parameters", 2); + write.AddContainer("MTime", "Events", kFALSE); + 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("ThetaSquared", "Events"); + write.AddContainer("Disp", "Events"); + write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); + write.AddContainer("MMcEvt", "Events", kFALSE); + write.AddContainer("IncidentAngle", "Events", kFALSE); + write.AddContainer("MPointingPos", "Events", kFALSE); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + + MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0", + "MuonHistCut"); + fillmhn2.SetFilter(&fmuonhn); + + 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; + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kNoCut); + cuts.SetVariables(param); + + + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&apply); + tlist.AddToList(&clean); + tlist.AddToList(&hcalc); + tlist.AddToList(&cuts); + 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 (!txtfile.IsNull()) + tlist.AddToList(&writeascii); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + d->SetTitle(displaytitle, kFALSE); + d->SaveAs(displayfile); + + return 0; +} +