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