Index: trunk/MagicSoft/Mars/mtemp/mpisa/Changelog_pisa =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpisa/Changelog_pisa (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mpisa/Changelog_pisa (revision 4055) @@ -0,0 +1,37 @@ +Please do not write behind the end of line tag so that it is possible +to send a Changelog file to a line printer directly. Thanks. + +For your entries please follow this style: + +_yyyy/mm/dd:_[Your Name] +_ +___*_filename: +_____-_description +_ +___*_filename2, filename3: +_____-_description +_ +_ +_ + +While an underscore is a placeholder for a white-space or an empty line. + + -*-*- END OF LINE -*-*- + 2004/05/12: Antonio Stamerra + + * mtemp/macros + - added macros: + - first_ana.C given a file with hillas analysis produces ON-OFF plots + - production.C produces calibrated data (photons) from raw data + - HillasAnal.C makes standard Hillas analysis + + * mtemp/macros + - added this directory + + * added this file Changelog_pisa + - test + + 2004/05/11: Antonio Stamerra + + * mtemp/mpisa + - added directory mpisa in mtemp Index: trunk/MagicSoft/Mars/mtemp/mpisa/macros/HillasAnal.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpisa/macros/HillasAnal.C (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mpisa/macros/HillasAnal.C (revision 4055) @@ -0,0 +1,191 @@ +/* ======================================================================== *\ +! +! * +! * 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, 12/2000 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////// +// +// MagicHillas.C +// ============= +// +// This is a demonstration program which calculates image (Hillas +) +// parameters using as input a Merpp output file (raw data). +// All parameters are written to an output file called hillas.root. Also +// filles histograms are displayed. +// For the calculation an arbitrary signal extractor (MCerPhotAnal2/Calc) +// is used. +// +/////////////////////////////////////////////////////////////////////////// + +void HillasAnal(const char *filename="output_test.root") +{ + // + // 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); + + // + // Use this if you want to change the binning of one of + // the histograms. You can use: + // BinningConc, BinningConc1, BinningAsym, BinningM3Long, + // BinningM3Trans, BinningWidth, BinningLength, BinningDist, + // BinningHeadTail, BinningAlpha, BinningSize, BinningDelta, + // BinningPixels and BinningCamera + // + // For more information see MBinning and the corresponding + // histograms + // + // MBinning binsalpha("BinningAlpha"); + // binsalpha.SetEdges(90, 0, 90); // 90 bins from 0 to 90 deg + // plist.AddToList(&binsalpha); + + // MBinning binssize("BinningSize"); + // binssize.SetEdgesLog(50, 1, 1e7); + // plist.AddToList(&binssize); + + // + // Craete the object which hlods the source positions in the camera + // plain in respect to which the image parameters will be calculated. + // For real data the containers will be filled by a task. + // + MSrcPosCam source; + source.SetReadyToSave(); + plist.AddToList(&source); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + // The first argument is the tree you want to read. + // Events: Cosmic ray events + // PedEvents: Pedestal Events + // CalEvents: Calibration Events + // + MReadMarsFile read("Events", filename); + read.DisableAutoScheme(); + + // Setup a task which makes sure that all used arrays have + // the correct size + MGeomApply geomapl; + + // tasks used if MC files are detected to calculate pedestals + // MMcPedestalCopy pcopy; + // MMcPedestalNSBAdd pnsb; + + // calculate the signal in a very simple way + // for real adat files use MCerPhotAnal2 instead which also + // calculates the pedestal. + // MCerPhotCalc ncalc; + // + // Alternative photon calculation: + // Example: use only 2nd to 6th FADC slices for photon calculation: + // + // const Float_t x[15]={0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + // const TArrayF w(15,(Float_t*)x); + // ncalc.SetWeights(w); + // + + // setup image cleaning and blind pixel treatment + MImgCleanStd clean; + MBlindPixelCalc blind; + + // + // Instead of unmapping the pixels you can also (The use of this + // class is deprecated, it will be replaced by MBadPixels*) + // + // blind.SetUseInterpolation(); + // blind.SetUseCetralPixel(); + // + + // setup tasks to calculate image parameters + MHillasCalc hcalc; + MHillasSrcCalc csrc1; + + // setup tasks to fill histograms + MFillH hfill1("MHHillas", "MHillas"); + MFillH hfill2("MHHillasExt"); + MFillH hfill3("MHStarMap", "MHillas"); + MFillH hfill4("HistExtSource [MHHillasExt]", "MHillasSrc"); + MFillH hfill5("HistSource [MHHillasSrc]", "MHillasSrc"); + //MFillH hfill6("MHNewImagePar"); + + // setup task to write containers to a file + MWriteRootFile write("hillas_Mrk421.root"); + write.AddContainer("MHillas"); + write.AddContainer("MHillasExt"); + write.AddContainer("MHillasSrc"); + write.AddContainer("MRawEvtHeader"); + write.AddContainer("MRawRunHeader"); + //write.AddContainer("HistExtSource"); + //write.AddContainer("MHNewImagePar"); + + // Setup the contents of zour tasklist + tlist.AddToList(&read); + tlist.AddToList(&geomapl); + //tlist.AddToList(&pcopy); + //tlist.AddToList(&pnsb); + //tlist.AddToList(&ncalc); + tlist.AddToList(&clean); + tlist.AddToList(&blind); + + tlist.AddToList(&hcalc); + tlist.AddToList(&csrc1); + + tlist.AddToList(&hfill1); + tlist.AddToList(&hfill2); + tlist.AddToList(&hfill3); + tlist.AddToList(&hfill4); + tlist.AddToList(&hfill5); + //tlist.AddToList(&hfill6); + tlist.AddToList(&write); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute your analysis + // + MProgressBar bar; + evtloop.SetProgressBar(&bar); + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // + // After the analysis is finished we can display the histograms + // + plist.FindObject("MHHillas")->DrawClone(); + plist.FindObject("MHHillasExt")->DrawClone(); + plist.FindObject("MHStarMap")->DrawClone(); + plist.FindObject("HistSource")->DrawClone(); + plist.FindObject("HistExtSource")->DrawClone(); + //plist.FindObject("MHNewImagePar")->DrawClone(); +} + Index: trunk/MagicSoft/Mars/mtemp/mpisa/macros/first_ana.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpisa/macros/first_ana.C (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mpisa/macros/first_ana.C (revision 4055) @@ -0,0 +1,125 @@ +/* ======================================================================== *\ +! +! * +! * 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): Alessio Piccioli, 04/2004 +! Author(s): Antonio Stamerra, 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////// +// +// first_ana.C +// ============= +// +// This macro produces a ON/OFF plot using the alpha HillasSrc branch +// taken from an analysis file produced by the AnalisiHillas.C macro +// The significance of excess events is also computed (simple calculation, +// no Li-Ma) +// +/////////////////////////////////////////////////////////////////////////// + + +Float_t degTomm ( Float_t deg ) { + return deg / 57.29577951 * 17000 ; +} + +Int_t first_ana ( Char_t *fNameOn = "crab.root" , Char_t *fNameOff = "offcrab.root" , Float_t minsize = 1000 , Float_t minWdeg = 0.00 , Float_t maxWdeg = 0.125 , Float_t minLdeg = 0.0 , Float_t maxLdeg = 0.26 , Float_t minDdeg = 0.2 , Float_t maxDdeg = 0.8 ) { + + TFile *fileOn = new TFile ( fNameOn , "READ" ) ; + TFile *fileOff = new TFile ( fNameOff , "READ" ) ; + + TTree *treeOn = (TTree *) fileOn -> Get ( "Parameters" ) ; + TTree *treeOff = (TTree *) fileOff -> Get ( "Parameters" ) ; + + TString title ; + + title = "Abs(Alpha) plot ON-Source" ; + TH1D *absalphaOn = new TH1D ( "absalphaOn" , title , 9 , 0 , 90 ) ; + + title = "Abs(Alpha) plot ON-Source" ; + TH1D *absalphaOff = new TH1D ( "absalphaOff" , title , 9 , 0 , 90 ) ; + + TString scut ; + scut = "MHillas.fSize>" ; + scut += minsize ; + cout << "CUT ON SIZE: " << scut << endl ; + TCut Scut = scut ; + + scut = "" ; + scut += degTomm(minWdeg) ; + scut += " < MHillas.fWidth && MHillas.fWidth < " ; + scut += degTomm(maxWdeg) ; + cout << "CUT ON WIDTH: " << scut << endl ; + TCut Wcut = scut ; + + scut = "" ; + scut += degTomm(minLdeg) ; + scut += " < MHillas.fLength && MHillas.fLength < " ; + scut += degTomm(maxLdeg) ; + cout << "CUT ON LENGTH: " << scut << endl ; + TCut Lcut = scut ; + + scut = "" ; + scut += degTomm(minDdeg) ; + scut += " < MHillasSrc.fDist && MHillasSrc.fDist < " ; + scut += degTomm(maxDdeg) ; + cout << "CUT ON DIST: " << scut << endl ; + TCut Dcut = scut ; + + treeOn -> Draw ( "abs(MHillasSrc.fAlpha)>>absalphaOn" , Scut && Wcut && Lcut && Dcut ) ; + treeOff -> Draw ( "abs(MHillasSrc.fAlpha)>>absalphaOff" , Scut && Wcut && Lcut && Dcut ) ; + Double_t alpha = absalphaOn -> Integral ( 4 , 9 ) / absalphaOff -> Integral ( 4 , 9 ) ; + Double_t evoff = absalphaOff -> Integral ( 1 , 3 ) ; + + absalphaOff -> Scale ( alpha ) ; + + Double_t fakeMax = absalphaOn -> GetMaximum ( ) ; + fakeMax *= 1.1 ; + + Double_t fakeMin = absalphaOn -> GetMinimum ( ) ; + fakeMin *= 0.8 ; + + TH2D *fake = new TH2D ( "fake" , "Mrk421 Sample B (ON/OFF) 15/02/2004" , 100 , 0 , 90 , 100 , fakeMin , fakeMax ) ; + fake -> GetXaxis ( ) -> SetTitle ( "Parameter Alpha (deg)" ) ; + fake -> GetYaxis ( ) -> SetTitle ( "NUmber of Events" ) ; + fake -> GetYaxis ( ) -> SetTitleOffset ( 1.3 ) ; + fake -> SetStats ( 0 ) ; + + TCanvas *c = new TCanvas ( "c" , "Alpha Canvas" , 800 , 600 ) ; + c -> SetGridx ( ) ; + c -> SetGridy ( ) ; + + absalphaOn -> SetMarkerStyle ( 22 ) ; + absalphaOn -> SetMarkerSize ( 2.7 ) ; + absalphaOff -> SetFillStyle ( 1001 ) ; + absalphaOff -> SetFillColor ( 3 ) ; + + fake -> Draw ( ) ; + absalphaOff -> Draw ( "SAME" ) ; + absalphaOn -> Draw ( "SAMEE1" ) ; + + Double_t evon = absalphaOn -> Integral ( 1 , 3 ) ; + + cout << endl << "EntriesON:" << absalphaOn -> GetEntries ( ) ; + cout << endl << "EntriesOFF:" << absalphaOff -> GetEntries ( ) ; + cout << endl << "ON:" << evon << " OFF:" << evoff << " EXCESS:" << evon-evoff*alpha << " SIGMA:" << (evon-evoff*alpha)/TMath::Sqrt(evoff)/alpha << endl << endl ; + +} + Index: trunk/MagicSoft/Mars/mtemp/mpisa/macros/production.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpisa/macros/production.C (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mpisa/macros/production.C (revision 4055) @@ -0,0 +1,295 @@ +/* ======================================================================== *\ +! +! * +! * 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): Hendrik Bartko, 03/2004 +! Markus Gaug, 03/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////// +// +// production.C +// ============ +// +// This macro takes calibration runs and pedestal runs +// to set pedestals and calibration factors +// Then raw data files are read and calibrated. +// The output file contains the calibrated data, pedestal subtracted, +// with the signal (in photons) for each pixels +// +// input: calibration file(s) +// pedestal file(s) +// raw data files +// +// output: root files with following containers (branches) +// - MCerPhotEvt +// - MPedestalCam +// - MCalibrationRelTimeCam +// - MCerPhotEvt +// - MRawEvtHeader +// - MBadPixelsCam +// - MPedPhotCam +// +/////////////////////////////////////////////////////////////////////////// + +// +// Set path for files and run numbers +// +const TString defpath = "/data1/earth/magic/data/Period014/rootdata/2004_02_15/"; +const TString defrout = "output_Mrk421_20040215.root"; + +const Int_t defpedr [] = {17284}; +const Int_t defcalr [] = {17285}; +const Int_t defdatar[] = {17286,17287,17288,17289,17290,17291,17292,17293, + 17294,17295,17296,17297,17298,17299,17300,17301, + 17302,17303,17304,17305,17306,17307,17308,17309, + 17310,17311,17312,17313,17314,17315,17316,17317, + 17318,17319,17320,17321,17322,17323,17324,17325, + 17326,17327,17328,17329,17330,17331,17332,17333, + 17334,17335,17336,17337,17338,17339,17340,17341, + 17342,17343,17344,17345,17346,17347,17348,17349, + 17350,17351,17352,17353,17354,17355,17356,17357, + 17358,17359,17360,17361,17362,17363,17364,17365, + 17366,17367,17368,17369,17370,17371,17372,17373, + 17374,17375,17376,17377,17378,17379,17380,17381, + 17382,17383,17384,17385,17386,17387,17388,17389, + 17390,17391,17392,17393,17394,17395,17396,17397, + 17398,17399,17400,17401,17402,17403,17404,17405, + 17406,17407,17408,17409,17410,17411,17412,17413, + 17414,17415,17416,17417,17418,17419,17420,17421, + 17422,17423,17424,17425,17426,17427,17428,17429, + 17430,17431,17432,17433,17434,17435,17436,17437, + 17438,17439,17440,17441,17442,17443,17444,17445, + 17446,17447,17448,17449,17450,17451,17452,17453, + 17454,17455,17456}; + + +void production(const TString inpath=defpath, + const Int_t psize=1, const Int_t pedruns[]=defpedr, + const Int_t csize=1, const Int_t calruns[]=defcalr, + const Int_t dsize=171, const Int_t dataruns[]=defdatar, + const TString resname=defrout) + +{ + + MExtractSlidingWindow extractor; + + MRunIter pruns; + MRunIter cruns; + MRunIter druns; + + for (Int_t i=0;iSetUpdateTime(3000); + display->Resize(850,700); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + /************************************/ + /* FIRST LOOP: PEDESTAL COMPUTATION */ + /************************************/ + + MParList plist1; + MTaskList tlist1; + plist1.AddToList(&tlist1); + + // containers + MPedestalCam pedcam; + MBadPixelsCam badcam; + // + // for excluding pixels from the beginning: + // + // badcam.AsciiRead("badpixels.dat"); + + + plist1.AddToList(&pedcam); + plist1.AddToList(&badcam); + + //tasks + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + static_cast(read).AddFiles(pruns); + + MGeomApply geomapl; + MPedCalcPedRun pedcalc; + MGeomCamMagic geomcam; + + tlist1.AddToList(&read); + tlist1.AddToList(&geomapl); + tlist1.AddToList(&pedcalc); + + // Create and execute the event looper + MEvtLoop pedloop; + pedloop.SetParList(&plist1); + pedloop.SetDisplay(display); + + cout << "*************************" << endl; + cout << "** COMPUTING PEDESTALS **" << endl; + cout << "*************************" << endl; + + if (!pedloop.Eventloop()) + return; + + tlist1.PrintStatistics(); + + // + // Now the short version: + // + // + // Now setup the new tasks for the calibration: + // --------------------------------------------------- + // + MJCalibration calloop; + calloop.SetInput(&cruns); + // calloop.SetFullDisplay(); + // + calloop.SetExtractor(&extractor); + // + // Set the corr. cams: + // + calloop.SetBadPixels(badcam); + // + // The next two commands are for the display: + // + calloop.SetDisplay(display); + + // + // Apply rel. time calibration: + // + calloop.SetRelTimeCalibration(); + // + // Use as arrival time extractor MArrivalTimeCalc2: + // + calloop.SetArrivalTimeLevel(2); + + // + // Do the event-loop: + // + cout << "***************************" << endl; + cout << "** COMPUTING CALIBRATION **" << endl; + cout << "***************************" << endl; + + if (!calloop.Process(pedcam)) + return; + + badcam.Print(); + + MBadPixelsCam &badbad = calloop.GetBadPixels(); + MCalibrationChargeCam &calcam = calloop.GetCalibrationCam(); + MCalibrationRelTimeCam &timecam = calloop.GetRelTimeCam(); + MCalibrationQECam &qecam = calloop.GetQECam(); + + badbad.Print(); + + /************************************************************************/ + /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */ + /************************************************************************/ + + // Create an empty Parameter List and an empty Task List + MParList plist3; + MTaskList tlist3; + plist3.AddToList(&tlist3); + + // containers + MCerPhotEvt photevt; + MPedPhotCam pedphotcam; + MSrcPosCam srccam; + MRawRunHeader runhead; + MExtractedSignalCam sigcam; + + plist3.AddToList(&geomcam ); + plist3.AddToList(&pedcam ); + plist3.AddToList(&calcam ); + plist3.AddToList(&qecam ); + plist3.AddToList(&badbad ); + plist3.AddToList(&timecam ); + plist3.AddToList(&sigcam ); + plist3.AddToList(&photevt); + plist3.AddToList(&pedphotcam); + plist3.AddToList(&srccam); + plist3.AddToList(&runhead); + + //tasks + MReadMarsFile read3("Events"); + read3.DisableAutoScheme(); + static_cast(read3).AddFiles(druns); + + MArrivalTimeCalc2 timecalc; + MCalibrateData photcalc; + photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate + // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data + // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data + + tlist3.AddToList(&read3); + tlist3.AddToList(&geomapl); + tlist3.AddToList(&extractor); + tlist3.AddToList(&timecalc); + // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data + tlist3.AddToList(&photcalc); + // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data + + MWriteRootFile write(resname); + + write.AddContainer("MGeomCam" , "RunHeaders"); + write.AddContainer("MRawRunHeader" , "RunHeaders"); + write.AddContainer("MSrcPosCam" , "RunHeaders"); + write.AddContainer("MCalibrationChargeCam" , "RunHeaders"); + write.AddContainer("MCalibrationQECam" , "RunHeaders"); + // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!! + write.AddContainer("MPedestalCam" , "RunHeaders"); + write.AddContainer("MCalibrationRelTimeCam", "RunHeaders"); + + write.AddContainer("MCerPhotEvt" , "Events"); + write.AddContainer("MRawEvtHeader" , "Events"); + write.AddContainer("MBadPixelsCam" , "Events"); + write.AddContainer("MPedPhotCam" , "Events"); + + tlist3.AddToList(&write); + + // Create and execute eventloop + MEvtLoop evtloop3; + evtloop3.SetParList(&plist3); + + cout << "*************************************************************" << endl; + cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; + cout << "*************************************************************" << endl; + + if (!evtloop3.Eventloop()) + return; + tlist3.PrintStatistics(); + +} + + +