/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Abelardo Moralejo 1/2004 ! Thomas Bretz, 5/2002 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // STARMCSTEREO - STandard Analysis and Reconstruction (for MC stereo files) // // This macro is the standard converter to convert raw data from stereo // camera simulation into image parameters // ///////////////////////////////////////////////////////////////////////////// void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 3) { Int_t CT[2] = {ct1, ct2}; Int_t NCTs = sizeof(CT)/sizeof(CT[0]); // ------------- user change ----------------- Char_t* AnalysisFilename = "G_XX_*_w0.root"; // File to be analyzed Char_t* OutFileTag = "gammas"; // Output file tag Float_t CleanLev[2] = {4., 3.}; // Tail cuts for image analysis Int_t BinsHigh[2] = {0, 5}; // First and last FADC bin of the range to be integrated, Int_t BinsLow[2] = {0, 5}; // for high and low gain respectively. // ------------------------------------------- // // This is a demonstration program which calculates the image // parameters from a Magic Monte Carlo root file. Units of Size are // for the moment, FADC counts. // // // Create a empty Parameter List and an empty Task List // The tasklist is identified in the eventloop by its name // MParList plist; MTaskList tlist; plist.AddToList(&tlist); TString s = " MSrcPosCam;"; s += CT[0]; MSrcPosCam src1(s); src1.SetReadyToSave(); plist.AddToList(&src1); if (NCTs==2) { s = " MSrcPosCam;"; s += CT[1]; MSrcPosCam src2(s); src2.SetReadyToSave(); plist.AddToList(&src2); } // // Now setup the tasks and tasklist: // --------------------------------- // MReadMarsFile read("Events"); read.DisableAutoScheme(); read.AddFile(AnalysisFilename); MGeomApply* apply = new MGeomApply[NCTs]; MMcPedestalCopy* pcopy = new MMcPedestalCopy[NCTs]; MExtractSignal* sigextract = new MExtractSignal[NCTs]; MMcCalibrationUpdate* mccalibupdate = new MMcCalibrationUpdate[NCTs]; MCalibrate* calib = new MCalibrate[NCTs]; MImgCleanStd** clean = new MImgCleanStd*[NCTs]; MHillasCalc* hcalc = new MHillasCalc[NCTs]; MHillasSrcCalc* scalc = new MHillasSrcCalc[NCTs]; TString outfile = "star_"; outfile += CT[0]; if (NCTs==2) { outfile += "_"; outfile += CT[1]; } // // We have two output files (will be later train and test sampls for random forest) // outfile += "_"; outfile += OutFileTag; outfile += "_train.root"; MWriteRootFile write1(outfile); outfile = "star_"; outfile += CT[0]; if (NCTs==2) { outfile += "_"; outfile += CT[1]; } outfile += "_"; outfile += OutFileTag; outfile += "_test.root"; MWriteRootFile write2(outfile); for (Int_t i = 0; i < NCTs; i++) { apply[i]->SetSerialNumber(CT[i]); pcopy[i]->SetSerialNumber(CT[i]); sigextract[i]->SetSerialNumber(CT[i]); sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]); mccalibupdate[i]->SetSerialNumber(CT[i]); calib[i]->SetSerialNumber(CT[i]); clean[i] = new MImgCleanStd(CleanLev[0], CleanLev[1]); clean[i]->SetSerialNumber(CT[i]); hcalc[i]->SetSerialNumber(CT[i]); scalc[i]->SetSerialNumber(CT[i]); write1.SetSerialNumber(CT[i]); write2.SetSerialNumber(CT[i]); write1.AddContainer(write1.AddSerialNumber("MMcEvt"), "Events"); write1.AddContainer(write1.AddSerialNumber("MHillas"), "Events"); write1.AddContainer(write1.AddSerialNumber("MHillasExt"), "Events"); write1.AddContainer(write1.AddSerialNumber("MHillasSrc"), "Events"); write1.AddContainer(write1.AddSerialNumber("MNewImagePar"), "Events"); write1.AddContainer(write1.AddSerialNumber("MSrcPosCam"), "RunHeaders"); write2.AddContainer(write2.AddSerialNumber("MMcEvt"), "Events"); write2.AddContainer(write2.AddSerialNumber("MHillas"), "Events"); write2.AddContainer(write2.AddSerialNumber("MHillasExt"), "Events"); write2.AddContainer(write2.AddSerialNumber("MHillasSrc"), "Events"); write2.AddContainer(write2.AddSerialNumber("MNewImagePar"), "Events"); write2.AddContainer(write2.AddSerialNumber("MSrcPosCam"), "RunHeaders"); } if (NCTs==2) { write1.AddContainer("MStereoPar", "Events"); write2.AddContainer("MStereoPar", "Events"); } write1.AddContainer("MRawRunHeader", "RunHeaders"); write1.AddContainer("MMcRunHeader", "RunHeaders"); write2.AddContainer("MRawRunHeader", "RunHeaders"); write2.AddContainer("MMcRunHeader", "RunHeaders"); tlist.AddToList(&read); for (i = 0; i < NCTs; i++) { tlist.AddToList(&(apply[i])); tlist.AddToList(&(pcopy[i])); tlist.AddToList(&(sigextract[i])); tlist.AddToList(&(mccalibupdate[i])); tlist.AddToList(&(calib[i])); tlist.AddToList(clean[i]); tlist.AddToList(&(hcalc[i])); tlist.AddToList(&(scalc[i])); } MStereoCalc stereocalc; stereocalc.SetCTids(CT[0],CT[1]); // // FIXME: telescope coordinates must be introduced by the user, since // they are not available nor in the camera file, neither in the reflector // output. // if (CT[0] == 1) stereocalc.SetCT1coor(0.,0.); else if (CT[0] == 2) stereocalc.SetCT1coor(0.,0.); else if (CT[0] == 3) stereocalc.SetCT1coor(60.,60.); // in meters else if (CT[0] == 4) stereocalc.SetCT1coor(60.,-60.); else if (CT[0] == 5) stereocalc.SetCT1coor(-60.,60.); else if (CT[0] == 6) stereocalc.SetCT1coor(-60.,-60.); else { cout << "Unknown CT id!" << endl; exit; } if (NCTs==2) { if (CT[1] == 1) stereocalc.SetCT2coor(0.,0.); else if (CT[1] == 2) stereocalc.SetCT2coor(0.,0.); else if (CT[1] == 3) stereocalc.SetCT2coor(60.,60.); // in meters else if (CT[1] == 4) stereocalc.SetCT2coor(60.,-60.); else if (CT[1] == 5) stereocalc.SetCT2coor(-60.,60.); else if (CT[1] == 6) stereocalc.SetCT2coor(-60.,-60.); else { cout << "Unknown CT id!" << endl; exit; } tlist.AddToList(&stereocalc); } MF filter1("MMcEvt;1.fEvtNumber<2501"); MF filter2("MMcEvt;1.fEvtNumber>2500"); // // ^^^^ Filters to divide output in two: test and train samples. // FIXME: It is better to separate the events in odd and even // event numbers (it is independent of the number of events in // the file), but it is not yet possible because we cannot use // the modulus operator(%) in the filter yet. // write1.SetFilter (&filter1); write2.SetFilter (&filter2); tlist.AddToList(&filter1); tlist.AddToList(&write1); tlist.AddToList(&filter2); tlist.AddToList(&write2); // // Create and set up the eventloop // MProgressBar bar; MEvtLoop evtloop; evtloop.SetProgressBar(&bar); evtloop.SetParList(&plist); // // Execute your analysis // if (!evtloop.Eventloop()) return; for (Int_t i= 0; i < NCTs; i++ ) delete clean[i]; tlist.PrintStatistics(); }