| 1 | 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")
|
|---|
| 2 | {
|
|---|
| 3 | double deltat = 17.5;
|
|---|
| 4 | // lvl1 = 2.5;
|
|---|
| 5 | // lvl2 = 0.5;
|
|---|
| 6 |
|
|---|
| 7 | // The sequence file which defines the files for the analysis
|
|---|
| 8 | MSequence seq(seqfile);
|
|---|
| 9 | if (!seq.IsValid())
|
|---|
| 10 | {
|
|---|
| 11 | gLog << "ERROR - Sequence invalid!" << endl;
|
|---|
| 12 | return 1;
|
|---|
| 13 | }
|
|---|
| 14 |
|
|---|
| 15 | // ------------------------------------------------------
|
|---|
| 16 |
|
|---|
| 17 | gLog.Separator("Star");
|
|---|
| 18 | gLog << "Calculate image parameters of sequence ";
|
|---|
| 19 | gLog << seq.GetFileName() << endl;
|
|---|
| 20 | gLog << endl;
|
|---|
| 21 |
|
|---|
| 22 | // ------------------------------------------------------
|
|---|
| 23 |
|
|---|
| 24 | gLog << "Inpath: " << inpath << endl;
|
|---|
| 25 | gLog << "Outpath: " << outpath << endl;
|
|---|
| 26 |
|
|---|
| 27 | const TString rule(Form("s/([0-9]+_[0-9]+)_C[.]root?$/%s\\/$1_I.root/",
|
|---|
| 28 | MJob::Esc(outpath).Data()));
|
|---|
| 29 | gLog << "Rule: " << rule << endl;
|
|---|
| 30 |
|
|---|
| 31 | MDirIter iter;
|
|---|
| 32 | if (seq.GetRuns(iter, MSequence::kFactCal, inpath)<=0)
|
|---|
| 33 | {
|
|---|
| 34 | gLog << "ERROR - Sequence valid but without files." << endl;
|
|---|
| 35 | return 2;
|
|---|
| 36 | }
|
|---|
| 37 |
|
|---|
| 38 | iter.Print();
|
|---|
| 39 |
|
|---|
| 40 | gLog.Separator();
|
|---|
| 41 |
|
|---|
| 42 | // --------------------------------------------------------
|
|---|
| 43 |
|
|---|
| 44 | MBinning bins1( 36, -90, 90, "BinningAlpha");
|
|---|
| 45 | MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
|
|---|
| 46 | MBinning bins4( 25, 0, 2.5, "BinningDist");
|
|---|
| 47 |
|
|---|
| 48 | MBinning binsM1(100, 0, 5, "BinningMuonRadius");
|
|---|
| 49 | MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation");
|
|---|
| 50 | MBinning binsM3(150, 0, 60, "BinningMuonTime");
|
|---|
| 51 | MBinning binsM4(300, 0, 30, "BinningMuonTimeRms");
|
|---|
| 52 | MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi");
|
|---|
| 53 | MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth");
|
|---|
| 54 | MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log");
|
|---|
| 55 | MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma");
|
|---|
| 56 |
|
|---|
| 57 | MBinning binsM9(100, 0, 5, "BinningRadius");
|
|---|
| 58 |
|
|---|
| 59 | // ---------------------------------------------------------
|
|---|
| 60 |
|
|---|
| 61 | // Filter to start muon analysis
|
|---|
| 62 | //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
|
|---|
| 63 | // Filter to calculate further muon parameters
|
|---|
| 64 | // MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&"
|
|---|
| 65 | // "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
|
|---|
| 66 | MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
|
|---|
| 67 | "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
|
|---|
| 68 | "MuonSearchCut");
|
|---|
| 69 |
|
|---|
| 70 | // Filter to fill the MHMuonPar
|
|---|
| 71 | MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
|
|---|
| 72 | "MuonFinalCut");
|
|---|
| 73 | // Filter to write Muons to Muon tree
|
|---|
| 74 | //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
|
|---|
| 75 |
|
|---|
| 76 | // ---------------------------------------------------------
|
|---|
| 77 |
|
|---|
| 78 | MStatusDisplay *d = new MStatusDisplay;
|
|---|
| 79 |
|
|---|
| 80 | MTaskList tlist;
|
|---|
| 81 |
|
|---|
| 82 | MParList plist2;
|
|---|
| 83 | plist2.AddToList(&tlist);
|
|---|
| 84 | plist2.AddToList(&bins1);
|
|---|
| 85 | plist2.AddToList(&bins3);
|
|---|
| 86 | plist2.AddToList(&bins4);
|
|---|
| 87 | plist2.AddToList(&binsM1);
|
|---|
| 88 | plist2.AddToList(&binsM2);
|
|---|
| 89 | plist2.AddToList(&binsM3);
|
|---|
| 90 | plist2.AddToList(&binsM4);
|
|---|
| 91 | plist2.AddToList(&binsM5);
|
|---|
| 92 | plist2.AddToList(&binsM6);
|
|---|
| 93 | plist2.AddToList(&binsM7);
|
|---|
| 94 | plist2.AddToList(&binsM8);
|
|---|
| 95 | plist2.AddToList(&binsM9);
|
|---|
| 96 |
|
|---|
| 97 | MMuonSetup muonsetup;
|
|---|
| 98 | plist2.AddToList(&muonsetup);
|
|---|
| 99 |
|
|---|
| 100 | MEvtLoop loop;
|
|---|
| 101 | loop.SetDisplay(d);
|
|---|
| 102 | loop.SetParList(&plist2);
|
|---|
| 103 |
|
|---|
| 104 | MReadMarsFile read("Events");
|
|---|
| 105 | read.DisableAutoScheme();
|
|---|
| 106 | read.AddFiles(iter);
|
|---|
| 107 |
|
|---|
| 108 | MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
|
|---|
| 109 |
|
|---|
| 110 | MGeomApply apply;
|
|---|
| 111 |
|
|---|
| 112 | MImgCleanStd clean(lvl1, lvl2);
|
|---|
| 113 | clean.SetMethod(MImgCleanStd::kAbsolute);
|
|---|
| 114 | clean.SetTimeLvl1(deltat);//17.5/*1.3*/);
|
|---|
| 115 | clean.SetTimeLvl2(deltat);//17.5/*1.3*/);
|
|---|
| 116 | clean.SetPostCleanType(3);
|
|---|
| 117 |
|
|---|
| 118 | MHillasCalc hcalc;
|
|---|
| 119 | hcalc.Disable(MHillasCalc::kCalcConc);
|
|---|
| 120 |
|
|---|
| 121 | MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
|
|---|
| 122 | hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
|
|---|
| 123 | hrate.SetName("Rate");
|
|---|
| 124 | hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
|
|---|
| 125 | hrate.InitLabels(MH3::kLabelsX);
|
|---|
| 126 | hrate.GetHist().SetMinimum(0);
|
|---|
| 127 |
|
|---|
| 128 | MFillH fillR(&hrate, "", "FillRate");
|
|---|
| 129 |
|
|---|
| 130 | // ------------------ Setup histograms and fill tasks ----------------
|
|---|
| 131 | MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
|
|---|
| 132 | MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
|
|---|
| 133 | MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
|
|---|
| 134 | evtC1.SetErrorSpread(kFALSE);
|
|---|
| 135 | evtC2.SetErrorSpread(kFALSE);
|
|---|
| 136 | evtC2.SetThreshold(0);
|
|---|
| 137 |
|
|---|
| 138 | MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
|
|---|
| 139 | MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
|
|---|
| 140 | MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
|
|---|
| 141 |
|
|---|
| 142 | MFillH fillD1("MHHillas", "MHillas", "FillHillas");
|
|---|
| 143 | MFillH fillD2("MHHillasExt", "", "FillHillasExt");
|
|---|
| 144 | MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
|
|---|
| 145 | MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
|
|---|
| 146 | MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
|
|---|
| 147 | MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
|
|---|
| 148 |
|
|---|
| 149 | MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
|
|---|
| 150 | halpha.SetName("AvsD;Dist;Alpha");
|
|---|
| 151 |
|
|---|
| 152 | MFillH filla(&halpha);
|
|---|
| 153 | filla.SetNameTab("AvsD");
|
|---|
| 154 | filla.SetDrawOption("colz");
|
|---|
| 155 |
|
|---|
| 156 | // ----------------------- Muon Analysis ----------------------
|
|---|
| 157 | //writem.SetFilter(&fmuon4);
|
|---|
| 158 |
|
|---|
| 159 | MHn mhn1;
|
|---|
| 160 | mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
|
|---|
| 161 | mhn1.InitName("MuonRadius");
|
|---|
| 162 | mhn1.InitTitle("MuonRadius");
|
|---|
| 163 |
|
|---|
| 164 | mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
|
|---|
| 165 | mhn1.InitName("MuonDeviation");
|
|---|
| 166 | mhn1.InitTitle("MuonDeviation");
|
|---|
| 167 |
|
|---|
| 168 | mhn1.AddHist("MMuonSearchPar.fTime");
|
|---|
| 169 | mhn1.InitName("MuonTime");
|
|---|
| 170 | mhn1.InitTitle("MuonTime");
|
|---|
| 171 |
|
|---|
| 172 | mhn1.AddHist("MMuonSearchPar.fTimeRms");
|
|---|
| 173 | mhn1.InitName("MuonTimeRms");
|
|---|
| 174 | mhn1.InitTitle("MuonTimeRms");
|
|---|
| 175 |
|
|---|
| 176 | MHn mhn3;
|
|---|
| 177 | mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
|
|---|
| 178 | mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
|
|---|
| 179 | mhn3.InitTitle("MuonRadius");
|
|---|
| 180 | mhn3.SetDrawOption("colz");
|
|---|
| 181 |
|
|---|
| 182 | MHn mhn2;
|
|---|
| 183 | mhn2.AddHist("MMuonCalibPar.fArcPhi");
|
|---|
| 184 | mhn2.InitName("MuonArcPhi");
|
|---|
| 185 | mhn2.InitTitle("MuonArcPhi");
|
|---|
| 186 |
|
|---|
| 187 | mhn2.AddHist("MMuonCalibPar.fArcWidth");
|
|---|
| 188 | mhn2.InitName("MuonArcWidth");
|
|---|
| 189 | mhn2.InitTitle("MuonArcWidth");
|
|---|
| 190 |
|
|---|
| 191 | mhn2.AddHist("MMuonCalibPar.fMuonSize");
|
|---|
| 192 | mhn2.InitName("MuonSize");
|
|---|
| 193 | mhn2.InitTitle("MuonSize");
|
|---|
| 194 |
|
|---|
| 195 | mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
|
|---|
| 196 | mhn2.InitName("MuonRelTimeSigma");
|
|---|
| 197 | mhn2.InitTitle("MuonRelTimeSigma");
|
|---|
| 198 |
|
|---|
| 199 | MHn mhn4;
|
|---|
| 200 | mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
|
|---|
| 201 | mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
|
|---|
| 202 | mhn4.InitTitle("MuonArcPhi");
|
|---|
| 203 | mhn4.SetDrawOption("colz");
|
|---|
| 204 |
|
|---|
| 205 | MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
|
|---|
| 206 | MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
|
|---|
| 207 | MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
|
|---|
| 208 | MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
|
|---|
| 209 | fillmhn1.SetNameTab("MuonS");
|
|---|
| 210 | fillmhn3.SetNameTab("DevVsRad");
|
|---|
| 211 | fillmhn2.SetNameTab("MuonC");
|
|---|
| 212 | fillmhn4.SetNameTab("WidthVsPhi");
|
|---|
| 213 | fillmhn2.SetFilter(&fmuon2);
|
|---|
| 214 | fillmhn4.SetFilter(&fmuon2);
|
|---|
| 215 |
|
|---|
| 216 | MMuonSearchParCalc muscalc;
|
|---|
| 217 | //muscalc.SetFilter(&fmuon1);
|
|---|
| 218 |
|
|---|
| 219 | MMuonCalibParCalc mcalc;
|
|---|
| 220 | mcalc.SetFilter(&fmuon2);
|
|---|
| 221 |
|
|---|
| 222 | MFillH fillmuon("MHSingleMuon", "", "FillMuon");
|
|---|
| 223 | MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
|
|---|
| 224 | fillmuon.SetFilter(&fmuon2);
|
|---|
| 225 | fillmpar.SetFilter(&fmuon3);
|
|---|
| 226 | fillmuon.SetBit(MFillH::kDoNotDisplay);
|
|---|
| 227 |
|
|---|
| 228 | // ---------------------------------------------------------------
|
|---|
| 229 |
|
|---|
| 230 | MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); // EffectiveOnTime
|
|---|
| 231 | write.AddContainer("MTime", "Events");
|
|---|
| 232 | write.AddContainer("MHillas", "Events");
|
|---|
| 233 | write.AddContainer("MHillasExt", "Events");
|
|---|
| 234 | write.AddContainer("MHillasSrc", "Events");
|
|---|
| 235 | write.AddContainer("MImagePar", "Events");
|
|---|
| 236 | write.AddContainer("MNewImagePar", "Events");
|
|---|
| 237 | write.AddContainer("MRawEvtHeader", "Events");
|
|---|
| 238 | write.AddContainer("MRawRunHeader", "RunHeaders");
|
|---|
| 239 | write.AddContainer("MGeomCam", "RunHeaders");
|
|---|
| 240 |
|
|---|
| 241 | MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
|
|---|
| 242 | "MuonHistCut");
|
|---|
| 243 | fillmhn2.SetFilter(&fmuonhn);
|
|---|
| 244 |
|
|---|
| 245 | tlist.AddToList(&read);
|
|---|
| 246 | tlist.AddToList(&cont);
|
|---|
| 247 | tlist.AddToList(&apply);
|
|---|
| 248 | tlist.AddToList(&clean);
|
|---|
| 249 | tlist.AddToList(&hcalc);
|
|---|
| 250 | tlist.AddToList(&fillC1);
|
|---|
| 251 | tlist.AddToList(&fillC2);
|
|---|
| 252 | tlist.AddToList(&fillC3);
|
|---|
| 253 | tlist.AddToList(&fillR);
|
|---|
| 254 | tlist.AddToList(&fillD1);
|
|---|
| 255 | tlist.AddToList(&fillD2);
|
|---|
| 256 | tlist.AddToList(&fillD3);
|
|---|
| 257 | tlist.AddToList(&fillD4);
|
|---|
| 258 | tlist.AddToList(&fillD5);
|
|---|
| 259 | tlist.AddToList(&fillD6);
|
|---|
| 260 | tlist.AddToList(&filla);
|
|---|
| 261 | //tlist.AddToList(&fmuon1);
|
|---|
| 262 | tlist.AddToList(&muscalc);
|
|---|
| 263 | tlist.AddToList(&fillmhn1);
|
|---|
| 264 | tlist.AddToList(&fillmhn3);
|
|---|
| 265 | tlist.AddToList(&fmuon2);
|
|---|
| 266 | tlist.AddToList(&fillmuon);
|
|---|
| 267 | tlist.AddToList(&mcalc);
|
|---|
| 268 | tlist.AddToList(&fmuonhn);
|
|---|
| 269 | tlist.AddToList(&fillmhn2);
|
|---|
| 270 | tlist.AddToList(&fillmhn4);
|
|---|
| 271 | tlist.AddToList(&fmuon3);
|
|---|
| 272 | tlist.AddToList(&fillmpar);
|
|---|
| 273 | //tlist.AddToList(&fmuon4);
|
|---|
| 274 | //tlist.AddToList(&writem);
|
|---|
| 275 | tlist.AddToList(&write);
|
|---|
| 276 |
|
|---|
| 277 | if (!loop.Eventloop())
|
|---|
| 278 | return 3;
|
|---|
| 279 |
|
|---|
| 280 | if (!loop.GetDisplay())
|
|---|
| 281 | return 4;
|
|---|
| 282 |
|
|---|
| 283 | TString title = "-- Image parameters #";
|
|---|
| 284 | title += seq.GetSequence();
|
|---|
| 285 | title += " --";
|
|---|
| 286 | d->SetTitle(title, kFALSE);
|
|---|
| 287 |
|
|---|
| 288 | TString path;
|
|---|
| 289 | path += Form("%s/20%06d_%03d-images.root", outpath,
|
|---|
| 290 | seq.GetSequence()/1000, seq.GetSequence()%1000);
|
|---|
| 291 |
|
|---|
| 292 | d->SaveAs(path);
|
|---|
| 293 |
|
|---|
| 294 | return 0;
|
|---|
| 295 | }
|
|---|
