Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C (revision 6045) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C (revision 6045) @@ -0,0 +1,328 @@ +/* ======================================================================== *\ + ! + ! * + ! * 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): Thomas Bretz 5/2002 + ! Abelardo Moralejo 1/2004 + ! Eva Domingo 1/2005 + ! + ! Copyright: MAGIC Software Development, 2000-2004 + ! + ! + \* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// STARMC - STandard Analysis and Reconstruction (MC example) +// +// This macro is a version of the standard converter to convert raw data +// into image parameters, made to show how to run analysis on MC files. +// +// (1/2005):: Spline extractor added and additional containers written out +// +///////////////////////////////////////////////////////////////////////////// + +#include "MImgCleanStd.h" + +void makeHillasMC() +{ + // + // This is a demonstration program which calculates the image + // parameters from Magic Monte Carlo files (output of camera). + + + TString* CalibrationFilename; + TString* AnalysisFilename; + TString* OutFilename1; + TString* OutFilename2; + + // ------------- USER CHANGE ----------------- + // Comment line starting "CalibrationFileName" to disable calibration. + // In that case the units of the MHillas.fSize parameter will be ADC counts + // (rather, equivalent ADC counts in inner pixels, since we correct for the + // possible differences in gain of outer pixels) + + // File to be used in the calibration (must be a camera file without added noise) + CalibrationFilename = new TString("/discpepe/root_0.73mirror/wuerzburg/gammas_nonoise/Gamma_zbin0_0_7_1000to1009_w0.root"); + // File to be analyzed + AnalysisFilename = new TString("/discpepe/root_0.73mirror/wuerzburg/gammas/Gamma_zbin*.root"); + + // Change output file names as desired. + // If you want only one output (not dividing the events in Train and Test samples), + // comment the initialization of OutFilename2. + OutFilename1 = new TString("/mnt/users/domingo/MAGIC/Disp05/All_gammas_0.73mirror_zbin0to12_cleaned_4035_Spline.root"); + // OutFilename1 = new TString("star_train.root"); // Output file name 1 (test) + // OutFilename2 = new TString("star_test.root"); // Output file name 2 (train) + + // Fraction of events (taken at random) which one wants to process from the + // file to be analyzed (useful to make smaller files if starting sample is too large). + Float_t accepted_fraction = 1.; + + + // ------------- USER CHANGE ----------------- + // Cleaning levels + Float_t CleanLev[2] = {4.0, 3.5}; // Tail cuts for image analysis + + // Signal extractor + // (default=Spline, other extraction methods can be used) + const Int_t wsize=2; + MExtractor* sigextract = new MExtractTimeAndChargeSpline(); + ((MExtractTimeAndChargeSpline*)sigextract)->SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum); + ((MExtractTimeAndChargeSpline*)sigextract)->SetChargeType(MExtractTimeAndChargeSpline::kIntegral); + ((MExtractTimeAndChargeSpline*)sigextract)->SetRiseTime((Float_t)wsize*0.25); + ((MExtractTimeAndChargeSpline*)sigextract)->SetFallTime((Float_t)wsize*0.75); + // Define FADC slices to be integrated in high and low gain: + sigextract->SetRange(1, 11, 2, 12); + // Defines when to switch to low gain + sigextract->SetSaturationLimit(240); + + + // ------------------------------------------- + // + // 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); + + // MGeomCamMagic geomcam; + + // MSrcPosCam src; + // WOBBLE MODE!! + // src.SetX(120.); // units: mm + // src.SetXY(0.,0.4/geomcam.GetConvMm2Deg()); + // src.SetXY(0.,0.); + // src.SetReadyToSave(); + + MBadPixelsCam badpix; + + // plist.AddToList(&geomcam); + // plist.AddToList(&src); + plist.AddToList(&badpix); + + + // Now setup the tasks and tasklist: + // --------------------------------- + // + MReadMarsFile read("Events"); + if (CalibrationFilename) + read.AddFile(CalibrationFilename->Data()); + read.DisableAutoScheme(); + + MGeomApply geom; + // Reads in geometry from MC file and sets the right sizes for + // several parameter containers. + + MMcPedestalCopy pcopy; + // Copies pedestal data from the MC file run fadc header to the MPedestalCam container. + + MPointingPosCalc pointcalc; + // Creates MPointingPos object and fill it with the telescope orientation + // information taken from MMcEvt. + + MMcCalibrationUpdate mccalibupdate; + + MCalibrateData calib; + // Transforms signals from ADC counts into photons. In the first + // loop it applies a "dummy" calibration supplied by MMcCalibrationUpdate, just + // to equalize inner and outer pixels. At the end of the first loop, in the + // PostProcess of MMcCalibrationCalc (see below) the true calibration constants + // are calculated. + calib.SetCalibrationMode(MCalibrateData::kFfactor); + + // MBlindPixelCalc blind; + // blind.SetUseInterpolation(); + + MMcCalibrationCalc mccalibcalc; + // Calculates calibration constants to convert from ADC counts to photons. + + MImgCleanStd clean(CleanLev[0], CleanLev[1]); + // Applies tail cuts to image. Since the calibration is performed on + // noiseless camera files, the precise values of the cleaning levels + // are unimportant for the calibration file processing (in any case, + // only pixels without any C-photon will be rejected). + clean.SetCleanRings(1); + // clean.SetRemoveSingle(kFALSE); + // clean.SetMethod(MImgCleanStd::kDemocratic); + + MHillasCalc hcalc; + // Calculates Hillas parameters. + + tlist.AddToList(&read); + tlist.AddToList(&geom); + tlist.AddToList(&pcopy); + tlist.AddToList(&pointcalc); + tlist.AddToList(sigextract); + tlist.AddToList(&mccalibupdate); + tlist.AddToList(&calib); + tlist.AddToList(&clean); + // tlist.AddToList(&blind); + tlist.AddToList(&hcalc); + tlist.AddToList(&mccalibcalc); + + + // + // Open output files: + // + + MWriteRootFile write1(OutFilename1->Data()); // Writes output1 + write1.AddContainer("MRawRunHeader", "RunHeaders"); + write1.AddContainer("MMcRunHeader", "RunHeaders"); + write1.AddContainer("MSrcPosCam", "RunHeaders"); + write1.AddContainer("MGeomCam", "RunHeaders"); + write1.AddContainer("MMcConfigRunHeader", "RunHeaders"); + write1.AddContainer("MMcCorsikaRunHeader", "RunHeaders"); + write1.AddContainer("MMcFadcHeader", "RunHeaders"); + write1.AddContainer("MMcTrigHeader", "RunHeaders"); + + write1.AddContainer("MRawEvtHeader", "Events"); + write1.AddContainer("MMcEvt", "Events"); + write1.AddContainer("MHillas", "Events"); + write1.AddContainer("MHillasExt", "Events"); + write1.AddContainer("MHillasSrc", "Events"); + write1.AddContainer("MImagePar", "Events"); + write1.AddContainer("MNewImagePar", "Events"); + write1.AddContainer("MConcentration", "Events"); + write1.AddContainer("MIslands", "Events"); + write1.AddContainer("MTopology", "Events"); + write1.AddContainer("MPointingPos", "Events"); + //write1.AddContainer("MArrivalTimeCam", "Events"); + //write1.AddContainer("MCerPhotEvt", "Events"); + + if (OutFilename2) + { + MWriteRootFile write2(OutFilename2->Data()); // Writes output2 + write2.AddContainer("MRawRunHeader", "RunHeaders"); + write2.AddContainer("MMcRunHeader", "RunHeaders"); + write2.AddContainer("MSrcPosCam", "RunHeaders"); + write2.AddContainer("MGeomCam", "RunHeaders"); + write2.AddContainer("MMcConfigRunHeader", "RunHeaders"); + write2.AddContainer("MMcCorsikaRunHeader", "RunHeaders"); + write2.AddContainer("MMcFadcHeader", "RunHeaders"); + write2.AddContainer("MMcTrigHeader", "RunHeaders"); + + write2.AddContainer("MRawEvtHeader", "Events"); + write2.AddContainer("MMcEvt", "Events"); + write2.AddContainer("MHillas", "Events"); + write2.AddContainer("MHillasExt", "Events"); + write2.AddContainer("MHillasSrc", "Events"); + write2.AddContainer("MImagePar", "Events"); + write2.AddContainer("MNewImagePar", "Events"); + write2.AddContainer("MConcentration", "Events"); + write2.AddContainer("MIslands", "Events"); + write2.AddContainer("MTopology", "Events"); + write2.AddContainer("MPointingPos", "Events"); + //write2.AddContainer("MArrivalTimeCam", "Events"); + //write2.AddContainer("MCerPhotEvt", "Events"); + + + // Divide output in train and test samples, using the event number + // (odd/even) to achieve otherwise unbiased event samples: + MF filter1("{MMcEvt.fEvtNumber%2}>0.5"); + MF filter2("{MMcEvt.fEvtNumber%2}<0.5"); + + write1.SetFilter(&filter1); + write2.SetFilter(&filter2); + } + + + // + // FIRST LOOP: Calibration loop (with no noise file) + // + MProgressBar bar; + bar.SetWindowName("Calibrating..."); + + MEvtLoop evtloop; + evtloop.SetProgressBar(&bar); + evtloop.SetParList(&plist); + + if (CalibrationFilename) + { + if (!evtloop.Eventloop()) + return; + mccalibcalc.GetHistADC2PhotEl()->Write(); + mccalibcalc.GetHistPhot2PhotEl()->Write(); + // Writes out the histograms used for calibration. + } + + + // + // SECOND LOOP: Analysis loop (process file with noise) + // + bar.SetWindowName("Analyzing..."); + + MIslands isl; + MArrivalTimeCam timecam; + MTopology topology; + plist.AddToList(&isl); + plist.AddToList(&timecam); + plist.AddToList(&topology); + + MArrivalTimeCalc2 timecalc; + // Calculates the arrival time as the mean time of the fWindowSize + // time slices which have the highest integral content. + // MArrivalTimeCalc timecalc; + // Calculates the arrival times of photons. Returns the absolute + // maximum of the spline that interpolates the FADC slices. + + MIslandsCalc islcalc; + islcalc.SetOutputName("MIslands"); + islcalc.SetAlgorithm(1); + // MIslandsClean islclean(40); + // islclean.SetInputName("MIslands"); + // islclean.SetMethod(1); + + MTopologyCalc topcalc; + + // Change the read task by another one which reads the file we want to analyze: + MReadMarsFile read2("Events"); + read2.AddFile(AnalysisFilename->Data()); + read2.DisableAutoScheme(); + tlist.AddToListBefore(&read2, &read); + tlist.RemoveFromList(&read); + + // Add new tasks (Islands and Topology calculations) + tlist.AddToListBefore(&timecalc,&mccalibupdate); + tlist.AddToListBefore(&islcalc,&hcalc); + // tlist.AddToListBefore(&islclean,&hcalc); + tlist.AddToListBefore(&topcalc,&hcalc); + + // Analyze only the desired fraction of events, skip the rest: + MFEventSelector eventselector; + Float_t rejected_fraction = 1. - accepted_fraction; + eventselector.SetSelectionRatio(rejected_fraction); + MContinue skip(&eventselector); + tlist.AddToListBefore(&skip, sigextract); + + // Remove calibration task from list. + tlist.RemoveFromList(&mccalibcalc); + + // Add tasks to write output: + if (OutFilename2) + { + tlist.AddToList(&filter1); + tlist.AddToList(&filter2); + tlist.AddToList(&write2); + } + tlist.AddToList(&write1); + + + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); +}