/* ======================================================================== *\ ! ! * ! * 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 // ///////////////////////////////////////////////////////////////////////////// // // User change. // Float_t ctx[7] = {0., 0., 0., 0., 0., 0., 0.}; Float_t cty[7] = {-70., -40., -30., 30., 50., 60., 70.}; // in meters // // FIXME: unfortunately present version of reflector was not prepared for // stereo configurations and keeps no track of CT position. So the positions // must be set above by the user, making sure that they correspond to the // files one is analysing. // void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 5) { if (ct1 > sizeof(ctx)/sizeof(ctx[0]) || ct2 > sizeof(ctx)/sizeof(ctx[0]) ) { cout << endl << "Wrong CT id number!" << endl; return; } Int_t CT[2] = {ct1, ct2}; // Only 2-telescope analysis for the moment Int_t NCTs = sizeof(CT)/sizeof(CT[0]); // ------------- user change ----------------- Char_t* AnalysisFilename = "gam-yXX-00001.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] = {5, 9}; // First and last FADC bin of the range to be integrated, Int_t BinsLow[2] = {5, 9}; // 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); MSrcPosCam src[NCTs]; MBadPixelsCam badpix[NCTs]; for (Int_t ict = 0; ict < NCTs; ict++) { TString s = "MSrcPosCam;"; s += CT[ict]; src[ict].SetName(s); src[ict].SetReadyToSave(); plist.AddToList(&(src[ict])); TString b = "MBadPixelsCam;"; b += CT[ict]; badpix[ict].SetName(b); badpix[ict].SetReadyToSave(); plist.AddToList(&(badpix[ict])); } // // 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. // stereocalc.SetCT1coor(ctx[CT[0]-1],cty[CT[0]-1]); stereocalc.SetCT2coor(ctx[CT[1]-1],cty[CT[1]-1]); tlist.AddToList(&stereocalc); MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5"); MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5"); // // ^^^^ Filters to divide output in two: test and train samples. // 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(); }