Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/callcleaning.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/callcleaning.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/callcleaning.C (revision 4088) @@ -0,0 +1,56 @@ +/* ======================================================================== *\ +! +! * +! * 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 Cherenkov 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): Daniel Mazin, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +// Daniel Mazin 14.05.2004 mazin@mppmu.mpg.de +// ********************************************************************************** +// this macro is used to make an image cleaning and calculate the Hillas parameters +// input: calibrated data (one or more files using wild card) +// output: hillas parameter +// ********************************************************************************** + +#include "cleaninghillas.C" +void callcleaning() +{ + gROOT->Reset(); + + TString typeInput ="ON"; // -> please specify "ON" or "OFF" + + const char *onfile="1*.ON"; + const char *offfile="1*.OFF"; + + + TString sourcename = "Crab"; + TString inPath; // directory where to find calibrated files + TString outPath; // directory where to write root-file with hillas parameters + + + inPath = "/.magic/magicserv01/scratch/David/CalibratedData/Crab/2004_01_27/"; + + outPath = "~mazin/data/Crab/2004_01_27/"; + + cleaninghillas(onfile, offfile, typeInput, inPath, outPath, sourcename); + + cout << " finished, back from callcleaning " << endl; + cout << " Hillas file " << outPath << sourcename << "Hillas" << typeInput << ".root is created " << endl; + +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/cleaninghillas.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/cleaninghillas.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/cleaninghillas.C (revision 4088) @@ -0,0 +1,230 @@ +/* ======================================================================== *\ +! +! * +! * 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 Cherenkov 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): Daniel Mazin, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +// Daniel Mazin 14.05.2004 mazin@mppmu.mpg.de +//********************************************************** +// this macro is a short form of the macro ONOFFAnalysis.C +// three jobs are executed: +// 1) image cleaning +// 2) calculation of hillas parameters +// 3) writing the root file +//********************************************************** + +gROOT->Reset(); + +//************************************************************************ +void cleaninghillas(const char *onfile, + const char *offfile, + TString typeInput, + TString inPath, + TString outPath, + TString sourcename) +{ + + TString fileON = inPath; + fileON += onfile; + fileON += ".root"; + + TString fileOFF = inPath; + fileOFF += offfile; + fileOFF += ".root"; + + // name of input root file + if (typeInput == "ON") + TString filenamein(fileON); + else if (typeInput == "OFF") + TString filenamein(fileOFF); + + // name of output root file + if (typeInput == "ON") + TString file(onfile); + else if (typeInput == "OFF") + TString file(offfile); + + TString outNameImage = outPath; + outNameImage += sourcename; + outNameImage += "Hillas"; + outNameImage += typeInput; + outNameImage += ".root"; + + + // set up the task list: + //----------------------------------------------------------- + MTaskList tliston; + MParList pliston; + + MObservatory observ; + MPointingPos mpos; + pliston.AddToList(&mpos); + + MMcEvt mmc; + pliston.AddToList(&mmc); + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenamein); + read.DisableAutoScheme(); + + MGeomApply apply; + + // NOTE: MMcEvt is used temporary. MPointingPos should be used as soon as synchronisation is working + // + // the following task is used to fill in a MMcEvt container with Theta and Phi values + // for each event using an ascii file (which is generated from central control report file) + // PLEASE NOTE: only needed if a derotation around the center of the camera is planned + // to be performed afterwords + MSourcePosfromStarPos sourcefromstar; + +// if (typeInput == "ON") +// { + sourcefromstar.AddFile("~mazin/mars/MarsCVS/Mars/data/positions2103.txt", 0); +// } +// else if (typeInput == "OFF") +// { +// sourcefromstar.AddFile("~mazin/data/dataMrk421_15feb/positionsFeb04_0", 0); +// } + + + +// ----------------------------------------------------- +// image cleaning ------------------------------------- +// ----------------------------------------------------- + MSigmabarCalc sigbar; + + MBlindPixelsCalc2 blind; + blind.SetUseInterpolation(); + blind.SetCheckPedestalRms(); + blind.SetName("BlindAfterPadding"); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + +// MImgCleanStd clean(3.0, 2.5); // Crab 27th old calibr (Hendrik) + MImgCleanStd clean(4.0, 3.0); +// clean.SetMethod(MImgCleanStd::kDemocratic); + clean.SetCleanRings(1); + + +// ----------------------------------------------------- +// calculation of image parameters --------------------- +// ----------------------------------------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + MHillasCalc hcalc; + hcalc.SetNameHillas(fHilName); + hcalc.SetNameHillasExt(fHilNameExt); + hcalc.SetNameNewImgPar(fImgParName); + + MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc); + hsrccalc.SetInput(fHilName); + + MFSelStandard selstandard(fHilNameSrc); + selstandard.SetHillasName(fHilName); + selstandard.SetImgParName(fImgParName); +// selstandard.SetCuts(usedpixelsmax, corepixelsmin, sizemin, distmin, distmax, lengthmin, widthmin); + selstandard.SetCuts(200, 5, 200, 0.0, 5.0, 0.07, 0.05); + MContinue contstandard(&selstandard); + contstandard.SetName("SelStandard"); + +// ----------------------------------------------------- +// write root file +// ----------------------------------------------------- + MWriteRootFile write(outNameImage); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MTime", "Events"); + write.AddContainer("MPointingPos", "Events"); + //write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + pliston.AddToList(&observ); + + pliston.AddToList(&source); + + //***************************** + // entries in MTaskList + + + tliston.AddToList(&read); + tliston.AddToList(&apply); + tliston.AddToList(&sourcefromstar); + + + tliston.AddToList(&sigbar); + tliston.AddToList(&blind); + tliston.AddToList(&contbasic); + + tliston.AddToList(&clean); + + tliston.AddToList(&hcalc); + tliston.AddToList(&hsrccalc); + tliston.AddToList(&contstandard); + tliston.AddToList(&write); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + //evtloop.ReadEnv(env, "", printEnv); + evtloop.SetProgressBar(&bar); + // evtloop.Write(); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + +} + + + Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/falsesourcemethod.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/falsesourcemethod.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/falsesourcemethod.C (revision 4088) @@ -0,0 +1,171 @@ +/* ======================================================================== *\ +! +! * +! * 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 Cherenkov 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): Daniel Mazin, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + + +// Daniel Mazin 14.05.2004 mazin@mppmu.mpg.de +// ********************************************************************************** +// this macro is used to produce false source plots +// either using one sample only (usually used for OFF, set TYPEOPTION = kFALSE) +// or in addition using the background estimation from OFF (set TYPEOPTION = kTRUE). +// The later option can be used ONLY if false source method for the OFF data +// has taken place before. +// input: hillas parameter file +// output: several 2D plots in a root file + +// ascii file with parameters like significance, Nexcess etc. +// ********************************************************************************** + + +// you may want to change these values: +#define STEPGRID 0.05 // [deg] +#define MINXGRID -1. // [deg] +#define MINYGRID -1. // [deg] +#define MAXXGRID 1. // [deg] +#define MAXYGRID 1. // [deg] + +#define ALPHAMAX 10. // [deg], ON region in the alpha plot, maybe 5 deg is better + // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible + +// if you wish to rotate set to kTRUE, else kFALSE +const Bool_t ROTOPTION = kFALSE; +const Bool_t TYPEOPTION = kFALSE; // kFALSE = estimate background from the same sample + // kTRUE = estimate background also from OFF data +const Bool_t CUTOPTION = kTRUE; // kTRUE = dynamical cuts + // kFALSE = statical cuts + +/* ****************************************************** */ + /* static cuts */ +#define LENGTHMIN 0.10 // deg +#define LENGTHMAX 0.32 // deg +#define WIDTHMIN 0.06 // deg +#define WIDTHMAX 0.15 // deg +#define SIZEMIN 2000. // photons +#define DISTMIN 0.6 // deg +#define DISTMAX 1.3 // deg +#define LEAKMAX 0.25 // deg +/* ****************************************************** */ + +/* ****************************************************** */ + /* dynamical cuts Crab 27th Jan 2004 */ +#define LENGTHMINParA 0.136 // deg +#define LENGTHMINParB 0.036 // +#define LENGTHMINParC -0.0038 // +#define LENGTHMAXParA 0.332 // deg +#define LENGTHMAXParB 0.037 // +#define LENGTHMAXParC 0.0261 // +#define WIDTHMINParA 0.063 // deg +#define WIDTHMINParB 0.013 // +#define WIDTHMINParC 0.0003 // +#define WIDTHMAXParA 0.123 // deg +#define WIDTHMAXParB 0.019 // +#define WIDTHMAXParC 0.0005 // +#define DISTMINParA 0.6 // deg +#define DISTMINParB 0.059 // +#define DISTMINParC 0. // +#define DISTMAXParA 1.25 // deg +#define DISTMAXParB 0.059 // +#define DISTMAXParC 0. // +/* ****************************************************** */ + + /* dynamical cuts Mrk 421 and Crab 15th Feb 2004 */ + +/* +#define LENGTHMINParA 0.12 // deg +#define LENGTHMINParB 0.034 // +#define LENGTHMINParC 0. // +#define LENGTHMAXParA 0.32 // deg +#define LENGTHMAXParB 0.034 // +#define LENGTHMAXParC 0. // +#define WIDTHMINParA 0.055 // deg +#define WIDTHMINParB 0.013 // +#define WIDTHMINParC 0.0 // +#define WIDTHMAXParA 0.12 // deg +#define WIDTHMAXParB 0.013 // +#define WIDTHMAXParC 0.0 // +#define DISTMINParA 0.6 // deg +#define DISTMINParB 0.059 // +#define DISTMINParC 0. // +#define DISTMAXParA 1.25 // deg +#define DISTMAXParB 0.059 // +#define DISTMAXParC 0. // +*/ +/* ****************************************************** */ + +// please do not change these values: +#define NUMSTEPS (int) ((MAXXGRID - MINXGRID) / STEPGRID + 1.5) +#define GRIDBINS (int) (NUMSTEPS * NUMSTEPS + 0.5) + + +#include "mtools.C" +#include // needed histograms + +// this helps me to check if the binning is ok +#define TOLERANCE 1e-3 + +#include "hillasfileread.C" +#include "gridloop.C" + +void falsesourcemethod() +{ + + if (TYPEOPTION == kTRUE) // output of OFF data sample +// DATPARAMIN = "data/ressignOffMrk3004.dat"; // 3000 + DATPARAMIN = "data/ressignOffCrab1502_1005.dat"; // 2000 + + + HILLASFILENAME = "~/data/Crab/2004_01_27/CrabHillasOFF.root"; + DATPARAMOUT = "data/paramOffCrab2701.dat"; + ROOTPLOTNAME = "/.magic/magicserv01/scratch/Daniel/plots/RootPlots/Crab/2004_01_27/CrabOffSkyPlots.root"; + hillasfileread(); + gridloop(); + +/* +// use this loop for different subsamples + for(Int_t i=1; i<10; i++) + { + HILLASFILENAME = "~/data/Crab15022004Update/sample"; + DATPARAMOUT = "data/ressignOnOffCrab1502Samp"; + ROOTPLOTNAME = "/.magic/magicserv01/scratch/Daniel/plots/RootPlots/Crab15022004Update/Crab1502Sample"; + + HILLASFILENAME += i; + HILLASFILENAME += "/"; + HILLASFILENAME += "CrabHillasON"; + HILLASFILENAME += ".root"; + + DATPARAMOUT += i; + DATPARAMOUT += "_AlphaCut"; + DATPARAMOUT += TMath::Nint(ALPHAMAX); + DATPARAMOUT += "deg"; + DATPARAMOUT += ".dat"; + + ROOTPLOTNAME += i; + ROOTPLOTNAME += "_AlphaCut"; + ROOTPLOTNAME += TMath::Nint(ALPHAMAX); + ROOTPLOTNAME += "deg"; + ROOTPLOTNAME += ".root"; + + + hillasfileread(); + gridloop(); + } +*/ +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/gridloop.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/gridloop.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/gridloop.C (revision 4088) @@ -0,0 +1,432 @@ +/* ======================================================================== *\ +! +! * +! * 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 Cherenkov 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): Daniel Mazin, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +// Daniel Mazin 17.05.2004 mazin@mppmu.mpg.de +// ********************************************************************************** +// this macro is used to produce false source plots +// and it is a part from the macro falsesourcemethod.C. +// See falsesourcemethod.C for more explanation +// input: hillas parameters as vectors +// output: several 2D plots in a root file + +// ascii file with parameters like significance, Nexcess etc. +// ********************************************************************************** + + +void gridloop() +{ + + printf("inside of gridloop() ... \n"); + +/* emty histograms */ + + histskyplot.Reset(); + histNexOnOff.Reset(); + histChi2Off.Reset(); + histskyplotOnOff.Reset(); + histNex.Reset(); + histChi2.Reset(); + histskyplot4.Reset(); + histskyLiMa.Reset(); + histsign.Reset(); + histsign4.Reset(); + histBerlin.Reset(); + +/* end emty histograms */ + + + TH1F histalphagrid; + histalphagrid.SetName("Alpha"); + histalphagrid.SetTitle("Alpha"); + histalphagrid.SetXTitle("Alpha (deg)"); + histalphagrid.SetYTitle("Counts"); + histalphagrid.SetDirectory(NULL); + histalphagrid.SetFillStyle(4000); + histalphagrid.UseCurrentStyle(); + + bins.SetEdges(18, 0.0, 90.); + bins.Apply(histalphagrid); + + Double_t xcorr, ycorr, xtest, ytest, xtilde, ytilde, fconc1; + + const Float_t binwidth1 = histalphagrid.GetBinWidth(1); + if (TMath::Abs( ALPHAMAX/binwidth1 - TMath::Nint(ALPHAMAX/binwidth1)) > TOLERANCE) + { + cout << "wrong binning for histogram: histalphagrid \n SORRY ... ABORT" << endl; + exit(1); + } + const Int_t numbbinOnMax = TMath::Nint(ALPHAMAX / binwidth1); // number bins in ON region + + /* fit functions */ + TF1 * fitsig = new TF1("fsig", "[0]", 0., ALPHAMAX); + fitsig->SetLineColor(4); + TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", 30., 90.); + fitbgpar->SetLineColor(2); + + TF1 * fitbg4 = new TF1("fbg4", "[0]*x*x*x*x + [1]*x*x + [2]", 30., 90.); + fitbg4->SetLineColor(3); + /* end fit functions */ + + Float_t Nex, Non, Noff, Sign; + Float_t xsource, ysource, xsournew, ysournew; + + Double_t Nex4, Noff4, Sign4; + Double_t alpha, SignLiMa; + + Double_t chisquarefit; + + Double_t integon, integoff, normf, NoffOFF, SignOnOff, NexOnOff; + Float_t logsize, lgsize, lgsize2, tanbeta, beta; + const Float_t LOG3000 = log(3000.); + + + + Double_t chi2par[GRIDBINS], apar[GRIDBINS], bpar[GRIDBINS]; // fit parameters from OFF sample + Float_t offbin1[GRIDBINS], offbin2[GRIDBINS], offbin3[GRIDBINS]; // Number of events in first 3 bins + Double_t xpos, ypos, binwidth2; + Int_t numgrid; + + + Int_t numbinsoff = 12; // FIXME: 12 should be changed according to the binning + Int_t nbinx, nbiny, gridpoint = 0; + + Char_t psname[120], titelname[120], plot1name[120], plot2name[120], stringsig[120]; + + FILE * fp; + + if (TYPEOPTION == kTRUE) // read fit parameters from the OFF sample + { + Char_t * datin = DATPARAMIN; + + fp = fopen(datin, "r"); + + for (Int_t i=0; i 90.) alphapar(k) -= 180.; + if(alphapar(k) < -90.) alphapar(k) += 180.; + +// histalpha.Fill(alphapar(k),1.); +// histdist.Fill( distpar(k),1.); + + if(CUTOPTION==kFALSE) + { +/* ****************************************************** */ + /* static cuts */ + if (lengthpar(k) > LENGTHMIN && lengthpar(k) < LENGTHMAX) + if (widthpar(k) > WIDTHMIN && widthpar(k) < WIDTHMAX) + if (sizepar(k) > SIZEMIN) + if (distpar(k) > DISTMIN && distpar(k) < DISTMAX) + if (leak1(k) < LEAKMAX) +/* ****************************************************** */ + { + alphapar(k) = TMath::Abs(alphapar(k)); +//printf("alpha passed the cuts: %5.2f \n", alphapar(k)); +//printf("histalphagrid: %5.2f \n", histalphagrid.GetNbinsX()); + histalphagrid.Fill(alphapar(k),1.); + + } + + } +/* ****************************************************** */ +/* dynamical cuts */ + else // if(CUTOPTION==kTRUE) + { + logsize = log(sizepar(k)); + lgsize = logsize-LOG3000; + lgsize2 = lgsize*lgsize; + if ( sizepar(k) > SIZEMIN ) + if ( leak1(k) < LEAKMAX ) + if ( lengthpar(k) > (LENGTHMINParA + LENGTHMINParB*lgsize + LENGTHMINParC*lgsize2) && + lengthpar(k) < (LENGTHMAXParA + LENGTHMAXParB*lgsize + LENGTHMAXParC*lgsize2)) + if ( widthpar(k) > (WIDTHMINParA + WIDTHMINParB*lgsize + WIDTHMINParC*lgsize2) && + widthpar(k) < (WIDTHMAXParA + WIDTHMAXParB*lgsize + WIDTHMAXParC*lgsize2) ) + if ( distpar(k) > (DISTMINParA + DISTMINParB*lgsize + DISTMINParC*lgsize2) && + distpar(k) < (DISTMAXParA + DISTMAXParB*lgsize + DISTMAXParC*lgsize2) ) + { + alphapar(k) = TMath::Abs(alphapar(k)); + histalphagrid.Fill(alphapar(k),1.); + if(alphapar(k) < 8.) histBerlin.Fill(xsournew, ysournew); + } + } + + } +// histalphagrid.Fit("fsig","NR"); +//printf("content bin 1: %4.2f\n", histalphagrid.GetBinContent(1)); + histalphagrid.Fit(fitsig,"NR"); + histalphagrid.Fit("fbg4","+NR"); + histalphagrid.Fit("fbgpar","+NR"); + +// cout << "OK" << endl; +// calculate significance: + Non = 0.; + for(Int_t i=1; i<=numbbinOnMax;i++) Non += histalphagrid.GetBinContent(i); +// Non = (fitsig->GetParameter(0)) * numbbinOnMax; + +// histalphagrid.Draw(); + + chisquarefit = fitbgpar->GetChisquare(); + Noff = (1./3. * (fitbgpar->GetParameter(0)) * pow(ALPHAMAX,3.) + + (fitbgpar->GetParameter(1)) * ALPHAMAX) / binwidth1; + + Nex = Non - Noff; + + Noff4 = (1./5. * fitbg4->GetParameter(0) * pow(ALPHAMAX,5.) + + 1./3. * (fitbg4->GetParameter(1)) * pow(ALPHAMAX,3.) + + fitbg4->GetParameter(2) * ALPHAMAX ) / binwidth1; + + + Nex4 = Non - Noff4; + +// calculate significance: + +// Sign4 = Nex4 / sqrt(Nex4 + 2.* Noff4); + Sign4 = LiMa17(Non,Noff4,1.); +// Sign = Nex / sqrt(Nex + 2.* Noff); + Sign = LiMa17(Non,Noff,1.); + + /* calculate Noff from the OFF data */ + if (TYPEOPTION == kTRUE) + { + // ON: + + integon = 0.; // number of events between 30 and 90 degrees + + for (Int_t ik = 0; ik < numbinsoff; ik++) + { + integon += histalphagrid.GetBinContent(7+ik); // FIXME: 7 should be changed according to the binning + } + // OFF: + + integoff = ((1./3. * apar[gridpoint] * pow(90.,3.) + bpar[gridpoint] * 90.) - + (1./3. * apar[gridpoint] * pow(30.,3.) + bpar[gridpoint] * 30.)) / binwidth2; + + normf = integoff / integon; + + NoffOFF = (1./3. * apar[gridpoint] * pow(ALPHAMAX,3.) + + (bpar[gridpoint] * ALPHAMAX)) / binwidth2 / normf; + + NexOnOff = Non - NoffOFF; + + SignOnOff = NexOnOff / sqrt(NexOnOff + 2.* NoffOFF); + + // calculate according to Li Ma: + SignLiMa = LiMa17(Non,NoffOFF*normf,1./normf); + + } + + cout << endl << " xsource: " << xsource << " ysource: " << ysource << endl; + + cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl; + + if (TYPEOPTION == kTRUE) + cout << " ON and OFF : Non : " << Non << " Noff (estimated) : " << NoffOFF << " Nex : " << NexOnOff << endl; + + cout << endl; + cout << " significance (parab): " << Sign << " sigma" << endl; + cout << " significance (4th order): " << Sign4 << " sigma" << endl; + if (TYPEOPTION == kTRUE) + { + cout << " significance (using ON and OFF): " << SignOnOff << " sigma" << endl; + cout << " significance (LiMa 17): " << SignLiMa << " sigma" << endl; + } + + + if (TYPEOPTION == kTRUE) + fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf \n", + xsource, ysource, Sign, Non, Noff, Nex, SignLiMa, + SignOnOff, NoffOFF, NexOnOff, chi2par[gridpoint]); + else + fprintf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", + gridpoint+1, xsource, ysource, + fitbgpar->GetParameter(0), fitbgpar->GetParameter(1), + fitbgpar->GetChisquare(), histalphagrid.GetBinWidth(1), + histalphagrid.GetBinContent(1), histalphagrid.GetBinContent(2), + histalphagrid.GetBinContent(3)); + + + histskyplot.SetBinContent(gridx+1, gridy+1, Sign); + histNex.SetBinContent(gridx+1, gridy+1, Nex); + histChi2.SetBinContent(gridx+1, gridy+1, fitbgpar->GetChisquare() / 10.); + + if (TYPEOPTION == kTRUE) histskyplotOnOff.SetBinContent(gridx+1, gridy+1, SignOnOff); + if (TYPEOPTION == kTRUE) histNexOnOff.SetBinContent(gridx+1, gridy+1, NexOnOff); + if (TYPEOPTION == kTRUE) histChi2Off.SetBinContent(gridx+1, gridy+1, chi2par[gridpoint] / 16.); + + if (TYPEOPTION == kTRUE) histskyLiMa.SetBinContent(gridx+1, gridy+1, SignLiMa); + histskyplot4.SetBinContent(gridx+1, gridy+1, Sign4); + histsign.Fill(Sign,1.); + histsign4.Fill(Sign4,1.); +// break; + + if (TYPEOPTION == kTRUE) + { + TF1 * bgoff = new TF1("bgoff", parabfunc, 0., 90., 3); + bgoff->SetParameters(apar[gridpoint], bpar[gridpoint], normf); + bgoff->FixParameter(0, apar[gridpoint]); + bgoff->FixParameter(1, bpar[gridpoint]); + bgoff->FixParameter(2, normf); + bgoff->SetLineColor(8); + } + +// sprintf(stringsig,"S = %.2f sigma", Sign); + if (TYPEOPTION == kTRUE) sprintf(stringsig,"S = %.2f sigma", SignOnOff); +/* + sprintf(psname,"/.magic/magicserv01/scratch/Daniel/plots/RootPlots/040215Mrk421_misp/all/alphaX%.2fY%.2fON.root", xsource, ysource); + sprintf(titelname,"alpha plot. assumed source position: x = %.2f y = %.2f", xsource, ysource); + TCanvas * cplot = new TCanvas(); + histalphagrid.SetTitle(titelname); + histalphagrid.Draw(); + bgoff->Draw("same"); + leg = new TLegend(0.1,0.15,0.52,0.35); + leg->Draw(); + leg->AddEntry(fitsig,"fit for ON region (0-20 deg)","l"); + leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l"); + leg->AddEntry(bgoff,"fit from OFF data (0-90 deg)","l"); + leg->SetHeader("Legend"); + leg->SetFillColor(19); + leg->Draw(); + + text = new TPaveText(50,20,90,80); + text->AddText(0.5, 0.8, "calculated significance"); + text->AddText(0.5, 0.5, "using ON and OFF:"); + text->AddText(1, 0.2, stringsig); + text->SetTextSize(0.04); + text->Draw(); + + cplot->Modified(); + cplot->Update(); + + cplot->SaveAs(psname); + delete cplot; +*/ + if (TYPEOPTION == kTRUE) delete bgoff; + + gridpoint++; + + } + } + + /* close file */ + fclose(fp); + + + TFile rootfile(ROOTPLOTNAME, "RECREATE", + "sky plots in all variations"); + + +// cout << " end of one of the grids 1 " << endl; + + histskyplot.Write(); + + histNex.Write(); + + histChi2.Write(); + + histsign.Write(); + + histskyplot4.Write(); + + histsign4.Write(); + + histBerlin.Write(); + + if (TYPEOPTION == kTRUE) + { + histNexOnOff.Write(); + histChi2Off.Write(); + histskyplotOnOff.Write(); + histskyLiMa.Write(); + } + + rootfile.Close(); + + delete fitsig; + delete fitbgpar; + delete fitbg4; + + cout << " end of one of the grids " << endl; +} + Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/hillasfileread.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/hillasfileread.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/hillasfileread.C (revision 4088) @@ -0,0 +1,178 @@ +/* ======================================================================== *\ +! +! * +! * 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 Cherenkov 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): Daniel Mazin, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +// Daniel Mazin 17.05.2004 mazin@mppmu.mpg.de +// ********************************************************************************** +// this macro is a part from the macro falsesourcemethod.C +// This macro reads mars file with hillas parameters etc. and puts the values into vectors +// input: hillas parameter file +// output: vectors with the same information +// ********************************************************************************** + + +void hillasfileread() +{ + lengthpar.ResizeTo(MAXNUMEVENTS); + widthpar.ResizeTo(MAXNUMEVENTS); + sizepar.ResizeTo(MAXNUMEVENTS); + xpar.ResizeTo(MAXNUMEVENTS); + ypar.ResizeTo(MAXNUMEVENTS); + deltapar.ResizeTo(MAXNUMEVENTS); + alphapar.ResizeTo(MAXNUMEVENTS); + distpar.ResizeTo(MAXNUMEVENTS); + cosgam.ResizeTo(MAXNUMEVENTS); + singam.ResizeTo(MAXNUMEVENTS); + conc1.ResizeTo(MAXNUMEVENTS); + leak1.ResizeTo(MAXNUMEVENTS); + + lengthpar.Zero(); + widthpar.Zero(); + sizepar.Zero(); + xpar.Zero(); + ypar.Zero(); + deltapar.Zero(); + alphapar.Zero(); + distpar.Zero(); + cosgam.Zero(); + singam.Zero(); + conc1.Zero(); + leak1.Zero(); + + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + + read.AddFile(HILLASFILENAME); + + MHillas mhillas; + plist.AddToList(&mhillas); + + MHillasSrc mhillassrc; + plist.AddToList(&mhillassrc); + + MNewImagePar mnewimpar; + plist.AddToList(&mnewimpar); + + MGeomCamMagic cam; + plist.AddToList(&cam); + + MMcEvt mcevt; + plist.AddToList(&mcevt); + + MObservatory observ; + plist.AddToList(&observ); + +// printf("observ: %f\n",observ); + + MHStarMap starmap; +// plist.AddToList(&starmap); + starmap.ReInit(&plist); +// MRawRunHeader header; +// plist.AddToList(&header); + + tlist.AddToList(&read); + + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + Float_t fMm2Deg; + +cout << " file name: " << HILLASFILENAME << endl; + + if (!tlist.PreProcess(&plist)) + return; + + fMm2Deg = cam->GetConvMm2Deg(); + Int_t event = 0; + Int_t filenumber = 0; + Float_t theta, phi; + + imnum = 0; + + while (tlist.Process()) + { + event++; +//cout << event << endl; + + if (mhillas->GetLength() != -1.) + { + + + if (ROTOPTION == kTRUE) + { + theta = mcevt->GetTelescopeTheta(); + phi = mcevt->GetTelescopePhi(); +// printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415); + observ.RotationAngle(theta, phi, singam(imnum), cosgam(imnum)); + } + + lengthpar(imnum) = (mhillas->GetLength()) * fMm2Deg; + widthpar(imnum) = (mhillas->GetWidth())*fMm2Deg; + sizepar(imnum) = mhillas->GetSize(); + xpar(imnum) = (mhillas->GetMeanX())*fMm2Deg; + ypar(imnum) = (mhillas->GetMeanY())*fMm2Deg; + deltapar(imnum) = mhillas->GetDelta(); + conc1(imnum) = (mnewimpar->GetConc1()); + leak1(imnum) = mnewimpar->GetLeakage1(); + + imnum++; + + } + else filenumber++; + } + + if (!tlist.PostProcess()) + return; + + printf("imnum = %d\n",imnum); + + lengthpar.ResizeTo(imnum); + widthpar.ResizeTo(imnum); + sizepar.ResizeTo(imnum); + xpar.ResizeTo(imnum); + ypar.ResizeTo(imnum); + deltapar.ResizeTo(imnum); + alphapar.ResizeTo(imnum); + distpar.ResizeTo(imnum); + cosgam.ResizeTo(imnum); + singam.ResizeTo(imnum); + conc1.ResizeTo(imnum); + leak1.ResizeTo(imnum); + + + cout << " conversion factor is: " << fMm2Deg << endl; + cout << " events read in from file : " << event << endl; + cout << " runs found in the file : " << filenumber << endl; + +} + Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/mtools.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/mtools.C (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/mtools.C (revision 4088) @@ -0,0 +1,19 @@ +Double_t LiMa17(Double_t Non, Double_t Noff, Double_t alpha) +{ + Double_t term1, term2, signif; + term1 = ((1.+alpha) / alpha) * (Non / (Non + Noff)); + term2 = (1.+alpha) * (Noff / (Non + Noff)); + + signif = sqrt(2.* (Non*log(term1) + Noff*log(term2))); + + if (Non>(alpha*Noff)) return signif; + else return -signif; + +} + +Double_t parabfunc(Double_t *x, Double_t *par) +{ + Double_t result=0.; + if (par[2]!=0) result = (par[0] * x[0] * x[0] + par[1]) / par[2]; + return result; +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/skyplot.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/skyplot.h (revision 4088) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/skyplot.h (revision 4088) @@ -0,0 +1,201 @@ +// temporal solution: +#define MAXNUMEVENTS 500000 +#define LOG3000 log(3000.) + +TString HILLASFILENAME; +TString DATPARAMIN; +TString DATPARAMOUT; +TString ROOTPLOTNAME; +Int_t imnum = 0; + +TVectorD lengthpar(MAXNUMEVENTS); +TVectorD widthpar(MAXNUMEVENTS); +TVectorD distpar(MAXNUMEVENTS); +TVectorD alphapar(MAXNUMEVENTS); +TVectorD sizepar(MAXNUMEVENTS); +TVectorD deltapar(MAXNUMEVENTS); +TVectorD xpar(MAXNUMEVENTS); +TVectorD ypar(MAXNUMEVENTS); +TVectorD conc1(MAXNUMEVENTS); +TVectorD leak1(MAXNUMEVENTS); +TVectorD cosgam(MAXNUMEVENTS); +TVectorD singam(MAXNUMEVENTS); + +// the alpha histogram +MBinning bins; + +//TH1F histalphagrid; +//histalphagrid.SetName("Alpha"); +//histalphagrid.SetTitle("Alpha"); +//histalphagrid.SetXTitle("Alpha [deg]"); +//histalphagrid.SetYTitle("Counts"); +//histalphagrid.SetDirectory(NULL); +//histalphagrid.SetFillStyle(4000); +//histalphagrid.UseCurrentStyle(); + +bins.SetEdges(18, 0.0, 90.); +bins.Apply(histalphagrid); + + + +// possible histograms: +TH1F histlength; +TH1F histwidth; +TH1F histsize; +TH1F histalpha; +TH1F histdist; +TH1F histmeanx; +TH1F histmeany; +TH1F histalphafinal; + +TH2F hist2xy("CoG1","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8); +TH1F histLoverS; +TH1F aftercuthistlength; +TH1F aftercuthistwidth; +TH1F aftercuthistsize; +TH1F aftercuthistalpha; +TH1F aftercuthistdist; +TH1F aftercuthistmeanx; +TH1F aftercuthistmeany; +TH2F aftercuthist2xy("CoG2","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8); +TH2F hist2lw("Length-Width", "correlation Length-Width", 100, 0.0, 1.0, 100, 0.0, 0.5); +TH2F hist2lalpha("Length-Alpha", "correlation Length-Alpha", 100, 0.0, 1.0, 100, -100., 100.); +TH2F hist2ldist("Length-Dist","correlation Length-Dist", 100, 0.0, 1.0, 100, 0.0, 1.6); +TH2F hist2walpha("Width-Alpha","correlation Width-Alpha", 100, 0.0, 0.5, 100, -100., 100.); +TH2F hist2wdist("Width-Dist","correlation Width-Dist", 100, 0.0, 0.5, 100, 0.0, 1.6); +TH2F hist2alphadist("Alpha-Dist","correlation Alpha-Dist", 100, -100., 100, 100, 0.2, 1.4); + + +TH2F histskyplot("Sky Plot 1","Sky Plot of significance (2nd order fit)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histskyplot.SetXTitle("-RA [deg]"); +//histskyplot.SetYTitle("DEC [deg]"); +histskyplot.SetDirectory(NULL); +histskyplot.SetFillStyle(4000); +histskyplot.UseCurrentStyle(); + +TH2F histNexOnOff("Sky Plot 2","Sky Plot of Number of excess events (using ON and OFF data)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histNexOnOff.SetXTitle("-RA [deg]"); +//histNexOnOff.SetYTitle("DEC [deg]"); +histNexOnOff.SetDirectory(NULL); +histNexOnOff.SetFillStyle(4000); +histNexOnOff.UseCurrentStyle(); + +TH2F histChi2Off("Sky Plot 3","Sky Plot of Chi2/ndf (2nd order fit, OFF data)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histChi2Off.SetXTitle("-RA [deg]"); +//histChi2Off.SetYTitle("DEC [deg]"); +histChi2Off.SetDirectory(NULL); +histChi2Off.SetFillStyle(4000); +histChi2Off.UseCurrentStyle(); + +TH2F histskyplotOnOff("Sky Plot 4","Sky Plot of significance (using ON and OFF data)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histskyplotOnOff.SetXTitle("-RA [deg]"); +//histskyplotOnOff.SetYTitle("DEC [deg]"); +histskyplotOnOff.SetDirectory(NULL); +histskyplotOnOff.SetFillStyle(4000); +histskyplotOnOff.UseCurrentStyle(); + +TH2F histNex("Sky Plot 5","Sky Plot of Number of excess events (2nd order fit)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histNex.SetXTitle("-RA [deg]"); +//histNex.SetYTitle("DEC [deg]"); +histNex.SetDirectory(NULL); +histNex.SetFillStyle(4000); +histNex.UseCurrentStyle(); + +TH2F histChi2("Sky Plot 6","Sky Plot of Chi2/ndf (2nd order fit)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histChi2.SetXTitle("-RA [deg]"); +//histChi2.SetYTitle("DEC [deg]"); +histChi2.SetDirectory(NULL); +histChi2.SetFillStyle(4000); +histChi2.UseCurrentStyle(); + +TH2F histskyplot4("Sky Plot 7","Sky Plot of significance (4th order fit)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histskyplot4.SetXTitle("-RA [deg]"); +//histskyplot4.SetYTitle("DEC [deg]"); +histskyplot4.SetDirectory(NULL); +histskyplot4.SetFillStyle(4000); +histskyplot4.UseCurrentStyle(); + +TH2F histskyLiMa("Sky Plot 9","Sky Plot of significance (ON - OFF LiMa17)", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histskyLiMa.SetXTitle("-RA [deg]"); +//histskyLiMa.SetYTitle("DEC [deg]"); +histskyLiMa.SetDirectory(NULL); +histskyLiMa.SetFillStyle(4000); +histskyLiMa.UseCurrentStyle(); + +TH2F histBerlin("Sky Plot Berlin","Sky Plot of Events below alpha 8 degrees", + NUMSTEPS, MINXGRID - 0.5*STEPGRID, MAXXGRID + 0.5*STEPGRID, + NUMSTEPS, MINYGRID - 0.5*STEPGRID, MAXYGRID + 0.5*STEPGRID); +//histskyLiMa.SetXTitle("-RA [deg]"); +//histskyLiMa.SetYTitle("DEC [deg]"); +histBerlin.SetDirectory(NULL); +histBerlin.SetFillStyle(4000); +histBerlin.UseCurrentStyle(); + +TH1F histsign; +histsign.SetName("Significance 1"); +histsign.SetTitle("Significance (2nd order fit)"); +histsign.SetXTitle("significance"); +histsign.SetYTitle("Counts"); +histsign.SetDirectory(NULL); +histsign.SetFillStyle(4000); +histsign.UseCurrentStyle(); + + bins.SetEdges(20, -5., 8.); + bins.Apply(histsign); + +TH1F histsign4; +histsign4.SetName("Significance 2"); +histsign4.SetTitle("Significance (4th order fit)"); +histsign4.SetXTitle("significance"); +histsign4.SetYTitle("Counts"); +histsign4.SetDirectory(NULL); +histsign4.SetFillStyle(4000); +histsign4.UseCurrentStyle(); + + bins.SetEdges(20, -5., 8.); + bins.Apply(histsign4); + + +TH1F histalphagrid; +histalphagrid.SetName("Alpha"); +histalphagrid.SetTitle("Alpha"); +//histalphagrid.SetXTitle("Alpha [deg]"); +histalphagrid.SetYTitle("Counts"); +histalphagrid.SetDirectory(NULL); +histalphagrid.SetFillStyle(4000); +histalphagrid.UseCurrentStyle(); + + bins.SetEdges(18, 0.0, 90.); + bins.Apply(histalphagrid); + +/* +#define n 100 +Double_t binsize[n]; +Float_t nmin = 100.; +Float_t nmax = 1e7; + +for(Int_t i=0; i