Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2907) +++ trunk/MagicSoft/Mars/Changelog (revision 2908) @@ -6,4 +6,8 @@ 2004/01/24: Abelardo Moralejo + + * macros/starmcstereo.C + - Added. Example of how to run the analysis chain for MC files + containing simulation of stereo systems of 2 telescopes. * mcalib/MCalibrate.cc Index: trunk/MagicSoft/Mars/macros/starmcstereo.C =================================================================== --- trunk/MagicSoft/Mars/macros/starmcstereo.C (revision 2908) +++ trunk/MagicSoft/Mars/macros/starmcstereo.C (revision 2908) @@ -0,0 +1,284 @@ +/* ======================================================================== *\ +! +! * +! * 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(); +} Index: trunk/MagicSoft/Mars/mcalib/MCalibrate.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 2907) +++ trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 2908) @@ -144,7 +144,13 @@ else { - signal = sig.GetExtractedSignalHiGain(); + if (sig.GetExtractedSignalHiGain() > 9999.) + { + signal = 0.; + signalErr = 0.; + } + else + signal = sig.GetExtractedSignalHiGain(); } - + // Float_t calibrationConversionFactor = pix.GetMeanConversionFFactorMethod(); const Float_t calibrationConversionFactor = pix.GetMeanConversionBlindPixelMethod(); @@ -162,5 +168,5 @@ MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixid, nphot, nphotErr); - + if (sig.GetNumLoGainSaturated() > 0) cpix->SetPixelSaturated();